M E C H AN I K A TEORETYCZNA I STOSOWANA ,. 3, 12 (1974)
ZASTOSOWANIE ITERACJI SEIDLA W METODZIE ELEMENTÓW SKOŃ CZONYCH NA PRZYKŁADZIE OBLICZEŃ STATYCZNYCH PŁYT
[JERZ Y G O Ł A Ś , Z YG M U N T K A S P E R S K I , AN N A P E E R - K A S P E R S K A , JERZ Y M A K O W S K I (OP OLE)
1. Wprowadzenie
W metodzie elementów skoń czonych, jak wiadomo, zasadniczymi czynnikami każ dego program u są nastę pują ce podprogram y: podprogram do budowania macierzy sztywnoś ci poszczególnych elementów, podprogram do budowania ukł adu algebraicznych równań liniowych oraz podprogram do rozwią zywania ukł adu równań. Zbież noś ć wyników obli-czeń numerycznych do rozwią zań dokł adnych zależy gł ównie od wyboru kształ tu elementu, postaci funkcji przemieszczeń oraz gę stoś ci podział u ustroju cią gł ego na elementy. Zwię k-"szenie liczby pun któw wę zł owych co prawda lepiej aproksymuje ustrój cią gł y, jednak pro-wadzi do duż ego ukł adu równań i zasadniczo limituje moż liwoś ć stosowania maszyn cyfrowych do obliczeń. Pojawia się zatem problem opracowywania programów wykorzystują -cych jak najoptymalniej pojemność pamię ci maszyny cyfrowej. D użą efektywność daje zastosowanie iteracji Seidla do rozwią zywania ukł adu równań. Pozwala ona na wielo-stopniowe budowanie ukł adu równań, bez pamię tania cał ej macierzy sztywnoś ci ukł adu, oraz na etapowy charakter obliczeń, a pon adto pozwala na peł ne wykorzystanie wszyst-kich wł asnoś ci macierzy sztywnoś ci, takich jak: pasmowoś ć, symetria oraz regularny rozkł ad wyrazów zerowych wystę pują cych w paś mie. A zatem, każ de równanie zapamię -tywane jest w postaci zwartej, z pominię ciem wyrazów zerowych.
Zasadnicze procedury program u wykorzystują cego wyż ej wymienione moż liwoś ci przedstawione zostaną na przykł adzie pł yt izotropowych, przyjmują c prostoką tny element skoń czony o 12 param etrach wę zł owych [1], W przypadku przyjmowania wię kszej liczby stopni swobody dla elementu prostoką tnego, np. 16, czy też jak dla pł yty z ż ebrami jedno-stronnymi i pł yty trójwarstwowej minimum 20, podstawowa istota budowania odpowied-nich podprogram ów nie ulegnie zmianie. Jedynie wymiary podmacierzy sztywnoś ci ele-mentu, zamiast wymiarów (3 x 3), przyjmą odpowiednio wielkoś ci ( 4x4) i (5 x 5) oraz wię ksza bę dzie liczba warunków brzegowych, [2], [3].
Przedstawiony algorytm obliczeń numerycznych swą ogólnoś cią obejmuje analizę statyczną pł yt dowolnie obcią ż onych, pł yt o kształ tach zł oż onych z prostoką tów, pł yt o do-wolnych warunkach brzegowych i podpartych w swoim obszarze oraz pł yt osł abionych otworami prostoką tnym i. Sztywność gię tna pł yt może być zmienna, lecz stał a w obrę bie elementu.
N iektóre elementy przedstawionego algorytmu mogą być rozszerzone n a inne zagadnie-nia metody elementów skoń czonych.
342 J. G OŁ AŚ, Z . KASPERSKI, A. PEER- KASPERSKA, J. MAKOWSKI
2. Algorytm programu
2.1. Dyskretyzacja ukł adu. Rysunek 1 przedstawia prostoką tny obszar pł yty, który po-dzielono na z = (n — 1) (m—l) elementów skoń czonych o dowolnych wymiarach oraz ilustruje przyję tą numerację punktów wę zł owych i poszczególnych elementów. Przyję ty na rysunku kształ t nie ogranicza zastosowań programu do pł yt prostoką tnych. W przy-padku pł yt o innych kształ tach (ale dają cych się podzielić n a elementy prostoką tne) od-powiednim elementom należy n adać bardzo mał ą grubość w porównaniu do gruboś ci
(1) k 3 -- ca y (m- l)n'l 2n.1 n*l
©
i h, nrt)n>2 GE£> n*2 I la,„a
3® j.1II
1 i1
fl'4 4 ! 1 1 II ln- 1 OfT) n- l na- l mn (m- i)n QmD (ED n- l nK,
X Rys. 1pł yty rzeczywistej. D opeł n ien ie obszaru pł yty do p ro st o ką ta pozwala zach ować cią gł ość numeracji, a t ym sam ym moż liwe jest cał kowite zautom atyzowan ie obliczeń, przy prostym sposobie przygotowywan ia dan ych . Jedyn ym i wielkoś ciami okreś lają cymi podział pł yty n a elementy są: liczby n, m (liczby p u n kt ó w wę zł owych znajdują cych się odpowiedn io n a poziomej i pion owej krawę dzi), wym iary aubj (i = 1, 2, ...,n; j = 1, 2, ..., ni), szerokość i wysokość i- tego pion owego i / - tego poziom ego pasa podział u pł yty, grubość tk(k =
= 1,2, ..., z) poszczególnych elem en tów.
2.2. Sposób tworzenia macierzy sztywnoś ci ukiadu. Warunki równowagi dla poszczególnych punktów wę zł owych podział u konstrukcji prowadzą do ukł adu równań
(2- 2.1)
gdzie [K]
{<5}
{F}
macierz sztywnoś ci ukł adu,
wektor uogólnionych przemieszczeń, wektor obcią ż eń.
ZASTOSOWANIE ITERACJI SEIDLA 343
M acierz sztywnoś ci ukł adu [K] buduje się z odpowiednio skł adanych i sumowanych podmacierzy kpq, bę dą cych elementami macierzy sztywnoś ci poszczególnych elementów. D la pł yty izotropowej macierz sztywnoś ci elementu prostoką tnego (rys. 2) przyjmuje postać (2.2.2) Rys. 2 gdzie elementy lciP są podmacierzami kwadratowymi o wymiarach (3 x 3). Podmacierze te wyznaczono w postaci jawnej i przykł adowo wynoszą : (2.2.3) f PP -m k — Ką ą
l- 2v
Wab
(1 +4v)
10a
l+4v
106
" (7- 2w)
10a6(1- 0
10a
(1 +4iO
106
[10a
106
I 1 1 1 1 _ ,6
a
3a
b
2b
b a a3 + b3a
b
2b
a
2b
2a
3a
b
2b
2a
1a
1 63(1 +4v)
10a
46(1 - v)15a
a (1 —v) b3 10a 6(1 — v)15a
r\ (1 +4v) a 10a ' b246(1 - 0
Ą a ť15a ' 36 ab
z4a
36
4
1+4
106
4a(l
" + b + a2 - v) , 46 156 ' 3a_a
b
2la
36
1+4J»106
3(1 - 0156
(1 +4v)
106 'o
4a{l- v)156
- 6 ~a
246
1 3a>
6 ~2a
226
3a
gdzie 2a, 26, t oznaczają odpowiedn io: dł ugoś ć, szerokość i grubość elementu prostoką tne-go pł yty.
344 J. G OŁAŚ, Z . KASPERSKI, A. PEER- KASPERSKA, J. MAKOWSKI
W dalszym cią gu, gdziekolwiek bę dzie mowa o wierszu bą dź kolumnie macierzy [K\ , to bę dziemy rozumieli pod tym odpowiednią trójkę wierszy lub kolumn, jak to pokazano na rys. 3. Przez i oznaczono numer pun ktu wę zł owego podział u pł yty, i = I, 2, ...,nxm.
-• i : : — -r: _
-1
4-—
if
- ,
:
T
- liJ
-11 i i ~ f. i I i -i -_ h . - ' i . ii ~ i - li —I i U -:ih_ ) - i - J+r , i J l -— ... • V 1 • i- :-Objaś nienie zapisu. Zapis (2.2.4) Rys. 3d- k
pq- * (u,v) : (h)
bę dzie oznaczał : należy wzią ć wymiary a, b, t elementu skoń czoneg o o numerze d, utwo-rzyć macierz kM, dodać ją do wyrazów macierzy [K] leż ą cych w wierszu u i w kolumnie v.
Z uwagi n a symetrię macierzy [K], bę dziemy budowali tylko jej wyrazy leż ą ce n a i nad przeką tną , dlatego u < v. Symbol (Ii) bę dzie wyjaś niony w dalszej czę ś ci pracy.
W zależ noś ci od poł oż enia punktów wę zł owych podział u pł yty rozróż nia się nastę -pują ce przypadki :
Przypadek I
i — 1 (wę zeł leż ą cy w lewym dolnym naroż u)
1—fc
w- > (l,l) :(1),
1 - * H - * ( 1 . « + 1) :( 4 ) , 1 — kpr —y ( 1, n+2) : (5),
1- *:„ - > (1,2) :(2).
Przypadek II
i — 2,3, ...,rt— 1 (wę zł y leż ą ce n a dolnej poziomej krawę dzi pł yty, nie liczą c wę zł ów
naroż nych) i- l i—l i- l i i i K sq ^ — k„ -»• . / , . — k pp ~* — kps - *• (i > (', (l> (i, (i, (i, (i, rt+i—
n+i)
i) i) rt + 1) ; + i ) O : (3), :(4), :( 1) , :( 1) . :(4), 1) : (5), :( 2 ) .Z ASTOSOWAN IE I TER AC JI SE I D LA 345
Przypadek III
i ~ n (wę zeł leż ą cy w prawym dolnym naroż u)
rt- 1— fc,s- > («, 2/ 7- 1) :( 3 ) ,
n- \ — / cs, - > (n,2n) :( 4) ,
«—1 — kss ->• ( », ») : (1). Przypadek IV
i = « + l , 2 « + l , ..., ( m —2) « + l (wę zł y leż ą ce n a lewej krawę dzi nie liczą c wę zł
ów na-roż n ych ). O góln ie: i =jn + l, gdzie./ = 1, 2, 3, . . . , m —2
(. / - 1 )«- 7 + 2 — V - + (/, i +1) : (2),
y« - y + 1 — / <pg - > (f, i+ «) : (4),
jn- j+l - ,kpr - + (/ , J+ W+ l) : (5),
Przypadek V
jn+k, gdzie fc = 2, 3, ..., « — 1; y = 1, 2, ..., m - 2 (wę zł y leż ą ce wewną trz obszaru
pł yty) U- l) ( r t - 1 ) + / c ~ 1 — /c,.r - > (/, 0 : (1), ( j- 1 ) ( «- 1 ) + / c — / c,r -*• (/, j + 1) : (2), j(n- l)+k- l - ksq - . (/, / + » - ! ) : (3), , / («- 1) +k- 1 — /c.„. - > (z, / + «) : (4), j(n- \ )+k - \ —kss - > ( / , / ) : ( 1 ) ,
./ (« - 1 ) + A — fc
pe- » (/ , i + «) : (4),
: ( 5),
Przypadek VIi = (j+l)n, gdziej = 1, 2, ..., m —2 (wę zł y leż ą ce na prawej krawę dzi nie liczą c wę zł ów
naroż nych)
./ (rt- 1) — fcrr- * (/ ,;) - .(I ).
C/ + \ )(n- 1) - / cs, - . 0\ i+tt- 1) • (3),
Przypadek VII
i = (m~l)n + l (wę zeł leż ą cy w lewym górnym naroż u)
346 J. G OLAS, Z. KASPERSKI, A. PEER- KASPERSK:A, J. MAKOWSKI
Przypadek VIII
i= (m- l)n+k, gdzie k = 2, 3, ..., n—l (wę zł y leż ą ce n a górnej krawę dzi, nie liczą c
wę zł ów naroż nych)
- l —*: „ - > {i, i) : (1),
: (2).
Przypadek IX
i — mxn (wę zeł leż ą cy w prawym górnym naroż u)
M acierz [AT] jest macierzą pasmową o duż ej liczbie wyrazów zerowych. P rzykł adowo: /- ty wiersz macierzy [K] m a postać przedstawioną na rys. 4, gdzie pola zakreskowane oznaczają wyrazy niezerowe. Przy duż
ych n liczba wyrazów zerowych wzrasta bardzo szyb-\
• I S
L+n- 1 w; ' • Rys. 4 I \ 2 5 4- S \ Rys. 5ko. Łatwo zauważ yć, że maksymalna liczba wyrazów niezerowych w każ dym wierszu [K] (po uwzglę dnieniu symetrii) wynosi 15. D latego też, każ dy wiersz zapamię tywany jest w tablicach W l, W 2, W 2> o wymiarach [1:15], [1:14] [1:13] (rys. 5), daje to oszczę dność pamię ci maszyny cyfrowej, a pon adto skraca czas obliczeń, ponieważ n a elementach ze-rowych nie wykonuje się ż adnych operacji arytmetycznych.
Zapamię tanie każ dego wiersza [K] w postaci zwartej powoduje zmianę indeksów poszczególnych wyrazów. U jmuje to liczba h (zob. 2.2.4), oznaczają ca numer pozycji, do której ś cią gnię to (przesunię to) odpowiedni wyraz / - tego wiersza. Wyrazem tym może być jedna, bą dź suma dwóch i wię cej podmacierzy (2.2.2).
2.3. Algorytm rozwią zywania ukł adu równań liniowych z uwzglę dnieniem warunków brzegowych. Rozwią zanie ukł adu równań (2.2.1), gdzie
Klt Kit K- rt Kr, d2
F
r Ftt m 3mxn jest liczbą stopni swobody ukł adu, otrzymuje się m etodą iteracyjną Seidla z
czynnikiem nadrelaksacji co. Odpowiedni wiersz macierzy [K] umieszczony jest w tabli-cach W l, W 2, W 3.
Z ASTOSOWAN I E ITERAC II SE I D LA 347
Proces iteracyjny w metodzie Seidla odbywa się cyklicznie. Wpierw obiera się dowolne przybliż enie d0 = [<5£ ,6%, ..., <5Ó]
T
, a nastę pnie w C/4- 1) kroku iteracyjnym oblicza się wektor <5J+ 1 wedł ug wzoru
r- l I
(2.3.1) a 5 + i - . a $ - - j ( gdzie Jfy są wyrazami [K].
Czynnik nadrelaksacji przyś pieszają cy zbież ność rozwią zania powinien speł niać nierów-ność 1 < co < 2, [4].
Ze wzoru (2.3.1) wynika, że w przeciwień stwie do metod dokł adnych rozwią zywania ukł adu równań, elementy macierzy sztywnoś ci ukł adu [K] nie zmieniają się. Pozwala to, wykorzystać wszystkie wł asnoś ci tej macierzy, jak: symetrię, pasmowość oraz regularny rozkł ad duż ej iloś ci elementów zerowych wystę pują cych w paś mie. D o obliczenia przybli-ż eń r- tej niewiadomej potrzebny jest tylko r- ty wiersz macierzy [K],
Tok postę powania poszczególnych etapów przedstawionej metody w programie na maszynę cyfrową jest nastę pują cy (zob. schemat blokowy, rys. 6):
a) dla dowolnego punktu wę zł owego i tworzy się tablice W\ , W2, W3 wedł ug sposobu podanego w punkcie 2.2,
b) wyznacza się odpowiednią trójkę prawych stron równania (2.2.1),
c) sprawdza się, czy dany wę zeł jest wę zł em brzegowym, jeś li tak, to jakiego rodzaju warunek brzegowy wystę puje tam . W przypadku n p. cienkich pł yt izotropowych rozróż nia się nastę pują ce warunki brzegowe:
I (<5r = 0, <5r+ 1 = 0, <5p+2 = = - O), I I (0P = O, ór + 1 = 0, <Y+2 * 0), ( 1 3 '2 ) I I I (<T - 0, dr+1 * 0, <5r+ 2 = 0), IV (<Sr = 0, dr+i ± 0, ór+2 jk 0), gdzie: dr oznacza ugię cie, zaś ór + 1 i dr+z
oznaczają ką ty obrotów kolejno wzglę dem osi równoległ ych do osi x i y.
Warunki brzegowe do pamię ci maszyny cyfrowej wprowadza się podając na taś mie kolejno: rodzaj warunku brzegowego, odpowiednią ich liczbę oraz numery wę zł ów, w którym dany warunek wystę puje. Wówczas w każ dym kroku iteracyjnym za odpowiednią niewiadomą podstawia się wartość zerową;
d) oblicza się <ST wedł ug wzoru (2.3.1) biorąc za krq odpowiedni wyraz tablic W\ , W2, W3>. Przy czym pamię tać należy równocześ
nie o tym, że w tablicach W wyrazy nie ozna-czają kolejnych wyrazów r- tego wiersza macierzy [K], lecz są to tylko niezerowe elementy tego wiersza. W zwią zku z tym niewiadome <5j+1
oraz d] we wzorze (2.3.1) muszą wystę-pować w odpowiedniej kolejnoś ci.
Po wykonaniu iteracji dla wszystkich niewiadomych ukł adu bada się wartość wyraż enia (2.3.3) .
Wielkość powyż szego wyraż enia charakteryzuje szybkość zbież noś ci do rozwią zania dokł adnego oraz okreś la w pewnym sensie wielkość bł ę du bezwzglę dnego, Jeś li u ^ e, gdzie e jest zadaną mał ą wielkoś cią, wówczas obliczanie przemieszczeń koń czy się, przy
348 J. G OLAS, Z . KASPERSKI, A. PEER- KASPERSKA, J. MAKOWSKI
czym przemieszczenia przyjmują takie wielkoś ci, jakie otrzymano w ostatnim kroku iteracji. Jeś li natomiast u > e, to zaczyna się tok postę powania od pun ktu a) do d), czyli wykonuje się nastę pny krok iteracji aż do otrzymania u ^ s.
Po rozwią zaniu ukł adu równań (2.2.1), czyli po wyliczeniu przemieszczeń <5, oblicza się wielkoś ci sił wewnę trznych, co nie przedstawia już wię kszych trudnoś ci i nie bę dzie omówione w niniejszej pracy.
JTAHT 1 Czytanie warunków brzegowych
Tworzenie uporzą dkowanego zbioru Z,
którego elementu sq numerami nienia-- domych równycti 0( warunki brzegowe)
Czytanie • • n,m,9,Ę ,Q
a; bj ł , (wymiary i gruboii elementu)
£,d, o(dokł .przijbl'.począ tkcme,ai/ miinak,
(nadanie m
UtNorzenie i- fego wiersza macierzy
ukł adu KO'Fdla odppw. wę zł a i
podstawienie pod tablice wi wz W Czu wezetiestorzeiiomii tak me Nadanie wartoś ci zerowych odpowiednim niewiadomym
Uizz
tak keu (V ? me tak Bruk na rn.ą njocz 1 oraz nitorze I nie 'ke<j[Q)\ tak nie me tak Obliczanie i wydruk sił wewnę trznych1
STOPRys. 6. Ogólny schemat blokowy programu zrealizowanego na maszynie cyfrowej «ODRA 1204» Zapis k e y (0 oznacza pytanie, czy został wciś nię ty klawisz o numerze i dla maszyny cyfrowej «ODRA 1204»
Z ASTOSOWAN I E I TER AC JI SE I D LA 349
3. U wagi koń cowe
Autorzy są dzą , że przedstawiona metoda iteracji może być wykorzystana wszę dzie tam, gdzie otrzymuje się dużą liczbę równ ań , których rozwią zanie metodą dokł adną np. metodą eliminacji G aussa przekracza moż liwoś ci pojemnoś ci pamię ci maszyny cyfrowej.
]27 120 115 106 99 02 85 76 71 54 57 50 43 J6 29 22 15 B 1
1 ł
28 86 ff 2 29 j 8' 50 116 109 102 95 83 61 • li 57 50 53 39 32 Z5 78 11 4 Sl 89 47 J i * ł ł ł i 52 155 126 119 112 105 W K 91 u 77 70 53 50 U 49 44 35 n 21 U S 7 ' j ,i X _»»_4
Rys. 7Wadą metody jest to, że oblicza się wielokrotnie te same wielkoś ci tzn. powtarza się zawsze tok postę powania od pun ktu a) aż do d). Jest to jednak konieczne ze wzglę du na ograniczone pojemnoś ci maszyn cyfrowych. Czas obliczeń zależy od gę stoś ci podział u
350 J. GOŁAŚ, Z . KASPERSKI, A. PEER- KASPERSKA, J. MAKOWSKI
obszaru pł yty na elementy i od odpowiedniego doboru czynnika nadrelaksacji co. Przykł a-dowo dla okoł o 1000 stopni swobody i co = 1,8 czas obliczeń wynosił 2,5 godziny.
M oż na pon adto wysuną ć nastę pują ce wnioski:
1. M etody iteracyjne posiadają waż ną wł asność tzw. autokorekcji.
2. Przedstawiony program moż na rozszerzyć n a zagadnienia statyki pł yt trójwarstwo-wych, pł yt z ż ebrami jednostronnymi itp.
3. Modyfikacja wektora obcią ż eń w każ dym kroku iteracyjnym, przez wprowadzenie czynnika nadrelaksacji znacznie skraca czas obliczeń na maszynie cyfrowej.
4. Przykł ad liczbowy
Opracowany program sprawdzono n a szeregu przykł adów liczbowych. D la pł yty trój-polowej, (rys. 7), przy podziale n a 108 elementów skoń czonyc h otrzymane wyniki mo-mentów zginają cych w porównaniu z wynikami znanymi z literatury [5] róż nią się o okoł o 0,5%- s- 1,5% (rys. 8).
Rys. 8. Wykres momentów zginają cych M w przekroju A •— A; O — wyniki znane z literatury, ... A • • • — wyniki otrzymane
Literatura cytowana w tekś cie
1. O. C. ZtENKrEWicz, Metoda elementów skoń czonych, Arkady, Warszawa (1972).
2. R. GANOWICZ, J. GOŁAŚ, Pewne problemy teorii pł yt z ż ebrami jednostronnymi, Rozpr. Inż ., 3,15, (1967). 3. J. GOŁAŚ, Pewne rozwią zania dla koł owych pł yt trójwarstwowych obcią ż onych na krawę dzi, Rozpr. Inż.
19, 2, (1971).
4. A. RALSTON, Wstę p do analizy numerycznej, Warszawa 1971.
5. Poradnik inż yniera i technika budowlanego, Praca zbiorowa, Arkady, Warszawa 1968.
P e 3 IO M e
n P H M E H E H H E H TEPA1JH H 3EH,Ii;jL$I K M ETOD Y K OH E ^H LI X 3J I E M E H T 0B H A ITPHMEPE CTATH *IECKH X P AC TE TOB I U I H T
B pa6oie paccMOTpena nporpaiwiwa craTiraecKoro pac<jeia3 MeTOflOM Kone^mbix 9JieMeHTOB3
npHMO-yronBHbix rniiiT co CKatmoo6pa3HO iweiwnomeHCH JKCCTKOCTLMJ n p n npoH3BOJiŁHbix Harpy3Kax H npoH 3-BOJIŁHŁK KpaeBwx ycnoBHHX. PenienH
e BŁiTeKaiomeii H3 nponeflypŁi flHCKpeiH3aiiHHj CHCTeiwbi ypaBHe-Z ASTOSOWAN IE I TE R AC JI SE I D LA 351
HHH noJiy^eH O c npHMeHeHHeM HTepau,noHHoro iweTofla 3eHflJiH co CBepxpejiaKcaqneii. ripH MeireH ue MeTofla 3eH,n;juł flano BO3MOWHOCTB M H ororayneH ^aToro nocTpoemiH cHCTeMM ypaBtieHHH n p n oflHOBpe-MeHHOM ee peineH H H . TIpoi- paiwMa n posepeH a Ha pnfle iH CJiemn.ix npninepoB. PaccMaTpuBaeMbiii MeTOfl MOJKeT HcnoJib3OBaTbcH fljiH n p yr u x sansn pemaeMLix meTOfloiw KOHe^HLix ajieM emoB.
S u m m a r y
APPLICATION OF SEID EL'S ITERATION TO TH E F IN ITE ELEMEN T METH OD , ON TH E SAMPLE OF STATIC PLATE ANALYSIS The paper presents a program of the finite element method applied to the static analysis of rectangular plates with jumptype variable rigidity, under arbitrary boundary conditions and loads. The set of equations resulting from the discretization process is solved by the iterational Seidel method with super- relaxation. Application of the Seidel method enables the multi- step construction of the set of equations and its si-multaneous solution. The program is verified on several numerical examples. The procedure discussed may be applied to other finite element problems.
WYŻ SZA SZKOŁA IN Ż YN IERSKA, O P O L E