Content extract
TERVEZÉS-INFORMATIKAI FÜZETEK SZABÓ TAMÁS ALKALMAZOTT MECHANIKA VÉGESELEM-MÓDSZER ALKALMAZÁSAI MISKOLCI EGYETEM 2003 2 Készült az Oktatási Minisztérium megbízásából. Szerződéskötő Hatóság: VÁTI Területfejlesztési Igazgatóság (NARD) Főprojekt címe: Phare HU0008-02 ESZA-típusú kísérleti projekt a képzésből a munka világába történő átmenet támogatására Alprojekt címe: Phare HU0008-02-04 A felnőttoktatás és az élethosszíg tartó tanulás lehetőségeinek javítása Projekt címe: Phare HU0008-02-04-0005 Moduláris Tervezés-informatika tanfolyam műszakiaknak. Projektvezető: Dr. Takács György egyetemi docens Lektorálta: Dr. Nándori Frigyes egyetemi docens Miskolc-Egyetemváros, 2003 november. 3 4 TARTALOMJEGYZÉK 1 RUGALMASSÁGTANI ISMERETEK ÖSSZEFOGLALÁSA 1D-S FELADAT ESETÉN 7 1.1 1.2 1.3 1.4 Rugalmas peremérték-feladat és megoldása A virtuális munkaelv variációs alakja egydimenziós esetben
Teljes potenciális energia minimum elve Ritz-módszer 1D-s feladatra 8 14 16 19 2 LOKÁLIS APPROXIMÁCIÓ ELVE, VÉGESELEMES DISZKRETIZÁCIÓ EGYDIMENZIÓS FELADATRA 25 2.1 2.11 2.12 2.2 Húzott-nyomott rúdelem Szerkezeti mátrixok A csomóponti elmozdulások meghatározása. Rácsos szerkezet vizsgálata húzott-nyomott rúdelemekkel 3 IZOPARAMETRIKUS ELEMCSALÁD 3.1 3.2 3.21 3.22 3.23 3.24 3.25 3.26 3.27 28 31 33 37 43 1D-s húzott-nyomott rúdelem 45 Síkbeli feladatok 49 Általánosított síkfeszültségi állapot( ÁSF, Tárcsafeladat, Plane stress problem): 49 Síkalakváltozási állapot (Plane strain problem): 50 Forgás v. tengelyszimmetrikus állapot (Axi-symmetric problem): 51 Rugalmas síkbeli (ÁSF) peremérték feladat kitűzése 52 A teljes potenciális energia 54 Síkbeli négy csomópontú izoparametrikus elem 55 Példa torzult síkbeli elem elfajuló leképzésére: 60 4 DINAMIKAI FELADATOK VIZSGÁLATA VÉGESELEM-MÓDSZERREL 62 4.1 4.2 4.3 4.4 Végeselem
mátrixok és a tehervektor: Diszkretizált mozgásegyenlet Állandósult harmonikus gerjesztés Sajátrezgések 64 66 67 68 5 HŐFESZÜLTSÉGEK SZÁMÍTÁSA 70 5 5.1 5.2 Hőmérséklet hatásának figyelembevétele 1D-s feladatnál: Hőmérséklet hatásának figyelembevétele 2D-s feladatnál: 71 73 6 AZ IDEAS PROGRAMRENDSZER ALKALMAZÁSA VÉGESELEMES FELADATOK MEGOLDÁSÁRA 75 6.1 6.11 6.2 6.21 6.3 6.31 76 77 78 79 80 81 Rúdszerkezet vizsgálata Rúdfeladat megoldásához alkalmazott Ideas utasitások Tárcsafeladat Tárcsafeladat megoldásához alkalmazott Ideas utasitások Térbeli feladat Térbeli feladat megoldásához alkalmazott Ideas utasitások 6 1 RUGALMASSÁGTANI ISMERETEK ÖSSZEFOGLALÁSA 1D-S FELADAT ESETÉN A Végeselem-módszer alkalmazásai előadás vázlat két tankönyv Páczelt I.: VÉGESELEM MÓDSZER A MÉRNÖKI GYAKORLATBAN, I. KÖTET, (MISKOLCI EGYETEMI KIADÓ, 1999.) és Szabó B, Babuska I,: Finite Elelement Analysis, (John Wiley
& Sons, INC. New York, 1991) valamint a szerző oktatási és kutatási tapasztalatai alapján készült. 7 1.1 RUGALMAS PEREMÉRTÉK-FELADAT ÉS MEGOLDÁSA Az 1.a ábrán egy l hosszúságú és A keresztmetszetű prizmatikus rúd látható A rúd rúdirányú önsúlyával és a véglapon megoszló erővel terhelt A térfogaton egyenletesen megoszló ( ρ g ) gravitációs erő és az egyik véglapját normálirányú ( f l ) megoszló erő rúdirányba terheli. A másik véglap befalazott Az ilyen típusú problémát húzott-nyomott rúdfeladatnak is szokás nevezni. A rúd anyaga homogén, izotróp és lineárisan rugalmas (E a Young féle modulus). A feladat megoldása során feltételezzük, hogy a rúd tetszőleges keresztmetszetében csak rúdirányú egyenletesen megoszló normálfeszültség ébred Az ilyen feladatot mechanikai szempontból egydimenziós feladatnak tekintjük, az eredeti feladattal egyenértékű egydimenziós rúdmodellt az 1.b ábra szemlélteti
dV=A dz Au Ap Ap r f r ρg Au E, A r F r r p = Aρg z dz z dz l l 1.a ábra A húzott-nyomott rúd térbeli rajza 1b ábra A rúd feladat egydimenziós modellje Az ábrákon alkalmazott jelölések: A – a keresztmetszet területe, ρ - a sűrűség, g – a rúdirányú gravitációs gyorsulás, r r p = p z e z - vonal menti terhelés intenzitás, r r f = f z e z - felületen megoszló terhelés intenzitása, r r r F = Fz e z = Af - véglapot terhelő megoszló erő eredője, l – a rúd hossza, dV – az elemi rúdtérfogat, dz - az elemi rúdhossz, 8 Ap - a terhelt felület, dinamikai perem, Au - az elmozdulás előírását tartalmazó felület, kinematikai perem. A vázolt feladat esetén keressük a w(z) elmozdulást, mint a hely függvényét. Mivel az elmozdulás függvény a rúd összes pontjának elmozdulását magába foglalja, ezért szokás elmozdulási mezőnek is nevezni. Az elmozdulás meghatározásához a szilárdságtanban tanult egyenleteket
alkalmazzuk. Jelen esetben a rugalmas peremérték feladat egyenletei: Kinematikai (geometriai) egyenlet: εz = dw dz 0< z<l, (1) ahol ε z a z irányú fajlagos nyúlás. Anyagegyenlet (egyszerű Hooke-törvény): σ z = Eε z N = A σ z = AEε z ⇒ 0< z<l, (2) ahol σ z a rúdirányú normálfeszültség, N a rúderő. Egyensúlyi egyenlet: Az egyenlet származtatásához tekintsük a dz elemi hosszúságú rúd egyensúlyát r pz N(z) N(z+dz) dz 2. ábra: Az elemi hosszúságú rúdelem egyensúlya Egyensúly esetén a z irányú vetületi egyenlet: − N ( z ) + N ( z + dz ) + p z dz = 0 . Az N(z+dz) rúderőt sorba fejtve N ( z + dz ) = N ( z ) + dN dz + magasabb rendű dz tagok és visszahelyettesítve a vetületi egyenletbe kapjuk az egyensúlyi egyenletet dN + pz = 0 dz 9 0< z<l. (3) Az (1)-(3) egyenletek tulajdonképpen egy differenciál-egyenletrendszert alkotnak, amelyhez két peremfeltétel írható fel esetünkben.
Kinematikai peremfeltétel: Az 1.a-b ábrák alapján láthatjuk, hogy a befalazásnál az Au felületen az elmozdulás adott, értéke: w(0)=0. (4) Dinamikai peremfeltétel: Az 1.a-b ábrák alapján láthatjuk, hogy a szabad rúdvégen az Ap felületen a terhelés adott, így ott a rúderő ismert, értéke: N(l)= Fz . (5) Az (1)-(5) rugalmas peremérték feladatban tulajdonképpen három ismeretlen menynyiség ( w, ε z , N ) fordul elő. A 3 ismeretlen meghatározásához rendelkezésünkre áll 3 mezőegyenlet ((1) – (3)) és a kinematikai valamint a dinamikai peremfeltétel ((4), (5)) A rugalmas peremérték feladat ezen alakját a probléma erős alakjának is szokás nevezni, mert pontszerűen került megfogalmazásra. Az (1)-(5)-el megadott peremérték-feladat analitikus megoldását szokás tényleges megoldásnak is nevezni. (Megjegyezzük, hogy egy általános térbeli rugalmas peremérték-feladat esetén az ismeretlenek száma 15 és a skalár egyenletek száma
szintén 15, valamint rendelkezésünkre áll még a kinematikai és a dinamikai peremfeltétel vektorikus alakban.) Tényleges megoldás: Az (1)-(5)-tel adott rugalmas permérték feladat megoldása ebben az egyszerű esetben analitikusan is előállítható. Helyettesítsük (1)-t (2)-be N = AE dw . dz Ezután a (6) kifejezést a (3)-as egyensúlyi egyenletbe helyettesítjük 10 (6) dw d AE dz + pz = 0 . dz (7) Mivel a rúd prizmatikus az AE konstansok kiemelhetők zárójel elé és kapjuk az elmozdulásra vonatkozó alapegyenletet AE d 2w + pz = 0 . dz 2 (8) A (8) egyenlet átrendezése után kapott egyenletet integráljuk kétszer egymás után z- szerint p d 2w =− z 2 dz AE / ∫ dz p dw = − z z + C1 dz AE w( z ) = − (9) / ∫ dz (10) pz 2 z + C1 z + C 2 2 AE (11) A (11) megoldásban szereplő C1 ,C 2 konstansokat a (4) és (5) peremfeltételekből határozhatjuk meg. Vagyis a (4)-s kinematikai peremfeltételt
segítségével megkapjuk a második konstanst: w(0 ) = 0 = − pz 0 + C1 0 + C2 2 AE ⇒ C2 = 0 . Az (5)-s dinamikai peremfeltétel szolgáltatja az első konstanst (10) felhasználásával N (l ) = AE dw (l ) = AE − p z l + C1 = Fz . dz AE (11) A (11)-ből az első konstans: C1 = Fz + p z l . AE Végül az analitikus vagy tényleges megoldás: elmozdulás rúderő w( z ) = − p z 2 Fz + p z l z, z + 2 AE AE N (z ) = p z (l − z ) + Fz . (12) (13) Habár ebben az egyszerű esetben ismert a (12) és (13) analitikus megoldás, de egy bonyolult térbeli feladat esetén általában nem tudjuk előállítani. Ekkor csak közelítő meg- 11 oldással tudunk szolgálni. A továbbiakban majd kitérünk ezen egyszerű egydimenziós feladathoz történő közelítő megoldás előállítására Közelítő megoldással szemben bizonyos elvárásokkal élünk mind az elmozdulás és mind a feszültség vonatkozásában. Ezért bevezetünk két
definíciót Kinematikailag lehetséges elmozdulásmező ( w* ( z ) ) Definició: Kinematikailag lehetséges az az elmozdulásmező ( w* ( z ) ), amely folytonos és elegendően sokszor differenciálható, valamint kielégíti a kinematikai peremfeltételt. Esetünkben a definíció tömören: ε* = dw* , dz w* (0) = 0 . Elegendő az egyszeri differenciálhatóságot megkívánni. Nyilvánvaló, hogy a tényleges megoldás eleget tesz a kinematikai lehetséges elmozdulásmező definíciónak. Statikailag lehetséges feszültségmező Definició: Statikailag lehetséges az a feszültségmező, amely kielégíti az egyensúlyi egyenletet és a dinamikai peremfeltételt. Esetünkben a feszültségi mező az N rúd erőnek felel meg, a definíció pedig tömören az alábbiak szerint írható: dN + pz = 0 , dz N (l)= Fz . Az olyan közelítőmegoldást előállító módszert, amely a kinematikailag lehetséges elmozdulásra alapozott elmozdulási módszernek nevezzük. Ha a
közelítő módszer a statikailag lehetséges feszültségre alapozva állítja elő a megoldást, akkor erőmódszernek nevezzük A gyakorlatban legelterjedtebb az elmozdulás módszerre alapozott eljárás Az erőmódszer alkalmazása általában lényegesen nehezebb, mint az elmozdulás módszer. Az utóbbi években egyre több kutatás foglalkozik a két mező együttes alkalmazásával, amelyet a vegyes mezők módszerének nevezzük. Jelen kurzus során az elmozdulási módszerrel fogunk foglalkozni. A közelítő megoldás előállításához olyan hiba elvre van szükségünk, amely egy adott függvénytérből a 12 lehető legjobb közelítéssel szolgál a megoldás vonatkozásában. A bemutatásra kerülő két hibaelv a pontszerű megfogalmazással szemben integrál értelmű megfogalmazást alkalmaz. Ezen két bemutatásra kerülő elv a virtuális munka elv variációs alakja és a potenciális energia minimum elv. 13 1.2 A VIRTUÁLIS MUNKAELV VARIÁCIÓS
ALAKJA EGYDIMENZIÓS ESETBEN A virtuális munkaelv származtatása előtt szükségünk van az esetünkre vonatkozó két definíció megfogalmazására. Az alábbi definíciók általános esetre megtalálhatók pl Béda Gy.-Kozák I Rugalmas testek mechanikája című könyvben is Virtuális elmozdulásmező: Két különböző kinematikailag lehetséges elmozdulásmező különbségét virtuális elmozdulásnak nevezzük, és δw( z ) módon jelöljük ( δw( z ) = w1* (z ) − w2 ( z ) ). Elmozdulásmező variációja: Ha a virtuális elmozdulásmezőt úgy számítjuk, mint a tény- leges elmozdulásmező és annak elég kis környezetében lévő kinematikailag lehetséges elmozdulásmező különbségét, akkor ezt az elmozdulásmező variációjának nevezzük és δw( z ) módon jelöljük . ( δw(z ) = w1* ( z ) − w( z ) , lásd a 3. ábrát ) Megjegyzés: A virtuális elmozdulásmező és az elmozdulás variációja megegyezik, ha a tényleges mező közvetlen
környezetében a tényleges elmozdulásmezőre vonatkozó virtuális elmozdulásmező elegendően kicsiny. w tényleges w( z ) δw( z ) kinematikailag lehetséges közelítő w* ( z ) l z 3. ábra: Az elmozdulásmező variációja 14 Nyilvánvaló, hogy az elmozdulásmező variációja az egy tetszőleges folytonos és elegendően sokszor (esetünkben egyszer) differenciálható függvény és a kinematikai peremen az értéke zérus. Tömören írva: δε = d (δw) , dz δw(0 ) = 0 . Szorozzuk meg a (3) egyensúlyi egyenletet δw -el és integráljuk a 0-tól l-ig: l ∫ δw ⋅ 0 l dN dz + ∫ δw ⋅ p z dz = 0 dz 0 u u = (14) v d (δw) dz v=N A parciális integrálást elvégezve (vesszőtlenek szorzata mínusz a kiszámítottak szorzatának integrálja) kapjuk, hogy d (δw) δw ⋅ N | − ∫ ⋅ N dz + ∫ δw ⋅ p z dz = 0 dz 0 0 l l l 0 (15) Figyelembe véve, hogy δw(0 ) = 0 és δw(l ) = δwl , valamint alkalmazva az (5) dinamikai peremfeltételt,
átrendezés után kapjuk a virtuális munka elv variációs alakját: d (δw) ∫0 dz ⋅ N dz = ∫0 δw ⋅ p z dz + Fz ⋅ δwl l l (16) A (16) virtuális munka elv variációs alakja azt fejezi ki, hogy a tényleges megoldásnál a belső erők virtuális munkája megegyezik a külső erők virtuális munkájával. A (2) anyagegyenlet behelyettesítése után megkapjuk a rugalmas peremérték feladatunk gyengealakját: d (δw) dw ∫0 dz ⋅ AE ⋅ dz dz = ∫0 δw ⋅ p z dz + Fz ⋅ δwl l l (17) A (17) egyenlet igen fontos szerepet tölt be a közelítő megoldások előállítása során. Többek között ezt az alakot alkalmazzuk a végeselmes formulák előállítására is. 15 1.3 TELJES POTENCIÁLIS ENERGIA MINIMUM ELVE A teljes potenciális energia az alakváltozási energia és a külső erők potenciáljának összege. A külső erők potenciálja helyett szokás annak –1-szeresét a külső erők virtuális munkájának is nevezni. Így a telje
potenciális energia: Π p = U −W , (18) ahol U az alakváltozási energia, W a külső erők virtuálismunkája. Ez a kifejezés tulajdonképpen egy funkcionál A matematikában azokat az operátorokat, amelyeknek az értékkészlete valós számhalmaz funkcionáloknak nevezzük A konkrét (1)-(5) peremérték feladathoz rendelt funkcionál alakja: Π p (w ) = 2 1 dw AE ⋅ dz − ∫ w ⋅ p z dz − Fz wl . ∫ 20 dz 0 l l (19) Ha (19)-t úgy tekintjük, mint a w(z) tényleges megoldásra felír teljes potenciális energia akkor felvetődik a kérdés, hogy ehhez képest milyen potenciális energiát szolgáltat egy w* ( z ) kinematikailag lehetséges közelítő elmozdulásmező 2 l l dw* 1 dz − ∫ w* ⋅ p z dz − Fz wl . Π w = ∫ AE ⋅ 20 dz 0 * p ( ) * (20) Amint az a 3. ábrán is látható a közelítő mező felírható a tényleges és az elmozdulásmező variációja összegeként: w* ( z ) = w( z ) +
δw( z ) (21) Így a (20) potenciális energiába behelyettesítve (21)-t Π *p (w + δw) = 1 d (w + δ w ) AE ⋅ dz − ∫ (w + δw) ⋅ p z dz − Fz (wl + δwl ) , (22) ∫ 20 dz 0 l 2 l majd célszerűen átrendezve, az alábbi kifejezést kapjuk: 16 2 1 dw Π (w + δw) = ∫ AE ⋅ dz − ∫ w ⋅ p z dz − Fz ⋅ wl + 20 dz 0 l l * p d (δw) dw +∫ ⋅ AE ⋅ dz − ∫ δw ⋅ p z dz − Fz ⋅ δwl + dz dz 0 0 l l (23) 1 d (δw) AE ⋅ dz ∫ 20 dz 2 l + Ha megfigyeljük (23)-t láthatjuk, hogy az első sor megegyezik a (19) tényleges megoldás teljes potenciális energia kifejezés jobboldalával, a második sor a (17) virtuális munkaelv variációs alakjának nullára rendezett alakját szolgáltatja, a harmadik sor pedig egy kvadratikus kifejezés integrálja, ami biztos, hogy nagyobb mint nulla. Vagyis a teljes potenciális energia minimummal rendelkezik a tényleges
megoldásnál. A (23) második sora tulajdonképpen a teljes potenciális energia első variációját l l d (δw) dw δΠ p = ∫ ⋅ AE ⋅ dz − ∫ δw ⋅ p z dz − Fz ⋅ δwl , dz dz 0 0 (24) míg a harmadik sor a második variációját szolgáltatja δ 2Π p = 1 d (δw) AE ⋅ dz ∫ 20 dz l 2 , (25) ez valójában a w* közelítő megoldás hibájából származó alakváltozási energia. Első variáció :Egy funkcionál első variációját képezhetjük az ún. iránymenti de- rivált segítségével is δΠ p = d = dε d dε ε =0 ε =0 Π p (w + εδw) = 2 l 1 l d (w + εδw) (w + εδw) ⋅ p z dz − Fz (wl + εδwl ) dz − AE ⋅ ∫ ∫ dz 2 0 0 . (26) Első lépésben végezzük el a deriválást ε paraméter szerint l l d (δw) d (w + εδw) δΠ p = ∫ dz − ∫ δw ⋅ pdz − Fz ⋅ δwl ⋅ AE ⋅ dz 0 0 dz ε =0 (27) majd az ε =
0 helyettesítést és megkapjuk (24)-t: d (δw) dw δΠ p = ∫ ⋅ AE ⋅ dz − ∫ δw ⋅ p z dz − Fz ⋅ δwl . dz dz 0 0 l l 17 (28) A variáció számítás elméletéből ismert, hogy egy funkcionál minimuma létezésének szükséges feltétele, hogy az első variációja legyen zérus, azaz d (δw) dw ⋅ AE ⋅ dz − ∫ δw ⋅ p z dz − Fz ⋅ δwl , dz dz 0 0 l l δΠ p = 0 = ∫ (29) ami tulajdonképpen megegyezik a (17) virtuális munka elve variációs alakjával. Ez azt jelenti, hogy a két hibaelv rugalmas peremérték feladatokra nézve egyenértékű. Példa a minimum elvre: Egy k merevségű rugót F erő terheli. Határozzuk meg a rúgó z irányú w megnyúlását a teljes potenciális energia minimum elvének felhasználásával! F z k A rugalmas rendszer potenciális energiája: Π p (w ) = 1 k ⋅ w2 − F ⋅ w . 2 (30) Most a teljes potenciális energia kifejezése egy közönséges valós függvény. A potenciális energia minimum
elv értelmében a függvénynek keressük a minimumát A létezésének szükséges feltétele, hogy első deriváltja zérus legyen: [ ] d Π p (w ) = 0 = k ⋅ w − F . dw (31) A (31) egyenletből pedig átalakítás után megkapjuk a jól ismert összefüggést az elmozdulásra: w= F k (32) 18 1.4 RITZ-MÓDSZER 1D-S FELADATRA Legyen valamely rugalmas peremérték feladat kinematikailag lehetséges elmozdulása az alábbi alakban adott w* ( z ) = C0 + C1 z + C2 z 2 + L , (33) ahol C0 , C1 , C2 , L ismeretlen paraméterek. Ha ezt a közelítő megoldást a (20) funkcionálba helyettesítjük, akkor a teljes potenciális energia tulajdonképpen a C0 , C1 , C2 , L paraméterek többváltozós függvényévé válik Π *p (C0 , C1 , C 2 , L) . A teljes potenciális energia minimum elv alkalmazása ebben az esetben azt jelenti, hogy egy többváltozós függvény szélsőértékét keressük a szükséges feltétel alkalmazásával: min Π (C0 , C1 , C 2 ,L) ⇒ * p
∂Π *p ∂Ci =0 (i = 0,1, 2,L) . (34) (34) alkalmazása az ismeretlen paraméterek számával megegyező egyenletrendszert eredményez azok meghatározásához. Példa Ritz-módszerre (lineáris approximáció): Tekintsük ismét az 1.b ábrán adott feladatot Keressük a feladat közelítő megoldását az alábbi alakban w* ( z ) = C0 + C1 z (35) Ellenőrizzük, hogy a (35) feltételezett megoldás kinematikailag lehetséges-e. Kinematikai peremfeltétel w* (0) = 0 = C0 + C1 0 , (36) csak akkor teljesül, ha C0 = 0 , azaz w* ( z ) = C1 z és a deriválhatóság feltétele 19 (37) ε* = dw* = C1 , dz (38) is teljesül tehát a (37) alakja az elmozdulásnak kinematikailag lehetséges. A (37), (38)-t behelyettesítve (20)-ba kapjuk ( ) Π *p w = Π p (C1 ) = l l 1 2 AE ⋅ (C1 ) dz − ∫ C1 z ⋅ p z dz − Fz C1l ∫ 20 0 (39) A (34) minimum feltétel szerint: min Π (C1 ) * p ∂Π *p ⇒ ∂C1 l l 0 0 = 0 = ∫ AE ⋅ C1dz − ∫ z ⋅ p z dz −
Fz l. (40) A kijelölt integrálást elvégezve és az ismeretlen paramétert kifejezve l2 0 = AEC1l − p z − Fz l 2 ⇒ C1 = Fl + p z AE l 2. (41) felírhatjuk a közelítő megoldást: Fl + p z l 2 z, elmozdulás w* ( z ) = rúderő l N * ( z ) = Fl + p z . 2 AE (42) (43) Ha összevetjük a (42), (43) közelítő megoldásokat a megfelelő (12), (13) analitikus megoldásokkal akkor valóban mind a két esetben lényeges eltérés mutatkozik. A megoldás hibája: Közelítő megoldás minősítéséhez definiálnunk kell a hiba fogalmát: def e( z ) = wanalitikus ( z ) − w * ( z ) (44) ahol e(z ) a hiba és z-nek a függvénye. Gyakorlati számításokra a hiba energia normáját szokás alkalmazni, amely tulajdonképpen a hibafüggvényből számolt alakváltozási energia - a (25) alatti második variáció - négyzetgyöke: l e = 2 ( ) 1 de AE dz = Π *p w − Π p (w) . ∫ 20 dz A hiba mértéke lineáris approximáció
esetén: 20 (45) def e( z ) = wanalitikus ( z ) − w* ( z ) = l Fl + p z p pl p z 2 Fz + p z l 2 = − z = − z z2 + z z z + z − 2 AE 2 AE AE 2 AE AE (46) A hiba függvény deriváltja: pl p p l de = − z z + z = z − z + dz AE 2 AE AE 2 (47) Az energia normával értelmezett hiba: l 1 1 ( pz ) de e = AE dz = − z + dz = ∫ ∫ 20 2 0 AE 2 dz 2 l z l 1 ( pz ) z − l+z = 2 AE 3 2 4 0 2 = 3 2 2 2 2 l l ( pz ) 2 (48) 3 l AE 24 Joggal vetődik fel a kérdés, hogy hogyan változik a megoldás pontossága az approximáció fokának növelésekor. A következő példában az elmozdulást kvadratikus függvénnyel approximáljuk. Példa Ritz-módszerre (kvadratikus approximáció): Tekintsük ismét 1.b ábrán adott feladatot Keressük a feladat közelítő megoldását az alábbi alakban w* ( z )
= C0 + C1 + C2 z 2 (49) Ellenőrizzük, hogy a (49) feltételezett megoldás kinematikailag lehetséges-e. Kinematikai peremfeltétel w* (0) = 0 = C0 + C1 0 + C 2 0 , (50) csak akkor teljesül, ha C0 = 0 , azaz w* ( z ) = C1 z + C 2 z 2 és a deriválhatóság feltétele 21 (51) ε* = dw* = C1 + 2C 2 z , dz (52) is teljesül tehát a (49) alakja az elmozdulásnak kinematikailag lehetséges. A (51), (52)-t behelyettesítve (20)-ba kapjuk, hogy Π *p (w ) = Π p (C1 , C 2 ) = 1 2 AE ⋅ (C1 + 2C 2 z ) dz − ∫ (C1 z + C 2 z 2 )⋅ p z dz − Fz (C1l + C 2 l 2 ) ∫ 20 0 (53) l l A (34) minimum feltétel szerint: min Π *p (C1 , C2 ) ∂Π *p ∂C1 l l = 0 = ∫ AE ⋅ (C1 + 2C 2 z )dz − ∫ z ⋅ p z dz − Fz l. 0 ∂Π *p ∂C2 ⇒ (54) 0 l ( ) l = 0 = ∫ AE ⋅ 2C1 z + 4C 2 z dz − ∫ z 2 ⋅ p z dz − Fz l 2 . 2 0 (55) 0 A kijelölt integrálást elvégezve és az ismeretlen paramétereket kifejezve ( ) 0 = AE C1l + C2l 2 − p z l2 − Fz
l 2 4 l3 0 = AE C1l 2 + C 2 l 3 − p z − Fz l 2 3 3 . . (56) (57) C1 = Fz + p z l , AE C2 = − pz , 2 AE (58) felírhatjuk a közelítő megoldást: elmozdulás rúderő w* ( z ) = Fz + p z l p z − z z2, AE 2AE N * ( z ) = Fz + p z l − p z z. (59) (60) Ha összevetjük a (54), (66) közelítő megoldásokat a megfelelő (12), (13) analitikus megoldásokkal, akkor láthatjuk, hogy most vissza kaptuk az analitikus megoldást. 22 Amint azt a két különböző approximáció alkalmazásánál láttuk a közelítő megoldás pontosságának növelése egyik lehetséges útja az approximáció fokának növelése. A továbbiakban egy másik utat fogunk tanulmányozni, amikor is a tartományt résztartományokra bontjuk, és ezeken a résztartományokon az ismeretlen elmozdulásmezőt külön-külön lokálisan közelítjük. Ezt abban a reményben teszzük, hogy a számítás pontossága a
résztartományok számának növelésével szintén növelhető Ezeket a résztartományokat véges méretű elemeknek, tömören végeselemeknek fogjuk nevezni A résztartományok (elemek) határain csomópontokat jelölünk ki és az approximációt a közelitendő mező csomóponti értékein keresztül fejezzük ki. A módszer előnye, hogy könnyen programozható és így bonyolult szerkezetek nagypontosságú elemzésére nyílik lehetőség. Megoldások szemléltetése: w Analitikus megoldás (Kvadratikus közelítés) Lineáris közelítés z A megoldások szemléltetése: 23 N Analitikus megoldás (Kvadratikus közelítés) Lineáris közelítés z 4. ábra A számolt rúderők 24 2 LOKÁLIS APPROXIMÁCIÓ ELVE, VÉGESELEMES DISZKRETIZÁCIÓ EGYDIMENZIÓS FELADATRA A vizsgálatainkat továbbra is az (1)-(5) alatt definiált peremérték feladatra végezzük. Az L hosszúságú elemek határait csomópontok jelölik. pz L Fz a. L w w3 w2 b. 1 1 2
2 3 z b. w w2 c. 1 1 2 z w w2 w3 d. 2 2 5. ábra 25 3 z A korábban 1.b ábrán vázolt tartományt most gondolatban két egyenlő hosszúságú résztartományra, azaz végeselemre bontjuk és az elmozdulásmezőt az elemeken külön-külön approximáljuk Elemenként külön-külön lineárisan közelítjük az ismeretlen elmozdulásmezőt (5.c és d ábrák) és gondoskodunk ezek illesztéséről is. (A b ábrán folytonos vonal jelöli az egzakt megoldást, szaggatott a közelítést) Ez a közelítés felépíthető csomópontokhoz rendelt approximációs függvények segítségével is. Approximáció csomópontokhoz rendelt függvényekkel: w* w3 w2 a. 1 1 2 2 3 z N1 w1 b. 1 2 3 z N2 w2 1 1 c. 1 2 3 z N3 w3 1 1 d. 1 2 3 6. ábra 26 z Az a. ábrán vázolt közelítő függvény (folytonos vonal) felépíthető szakaszonként lineáris függvények (szaggatott vonal) összegeként, azaz a csomópontokhoz rendelt N i
(i=1,2,3) alakfüggvények lineáris kombinációjaként (lásd 6. b,c,d ábrákat) A csomóponti függvények lineáris kombinációja szolgáltatja az approximációt: 3 w* ( z ) = ∑ N i (z ) ⋅ wi (megjegyezzük, hogy w1 = 0 ). i =1 27 (61) 2.1 HÚZOTT-NYOMOTT RÚDELEM Az 1. fejezetben láthattuk, hogy közelítő megoldás keresése során (Ritz-módszernél) a szerkezet teljes potenciális energiáját kellet felírni. Hasonlóan kell eljárni, amikor a vizsgált tartományt résztartományokra, azaz végeselemekre bontjuk A szerkezet teljes potenciális energiája az egyes elemeken számolt potenciális energiák összegeként állítható elő a koncentrált erő munkájával együtt: 2 Π p = ∑ Π ep − Fz w3 . (62) e =1 ahol (62) az e index a végeselem sorszámát jelöli, itt és a továbbiakban közelítő megoldásról fogunk beszélni, de a korábban alkalmazott * jelölést a jobb felső indexben elhagyjuk. Az elemen számolt potenciális energia L L
dwe 1 dwe dζ − ∫ we ⋅ p z dζ ⋅ AE ⋅ Π (wi , w j ) = ∫ 2 0 dζ dζ 0 e p (63) Tekintsük a 5.d ábrán és a 7 ábrán is vázolt 2-es végeslemen az elmozdulás approximációt w w3 w2 z 2 L 3 ζ 7. ábra: Approximáció 2-s elemen A 7. ábrán az elemhez kötötten bevezettünk egy új ζ koordinátát oly módon, hogy az origója egybe essen az elem egyik végpontjával és iránya megegyezik az eredeti z iránynyal. Írjuk fel a ζ koordináta segítségével kifejezve az elmozdulásmezőt az elem mentén: 28 w3 − w2 ζ, L w 2 (ζ ) = w2 + (64) ahol a w jobb felső indexében a 2-es szám az elem sorszámára utal. Rendezzük át (64)-t a csomóponti elmozdulások szerint, majd sor és oszlop mátrixokkal is kifejezve ζ ζ w2 ζ ζ ⋅ = [w2 w 2 (ζ ) = 1 − w2 + w3 = 1 − L L L L w3 ζ
1 − L . w3 ] ⋅ ζ L (65) A 2-es elem approximációja alapján felírhatjuk egy általános e sorszámú és i, j csomópontú elem közelítését is ζ ζ wi ζ ζ w e (ζ ) = 1 − wi + w j = 1 − ⋅ = wi L L L L w j [ ζ 1 − w j ⋅ L . (66) ζ L ] A (66)-s képlet alkalmas az e=1-es sorszámú elem elmozdulásának a leírására is a megfelelő i=1és j=2 csomóponti elmozdulás behelyettesítésével. A (66) elmozdulásmező ismeretében az (1) egyenlet segítségével számolhatjuk az alakváltozást, a deriválást értelemszerűen ζ szerint hajtjuk végre ε ze (ζ ) = dw (ζ ) = dζ e w j − wi L −1 = L 1 wi ⋅ = wi L w j [ − 1 wj ⋅ L 1 L ] (67) Alkalmazva a (2)
anyagtörvényt meghatározható az elemen a rúderő is −1 N e (ζ ) = AEε ze (ζ ) = AE L 1 wi ⋅ = wi L w j [ − 1 w j ⋅ L AE 1 L ] (68) Most már felírhatjuk az e-ik elem potenciális energiáját is (66) és (67) segítségével, (63) felhasználásával 29 − 1 −1 w j L AE 1 L L [ L 1 wi dζ − ∫ wi L w j 0 ] 1 Π ep (wi , w j ) = ∫ wi 20 [ L ζ 1 − L p dζ . w j z ζ L ] (69) A (69) első integrálja az elem alakváltozási energiája U e , a második integrál a megoszló erőrendszer munkája W e . A csomóponti paraméterek az integrálás szempontjából konstansnak tekinthetők ezért kiemelhetjük az integrál jel elé [ 1 U e = wi 2 − 1 −1 w j ∫ L AE 1 0 L
L 1 wi 1 dζ = wi L w j 2 ] [ L AE 2 wj ∫ L AE 0 − L2 ] L AE L2 dζ wi AE w j L2 − (70) A (70)-ben lévő mátrix egy elemének integrálja L AE AE ∫0 L2 dζ = L2 ζ L = 0 AE . L (71) Visszahelyettesítve (71)-t (70)-be Ue = AE wj L AE − L [ ] 1 wi 2 − AE L wi = 1 q eT K e q e , AE w j 2 L (72) ahol a 2x2-es mátrixot az e elem merevségi mátrixnak nevezzük és K e -vel jelöljük, a függőleges 2x1-es oszlopvektort és a vízszintes 1x2-es sorvektort az e elem csomóponti vektorának nevezzük és q e , q eT -vel jelöljük, T a transzponálás jele. A (69)-ben szereplő második integrál a külső erő munkája ez elemen W e [ W e = wi ζ 1 − L p dζ w j ∫ z ζ 0 L ] L A (73) oszlopvektor elemeinek integrálása
után 30 (73) ζ2 ζ ∫0 1 − L p z dζ = ζ − 2 L p z L L ζ ∫ L p dζ z = 0 L ζ2 = pz 2L 0 L pz L , 2 = 0 (74) pz L , 2 (75) kapjuk, hogy pz L 2 eT e wj = q fp , pz L 2 [ ] W e = wi (76) ahol a 2x1-es oszlopvektor az elem tehervektora és f pe -vel jelöljük. Végül is egy általános e elem potenciális energiája Π ep = [ 1 wi 2 AE wj L AE − L ] AE L wi − w i AE w j L − [ pz L wj 2 . p L z 2 ] (77) 2.11 SZERKEZETI MÁTRIXOK Az elemek potenciális energiájának ismeretében felírhatjuk a szerkezet potenciális energiáját 2 Π p (w1 , w2 , w3 ) = ∑ Π ep − Fz w3 = e =1 1 = [w1 2 1 + [w2 2 AE w2 ] L AE − L AE w3 ] L AE − L AE L w1 − [w 1 AE w2
L − AE L w2 − [w 2 AE w3 L − 31 pz L w2 ] 2 + pL z 2 pz L w3 ] 2 − Fz w3 . p L z 2 (78) A szerkezet teljes potenciális energiája tömörebben is felírható, hiszen a szomszédos elemek közös csomópontjában az elmozdulás megegyezik – jelen esetben w2 – így értelemszerűen csak egyszer szerepeltetjük a kifejezésben, ezáltal biztosítjuk az elemek illesztését Π p (w1 , w2 , w3 ) = 1 [w1 2 w2 AE L AE w3 ]− L 0 0 w1 AE w2 − [w1 − L AE w3 L AE L AE 2 L AE − L − w2 pz L 2 pz L w3 ] pz L + F z 2 (79) Az elemek illesztésén túl van egy másik fontos feltétel - amit még teljesítenünk kell - a kinematikai peremfeltétel. Ez azt
jelenti, hogy a befalazásnál lévő csomópontban gon- doskodnunk kell arról, hogy w1 = 0 legyen Π p (w2 , w3 ) = 1 [0 w2 2 AE L AE w3 ]− L 0 AE L AE 2 L AE − L − 0 0 AE w2 − [0 w2 − L AE w3 L pz L 2 p z L . (80) w3 ] pz L + F z 2 Figyeljük meg, hogy 0-val szorozzuk a merevségi mátrix első sorát és oszlopát, a tehervektor vonatkozásában meg csak az első elemet. Ezért ezen sor és oszlop, valamint elem a szerkezeti mátrixból és vektorból elhagyható, 1 Π p (w2 , w3 ) = [w2 2 ahol AE 2 w3 ] L AE − L AE L w2 − [w 2 AE w3 L − pz L , w3 ] p z L + Fz 2 (81) a 2x2-as mátrixot szerkezeti merevségi mátrixnak nevezzük és K -val jelöljük, a csomóponti elmozdulásokat
tartalmazó függőleges 2x1-es oszlop- és a vízszintes 1x2-es sorvektort szerkezeti csomóponti vektornak nevezzük és q, q T -vel jelöljük, végül a 2x1-es 32 oszlopvektort, amely a szerkezet terhelését tartalmazza szerkezeti tehervektornak nevezzük és f -vel jelöljük. E jelölések bevezetése után (81) a szerkezeti mátrixokkal is felírható Π p (q ) = 1 T q Kq − q T f . 2 (82) 2.12 A CSOMÓPONTI ELMOZDULÁSOK MEGHATÁROZÁSA A (81) szerkezet teljes potenciális energiája a csomóponti elmozdulási paraméterek kétváltozós függvénye. A potenciális energia minimum elv értelmében keressük ennek a többváltozós függvénynek a minimumát, a minimum létezésének szükséges feltétele a Ritzmódszernél is bemutattak szerint (lásd (34)-t) min Π p (w2 , w3 ) ⇒ 0= ∂Π p (w2 , w3 ) ∂w2 , 0= ∂Π p (w2 , w3 ) ∂w3 . (83) ugyanez tömörebben is felírható min Π p (q ) ⇒ 0= ∂Π p (q ) ∂q (84) Először (84)-s alatti
műveleteket mutatjuk be 0= ∂Π p (q ) ∂q 1 ∂ q T Kq − q T f 2 = Kq − f , = ∂q (85) átrendezés után egy lineáris egyenletszert kapunk Kq = f azaz 33 (86) AE w p z L 2 L = . AE p z L w3 + Fz L 2 AE 2 L AE − L − (87) Ez a lineáris algebrai egyenletrendszer két egyenletet és két ismeretlent tartalmaz, mivel az együtthatók determinánsa nyilvánvalóan nem nulla biztos, hogy megoldható a csomóponti elmozdulás paraméterekre: w2 = 3 p z L2 Fz L + , 2 AE AE w3 = 2 p z L2 FL +2 z . AE AE (88) A rúderőket (68) felhasználásával állíthatjuk elő: −1 N (ζ ) = AE L 1 −1 N 2 (ζ ) = AE L = 3 p z L + Fz 2 1 w1 −1 ⋅ = AE L w2 L 0 1 3 p L2 F L ⋅ z z L 2 AE + AE
1 w2 −1 ⋅ = AE L w3 L 3 p z L2 Fz L + AE 1 1 2 AE ⋅ = p z L + Fz L p z L2 Fz L 2 +2 2 AE AE N N1 Analitikus megoldás N2 1 2 8. ábra: A rúderő 34 z A (83) alatti műveletek részletezése: Először végezzük el (81)-ben kijelölt szorzásokat Π p (w2 , w3 ) = [ ] 1 AE pL 2 2 2(w2 ) − 2 w2 w3 + (w3 ) − w2 p z L − w3 z + Fz , (89) 2 L 2 és helyettesítsük be (83)-két egyenletébe [ ] 1 AE p L 2 2 2(w2 ) − 2 w2 w3 + (w3 ) − w2 p z L − w3 z + Fz ∂ ∂Π p (w2 , w3 ) 2 L 2 0= = = ∂w2 ∂w2 = 1 AE [4(w2 ) − 2w3 ] − p z L = = AE [2(w2 ) − w3 ] − p z L , L 2 L [ (90) ] 1 AE pL 2 2 2(w2 ) − 2 w2 w3 + (w3 ) − w2 p z L − w3 z + Fz ∂ ∂Π p (w2 , w3 ) 2 L 2 0= = =
∂w3 ∂w3 = 1 AE [− 2w2 + 2w3 ] − p z L + Fz = AE [− w2 + w3 ] − p z L + Fz . 2 L 2 2 L (91) Ha megnézzük a (90)-ben kapott eredményt akkor látjuk, hogy (87) első sorát kaptuk vissza, hasonlóan (91) alatti eredmény (87) második sorával egyezik meg. A végeselem-módszer rövid összefoglalása: A végeselem-módszer alapgondolata, hogy a vizsgált szerkezetet gondolatban véges számú részre, elemre bontjuk, s a keresett megoldást elemenként külön - külön közelítjük. Az elemek valóságos kapcsolódásának megfelelően az elemeket illesztjük egymáshoz Erre szolgálnak az elemek határain kijelölt kapcsolódási pontok, vagy csomópontok, illetve azok elmozdulásai. Így a teljes szerkezetre érvényes közelítést kapunk, amely már csak a csomópontok jellemző elmozdulásait tartalmazza. Ebből kiindulva kiszámítható a teljes szerkezet alakváltozási energiája a csomóponti elmozdulások
függvényében. Hasonlóan 35 eljárva, a szerkezetre működő terhelések munkáját is kifejezhetjük a csomópontok elmozdulásaival. A szerkezet mozgását korlátozó kényszereket is csomópontokra vonatkozó kinematikai előírásokkal vesszük figyelembe Ez azt jelenti, hogy a megfelelő csomópont teljes, vagy adott irányú elmozdulását meggátoljuk. Energetikai megfontolásokból származtatható a közelítésben felvett összes csomóponti paraméter (elmozdulási koordináták) kiszámítására szolgáló egyenletrendszer. Lineárisan rugalmas szerkezet statikus terhelése mellett ez nagyméretű lineáris egyenletrendszer, amely a szerkezet egyensúlyát fejezi ki. Az egyenletrendszer megoldása után - a csomóponti paraméterekből kiindulva -, meghatározható tetszőleges szerkezeti elem szilárdságtani állapota, azaz tetszőleges pontban megkapjuk az elmozdulási, alakváltozási és feszültségi állapot jellemzőit Az eddigi végeselemes
számításnál kihasználtuk, hogy a z és ζ koordináta tenge- lyek egyirányúak. A továbbiakban olyan szerkezetet fogunk vizsgálni, ahol a rúdelemek szöget zárnak be egymással. 36 2.2 RÁCSOS SZERKEZET VIZSGÁLATA HÚZOTT-NYOMOTT RÚDELEMEKKEL Y 1 ζ2 Fy e=1 ζ1 3 1m e=2 2 α2 2m Z 9. ábra: A szerkezet geometriája: L1 = 1 m, L2 = 1 + 2 2 = 5 m, A1 = 10 −4 m 2 , A2 = 2 ⋅10 −4 m 2 , α1 = 0 , sin α 2 = 1 2 , cosα 2 = , 5 5 Anyagjellemző: mindkét rúd acélból van E = 2 ⋅ 1011 N Terhelés: Fy = 100 N . m2 Jelölje rendre V,W az Y, Z irányú elmozdulást. A kinematikai peremfeltétel: V1 = 0, W1 = 0, V2 = 0, W2 = 0 . A megoldás során úgy mint korábban most is a szerkezet teljes potenciális energiájának előállítása útján történik 2 Π p = ∑ Π ep − FyV3 . e =1 37 (92) Az elemeken most megoszló erő nem működik, ezért az elemek potenciális energiája egyenlő az alakváltozási energiával. Felhasználva
a lokális koordinátarendszerben felírt (77) kifejezést a két rúdra az energia Π1p = Π 2p = [ 1 2 w2 2 [ 1 1 w1 2 A1 E L 1 w3 1 − A1 E L1 ] A2 E L 2 w3 2 − A2 E L2 A1 E 1 L1 w1 , A1 E 1 w3 L1 − A2 E 2 L2 w2 2 , A2 E w3 L2 (93) − ] (94) ahol a két különböző elemen a megfelelő lokális koordináta rendszerben (bal felső index jelöli) értelmezett 1 w3 ≠ 2 w3 elmozdulás nem azonos. Hiszen két különböző irányú ζ 1 ,ζ 2 tengelyhez tartoznak. Ezt a problémát csak úgy tudjuk feloldani, ha közös koordinátarendszert alkalmazunk, azaz a globális Y,Z rendszert Y V w W α Z 10. ábra: az elmozdulás transzformációja A 10. ábrán látható geometriai viszonyok alapján a lokális és globális elmozdulások közötti kapcsolat w = V sin α + W cosα = [sin α sin α V cosα ]
= [V W ] . cosα W 38 (95) Egy általános e elem elmozdulásának a transzformációja i,j csomóponti elmozdulásokkal kifejezve wi sin α e e = w j 0 cosα e 0 e 0 sin α e Vi 0 Wi cosα e V j W j (96) Az e=1 elem elmozdulásának transzformációja ( α 1 = 0 ): V1 w1 0 1 0 0 W1 1 = w3 0 0 0 1 V3 W 3 1 1 w2 5 2 = w3 0 2 (97) 2 5 0 0 1 5 V2 0 W 2 2 V3 5 W3 (98) Fejezzük ki (93) és (94) energiákat a globális elmozdulásokkal Π1p = Π 2p = 1 [V1 W1 V3 2 1 [V2 W2 V3 2 W3 ] 0 1 W3 ] 0 0 1 5 2 5 0 0 0 A E 1 0 L1 0 − A1 E 1 L1 0
AE 0 2 L2 1 A2 E − 5 L2 2 5 V1 A1 E L1 0 1 0 0 W1 , (99) A1 E 0 0 0 1 V3 W L1 3 − A2 E 1 L2 5 A2 E 0 L2 − 2 5 0 0 1 5 , és a kijelölt szorzást elvégezve 39 V2 0 W 2 V 2 3 5 W3 (100) Π1p = 1 [V1 W1 V3 2 1 [V2 W2 V3 2 Π 2p = 0 0 A1 E 0 L1 W3 ] 0 0 A 0 − 1 E L1 A2 E 5L 2 A 2 2E 5 L2 W3 ] AE − 2 5 L2 2 A2 E − 5 L2 2 A2 E 5 L2 4 A2 E 5 L2 2A E − 2 5 L2 4 A2 E − 5 L2 0 A1 E V1 0 − L1 W1 , 0 0 V3 A1 E 0 W3 L1 0 A2 E 5 L2 2A E − 2 5 L2 A2 E 5 L2 2 A2 E 5 L2 − (101) 2 A2 E 5 L2 4 A2 E V2 − 5 L2 W2 (102) 2 A2 E
V3 5 L2 W3 4 A2 E 5 L2 − a globális koordinátákra Πp = 1 [V1 W1 V2 2 0 0 A1 E 0 L1 0 0 0 0 0 0 A1 E 0 − L 1 W2 V3 W3 ] 0 0 0 0 0 0 A2 E 5 L2 2 A2 E 5 L2 AE − 2 5 L2 2A E − 2 5 L2 2 A2 E 5 L2 4 A2 E 5 L2 2A E − 2 5 L2 4A E − 2 5 L2 A2 E 5 L2 2 A2 E − 5 L2 A2 E 5 L2 2 A2 E 5 L2 − 0 V1 A1 E − L1 W1 2 A2 E − V 5 L2 2 A2 E − FyV3 − 5 L2 W2 2 A2 E 5 L2 V3 A1 E 4 A2 E + 5 L2 W3 L1 (103) Ha figyelembe veszzük a kinematikai peremfeltételt, hogy V1 = 0, W1 = 0, V2 = 0, W2 = 0 vagyis az első négy sor és oszlop nullával szorzódik, akkor az energia kifejezés lényegesen egyszerűsödik 40 Πp = 1 [V3 2 A2 E 5L 2 W3 ] 2 A2 E 5L 2
2 A2 E V 5L2 F 3 − [V3 W3 ] y , A1 E 4 A2 E W3 0 + 5 L2 L1 (104) A (104) kifejezés tömörebben is írható Π p (q ) = 1 T q Kq − q T f . 2 (105) Alkalmazva a potenciális energia minimum elvét min Π p (q ) ⇒ 0= ∂Π p (q ) ∂q (106) a deriválást végrehajtva 1 ∂ q T Kq − q T f ∂Π p (q ) 2 = Kq − f , 0= = ∂q ∂q (107) átrendezés után a jól ismert lineáris egyenletszert kapjuk Kq = f . (108) részletesen felírva 2 A2 E V3 Fy 5L2 = A1 E 4 A2 E W + 0 5L2 3 L1 A2 E 5L 2 2 A2 E 5L 2 (109) a megoldás a globális koordinátákra V3 = −2 Fy L1 5 A1 L2 100 ⋅1 5 ⋅10 −4 5 5 5 − m, − − 2 = −2 −4 − 2 = 10 −5 2 +
11 −4 A1 E 2 A2 L1 4 ⋅ ⋅ ⋅ ⋅ 10 2 10 2 2 10 W3 = −2 Fy L1 A1 E = −2 100 ⋅1 = −10 −5 m. 11 10 ⋅ 2 ⋅10 −4 A megoldás a lokális rendszerben (97) és (98) alapján 41 0 0 0 1 w1 0 1 0 0 5 5 − 5 1 = + = − 10 −5 w3 0 0 0 1 10 2 4 − 10 −5 1 2 w2 5 2 = w3 0 2 5 0 0 1 5 0 0 0 0 10 −5 2 + 5 5 = 5 ⋅ 10 −5 2 4 4 5 − 10 −5 (110) (111) A rúderők kiszámítása: (68) felhasználásával történik −1 N = A1 Eε ζ (ζ ) = A1 E L1 1 1 1 1 w1 ⋅ = L1 1 w3 − 1 1 0 = 10 −4 ⋅ 2 ⋅1011 ⋅ =
−200 N . −5 1 1 − 10 −1 N = A2 Eε ζ (ζ ) = A2 E L2 2 2 −1 = 2 ⋅10 − 4 ⋅ 2 ⋅1011 5 2 1 w2 = ⋅ L2 2 w3 1 0 ⋅ 5 −5 = 100 5 N 5 4 10 (112) (113) A végeselem-módszerrel itt előállított megoldás most azonos az analitikus megoldással. 42 3 IZOPARAMETRIKUS ELEMCSALÁD A kereskedelmi szoftverekben leggyakrabban ún. izoparametrikus elemeket alkalmaznak Az „izoparametrikus” jelző abból adódik, hogy az ilyen elem típus egyik jellemzője, hogy a geometria leképzésére alkalmazott paraméterek száma azonos az ismeretlen elmozdulásmező közelítésére felvett paraméterek számával. Az elterjedését az segítette elő, hogy az elem merevségi mátrixának és tehervektorának előállításakor a numerikus integrálás könnyen végrehajtható. A leggyakrabban alkalmazott izoparametrikus elemek: 1D-s elemek 2 és 3
csomópontú rúdelemek 2D-s elemek 11. ábra 3 és 6 csomópontú háromszög, 4 és 8 csomópontú négyszög síkbeli elemek, tengelyszimmetrikus elemek, héjelemek 3D-s elemek 4 és 6 csomópontú tetraéder elemek 12. ábra 8 és 20 csomópontú hexahedron elemek Az itt ábrázolt egyenes élű elemek - csak az élek végpontjain tüntetettünk fel csomópontot - lineáris közelítést tartalmaznak, míg a görbült élű elemek - amelyeknél az ol- 43 dalfelezőnél is feltüntettünk csomópontot - kvadratikus közelítést alkalmaznak mind a geometria és mind az elmozdulásmezőre nézve. A kvadratikus elemekkel tehát pontosabb megoldás nyerhető, hiszen a vizsgált tartomány geometriáját jobban megközelíthetjük és az ismeretlen elmozdulásmező is magasabb fokú függvénnyel approximáljuk. A továbbiakban részletesen megmutatjuk a rúdelem és a síkbeli négy csomópontú elem merevségi mátrixának és tehervektorának származtatását. 44 3.1
1D-S HÚZOTT-NYOMOTT RÚDELEM Vizsgáljuk ismét az 2. fejezetben is bemutatott húzott-nyomott rúdfeladatot (Vesd össze za 1. ábrát és 13 ábrát) pz Fz L L 13. ábra Tekintsük most csak a második elemet, mint az általános i,j csomópont párral adott e-dik elemet! y i e=2 P’(z) j z zi wi zj L wj z P’(z) zj b zi -1 0 P( ξ ) 1 ξ 14. ábra Geometriai leképzés z ⇔ ξ koordináták között A 14. ábrán a rúdelemhez egy lokális ún természetes ξ koordinátatengelyt kötöttünk Keressük a természetes ξ koordinátájú pont és a hozzátartozó pont globális z koordi- 45 nátája közötti kapcsolatot, azaz a leképző függvényt. A ξ tengelyre merőlegesen felmérjük a csomópontok z koordinátáját és egyenessel összekötve megkapjuk a leképzés függvény képét. Egy tetszőleges P( ξ ) pontból függőlegesen felvetítve megkapjuk a hozzárendelt P’(z) képét. A leképző függvény m meredeksége és b metszéspontja
alapján könnyen felírhatjuk az egyenes egyenletét, amelyet célszerűen átrendezünk z (ξ ) = z j − zi z j + zi 1 (1 − ξ )zi + 1 (1 + ξ )z j , (3.1) 2 2 2 2 1 1 ahol a csomóponti koordináták együtthatói ( N1 (ξ ) = (1 − ξ ) és N 2 (ξ ) = (1 + ξ ) ) az ún. 2 2 alakfüggvényeknek. N1 (ξ ) = ξ+ = 1 (1 − ξ ) 2 N 2 (ξ ) = 1 -1 1 (1 + ξ ) 2 1 0 1 ξ -1 0 1 ξ 15. ábra A 2 csomópontú izoparametrikus elem alakfüggvényei Az elmozdulásmező közelítését, ezen két alakfüggvény és a wi , w j csomóponti elmozdulások segítségével fogjuk közelíteni w 1 1 1 1 (1 + ξ ) i = N e (ξ )q e (3.2) w e (ξ ) = (1 − ξ ) wi + (1 + ξ ) w j = (1 − ξ ) 2 2 2 2 w j Az elmozdulás ismeretében az alakváltozás előállítható a láncszabály alkalmazásával ε ze (ξ ) = dw e (ξ ) dw e (ξ ) dξ = , dz dξ dz (3.3) ahol az első tag ξ -szerinti deriválása (3.2) behelyettesítése után
végrehajtható, a második tag közvetlenül nem, de az inverze (3.1) ismeretében képezhető dz (ξ ) z j − z i L = = . dξ 2 2 Ezután visszahelyettesítve (3.2)-t és (34) reciprokát (33)-ba kapjuk, hogy 46 (3.4) ε ze (ξ ) = 1 wi 1 = − 2 w j L 2 1 − L 2 1 − L wj 1 L 1 wi = wi L w j [ ] (3.5) Az alakváltozás ismeretében a feszültség az egyszerű Hooke-törvény alapján számítható 1 1 wi σ e (ξ ) = Eε e (ξ ) = E − (3.6) L L w j Célunk, hogy előállítsuk az elem potenciális energiáját (w , w ) = 1 ∫ dw L dwe ⋅ AE ⋅ dz − ∫ we ⋅ p z dz 2 0 dz dz 0 L Π e p i j e (3.7) ahol az első integrálból származtatható az elem merevségi mátrixa, a másodikból az elem tehervektora. Az
integrálást most nem z-szerint hanem ξ szerint hajtjuk végre, a dz előállításához felhasználjuk (34)-et dz = dz (ξ ) L dξ = dξ . dξ 2 (3.8) A merevségi mátrix: 1 − L 1 1 L wi wj ∫ AE − dξ 1 L L 2 w j ξ = −1 L4 44244444 4 14 3 dwe 1 dwe 1 dz = wi ⋅ AE ⋅ ∫ 2 0 dz 2 dz [ L ] 1 Ke (3.9) ahol AE K e = ∫ 2L AE −1 − 2L AE AE 2 L dξ = L AE AE − 2L L − 1 AE L AE L − (3.10) Konstans megoszló terhelést feltételezve a tehervektor: L ∫w e [ ⋅ p z dz = wi 0 1 2 (1 − ξ ) L wj ∫ p z dξ = wi 1 2 ξ = −1 (1 + ξ ) 24424 443 14 ] [ 1 fe ahol az egyes elemek integráljai 47 pz L wj 2 p L z 2 ] (3.11) 1 [ ] = pL 2
[ ] = pL 2 ∫ 2 (1 − ξ ) pL pL ξ −ξ 2 dξ = 2 4 ∫ 2 (1 + ξ ) pL pL ξ −ξ 2 dξ = 2 4 −1 1 1 1 −1 1 −1 1 −1 Végül az elem potenciális energiája: Π ep = [ 1 wi 2 AE wj L AE − L ] AE L wi − w i AE w j L − [ pz L wj 2 . pz L 2 ] (3.12) amely a 2. fejezetbeli (70) alakkal egyezik Így a további lépések az alkalmazás tekintetében azonosan hajtható végre, mint a 2 és 3 fejezetben ezért a megoldás lépéseit itt már nem ismételjük meg. 48 3.2 SÍKBELI FELADATOK A valóságban térbeli mechanikai problémák egy része bizonyos feltételek esetén visszavezethetők síkbeli, pontosabban 2 dimenziós (2D-s) feladatokra. Ezen 2D-s feladatok közül a következő három hasonlóan tárgyalható: - általánosított síkfeszültségi állapotú feladat, azaz tárcsafeladat, - síkalakváltozási feladat, -
tengelyszimmetrikus feladat. 3.21 ÁLTALÁNOSÍTOTT SÍKFESZÜLTSÉGI ÁLLAPOT( ÁSF, TÁRCSAFELADAT, PLANE STRESS PROBLEM): y p vi i ui x 16. ábra Általánosított síkfeszültségi feladat Lemezszerű szerkezet sajátsíkjában terhelt. A feszültségek valójában a b falvastagság mentén képezett átlagértékek, de ezt külön nem jelöljük A feszültségi tenzor és a független elemekből képzett feszültségi vektor: 49 σ x τ xy T = τ yx σ y 0 0 0 0 , 0 σ x σ = σ y . τ xy (3.13) Hasonló alakot ölt az alakváltozási tenzor és a független elemekből képzett alakváltozási vektor: εx 1 A = γ yx 2 0 1 γ xy 2 εy 0 0 0, εz εx ε = ε y , γ xy (3.14) Az alakváltozási vektorban ε z mennyiséget azért nem tüntettük fel, mert az alakváltozási energiába sem játszik
szerepet, hiszen a feszültségi párja zérus. A feladat jellemzője, hogy a végeselem háló csomópontjaiban csak x és y irányú elmozdulás ismeretlen paraméterekről beszélünk, valamint ennek megfelelő irányú erők működtethetők. 3.22 SÍKALAKVÁLTOZÁSI ÁLLAPOT (PLANE STRAIN PROBLEM): y vi ui x 17. ábra Egy folyómentén épített gát keresztmetszete A keresztmetszett síkjára merőlegesen végtelen hosszúnak tekintett test bármelyik keresztmetszetében ugyan olyan alakváltozási és feszültségi állapot ébred. A modell egység- 50 nyi vastagságú metszetet vizsgál. Az alakváltozási tenzor és a független elemekből képzett alakváltozási vektor: εx 1 A = γ yx 2 0 0 εx ε y 0 , ε = ε y . (3.14) γ xy 0 0 Hasonló alakot ölt a feszültségi tenzor és a független elemekből képzett feszültségi vektor: A feszültségi vektorban σ z 1 γ
xy 2 σ x σ x τ xy 0 σ = σ y . (3.15) T = τ yx σ y 0 , τ xy 0 0 σ z mennyiséget azért nem tüntettük fel, mert az alakváltozási energiába sem játszik szerepet, hiszen az alakváltozási párja zérus. A feladat jellemzője, hogy a végeselem háló csomópontjaiban csak x és y irányú elmozdulás ismeretlen paraméterekről beszélünk, valamint ennek megfelelő irányú erők működtethetők. 3.23 FORGÁS V. TENGELYSZIMMETRIKUS ÁLLAPOT (AXI-SYMMETRIC PROBLEM): z wi ui p r 18. ábra Egy csavar tengelyszimmetrikus terhelése 51 A forgásszimmetrikus test terhelése is forgásszimmetrikus és bármelyik meridián metszetében ugyan olyan alakváltozási és feszültségi állapot ébred. Az alakváltozási tenzor és a független elemekből képzett alakváltozási vektor: 1 εr γ rz 0 εr εϕ 2 A= 0 , (3.14) ε 0 ε = . ϕ ε
z 1 γ γ rz 0 ε zr z 2 Hasonló alakot ölt a feszültségi tenzor és a független elemekből képzett feszültségi vektor: σ r σ ϕ σ= . (3.15) σ ϕ 0 , σz 0 σ z τ xz A feladat jellemzője, hogy a végeselem háló csomópontjaiban csak r és z irányú σ r T = 0 τ zr τ rz 0 elmozdulás ismeretlen paraméterekről beszélünk, valamint ennek megfelelő irányú erők működtethetők. A három feladat végeselem vizsgálata azért nagyon hasonló, mert a csomóponti elmozdulásnak csak síkba eső koordinátája fordul elő. A továbbiakban részletesen csak a síkfeszültségi állapotú végeselem előállítását részletezzük. 3.24 RUGALMAS SÍKBELI (ÁSF) PEREMÉRTÉK FELADAT KITŰZÉSE y Ap p n V Au r r r ρk x 19. ábra 52 Elmozdulásmező (elsődleges ismeretlen): r r r r u * ( r ) = u ( x , y )e x + v ( x, y ) e y , u ( x, y
) u= , v( x, y ) (3.16) εx ε = ε y γ xy (3.17) σ x σ = σ y , τ xy (3.18) ahol a * a feladat síkbeli jellegére utal. Az átlagos alakváltozások síkbeli része: ( ) r 1 r A = u * o ∇ + ∇ o u , 2 * ahol ∇ * = ν (ε x + ε y ) ∂ r ∂ r ex + e y , és a síkra merőleges átlagos nyúlás ε z = − ∂y ∂x 1 −ν Az átlagos feszültségek tenzora (Hooke-törvény): ν (ε x + ε y ) * T = 2G ( A + I ), 1 −ν * * az anyagtörvény mátrixos formában is felírható: 1 ν σ x 0 ε x E ν 1 0 ε y σ = Dε = σ y = 2 1 −ν 1 −ν τ xy 0 0 γ xy 2 144424443 (3.19) D Egyensúlyi egyenlet: k x k= . k y ahol k pl. a gravitációs vagy a forgásból származó gyorsulás vektora r r T* o ∇ + ρk = 0 , (3.20)
Kinematikai peremfeltétel: ( ) ( ) r r r r u * r = u0 r r r * ∈ Au 53 (3.21) Dinamikai peremfeltétel: px p= . (3.22) py A síkbeli feladat összesen 8 ismeretlen mezőt tartalmaz, ezek egyértelmű meghatározásához 8 skalár egyenlet (részben parciális differenciálegyenlet) és a megfelelő peremfeltételek állnak rendelkezésre. Az olyan megoldást amely (316)-(322)-nek eleget tesz analitikus megoldásnak, egzakt megoldásnak vagy tényleges megoldásnak nevezzük Azonban a következőkben is közelítő megoldást keresünk a potenciális energia mi- r r T* ⋅ n = p r r * ∈ Ap , nimum elv felhasználásával. 3.25 A TELJES POTENCIÁLIS ENERGIA A 4.24 pont alatt bevezetett mennyiségekkel a rugalmas síkbeli feladatra a teljes potenciális energia az alábbi alakban írható: Π p (u ) = 1 T ε DεdV − ∫ u T ρkdV − ∫ u T pdA ∫ 2V V Ap (3.23) Végeselem-módszer alkalmazásakor az elemekre bontott tartományokon
lokálisan approximált elmozdulással fejezzük ki a potenciális energiát: ne { ( )} Π p = ∑ Π ep u e (3.24) e =1 ahol az elemenkénti approximációval kifejezett potenciális energia ( ) Π ep u e = 1 eT e ε D ε dV − ∫ u eT ρk e dV − ∫ u eT p e dA . ∫ 2 Ve Ve Ae (3.25) p Végül a potenciális energia minimum elvből határozható meg az elmozdulásmező ismeretlen paraméterei (csomóponti elmozdulások). A következőkben bemutatjuk egy elem potenciális energiájának az előállítását, azaz a merevségi mátrix és tehervektor származtatását 54 3.26 SÍKBELI NÉGY CSOMÓPONTÚ IZOPARAMETRIKUS ELEM y η 4 η 4 +1 3 -1 P +1 -1 2 v3 3 u3 P’ ξ ξ leképzés 1 1 2 x 20. ábra A globális (x,y) és a természetes ( ξ , η ) lokális koordinátarendszer közötti leképzés A leképzés előnye elsősorban abban jelentkezik, hogy a természetes koordinátarendszer feletti integrálásra létezik könnyen programozható
numerikus algoritmus. A geometria leképzése: 4 4 x(ξ , η ) = ∑ N i (ξ , η )xi y (ξ , η ) = ∑ N i (ξ , η ) yi i =1 (3.26) i =1 ahol az alakfüggvények N1 (ξ ,η ) = 1 (1 − ξ )(1 − η ) , 4 1 N 3 (ξ ,η ) = (1 + ξ )(1 + η ) , 4 Az elmozdulás közelítése: N 2 (ξ ,η ) = 1 (1 + ξ )(1 − η ) , 4 1 N 4 (ξ ,η ) = (1 − ξ )(1 + η ) , 4 Az izoparametrikus elnevezésből következően ugyan azokat az alakfüggvényeket alkalmazzuk az elmozdulási mező közelítésére is, mint a geometria leképzésére 4 u (ξ , η ) = ∑ N i (ξ , η )ui , i =1 4 v(ξ , η ) = ∑ N i (ξ , η )vi . (3.27) i =1 A (3.27) felírható mátrixos alakban is N u e (ξ ,η ) = Nq e = 1 0 0 N1 N2 0 55 0 N2 N3 0 0 N3 N4 0 u1 v 1 u 2 0 v2 (3.28) N 4 u3 v3 u 4 v4 Az Ni alakfüggvényeket szokás interpolációs függvényeknek is nevezni,
szemleltetésüket axonometrikus ábrán mutatjuk meg. Az alakfüggvények: 1 1 1 1 N 3 (ξ ,η ) 21. ábra Az alakfüggvények N1 (ξ ,η ) N 2 (ξ ,η ) N 4 (ξ ,η ) u4 u3 u2 u1 22. ábra: 4 Az u (ξ ,η ) = ∑ N i (ξ ,η )ui szemléltetése i =1 Az alakváltozási jellemzők számítása: ∂ 4 ∂u N i (ξ ,η )ui ∑ ∂x i =1 ε x ∂x ∂ 4 ∂v e = ε = ε y = ∑ N i (ξ ,η )vi ∂y ∂y i =1 γ xy ∂u ∂v ∂ 4 ∂ 4 + ∑ N i (ξ ,η )ui + ∑ N i (ξ ,η )vi ∂x i =1 ∂y ∂x ∂y i =1 (3.29) Felírhatjuk (3.29)–et az alakfüggvények x,y szerinti deriválását alsóindexben történő jelölésének bevezetésével mátrixos alakban is 56 N1 x ε = B(ξ , η )q = 0 N1 y e e 0 N1 y N1 x N2x 0 N2y N3x 0 N2y 0 N3 y 0 N3y N2x N3x N4x 0 N4y
u1 v 1 u 0 2 v N 4 y 2 (3.30) u N 4 x 3 v3 u 4 v4 ahol B az úgynevezett elmozdulás-alakváltozási mátrix. Láthatjuk, hogy(3.29) és (330)-ban az alakfüggvények x,y szerinti deriváltjai szerepelnek, ez a deriválás a láncszabály alkalmazásával hajtható végre ∂N i (ξ ,η ) ∂N i (ξ ,η ) ∂ξ + ∂N i (ξ ,η ) ∂η ∂ξ ∂η ∂N i (ξ ,η ) ∂x ∂ξ ∂x ∂η ∂x ∂x ∂x ∂ξ ∂N (ξ ,η ) = ∂N (ξ ,η ) ∂ξ ∂N (ξ ,η ) ∂η = ∂ξ ∂η ∂N (ξ ,η ) i i i + i y y ∂ ∂ y ∂ ∂y ∂η ∂y 14243 ∂η ∂ξ (3.31) J −1 Közvetlenül a Jacobi mátrix inverze nem állítható elő, de a Jacobi mátrix igen ∂x ∂ξ J (ξ ,η ) =
∂x ∂η ∂y ∂ 4 ∑ N i (ξ ,η )xi ∂ξ ∂ξ i =1 = ∂y ∂ 4 N i (ξ ,η )xi ∂η ∂η ∑ i =1 ∂ 4 N i (ξ ,η ) yi ∑ ∂ξ i =1 . ∂ 4 ∑ N i (ξ ,η )yi ∂η i =1 (3.32) A Jacobi mátrix ismeretében az inverz képzést alkalmazzuk a numerikus számításoknál J −1 (ξ ,η ) = adj (J ) . det(J ) (3.33) A Jacobi mátrix determinánsa igen fontos szerepet játszik az elem leképzésének ellenőrzése során is, ha az értéke zérus, akkor a leképzés szinguláris. A gyakorlati számítás során azt érzékeljük, hogy különböző pontokban előjelet vált ilyenkor az elem nagyon eltorzult alakot mutat. 57 A feszültség származtatása: σ e (ξ ,η ) = Dε(ξ ,η ) = DB e (ξ ,η )q e (3.34) Az elem potenciális energiája (3.25) az alakváltozási energiából és a külső erők potenciáljából (vagy külsőerők munkájából) áll Nézzük először az
alakváltozási energiát : 1 1 1 eT e 1 1 ε D ε dV = q eT ∫ B T D B e bdxdyq e = q eT ∫ ∫ B T D B e b det(J )dξdηq e ∫ 2 Ve 2 2 −1 −1 Ae (3.35) Az integrálban a dV térfogat elem tárcsa b vastagsága és az újváltozóra történő átalakítás miatt a Jacobi determináns segítségével van felírva: dV = bdxdy = b det( J )dξdη A térfogati erő munkája: 1 1 eT e eT T e eT T e ∫ u ρk dV = q ∫ N ρk bdxdy = q ∫ ∫ N ρk b det(J )dξdη Ve (3.36) −1 −1 Ae A megoszló terhelés munkája: η y p dΓ 3 dy dx 4 1 2 dΓ = (dx ) 2 2 + (dy ) 2 2 dx dy = + dξ dξ dξ x 23. ábra 2 2 dx dy ∫e u p dA = q ∫Γ N (ξ ,η = 1) p bdΓ = q −∫1 N (ξ ,η = 1) p b dξ + dξ dξ Ap 1 eT e eT T e eT T e (3.37) 58 0 ahol p = most csak függőleges irányú megoszló terhelést tartalmazza.
py A (3.35) és (436) felületi integrálokat tartalmaznak a kétegység élű négyzet felett, míg (3.37) csupán egy vonalintegrál a [-1,1] tartományon Mind a kéttípusú integrálra alkalmazható az ún Gauss-féle kvadratúra Numerikus integrálás, (Gauss kvadratúra): α (ξ ) ξ -1 ξ1 ξ2 1 24. ábra Vonalintegrál 1 NG −1 i =1 ∫α (ξ ) dξ = ∑α (ξ i ) Wi ahol ξ i ,Wi a Gauss koordináták és a megfelelő Gauss súlyok. Polinomok esetén (2*NG-1) rendig pontos értéket szolgáltat. 4.1 Táblázat: Gauss koordináták és Gauss súlyok ξi NG 1 2 3 Wi 0 -0.57735 02691 89626 0.57735 02691 89626 -0.77459 66692 41483 0.0 0.77459 66692 41483 2.0 1.0 1.0 0.55555 55555 55555 0.88888 88888 88888 0.55555 55555 55555 Felületi integrál: 1 1 NG NG −1 −1 i =1 j =1 ∫ ∫ β (ξ , µ ) dξdη = ∑∑ β (ξ , µ ) Wi W j 59 Merevségi mátrix: Ke = T e T e ∫ ∫ B D B b det(J )dξdη = ∑∑ B (ξ i ,η j )D B (ξ i ,η j )b det(J
(ξ i ,η j )) WiW j 1 1 NG NG −1 −1 i =1 j =1 (3.38) Tehervektor a térfogati terhelésből: f ρ = ∫ ∫ N T ρk e b det(J )dξdη = ∑∑ N T (ξ i ,η j )ρk e (ξ i ,η j )b det (J (ξ i ,η j ))WiW j 1 1 NG NG −1 −1 i =1 j =1 e (3.39) Tehervektor a felületi terhelésből 2 2 NG dx dy f = ∫ N (ξ ,η = 1) p b + dξ = ∑ N T (ξ i ,η = 1) p e b i =1 dξ dξ −1 1 e p T e 2 2 dx dy + Wi dξ dξ (3.40) 3.27 PÉLDA TORZULT SÍKBELI ELEM ELFAJULÓ LEKÉPZÉSÉRE: η = −ξ + 2 3 4 3 1 2 25. ábra:Torzult elem Az ábrán vázolt elem sarok pontjainak koordinátái adottak: 60 Csomópont 1 x y 0 0 2 4 0 3 1 1 4 0 4 Az elem leképzése: 4 x(ξ ,η ) = ∑ N i (ξ ,η )xi = 1 (1 + ξ )(1 − η )4 + 1 (1 + ξ )(1 + η )1 = 1 (1 + ξ )(5 − 3η ) 4 4 4 i =1 4 1 1 1 y (ξ ,η ) = ∑ N i (ξ ,η ) yi = (1 +
ξ )(1 + η )1 + (1 − ξ )(1 + η )4 = (1 + η )(5 − 3ξ ) 4 4 4 i =1 A Jacobi mátrix ∂x ∂y (5 − 3η ) − 3(1 + η ) ∂ξ ∂ξ 4 4 J (ξ ,η ) = = ( ) ( ) 5 − 3ξ ∂x ∂y − 3 1 + ξ 4 4 ∂η ∂η A determináns: (5 − 3η ) − 3(1 + η ) 25 − 15ξ − 15η + 9ξη 9 + 9ξ + 9η + 9ξη 4 4 det[J (ξ ,η )] = 0 = = − = − 3(1 + ξ ) (5 − 3ξ ) 4 4 4 4 16 − 24ξ − 24η = = 4 − 6ξ − 6η = 0 4 Azon pontok mértani helye ahol a leképzés szinguláris: 61 η = −ξ + 2 3 4 DINAMIKAI FELADATOK VIZSGÁLATA VÉGESELEMMÓDSZERREL A vizsgálatainkat 1D-s esetre mutatjuk be. Az 1 Fejezetben ismertetett statikai peremérték feladatot (1.a és 1b ábrák) úgy bővítjük ki, hogy figyelembe vesszük a tömegerőket és az erőhatások időbeli lefutását is. Így egy dinamikai feladatot kapunk, amely kezdeti érték feladatként definiálható. Ez azt jelenti,
hogy a (3) egyensúlyi egyenlet jobboldalán az impulzus tétel alapján megjelenik az inercia erő. A többi egyenlet megismétlődik oly módon, hogy az idő implicite megjelenik. Erős megfogalmazás Kinematikai (geometriai) egyenlet: ε z (z, t ) = dw(z , t ) dz 0< z<l, 0 < t < t∞ (4.1) ahol t az idő, és t ∞ a vizsgálat időtartama. Anyagegyenlet (egyszerű Hooke-törvény): σ z = Eε z ⇒ N ( z , t ) = A σ z = AEε z 0 < z < l , 0 < t < t∞ (4.2) Impulzus tétel: dN ( z , t ) &&( z , t ) + p z ( z , t ) = Aρw dz 0< z<l. 0 < t < t∞ (4.3) && a rúdirányú gyorsulás. ahol w Kinematikai peremfeltétel: w(0,t)=0. Dinamikai peremfeltétel: 62 0 < t < t∞ (4.4) N(l,t)= Fz (t). Kezdeti feltételek: 0 < t < t∞ w( z ,0) = w0 ( z ) w& ( z ,0) = w& 0 ( z ) t =0 t =0 (4.5) (4.6) (4.7) A gyengealak származtatása valójában megegyezik az 1. Fejezetben leírtakkal A
d’Alambert elv értelmében az inercia erő –1 szeresét mint külső erőt veszzük figyelembe és így a p z megoszló erővel együtt szerepeltetjük dN ( z , t ) &&( z , t )] = 0 . + [ p z ( z , t ) − Aρw dz Hangsúlyozzuk, hogy az elmozdulás variációja δw( z ) időtől független mennyiség. (4.8) Végül is a (17) gyengealak egyszerűen átírható d (δw) dw ∫0 dz ⋅ AE ⋅ dz dz = ∫0 δw ⋅ [ p z − Aρw&&]dz + Fz ⋅ δwl l l (4.9) Gyenge megfogalmazás Gyengealak: d (δw) dw ∫0 dz ⋅ AE ⋅ dz dz + ∫0 δw ⋅ Aρw&&dz = ∫0 δw ⋅ p z dz + Fz ⋅ δwl 0 < t < t∞ l l l (4.9) Kezdeti feltételek: l l ∫ δw ⋅ Aρw(z,0)dz = ∫ δw ⋅ Aρw0 (z )dz 0 0 l l 0 0 ∫ δw ⋅ Aρw& (z,0)dz = ∫ δw ⋅ Aρw& (z )dz 0 63 t =0 (4.10) t =0 (4.11) 4.1 VÉGESELEM MÁTRIXOK ÉS A TEHERVEKTOR: Az első fejezetben hangsúlyoztuk, hogy a virtuális munka elv rugalmas feladatok esetén
egyenértékű a teljes potenciális energia minimum elvvel. Ez azt jelenti, hogy a végelemmódszer során előálló merevségi mátrix és tehervektor végül is független attól, hogy melyik elvből származtatjuk Így pl (49) baloldalán az első integrálból áll elő a merevségi mátrix és a jobboldalon szintén az első integrálból, pedig a tehervektor származtatható. Egy szerkezet esetén a vizsgált tartományt elemekre bontjuk és előállítjuk elemenként a megfelelő integrál mennyiségeket. Egy általános húzott-nyomott rúdelemen ((66) alapján) az elmozdulás approximáció ζ ζ wi ζ ζ (4.12) ⋅ , we (ζ , t ) = 1 − wi (t ) + w j (t ) = 1 − L L L L w j a gyorsulás is hasonlóan írható fel &&i ζ ζ w ζ ζ && e (ζ , t ) = 1 − w &&i (t ) + w && j (t ) = 1 − ⋅ , w
&& j L L L L w (4.13) az elmozdulás variációja azonban nem függ az időtől ζ ζ δw e (ζ ) = 1 − δwi + δw j = [δwi L L Az elem merevségi mátrixa: 64 ζ 1 − δw j ⋅ L . ζ L ] (4.14) L ∫ 0 d (δw) dw ⋅ AE ⋅ ζ dζ = dζ dz − 1 − 1 1 wi = ∫ δwi δw j L AE w dζ = 1 L L j 0 L AE AE − L L wi = δq eT K e q e = δwi δw j AE AE w j − L L L [ ] [ (4.15) ] Az elem tehervektora: ∫ δw ⋅ pz dζ == ∫ [δwi L L 0 0 ζ 1 − δw j L p z dζ = δwi ζ L ] [ pz L 2 = δq eT f e (4.16) δw j p p L z 2 ] Az elem tömegmátrixa: ζ
&&i 1 − L ζ ζ w Aρ 1 − & & ⋅ = ⋅ w A w d w w ⋅ dζ = δ ρ ζ δ δ j ∫0 ∫0 i && j ζ L L w L AρL AρL && 6 wi = δq eT M e q && e = δwi δw j 3 AρL AρL w && j 3 6 L L [ [ ] ] 65 (4.17) 4.2 DISZKRETIZÁLT MOZGÁSEGYENLET Egy összetett szerkezetnél, mint a teljes potenciális energiánál összegezni kell a megfelelő mennyiségeket, végül is a szerkezetre nézve elő kell állítani a merevségi mátrixot, a tömegmátrixot és a tehervektort, azaz (4.9) átírható && + δq T Kq = δq T f (t ) δq T M q (4.18) A megfelelő kinematikai peremfeltételek figyelembevétele után a tetszőleges variációból következően az alábbi mozgásegyenletet kapjuk: && + Kq = f (t ) Mq 66 (4.19) 4.3
ÁLLANDÓSULT HARMONIKUS GERJESZTÉS A (4.19) egyenlet jobboldalán a gerjesztő erő időben harmonikus függvény szerint változik f (t ) = f 0 cos ωt . (4.20) Feltételezzük, hogy a gerjesztő erőhatás tartósan működik, vagyis (4.19) mozgásegyenlet homogén része a valóságban mindig jelenlévő csillapítás következtében lecseng Keressük (4.19) mozgásegyenlet partikuláris megoldását a (420) zavaró függvény alakjában: q(t ) = q 0 cos ωt . (4.21) &&(t ) = −ω 2q 0 cos ωt . q (4.22) Időben kétszer deriválva Visszahelyettesítve (4.20)-(421)-et (419)-be − ω 2 Mq 0 cos ωt + Kq 0 cos ωt = f 0 cos ωt (4.23) (K − ω M )q (4.24) átrendezés után kapjuk, hogy 2 0 cos ωt = f 0 cos ωt . A cos ωt diszkrét pontok kivételével különbözik zérustól, ezért egyszerűsíthetünk vele és egy lineáris algebrai egyenletrendszert kapunk az elmozdulás amplitúdóra. (K − ω M )q = f , (4.25) ahol az együttható mátrixot (K
− ω M ) = K szokás dinamikai merevségi mátrixnak is 2 0 0 2 D nevezni. Az elmozdulás amplitúdó nem csak a gerjesztő erő amplitúdójától, hanem a gerjesztés ω körfrekvenciájától is függ ( q 0 = K − ω 2M ) 67 −1 f0 . (4.26) 4.4 SAJÁTREZGÉSEK Tekintsünk el a (4.19) differenciálegyenletet jobboldalán lévő gerjesztő erőhatástól && + Kq = 0 Mq (4.27) az így kapott (4.27) homogén differenciálegyenlet megoldását időben harmonikus függvény alakjában keressük q(t ) = φ cosαt . (4.28) &&(t ) = −α 2 φ cosαt . q (4.29) Idő szerint kétszer deriválva Visszahelyettesítve (4.28)-(429)-et (427)-be − α 2 Mφ cosαt + Kφ cosαt = 0 (4.30) (K − α M )φ cosαt = 0 . (4.31) átrendezés után kapjuk, hogy 2 A cosαt diszkrét pontok kivételével különbözik zérustól, ezért egyszerűsíthetünk vele és egy lineáris homogén algebrai egyenletrendszert kapunk az elmozdulás amplitúdóra (K
− α M )φ = 0 . 2 (4.32) (4.32) egy sajátértékfeladat, ahol az α ún sajátkörfrekvencia nem ismert A homogén lineáris algebrai egyenletrendszernek akkor van triviálistól különböző megoldása, ha az együttható mátrix determinánsa zérus ( ) ( ) det K − α 2 M = 0 = polinom α 2 . (4.33) A (4.33)-ban kifejtett determináns α 2 -re nézve egy n-ed fokú polinomot ad, amelyet karakterisztikus polinomnak nevezünk ( K és M mérete =nxn) Bizonyítható, hogy 68 (4.33) gyökei α 2 -re nézve valósak és egyenlő vagy nagyobbak zérusnál Végül is (432) és (4.33) megoldásai n számú sajátértékpárt szolgáltat ( α 12 , φ1 ), ( α 22 , φ 2 ), ( α n2 , φ n ) ahol φ i az i-edik sajátkörfrekvenciához tartozó sajátvektor. A sajátérték feladatot gyakran az alábbi alakban szokás felírni Kφ i = α i2 Mφ i . A gyakorlatban azért fontos a sajátkörfrekvencia ismerete, mert ezen a frekvencián gerjesztve a rezgőrendszert az
elmozdulás amplitúdó rendkívül megnő, csillapítás mentes esetben elvileg tart a végtelenbe. A sajátvektor vagy más néven rezgéskép ismerete segíti a tervezőt abban, hogy esetleges pótlólagos megfogásokat alkalmazva a nagy amplitúdójú pontokban a rendszer eredményesen elhangolható legyen egy nem kívánt sajátfrekvenciáról. 69 5 HŐFESZÜLTSÉGEK SZÁMÍTÁSA Ha egy test hőtágulását nem korlátozzuk, akkor a konstans vagy a hely lineáris függvénye szerinti hőmérséklet eloszlásból nem származik hőfeszültség. Ellenkező esetben, ha egy test hőtágulását korlátozzuk, akkor a testben feszültség ébred még konstans hőmérséklet eloszlás esetén is. Az első megállapításból következik, hogy a feszültség számításánál figyelembe kell venni a hőtágulás okozta alakváltozást. 70 5.1 HŐMÉRSÉKLET HATÁSÁNAK FIGYELEMBEVÉTELE 1D-S FELADATNÁL: Húzott nyomott rudak esetén a hőtágulással módosított
egyszerűsített Hooke-törvény: σ z = E (ε z − α∆T ) , (5.1) ahol α a fajlagos hőtágulási együttható, ∆T a hőmérséklet növekmény. Megszorozva (51) mindkét oldalát a rúd A keresztmetszetével kapjuk a rúderőre vonatkozó anyagtörvényt dw − α∆T . N = AE (ε z − α∆T ) = AE dz (5.2) Az (5.1) és (52) képletekben a tényleges alakváltozásból levontuk a hőtágulás okozta fajlagos nyúlást, ezért a potenciális energia felírásánál is hasonlóan járunk el: Π ep (wi , w j ) = L L dw e 1 dw e − ∆ ⋅ ⋅ − ∆ α T AE α T d ζ − w e ⋅ p z dζ ∫ dζ 2 ∫0 dζ 0 (5.3) elvégezve a kijelölt szorzást L L L dw e dw e 1 dw e e ⋅ AE ⋅ dζ − ∫ w ⋅ p z dζ − ∫ ⋅ AE ⋅ α∆T dζ Π (wi , w j ) = ∫ 2 0 dζ dζ dζ 0 0 e p L + 1 2 AE
⋅ (α∆T ) dζ ∫ 20 (5.4) Az utolsótag (5.4)-ben nem függ a w elmozdulástól, így az a variációképzés szempontjából konstansnak tekinthető, ezért el is hagyjuk: 71 L L L dw e dw e 1 dw e e ⋅ AE ⋅ dζ − ∫ w ⋅ p z dζ − ∫ ⋅ AE ⋅α∆T dζ Π (wi , w j ) = ∫ dζ 2 0 dζ dζ 0 0 e p (5.5) A korábban a 2. fejezetben ismertetett potenciális energia most egy plusz taggal bővült, ebből a tagból a végeselemre nézve egy ún. hőtehervektor származtatható Amint azt láttuk (67) alatt az elmozdulás deriváltja a következő módon számolható dw (ζ ) = wi dζ [ e − 1 wj L 1 L ] (5.6) ezt behelyettesítve (5.5) utolsó integráljába és elvégezve a kijelölt műveletet megkapjuk a hőtehervektort dw e ∫0 dζ ⋅ AE ⋅α∆T dζ = wi [ L − 1 w j ∫ L
⋅ AE ⋅ α∆T dζ = wi 1 0 L ] [ L − AEα∆T wj AEα∆T ] (5.7) Végül is az elem (77) potenciális energiája kissé módosul a terhelési vektorban Π ep = [ 1 wi 2 AE wj L AE − L ] AE L wi − w i AE w j L − [ 72 w j ] pz L − AEα∆T 2 . pz L + AEα∆T 2 (5.8) 5.2 HŐMÉRSÉKLET HATÁSÁNAK FIGYELEMBEVÉTELE 2D-S FELADATNÁL: Az 1D-s feladathoz hasonlóan módosul a (3.19) Hooke-törvény, 2D-s feladatoknál is a fajlagos nyúlásokból le kell vonni a hőtágulás okozta fajlagos nyúlást: 1 ν 0 ε x α∆T E ν 1 0 ε y − α∆T σ = D(ε − ε 0 ) = 2 1 −ν 1 −ν γ xy 0 0 0 2 (5.9) Az elem
(3.25) potenciális energiáját is hasonlóan módosítjuk ( ) Π ep u e = ( ) ( ) 1 ε eT − ε T0 D ε e − ε 0 dV − ∫ u eT ρk e dV − ∫ u eT p e dA 2 V∫e Ve Aep (5.10) a kijelölt szorzást elvégezve ( ) Π ep u e = 1 eT e 1 ε D ε dV − ∫ u eT ρk e dV − ∫ u eT p e dA − ∫ ε eT D ε 0 dV + ∫ ε T0 D ε 0 dV (5.11) ∫ 2 Ve 2 Ve Ve Aep Ve az utolsó integrál most is független az elmozdulástól ezért a variációképzés szempontjából ez a konstansak tekinthető ezért el is hagyjuk Π ep (u e ) = 1 eT e ε D ε dV − ∫ u eT ρk e dV − ∫ u eT p e dA − ∫ ε eT D ε 0 dV . (512) ∫ 2 Ve Ve Ae Ve p Az (5.12) utolsó integráljából most is egy ún hőtehervektor származtatható felhasználva (4-30)-at ∫ε Ve 1 1 eT D ε 0 dV = q eT ∫ B (ξ ,η )D ε bdxdy = q ∫ ∫ B (ξ ,η ) D ε b det(J )dξdη = q T eT T 0 Ve 0 eT fαe −1 −1 ahol fαe az elem hőteher vektora, numerikusan a Gauss kvadratúrával
számolható ki. 73 fαe = T T ∫ ∫ B (ξ ,η ) D ε 0b det(J )dξdη = ∑∑ B (ξ i ,η j ) D ε 0 (ξ i ,η j )b det(J (ξ i ,η j ))WiW j 1 1 NG NG −1 −1 i =1 j =1 A hőtágulásból származó fajlagos nyúlás is lehet a hely függvénye, ha a hőmérséklet is függ a helytől 4 α ∑ N i (ξ ,η )∆Ti α∆T (ξ ,η ) i =41 ε 0 (ξ ,η ) = α∆T (ξ ,η ) = α ∑ N i (ξ ,η )∆Ti i =1 0 0 74 6 AZ IDEAS PROGRAMRENDSZER ALKALMAZÁSA VÉGESELEMES FELADATOK MEGOLDÁSÁRA Ebben a fejezetben három példát mutatunk be a végeselem-módszer alkalmazására. Az első feladat 1D-s elemet azaz rúdelemet alkalmaz a modellezéshez. A rudak nemcsak húzva hanem hajlítva és nyírva is vannak Ennek megfelelően a csomópontokban a három koordináta irányú elmozdulási paraméter mellett a három kordináta tengely irányú szögelfordulás is megjelenik
mint ismeretlen Ekkor a csomópontokban, a koordináta irányokban az erőkön kívül erőpárok (nyomaték) is működtethető. A második feladat 2D-s elemeket alkalmaz a tárcsfeladat, azaz az általánosított síkfeszültségi állapotú peremérték probléma modellezéséhez. Csomópontonként ebben az esetben a két koordináta irányú elmozdúlás paraméter jelentkezik, mint ismeretlen. Ennél a feladatnál a csomópontokban csak síkbeli erőhatások alkalmazhatók, de nyomatékok nem Harmadik esetben 3D-s elemekkel modellezzük a kitűzött peremérték feladatot. Itt is a csomópontokban csak koordináta irányú elmozdulások az ismeretlenek, így terhelés csak erő jellegű lehet. 75 6.1 RÚDSZERKEZET VIZSGÁLATA A szerkezet geometriája A körszelvényű rudak (D100/d90, D50/d40), F=1000N, a megtámasztás csuklós és görgős. 76 6.11 RÚDFELADAT MEGOLDÁSÁHOZ ALKALMAZOTT IDEAS UTASITÁSOK Ideas Open GL / Model file name: /Application: /Task : lev01
Simulation Master modeler Option /Units / mm(newton) MASTER MODELLER: B(2,3) Wokplace Apperarence (-2000, -2000) ( 2000, 2000) C(1,1) Refresh A(2,1) Line A(4,1) Dimension B(2,1) Modify . BEAM SECTION A(1,3) Circular beam A(5,2) Directory . MESHING A(1,1) B(1,1) B(1,3) . A(2,2) (100/90) CSO1 Define beam mesh/ Element length / Beam options / / @@ Labels (Nodes, Elements) (Label size) (1000) (CSO1/CSO2) Coincident Nodes BOUNDARY CONDITION A(4,2) Displ. Restraint/ A(2,1) . (50/40) (CSO2) Pin (z) Roller (x-trans, z-rotation) (-2500) Force/ Y-Force MODEL SOLUTION A(1,2) Create/ A(2,1) Solve . POSTPROCESSING A(1,1) A(1,2) A(2,1) Display Options (Done) 77 Element forces/ 6.2 TÁRCSAFELADAT A szerkezet geometriája A szerkezet terhelése Total force 1500N, merev megfogás a furat mentén, vastagság 2mm. 78 6.21 TÁRCSAFELADAT MEGOLDÁSÁHOZ ALKALMAZOTT IDEAS UTASITÁSOK Ideas Open GL Option / Model file name: /Application: /Task : lev02 Simulation Master modeler
/Units / mm(newton) MASTER MODELLER: B(2,3) Wokplace Apperarence C(1,1) A(2,1) A(4,1) B(2,1) . A(4,2) . A(3,1) A(5,1) A(4,3) (-100, -100) ( 100, 100) Refresh Line Dimension Modify Filet (R5) Circle Surface by boundary Trim at curve BOUNDARY CONDITION A(4,2) Displacement Restraint A(2,1) Total force / Inplane (1500N) MESHING A(7,2) A(1,1) Surface thickness (2) Def. Shell mesh /Length (1) /Element family (Plane stress) @@ (megtekintes) MODEL SOLUTION A(1,2) Create A(2,1) Solve . POSTPROCESSING A(1,1) A(1,2) A(2,1) Display (Done) 79 6.3 TÉRBELI FELADAT A tengelycsonk geometriája A tengely 10000N megoszló erővel terhelt, a hátsó síklap mereven megfogott 80 6.31 TÉRBELI FELADAT MEGOLDÁSÁHOZ ALKALMAZOTT UTASITÁSOK Ideas Open GL / Model file name: /Application: /Task : lev01 Simulation Master modeler Option /Units / mm(newton) MASTER MODELLER: B(2,3) Wokplace Apperarence (-100, -100) ( 100, 100) C(1,1) Refresh A(2,1) Rectangle by two corners C(3,2)
Axiális nézet A(5,1) A(1,1) A(3,1) A(5,1) A(5,2) Extrude (100) Sketch in place Circle (D50) Extrude (200) Filet (5) A(2,1) C(3,2) A(5,1) Line Axiális nézet Extrude (300) Protrude BOUNDARY CONDITION A(4,2) Displacement Restraint A(2,1) Total force / Vector Inplane (10000N) MESHING A(1,1) Solid mesh /Length (10) /Element family (Linear) @@ (megtekintes) MODEL SOLUTION A(1,2) Create /option / Iterative A(2,1) Solve . POSTPROCESSING A(1,1) A(1,2) A(2,1) Display (Done) 81 IDEAS