Univerza v Ljubljani Fakulteta za gradbeništvo ii geodezijo PODIPLOMSKI ŠTVDIJ GRADBENIŠTVA KONSTRUKCIJSKA SMER DOKTORSKI ŠTUDU Kandidat: Matija Gams, univ.dipl.inž.grad. GEOMETRIJSKO TOČNA DINAMIČNA ANALIZA ELASTIČNIH IN ARMIRANOBETONSKIH OKVIRNIH KONSTRUKCIJ Doktorska disertacija štev.: 178 GEOMETRICALLY EXACT DYNAMIC ANALYSIS OF ELASTIC AND REINFORCED CONCRETE FRAME STRUCTURES Doctoral thesis "No.: 178 Temo doktorske disertacije je odobril Senat Univerze v Ljubljani ha svoji 7. seji dne 27. junija 2006 in imenoval mentorja prof.dr. Mirana Sajeta in somentorja izr.prof.dr. Igorja Planinca Ljubljana, 14. marec 2008 Y Vili jh » * i f i iTi] r iS i i-M I il) I ll I yi|^ui^ Univerza v Ljubljani Fakulteta za gradbenišovo m iiitti, rn.iTn Ini ili i ilii_n| liiTîiimiiiîi Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi prof.dr. Miran Saje izr.prof.dr. Igor Planine doc.dr. Marko Kegl, UM FS doc.dr. Gordan Jelenić, Građevinski fakultet, Sveučilište u Rijeci, Hrvaška je imenoval Senat Fakultete za gradbeništvo in geodezijo na 9. redni seji dne 19. aprila 2006. Komisijo za oceno doktorske disertacije v sestavi izr.prof.dr. Jože Korelc doc.dr. Marko Kegl, UM FS doc.dr. Gordan Jelenić, Građevinski fakultet, Sveučilište u Rijeci, Hrvaška je imenoval Senat Fakultete za gradbeništvo in geodezijo na 14. redni seji dne 30. januarja 2008. Komisijo za zagovor doktorske disertacije v sestavi prof.dr. Bojan Majes, dekan, predsednik prof.dr. Miran Saje izr.prof.dr. Igor Planine izr.prof.dr. Jože Korelc doc.dr. Marko Kegl, UM FS doc.dr. Gordan Jelenić, Građevinski fakultet, Sveučilište u Rijeci, Hrvaška je imenoval Senat Fakultete za gradbeništvo in geodezijo na 15. redni seji dne^-^jj^ 27. februarja 2008 ^^-----^ Univerza v Ljubljani Fakuketaza jsr* m 11 t Tri r 1 *j mu li i t il i l ih j_ij_ iAJüLLLLl IZJAVA O AVTORSTVU Podpisani MATIJA GAMS, univ.dipl.inž.grad., izjavljam, da sem avtor doktorske disertacije z naslovom: »GEOMETRIJSKO TOČNA DINAMIČNA ANALIZA ELASTIČNIH IN ARMIRANOBETONSKIH OKVIRNIH KONSTRUKCIJ«. Ljubljana, H.marec 2008 (podpii Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij V Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. BIBLIOGRAFSKO-DOKUMENTACIJSKA STRAN IN IZVLEČEK UDK: 531.3:624.012.45(043.2) Avtor: Matija Gams Mentor: prof. dr. Miran Saje Somentor: izr. prof. dr. Igor Planine Naslov: Geometrijsko točna dinamična analiza elastičnih in armiranobetonskih okvir- nih konstrukcij Obseg in oprema: 194 str., 9 pregi., 37 si., 267 en. Ključne besede: točna kinematika, Reissnerjev nosilec, dinamika, časovna integracija, optimizacija dinamičnih sistemov, armirani beton Izvleček V disertaciji se ukvarjamo z dinamiko kinematično točnih ravninskih elastičnih in armiranobetonskih konstrukcij. Osrednja tema je upoštevanje točnih (Reissnerjevih) kinematičnih zvez v problemih, povezanih z dinamiko ravninskih nosilcev z upoštevanjem osnih, strižnih in upogibnih deformacij. Znotraj te (širše) teme je disertacija razdeljena na štiri podteme. V sklopu prve rešujemo problem dinamike elastičnih nosilcev na način, kjer kot edine neznane funkcije upoštevamo deformacijske količine. Izkaže se, da direkten pristop numerično ni najbolj učinkovit, saj v njem nastopa gnezdenje integralov in zato razvijemo izboljšan numerični algoritem za reševanje takih problemov. Lastnosti osnovne in izboljšane formulacije preverimo na numeričnih primerih. V nadaljevanju disertacije se lotimo vprašanja časovne integracije, pri čemer se omejimo na končne elemente, v katerih so interpolirane le kinematične količine (pomiki in zasuki). V sklopu te teme predstavimo klasični pristop k integraciji, zatem pa temeljito preučimo ti. energijski pristop. Energijski pristop izhaja iz pogoja ohranjanja mehanske energije, ki je kriterij stabilnega reševanja. Med analizo energijskih integratorjev izpeljemo še nov časovni integrator, ki temelji na upoštevanju šibkih (časovno odvajanih) kinematičnih enačb. Za vse analizirane integratorje analitično in numerično preverimo ohranjevalne lastnosti energije in ostalih količin gibanja. Znotraj te teme v integratorje vgradimo še numerično dušenje, ki se lahko s specialnim novim algoritmom samodejno vklaplja in izklaplja. Pokaže se, da tako dušenje manj prizadene nihanje v osnovnih nihajnih oblikah, kar preverimo tudi z numeričnimi primeri. Tretja tema disertacije se ukvarja z optimizacijo elastičnih nosilcev v dinamiki. Prispevek v sklopu te teme predstavlja razvoj in aplikacija občutljivostne analize za novo razviti šibki časovni integrator. Razvito orodje omogoča reševanje izrednega nabora problemov iz optimizacije dinamičnih sistemov, kjer se lahko optimira tako oblika in odpornost konstrukcij na dinamično obtežbo, kot tudi operacije in geometrija strojnih manipulatorjev. V zaključnem delu disertacije se ukvarjamo s kinematično točno dinamiko armiranobetonskih konstrukcij. Razviti računski postopek verficiramo s primerjavo s programom Opensees in validiramo na eksperimentalnih rezultatih. Razvoj večine numeričnih postopkov v disertaciji temelji na programskih sistemih AceGen in AceFEM, ki omogočata razvoj programske kode končnih elementov v simbolnem zapisu in optimizacijo generirane kode ter zmogljivo vsestransko okolje za preverjanje in uporabo analiz s končnimi elementi. VI Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION UDC 531.3:624.012.45(043.2) Author: Matija Gams Supervisor: prof. dr. Miran Saje Co-supervisor: Assoc. prof. dr. Igor Planinc Title: Geometrically exact dynamic analysis of elastic and reinforced concrete fra- me structures Notes: 194 p., 9 tab., 37 fig., 267 eq. Key words: exact kinematics, Reissner’s beam, dynamics, time integration, optimization of dynamic systems, reinforced concrete Summary The dissertation deals with dynamics of geometrically exact elastic and reinforced concrete planar structures. The main topic is the treatment of exact (Reissner’s) kinematic equations in problems concerning dynamics of planar beams. The kinematic equations include axial, bending as well as shear strains. The dissertation is divided into four themes. The first theme concerns solving the problem of dynamics of elastic planar beams by taking the strains as the only unknown functions. It turns out that the direct approach to solving problems, defined in this way, is numerically inefficient as it involves multiple nesting of integrals. An enhancement of the approach is therefore developed and tested on numerical examples. In the next theme we discuss the issue of time integration in kinematically exact dynamics. During the analysis we restrict ourselves to displacement and rotation based finite elements, which is the standard approach to solving problems in mechanics. The classical approach is presented first followed by the detailed energy approach. As a part of the analysis of energy conservation based integrators, a new time integration scheme is developed. It is based on the use of kinematic equations, which are differentiated with respect to time. Conservational properties of all analysed integration schemes are analytically and numerically tested. The schemes have been supplemented with numerical dissipation, which can be arbitrarily turned on or off according to a special algorithm in order to affect the lower modes of response as little as possible. This is also verified by the numerical examples. The third theme of the dissertation deals with the optimization of kinematically exact elastic structures. This part of the dissertation contributes to the development and the application of the sensitivity analysis for the newly developed time integration scheme. A wide array of problems dealing with optimization of dynamical systems can be solved. These problems include the optimization of shape and resistance of structures as well as the loading regimes and the optimal shape of mechanical manipulators. In the final part of the dissertation, the dynamics of kinematically exact reinforced concrete structures is discussed. A numerical procedure is developed, verified by the Opensees program and validated by experimental results. The majority of the numerical procedures presented in the dissertation have been developed in AceGen and AceFEM computer programs, through the symbolic programming of the finite element computer code and the expression optimization. These programs have been found to offer a versatile environment for testing and using finite element based analyses. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij VII Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Zahvale Od srca se zahvaljujem mentorju prof. Miranu Sajetu in izr prof. Igorju Planincu za strokovno vodstvo pri pripravi doktorske disertacije. Vajina predanost znanosti je vredna najvišjega spoštovanja! Posebna zahvala gre moji ženi Ivani, ki mije tekom študija neomajano stala ob strani. Za podporo med študijem se iskreno zahvaljujem svojim staršem. Hvaležen sem tudi vsem sodelavcem na Katedri za mehaniko za prijetno in prijateljsko vzdušje. Med njimi moram izpostaviti sodelavce, s katerimi sem si v času študija delil pisarno: Bojana Časa, Sebast-jana Bratino, Nano Krauberger, Tomaža Hozjana in Mojco Markovič. Za prijetne znastvene pogovore, ki so pripomogli pri nastanku tega dela sem hvaležen več ljudem, kot jih je možno napisati na to stran. Med njimi moram izpostaviti Miho Kramarja in Mateja Rozmana, ki sta mi pomagala bolj, kot bi od njiju lahko pričakoval. VIII Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. KAZALO VSEBINE 1 UVOD 1 1.1 Predstavitev problema in pregled stanja na področju dinamike ravninskih okvirnih kon- 2 OPIS PROBLEMA 9 2.3 Osnovni izreki in definicije................................. 11 2.4 Vodilne/glavne enačbe ................................... 13 2.5 Analiza zapisanih enačb .................................. 15 3 DINAMIKA LINEARNO ELASTIČNIH LINIJSKIH KONČNIH ELEMENTOV Z INTERPOLACIJO DEFORMACIJSKIH KOLIČIN 17 3.2 Metoda uteženih ostankov ................................. 19 3.3 Galerkinova metoda končnih elementov.......................... 20 3.4 Numerični vidik formulacije................................ 21 3.7 Izboljšava elementov z interpolacijo deformacij...................... 23 4 DINAMIKA LINEARNO ELASTIČNIH LINIJSKIH KONČNIH ELEMENTOV Z INTERPOLACIJO KINEMATIČNIH KOLIČIN 30 Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij IX Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. 4.2 Klasicˇen pristop ....................................... 30 4.3 Energijsko konzervativne sheme .............................. 34 4.3.1 Deformacijska energija ............................... 35 4.3.1.1 Algoritmi s krepkimi kinematicˇnimi vezmi ............... 36 4.3.1.2 Algoritmi s sˇibkimi kinematicˇnimi vezmi ................ 38 4.3.2 Kineticˇna energija ................................. 41 4.3.3 Obtezˇba in potencial zunanjih sil .......................... 43 4.3.4 Povzetek enacˇb ................................... 44 4.3.5 Ohranjanje gibalne kolicˇine ............................ 44 4.3.6 Ohranjanje vrtilne kolicˇine ............................. 45 4.4 Algoritmi .......................................... 47 4.5 Numericˇni primer ...................................... 49 4.6 Sklepi ............................................ 52 4.7 Dusˇenje ........................................... 53 4.8 Numericˇni primer ...................................... 56 4.9 Sklepi ............................................ 61 5 OPTIMIZACIJA V DINAMIKI LINEARNO ELASTICˇNIH LINIJSKIH NOSILCEV 62 5.1 Uvod ............................................ 62 5.2 Optimizacijski problem ................................... 63 5.3 Projektne spremenljivke .................................. 64 5.4 Obcˇutljivostna analiza .................................... 64 5.5 Optimizacijski algoritem .................................. 66 5.6 Numericˇni primer ...................................... 66 5.7 Sklepi ............................................ 69 6 DINAMIKA ARMIRANOBETONSKIH LINIJSKIH KONSTRUKCIJ 70 6.1 Uvod ............................................ 70 6.2 Materialni modeli ...................................... 70 6.2.1 Jeklo ........................................ 71 6.2.2 Beton ........................................ 73 6.3 Integracija po prerezu .................................... 76 6.4 Numericˇna primera ..................................... 78 6.4.1 Pushover analiza (Prohitech) ............................ 78 X Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 6.4.2 Dinamična analiza montažne hale......................... 84 7 KNJIŽNICA KONČNIH ELEMENTOV 92 7.1.1 Kinematične teorije................................. 92 7.1.2 Stopnje interpolacije in integracije......................... 94 7.1.3 Materialni modeli.................................. 94 7.1.4 Časovni integratorji................................. 94 7.1.5 Vrsta analize in oznake elementov......................... 94 7.1.6 Uporabne funkcije................................. 94 7.2 Knjižnica končnih elementov................................ 95 7.2.1 Končni elementi za statično analizo........................ 95 7.2.2 Končni elementi za dinamično analizo....................... 96 8 ZAKLJUČKI 97 8.1 Prispevki k znanosti..................................... 99 VIRI 100 PRILOGE 109 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij XI Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. KAZALO SLIK Slika 2.1 Geometrija nosilca, krajevni vektorji in vektor pomikov, levo; linijska in točkovna obtežba, sredina; smeri rezultant napetosti ter prečni prerez, desno......... 10 Slika 3.1 Geometrija problema, levo; zasuk v vpetišču v odvisnosti od časa, desno...... 21 Slika 3.2 Odziv sistema in primerjava z rezultati Sima in Vu-Quoca (1986).......... 22 Slika 3.3 Računski model konzole z obtežbo, levo; poveš prostega konca (vb) v odvisnosti Slika 4.1 Geometrijski podatki o primeru, levo; odvisnost obtežbe od časa, desno....... 49 Slika 4.3 Shematski prikaz primerov ob vklopu in izklopu dušenja............... 54 Slika 4.4 Geometrija in obtežba Bauchaujevega nihala..................... 55 Slika 4.5 Bauchaujevo nihalo: vodoravni in navpični pomik v točki B............. 56 Slika 4.6 Bauchaujevo nihalo: povečava vodoravnih in navpičnih pomikov v točki B s slike Slika 4.7 Bauchaujevo nihalo: disipacija energije z različnimi stopnjami dušenja....... 58 Slika 4.8 Bauchaujevo nihalo: temna območja so dušena.................... 59 Slika 4.9 Bauchaujevo nihalo: notranje sile........................... 59 Slika 4.10 Bauchaujevo nihalo: frekvenčni sestav odziva notranjih sil.............. 60 Slika 5.1 Tipična uporaba projektnih spremenljivk: klasične, levo; oblike, sredina; vzbu- Slika 5.2 Računski model okvirja. C P je oznaka za kontrolno točko, DE pa za projektno telo. 66 Slika 5.3 Pomik podpor v odvisnosti od časa.......................... 67 Slika 5.4 Pomiki na sredini višine konstrukcije s črtkano črto in na vrhu konstrukcije s polno črto po prvi fazi optimizacije, levo in po končani optimizaciji, desno......... 68 Slika 5.5 Oblika konstrukcije po prvi fazi optimizacije, levo in po končani optimizaciji, desno. 68 XII Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 6.1 Shematski prikaz določitve točke To.......................... 71 Slika 6.2 Zgodovina deformacije v odvisnosti od 'časa’, levo. Pripadajoča napetost, desno. . 73 Slika 6.3 Zgodovina deformacije v odvisnosti od 'časa’, levo. Pripadajoča napetost, desno. . 75 Slika 6.4 Prerez, skrajno levo. Razdelitev prereza na lamele in njihovo številčenje, sredina. Ilustracija treh možnih integracij po prerezu, desno.................. 77 Slika 6.5 Fotografija celotne konstrukcija z označenim analiziranim delom.......... 78 Slika 6.6 Fotografija (levo) in načrt (desno) obravnavanega dela konstrukcije v vzdolžni smeri. 78 Slika 6.7 Načrt konstrukcije v prečni smeri, levo. Računski model, desno........... 79 Slika 6.8 Armaturni načrt konstrukcije z označenimi območji prečnih prerezov in sodelujočimi širinami................................. 80 Slika 6.9 Prečniprerezi.............................. 81 Slika 6.10 Rezultat monotone pushover analize. . ................ 82 Slika 6.11 Rezultat ciklicˇne pushover analize ................... 83 Slika 6.12 Geometrija obravnavane konstrukcije; s puščicami sta označeni mesti nanosa obtežbe. 84 Slika 6.13 Prečni prerez stebra z vzdolžno in strižno armaturo.......... 85 Slika 6.14 Računski model konstrukcije..................... 85 Slika 6.15 Normaliziran akcelerogram. . .................... 86 Slika 6.16 Odziv na obremenitev skalirano z 0.14 g. . ............. 88 Slika 6.17 Odziv na obremenitev skalirano z 0.35 g. . ............. 89 Slika 6.18 Odziv na obremenitev skalirano z 0.525 g ............... 90 Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij XIII Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. LIST OF FIGURES Figure 2.1 Geometry of the beam, position vectors and vector of displacements, left; distributed and point loads, middle; stress resultants and cross-section, right. . . . . 10 Figure 3.1 Geometry of the problem, left; rotation at pinned end vs time, right. . ....... 21 Figure 3.2 Response of the system and comparison to results of Simo and Vu-Quoc (1986). . 22 Figure 3.3 Cantilever beam model and load, left; time-history response for tip deflection vB, right .......................................... 26 Figure 4.1 Geometrical data, left; loads vs time, right ...................... 49 Figure 4.2 Free flight of spaghetti. . .............................. 50 Figure 4.3 Schematic representations of cases when damping is engaged and disengaged. . . 54 Figure 4.4 Geometry and loading of the Bauchau swinging pendulum. . ........... 55 Figure 4.5 Bauchau swinging pendulum: horizontal and vertical displacements at point B. . 56 Figure 4.6 Bauchau swinging pendulum: magnification of the horizontal and vertical displacements at point B from Fig. 4.5. . ....................... 57 Figure 4.7 Bauchau swinging pendulum: energy dissipation with different levels of damping. 58 Figure 4.8 Bauchau swinging pendulum: dark areas mark presence of damping ........ 59 Figure 4.9 Bauchau swinging pendulum: internal forces. . .................. 59 Figure 4.10 Bauchau swinging pendulum: frequency content of internal forces ......... 60 Figure 5.1 Typical applications of design variables: classical, left; shape, center; actuation, right .......................................... 64 Figure 5.2 The model of the frame. CP are control points and DE design elements. . . . . . 66 Figure 5.3 Movement of supports vs time. . .......................... 67 Figure 5.4 Displacements at mid-height and at the top of the structure are plotted with dashed and full line, respectively. The results of the first phase are on the left hand side and the final results on the right hand side of the figure. . ............. 68 Figure 5.5 The shape of the structure after the first phase of optimization is on the left hand side and the final shape on the right hand side of the figure. . ........... 68 XIV Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Figure 6.1 Schematical representation of locating point T0. . ................. 71 Figure 6.2 Strain history in the fibre, left. Related stress, right. . ............... 73 Figure 6.3 Strain history in the fibre, left. Corresponding stress, right. . ........... 75 Figure 6.4 Cross-section, left. Division and enumeration of the cross-section into layers, middle. Illustrations of some of the possible integration schemes, right. . ..... 77 Figure 6.5 Photography of the entire structure. The analysed section is framed. . ...... 78 Figure 6.6 Photography (left) and plan (right) of the analysed section of the structure in longitudinal direction ................................. 78 Figure 6.7 Schematics of the strucutre in transversal direction, left. Mathematical model of the structure, right. . ................................ 79 Figure 6.8 Plan of reinforcement with marked cross-sections and effective widths ....... 80 Figure 6.9 Cross-sections. . ................................... 81 Figure 6.10 The results of the monotone pushover analysis. . .................. 82 Figure 6.11 The results of the cyclic pushover analysis ...... ................ 83 Figure 6.12 The geometry of the structure; arrows depict locations and directions of the loads. 84 Figure 6.13 Cross-section of the beam with longitudinal and shear reinforcement. . . ..... 85 Figure 6.14 Numerical model of the structure. . ......................... 85 Figure 6.15 Normalized accelerogram. . ............................ 86 Figure 6.16 Response to loads scaled to 0.14 g .......................... 88 Figure 6.17 Response to loads scaled to 0.35 g .......................... 89 Figure 6.18 Response to loads scaled to 0.525 g. . ....................... 90 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij XV Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. KAZALO PREGLEDNIC Preglednica 3.1 Izboljsˇana formulacija. . ............................ 26 Preglednica 3.2 Osnovna verzija formulacije. . ........................ 27 Preglednica 4.1 Algoritmi iz Newmarkove druzˇine integratorjev. . .............. 31 Preglednica 4.2 Parametri algoritmov. . ............................ 43 Preglednica 4.3 Sˇpaget: rezultati racˇuna ............................. 50 Preglednica 6.1 Parametri materialnega modela za jeklo. . .................. 70 Preglednica 6.2 Parametri materialnega modela za beton. . .................. 73 Preglednica 7.1 Knjizˇnica koncˇnih elementov za staticˇno analizo. . .............. 94 Preglednica 7.2 Knjizˇnica koncˇnih elementov za dinamicˇno analizo. . ............ 95 XVI Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. LIST OF TABLES Table 3.1 Numerically improved formulation. . ........................ 26 Table 3.2 Basic version of the formulation. . .......................... 27 Table 4.1 Algorithms of the Newmark family of integrators. . ................. 31 Table 4.2 Algorithm parameters. . ............................... 43 Table 4.3 Spaghetti: results. . .................................. 50 Table 6.1 Parameters of the steel material model. . ...................... 70 Table 6.2 Parameters of the concrete material model. . .................... 73 Table 7.1 Finite element library for static analysis. . ...................... 94 Table 7.2 Finite element library for dynamic analysis ...................... 95 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij 1 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 1 UVOD 1.1 Predstavitev problema in pregled stanja na podrocˇju dinamike ravninskih okvirnih konstrukcij Linijski nosilci so osnovni gradniki matematičnih modelov gradbenih konstrukcij ter mehanizmov in manipulatorjev v strojni, letalski in vesoljski industriji. Temeljito razumevanje njihovega mehanskega obnašanja je zato pogoj uspešne analize. Mehaniko nosilcev lahko opišemo s štirimi neodvisnimi pojavi: z ravnotežjem, s kinematiko, s konstitucijo in s konsistenčnimi pogoji. V gradbenih konstrukcijah ima od teh pojavov na odziv nedvomno najbolj pomemben vpliv konstitucija oz. materialni zakoni. Parametri materialnih zakonov so določeni z eksperimenti in so kot taki podatek v analizi. Velik vpliv na odziv ima tudi kinematika. Kinematične teorije lahko razdelimo na tri skupine: na linearno teorijo, na točno teorijo in na teorije med linearno in točno. Za nekatere slednje se uporablja izraz ‘teorije drugega reda’. Točna kinematična teorija (tudi ‘geometrijsko točna’ in ‘Reissnerjeva’) ima izmed omenjenih najširše območje veljavnosti, ki pa ni neomejeno. Njene najbolj značilne lastnosti so upoštevanje obtežbe na deformirani obliki nosilca, povezanost upogibnega in osnega deformiranja in upoštevanje striga. Linearna teorija ima ravno nasprotne lastnosti - obtežbo upošteva na začetni legi in ne upošteva povezanosti upogibnega in osnega deformiranja. Zato je njeno območje veljavnosti najožje. Edina prednost linearne kinematične teorije je njena preprostost in linearna forma, ki vedno omogoča analitično reševanje, če je material linearno elastičen. Teoriji sta si v resnici bolj blizu, kot bi lahko ocenili iz diametralno nasprotnih lastnosti, saj linearno teorijo dobimo z linearizacijo točne teorije. Konsistenčni pogoji so pogoji, s katerimi zahtevamo enakost konstitucijskih in ravnotežnih notranjih sil. Pri večini formulacij teh pogojev v eksplicitni obliki ni, v določenih pa so in imajo vpliv na njihove lastnosti, toda več o njih in njihovem vplivu na lastnosti v nadaljevanju. Ravnotežje namenoma omenjamo zadnje, saj je za gradbenike ravnotežje osnovni pojem in vsakdanji kruh. Ravnotežje v statiki pomeni ničelnost rezultant sil in momentov. Če v dinamiki produkte pospeškov in mas nadomestimo z linijskimi obtežbami, si lahko ravnotežje v dinamičnem smislu predstavljamo na enak način. V disertaciji bomo pokazali, da tako dinamično ravnotežje ni edino možno in bomo predstavili alternativno, t.i. ‘energijsko dinamično ravnotežje’. Za formuliranje in reševanje problemov moramo te pojave združiti. To bomo naredili na različne načine, vedno pa bomo uporabili kinematično točno teorijo, ki je zato rdeča nit disertacije, skupna vsem problemom, ki jih bomo obravnavali. Najprej se bomo dotaknili fundamentalnega vprašanja izbire neznanih funkcij, zatem se bomo ukvarjali s problemom časovne integracije in takrat se bomo srečali z energijskim ravnotežjem. V sklopu obravnavanja časovne stabilnosti se bomo ukvarjali tudi z vprašanjem 2 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. najbolj učinkovitega numeričnega dušenja. Posebno poglavje bo namenjeno optimizaciji v dinamiki geometrijsko točnih elastičnih nosilcev. Na koncu bomo točno kinematično teorijo uporabili še pri modeliranju armiranobetonskih konstrukcij v dinamiki. V splošnem lahko med formulacijami za dinamiko nosilcev zasledimo dva fundamentalno različna opisa gibanja dinamičnega sistema. Prvi pristop je t.i. pristop premikajočega se ogrodja oz. ‘floating frame’, kjer so kinematične spremenljivke sistema določene glede na toge skelete konstrukcij, ki opišejo togi del gibanja in se premikajo skupaj s podajnim okvirjem. Prednost tega pristopa je, da so pomiki in zasuki togih okvirjev lahko poljubno veliki, deformacijske količine pa so običajno majhne. Zato ima izraz za potencialno energijo nosilca preprosto obliko. Slabost pristopa je, da so, kljub majhnim deformacijam glede na togi skelet, enačbe izrazito nelinearne in da so vztrajnostni členi v posameznih smereh močno povezani - glej Crisfield in Shi (1994), Hsiao in Jang (1991) in Kane, Likins in Levinson (1983). Včasih ta pristop imenujemo korotacijski oz. ‘co-rotational’. Drugi pristop, ki se je bolj uveljavil in ga privzemamo tudi v tem delu, se imenuje vztrajnostni oz. ‘inertial frame approach’. Pri njem se kinematične spremenljivke določi glede na fiksni koordinatni sistem in zato so izrazi za vztrajnostne člene preprosti. Nelinearnost enačb in problema pa se pojavi pri zapisu togostnega dela enačb oz. deformacijske energije nosilca. Različnih formulacij za opis dinamike ravninskih elastičnih nosilcev s točnimi ali nelinearnimi geometrijskimi zvezami je več. Formulacije, ki interpolirajo pomike in zasuke, so predstavili Crisfield in Shi (1994), Crisfield, Galvanetto in Jelenič (1997), Hsiao in Jang (1991), Ibrahimbegovič in Mamouri (1999), Sansour, Sansour in Wriggers (1996), Simo in Vu-Quoc (1986), Simo, Tarnow in Doblare (1995), formulacije z interpolacijo absolutnih koordinat v vozliščih pa von Dombrowski (2003), Dufva, Sopanen in Mikkola (2005), Omar in Shabana (2001), Shabana (1997) in Shabana (1998). Formulacije z interpolacijo absolutnih vozliščnih koordinat se v izvirniku imenujejo ‘ANCF’ oz. ‘absolute nodal coordinate formulation’. Jelenič (1990) in Vratanar (1998) sta predstavila formulacije, v katerih se interpolirajo izključno zasuki. Prva tema, s katero se v disertaciji spoprimemo, je vprašanje izbire neznanih količin. Ker je možnosti veliko, se osredotočimo na možnost upoštevanja deformacijskih količin kot osnovnih neznank problema. Ta možnost je zanimiva zato, ker ima tako formulirani problem posebno obliko enačb, ki prinese nekatere ugodne lastnosti. Posebnost oblike enačb je v tem, da v njih eksplicitno nastopajo konsistenčni pogoji. Pomembnejše lastnosti formulacij s tako obliko enačb so odpornost na blokiranje, popolna konsistentnost količin in velika natančnost. Podobno obliko enačb imajo tudi formulacije, ki za neznane funkcije upoštevajo le zasuke. Taka je Jeleničeva (1990) formulacija za dinamično analizo hiperelastičnih nosilcev in Vratanarjeva (1998) dinamična analiza idealno elastičnih nosilcev z drugo časovno integracijsko shemo in s kolokacijsko metodo. Izključno na zasukih temelji tudi statična analiza ravninskih okvirjev Sajeta in sodel. (1990, 1991, 1994, 1997) in Jeleniča in Sajeta (1994, 1995). Zelo podobno obliko enačb imajo tudi formulacije za analizo statičnih problemov, ki temeljijo na deformacijah. Prvi je take elemente izpeljal Planine (1998), ko je analiziral problem stabilnosti elastičnih in elastoplastičnih ravninskih nosilcev. Končne elemente z interpolacijo deformacijskih količin je uporabil tudi Bratina s sodelavci (2003, 2004) pri požarni analizi ravninskih armiranobetonskih konstrukcij in Čas (2004) pri analizi kom-pozitnih nosilcev z zdrsom med sloji. Zupan (2001, 2003) in Zupan in Saje (2003) sta jih uporabila pri modeliranju prostorskih nosilcev, Hozjan, Turk in Srpčič (2007) pa pri analizi jeklenih konstrukcij med požarom. Schnabl (2007) je z njimi modeliral kompozitne nosilce med požarom z upoštevanjem oglenenja lesa. Krauberger s sodel. (2007) je podobne enačbe s konsistenčnimi pogoji uporabila kot izhodišče za določitev uklonskih sil armiranobetonskih stebrov. Formulacija, ki jo zgradimo okoli interpolacije deformacijskih količin, uporablja vztrajnostni pristop in, Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 3 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. podobno kot deformacijske formulacije za statične primere, eksplicitno zadošča konsistenčnim enačbam. Formulacija je relativno zapletena, dobljeni sistem enačb pa ni samo diferencialen, temveč diferencialno-algebrajski. S tem so povezane težave pri časovni integraciji. V drugem delu disertacije se odmaknemo od elementov z interpolacijo deformacij, saj za osnovne neznane količine privzamemo kinematične spremenljivke in se spoprimemo s temo časovne integracije. Ključni problem časovne integracije v nelinearni dinamiki je stabilnost. Kuhi in Crisfield (1999) celo trdita, da je stabilnost bolj pomembna od natančnosti. Kljub temu pa se ukvarjamo z metodami z vsaj kvadratično konvergenco, saj nižje stopnje konvergence niso sprejemljive. Matematični model gradbenih konstrukcij v geometrijsko točni teoriji nosilcev je ponavadi predstavljen s togim sistemom diferencialnih enačb. Togost enačb je posledica za rede velikosti različnih togosti posameznih delov konstrukcije ali pa bistveno različnih osnih in upogibnih togosti nosilca. Togost lahko povzročijo tudi nenadne spremembe v sistemu, ki so lahko posledica hitre spremembe obtežbe. Dodaten izvor togosti je uporaba metode Lagrangevih množiteljev za upoštevanje kinematičnih vezi, saj postane sistem enačb zaradi algebrajskih vezi diferencialno-algebrajski. Algebrajskim enačbam pripadajo neskončne lastne vrednosti (Gams, Planine in Saje (2007c)), zato so inherentno togi. K integraciji takih problemov so prispevali Bauchau, Damilano in Theron (1995), Bauchau in Theron (1996a, b), Bauchau, Bottasso in Trainelli (2003), Betsch in Steinmann (2003), Bottasso, Bauchau in Choi (2002), Bottasso in Borri (1998), Goicolea in Orden (2002) in Ibrahimbegovič s sodel. (2000). Zelo pogosto so togi tudi podajni mehanizmi (angl. ‘flexible multibody systems’), kjer so podajni nosilci med sabo povezani z različnimi vrstami vezi oz. členkov (Bauchau in Theron (1996a)). Tipična značilnost togih sistemov je pojav visokih frekvenc umetnega izvora. Sposobnost kontroliranja teh frekvenc je zagotovilo kakovostne in časovno stabilne metode. Integriranje togih enačb je zelo težavno, saj je problem časovne nestabilnosti pri njih najbolj izrazit. Klasične integracijske sheme, ki so bile razvite za linearne sisteme, kot je npr. Newmarkova družina in-tegratorjev in a integratorji, imajo pri reševanju nelinearnih togih diferencialnih enačb velike probleme, kar so nazorno pokazali Crisfield, Galvanetto in Jelenič (1997), Géradin in Cardona (2001), Ibrahimbegovič in Mamouri (1999), Jelenič in Crisfield (2001), Stander in Stein (1995) in Erlicher, Bonaventura in Bursi (2002). Stabilnost metode lahko izboljšamo na več načinov. En izmed pristopov je, da v formulacijo vpeljemo numerično dušenje, ki disipira energijo višjih oblik nihanja. Osnovni namen numeričnega dušenja je iz odziva izločiti višje frekvence, toda z dušenjem se pogosto izboljša tudi stabilnost reševanja. Zelo razširjena je disipativna shema Hilberja, Hughesa in Taylorja (1977) oz. metoda HHT in njena posplo-šitev, HHT-a; metoda Chunga in Hulberta (1993). Ti shemi sta zelo učinkoviti pri reševanju linearnih sistemov, saj sta tedaj brezpogojno stabilni, pri reševanju nelinearnih sistemov pa po poročanju Géradina in Cardone (2001) pogosto zadovoljivo delujeta, vedno pa ne. Bauchau, Damilano in Theron (1995) in Crisfield, Galvanetto in Jelenič (1997) so na praktičnih primerih pokazali, da lahko pride do rasti mehanske energije tudi v primeru uporabe teh shem (z dušenjem), kar neizogibno vodi v nestabilnost. Numerično dušenje so na nekoliko drugačen način implementirali Kuhi in Crisfield (1999) na prostorskem paličju in Kuhi in Ramm (1999) na lupinah in ta pristop poimenovali posplošeni energijsko količinski pristop oz. ‘generalized energy momentum’. Njihov pristop je kombinacija pristopa posplošenih a metod in pristopa z ohranjanjem energije in količin Sima, Tarnowa in Doblareja (1995). Nezvezni oz. ‘time-diseontinous’ Galerkinov pristop za dušenje v svojih delih uporabijo Bauchau in Theron (1996a, b), Bottasso in Borri (1998) in Bottasso, Bauchau in Choi (2002). Armero in Romero 4 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. (2001a, b) vpeljeta disipacijo energije z diferenčno metodo. Goicolea in Orden (2002) za dušenje pri reševanju nelinearnega Hamiltonskega sistema uporabita koncept diskretnih odvodov in kazensko metodo (angl. ‘penalty method’) za zadoščanje vezi. Poročata o izredni učinkovitosti in robustnosti svojega pristopa. Romero in Armero (2002) vpeljeta dušenje v formulacijo prek spremenjenih konstitucijskih enačb. Enak pristop uporabita Ibrahimbegovič in Mamouri (2002) na prostorskih nosilcih. Za reševanje togih enačb gibanja v splošnem Hamiltonskem kontekstu s časovno zvezno (angl. ‘tirne-continous’) Galerkinovo metodo so svoje pristope predstavili tudi Betsch in Steinmann (2000) in Bui (2003). Pristop Armera in Petöcza (1999) temelji na modificiranem računu deformacij. Vse omenjene sheme so implicitne in dokazano disipirajo energijo nelinearnih sistemov, kar je po Belytschku in Schoe-berlu (1975) formalen dokaz stabilnosti sheme. Bottasso, Borri in Trainelli (2001) so pokazali, da lahko disipacijo za nelinearne sisteme vpeljejo tudi v eksplicitne integratorje družine Runge-Kutta. V pričujoči disertaciji je razvita nova shema z dušenjem, ki deluje na deformacijah - podobno kot postopek Armera in Petöcza (1999), le daje nadgrajena z algoritmom, ki po potrebi vklaplja in izklaplja dušenje. Slednje tako deluje le na mestih in pri časih, ko je to potrebno. Analiza frekvenčnega odziva pokaže, da s tem načinom disipiramo manj osnovnih frekvenc. Drug pristop za izboljšanje stabilnosti je vpeljava dodatnih vezi po metodi Lagrangevih množiteljev, s katerimi se eksplicitno zagotavlja ohranjanje energije in količin gibanja. Po poročanju Kuhla in Ramma (1996) rezultati kažejo, da take metode popolnoma ohranjajo zahtevane količine, vendar vseeno niso stabilne. Tretji pristop za izboljšanje stabilnosti reševanja so t.i. konzervativne sheme, ki algoritemsko ohranjajo energijo in količine gibanja. Algoritemsko ohranjanje pomeni, da rešitev enačb problema avtomatično zadosti tudi pogojem ohranjanja količin. Take sheme so brezpogojno stabilne (Hughes (1992)), vendar se lahko v odzivu pojavijo nerealistično visoke frekvence, kar lahko po poročanju Bauchauja, Damilana in Therona (1995) vpliva na hitrost konvergence in uspešnost Newton-Raphsonovega postopka reševanja sistema enačb. Pionirska dela na področju energijsko konzervativnih shem so predstavili Simo, Tarnow in Wong (1992), Simo, Tarnow in Doblare (1995) in Stander in Stein (1996). Ibrahimbegovič s sodel. (1999, 2002) je razvil podobno shemo kot Simo, Tarnow in Doblare (1995), vendar z drugačno interpolacijo zasukov. Tej razliki bomo namenili nekaj besed tudi v nadaljevanju. Energijske sheme uporabljajo t.i. ‘midpoint’ pravila, ki enačbe izvrednotijo pri sredinskem času, vendar jih lahko v primeru nelinearnih zvez interpretiramo na več načinov, kot je nazorno pokazal Graham (2001). Jelenič in Crisfield (2001) sta implementirala energijsko shemo za prostorske nosilce, v kateri členke modelirata s ti. ‘master-slave’ pristopom. Crisfield in Shi (1996) sta izpeljala energijski integrator za korotacijsko formulacijo ravninskegapaličja. Sansour, Sansour in Wriggers (1996) so uporabili midpoint formulacijo za ravninske nosilce z uporabo Cosseratove teorije nosilcev. Znaten del disertacije se ukvarja z analizo energijskih shem za časovno integracijo in v sklopu teh analiz razvijemo novo energijsko konzervativno shemo, za katero bomo pokazali tudi ohranjanje gibalne in vrtilne količine. Posebnost formulacije je, da se upošteva kinematične enačbe v časovno odvajani obliki in da se deformacije računa na inkrementni način. Tip formulacije in interpolacije sledi delu Sima in Vu-Quoca (1986), kjer se privzame vztrajnostni pristop in interpolira vse kinematične spremenljivke, končne enačbe pa so diferencialne. Prvič so idejo o inkrementnem računu deformacij in upoštevanju kinematičnih enačb v šibki obliki predstavili Sansour, Wriggers in Sansour (1997) za lupine. Kasneje je Sansour s sodel. (2002) ta pristop apliciral na klasično in kaotično gibanje lupin in dinamično stabilnost lupin. Sansour, Wriggers in San- Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 5 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. sour (2004) so predstavili posplosˇeno teorijo s takim pristopom, ki nacˇelno velja za poljubne konstrukcije in 3D kontinuum. V disertaciji se bezˇno dotaknemo tudi optimizacije. V splosˇnem predstavlja optimizacija mozˇnost zadostiti kriterijem naloge ceneje, hitreje, boljsˇe in bolj ucˇinkovito. Potencial mozˇnosti in razpon optimizacije je za podajne nosilce izredno velik – od lazˇjih, cenejsˇih in bolj odpornih konstrukcij, do manipulatorjev in robotskih rok, ki izkorisˇcˇajo podajnost mehanizma, da nalogo opravijo hitreje in z manjsˇo porabo energije. Zacˇetki poglobljenih analiz optimizacije linearnih dinamicˇnih sistemov segajo v sedemdeseta leta prej-sˇnjega stoletja, kot porocˇajo Kang, Park in Aurora (2006), ki v svojem cˇlanku temeljito opisˇejo razvoj gradientnih metod optimizacije na podrocˇju dinamike linearnih sistemov. Le stezˇka pa najdemo dela, ki bi se ukvarjala z optimizacijo kinematicˇno ali pa materialno nelinearnih problemov. Glavni razlog je relativno zahteven racˇun odziva konstrukcije, ki je hkrati tudi motiv za uporabo gradientnih optimizacijskih metod, kot navajajo Cardoso in Arora (1992) in Arora in Dutta (1997). Stupkiewicz (2001) v svojem prispevku pokazˇe, kako dobiti priblizˇne vrednosti za obcˇutljivostno analizo v nelinearnih dinamicˇnih problemih. Kulkarni in Noor (1995) opisˇeta obcˇutljivostno analizo ravninskih viskoplasticˇnih konstrukcij glede na materialne parametre. Cho in Choi (2000a, b) v svojih delih predstavita postopek obcˇutljivostne analize in gradientne optimizacije v primeru optimizacije prerezov in oblike konstrukcije. Optimizacijo ciklicˇno obremenjenih elastoplasticˇnih konstrukcij opisujejo tudi Sousa, Cardoso in Valido (1997). Pedersen (2004) v svojem prispevku opisˇe optimizacijo topologije elastoplasticˇnih ravninskih nosilcev. V sklopu optimizacije podajnih nosilcev z novo razvitim integratorjem za podajne elasticˇne nosilce s sˇibkimi kinematicˇnimi vezmi in inkrementnim nacˇinom racˇuna deformacij izpeljemo in implementiramo obcˇutljivostno analizo, ki je ena izmed kljucˇnih komponent optimizacije. Za resˇevanje optimizacijskega problema uporabimo optimizacijski algoritem, ki so ga razvili Kegl, Butinar in Kegl (2002) in je implementiran v racˇunalnisˇkem programu ‘iGO’ (M. Kegl, 2007). Tukaj predstavljamo le posebnosti, s katerimi se srecˇamo ob optimiranju dinamicˇnih problemov in posebnosti integratorja in elementa, ki jih moramo pri tem uposˇtevati. Ko kinematicˇno tocˇno teorijo ravninskih nosilcev nadgradimo z materialnim modelom betona po Popo-vicsu (1973), ki razbremenjevanje v tlacˇni coni uposˇteva po Karsanu in Jirsi (1963) in z materialnim modelom jekla po Menegottu in Pintu (1973), se nam odpre mozˇnost modeliranja armiranobetonskih konstrukcij in testiranja cˇasovno integracijskih shem na materialno nelinearnih problemih. Dinamicˇne analize armiranobetonskih konstrukcij z uposˇtevanjem prej omenjenih materialnih modelov so dolga leta temeljile na geometrijsko linearni teoriji, saj so bili znanstveniki bolj osredotocˇeni na modeliranje razlicˇnih materialnih pojavov. Analize, ki bi hkrati uposˇtevale geometrijsko in materialno nelinearnost, so zato redke. Najpogosteje uporabljeni nacˇini za uposˇtevanje geometrijske nelinearnosti v teoriji konstrukcij so t.i. P-? metode ali pa t.i. teorije drugega reda, ki niso veljavne v primeru velikih deformacij. Razvoj materialno nelinearnih modelov v teoriji konstrukcij se je zacˇel s koncˇnimi elementi s koncentrirano plasticˇnostjo (plasticˇni cˇlenki), kot v svojih preglednih delih navajajo Taucer, Spacone in Filippou (1991) ter Spacone, Taucer in Filippou (1996). Taki elementi se v praksi sˇe danes veliko uporabljajo, predvsem na racˇun svoje zanesljivosti in relativne preprostosti, vendar pa so sposobni opisati le nekaj osnovnih pojavov v armiranobetonskih konstrukcijah. Alternativa tem elementov so elementi s t.i. porazdeljeno plasticˇnostjo oziroma elementi z integracijo po t.i. vlaknih (angl. ‘fiber elements’). Vlakna predstavljajo materialne niti, iz katerih je sestavljen precˇni prerez. Vsaka nit ima svoj neodvisen ma- 6 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. terialni zakon in z integracijo napetosti po materialnih nitih (vlaknih) dobimo konstitucijske notranje sile v prečnem prerezu. Primerjavo nekaterih elementov s točkovno in porazdeljeno plastičnostjo lahko najdemo v članku Coma, De Stefana in Ramascota (2003), kjer je ugotovljeno, da dajo elementi z integracijo po vlaknih bistveno boljše rezultate kot elementi s točkovno plastičnostjo, saj so v nasprotju z elementi s točkovno plastičnostjo inherentno sposobni upoštevati interakcijo osne sile in upogibnega momenta. Elementi z integracijo po vlaknih, osnovani na klasični metodi končnih elementov, kjer se interpolirajo kinematične količine, imajo v statičnih problemih po poročanju Taucerja, Spaconeja in Fillipouja (1991) težave pri modeliranju mehčanja materiala. Mehčanje je upoštevanje nosilnosti po preseženi maksimalni nosilnosti in je zelo značilno za armiranobetonske konstrukcije. Takrat postanejo elementi numerično zelo občutljivi, saj raziskovalci poročajo o problemih z nestabilnostjo reševanja, predvsem v povezavi z neustrezno interpolacijo, ki ne uspe zajeti koncentracije ukrivljenosti v plastičnih členkih. Za pristop, ki interpolira kinematične količine, se je v krogih, ki se ukvarjajo z armiranobetonskimi konstrukcijami, uveljavil izraz ‘togostni’. Alternativa temu pristopu je ‘podajnostni’ pristop, ki temelji na interpolaciji sil in se je zelo uveljavil, predvsem po objavi članka Spaconeja, Taucerja in Filippouja (1996). V njem so namreč rešene ključne ovire za vgradnjo podajnostnih elementov v standardne programe za analizo po metodi končnih elementov. Podajnostno formulacijo sta Neuenhofer in Filippou (1998) razširila na geometrijsko nelinearnost z upoštevanjem interakcije osnega in upogibnega deformiranja, vendar pa morajo deformacije v tej formulaciji ostati majhne, čeprav velikost pomikov ni omejena (teorija drugega reda). Njuna formulacija je osnovana na interpolaciji ukrivljenosti po dolžini nosilca. Petrangeli in Ciampi (1997) v svojem članku ugotavljata prednosti podajnostnega pristopa v primerjavi s togostnim pristopom in mešanimi pristopi, vendar le za geometrijsko linearno teorijo. Ademar in White (2005) opravita podobno analizo, toda uporabita teorijo drugega reda, katere veljavnost omejita na srednje velike rotacije in majhne deformacije, ali pa P-A metodo. Njun zaključek je, da je mešana formulacija najboljša za materialno in geometrijsko nelinearne probleme. Sivaselvan in Reinhorn (2002) za dinamično analizo armiranobetonskih konstrukcij uporabita podajnostni pristop, toda namesto običajnih končnih elementov sestavljata makro elemente, ki potekajo med fizičnimi vozlišči konstrukcije. Uporabljene geometrijske zveze so točne, toda uporabljena je korotacij-ska formulacija. S splošnim programom za reševanje diferencialno algebrajskih enačb lahko analizirata konstrukcije do porušitve. Članek Hjelmstada in Tacirogluja (2003) poskuša sistematično primerjati geometrijsko (teorija drugega reda) in materialno nelinearne formulacije, kjer je posebna pozornost namenjena računski zahtevnosti. Primerjane so tri variacijsko konsistentne formulacije, ki jih splošno lahko označimo kot eno, dvo in tripoljske oz. kot klasični pristop (interpolacije in variacije pomikov), Hellinger-Reissnerjev pristop (interpolacije in variacije pomikov in napetosti) ter mešani pristop, ki dopušča interpolacijo in variacije tako pomikov kot deformacij in napetosti. Zaključek v članku je, da čistega zmagovalca med formulacijami ni, le določene formulacije so bolj primerne za določene vrste problemov. Lokalizacijo deformacij v statični analizi armiranobetonskih okvirjev so učinkovito opisali Bratina, Saje in Planine (2004). Njihova formulacija temelji na končnih elementih z interpolacijo deformacij, s kratkimi končnimi elementi na mestih pričakovanih koncentracij deformacij. Posebnost kratkih elementov je, da se lahko mehčajo. Z večanjem obtežbe preko maksimalne nosilnosti se v kratkih elementih povečajo deformacije in zmanjša odpornost. Ostali elementi ne morejo preiti v mehčanje in se lahko le razbremenijo. Lokalizacijo s končnimi elementi z interpolacijo kinematičnih količin in Timošenkovo teorijo nosilcev sta Ehrlich in Armero (2005) modelirala z nezveznostjo pomikov v območju plastičnih členkov. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 7 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 1.2 Vsebina dela Delo je, poleg uvodnega poglavja, sestavljeno še iz sedmih poglavij, ki obravnavajo različne teme iz dinamike ravninskih nosilcev. V drugem poglavju ‘Opis problema’ zberemo vse izreke, definicije in količine, potrebne za opis dinamike elastičnih ravninskih nosilcev in predstavimo osnovno obliko enačb problema. Enačbe in neznanke analiziramo ter predstavimo možnosti za izbiro osnovnih neznanih količin problema. Izbira osnovnih neznanih funkcij je namreč ključno in izhodiščno vprašanje, ponuja pa izredno velik nabor možnosti. V tretjem poglavju ‘Dinamika linearno elastičnih linijskih končnih elementov z interpolacijo deformacij-skih količin’ izpeljemo in rešimo dinamični problem ravninskih elastičnih nosilcev, v katerem so edine neznane funkcije deformacije. Izkaže se, da so izpeljani elementi neobčutljivi na blokiranje ter zelo natančni. Žal pa imajo v osnovni obliki formulacije dve izraziti pomanjkljivosti: (i) v računskem smislu so izredno neučinkoviti in (ii) v kombinaciji z Newmarkovo časovno integracijsko shemo so nestabilni. Zato v nadaljevanju tega poglavja z modifikacijo metode drastično izboljšamo numerično učinkovitost, z izbiro druge časovno integracijske sheme pa občutno izboljšamo problem časovne stabilnosti, vendar ga žal ne odpravimo popolnoma. Osnovno in izboljšano formulacijo preverimo in primerjamo na praktičnih računskih primerih. Rezultati računa potrdijo teoretične ugotovitve. V obsežnem četrtem poglavju ‘Dinamika linearno elastičnih linijskih končnih elementov z interpolacijo kinematičnih količin’ rešujemo problem tako, da obravnavamo pomike in zasuk kot osnovne neznanke problema. Tak pristop je sicer najbolj razširjen v komercialnih programih in tudi med raziskovalci. V prvem delu poglavja predstavimo klasični pristop ter njegovo aplikacijo na nelinearne probleme, zatem pa pozornost namenimo t.i. energijskemu pristopu. Posebnost tega pristopa je, da integratorjev ne izpeljujemo iz osnovnih enačb, temveč iz izreka o ohranitvi mehanske energije. Veljavnost tega izreka je namreč potrebni in zadostni pogoj za brezpogojno stabilnost. Iz energijskega izhodišča je možno izpeljati več različnih integratorjev in to tudi storimo, dodatno pa izpeljemo še nov integrator za ravninske nosilce, ki prav tako zadošča pogoju stabilnosti. Ohranjevalne lastnosti integratorjev analitično dokažemo, zatem pa preverimo tudi z računskim primerom. Potrdijo se vse teoretične ugotovitve: ohranjanje energije in količin je neodvisno od velikosti časovnega koraka in stopnje interpolacije. V zaključnem delu tega poglavja na integratorje apliciramo še dušenje, ki deluje na deformacijah. Posebnosti in novosti pristopa sta dve: (i) dušenje je odvisno od trigonometrične funkcije arctan in ne od linearne funkcije, kot običajno, ter (ii) v dušenje vključimo algoritem, ki z relativno malim številom dodatnih parametrov dobro ugotovi, ali je dušenje potrebno ali ne. S tem, ko dušenje apliciramo le na mestih, kjer je to potrebno in le pri časih, ko je to potrebno, dosežemo, da se udušijo le najvišje frekvence, odziv v osnovnih nihajnih oblikah pa ostane relativno nedotaknjen. Sledi peto poglavje ‘Optimizacija v dinamiki linearno elastičnih linijskih nosilcev’, kjer predstavljamo osnovne koncepte optimizacije prerezov, oblike konstrukcije in premikov podpor. Te koncepte apliciramo na novo razvito formulacijo s šibkim integratorjem in povprečno rotacijsko matriko, za katero predstavimo tudi občutljivostno analizo. Na zanimivem primeru prikažemo možnost optimizacije oblike, v prilogi pa optimizacijo premikanja podpor. Optimizacija premikanja podpor spada v skupino kontrolnih problemov. V šestem poglavju ‘Dinamika armiranobetonskih linijskih konstrukcij’ najprej predstavimo materialne modele za dinamično (ciklično) obremenjevanje betona in jekla. Te modele zatem vgradimo v računalniški program v kombinaciji z različnimi časovno integracijskimi shemami. Materialni model validiramo in verificiramo s primerjavo eksperimentalnega ‘pushover’ testa in numerične simulacije armiranobetonske konstrukcije, zatem pa še z dinamičnim testom preproste konstrukcije. 8 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. V sedmem poglavju opišemo knjižnico končnih elementov, ki so osnovani na interpolaciji kinematičnih količin. Elementi so razdeljeni v dve knjižnici: eno za elemente za statično analizo ‘2DBSta’ in drugo za elemente za dinamično analizo ‘2DBDyn’. V knjižnici ‘2DBSta’ so zbrani končni elementi za statično analizo elastičnih, jeklenih in armiranobetonskih konstrukcij in imajo vgrajene tri možne kinematične teorije - linearno, teorijo drugega reda in točno Reissnerjevo teorijo. V knjižnici so na voljo elementi s poljubnim materialnim modelom v kombinaciji s štirimi različnimi stopnjami interpolacije, od linearne do četrtega reda, in poljubno od treh kinematičnih teorij. Skupno je v knjižnici 36 končnih elementov. V knjižnici ‘2DBDyn’ so elementi za dinamično analizo, vendar le v kombinaciji z linearno in Reissnerjevo kinematično teorijo. Na voljo so v kombinaciji s tremi materialnimi modeli in redi interpolacij od linearnega do četrtega. Elementi so izdelani za sedem različnih integratorjev: Newmarkov integrator (Newmark (1959)), integrator Huberta, Hughesa in Taylorja (1977), integrator Standerja in Steina (1996), integrator Sima, Tarnowa in Doblareja (1995), integrator Ibrahimbegovića in Mamourija (1999) ter dva nova integratorja, ki sta razvita v tej disertaciji. Imenujemo ju šibki integrator s povprečno rotacijsko matriko in šibki z midpoint rotacijsko matriko. Končni elementi so razviti v programskem okolju AceGen (Korelc (2006)). To okolje omogoča avtomatsko generacijo programske kode za končne elemente in simbolno odvajanje. Vsi elementi so pripravljeni za delovanje v okolju AceFEM z naborom spremljajočih funkcij za analizo rezultatov in vmesnih rezultatov. Oba programa sta sicer komercialna, toda marsikateri primer je možno izračunati z demo verzijo, ki omejuje le število končnih elementov in računskih korakov v analizi. Demo verziji sta dostopni na internetnem naslovu www.fgg.uni-lj.si/Symech. Knjižnici končnih elementov in spremljajoča dokumentacija so javno dostopni na internetnem naslovu www.km.fgg.uni-lj.si/Matija/index.htm. Dostop do njih je omočen tudi skozi izmenjevalnik končnih elementov AceShare znotraj programskega sistema AceFEM. V osmem, sklepnem poglavju so povzeti zaključki disertacije in opisani prispevki k znanosti. Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij 9 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 2 OPIS PROBLEMA 2.1 Uvod Problemi v gradbeništvu so definirani z geometrijo konstrukcije, gradbenim materialom in obtežbo. Vsaka od teh kategorij je ponavadi zelo zapletena, tako da tudi navidezno preproste konstrukcije kmalu postanejo matematično tako zapletene, da analitične rešitve ne obstajajo. Ko problem obravnavamo dinamično, je primerov, za katere obstajajo analitične rešitve, še manj, saj so znane le za kombinacijo geometrijske linearnosti, linearno elastičnega materiala, matematično regularnih oblik (npr. nosilec, cilinder, ipd.) in matematično preprostih oblik obtežb, kot so konstantna, linearna ali periodična po času (npr. Paz in Leigh (2004)). Ob očitnem pomanjkanju analitičnih rešitev se gradbeniki pri projektiranju gradbenih konstrukcij praviloma poslužujemo v računalniške programe vgrajenih numeričnih metod. Besedilo obravnava dinamiko elastičnih, jeklenih in armiranobetonskih linijskih konstrukcij, ki jih lahko modeliramo z v začetni legi ravnimi, ravninskimi nosilci z neomejenimi pomiki in zasuki. Da lahko problem matematično opišemo, najprej vpeljemo osnovne količine, s katerimi je problem opisan, zatem naštejemo osnovne izreke in definicije in na koncu zapišemo vodilne enačbe problema. 2.2 Osnovne količine Osnovne količine, s katerimi opišemo geometrijo nosilca, so prikazane na sliki 2.1. X in Y sta osi globalnega koordinatnega sistema, y in z pa osi lokalnega koordinatnega sistema prereza. Kinematične spremenljivke. S tem terminom označujemo pomik u(s, t) in zasuk (p(s, t) materialne točke težiščne osi nosilca, določene z ločno dolžino s točke v nedeformirani legi. Obe količini sta funkciji ločne dolžine s in časa t. Vektor pomikov u(s, t) ima komponenti v globalnih smereh X in Y, ki ju označujemo z u(s, t) in v(s, t). Globalni koordinatni sistem je označen na sliki 2.1. Zasuk (p(s, t) predstavlja zasuk prečnega prereza nosilca. Kadar strižne deformacije niso prisotne, predstavlja zasuk ip(s, t) zasuk težiščne osi (Vratanar in Saje (1999)). V nadaljevanju besedila zapis funkcijske odvisnosti od časa in kraja pogosto opuščamo. Deformacijske spremenljivke. Definicije deformacij je podal Reissner (1972). Podane so glede na težiščno os in jih označujemo z s(s, t) in n(s, t), kjer sta komponenti vektorja e(s, t) = [e(s, t),7(s, t)] (2.1) osna deformacija e(s,t) in strižna deformacija j(s,t). Osna deformacija deluje v smeri normale na prerez. Strižna deformacija 7 predstavlja zasuk prečnega prereza glede na normalo na deformirano težiščno os nosilca. 10 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 2.1: Geometrija nosilca, krajevni vektorji in vektor pomikov, levo; linijska in tocˇkovna obtezˇba, sredina; smeri rezultant napetosti ter precˇni prerez, desno. Figure 2.1: Geometry of the beam, position vectors and vector of displacements, left; distributed and point loads, middle; stress resultants and cross-section, right. Zaradi ravninske narave problema so deformacije simetrične glede na lokalno os z in neodvisne od lokacije v smeri lokalne osi y (slika 2.1). Konstitucijske spremenljivke. Konstitucijske spremenljivke so rezultante napetosti po prerezu. Označujemo jih z Ne (s,t) in Me (5, t). V vektorju Ne (5, t) = [iVc(s, t), Qc{s, t)]T se ‘skrivata’ konstitu-cijska osna in prečna sila, Mç (s, t) je konstitucijski upogibni moment. Računamo jih s konstitucijskimi enačbami, ki so zapisane v nadaljevanju teksta. Ravnotežne spremenljivke. Ravnotežne spremenljivke oz. ravnotežne sile v prečnem prerezu so N(s,L) in M(s,L). Vektor N(s,L) = [7V(s,L), Q(s,t)] je sestavljen iz ravnotežne osne in ravnotežne prečne sile. Ravnotežne spremenljivke dobimo tako, da konstrukcijo na izbranem mestu prerežemo in na podlagi dinamičnih ravnotežnih enačb (od tod njihovo ime) določimo neznane sile N(s, t) in M(s, t) v prečnem prerezu. Osnovni opis lege in prereza. Lokacija poljubne materialne točke na težiščni osi v začetnem stanju ravnega, a za kot y?0 zasukanega nosilca, je določena s krajevnim vektorjem tq: ro(s) = X(s) = Y (s) Xq + S COS (Pq Yq + s sin (fQ (2.2) kjer je s ločna dolžina nedeformiranega nosilca, Xq in 1q pa sta koordinati pri s = 0 oz. na začetku nosilca. Ta materialna točka na osi deformiranega nosilca je določena s krajevnim vektorjem r(s,t): r(s, t) = rQ + 11(5, t) = X + u(s,t) Y + v(s, t) (2.3) Z A označujemo ploščino prečnega prereza, z As ploščino strižnega prečnega prereza in z J vztrajnostni moment prečnega prereza nosilca glede na njegovo težišče. Nosilec je v osi obtežen z linijsko obtežbo q Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 11 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. in linijsko momentno obtežbo mz', obe sta merjeni na enoto dolžine nedeformiranega nosilca. Komponenta q v globalni smeri X je qx, komponenta v globalni smeri Y pa qy. A je rotacijska matrika, ki pri matričnem množenju z leve zasuka vektor v protiurnem smislu; zapiše se: . r coso? — sin ip 1 A = . (2.4) simp cos (p Rotacijska matrika A je ortogonalna matrika, zato med drugim velja A-1 = AT, kar bomo v nadaljevanju pridoma izkoriščali. Dodatno velja še: A' = WAy/, (2.5) kjer smo z apostrofom (•) = d (•) /ds označili odvajanje po ločni koordinati s; W je antisimetrična matrika: -tir \ 0 — 1 1 VV = . . . (2.6) 1 U Omenimo še, daje množenje matrik A in W komutativno, tudi če je katerakoli od matrik transponirana. Ko vpeljemo še oznaki i? in G za elastični in strižni elastični modul, lahko togostne in vztrajnostne lastnosti prečnega prereza zapišemo z matrikama I" E A 0 1 0 GAs E A 0 CE = ) (2-7) 0 GAs Ap 0 kjer je I identitetna matrika velikosti 2x2. r Ap o i , . .. CP = r. ^ = ^P i) (2-0) o Ap 2.3 Osnovni izreki in definicije Če ne upoštevamo termomehanskih zakonov, ki niso predmet disertacije, za obravnavani problem veljajo štirje osnovni izreki na katere se bomo v nadaljevanju sklicevali in jih zato na tem mestu zapišemo. Izrek o ohranitvi mase pravi, da se masa poljubnega dela konstrukcije, če mu nič mase ne dodamo in nič ne odvzamemo, ne spreminja. Iz tega sledi, daje p dV = Po dVo, (2-9) kjer sta p in dV gostota materiala in infinitezimalni delec volumna v trenutni (deformirani) legi ter p0 in dVo isti količini v začetni (nedeformirani) legi. Izrek o ohranitvi energije pravi, da se mehanska energija sistema ohranja, če na sistem deluje le konzervativna obtežba, lahko pa se spreminja njena oblika (npr. iz kinetične v deformacijsko). Mehanska energija sistema se lahko ohranja tudi, če je sistem podprt. Izrek o gibalni kolicˇini pravi, daje odvod gibalne količine p sistema (telesa, nosilca) v vsakem trenutku enak rezultanti zunanjih sil Urez, ki delujejo nanj: p = Rrez- (2.10) 12 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Gibalna količina je definirana kot L p= / Apuds. (2.11) o Za poljuben algoritem, ki želi ohranjati gibalno količino neobremenjenega sistema, mora veljati Ap = 0. Ker iz tega pogoja sledi p = 0, lahko iz enačbe oz. izreka o gibalni količini (2.10) razberemo, da mora biti rezultanta zunanjih sil enaka nič, če želimo imeti ohranjanje gibalne količine. Ker podpore na konstrukcijo predstavljajo obtežbo, je očitno, da se v primeru podpor gibalna količina v splošnem ne ohranja. V primeru ohranjanja energije so podpore lahko prisotne, saj na mestih podpor ni pomikov in s tem je potencial zunanjih sil na teh mestih enak nič. Izrek o vrtilni kolicˇini pravi, daje odvod vrtilne količine L sistema (telesa, nosilca) v vsakem trenutku enak rezultanti zunanjih momentov Mrez, ki delujejo nanj: L = Mrez. (2.12) Vrtilna količina sistema je definirana kot L L = I (r x p + J p (p) ds. (2.13) o Vektorski produkt v enačbi (2.13) je zapisan za dvorazsežna vektorja in je definiran v naslednjem razdelku. Če želimo doseči ohranjanje vrtilne količine konstrukcije, ta ne sme biti podprta, kar lahko pokažemo s podobno izpeljavo kot pri gibalni količini. Definicija vektorskega produkta v ravnini. Matematična operacija vektorskega množenja je definirana v trirazsežnem prostoru. Zato vpeljemo definicijo vektorskega produkta za množenje dvorazsežnih vektorjev, ki jo bomo uporabljali v tem delu. Definicijo vpeljemo s štirimi generičnimi simboli a, b, c in d : [al Tel Ta 6] , x , = det = a d — b c. (2.14) o d c d Vektorsko množenje vrne skalar, ki pove velikost (dolžino) vektorja, pravokotnega na ravnino vektorjev množiteljev. Ker vektorja množitelja ležita v ravnini (X, Y), je smer pravokotnega vektorja v globalni smeri Z. Vektorsko množenje pa lahko zapišemo tudi z vpeljavo prej opisane matrike W: [al [ c 1 [al ^-,tT [ c 1 , x , = , W , b d b " W • (2.15) d V besedilu po potrebi uporabljamo oba zapisa - vektorskega in matričnega. Deformacijska energija ravninskega nosilca je definirana kot: 1 fL 2 2\ 1 fL ildef = - (EA e + GAs 7 + E J n ) ds = -i (e ce L + E J n ) ds. (2.16) 2 7o 2 0 Potencial zunanjih sil je definiran kot: L nzun = — S U* — / (q u + rriiLp) ds, (2.17) o kjer je S = [Si, S2, Ss, S4, S$, Sq] vektor točkovnih sil v krajiščih nosilca (pri s = 0 in s = L), U* pa vektor pomikov v krajiščih nosilca, ki ima enako obliko in smer kot vektor točkovnih sil. Potencial zunanjih sil je po definiciji negativno delo. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 13 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Potencialna energija je vsota deformacijske energije in potenciala zunanjih sil: ripot = ndef + nzun. (2.18) Kinetična energija je definirana z naslednjim izrazom: 1 fL 1 fL Hon = — (Ap u + Ap v + E J (p ) ds = — \ (u Cp u + E J (p ) ds. (2.19) 2 7o 2 0 Mehanska energija oziroma totalna energija je definirana kot vsota kinetične in potencialne energije: n = nkin + npot. (2.20) 2.4 Vodilne/glavne enačbe Dinamične ravnotežne enačbe. Če si produkt mase in pospeška predstavljamo kot negativno linijsko obtežbo, lahko dinamične ravnotežne enačbe interpretiramo kot ravnotežne enačbe, t.j. kot zahtevo, da je vsota vseh sil in vsota vseh momentov na sistem enaka nič. Za konstrukcije lahko ravnotežne enačbe zapišemo na dva različna načina: lahko jih zapišemo v obliki, ki velja v vsaki točki na nosilcu. Taki obliki pravimo krepka. Lahko pa jih zapišemo v obliki, ki zahteva, daje ravnotežje zadoščeno v nekem povprečnem smislu. To je t.i. šibka oblika. Krepke ravnotežne enačbe se zapišejo za posamezno točko konstrukcije in imajo obliko diferencialnih enačb: (AN)' + q =Apü, (2.21) M' + r' x AN+mz = J pip. Z nadpisnimi pikami označujemo odvajanje po času t. Primer šibke oblike ravnotežnih enačb je princip virtualnega dela, ki se glasi: L / [N 5e + M on — (q — ücp) öu — (niz — Jpfi) Sip\ ds = 0. (2.22) o V enačbi (2.22) smo z 5 označili variacijo. V analitičnem smislu sta krepka in šibka oblika enačb enakovredni: če velja šibka, velja tudi krepka in obratno. Ko količine v teh enačbah diskretiziramo, pa ta enakost ne velja več. Dodatno lahko ugotovimo, da princip virtualnega dela velja v integralskem smislu, krepke ravnotežne enačbe pa eksplicitno v vsaki točki. Tretja, iz numeričnega vidika najbolj pomembna razlika je, da v krepki obliki enačb (2.21) nastopajo za eno stopnjo višji odvodi kot v šibki (2.22). To implicira zahtevo zveznosti prvih odvodov neznanih količin (C\ zveznost) in ne samo zveznosti neznanih količin samih (Co zveznost). Zienkiewicz, Taylor in Zhu (2005) navajajo, da se pogosto zgodi, da formulacije s šibkimi ravnotežnimi enačbami (z implicirano Cq zveznostjo) dajejo bolj ‘fizikalno realistične’ rezultate. Matrični produkt AN, ki nastopa v enačbi (2.21), bomo včasih nadomestili z oznako R. R je vektor projekcij rezultant napetosti na globalni smeri X in Y (shematski prikaz je tudi na sliki 2.1, desno zgoraj): AN = R = [Kx, -fvrj • {2"^5) 14 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Kinematične enačbe. Naše izhodišče za opis zveze (vezi) med kinematičnimi in deformacijskimi spremenljivkami so Reissnerjeve enačbe (Reissner (1972)). Zapišemo jih v naslednji obliki: oziroma v komponentni obliki: e 7 L = A (u + roj + Ci = A r +ci, [u + X ) cos (p + (v + Y ) sin ^ - 1 - [u + X ) sin (/? + (f + r J cos tL> (2.24) (2.25) Smiselno je zapisati še obratne zveze, kjer so pomiki izraženi z deformacijami, saj se bomo v nadaljevanju nanje sklicevali: r' = A (e - ci), (2.26) (p' = K. Enačbe (2.24) oz. (2.25) so izpeljane in veljavne ob naslednjih predpostavkah: • pomiki in deformacije so po velikosti neomejeni1; • prečni prerez ostane po deformiranju nosilca raven in nedeformiran (tog); • prečni prerez po deformiranju nosilca ni več nujno pravokoten na os nosilca; • enačbe upoštevajo vpliv osnih, strižnih in upogibnih deformacij na pomike; • enačbe veljajo za ravne in ukrivljene nosilce. Ker so enačbe zapisane v ravnini, velja še: • prečni prerez je poljuben, toda simetričen glede na ravnino (X, Y); • materialne lastnosti se po prerezu spreminjajo poljubno, toda simetrično glede na ravnino (X, Y) ; • gibanje izven ravnine (X, Y) ni možno. Konstitucijske enačbe. S konstitucijskimi enačbami iz znane deformacije in materialnega zakona v točki osi nosilca določimo konstitucijske notranje sile v prerezu. Zveze matematično zapišemo Nc = <7nn(e) dA, J A Qc = / Crai (7) dA, J A Mc = / (Tnn(e) z d A. (2.27) A ann in a nt označujeta napetosti v ravnini prereza, določeni z normalo en (slika 2.1). Napetost ann deluje v smeri normale, napetost ant pa v smeri tangente na ravnino. To sta edini od nič različni napetosti in sta funkciji le deformacij. 'Edina omejitev je fizikalna, ki preprečuje, da bi se material stisnil v točko: e > —1 _ Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 15 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Napetosti v geometrijsko nelinearnih problemih v splošnem niso trivialna reč, saj smo dolžni uporabljati take napetosti, ki so fizikalno komplementarne deformacijam. Fizikalna komplementarnost pomeni, da je produkt deformacij in napetosti enak deformacijski energiji. Reissnerjeve deformacije so energijsko komplementarne 2. Piola-Kirchhoffovim napetostim. Ker za obravnavane materiale poznamo materialni zakon na nivoju 2. Piola-Kirchhoffovih napetosti ann in Reissnerjevih deformacij e, so te zveze izhodišče za račun konstitucijskih spremenljivk. Za odnos med strižnimi deformacijami in napetostmi predpostavimo linearno elastični zakon. Bolj podrobno bomo funkcije (0) — L n ds = 0, T (3.12) R(0)+ [Si, S2J = 0, M(0) + S3 = 0, fL •• j r o o iT r. R(0) + ln (Apu — q) as — 04, 05 = U, /~v iL ["/ \TtttTtvt "1 7 rf m(U) + In (e — ci)W JN— m z + Joti1 as — 06 = 0. Ob teh enačbah moramo upoštevati tudi zveze, podane v enačbah (3.8) in (3.9), ki pa so le predpisi za račun pomikov in sil v odvisnosti od deformacij. Edine neznane funkcije koordinate s, ki nastopajo v enačbah, so deformacije e in k; v enačbah pa nastopa še devet neznanih parametrov: u(0), (0), R(0), M(0), u(L), ip(L). Vse omenjene količine so funkcije časa t. Enačbe rešujemo z metodo končnih elementov, kot je opisano v naslednjem poglavju. 3.2 Metoda uteženih ostankov Z metodo uteženih ostankov lahko enačbe problema (3.12) zapišemo v integralski obliki, kar ima, kot navajajo Zienkiewicz, Taylor in Zhu (2005), določene prednosti. V skladu s to metodo sistem enačb (3.12) pomnožimo s poljubnimi utežmi, ki jih označimo z H in 0, ter integriramo po območju veljavnosti enačbe. Uteži H = H(s) so funkcije kraja, uteži 0 pa so skalarji. Ker sta le prvi dve enačbi veljavni na območju s G [0, L], integriramo le njiju. Ostale enačbe so veljavne le v enem od krajišč nosilca, zato integriranje ni potrebno. Tako dobimo: Jo (.-N — CE L) ai ds = 0, rL j-^ \ —1 L, (M — tiJ k) Lio ds = 0, U ' ^ I rL { 4 / \ I \ \ r^ u(L) — u(0)— |n A (e — Ci) — rn) ds ) fe)i = 0, Uv *-> u L (tp(L) — ü?(0) — Cf k ds) 0o = 0, ' V T u ^ ( T\t (3.13) (R(0)+[Si, 6*2] j ©3 = 0, (M(0) + 6*3) ©4 = 0, [ tt(L) — R(0) — |n (Apu — q) ds) Ö5 = 0, U ' ( r\\ rL \, \TiirTi\T t •• 1 \ r\ [ M(L) — M(U) — L (e — ci) W JN— m z + J p(p\ ds) Hr = 0. Z izbiro različnih utežnih funkcij dobi metoda uteženih ostankov različna imena. Če vpeljemo za utežne funkcije enake nastavke kot za neznane funkcije, se metoda imenuje Galerkinova metoda. Če vpeljemo Diracove delta funkcije, je metoda kolokacijska. 20 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 3.3 Galerkinova metoda končnih elementov Interpolacija po kraju. Za interpolacijo neznanih funkcij uporabimo Lagrangeve polinome stopnje nap, ki jih označujemo s Pi, i = 1,..., mtp : r e 1 v^ ti r Li 1 = > Fi K ^—' Ki i=l (3.14) Lj in Ki označujeta diskretne vrednosti deformacij v interpolacijskih točkah, ki jih ekvidistantno razporedimo po dolžini elementa: Si = (i — 1)——j. Ker želimo izpeljati Galerkinovo metodo končnih elementov, izberemo enak nastavek tudi za interpolacijo uteži: — = > K- „ -2 ^ -2j (3.15) Uteži ®i so konstante in ne funkcije, zato jih ni potrebno aproksimirati. Ko interpolacijska nastavka (3.14) in (3.15) vstavimo v enačbe (3.13), dobimo: J0 (N - cE }^i=i pi Li) l^j=i pj ai,j ds = °> J0 (M — L/J 2^=1 ^i Ki) Z^7=l ^J ^2,J "S = U, ( u(L) — u(0)— L (A (e u v ci ds ©i = 0, 3 = 0, (M(0) + 6*3) 04 = 0, T-» / T-» fL / a 7 \ y^. ril L) — K(0) — ln (Apu — q) as) Ö5 = 0, v (jvr 1 li/r/ r\\ rL Vi \TtttT t ••! 7 \ r\ MIL) — M(U) — ln e - Ci VV N-mz + Jaw as) 06 = 0. (3.16) Ker so uteži Hij, ^2,j in ®i poljubne vrednosti, morajo biti členi, s katerimi so te uteži pomnožene, enaki nič, da bo enačbam (3.16) identično zadoščeno. Tako dobimo sistem 3 n\tv + 9 enačb za 3 n\tv + 9 neznank: rL /tvt \—V^itD \ 7 Jo (-N ~~ ce l^i=i Pi Li) Pj as= 0, j = 1,..., nap fL 171 /r 771 t v^raitp \ A n J0 [M. — t.J z^i=i Pi Ki) Pj as = U, j = 1,..., nitp fL I » / \ / \ 7 u(L) — u(0) — /n A ( L — Ci ) — rn os = 0 U v x/ u 0 0 (3.17) tL>(L) — o?(0) — L n ds = 0 ' T' U (0j+ [ol, O2J = Ö M(0) + 5*3 = 0 T-» / X-» r^ 7 KlL) — K(0) — ln (Apu — q) OS = 0 v v Ur 7iT/r\ 71 /rIr\\ fL Ti \TtttTivt ••! 7 Mit) — Milj) — /n (e — Ci ) VV JN— mz + J pipi as = 0 Zapisani sistem enačb je semi-diskreten sistem diferencialno-algebrajskih enačb. Pred reševanjem moramo vpeljati še diskretizacijo po času in s tem sistem prevesti na algebrajskega. Opomnimo, da dobljeni sistem enačb nima standardne oblike M T + K(T)T + F(T) = 0, kjer bi M predstavljal masno matriko, K togostno matriko, F residual in T vektor neznanih količin. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 21 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Interpolacija in integracija po času. Za interpolacijo po času uporabimo Newmarkovo (Newmark (1959)) časovno integracijsko shemo (o časovni integraciji bo v nadaljevanju še veliko govora, predvsem v naslednjem poglavju). V enačbah nastopajo prvi in drugi časovni odvodi količin Lj, «j, u(0), (0), ki jih po času interpoliamo z Newmarkovimi diskretizacijskimi nastavki: Li,n+l = Li,n H----1~((1 — 2/3t) Litn + 2ßt Litn+l) Î = 1,¦ ¦ ¦ , ^itp Ki,n+1 = Ki,n H----1~((1 — 2/3+) liitn + 2/3+ Hi,n+l) Î = 1,¦ ¦ ¦ , ^itp u(0)ra, 1 = u(0)ra H—1~((1 ~~ 2/3+) u(0)jn + 2/3+ üj)n_|_i) )ra+l) (3.18) v(0)n+i = (0)„ H—1-((1 — 2/3+) (o(0)„ + 2/3+ ó?(0)r Li,n+l = Li,n + At((l — 7+) Litn + 7+ Li,n+l ) I = 1,. . . , Tlüp Ki,n+1 = Ki,n + At((l — 7+) Kj>n + 7+ K^n+l) 2 = 1,..., «itp U(0)ra , 1 = U(0)in + At((l — 7+) Üi>n + 7+ Üi)n-|-i) ^(0)n+i = (0)ra + At((l — 7+) (0)n + 7+ (p(0)n-\-i) (3.19) Indeks (»)raoznačuje znane (skonvergirane) vrednosti pri času tn, indeks (•)n+1 pa vrednosti pri času tn+i = tn + At, kjer je At velikost časovnega koraka. ßt in 7+ sta parametra Newmarkove časovno integracijske sheme. Pri vrednostih ßt = 1/4 in 7+ = 1/2 dobimo znano in pogosto uporabljano trapezno pravilo s konstantnimi pospeški. Enačbi (3.18) in (3.19) se pri implementaciji uporabljata v obrnjeni obliki, saj iz enačbe (3.18) iz znanih vrednosti (•)n+1 računamo pospeške (•)n+1, ki jih nato vstavimo v enačbo (3.19), da dobimo hitrosti (*)n+l- 3.4 Numerični vidik formulacije Ko poznamo diskretne vrednosti deformacij v integracijskih točkah Li}U+\ in K^n+i, lahko izvrednotimo vse ostale količine. Problem je v tem, da ti računi niso preprosti, saj zahtevajo izvrednotenje večkrat vgnezdenih integralov. Denimo, da želimo izračunati vrednosti projekcij sil R v i-ti točki. Po enačbi (3.9) je za to potrebno integrirati pospešek ü, ki ga dobimo z dvakratnim odvajanjem pomika u iz enačbe (3.8) po času : ril Si) = K(0) + ln (Apu — q) as = = R(0) + T/i Apü(0)+-772 fn (A (e — ci) — rn ds — q ds (3.20) Ker je potrebno vse integrale izvrednotiti numerično, vpeljava dvojne integracije drastično poveča število računskih operacij, saj je potrebno za izvrednotenje v vsaki integracijski točki zunanjega integrala izvrednotiti celoten notranji integral. Število potrebnih računskih operacij za izvrednotenje enojnega integrala, Tijtg, naraste na n^g + nft operacij za izvrednotenje integrala z enim vgnezdenim integralom. Pri izvred-notenju izraza v enačbi (3.20) se dejansko pojavi še eno vgnezdenje, ker v rotacijski matriki A nastopajo zasuki, ki se jih računa z integracijo upogibnih deformacij po enačbi (3.8). V seštevku je najbolj zapleten člen v togostni matriki kar petkrat gnezden. To je računsko izredno potratno in algoritem je v takšni obliki nekonkurenčen drugim formulacijam. Pomembno je tudi vprašanje stopnje integracije n^g, saj v izrazu nastopajo kombinacije trigonometričnih funkcij, za katere ne moremo vnaprej določiti zadostne stopnje integracije in ne preostane nam drugega, kot da jo določimo s poskušanjem. 22 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 3.5 Numerični primer Zasuk robotske roke. Račun obravnava t.i. dvig podajne robotske roke, ki sta ga za testiranje zmogljivosti in natančnosti formulacij geometrijsko točnih nosilcev predlagala Simo in Vu-Quoc (1986). 10 m dolga podajna ravna robotska roka leži na tleh, zatem pa se jo v vpetem vozlišču zasuka okoli navpične osi za 1.5 radianov v 5 sekundah, kot je prikazano na sliki 3.1. Materialni parametri primera so: E A = GAS = 10000 N, EI = 500 N/cm2, (3.21) Ap = 1 kg/cm, J p = 10 kg cm, zasuk v vpetišču v odvisnosti od časa pa je prikazan na sliki 3.1, desno. Slika 3.1: Geometrija problema, levo; zasuk v vpetišču v odvisnosti od časa, desno. Figure 3.1: Geometry of the problem, left; rotation at pinned end vs time, right. Problem rešimo s štirimi končnimi elementi s konstantno interpolacijo deformacij in 5-točkovno Lo-battovo integracijo, za časovno integracijo uporabimo trapezno pravilo ( ßt = 1/4 in jt = 1/2) ter konstantni časovni korak ?t = 0.5 s. Rezultati so prikazani na sliki 3.2, kjer so za primerjavo vrisani tudi rezultati, ki sta jih dobila Simo in Vu-Quoc (1986) s končnimi elementi z linearno interpolacijo pomikov ter s selektivno reducirano integracijo zaradi preprečitve blokiranja. Opazimo lahko, da se rezultati obeh formulacij dobro ujemajo. Do blokiranja deformacijskih elementov ni prišlo. Več numeričnih testov te formulacije lahko bralec najde v članku Gamsa s sodel. (2007a) v prilogi A. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 23 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 3.2: Odziv sistema in primerjava z rezultati Sima in Vu-Quoca (1986). Figure 3.2: Response of the system and comparison to results of Simo and Vu-Quoc (1986). 3.6 Sklepi Izpeljali in implementirali smo končne elemente, ki temeljijo na interpolaciji deformacijskih količin kot edinih neznanih funkcijah problema. V formulaciji so uporabljene točne kinematične zveze, v katerih upoštevamo vpliv osnih, upogibnih in strižnih deformacij. Na numeričnih primerih seje pokazala velika natančnost teh elementov tudi v primerih mrež z majhnim številom končnih elementov. Ker elementi ne blokirajo, ni potrebno uporabiti nobenih posebnih ukrepov za odpravljanje blokiranja, kar je prednost. V računih seje pokazalo, da kombinacija Newmarkove časovne integracijske sheme in elementov z interpolacijo deformacij ni brezpogojno stabilna. To se pokaže z rastjo mehanske energije sicer konzervativne konstrukcije ter končno v izgubi konvergence Newton-Raphsonovega postopka reševanja sistema nelinearnih algebrajskih enačb. Kot zadnje pa omenimo, daje formulacija s stališča računskih časov inferiorna v primerjavi z drugimi formulacijami zaradi gnezdenja integralov. Zadnji dve hibi formulacije poskušamo odpraviti v naslednjem razdelku. Problem časovne nestabilnosti občutno zmanjšamo, ne pa tudi odpravimo, z uporabo drugačne časovne integracijske sheme, numerično zahtevnost pa bistveno zmanjšamo z drugačnim pristopom k reševanju. 3.7 Izboljšava elementov z interpolacijo deformacij Izboljšava numeričnega postopka. S posebnim ukrepom se je možno izogniti gnezdenju integralov v izrazih. Ideja, ki nas pripelje do tega ukrepa, posnema ‘filozofijo’ numerične integracije. Pri numerični integraciji integrand interpoliramo z interpolacijskimi funkcijami P* kot linearno kombinacijo interpolacijskih funkcij in vrednosti integranda. Ta interpolacija, zapisana na poljubni generični 24 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. funkciji /(L), se glasi: »dtp /(C) = / f(L,i)Pi(0- (3.22) i=\ S L smo označili neodvisno spremenljivko, s ^ pa njene vrednosti v izbranih (zaenkrat poljubnih) interpolacijskih točkah. Opozorimo, daje na desni strani enačbe (3.22) le P*(L) funkcija koordinate L. Pri integriranju enačbe (3.22) po spremenljivki L na območju L G [—1,1], /(Ci) predstavljajo konstante. Integracijo interpolacijskih funkcij P*(0 opravimo le enkrat, ker se ne spreminjajo in njihove integrale označimo z Wi. Tako se integral enačbe (3.22) zapiše: ç\ ltp / /(C) dL, = y u>i /(CJ- (3.23) ?—i S pravilno izbiro interpolacijskih funkcij in interpolacijskih točk lahko izboljšamo natančnost integracije. Za integracijo polinomov je optimalna izbira Gaussova integracija, ki pri danem številu interpolacijskih točk mtp točno zintegrira polinom najvišje možne stopnje. Vendar pa pri tej integraciji ni uteži pri koordinatah L = ±1, kar si pri materialno nelinearnih nosilcih želimo. Dobra izbira integracije, ki vsebuje integracijske točke tudi pri L = ±1, je ti. Gauss-Lobattova integracija oz. skrajšano Lobattova integracija. Vsa sodobna matematična okolja (Mathematica, Matlab) uporabljajo adaptivno določevanje stopnje integracije, da lahko izračunajo integrale do zahtevane natančnosti. V metodi končnih elementov takih prijemov praviloma ne uporabljamo. V primeru naše formulacije se zgoraj opisana ideja numerične integracije implementira na naslednji način: začnemo z interpoliranjem odvodov pomikov iz enačbe (3.1): »dtp u'(s) = > Pi (s) (A (e — ci) — rn).. (3.24) i=\ Ta interpolacija ni običajna, saj ne interpoliramo direktno osnovnih neznanih funkcij kot npr. v enačbah (3.14) ali (3.15), temveč iz njih izpeljane izraze. Z integriranjem enačbe (3.24) dobimo pomike: «itp wis) = u(0) + > Pi(s) (A (e — ci) — rn).. (3.25) i=\ Z oznako Pi{s) = J0 Pi at, za i = 1,..., nnp smo označili integrale Lagrangevih polinomov. Zasuke dobimo direktno z integracijo zadnje izmed enačb (3.1): »«p tL>(s) = (0) + 2, Pi(s) Ki- (3.26) i=\ Nadalje interpoliramo tudi odvod ravnotežnega momenta, ki je zapisan v ravnotežnih enačbah (3.3): »dtp ( l M'=y Pi (s) ((e — ci) W N—niz + Jp&) • • (3.27) i=1 Po integraciji dobimo nitp M = M(0) + > Pj(s) ((e — Ci) W N—niz + Jp& • • (3.28) i=\ Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 25 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Z enačbami (3.25), (3.26) in (3.28) dobimo pomike, zasuke in momente v poljubni točki samo z integriranjem Lagrangevih polinomov, kar je potrebno storiti le enkrat. Interpolacije, ki jih opisujemo v tem poglavju, so le nastavki za numerično izvrednotenje integralov in jih ne smemo šteti za interpolacijo v smislu metode končnih elementov. Analogno kot v primeru pomikov, zasuka in momenta postopamo pri računu pospeškov, ki so vgnezdeni pri računu notranjih sil R in M: Ü(s) = Ü(0) + ^i=l Pi(s) Tip (-A- (e — cl) — ro) ' (p(s) = (p(0) + ^i=i 1Pi(s) Ki- (3.29) Po integraciji dobimo: (3.30) J0 Ü(s) ds = Ü(0) S + ^i=l nPi(s) AFI (-A- (e — ci) ~~ ro) • J0 (p(s) ds = (p(0) S + ^i=l Pi(s) Kii TT T~t ( rS TČ T~t i i/- ^ kjer smo z FAs) = /n /n F; dri at označili druge integrale Lagrangevih polinomov. ^ U U ^ ' ^ Kolokacijska metoda končnih elementov. Zaradi želje, da še bolj zmanjšamo število integralov po nosilcu, vrsto metode spremenimo v kolokacijsko. Motiv je tudi v tem, da želimo ves čas delati le z vrednostmi v interpolacijskih točkah, ki bi se ujemale s kolokacijskimi točkami, kar bi še nadalje poenostavilo numerični postopek. Interpolacijske in posledično tudi kolokacijske točke so še vedno ekvidistantno razporejene po nosilcu. Nastavki za interpolacijo deformacij ostanejo enaki kot prej, uteži H pa sedaj aproksimiramo z Diracovimi funkcijami: r Hi i r Eij i ,_, = Oj{s) , za j = 1,..., nitp (3.31) '—2 '—'2j kjer je 5i Diracova funkcija: oo: sicer r / °°; ce s = Sj djls) = . (3.32) 0; in zanjo velja J0 oj(sj as = 1 za j = 1,..., n^p. S s.,- za j = 1,..., n^p, smo označili koordinate kolokacijskih točk. Zaradi ujemanja kolokacijskih in interpolacijskih točk pri i = j velja Sj = Sj. Enačbe (3.13) z vpeljavo utežnih funkcij (3.31) postanejo: (N — cg Si=i -P» Li) ^i,i = 0) J = 1) • • • ; ^itp (M — _BJ ^i=i -P» Ki) ^2j =0, j = 1,..., nitp ( u(Lj — u(0j— J0 ^A (e — cij — r0J as j fe)i = 0, ((0) — L k ds = 0 (3.34) (0j+ [ol, 02J = 0 M(0) + 5*3 = 0 T-» T-» r^ 7 r\ tt(L) — K(0) — /n (A/9U — q) as = U v U ' n\ rL ["/ \TtttTtvt •• 1 7 M(L) — M(U) — ln (e — ci ) W JN— m z + Joti1 as = 0 Zopet smo dobili sistem 3 n;tp + 9 enačb za 3 nap + 9 neznank, ki pa nima standardne oblike. Omenimo še to, da smo na nivoju končnega elementa iz formulacije izločili vsakršno integracijo, le Lagrangeve polinome je potrebno že pred začetkom analize dvakrat zaporedoma integrirati. Količine Pi in Pi so odvisne od dolžine elementov in zato jih je potrebno pripraviti za vsak element posebej. Časovna integracija. Za interpolacijo po času tokrat uporabimo t.i. sredinsko oz. midpoint časovno integracijsko shemo (npr. Stander in Stein (1995)). O midpoint integraciji bomo podrobno govorili v naslednjem poglavju, zato naj omenimo le, da ob pravilni uporabi midpoint časovne integracije postanejo nelinearni končni elementi, osnovani na pomikih, brezpogojno stabilni, ker zadoščajo kriteriju ohranjanja energije. Preveriti smo želeli, če ta ugodna lastnost velja tudi za deformacijske končne elemente. Pri teh elementih moramo poleg časovnih odvodov, ki jih računamo po sredinskem pravilu, zadoščati tudi enačbam dinamičnega ravnotežja pri sredinskem času tm = tn + At/2. Pri vrednotenju količin pri sredinskem času pa je potrebno skoraj praviloma uporabiti aritmetično povprečje. Diskretizacijski nastavki so: u(0)TO = 2 vLi,n i L-i,n-\-l) i \ (Hi,n + Ki,n+l) , U(°)i,n + U(°)i,n+1 (3.35) (0)i,ra+l Li,m = ÂI (Li,n+l — Li,n) , V ) Y>(0) Ai U(°)i,n+1 - U(°)i,n Y>(0)*,n+l - Y>(0)*,n (3.36) i __ 2 /_ __ _ ~i,m Ai2 V L,^+1 ^i,n At ki Ai2 V^^ro+l ^i,n ^L l^i ü(0) 2 Ai2 2 Ai2 u(0)in_|_1 — u(0)in — At u(0)r ^(0)^n,1 — ^(0)^n — At (f(0) (3.37) _ _ i _ _ _ Indeks (»)TOoznačuje vrednosti pri času tm. K tem nastavkom moramo dodati še nastavek za popravek Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 27 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. hitrosti: ei,n+l = Ät (ehn+l — Li,n) ~ &i At K ijit U(°)n+1 = h (U(°)i,n+1 - U(°)i,n) ~ U(°)n' ^(0)n+i = s (^(0)^^+1 - Y>(0)ijTl) - Y?(0)n. (3.38) 3.8 Numerični primer Nihanje konzole. V tem računskem primeru na ‘preprostem’ primeru konzole primerjamo originalno in izboljšano formulacijo z interpolacijo deformacij. Geometrijski podatki konzole so namerno izbrani tako, da bi formulacija, osnovana na pomikih, blokirala (Rong in Lu, 2003). Slika 3.3: Računski model konzole z obtežbo, levo; poveš prostega konca (vb) v odvisnosti od časa, desno. Figure 3.3: Cantilever beam model and load, left; time-history response for tip deflection vb, right. Konzola (slika 3.3, levo) je na prostem koncu obremenjena s silo 250 N, ki jo hipoma nanesemo na konzolo. Prečni prerez konzole je poln pravokotnik. Po Rongu in Luju (2003) velja, da formulacija z interpolacijo pomikov blokira, če velja log10(L/iy) > 1, kjer je iy = ^ Iy/A vztrajnostni polmer prereza. V skladu s tem pogojem smo izbrali razmerje log10(L/iy) = 2, dolžino konzole L = 1 m in širino prečnega prereza b = 0.01 m. Višino prečnega prereza h smo dobili z računom h = L^/VŽ/IO2 o 0.03464 m. Elastični modul materiala je E = 2.1 x 1011 N/m2, kar ustreza jeklu. Za račun odziva smo uporabili časovni korak At = 0.0005 s, odziv pa smo računali do časa 0.1 s. Poveš prostega konca konzole v odvisnosti od časa je prikazan na sliki 3.3, desno. Primer smo reševali z osnovno (Galerkinova metoda, konstantna interpolacija deformacij, Newmarkova časovna integracijska shema) in izboljšano metodo (kolokacijska metoda, več stopenj interpolacije deformacij, midpoint časovna integracijska shema) ter rezultate strnili v preglednicah 3.1 in 3.2. Omenimo še, da smo pri osnovni verziji uporabljali Lobattovo integracijo 5. stopnje, v izboljšani formulaciji pa klasična numerična integracija ni potrebna. V preglednicah je z oznako ‘FEx’ označeno število končnih elementov in uporabljena stopnja interpolacije; 2 FE3 tako pomeni analizo, kjer je bila konstrukcija modelirana z dvema končnima elementoma s kubično interpolacijo deformacij. Z ‘DOF’ je označeno število prostostnih stopenj, relativno napako pa 28 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Preglednica 3.1: Izboljsˇana formulacija. Table 3.1: Numerically improved formulation. v —— 0.05 s L = 0.1 s FE DOF v (L) [m] Rel. nap. f (L) [m] Rel. nap. Max. v(L) [m] Rač. čas [s] 1 FE2 15 0.02296 4 v 10-02 0.00397 4 v 10-01 0.02405 4.0 1 FE3 18 0.02168 Owi n-02 0.00790 1 v 10-01 0.02400 4.7 1 FE4 21 0.02212 1 v 10-04 0.00708 1 v 10-02 0.02360 5.3 2 FE2 30 0.02177 9 v 10-02 0.00687 9 v 1 O-02 0.02404 6.5 2 FE3 36 0.02214 8 v 10-04 0.00697 5 v 10-03 0.02368 7.7 2 FE4 42 0.02216 9 v 10-03 0.00700 4 v 1 O-04 0.02372 9.3 4 FE2 60 0.02208 9 v 10-03 0.00705 7 v 10-03 0.02382 11.5 4 FE3 72 0.02213 4 v 10-04 0.00700 0 y 1 n -04 0.02375 14.3 4 FE4 84 0.02212 r y 1 n-05 0.00700 9 v 10-05 0.02377 17.1 20 FE3 360 0.02212 0.00700 0.02377 64.4 Preglednica 3.2: Osnovna verzija formulacije. Table 3.2: Basic version of the formulation. v —— 0.05 s L = 0.1 s FE DOF v (L) [m] Rel. nap. v (L) [m] Rel. nap. Max. v(L) [m] Rač. čas [s] 1 FE0 9 0.00672 7 v 1 O-01 0.01688 1 x 100 0.01806 23.6 2 FE0 18 0.02208 9 v 1 O-03 0.00039 9 x 10-01 0.02256 41.8 3 FE0 27 0.02209 1 v 1 O-03 0.00317 r w 1 n-01 0.02338 60.9 4 FE0 36 0.02231 9 v 1 O-03 0.00506 3 V 1 O-01 0.02357 77.0 5 FE0 45 0.02177 9 v 1 O-02 0.00540 9 V 1 O-01 0.02376 98.4 6 FE0 54 0.02181 1 v 1 O-02 0.00610 1 v 1 O-01 0.02385 115.8 8 FE0 72 0.02200 5 v 1 O-03 0.00684 9 V 1 O-02 0.02393 148.7 10 FE0 90 0.02205 0 y 1 n -03 0.00699 0 y 1 n-05 0.02388 176.9 40 FE0 360 0.02211 0.00699 0.02379 743.0 smo računali glede na najbolj točno rešitev, ki smo jo v osnovni verziji formulacije dobili z analizo 40 FE0 in v izboljšani verziji z 20 FE3. Na podlagi rezultatov v preglednicah 3.1 in 3.2 lahko ugotovimo: (i) obe formulaciji konvergirata k istim rešitvam, s čimer se potrdi (validira) numerični postopek; (ii) hitrost konvergence rešitve narašča z naraščajočim redom interpolacije; (iii) natančne rešitve dobimo z zelo majhnim številom končnih elementov; (iv) blokiranja ni in (v) izboljšana formulacija je približno 10-krat hitrejša kot osnovna, kar predstavlja občuten napredek v učinkovitosti formulacije. Dodatni računski primeri so prikazani v članku Gamsa, Planinca in Sajeta (2007b), ki je dodan kot priloga B. 3.9 Sklepi Razvili in implementirali smo metodo končnih elementov za dinamiko ravninskih elastičnih nosilcev z interpolacijo deformacijskih količin. V formulaciji zato inherentno nastopajo konsistenčne enačbe v eksplicitni obliki, česar, kot bomo videli v nadaljevanju, pri elementih z interpolacijo pomikov ni. Bistvena Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 29 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. posledica eksplicitnega zadoščanja konstitucijskih enačb je, da elementi ne blokirajo. Dodatno se od deformacijskih elementov pričakuje, da v primeru nelinearnih materialnih modelov bolj konsistentno računajo notranje sile, saj so prek osnovnih spremenljivk - deformacij podane z zveznimi funkcijami in ni pojava ‘nenatančnega’ računa odvodov pomikov in zasukov, ki jih potrebujemo za račun notranjih sil. S tem se izognemo ekstrapolaciji notranjih sil iz ti. Barlowih točk (Barlow (1976)). Elementi poleg odpornosti na blokiranje izkazujejo tudi veliko natančnost, saj potrebujemo le malo končnih elementov za zelo natančne rezultate. Nekoliko drugačen je tudi pristop k integraciji v končnih elementih, saj za numerično integracijo v osnovni verziji nimamo analitičnih dokazov o minimalni potrebni stopnji integracije, zato jo določimo s poskušanjem; v izboljšani verziji pa numerična integracija po končnem elementu sploh ni potrebna. Pokazalo se je, da Newmarkova časovna integracijska shema ni primerna za integracijo pri deformacij-skih elementih. Veliko boljša izbira je t.i. midpoint časovna integracijska shema, ker je bolj stabilna. Žal tudi ta shema v kombinaciji z elementi z interpolacijo deformacij ni brezpogojno stabilna. Vse kaže, da bi bilo za brezpogojno stabilnost potrebno izpeljati posebej prilagojeno časovno integracijsko shemo; ta trditev je nekoliko izven konteksta tega poglavja, zato jo bomo podrobneje utemeljili v naslednjem poglavju. Lahko pa trdimo, da problem brezpogojno stabilnega integratorja za dinamiko deformacijskih elementov še ni rešen. Ker v osnovni obliki formulacije nastopajo izrazi z vgnezdenimi integrali, ki drastično upočasnijo računski postopek in naredijo formulacijo nekonkurenčno, smo razvili numerični postopek, ki se temu problemu izogne. Novi postopek temelji na numerični interpolaciji integrandov vgnezdenih funkcij in prenosu integracije iz osnovnih izrazov na interpolacijske funkcije. Prednosti formulacije z interpolacijo deformacij do določene mere zasenči zapletenost teh elementov in zaenkrat še nerešen problem ustrezne časovne integracije. 30 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. ___ ___ ___ v 4 DINAMIKA LINEARNO ELASTIČNIH v ___ ___ ___ ___ LINIJSKIH KONČNIH ELEMENTOV Z ______ ___ ___ v v INTERPOLACIJO KINEMATICNIH KOLIČIN 4.1 Uvod V tem poglavju se ukvarjamo s časovnimi integratorji, pri tem pa se omejimo na končne elemente, ki za osnovne neznanke privzamejo kinematične količine: pomik u in zasuk p. Razvili bomo dve varianti novega časovnega integratorja za analizo ravninskih nosilcev in ju primerjali z drugimi, standardnimi, kot sta Newmarkov (1959) ali pa HHT-a (1977) ter energijsko konzervativnimi, kot so integrator Standerja in Steina (1996), integrator Sima, Tarnowa in Doblareja (1995) in integrator Ibrahimbegovića in Mamourija (1999). Pri navajanju že razvitih energijsko konzervativnih integratorjev se bomo v veliki meri naslonili na delo Menica in Crisfielda (2002), ki sta uspela večino energijsko konzervativnih integratorjev zapisati na enoten način. Ta zapis smo tu primorani v prilagojeni obliki do določene mere reproducirati, da mu lahko dodamo novo razvita integratorja. Ponovili bomo tudi njun dokaz o ohranitvi gibalne in vrtilne količine, saj se da ta dokaz direktno aplicirati na novo razvita integratorja. Poleg tega bomo v tem poglavju razvili in predstavili tudi novo shemo za umetno dušenje višjih oblik nihanja. Temeljna novost je adaptivno vklapljanje in izklapljanje dušenja v odvisnosti od prisotnosti visokih frekvenc z računsko nezahtevnim algoritmom. V primerjavi z drugimi tu predstavljena shema dušenja uduši občutno manjši del osnovnih nihajnih oblik. 4.2 Klasičen pristop Princip virtualnega dela. Osnovne enačbe zapišemo s pomočjo principa virtualnega dela (PVD): L / [N fe + Mon — (q — ücp) öu — (niz — Jp'f) S e(u, N(s(u, (p)), M(n( PVD. (4.2) Galerkinova metoda. Neznanke u,

Pi(s) -, , (4.3) ip ^—' ®i u v—v, Uj = > Pi(S) -, f i=1 rei raitP r cxT i du v^ Tw \ oU,- r = / * i\S) r^ • (4-4) dip ^—' o /»(s) rj = > r^(s)(roi + Uj) = > Pi(s) ,, T. . (4.5) ^—' ^—' ^—' jo + Si sin r', (4.6) Ko v princip virtualnega dela (4.1) vstavimo variacije deformacij (4.6), dobimo: L / [NA öu' + N W A r' öip + M öpJ — (q — üc^) öu — (niz — Jp&) öip\ ds — S $U* = 0. o (4.7) V dobljeno enačbo (4.7) vstavimo interpolacijske nastavke (4.4), kar da r 5^7=i NT AT P'-öXJ j + (NrWrArr' )Pjö$>j + M P'- ö&j o — (q — ücpj PjöXJj — (niz — Jp0) Pjö&j] ds — S ÖXJ* = 0. V skladu z osnovnim izrekom variacijskega računa morajo biti členi pri variacijah enaki nič in tako dobimo sistem 3 nnp zveznih diferencialnih enačb 2. reda po času: g = NAP' — (q — ücp) Pj \ ds = 0, j = 1,..., nnp o L r -, / M P'- + N W A r' Pj — (niz — Jpfi) Pj ds = 0, j = 1,..., n\tp o = 0 (4.9) s pripadajočimi naravnimi in/ali bistvenimi robnimi pogoji, ki so zapisani v enačbah (3.5) in (3.6) in z začetnimi pogoji iz enačb (3.7). 1Y primeru ukrivljenih nosilcev ta trditev drži le, če je interpolacija pomikov sposobna točno opisati začetno obliko ukrivljenega nosilca. Enačba (4.5) velja le za ravne nosilce. 32 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Edini neznani funkciji, ki nastopata v enačbi (4.9), sta pomik in zasuk, ki ju interpoliramo po nastavkih (4.3). Neznane ostanejo še diskretne vrednosti pomikov in zasukov v interpolacijskih točkah kot diskretne neznanke problema. Enačbe so semi-diskretne, saj v njih nastopajo tudi časovni odvodi. Enačbe (4.9) se da zapisati kot vsoto treh delov, ki pripadajo notranjim silam gnot, vztrajnostnim silam gvzt in zunanjim silam gzun : g = gnot + gvzt + gzun = 0. (4.10) Tu so posamezni členi definirani tako: gnot NAP' ds, j = 1,..., nitp ( M P'- + N W A r' Pj j ds, j = 1,..., n;tp o (4.11) gvzt ücpPj ds, j = 1,..., nitp gzun J pip P j ds, j = 1,..., nitp J o L / ^Pj ds, j = 1,..., riitp J o L / mzPj ds, j = 1,..., n;tp o (4.12) (4.13) Časovna diskretizacija. Klasična pristopa k diskretizaciji časovnih odvodov v enačbi (4.9) sta uporaba diferenčne metode (npr. Graham (2001)) ali metode uteženih residualov (npr. Zienkiewicz, Taylor in Zhu (2005)). Najbolj razširjena metoda oz. časovni integrator je t.i. družina Newmarkovih integratorjev (Newmark (1959)), ki temelji na diskretizacijskih nastavkih: Un+l = Ura + Un At + At2 [(g — ßt) Üra + /3tÜra-|-i] ^ra+1 = ^Pn + ftn^-t + ^2 [Q — /?t) (pn + ßt&n+l] (4.14) (4.15) Ün+l = Üra + [(1 — 7t) At Üra + 7tÜra-|-i] ftn+l = ftn^ [(1 — 7t) At ipn + 7t n_|_i] Z različnimi vrednostmi parametrov ßt in jt dobimo različne algoritme oz. člane družine. Nekaj najbolj uporabljanih kombinacij je prikazanih v preglednici 4.1. Preglednica 4.1: Algoritmi iz Newmarkove družine integratorjev. Table 4.1: Algorithms of the Newmark family of integrators. ime ßt It komentar Trapezno pravilo Linearni pospeški Fox-Goodwin Central difference 0 4 2 1 1 6 2 1 1 12 2 0 1 2 implicitna, brezpogojno stabilna implicitna, brezpogojno stabilna implicitna, brezpogojno stabilna eksplicitna, pogojno stabilna Newmarkovo družino integratorjev lahko razširimo na t.i. družino a integratorjev tako, da določene člene v enačbah (4.9) vrednotimo pri različnih časih znotraj intervala od tn do Ln+i. Če vpeljemo oznaki o _ L _ _ Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 33 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. anot in avzt za skalami vrednosti med 0 in 1 in definiramo časa tanot = tn + anot At ter tavzt = tn+ avzt At, lahko enačbe (4.9) zapišemo (Chung in Hulbert (1993)) üavztcp) Pj\ ds = 0, j = 1,...,nitp {mz,anot ~ Jp(pavzt) Pj = ^' i = 1' • • • ' nitp (4.16) Z indeksoma (•)Q, = (•) (t«not) in (•)Q, = (•) (t«vzt) označujemo vrednosti pri časih tanot in t«vzt, ki se računajo po aproksimacijskih pravilih: (•) saJ imamo na voljo tri možnosti (Graham (2001)): gnot,anot v C^notJ gnot,ro ~r ttnot gnot,ro+l? v*"-!"/ gnot,anot = gnot((l ~ «not) Sn + anot Ln+l) , (4.19) gnot,anot = gnot (e ((1 ~ «not) Ura + anot Ura+i)) . (4.20) Varianta (4.18) je osnovna interpretacija v smislu a metod, varianta (4.19) je energijska (Energy-Momentum), varianta (4.20) pa simplektična (Symplectic-Momentum). Sredinska oz. ‘midpoint’ pravila. Če v vseh treh variantah a-metode za nelinearne probleme privza-memo sredinsko shemo, imamo anot = 2 m avzt = 2- ^° vPeiJemo oznako m za količine, ki so izvrednotene pri sredinskem času po enačbi (4.17), lahko enačbe (4.18)-(4.20) zapišemo: 1 gnot,m T lgnot,n ~r gnot,ro+l/ ; v^'^*-) 2 gnot,m = gnot (^m) , (4.22) gnot,m = gnot (e (um)) (4.23) g = NTAT P'- (q«n t / T T T t\ Manot P i + (N W ^ Y)n Pj - U in dobimo midpoint predstavnike shem (4.18)–(4.20). 34 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Z energijsko interpretacijo (4.22) lahko g iz enačbe (4.16) zapišemo tako: g = N^ A^ P'- — (q^ — ümCp) Pj\ ds = 0, j = 1,...,n\tp o / Mm P j + N^j W A^j r'm Pj — (niz,m — J P (fim) Pj\ ds = 0, j = 1,..., mtp O (4.24) Z upoštevanjem trapeznega pravila (preglednica 4.1) postanejo nastavki (4.14) in (4.15) ura+i = un + uraAt + 4f- üm, 2 (4-25) Un+1 = Ura + At ÜTO, (4.26) ') ds. o Ko vpeljemo interpolacijo inkrementov pomikov in zasuka (enačbe (4.3)) = Ei=l Pi(s) (4.50) (4.51) Au Aip Au' Aip' zn ms AUi A$i AUi A$i (4.52) dobimo Aildef = ^itp r -i T EAUj f A$i o B* NW A„ r'mPi + MmP ds, AmNmP( 0 B*Nm^ m m t kjer (•) označuje skalarno množenje in določimo prispevek notranjih sil k residualu v j-tem vozlišču: i] (4.53) &j,not r -^o- • i L r &7,not / 71 /T == M n / y j, not O A N P' B* N^WTA^ rfmPj + MmP *] ds. yj',not q 777, m m j ,hl j ¦*¦' -'*-' (4.54) Z oznako gj,not smo označili del residuala notranjih sil v j-tem vozlišču zaradi sil, z gj>not pa del zaradi momentov. Enačbo (4.54) lahko zapišemo z vektorskim produktom namesto z antisimetrično matriko W: &j,not N i L r g jnot / ]\/f ' —— 9'j,not o A N P' L i/ «S- 0 ~B* r'm X (AmNm) Pj + MmP Oba izraza, enačba (4.54) in enačba (4.55), sta matematično enakovredna. (4.55) Algoritem Sima, Tarnowa in Doblareja. Posebnost algoritma Sima, Tarnowa in Doblareja je, da pred interpoliranjem inkrementov pomikov in zasuka v enacˇbi (4.49) vpelje skaliranje inkrementa zasuka: Aips = tan A

tra+i izraža z vrednostjo njenega odvoda pri vmesnem času: ^n-\-i A (•) = / (•) (t) dt = (•)«(„) At. (4.66) Gornji rezultat moramo razumeti tako, da za nek določen a(n) G [tn, tra+i] velja, da je prirastek na časovnem intervalu enak odvodu funkcije (nekje) na tem intervalu, pomnožen z velikostjo časovnega intervala. a(n) je neznanka, ki se spreminja s časovnim korakom - je funkcija časovnega koraka n. Te neznanke se znebimo, če predpostavimo, daje a(n) = tn + -L¦ = tm: A (•) ~ (•) (tm)At. (4.67) To pravilo uporabimo pri integraciji enačb (4.65): ^n-\-l ^n-\-l e dt = I IjTl *^Ti ^n-\-i ^n-\-l r -i e dt = Ar' + A u' dt, L>T v TI ^n-\-i ^n-\-i / k dt = I (p' dt h T *->T . . (4.68) ^n+1 T"n-\-l in dobimo As« A r +A u\(t m )At, (4.69) An ~ (*m) ~ Vm = 2 { u'(tm)At = Au', P'(tm) =^r => P'(tm)At = Au?', (4.72) V(tm) = at => V(tm)At = Alo in enačbo (4.70) lahko zapišemo kot Ae = W A (tm) v'rn ^-^P + A (tm) All', (4.73) An = Aip . Diskretizirati moramo le še rotacijsko matriko, kar lahko storimo na dva različna načina, zato ločimo dve varianti šibkih integratorjev. 40 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Varianta z midpoint rotacijsko matriko. V duhu aproksimacijskega nastavka (4.71) zasuke v rotacijski matriki nadomestimo s povprečnimi: * * / r COS(fim — sin') ds. (4.77) o Ko izraze iz enačbe (4.75) vstavimo v (4.40), dobimo L 0 V tej enačbi interpoliramo inkremente pomikov in zasukov po nastavkih (4.61) ter dobimo n itp T ATT v^ AUi f A(tfim) NmP/ AÜHef = > a ;r ' TV-I-T-I-I7-T * T/ / ti i/ • (4.78) ć-^ A$i o ^"m** -^ vPm) rm^i + MmP,\ Ker velja 1 A(tL> ) =-------tAm, (4.79) cos tt1 v enačbi (4.78) zaradi konformnosti z ostalimi integratorji nadomestimo rotacijsko matriko A(ipm) s tem izrazom. Ko vpeljemo še oznako 1 A =-------r, (4.80) cos ^ 2 lahko enačbo (4.78) zapišemo v obliki P T ATT 1 L A * A TVT t/ AiJi / A Am\\mPi ^^i Jo A* N^WTA^ r'mPi + MmP[ p ATT lj EAUj A$i o ATTi ^tmm1 v ,. n1, Alldef = ' TT T • O+.ol) Iz nje določimo prispevek notranjih sil k residualu v j-tem vozlišču: K7not A* ATONTOP' gj,not = M = /l* tvtTtst-TaT / d i j\/T D/ "S" (4-82) Wg j not L ^4* ATONTOP' 5j,not o ^* ^mW ^ Če uporabimo vektorski produkt namesto matrike W, dobimo &?',not 9j,not o —^4* rTO x (ATONTO) Pj + MmPi M = A* I „ (K TVT N D I A/f D/ "S- (4-83) 9 j,not o — -^ rm X (•^¦m^"mj -* J + ^m,Pj Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 41 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Varianta s trapezno rotacijsko matriko. Rotacijsko matriko A(tm) lahko diskretiziramo tudi z nastavkom (4.71) oziroma (4.45), kakor daje A samostojna neznanka: 1 A(tm) ~ — (Ara + Ara+i) = Am. (4.84) 2 V tem primeru enačba (4.73) postane Ae = W A^j r'm Aip + A^ Au', (4.85) An = Aip . To je zopet pravilo, ki narekuje račun deformacij z inkrementi in je prikazano v enačbi (4.76). Ko izraze iz enačbe (4.85) vstavimo v (4.40), dobimo L Aildef = / (Nto W Am r'm Aip + NTO Am Au' + Mm AtL>') ds. (4.86) o Po nastavkih (4.61) interpoliramo inkremente pomikov in zasukov ter dobimo P T ATT 1 L r A TVT Ti/ 1 EAUj / Am \\m^i A$i o ^roW^A^ r'm,Pi + ^m,P( ltp ATT L ATT V^ AUj / ^^l; Allfjef = / ' TT T • (4.87) Iz te enačbe določimo prispevek notranjih sil k residualu v j-tem vozlišču: r N~ 1 L r a "XT d/ g7-not / J\m\Smr- gj>ot = m = TTT / d/ ds- (4M) Ce uporabimo vektorski produkt namesto matrike W, dobimo &7,not / 9j,not Jo —rm X (AmNm) Pj + MmP: wo- a i\t p' Ki not = M = I I * Â T-i/ «S- (4.89) 9j,not o —rm x (AmNm)ij- + MmF- 4.3.2 Kinetična energija Razlika kinetičnih energij je 1 fL t-T . -T . t 2 2 ^ Allkin = — (ura+1c/3ura+i — UyjCpUyj + EJ(pn+i — EJ(pn ds. (4.90) Razlike kvadratov v enačbi (4.90) razvijemo podobno kot v enačbi (4.36): 1 /Vi2 ,\2\ Un + l+Un ( • • • A • — I 1] -, --- 1] I = ------—-------- W,n i 1 --- U ZÄH = XWrnLXt. . At . f * ' (4.95) /Ó Vn + l — Pn A

rmj x Ap (un+i — un) Fj as — J p [fin+i — î^itr. T T 1LP r L L 1 AL = y Wm,j x / At gj,vzt ds — At 9j,vzt ds . (4.125) •_i o o - 0 o Ngj,vzt nadomestimo z wgj)Zun -N gj>ot in M9j,vzt z M9j,zun ~M g j,not ter dobimo 0 AL = 5^7=1 r"i,i x ^ ( "I™ Pj ~ A* ATONTOP' j ds— o — / At I B* r'm x (ATONTO) Pj + C* MmP'- \ ds . o (4.126) Pri členu C*MmP'- upoštevamo zvezo (4.119), pri členu B* r'm x (ATONTO) Pj pa zvezo (4.115), izpostavimo At ter člene preuredimo in dobimo n itp \v v 7 =1 1m,J ^ L L i + rmj x I A* AmNmP- ds — B* r'm x (ATONTO) ds 0 0 AL = At 5^7=i r"i,i x "I™ L,- ds + (4.127) V zadnjem koraku izpeljave obtežne člene zanemarimo, saj zahtevamo ohranjanje v odsotnosti obtežbe, interpoliamo r'm in diskretne vrednosti krajevnega vektorja postavimo pred integral: izp r i AL = At y rm,j x A* AmNmPJ ds — rmj x / I3*ATONTO P- ds . (4.128) 7_i o o 3 = Enačbo preuredimo in dobimo: AL = At y rm,j x / (-A* ~~ B*) ATONTOP' ds. (4.129) 7—1 ^ Potrebni in zadostni pogoj za ohranjanje vrtilne količine je, da na (nepodprti) sistem ne deluje obtežba in da velja ^* _ ß* (4.130) Temu pogoju zadoščajo vsi algoritmi v preglednici 4.2 razen algoritma Standerja in Steina. Dokaz je povzet po delu Jeleniča in Crisfielda (2002). Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 47 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 4.4 Algoritmi Algoritmi posameznih integratorjev zahtevajo upoštevanje določenih specifičnih pravil v računu. V tem poglavju pravila še enkrat predstavimo in poudarimo razlike med integratorji. Za integratorje vpeljemo naslednje oznake: Newmarkov integrator (New), integrator Hilberja, Hughesa in Taylorja (HHT), integrator Standerja in Steina (StaSt), integrator Sima, Tarnowa in Doblareja (STD), integrator s šibkimi kinematičnimi vezmi in povprečno rotacijsko matriko (Wpov), integrator s šibkimi kine-matičnimi vezmi in sredinsko rotacijsko matriko (Wmid). Kot zadnjega navedemo integrator Ibrahimbe-goviča in Mamourija (IbMa), ki ga do sedaj nismo izpostavljali. Integrator IbMa je namreč v skoraj vseh pogledih identičen integratorju STD z drobno, a ključno razliko v numerični implementaciji, ki jo bomo prikazali v naslednjih vrsticah. 1. Začetek algoritma t = 0: Newmark, HHT: Iz znanih začetnih pogojev (3.7) se izračuna Uo, 3>o iz enačbe (4.29). StaSt, STD, Wpov, Wmid, IbMa: Znane so vrednosti začetnih pogojev (3.7). 2. Račun pri času tn+i\ Znane so vrednosti neznank pri času tn: Newmark, HHT: Ura, 3>ra, Ura, 4>ra, Ura, 3>ra. StaSt, STD, Wpov, Wmid, IbMa: Ura, 3>ra, Ura, 4>ra. > Vstop v Newtonovo iteracijsko zanko. a. Inicializacija: Ura+i = Ura, *ra+i = *„. b. Račun hitrosti in pospeškov: Newmark, HHT: tjn+i, $ra-|-i po enačbah (4.14). Un-|_i, $n+i po enačbah (4.15). StaSt, Wpov, Wmid: Un-|-i, $n+i po enačbah (4.97). UTO, 3>TO iz enačb (4.98). STD, IbMa: Ura-|-i po enačbah (4.97). ^n+i po enačbah (4.103). UTO, 3>TO po enačbah (4.95). 48 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. • Vstop v integracijsko zanko z reducirano integracijo. c. Interpolacija: Newmark, HHT, StaSt, Wpov, IbMa: un+i) Vn+ii un> Vn> um, Au, Aip, po enačbah (4.3) in (4.52). Wmid: ura+i) ^ri+i) un, <ßn> um, fmi Au, AtL>, P° enačbah (4.3) in (4.52). STD : un, um, Au, po enačbah (4.3) in (4.52). Aips^n+i, po enačbi (4.61) - interpolira se skalirane zasuke! Aipn+l, iz enačbe (4.56): AtL>ra+1 = 2 arctan^n+i (v integracijskih točkah). u«+i> iz inkrementov: ura+i = un + Au, ipn+1 =(fin + Aip. d. Račun deformacij: Newmark, StaSt, STD, IbMa: Ln+i, Kn+i po enačbah (2.24). HHT (račun pri tn + at At): u(tn + at At), cp(tn + at At) po enačbah (4.17). s(tn + at At), n(tn + at At) po enačbah (2.24). Wmid: Ln+i, Hn+i po enačbah (4.75) in (4.76). Wpov: Ln+i, Hn+i po enačbah (4.85) in (4.76). e. Račun notranjih (konstitucijskih) sil: Newmark, HHT: Nn-|_i, Mn-\-i po enačbah (2.30) in (2.28). StaSt, STD, Wpov, Wmid, IbMa: Lm in nm po enačbah (4.38). NTO, Mm po enačbah (4.39) in (2.30). f. Zapis ‘togostnega’ dela enačb gnoi> Togostnim silam pripadajoči del enačb (4.111): A* AmNP' fL A* ATONTOP' Jo B* N^WTA^ r'mPj + C* MmP' gnot = T T T i ^S' J = 1) • • • )nitp- B* N„W A^j rmPj + C* MmP- Tangentna matrika (odvod vseh enačb gnot po vseh neznankah \Jn+i in ^n+i): ^gnot - določimo s programom AceGen3. 2Interpolacija in racˇun zasukov je edina razlika med integratorjema STD in IbMa. 3Program AceGen (Korelc (1997)) je program za avtomatsko odvajanje in generacijo programske kode. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 49 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. g. Račun količin z vrednostmi pri tn-^\: Deformacijska energija ü^ef po enačbi (2.16). • Izstop integracijske zanke z reducirano integracijo. • Vstop v integracijsko zanko s polno integracijo. c. Interpolacija: Newmark, HHT: üri+i, <ßn+i P° interpolacijskem nastavku (4.3) na vrednostih Ura-|-i in $n_|_i. StaSt, STD, Wpov, Wmid, IbMa: üm, (pm po interpolacijskem nastavku (4.3) na vrednostih UTO in $TO. f. Zapis ‘vztrajnostnega’ delagvzt: Vztrajnostnim silam pripadajoč del enačb (4.111): 1 fL Ap (ura+i — ura) Pj gvzt = -TT „ T / . . , ds. ^¦t o tij [ Ponavljaj dokler \\5Y\\ > predpisana natančnost. Edina razlika med integratorjema STD in IbMa je pri interpolaciji oz. računu zasukov. V algoritmu IbMa se namreč interpolira zasuke oz. inkremente zasukov, medtem ko se v algoritmu STD interpolira skalirane inkremente zasukov. Posledica je, da integrator IbMa ne uspe enako dobro ohranjati mehanske energije kot ostali energijski integratorji. 4.5 Numerični primer Leteči špaget. Računski primer ti. letečega špageta, ki sta ga predlagala Simo in Vu-Quoc (1986) in kasneje reševala tudi Hsiao in Jang (1991), je idealen za preverjanje konzervativnih lastnosti integrator-jev, saj je sistem nepodprt in bi moral v odsotnosti obtežbe ohranjati vse količine gibanja: mehansko energijo ter gibalno in vrtilno količino. Primer obravnava zelo podajen nosilec, ki se močno upogibno deformira in v tem pogledu spominja na let kuhanega špageta, od koder izvira tudi ime tega primera. 50 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 4.1: Geometrijski podatki o primeru, levo; odvisnost obtezˇbe od cˇasa, desno. Figure 4.1: Geometrical data, left; loads vs time, right. V začetnem stanju nosilec miruje v poševni legi, kot je razvidno s slike 4.1, levo. Na spodnjem robuje obremenjen s točkovno silo in momentom, ki na nosilec delujeta 2.5 sekund, zatem pa prenehata delovati, kot prikazujeta grafa na desni strani slike 4.1. Materialni podatki so: E A = GAS = 10000 N, E J = 100 Nm , Ap = lkgm- , J p = lOkgm. Let špageta, ki je prikazan na sliki 4.2, je izračunan z integratorjem Standerja in Steina, z mrežo 10-ih končnih elementov s kvadratično interpolacijo in s časovnim korakom At = 0.01 s. Zaporedje desetih deformiranih leg je prikazano na vsakih 1.25 s. Analizo z različnimi integratorji izvršimo do računskega časa 100 s ter opazujemo ‘kvaliteto’ ohranjanja količin gibanja. Kvaliteto ohranjanja ocenimo z razliko med maksimalno in minimalno vrednostjo količine na intervalu 2.5 do 100 s. Nihanje mehanske energije označimo z Ali, nihanje komponent gibalne količine v smereh X in Y z Apx in Apy, nihanje vrtilne količine pa z ALZ. New-markov integrator uporablja trapezno pravilo. Kot kriterij konvergence Newton-Raphsonovega postopka je upoštevana absolutna norma prirastkov 1 x 10-12. Da bi tudi z računskim primerom potrdili veljavnost dokazov za različne velikosti časovnih korakov, smo vsako analizo opravili z dvema različnima velikostima časovnih korakov. Večji od uporabljenih časovnih korakov je dolg 0.1 s, kar pomeni, da za analizo potrebujemo 1000 časovnih korakov; manjši časovni korak je 0.01 s in z njim potrebujemo 10000 časovnih korakov. Neodvisnost od stopnje interpolacije po kraju preverimo z uporabo dveh različnih stopenj integracije: z linearno interpolacijo, ki jo označimo z ‘II’ in z interpolacijo četrte stopnje, ki jo označimo z ‘14’. Rezultati so zbrani v preglednici 4.3. Ob analizi rezultatov v preglednici 4.3 ugotovimo, da mehansko energijo z natančnostjo, kije primerljiva ostanku (residualu), ohranjajo vsi integratorji razen Newmarkovega s trapeznim pravilom (da New-markova shema s trapeznim pravilom ne ohranja mehanske energije, sta pokazala že Crisfield in Shi (1996)) in integrator Ibrahimbegovića in Mamourija. Slabo ohranjanje energije integratorja IbMa je posledica direktne interpolacije zasukov v kombinaciji z enačbami, ki pripadajo interpolaciji skaliranih zasukov. Integratorji, ki energijo ohranjajo, imajo približno enako kvaliteto ohranjanja energije ne glede na velikost časovnega koraka in stopnjo interpolacije. Newmarkov časovni integrator s trapeznim pravilom je energijo praviloma ohranjal vsaj za dva reda velikosti bolje, če smo uporabili (za red velikosti) manjši časovni korak. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 51 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 4.2: Let sˇpageta. Figure 4.2: Free flight of spaghetti. Preglednica 4.3: Špaget: rezultati računa. Table 4.3: Spaghetti: results. StaSt Alg. ?t [s] Int. ?? [J] ?px [kg m/s] ?py [kg m/s] ?L z [Nms] 0.1 0.01 II 14 II 14 1.2 x 10 2.5 x 10 4.1 x 10-6 7.5 x 10-7 1 2 x 10-11 6 9 x 10-12 4 7 x 10-10 2 1x10-10 2 9 y 10-13 3 8 y 10- 13 1 1 v 10-11 75x10-12 1.3 x 10 1.4 x 10-7 7.8 x 10 3.8 x 10 2 2x10-11 5 8 y 10-12 9 1 x 10-10 4 7x10-10 3 9 y 10- 13 3 2 y 10-13 2 1x10-11 75x10-12 4.1 x 10-1 3.6 x 10 4.0 x 10-3 5.5 x 10-3 3.1 x 10 3.3 x 10-8 4.5 x 10 3.3 x 10 STD 0.1 0.01 II 14 II 14 Wpov 0.1 0.01 II 14 II 14 1.2 x 10-7 2.2 x 10 4.0 x 10-6 7.5 x 10-7 14x10-11 11x10-11 41x10-10 1 8 y 10- 10 5 7 y 10-13 2 5 y 10-13 1 3 y 10-11 9 2 x 10-12 2.6 x 10-8 4.5 x 10 1.1 x 10-6 2.0 x 10-7 Wmid 0.1 0.01 II 14 II 14 1.2 x 10 2.1 x 10-7 4.1 x 10-6 7.5 x 10 12x10-11 6 2 x 10-12 3 5 y 10-10 3 1 y 10- 10 66 x 10-13 2 8 y 10- 13 7 2x10-12 9 5 y 10-12 1.4 x 100 1.3 x 10 1.0 x 10 2.6 x 10-2 2 0x10-11 12x10-11 1.2 x 10 2 9 x 10-10 3 0 y 10-13 2 7x10-13 14x10-11 8 4 x 10-12 2.7 x 10 4.4 x 10-8 1.1 x 10-6 2.0 x 10 2.9 x 10-8 3.3 x 10 4.8 x 10 2.5 x 10-9 IbMa 0.1 0.01 II 14 II 14 8.5 x 10 1.1 x 10 6.8 x 10-2 1.8 x 10 2 5 y 10- 11 2 0x10-11 1.0 x 10-9 8 5 y 10-10 9 4 x 10-13 3 7 y 10- 13 2 6 x 10-11 1 3 y 10- 11 1.7 x 10 2.2 x 10 1.4 x 10-2 2.2 x 10 New 0.1 0.01 II 14 II 14 52 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Gibalno količino po pričakovanju ohranjajo vsi integratorji, ne glede na velikost časovnega koraka in stopnjo interpolacije. Vrtilno količino, zopet ne glede na velikost časovnega koraka in stopnje interpolacije, ohranjajo vsi integratorji razen integratorja Standerja in Steina, kar smo pokazali že v teoretičnem delu, in Newmarkovega integratorja s trapeznim pravilom (za slednjega so to pokazali že Simo, Tarnow in Wong (1992)). 4.6 Sklepi Pokazali smo, da lahko integratorje izpeljemo iz dveh vsebinsko različnih izhodišč. Klasični pristop zapiše ravnotežne enačbe oz. princip virtualnega dela, v njih upošteva kinematične in konstitucijske enačbe ter z diskretizacijo po kraju in času izpelje metodo. Vse enačbe torej držijo v svoji osnovni obliki, zdi pa se, da na tak način ne moremo doseči ohranjanja vseh količin gibanja, kar je posledica oblike enačb samih. Po drugi strani pa ‘energijski’ pristop v osnovi izhaja iz zahteve po ohranitvi mehanske energije. Osnovna oz. izhodiščna enačba je torej bistveno drugačna. V nadaljevanju moramo z upoštevanjem ostalih enačb problema izhodiščno enačbo predelati v tako obliko, da iz enačb oz. residuala izpostavimo inkremente neznank problema. Postopkov za dosego tega cilja je več. Ob izpeljavi v vsakem postopku vpeljemo določene aproksimacije, s katerimi si pomagamo do končnega cilja. Ne glede na izbrano pot so končne enačbe de facto izpeljane iz izhodiščne enačbe. Kot smo opazili pri izpeljavi energijskih integratorjev, gre pri šibkih algoritmih in pri algoritmu Standerja in Steina za enako časovno diskretizacijo kot pri Newmarku s trapeznim pravilom. Gre torej za konstantno interpolacijo pospeškov po času, linearno interpolacijo hitrosti po času in kvadratično interpolacijo pomikov po času. Interpolacija zasukov po času v primeru integratorja Sima, Tarnowa in Doblareja je nekoliko drugačna in ni primerljiva z interpolacijo po Newmarkovi shemi. Izredno zanimiva je ugotovitev, daje izpeljana oblika ravnotežnih enačb, zapisana v enačbi (4.111) in v pripadajoči preglednici 4.2, v skoraj vseh primerih - razen v primeru ‘šibkega integratorja s povprečno rotacijsko matriko’ - enaka osnovni obliki ravnotežnih enačb (4.24) le v limiti A? —> 0, ko členi A*, B* in C* limitirajo k 1. Za A? ^ 0 se klasično ravnotežje razlikuje od izpeljane oblike ravnotežja, ki ga zato imenujemo ‘energijsko ravnotežje’. Če klasično ravnotežje ustreza splošni predstavi ravnotežja, kjer je vsota vseh sil in vsota vseh momentov enaka nič, se lahko upravičeno vprašamo, kakšno je ‘energijsko’ ravnotežje. Energijsko ravnotežje je tako ravnotežje, ki zagotavlja ohranjanje mehanske energije in s tem potrebni pogoj za stabilnost časovnega reševanja. Poleg tega v limiti A? —> 0 zagotavlja veljavnost klasičnega ravnotežja, za končno velik A? pa ne. Zato vsota vseh sil in vsota vseh momentov v splošnem ne bosta enaka nič, kar je v gradbeništvu malodane krivoverstvo, toda taka je cena stabilnosti reševanja pri omenjenih integratorjih. Ohranjanje mehanske energije smo dokazali za vsako velikost časovnega koraka. To je zelo pomembno, saj nam energijsko ravnotežje zagotavlja ohranjanje mehanske energije popolnoma neodvisno od velikosti časovnega koraka. To je v nasprotju s klasičnimi integratorji, kjer ponekod obstaja ‘ljudsko’ prepričanje, da manjši časovni korak izboljša stabilnost in ohranjanje mehanske energije. V našem računskem primeru se je pokazalo, da se ohranjanje mehanske energije z manjšanjem časovnega koraka izboljšuje, toda ta računski primer je netog. Za toge sisteme je znano (Stander in Stein (1996)), da manjšanje časovnega koraka lahko pri nekonzervativnih integratorjih pripelje do še bolj zgodnjega pojava časovne nestabilnosti. Ohranjanje energije smo dokazali neodvisno od stopnje in vrste interpolacije Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 53 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. po času, kar smo preverili tudi na računskem primeru. V tem poglavju smo teoretično in računsko pokazali tudi, da vsi obravnavani integratorji ohranjajo gibalno količino. Izpeljali smo tudi pogoj za ohranjanje vrtilne količine in računski primer je izpeljani pogoj potrdil. Velikost časovnega koraka je brezpredmetna tudi pri ohranjanju gibalne in vrtilne količine. Pokazali smo tudi, da lahko izpeljane algoritme ločimo na dve skupini glede na obliko kinematičnih enačb: na algoritme s šibkimi kinematičnimi vezmi in algoritme s krepkimi kinematičnimi vezmi. Če v izpeljavi uporabimo krepke kinematične vezi, nimamo možnosti, da bi se energijsko in klasično ravnotežje ujemala. Nasprotno pa, če uporabimo šibke kinematične vezi in rotacijsko matriko diskretiziramo tako, da pri sredinskem času vzamemo njeno povprečno vrednost, se bosta obliki energijskega in klasičnega ravnotežja ujemali. Ne glede na to, kako obrnemo problem, se moramo za ohranjanje energije nečemu odpovedati. V primeru integratorjev s krepkimi kinematičnimi vezmi se odpovemo klasični obliki ravnotežja, ki zagotavlja ničelnost rezultant. Zadoščeno pa je energijskemu ravnotežju, ki pa je le malo drugačno od klasičnega. V primeru algoritmov s šibkimi kinematičnimi vezmi se odpovemo krepki obliki kinematičnih enačb. V analitičnem smislu sta krepka in šibka oblika kinematičnih enačb enakovredni, vendar z diskretizacijo enakost izničimo. Zato krepke Reissnerjeve kinematične enačbe niso izpolnjene. Za ohranjanje mehanske energije smo se torej odpovedali točnim krepkim kinematičnim enačbam, izpeljano energijsko ravnotežje pa je enako klasičnemu, kar je prednost. Integratorji z ohranjanjem količin nedvomno izboljšajo stabilnost v računskih postopkih, vendar zaradi oblik enačb problema od nas zahtevajo, da se nečemu odpovemo: bodisi ravnotežju v klasičnem smislu ali pa kinematičnim enačbam v krepki obliki. 4.7 Dušenje Uvod. Odziv konstrukcije na splošno obtežbo je sestavljen iz neskončnega števila njenih nihajnih oblik. Te imajo v odzivu različne stopnje participacije, ki praviloma z naraščajočimi frekvencami upadajo, z inženirskega stališča pa hitro postanejo zanemarljivo majhne, saj za zadovoljivo natančnost opisa odziva večinoma zadostuje že nekaj prvih lastnih oblik nihanja. Diskreten matematični model konstrukcije ima končno število prostostnih stopenj in enako število nihajnih oblik. Vsaki nihajni obliki pripadajo svoje frekvence in nihajni časi. Za vsako nihajno obliko, ki jo želimo ‘zadovoljivo’ upoštevati v odzivu, moramo pri numerični časovni integraciji predvideti približno deset časovnih korakov za en nihajni čas, da obliko zadovoljivo natančno zajamemo v odzivu (Bathe (1996)). Nihajni časi relativno hitro padajo, tako da smo pri integraciji vedno v položaju, da celo množico nihajnih oblik računamo s prevelikim časovnim korakom. Zaradi prevelikega časovnega koraka take nihajne oblike povzročajo numerične motnje, ki se kažejo kot šum, ki lahko postane moteč. Dodatno je mreža končnih elementov v krajevnem smislu praviloma pregroba, da bi zadovoljivo modelirala deformiranje v višjih nihajnih oblikah in to je dodaten vir numeričnih motenj v analizi. Te vire napak in motenj v problemu omilimo oz. odstranimo z umetnim, numerično povzročenim dušenjem. Osnovni namen takega dušenja je izločiti višje nihajne oblike iz odziva. Ključen problem je, kako doseči, da bo umetno dušenje vplivalo samo na nihajne oblike z visokimi frekvencami. V nadaljevanju besedila predstavimo enega od možnih pristopov, ki temelji na ideji, da dušenje deluje le, kadar in kjer je potrebno, sicer pa ne, ter s tem minimizira vpliv numeričnega dušenja na osnovne nihajne oblike. Sheme z numeričnim dušenjem. Poudariti moramo, da se tu ukvarjamo z ‘numeričnim’ dušenjem, ki 54 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. z realnim dušenjem materiala nima povezave. Dušimo nenatančno izračunane višje nihajne oblike, ki so posledica napak računske metode. Metod z numeričnim dušenjem je moč v literaturi zaslediti kar nekaj. V linearni dinamiki sta dušenje in stabilnost metode temeljito raziskana, kot poroča Fung (2003) v svojem preglednem članku, težišče raziskav pa se je premaknilo na izboljšanje natančnosti (reda) in-tegratorjev. V nelinearni dinamiki pa razvoj še vedno poteka in v uporabi je več vrst integratorjev z numeričnim dušenjem. Nekatere, kot so Huberta, Hughesa in Taylorja (1977) ali pa Chung - Hulbertova (1993) posplošena a-metoda, smo že omenili, v novejši literaturi pa se pojavljajo še druge, npr. metode Ibrahimbegoviča in Mamourija (2002), Armera in Petöcza (1999) in Armera in Romera (2003). Sheme, ki temeljijo na posplošeni a-metodi oz. ki so bile razvite prvenstveno za linearne sisteme, ne zagotavljajo disipacije energije, kar so nazorno pokazali Crisfield, Galvanetto in Jelenič (1997) in Erlicher, Bonaventura in Bursi (2002). Koncept in enačbe. Dušenje apliciramo le na deformacije na podoben način kot Armero in Petöcz (1999), Armero in Romero (2003) in Sansour, Wriggers in Sansour (2003): (4.131) kjer smo z e*m in n*m označili deformacije z dušenjem, z a^am faktor za skaliranje dušenja in z d zaenkrat nedefinirano funkcijo, ki ohranja predznak. Ko em in km zamenjamo z izrazi (glej enačbo (4.94)) i. Ae f^YTI —— Ai (4.132) (4.133) in At nesemo v funkcijo d, ki zato postane funkcija d, dobimo: m ^m ~r C^dam & v J > Razlika energije Ali v časovnem koraku je (enačbe (4.33), (4.37) in (4.102)) L An = / {s^ceAs + E J Km An + iimCpAu + EJ(pmA(ps) ds. (4.134) 0 V tej enačbi deformacije em in Km nadomestimo z em in nm iz enačbe (4.133) ter dobimo L An = I (cßLm) As + E J K,m Ak, + iimCpAu + EJ(pmA(p J ds. (4.135) o o Deformacije razvijemo in dobimo L L r -, An* = {^Eem) As + E J K,m Ak, + iimCpAu + EJ(pmA(p ds + o L r , \ T / \ T ~\ + / «dam fcE d(Ae)J Ae + adam (cEd(An)) An\ ds. (4.136) Prva vrstica te enačbe je enaka razliki mehanskih energij An, kakršno smo imeli že doslej, drugo vrstico pa razvijemo tako, da združimo ( cE d (AeTO) j Ae v funkcijo d* (Ae2) in I cE d (An) j An v d* (Ak2): L Au* = An + / [cKdamd* (Ae ) + «dam d* (Ak )] ds. (4.137) o Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 55 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Prispevek integralskega člena v tej enačbi disipira energijo za vsak cetani > 0, če je le d* taka funkcija, ki ohranja predznak argumenta. Algoritem. Kot smo že omenili, želimo dušenje aplicirati le tam, kjer je to potrebno in še to le takrat, ko je to pomembno. Algoritem temelji na opazovanju dogajanja v integracijskih točkah in vklapljanju ter izklapljanju dušenja v teh točkah. Integracijskih točk v konstrukciji pa je lahko relativno veliko, zato mora biti algoritem za zaznavo potrebe po dušenju karseda preprost. Predlagani algoritem temelji na opazovanju velikosti in predznaka sprememb (inkrementov) deformacij v zadnjih nhis korakih. Algoritem se začne brez dušenja. Ko oziroma če zaporedje predznakov opazovanih inkrementov deformacij v neki točki tvori popolno žagasto zaporedje pozitivnih (+1) in negativnih (—1) vrednosti, se dušenje v tej točki vklopi. Dodatno zahtevamo, da mora biti za vklop dušenja nihanje zadosti veliko (večje od predpisanega praga t), kar testiramo tako, da pogledamo največjo absolutno vrednost opazovanih inkrementov deformacij. V točkah, kjer je prisotno dušenje, opazujemo, ali je nihanje zadosti majhno, da bi dušenje lahko izklopili. Kriterij za izklop dušenja je enakost predznakov inkrementov deformacij v zadnjih n^ korakih. Za izklop zadostuje tudi, da so največje vrednosti inkrementov deformacij manjše od predpisanega praga r. Po izklopu dušenja se v vsakem koraku zopet preveri, ali je morda potrebno dušenje ponovno vklopiti. Značilen primer, ko se dušenje vklopi in izklopi, je shematično prikazan na sliki 4.3. Slika 4.3: Shematski prikaz primerov ob vklopu in izklopu dušenja. Figure 4.3: Schematic representations of cases when damping is engaged and disengaged. Prag vrednosti r, kije pomemben pri vklopu in izklopu dušenja, je tesno povezan z vrednostmi notranjih sil, saj je ključen za odločitev o tem, kaj predstavlja šum in kaj ne. Odločitev o tem, kako velike sile nam predstavljajo šum, mora sprejeti uporabnik sam pred analizo. Na podlagi izbrane velikosti sile, ki predstavlja šum, se dejansko vrednost praga dobi s kvocientom te vrednosti in pripadajoče togosti prereza. Število opazovanih korakov nhis, ki je teoretično lahko različno pri vklopu in izklopu, je neke vrste parameter občutljivosti; manjši ko je, hitreje se bo algoritem odločil, da gre za primer, ko je potrebno dušenje vklopiti ali izklopiti. Na sliki 4.3 je prikazan primer, ko je vrednost n^s = 4 enaka za vklop in izklop. Zgornja meja teoretično ni omejena, spodnja pa je dva. Praktični računi so pokazali, da so smiselne vrednosti od 4 do 10. Do sedaj smo od oblike funkcije dušenja zahtevali le, da ohranja predznak. Izbira je torej še povsem odprta. Praktično imamo na voljo le tri konceptualno različne izbire: (i) polinom lihe stopnje; (ii) linearno funkcijo in (iii) navzgor omejene funkcije , kakršna je funkcija arctan. Načelno se zdi najbolj primerna liha polinomska funkcija, saj bi količina dušenja z rastjo hitrosti naraščala s potenco. 56 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Računski primeri so pokazali, da so funkcije takega tipa preveč agresivne in/ali preveč udušijo odziv, ali pa ob vklopu dušenja povzročijo tak numerični šok za metodo, da nastopijo težave s konvergenco Newton-Raphsonovega postopka. Linearna funkcija se zdi primerna, ker je preprosta. Kljub temu se je na praktičnih primerih pokazalo, da je tretja vrsta funkcije dušenja bolj primerna. Navzgor omejena funkcija dušenja, kakršna je funkcija arctan, namreč omeji količino dušenja v posameznem koraku analize in s tem omogoča bolj gladek prehod med fazami z dušenjem in brez in s tem v največji možni meri ohranja stabilnost Newton-Raphsonovega postopka reševanja enačb. Funkcijo še bolj zgladimo, če vpeljemo dodaten parameter /?dam: €rn) — arctan----- ) . (4.138) A^dam 4.8 Numerični primer Bauchaujevo nihalo. Računski primer je poimenovan po njegovem prvem avtorju inje bil predstavljen v članku Bauchauja, Damilana in Therona (1995). Kasneje so ga preračunavali tudi drugi avtorji: Bauchau in Theron (1996), Ibrahimbegovič in Mamouri (1999), Ibrahimbegovič in sodel. (2000) in Ibrahimbe-govič in Mamouri (2002). Posebnost primera je, daje pripadajoči sistem vodilnih enačb izredno tog. Nihalo je sestavljeno iz podajnega vodoravnega nosilca, ki je členkasto vpet v dve togi palici v točkah A in E. Sredi podajnega nosilca se nahaja točkovna masa m = 0.5 kg. Točka B je opazovališče. Leva palica je navpična, desna pa pod kotom 7r/4 od vodoravnice v protiurnem smislu, kot je prikazano na sliki 4.4, levo. Togi palici sta brez mase (p = 0 kg/m3), njuno togost pa modeliramo z veliko vrednostjo elastičnega modula, ki jo po avtorjevih (Bauchau) navodilih dobimo tako, da osnovni elastični modul pomnožimo s faktorjem 10. Strižne napetosti v računu zanemarimo tako, da privzamemo veliko vrednost strižnega modula {G = 100 E). V začetnem stanju sistem miruje, zatem pa nanj prične delovati točkovna sila S, ki ima prijemališče na mestu mase. Sila S linearno narašča od vrednosti 0 v začetnem stanju do maksimalne vrednosti 2 N pri času 0.128 s. Zatem linearno pada do časa 0.256 s, ko preneha delovati. Od tega trenutka dalje je sistem prepuščen prostemu nihanju. Velikost sile S v odvisnosti od časa je prikazana na sliki 4.4, desno. Slika 4.4: Geometrija in obtezˇba Bauchaujevega nihala. Figure 4.4: Geometry and loading of the Bauchau swinging pendulum. V analizi uporabimo enako velikost časovnega koraka, kot ga priporoča literatura (At = 0.0005 s). Po- Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 57 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. dajni nosilec modeliramo z osmimi elastičnimi končnimi elementi, vsako togo palico pa s po štirimi končnimi elementi s kvadratično interpolacijo pomikov in zasuka. Računamo s šibkim časovnim inte-gratorjem s sredinsko rotacijsko matriko Wmid. Materialni in geometrijski podatki primera so: E = 73 x 10 N/m , L = 0.72 m, A = 0.005 x 0.001 m , 1 3 4 S J = — x 0.005 x 0.001 m , p = 2700 kg/m , ^nosilca — 0.00972 kg. 12 Na začetku obe palici nihata v nasprotni smeri urinega kazalca, zatem pa pri času t o 0.64 s pride do nenadne spremembe smeri nihanja in desna palica začne nihati v smeri urinega kazalca. Hkrati pride do hipne spremembe horizontalne komponente hitrosti gibanja mase m. Ker predstavlja točkovna masa bistveni del mase sistema, saj je približno 50-krat večja od mase nosilca, dogodek na sistem deluje kot nekakšen trk. V konstrukciji to povzroči nenadno in izrazito spremembo v osnem nihanju nosilca. Kljub temu pa sistem nadaljuje pretežno nihajoče gibanje okoli točk 0\ in O^- Slika 4.5: Bauchaujevo nihalo: vodoravni in navpični pomik v točki B. Figure 4.5: Bauchau swinging pendulum: horizontal and vertical displacements at point B. Primer analiziramo s tremi različnimi vrednostmi parametra o^am- (i) brez dušenja («dam = 0), (ii) z majhno stopnjo dušenja («dam = 0.01, /?dam = 100, nhis — 5) in (iii) s srednjo stopnjo dušenja (<^dam — 0.1, /?dam = 100, Ti^is = 5). Mejna vrednost šuma znaša 0.01 N za osne in strižne sile ter 0.00001 Nm za upogibne momente. Na sliki 4.5 so prikazani vodoravni in navpični pomiki v točki B. Razlike med analizami z različnimi stopnjami dušenja na tej sliki niso vidne, zato je na sliki 4.6 prikazana povečava. Na povečavi lahko vidimo, da rešitev z majhno stopnjo dušenja le malenkostno zgladi nedušeno rešitev. Rešitev s srednjim dušenjem predstavlja bistveno bolj gladek odziv. Odziv glavnih nihajnih oblik je navidez ohranjen, tako s stališča amplitud kot tudi nihajnih časov. Primerjava količine disipirane energije je prikazana na sliki 4.7. Očitno je, da se disipira le majhen delež celotne energije. Hkrati pa je razviden trend disipacije, ki se po začetnem skoku zvezno umiri in praktično preneha delovati. Informacijo o tem, kdaj in kje je dušenje delovalo, nam daje slika 4.8, kjer so področja z dušenjem potemnjena. Na abscisno os je projecirana celotna konstrukcija, na ordinato pa čas. Opazimo lahko, da se večina dušenja aktivira v podajnem nosilcu. V togih palicah nastopajo le osne deformacije, na katerih se tudi aktivira dušenje. Zanimivo je, da večina dušenja nasploh odpade na osne deformacije in da dušenja upogibnih deformacij skorajda ni. Nekoliko presenetljivo pa je dejstvo, da, čeprav se količina dušenja zmanjšuje, še zdaleč ne pride do tega, da bi dušenje popolnoma prenehalo delovati. Prav tako ne 58 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 4.6: Bauchaujevo nihalo: povecˇava vodoravnih in navpicˇnih pomikov v tocˇki B s slike 4.5. Figure 4.6: Bauchau swinging pendulum: magnification of the horizontal and vertical displacements at point B from Fig. 4.5. Slika 4.7: Bauchaujevo nihalo: disipacija energije z razlicˇnimi stopnjami dusˇenja. Figure 4.7: Bauchau swinging pendulum: energy dissipation with different levels of damping. kaže, da bi se dušenje kdaj v celoti izklopilo. Na podlagi disipacije energije (slika 4.7), ki se zelo hitro umiri, lahko sklepamo, da se dušenje ves čas vklaplja, vendar pa se njegova intenziteta zelo zmanjšuje. Bistveno bolj kot informacija o tem, kdaj in kje se je vklopilo dušenje, je pomemben vpliv dušenja na pomike in notranje sile. Slednje so prikazane na sliki 4.9. Notranje sile so prikazane v prvi Gaussovi integracijski točki levo od točke opazovališča B. Nihanje v izredno visokih frekvencah, ki je razvidno iz slik notranjih sil, je očitno posledica prisotnosti nihajnih oblik z nihajnimi časi, ki so bistveno manjši od računskega časovnega koraka, zato je v vsakem primeru smiselno prispevke teh oblik udušiti. Opozoriti pa moramo, da so to nihanja deformacij in ne osnovnih neznank (pomikov in zasuka) in da je gladkost deformacij garancija za gladek potek pomikov, obratno pa ni nujno res. Tiste notranje sile, ki jih dobimo z analizo brez dušenja, izkazujejo izredno velike Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 59 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 4.8: Bauchaujevo nihalo: temna obmocˇja so dusˇena. Figure 4.8: Bauchau swinging pendulum: dark areas mark presence of damping. amplitude nihanj osne sile, in te se zdijo v tej analizi najbolj problematične. Nizka stopnja dušenja skoraj popolnoma odpravi težave z oscilacijo osne sile, blagodejen učinek pa se pozna tudi na prečnih silah in upogibnih momentih. Z višanjem stopnje dušenja se vpliv višjih nihajnih oblik zmanjšuje, kar je pričakovano. Na koncu prikazujemo vpliv dušenja na frekvenčni spekter notranjih sil, ki ga dobimo z Matlabovo funkcijo FFT (‘Fast Fourier Transform’). Analiza je prikazana na sliki 4.10 in prikazuje veliko prisotnost visokih frekvenc v nedušenem odzivu in drastično izboljšanje z dušenjem. Ključno pri tej analizi pa je, da dušenje ne prizadene nizkih rrekvenc. Izjema je vpliv močnejšega dušenja na stnzne derormacije, ki so v tem primeru relativno nepomembne. Analiza frekvenčnega spektra pomikov in zasuka pokaže, da dušenje (v tem primeru) ne vpliva na frekvenčni sestav odziva pomikov in zasuka. 60 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 4.9: Bauchaujevo nihalo: notranje sile. Figure 4.9: Bauchau swinging pendulum: internal forces. Slika 4.10: Bauchaujevo nihalo: frekvencˇni sestav odziva notranjih sil. Figure 4.10: Bauchau swinging pendulum: frequency content of internal forces. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 61 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 4.9 Sklepi Formulacijo geometrijskih nosilcev, vpeljano v razdelkih 4.1-4.6, smo razširili z algoritmom za dodajanje umetnega dušenja na deformacijah. Za dušenje smo uporabili trigonometrično funkcijo arctan, ki omogoča večjo robustnost v primeru nenadnih skokov hitrosti. Dušenje smo nadgradili s pametnim algoritmom, ki sam določa kraj in trajanje dušenja s tem, da ga po potrebi vklaplja in izklaplja na nivoju integracijske točke. Algoritem je računsko zelo nezahteven, saj mora za odločitev poznati le nekaj preteklih predznakov inkrementov deformacij in največjo vrednost inkrementov deformacij v opazovanem obdobju. Zaznavanje prekomernega nihanja v visokih frekvencah poteka na nivoju deformacij, kar se je izkazalo kot bolj primerno od opazovanja nihanja pomikov in zasukov. Ključni test algoritma je količina udušenih osnovnih in višjih nihajnih oblik. Na računskem primeru se je pokazalo, da je algoritem sposoben udušiti višje nihajne oblike, ne da bi pomembneje vplival na osnovne oblike. Pomanjkljivost algoritma je, da faktorja dušenja adam ne moremo povezati s fizikalnim parametrom kritičnega dušenja, kar je sicer splošna pomanjkljivost shem z numeričnim dušenjem. Optimalno vrednost parametra je potrebno poiskati s poskušanjem za vsak primer posebej. Ker algoritem ne temelji na teoretičnem dokazu, temveč bolj na izkušnjah in opazovanjih, ga imenujemo ‘hevrističen’. Koncept in algoritem pa je možno aplicirati na poljubno formulacijo. V računskih primerih smo prikazali le en primer, več jih je moč najti v priloženem (priloga D) članku Gamsa, Planinca in Sajeta (2007d). 62 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 5 OPTIMIZACIJA V DINAMIKI LINEARNO ELASTICˇ NIH LINIJSKIH NOSILCEV 5.1 Uvod Optimizacijo dinamičnih sistemov lahko opišemo kot iskanje take kombinacije geometrijskih, materialnih in obtežnih podatkov, da je določena lastnost sistema v nekem specifičnem času minimalna ali maksimalna. Iskanje je lahko podvrženo različnim omejitvam. Če znamo rešiti tak problem, lahko stvari naredimo ceneje, hitreje, z manj materiala in z manjšo porabo energije. V gradbeništvu lahko na primer naredimo lažje ali bolje oblikovane nosilne konstrukcije, ki so bolj odporne na dinamične obremenitve, kot sta veter in potres. Ovira širši uporabi optimizacije je njena težavnost, saj so računi relativno zapleteni - po eni strani zaradi zapletene matematične formulacije problema, po drugi pa zaradi numerične zahtevnosti. Nemalokrat se namreč zgodi, da optimizacijski algoritem ne najde sprejemljive rešitve. Za reševanje optimizacijskih problemov obstaja nekaj konceptualno različnih pristopov. Značilen primer so genetski oz. evolucijski algoritmi, ki iz začetne populacije s postopki, ki oponašajo naravno evolucijo (dedovanje, mutacije, ipd.), določijo optimalne primerke. Vendar pa se tovrstni algoritmi v gradbeništvu niso izkazali oz. uveljavili, ker so bistveno počasnejši od njihovih tekmecev - gradientnih algoritmov. Gradientni algoritmi temeljijo na znanih odvodih po projektnih spremenljivkah, s katerimi se v smislu inkrementne gradientne metode približujemo optimalni konstrukciji. V tem delu se ukvarjamo le z gra-dientnimi algoritmi. V postopku optimizacije se oblika in lastnosti konstrukcije spreminjajo na nepredvidljiv način. Zato je možno, da se bo tekom tega postopka količina uporabljenega materiala zmanjšala do te mere, da bo konstrukcija postala bolj podajna. Podajna konstrukcija pa lahko v analizi ‘uide’ iz območja veljavnosti linearne teorije. Zato je potrebno upoštevati kinematično nelinearne zveze, kot so na primer Reissnerjeve enačbe, ki imajo bistveno večje območje veljavnosti, saj velikosti deformacij in pomikov niso omejene. Dodatno se lahko zgodi, da postane sistem enačb konstrukcije v tem postopku ‘tog’. Kot vemo že iz prejšnjih poglavij, je numerično reševanje togih sistemov problematično, zato je potrebno v analizi uporabljati stabilne integracijske sheme, kakršne so na primer energijske integracijske sheme. Optimizacija je samostojno raziskovalno področje. V tem poglavju ne izboljšujemo optimizacijskega postopka, temveč je prispevek le aplikacija optimizacije na novo razviti šibki integracijski shemi s povprečno rotacijsko matriko s poudarkom na obravnavanju posebnosti formulacije pri občutljivostni analizi. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 63 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 5.2 Optimizacijski problem Matematična formulacija. V matematični zapis enačb problema vpeljemo ti. projektne spremenljivke b, s katerimi opišemo vse spremenljive lastnosti konstrukcije. Cilj je določiti take vrednosti teh spremenljivk, da bo zadoščeno optimizacijskemu problemu. Slednji v matematičnem zapisu vsebuje enačbe, neenačbe in zahtevo po minimumu ali maksimumu neke funkcije (tako definirani problemi spadajo v področje nelinearnega matematičnega programiranja): min/o, (5.1) fi < 0, i = 1,..., npog, (5.2) gm = 0. (5.3) Tu smo s /o označili skalarno namensko funkcijo, s fi omejitvene pogoje in z %>0g število omejitvenih pogojev. Namenska funkcija /o in omejitveni pogoji fi so lahko definirani pri poljubnih časih. Enačba (5.3) predstavlja vodilne enačbe mehanskega problema. Rešitev iščemo na intervalu b < b < b , (5.4) kjer z b in b označujemo spodnjo in zgornjo mejo projektnih spremenljivk. Analiza odziva konstrukcije je reševanje mehanskega problema, določenega z enačbami (5.3) oziroma (4.111) in s pripadajočimi robnimi in začetnimi pogoji (3.5), (3.6) in (3.7). Rešitev so diskretne spremenljivke T na časovnem intervalu t G [0, Lend]> kjer je tend čas, do katerega poteka analiza. Kot smo že omenili v prejšnjih poglavjih, je T = [U, $] vektor vseh diskretnih neznank. Ker so lastnosti konstrukcije sedaj odvisne od projektnih spremenljivk b, je tudi odziv T = T(t, b) funkcija projektnih spremenljivk. Ker so namenska funkcija in omejitveni pogoji odvedljivi po projektnih spremenljivkah, lahko opti-mizacijski problem (5.1)—(5.2) rešujemo s konveksnimi gradientnimi metodami matematičnega programiranja. V tem delu uporabimo reševanje s postopkom, ki so ga razvili Kegl, Butinar in Kegl (2002) in Kegl in Brank (2006). Namenska funkcija. Izbira namenske funkcije je odvisna od tega, kaj želimo minimizirati. V gradbeništvu se pogosto odločimo, da je to volumen materiala, lahko pa minimiziramo tudi druge količine. Na primer, lahko zahtevamo, da sistem pri predpisanem času obmiruje ali pa minimiziramo deformacij-sko energijo, kar sili konstrukcijo v drugačne mehanizme prenašanja obremenitev, npr. z osnimi silami namesto z upogibnimi momenti. Možno je tudi minimiziranje posameznih komponent neznank (pomika ali zasuka). Problem namenske funkcije je, da je skalarna in da jo lahko zapišemo le pri enem specifičnem času - minimiziramo lahko le eno funkcijo. Pogosto pa želimo minimizirati več količin hkrati, ali pa eno količino pri različnih časih. Za take primere lahko sestavimo namensko funkcijo na več načinov, od katerih je najbolj preprost ta, da količine, ki jih želimo minimizirati, seštejemo ali množimo in če so te količine podobnega velikostnega reda, pristop zadovoljivo deluje. Bolj robusten pa je pristop, ki za namensko funkcijo vpelje novo projektno spremenljivko 6new in jo preko omejitvenih pogojev poveže s količinami, ki jih želimo minimizirati. Če želimo npr. minimizirati funkciji /o,i in /0,2, lahko ta pristop matematično zapišemo tako: /0,1 ^ ^new /o,2 < &new (5.5) /o = &new 64 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. V splosˇnem lahko tako minimiziramo poljubno sˇtevilo funkcij. 5.3 Projektne spremenljivke Za opis parametrov uporabljamo tri vrste projektnih spremenljivk, ki jih lahko klasificiramo kot: klasicˇne, posredne in projektne spremenljivke vzbujanja. Shematsko je uporaba teh spremenljivk prikazana na sliki 5.1. Klasicˇne projektne spremenljivke. K njim sˇtejemo vse projektne spremenljivke, s katerimi neposredno opisujemo posamezne parametre v metodi koncˇnih elementov. Tipicˇen primer takih spremenljivk so dimenzije precˇnega prereza (sizing optimization), uporablja pa se jih tudi pri optimizaciji oblike (shape optimization). Pri optimizaciji oblike se take projektne spremenljivke priredi koordinatam posameznega vozlisˇcˇa. Tak pristop se zdi najbolj naraven, saj zaradi neposrednosti ni potrebno vpeljevati dodatnih spremenljivk in se ga posluzˇuje veliko avtorjev (Pedersen in Nielsen (2003), Stocki in sodelavci (2001), Wang, Zhang in Jiang (2002a, b, 2004)), vendar pa ni najbolj primeren, saj pogosto vrne negladke oblike (Harl in Kegl (2005)). Posredne projektne spremenljivke oblike. S tem terminom oznacˇujemo vse projektne spremenljivke, s katerimi preko projektnih teles (design bodies) opisˇemo obliko konstrukcije. Razlika v primerjavi s klasicˇnimi spremenljivkami oblike je v tem, da projektne spremenljivke tokrat niso vezane neposredno na podatke v metodi koncˇnih elementov, temvecˇ na projektna telesa. Ta predstavljajo del ali pa celotno konstrukcijo, na katero je vezana mrezˇa koncˇnih elementov. Bistvena prednost tehnike s projektnimi telesi in projektnimi spremenljivkami oblike je, da obliko konstrukcije opisˇemo z relativno majhnim naborom spremenljivk in da so rezultati boljsˇi, saj prevzamejo ugodne lastnosti projektnih elementov. V tem delu za opis projektnih teles uporabimo Be´zierjeva projektna telesa, ki se po porocˇanju raziskovalcev dobro obnesejo (Harl in Kegl (2005), Kegl (2000, 2005), Kegl in Brank (2006)). Omogocˇajo namrecˇ izreden nabor oblik, dajejo gladke oblike in so sposobna popolno opisati osnovne geometrijske oblike (Kegl in Brank (2006)). Projektne spremenljivke vzbujanja. Za opis dinamicˇne obtezˇbe v odvisnosti od cˇasa uporabimo projektne spremenljivke vzbujanja. Ker se s takimi projektnimi spremenljivkami lahko opisˇe bistvene robne pogoje, ki v enacˇbe ne vnasˇajo novih neznank, je pristop nekoliko poseben. Podrobneje je razlozˇen v prilogi E. Optimizacija obtezˇbe je zelo pomembna pri problemih krmiljenja – npr. za dolocˇitev funkcij zasukov podajnih robotskih rok in manipulatorjev. Tipicˇen primer uporabe je tudi minimiziranje cˇasa dolocˇene aktivnosti. 5.4 Obcˇutljivostna analiza Obcˇutljivostna analiza (sensitivity analysis) je tisti del optimizacijskega postopka, kjer se racˇuna odvode po projektnih spremenljivkah. Ta del je praviloma zelo zahteven in mora uposˇtevati posebnosti formulacij. V tem poglavju se omejimo na sˇibki algoritem s povprecˇno rotacijsko matriko, ki ima v smislu numericˇnega racˇuna pomembne posebnosti, ki jih moramo uposˇtevati pri racˇunu odvodov po projektnih spremenljivkah. Obcˇutljivostna analiza je racˇun odvodov namenske funkcije in omejitvenih pogojev po projektnih spre- Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 65 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. I2 in grafično prikazana na sliki 5.3. V računskem modelu smo za vsak steber in vsako prečko uporabili po en končni element s kvadratično interpolacijo in šibki časovni integrator s povprečno rotacijsko matriko. Časovni korak At = 0.02 s je konstanten. Projektne spremenljivke so uporabljene za opis projektnih elementov, s katerimi modeliramo po višini sprejemljivo širino okvirja. Projektni elementi so označeni s črtkano črto na levi strani slike 5.2. Vsak projektni element je definiran z osmimi kontrolnimi točkami (CP). Dovoljene vrednosti projektnih spremenljivk so take, da je lahko skupna širina okvirja od 12 do 24 m, višina spodnje etaže je med 3 in 6 m. Širina notranje ladje je lahko od 4 do 10 m. Z namensko funkcijo minimiziramo nihanje konstrukcije po prenehanju delovanja sunka. S spreminjanjem geometrije konstrukcije želimo doseči, da bi okvir po zaključku sunka obmiroval, če je to mogoče. Najbolj primerna funkcija za zagotavljanje mirovanja je mehanska energija, saj ničelna mehanska energija zagotavlja absolutno mirovanje sistema. 1Racˇun odziva konstrukcije, vrednotenje namenske funkcije in omejitvenih pogojev ter obcˇutljivostno analiza smo vgradili v program EMS (Kegl in sodelavci). 2Uporabljen je bil optimizacijski program ‘iGO’. Avtor programa je M. Kegl. Algoritem programa je opisan v cˇlanku Kegl, Butinar in Kegl (2002). Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 67 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 5.2: Racˇunski model okvirja. CP je oznaka za kontrolno tocˇko, DE pa za projektno telo. Figure 5.2: The model of the frame. CP are control points and DE design elements. Slika 5.3: Pomik podpor v odvisnosti od časa. Figure 5.3: Movement of supports vs time. Izkaže se, da tako zastavljenega optimizacijskega problema ne moremo uspešno rešiti v eni potezi. Postopek zato razbijemo na dve fazi: v prvi fazi minimiziramo vodoravne pomike na sredini vrhnje etaže pri časih 1 s in 1.2 s. Ko najdemo optimum tako zastavljenega problema, dodatno minimiziramo še pomike pri času 1.4 s. Za minimiziranje količin pri več časih uporabimo pristop, kije opisan v enačbi (5.5). Odziv konstrukcije po tej fazi optimizacije je prikazan na grafu na levi polovici slike 5.4. Kot je tam razvidno, so pomiki na sredini višine konstrukcije še vedno občutni, medtem ko se konstrukcija pri vrhu po sunku skoraj popolnoma umiri. 68 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. V drugi fazi smo za izhodišče vzeli optimirano obliko konstrukcije iz prve faze in kot namensko funkcijo določili mehansko energijo. Postopek je v tem primeru uspešno skonvergiral in konstrukcija po sunku skoraj popolnoma miruje, kot kaže graf na desni strani slike 5.4. Oblika konstrukcije po prvi fazi računa in končna, optimalna oblika konstrukcije sta prikazani na sliki 5.5. Več primerov z zanimivimi izbori namenskih funkcij in obtežb je prikazanih v prilogi E. Slika 5.4: Pomiki na sredini visˇine konstrukcije s cˇrtkano cˇrto in na vrhu konstrukcije s polno cˇrto po prvi fazi optimizacije, levo in po koncˇani optimizaciji, desno. Figure 5.4: Displacements at mid-height and at the top of the structure are plotted with dashed and full line, respectively. The results of the first phase are on the left hand side and the final results on the right hand side of the figure. Slika 5.5: Oblika konstrukcije po prvi fazi optimizacije, levo in po koncˇani optimizaciji, desno. Figure 5.5: The shape of the structure after the first phase of optimization is on the left hand side and the final shape on the right hand side of the figure. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 69 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 5.7 Sklepi V tem poglavju smo predstavili postopek optimiranja dinamičnih sistemov. Osredotočili smo se na opis občutljivostne analize in njeno aplikacijo na novo razviti časovni integrator za ravninske nosilce - šibki integrator s povprečno vrednostjo rotacijske matrike. Z vpeljavo treh vrst projektnih spremenljivk: klasičnih, posrednih oblikovnih in spremenljivk vzbujanja smo sposobni modelirati izredno velik nabor inženirskih problemov. V računskem primeru smo demonstrirali le uporabo posrednih projektnih spremenljivk oblike, v prilogi E pa lahko najdemo še zanimive primere optimizacije robotskih manipula-torjev. Bézierjeva telesa so se izkazala kot primerna za opis oblike v optimizaciji dinamičnih problemov. Pri optimizacijskih problemih pogosto želimo optimirati več kot eno količino. Splošnega pristopa k optimizaciji večih funkcij ni, pristop z uvedbo nove projektne spremenljivke, ki smo ga prikazali, pa je bil uspešen. Dodaten problem v dinamiki je, kako upoštevati časovno odvisne pogoje, še posebej, če so podani na časovnem intervalu. Literatura in naši primeri so pokazali, da lahko take pogoje nadomestimo z manjšim številom pogojev pri specifičnih (diskretnih) časih. Pokazalo seje, da v konkretnem računskem primeru optimizacije nismo bili sposobni uspešno opraviti v enem koraku inje bilo potrebno postopek razbiti na dva koraka. Optimizator tipično ne more konvergirati v optimum, kadar so neugodno izbrane začetne projektne spremenljivke, saj se optimizacija ujame v nesprejemljive lokalne minimume, v katerih so lahko celo kršeni omejitveni pogoji. Vendar seje pokazalo, da ko se rešitev zadosti približa optimumu, postopek zelo hitro, v nekaj iteracijah skonvergira. Da bi našli ugodna izhodišča za uspešno optimizacijo, je potrebno nekoliko izvirnosti v strategiji reševanja in pripravljenost poskušati. Izkušnje tako ostajajo pomemben del uspešne optimizacije. Nazadnje moramo poudariti koristnost programa za simbolično računanje odvodov in avtomatsko generacijo kode AceGen, s katerim smo lahko pripravili podatke za občutljivostno analizo nekajkrat hitreje, kot če bi to opravili ‘peš'. 70 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 6 DINAMIKA ARMIRANOBETONSKIH LINIJSKIH KONSTRUKCIJ 6.1 Uvod Z upoštevanjem realnih lastnosti materialov v dinamični analizi lahko izboljšamo natančnost opisa obnašanja gradbenih konstrukcij. Delo, ki ga predstavljamo v pričujočem poglavju, je prispevek na tem področju, saj poleg točnih geometrijskih zvez upošteva tudi realne lastnosti armiranega betona. Klasične analize armiranobetonskih konstrukcij upoštevajo bodisi linearne kinematične zveze, P — A metodo ali različne poenostavljene verzije teorij drugega reda. V tem poglavju ugotavljamo, če analize s točnimi geometrijskimi zvezami bolj realno opišejo obnašanje armiranobetonskih konstrukcij pod vplivom dinamične obtežbe in kako pomembna je izbira časovne integracije. Pri velikih obremenitvah se v konstrukciji oblikujejo območja, kjer se poškodbe ‘skoncentrirajo’ in jih imenujemo plastični členki. Taka območja se praviloma pojavijo na mestih vpetja stebrov in prečk. S pojavom plastičnega členka v elementu se točke v okolici členka razbremenijo, deformacije in poškodbe pa se bolj in bolj skoncentrirajo v plastičnem členku. S stališča implementacije materialnega modela v formulacijo ni večjih sprememb v algoritmu računa, kakor je bil prikazan v poglavju ‘4.4 Numerični algoritem’. Razlika je le v tem, da se tokrat konstitucijski sili Nc in Mc računata z numerično integracijo po prerezu. Ta dodatni postopek je podrobneje opisan v nadaljevanju. 6.2 Materialni modeli Racˇunski modeli za material so, tako kot vsi racˇunski modeli, idealizacija realnega materialnega stanja. Dobljeni so na podlagi standardiziranih tlacˇnih in nateznih preizkusov in podajajo enoosno odvisnost med napetostmi in deformacijami. V analizi jih definiramo za vsako posamezno materialno vlakno, iz katerih sestavimo precˇni prerez. Racˇunski materialni modeli, ki jih uporabljamo v tem delu, so po svoji matematicˇni obliki eksplicitni. To pomeni, da lahko pri znani deformaciji in zgodovini dogajanja eksplicitno izrazimo oz. izracˇunamo napetost v vlaknu. Kljub eksplicitni obliki pa so za uporabo relativno zapleteni, ker uposˇtevajo zgodovino napetostnega stanja. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 71 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 6.2.1 Jeklo Za matematično modeliranje jekla uporabimo materialni model, ki sta ga razvila Menegotto in Pinto (1973), njegova značilnost pa je, da z ukrivljenim prehodom iz elastičnega v plastično območje zmore dobro modelirati Bauschingerjev efekt. Računski model v osnovi upošteva samo kinematično utrjevanje, dodatno pa je vanj vgrajena tudi možnost upoštevanja izotropnega utrjevanja po Filippou, Popovu in Berteru (1983). V literaturi je ta model zelo cenjen, ker lahko z njim zelo natančno ujamemo eksperimentalne rezultate cikličnega obremenjevanja armaturnih palic. Parametri modela. So zbrani v preglednici 6.1. Preglednica 6.1: Parametri materialnega modela za jeklo. Table 6.1: Parameters of the steel material model. Parameter Ime Priporočena vrednost Eq Začetni elastični modul Ponavadi na voljo h Količnik utrjevanja Ponavadi na voljo fy Meja plastičnosti Ponavadi na voljo Ro Začetna ukrivljenost prehoda 10 do 20 cR\ Vpliv na ukrivljenost 0.925 Ci?2 0.15 a\, a2 Vpliv na izotropno utrjevanje o\ = 0 in Ü2 = 1 ni v tlaku izotropnega utrjevanja v tlaku as, (14 Vpliv na izotropno utrjevanje Os = 0 in Ü4 = 1 ni v nategu izotropnega utrjevanja v nategu Enačbe materialnega modela. Napetost v materialnem vlaknu je podana v eksplicitni odvisnosti od njegove deformacije in zgodovine obremenjevanja. Postopek upoštevanja zgodovine obremenjevanja razdeli opazovanje dogajanja v materialnem vlaknu na faze. Vsaka faza predstavlja režim monotonega spreminjanja deformacije. Prehod med fazami se sproži s spremembo predznaka prirastka deformacije Ae. Napetost v vlaknu pri dani deformaciji je določena z enačbo: a* = h e* + (1 — h) e* i /i i (f*\ ) kjer sta e* in a* definirana z zvezama in e* a* = L — Lr Lq — Lr a — ar (6.1) (6.2) (6.3) do — or Točki To (eo, /, Lo V tej enačbi je sy konstantna deformacija na meji plastičnosti: Ly = fy Eq (6.4) (6.5) (6.6) V enačbi (6.5) je Lpiimax največja (do tedaj) dosežena plastična deformacija, ki predstavlja v fazah s pozitivnimi inkrementi deformacije največjo (do tedaj) doseženo plastično deformacijo v nategu (smax)> in v fazah z negativnimi inkrementi deformacije največjo (do tedaj) doseženo plastično deformacijo v tlaku (čmin). Ti dve vrednosti shranjujemo za vsako vlakno. Začetna vrednost maksimalne plastične deformacije za neobremenjeno vlakno je LpZ,max (zaˇcetni) ( Ly, ˇce ?s > 0 — Ly, ˇce ?L < 0 (6.7) v nadaljevanju pa se določa glede na doseženo maksimalno plastično deformacijo: %>Z,max J max(S|/,Lmax), ˇce ?s > 0 min(—Ly, Lmin), ˇce ?s < 0 (6.8) Izotropno utrjevanje v tlaku upoštevamo tako (Filippou, Popov in Bertero (1983)), da presečišče poltraka in premice izračunamo po spremenjeni formuli: Lq = —fy at + h Eq Ly at — katerih pomen bomo opisali v nadaljevanju, v vrednostih napetosti ob koncu prejšnjega koraka apre in v razbremenitvenih tangentnih modulih v tlačni (i?Cc,raz) in natezni coni (EctiTdbZ). A) Ovojnica napetosti v tlaku. Sekantni modul Esec določimo po enačbi -^sec fc Glede na vrednost sekantnega modula je definirana pomožna spremenljivka r: T == 400, Esec > Eo 1 SeC — U Ec0-E sec Po vpeljavi spremenljivke r\ i] = ^c0 (6.17) (6.18) (6.19) z Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 75 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. lahko zapišemo enačbo za napetost r a = jcr]------------------- (6.20) (r — 1 + rf) in tangentni modul (r — 1) (1 — i] ) Ect = tcr-------------------------^. (6.21) p n (r — ^ -\- rir'\ ccu ^ / / B) Ovojnica napetosti v nategu. Po računu deformacije na meji plastičnosti v nategu, Jet eto = , (6.22) Eq izračunamo napetost in tangentni modul po enačbah: e Eq, eto < e < 0 CT = fa ßLtu-m } etu 2 „ ^ J = /2in/ /- • {626) 0.145 Ç +U.13Ç, Ç<2 Z razmerjem r^ j lahko določimo sečišče razbremenitve in abscise, ki je v koordinatnem sistemu e — a določeno s koordinatami (eend, 0): eend = rKj eco- (6.27) Z znanim sečiščem lahko izračunamo naklon razbremenitvene daljice: V^min ^endj ^ ß . Eco, (emm — eend) > w~ (*) cu, v nun e _ ßcQ \ j Eccraz = L^in Lnd mm x ™aaff • (6.28) Če je v enačbi (6.28) merodajen drugi pogoj (*), je potrebno deformacijo eend računati po pravilu: a Lend = črnin — J^ ¦ (6.29) emm predstavlja najmanjšo doseženo deformacijo, in si jo moramo, tako kot deformaciji eend in emax, za vsako vlakno shranjevati. Vlakno, ki ima določeno stopnjo poškodovanosti v tlaku (eend 7^ 0), pri ponovnem obremenjevanju do deformacije eend ne nudi odpora. Šele ko je ta deformacija presežena, vlakno nudi odpornost, kije proporcionalna tangentnemu modulu ECCti:az. D) Pravila za razbremenjevanje v natezni coni. Razbremenjevanje iz katerekoli točke v natezni coni je predvideno skozi koordinatno izhodišče: Ect^z = -¦ (6.30) e 76 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Obnašanje modela. Materialno vlakno obremenimo z deformacijo z naraščajočim žagaš tim vzorcem, ki je prikazan na levi strani slike 6.3. Pripadajoča napetost v betonskem vlaknu z materialnimi karakteristikami T7I o vx lO 2 p o w 7 2 p 0 w 6 2 A?o — o.2 x 10 N/m , jc = -d x 10 N/m , jct = ^x 10 N/m , Lco = -0.00214, e&r = -0.0069, e^c/ = 0.0007, (6.31) ß = 0.1 je prikazana na desni strani slike 6.3. Slika 6.3: Zgodovina deformacije v odvisnosti od ‘cˇasa’, levo. Pripadajocˇa napetost, desno. Figure 6.3: Strain history in the fibre, left. Corresponding stress, right. 6.3 Integracija po prerezu Pri postopku numerične integracije upoštevamo koncept, daje prečni prerez sestavljen iz vlaken. Vlakna so, kot smo omenili v uvodnem poglavju, enodimenzionalne materialne niti, ki potekajo po celotni dolžini elementa. Zaradi zvezne porazdeljenosti mase po nosilcu je v prerezu neskončno materialnih vlaken. Predpostavka je, da napetostno stanje v posameznem vlaknu nima neposrednega vpliva na sosednja vlakna. Zato lahko v vlaknu uporabimo enoosni materialni model. Zaradi natančnejše numerične integracije prerez nosilca razdelimo na podprereze, na vsakem podprerezu pa uporabimo ustrezno numerično metodo integracije. Uporabljamo Gaussovo ali Lobattovo integracijo. Takšna integracija zahteva podatke o deformaciji, napetosti in materialnem modelu v izbranem vlaknu. Te dobimo, brž ko so znani pomiki nosilca, vrsta materiala in njegov konstitucijski zakon. Pri ravninski analizi okvirjev imajo podpodročja obliko vzporednih lamel (slojev), v katerih se mehanske količine spreminjajo samo v prečni smeri na lamele, to je, le v odvisnosti od lokalne koordinate z in zato integracija po prerezu poteka dejansko samo po koordinati z - po višini prereza. V splošnem so poteki napetosti in materialnih modelov po prerezu pri času tn zvezne ali celo nezvezne funkcije. Z delitvijo prereza na lamele in vpeljavo numerične integracije funkcijo zamenjamo z množico diskretnih vrednosti v pripadajočih diskretnih točkah. Izbira točk je pomembna s stališča natančnosti numerične integracije. O vplivu izbire na natančnost smo govorili že v poglavju 3.7 in vse ugotovitve veljajo tudi v primeru integracije po prerezu. Lobattovi integraciji dajemo prednost pred Gaussovo, ker zajame tudi skrajne robove vlaken in z večanjem stopnje integracije izboljšujemo natančnost integrala. V analizi moramo upoštevati še dve posebnosti, ki vplivata na način delitve prereza na lamele. Prva Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 77 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. zadeva obliko prečnega prereza. Ta je pogosto nepravokotna; tak je T-prerez. V teh primerih moramo prerez razdeliti na podpodročja z zveznim robom. T-prerez moramo deliti na območji pasnice in stojine. Druga posebnost je oblika funkcij, kijih integriramo. Te so lahko take, da so na nekem območju prereza enake nič, na preostalem delu pa imajo od nič različne vrednosti. Meja med ničelnimi in neničelnimi vrednostmi pa je v vsakem časovnem koraku drugje. Zato je optimalno, da prerez razdelimo na več lamel, po katerih integriramo s poljubno stopnjo integracije. Z delitvijo na več lamel povečamo možnost, da bo meja lamele blizu meje med ničelnimi in neničelnimi vrednostmi, z uporabo integracije višje stopnje pa hkrati občutno izboljšamo natančnost integracije. Nekaj primerov različnih integracij je prikazanih na desni strani slike 6.4. Pri integraciji moramo upoštevati še lastnost betona, da ima s stremeni objeti del povečano tlačno trdnost fc in bistveno povečani deformaciji eco in ec\j. Karakteristike objetega betona izračunamo po Mander-jevem modelu (Mander, Priestley in Park (1988)). Objeti in neobjeti del betona ločeno razdelimo na lamele (glej sliko 6.4). Objeti del betona ‘razrežemo’ na niam lamel, kot je prikazano na osrednjem delu slike 6.4. Na podoben način ‘razrežemo’ tudi neobjeti del betona (zaščitno plast) in dobimo še dodatno nniam lamel neobjetega betona. Ko predpostavimo še vrsto in stopnjo integracije po lamelah, lahko prispevke posameznih lamel seštejemo in dobimo konstitucijsko osno silo Nq in konstitucijski moment Me- Nekaj možnih izbir integracije je prikazanih na desni strani slike 6.4. Število lamel oz. integracijskih točk po lamelah ima zelo pomemben vpliv na hitrost računa. V literaturi je ponavadi integracija po lamelah implementirana tako, da se za posamezno lamelo predpostavi, da imajo v njej vsa vlakna enako, konstantno napetost (npr. Taucer, Spacone in Filippou (1991)). Ta predpostavka je enakovredna enotocˇkovni Gaussovi integraciji (slika 6.4) in zahteva sorazmerno vecˇje sˇtevilo lamel po prerezu za enako natancˇnost kot vecˇtocˇkovna integracija. Slika 6.4: Prerez, skrajno levo. Razdelitev prereza na lamele in njihovo osˇtevilcˇenje, sredina. Ilustracija treh mozˇnih integracij po prerezu, desno. Figure 6.4: Cross-section, left. Division and enumeration of the cross-section into layers, middle. Illustrations of some of the possible integration schemes, right. Koordinate Zij integracijskih točk so določene glede na izhodišče koordinatnega sistema y — z, ki je 78 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. postavljeno v težišče prereza. Indeks i označuje lamelo, indeks j pa integracijsko točko. Višina lamele objetega betona je h\amti, njena širina je bij; za lamelo neobjetega betona veljajo oznake n/Mam,i in nbi. V implementaciji postopka je predpostavljeno, da se lahko širine lamel objetega betona po višini spreminjajo konstantno ali linearno, širine lamel neobjetega betona pa so konstantne. Prvi korak analize je določitev lastnosti prečnega prereza, ki jih lahko priročno določimo po enačbah (formulah), zbranih npr. v diplomi Markovičeve (2006). Te formule za račun težišča, ploščine in vztraj-nostnih momentov uporabljajo koordinate oglišč poligona, s katerim opišemo zunanjo konturo prereza. Vzdolžne deformacije e v integracijskih točkah izračunamo po enačbi €ij = e + k Zij (6.32) za vse lamele i = 1,..., n\am in vse integracijske točke po lamelah j = 1,..., ^itg,iam objetega betona. Analogno velja za neobjeti beton: n€ij = e + k nZij. (6.33) Z znanimi vzdolžnimi deformacijami lahko s pomočjo enačb (6.17)-(6.30) izračunamo pripadajoče napetosti ac(eij). Armaturne palice obravnavamo točkovno. Posamezne palice označimo z indeksom k = 1,..., npa\ in njihove ploščine z Apa\tk- Vzdolžno deformacijo v njih računamo analogno kot v betonskih vlaknih, iz nje pa z enačbami (6.1)—(6.15) določimo napetost 8/5 cm ob vpetju (do višine 90 cm) in z enakimi stremeni <3> 8/15 cm nad to višino. Prečni prerez stebra je prikazan na sliki 6.13. Zaščitna plast betona je 2.4 cm, statična višina pa 36 cm. Slika 6.13: Prečni prerez stebra z vzdolžno in strižno armaturo. Figure 6.13: Cross-section of the beam with longitudinal and shear reinforcement. Materialni parametri za beton so: fc = -55 MPa, Ecq = 37000 MPa, sco = -0.0018, ecu = -0.0035, (6.39) fct = 0 MPa, e tu = 0.00015, vc = 0.2, ß = 0.1, Vrednost elastičnega modula je določena po formuli Ec = 57000\/|/c|- Natezne trdnosti betona v 86 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. racˇunih nismo uposˇtevali, saj se je pokazalo, da uposˇtevanje natezne trdnosti vrne bistveno pretog odziv konstrukcije ob zacˇetku obremenjevanja. To dodatno utemeljujemo z dejstvom, da je bila konstrukcija do dolocˇene mere obremenjena pred samimi testi z namenom umeritve merilnih naprav in potisnih batov, kar je lahko povzrocˇilo razpokanje betona. Karakteristike objetega betona so dolocˇene po Manderjevem modelu objetega betona: ° Je 1fc — -70.5 MPa, ° Je 1Lcq = -0.006, ° Je 1Lcu = -0.019. (6.40) Materialne karakteristike armature so: fy = 55.5 MPa, i?o = 200 GPa, vs = 0.3, h = 0.01, Rq = 15, cR\ = 0.925, ci?2 = 0.15, CL\ = 0, Ü2 = 1, &3 = 0, a4 = 1. (6.41) Karakteristike so bile določene na podlagi specifične kvalitete vgrajenih materialov (beton C40/50 in jeklo B500H). Izmerjene vrednosti žal niso bile na voljo. Računski model in obremenitve. Ker je eksperiment pokazal, da se strešna konstrukcija obnaša kot toga diafragma in ker je konstrukcija obremenjena le v eni smeri, lahko uporabimo ravninski računski model. Konstrukcijo modeliramo z enim, 5 m visokim, togo vpetim stebrom, na vrhu katerega se nahaja šestina mase strešne konstrukcije (slika 6.14). V računu uporabimo 5 elementov s kvadratično interpolacijo pomikov vzdolž elementa in reducirano Lobattovo integracijo za togostne člene in polno za vztrajnostne. Pri integraciji po prerezu uporabimo tri lamele s 5 točkovno integracijo. V računu upoštevamo, da je edina masa točkovna masa na vrhu stebra. Namen primera je opazovati vpliv različnih časovnih integracijskih shem, zato primer rešimo s šestimi integracijskimi shemami. Slika 6.14: Računski model konstrukcije. Figure 6.14: Numerical model of the structure. V analizi je upoštevano 3% dušenje na prvi in drugi nihajni obliki; matrika dušenja je sestavljena kot linearna kombinacija masne in togostne matrike. Literatura (npr. Paz in Leigh (2004)) priporoča, da se za armiranobetonske konstrukcije tovrstno dušenje upošteva s 3 do 5% deležem kritičnega dušenja. Ker gre za psevdo-dinamični eksperiment, smo se odločili za spodnjo mejo s tega intervala. Konstrukcija je bila psevdo-dinamično obremenjena z akcelerogramom potresa v Črni gori, ki je bil umetno prilagojen na projektne pospeške 0.05 g, 0.14 g, 0.35 g in 0.525 g. Prva analiza (0.05 g) je služila le za preverjanje elastičnih lastnosti konstrukcije in je v računu ne upoštevamo. Ostale tri dinamične obtežne primere preračunamo z numeričnim postopkom. Pri eksperimentu je bila konstrukcija Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 87 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 6.15: Normaliziran akcelerogram. Figure 6.15: Normalized accelerogram. po dinamičnih preizkusih obremenjena še s kontroliranim pomikom do porušitve, vendar te faze v tem primeru ne obravnavamo. Akcelerogram je prikazan na sliki 6.15. Rezultati numeričnih analiz. Rezultati analiz problema z različnimi integratorji so praktično identični. Na grafu se razlik med rezultati ne da ločiti s prostim očesom niti v primeru, ko so v analizi uporabljeni zelo veliki časovni koraki, zato v nadaljevanju prikazujemo le rezultate, ki so dobljeni z Newmarkovim integratorjem s trapeznim pravilom. V analizi je uporabljen konstanten časovni korak ?t = 0.02 s. V sklopu primerjav z eksperimentalnimi rezultati primerjamo vodoravni pomik na vrhu konstrukcije u m in strižno silo ob vpetju Q. Rezultati so prikazani na slikah 6.16, 6.17 in 6.18 za obremenitve, skalirane na 0.14 g, 0.35 g in 0.525 g. Iz prvega grafa na sliki 6.16, ki prikazuje zvezo med Q in um, lahko razberemo, daje začetna togost konstrukcije precenjena kljub temu, da natezne trdnosti betona nismo upoštevali, zatem pa izkazuje celo nekoliko prenizko togost. Neujemanje z eksperimentom je najverjetneje posledica pojavov, ki v našem računskem modelu niso upoštevani. Dodatno lahko opazimo, daje računska histereza ‘preozka’. Histereza dobi ‘širino’ zaradi plastifikacije materiala, ki se v naši analizi ni pojavila. To je posledica dejstva, da smo v računu uporabili materialne karakteristike osnovane na projektni kakovosti vgrajenih materialov in ne na dejanskih izmerjenih karakteristikah. Na srednjem grafu je prikazan vodoravni pomik na vrhu konstrukcije v odvisnosti od časa. Opazimo lahko, da so pomiki konstrukcije nekoliko podcenjeni, predvsem v drugi polovici analize; računski nihajni časi analize se skoraj popolnoma ujamejo z eksperimentalnimi. Strižna sila je, tako kot pomiki, nekoliko podcenjena. Na sliki 6.17 je prikazan odziv konstrukcije na obremenitev, ki je skalirana na pospešek zemeljskih tal 0.35 g. Histereza je zopet preozka, sicer pa ne prekorači eksperimentalne ovojnice in kaže primerljivo togost. Pomiki in strižne sile se zelo dobro ujamejo z eksperimentalnimi, so le za malenkost podcenjeni. Primerjeva zadnje faza dinamičnega preizkušanja in numerične analize je prikazana na sliki 6.18. Tokrat se dobro ujamejo vse merjene količine, vključno s histerezo (zgornji graf), ki je do sedaj izkazovala največja razhajanja med numeričnimi in eksperimentalnimi rezultati. Primer je bil v osnovi mišljen kot primerjava časovnih integratorjev med seboj na primeru armiranobetonske konstrukcije. Izkazalo seje, da so vsi integratorji problem rešili enako kvalitetno celo v primerih analiz z velikimi časovnimi koraki, kjer bi se morale najprej pokazati razlike. Primer pa nam služi 88 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 6.16: Odziv na obremenitev skalirano z 0.14 g. Figure 6.16: Response to loads scaled to 0.14 g. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 89 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 6.17: Odziv na obremenitev skalirano z 0.35 g. Figure 6.17: Response to loads scaled to 0.35 g. 90 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Slika 6.18: Odziv na obremenitev skalirano z 0.525 g. Figure 6.18: Response to loads scaled to 0.525 g. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 91 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. tudi kot dodatna validacija racˇunskega postopka, saj se je pokazalo, da je ujemanje med racˇunskimi in eksperimentalnimi rezultati primerno. 6.5 Sklepi V tem poglavju smo izpeljali in aplicirali racˇunski postopek za geometrijsko tocˇno analizo ravninskih armiranobetonskih konstrukcij. Namen poglavja je bil raziskati vpliv razlicˇnih cˇasovno integracijskih shem na odziv takih konstrukcij. Na racˇunskih primerih realnih konstrukcij, ki so bili podprti z eksperimentalnimi rezultati, se je pokazalo, da je vpliv cˇasovnega integratorja pri analizi zanemarljiv. Delo predstavlja novost, saj v literaturi ni mocˇ zaslediti tako osnovane formulacije za geometrijsko tocˇno analizo armiranobetonskih konstrukcij. V sklopu numericˇnih primerov smo racˇunski postopek verificirali s primerjavo rezultatov z racˇunalnisˇkim programom Opensees in validirali na eksperimentalnih rezultatih ciklicˇne pushover analize in psevdo – dinamicˇne analize. 92 Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. __ v v __ __ __ __ 7 KNJIŽNICA KONČNIH ELEMENTOV 7.1 Uvod V sklopu dela z elementi z interpolacijo kinematičnih količin z različnimi materialnimi modeli je bila razvita obsežna knjižnica končnih elementov, ki so v tem poglavju našteti in na kratko opisani. Vsi elementi so dostopni na internetnem naslovu www.km.fgg.uni-lj.si/Matija/index.htm. Tam so na voljo rutine elementov {.dll) in izvorne kode elementov v programskem okolju AceGen. 7.1.1 Kinematične teorije V končnih elementih za statično analizo nastopajo tri kinematične teorije: linearna, teorija drugega reda in točna teorija. Enačbe točne teorije poznamo, enačbe linearne teorije in teorije drugega reda pa dobimo, če točno teorijo razvijemo v Taylorjevo vrsto okoli začetne (nedeformirane) lege. Reissnerjeve enačbe (2.25) so funkcije pomikov in jih ne moremo direktno razviti v Taylorjevo vrsto, zato si pomagamo s smernim odvodom (smerni odvod je posplošitev odvoda, ki predstavlja spremembo poljubne količine pri majhni spremembi parametra, od katerega je količina odvisna. Osnovna količina je popolnoma poljubna - lahko je funkcija, funkcional, matrika, itd.). Smerni odvod öf (a)1 splošne funkcije /(a) je definiran kot d öf (a) [h] = aç /(a+ç h), (7.1) kjer je çh sprememba argumenta a, ç pa skalar, s katerim je skaliran h. 5f(a) [h] označuje smerni odvod funkcije /(a) v smeri h. Oznako [h] v nadaljevanju besedila opustimo. Funkcijo / razvijemo v Taylorjevo vrsto po skalarnem parametru ç okrog ç = 0 in dobimo ^členov -i /(a+ç h) = /(a) + y — o /(a+ç h) ç . (7.2) *-^ nI n=l Preden Reissnerjeve kinematične enačbe razvijemo v Taylorjevo vrsto na opisani način, pripravimo prve in druge variacije deformacij (tokrat v komponentni obliki): e = [u + X ) cos (p + (v + Y ) sin (p — 1, = ou cos (p — [u + X ) simp dip + ov simp + (v + Y ) cosip dip, (7.3) o e = — 2 ou smtL> dip — (u + X) cost/? dip + 2 ov cosip dip — [v + Y) simp dip , 1 Smerni odvod in variacija sta označena z istim simbolom (S). Gre namreč za isti količini. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 93 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 7 = — [u + X ) srn

— [v + Y J sin tL> oy>, (7-4) 0 7 = — 2 ou cos ip o

«' sin tL>0 — 0 + 2 tL> f ' cos cp0 — ip2 (v + Y)' sin fL>0) , 72orK = 7LinK + 2 (—^ ^ w'cos ^o + Y>2 (w + X) sin ip0 — 2 (p v' sin fL>0 — 2 (f + Y)' cos tpo) , (7.10) Z oznako ‘2orK’ smo označili kinematično teorijo drugega reda. Točna teorija. Enačbe točne teorije (2.25) poznamo in jih označujemo z oznako ‘Reis’. 94 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 7.1.2 Stopnje interpolacije in integracije Vsi elementi so razviti za štiri interpolacijske stopnje. Najnižja stopnja interpolacije je linearna, kjer imamo le dve interpolacijski funkciji in dve vozlišči na elementu. Najvišja stopnja interpolacije je četrtega reda in ima pet vozlišč oz. interpolacijskih točk. Interpolacijske točke so vedno enakomerno razporejene po dolžini elementa. V programskem okolju AceGen/AceFEM so vsi elementi - topološko gledano - daljice, ki jih podajamo le z začetnim in končnim vozliščem. Ostala vozlišča so obravnavana kot notranja. Oznake interpolacijskih stopenj so: ‘II’ - linearna, '12’ - kvadratična, '13’ - kubična in '14’ - interpolacija četrtega reda. Integracija po dolžini končnega elementa ni poljubna, temveč je vnaprej določeno, da je stopnja integracije reducirana za togostni del enačb in polna za vztrajnostni del. Izbiramo lahko med Lobattovo in Gaussovo integracijo, vendar pa moramo odločitev eksplicitno navesti v programskem okolju AceFEM. Elementi z linearno interpolacijo ne morejo imeti Lobattove integracije prve stopnje, saj ta ne obstaja. 7.1.3 Materialni modeli V elementih lahko upoštevamo tri vrste materialnih modelov - elastičnega, jeklenega in betonskega. Vsi modeli so že bili podrobno opisani. Za njih vpeljemo oznake: ‘El’ - za elastičnega, ‘St’ - za jeklenega in ‘RC - za betonskega. 7.1.4 Cˇasovni integratorji Implementiranih je sedem časovnih integratorjev: Newmark (1959) - New; Hilber, Hughes in Taylor (1977) - HHT; Simo, Tarnow in Doblare (1995) - STD; Ibrahimbegovič in Mamouri (1999) - IbMa; Stander in Stein (1996) - StaSt; šibki s povprečno rotacijsko matriko - Wpov in šibki z midpoint rotacijsko matriko - Wmid. Časovni integratorji so razviti le za točno kinematično teorijo z izjemo New-markovega in HHT, ki sta narejena tudi za linearno kinematično teorijo. 7.1.5 Vrsta analize in oznake elementov Analiza je lahko dinamična (oznaka ‘Dyn’) ali statična (oznaka ‘Sta’). V primeru statične analize so oznake elementov sestavljene iz oznak: 2DB + Sta + oznaka materialnega modela + oznaka kinematične teorije + oznaka stopnje interpolacije. Oznaka ‘2DB’ pove, da gre za dvo-dimenzionalne nosilce (2 dimensional beams), tipična oznaka končnega elementa pa je recimo ‘2DBStaRCLinKI4’, ki označuje element za statično analizo armiranobetonskih konstrukcij z linearno kinematično teorijo in četrto stopnjo interpolacije neznanih količin po elementu. 7.1.6 Uporabne funkcije Za vse elemente so razvite priročne funkcije, ki omogočajo bolj podrobne analize in več možnosti modeliranja. Detajli uporabe funkcij so prikazani v demonstracijskih datotekah na spletu. Vgrajene so funkcije za • Dodajanje členkov: Možno je dodati translacijske členke v globalnih smereh X in Y in momentne členke. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 95 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. • Dodajanje točkovnih mas: V statični analizi te možnosti ni, saj tam maso modeliramo s točkovno silo. • Upoštevanje gravitacijskega pospeška: Določimo lahko poljubno vrednost gravitacijskega pospeška. • Obremenjevanje z linijskimi obtežbami: Podamo lahko konstantne ali linearno spreminjajoče se linijske obtežbe po dolžini posameznih elementov. Obtežba se integrira v reducirani zanki integracije, zato je potrebno biti pozoren, da se ne uporablja elementov s konstantno interpolacijo ('H’) in linearno spreminjajoče se linijske obtežbe. S podajanjem linijske obtežbe lahko modeliramo pospeške temeljnih tal. • Analiza lastnih vrednosti: V vsakem trenutku analize je v statični analizi možen račun lastnih vektorjev (oblik) in lastnih vrednosti in v dinamični analizi račun nihajnih oblik in nihajnih časov. 7.2 Knjižnica končnih elementov Uporaba naprednih orodij za generacijo programske kode končnih elementov nam je omogočila, da smo naredili veliko različnih elementov. Elementi so zbrani v preglednicah 7.1 in 7.2. 7.2.1 Končni elementi za statično analizo V preglednici 7.1 so zbrani elementi za statično analizo. Preglednica 7.1: Knjižnica končnih elementov za statično analizo. Table 7.1: Finite element library for static analysis. Material Stopnja int. Lin. teorija Teorija 2. reda Točna teorija Linearna 2DBStaElLinKIl 2DBStaE12orKI1 2DBStaElReisI1 Elastičen Kvadratična Kubična 2DBStaElKinKI2 2DBStaElLinKI3 2DBStaE12orKI2 2DBStaE12orKI3 2DBStaElReisI2 2DBStaElReisI3 4. reda 2DBStaElLinKI4 2DBStaE12orKI4 2DBStaElReisI4 Linearna 2DBStaStLinKIl 2DBStaSt2orKIl 2DBStaStReisIl Jeklen Kvadratična Kubična 2DBStaStKinKI2 2DBStaStLinKI3 2DBStaSt2orKI2 2DBStaSt2orKI3 2DBStaStReisI2 2DBStaStReisI3 4. reda 2DBStaStLinKI4 2DBStaSt2orKI4 2DBStaStReisI4 Linearna 2DBStaRCLinKI1 2DBStaRC2orKI1 2DBStaRCReisIl AB Kvadratična 2DBStaRCKinKI2 2DBStaRC2orKI2 2DBStaRCReisI2 Kubična 2DBStaRCLinKI3 2DBStaRC2orKI3 2DBStaRCReisI3 4. reda 2DBStaRCLinKI4 2DBStaRC2orKI4 2DBStaRCReisI4 96 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 7.2.2 Koncm elementi za dinamično analizo V preglednici 7.2 so zbrani elementi za dinamično analizo. Preglednica 7.2: Knjižnica končnih elementov za dinamično analizo. Table 7.2: Finite element library for dynamic analysis. Alg. Stopnja int. Elastičen mat. Jeklen mat. AB Linearna 2DBDynElReisI1New 2DBDynStReisI1New 2DBDynRCReisI1New Kvadratična 2DBDynElReisI2New 2DBDynStReisI2New 2DBDynRCReisI2New Kubična 2DBDynElReisI3New 2DBDynStReisI3New 2DBDynRCReisI3New 4. reda 2DBDynElReisI4New 2DBDynStReisI4New 2DBDynRCReisI4New Linearna 2DBDynElReisI1HHT 2DBDynStReisI1HHT 2DBDynRCReisIlHHT Kvadratična 2DBDynElReisI2HHT 2DBDynStReisI2HHT 2DBDynRCReisI2HHT Kubična 2DBDynElReisI3HHT 2DBDynStReisI3HHT 2DBDynRCReisI3HHT 4. reda 2DBDynElReisI4HHT 2DBDynStReisI4HHT 2DBDynRCReisI4HHT Linearna 2DBDynElReisIlStaSt 2DBDynStReisIlStaSt 2DBDynRCReisIlStaSt Kvadratična 2DBDynElReisI2StaSt 2DBDynStReisI2StaSt 2DBDynRCReisI2StaSt Kubična 2DBDynElReisI3StaSt 2DBDynStReisI3StaSt 2DBDynRCReisI3StaSt 4. reda 2DBDynElReisI4StaSt 2DBDynStReisI4StaSt 2DBDynRCReisI4StaSt Linearna 2DBDynElReisIlSTD 2DBDynStReisIlSTD 2DBDynRCReisIlSTD Kvadratična 2DBDynElReisI2STD 2DBDynStReisI2STD 2DBDynRCReisI2STD Kubična 2DBDynElReisI3STD 2DBDynStReisI3STD 2DBDynRCReisI3STD 4. reda 2DBDynElReisI4STD 2DBDynStReisI4STD 2DBDynRCReisI4STD Linearna 2DBDynElReisI1IbMa 2DBDynStReisI1IbMa 2DBDynRCReisI1IbMa Kvadratična 2DBDynElReisI2IbMa 2DBDynStReisI2IbMa 2DBDynRCReisI2IbMa Kubična 2DBDynElReisI3IbMa 2DBDynStReisI3IbMa 2DBDynRCReisI3IbMa 4. reda 2DBDynElReisI4IbMa 2DBDynStReisI4IbMa 2DBDynRCReisI4IbMa Linearna 2DBDynElReisI1Wavg 2DBDynStReisI1Wavg 2DBDynRCReisil Wavg Wavg Kvadratična 2DBDynElReisI2Wavg 2DBDynStReisI2Wavg 2DBDynRCReisI2Wavg Kubična 2DBDynElReisI3Wavg 2DBDynStReisI3Wavg 2DBDynRCReisI3Wavg 4. reda 2DBDynElReisI4Wavg 2DBDynStReisI4Wavg 2DBDynRCReisI4Wavg Linearna 2DBDynElReisI1Wmid 2DBDynStReisIlWmid 2DBDynRCReisI1Wmid Kvadratična 2DBDynElReisI2Wmid 2DBDynStReisI2Wmid 2DBDynRCReisI2Wmid Kubična 2DBDynElReisI3Wmid 2DBDynStReisI3Wmid 2DBDynRCReisI3Wmid 4. reda 2DBDynElReisI4Wmid 2DBDynStReisI4Wmid 2DBDynRCReisI4Wmid Gams, M. 2008. Geom. točna dinamična analiza elastičnih in AB okvirnih konstrukcij 97 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. 8 ZAKLJUČKI V disertaciji smo se ukvarjali s štirimi področji mehanike kinematično točnih ravninskih nosilcev: (i) z izpeljavo in aplikacijo končnih elementov z interpolacijo osne, strižne in upogibne deformacije kot osnovnih neznanih količin v dinamiki nosilcev; (ii) s časovno integracijskimi shemami v dinamiki končnih elementov z interpolacijo pomikov s poudarkom na integratorjih, ki ohranjajo mehansko energijo in količine gibanja. Pri tem smo razvili tudi novo časovno integracijsko shemo za geometrijsko točne ravninske nosilce; (iii) z optimizacijo v dinamiki kinematično točnih elastičnih nosilcev in (iv) z uporabo kinematično točnih nosilcev v dinamiki armiranobetonskih konstrukcij. V sklopu numeričnih in teoretičnih raziskav dinamike nosilcev z interpolacijo deformacijskih količin smo prišli do naslednjih ugotovitev: • Razvili in predstavili smo osnove pristopa k reševanju tako zastavljenega problema in ga preverili z računskimi primeri. V njem smo upoštevali osne, strižne in upogibne deformacije. Ugotovili smo, da tak pristop ne blokira in daje zelo natančne rezultate, vendar pa je s stališča računskega časa zaradi gnezdenja integralov nekonkurenčen in časovno nestabilen pri uporabi Newmarkove časovno integracijske sheme. • Učinkovitost osnovne formulacije je možno numerično izboljšati s posebnim pristopom k nu-merični integraciji vgnezdenih funkcij. Z izboljšavo se računski čas v primerjavi z osnovno formulacijo skrajša približno za faktor 10. Izboljšana formulacija ima podobno število numeričnih operacij na iteracijo Newton-Raphsonovega postopka kot standardne formulacije s pomiki in je v tem pogledu povsem konkurenčna drugim formulacijam za geometrijsko točno dinamično analizo. • Izboljšana formulacija je zastavljena tako, da v njej ni več eksplicitne integracije po dolžini elementa pri vsakem računu integralov. V primerjavi s prvotno formulacijo, kjer je stopnjo integracije potrebno določiti s poskušanjem, je v izboljšani formulaciji stopnja integracije že vsebovana v formulaciji in ni uporabniški podatek, kar je prednost. • Časovno stabilnost na deformacijah osnovane formulacije je moč izboljšati z izbiro ustrezne časovne integracijske sheme, ki pa bi morala biti brezpogojno stabilna. Izpeljava take sheme je še odprto vprašanje. • Prisotnost konsistenčnih enačb v formulaciji nam zagotavlja konsistentnost količin in odpornost na blokiranje. Ker so osnovne spremenljivke deformacije, lahko deformacije in napetosti v poljubni točki nosilca računamo z enako natančnostjo. To je drugače kot pri formulacijah s pomiki, kjer je deformacije potrebno računati z ekstrapolacijo deformacij, izračunanih v Barlowih točkah (Barlow (1976)). 98 Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. V sklopu raziskav časovnih integracij za končne elemente, ki temeljijo na interpolaciji kinematičnih količin, ugotavljamo: • Energijske integratorje je moč izpeljati iz pogoja o ohranitvi energije, kar sta pokazala že Jelenič in Crisfield (2002). Njuna izpeljava velja za formulacije, ki uporabljajo kinematične enačbe v osnovni (krepki) obliki. V disertaciji smo pokazali, da lahko isto izhodišče služi za izpeljavo integratorjev, ki temeljijo na kinematičnih enačbah v šibki obliki. • V pristopu, ki izhaja iz zahteve po ohranitvi energije, ravnotežne enačbe niso osnovne, temveč izpeljane in imajo v splošnem drugačno obliko kot klasične. Izpeljano ravnotežje imenujemo ‘energijsko ravnotežje’ in je klasičnemu enako le v limiti, ko gre inkrement zasukov proti nič (?ip —> 0). • Z ustrezno interpolacijo rotacijske matrike lahko dosežemo, da je energijsko ravnotežje enako klasičnemu, toda le v primeru integratorja, ki uporablja šibke kinematične enačbe. • Novo razvita integratorja, ki temeljita na uporabi šibkih kinematičnih enačb, ohranjata vse količine gibanja (mehansko energijo, gibalno in vrtilno količino). Ohranjanje količin se da analitično dokazati, bilo je preverjeno tudi z numeričnim primerom. • Za vse analizirane integratorje je kvaliteta ohranjanja količin gibanja, če se te ohranjajo, neodvisna od velikosti časovnega koraka in stopnje krajevne interpolacije. • Integratoci z ohranjevalnimi lastnostmi so veliko bolj stabilni od klasičnih, nekonzervativnih, kar je posebej očitno v primeru togih enačb. • Klasični model dušenja deformacij po Armem in Petöczu lahko nadgradimo s pametnim algoritmom, ki vklaplja in izklaplja dušenje po potrebi. Identifikacija situacij za vklop in izklop dušenja poteka na osnovi spreminjanja predznakov sprememb deformacij in prekoračitve praga velikosti. S tem se zmanjša količina disipirane energije in, kar je prav tako pomembno, vpliv na osnovne nihajne oblike. Pri aplikaciji optimizacijskih postopkov na dinamiko kinematično točnih nosilcev smo ugotovili: • Vpeljava posrednih projektnih spremenljivk za opis oblike preko Bézierjevih teles, ki je tako uspešna pri optimizaciji statičnih problemov, se je obnesla tudi v dinamiki. • Potrdilo se je znano dejstvo iz teorije optimizacije, da lahko pogoje, ki so podani na zveznem območju, skoraj enakovredno nadomestimo z manjšim številom diskretnih pogojev. • Optimizacijski postopek v zahtevnejših primerih ni bil sposoben skonvergirati k sprejemljivi rešitvi že v prvem poskusu. Optimizacijo pa smo lahko vedno uspešno izpeljali do konca, če smo problem prevedli na več preprostejših faz. Izkušnje in pripravljenost poskušati ostajata del uspešne optimizacije. • Predstavili smo algoritem za občutljivostno analizo v primeru uporabe šibkega časovnega integratorja s povprečno rotacijsko matriko. Reševanje optimizacijskih problemov je potrdilo velik potencial take formulacije pri optimiranju dinamičnih problemov. Ob izpeljavi in aplikaciji numeričnega postopka za dinamično analizo kinematično točnih ravninskih armiranobetonskih konstrukcij smo ugotovili: Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 99 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. • Formulacija predstavlja novost na tem področju, saj so formulacije s točnimi kinematičnimi zvezami za analizo armiranobetonskih konstrukcij redke. • V analizah smo uporabili več različnih konzervativnih shem, toda pokazalo se je, da so za obravnavane primere enakovredne. S tega stališča ugotavljamo, daje izbira algoritma za račun obravnavanih AB konstrukcij drugotnega pomena. Sklep ne velja za splošne AB konstrukcije. 8.1 Prispevki k znanosti Pomembnejše izvirne prispevke k tehničnim znanostim predstavljajo naslednji novi numerični postopki v kinematično točni dinamični analizi elastičnih in armiranobetonskih konstrukcij, izpeljani v disertaciji: • Izpeljali in implementirali smo novo formulacijo za numerično reševanje enačb dinamike kinematično točnih ravninskih nosilcev, ki za osnovne neznane količine uporablja osno, strižno in upo-gibno deformacijo. Ključne lastnosti formulacije so odpornost na blokiranje, natančnost, konsistentnost napetostnih rezultant v prerezu, diferencialno-algebrajski sistem vodilnih enačb in pojav gnezdenih integralov. • Izpeljali smo numerični postopek s katerim se izognemo zaporednemu računu notranjih integralov v primeru vgnezdenih integralov. Ideja postopka temelji na interpolaciji argumentov gnezdenih integralov, s čimer se občutno zmanjša število računskih operacij pri vrednotenju integralov. Postopek smo uporabili v formulaciji z interpolacijo deformacij in dosegli občutno skrajšanje računskih časov. • Razvili smo nov časovni integrator, ki za kinematično točne, linearno elastične ravninske nosilce, ki točno ohranja mehansko energijo in količine gibanja in to neodvisno od velikosti časovnega koraka in krajevne interpolacije. Posebnost integratorja je uporaba kinematičnih enačb v šibki obliki. Za integrator smo ohranjanje količin dokazali analitično in preverili na numeričnih primerih. • Integrator smo dopolnili tudi z umetnim dušenjem za stabilizacijo rešitve, ki s selektivnim vklap-ljanjem in izklapljanjem numeričnega dušenja dobro uduši visoke oblike nihanja, medtem ko skoraj v celoti ohranja vpliv osnovnih nihajnih oblik. • Izpeljali in implementirali smo tudi občutljivostno analizo v dinamiki kinematično točnih elastičnih konstrukcij s šibkim časovnim integratorjem s povprečno rotacijsko matriko. Z razvitim programskim orodjem smo nato optimirali zanimive probleme iz dinamike konstrukcij in strojnih manipulator] ev. • Izpeljane in implementirane formulacije smo nato priredili tudi za račun kinematično točnih armiranobetonskih konstrukcij. 100 Gams, M. 2008. Geom. tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. VIRI Ademar, B.N., White, D.W. 2005. Displacement, flexibility, and mixed beam-column finite element formulations for distributed plasticity analysis. J. Struct. Eng. 131, 12: 1811–1819. Armero, F., Peto˝cz, E. 1999. A new dissipative time-stepping algorithm for frictional contact problems: formulation and analysis. Comput. Methods Appl. Mech. Eng. 179: 151–178. Armero, F., Romero, I. 2001a. On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. Part I: Low-order methods for two model problems and nonlinear elasto-dynamics. Comput. Methods Appl. Mech. Eng. 190: 2603–2649. Armero, F., Romero, I. 2001b. On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. Part II: Second-order methods. Comput. Methods Appl. Mech. Eng. 190, 6783–6824. Armero, F., Romero, I. 2003. Energy-dissipative momentum-conserving time-stepping algorithms for the dynamics of nonlinear Cosserat rods. Comput. Mech. 31: 3–26. Arora, J.S., Dutta, A. 1997. Explicit and implicit methods for design sensitivity analysis of a nonlinear structures under dynamic loads. Appl. Mech. Rev. Trans. ASME 50: S11–S19. Barlow, J. 1976. Optimal stress locations in finite element models. Int. J. Numer. Methods Eng. 10: 243–251. Bathe, K.J. 1996. Finite element procedures. New Jersey. Prentice-Hall, Inc: 1037 str. Bauchau, O.A., Damilano, G., Theron, N.J. 1995. Numerical integration of non-linear elastic multi-body systems. Int. J. Numer. Methods Eng. 38: 2737–2751. Bauchau, O.A., Theron, N.J. 1996a. Energy decaying scheme for nonlinear elastic multi-body systems. Comp. Struct. 59: 317–331. Bauchau, O.A., Theron, N.J. 1996b. Energy decaying scheme for nonlinear beam model. Comput. Methods Appl. Mech. Eng. 134: 37–56. Bauchau, O.A., Bottasso, C.L., Trainelli, L. 2003. Robust integration schemes for flexible multibody systems. Comput. Methods Appl. Mech. Eng. 192: 395–420. Belytschko, T., Schoeberle, D.F. 1975. On the unconditional stability of an implicit algorithm for nonlinear structural dynamics. J. Appl. Mech. 42: 865–869. Betsch, P., Steinmann, P. 2000. Inherently energy conserving time finite elements for classical mechanics. J. Comput. Phys. 160: 88–116. Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij 101 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Betsch, P., Steinmann, P. 2003. Constrained dynamics of geometrically exact beams. Comput. Mech. 31: 49–59. Bottasso, C.L., Borri, M. 1998. Integrating finite rotations. Comput. Methods Appl. Mech. Eng. 164: 307–331. Bottasso, C.L., Borri, M., Trainelli, L. 2001. Integration of elastic multibody systems by invariant conserving/dissipating algorithms. Part II: Numerical schemes and applications. Comput. Methods Appl. Mech. Eng. 190: 3701–3733. Bottasso, C.L., Bauchau, O.A., Choi, J.Y. 2002. An energy decaying scheme for nonlinear dynamics of shells. Comput. Methods Appl. Mech. Eng. 191: 3099–3121. Bratina, S. 2003. Odziv armiranobetonskih linijskih konstrukcij na pozˇarno obtezˇbo. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo: 159 str. Bratina, S., Saje, M., Planinc, I. 2004. On materially and geometrically non-linear analysis of reinforced concrete planar frames. Int. J. Solids Struct. 41: 7181–7207. Bui, Q.V. 2003. Energy dissipative time finite elements for classical mechanics. Comput. Methods Appl. Mech. Eng. 192: 2925–2947. Cardoso, J.B., Arora, J.S. 1992. Design sensitivity analysis of nonlinear dynamic response of structural and mechanical systems. Struct. Optim. 4: 37–46. Cho, S., Choi, K.K. 2000a. Design sensitivity analysis and optimization of non-linear transient dynamics. Part I – sizing design. Int. J. Numer. Methods Eng. 48: 351–373. Cho, S., Choi, K.K. 2000b. Design sensitivity analysis and optimization of non-linear transient dynamics. Part II – configuration design. Int. J. Numer. Methods Eng. 48: 375–399. Chung, J., Hulbert, G.M. 1993. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-? method. ASME J. Appl. Mech. 60: 371–375. Como, M., De Stefano, M., Ramasco, R. 2003. Effects of column axial force – bending moment interaction on inelastic seismic response of steel frames. Earth. Eng. Struct. Dyn. 32: 1833–1852. Corte, G.D., Barecchia, E., Mazzolani, F.M. 2006. Seismic upgrading of RC buildings by FRP: Full-scale tests of a real structure. J. Mat. Civil Engrg. ASCE 18, 5: 659–669. Crisfield, M.A., Shi, J. 1994. A co-rotational element/time-integration strategy for non-linear dynamics. Int. J. Numer. Methods Eng. 37: 1897–1913. Crisfield, M.A., Shi, J. 1996. An energy conserving co-rotational procedure for non-linear dynamics with finite elements. Nonlinear Dynamics 9: 37–52. Crisfield, M.A., Galvanetto, U., Jelenic´, G. 1997. Dynamics of 3-D co-rotational beams. Comput. Mech. 20: 507–519. Cˇ as, B. 2004. Nelinearna analiza kompozitnih nosilcev z uposˇtevanjem zdrsa med sloji. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo: 136 str. Della Corte, G., Barecchia, E., Mazzolani, F.M. 2005. Seismic upgrading of RC buildings by FRP: full-scale tests of a real structure. J. Mater. Civ. Eng. 18, 5: 659–669. Dombrowski, S. von. 2002. Analysis of large flexible body deformation in multibody systems using absolute coordinates. Multibody System Dynamics 8: 409–432. 102 Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Dufva, K.E., Sopanen, J.T., Mikkola, A.M. 2005. A two-dimensional shear deformable beam element based on the absolute nodal coordinate formulation. J. Sound Vibration 280: 719–738. Erlicher, S., Bonaventura, L., Bursi, O.S. 2002. The analysis of the generalized-? method for the nonlinear dynamic problems. Comput. Mech. 28: 83–104. Ehrlich, D., Armero, F. 2005. Finite element methods for the analysis of softening plastic hinges in beams and frames. Comput. Mech. 35: 237–264. Fajfar, P. 1984. Dinamika gradbenih konstrukcij. Ljubljana, Univerza Edvarda Kardelja, FAGG: 550 str. Filippou, F.C., Popov, E.P., Bertero, V.V. 1983. Effects of bond deterioration on hysteretic behavior of reinforced concrete joints. EERC Report 83-19. Berkeley, Earthquake engineering research center Berkeley. Fung, T.C. 2003. Numerical dissipation in time-step integration algorithms for structural dynamic analysis. Prog. Struct. Engrg. Mater. 5: 167–180. Gams, M., Saje, M., Srpcˇicˇ, S., Planinc, I. 2007a. Finite element dynamic analysis of geometrically exact beams. Comput. Struct. 85, 17/18: 1409–1419. Gams, M., Planinc, I., Saje, M. 2007b. The strain-based finite elements in multibody dynamics. J. Sound Vibration 305, 1–2: 194–210. Gams, M., Planinc, I., Saje, M. 2007c. Energy conserving time integration scheme for geometrically exact beams. Comput. Methods Appl. Mech. Engrg 196, 17–20: 2117–2129. Gams, M., Planinc, I., Saje, M. 2007d. A heuristic viscosity-type dissipation for high frequency oscillation damping in time integration algorithms. Comput. Mech. 41, 1: 17–29. Ge´radin, M., Cardona, A. 2001. Flexible Multibody Dynamics. A Finite Element Approach. England, John Wiley and Sons Ltd: 340 str. Goicolea, J.M., Orden, J.C.G. 2002. Quadratic and higher-order constraints in energy-conserving formulations of flexible multibody systems. Multibody Sys. Dyn. 7: 3–29. Graham, E. 2001. A survey of current time-stepping algorithms used in non-linear, dynamic structural analysis. Internal report. London, Department of Aeronautics, Imperial College of Science, Technology & Medicine. Harl, B., Kegl, M. 2005. Optimization of space trusses. J. Mech. Engrg. 51, 9: 570–588. Hilber, H.M., Hughes, T.J.R., Taylor, R.L. 1977. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Eng. Struct. Dyn. 5: 282–292. Hjelmstad, K.D., Taciroglu, E. 2003. Mixed variational methods for finite element analysis of geometrically non-linear, inelastic Bernoulli-Euler beams. Commun. Num. Methods Engrg. 19: 809–832. Hozjan, T., Turk, G., Srpcˇicˇ, S. 2007. Fire analysis of steel frames with the use of artifical neural networks. J. Constr. Steel Res. 63, 10: 1396–1403. Hsiao, K.M., Jang, J.Y. 1991. Dynamic analysis of planar flexible mechanisms by co-rotational formulation. Comput. Methods Appl. Mech. Engrg 87: 1–14. Hughes, T.J.R. 1992. The Finite Element Method, Prentice Hall, Englewood Cliffs, N.J. Ibrahimbegovic´, A., Mamouri, S. 1999. Nonlinear dynamics of flexible beams in planar motion: formulation and time-stepping scheme for stiff problems. Comput. Struct. 70: 1–22. Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij 103 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Ibrahimbegovic´, A., Mamouri, S., Taylor, R.L., Chen, A.J. 2000. Finite element method in dynamics of flexible multibody systems: Modeling of holonomic constraints and energy conserving integration schemes. Multibody Sys. Dyn. 4: 195–223. Ibrahimbegovic´, A., Mamouri, S. 2002. Energy conserving/decaying implicit time-stepping scheme for three-dimensional beams undergoing finite rotations. Comput. Methods Appl. Mech. Eng. 191: 4241– 4258. Iura, M., Atluri, S.N. 1995. Dynamic analysis of planar flexible beams with finite rotations by using inertial and rotating frames. Comput. Struct. 55: 453–462. Jelenic´, G. 1990. Dinamika ravninskih hiperelasticˇnih nosilcev pri neomejenih deformacijah. Magistrska naloga. Ljubljana, Univerza v Ljubljani, Fakulteta za arhitekturo, gradbenisˇtvo in geodezijo. Jelenic´, G., Saje, M. 1994. Finite deformations of linear elastic space beams. Z. angew. Math. Mech. 74, 4:298–300. Jelenic´, G., Saje, M. 1995. A kinematically exact space finite strain beam model – finite element formulation by generalized virtual work principle. Comput. Methods Appl. Mech. Eng. 120:131–161. Jelenic´, G., Crisfield, M.A. 2002. Stability and convergence characteristics of conserving algorithms for dynamics of 3D rods. Internal report, Department of Aeronautics, Imperial College of Science, Technology & Medicine, London, UK. Jelenic´, G., Crisfield, M.A. 2001. Dynamic analysis of 3D beams with joints in presence of large rotations. Comput. Methods Appl. Mech. Eng. 190: 4195–4230. Kane, T.R., Likins, P.W., Levinson, D.A. 1983. Spacecraft Dynamics. McGraw-Hill, New York. Kang, B.S., Park, G.J., Arora, J.S. 2006. A review of optimization of structures subjected to transient loads. Struct. Multidisc. Optim. 31, 2: 81–95. Kante, P. 2005. Potresna ranljivost armiranobetonskih sten, Doktorska disertacija, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Karsan, I.D., Jirsa, J.O. 1969. Behavior of concrete under compressive loading. J. Struct. Division ASCE 95, ST12: 2543–2563. Kegl, M. 2000. Shape optimal design of structures: an efficient shape representation concept. Int. J. Numer. Methods Eng. 49: 1571–1588. Kegl, M., Butinar, B., Kegl, B. 2002. An efficient gradient-based optimization algorithm for mechanical systems. Comm. Numer. Meth. Engrg. 18: 363–371. Kegl, M. 2005. Parameterization based shape optimization: theory and practical implementation aspects. Eng. Comput. 22, 5–6: 646–663. Kegl, M., Brank, B. 2006. Shape optimization of truss-stiffened shell structures with variable thickness. Comput. Methods Appl. Mech. Eng. 195: 2611–2634. Kegl, M. 2007. Program ‘iGO’. http://fs-server.uni-mb.si/si/inst/im/loms/KeglM-si.html. Korelc, J. 1997. Automatic generation of finite-element code by simultaneous optimization of expressions. Theor. Comput. Sci. 187: 231–248. Korelc, J. 2007. AceGen, AceFEM and AceShare. http://www.fgg.uni-lj.si/Symech/Flyer.pdf 104 Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Krauberger, N., Saje, M., Planinc, I., Bratina, S. 2007. Exact buckling load of a restrained RC column. Struct. Eng. Mech. 27, 3: 293–310. Kuhl, D., Ramm, E. 1996. Constraint energy momentum algorithm and its application to non-linear dynamics of shells. Comput. Methods Appl. Mech. Eng. 136: 293–315. Kuhl, D., Crisfield, M.A. 1999. Energy-conserving and decaying algorithms in non-linear structural dynamics. Int. J. Numer. Methods Eng. 45: 569–599. Kuhl, D., Ramm, E, 1999. Generalized energy-momentum method for non-linear adaptive shell dynamics. Int. J. Numer. Methods Engrg. 178: 343–366. Kulkarni, M., Noor, A.K. 1995. Sensitivity analysis of the nonlinear dynamic viscoplastic response of 2-D structures with respect to material parameters. Int. J. Numer. Methods Eng. 38: 183–198. Mander, J.B., Priestley, M.J.N., Park, R. 1988. Theoretical stress-strain model for confined concrete. J. Struct. Eng. ASCE 114, 8: 1804–1825. Markovicˇ, M. 2006. Mejna nosilnost AB precˇnih prerezov pri dvojno ekscentricˇni osni obremenitvi. Diplomska naloga. Ljubljana. Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Menegotto, M., Pinto, P.E. 1973. Method of analysis for cyclically loaded reinforced concrete plane frames including changes in geometry and non-elastic behavior of elements under combined normal force and bending, Proceedings, IABSE Symposium on resistance and ultimate deformability of structures acted on by well defined repeated loads: 15–22. Nuenhofer, A., Filippou, F.C. 1998. Geometrically nonlinear flexibility-based frame finite element. J. Struct. Eng. 124, 6: 704–711. Newmark, N.M. 1959. A method of computation for structural dynamics. ASCE J. Engrg. Mech. Division 85: 67–94. Omar, M.A., Shabana, A.A. 2001. A two-dimensional shear deformable beam for large rotation and deformation problems. J. Sound Vibration 243, 3: 565–576. Opensees. 2006. Open System for Earthquake Engineering Simulation. http://opensees.berkeley.edu/index.php. Paz, M., Leigh, W. 2004. Structural Dynamics. New York, Springer, 5th ed. Pedersen, C.B.W. 2004. Crashworthiness design of transient frame structures using topology optimization. Comput. Methods Appl. Mech. Eng. 193, 6–8: 653–678. Pedersen, N.L., Nielsen, A.K., 2003. Optimization of practical trusses with constraints on eigenfrequen-cies, displacements, stresses and buckling. Struct. Multidisc. Optim. 25, 5–6: 436–445. Petrangeli, M., Ciampi, V. 1997. Equilibrium based iterative solutions for the non-linear beam problem. Int. J. Numer. Methods Engrg. 40: 423–437. Planinc, I. 1998. Racˇun kriticˇnih tocˇk konstrukcij s kvadraticˇno konvergentnimi metodami. Doktorska disertacija, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Planinc, I., Saje, M., Cˇ as, B. 2001. On local stability condition in planar beam finite element. Struct. Eng. Mech. 12: 507–526. Popovics, S. 1973. A numerical approach to the complete stress strain curve for concrete. Cement and concrete research 3, 5: 583–599. Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij 105 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Prohitech 2004. Project n. INCO-CT-2004-509119. Earthquake protection of historical buildings by reversible mixed technologies. Project acronym: Prohitech. WP8: NUMERICAL ANALYSES – Preliminary study. Reissner, E. 1972. On one-dimensional finite-strain beam theory: the plane problem. J. Appl. Math. Physics (ZAMP) 23: 795–804. Romero, I., Armero, F. 2002. Numerical integration of the stiff dynamics of geometrically exact shells: an energy-dissipative momentum-conserving scheme. Int. J. Numer. Methods Engrg. 54: 1043–1086. Rong, T.-Y., Lu, A.-Q. 2003. Generalized mixed variational principles and solutions of ill-conditioned problems in computational mechanics. Part II: Shear locking. Comp. Methods Appl. Mech. Engrg. 192: 4981–5000. Rozman, M., Fajfar, P. 2006. Earthquake protection of historical buildings by reversible mixed technologies. Project report. WP9: numerical analyses. Preliminary study. Saje, M. 1990. A variational principle for finite planar deformation of straight slender elastic beams. Int. J. Solids Struct. 26, 8: 887-900. Saje, M. 1991. Finite-element formulation of finite planar deformation of curved elastic beams. Comput. Struct. 39, 3–4: 327-337. Saje, M., Jelenic´, G. 1994. Finite element formulation of hyperelastic plane frames subjected to nonconservative loads. Comput. Struct. 50, 2: 177–189. Saje, M., Planinc, I., Turk, G., Vratanar, B. 1997. A kinematically exact finite element formulation of planar elastic-plastic beam. Comp. Methods Appl. Mech. Engrg. 144: 125–151. Sansour, C., Sansour, J., Wriggers, P. 1996. A finite element approach to the chaotic motion of geometrically exact rods undergoing in-plane deformations. Nonlinear Dyn. 11: 189–212. Sansour, C., Wriggers, P., Sansour, J. 1997. Nonlinear dynamics of shells: Theory, finite element formulation and integration schemes. Nonlinear Dyn. 13: 279–305. Sansour, C., Wagner, W., Wriggers, P., Sansour, J. 2002. An energy-momentum integration scheme and enhanced strain finite elements for the non-linear dynamics of shells. Int. J. Non-Lin. Mech. 37: 951–966. Sansour, C., Wriggers, P., Sansour, J. 2003. A finite element post-processed Galerkin method for dimensional reduction in the non-linear dynamics of solids. Application to shells. Comput. Mech. 32: 104–114. Sansour, C., Wriggers, P., Sansour, J. 2004. On the design of energy-momentum integration schemes for arbitrary continuum formulations. Applications to classical and chaotic motion of shells. Int. J. Numer. Methods Eng. 60: 2419–2440. Schnabl, S. 2007. Mehanska in pozˇarna analiza kompozitnih nosilcev. Doktorska disertacija. Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Shabana, A.A. 1997. Definition of the slopes and the finite element absolute nodal coordinate formulation. Multibody System Dynamics 1, 3: 339–348. Shabana, A.A. 1998. Dynamics of multibody systems. Second ed. Cambridge. Cambridge University Press: 384 str. 106 Gams, M. 2008. Geometrijsko tocˇna dinamicˇna analiza elasticˇnih in AB okvirnih konstrukcij Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbenisˇtvo in geodezijo, Konstrukcijska smer. Simo, J.C., Vu-Quoc, L. 1986. On the dynamics of flexible beams under large overall motions – the plane case: Part I and Part II. ASME J. Appl. Mech. 53: 849–863. Simo, J.C., Tarnow, N., Wong, K.K. 1992. Exact energy-momentum conserving algorithms and sym-plectic schemes for non-linear dynamics. Comput. Methods Appl. Mech. Eng. 100: 63–116. Simo, J.C., Tarnow, N., Doblare, M. 1995. Non-linear dynamics of three-dimensional rods: exact energy and momentum conserving algorithms. Int. J. Numer. Methods Eng. 38: 1431–1473. Sivaselvan, M.V., Reinhorn, A.M. 2002. Collapse analysis: Large inelastic deformations analysis of planar frames. J. Struct. Eng. 128, 12: 1575–1583. Sousa, L.G., Cardoso, J.B., Valido, A.J. 1997. Optimal cross-section and configuration design of elastic– plastic structures subject to dynamic cyclic loading. Struct. Optim. 13: 112–118. Spacone, E., Filippou, F.C., Taucer, F.F. 1996. Fibre beam-column model for non-linear analysis of R/C frames: Part I. Formulation, Part II. Applications. Earth. Eng. Struct. Dyn. 25: 711–742. Stander, N., Stein, E. 1996. An energy-conserving planar finite beam element for dynamics of flexible mechanisms. Engrg. Comput. 13, 6: 60–85. Stocki, R., Kolandek, K., Jendo, S., Kleiber, M. 2001. Study on discrete optimization techniques in reliability-based optimization of truss structures. Comput. Struct. 79, 22–25: 2235–2247. Stupkiewicz, S. 2001. Approximate response sensitivities for nonlinear problems in explicit dynamic formulation. Struct. Multidisc. Optim. 21: 283–291. Taucer, F.F., Spacone, E., Filippou, F.C. 1991. A fiber beam-column element for seismic response analysis of reinforced concrete structures. EERC Report 91-17, Earthquake engineering research center Berkeley. Toniolo, G. (coordinator) 2007. Final report of the EU Research Project: Seismic Behaviour of Precast Concrete Structures with respect to EC8 (Contract No. G6RD-CT-2002-00857). Vratanar, B. 1998. Geometrijsko tocˇna dinamika ravninskih nosilcev. Doktorska disertacija. Lubljana, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Vratanar, B., Saje, M. 1999. A consistent equilibrium in a cross-section of an elastic-plastic beam. Int. J. Solids Struct. 36, 2: 311–337. Wang, D., Zhang, W.H., Jiang, J.S. 2002a. Combined shape and sizing optimization of truss structures. Comput. Mech. 29, 4–5: 307–312. Wang, D., Zhang, W.H., Jiang, J.S. 2002b. Truss shape optimization with multiple displacement constraints. Comput. Methods Appl. Mech. Eng. 191, 33: 3597–3612. Wang, D., Zhang, W.H., Jiang, J.S. 2004. Truss optimization of shape and sizing with frequency constraints. AIAA J. 42, 3: 622–630. Wood, W.L., Bossak, M., Zienkiewicz, O.C. 1980. An alpha modification of the Newmark’s method. Int. J. Numer. Methods Engrg 15: 1562–1566. Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z. 2005. The Finite Element Method: Its Basis and Fundamentals, Sixth edition, Elsevier Butterworth-Heinemann. Zupan, D. 2001. Teorija prostorskih nosilcev, osnovana na ukrivljenostih. Magistrska naloga. Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Gams, M. 2008. Geometrijsko točna dinamična analiza elastičnih in AB okvirnih konstrukcij 107 Doktorska disertacija. Ljubljana, UL, Fakulteta za gradbeništvo in geodezijo, Konstrukcijska smer. Zupan, D. 2003. Rotacijsko invariantne deformacijske kolicˇine v geometrijsko tocˇni teoriji prostorskih nosilcev. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbenisˇtvo in geodezijo. Zupan, D., Saje, M. 2003. A new finite element formulation of three-dimensional beam theory based on interpolation of curvature. Comput. Model. Eng. Sci. 4, 2: 301–318. PRILOGE PRILOGA A Gams, M., Saje, M., Srpcˇicˇ, S., Planinc, I. 2007a. Finite element dynamic analysis of geometrically exact beams. Computers and structures 85, 17/18: 1409–1419. PRILOGA B Gams, M., Planinc, I., Saje, M. 2007b. The strain-based finite elements in multibody dynamics. Journal of sound and vibration 305, 1-2: 194–210. PRILOGA C Gams, M., Planinc, I., Saje, M. 2007c. Energy conserving time integration scheme for geometrically exact beams. Computer methods in applied mechanics and engineering 196, 17–20: 2117–2129. PRILOGA D Gams, M., Planinc, I., Saje, M. 2007d. A heuristic viscosity-type dissipation for high frequency oscillation damping in time integration algorithms. Computational mechanics 41, 1: 17–29. PRILOGA E Gams, M., Saje, M., Planinc, I., Kegl, M. 2007. Optimization of geometrically exact planar beams in dynamics. Neobjavljeno.