MKE za dvodimenzionalnu i trodimenzionalnu elastičnost
U ovoj glavi pokazano je kako se problem dvodimenzionalne i trodimenzionalne elastičnosti
može numerički riješiti metodom konačnih elemenata. Predstavljen je matrični koncept MKE, osnovni elementi koji se nalaze i u bibliotekama softverskih paketa koji rješavaju problem elastičnosti metodom konačnih elemenata, kao i način formiranja globalne matrice krutosti sistema.
5.1 Matrica krutosti pravougaonog KE sa četiri čvora Jedan od pristupa da se numerički r ješi iješi problem određivanja napona i deformacija tanke ploče na slici 5.1(a) čije analitičko rješenje nije poznato je da se traži približno rješenje pretpostavljajući polja pomjeranja na cijelom domenu koje će približno zadovoljiti jednačine ravnoteže. Drugi pristup je da se ploča podijeli na ograničen broj poddomena (konačnih elemenata, slika 5.1(b)) i da se polje pomjeranja na cijelom domenu aproksimira serijom jednostavnijih funkcija funkcija na poddomenima. poddomenima.
131
F 1
F 1
F 2
F 2
(b)
(a)
Slika 5.1 Tanka ravna ploča pod dejstvom opterećenja (a) i podjela
geometrijskog domena na mrežu pravougaonih KE (b)
Na slici 5.2(a) prikazan je jedan pravougaoni KE kojim je diskretizovan domen ploče. Konačni
element je preko zajedničkih čvorova vezan za ostale KE. Ut icaj susjednih konačnih elemenata zamjenjuje se silama u čvorovima, slika 5.2(b). U čvorovima KE može djelovati i vanjsko opterećenje.
D
A
C
B
(a)
D
C
B
A
(b)
Slika 5.2 Pravougaoni konačni element (a) i sile u čvorovima
konačnog elementa (b)
Da bi se opisalo ponašanje cijelog tijela pod dejstvom opterećenja, izvest će se prvo
jednačine koje opisuju ponašanje KE kojima je tijelo dis kretizovano pod dejstvom opterećenja. Vez a između sila koje djeluju u vrhovima KE i pomjeranja izvest će se uz pomoć principa o minimumu ukupne potencijalne energije sistema. Radi jednostavnosti, usvojit će se da su dimenzije KE , a poslije će se dobijeni rezultati generalisati za KE proizvoljnih dimenzija. Neka je ishodište koordinatnog sistema u geometri j jskom središtu KE. Da bi se procijenila vrijednost deformacionog rada za konačni element potrebno je
132
pretpostaviti polje pomjeranja na domenu konačnog elementa. Pretpostavit će se da se polje pomjeranja
i
u
i
pravcu mijenja prema sljedećim jednačinama:
Osam konstanti
i
) u prethodnim jednačinama mogu se izraziti preko
(
komponenti vektora pomjeranja ,
(5.1)
,
i
i
u čvorovima KE. Naprimjer, za čvor A vrijedi: i uvrštavanjem ovih vrijednosti u jednačine
(5.1) slijede dvije jednačine,
(5.2)
Na isti način se za tri preostala čvora mogu napisati po dvije jednačine tako da se iz osam jednačina mogu odrediti osam nepoznatih konstanti i . Nakon određivanja ovih osam konstanti komponente vektora pomjeranja pomjeranja u čvorovima na sljedeći način:
i
mogu se na domenu KE izraziti u funkciji
Funkcije
(5.3)
nazivaju se interpolacione funkcije jer služe za interpolaciju
vrijednosti polja pomjeranja unutar KE. Osobina interpolacione funkcije
jednaka nuli u svim čvorovima osim u čvoru koordinate tačke A
jedan, dok u tačkama B
gdje je jednaka jedinici. Naprimjer, za
interpolaciona funkcija je ,C
,iD
je takva da je ima vrijednost
vrijednost ove funkcije je jednaka nuli.
Na osnovu jednačina (5.3) za polja pomjeranja, dilatacije i klizanja su:
(5.4)
(5.5) 133
(5.6)
Potencijalna energija deformacije određena jednačinom (2.26) svodi se za slučaj ravnog naponskog stanja, za koji vrijedi , na sljedeći izraz:
(5.7)
Uvrštavanjem u prethodnu jednačinu izraza za napone, koji su z a slučaj linearno elastičnog tijela i ravnog naponskog stanja određeni jednačinama (2.3 9), potencijalna energija deformacije može se izraziti samo u funkciji dilatacija i klizanja:
(5.8)
U skladu sa jednačinom (2.3 7) izvod potencijalne energije deformacije po pomjeranju jednak je sili koja djeluje u čvoru A,
(5.9)
Na osnovu izraza (5.8) slijedi:
(5.10)
134
Na osnovu izraza (5.4)-(5.6) iz izraza (5.10) slijedi:
Imajući u vidu da je jednačine (5.11) dobija se
, gdje je
(5.11)
debljina ploče (odnosno KE), nakon integracije
(5.12)
Ponavljanjem prethodnog postupka, parcijalni izvod potencijalne energije deformacije po pomjeranju
je
(5.13)
odnosno
(5.14)
Ponavljanjem postupka za preostala tri čvora KE dobija se dodatnih šest jednačina koje se zajedno sa jednačinama (5.9) i (5.14) mogu pisati u matričnom obliku na sljedeći način:
135
C
Neupisani elementi matrice su simetrični u odnosu na glavnu dijagonalu
(5.15)
gdje je
matrica krutosti konačnog elementa u odnosu na definisani
i
koordinatni sistem prikazan na slici 5.2(b), a koja daje vezu između pomjeranja čvorova KE i
sila koje djeluju u tim čvorovima. Primjena pravouganog četverougaonog KE biće prikazana na sljedećem primjeru.
Primjer 5.1
Potrebno je izračunati pomjeranje napadne tačke sile intenziteta
tanku konzolnu ploču debljine Modul elastičnosti materijala ploče
koja djeluje na
čija geometrija i opterećenje su dati na slici 5.3(a). i Poissonov koeficijent
.
Geometrijski domen ploče diskretizovat će se sa 4 pravougaona konačna elementa koji imaju ukupno 9 čvorova, slika 5.3(b). Konačni elementi su međusobno u vezi preko čvorova, a čvorovima 1, 2 i 3 je spriječeno pomjeranje kako bi se time aproksimirao graničn i uslov po pomjeranjima na ukli ještenom dijelu konzole. Da bi se primijenio princip o minimumu ukupne potencijalne energije sistema, potencijalna energija deformacije će se računati kao zbir potencijalnih energija deformacije konačnih elemenata , ,
odakle za parcijalni izvod po pomjeranju
136
vrijedi:
(5.16)
F
F
3
6
2 4m
4m
4
1 1
8
5
2
9
3 4
7
(b)
(a)
Slika 5.3 Dimenzije i opterećenje konzolne ploče (a) i diskretizacija
ploče mrežom konačnih elemenata (b) Svi KE na domenu ploče su jednakih dimenzija kao i KE na slici 5.2 tako da se može iskoristiti matrica krutosti u sistemu jednačina (5.15) . Kada se u sistem jednačina (5.15) uvrst e
vrijednosti debljine ploče, modula elastičnosti i Po issonovog koeficijenta ovaj sistem se može pisati u sljedećem oblik:
5.4
1.95
-3.3
-0.15
-2.7
-1.95
0.6
0.15
1.95
5.4
0.15
0.6
-1.95
-2.7
-0.15
-3.3
-3.3
0.15
5.4
-1.95
0.6
-0.15
-2.7
1.95
-0.15
0.6
-1.95
5.4
0.15
-3.3
1.95
-2.7
-2.7
-1.95
0.6
0.15
5.4
1.95
-3.3
-0.15
-1.95
-2.7
-0.15
-3.3
1.95
5.4
0.15
0.6
0.6
-0.15
-2.7
1.95
-3.3
0.15
5.4
-1.95
0.15
-3.3
1.95
-2.7
-0.15
0.6
-1.95
5.4
(5.17)
Da bi se iskoristio prethodni sistem jednačina potrebno je da se za svaki konačni element
identifikuju odgovarajući brojevi (globalnih) čvorova sa slovima koja označavaju čvorove u prethodnim jednačinama. Naprimjer, na osnovu slike 5.4 za konačni element broj 1 (lokalnom) čvoru A odgovara (globalni) čvor broj 1, čvoru B odgovara čvor 4, itd. Na osnovu veze između globalnih i lokalnih čvorova sistem jednačina (5.17) za KE 1 ima sljedeći oblik:
137
3
6
9 D
2
C
2
5
4 5
2
1
1
8 A
3
1
B
1
4
7
4
Slika 5.4 Veza između lokalnih i globalnih čvorova mreže
5.4
1.95
-3.3
-0.15
-2.7
-1.95
0.6
0.15
1.95
5.4
0.15
0.6
-1.95
-2.7
-0.15
-3.3
-3.3
0.15
5.4
-1.95
0.6
-0.15
-2.7
1.95
-0.15
0.6
-1.95
5.4
0.15
-3.3
1.95
-2.7
-2.7
-1.95
0.6
0.15
5.4
1.95
-3.3
-0.15
-1.95
-2.7
-0.15
-3.3
1.95
5.4
0.15
0.6
0.6
-0.15
-2.7
1.95
-3.3
0.15
5.4
-1.95
0.15
-3.3
1.95
-2.7
-0.15
0.6
-1.95
5.4
(5.18)
Za KE 2 sistem jednačina (5.17) ima oblik:
5.4
1.95
-3.3
-0.15
-2.7
-1.95
0.6
0.15
1.95
5.4
0.15
0.6
-1.95
-2.7
-0.15
-3.3
-3.3
0.15
5.4
-1.95
0.6
-0.15
-2.7
1.95
-0.15
0.6
-1.95
5.4
0.15
-3.3
1.95
-2.7
-2.7
-1.95
0.6
0.15
5.4
1.95
-3.3
-0.15
-1.95
-2.7
-0.15
-3.3
1.95
5.4
0.15
0.6
0.6
-0.15
-2.7
1.95
-3.3
0.15
5.4
-1.95
0.15
-3.3
1.95
-2.7
-0.15
0.6
-1.95
5.4
(5.19)
Za 3 sistem jednačina (5.17) ima oblik:
138
5.4
1.95
-3.3
-0.15
-2.7
-1.95
0.6
0.15
1.95
5.4
0.15
0.6
-1.95
-2.7
-0.15
-3.3
-3.3
0.15
5.4
-1.95
0.6
-0.15
-2.7
1.95
-0.15
0.6
-1.95
5.4
0.15
-3.3
1.95
-2.7
-2.7
-1.95
0.6
0.15
5.4
1.95
-3.3
-0.15
-1.95
-2.7
-0.15
-3.3
1.95
5.4
0.15
0.6
0.6
-0.15
-2.7
1.95
-3.3
0.15
5.4
-1.95
0.15
-3.3
1.95
-2.7
-0.15
0.6
-1.95
5.4
(5.20)
Za KE 4 sistem jednačina (5.17) ima oblik:
5.4
1.95
-3.3
-0.15
-2.7
-1.95
0.6
0.15
1.95
5.4
0.15
0.6
-1.95
-2.7
-0.15
-3.3
-3.3
0.15
5.4
-1.95
0.6
-0.15
-2.7
1.95
-0.15
0.6
-1.95
5.4
0.15
-3.3
1.95
-2.7
-2.7
-1.95
0.6
0.15
5.4
1.95
-3.3
-0.15
-1.95
-2.7
-0.15
-3.3
1.95
5.4
0.15
0.6
0.6
-0.15
-2.7
1.95
-3.3
0.15
5.4
-1.95
0.15
-3.3
1.95
-2.7
-0.15
0.6
-1.95
5.4
(5.21)
Treba primijetiti da su u sistemu jednačina (5.18) do (5.21) u matricama vektor kolonama pomjeranja i opterećenja komponente ovih vektora u globalnom sistemu. Obzirom da su lokalne koordinate i koje su vezane za konačni element paralelne osama i globalnog koordinatnog sistema matrica transformacija je jedinična matric a.
U skladu sa principom o minimumu ukupne potencijalne energije, ili Castiglianovoj teoremi, izvod ukupne potencijalne energije sistema po pomjeranju jednak je sili u pravcu tog
pomjeranja. Za usvojenu diskretizaciju u datom problemu postoji 9 čvorova i u svakom čvoru dva stepena slobode, tako da se mogu pisati za svaki čvor sljedeć i izrazi:
(5.22)
, što ukupno čini 18 jednačina sa 18 nepoznatih komponen ti vektora pomjeranja u 9 čvorova. gdje je
Za formiranje sistema jednačina (5. 22) iskoristit će se matrice krutosti pojedinih KE. Naprimjer, izraz u prvoj jednačini u sistemu jednačina (5.18) sa lijeve strane predstavlja parcijalni izvod potencijalne energije deformacija po pomjeranju
,
. Ovaj izraz čini
doprinos u prvoj jednačini sistema jednačina (5. 22) za . Sistem jednačina (5. 22) pisan u matričnoj formi dat je izrazom (5. 23) a doprinos iz prve od jednačina sistema (5.18) unesen je u izraz (5.23) vodeći računa da elementi matrice krutosti KE 1 koji se množe sa odgovarajućom komponentom vektora pomjeranja dodaju elementima globalne matrice krutosti koji se množe sa istom komponent om vektora pomjeranja u globalnoj matrici krutosti. Naprimjer, element u matrici krutosti u izrazu (5.18) koji se nalazi u trećoj vrsti i sedmoj koloni matrice krutosti iznosi . Izraz sa lijeve strane u trećoj vrsti
sistema jednačina (5.18) odnosi se na parcijalni izvod
i ovaj element treba dodati
139
. Element množi se u sistemu jednačina (5.18) sa komponentom pomjeranja Ovom komponentom se u globalnoj matrici krutosti množe elementi u desetoj vrsti tako da se element dodaje elementu . Na isti način se i ostali elementi matrice krutosti KE 1 sabiraju sa globalnoj matrici krutosti u sedmu vrstu koja se odnosi na jednačinu
.
odgovarajućim elementima globalne matrice krutosti. 5.4
1.95
0.6
0.15
-3.3
-0.2
-2.7
-1.95
1.95
5 .4
-0.2
-3.3
0.15
0 .6
-2
-2.7
0
-0.2
5 .4
-2
-2.7
1.95
-3.3
0.15
0.15
-3.3
-2
5 .4
1.95
-2.7
-0.2
0.6
-3.3
0.15
-2.7
1.95
5.4
-2
0.6
-0.15
-0.2
0 .6
1.95
-2.7
-2
5 .4
0.15
-3.3
-2.7
-2
-3.3
-0.2
0.6
0.15
5 .4
1.95
-2
-2.7
0.15
0.6
-0.2
-3.3
1.95
5.4
=
(5.23)
Sabiranjem matrica krurosti ostala tri KE dobija se sljedeći sistem jednačina:
5.4
1.95
0.6
0.15
0
0
-3.3
-0.15
-2.7
-1.95
0
0
0
0
0
0
0
0
1.95
5.4
-0.15
-3.3
0
0
0.15
0.6
-1.95
-2.7
0
0
0
0
0
0
0
0
0
-0.15
10.8
0
0.6
0.15
- 2.7
1.95
- 6.6
0
- 2.7
-1.95
0
0
0
0
0
0
0.15
-3.3
0
10.8
-0.15
-3.3
1.95
-2.7
0
1.2
-1.95
-2.7
0
0
0
0
0
0
0
0
0.6
-0.15
5.4
-1.95
0
0
-2.7
1.95
-3.3
0.15
0
0
0
0
0
0
0
0
0.15
-3.3
-1.95
5.4
0
0
1.95
-2.7
-0.15
0.6
0
0
0
0
0
0
-3.3
0.15
-2.7
1.95
0
0
10.8
0
1.2
0
0
0
-3.3
-0.15
-2.7
-1.95
0
0
-0.15
0.6
1.95
- 2.7
0
0
0
10.8
0
-6.6
0
0
0 .15
0.6
-1.95
- 2.7
0
0
-2.7
-1.95
-6.6
0
-2.7
1.95
1.2
0
21.6
0
1.2
0
-2.7
1.95
-6.6
0
-2.7
-1.95
-1.95
-2.7
0
1.2
1.95
-2.7
0
-6.6
0
21.6
0
-6.6
1.95
-2.7
0
1.2
-1.95
-2.7
0
0
-2.7
-1.95
-3.3
-0.15
0
0
1.2
0
10.8
0
0
0
-2.7
1.95
-3.3
0.15
0
0
-1.95
-2.7
0.15
0.6
0
0
0
-6.6
0
10.8
0
0
1.95
-2.7
-0.15
0.6
0
0
0
0
0
0
-3.3
0.15
-2.7
1.95
0
0
5.4
-1.95
0.6
-0.15
0
0
0
0
0
0
0
0
-0.15
0.6
1.95
-2.7
0
0
-1.95
5.4
0.15
-3.3
0
0
0
0
0
0
0
0
-2.7
-1.95
-6.6
0
-2.7
1.95
0.6
0.15
10.8
0
0.6
-0.15
0
0
0
0
0
0
-1.95
-2.7
0
1.2
1.95
-2.7
-0.15
-3.3
0
10.8
0.15
-3.3
0
0
0
0
0
0
0
0
-2.7
-1.95
-3.3
-0.15
0
0
0.6
0.15
5.4
1.95
0
0
0
0
0
0
0
0
-1.95
-2.7
0.15
0.6
0
0
-0.15
-3.3
1.95
5.4
=
(5.24) 140
Sistem jednačina (5.24) može se riješiti nakon primjene geometrijskih graničnih uslova.
Čvorovi 1, 2 i 3 nalaze se na mjestu uklještenja konzole tako da su komponente vektora pomjeranja . Reakcije veze u ovim čvorovima, to jest, sile su nepoznate i mogu se odrediti iz prvih šest jednačina nakon što su poznati vektori pomjeranja u čvorovima. Vanjske sile u čvorovima 4 9 su jednake nuli osim komponente sile (slika 5.3). Primjenom pomenutih graničnih uslova preostaje da se riješi sljedeći sistem jednačina (koji je markiran u izrazu (5.24)):
10.8
0
1.2
0
0
0
-3.3
-0.15
-2.7
-1.95
0
0
0
10.8
0
- 6.6
0
0
0.15
0.6
-1.95
-2.7
0
0
1.2
0
21.6
0
1.2
0
-2.7
1.95
-6.6
0
-2.7
-1.95
0
-6.6
0
21.6
0
-6.6
1.95
-2.7
0
1.2
-1.95
-2.7
0
0
1.2
0
10.8
0
0
0
-2.7
1.95
-3.3
0.15
0
0
0
- 6.6
0
10.8
0
0
1.95
- 2.7
-0.15
0.6
-3.3
0.15
-2.7
1.95
0
0
5.4
-1.95
0.6
-0.15
0
0
-0.15
0.6
1.95
-2.7
0
0
-1.95
5.4
0.15
-3.3
0
0
-2.7
-1.95
-6.6
0
-2.7
1.95
0.6
0.15
10.8
0
0.6
-0.15
-1.95
-2.7
0
1.2
1.95
-2.7
-0.15
-3.3
0
10.8
0.15
-3.3
0
0
-2.7
-1.95
-3.3
-0.15
0
0
0.6
0 .15
5.4
1 .95
0
0
-1.95
-2.7
0 .15
0.6
0
0
-0.15
-3.3
1 .95
5.4
=
(5.25)
Rješavanjem prethodnog sistema jednačina dobijaju se komponente vektora pomjeranja u čvorovima:
Reakcije veze u čvorovima 1, 2 i 3 na mjestu uklještenja slijede iz prvih 6 jednačina sistema (5.24):
141
Nakon što su poznate komponente vektora pomjeranja, iz jednačina (5.4) – (5.6) mogu se odrediti komponente deformacija, odnodno naponi iz konstitutivnih relacija.
Usvojena diskretizacija prostornog domena na samo četiri KE nije dovoljna da k valitetno opiše polje pomjeranja i napona u ploči. Da bi se postigli tačniji numerički rezultati potrebno je geometrijski domen ploče podijeliti na znatno veći broj konačnih elemenata. U prethodnom primjeru moglo se vidjeti da se mnoge matematske operacije, kao i postupak formiranja globalne matrice krutosti na osnovu matrica krutosti KE ponavljaju . Ovo čini primjenu MKE pogodnom za programiranje nekim od računarskih programskih jezika.
Posebno je važno da se postupak formiranja osnovnih jednačina za primjenu MKE može iskazati u matričnoj formi što također predstavlja prednost prilikom programiranja . U narednom poglavlju date su osnovne jednačine MKE u matričnoj formi.
5.2 Matrični koncept MKE Izraz (2.28) za ukupnu potencijalnu energiju sistema može se napisati u sljedećem matričnom obliku:
gdje
su
matrica
komponenti
matrica komponenti deformacija,
komponenti vektora zapreminskih sila,
površinskih sila, i
(5.26)
napona, matrica
matrica komponenti vektora
matrica komponenti vektora pomjeranja.
Konstitutivne relacije (2.16) mogu se izraziti u matričnom obliku
gdje je
142
matrica elastičnosti
(5.27)
(5.28)
simetrično
Imajući u vidu da je transponovan proizvod dvije matrice jednak proizvodu transponovanih , i da za simetričnu matrica koje se množe obrnutim redom, to jest, matricu vrijedi relacija
, izraz (5.26) može se pisati u sljedećem obliku:
Komponente matrice deformacija
(5.29)
definisane su preko komponenti vektora pomjeranja
jednačinama (2.4) – (2.7) i mogu se u matričnom obliku pisati na sljedeći način:
gdje je
matrica operatora:
(5.30)
(5.31)
Komponente vektora pomjeranja na domenu KE opisuju se izrazom
gdje su:
(5.32)
matrica interpolacionih funkcija
143
sa komponentama
(5.33)
matrica kolona diskretnih vrijednosti zavisnih varijabli
,
(stepeni slobode) u čvorovima KE
(5.34)
je broj čvorova KE. Naprimjer, za ravanski četverougaoni KE sa dva stepena slobode . Primjer komponenti matrice kretanja po čvoru interpolacionih funkcija za dvodimenzionalni problem može se vidjeti u sistemu a
jednačina (5.3) gdje je, naprimjer,
.
, ili
Uvršavanjem izraza (5.30) i (5.32) u jednačinu (5.29) dobija se sljedeći izraz:
(5.35)
(5.36)
odnosno
gdje su
matrica deformacija,
,a
matrica krutosti KE,
(5.37)
Minimiziranje izraza (5.36) u odnosu na komponente vektora pomjeranja u čvorovima slijedi:
144
(5.38)
gdje se matrica vektor kolona komponenti sila u čvorovima KE računa iz izraza
(5.39)
5.3 Ravanski četverougaoni linearni KE Na slici 5.1(b) prikazana je diskretizacija geometrije ravanskog problema sa pravougaonim
KE. Očigledno je da sa ovakvim elementima nije dobro aproksimirana geometrija. Jedan od načina da se bolje aproksimira geometrija je da se smanji veličina KE. Efikasniji pristup je da se koriste ravanski četverougaoni KE kao što je to prikazano na slici 5.5 (b) gdje je očigledno da je sa ovim KE mnogo bolja aproksimacija geometrije problema.
F 1
F 1
F 2
(a)
F 2
(b)
Slika 5.5 Tanka ravna ploča pod dejstvom opterećenja (a) i podjela
geometrijskog domena na mrežu četverougaonih KE sa četiri čvora
Na slici 5.6 prikazan je četverougaoni konačni element sa četiri čvora. Jedan od načina da se izvede matrica krutosti ovog elementa je da se izvrši preslikavanje geomtrije KE iz ravni na ravan u jednostavniji oblik. Preslikavanje određeno jednačinama
(5.40)
145
preslikava geometrijski domen četverougaonog KE sa četiri čvora u ravni na pravougaoni (kvadratni) element sa četiri čvora u ravni sa granicama domena ograničenim pravama
i
, slika 5.6. Koordinate
i
3
4
nazivaju se i prirodnim koordinatama.
4
3
1
2
1
2
Slika 5.6 Preslikavanje geometrijskog domena četverougaonog elementa na pravougaoni element
Funkcije
(5.41)
imaju osobinu da su jednake jedinici u jednom od čvorova elementa, dok su u ostalim čvorovima jednake nuli. Polje pomjeranja na domenu KE može se opisati istim interpolacion im funkcijama:
146
(5.42)
Da bi se izračunala matrica krutosti
data izrazom (5.37) potrebno je izvesti integraciju izraza u
prirodnim koordinatama. Za dvodimezionalni KE debljine
gdje je
vrijedi relacija
determinanta Jacobianove matrice
Izraz (5.37) sada se može pisati :
Da bi se mogla izvesti integracija u izrazu (5.44) potrebno je matricu deformacija izraziti u funkciji prirodnih koordinata interpolacionih funkcija po
,
(5.43)
(5.44)
. U matrici deformacija se nalaze parcijalni izvodi
koordinatama (vidjeti izraz 5.31). Parcijalni izvodi
interpolacionih funkcija se mogu izraziti u funkciji prirodnih koordinata na sljede ći način:
gdje se iz sistema jednačina (5.45) mogu izračunati
(5.45)
i .
Podintegralne funkcije u jednačini (5.44) mogu biti komplikovane te se iz toga razloga vrši numerička integracija u ovom jednačinama od koje je najzustupljenija Gaussova metoda integracije (ili Gaussove kvadraturne formule).
Prilikom izvođenja matrice krutosti elementa korištene su iste interpolacione funkcije za preslikavanje geometrijskog domena KE i opis polja pomjeranja na domenu KE. Ovakvi elementi kod kojih se koriste iste interpolacione funkcije za preskivanja geometrijskog domena KE i opis polja pomjeranja na KE zovu se izoparametarski KE.
Ravanski trougaoni KE sa tri čvora , koji se često nalazi u standardnim bibliotekama MKE softverskih paketa, prikazan je na slici 5.7. Polje pomjeranja za ovaj element opisano je
sljedećim linearnim funkcijama :
(5.46)
147
3
Dilatacije i klizanja na domenu ovog
elementa su konstantni,
1
, i
Dakle,
2
za
ovaj
element
, .
postoji
diskontinuitet komponenti deformacija na granicama KE. Trougaoni KE sa tri čvora daje u opštem slučaju lošije rezultate od
četverougaonog elementa sa četiri čvora i
Slika 5.7 Ravanski KE sa tri čvora
rijetko se koristi.
Da bi se što kvalitetnije opisala deformacija duž KE uvode se i KE kod kojih je promjena deformacije u elementu kvadratna, kubna ili još višeg reda. Kod ovih elemenata red polinoma kojim se aproksimira polje pomjeranja takođ e raste, a time i broj stepeni slobode kretanja elementa, odnosno matrica krutosti elementa. U narednom poglavlju dat je opis
načešće korištenih ravanskih KE višeg reda.
5.4 Ravanski konačni elementi višeg reda U slučaju kada su granice geometrijskog domena problema koji se rješava zakrivljene linije ili površine, KE ograničen pravim linijama ne može da dobro opiše geometriju problema. Iz tog razloga poželjno je imati i KE koji su sposobni da dobro prate krivolinijske konture domena. Jedan od takvih KE sa devet čvorova koji pored čvorova u vrhovima ima i čvorove na sredinama strana i jedan čvor u sredini elementa prikazan je na slici 5.8.
4
3
7
8
9
1
6
2
5
7
3
4
6
9
8 1
2 5
Slika 5.8 Preslikavanje geomtrijskog domena četverougaonog KE
sa devet čvorova na pravougaoni domen 148
Konačni element ima u svakom čvoru po dva translatorna stepena slobode kretanja, što znači da ima ukupno 18 stepeni slobode. Interpolacione funkcije za KE sa devet čvorova su
(5.47)
dok je geometrija elementa opisana preko prirodnih koordinata jednačinama
(5.48)
odnosno, polje pomjeranja je za slučaj izoparametarskog elementa op isano pomoću istih funkcija
gdje su
i
koordinate čvorova u globalnom koordinatnom sistemu, a
čvorova u pravcu osa globalnog koordinatnog sistema.
(5.49)
i
pomjeranja
Važno je primijetiti da ovaj KE osim što bolje aproksimira krivolinijske granice domena aproksimira i polje pomjeranja polinomima višeg reda što daje prednost u odnosu na četverougaoni KE sa četiri čvora. Ovaj element, kao i ostali KE sa interpolacionim funkcijama višeg reda, daje tačnije rezultate od line arnih elemenata. Ako bi se tražio nedostatak ovih elemenata u odnosu na linea rne elemente to je što su matrice krutosti mnogo veće kod ovakvih elemenata tako da u opštem slučaju zahtjevaju i više računarskog vremena. Ipak , u opštem slučaju, ovi elementi su mnogo efektivniji od linearnih elemenata. Elementi višeg reda koji su često u upotrebi su i pravougaoni KE sa 8 čvorova i trougaoni KE sa šest čvorova, slika 5.9. KE sa osam čvorova je posebno pogodan u slučaju kada su prisutne deformacije usljed savijanja. 149
4
7
3
5
3
2
8
6
6
4 1
2 5
1 (b)
(a)
Slika 5.9 Ravanski četverougaoni KE sa osam (a) i trougaoni KE sa
šest čvorova (b) Na sljedećem primjeru može se procijeniti i tačnost pojedinih KE.
Primjer 5.2
U primjeru je analizirana brzina konvergencije pojedinih KE ka tačn om rješenju problema datog na slici 5.10. Problem je analiziran kao model ravnog naponskog stanja sa trougaonim
KE sa tri čvora, trougaonim KE sa šest čvorova, i četverougaonim KE sa 8 čvorova .
8m
2m
debljina ploče je 0.02 m
Slika 5.10 Konzolna tanka ploča
Na slici 5.11(a) prikazan je numerički proračun (ADINA softverom) maksimalnog pomjeranja
tačke na slobodnom kraju ploče za uniformne mreže različitih stepeni slobode kretanja i različite vrste KE. Na slici se vidi da oba KE višeg reda , trougaoni KE sa 6 i čeverougaoni KE sa 8 čvorova, daju veoma dobre rezultate i sa grubim mrežama. Trougaoni KE sa tri čvora pokazuje veoma sporu konvergenciju prema tačnom rješenju. Na slici 5.11(b) date su vrijedno sti greške numeričkog proračuna maksimalnog pomjeranja za
mreže sa različitim brojem stepeni slobode kretanja. Oba elementa višeg reda daju grešku manju od 2.5% već za mreže sa 20 čvorova (40 stepeni slobode kretanja), dok trougaoni KE za isti broj čvorova pravi grešku oko 45%, a potrebno je preko 500 čvorova da bi greška pala ispod 2.5%. 150
3.00
50
2.50
40 35
2.00
) m (
4 -
Trougaoni KE sa 3 čvora Trougaoni KE sa 6 čvorova pravougaoni KE sa 8 čvorova
45
) 30 % ( a25 k š 20 e r G15
1.50
0 1 v
1.00 Trougaoni KE sa 3 čvora
10
Trougaoni KE sa 6 čvorova
0.50
Četverougaoni KE sa 8 čvorova
0.00 0
100 Broj stepeni slobode
200
5 0 0
40
80
120
160
200
Broj stepeni slobode
(b)
(a)
Slika 5.11 Numerički rezultati za pomjeranje (a) i greška
numeričkog proračuna (b) za različite KE
5.5 Osnosimetrični konačni elementi U poglavlju 2.8.3 date su jednačinama (2.43) veze između komponentnih deformacija i komponenti vektora pomjeranja i jednačinama (2.44) konstitutivne relacije za slučaj
osnosimetričnog problema.
O
(a)
(b)
Osnosimetrični problem sa mrežom KE (a) i osnosimetrični KE (b) Slika
5.12
151
Za izvođenje matrice krutosti osnosimetričnih KE potrebno je pretpostaviti polje pomjeranja na domenu KE. Polje pomjeranja je funkcija
i
koordinata,
i
izražava se u funkciji stepeni slobode pomjeranja u čvorovima KE. Naprimjer, za slučaj pravougaonih KE na slici 5.12(a) i 5.12(b) u svakom čvoru KE ima dva stepena slobode (naprimjer, za čvor broj 3 stepeni slobode su i ). Konačni elementi koji se koriste za problem ravnog naponskog i ravnog deformacionog sta nja koriste se i za osnosimetrične probleme. Matrica krutosti KE formira se na osnovu izraza (5.37) gdje za osnosimetrične KE vrijedi
:
Za osnosimetrične KE također postoji veličina sila u čvorovima:
(5.50)
u izrazu (5.39) za vektor kolone komponenti
(5.51)
tako da se ova veličina obično izostavlja pri formiranju matrice krutosti i matrice vektor kolone opterećenja jer se u izrazu (5.38) krati.
Primjer 5.3
Cilindar na slici 5.13(a) nalazi se pod unutrašnjim pritiskom od
vanjski poluprečnik cilindra je Ravnima cilindra
i
i
. Unutrašnji i
. Visina cilindra je
spriječeno je kretanje u pravcu
.
ose. Materijal cilindra ima
modul elastičnosti i Poissonov koeficijent . Potrebno je numerički izračunati polje pomjeranja i napon i uporediti numeričke rezultate sa analitičkim rješenjem. Analitičko rješenje za polje pomjeranja i napona je (Timoshenko i Goodier , 1970):
152
(5.52)
O
(a)
(b)
(c)
Slika 5.13 Cilindar pod pritiskom (a), model osnosimetričnog problema (b) i model ravnog deformacionog stanja (c)
Na slici 5.13(b) prikazan je osnosimetričan model problema. Simbolima klizača označeni su granični uslovi na mjestima gdje je spriječeno po mjeranje u pravcu ose. Domen je podijeljen u radijalnom pravcu sa 5 osnosimetričnih KE sa osam čvrorova. Podjela dom ena u
pravcu je nebitna jer su u tom pravcu sve varijable konstantne.
Problem se može riješiti i kao problem ravnog deformacionog stanja sa geometrijskim domenom kao na slici 5.13(c), jer su dilatacije u pravcu
ose jednake nuli, slika 5.13(c).
Obzirom na simetriju problema modelirana je samo četvrtina domena (slika 5.13(c)) , mada je bilo moguće modelirati i bilo koji isječak domen a. Ovaj model imao je u radijalnom i cirkularnom pravcu po 5 KE sa osam čvorova. Na slikama 5.14(a) do 5.14(d) prikazani su uporedo numerički i analitički rezultati proračuna
polja radijalnog pomjeranja (cirkularnog) napona
, normalnog napona u radijalnom pravcu
, i normalnog napona u pravcu ose simetrije
, normalnog
. Na slikama se vidi
da oba numerička modela daju sa izabranom gustoć om mreže gotovo identične rezultate visoke tačnosti.
153
r
0.1
1.2
0.12
(m)
0.14
0.16
0.18
0.2
0 1 -2
) a P -4 M (
0.8
) m ( r
r
0.6
u
5
0 10.4
-8
Analitičko rješenje MKE osnosimetrični KE
0.2
0.14 r
0.16
0.18
MKE osnosimetrični KE MKE ravno def. stanje
0 0.12
Analitičko rješenje
-10
MKE ravno def. stanje
0.1
-6
-12
0.2
(m)
(a)
(b)
18
4
16 14
3
12
) a P 10 M (
8
) a P M (
6
2
z
Analitičko rješenje MKE osnosimetrični KE
4 2
1
Analitičko rješenje MKE osnosimetrični KE
MKE ravno def. stanje
MKE ravno def. stanje
0
0 0.1
0.12
0.14 0.16 r (m)
0.18
0.2
0.1
0.12
0.14
0.16 r (m)
0.18
0.2
(d)
(c)
Slika 5.14 Numerički i analitički rezultati proračuna radijalnog pomjeranja (a), normalnog napona u radijalnom pravcu (b), normalnog (cirkularnog) napona (c), i normalnog napona u pravcu ose simetrije (d)
5.6 Trodimenzionalni konačni elementi Na slici 5.15 prikazani su neki od standardnih trodimenzionalnih KE. Na slikama 5.15(a) i 5.15(b) su heksaedarski KE sa 8 i 20 čvorova, a tetraedarski KE sa 4 i 10 čvorova su prikazani na slikama 5.15(c) i 5.15(d). Elementi imaju tri stepena slobode kretanja po čvoru što znači 154
da, naprimjer, heksaedarski element sa 20 čvorova ima 60 stepeni slobode kretanja koliki je i red matrice krutosti ovog elementa.
(a)
(c)
(b)
(d)
Slika 5.15 Heksaedarski KE sa 8 čvorova (a) i 20 čvorova (b), i
tetraedarski KE sa 4 čvora (c) i deset čvorova (d)
Za izoparametarski KE sa 8 čvorova (slika 5.1 6) interpolacione funkcije izražene preko prirodnih koordinata , i se mogu pisati u sljedećem obliku:
(5.53)
a za opis geometrije KE i polja pomjeranja vrijede relacije:
155
(5.54)
(5.55)
2 3
3 1 6
7
7 5
8
1
4
4
2
8
6 5
Slika 5.16 Preslikavanje geometrijskog domena heksaedarskog KE
sa osam čvorova i oznake čvorova Kod primjene trodimenzionalnih KE vrijedi isto pravilo kao i kod ravanskih KE. Konačni elementi višeg reda su tačniji od linearnih KE. Konačni element sa 27 i 20 čvorova daj e rezultate visoke tačnosti mada analiza ovim elementom zahtijeva i veće računarsko vrijeme. Ovi elementi daju također i najbolje rezultate u slučaju kada su pravougaoni . Heksaedarski element sa 20 čvorova se preporučuje i za slučaj trodimenzionalnih tijela sa stijenkama izloženih savijanju. Standardni heksaedarski KE sa 8 čvorova i tetraedarski KE sa 4 čvora se ne preporučuju za dijelove konstrukcije gdje su dominantni efekti savija nja.
5.7 Kriteriji konvergencije MKE rješenja Kao i kod svakog numeričkog postupka postavlja se pitanje, da li s povećanjem stepena diskretizacije, to jest, u ovom slučaju sa smanjivanjem veličine KE, numeričko rješenje teži
156
tačnom (analitičkom) rješenju. Konvergencija rješenja zavisi od vrste konačnog elementa koji se koristi, odnosno od polinoma kojim se aproksimira polje nezavisno promjenjive.
Da bi rješenje varijacionom fomulacijom MKE bilo konvergentno dovoljno je da bud u ispunjeni tzv. kriteriji kompatibilnosti i kompletnosti. U primjeni su i neki KE koji ne zadovoljavaju pomenute kriterije konvergencije. Iako konvergencija s njima nije
zagarantovana često se postižu rezultati visoke tačnosti, koji brzo konvergiraju (obično ne monotono) prema analitičkom rješenju. Kriterij kompatibilnosti podrazumijeva da polje zavisne varijable, kao i izvodi ove varijable do
za jedan red manji od izvoda varijable koji se pojavljuju u integralnoj formulaciji jednačina elementa, unutar e lementa i duž njegovih granica koje dijeli sa susjednim elementima moraju biti neprekidni. Prilkom izvođenja jednačina ravnoteže za konačne elemente zavisne varijable su bil e
komponentna pomjeranja, a najveći izvodi komponenti pomjeranja bili su prvi parcijalni izvodi ovih pomjeranja koji su se javljali u komponentama matrice deformacije (dilatacije i klizanja). Prema tome, u skladu sa kriterijom kompatibilnosti potrebno je da samo polje pomjeranja unutar domena elementa i na njegovim granicam bude neprekidno.
Prilikom izvođenja jednačina konačnih elemenata za štap i gredu i formiranja matrice krutosti sistema, usvojeno je da dva susjedna konačna elementa d i jele zajednički čvor tako da je uslov kompatibilnosti bio automatski zadovoljen.
Za četverougaoni KE prikazan na slici 5.6 polje pomjeranja između čvorova 1 i 2 za : komponentu dobije se na osnovu prve od jednačina (5.42) za
(5.56)
Iz prethodne jednačine se vidi da se polje pomjeranja na granici KE određenoj sa čvorovima 1 i 2 mijenja linearno od vrijednosti
do
. Obzirom da je polje pomjeranja i susjednog KE
koji dijeli čvorove 1 i 2 opisano istom funkcijom uslov kompatibilnos ti je automatski zadovoljen. Kriterij kompletnosti podrazumijeva da pretpostavljeno polje zavisne varijable kao i izvodi ove varijable do reda jednakog izvodu varijable koji se pojavljuju u integralnoj formulaciji
jednačina elementa mora ju osigurati i mogućnost da ova varijabla i njezini izvodi budu konstantni na elementu kada veličina elementa teži nuli. U slučaju pomjeranja kao zavisne varijable , kriterij kompletnosti se svodi na zahtjev da pretpostavljeno polje pomjeranja mora osigurati da se ne događa deformacija elementa ako tijelo ne doživljava deformaciju pri kretanj u (kreće se kao kruto tijelo), i da omogućuje konstantnu deformaciju na polju KE. Naprimjer, kriterij kompletnosti za slučaj ravanskog 157
trougaonog KE prikazanog na slici 5.7 čije polje pomjeranj a je određeno jednačinama (5.46), je ispunjen jer koeficijenti i osiguravaju da se može opisati kretanje elementa kao krutog tijela, dok koeficijenti i osiguravaju mogućnost konstantne deformacije
(
158
.
