Ghostscript wrapper for D:\Digitalizacja\MTS80_t18z1_4_PDF_artyku³y\mts80_t18z2.pdf M E C H AN I K A TEORETYCZNA I  STOSOWANA 2,  18 (1980) N IELIN IOWA,  STATYCZN A  I  D YN AM ICZN A  AN ALIZA  P OWŁ OK  M E T O D Ą ELEM EN TÓW SKOŃ C Z ON YC H1* M ICH AŁ  K  L  E I B  E R  (WARSZAWA) 1,  Uwagi wstę pne W  wyniku  znacznego  w  ostatnim  dziesię cioleciu  rozwoju  metody  elementów  skoń- czonych  stał o  się   obecnie  moż liwe  uzyskiwanie  zadawalają cych  projektanta  (tj.  dosta- tecznie  dokł adnych)  rozwią zań  nawet  dla  trudnych  zagadnień  brzegowych  w  zakresie liniowej  statyki  powł ok.  Postę pują ce  udoskonalenia  odpowiednich,  programów  oblicze- niowych  uczynił y  z  rozwią zywania  prawie  wszystkich  zagadnień  dotyczą cych  liniowej teorii  powł ok  standardową   czynnoś ć,  stosunkową   .prostą   do  bezpoś redniego  i  natych- miastowego  wykonania.  Pomimo,  iż  dla  pewnej  klasy  zagadnień  (jak  n p.  obcią ż enia skupione,  koncentracja  naprę ż eń,  niecią głe  warunki  brzegowe)  problemy  wyboru  typu elementu czy też sposobu  dyskretyzacji  konstrukcji  wymagają   dalszych  studiów  i  ulepszeń, moż na  w  zasadzie  stwierdzić,  że  istnieją ce  procedury  numeryczne  pozwalają   na  obli- czanie  w  ekonomicznie  uzasadniony  sposób  liniowo- sprę ż ystych  powł ok  z  dostateczną dla  celów  praktycznych  dokł adnoś cią.  , W  dziedzinie  obliczeń  konstrukcji  powł okowych  pracują cych  (lub  przewidzianych do  pracy  np.  w  warunkach  awaryjnych)  w  zakresie  geometrycznie  i  materiał owo  nieli- niowym  sytuacja  jest  daleko  bardziej  skomplikowana.2*  N iezwykł a  zawił ość  analityczna (czy nawet, niekiedy, poję ciowa)  zagadnień takich jak:  redystrybucje  naprę ż eń  wynikają ce z lokalnych uplastycznień  materiał u, lokalne  odcią ż enia  w  obszarach uprzednio  uplastycz- nionych,  róż ne  rodzaje  utraty  statecznoś ci  czy  wzajemny  wpł yw  na  siebie  (sprzę ż enie) obu  przyczyn  nieliniowoś ci  (nieliniowej  charakterystyki  materiał u  i  duż ych  przemiesz- czeń)  był y  powodem,  że  dotychczas  podejmowano  tylko  niewiele  prób  peł nego,  nielinio- wego  opisu  zachowania  się   powł ok.  W  szczególnoś ci  rozwią zanie  zagadnień  utraty  sta- tecznoś ci  nie jest, i  nie może  być,  proste i  natychmiastowe.  Podejś cie  do  badan ia  statecz- noś ci  konstrukcji  róż ni się   bowiem  od przypadku  do przypadku  i nie jest  moż liwe  zbudo- wanie jednego programu zdolnego do przeanalizowania  wszystkich  problemów  zwią zanych. z  poszukiwaniem  osobliwoś ci  w  rozwią zaniach  dla  powł ok  o  dowolnym  kształ cie.  Jeś li chodzi  o  problematykę   plastycznoś ci  to  najwię ksze  kł opoty sprawia  tu  wybór  odpowied- nich  elementów  skoń czonych  zdolnych  do  zadawalają cego  opisu  rozprzestrzeniania  się stref  uplastycznionych  zarówno  wzdł uż  powierzchni  jak  i  po  gruboś ci  powł oki. 1 5  Uwagi  w oparciu o doś wiadczenia  uzyskane przy  opracowywaniu  programów  SH E LAX,  SH E L I N , SH E N ON ,  LARSTRAN ,  D YN AX. 2 )  Szerokie  omówienie  róż nych  sytuacji  praktycznych,  wymagają cych  obliczenia  konstrukcji  powł o- kowych  w  oparciu  o  teorie  nieliniowe  przedstawiono  w  [1]. 222  M .  KLEIDER Aby  wię c  podją ć  t ru d  dokł adn ej,  nieliniowej  analizy  powł oki  potrzeba  uprzedn io d o brze zrozum ieć i przewidzieć  róż ne skom plikowane  zjawiska  zachodzą ce w  kon strukcjach po wł o ko wych  i  stosown ie  do  tego  zdecydować  się   n a  taką   lub  inną   technikę   numeryczną . N a we t  dyspon ują c  gotowym  program em  nieliniowej  analizy  powł ok  uż ytkownik  bę dzie m usiał   podejm ować  decyzje  dotyczą ce  uwzglę dniania  jedn ych  i  pomijania  innych  efektów n ielin iowych  w  toku  kon kretn ych  obliczeń  num erycznych.  U wzglę dnienie  wszystkich przyczyn  ewen tualn ego  nieliniowego  zachowan ia  się   powł oki  powodować  bę dzie  bowiem prawie  zawsze  olbrzym ie  kom plikacje  numeryczne  (konieczność  szerokiego  korzystan ia z  pam ię ci  zewn ę trzn ych  maszyny  cyfrowej,  kł opoty ze  zbież noś cią   procedur  iteracyjnych, ewen tualn a  kon ieczn ość  pracy  z  wykorzystaniem  podwójnej  precyzji,  wielokrotn e  zwię k- szenie  się   kosztu  obliczeń ). Z  rozważ ań  powyż szych  wyn ika  pierwsza  ogólna  konkluzja  (n astę pne podam y  w  dal- szym  t o ku  p r a c y) : Stosują c  m etody  n um eryczn e  w  nieliniowych  zagadnieniach  powł ok  szczególną uwagę   przypisać  należy  grun town em u  zrozum ieniu  fizycznych  podstaw  odnoś nej problem at yki.  Z n aczą ce  wyniki  i  - ich  odpowiednią   in terpretację   uzyskać  m oż na wył ą cznie  w  oparciu  o  taką   wiedzę . W  szczególnoś ci  go dn e  polecenia  wydaje  się   przestudiowanie  przed  przystą pieniem d o  ko n kret n yc h  obliczeń  istnieją cych  rozwią zań  analitycznych  i  numerycznych  dotyczą - cych  u t r a t y  stateczn oś ci  i  zach owan ia  pokrytycznego  róż nych  prostych  kon strukcji  sprę - ż ystych  i  n iesprę ż ystych  takich  ja k  kolum n y,  ukł ady  prę tów,  prostsze  przypadki  pł yt. 2.  Podstawy  nieliniowej  analizy  powł ok E lem en ty  skoń czone  stosowan e  dotychczas  do  analizy  powł ok  podzielić  m oż na  n a pię ć  zasadn iczych  gr u p : a)  elem en ty  pł askie,  otrzym ywan e  n a  drodze  prostej  superpozycji  elementów  tarczowych i  pł ytowych ,  •   . b)  elem en ty  zakrzywion e,  otrzymywane  n a  bazie  róż nych  teorii  powł ok  (najczę ś ciej  n a bazie  teorii  cien kich  powł ok  m ał o  wyniosł ych), c)  trójwym iarowe  elem en ty'izoparam etryczn e,  dla  których  wprowadzon o  w  dyskretyzo- wan ej  postaci  pewn e  zał oż enia  teorii  powł ok, d)  „ n ie fo r m a ln e "  elem enty  o p art e  n a  tzw.  metodzie  „ physical  lu m pin g", e)  elem en ty  d o  an alizy  zagadn ień  specjalnych,  wś ród  których  najwię kszą   grupę   stanowią elem en ty  d o  an alizy  zagadn ień  o  osiowosymetrycznej  geom etrii  z  osiowo-   lub  nie- osiowosym etryczn ym  obcią ż eniem. Ż a d na  z  wym ien ion ych  wyż ej  grup  elementów  (a) -  (d)  nie  został a w  praktyce  uzn an a za  zdecydowan ie  lepszą   od  pozostał ych .  Elementy  „ dokł adn iejsze"  dają   lepsze  wyniki kosztem  duż ego  n a kł a d u  pracy  przy  opracowywaniu  dan ych  wejś ciowych,  skom pliko- wan ia  algorytm ów  n um eryczn ych  oraz  trudn oś ci  z  interpretacją   wyników.  Stosowanie elem en tów  prostszych  ( n p . pł askich) usuwa  wprawdzie  wię kszość  z  tych  trudn oś ci,  powo- duje  jed n a kże  róż ne  kom plikacje  szczególnie  przy  opisie  powł ok  silnie  zakrzywionych o raz  w  zagadn ien iach wraż liwych  n a  imperfekcje  geometrii. W  ogólnoś ci  stwierdzić  należ y, ż e: AN ALIZ A  POWŁOK  223 W  obecnej  fazie  rozwoju  metody  elementów  skoń czonych  wydaje  się   celowe  bu- dowanie  programów  w  oparciu  o  wszystkie  wymienione  wyż ej  rodzaje  elementów. Ostateczna decyzja  co do ich uż ycia  w konkretnym zagadnieniu spoczywać  powinna w  rę ku  uż ytkownika  i  zapadać  w  oparciu  o  zwię kszają cą   się   stale  liczbę   testów i  porównań  publikowanych  w  literaturze  ś wiatowej. Róż ne  rodzaje  elementów  wykorzystanych  w  programach  opracowanych  w  IP P T omówimy  w  rozdz.  4,  próbują c  scharakteryzować  szerzej  ich  wady  i  zalety.  Przedtem jednak  przedstawimy  pewne  podstawowe  aspekty  nieliniowego  opisu  powł ok  obowią - zują ce  niezależ nie  od  rodzaju  zastosowanego  elementu. Formuł ują c  problemy  z  zakresu  nieliniowej  mechaniki  powł ok  stosujemy  zazwyczaj przyrostowy  opis  zagadnień,  brzegowo- począ tkowych.  D la  materiał ów  sprę ż ystych  ana- lizowanych  w  zakresie  geometrycznie  nieliniowym  opis  taki  jest  jedynie  jednym  z  moż- liwych  podejść  do  rozwią zywania  zagadnień  brzegowo- począ tkowych  podczas  gdy  dla materiał ów  charakteryzowanych  przyrostowym  równaniem  konstytutywnym  jest  on koniecznoś cią. Celem naszym  bę dzie  okreś lenie  poł oż enia analizowanego  ciał a  S8 i istnieją cego  w  tym ciele  stanu  naprę ż enia  w  kolejnych,  dyskretnych  chwilach  czasu  t 0 ,  t o +At,  t o +2At,  ... gdzie  At  jest  przyrostem  czasu  zaś  t 0   począ tkową   chwilą   analizy.  Zakł adać  bę dziemy, że znamy pewne rozwią zanie  dla wszystkich  kroków  po czasie z przedział u  ft 0 ,  t]  (wł ą cznie z  chwilą , t,  którą   nazywać  bę dziemy  chwilą   aktualną ), poszukiwać  zaś  bę dziemy  rozwią - zania  w  chwili  t+At.  Sformuł owany  w  ten  sposób  problem  jest  problemem  typowym (powtarzalnym)  i  zastosowany  do  kolejnych  chwil  czasu  umoż liwia  znalezienie  cał ej drogi  poł oż eń  równowagi  ciał a. Istotną   rolę   przy  ocenie  efektywnoś ci  numerycznej,  geometrycznie  nieliniowej  analizie powł ok  odgrywają   ukł ady  współ rzę dnych, w  których  wyraż ane  są   poszczególne  równania teorii. W zagadnieniach geometrycznie  liniowych  bezdyskusyjny  jest  wybór  nieruchomego • układu  odniesienia.  Wykorzystywane  w  standardowych  algorytmach  metody  elementów skoń czonych  lokalne  (dla  każ dego  elementu), ale  stał e  w  czasie  ukł ady  współ rzę dnych służą   jedynie  wygodzie  algebraicznych  przekształ ceń  na  poziomie  elementu  i  nie  są   wa- runkowane  kinematyką   procesu  deformacji.  Odmienną   sytuację   napotykamy  w  zagad- nieniach geometrycznie  nieliniowych,  w  których  duże przemieszczenia w n aturaln y  sposób podsuwają   jeden  z  dwu. nastę pują cych  lokalnych  (dla  każ dego  elementu)  opisów: a)  opis  w  stał ym  ukł adzie  począ tkowym, b)  opis w ukł adzie konwekcyjnym  (lub w ukł adzie bę dą cym jego modyfikacją ,  n p . w ukł a- dzie  współ obrotowym). Aby  to  bliż ej  wyjaś nić  rozpatrzmy  duże  przemieszczenia  pł askiego  elementu  powł o- kowego  pokazane  na  rys.  la.  Stwierdzamy,  że  aby  w  peł ni  zachować  korzyś ci  pł yną ce z  dwuwymiarowoś ci  opisu  tych  elementów  moż emy: a)  sprowadzić  wszystkie  rozważ ania  do  konfiguracji  pierwotnej  °C  (stał ego  w  czasie ukł adu  współ rzę dnych  począ tkowych),  w  której  elementy  zachowują   swój  pierwotny kształ t,  . b)  wprowadzić  konwekcyjny  ukł ad  współ rzę dnych; ze wzglę dów  rachunkowych  wygodnie jest. ograniczyć  się   do  ukł adów  zachowują cych  swą   ortogonalność  i  prostoliniowoś ć, a  wię c  to  tzw.  ukł adów  współ obrotowych. 224 M .  KLEIBER Podejś cie  (a)  nazywać  bę dziemy  stacjonarnym  opisem  Lagrange'a  zaś  podejś cie  (b) uaktualnionym  opisem  Lagrange'a. Rozważ ania  powyż sze  są   przekonywują ce  w  przypadku  przedstawionego  na  rys.  la najprostszego,  pł askiego  elementu  powł okowego.  Jednakż e,  jak  się   okazuje,  również dla  elementów  zakrzywionych  wygodnie  jest  mieć  ortogonalny  ukł ad. współ rzę dnych b) Rys.  1 * prostoliniowych  zwią zany  z  elementem w  sposób  pokazany  na  rys.  lb,  przy  czym  należy oczywiś cie  odpowiednio  zdefiniować,  co  rozumie  się   pod  słowem  „współ obrotowy". Obroty  róż nych  czą stek  elementu  są   bowiem  róż ne  i  należy  wybrać  pewien  uś redniony, „ globalny"  obrót  elementu, definiowany  najwygodniej  przez  ruch  wę zł ów. Omawiany  ukł ad współ rzę dnych  lokalnych jest  oczywiś cie  tylko jednym  z kilku  ukł a- dów współ rzę dnych, które wprowadzić  należy analizują c  konstrukcje  powł okowe. W szcze- gólnoś ci  zawsze  wprowadzić  musimy: a)  stał y w  czasie  globalny  ukł ad współ rzę dnych  traktowany  jako  ukł ad odniesienia  przy definiowaniu  danych  wejś ciowych, b)  wspólny  dla  wszystkich  elementów  ukł ad  współ rzę dnych  umoż liwiają cy  „ zł oż enie" konstrukcji  w  cał oś ć. Przykł ady  definiowania  wszystkich  wymienionych  wyż ej  ukł adów  współ rzę dnych podamy  w  rozdz.  4.  N a  podstawie  powyż szych  uwag  stwierdzamy,  ż e: Przystę pując  do  opracowywania  programu  nieliniowej  analizy  powł ok  należy zawsze  dokonać  ś wiadomego  wyboru  konfiguracji  odniesienia  i  odpowiednich ukł adów  współ rzę dnych  kierują c  się   przy  tym: a)  wygodą   (kosztem)  obliczeń, b)  fizycznym  charakterem  równań  (np.  postulatem  przestrzennej  addytywnoś ci miar  odkształ ceń  sprę ż ystych  i  plastycznych). We  wszystkich  omawianych programach za podstawę  przyję to  klasyczną   teorię  powł ok Love'a- Kirchhoffa  oraz  klasyczne  sformuł owanie  teorii  sprę ż ysto- plastycznoś ci  Prandtla- Reussa  z  izotropowym  wzmocnieniem. Do opisu  kinematyki  procesu deformacji  stosowa- no  uaktualniany  opis  Lagrange'a  wykorzystują c  współ obrotowy,  lokalny  ukł ad współ - rzę dnych  kartezjań skich.  Konsekwencją   tego  podejś cia  był o  postulowanie  równania AN ALIZA  POWŁOK  225 konstytutywnego  w  postaci  zwią zku  mię dzy  Il- im  tensorem  przyrostu  naprę ż enia  Pioli- Kirchhoffa  odniesionym  do  aktualnej  konfiguracji  ciał a  a  przyrostem  odkształ cenia G reena zdefiniowanym  na tej  konfiguracji.2'  Takie definicje  podstawowych  niewiadomych wymagał y  przeprowadzania  odpowiedniego  procesu  akumulacji  poszczególnych  wielkoś ci przyrostowych  oraz  transformacji  wielkoś ci  globalnych  przy  koń cu  każ dego  kroku, por.  [2],  [3]. Okreś lenie  macierzy  konstytutywnej  ł ą czą cej  przyrosty  uogólnionych  naprę ż eń  (tj. sił wewnę trznych)  z  przyrostami  uogólnionych  odkształ ceń  (tj.  odkształ ceń  i  krzywizn  po- wierzchni  ś rodkowej  powł oki) wymagał o,  odmiennie niż  to jest  w  przypadku  sprę ż ystym, stosowania  metod  cał kowania  numerycznego.  Spowodowane  to  jest  skomplikowanym, nieliniowym  charakterem  rozkł adu naprę ż eń  po  gruboś ci  powł oki  i  zależ noś cią   macierzy konstytutywnej  materiał u  od  aktualnego  stanu  naprę ż enia.  Ponieważ  w  przekrojach, w  których  dominują cym  typem  deformacji  jest  zginanie  sztywność  materiał u  konstrukcji w  zakresie  sprę ż ysto- plastycznym  zmienia  się   po  gruboś ci  bardzo  znacznie, powodował o to  konieczność  uwzglę dniania  w  trakcie  procesu  cał kowania  numerycznego  duż ej  liczby punktów  cał kowania (10,  12 i wię cej  punktów). W  opracowywanych  programach  przyję to prosty  sposób  cał kowania  numerycznego  po  gruboś ci  powł oki  mają cy  fizyczną   interpre- tację   w  postaci  analizy  powł oki  podzielonej  na  warstwy. Aby  zapewnić  zbież ność  cią gów  rozwią zań  przybliż onych  otrzymywanych  w  oparciu o  metodę   elementów  skoń czonych  do  rozwią zania  dokł adnego  funkcje  aproksymują ce stan  przemieszczenia  speł niać muszą   trzy  zasadnicze  warunki: 1.  dostateczną   gł adkość  poprzez  brzegi  mię dzyelementowe, 2.  moż liwość  modelowania  jednorodnych  stanów  odkształ cenia, 3.  wytwarzanie  zerowych  odkształ ceń przy  sztywnych  ruchach  elementów. Warunki  te  są   identyczne  jak  w  analizie  liniowej.  Uwagi  dotyczą ce  ich  speł nienia przez  róż ne  elementy - zamieś cimy  w  trakcie  omawiania  poszczególnych,  zrealizowanych programów. N astę pną   uwagę   poś wię cimy  istotnej  sprawie  „ skł adania" pł askich elementów powł o- kowych  ze  sobą   w  celu  utworzenia  podstawowego  ukł adu  równań  konstrukcji.  Jako wspólny  ukł ad  współ rzę dnych  sł uż ą cy  temu  celowi  przyjmuje  się   zazwyczaj: a)  globalny  ukł ad współ rzę dnych  kartezjań skich, b)  ukł ad współ rzę dnych krzywoliniowych  na  powierzchni  ś rodkowej  powł oki, przy  czym przeważ nie  ukł ad ten definiowany  jest  poprzez  pł aszczyzny  styczne  do  powł oki  przyj- mowane jako  pewne  uś rednione  pł aszczyzny  są siadują cych  ze  sobą   elementów. Jak  podkreś laliś my  poprzednio,  poszczególnym  wę zł om  konstrukcji  powł okowej przypisuje  się  zazwyczaj  pię ć stopni  swobody, w  taki  sposób,  że wę zły  te  nie mają   sztyw- noś ci  wzglę dem  obrotu  wokół   normalnej  do  powł oki.  Aby  „ zł oż yć"  pł askie  elementy powł okowe  ze  sobą   w  globalnym,  kartezjariskim  ukł adzie współ rzę dnych,  którego  orien- tacja przestrzenna może być znacznie róż na od orientacji  ukł adów lokalnych w elementach, czę sto  konieczne jest  wzię cie  pod  uwagę   wszystkich  sześ ciu  stopni  swobody  wę zł a.  Spo- wodowane jest to koniecznoś cią   posiadania nieosobliwej  transformacji  od ukł adu  lokalnego 2 )  D efinicje  tych  wielkoś ci, a  także  uzasadnienie  koniecznoś ci  ich  uż ycia  przy  uaktualn ion ym  opisie Lagrange'a  podan o  w  pracy  [2]. 226  M .  KLEIBER do  ukł adu  globalnego  w  przypadku  koplanarnoś ci  elementów  są siadują cych  z  danym wę zł em.  W  praktyce  dokonuje  się   tej  ransformacji  przypisują c  rozpatrywanemu  wę zł owi pewną   fikcyjną   sztywność  wzglę dem  obrotu  wokół   normalnej,  stwierdzono  bowiem  n a drodze  numerycznych  doś wiadczeń  zadawalają cą   zgodność  wyników  otrzymywanych nawet  dla  znacznie  róż nią cych  się   od  siebie  wartoś ci  tej  fikcyjnej  sztywnoś ci. Z n aczn ie  proś ciej  przedstawia  się   sytuacja  w  przypadku  korzystania  w  celu  utworzenia globalnej  macierzy  sztywnoś ci  z  ukł adu  współ rzę dnych  powierzchniowych.  Przyjmują c za  pł aszczyznę   styczną   do  powł oki  pewną   pł aszczyznę   powstał ą   z  uś rednienia  pł aszczyzn elementów  są siadują cych  z danym  wę zł em  dokonać moż na  transformacji  mię dzy  poszcze- gólnymi  grupam i  pię ciu  stopni  swobody  wyraż onych  w  ukł adach  lokalnych  a  pię cioma stopniam i  swobody  wyraż onymi  we  wspólnym  ukł adzie powierzchniowym,  wpł yw  bowiem szóstego  stopn ia  swobody  na  tę  transformację   uznać moż na za pomijainie  mał y. Prowadzi to  do  znacznych  oszczę dnoś ci  czasu  obliczeń 35  ma jedn ak  wadę   w  postaci  trudnoś ci  przy analizowaniu  powł ok  o  gwał townych  zmianach kształ tu  (np. analiza  ż ebra  traktowanego jako  zwykł y  element  powł okowy).  D alszymi  zaletami  takiego  podejś cia  są : —  otrzymywanie  skł adowych  stanu  pomieszczenia  i  obrotu  w  ukł adzie  współ rzę dnych ś ciś le  zwią zanym  z  powł oką , —  brak  jakichkolwiek  kł opotów  w  przypadku  koplanarnoś ci  wszystkich  elementów są siadują cych  z  danym  wę zł em. Jak  wiadom o,  bardzo  szeroką   klasę   nieliniowych  zagadnień  mechaniki  continuum, w  tym  również  zagadnienia  nieliniowej  mechaniki powł ok,  opisać  moż na w  ramach prze- mieszczeniowego  modelu  metody  elementów  skoń czonych  ukł adem  równań  algebraicz- nych  postaci (1) gdzie / 2"\   Jf  _  irkon st.  ,  - g- pocz.priem.  i  g- pocz.napr.  i  gn iekon s.obc. p r z y  c z ym  w p r o w a d z o n o  t e  o z n a c z e n i a : K  —  gl o b a l n a  m a c i e r z  szt ywn o ś ci  ko n st r u kc ji, • gkonst  —m a c i e r z  szt ywn o ś ci  k o n st yt u t ywn e j, jjpocz.przem.  —m a c i e r z  p o c z ą t k o wych  p r ze m ie szc ze ń , jpocz.n apr.  —m a c i e r z  p o c z ą t k o wych  n a p r ę ż e ń, jcniekons.obc.,—  m a c i e r z  o b c ią ż eń  n ie ko n se r wa t ywn yc h , C  —  m a c i e r z  t ł u m i e n i a , M   —  m a c i e r z  m a s,  •   '  , R t   —  we k t o r  o bc ią ż e n ia  ze wn ę t r zn e go  wę zł ów  (w  ch wili  t), r t   —  we k t o r  u o gó l n i o n yc h  p r ze m ie szc ze ń  wę zł ów  ( w  ch wili  t)., Jt- at  —we kt o r  „ wewnę trznych"  sił   wę zł owych  odpowiadają cy  poprzedniemu krokowi  analizy. 3 )  U wzglę dnienie  6- go  stopnia  swobody powoduje  wzrost czasu  rozwią zywania  podstawowego  ukł adu / 6 \ 3 równ ań  okoio  I —I  w 1.728  razy. AN ALIZA  POWŁOK  227 Równanie  (1),  zależ nie  od  geometrii  i  materiał u  konstrukcji,  rodzaju  obcią ż enia,  zasto- sowanego  opisu  (stacjonarny  lub  uaktualniony)  oraz  wybranego  ukł adu  współ rzę dnych, doprowadzić  moż na  do  jawnej  postaci,  dla  której  wybrać  należy  odpowiednią   m etodę rozwią zywania. 3.  Uwagi  na  temat  algorytmów  rozwią zywania  nieliniowych  zagadnień  powł okowych M etody rozwią zywania  statycznych  i dynamicznych zagadnień  powł okowych  nie róż nią się   zasadniczo  od  metod  rozwią zywania  odpowiednich  problemów  mechaniki con tin uum . Poniż ej  dokonamy  krótkiego  przeglą du  tych  metod  zaś  do  klasyfikacji  tej  nawią ż emy w  rozdz.  4  wskazują c  na  podejś cia  zastosowane  w  zrealizowanych  program ach. (A)  Zagadnienia  statyczne Pomijają c  czł ony  dynamiczne  zapiszmy  równanie  (1)  w  postaci (3)  (K L  + K 2 ) Ar  =  AR- J* gdzie  Ki  jest  czę ś cią   macierzy  sztywnoś ci  K  stał ą   w  trakcie  procesu  deformacji,  K 2  = s  K—Kx,  zli? jest  przyrostem  obcią ż enia  zewnę trznego  na  danym  kroku  zaś  / *  odpo- wiednio  zmodyfikowanym  wektorem  wypadkowych  sił   wę zł owych  wyraż ają cych  fakt niezrównoważ enia  wę zł ów.  N iezrównoważ enie  t o  jest  wynikiem  bł ę dów  przybliż eń  n u- merycznych  powstał ych  na  poprzednich  krokach  analizy.  M acierz  K 2  jest  funkcją   aktual- nego  stanu  powł oki  co  symbolicznie  zapiszemy  jako  K 2   — K2(ff, r,  ...). Istnieją   dwie  podstawowe  metody  rozwią zywania  zagadnienia  przyrostowego  opi- sanego  równaniem  (3),  które  scharakteryzujemy  krótko  poniż ej. (Al)  Metoda zmiennej  sztywnoś ci.  M etoda  ta  polega  n a  rozwią zywaniu  równ an ia  (3) dla  kolejnych  przyrostów  obcią ż enia  uwzglę dniając  zmiany  macierzy  sztywnoś ci  odpo- wiadają ce  rozprzestrzenianiu  się   stref  plastycznych  w  analizowanej  powł oce. P o  wykon an iu obliczeń  danego  kroku  należy  akiimulować  otrzymane wielkoś ci  przyrostowe  i  n a  podsta- wie  obliczonych  wielkoś ci  globalnych  charakteryzują cych  proces  (naprę ż enie,  odkształ - cenie  plastyczne,  parametry  wzmocnienia)  okreś lić  nową   macierz  sztywnoś ci  powł oki. M etoda  powyż sza,  ze  wzglę du  n a  swą   prostotę   i  ł atwość  interpretacji  poszczególnych operacji  algebraicznych  w  ram ach  mechaniki,  jest  szeroko  stosowana  we  wszystkich prawie  programach  dotyczą cych  dyskretyzowanych  m etod  w  nieliniowej  analizie  powł ok. Rozpatrzmy  dowolny  / - ty  krok  analizy  prowadzonej  w  oparciu  o  równ an ie  (3),  i  = =   1, 2,  . . . , «.  Wprowadzimy,  dla  dowolnej  funkcji  /   charakteryzują cej  rozpatrywany proces,  symbol  / ( i )  na  oznaczenie  wartoś ci  funkcji  /   odpowiadają cej  począ tkowi  z- tego kroku,  zaś  symbol  'Afi''i+1>  m a  oznaczenie  wartoś ci  przyrostu  fu n k c ji/ n a  tym  kroku. W  metodzie  zmiennej  sztywnoś ci  stosujemy  algorytm  postę powania  dan y  równ an iam i (4) 228  M.  KLEIBER M etoda  zmiennej  sztywnoś ci  nie jest  wolna  od  wad.  Zaliczyć  do  nich  należy  przede wszystkim: 1.  konieczność wielokrotnego  (na każ dym  kolejnym  kroku) rozwią zywania  podstawowego ukł adu  równań  z  nową ,  uaktualnioną   macierzą   sztywnoś ci, 2.  mał a  efektywność  metody w zagadnieniach prowadzą cych  do niesymetrycznej  macierzy sztywnoś ci  (np.  przy  uwzglę dnieniu  niestowarzyszonych  praw  pł ynię cia  plastycznego materiał u  powł oki  czy  też  obcią ż eń  zewnę trznych  typu  ś ledzą cego), 3.  trudne  do  oszacowania  bł ę dy  „odchodzenia"  od  rozwią zania  dokł adnego  wraz  ze wzrostem  liczby  przyrostów, 4.  trudnoś ci  zwią zane  ze  zł ym  uwarunkowaniem  globalnej  macierzy  sztywnoś ci  w  oto- czeniu  punktów  krytycznych. Istnieje  wiele metod poprawiania  dokł adnoś ci rozwią zań  uzyskiwanych  metodą  zmien- nej  sztywnoś ci.  Jedną   z  nich, tzw.  metodę  sił  korekcyjnych,  opisać  moż na  formalnie  za- leż noś ciami = 0 , (5) Stosowanie  metody  zmiennej  sztywnoś ci  uwzglę dniają cej  sił y  korekcyjne  wymaga wykonania  jedynie  niewielkiej  liczby  dodatkowych  operacji  mnoż enia,  potrzebnych  do okreś lenia  wektora  danego  wzorem  ( 5) 6.  N ie  wymaga  ono  natomiast  rozwią zywania nowego  ukł adu  równań  ani  też  stosowania  jakichkolwiek  procedur  iteracyjnych. (A2) Metoda począ tkowych  obcią ż eń.  Metoda  ta  ma  charakter  iteracyjny;  ft- tą  iterację   pa ż - tym  kroku  analizy  przyrostowej  opisać  moż na  w  sposób  nastę pują cy: (6) gdzie  • »;  jest  zadaną   dokł adnoś cią   procesu  iteracyjnego.  Proces  powyż szy  kontynuujemy aż  do  momentu, w  którym  dwa  kolejno  obliczone  wektory  wę złoyfych  obcią ż eń  począ t- kowych  7*1**1,  j*o*+D  róż nić  się   bę dą   dostatecznie  mał o  od  siebie. N a każ dym  kroku analizy proces  iteracyjny  rozpoczą ć moż na zakł adają c  ^ |r ( M + 1 : 0 )  =  0 lub, co jest  bardziej e f e k t y w n e ,  A r < ' - i + 1 : O i  =   A ^ 1 ^ AN ALIZ A  POWŁ OK  229 Zasadniczą  zaletą  metody  obcią ż eń  począ tkowych  jest  konieczność  tylko  jedn okrot- nego  odwrócenia  macierzy  sztywnoś ci  (która  jest  stał a  w  trakcie  procesu  deformacji), wada  zaś  natomiast  czę sto  napotykane  kł opoty  ze  zbież noś cią  procedury  iteracyjnej. N ie istnieją  bowiem  w  przypadku  powł ok ż adne ś cisłe kryteria  dotyczą ce  tej  zbież noś ci*). Bardzo  istotną  sprawą  przy  stosowaniu  metody  począ tkowych  obcią ż eń  jest  sposób kolejnego  okreś lania  sil  wę zł owych,  bowiem  w  praktyce  nie  korzysta  się  tu  bezpoś rednio ze wzoru  (6)! wprowadzonego  powyż ej  jedynie  w  celu  uzyskania  wię kszej  zwartoś ci  opisu omawianego  algorytmu  obliczeń.  Dwie,  zasadniczo  róż ne  od  siebie,  metody  okreś lenia dodatkowych  sił  wę zł owych  noszą  w  literaturze  nazwy:  metoda  począ tkowych  naprę ż eń i metoda począ tkowych  odkształ ceń. W  praktyce  okazuje  się, że zbież ność  metody począ t- kowych  naprę ż eń jest  znacznie  sł absza  od  zbież noś ci  metody  począ tkowych  odkształ ceń i wł aś nie ta druga metoda stosowana jest powszechnie w programach  sprę ż ysto- plastycznej analizy  powł ok. (B)  Zagadnienie dynamiczne Podstawowe  algorytmy  rozwią zywania  ukł adu  równań  (1)  omówimy  na  przykł adzie prostszego  przypadku  dynamicznej  analizy  zagadnień  liniowo- sprę ź ystych  opisanych równaniem (7)  Mr t   + Cr t   + Kr t   =   £ t . Procedura  tworzenia  globalnych  macierzy  mas  i  tł umienia jest  analogiczna  do  procedury budowy  globalnej  macierzy  sztywnoś ci  w  zagadnieniach  statycznych35.  Polega  ona  na bezpoś rednim  sumowaniu  odpowiednich  macierzy  poszczególnych  elementów, utworzo- nych w lokalnych  ukł adach współ rzę dnych i przetransformowanych  do ukł adu  globalnego. Podstawy  konstrukcji  elementowych  macierzy  tł umienia  są  bardzo  sł abo  zbadane co  wynika  z  ogólnego  braku  informacji  na  temat  mechanizmów  pochł aniania  energii w  skomplikowanych  ukł adach  konstrukcyjnych.  Zamiast  więc  prób  bardziej  ogólnego definiowania  tej  macierzy  ustala  się  zazwyczaj  pewne  parametry  tł umienia  dla  poszcze- gólnych  postaci  drgań  powł oki  w  oparciu  o  doś wiadczenia  uzyskane  przy  analizie  po- dobnych  ukł adów.  U stalone  w  ten  sposób  parametry  tł umienia  modalnego  mogą  być bezpoś rednio  wykorzystane  do  analizy  o  ile  stosujemy  przy  rozwią zywaniu  tzw.  metodę superpozycji  modalnej  lub  też  przetransponowane  do  postaci  równoważ nej  macierzy tł umienia  C  w  przypadku  bezpoś redniego  cał kowania  równań  ruchu.. D la  wszystkich  elementów  powł okowych  istnieje  moż liwość  budowania  tzw.  kon- syŝ entnej  macierzy  mas  w  oparciu  o  te  same  zał oż enia  aproksymują ce,  które  przyję to przy  konstrukcji  macierzy sztywnoś ci.  Konsystentną  macierz mas  okreś la  się  n a  podstawie równania,  które  symbolicznie  napisać  moż na  jako (8)  Ji  - - 4 )  K ryt eria  t akie  istn ieją  t ylko  d la  n iekt ó r yc h ,  n ajpro st szych  za ga d n ie ń  z  za kr e su  sp r ę ż yst o - p la st yc z- n o ś ci,  n p .  zagad n ień  p ł askiego  st a n u  n a p r ę ż en ia  w  m a t e r ia le  ze  wzm o c n ien iem  i z o t r o p o wym . 5 )  P eł ną  an alo gię  wykazuje  r ó wn ież  zagad n ien ie  szó st ego  ( o b r ó t  wo kó ł   n o r m a ln e j  d o  p o wł o ki)  st o p n i a swo bo d y  w  wę ź le. 6  Mech.  Teoret.  i  Slos.  2/80 230  M .  KLEIBER. gdzie  Jł   jest  macierzą   m as  rozpatrywan ego  elementu,  Q gę stoś cią  masy  m ateriał u zaś 0 wektorem  odpowiedn ich  funkcji  kształ tu.  N a podstawie  (8) stwierdzamy,  że dla  funkcji kszt ał t u  bę dą cej  wielom ian em  n p .  3- go  stopn ia  przy  okreś lan iu  m as  cał kowan iu  podlega wielom ian  st o p n ia  6- go.  W przypadku  korzystan ia  z  procedur  cał kowan ia  num erycznego prowadzić  t o m oże  d o kon ieczn oś ci  uż ycia  znacznej  liczby  pun któw  cał kowan ia  czynią c z  takiego  po st ę po wan ia  m etodę   bardzo  kł opotliwą   i  kosztowną   (w  porówn an iu  z  n u- m eryczn ym  cał kowan iem  potrzebn ym  d o  okreś lenia  macierzy  sztywnoś ci  tego  samego elem en t u ) . Altern atywn ą   procedurą  jest tworzen ie  uproszczonej, diagonalnej  macierzy  m as w opar- ciu  o  zał oż en ie  istn ien ia  skon cen trowan ych  m as  w pu n kt ach  wę zł owych.  Z e wzglę du n a kł o p o t y  n um eryczn e  p rzy  okreś lan iu  macierzy  konsystentnej  oraz  zupeł nie  zadowalają ce wyniki  otrzym ywan e  przy  korzystan iu  z  macierzy  diagonalnej  w  wielu  program ach dy- n am iczn ej  an alizy  powł ok  stosuje  się  to  drugie  podejś cie.  D iagon aln a  macierz  m as za- pewn ia  bowiem  dużą   efektywność  cał kowan ia wzglę dem  czasu  równ an ia  (7) n awet w przy- p a d ku  zn aczn ej  liczby  stopn i  swobody  an alizowan ego  zakł adu. P on iż ej  sch arakteryzujem y  kr ó t ko  dwa  podstawowe  podejś cia  stosowan e  d o  przy- bliż on ego  cał kowan ia  ró wn ań  ruch u  (7), za które  uważa  się  grupę   m etod  bezpoś redniego cał ko wan ia  oraz  m et o d ę   superpozycji  m odaln ej. (Bl)  Metody  bezpoś redniego  całkowania. W  m etodach  tych  równ an ie  (7)  cał kowan e  jest wzglę dem  czasu  m etodą   „ kro k  p o  kr o ku ",  zaś  okreś lenie  „ bezpoś redn ie  cał kowan ie" ozn acza,  że  przed  przystą pien iem  do  cał kowan ia  równ an ia  tego  nie podajem y  ż adnej tran sform acji.  M et o d y  bezpoś redn iego  cał kowania  oparte  są   n a  dwu  podstawowych zał oż en iach.  P o pierwsze  zakł adam y, że równ an ie  (7) m a być speł niane tylko  w  wybra- n ych  chwilach  czasu  t 0 ,  t o +At,  ł o +2At,  . . . , t 1   nie  zaś  w cał ym przedziale  czasu  [t 0 ,  f j . P o  drugie,  zakł adać  bę dziemy  z  góry  pewną   zm ienność  przemieszczeń,  prę dkoś ci  i  przy- spieszeń  w  zapatrywan ym ,  typowym  przedziale  czasu  [t, t+At].  P onieważ  zaś równ an ie (7)  jest  wektorowym  równ an iem  róż niczkowym  o  stał ych  współ czynnikach,  d o  wyra- ż en ia  prę dkoś ci  i  przyspieszeń  w  funkcji  przemieszczeń  zastosować  m oż na  w  zasadzie dowoln ą   aproksym ację   róż nicową   poch odn ych  czasowych.  Jak się   jedn ak  okazuje,  tylko kilka  wyraż eń  prowadzi  do dostateczn ie  efektywnych  algorytm ów. W  m etodzie  róż n ic  cen traln ych  przyjmujemy h = co  w  poł ą czen iu  z  (7) prowadzi  d o c 1 \ h- Af[(At) 2   2At\ Wzór  (11)  przedstawia  ukł ad  równ ań  liniowych, z którego  znają c  r t , r t - ń t   oraz  R t   wyliczyć m o ż na  r, +At .  P ro ced u rę   powyż szą   nazywamy  metodą   cał kowan ia  jawn ego ;  zauważ my, że  p r o c ed u r a t a n ie wym aga  odwracan ia  macierzy  sztywnoś ci  K co, szczególnie  gdy  m am y d o czyn ien ia z diagon aln ym i  m acierzam i m as i tł um ienia, jest jej  olbrzymią   zaletą .  D odajm y AN ALIZA  POWŁOK  231 po n ad t o ,  że pon ieważ  do obliczenia  r t+At   pot rzebn a jest  zn ajom ość  r t   i r,_ M ,  dla t  = t 0 zachodzi  konieczność  opracowan ia  pewnej  procedury  „ st a r t o wej",  kt ó r a  w  celu  okreś- lenia  r to +4t  nie  korzystał aby  z wielkoś ci  r ta - M .  I n n ą , zn aczn ie  poważ niejszą   wadą   m etody róż n ic  centralnych jest  konieczność  doboru  dł ugoś ci  kro ku  po czasie  At  w  t aki  sposób, aby  był a  on a  mniejsza  od  pewn ego  czasu  krytycznego  A t kr ,  zależ n ego  od  wł asn oś ci  cał ego ukł adu.  W  przeciwnym  bowiem  wypadku  otrzym ywan e  rozwią zanie  cechuje  brak  sta- bilnoś ci.  Z achodzić  powin n o T gdzie  T „ jest  najmniejszym  okresem  drgań  wł asnych  u kł ad u .  Waru n ek  ten  je st  bard zo silnym  ograniczeniem  dł ugoś ci  kro ku  i, m im o iż w  pewnych  sytuacjach  n ie m usi  być  o n w  peł ni  przestrzegan y,  stwarza  poważ ne  przeszkody  w  efektywnym  wykorzystan iu  o m a- wianego  algorytm u. D rugą   m etodę , którą   przedstawim y  w niniejszym  opracowan iu, jest m et o d a cał ko wan ia niejawnego  typu  Wilsona,  zwan a  także  metodą   p aram et ru  &. W m etodzie  tej zakł ad am y (12)  '  r t +i  =   % gdzie  T e  [0, • d- At],  • &  >  1, co prowadzi  do zależ noś ci (13)  %  ( } P odstawiają c  (13)  i  (14)  d o równ an ia  ruchu  wypisan ego  dla chwili  t+&At  otrzym ujem y równ an ie, z którego  wyliczyć m o ż na f, + ^ r  a n astę pn ie otrzym ać r,+jt,  h+At»>"t+A t •  W algo- rytm ie  tym n ie zachodzi  p o t rzeba  opracowywan ia  dodatkowej  pro cedu ry  „ st a r t o we j", bowiem  przemieszczenia,  prę dkoś ci  i przyspieszenia  w chwili  t+At  wyraż one  są  w  funkcji wielkoś ci  zdefiniowanych  jedyn ie  dla chwili  t.  M et o d a  Wilson a  wym aga  trian gularyzacji macierzy  sztywnoś ci  dla każ dej  dyskretnej  chwili  czasu.  M et o d a  t a jest  bezwaru n ko wo stabiln a  dla • & >  1,37; powszechn ie  przyjmuje  się  & — 1,4. (B2) Metoda superpozycji  modalnej.  M etoda  t a  polega  n a  przetran sform owan iu  ró wn an ia ruch u  (7) do prostej,  rozprzę ż on ej  p o st a c i6 ) (15)  x(t)  + i22x(t)  -   P(t), otrzym anej  za  pom ocą   m acierzy  transform acji  utworzon ej  z  rozwią zań  0 L ,  & 2 ,  ...,  &„ zagadn ien ia  n a wartoś ci  i  wektory  wł asne  postaci (16)  K * =  w2M«f>, przy  czym  Q2  jest  diagon aln ą   macierzą   o  wyrazach  bę dą cych  kwad rat am i  kolejn ych wartoś ci  wł asnych  oĄ ,m\ ,  ...,, z}  (a wł aś ciwie  ukł ad  {r, z}  na wybranej  poł udniowej pł aszczyź nie przekroju cp  =   const.). Ze wzglę du na brak moż liwoś ci  obrotu wokół  normalnej do powł oki w  przypadku  powł ok  osiowosymetrycznych  nie  wystę pują  opisane  w  rozdz.  3  problemy z  wprowadzaniem  fikcyjnej  sztywnoś ci  odpowiadają cej  temu  stopniowi  swobody. 234  M .  K L E I BE R Program  napisany  jest  w ję zyku  FORTRAN   IV  i  skł ada  się   z  okoł o  1400  wyraż eń (kart  perforowanych);  w  trakcie  jego  opracowywania  korzystano  z  maszyny  cyfrowej CYBER  72,  na  której  też  policzono  wszystkie  dotychczasowe  przykł ady. Przykł ad  1 Z a  pomocą   programu  SHELAX  przeprowadzono  analizę   czaszy  kulistej  poddanej dział aniu  sił y  P  przył oż onej  w  wierzchoł ku  czaszy,  rys.  3 7 ) .  Omawianie  tego  problemu rozpoczniemy  od  przytoczenia  wyników  dotyczą cych  sprę ż ystej  czaszy  utwierdzonej  na obwodzie.  D la  przyrostu  obcią ż enia  AP  =   lib  stwierdzono  pewną   rozbież ność  wyników w  porównaniu  z  rozwią zaniem  ś cisł ym  podanym  w  [9], szczególnie  w  ozę ś ci  ś rodkowej cał ego  zakresu  analizowanych  deformacji.  Wią że  się   to  z  najmniejszą   sztywnoś cią   kon- strukcji  w  tym  zakresie.  W  celu zwię kszenia  dokł adnoś ci obliczeń zmniejszono  pię ciokrot- nie przyrost  obcią ż enia  co  wpł ynę ło  ńa znaczną   poprawę   dokł adnoś ci rozwią zania.  Oma- wianą   czaszę   poddano  nastę pnie  analizie  w  zakresie  sprę ż ysto- plastycznym.  Odpowiedni wykres  obcią ż enie- przemieszczenie  podano  na  rys.  4,  na  którym  zilustrowano  również charakter  powstają cych  stref  plastycznych  (naniesionych na nieodkształ coną   konfigurację powł oki). Zmieniają c  warunki  brzegowe  powł oki  z  utwierdzenia  na  przegubowe  podparcie nieprzesuwne  obserwujemy  jakoś ciową   róż nicę   w  otrzymywanym  rozwią zaniu.  Umoż li- wiliś my  bowiem  w  ten  sposób  wystą pienie  zjawiska  globalnej  utraty  statecznoś ci  czaszy w  postaci  tzw.  przeskoku.  Wartoś ci  obcią ż enia  krytycznego  dla  powł oki sprę ż ystej  i sprę - ż ysto- plastycznej,  jak  również  odpowiednie  kształ ty powł oki  w  chwili  przeskoku  podano n a  rys.  5. (B)  program  SH E N O N 8) — statyczna  analiza  duż ych  przemieszczeń  i  statecznoś ci cienkich  powł ok  dowolnego  kształ tu,  [1],  [4],  [5], [6]. Program  SH EN ON   umoż liwia  analizę   w  wyż ej  wymienionym  zakresie  dowolnych powł ok  sprę ż ystych  oraz  sprę ż ysto- lepkoplastycznych.  Wykorzystują c  stacjonarne  wł as- noś ci  rozwią zań  sprę ż ysto- lepkoplastycznych  program  umoż liwia  również  otrzymywanie rozwią zań  sprę ż ysto- plastycznych.  W  programie  uwzglę dniono  nastę pują ce  powł okowe elementy  skoń czone  typu  przemieszczeniowego: a)  pł aski  element  trójką tny  charakteryzowany  przez: (al)  stan  membranowy —  element z  liniową   funkcją   kształ tu, (a2)  stan  zgię ciowy —  zgodny  element  pł ytowy  zaproponowany  w  [7], b)  niepł aski  element  czworoką tny  charakteryzowany  przez: (bl)  stan  membranowy — ukł ad  czterech  pł askich  elementów  trójką tnych  z  kwadra- tową   funkcją   kształ tu  poddany  wię zom  kinematycznym  prowadzą cym  do  linio- woś ci  przemieszczeń  wzdł uż  czterech  zewnę trznych  boków  czworoką ta, (b2)  stan  zgię ciowy —  ukł ad  czterech  elementów  pł ytowych  wymienionych  w  (a2). 7 )  Ab y  u m o ż liwić  p r zep r o wa d zen ie  p o r ó wn a ń  o t rzym an yc h  wyn ikó w  z  ro zwią zan iami  o t r zym a n ym i p r z e z  i n n yc h  a u t o r ó w,  w  p r zykł a d zie  t ym  ( i  p a r u  in n ych ) przyję to  u kł a d jed n o st ek  st o so wan y  powszech n ie w  U S A. 8 )  O p is  p r o gr a m u  S H E L I N ,  [4], bę d ą c ego  lin iową   wersją   an alizy  p o wł o k  n a  bazie  kt ó rej  o p r a c o wa n o wprsję   n ielin io wą   za t yt u ł o wa n ą   S H E N O N ,  zawart y  jest  a u t o m a t yc zn ie  ( jako  czę ść  o pisu  p r o gr a m u  S H E - N O N ) '  w  p o n iż szym  t ekś cie. o  Metoda  zmiennej  sztywnoś ci,  AP=1 Ib •   Metoda  zmiennej  sztywnoś ci  ,&P- 1/ 5 I b —'  rozw.  dokT. \vg  [ 9 3 Czas  obliczeń  ICDC 6600): CP~  750 sek  przy 120  krokach H=0.08589 1  P(lb) 90 80 70 6 0 ' 50 40 30 20 10 0 E  = 10 7  psi ugię cie  ś r od ka  powł oki  (in) Rys.  3 P=3S Ib rozwiazqnie  sprę ż yste  P= 2 7 l b = 1  Ib) P= 20lb rozwiqzanie  plastyczne (Przyrost  obc  AP=1lb ) Zmiany  kształ tu  powierzchnil ś rodkowej  powłoki E= 1 0 7 p s i )>  = 0.3  p si S_- ?.in« psi 10  war st w 2 0  e l e m e n t ó w P=5  Ib 0,05  0,10  0,15  p ionowe  p r zem . wier zc ho ł k a  ( in) Rys.  4 Rozwój  stref  up last yc znio nyc h (konfig urac ja  nieoc lkszt alc ona) [235] AN ALIZA  POWŁOK  237 Przy opracowywaniu  programu zdecydowano  się  na element zł oż ony z  czterech pł askich elementów  trójką tnych  ze  wzglę du  n a: 1.  brak  w  naszych  krajowych  warunkach  doś wiadczenia  w  korzystaniu  ze  skompliko- wanych  programów  numerycznych  w  ogóle,  a  z  programów  dotyczą cych  powł ok w  szczególnoś ci, 2.  konieczność  minimalizacji  wykorzystywanego  obszaru  pamię ci  wewnę trznej  maszyny, 3.  prostotę   i  szybkość  przygotowywania  danych  wejś ciowych, 4.  wzglę dnie  dobrą  aproksymację   rzeczywistej  geometrii powł oki w przypadku niepł askiego elementu  czworoką tnego  przy  zachowaniu  niewielkiej  liczby  danych  wejś ciowych, 5.  przejrzystość  programu  i  ł atwość jego  rozbudowywania, 6.  ł atwą   interpretację   otrzymywanych  wyników, 7.  moż liwość  (prawie  we  wszystkich  przypadkach)  redukcji  bł ę dów  wprowadzonych w trakcie idealizacji  powierzchni zakrzywionej  elementami pł askimi poprzez  zagę szcza- nie  siatki  elementów. Przeprowadzone  testy  wskazują   na  znacznie  lepszą   dokł adność osią ganą   przy  uż yciu elementu  czworoką tnego  dlatego  też  elementy  trójką tne  zaleca  się   uż ywać  wył ą cznie jako  uzupeł nienie siatki  elementów  czworoką tnych. Macierze  sztywnoś ci  elementów  wyprowadzono  korzystają c  z  bardzo  wygodnych tzw.  polowych  współ rzę dnych  wprowadzonych  w  obszarze  każ dego  z  trójką tów.  Odpo- wiednie  funkcje  kształ tu  oraz  otrzymane  na  ich  podstawie  macierze  sztywnoś ci  podano w  [4],  [5]. Opracowany  element czworoką tny  posiada  37  stopni  swobody,  z  których  jed- nakże tylko  20  ma charakter  stopni zewnę trznych.  Reszta, tj.  17 jest  eliminowana  w pro- cesie  statycznej  kondensacji  na  szczeblu  elementu  zapewniają c  bardzo  korzystną ,  nie- wielką   szerokość  pasma  globalnej  macierzy  sztywnoś ci. Omawiany  element  czworoką tny  jest  elementem  zgodnym  tylko  w  przypadku  ko- planarnoś ci  wszystkich  swoich  trójką tnych  podelementów;  brak  zgodnoś ci  w  przypadku niepł askim nie spowodował  w dotychczasowych  testach zauważ alnych  problemów ze zbież- noś cią   cią gów  rozwią zań  przybliż onych.  Opracowana w  programie macierz  począ tkowych naprę ż eń jest  tzw.  macierzą   niekonsystentną   tzn.  opartą   na nieco uproszczonym ukł adzie aproksymują cych  funkcji  kształ tu,  por.  [5],  U proś ciło  to  znacznie  procedurę   jawnego otrzymywania  macierzy  począ tkowych  naprę ż eń  nie prowadzą c  (w  dotychczasowych  .za- stosowaniach)  do ż adnych  niekorzystnych  wł asnoś ci takiego  sformuł owania. N umeryczne cał kowanie potrzebne przy okreś laniu tej macierzy przeprowadza się  w programie wykorzy- stują c  wzory  P.  C. H ammera,  O. P . Marlowe'a  i  A. H . Strouda  podane  np.w  [8], s.  421. „ Skł adanie"  macierzy  poszczególnych  elementów w  globalną   macierz ukł adu  odbywać się   może  w  programie  alternatywnie  w  oparciu  o : a)  globalny  ukł ad  współ rzę dnych  prostoką tnych  dla  translacyjnych  stopni  swobody (3  skł adowe) i  powierzchniowy  ukł ad  współ rzę dnych  dla  obrotowych  stopni  swobody (2  skł adowe), b)  powierzchniowy  ukł ad  współ rzę dnych  dla  wszystkich  pię ciu  skł adowych  stanu  uogól- nionych  przemieszczeń. Zgodnie z poprzednimi uwagami  ż adna z  tych metod nie prowadzi  do  trudnoś ci zwią - zanych  z  brakiem  w  wę ź le  szóstego  stopnia  swobody  odpowiadają cego  obrotowi  wokół normalnej  do powł oki. 2 3 8  M .  K.LE1BER Podstawową   metodą   rozwią zywania  ukł adu równań  zastosowaną   w  programie SHE- N O N  jest  metoda począ tkowych  obcią ż eń  w  odniesieniu  do nieliniowoś ci  typu  fizycznego oraz  metoda  zmiennej  sztywnoś ci  w  odniesieniu  do  nieliniowoś ci  typu  duż ych  przemiesz- czeń. Program  SH EN ON   pozwala  na  znaczne  uproszczenie  i  zmniejszenie  liczby  danych wejś ciowych  poprzez  wykorzystanie  szeregu  podprogramów  generacyjnych.  W  szczegól- noś ci  program  umoż liwia: a)  generację   współ rzę dnych  wę zł ów  dla  pię ciu  czę sto  spotykanych  typów  powierzchni, b)  generację   cosinusów  kierunkowych  współ rzę dnych  powierzchniowych, c)  generację   numerów  wę zł ów, d)  generację   danych  materiał owych, e)  generację   przemieszczeniowych  warunków  brzegowych, f)  generację   obcią ż eń  zewnę trznych. Program  SH EN ON   napisany jest  w ję zyku  FORTRAN  IV  i  skł ada  się   z  okoł o 3.000 wyraż eń  (kart perforowanych);  w trakcie jego  uruchamiania korzystano  z maszyny  cyfro- wej  CYBER  72. Przykł ad  2  .  .  . Za  pomocą   programu  SH EN ON   1  przeprowadzono  analizę   duż ych  ugię ć  wycinka sfery  kulistej  obcią ż onego  siłą   skupioną   P,  rys.  6.  N a  wykresie  przemieszczenie  ś rodka powł oki —  obcią ż enie  zewnę trzne  porównano  wyniki  otrzymane  przez  róż nych  autorów. (C)  LARSTRAN  —  analiza  duż ych  odkształ ceń  konstrukcji  sprę ż ystych,  sprę ż ysto" plastycznych,  sprę ż ysto- lepkoplastycznych. System  LARSTRAN  jest duż ym systemem nieliniowej, dynamicznej analizy  konstrukcji opracowanym  w  Instytucie  Statyki  i  Dynamiki  Uniwersytetu  w  Stuttgarcie,  [13],  [14]. Istnieją ce  w  tym  systemie  moż liwoś ci  analizy  powł ok  opisane  został y  w  pracach  [15]', [16],  [17]. Przedstawiony  w  tych pracach element TRU M P może być  stosowany  w bardzo wielu  zagadnieniach  praktycznych  aczkolwiek  autorzy  z  naciskiem  podkreś lają   jego przybliż ony  charakter i,  zwią zaną   z  tym, jego jedynie  „inż ynierską"  dokł adnoś ć.  Mówią c ogólnie,  proponują c  element TR U M P  kierowano  się   przekonaniem, że  obniż ona  dokł ad- ność rozwią zań  rekompensowana  bę dzie  niż szym  kosztem  obliczeń  i  mniejszym  zapotrze- bowaniem n a pamię ć komputera, co w skomplikowanych  problemach nieliniowej mechaniki powł ok  może  być  sprawą   na  wagę   moż liwoś ci  uzyskania  jakiejkolwiek  uż ytecznej  infor- macji.  Element TR U M P  zaliczyć  trzeba  do klasy  elementów  „nieformalnych", jego  przy- datność  bę dzie musiał a zostać  w  przyszł oś ci  potwierdzona  na  drodze  róż norodnych  eks- perymentów  numerycznych. Omawiany  element  skoń czony  jest  trójką tnym,  pł askim  elementem  powł okowym z  trzema  punktami  wę zł owymi,  z których  każ dy  posiada  sześć  stopni  swobody  w postaci trzech  przemieszczeń  i  trzech  obrotów.  Zgodnie  z  poprzednimi  uwagami  konieczność wprowadzenia  sztywnoś ci  elementu wzglę dem  obrotu wokół  normalnej do powł oki wynika ze  skł adania  macierzy  sztywnoś ci  w  globalnym  ukł adzie  współ rzę dnych;  sztywność  ta ma  charakter  czysto  „ numeryczny". N aturalny opis  stanu deformacji  elementu powoduje, że  element  posiada  pozornie jedynie  12  stopni  swobody  ( 3x6  stopni  swobody  wę złów minus  6  stopni  swobody  odpowiadają cych  sztywnym  ruchom  elementu  w  przestrzeni trójwymiarowej).  Zastosowanie  koncepcji  elementu  warstwowego  do  analizy  powł ok AN ALIZA  POWŁOK 239 niesprę ż ystych  odpowiada  metodzie  omówionej  na  przykł adach  programu  SH ELAX; odpowiednio budowana  macierz sztywnoś ci  elementu uwzglę dnia  jego  wł asnoś ci  membra- nowe  i  zgię ciowe  zaś  struktura  warstwowa  umoż liwia  opis  zmian  wł asnoś ci  materiał u „ po  gruboś ci"  powł oki: a  = 61,8031"  (156.99 cm) h  (grubość  powłoki) 3.9154" I  9.95cm) Rc  =100.0"  (254,0cm) E  (moduł  Younga)=  105psi y(wsp.Poissona)=0,3. Leicester (szeregil,[111 Thomas  [ 1 2 ] Boland  [ 1 0 ] SHENON  1 O  0,2  0,4  0,6  0,8  1.0  1.2  1,4  1,6  V Rys.  6 W  programie  zał oż ono,  że  poszczególne  warstwy  elementu  TR U M P  odpowiadają najprostszemu  elementowi  pł askiego  stanu  naprę ż enia  charakteryzowanemu  stał ym- rozkł adem  naprę ż enia. Ponieważ jednak  element zł oż ony z  warstw  o  stał ym  naprę ż eniu nie jest w stanie opisać deformacji  elementu typu antysymetrycznego  zginania w  programie wykorzystano  dodatkowo  analogię   typu  pł yta trójką tna  — ukł ad  trzech  zginanych  belek. W  trakcie  statycznej  analizy  konstrukcji  niesprę ż ystych  istnieje  w  systemie  LARS- TRAN   moż liwość  wykorzystywania  wszystkich  omówionych  w rozdz.  3 procedur  analizy nieliniowej  tj.  metody  zmiennej  sztywnoś ci  oraz  obu  wersji  metody  począ tkowych,  ob- cią ż eń.  Element  TRU M P  jest  obecnie  wykorzystywany  także  do  zagadnień  dynamiki powł ok.  . ' • . • .. 240 M .  KLEIBER Przykł ad  3 Jednym  z  testów  elementu TRU M P  był a analiza  powł oki omówionej  w Przykł adzie 1, [18].  Wycinek  w  postaci  1- go  radiana  kopuł y poddano  dyskretyzacji  za  pomocą   36 ele- mentów  TR U M P .  Sprę ż ysto- plastyczne  rozwią zanie  O- A- B- C,  rys.  7,  otrzymane  przy AP  =   0,5/ 6  wykazał o  dość  dobrą   zgodność  z  rozwią zaniem  otrzymanym  za  pomocą ,  H  0,08589" [ Ib ] . 50 40 30 20 10 przyrostowe  rozw.spc leptoplast.,LARSTRAN 90 przyrostów AP=0,5lb Co B, A rozw. spr.- plabt„ lASTfWi 90 przyrostów  AP= 0,5  Ib vrozw.  spr- piast, SHELAX 90przyrostowAP = 0,5 Ib 0,05 0,10 1/ sec,  E - 107 Rys.  7 0,15 6 0 =2 - 1 0 4 p si. SH ELAX- u.  W  celu  wykorzystania  wł asnoś ci  rozwią zań  sprę ż ysto- lepkoplastycznych opisanych  w  [2],  [17],  [18] przył oż ono nastę pnie w jednym  kroku  obcią ż enie  odpowiada- ją ce  punktowi  A o   i  nastę pnie, w  procesie  pł ynię cia lepkoplastycznego  pod  stał ym obcią - ż eniem  otrzymano  leż ą cy  na  odpowiedniej  krzywej  punkt  A.  Podobnie  dla  obcią ż eń odpowiadają cych  punktom  B o   i  C o  otrzymano  punkty  B  i  C,  a  nastę pnie rozpatrzono także  inne  historie  obcią ż enia  (np. 0~A o —B x —B~  C 2  — C), wszystkie  wykazują ce  zgod- ność  z  odpowiednim  rozwią zaniem  sprę ż ysto- plastycznym.  Liczba  iteracji  konieczna do  uzyskania  odpowiedniej  zbież noś ci  procesu  numerycznego  opisują cego  pł ynię cie lepkoplastyczne  wahał a  się   od  37  (w  przypadku  pł ynię cia  C o  - >•   C)  do  8  (w  przypadku pł ynię cia  C 3  -> C). N a rys.  8 przedstawiono  rozwią zanie  tego samego  zadania dla róż nych wartoś ci  współ czynnika  lepkoś ci;  otrzymane  krzywe  wykazują   tendencję   do  asympto- tycznego  zbliż ania  się   do  rozwią zania  sprę ż ysto- plastycznego  (y  =  oo). AN ALIZA  POWŁOK 241 50  - 4 0 - 3 0 - 2 0 - 10 rozwią zania  przyrostowe 90x0,5 lb 1 1 1 \ / \ \   rozw. spn piast. (y=oo) \   \ rozw.spr lepkoplast. \   \   T- 1500  1/sec. \   \ rozw. spn lepkoplost. \   7=1000  1/sec: \ razw. spr lepkoplast. T=500  1/sec, 0,05 0,10  0,15 Rys.  8 0,20 in] (D )  WH AM  — nieliniowa  dynamiczna  analiza  ciał   i  konstrukcji System  WH AM   jest  duż ym  systemem  nieliniowej,  dynamicznej  analizy  konstrukcji opracowanym  na  U niwersytecie  N orthwestern  w  Chicago,  U SA  przez  T.  Belytschke i  B.  Mullena.  Istnieją cy  w  tym  systemie  element  powł oki  osiowosymetrycznej  opisano w  [19]. Element ten jest  podobny  do  elementu istnieją cego  w  programie  SH ELAX.  Jego tworzą ca  opisywana  jest  równaniem,  por.  (19) zaś  funkcje  aproksymują ce  stan  przemieszczenia  przyję to  w  postaci,  por.  (20) Rys.  9 242 M.  KLEIBER F unkcje  te  opisują   stan  przemieszczenia  w  elemencie  skoń czonym  o  sześ ciu,  stopniach swobody  i  dwu  wę zł ach zewnę trznych,  rys.  9. Zmianę   wł asnoś ci  materiał u  wzdł uż  gruboś ci  powł oki  uwzglę dniono  wprowadzają c model  warstwy.  D o  cał kowania  równań  ruchu  zastosowano  metodę   róż nic centralnych. Przykł ad  4 Moż liwoś ci  obliczeniowe  programu  WHAM   zilustrujemy  na  przykł adzie  analizy utwierdzonej  n a  obwodzie  powł oki  sferycznej  obcią ż onej  impulsem  ciś nienia,  rys.  10. N a  rysunku  tym  pokazano  zmianę   ugię cia  ś rodka  kopuł y  w  czasie  dla  przypadków: sprę ż ystego  i  sprę ż ysto- plastycznego. R- h . d » E = 2227  i 0.4  in 26.67 ° 10.5- 10 0 , ^2 4 ' i r j 3 EP-•   2.UWS n S  % J lb/ m2 ?=> 2.45»10 ,-4  lbsec ł - Q06 m at er ial  sprę ż ysty czas (msec^ material sprę ż ysto- plaatyszny Rys.  10 (E)  D YN AX  — dynamiczna  analiza  osiowosymetrycznych  ciał   i  powł ok  sprę ż ystych poddanych  dowolnemu  obcią ż eniu. P rogram  D YN AX,  [20], stanowi  pierwszą   fazę  pracy nad programem  umoż liwiają cym dynamiczną   analizę  róż norodnie obcią ż onych  konstrukcji  o osiowosymetrycznej  geometrii w  zakresie  fizycznej  i  geometrycznej  nieliniowoś ci. W  odróż nieniu od programu  SHELAX i  powł okowej  czę ś ci  programu  WH AM   program  D YN AX  dopuszcza  powstawanie  do- wolnego  (a  nie  tylko  osiowosymetrycznego)  stanu  odkształ cenia  w  analizowanych  kon- strukcjach).  Aby  umoż liwić  taką   ogólność  analizy  zmienne wokół  osi  symetrii  obcią ż enie przedstawiane jest  w  postaci  szeregów  F ouriera. N a każ dym  kroku  analizy  otrzymujemy rozwią zanie  dla  każ dego  wyrazu  rozwinię cia  oddzielnie  i  sumujemy  otrzymane  rozwią - zania  otrzymują c  wynik  odpowiadają cy  cał emu  obcią ż eniu. AN ALIZA  POWŁOK 243 P rogram  dysponuje  obecnie  trzem a  typam i  elem en tów: 1.  elem ent  powł okowy  identyczny  do  wykorzystywan ego  w  program ie  SH E L AX, 2.  „ t ró jką t n y"  elem en t  ciał a  osiowosym etrycznego  (grubej  powł oki)  z  lin iową   funkcją kształ tu  w  pł aszczyź nie  poł udn ikowej, 3.  „ czworoką t n y"  elem en t  ciał a  osiowosymetrycznego  zł oż ony  z  czterech  elem en tów „ t ró jką t n ych"  bez  wę zła  ś rodkowego  elim in owan ego  w  procesie1  statyczn ej  ko n d en - sacji. D la  wszystkich  elem en tów  w  program ie  istnieją   moż liwoś ci  obliczania  kon systen tn ych macierzy  m as.  Współ czyn n iki  leż ą ce  n a  diagonali  w  wygodniejszej  w  an alizie  diagon aln ej macierzy  m as  otrzym ywan e  są   poprzez  sum owan ie  wszystkich  współ czyn n ików  wystę pu- ją cych  w  poszczególnych  wierszach  macierzy  kon systen tn ej. M acierz  tł um ien ia  przyjm owan a  jest  w  postaci  proporcjon aln ej  jako C  =   aM + / SK gdzie  współ czynniki  KfleHBi  ocHOBHbie  acneKTBi  H ejiH H eftH oro  aHaxcn3a  c  npHMeHeHKeiw 3 B M T O H K H X  o6onoM;eK  npon3B0JibH 0H   cj>opMM.  y*iHTbiBaeTCH   H3nqecKaH   H   reoM eTpnqecK H OC TŁ  aw,ann-   n peflcraBJieH Bi  O C H O BH Ł K  amropHTMBij  npH MeiraeMBie  iip n  nodaH OBKe H  flH KaM H ^ecKH X 3a, ą aq.  P aSoTa  n pojuun ocTpupoBaH H aH   MH orc- mcneH H biMu  npnM epai«n  p ac i&ro B, Bbi- n o jm e H H t ix  c  noM ombM   n porpaM M   Ha  3 B M 3  pa3pa6oT am ibix S u m m a r y N O N LI N E AR  STATIC  AN D  D YN AM I C  AN ALYSIS  OF  SH ELLS  BY  TH E F I N I TE ELEM EN T  M ETH OD ,  I n the paper  some  fundamental  aspects  of  numerical  nonlinear  analysis  of'free- form  thin  shells are discussed. Both material an d geometrical  nonlinearities  are taken into account. The basic numerical methods are  described  for  nonlinear  static  and  dynamic  problems.  A  number  of  examples  is  included. I P P T  PAN Pracą   został a  zł oż ona  w  Redakcji  dnia 20  marca 1979 roku