Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo l^W.TM h iTFfi 1 ill iTTl In ill 11 ill m| I TTTTltTTlTTTTi PODIPLOMSKI ŠTUDIJ GRADBENIŠTVA DOKTORSKI ŠTUDIJ Kandidatka: EVA ZUPAN, prof. matem., specialistka DINAMIKA PROSTORSKIH NOSILCEV S KVATERNIONSKO PARAMETRIZACIJO ROTACIJ Doktorska disertacija štev.: 201 DYNAMICS OF SPATIAL BEAM STRUCTURES USING QUATERNION PARAMETRIZATION OF ROTATIONS Doctoral thesis No.: 201 Temo doktorske disertacije je odobrila Komisija za doktorski študij, po pooblastilu s 25. seje Senata Univerze v Ljubljani, dne 24. junija 2008 in imenovala mentorja prof.dr. Mirana Sajeta in somentorja izr.prof.dr. Igorja Planinca. Ljubljana, 7. maj 2010 Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo l^W.TM h iTFfi 1 ill iTTl In iih i ill m| n;TTiiTii:ii!i Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi prof.dr. Miran Saje izr.prof.dr. Igor Planine izr.prof.dr. Gordan Jelenič (Univerza na Rijeki, GF) prof.dr. Miha Boltežar (UL, FS) je imenoval Senat Fakultete za gradbeništvo in geodezijo na 16. redni seji dne 26. marca 2008. Komisijo za oceno doktorske disertacije v sestavi prof.dr. Goran Turk izr.prof.dr. Gordan Jelenič (Univerza na Rijeki, GF) prof.dr. Miha Boltežar (UL, FS) je imenoval Senat Fakultete za gradbeništvo in geodezijo na 8. redni seji dne 24. februarja 2010. Komisijo za zagovor doktorske disertacije v sestavi prof.dr. Matjaž Mikoš, dekan UL FGG prof.dr. Goran Turk izr.prof.dr. Gordan Jelenič (Univerza na Rijeki, GF) prof.dr. Miha Boltežar (UL, FS) je imenoval Senat Fakultete za gradbeništvo in geodezijo na 10. redni seji dne 21. aprila 2010. Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo l^W.TM h iTFfi 1 ill iTTl In iih i ill m| n;TTiiTii:ii!i IZJAVA O AVTORSTVU Podpisana EVA ZUPAN, prof. matem., specialistka, izjavljam, da sem avtorica doktorske disertacije z naslovom: »DINAMIKA PROSTORSKIH NOSILCEV S KVATERNIONSKO PARAMETRIZACIJO ROTACIJ« Ljubljana, 7. maj 2010 (podpis) BIBLIOGRAFSKO-DOKUMENTACIJSKA STRAN IN IZVLEČEK UDK Avtor: Mentor: Somentor: Naslov: Obseg in oprema: Ključne besede: 512:624.074(043.3) Eva Zupan prof.dr. Miran Saje izr.prof.dr. Igor Planine Dinamika prostorskih nosilcev s kvaternionsko parametrizacijo rotacij 163 str., 12 pregl., 36 sl., 400 en. Reissnerjev nosilec, prostorski nosilci, rotacija, rotacijski kvaternion, kvater-nionska algebra, statika, dinamika, potujoci delec Izvlecek V doktorski disertaciji izpeljemo formulacijo za reševanje enacb prostorskih nosilcev po kinematicno tocni Reissnerjevi teoriji, osnovano na kvaternionski parametrizaciji rotacij in pripadajocimi operacijami iz kvaternionske algebre. Tak zapis enacb ohranimo tudi pri izpeljavi linearizirane oblike in v numericni implementaciji. Predstavljena numericna metoda je osnovana na metodi koncnih elementov. Osnovne neznanke so kinematicne kolicine - pomiki in rotacijski kvaternioni. Enacbe po kraju diskretiziramo s kolokacijsko metodo, neznanke pa interpoliramo. Tocke interpolacije in kolokacije poenotimo; v notranjih tockah zahtevamo ujemanje odvodov konstitucijskih notranjih in ravnoteznih sil, v robnih tockah pa zadošcamo ravnoteznim pogojem. Druzina izpeljanih koncnih elementov omogoca poljubno zacetno lego, obremenjevanje s staticno in z dinamicno obtezbo in vgradnjo dokaj splošnega nelinearnega materiala. Diskretizirane nelinearne enacbe za staticno analizo rešimo z Newtonovo metodo. Pomemben algoritem tega reševanja je konsistenten postopek za dodajanje linearnih popravkov rotacijskih kvaterni-onov in ostalih neaditivnih kolicšin. Predstavljena numericšna metoda nima strizšnega blokiranja, izbrana interpolacija rotacijskih kvaternionov pa zagotavlja objektivnost deformacij. (Časovno integracijo di-namicnih enacb izvedemo na dva nacina: s klasicno metodo Runge-Kutta in s kombinacijo casovne diskretizacije enacšb s posplosšeno metodo Newmark, prirejeno za rotacijske kvaternione, in Newtonove metode za reševanje dobljenih nelinearnih algebrajskih enacb. Prva metoda casovne integracije ne zahteva linearizacije enacšb, potrebno je le prevesti sistem diferencialnih enacšb prostorskega nosilca drugega reda na sistem diferencialnih enacb prvega reda, medtem ko druga metoda upošteva neaditivno naravo rotacij in ostalih z rotacijami povezanih kolicin ter je ze namenjena reševanju diferencialnih enacb drugega reda. Pri numericni formulaciji za dinamicno analizo prostorskih nosilcev z Newmarkovo integracijo po casu dodatno obravnavamo vpliv delca z maso, ki se giblje vzdolz tezišcne osi nosilca. Gibalno enacbo delca rešujemo socasno in povezano z enacbami za dinamicno analizo nosilcev. Pri tem lahko upoštevamo poljubno zacetno geometrijo in poljubne robne in zacetne pogoje delca in konstrukcije. BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION Title: UDC Author: Supervisor: Co-supervisor: 512:624.074(043.3) Eva Zupan Prof. Miran Saje Asoc. Prof. Igor Planinc Dynamics of spatial beam structures using quaternion parametrization of rotations Notes: Key words: 163 p., 12 tab., 36 fig., 400 eq. Reissner's beam, space beam, rotation, rotational quaternion, quaternion algebra, statics, dynamics, moving mass Summary In the present thesis a new formulation of the three dimensional geometrically exact Reissner beam theory based on quaternion algebra is presented. Rotations are fully replaced by rotational quaternions. The same approach is used for both linearization and the numerical implementation. The finite element method for solving nonlinear differential equations is presented where the kinematic quantities displacements and rotational quaternions are choosen as the primary unknown functions. Equations are discretizied by collocation mehod and primary unknowns are interpolated. The interpolation and the collocation points coincide. At internal nodes the formulation satisfies the weak equality of equilibrium and constitutive internal forces, while the equilibrium equations are satisfied at the edges of the beam. The present family of finite elements for static and dynamic analysis allows arbitrary initial geometry and a general form of the constitutive law describing the material of the beam. The static equations are solved using the Newton iteration procedure for solving nonlinear equations. For the update procedure of the rotational quaternion and other non-aditive quantities special formulae are derived. Proposed numerical method does not suffer from shear locking phenomenon and the choosen interpolation of rotational quaternions is objective with respect to conservation of strains with rigid body motion. For time integration two different approaches are presented: the classic Runge-Kutta method and the combination of the extension of the Newmark discretization algoritm to the quaternion parametrization of rotation and the Newton iteration method. For both approaches the algebraic boundary condition are modified to the differential form with respect to time. The first integration method does not require the linearized equations however it has been developed for the sistem of first-order differential equations and the beam equations must be written this way while the second approach properly considers the nonlinearity of the rotations and related quantities and is directly applicable for solving second-order differential equations. Finally the moving mass problem is studied. Mass moves along the centroid line of a beam according to the motion equation of the moving mass which is coupled with the system for governing equations of spatial beams under dynamic loading. The geometry of the beam as well as the boudary and initial conditions of the beam and the mass are arbitrary. Zahvale Na mojo odločitev, da nadaljujem s študijem in na končne rezultate mojega dela je vplivalo mnogo ljudi. Najprej se zelim zahvaliti mojemu mentorju Miranu za zaupanje, moznost in odprtost za nove ideje ter nenazadnje za strokovno pomoč. Ostalim kolegom na katedri in podiplomskim študentom sem hvalezna za pozitivno delavno vzdušje. Hani, Paskalu in Jonatanu sem dolzna vso priznanje za razvoj moje pridnosti, kije temeljni pogoj za uspešno raziskovalno delo. Frančiški in Janezu se zahvaljujem za pomoč pri skrbi za otroke. Najbolj pa je na moje odločitve in delo vplival Dejan; za vso spodbudo, strokovne razgovore in potrpezšljivost se mu iskreno zahvaljujem. Kazalo 1 Uvod 1 1.1 Pregled stanja................................................................................1 1.2 Opis dela....................................................................................3 2 Definicija linijskega nosilca 5 2.1 Začetna lega..................................................................................5 2.2 Trenutna lega................................................................................6 3 Prostorske rotacije 9 3.1 Definicija rotacij in njihove osnovne lastnosti..............................................9 3.1.1 Referenčna in pomična baza........................................................10 3.2 Parametrizacija rotacij z rotacijskim vektorjem ............................................12 3.3 Odvodi in variacije rotacijskih kolicin......................................................13 4 Parametrizacija rotacij z rotacijskim kvaternionom 16 4.1 Kvaternionska algebra in rotacijski kvaternion..............................................16 4.2 Matricni zapis kvaternionov in operatorjev..................................................20 4.2.1 Referencna in pomicna baza kvaternionskega prostora............................22 4.2.2 Zapisi rotacijskih kvaternionov v razlicnih bazah..................................23 4.3 Odvod in variacija ............................................................................24 5 Enacbe prostorskega nosilca za statiko 29 5.1 Enacbe s parametrizacijo rotacije z rotacijskim vektorjem..................................29 5.2 Enacbe s kvaternionsko parametrizacijo rotacij ............................................32 5.2.1 Izbira in priprava glavnih in robnih enacb..........................................33 5.3 Reševanje staticnih kvaternionskih enacb ..................................................34 5.3.1 Diskretizacija rešitve in enacb......................................................35 5.3.2 Newtonova iterativna metoda........................................................36 5.3.3 Priprave na linearizačijo enačb: linearizačija posameznih količin..................37 5.3.4 Linearizačija enačb nosilča..........................................................42 5.3.5 Numerična integračija..............................................................43 5.3.6 Postopek dodajanja linearnih popravkov in popravljanje ostalih količin ..........44 5.4 Numerični testi..............................................................................50 6 Enačbe prostorskega nosilca za dinamiko 60 6.1 Gibalna količina in njen odvod po času......................................................60 6.2 Vrtilna količina in njen odvod po času......................................................61 6.3 Dopolnitev enačb prostorskega nosilča za dinamiko........................................69 6.3.1 Izrek o gibalni količšini na delu nosilča ..............................................71 6.3.2 Izrek o vrtilni količšini na delu nosilča ..............................................72 6.4 Kvaternionska oblika enačšb ..................................................................75 6.4.1 Osnovne neznanke ..................................................................76 6.5 Uvod v numerično reševanje enačb..........................................................78 6.6 Diskretizačija rešitve in enačb ..............................................................79 6.7 Izbira osnovnih enačšb in integračije po čšasu ................................................80 6.8 Naravni robni pogoji v diferenčialni obliki..................................................80 7 ReSevanje dinamičnih enaCb z metodami druZine Runge-Kutta 83 7.1 Metode druzine Runge-Kutta................................................................83 7.2 Priprava glavnih in robnih enačšb ............................................................85 7.3 Numeričšni testi ..............................................................................88 8 ReSevanje dinamičnih enačb s posplošeno metodo Newmark 96 8.1 Izpeljava posplošene metode Newmark za kvaternionske enačbe..........................96 8.1.1 Začetni priblizki za nov časovni korak: prediktor.................102 8.2 Priprava enačb za numerični račun.............................103 8.3 Linearizačija hitrosti in pospeškov.............................104 8.4 Linearizačija posameznih členov enačb ..........................106 8.5 Numerična integračija po kraju...............................108 8.6 Lokalna napaka in prilagajanje časovnega koraka metode.................108 8.7 Numerični testi.......................................109 9 Potujoci delec 114 9.1 Izbira krivocrtne baze in ukrivljenost tezišcne osi.....................114 9.2 Hitrost in pospešek delca..................................116 9.3 Gibalne enacbe delca....................................117 9.4 Upoštevanje vpliva gibajocega se delca na nosilec.....................119 9.5 Numericno reševanje....................................121 9.6 Linearizacija ................................................................................122 9.6.1 Linearizacija krajevnega vektorja in njegovih odvodov..............122 9.6.2 Linearizacija baznih vektorjev...........................123 9.6.3 Linearizacija ukrivljenosti.............................125 9.6.4 Linearizacija Kontaktnih sil.............................125 9.6.5 Linearizacija gibalne enacbe delca.........................126 9.7 Dolocanje lege delca v konstrukciji.............................129 9.8 Numericni primeri .....................................129 9.8.1 Delec na klancu...................................130 9.8.2 Delec na prostolezšecšem nosilcu ....................................................131 9.8.3 Delec na previsnem nosilcu ........................................................135 9.8.4 Vodni tobogani...................................136 10 Zakljucek 141 11 Povzetek 143 12 Summary 146 Literatura 148 Priloga 154 Dodatki 155 A Odvod, integral in variacija ..................................................................155 A1 Odvod in variacija operatorja............................155 A2 Odvod in variacija vektorske funkcije skalarnega argumenta...........157 B Eksponentna preslikava kvaternionskega argumenta....................157 B1 Eksponentna oblika rotacijskega kvaterniona...................158 B2 Odvod in variacija rotacijskega kvaterniona....................159 B3 Inverz in njegova linearizacija...........................160 C Prevedba kvaternionske rotacijske matrike v klasicno Rodriguesovo formulo.......161 Seznam slik 2.1 Nosilec v zacetni in trenutni legi............................................................7 2.1 The initial and current position of the beam................................................7 5.1 Geometrijska interpretacija Newtonove iteracije............................................37 5.1 Geometrical interpretation of the Newton iteration ........................................37 5.2 Konzola s precno silo v prostem krajišcu ..................................................51 5.2 Cantilever under a free-end vertical force ..................................................51 5.3 Precni pomik prostega krajišca v odvisnosti od parametra GL2/Eh2......................51 5.3 Vertical displacement at the free-end vs parameter GL2/Eh2 ............................51 5.4 Primerjava precnega pomika prostega krajišca v odvisnosti od parametra GL2/Eh2 s klasicnim Timošenkovim nosilcem s striznim blokiranjem ................................52 5.4 Vertical displacement at the free-end vs parameter GL2/Eh2; comparison with the classical Timoshenko beam with shear locking problem ........................................52 5.5 Konzola z momentom ......................................................................52 5.5 Cantilever under free-end moment..........................................................52 5.6 Zviti nosilec ................................................................................53 5.6 Twisted beam................................................................................53 5.7 Bocna zvrnitev pravokotnega okvirja ......................................................55 5.7 Lateral buckling of a right-angle frame ....................................................55 5.8 Ukrivljena konzola..........................................................................56 5.8 Cantilever 45° bend ........................................................................56 5.9 Konzola, upognjena v spiralo ..............................................................57 5.9 Cantilever bent to the helical form..........................................................57 5.10 Koncna deformirana oblika konzole, upognjene v spiralo..................................58 5.10 Cantilever: deformed shape at the final load step ..........................................58 5.11 Konzola, upognjena v spiralo: odvisnost pomik-sila........................................58 5.11 Load-displacement curve for the helical beam....................... 58 6.1 ObteZba na delu nosilca .................................. 71 6.1 Loading at infinitesimal part of a beam .......................... 71 7.1 Geometrija in obteZba upogibnega nosilca......................... 89 7.1 Force driven flexible beam: geometry and loading data.................. 89 7.2 Pomiki prostega krajišca v odvisnosti od casa za primer upogibnega nosilca. 10 linearnih elementov; primerjava z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) . 90 7.2 Free-end displacement of force driven flexible beam. 10 linear elements; comparison with Ibrahimbegovic and al Mikdad (1998) (dashed-line)................. 90 7.3 Pomiki prostega krajišca v odvisnosti od casa za primer upogibnega nosilca. 2 linearna elementa; primerjava z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) . . 90 7.3 Free-end displacement of force driven flexible beam. 2 linear elements; comparison with Ibrahimbegovic and al Mikdad (dashed-line)........................ 90 7.4 Geometrija in obtezba kolenaste konzole ......................... 91 7.4 Right-angle cantilever beam: geometry and loading data ................. 91 7.5 Odvisnost pomikov od casa za koleno (A) in za prosto krajišce (B) kolenaste konzole; sive crte prikazujejo rezultate Sima in Vu-Quoca (1988) za primer dveh in desetih elementov ........................................... 92 7.5 The right-angle cantilever beam: time histories of elbow (A) and tip (B) displacements; comparison with Simo and Vu-Quoc (1988) (grey lines) employing 2 and 10 elements . 92 7.6 Geometrija in obtezba nepodprtega nosilca ........................ 93 7.6 Free-free flexible beam: geometry and loading data.................... 93 7.7 Ravninske projekcije poteka gibanja nepodprtega nosilca do casa t = 5 ......... 94 7.7 Projections of the deformed shapes on the coordinate planes for the free-free flexible beam to time t = 5 ..................................... 94 7.8 Aksonometricni prikaz zaporedja deformiranih leg nepodprtega nosilca pri znacilnih casih na casovnem obmocju [2.5,15.5] .......................... 95 7.8 Free-free flexible beam. Perspective view of deformed shapes on time interval [2.5,15.5] 95 8.1 Pomiki prostega krajišca v odvisnosti od casa za primer spiralnega gibanja upogibnega nosilca z uporabo metode Newmark; 10 linearnih elementov; primerjava z rešitvami po metodi Runge-Kutta in z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) . 110 8.1 Free-end displacement of force driven flexible beam using the Newmark time integration method; 10 linear elements; comparison with results obtained by the Runge-Kutta time integration method and with results by Ibrahimbegovic and al Mikdad (1998) (dashed-line) .............................................110 8.2 Pomiki prostega krajišca v odvisnosti od casa za primer spiralnega gibanja upogibnega nosilca z uporabo metode Newmark; 2 linearna elementa; primerjava z rešitvami po metodi Runge-Kutta in z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) . 110 8.2 Free-end displacement of force driven flexible beam using the Newmark time integration method; 2 linear elements; comparison with results obtained by the Runge-Kutta time integration method and with results by Ibrahimbegovic and al Mikdad (1998) (dashed- line).............................................110 8.3 Pomiki kolena (tocka A) in prostega krajišca (tocka B) v odvisnosti od casa po metodi Newmark za primer kolenaste konzole; primerjava z resšitvami po metodi Runge-Kutta in z rešitvijo iz literature [Simo, Vu-Quoc, 1988] .....................111 8.3 The right-angle cantilever beam: time histories of elbow (A) and tip (B) displacements employing the Newmark integration scheme and comparison with results obtained by the Runga-Kutte method and [Simo, Vu-Quoc, 1988].....................111 8.4 Ravninske projekcije poteka gibanja nepodprtega nosilca do casa t = 5 .........112 8.4 Projections of the deformed shapes on the coordinate planes for the free-free flexible beam to time t = 5.....................................112 8.5 Aksonometricni prikaz zaporedja deformiranih leg nepodprtega nosilca pri znacilnih casih na casovnem obmocju [2.5,15.5] ..........................113 8.5 Free-free flexible beam: perspective view of deformed shapes on time interval [2.5,15.5] 113 9.1 (a) Sistem sil na nosilcu; (b) sistem sil na delcu......................118 9.1 (a) Forces acting on a beam; (b) forces on a particle....................118 9.2 Zacetna lega za primer delca na klancu ..........................130 9.2 Particle on a slope: the initial configuration ........................ 130 9.3 Zacšetna lega za primer delca na prostolezšecšem nosilcu .................. 131 9.3 Simply supported beam: the initial configuration ..................... 131 9.4 Primerjava precnih pomikov sredine prostolezecega nosilca s potujoco silo in potujocim delcem; vijolicna pikcasta crta prikazuje numericno rešitev ob upoštevanju zacetno deformiranega nosilca, v ostalih primerih je privzet zacšetno raven nosilec ......... 133 9.4 Lateral displacements at the mid point of simply supported beam when moving particle/load is applied; the purple dotted-line represents the numerical solution considering initially deformed shape due to gravity, other numerical results consider initially straight beam ............................................ 133 9.5 Primerjava rezultatov za drugi primer prostolezecega nosilca; zelena crtkana crta prikazuje rezultate v primeru, ko upoštevamo zacetno deformirano lego nosilca, pri ostalih rezultatih predpostavljamo zacšetno raven nosilec ..................... 134 9.5 The lateral displacements at the mid point of simply supported beam (case 2) when moving particle/load is applied; the green dashed-line represents the numerical solution considering initially deformed shape due to gravity, other numerical results consider initially straight beam ..................................... 134 9.6 Primerjava rezultatov za primer prehoda delča prek previsnega nosilča..........135 9.6 Moving partičle on čantilever beam: čomparison of approximative solutions.......135 9.7 Prečni prerez visokih vodnih toboganov; geometrijski podatki ravnega tobogana .... 137 9.7 Cross-sečtion of high water slides; geometričal data for straight water slide .......137 9.8 Potek osnih deformačij po obodu prereza na deformirani legi za primer statične analize ravnega tobogana; pomiki konstrukčije so pomnozeni s faktorjem 30, barvna skala označuje deformačije v pročentih .............................138 9.8 Statič analysis of straight water slide: axial deformations on the deformed čonfiguration; displačements are multiplied by fačtor 30, deformations are multiplied by fačtor 100 (in per čent) ..........................................138 9.9 Ogrinjača poteka osnih deformačij na obodu prereza na deformirani legi za primer dinamične analize ravnega tobogana; pomiki konstrukčije so pomnozeni s faktorjem 30, barvna skala označuje deformačije v pročentih ......................138 9.9 Dynamič analysis of straight water slide: evelope of axial deformations on the deformed čonfiguration; displačements are multiplied by fačtor 30, deformations are multiplied by fačtor 100 (in per čent)...................................138 9.10 Začetna geometrija ukrivljenega tobogana.........................139 9.10 Initially čurved water slide: the initial čonfiguration ...................139 Seznam preglednic 4.1 Zapisi rotacijskih kvaternionov v pripadajocih bazah ......................................23 4.1 Rotational quaternions with respect to natural bases........................................23 5.1 Popravki rotacijskega vektorja in njegovih odvodov........................................46 5.1 Update of rotational vector and its derivatives..............................................46 5.2 Pomiki in zasuki prostega krajišca konzole z momentno obtezbo..........................53 5.2 Free-end displacements and rotations of cantilever beam under free-end moment . . . . 53 5.3 Pomiki prostega krajišca po linearni teoriji za zviti nosilec ................................54 5.3 Free-end displacements of twisted beam; linearized theory................................54 5.4 Kriticni upogibni moment Mcr prostolezece podprtega pravokotnega okvirja ............55 5.4 Critical moment Mcr of the simply supported right-angle frame ..........................55 5.5 Komponente r\,r2, r3 krajevnega vektorja prostega krajišca pri upogibu ukrivljene kon- zole ..........................................................................................56 5.5 Components r\,r2, r3 of the free-end position vector of the cantilever 45° bend under out-of-plane force ..........................................................................56 5.6 Primerjava pomikov prostega krajišca v odvisnosti od števila korakov nalaganja obtezbe za primer ukrivljene konzole................................................................57 5.6 Cantilever 45° bend: free-end displacements as a function of the number of load steps . 57 8.1 Posplošena implicitna Newmarkova shema za rotacijski kvaternion ...........100 8.1 Generalized Newmark implicit time-stepping algorithm for rotational quaternion .... 100 8.2 Aproksimacije za popravljanje dinamicnih kolicin ....................102 8.2 Approximate formulae for dynamic quantities.......................102 9.1 Najvišje vrednosti deformacij in pomikov za primer ravnega tobogana..........137 9.1 Straight water slide: maximum values of deformations and displacements........137 9.2 Najvišje vrednosti deformacij in pomikov za primer ukrivljenega tobogana .......140 9.2 Curved water slide: maximum values of deformations and displacements ................140 9.3 Ukrivljen tobogan: pomiki, hitrosti in pospesški delča za primera deformabilnega in to- gegatobogana ....................................... 140 9.3 Displačements, veločities and aččelerations for a deformable and rigid čurved water slide 140 Seznam opomb Opomba 1 ................................................................................... 6 Opomba 2 ................................................................................... 8 Opomba 3 ................................................................................... 8 Opomba 4 .................................................................................. 10 Opomba 5 .................................................................................. 12 Opomba 6 .................................................................................. 24 Opomba 7 .................................................................................. 28 Opomba 8 .................................................................................. 33 Opomba 9 .................................................................................. 35 Opomba 10 ................................................................................40 Opomba 11 ................................................................................ 45 Opomba 12 ................................................................................47 Opomba 13 ................................................................................ 66 Opomba 14 ................................................................................69 Opomba 15 ................................................................................78 Opomba 16 ................................................................................88 Opomba 17 ................................................................................92 Opomba 18 ...............................................................................124 Opomba 19 ...............................................................................160 1 Uvod 1.1 Pregled stanja Konstrukčije v gradbeništvu in na drugih inzenirskih področjih sestavljajo trirazsezni objekti, ki imajo pogosto eno ali dve dimenziji bistveno večji od preostalih. Take konstrukčijske dele imenujemo ploskovni (plošče, lupine) in linijski (nosilči) elementi. Kadar pri matematičnem modeliranju takih elementov upoštevamo dimenzijske posebnosti tako, da neznanke redučiramo na osnovno ploskev lupine oziroma na teziščno os nosilča, se prostorska in časovna zahtevnost numeričnih izračunov bistveno zmanjšata. Razvoj numeričnih metod na teh področjih je še vedno zivahen, saj tako poenostavljeni modeli v nasprotju s polnimi prostorskimi elementi prinašajo učinkovite in še vedno natančne numerične algoritme. Za dovolj natančno analizo konstrukčij pa je pomembno v model zajeti geometrijsko in materialno nelinearnost ter tudi druge dejavnike, ki močno vplivajo na dinamični odziv konstrukčij. Eden takih je vpliv premikajočših se teles po konstrukčiji. V primeru lahkih konstrukčij, kot so polimerni visoki vodni tobogani, ali v primeru zelo visokih hitrosti teles, na primer prehod hitrega vlaka prek jeklenega zelezniškega mostu, lahko premikajoče se telo povzroči dinamični odziv konstrukčije, ki ga s preprostimi analizami ne moremo ustrezno opisati in ga lahko zato močno podčenimo. Na področju teorij linijskih nosilčev večina raziskovalčev uporablja enega od znanih linijskih modelov nosilča. Preprost Euler-Bernoullijev model, ki zanemari strizne deformačije, predvideva linearno kine-matiko nosilča. Timošenkov model linijskega nosilča vključuje tudi strizne deformačije, vendar je v osnovi prav tako kinematično linearen, s posplošitvami pa tudi višjih redov. Reissnerjev model nosilča je kombinačija točšne (nelinearne) kinematike in Timosšenkovega modela nosilča. Mnoge sodobne metode nosilčev izhajajo iz točnih kinematičnih enačb, razlikujejo pa se v numerični implementačiji. Prvi, ki je Reisnerjev model prostorskega nosilča uspešno numerično implementiral, je bil Simo (1995,1996). Kasneje so podoben pristop uporabili mnogi avtorji, tako za statiko kot za dinamiko. S tem področjem so se ukvarjali: Simo in sodelavči (1988, 1991, 1995), Jelenič in Saje (1995), Bottasso in Borri (1997), Ibrahimbegovič in al Mikdad (1998), Jelenič in Crisfield (1999a), Romero in Armero (2002), D. Zupan in Saje (2003), Mata in sodelavči (2007) ter Ghosh in Roy (2008). Vse omenjene formulačije izhajajo iz Reissnerjeve teorije prostorskih nosilčev [Reissner, 1981] in večinoma za osnovne neznanke problema uporabijo kinematične količine: pomike in rotačijske vektorje, razen Jelenič in Saje, ki uporabita zgolj rotačijske vektorje ter D. Zupan in Saje, ki za osnovne neznanke izbereta deformačijske količine. Rota-čijski vektorji so vektorska parametrizačija rotačijskega operatorja. Zaradi narave mehanskih problemov se rotačijam tezko izognemo. Zaradi zahtevnosti matematičnih struktur, kijih tvorijo rotačijski operatorji in njihove vektorske parametrizačije, se s problematiko uspešne parametrizačije zasukov ukvarjajo številni avtorji s širšega področja mehanike, na primer [Argyris, 1982; Atluri, Cazzani, 1995; Geradin, Rixen, 1995; Ibrahimbegovič et al., 1995; Zupan, 2003]. Omenimo, da zelo uveljavljena parametrizačija prostorskih rotačij z rotačijskim vektorjem ni enolična in vsebuje singularnost. Avtorja MčRobie in Lasenby (1999) sta v novejšem casu na podrocju dinamike prostorskih nosilcev predstavila nesingularno parametrizacijo rotacij z uporabo Cliffordovih algeber. Pri tem se ne omejita le na rotacije v trirazseznem prostoru, temvec njuna parametrizacija velja za rotacije v Evklidskem prostoru poljubnih dimenzij. Pri omejitvi na trirazsezni prostor je parametrizacija s Cliffordovo algebro enakovredna parametrizaciji z algebro kvaternionov. Rotacijski kvaternioni, glej [Poreous, 1995; Ward, 1997], v nobeni od prej naštetih teorij prostorskih nosilcev ne nastopajo kot osnovna spremenljivka, se pa uporabljajo za stabilizacijo numericnih procesov v nekaterih podalgoritmih na primer Spurrierov algoritem [Spurrier, 1978], metoda za ohranjanje ortogonalnosti rotacijske matrike [Simo, Wong, 1991] in za posredno interpolacijo rotacijskih vektorjev - sfericno interpolacijo kvaternionov (SLERP, glej [Shoemake, 1985]) uporabita Ghosh in Roy (2008), da dosezeta deformacijsko objektiven element, ki temelji na vektorski parametrizaciji rotacij; podoben postopek interpolacije primerja Romero (2004) z drugimi pristopi. Nasprotno pa se rotacijski kvaternioni kot osnovna neznanka problema in parametrizacija prostorskih rotacij pogosto uporabljajo na drugih raziskovalnih podrocjih, kjer so sistemi enacb preprostejši, na primer na podrocju racunalniške grafike [Shoemake, 1985], molekularne dinamike [Martys, Mountain, 1999] in na podrocju rotacije togih tridimenzionalnih teles v prostoru [Johnson et al., 2007]. Za casovno integracijo sistema diferencialnih enacb, ki opisujejo dinamiko prostorskih nosilcev, uporabljajo avtorji Simo in sodelavci (1988, 1991, 1995), Bottasso in Borri (1997), Ibrahimbegovic in al Mikdad (1998), Jelenic in Crisfield (1999a) ter Romero in Armero (2002) sredinske ali 'midpoint' integratorje (implicitna Eulerjeva metoda casovne integracije), ter osnovne ali modificirane Newmarkove metode za integracijo sistema diferencialnih enacb drugega reda. Bottasso samostojno (1997) in skupaj z Borrijem (1998) uporabita tudi sicer splošno uveljavljene metode druzine Runge-Kutta (glej npr. [Evans, 1995; Gerald, Wheatley, 1994]). Zivahen razvoj na podrocju dinamike prostorskih nosilcev poteka na podocju integratorjev, ki ohranjajo ali celo reducirajo invariante dinamicnih sistemov, to je energijo in/ali gibalno in vrtilno kolicino, saj je to eden od nacinov za stabilizacijo numericnega procesa na daljšem casovnem intervalu. Raziskave na to temo najdemo v [Bauchau et al., 1995; Bauchau, Theron, 1996; Bottasso, Borri, 1997; Gams et al. 2007; Simo, Wong, 1991; Simo et al., 1995] in drugih. Integratorji, ki ohranjajo invariante dinamicnih sistemov, so omejeni na metode nizkega reda za reševanje sistemov diferencialnih enacb, kar posledicno ne zagotav-lja visoke natancšnosti elementov. Neaditivnost, neenolicšnost ter sigularnost parametriza-cije prostorskih rotacij z rotacijskim vektorjem in tudi drugi še ne dovolj raziskani vzroki na podrocju dinamike prostorskih nosilcev povzrocajo številne probleme numericnih formulacij, kot so blokiranje koncnih elementov, nestabilnost racunskih postopkov pri uporabi metod višjega reda za reševanje diferencialnih enacb, odvisnost od poti nalaganja obtezbe, neobjektivnost deformacij klasicnih elementov dinamike prostorskih nosilcev (po definiciji Crisfielda in Jelenica (1991)). Objektivnost sta rešila Crisfield in Jelenic [Crisfield, Jelenic, 1999; Jelenic, Crisfield, 1999a] z zelo zahtevno interpolacijo zasukov. Objektivna je tudi formulacija, kot stajo predstavila Bottasso in Borri in uporaba linearne sfericne interpolacije za rotacijske kvaternione, ki stajo v teoriji nosilcev prva implementirala Ghosh in Roy (2008). Vsi trije omenjeni pristopi, ki ohranjajo objektivnost deformacij, upoštevajo naravo rotacij in so medsebojno enakovredni, kot navaja [Romero, 2004]. Zahtevnost implementacije take interpolacije je bistveno nizja ob izbiri rotacijskega kvaterniona za parametrizacijo rotacije in SLERPa. Objektivnost deformacij pa lahko dosezemo tudi prek interpolacije vseh komponent rotacijske matrike, ortogonalnost rotacijskih matrik pa upoštevamo z ortogonalizacijo z interpolacijo dobljenih (v splošnem ne-ortogonalnih) matrik; enakovredno parametrizacijo, ki ohranja objektivnost deformacij, ob socasni manjši casovni in prostorski zahtevnosti, dosezšemo tudi s kombinacijo neodvisne interpolacije sštirih kvaternionskih parametrov in z normalizacijo rotacijskega kvaterniona v vmesnih, interpoliranih tockah [Romero, 2004]. Reševanje problema gibanja delca po nosilcu je staro ze vec kot 150 let, izvira pa iz zelje po analizi zelezniških mostov. Iz starejših virov omenimo le najpomembnejša. G. G. Stokes je ze leta 1849 raz- pravljal o reševanju diferencialnih enacb, ki vplivajo na prelom zelezniških mostov, S. Timošenko pa je v clanku iz leta 1927 racunal vibracije mostov. Objave v mednarodno priznanih revijah v zadnjem desetletju kazejo, da se nekaj skupin raziskovalcev po svetu še vedno ukvarja s to problematiko, vendar le na preprostih linijskih modelih ravninskih nosilcev. Najpogosteje za model nosilca uporabijo geometrijsko linearni Euler-Bernoullijev ravninski model nosilca s preprostim podpiranjem (prostolezeci nosilec), kot na primer Esmailzadeh in Ghorashi (1995), ki obravnavata lokalno porazdeljeno maso delca; Michaltsos in soavtorji (1996), ki obravnavajo pomik delca po prostolezecem nosilcu s konstantno hitrostjo; Wang (1998), ki analizira pospešeno gibanje delca po nosilcu zaradi vlecne/zaviralne sile; Siddiqui in sodelavci (2000) z vzmetjo modelirajo silo, ki povzroca pomik delca in je sorazmerna pomiku delca; Park in Youm (2001) modelirata pomik vozicka s predpisano hitrostjo. Poenostavljni Euler-Bernoullijev model ravninskega nosilca - sistem togih palic in vzmeti - z vec moznimi podpiranji nosilca pa za analizo prehoda delca s konstantno hitrostjo predstavita Mofid in Akin (1996). Posplošitev na nelinearni Euler-Bernoullijev ravninski linijski nosilec predlagajo Siddiqui in sodelavci (2003), geometrijsko nelinearnost do cetrtega reda pa prek Hamiltonovega principa upoštevajo še Xu in sodelavci (1997). Nekateri raziskovalci obravnavajo vpliv premikajocega se delca tudi na linearnem Timošenkovem ravninskem nosilcu s preprostim podpiranjem, na primer Esmailzadeh in Gorashi (1997) z lokalno porazdeljeno masa delca, Lee (1996a,1996b) pa obravnava delec s konstantno hitrostjo oziroma pospeškom. Yavari s soavtorji (2002) obravnavajo prehod delca s konstantno hitrostjo in z vec moznostmi podpiranja na poenostavljenem Timošenkovem modelu nosilca (sistem togih palic in dvojnih vzmeti). Wu, Whittaker in Cartmell (2000) so opozarjali na pomanjkanje razvoja s podrocja analize prostorskih nosilcev ob obtezbi z gi-bajoco silo ali delcem, ceprav se v praksi pojavlja potreba po takih analizah. Zato so v [Wu et al., 2000] in [Wu et al., 2001] predstavili moznost vgraditve poenostavljenega modela gibanja sile in delca vzdolz prostorskega nosilca v komercialne programe. 1.2 Opis dela V delu predstavljam novo druzino koncnih elementov za staticno in dinamicno analizo prostorskih linijskih konstrukcij, osnovano na kinematicno tocni Reissner-Simovi teoriji prostorskih nosilcev. Linijsko konstrukcijo modeliramo s tezišcno osjo in druzino precnih prerezov. Izbrane enacbe omogocajo poljubno zacetno ukrivljeno in zvito obliko in nepravokoten kot med tezišcno osjo in prerezi nosilca. Kinematicno tocne formulacije temeljijo na tocnih kinematicnih enacbah. Kinematicne kolicine, pomike tezišcne osi in zasuke precnih prerezov, izberemo za osnovne neznanke problema. Zasuk je natancno dolocen z rotacijskim operatorjem, ki ga parametriziramo z rotacijskim kvaternionom. Zaradi take izbire parametrizacije rotacije sistem vodilnih enacb prostorskih nosilcev v celoti zapišemo z operacijami in elementi kvaternionske algebre in se tako popolnoma izognemo uporabi rotacijske matrike in rotacijskega vektorja. Tak zapis uporabimo tudi za numericno implementacijo. Prepletanje rotacijske matrike in rotacijskih parametrov je sicer znacilno za formulacije, osnovane na parametrizaciji rotacij s tremi skalarnimi kolicinami. S prevedbo enacb v kvaternionski zapis se tej dvojnosti popolnoma izognemo in povecšamo racšunsko ucšinkovitost posameznih izrazov. V delu predstavim osnove rotacij in parametri-zacijo rotacij z rotacijskim vektorjem; podrobno obravnavamo algebro kvaternionov in izpeljemo povezave med kvaternionsko in ustaljeno parametrizacijo rotacij z rotacijskim vektorjem. Enacbe rešujemo po metodi koncnih elementov, pri tem za osnovne neznanke izberemo pomike in rotacijske kvaternione. Neznanke diskretiziramo in interpoliramo s polinomi poljubne stopnje; enacbe diskretiziramo z metodo kolokacije. Z izbiro izoparametricne interpolacije za vse štiri komponente rotacijskega kvaterniona ob uporabi normalizacije v vmesnih, interpoliranih tockah, dobimo element, ki ohranja objektivnost deformacij. Tocke diskretizacije in kolokacije poenotimo. V notranjih kolokacijskih tockah zadošcamo šibkim konsistencnim enacbam, kar pomeni enakost odvodov rezultantnih ravnoteznih in materialnih sil in momentov v prerezu; na robu zadošcamo ravnoteznim enacbam. Sistem diskretnih enacb za staticno analizo linijskih konstrukcij rešujemo z Newtonovo metodo. V delu predstavimo postopek za konsistentno upoštevanje linearnih iterativnih popravkov rotacijskih kvaternionov in ostalih rotacijskih kolicin. Rotacijski kvaternioni so prav tako kot rotacije neaditivne kolicine, popravek, dobljen z Newtovovo ite-racijo, je element tangentnega prostora in je v splosšnem neenotski kvaternion. S pravilnim uposštevanjem lastnosti kolicin, ki so povezane z zasuki, izpeljemo postopek za multiplikativen popravek, kije rotacijski kvaternion, in ga lahko dodamo predhodnjemu rotacijskemu kvaternionu z ustreznim kvaternionskim produktom. Linearizacija enacšb, ki vsebujejo rotacije, je ob uporabi rotacijskega kvaterniona bistveno manj zahtevna, kot pri uporabi vektorske parametrizacije rotacij, saj je v kvaternionski algebri rotacija tri-razsezšnega prostora enakovredna produktu dveh linearnih operatorjev realnega sštirirazsezšnega prostora. S tem se popolnoma izognemo Liejevi grupi in smernim odvodom v nelinearnih prostorih. Ustreznost racunskega modela za staticno analizo konstrukcij prikazujemo na znacilnih primerih iz literature. S primeri prikazemo splošnost (poljubna zacetna geomerija, poljubna stopnja interpolacije), natancnost in ucinkovitost predstavljenega algoritma. Numericne simulacije kazejo tudi neobcutljivost postopka na strizšno blokiranje. Pri dinamicni analizi ohranimo model diskretizacije po kraju, uporabimo pa dva zelo razlicna postopka integracije po cšasu. Za potrebe dinamicšne analize robne ravnotezšne enacšbe preoblikujemo v diferencialno obliko tako, da v njih nastopajo osnovne spremenljivke odvajane po casu. S tem se izognemo nu-mericno obcutljivemu reševanju diferencialno-algbrajskega sistema enacb. Za prvi postopek integracije po casu izberemo dve uveljavljeni metodi druzine Runge-Kutta, kot sta implementirani v programskem okolju Matlab. Izbrani pristop za prostorske nosilce nacelno ni ustrezen, ker so klasicne metode Runge-Kutta namenjene reševanju diferencialnih enacb v aditivnih konfiguracijskih prostorih. Kljub temu se za izpeljani model izkazeta ti dve metodi kot dovolj natancni tudi pri zahtevnejših prostorskih primerih. Vzrok temu je ucšinkovitost in robustnost zapisa enacšb v kvaternionski algebri. Prednost tega pristopa je preprostost, saj diskretizacija kolicšin po cšasu in linearizacija enacšb nista potrebni, poleg tega pa metodi prilagajata korak glede na zahtevano lokalno natancšnost. Pri drugem nacšinu integracije enacšb po cšasu pokazšemo, da lahko za predstavljeni koncšni element priredimo tudi integrator, kot sta ga razvila Simo in Vu-Quoc (1988) za klasicne kinematicne spremenljivke. Integrator izhaja iz metode Newmark in je prirejen za uporabo na rotacijskih vektorjih. Integrator prilagodimo izbrani parametrizaciji rotacij z rotacijskim kvaternionom in metodo dopolnimo s preverjanjem lokalne napake. Z obema metodama izracunamo referencne primere iz literature. Primerjava rezultatov po obeh metodah in z drugimi avtorji kaze, da izbira casovnega integratorja mocno vpliva tako na natancnost rezultatov kot na potek reševanja (casovna in prostorska zahtevnost), poleg casovnega integratorja pa ima na numericno reševanje di-namicnih enacb bistven vpliv tudi izbira koncnih elementov. Ob koncu obravnavamo potovanje delca prek prostorskega nosilca in analiziramo njegov odziv. Delec drsi po tezšisšcšni osi nosilca skladno z gibalno enacšbo delca in ne s predpisano hitrostjo ali pospesškom. Rešujemo torej povezan problem, kjer hkrati rešujemo nelinearne parcialne diferencialne enacbe nosilca in delca. Izbranim osnovnim kinematicnim neznankam kot novo neznanko dodamo dolzino opravljene poti delca. Pri tem upoštevamo tocno parametrizacijo tezišcne osi, saj parameter tezišcne osi v deformirani legi ni naravni parameter, in vplive zapisšemo v krivocšrtnem Frenetovem koordinatnem sistemu. Za modeliranje vpliva trenja uporabimo Coulombov zakon. Nosilci, prek katerih potuje delec, imajo poljubne podpore in zacetno obliko. Postopek zaradi primerjave z analiticnimi rezultati in drugimi avtorji poenostavimo še za primer prehoda sile s konstantno hitrostjo. Rezultati kazejo na pomen izbire ustreznega modela, saj lahko s poenostavljanjem modela bistveno podcenimo dinamicen odziv konstrukcije. 2 Definicija linijskega nosilca Nosileč pred začetkom deformiranja ob času t = 0 zavzema začetno lego, ob poljubnem času t > 0 pa trenutno lego. Večina avtorjev prostorskih nosilčev je v svojih začetnih delih izhajala iz nedeformirane začetne lege, na primer Simo samostojno in sodelavči (1985, 1986, 1988), Jelenič in Saje (1995), Jelenič in Crisfield (1999a), Ibrahimbegovič in Mamouri (1999), Romero in Armero (2002), to pomeni, daje teziščna os nosilča daljiča, prerezi pa so vzporedni ter pravokotni na teziščno os nosilča. Nekateri med njimi navajajo načeln postopek, ki omogoča splošno (ukrivljeno, deformirano) začetno lego, vendar je potrebno deformirano začšetno lego obvezno preslikati v nedeformirano, nato pa nadaljevati v trenutno (v statiki v končno) lego. Večina omenjenih avtorjev v kasnejših delih preide na ukrivljeno začetno lego, na primer Simo v delih [Simo et al., 1995], [Simo, Wong, 1991], Crisfield in Jelenič (1999) ter Ibrahimbegovič samostojno in s sodelavči (1995, 1997, 1998). Opis splošne začetne lege preko izbire referenčne ravnine, ki vsebuje referenčen prerez, najdemo v delu Ritto-Corree in Camotima (2002) ter v delih Zupana in Sajeta (2003, 2004, 2006), in predstavlja naravnejši pristop k opisu prostorskih nosilčev. Izhajamo namreč iz poljubne ukrivljene začetne lege in nadaljujemo neposredno v trenutno deformirano lego, izbrani referenčni prerez pa sluzi za opis obeh. Podrobneje opisšimo začšetno in trenutno lego nosilča. 2.1 Začetna lega Naj bo IR3 realni vektorski prostor z izhodiščem O in mirujočo ortonormirano desnosučno bazo ( 0. Lahko jo opišemo na dva nacina. 1. Trenutno lego nosilca v casu t izpeljemo kot posplošitev zacetne lege (2.1) tako, da dodamo casovni parameter. Krajevni vektor poljubne tocke D (x, 6, 6) v trenutni legi nosilca v tem smislu zapisšemo kot (x, ^2,^3,t) = p(x, t) + R (x, t) p® (6,6) , (2.2) kjer r (x, t) za x G [0, L] doloca gladko prostorsko krivuljo r (t), ki opiše deformirano težiščno os nosilca v trenutni legi pri cčasu t, r(t) = {r (x, t), x G [0, L]} . R (x, t) pa so rotacije R (x, t) : IR3 ^ RR3, ki referencne precne prereze zasucejo v prečne prereze nosilca v trenutni legi R (x, t) : Ar 1—^ A (x, t). trenutna lega, t>0 Slika 2.1: Nosilec v zacetni in trenutni legi Figure 2.1: The initial and current position of the beam 2. Drug pristop zapisa trenutne lege nosilca izhaja iz zacetne lege nosilca. Tocke tezišcne osi v trenutni legi izrazimo s pomikom u (x, t) tock osi iz zacetne lege f (x, t) = f[0] (x) + f (x, t), (2.3) skupaj pa opisujejo tečiččno os v trenutni legi rt rt = jr[0] (x) + f (x, t), x G [0, L] j . Druzino precnih prerezov v trenutni legi dobimo z druzino rotacij Z (x, t) precnih prerezov iz zacetne lege: Z (x, t) : IR3 ^ IR3, Z (x, t) : A[0] (x) i—► A (x, t) . Potem lahko krajevni vektor poljubne tocke D (x, {2, {3) nosilca v trenutni legi zapišemo tako: r D (x, {2, {3, t) = f[0] (x) + f (x, t) + Z (x, t) pD[0l ({2, {3) = f[0] (x) + f (x, t) + Z (x, t) R[0] (x) pD ({2, {3). (2.4) Ker sta opisa (2.2) in (2.4) enakovredna, saj za krajevne vektorje velja aditivnost in ker je kompozitum dveh rotacij ponovno rotacija, sledi R (x, t) = Z (x, t) R[0] (x). Linijski prostorski nosilec v trenutni legi v času t zapišemo kot mnozico tock b (t) = {rD (x,6,6,t), x e [0,L] in (&,&) e Ar} . Opomba 2 Tečiččne osi nosilcev so v vseh legah parametrizirane z istim parametrom krivulje, kije pri t = 0 naravni parameter krivulje T[0\ ki opisuje os nosilca v začetni legi: x G [0, L]. Dolčina krivulje r (t) za t > 0 je v spločnem različna od L, torej x ni njen naravni parameter. Opomba 3 Ker prečne prereze v trenutni legi dobimo iz začetnih prečnih prerezov A (x) z rotacijo Z (x, t), se ohrani ravnost (planarnost) prečnih prerezov nosilca. S tem privzemamo znano Bernou-llijevo hipotezo o ravnih prečnih prerezih. Se vedno je dopuščen zasuk normale prereza glede na tangento osi nosilca. 3 Prostorske rotacije V tem poglavju najprej predstavljamo nekaj osnov parametrizačije rotačij z rotačijskim vektorjem, ki ga potrebujemo predvsem za primerjavo z drugimi avtorji. Ker ta parametrizačija ni bistvena za pričujoče delo, je predstavitev rotačijskega vektorja in z njim povezanimi matričnimi operatorji kratka in skopa, bralča pa naprošam, da si po potrebi manjkajoče poišče v premnogih delih o rotačijah, na primer v delih avtorjev Argyris (1982), Agyris in Poterasu (1993), Atluri in Cazzani (1995), Crisfield (1997), Geradin in Rixen (1995) ter Zupan (2003). 3.1 Definicija rotacij in njihove osnovne lastnosti Beseda rotačija ima ze brez matematičnega predznanja nek pomen. V inzenirskih besedilih se skupaj z rotačijo pogosto vpelje še geometrijsko nazorna pojma os rotačije in velikost zasuka rotačije, kar pa ne zadošča za matematično definičijo te operačije. Obstaja več načinov, kako enolično definiramo rotačijo; omenimo vsaj enega. Definicija 1 (Prostorska) rotacija R je linearna preslikava iz vektorskega prostora IR3 v vektorski prostor IR3, ki poljubno desnosučno ortonormirano bazo vektorskega prostora IR3 preslika v desnosučno ortonormirano bazo vektorskega prostora IR3. Kar naenkrat smo vpletli vrsto matematičnih pojmov: linearna preslikava, vektorski prostor, desnosučna ortonormirana baza. Teh osnov v pričujočem delu ne obravnavamo, saj so opisane v večini učbenikov, namenjenih dodiplomskim študentom tehničnih fakultet, na primer v [Krizanič, 1993] in [Krizanič, 1990]. Rotačijaje na prvi pogled 'lepa' preslikava, saj ima čelo vrsto prijetnih lastnosti: • kompozitum dveh rotačij je ponovno rotačija; • premiče preslika v premiče, ravnine v ravnine in like v skladne like; • ohranja orientačijo, dolzine, kote, ploščine likov in prostornine teles. Vendar se izkaze, da je rotačija zelo zahtevna preslikava. Predstavlja osnovno znanje pri formulačiji Reissner-Simovega modela prostorskega nosilča. Za konkretno rabo pa abstraktni, brezkoordinatni opis rotačije ni najbolj primeren, zato vpeljemo tudi baze trirazseznega evklidskega prostora in pripadajoče rotačijske matrike, o čšemer govorimo v naslednjem razdelku. 3.1.1 Referenčna in pomična baza Standardni mehanski opis gibanja delca po prostoru vsebuje zapis v dveh bazah. Referenčno prostorsko bazo Ba = (pi, g2, p3) trirazseznega evklidskega vektorskega prostora z izhodišcem v mirujoci tocki O smo spoznali ze v poglavju 2.1, slika 2.1. Drugo, prav tako desnosucno ortonormirano bazo Bg = ^Gi, G2, G3j, izberemo tako, dajo pripnemo v tezišcno tocko opazovanega prereza in premikamo ter vrtimo skupaj z opazovanim prerezom. Preslikava med dvema desnosucnima ortonormiranima bazama je prav rotacija: R :(Pi, P2,93) (£i,G2,G3) , kjer je Gi = Epi, za i = 1, 2, 3. Pri opisu dinamicnega deformiranja nosilca potrebujemo tako bazo za vsako tocško vzdolzš osi nosilca v vsakem cšasu, j (Gi (x, t), G2 (x, t), G3 (x, t)) , za t > 0 in x G [0, L]} , pri cemer velja Gi (x, t) = R (x, t) pi, za i G 1,2,3. (3.1) Posamicno bazo ^Gi (x, t), G2 (x, t), G3 (x, t) j pri izbranih x G [0, L] in t > 0 imenujemo materialna ali pomična baza, pripadajoce koordinatne osi pa oznacimo z malimi tiskanimi crkami x, y, z. Opomba 4 Zaradi boljče preglednosti zapis odvisnosti vektorjev materialne baze od x in t pogosto opustimo in pičemo kratko Gi, G2, G3. Vsak vektor lahko zapišemo v razlicnih bazah. V teoriji prostorskih nosilcev so nekatere vektorske spremenljivke izrazene glede na referencno, druge pa glede na pomicno bazo. V izogib nejasnostim z indeksom oznacujemo, v kateri bazi je zapisan vektor; na primer vektorju p = Agigi + afl2p2 + afl3p3 = AGiGi + QG2G2 + OG3G3 v bazi Bg pripada stolpec v bazi Bg pa stolpec a, = [ flgi ag2 ag3 ] T aG = [ AGi AG2 AG3 ] T Tudi bazne vektorje pomicne baze Gi lahko zapišemo v referencni bazi, le da pri komponentah indeks baze zaradi vecje preglednosti izpustimo Gi = Gii9i + G2i92 + G3i93 Gi,, = [ Gii G2i G3i ] . Do nejasnosti ne more priti, saj komponente baznih vektorjev pomicne baze glede na pomicno bazo ne potrebujejo crkovnih oznak g i,G = 1 0 0 G 2,G = G 3,G = 0 0 1 Matriki rg, ki pripada rotaciji R, v bazi Bg pripada komponentni zapis Rg — G11 G\2 G13 G21 G22 G23 G31 G32 G33 med komponentnima zapisoma vektorja ~a pa velja povezava ag = Rg aG. (3.2) (3.3) Rotacijska matrika rg v (3.2) je ortogonalna matrika z enotsko determinanto [Geradin, Rixen, 1995] R-1 = Rj in detRg = 1, (3.4) zato mnozica vseh rotacijskih matrik evklidskega vektorskega prostora v izbrani ortonormirani bazi tvori specialno ortogonalno grupo SO3 z notranjo operacijo komponiranja matrik. Lastnosti (3.4) obicajno zdruzšimo v zapis (3.5) RgRj = Rj Rg = 1, v katerem I oznacuje enotsko matriko 100 010 0 0 1 Med zapisoma vektorja a v obeh bazah tako velja še inverzna zveza T, aG = Rj ag. (3.6) Vidimo, da ima rotacija R dvojni pomen: 1. je preslikava, ki preslika poljuben vektor a v nov, zavrten vektor Ra. Primer: bazni vektor gg1, zapisan v referencni bazi 1 01g = 0 0 se z rotacijsko matriko Rg zasucše v bazni vektor G1 , rezultat preslikave je stolpec = G1,g ki pa je še vedno zapisan v isti (referencni) bazi; 2. je koordinatna transformacija med referencno in materialno bazo. Komponente poljubnega vektorja ~d, zapisanega v materialni bazi Bg, matrika Rg transformira v komponente tega vektorja v referencni bazi Bg. Na primer baznemu vektorju G1 v pomicni bazi pripada stolpec G11 G12 G13 1 G11 G21 G22 G23 0 = G21 G31 G32 G33 0 G31 g 1 ,G = 1 0 0 Ko ga transformiramo z matriko rg, dobimo po enacbi (3.3) zapis istega vektorja v referencni bazi: G11 G12 G13 1 G11 G21 G22 G23 0 = G21 = G1g G31 G32 G33 0 G31 Rotacijska matrika rg transformira tudi med matricnima zapisoma ag in Ag linearnega operatorja A: (3.7) ag = Rtagrg in AG = rg agrt . lGRg G = Rg AgR T Enacbo (3.7) smemo uporabiti tudi za transformacijo rotacijske matrike Rg v rotacijsko matriko Rg: Ob upoštevanju lastnosti (3.5) dobimo RG = Rg Rg rT . Rg = Rg = R. (3.8) S tem smo pokazali, da ima rotacijska matrika v obeh bazah, med katerima slika rotacija, enak kompo-nentni zapis. Zato indeks baze za rotacijsko matriko pogosto opušcamo. 3.2 Parametrizacija rotacij z rotacijskim vektorjem Vzemimo rotacijo za kot $ okoli osi, doloceni z enotskim vektorjem n. Rotacijski operator R in rotacijski vektor $ = $ n povezuje vektorska enacba ^ $ sin $ - $ 1 — cos $ -v . _ $. , $ Ra = a H--— $ x a H--^-$ x $ x a , za vsak vektor a. $ $2 (3.9) Enacba (3.9) je znana kot Rodriguesova formula, njeno izpeljavo pa najdemo na primer v delu [Geradin, Rixen, 1995] in v [Zupan, 2003]. Opomba 5 V enačbi (3.9) smo rotacijo R zapisali v odvisnosti od rotacijskega vektorja $. Odvisnost je nelinearna; kadar jo čelimo posebej poudariti, zapičemo R = R $ in rečemo, da je rotacije R parametrizirana z rotacijskim vektorjem $. Naj bosta $ in v neka rotacijska vektorja, ki preko enačbe (3.9) določata rotaciji R (v) in R . Rotacijo, ki pripada dvema zaporednima rotacijama izračunamo s kompozitumom operatorjev: R (v) R . Z uporabo (3.9) takoj vidimo, da ne velja aditivnost rotacijskih vektorjev v smislu komponiranja rotacij, saj je R (V + — = R (V) R ; rezultat je direktna posledica nelinearne povezave med rotacijo in parametrom. Zaporedno združevanje rotacijskih vektorjev v nov rotacijski vektor, ki pripada kompozitumu pripadajočih rotacij lahko formalno pičemo vendar se moramo zavedati, da je enak tistemu rotacijskemu vektorju, ki pripada rotacijski matriki kom-pozituma R ($) r($ . Opozorimo še na to, da postopek določitve rotacijskega vektorja, ki pripada rotacijski matriki ni enoličen (poleg vektorja $ © $ kompozitumu R ($) R ($ J ustreza tudi vsak koli- nearen vektor velikosti $ © $ + 2kn za vsako celo število k). Kljub temu pa obstaja numerično stabilen postopek za določitev rotacijskega vektorja iz rotacijske matrike velikostnega reda do 2n znan pod imenom Spurrierov algoritem [Spurrier, 1978]. Opis uporabe tega algoritma najdemo v [Simo, Vu-Quoc, 1986, pregl. 12]. Nelinearno odvisnost operatorja R od parametra $ pa ne smemo zamenjevati linearnim delovanjem operatorja R na argumentu $, torej z dejstvom, da je R linearna preslikava: R (a + V = Ra + R$ R (Xa) = XRa, za poljubna vektorja $,b G IR3 in skalar \ G IR. Vektorski produkt' x' vektorja $ s poljubnim vektorjem $ lahko alternativno zapišemo z antisimetricnim operatorjem 0 z enolicnim predpisom 0$ = $ x $, za vsak vektor a. Ker je R linearen operator argumenta $, lahko iz vektorske oblike Rodriguesove formule (3.9) izlocimo argument va in dobimo R = i + sin$ e + i-^ e2, (3.10) $ $2 kjer I oznacuje identicno preslikavo, Ta = ~a za vsak ~a. (Če rotacijski vektor razcepimo po vektorjih baze $ = $gl9l + $g2V2 + $g3V3, (3.11) potem operatorju © pripada v tej bazi antisimetricna matrika ?g2 1 (3.12) _ -$g2 $g1 0 rotacijski matriki pa zapis 0 -$g3 $g2 ©g = 3 g $ 0 —$gi L — $g2 $gi 0 R = i + ^ eg + ež. (3.13) Vektor $v imenujemo osni vektor operatorja e in pisšemo e = e $v 3.3 Odvodi in variacije rotacijskih količin Pripravili smo potrebne osnove, da lahko izpeljemo fizikalne kolicine ukrivljenost, kotna hitrost in pospešek, ki igrajo kljucno vlogo v teoriji prostorskih nosilcev, ter si ogledamo variacije rotacijskih kolicin. Pri tem se opiramo na osnovne definicije odvajanja in variacije vektorskih funkcij in operatorjev, zbrane v dodatku A. Pri izpeljavi vseh treh omenjenih kolicin izhajamo iz zveze (3.1) med referencno in pomicno bazo Gi (x,t) = R (x,t) $i za i = 1,2,3, (3.14) in njenim inverzom 9i = RT(x, t) Gi (x, t), (3.15) v katerem nastopa rotacija kot tista preslikava, ki slika med poljubnima dvema ortonormiranima bazama vektorskega prostora. Zvezo (3.14) odvajamo posebej po kraju in po casu. Odvoda operatorja R oznacimo z R' in E, pomenita pa krepka odvoda operatorja v skladu z definicijo 2 (glej dodatek A1): G = R'pi (3.16) Gi = R pi. (3.17) V teh izrazih 9 nadomestimo s pravilom (3.15) Gi' = R'Rt Gi (3.18) Gi = R Rt Gi. (3.19) Produkt operatorjev ERT = i (p) je antisimetricen, i (P) = -iT (p), kar pokazemo s preprosto izpeljavo EET = / jd (RET) = -d/, dt v ; dt ' E Rt + EE t = 0, i?Et = - (iRt)T . (3.20) Ker se dolzina vektorja Gi s casom ohranja, doloca antisimetricni operator i (p) hitrost vrtenja baznih vektorjev Gi s casom. Zato ga imenujemo operator kotne hitrosti. Ker je antisimetricen, ga lahko enolicno opišemo z osnim vektorjem p i (p) a = P x a, za vsak a. (3.21) Osni vektor P operatorja i (P) imenujemo kotna hitrost. Z enakim sklepanjem pokazemo, daje tudi operator i (K) = R'RT antisimetricen, Q (K) = -QT (K). Ker doloca hitrost spreminjanja smeri baznih vektorjev Gi pri spreminjanju lege tocke na osi nosilca, ga imenujemo operator (psevdo-) ukrivljenosti, njegov osni vektor p pa (psevdo-) ukrivljenost. Zanj velja Q (p) a = p x a, za vsak a. Ukrivljenost in kotna hitrost sta tisti kolicini, ki povezujeta vektorje pomicne baze z njenimi krajevnimi in cšasovnimi odvodi, Gi' = Q (K) Gi = P x Gi, (3.22) č?i = n (P) Gi = P x Gi (3.23) za vsak i = 1, 2, 3. Operator kotne hitrosti Q = iiRT ponovno odvajamo po casu in dobimo antisimetricen operator kotnega pospečka: Ji = E Rt + ii iiT. (3.24) Antisimetricnost operatorja (3.24) je preprosto pokazati, le izraz (3.20) ponovno odvajamo R rt + R RT = - (R Rt + R RT'T Osni vektor operatorja 11 imenujemo kotni pospeček in obicajno oznacujemo z L J I /v> _ /y>L 1 rv> _ /y>L J j ™ _ /v>L J I A[s] = ^ xo j ...{x xi-1) (x xl+1) ... xs j l ' f >1 JsA fjs]_ JsM fjs]_ JsM fjs]_ >P /y>L J _ /y>L J I I /y>L J _ /y>L J 1 | /y>L J _ /y>L J 1 | /y>L J _ rp ^n J ♦♦♦ \ ^*_i i \ ^* ^l+1 ' ' ^^ ^ za Gaussove točkex[s] z intervala [a, b], i = 0,1,..., s. V ucbenikih predstavljene Gaussove tocke in uteZi za integracijo na intervalu [-1,1] z linearno transformacijo prenesemo na ustrezno območje integracije, na primer L/2 1 M J /(x) = 4 J e + ^ ^ 4 jr /(^j), 0 -1 j=0 kjer je / integrabilna funkcija, Xj G (0, -f) so z linearno preslikavo x = -{ + - prirejene klasicne Gaussove integracijske tocke {j G (-1,1) in so Gaussove utezi. 5.3.6 Postopek dodajanja linearnih popravkov in popravljanje ostalih količin rp Rezultat Newtonove iteracije i + 1 so linearni popravki "ug and . V linearnih prostorih popravke v skladu s (5.51) kar prištejemo trenutni vrednosti priblizka iz iteracije i. Tako na primer postopamo pri popravljanju pomikov up+1,g = up,g + Tako popravljanje pa ni ustrezno pri rotacijskih kvaternionih. Kot smo omenili ze v poglavju 4.2.2, rotacijske kvaternione zaporednih rotacij zdruzujemo multiplikativno: v referencni bazi jih dodajamo z leve, enacba (4.36), v pomicni pa z desne, enacbi (4.39) in (4.40). Po dodajanju popravka pa mora rotacijski kvaternion ostati rotacijski, kar pomeni, da ohrani normo 1. To se zgodi le, kadar je tudi popravek rotacijski kvaternion. Toda linearni popravek ni rotacijski kvaternion, saj njegova velikost ni enaka ena in se z vecanjem števila iteracij praviloma priblizuje nic, torej postaja zelo majhna. Zato popravka ne moremo direktno dodajati. Iz poglavja 4.3 poznamo povezavo med linearnim popravkom rotacijskega kvaterniona in linearnim popravkom rotacijskega vektorja , enacba (4.57), iz definicije rotacijskega kvaterniona (4.22) pa iz znamo izracunati pripadajoč rotacijski kvaternion Afcg: "tfg = 0 + "tfg = 2"fcg O ^ (5.77) = |čtfg | Afcg =cos — + sin T. (5.78) A/cg imenujemo multiplikativnipopravek rotacijskega kvaterniona, saj ima vselej normo 1 in zato predstavlja pravi rotacijski kvaternion. Ker je izrazen glede na bazo Bg, ga multiplikativno dodamo z leve strani. Rotacijski kvaternion v odvisnosti od rotacijskega vektorja lahko zapišemo tudi z eksponentno preslikavo kvaternionskega argumenta Afcg =cos — + sin — = exp( "/j . (5.79) Dokaz najde bralec v dodatku B. Zanimivo je, da se eksponentna preslikava pojavi tudi pri popravljanju trenutne rotacije v vektorski parametrizaciji rotacij, le da gre pri rotacijskem vektorju za eksponentno preslikavo linearnega antisimetricnega operatorja z osnim vektorjem [Zupan, 2003]. Popravljanje p rotacijskih kvaternionov z linearnim popravkom "fcg ima torej obliko kp+1,g = exP ("kp O O O ^p,g. (5.80) Pravilnost izraza (5.80) za novi rotacijski kvaternion di+1;g potrdimo še z obratno analizo, v kateri li- •—p neariziramo desno stran v (5.80) in preverimo, ce je rezultat ravno 5kg. Po definiciji 3 (dodatek A1) je smerni odvod operatorja F v tocki a v smeri b enak DFa _d_ da F a + ab V nasšem a=0 vg vg vg ^ primeru išcemo smerni odvod di+1;g v tocki dig v smeri 5jg, le formalna vsota a + ab se v primeru multiplikativnih prostorov pacš izvede drugacše: D( kg+1>g; k i, 5dg da exp ( a5kPg o jp i,g i O Vig a=0 Ker skalar a nastopa direktno in ne samo posredno v argumentu eksponentne funkcije, je njen odvod po a kar obicajen odvod eksponentne funkcije kk ui+1,g k i,, 5dg exp ( a5jpg o jp i,g I O 5kg O ^g O ^Ig a=0 jg vg* v jg vg 5dg O di,g O V O di,g = 5dg Rezultat ponovno potrjuje pravilnost upoštevanja linearnega popravka rotacijskega kvaterniona po enacbi (5.80). Ob upoštevanju enacbe (5.79) lahko pravilo popravljanja (5.80) zapišemo krajše kk i+1,g avp O kg,g. (5.81) b Opomba 11 Podoben postopek za upočtevanje popravkov v Cliffordovi algebri sta izpeljala tudi McRo-bbie in Lasenby (1999). Abstrakten algebrajski zapis količin v Cliffordovi algebri je na prvi pogled sicer bistveno drugacčen od kvaternionske algebre, ki se ves cčas spogleduje z navadnim, sčtirirazsezčnim prostorom IR4, toda kvaternionska in čtirirazsečna Cliffordova algebra dejansko predstavljata isto algebrajsko strukturo [Ward, 1997]. Tehnično razliko med obema algebrama predstavlja definicija notranje operacije množenja - predznak skalarnega produkta, ki nastopa v kvaternionskem produktu, je negativen - le ta pa vpliva na ostale tehnične razlike. Vendar Lasenby in McRobbie (1999) le delo Simota in Vu-Quoca (1988), ki temelji na interpolaciji komponent rotacijskega vektorja, prevedeta v jezik algebre; pri tem pokačeta, da je linearizacija enačb prostorskega nosilca v Cliffordovi algebri precej preprostejča, kot pa v Liejevi grupi SO (3). Poleg osnovnih neznank in njihovih variacij moramo v vsakem koraku iteracije popravljati prve in druge odvode po x. Prek linearizacije interpolacijskih nastavkov smo ugotovili, da polacijo kot za osnovne neznanke uporabljamo tudi za njihove variacije, glej (5.57)-(5.58), pa zaradi aditivnosti tudi za odvode osnovnih spremenljivk N +1 ug (x) = £ L'g (x) ug (5.82) p=0 N +1 ug' (x) = £ L' (x) ug. (5.83) p=0 V primeru rotacijskih kvaternionov tak postopek ni ustrezen. Odvode kvaternionov popravljamo prek ukrivljenosti in njenega odvoda. Zato si najprej oglejmo, kako popravljamo ukrivljenost. Ukrivljenost K popravljamo v skladu z njeno 'aditivno' naravo v tocno doloceni bazi. Iz enacbe (4.48, levo) z uporabo tudi njihove enako inter-pri pomikih (5.81) izrazimo Ki+1,g: Ki+1,3 = 2 (Ag o o k* o Afcg* = (2Akg o +2Afcg o ky o k * ◦ Akg* = o Akg* + 2Akg o fc^ o o Akg* = AKg + Akg o Ki>g o Akg (5.84) AKg = 2Akg o Akg*. (5.85) Argument fci+1,g ukrivljenosti Ki+1g opozarja, da upoštevamo le tisti del celotne ukrivljenosti, ki je posledica dodane rotacije in ne celotne rotacije ki+1,g, kot na primer v enacbi (4.48). Argu- menta praviloma ne bomo vec pisali, razen takrat, ko bi utegnilo priti do zamenjave pojmov. V primeru uporabe zapisa v bazi Bg ne moremo ukrivljenosti direktno seštevati (5.84), temvec moramo prispevek ukrivljenosti iz predhodnega koraka ustrezno transformirati preden ga prištejemo popravku ukrivljenosti. Za izracun popravka ukrivljenosti AKg potrebujemo odvod popravka Akg, kije izpeljan v dodatku B, enacba (B.7) Akg (»g,«g) = a (f) £. Matrika A ima obliko neskoncne vrste (B.6), popravek rotacijskega vektorja 5$g in njegov odvod pa izracunamo v skladu z enacbo (4.57) in njenim odvodom; izraza sta predstavljena v preglednici 5.1. Preglednica 5.1: Popravki rotacijskega vektorja in njegovih odvodov Table 5.1: Update of rotational vector and its derivatives 5tfg = 25kg o ki,g R3 5^g = 25kg 25kg' o ki,g + 25kg o ki,g R3 - -// + 25kg o ki,g 5< = o ki,g + 45^ o ki,g R3 Popravek 5kg in njegov odvod 5kg smo ze predhodno izrazili z vozlišcnimi popravki, enacba (5.58). Ce enacbo (5.84) transformiramo v pomicno bazo , dobimo ki+1,Gi+i = ki*+1,g o AKg o ki+1,g + k*+1,g o Akg o ki,g o Akg* o ki+1,g = ki*+1,g o AKg o ki+1,g + k*g o Ki,g o ki,g = AKGi+i + Kj,Gi. (5.86) Ukrivljenost je torej aditivna kolicina, vendar le v primeru, ko za vsako nastopajoco ukrivljenost uporabimo ustrezno pomicno bazo. Izraz (5.86) uporabimo za izracun novega priblizka Ki+1,Gi+1, odvod rotacijskega kvaterniona ki+1,g pa izracunamo iz enacbe (4.48, levo), tako da le ukrivljenost Ki+1g zapišemo v bazi Bci+1: ki+1,g = 2Ki+1,g O ^»+1,g = 1gi+1,g O gi*+1,g O Ki+1,g O gi+1,g O gi*+1,g O fci+1,g = 1 gi+1,g O Ki+1,Gi+i O ?g01 * (5.87) Opomba 12 V enačbi (4.48) poleg celotnega kvaterniona rotacije g) + AoGi+i. Zadnja vrstica koncnega izraza (v + vi) se odšteje s clenoma (ii + iv). Odvod ukrivljenosti v pomicni bazi Kl'+1,Gi+i je enak vsoti prvih dveh vrstic koncnega izraza (v + vi) ter clenov (i) in (iii): ' Kl+1,Gi+i = = fl*+1,g o AKg o fl+1,g + -AKGi+i o Kl+1,Gi+i (fl+1,g) - -Kl+1,Gi+i (fl+1,g) o AKGi+i + fl*g o Kl',g o fl,g + 1 Kl,Gi o Kl,Gi (fl,g) - 1 Kl,Gi (fl,g) o Kl,Gi. Dvakrat uporabimo enacbo (4.52) in dobimo preprosto zvezo K i+1)Gi+1 = AK G i+i + K iU. (5.90) Tudi odvodi ukrivljenosti so v primernih pomicnih bazah 'aditivni'! Popravek ukrivljenosti v trenutni bazi izracunamo v skladu z enacbo (4.52) AKGi+i = j*+1,g O AKg O ji+1,g + 1 (0R (KGi+1 (ji+1,g)) - 0L (KGi+i (ji+1,g))) AKGi+i. +i - qi+1,g ◦ ^Kg ◦ 2 VK°i+i v^+lg;; - vl vKG;+i ^rn.g;;; ^KGi+i. Drugi odvod rotacijskega kvaterniona, di+1;g, izracunamo iz prvega odvoda ukrivljenosti, KGi+i. V enacšbi 1 i+i ~ "*i+1,g " '"t+1,g " "*i+1,g 1 2 v^R v'"Gi+i vit+1,gyy vl v'"Gi+i v^t+1,^ 'n+i,Gi+i Vt'+1,Gi+i = ji*+1,g O Ki+1,g O ji+1,g + ^ (0R (KGi+i (Vi+1,g)) - 0L (KGi+i (ji+1,g))) Ki+1,Gi Ki+1g nadomestimo z izrazom v/ v'' v* v' j * ' Ki+1,g = di+1,g O di+1,g + 2di+1,g O di+1,g '' ' ' = di+1,g O di+1,g — 2di+1,g O di+1,g O di+1,g O di+1,g _ v'' v * 1 v = ki+1,g O ki+1,g — 2 Ki+1,g O Ki+1,g v'' in izpostavimo di+1;g ter dobimo v'' ^ v v v' v* v di+1,g = 2 Ki+1,g ◦ Ki+1,g ◦ di+1,g + Vi+1,g O Ki+1,Gi+i O Vi+1,g O di+1,g - ji+1,g O 2 [0R (KGi+i (ji+1,g)) - 0L (KGi+i (ji+1,g))] Ki+1,Gi+i O ji*+1,g O ki+1,g. Dobljeni izraz še naprej preoblikujemo, dokler v njem ne nastopajo samo predhodno pripravljene koli-cine, na primer Ki+1,g O Ki+1,g O di+1,g = 2di+1,g O di*+1,g O 2di+1,g O di*+1,g O di+1,g = di+1,g O Ki+1,Gi+i O Ki+1,Gi+i in (0R (KGi+i (qi+1,g)) - 0L (KGi+i (qi+1,g))) Ki+1,Gi+i = Ki+1,Gi+i O Ki+1,Gi+i (Vi+1,g) - Ki+1,Gi+i (ji+1,g) O Ki+1,Gi+ Ki+1,Gi+i O Ki+1,Gi+i - Ki+1,Gi+i O Ki+1,Gi+i [0] v[0| v[0] Ki+1,Gi+i O KG[o] - KG[o] o ki+1,Gi+i = [0L (Ki+1,Gi+0 - 0R (Ki+1,Gi+0] V(3[o] ter dobimo 1 [0] di+1,g = 2 di+1'g O 7i+1,Gi+i O Vi+1,Gi+i + ji+1,g O Vi+1,Gi+i O jg - 2 ji+1,g O (Vi+1,Gi+i) - 0R (Vi+1,Gi+i)] VJ0!o] O ^ *. Zacetno ukrivljenost K[0] v prvi pomicni bazi BgM izracunamo v skladu z enacbo (4.48, desno) V[°] =2V[0]* oV[0]' VG[o] = 2qg O qg . 5.4 Numericni testi Predstavljeno metodo preverimo s primeri iz literature. Izbrani testi jasno kazejo na pravilnost formulacije in zmogljivosti implementacije, za katero velja, da • poleg zacetne ravne lege lahko podajamo tudi zacetno upogibno in torzijsko ukrivljeno lego nosilca; • so dopušceni poljubno veliki pomiki in zasuki vecji od 2n; • lahko natancnost analize povecujemo z vecanjem števila elementov in/ali z uporabo elementov višjega reda. Numericni rezultati so v celoti pridobljeni v programskem okolju Matlab [The MathWorks, 1999]. Pri izracunih se omejujemo na linearno elasticni materialni model, kjer imata linearna operatorja Cn in Cm iz enacb (6.51)-(6.52) v matricnem zapisu diagonalno matriko: CN = CM = EA1 0 0 0 GA2 0 0 0 GA3 GJ1 0 0 0 EJ2 0 0 0 EJ3 Y G = C N Y G, Kg = c M Kg. (5.91) (5.92) E in G sta elasticni in strizni modul materiala; A1 je plošcina prereza nosilca; J1 je tezišcni torzijski vztrajnostni moment prereza; A2 in A3 sta strizna prereza v smereh glavnih vztrajnostnih osi prereza G2 in G3 s pripadajocima tezišcnima vztrajnostnima momentoma prereza J2 in J3. Za namene primerjave z drugimi avtorji rezultate rotacijskih kvaternionov kg = k0 + k preracunamo v rotacijske vektorje po formuli k k (5.93) ti g = 2 arccos (ko) jn. k Osi koordinatnega sistema, ki ga dolocajo bazni vektorji referencne baze B z izhodišcem v tocki O, oznacujemo z velikimi tiskanimi crkami X, Y in Z. Prikaz neobčutljivosti elementa na strižno blokiranje. Klasicni elementi kazejo podcenjene pomike, kadar je strizni modul G zelo velik in kadar je višina prereza h zelo majhna v primerjavi z dolzino koncnega elementa (L/h je veliko število). Formalno lahko za klasicne elemente rešitev izrazimo v odvisnosti od parametra GL2/Eh2 (glej [Bathe, 1996]), kar utemeljuje te tezave. Kadar je pri reševanju enacb nosilca po metodi koncnih elementov kljub zgošcanju mreze pomik še naprej mocno podcenjen, govorimo o strižnem blokiranju elementa. Za predstavljeni koncni element se rešitev ne izraza direktno v odvisnosti od parametra GL2/Eh2, lahko pa numericno preverimo obcutljivost rešitve v odvisnosti od velikih vrednosti G in L/h. Neobcutljivost predstavljene metode koncnih elementov na strizno blokiranje ilustriramo s standardnim testom. Analiziramo ravno konzolo s precno tockovno silo F = 1 v prostem krajišcu, kot prikazuje slika 5.2. Ostali podatki o konzoli so: E = 107 G = 1013 L = 1 t = 0.1. Z Slika 5.2: Konzola s precno silo v prostem krajišcu Figure 5.2: Cantilever under a free-end vertical force Nerealno velik strizni modul G numericno doloca nosilec, ki ni sposoben striznega deformiranja. Ker zmanjševanje višine precnega prereza pri standardnih koncnih elementih hitreje privede do striznega blokiranja, višino h spreminjamo od (majhnih) 0.1 do (velikih vrednosti) 10. Konzolo modeliramo z enim, s petimi in z desetimi kvadraticnimi elementi (ena notranja tocka). Rezultate pomikov prostega krajišca v precni smeri zaradi nazornosti normiramo z referencno vrednostjo UZ,ref (numericna rešitev za mrezo 100 elementov). Spreminjanje napake pomika uz/«Z,ref v odvisnosti od strukturnega parametra GL2/Eh2 za razlicno goste mreze, ne = 1, 5,10, prikazujemo graficno v logaritemskem merilu na sliki 5.3. UZ /Uz,ref 1.15 1.1 1.05 1 0.95 0.9 4 ----n =1 ---n =5 n =10 7 „ 0 8 logJGL2/ Eh2) Slika 5.3: Precni pomik prostega krajišca v odvisnosti od parametra GL2/Eh2 Figure 5.3: Vertical displacement at the free-end vs parameter GL2/Eh2 5 6 Dobljeni rezultati so popolnoma neobcutljivi na spreminjanje strukturnega parametra za vse tri razlicno goste mreze koncnih elementov ter se z vecanjem števila elementov priblizujejo referencni rešitvi. Na sliki 5.4 podajamo primerjavo med našo rešitvijo, kjer nosilec modeliramo z enim elementom in rešitvijo s klasicnim Timošenkovim nosilcem, kije osnovan na pomikih [Bathe, 1996]. Test torej kaze, da uporabljena numericšna metoda nima strizšnega blokiranja. 2 0 -2 -4 tg lo -6 -8 -10 -naša rešitev — — Timošenko 5 6 o „ 7 2 2 logJGL /Eh ) 4 8 Slika 5.4: Primerjava precnega pomika prostega krajišca v odvisnosti od parametra GL2/Eh2 s klasicnim Timošenkovim nosilcem s striZnim blokiranjem Figure 5.4: Vertical displacement at the free-end vs parameter GL2/Eh2; comparison with the classical Timoshenko beam with shear locking problem Konzola z momentom. Vzemimo ravno konzolo dolZine L = 100 v smeri osi X in jo v prostem krajišcu obremenimo s tockovnim momentom My, slika 5.5. Nosilec se zaradi obtezbe deformira v ravnini (X, Z). Ostali materialni in geometrijski podatki so E = 2.1 ■ 104 G = 1.05 ■ 104 Ai = 20 A2 = A3 = 16 Ji = 6.4566 J2 = 1.6667 J3 = 666.66. O prerez: M2 x aG3(X) G2(X) Slika 5.5: Konzola z momentom Figure 5.5: Cantilever under free-end moment V preglednici 5.2 primerjamo pomika in zasuk prostega krajišca z analiticnimi rešitvami [Saje, Srpcic, 1986] za dva obtezna primera: My = 1 (majhna obtezba) in My = 100 (velika obtezba). V prvem obteznem primeru so pomiki in zasuki majhni, zato dosezemo dovolj dobro ujemanje z analiticno rešitvijo ze z izbiro enega elementa z dvema notranjima tockama. V drugem obteznem primeru so pomiki in zasuki veliki, pa vseeno ze dva elementa iste stopnje zadošcata za dobro ujemanje numericne resšitve z analiticšno. Preglednica 5.2: Pomiki in zasuki prostega krajišca konzole z momentno obtezbo Table 5.2: Free-end displacements and rotations of cantilever beam under free-end moment ne N u1 u3 §2 M2 = 1 1 1 0.00010 0.14286 0.00286 2 0.00014 0.14286 0.00286 2 1 0.00013 0.14286 0.00286 2 0.00014 0.14286 0.00286 tocna nelin. [Saje, Srpcic, 1986] 0.00014 0.14286 0.00286 M2 = 100 1 1 1.01867 14.23718 0.28571 2 1.35507 14.18880 0.28571 2 1 1.27107 14.20086 0.28571 2 1.35501 14.18880 0.28571 tocna nelin. [Saje, Srpcic, 1986] 1.35500 14.18880 0.28571 ne=število elementov, N=število notranjih kolokacijskih tock. Slika 5.6: Zviti nosilec Figure 5.6: Twisted beam Zviti nosilec. Za prikaz sposobnosti metode, da upošteva zasukane prereze v zacetnem stanju, smo izbrali standardni test za metodo koncnih elementov MacNeala in Harderja (1985). Nosilec je enostransko togo vpet, prosto krajišce pa je obtezeno z enotsko silo v precnih smereh Y in Z, kot kaze slika 5.6. Tezišcna os nosilca je ravna, pravokotni precni prerez pa se od zacetka do konca nosilca zasuka za kot 2n tako, da se os rotacije ujema s tezišcno osjo nosilca in da se velikost kota v odvisnosti od naravnega parametra x linearno spreminja od velikosti § = 0 pri x = 0 do § = 1 n pri x = L (slika 5.6). Tezišcna os nosilca v zacetni legi ostaja ves cas pravokotna na ravnine precnih prerezov. Zupan in Saje (2004) opozarjata, da taka zacetna geometrija nosilca pomeni, da so robovi nosilca v zacetni legi ukrivljeni. Zaradi pravokotnega prereza locimo dva obtezna primera za dve precni smeri nosilca, zaradi zvitosti pa se nosilec v obeh primerih deformira prostorsko in ne ravninsko. Ostali geometrijski in materialni podatki so: h = 1.1 t = 0.32 L = 12 E = 29 ■ 106 G = 11.885 ■ 106. Pomike prostega krajišca v smeri obteZbe po linearni teoriji (le prvi korak Newtonove iteracije) primerjamo s teoreticnimi rezultati MacNeala in Harderja (1985) in z analiticnimi rešitvami, dobljenimi iz linearizirane oblike enacb za prostorski nosilec po Reissner-Simovi teoriji nosilcev, glej preglednico 5.3. Stevilo uporabljenih elementov je bistveno nizje od števila priporocenih elementov (12) s strani avtorjev standardnega testa [MacNeal, Harder, 1985]. Numericni rezultati pomikov v prostem krajišcu nosilca kazejo veliko natancnost rezultatov ze pri majhnem številu elementov nizkega reda. Rezultati se ujemajo z analiticnimi na dve nenicelni decimalni mesti ze z uporabo enega samega elementa s petimi notranjimi kolokacijskimi tocškami ali z dvema elementoma s po dvema notranjima kolokacijskima tocškama. Za ujemanje na štiri znacilna decimalna mesta zadošca en sam element s sedmimi notranjimi tockami, dva elementa s štirimi internimi tockami in štirje elementi s po tremi internimi tockami. Primerljivo tocne rešitve najdemo še pri Zupanu in Sajetu (2004, 2006). Odstopanja od analiticne rešitve pa najdemo pri na primer Ibrahimbegovicu in Freyu (1993), ceprav uporabita 12 in 24 elementov. Teoreticni rezultati MacNeala in Harderja (1985) zaradi drugacšne teorije nosilcev malenkostno odstopajo od predstavljene analiticšne in numericšne resšitve. Preglednica 5.3: Pomiki prostega krajišca po linearni teoriji za zviti nosilec Table 5.3: Free-end displacements of twisted beam; linearized theory Fi F2 ne N U2 U3 U2 U3 1 3 0.005005 0.001584 0.001431 0.001677 5 0.005427 0.001725 0.001725 0.001754 6 0.005430 0.001720 0.001719 0.001749 7 0.005429 0.001719 0.001719 0.001749 2 2 0.005428 0.001728 0.001655 0.001711 3 0.005426 0.001716 0.001714 0.001747 4 0.005429 0.001719 0.001719 0.001750 4 1 0.005480 0.001795 0.001795 0.001755 2 0.005429 0.001719 0.001715 0.001747 3 0.005429 0.001719 0.001719 0.001750 analiticšno 0.005429 0.001719 0.001719 0.001750 teoreticšno [MacNeal, Harder, 1985] 0.005424 0.001754 ne=število elementov, N=število notranjih kolokacijskih tock. Bočna zvrnitev pravokotnega okvirja. S tem klasicnim primerom, ki gaje vpeljal ze Argyris in sodelavci (1978) in so ga študirali tudi mnogi drugi, analiziramo prostolezece podprt okvir oblike enako-krakega pravokotnega trikotnika s prerezom v obliki ozkega pravokotnika. Geometrijske karakteristike okvirja so prikazane na sliki 5.7. Preostale karakteristike okvirja, povzete po [Zupan, Saje, 2003], so: J1 = 2.16 Ai = 18 E = 71240 J2 = 0.54 A2 = 21.6 G = 27191 J3 = 1350 A3 = 21.6 L = 240. Slika 5.7: Bocna zvrnitev pravokotnega okvirja Figure 5.7: Lateral buckling of a right-angle frame Kljub temu, daje zacetna geometrija skupaj z momentno obtezbo v podprtih vozlišcih popolnoma ravninska, pa se zaradi izrazito visokega pravokotnega prereza ob narašcanju momentne obtezbe pri neki kriticni vrednosti momenta, Mcr, zgodi izklon iz ravnine, imenovan bocna zvrnitev. Primer dokazuje sposobnost metode, da zazna take vrste kriticnega stanja konstrukcije. Pojav bocne zvrnitve zaznamo s singularnostjo tangentne matrike tako, da se pricakovani kriticni vrednosti momentne obtezbe postopno priblizujemo. Upogibni moment Mcr smo dolocili za razlicno število elementov in razlicne stopnje elementov ter s tem pokazali ustrezno delovanje formulacije in racunalniškega programa za razlicno goste mrezše koncšnih elementov. Preglednica 5.4: Kriticni upogibni moment Mcr prostolezece podprtega pravokotnega okvirja Table 5.4: Critical moment Mcr of the simply supported right-angle frame ne = 1 ne = 2 ne = 4 ne = 6 ne = 8 ne = 10 N = 1 ±622.00 ±622.00 ±622.00 ±622.20 ±622.20 ±622.20 N = 3 ±622.00 ±622.43 ±622.24 ±622.23 ±622.22 N = 5 ±622.19 ±622.22 ±622.22 ±622.22 N = 7 ±622.22 ±622.22 [Jelenic, Saje, 1995] ±622.2 [Nour-Omid, Rankin, 1991] ±626.7 [Simo, Vu-Quoc, 1986] ±626 analiticna rešitev [Timoshenko, Gere, 1961] ±622.21 ne=število elementov, N=število notranjih kolokacijskih tock. V preglednici 5.4 primerjamo numericne rezultate z analiticnimi iz [Timoshenko, Gere, 1961] ter z nu-mericnimi rezultati drugih avtorjev. Na tri znacilna mesta natancno numericno rešitev dobimo ze z enim elementom nizke stopnje, vsekakor pa opazimo hitro konvergenco rešitve ob vecanju stopnje elementa. Popolnega ujemanja z analiticno rešitvijo [Timoshenko, Gere, 1961] pa ne moremo pricakovati, ker ta predvideva neskoncšne vrednosti za obe strizšni in za osno togosti ter za upogibno togost v ravnini kolena. Slika 5.8: Ukrivljena konzola Figure 5.8: Cantilever 45° bend Upogib ukrivljene konzole. Tudi ta primer, ki sta ga prva predstavila Bathe in Bolourchi (1979), sodi med klasicne testne primere predvsem zato, ker vkljucuje vse mogoce oblike deformiranja nosilca: upogib, strig, nateg in torzijo. Tezišcna os nosilca ima obliko osmine kroznega loka s polmerom 100. Nosilec je enostransko togo vpet, prosto krajišce paje obtezeno s silo, ki kaze ven iz ravnine. Loamo dva obtezna primera: F je 300 in 600. Prerez je enotski kvadrat. Geometrija nosilca je predstavljena na sliki 5.8, preostale geometrijske in materialne karakteristike pa so: h = 1 t = 1 R = 100 E = 107 G = E/2. Preglednica 5.5: Komponente r1,r2,r3 krajevnega vektorja prostega krajišca pri upogibu ukrivljene konzole Table 5.5: Components r1, r2, r3 of the free-end position vector of the cantilever 45° bend under out-of-plane force F = 300 F = 600 formulacija r1 r2 r3 r1 r2 r3 raven, N = 1 22.286 58.786 40.238 15.765 47.168 53.566 raven, N = 3 22.275 58.782 40.156 15.738 47.150 53.433 ukrivljen, N = 1 22.242 58.785 40.277 15.688 47.173 53.608 ukrivljen, N = 3 22.245 58.779 40.192 15.685 47.150 53.475 [Bathe, Bolourchi, 1979] 22.5 59.2 39.5 15.9 47.2 53.4 [Simo, Vu-Quoc, 1986] 22.33 58.84 40.08 15.79 47.23 53.37 [Cardona, Geradin, 1988] 22.14 58.64 40.35 15.55 47.04 53.50 [Crivelli, Felippa, 1993] 22.31 58.85 40.08 15.75 47.25 53.37 [Zupan, Saje, 2003] 22.28 58.78 40.16 15.74 47.15 53.43 število elementov=8, N=število notranjih kolokacijskih tock. Analiticne rešitve ni na voljo, zato rezultate primerjamo z numericnimi rešitvami drugih avtorjev. V preglednici 5.5 prikazujemo primerjave krajevnih vektorjev prostega krajišca. Vecina avtorjev je problem modelirala z osmimi ravnimi elementi, zato so primerjave rezultatov opravljene za mrezšo osmih elementov. Za primerjavo problem modeliramo še z osmimi ukrivljenimi elementi. Obtezbo v prvem obteznem primeru nanasšamo v sšestih, v drugem pa v dvanajstih enako velikih obtezšnih korakih. Zahtevano na-tancnost za zakljucek Newtonove iteracije je ^New = 10-7. Rezultati razlicnih avtorjev predstavljeni v preglednici 5.5 so primerljivi z našimi rezultati. Rezultati izracunani z mrezo osmih ravnih elementov s tremi internimi tockami se popolnoma ujemajo z rezultati iz [Zupan, Saje, 2003]. Med ravno in ukrivljeno modeliranimi nosilci ne opazimo vecjih razlik: pri vecji obtezbi so razlike rezultatov komaj opazne, za manjšo obtezbo pa so razlike rahlo vecje. Preglednica 5.6: Primerjava pomikov prostega krajišca v odvisnosti od števila korakov nalaganja obtezbe za primer ukrivljene konzole Table 5.6: Cantilever 45° bend: free-end displacements as a function of the number of load steps sštevilo korakov korakov Ui U2 U3 4 13.55148 -23.56059 53.43333 11 13.55149 -23.56057 53.43333 50 13.55149 -23.56057 53.43333 Prav ta primer Jelenic in Crisfield (1999a) uporabita za numericno testiranje odvisnosti od poti, lahko pa sluzi še za primerjavo konvergencnega polmera metod razlicnih avtorjev. V [Jelenic, Crisfield, 1999a] najdemo primerjavo odvisnosti od poti in potrebno število korakov nalaganja obtezbe za drugi (vecji) obtezni primer; primerjava vkljucuje številne koncne elemente za staticno analizo prostorskih nosilcev, vecini pomiki in rotacijski vektorji pomenijo osnovne spremenljivke problema, rotacije pa parametrizi-rajo z rotacijskim vektorjem. Primerjane metode za izracun rezultata potrebujejo tri korake nalaganja obtezbe. Za primerjavo izberemo mrezo osmih ravnih elementov s tremi internimi tockami. Konvergenco dosezšemo pri obtezševanju najmanj s sštirimi enako velikimi obtezšnimi koraki. V tabeli 5.6 prikazujemo rezultate pomikov ob obtezevanju konstrukcije s štirimi, enajstimi in petdesetimi enako velikimi obtezšnimi koraki. Primerjava nakazuje, da koncšni pomiki niso bistveno odvisni od sštevila in velikosti obteznih korakov. Glede na to, da lahko Jelenic in Saje (1995), z metodo kije osnovana na interpolaciji zgolj rotacijskih vektorjev ter Zupan in Saje (2003) z metodo, kije osnovana na interpolaciji deforma-cijskih kolicin, obtezbo nalozita v enem samem koraku lahko sklepamo, da izbira osnovnih neznank bistveno vpliva na konvergencni polmer metode. Slika 5.9: Konzola, upognjena v spiralo Figure 5.9: Cantilever bent to the helical form Konzola, upognjena v spiralo. Zadnji primer dokazuje sposobnost pricujoce formulacije in še posebej kvaternionske parametrizacije za primere, ko velikost zasukov pri deformaciji nosilca bistveno presega Slika 5.10: Koncna deformirana oblika konzole, upognjene v spiralo Figure 5.10: Cantilever: deformed shape at the final load step Slika 5.11: Konzola, upognjena v spiralo: odvisnost pomik-sila Figure 5.11: Load-displacement curve for the helical beam 2n, poleg tega pa se rotacije spreminjajo izrazito ciklicno. Primer ravne konzole, ki jo z momentom in s silo iz ravnine navijemo v spiralo, je prvi predstavil Ibrahimbegovic (1997). Geometrija konzole je prikazana na sliki 5.9, preostali geometrijski in materialni podatki pa so: GA2 = GA3 = EAi = 104 L = 10 EJ = EJ3 = GJi = 102. Obtezni primer vsebuje socasno nalaganje sile F = 50A in momenta M = 200nA, A G [0,1], v 1000 enakomernih korakih. Rezultat takega obtezevanja je konzola, upognjena v spiralo, ki jo sestavlja deset obrocev; oblika tezišcne osi pri polni obtezbi A = 1 je prikazana na sliki 5.10. Pomembno je poudariti, daje nosilec pri polni obtezbi upognjen v smeri, kije nasprotna smeri obtezne sile F. Rezultati so dobljeni z mrezšo 25 elementov s sedmimi notranjimi kolokacijskimi tocškami. Pomiki prostega krajišca v smeri obtezne sile u2 v odvisnosti od obteznega faktorja A so prikazani na sliki 5.11 in se ujemajo z rezultati iz [Ibrahimbegovic, 1997]. S tem smo pokazali, da kombinacija kvaternion-ske parametrizacije rotacij in predlagane numericne metode za reševanje enacb prostorskega nosilca po Reissner-Simovi teoriji omogoca analizo konstrukcij, pri katerih se razvijejo veliki zasuki z izrazitimi oscilacijami. 6 Enačbe prostorskega nosilca za dinamiko V tem poglavju enacbe prostorskega nosilca za staticno analizo, predstavljene v poglavju 5.1, razširimo na dinamicno analizo. Izpeljavo izvedemo v vektorski parametrizaciji rotacij, po primerjavi dobljenih enacb z drugimi avtorji pa enacbe preoblikujemo v kvaternionsko obliko in takšne tudi rešujemo. Oblika vecine enacb prostorskega nosilca za dinamiko je formalno enaka enacbam prostorskega nosilca za statiko, dodana pa je odvisnost od casa. Dinamicno ravnotezje telesa dolocata izreka o gibalni in vrtilni kolicšini, ki zahtevata, da je vsota vseh zunanjih sil enaka odvodu gibalne kolicšine po cšasu in da je vsota vseh momentov zunanjih sil enaka odvodu vrtilne kolicšine po cšasu: K je gibalna kolicina in Lo je vrtilna kolicina. Vrtilno kolicino in moment smo zapisali glede na iz-hodišce referencnega koordinatnega sistema O. Ker sta levi strani enacb staticnega (5.1)-(5.2) in di-namicnega (6.1)-(6.2) ravnotezja enaki, se v tem poglavju osredotocimo na izpeljavo gibalne in vrtilne kolicine in njunega odvoda po casu. 6.1 Gibalna količina in njen odvod po času Gibalna količina je integral produkta mase in hitrosti delca po obmocju nosilca (6.1) (6.2) (6.3) (6.4) (6.5) V Uporabimo razcep vektorja pD v referencni bazi, glej poglavje 3, p® £3) = 652 + (6.6) in integral po volumnu nosilca nadomestimo z integraloma po dolzini nosilca [0, L] in po površini prereza Ar. Upoštevamo še, da predstavlja os nosilca krivuljo, ki povezuje tezišca precnih prerezov in da sta zato statična momenta prereza glede na tezišce prereza enaka nic: J -2 dA = 0, Ar J -3 dA = 0. Ar Drugi clen integrala (6.5) je tako nic, saj se v njem pojavijo staticni momenti: J pQR (x,t) p? (-2,-3) dV V = J pQR (x, t) I y -2dAp2 + j -3dAp3 I dx = 0. 0 Ur Ar / Prvi clen izraza (6.5) je od -2 in -3 neodvisen, zato lahko pripadajoci integral po Ar izvrednotimo in dobimo preprosti integralski enacbi za gibalno kolicino in njen odvod: L K (t) = j pArp(x, t) dx, (6.7) 0 L K (t) = j pArp(x, t) dx. (6.8) 0 Upoštevali smo, da se pri izbrani parametrizaciji osi nosilca (naravni parameter osi nosilca v zacetni legi) meje integracije [0, L] ne spreminjajo, zato je odvod integrala kar integral odvoda integranda. Kadar potrebujemo gibalno kolicšino le za del nosilca, na primer od levega krajisšcša x = 0 do tocške x, velja x K (x,t) = J pArp(-,t) d-, 0 x K (x,t) = J pArp(-,t) d-. 0 6.2 Vrtilna kolicina in njen odvod po casu Vrtilna kolicina glede na izhodišce referencnega koordinatnega sistema je integral vektorskega produkta krajevnega vektorja in gibalne kolicine delca po obmocju nosilca (t) = J rD (x, -2,-3, t) X ppD (x, -2, -3, t) dV. V Uporabimo enacbi za krajevni vektor (2.2) in hitrost delca (6.4) in dobimo: = J (V + Rpr®) x p(r + QRpr®) dV. (6.9) V Zaradi vecje preglednosti zapisov smo opustili eksplicitni zapis odvisnosti od argumentov, zato še enkrat opomnimo, da so vse kolicine, razen p®, odvisne od casa t. Uporabimo razcep vektorja p® v referencni bazi (6.6), integral po volumnu nosilca nadomestimo z zaporednim integriranjem po dolzini nosilca [0, L] in po površini prerezov A ter iz integrala po površini prerezov izlocimo kolicine, ki so neodvisne od lege delca D v ravnini precnega prereza: = / (p + R (6p2 + 6p3)) x p p + QR (£252 + £353) dV V = p Arr x rdx + pr x QRp2 I / £2dA + pr x QRp3 / £3dA I dx 0 VA A L + J pR I p2 x r J £2dA + p3 x r J £3dA | dx 0 A A 33 L + E I] I I pRLi x ^Rpj y £i£jdA | dx. A i=2 j=2 ; L L Ker sta staticna momenta prereza glede na tezišce prereza enaka nic, so cleni tretje in cetrte vrstice zgornje enacbe enaki nic. Za preglednejši zapis zadnje, pete vrstice vpeljemo 'mehanske vztrajnostne momente prereza' z zacasno, nestandardno oznako = / p£i£j dA, i, j = 2,3 (6.10) Ar in dobimo, daje vrtilna kolicina glede na tezišce prereza enaka vsoti petih integralov vzdolz tezišcne osi: LO = J p Arr x rdx + £ £ J Ij (Rp x QRpj) dx. (6.11) Operator Q lahko zapišemo kot vektorski produkt (3.21) in nato uporabimo znano pravilo o dvojnem vektorskem produktu [Bronštejn, Semendjajev, 1988, str. 610], a x (pi x p =b (a • c) — c (a • b j , (6.12) ki velja za poljubne vektorje a, b in c. Upoštevamo še, daje tudi (Rgg1, Rf2, Rf3) ortonormirana baza in izpeljemo naslednji izraz za vrtilno kolicšino v vektorski obliki L L f o = J p Arf x rdx + J (/22 + I33) fdx (6.13) 0 0 3 3 L - E E / /ijRggi (Rgj' • f) dx i=2 j=2 0 V izreku o vrtilni kolicini nastopa odvod vrtilne kolicine po casu, ki ga izpeljemo z direktnim odvajanjem izraza za vrtilno kolicino (6.9): (f + QRpf) x p (r + IR^0) (6.14) V + (f + Rf®) x p (f + (QR)' f®) dV. Upoštevali smo, da se pri izbrani parametrizaciji nosilca (naravni parameter osi nosilca v zacetni legi, referencni prerez) meje integracije [0, L] x Ar ne spreminjajo, zato je odvod integrala kar integral odvoda integranda. Prva vrstica izraza (6.14) vsebuje vektorski produkt enakih vektorjev in je zato enaka nic, v drugi vrstici pa upoštevamo (QR)" = 11R + QQR in izraz za Lo razbijemo na vec clenov: LO = j pr x rdV (i) (6.15) V + J pr x 1 Rfr®dV (ii) v + J pr x QQRf®dV (iii) v + J pRf® x rdV (iv) v + J pRf® x 11 RfDdV (v) Vv + J pRf® x HRfDdV. (vi) V (Člene obravnavamo loceno. Ponovno uporabimo razcep f® po baznih vektorjih reference baze (6.6) in integral po volumnu nosilca nadomestimo z integraloma po obmocju [0, L] x Ar. Upoštevamo, da sta staticna momenta glede na tezišce prereza enaka nic. Edina kolicina, kije odvisna od £2 in £3, je f, zato lahko vse ostale kolicine izpostavimo iz integralskega znaka integracije po prerezu Ar. Prvi clen izraza (6.15) je neodvisen od £2 in £3, zato lahko integral po Ar izvrednotimo in dobimo (,)=/ PA (f x f) dx. L V drugem (ii), tretjem (iii) in cetrtem (iv) cienu izraza (6.15) nastopa pT le linearno, zato se v cienih pojavijo staticni momenti, ki so enaki nic: (ii) = J pr x 1R l J C2dAg2 + J &dAg3 | dx = 0, 0 Ur Ar (iii) = J pr x QQR l J C2dAg2 + J&dAg3 | dx = 0, 0 Ur Ar L (iv) = — J r x pRpjPdV 0 = —J r x pR U&dAg2 + J &dAg3 | dx = 0. 0 Ar Ar Pri spremembi zapisa petega clena izraza (6.15) uporabimo antisimetricnost operatorja , ki ga lahko zato zapišemo s pripadajocim osnim vektorjem uu: 1 a = uu x a in pravilo o dvojnem vektorskem produktu (6.12). Iz uporabe razcepa (6.6) sledi Rpr = 6R12+(3R93, Rpr • Rpr = C2 + C2. Ponovno uporabimo pomozno oznako Ij (6.10) in izraz (v) preuredimo v vektorsko obliko: L (v) = J (I22 + I33) ^dx 0 L — J Rg2 {122U • Rg2 + I23W • dx 0 L — Rg3 (123U • Rg2 +133^^ • Rg^ dx. 0 V šestem clenu izraza (6.15) najprej uporabimo prirejeno pravilo o dvojnem vektorskem produktu a x (b x ^c x dJJ = b ^a • ^c x d^ — ^c x d j ^a • b j , ki velja za poljubne vektorje a, b, c in d, nato pa ob upoštevanju razcepa (6.6) iz w x Rpr = C2w x Rxj2 + C3w x Rgg3 izpeljemo vektorsko obliko izraza (iv): L (iv) = - J (U • Rg2) [I22 (L x Rg2) + I23 (U x Rg3)] dx 0 L - J (U • Rg3) [I23 (L x Rg2) + I33 (u x Rg3)] dx. Sedaj lahko zapišemo odvod vrtilne kolicine v splošni vektorski obliki L L Lo = pAr (r x r) dx + (I22 + I33)Ldx 0 L -J RL2 0 L I22 U • Rg2 + I23 L • Rg3 dx dx 7 Rg3 V [I23 (L • R92) + I33 • Rg3 0 L - J (l • Rgg2) [I22 (L x Rgg2) + I23 (L x Rg3)] dx 0 L - J (l • Rgg3) [I23 (L x Rgg2) + I33 (L x Rg3)] dx. (6.16) Izraz (6.16) za Lo zelimo zapisati tudi v matricni obliki, saj pricakujemo krajši in preglednejši zapis. V ta namen vektorske kolicine zapišemo v referencni bazi (Li, g2,53): r = rigi + r2g2 + r3g3, Rg2 = G2 = G12gi + G22g2 + G32g3, Rg3 = G3 = G13gi + G23g2 + G33g3, L = Uigi + U2g2 + U3g3, č5 = U igi + U 2g2 + U 3g3 oziroma v zapisu s stolpci ri rg = r2 . r3 _ " Ui " U i ^g = U2 = U 2 . U3 . UU 3 Gi2 Gi3 G2g = G22 G3g = G23 G32 G33 (6.17) Opomba 13 Za poljuben vektor a z zapisom ag = [ a1 a2 a3 ]T v referenčni bazi velja, daje njegov odvod a' po x v referenčni bazi enolično določen z odvodi komponent: a'g = [ al a'2 a'3 ] .Podobno • T je odvod po času a določen s stolpcem ag = [ a 1 a 2 a 3 ] . Zato je izbira zapisa parcialnega odvoda u v referenčni bazi s parcialnimi odvodi komponent povsem ustrezna. Iz istih razlogov je vektor r' [ ]T •• [ ]T določen s stolpcem r'g = [ r' r2 r3 ] in vektor u s stolpcem ug = [ ii1 ii2 ii3 ] . Tak zapis pa ni ustrezen v primeru zapisa vektorja glede na pomično bazo, glej (4.51) in (4.53). Zapis operatorja R s pripadajoco matriko R v referencni bazi (gg1, u2, U3) Ze poznamo: G11 G12 G13 R = G21 G22 G23 G31 G32 G33 Delovanje operatorja R na referencni ravnini lahko opišemo z 'reducirano' pravokotno matriko R A = G12 G13 G22 G23 G32 G33 (6.18) V matricne enacbe vpeljemo tudi standarden zapis matrike mehanskih vztrajnostnih momentov (glej [Crisfield, 1997] ali [Saje, 1994]) (6.19) " J11 J12 J13 J = J21 J22 J23 J31 J32 J33 s komponentami J = / P fc2 + £2) dV V Jj = — J dV = Jji, V (6.20) (6.21) za paroma razlicne i, j, k = 1, 2, 3. Oznaka 'V' pod integralskim znakom pomeni integracijo po obmocju V, ki ga zavzemajo vse tocke telesa v zacetni legi. V primeru nosilcev imamo opravka samo z ravninskimi mehanskimi vztrajnostnimi momenti. Mehanski vztrajnostni moment telesa je definiran z integracijo po obmocju nedeformiranega telesa (6.20)-(6.21), v teoriji nosilcev pa nastopa vztraj-nostni moment prereza, kije dvorazsezen objekt. Za izracun vztrajnostnih momentov prereza vzamemo okoli izbranega x obmocje nosilca, dolocenega z odsekom [x — dx, x + d^] in izracunamo tezišcni vztrajnostni moment telesa, ki ga doloca vektor r0D (x,£2,£3), kjer parametri (x,£2,£3) pretecejo obmocje [x — dx,x + dx] x Ar. V limiti dx ^ 0 dobimo pripadajoci mehanski vztrajnostni moment prereza. Posebej si oglejmo mehanske vztrajnostne momente (6.21) za i = 1: ( x+dx/2 \ «1 jj1 = j1j = j p£j 6 dv = j p £j j™0 j Vr A = / P £j lim dA 0 Ar 2 V dx/2 , dx/2 dA / = / p £j ,lim . (x + dx/2)2 - (x - dx/2)2' dx—> 0 dA Ar = / p£7- lim (x ■ dx) dA = 0, J j dx^0 Ar (6.22) za j = 2, 3, torej so mehanski vztrajnostni momenti J21, J31, J12, J13 enaki nic. Oglejmo si še diagonalne clene matrike (6.19), na primer clen J22: / x+dx/2 \ J22 = y p (£2 + £2) dV = J p dxm0 J £2d£1 dA + J P£2dV (6.23) Vr Ar \ x-dx/2 / V(x)r = / p Ji^ ^2x2dx + dA + J p£fdV = 0 + /33. Analogno izpeljemo J33 = /22. Vztrajnostni moment J11 imenujemo polarni ali torzijski vztrajnostni moment precnega prereza in ga zato oznacujemo z Jp. Jp = J11 in izracšunamo kot J11 = J p (£22 + £2) = /22 + /33. V Zunajdiagonalne clene matrike mehanskih vztrajnostnih momentov J, Jj za i = j, imenujemo devia-cijski vztrajnostni momenti prereza. Zaradi lastnosti (6.22) dobi standardna matrika mehanskih vztrajno-stnih momentov blocšno-diagonalno obliko (6.24) Za lazjo izpeljavo vpeljemo še pomozno reducirano matriko vztrajnostnih momentov, ki jo po lastnosti (6.23) povezemo tudi s standardnimi vztrajnostnimi momenti: t A = ■ J11 0 0 J= 0 J22 J23 0 J23 J33 /22 /23 J33 — J23 /23 /33 — J23 J22 Spomnimo se tudi matrike Hg, ki pripada antisimetricnemu operatorju Q v referencni bazi "g= 0 -U3 U2 U3 0 -Ui -u2 U1 0 (6.25) Tako imamo vse pripravljeno, da enacbi (6.13) in (6.16) za vrtilno kolicino in njen odvod zapišemo v matricni obliki: L o P Ar S (r ) Ug + (I22 + I33) V g - RAlARTA^g dx L O = [pAr (rg x Ug) + (I22 + I33) V g ] dx (6.26) (6.27) - [RtIaRTaVg + HgRaIaRA"g] dx. Vrtilno kolicino in njen odvod smo zapisali v kompaktni matricni obliki, v kateri pa ne nastopa standardni zapis matrike vztrajnostnih momentov (6.19), niti standardna rotacijska matrika. Podcrtana clena v (6.26) zdruzimo tako, da upoštevamo za (I22 + I33) vg razširjen zapis R (I22 +133) 0 0 0 (I22 +133) 0 0 0 (I22 +133) Rt V prav tako pa smemo zadnji clen -RaIaRA"g razširiti do oblike, ki temelji na standardni rotacijski matriki R R 0 0 0 0 ^33 -^23 0 -123 I22 RT Vg. Po upoštevanju teh razširitev v (6.26) se pojavi standardna matrika vztrajnostnih momentov (6.19), z njo pa je standardni zapis vrtilne kolicine glede na koordinatno izhodišce referencnega koordinatnega sistema tak LO = [pArrg xUg + RJR vg] dx. (6.28) Prvi clen je posledica zapisa vrtilne kolicine glede na koordinatno izhodišce O referencnega koordinatnega sistema. Podobno izpeljemo tudi klasicni zapis odvoda vrtilne kolicine. Izhajamo iz zapisa (6.27), kjer prav tako I33) V g in -RaI aRAv g RJRt V g zdruzimo clena (I22 + I33) Vg in -RaIaRAvg v kompaktno obliko Preostalemu clenu -HgRAlARTvg prištejemo nicelni clen Hg (I22 + I33) vg = (I22 + I33) vg x vg = 0, rezultat pa je standardni zapis odvoda vrtilne kolicine: L LO = J [pAr (rg x Ug>+ RJRTg + "gRJRT"g] dx. 0 (6.29) L L L L Opomba 14 Matrika vztrajnostnih momentov ima obliko (6.24) le v ortonormirani bazi, kjer je prvi bazni vektor pravokoten na ravnino prereza in ima koordinatni sistem izhodišče v težišču prereza. Za referenčni prerez Ar je ustrezna baza referenčna baza prostora (gi, ~g2, #3), medtem ko je za prereze A (x,t) deformiranega nosilca, kijih dobimo iz referenčnega prereza z rotacijo R (x,t) : Ar ^ A (x,t), ustrezna lokalna baza, kije pripeta na prerez z izhodiččem na teč^čni osi; tako bazo smo če vpeljali pod imenom pomična baza (ci (x, t), G2 (x, t), G3 (x, t) j = (R (x, t) ~gi,R (x, t) ]g2, R (x, t) ~g3). Ker velja predpostavka o nespremenljivosti prerezov, imata Ar in A (x,t) v omenjenih bazah enako matriko mehanskih vztrajnostnih momentov Jg ( Ar) = JG(x,t) (A (x,t)) , kar velja za poljubna x G [0, L] in t > 0. V izrazih (6.28)-(6.29) nastopa konstantna matrika J mehanskih vztrajnostnih momentov referenčnega prereza Ar zapisana glede na bazo Bg / J = Jg (Ar); toda v enačbi nastopa kompozitumu matrik RJRT, kar pa v tem primeru ne predstavlja koordinatne transformacije, temveč 'rotacijo' matrike J (mehanske vztrajnostne momente prereza Ar transformiramo v mehanske vztrajnostne momente zavrtenega prereza A (x, t)). Na ta način dobimo matriko vztrajnostnih momentov prereza A (x,t), ki pa je če vedno zapisana v bazi Bg / Jg (A(x,t)) = RJRT. Enačbo (6.28) bi torej lahko zapisali L Lo = /M'rg x"g + J g (A (x,t)) Cg] dx 0 kjer se Jg spreminja s časom in v spločnem nima oblike (6.24). Kadar potrebujemo vrtilno kolicino in njen odvod le za del nosilca od £ = 0 do x, lahko zapišemo x Lo (x,t) = j [p Ar rg (£,t) xug (£,t)+ R (£,t) JRT (£,t) Cg (£,t)] d£ 0 x L o (x, t) = j pAr (rg (£, t) x Ug (£, t)) 0 x + f [R (£,t) JRt (£,t) c (£,t) 0 + Hg (£, t) R (£, t) JRt (£, t) Cg (£, t)] d£. 6.3 Dopolnitev enacb prostorskega nosilca za dinamiko Vse enacbe in robne ter zacetne pogoje za prostorski nosilec moramo zapisati v obliki, primerni za dinamiko. Izhajamo iz enacb, ki smo jih pripravili v poglavju 6. Enacbe, ki povezujejo kinematične in deformacijske količine (5.11)-(5.12). Te enacbe ostanejo nespremenjene, le kolicine so odvisne tudi od casa t. Tako je rg (x, t) = R (x, t) (yg (x, t) — CG (x, t)) ^g (x, t) = T-1 (tfg (x, t)) R (kg (x, t) — dG (x, t)) x rg (x, t) = rg (0, t) + J R (£, t) (7g (£, t) — cg (£, t)) d£ 0 x tfg (x, t) = tfg (0, t) + j T-1 (tfg (£, t)) (»G (£, t) — dG (£, t)) d£. (6.30) (6.31) (6.32) (6.33) Variacijski konstanti (glede na relativno variacijo) uin d sta v splošnem funkciji x in t. i, je metoda ekspličitna, sičer je impličitna. Po zapisu, ki gaje za metode Runge-Kutta vpeljal Butčher (1987), tako metodo predstavimo s strukturo ci an ai2 ■ ■ ais c2 fl2i a22 ■ ■ ■ a2s cs asi as2 ■ ■ ass bi b2 ■ ■ bs Naravno število s imenujemo stopnja metode. Ker funkčijo F v vsaki vrstiči sheme (7.2) izvrednotimo največ s-krat, nam s pomeni tudi mero za kompleksnost metode. Konstante aij, bi in ci so določene tako, daje napaka rešitve, oznacimo jo z O (hp+1), cim manjša oziroma daje red metode p cim vecji. Konstante dolocimo z izenacenjem koeficientov istoleznih clenov Taylorjeve vrste okoli tocke t za tocno y (t + h) in priblizno Y (t + h) rešitev. Ker so tako dobljene enacbe za proste konstante v splošnem nelinearne, obstaja vec metod Runge-Kutta iste stopnje. Posebej omenimo le dve metodi, vgrajeni v komercialni program Matlab [The MathWorks, 1999], kiju uporabljamo v nadaljevanju. Metoda, imenovana ode23tb, temelji na implicitni Runge-Kutta metodi stopnje 3 (pravzaprav na dveh metodah razlicnih redov, vendar govorimo le o eni metodi, saj dvojnost sluzi le za oceno napake), ki jo natanko doloca struktura (7.3) 0 0 0 0 T d d 0 1 w w d w w d 1—w 3w+1 d 3 3 3 kjer so: t = 2 — \/2, d = t/2 in w = \/2/4. (Četrta vrstica v shemi (7.3) doloca metodo reda 2, zadnja vrstica pa metodo reda 3; razlika nam sluzi za preverjanje lokalne napake, upoštevamo pa rezultat po metodi višjega reda. Metodo (7.3) lahko razbijemo na dva koraka in kot taka je v literaturi bolj znana pod imenom TR-BDF2 [Hosea, Shampine, 1996] (TR je kratica za trapezno pravilo, angl. trapezoidal rule in BDF2 je kratica za obratno diferencno metodo drugega reda, angl. backward differentiation formula). Druzina veckoracnih metod BDFje bila ze leta 1952 kot prva predlagana za reševanje togih diferencialnih enacb in z nekaterimi razširitvami še vedno velja za najbolj uporabno druzino metod za reševanje togih problemov (glej [Hairer, Wanner, 1991] za podrobnosti). Kot pove ze ime samo, je metoda TR-BDF2 dvokoracšna; prvi (eksplicitni) korak je trapezno pravilo, drugi (implicitni) korak pa obratna diferencšna metoda 2. reda. Za izracune v obeh korakih uporabimo isto iterativno matriko. Le na kratko predstavimo osnovne ideje metode. Vzemimo korak velikosti h in tocko t ter oznacimo tT = t + Th in t1 = t + h. Predpostavimo, da pri t poznamo tocno rešitev y (t) zacetne naloge (7.1). Po trapeznem pravilu velja h YT = y + T-(F (t, y) + F (tT, YT)). YT oznacuje priblizek za vrednost y (tT). Oznacimo še z Y1 priblizek za y (t1) ter z ZT in Z1 priblizka za hy' (tT) in hy' (t1), medtem ko nam z pomeni tocno vrednost za hy' (t). Za prvi eksplicitni korak izberemo priblizšek Z° = h F (t, y). Za zacetni priblizek prvega implicitnega koraka vzamemo kar Z° = z, in za vsak iterativen korak k implicitnega dela pri tT izracšunamo Yk = (y + d z)+ d Zk (2. vrstica sheme 7.3). (7.4) Nato z Newtonovo iteracijo rešujemo enacbo Z k — hF (tT, Y k) = 0 tako, da enacbo lineariziramo po Z k: r , dF dY k r ,, T I — h---T ~ j — hdJ. dYk dZk rk k dYk Pri tem smo uporabili aproksimacijo (7.4) za linearizacijo Yk: = dJ, J pa oznacuje numericno d Z k T d F aproksimirano Jacobijevo matriko , kije enaka za vse iteracijske korake k. Priblizna Newtonova iteracija dobi obliko (I — hdJ)AZk = hF tT,Yk) — ZT z k+1 = Z k+aZ k. Opozorimo, da se tako izracunani Z T razlikuje od hF (tT, Y T), kar je bistveno za predstavljeno metodo. Zacetni priblizek Zi implicitnega dela pri t1 je dobljen s Hermitovo interpolacijo skozi tocki t in tT: Z1 = (i.5 + V2) z + (2.5 + 2^2) Zt - (6 + 4.5^2) (Yt - y). Za vsak korak k nato izracunamo Y1 = (y + w z + w Z t )+ d Z k (3. vrstica sheme 7.3). (7.5) Z1 izracunamo iterativno po poenostavljeni Newtonovi metodi (I - hdJ) AZk = hF (ti, Yk) - Z1 Z k+i = Z k+aZ k, k dYk kjer smo uporabili aproksimacijo (7.5) za linearizacijo Yk: = dI, matrika J pa je enaka kot za dZ k implicitni del pri tT. Druga uporabljena metoda ode45 programa Matlab [The MathWorks, 1999] je standardna eksplicitna Runge-Kutta metoda stopnje 5 [Dormand, Prince, 1980]. Primerna je za natancno reševanje netogih diferencialnih enacb. Enolicno jo doloca struktura 0 1/5 3/10 4/5 8/9 1 1 1/5 3/40 9/40 44/45 -56/15 32/9 19372/6561 -25360/2187 -212/729 9017/3168 -355/33 46732/5247 49/176 -5103/18656 35/384 0 500/1113 125/192 -2187/6784 11/84 5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40 35/384 0 500/1113 125/192 -2187/6784 11/84 0 Zadnji dve vrsti dolocata metodi redov 4 in 5, razlika pa sluzi za oceno lokalne napake metode. 7.2 Priprava glavnih in robnih enacb Kot smo ze omenili, lahko metode druzine Runge-Kutta uporabimo le na sistemih navadnih diferencialnih enacb 1. reda. Za uporabo vecine metod, vgrajenih v okolje Matlab (ode23tb, ode45,...), moramo sistem enacšb podati v obliki M (t, y) y = F (t, y), (7.6) kjer je M masna matrika, F je vektor desnih strani sistema diferencialnih enacb, y pa je vektor neznank. Najugodneje je, ce je masna matrika konstantna, v splošnem pa je lahko odvisna od casa in neznank, ki so prav tako funkcije casa: M (t, y) = M (t, y (t)). V nadaljevanju bomo ugotovili, da masna matrika enacb prostorskega nosilca vedno zavzame najsplošnejšo obliko. Glavne enacbe (6.82)-(6.83) in pripadajoci naravni robni pogoji (6.86)-(6.89) sestavljajo sistem navadnih diferencialnih enacb 2. reda po casu. Zato moramo te enacbe najprej preoblikovati v dvakrat toliksšen sistem diferenčialnih enačšb 1. reda po čšasu. To storimo na standarden načšin. Za prve odvode neznank po času v točkah xi, i = 0,1,..., N + 1, uvedemo nove oznake vg (t)= ug (t) Al (t) = k (t), (7.7) (7.8) in jih dodamo naboru osnovnih neznank; tako je novi vektor vseh neznank y u0, v0, k, A0, <, < k, A0, u2, u oo .i „.i -i -i N+i „,N+i kN+i kN+i gg g' g' g gg , v k , A0 T (7.9) V momentni ravnotezni enačbi (6.71) kotno hitrost (kg in pospešek kg izrazimo z osnovnimi spremenljivkami in dobimo k6RK : k6S - 2k ◦ J (5* ◦ Ag + k* ◦ Ag ◦ Ag ◦ kJJ ◦ k - 4SfAg O k*) fk ◦ J fk* ◦ Ag ◦ k*l = k, kjer pomenita oznaki Ag = Ag o k0 + kg ◦ k0 ter Ag = Ag o k0 + 2Ag o k0 + kg o k0 prvi in drugi odvod kvaterniona čelotne rotačije (4.36). Glavne enačbe (6.82)-(6.83) nadomesti sistem diferenčialnih enačb 1. reda, zapisan glede na kolokačijske točke Xk, k = 1, 2,..., N: N+i N+i f ft : £ U g (t) Li (xfc )= £ vg (t) Li (xfc) i=0 i=0 N+i f i2 : pAr £ vg (t) Li (Xk) i=0 kk' ockk o kk *+kk oCkk o kk *+kk ocn ◦ kk *' N+i . . N+i f 2; : £ k (t) Pi (Xk) = £ Ag (t) Pi (Xk) i=0 i=0 JR3 + nk k22 : 20l (kk) 0R (kfc *) k (kfc *) 0R (^0) £ Ag (t) Pi (Xk) i=0 kk' O ckM O kk * + kk O cM o kk * + kk O ckM O kk *' + mk - S fkk oCN O kk* 1 kk' N+i - 2kk o ( j (2kk* o Ag o k0+kk o kk+kk* ◦ Ag ◦ Ag ◦ k i k k -jj k k kk - 4S (Ak o kk* kk O J (kk* O a; k ◦k k (7.10) (7.11) (7.12) (7.13) Enačbe sistema smo zapisali v obliki (7.6), primerni za uporabo v programskem okolju Matlab. Podobno zapišemo tudi robne enačbe (6.86)-(6.89). Odvisnost osnovnih spremenljivk od naravnega parametra x, ki preteče območje integračije po osi nosilča, bomo zaradi predolgih izrazov in nepreglednosti izpustili. Ohranimo pa zgornje indekse k za količine, izvrednotene v kolokačijskih točkah in L/2 i- k za kolicine, izvrednotene na sredini tezišcne osi nosilca. V levem krajišcu nosilca tako dobimo štiri diferencialne enacšbe 1. reda: N+i N+i k hii : £ uj (t) Li (0) = ^ vg (t) Li (0) (7.14) i=0 i=0 L/2 L/2 N+i ^ hi2 : pAr J Li (x) dx vg (t) = S0 + NL/2 + J ngdx i=0 0 i=0 0 L/2 L/2 = P0 + dL/2 O ^M/2 O dL/2* + f mg dx -j S (d O (Tn o d *) dg dx O C ot i i " ^g d 00 L/2 - 2 J d O (J (2d * O Ag O d0 + d0 O d0 + d * O Ag O Ag O d) ) od *dx 0 L/2 - 4 / d (Ag O d*' (dk O (j (dk* O Ag'' O d*' dx in še štiri v desnem krajišcu: N+i h32 : pAr E / Li (x) dx vg (t) = SL - N^/2 + ng dx i=0 N+i N+i i ig g i=0 i=0 L • * - 2 j d O ( J (2d O Ag O d0 + d0 O d0 + d O Ag O Ag O j) ) od dx L/2 (7.15) N+i • k N+i k h2i : E k (t) Pi (0) = E Ak (t) Pi (0) (7.16) i=0 i=0 N+i L/2 • h22 : E / 20l (d) (d*) j0l (d*) (d0) Pi (x) dxAg (t) (7.17) N+i N+i h3i : E ug (t) Li (L) = £ vg (t) Li (L) (7.18) i=0 i=0 L L (7.19) L/2 L/2 N+i h4i : £ k (t) Pi (L) = E Ag (t) Pi (L) (7.20) N +i h42 : E J 2^l (d) (d*) j 0l (d*) (d0) Pi (x) dxAg (t) (7.21) i=0 L/2 L L = PL - dL/2 o^M/2 O dL/2 * + f mg dx -j j(d o (Tn o j*) j, dx L/2 L/2 L L - 4 y S (Ag o j*) (jk ◦ (j (jk * o Ag)) o j*) L/2 Integrale v robnih enacbah (7.14)-(7.21) izvrednotimo numericno s klasicno Gaussovo kvadraturno metodo izbrane stopnje. Metodo smo predstavili v poglavju 5.3.5. Opomba 16 Metodi Runge-Kutta, predstavljeni v poglavju 7.1, sta razviti za integracijo v aditivnih kon-figuracijskih prostorih in pri popravljanju osnovnih neznank ne upoštevata multiplikativne narave rotacij. Zato so rotacijski kvaternioni po zaključenem popravljanju v splošnem neenotski. Da odstopanja od enotske norme niso prevelika do določene mere ureja preverjanje lokalne napake in posledično prilagajanje časovnega koraka. Poleg tega po zaključku vsakega časovnega koraka kvaternione, ki določajo rotacije normiramo, predvsem v izogib prehajanja v kompleksna čtevila v primeru, ko norma rotacijskega kvaterniona preseče vrednost ena. V primeru, ko v posamičnem časovnem koraku rotacijski kvaternioni bistveno ostopajo od enotske norme, moramo poseči po drugačnih ukrepih: čtirikomponentne enačbe reduciramo na trikomponentne (vektorski del kvaternionskih enačb), vektor vseh neznank pa zmanjčamo za prve komponente rotacijskih kvaternionov in njihovih odvodov. Rečitev skrčenega sistema po zaključku vsakega časovnega koraka ponovno razčirimo s prvimi komponentami rotacijskih kvaternionov po (4.23). 7.3 Numericni testi Predstavljeni postopek testiramo numericno s primerjanjem rezultatov znacilnih prostorskih primerov z rezultati literature. Numericni rezultati so v celoti pridobljeni z lastnim racunalniškim programom v programskem okolju Matlab [The MathWorks, 1999]. Pri izracunih se omejimo na linearno elasticen materialni model (5.91)-(5.92). Za vse primere uporabljamo metodo ode23tb programskega okolja Matlab 7.3.0 (R2006b), primerno za reševanje togih in zelo togih primerov. Le za primer nepodprtega nosilca v prostem gibanju za integracijo na daljšem casovnem intervalu uporabimo metodo višjega reda ode45. Predstavljeni so rezultati za pomike in rotacijske vektorje, ki smo jih za namene primerjave z drugimi avtorji preracunali iz rotacijskih kvaternionov z enacbo (5.93). Spiralno gibanje upogibnega nosilca. Primer prostorskega gibanja enostransko drsno-clenkasto vpetega nosilca, ki ga iz mirovanja s primerno silo in momentom sprozšimo v spiralno gibanje, sta predstavila Ibrahimbegovic in al Mikdad (1998). Za primer so znacilni zelo veliki zasuki (vecji od 2n) in izrazite globalna in lokalna nihanja. Geometrijski in materialni podatki so povzeti po [Ibrahimbegovic, al Mikdad, 1998]: L = 10 Ap = 1 3 = 104 EJ 20 0 0 0 10 0 0 0 10 Nosilec lezi v smeri osi X (oziroma v smeri <71), levo krajišce je podprto drsno in clenkasto tako, da sta dopušcena pomik v smeri Z in zasuk okoli osi Z, glej sliko 7.1. Zaradi takih robnih pogojev so nekatere komponente pomika in rotacijskega kvaterniona levega krajišca enake nic: ui = u0 = 0 kg — kg — 0. Slika 7.1: Geometrija in obteZba upogibnega nosilca Figure 7.1: Force driven flexible beam: geometry and loading data Spomnimo se, da zgornji indeks (0) oznacuje pogoj v levem krajišcu pri x = 0, spodnji indeksi pa dolocajo komponento ustreznega vektorja ali kvaterniona. (Čeprav je v levem krajišcu dopušcen zasuk prereza zgolj okoli ene osi, pa sta kar dve komponenti rotacijskega kvaterniona pri x = 0 brez pogoja, k0 in k0. Ker so tudi prvi odvodi po casu dodani seznamu osnovnih neznank, moramo v levem krajišcu podati tudi robne pogoje za hitrosti oziroma odvode rotacijskih kvaternionov po casu. V levem krajišcu dopustimo le hitrost v smeri Z in kotno hitrost (oziroma nenicelen odvod rotacijskega kvaterniona) okoli osi Z, torej v0 = v0 = 0 A0 = A0 = 0. Zacetne vrednosti vseh neznank so nic, razen seveda vrednosti prve komponente rotacijskih kvaternionov, ki so enake ena: [0] ,[0] ■u = u [0] [0] Al [0] [0] = k [0] [0] = (Ai0]) = (A =u [0] =k [0] [0] 0 ki01) * = i «n *=o o, za vse indekse i, ki dolocajo robni tocki x0 = 0 in xn+i = L, kolokacijske tocke x* = 1, 2,.., N in sredinsko tocko x = L. Gibanje nosilca sprozimo s socasno obtezbo v levem krajišcu s konstantnim momentom Mz = 80 okoli osi Z in s konstantno tockovno silo Fz = 0.05 Mz v smeri osi Z. Obe obtezbi prenehata delovati pri casu t = 2.5 (slika 7.1). Gibanje nosilca opazujemo do casa t = 50. Prav tako kot Ibrahimbegovic in al Mikdad smo uporabili mrezo desetih linearnih elementov. Podatki numericšnega racšuna o zahtevanih tolerancah so: ^abs 10 -3 ^rel [10 -3 10-3,10-1,10 -11 2 3 0 1 1 2 0 3 kjer se relativne tolerance zaporedoma nanašajo na pomike, rotacijske kvaternione, hitrosti in odvode rotacijskega kvaterniona. Za pomike in rotacijske kvaternione smo zahtevali višjo relativno natancnost Slika 7.2: Pomiki prostega krajišca v odvisnosti od casa za primer upogibnega nosilca. 10 linearnih elementov; primerjava z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) Figure 7.2: Free-end displacement of force driven flexible beam. 10 linear elements; comparison with Ibrahimbegovic and al Mikdad (1998) (dashed-line) Slika 7.3: Pomiki prostega krajišca v odvisnosti od casa za primer upogibnega nosilca. 2 linearna elementa; primerjava z rezultati Ibrahimbegovica in al Mikdada (1998) (crtkana crta) Figure 7.3: Free-end displacement of force driven flexible beam. 2 linear elements; comparison with Ibrahimbegovic and al Mikdad (dashed-line) izraCuna kot za njihove odvode, saj s tem pospešimo racun. Ibrahimbegovic in al Mikdad (1998) ne navajata podatkov o natancnosti racuna. Gibanje konstrukcije zajema drsenje vzdolz osi na dopušceni smeri Z in socasno vrtenje okoli osi Z za skoraj dva obrata. Pri tem nosilec niha v vseh treh smereh. Slika 7.2 spreminjanja pomikov prostega krajišca s casom ima podobno obliko, kot jo prikazujeta Ibrahimbegovic in al Mikdad, amplitude nihanja vseh treh vozlišcnih pomikov dosegajo iste vrednosti, nihajna casa pomikov ui in u2 pa sta nekoliko manjša. Opazne so drobne oscilacije na grafih vseh pomikov, medtem ko oscilacij za pomike u1 in u2 pri Ibrahimbegovicu in al Mikdadu (1998) ni opaziti. Ibrahimbegovic in al Mikdad sta uporabila konstanten casovni korak velikosti 0.5, kar znese 100 casovnih korakov, medtem ko je korak tu uporabljene metode prilagodljiv in znaša od 0.1302 do 0.0232, kar znese 576 casovnih korakov. Zaradi vecjega števila korakov in preverjanja ter omejevanja absolutne in relativne napake predvidevamo, daje naš rezultat tocnejši. Ker na dolzino casovnega koraka ne moremo vplivati, manj natancen izracun rezultatov dobimo z mrezo dveh elementov in z manj strogimi pogoji za absolutno in relativno napako: eabs = 10-2, £rei = [10-1,10-1, 10-1,10-1]. V tem primeru se grafi pomikov u1 in u2 prostega krajišca tako dobro ujemajo, da so odstopanja z graficno natancnostjo neopazna. Zares se nihajna casa pomikov u1 in u2 povecata in do graficne natancnosti ujameta z rezultatom iz literature, glej sliko 7.3. Ocitno pa model nosilca s samo dvema elementoma ni v zadostni meri sposoben zaznati drobnih lokalnih oscilacij pomikov v nobeni smeri. Kolenasta konzola z obtežbo pravokotno na ravnino. Primer pravokotno prepognjene konzole, katere geometrija je predstavljena na sliki 7.4, sta obravnavala Simo in Vu-Quoc (1988). Ostali podatki so: L = 10 GA = = 106 " 20 0 0 0 10 0 0 0 10 Ap = 1 EJ = GJ = 103 1 12 t B Slika 7.4: Geometrija in obtezba kolenaste konzole Figure 7.4: Right-angle cantilever beam: geometry and loading data Kolenasta konzola je v pregibu A obtezena s tockovno silo, ki deluje pravokotno na ravnino okvirja. Do casa t = 1 obtezba linearno narašca do vrednosti 50 in nato linearno pada do nicelne vrednosti pri casu t = 2. Po casu t = 2 je konstrukcija neobremenjena. Nadaljnje gibanje konzoleje prosto nihanje, ki vkljucuje znaten upogib in torzijo. Spremljamo ga vse do casa t = 30. Konstrukcija doseze velike pomike velikostnega reda dolzšin obeh nosilcev. «3 (A) - «3 (B) == [Simo, Vu-Quoc, 1988], 2 elementa [Simo, Vu-Quoc, 1988], 10 elementov 30 Slika 7.5: Odvisnost pomikov od casa za koleno (A) in za prosto krajišce (B) kolenaste konzole; sive crte prikazujejo rezultate Sima in Vu-Quoca (1988) za primer dveh in desetih elementov Figure 7.5: The right-angle cantilever beam: time histories of elbow (A) and tip (B) displacements; comparison with Simo and Vu-Quoc (1988) (grey lines) employing 2 and 10 elements Simo in Vu-Quoc racun izvedeta za dve mrezi, eno z dvema in drugo z desetimi kvadraticnimi elementi ter pri konstantnem casovnem koraku 0.25. Metoda ode23tb s tako majhnim številom elementov ne more doseci smiselno majhne relativne oziroma absolutne napake. (Če zahtevano natancnost sprostimo, racun sploh ni izvedljiv. Kombinacija mreze 12 elementov z dvema internima tockama (kubicen element) za vsak del konstrukcije dolzine 10 in numericno integracijo cetrtega reda (3 integracijske tocke) ter s pogojema za lokalno napako Sabs = 10-5 Srel = [10-5, 10-5, 10-1, 10-1] pa zmore analizo vse do casa t = 30. (Časovni grafi pomika kolena in prostega krajišca, prikazani na sliki 7.5, se do casa t = 20 dobro ujemajo z grafi iz literature [Simo, Vu-Quoc, 1988], kasneje pa so razlike izrazitejše. Vecje razlike po casu t = 20 lahko opazimo tudi ob primerjavi rezultatov razlicnih avtorjev iz literature, glej na primer še [Ibrahimbegovic, al Mikdad, 1998]. Velike razlike med razlicnimi numericnimi postopki reševanja raznih avtorjev se pokazejo v številu casovnih korakov in casovni ter prostorski zahtevnosti integratorja. Medtem ko [Ibrahimbegovic, al Mikdad, 1998] in [Simo, Vu-Quoc, 1988] uporabita le 126 oziroma 54 prostostnih stopenj, za koncni rezultat pa potrebujeta le 120 casovnih korakov, pa predstavljeni pristop zahteva 1535 casovnih korakov velikosti od 5.3350 ■ 10-5 do 0.04012, v vsakem koraku pa kar 504 prostostnih stopenj za pomike in rotacijske kvaternione. Vendar je v primerjavi potrebno upoštevati, da omenjena avtorja v vsakem koraku iterativno rešujeta sistem nelinearnih enacb, kar zahteva sestavljanje tangentne togostne matrike v vsakem iterativnem koraku, medtem ko pri uporabi metode ode23tb tangentno togostno matriko (oziroma le njen priblizek) v vsakem casovnem koraku izracunamo le enkrat (glej poglavje 7.1). Zaradi velikih pomikov se izjemno hitro spreminjajo tudi naklonski koti precnih prerezov nosilca. Zato pri uporabi metode ode23tb, kije zasnovana za integracijo v linearnih prostorih in ne ohranja narave rotacij oziroma enotskosti rotacijskih kvaternionov, prihaja v primeru upoštevanja vseh štirih kvaternionskih enacb do numericnega neskladja. Norma rotacijskih kvaternionov se zelo hitro oddalji od vrednosti 1, nakar uporabljena metoda ni vec sposobna zagotoviti zadostne natancnosti racuna. Zato smo za primer kolenaste konzole rešitev izracunali z redukcijo kvaternionskih enacb na tri komponente tako, da smo prve komponente kvaternionov k0 (xi, t) izlocili iz vektorja neznank in jih racunali posredno iz k (xi, t) z vezno enacbo (4.23). Opomba 17 Pri uporabi postopka redukcije štirih kvaternionskih enačb na tri (vektorski del) enačbe s sočasno redukcijo čtirih kvaternionskih neznank na tri (vektorski del) smo se izognili numeričnim tečavam, ki so posledica neupočtevanja narave rotacij pri popravljanju rotacijskih količin. Poleg omenjene redukcije, s katero smo dosegli zgolj ohranjanje lastnosti rotacijskega kvaterniona (4.23), bi lahko ohranjanje enotske norme rotacijskega kvaterniona dosegli tudi s popravljanjem uporabljenega inte-gratorja oziroma z uporabo nekega drugega integratorja, ki upočteva naravo rotacij. Popravljanje obstoječih komercialnih integratorjev ni namen pričujočega dela, medtem ko uporabo drugačnega integratorja izvedemo v nadaljevanju (glej poglavje 9). 8 t M3(t) F(t)=M2(t)/10 M3(t)/10=M2(t)/2 6 Slika 7.6: Geometrija in obtezba nepodprtega nosilca Figure 7.6: Free-free flexible beam: geometry and loading data Nepodprt nosilec v prostem gibanju. Zadnji primer obravnava gibanje popolnoma nepodprtega nosilca, ki ga po zacetni obtezbi prepustimo prostemu gibanju v prostoru. Podatki o obtezbi in geometriji nosilca so predstavljeni na sliki 7.6, preostali podatki pa so: L = 10 Ap = 1 = 104 EJ 10 0 0 0 10 0 0 0 10 EJ = GJ = 500 Prostorski primer nepodprtega nosilca v prostem gibanju sta študirala Ze Simo in Vu-Quoc (1988). Izbrali smo mreZo 10 linearnih elementov. Natancnost racunaje dolocena z zahtevami: ^abs = 10-3 £rel = [10-3, 10-3, 10-1, 10-1] . V zacetnem stanju lezi nosilec v ravnini, ki jo dolocata bazna vektorja gg1 in gg3. V desnem krajišcu je socasno obtezen s tockovno silo F1 v smeri gg1 in s tockovnima momentoma M2 in M3 okoli osi, dolocenih z gg2 in gg3. Vse obtezbe od nicelne vrednosti pri casu t = 0 dalje linearno narašcajo do casa t = 2.5, nato pa se linearno zmanjšujejo do casa t = 5, ko postanejo enake nic; prosto gibanje konstrukcije opazujemo vse do casa t = 15. Slika 7.7 prikazuje znacilne zaporedne deformirane lege nosilca do casa t = 5, projecirane na ravnini XZ in , slika 7.8 pa aksonometricni prikaz gibanja do casa t = 15. 8 6 Z 4 2 0 0 2468 10 X 12 14 16 18 8 6 Z 4 2 0 42 0 -2 -4 Y Slika 7.7: Ravninske projekcije poteka gibanja nepodprtega nosilca do casa t = 5 Figure 7.7: Projections of the deformed shapes on the coordinate planes for the free-free flexible beam to time t = 5 Rezultat gibanja je primerljiv z rezultatom iz literature [Simo, Vu-Quoc, 1988], prav tako je primerljiv velikostni red števila casovnih korakov; Simo in Vu-Quoc uporabita fiksen casovni korak velikosti 0.1, kar pomeni skupno 150 korakov, medtem ko je bilo s prilagodljivim korakom (velikosti med pribliZno 0.208 in 0.019) metode Runge-Kutta potrebnih skupno 190 casovnih korakov. Slika 7.8: Aksonometricni prikaz zaporedja deformiranih leg nepodprtega nosilca pri znacilnih casih na casovnem obmocju [2.5,15.5] Figure 7.8: Free-free flexible beam. Perspective view of deformed shapes on time interval [2.5,15.5] V [Simo et al., 1995] isti primer nepodprtega nosilca v prostem gibanju avtorji uporabijo za demonstracijo uspešnosti metode, ki pri konzervativnih primerih ohranja energijo, gibalno in vrtilno kolicino, kar je zadosten pogoj za dolgotrajno stabilnost numericnega racuna. Kadar za dolgotrajen racun uporabljamo metodo Runge-Kutta, je smiselno izbrati metodo cim višjega reda z visoko natancnostjo ob strogi zahtevi za lokalno napako, da preprecimo globalno kopicenje numericne napake racuna. Zato je za racun na dolgem casovnem intervalu [0, 500] primerna klasicna metoda Runge-Kutta reda 5, kije v Matlabu imenovana ode45. Za absolutno in relativno napako tokrat postavimo strozje zahteve: £abs = 10-5 £rel = [l0-5, 10-5, 10-2, 10-2] . Metoda višjega reda praviloma zazna tudi majhne oscilacije, zato postane ob strogi zahtevi za natancnost dolzina koraka razmeroma majhna in v obravnavanem primeru vecinoma velikostnega reda 10-4; zato je racun zelo dolgotrajen (v obravnavanem primeru okrog 36 ur na PC). Vendar racun tudi po vec kot dveh milijonih korakov ostaja stabilen (velikost koraka se po ustalitvi skoraj ne spreminja vec) in brez tezav doseze cas t = 500. 8 Reševanje dinamičnih enačb s posplošeno metodo Newmark Sodobne integracijske metode, kijih v zadnjih dveh desetletjih številni raziskovalci razvijajo na podrocju dinamike prostorskih nosilcev ob uporabi rotacijskega vektorja in z upoštevanjem nelinearne narave prostorskih rotacij, lahko uspešno priredimo tudi za izpeljano kvaternionsko obliko enacb. Med metodami je najpogosteje uporabljena druzina Newmarkovih metod, ki so namenjene prav reševanju sistema diferencialnih enacb 2. reda. Metode so predstavljene v mnogih ucbenikih, ki se podrobneje ukvarjajo z metodo koncnih elementov, na primer v [Bathe, 1996], [Belytschko et al., 2000]. (Časovni integrator posplošene Newmarkove metode, ki upošteva naravo rotacij, za parametrizacijo rotacij pa uporablja rotacijski vektor, sta razvila Simo in Vu-Quoc (1988). Naš cilj je, ta integrator prirediti tudi za kvaternionsko parametrizacijo rotacij. 8.1 Izpeljava posploSene metode Newmark za kvaternionske enačbe Sprva nakazemo kratko izpeljavo Newmarkove metode, kot je predstavljena in implementirana v [Simo, Vu-Quoc, 1988]. Izhajamo iz znacilnosti razvoja funkcije v Taylorjevo vrsto, povzetih po ucbe-niku [Krizanic, 1990]. Zvezno rešitev y = f (x) diferencialne enacbe drugega reda y" = F (x, y, y') razvijemo v Taylorjevo vrsto okrog tocke a y = f (x) = f (a) + (x - a) f' (a) + ^^ f " (a) + ... + (x--aP f(p) (a) + .... (8.1) 2 p. (Če vrsti (8.1) odvzamemo zacetne clene do vkljucno p-tega odvoda, preostanek vrste imenujemo ostanek Taylorjeve vrste in oznacimo s Sp OO Sp = £ ^^ f(i) (a). i=p+1 Za ostanek Taylorjeve vrste velja, da z vecanjem števila p konvergira proti nic, kar lahko zapišemo z razliko prave funkcijske vrednosti in njene aproksimacije s koncno Taylorjevo vrsto do vkljucno p-tega cšlena Sp = f (x) - >T fcj^ f(a) 0. i=1 (Če je funkcija f (x) v okolici tocke x (p + 1)-krat zvezno odvedljiva in ima odvod f (p+1) povsod omejen, potem je ostanek Taylorjeve vrste Sp enak vrednosti predhodnega (zadnjega odvzetega) clena vrste v neki tocki C, ki leZi med x in a [KriZanic, 1990, str. 502] Sp = (XPTi)~ f(P+1) (e) • Tako funkcijo lahko torej tocno zapišemo s koncno vrsto ( (p + 1) f (x) = £ f(i) (a) + f (p+1) (C). (8.2) i=1 Newmarkova integracijska shema za krajevne vektorje. Tocen razvoj funkcije v vrsto (8.2) uporabimo za zapis krajevnih vektorjev tezišcne osi nosilca tako, da izberemo p = 2: rg (tn + At) = rg (tn) + At r g (tn) + — rfl (C) • At2 .. T V izrazu nastopa prvi in drugi odvod krajevnega vektorja po casu; ker pomenita hitrost in pospešek tocke na osi, jima dodelimo oznaki Torej velja tocen izraz At2 rg (tn + At) = rg (tn) + At Vg (tn) + ag (C) Vrednosti C v splošnem ne poznamo, lezi pa med tn in tn + At. Predpostavimo, da je pospešek v tocki (casu) tn zvezna funkcija t. Potem je na dovolj majhni okolici tocke tn konstantna, narašcajoca ali padajoca funkcija casa. V primeru konstantnega pospeška je ustrezen prav vsak C £ [tn,tn + At]. V primeru narašcajocega ali padajocega pospeška pa lahko pospešek v tocki C zapišemo kot linearno kombinacijo robnih vrednosti ag (C) = (1 - b) ag (tn) + b ag (tn + At) (8.3) za nek skalar b G [0,1]. Formula (8.3) velja tudi za konstantno vrednost pospeška (1 - b) ag + b ag = ag = ag (C) • Izraz (8.3) je tocen za tako velikost At, daje pospešek na intervalu [tn, tn + At] monoton ali konstanten. Z vrsto (8.2) smemo aproksimirati tudi hitrosti. Za aproksimacijo hitrosti izberemo p = 1, skalar linearne kombinacije pospeškov pa oznacimo z y G [0,1]. Na koncu skalar b G [0,1] nadomestimo s polovicnim skalarjem P G [0,1], P = | in tako dobimo klasicno obliko Newmarkove aproksimacije pomikov in hitrosti pri casu tn+1 = tn + At: rgn+1] = rgn] + At vgn] + At2 2 - p) agn] + P agn+1] (8.4) vgn+1] = vgn] + At [ (1 - y) agn] + y agn+1^ , (8.5) za skalarja P G [0,1] in 7 G [0,1]. (Čas, na katerega se kolicina nanaša, je oznacen z zgornjim indeksom v oglatem oklepaju. Razlika rgn+1] — rgn] = rog + ugn+1] — ( rog + ugn]) je enaka spremembi pomika rg = vg rg = ag Au [n] u [n+1] — u [n] in je, zato ker so krajevni vektorji in pomiki aditivne vektorske količine, enaka [n] linearnemu delu spremembe pomika, Aug spremembe pomika zapišemo kot 5ugn] =At vgn] + At2 5u. n] Z upoštevanjem teh zvez v (8.4) linearni del 2 - * a [n] + *agn+1] (8.6) Newmarkova integracijska shema za rotacijske matrike, parametrizirane z rotacijskim vektorjem. Upoštevati moramo, da se sprememba rotačijske matrike iz časa pri stanju n v stanje n + 1 doda multiplikativno: R[n+1] = arN R[n] = exp 5©!n] i R[n] (8.7) kjer je 50gn] antisimetrična matrika z osnim vektorjem 5$^; zgornji indeks [n] pomeni popravek iz stanja pri času tn v čas tn+1. Količina 5$gn] ne predstavlja trenutnega iterativnega Newtonovega popravka, saj v enačbe še nismo vpeljali iterativnega pristopa popravljanja količin. Zaradi aditivnosti nekaterih z rotačijo povezanih količin, kadar so izrazene v pomični bazi, pri izpeljavi preidemo na pomično bazo; za podrobnosti glej disertačijo D. Zupana (2003), ali pa analogno izpeljavo aditivnosti za vektorje ukrivljenosti (5.86) iz poglavja 5.3. Tudi linearni del spremembe, 5©{n] oziroma 5${n], aproksimiramo v skladu z ze izpeljano Newmarkovo aproksimačijo za pomike (8.6) z vsoto prispevka kotne hitrosti v [n] predhodnem času uG in prispevka linearne kombinačije kotnih pospeškov v prejšnjem ter trenutnem [n] [n] [n+1] 1GI"+11' 5$ [n] G[n+11 At-Gti+At2 2 - * a [n] G[n] + * a[n+1] T * aG[n+1] Aproksimačija kotne hitrosti je enako oblikovana kot aproksimačija hitrosti (8.5): [n+1] G[n+11 Ui+At (1 - Y) aGt + Y aGl+++]11 [n+1] (8.8) (8.9) (Četudi skalarjev * in y ne določimo tako, da bi enačbi (8.8)-(8.9) predstavljali točna izraza za popravek rotačije in kotne hitrosti, pa predstavlja R[n+1] = exp ^50gn] (*)) R[n] lokalno aproksimačijo operatorja R[n+1] = exp ^50gn]) R[n] tretjega reda; dokaz najdemo v [Simo, Vu-Quoč, 1988]. Zaradi nelinearnosti diskretnih enačšb prostorskega nosilča za dinamiko jih je potrebno resševati iterativno. Izpeljane impličitne aproksimačije (8.4)-(8.9) se nanašajo na konstanten časovni korak. Z vpeljavo ločevanja med dvema zaporednima iterativnima stanjema i in i+1 v istem časovnem koraku pa izpeljemo ustrezne iterativne izraze, ki pa zavzamejo ekspličitno obliko. Kratko nakazšimo potek izpeljave le za kotne količine. Zaradi mnoziče različnih stanj je pomembno točno opredeliti količino in pripadajočo bazo, v kateri je izrazšena; zato v nadaljevanju sistematičšno navajamo pomene indeksnih oznak: MgL je količina (*) ob zaključku stanja v času tn, izrazeno v pomični bazi zaključnega stanja v čšasu tn, G[n]; (*)P [n] [n + 1] pomeni spremembo količine (*) od stanja v času tn do stanja v i-ti iteračiji ob času tn+1, izrazeno v predhodni iteračijski pomični bazi G, [n+1]. (*)P [n] [n+1] i+1,Gi++ pomeni spremembo količine (*) od stanja v času tn do stanja v i + 1-vi iteračiji ob času tn+1, izrazeno v (trenutni) pomični bazi G [n+1]; i+1 ; [n+1] [n + 1] G [n+1]. je kolicina (*) v i-ti iteraciji ob casu tn+1, izraZeno v predhodni iteracijski pomicni bazi (*)[n+1] [n+1, je kolicino (*) v i + 1-vi iteraciji ob casu tn+1, izrazeno v trenutni iteracijski pomicni i+1,Gi+1 bazi G [n+1] i+1 . Trenutne linearne popravke neznank (iz iteracije i v iteracijo i + 1 v casu tn+1) oznacujemo zgolj z indeksom baze, na primer in čtfg. Aproksimacijo spremembe rotacije iz casa tn do casa tn+1 za dve zaporedni iteracijski stanji i in i + 1 zapišemo v skladu z enacbo (8.8) in vpeljanimi oznakami: M tf [n] i,G. [n+1] M tf [n] i+1,G. [n+1] +1 At ^Gln]+At2 At ^G]n]+At2 2 - * 2 - * a a [n] G[n] [n] G[n] + * a [n+1] i,G [n+1] + * a [n+1] i+1,G [n+1] +1 (8.10) Enacbi medsebojno odštejemo in uredimo a [n+1] i+1,Gj++1] a[n+1] + a [n] + *At2 čtf [n] i+1,Gl++1] -čtf [n] i,G [n+1] (8.11) Enacba (8.11) predstavlja eksplicitni postopek za izracun priblizka kotnega pospeška iz iteracije i v itera- cijo i + 1 v casu tn+1. čtf [n] i+1,G [n+1] +1 predstavlja spremembo rotacijskega vektorja od cšasa tn do iteracije i + 1 v casu tn+1, izrazeno v trenutni bazi in ni enaka trenutnemu linearnemu popravku rotacijskega vektorja V i + 1-vi iteraciji ob casu tn+1 se izracuna iz predhodne spremembe rotacijskega vektorja [ n] od casa tn do i-te iteracije ob casu tn+1, izrazene v predhodni bazi, to je iz £tf[ ] [ +1], z multiplikativnim i,Gi + ] dodajanjem trenutne linearne spremembe £tf [n+1,, izrazene v trenutni bazi: G +1 ^tf[n] ,[n+1] = exp-1 i+1,G' +1 exp I čtf [n] i,G, [n+1] exp ^tfG[n + 1] Gi+1 (8.12) Kot je bilo omenjeno, trenutna linearna sprememba £tf [n+1, nima zgornjega indeksa. Oznaka exp G 1 +1 izrazu (8.12) je inverz eksponentne preslikave matricnega argumenta R in oznacuje postopek (in ne preslikave, saj postopek ni enolicen) dolocitve osnega vektorja tf antisimetricnega operatorja 0, ki pripada dani rotaciji R = exp(0). Ker pri implementaciji uporabljamo zgolj kvaternionsko parametrizacijo rotacije, se numericni izvedbi tega inverza podrobno ne posvetimo. Primer numericno stabilnega algoritma za inverz eksponentne preslikave matricnega argumenta je Spurrierov algoritem [Spurrier, 1978]; njegova uporaba je na podrocju tocnih teorij prostorskih nosilcev zelo razširjena, vpeljal pa jo je Simo z razlicnimi sodelavci (1986, 1988, 1991). Ob locevanju dveh zaporednih iterativnih stanj v casu tn+1 iz enacbe (8.9) sledita izraza za koncne kotne hitrosti u [n+1] [n+1] u i,G [n+1] i+1,Gi++1] u£L+At u£L+At n_ Y) a[n] + Ya[n+1] (1 Y) aG[n] + Y aiG[n + 1] (1 - Y) a" + ya£G 1 v Izraza odštejemo in uredimo [n+1] [n+1] [n+1] [n+1] Z uporabo enacbe (8.11) nadaljujemo s preoblikovanjem tako, da pospeške nadomestimo s spremembami zasukov: [n+1] i+1,G. [n+1] + 1 [n+1] [n+1] + Y pAt 5? [n] i+1,G. [n+1] +1 -5? [n] i,G, [n+1] (8.13) Dobili smo eksplicitno iterativno izrazavo kotne hitrosti v iteraciji i + 1 s kotno hitrostjo iz predhodne iteracije i. Newmarkova integracijska shema za rotacijski kvaternion. Koncno izpeljemo še Newmarkovo integracijsko shemo za rotacijske kolicine, izrazene s kvaternionsko parametrizacijo. Kolicine, izrazene z linearnim popravkom rotacijskega vektorja, moramo zapisati v odvisnosti od linearnega popravka rota- [n] cijskega kvaterniona, kar dosezemo z zvezo med obema (5.77). Najprej v izrazu (8.7) 5?g nadomestimo s gln], rotacijske matrike zamenjamo z rotacijskimi kvaternioni in dobimo ;[n+1] A?gn] o ?[n] = exP (- 5? gn] o g[n], za 5?gn] = 25?gn] o g[n]*. (8.14) Nato enacbo (8.8) z uporabo povezave med linearnima popravkoma rotacijskega vektorja in rotacijskega kvaterniona (4.57) preoblikujemo v kvaternionsko obliko: 2<3jG[II + 1] ◦ 5<3jG[n + 1] = G [n + 1] = At U G?[n] + At 5g[n] =- g[n]* n J A t r,[n] + At2 5gG[n + 1] = 2gG[n+1] 0\ At rG[n] + At ^ a [n] + P a [n+1] 2 P I ^G[n] + P aG[n+1] !_ ^ a [n] + P a [n+1] 2 P J ^G[n] + P ^G [n + 1] Enacba (8.9) ostane nespremenjena. Enacbe Newmarkove integracijske sheme, prirejene za rotacijski kvaternion, so zbrane v preglednici 8.1. Preglednica 8.1: Posplošena implicitna Newmarkova shema za rotacijski kvaternion Table 8.1: Generalized Newmark implicit time-stepping algorithm for rotational quaternion preme kolicine 5Un] = At 4n] + At2 v g , [n+1] g [n+1] g rin] + 5u [n] (2 - P) «gn] + p a [n+1] v + At /-I \ [n] , [n+1] (1 - y) ag]+y ag ] kotne kolicšine [n] 5gG] •~[n] 1 gG] o 2 Ata;GnL+At2 ggn+1]=exp 5ggn]o gjn]* w [n+1] G [n o gg d - p) a [n] + Pa[n+1] i- P aG[n + 1] [n]* X ^[n winL+At (1 - Y) ain[+++]1] + YUG[n| [n] [n] * [n] gG] o exP (gGwo 5gGin [n] 1 n + YaUL Preoblikujmo še eksplicitne iterativne izraze za popravljanje kotnih hitrosti in pospeškov iz iteracije i v iteracijo i + 1 v odvisnosti od osnovnih rotacijskih neznank problema 5gg. Najprej zapišimo kolicino G[n G[„+i] v odvisnosti od trenutnega popravka . Trenutni rotacijski kvaternion qj++1] ima v bazah Bg in B [n+i] enak komponentni zapis g[n+1] = g[n+1] Zapišemo ga lahko kot kompozitum prejšnjih stanj in popravkov na vec nacinov, na primer «:++4++i]=A«:+,c5++i]- (815) Upoštevali smo, da popravke rotacijskim kvaternionom v ustrezno izbranih bazah dodajamo z leve ali z desne (glej poglavje 4.2.2). (Če upoštevamo še, da ima kvaternion gtn] v bazah Bg in BG[n, enak zapis: lahko zapisšemo enakost £[nj = g[n] yg _ yG[n], Agg ◦ C+1] = ggn] ◦ Aq]+1,G[„+1]. (8.16) Iz (8.16) izrazimo č$[n] [ +1] v odvisnosti od rotacijskega kvaterniona in njegovega linearnega po- i+1>Gi+1i ] pravka: A?;+1)Gr++i] = ggn]*o Agg O ^ o *5]+1iGi++i, = exp-1 (ggn]* o Agg o ?!n+1]) (8.17) [n] -1 [n] [n+1] , = 2 exp-1 (ggn]*o Agg o gj^) . (8.18) Inverz eksponentne preslikave kvaternionskega argumenta exp-1 je opredeljen v dodatku B3. Izrazimo še linearno in multiplikativno spremembo med trenutnim rotacijskim kvaternionom gi++1] in rotacijskim kvaternionom gtn] zakljucenega stanja pri casu tn, torej med in Ag[+1, v odvisnosti od trenutne linearne in multiplikativne spremembe, čg in Aq. Najprej v enacbi (8.16) upoštevamo ggn] = [n] [n+1] [n+1] [n] «G[n] in yi+1,g = gi+1)(G[n + 1] ter izrazimo gg ] ?gn]* = g;nG+n]+*1] o Agf! [n+1]. (8.19) 'i 'i Nato iz enacbe (8.17) izrazimo čg[n] [n+1,, uporabimo (8.19) in dobimo i+1'Gi+1 [n] [n] 1 [n] [n+1] čgi+]1'G[n+1] = gg] o exp 1 (gg] o Aggo g- g ] (8.20) [n] 1 [n] [n+1] [n+1] = gg] o exp 1 Ag; G [n+1,o gi,g ] o Aggo gi,g ] = ggn] o exp-1 (Ag^[n+1, o Ag^n+1]) Agl+]1 G[n+1] = ^[n + 1] o AgG[n++1] • (8.21) Po primerjavi z enacbo (8.12) vidimo, daje popravljanje kolicine Ag[n] [n+1] v kvaternionski parame- i+1>Gi+1 M trizaciji rotacij povsem analogno popravljanju kolicine 5$+ G oziroma pripadajoce rotacijske matrike R 5$ [n] i+1,G . Rezultat se prav tako ujema z izpeljanim pravilom iz poglavja 4.2.2 o desnem dodajanju rotacijskih kvaternionov, kadar so ti zapisani v svojih koncnih bazah. S pomocjo izpeljanih izrazov (8.14), (8.18) in (8.20)-(8.21) lahko eksplicitne iterativne izraze za apro-ksimacijo hitrosti in pospeškov (8.11) in (8.13) direktno preuredimo v kvaternionski zapis. Eksplicitni iterativni izrazi za aproksimacijo hitrosti in pospeškov z uporabo kvaternionske parametrizacije rotacije so urejeni v preglednici 8.2. pomiki Preglednica 8.2: Aproksimacije za popravljanje dinamicnih kolicin Table 8.2: Approximate formulae for dynamic quantities zasuki exp o q[n+1] *) [n+1] i+1,g v[n+1] a [n+1] i+1,g [n+1] ] + 5u v[n+1] + 5u vi,g + 0At5u a[n+1] + 1 5u ai,g + ^At2 5u a[n+1] ai+1,g Agg 0 Li,g [n+1] 5-M 5 W^11 [n+1] i+1,«!++11 ' = q[n] 1 [n]* »-~-[n+1] gg ] o Ag„ o g[ ] n n [n+1] i,G [n+1] + Jr.a[n]* o + 0Atg 0 *g 5-[n] 5^[n] 5gi+1;Gin+^1] - 5gi,G!n+1] a [n+1] "i+1,G [n + 1] +1 a [n+1] + 2 a[n] * a, G [n+1] + ^At2 ag 0 5g [n] i+1,G [n+1] +1 — 5q [n] i,G, [n+1] 8.1.1 Zacetni približki za nov časovni korak: prediktor Po zakljuceni iteraciji pri casu tn potrebujemo zacetni priblizek za nov casovni korak tn+1. Rezultati predhodne zakljucene iteracije pri casu tn so po izkušnjah avtorjev clanka [Simo, Vu-Quoc, 1988] neprimerni, da jih direktno uporabimo kot zacetni priblizek pri naslednjem casovnem koraku tn+1, zato predlagajo, da na zacetku vsakega casovnega koraka predpostavimo nicelni prirastek pomikov in zasukov: rgn,+1] = rgn] in r0"+1] = R[n], enacbe za hitrosti in pospeške pa na podlagi teh predpostavk izpeljemo iz enacb za prehajanje med casovnimi koraki, preglednica 8.1. Obrazca za zacetne priblizke hitrosti ostaneta nespremenjena n [n+1] g,0 [n+1] G,0 = vgn] + At (1 — y) agn] + yag,o "(i — y ) a Gn]+ 7s Gn+1] [n+1] [n] G + At (8.22) (8.23) medtem ko sta obrazca za zacšetne priblizške pospesškov drugacšna _ [n+1] ag,o = 1 v" At p Vg —11 (2 —p ag"' S [n+1] = a G,0 _ 1 n [n] At p nG — i (2 - / )a G'1 (8.24) (8.25) Enacbe (8.22)-(8.25) imenujemo prediktor novega casovnega koraka. 8.2 Priprava enacb za numericni racun Pripravljeno imamo diskretizacijo neznank in glavnih enacb (6.82)-(6.83) in naravne robne pogoje (6.86)-(6.89). Skupaj tvorijo sistem 7 • (N + 2) navadnih nelinearnih diferencialnih enacb za prav toliko neznank: -p/k . J 1 • 5/k ocN (7 G, kG) o rk + 5k ocN (YG, «G) ◦ + ocN (Y G, KG) o 5 k ^k \ ^ k | -^k ^ 5 f ^.k \ ^ ^* k | -^k ^ \ ^ / i g + nk - pA ak = 0 k 52 • 5/k oc5m (y G, KG) o 5 *k + 5k ocM (7G, KG) o 5 *k + 5 oc5M (YG, KG) o 5 k .„k \ _ * k i -~-k _ 5 /„,k .„k \ _ * k k k /k JR3 R3 (8.26) (8.27) + mk-5 (qk oCn (7 G, kG O 5* M 5/k - q o J 5 o ak o 5 05 - H„ 50 J 5 o u„ o 5 o 5 = 0 j*k -—k <;k\\ L/2 L/2 h1 • S0 + NL/2 + J ng dx - J pAr ug dx = 0 f-1 . s -r N g T j ngdx - j pAr ugd 00 L/2 L/2 h 2 • P0 + ML/2 + I mg dx -I S (5 O Cn (y g, kg ) o 5*) 5' dx gg 00 L/2 L/2 - 5 o I J (5* o ag o 5)) 05 *dx - H g (5 o (J" u g) o 5*) dx = 0 (8.28) (8.29) L L h3 • SL - NL/2 + J ng dx - J pArug dx = 0 L/2 L/2 L L h4 • pL - ML/2 + mgdx - 5 (q o cn (Yg, kg) ◦ 5*) 5gdx L/2 L/2 (8.30) (8.31) - q o (J" (5* o ag o 5)) 05 *d£ - hhg (5 o (J" u g) o 5*) d£ = 5 L/2 L/2 Ko v diskreten sistem diferencialnih enacb gibanja nosilca vpeljemo aproksimacije dinamicnih kolicin - hitrosti, pospeške, kotne hitrosti in kotne pospeške, ki so zbrane v preglednici 8.2 - iz diferencialnih enacb nastanejo algebrajske enacbe, ki pa so še vedno nelinearne. Zato uporabimo Newtonovo shemo za reševanje nelinearnih enacb, ki smo jo opisali ze v poglavju 5.3.2. Newtonova metoda zahteva linearizacijo dobljenih nelinearnih algebrajskih enacb; kar zadeva linearizacije je bilo veliko dela storjenega ze v poglavju 5.3, kjer smo linearizirali staticne enacbe. Po primerjavi staticnih (5.41)-(5.46) in dinamicnih (8.26)-(8.31) enacb ugotovimo, da moramo pripraviti le še linearizacijo (kotnih) hitrosti in pospeškov in tistih clenov enacb, ki te kolicine vsebujejo. Obicajno dinamicnih clenov v linearizirani obliki ne dodamo tangentni togostni matriki K, temvec jih zdruzimo v novi matriki M. Tako enacba Newtonove iterativne metode (5.50) zavzame obliko K[n] + M[n]) % = -f[n] (8.32) L y [n+1] čy + y[n]. Matriko lineariziranih dinamicnih cienov imenujemo masna matrika, saj v najpreprostejšem primeru tockovnih mas M dejansko vsebuje njihove mase. Matrika K[n] pa se izraza enako kot pri staticni analizi nosilcev. 8.3 Linearizacija hitrosti in pospeškov Povsem preprosta je linearizacija premih hitrosti in pospeškov, izrazenih z enacbami iz preglednice 8.2: Y 1 = Mt^g čag = čug. (8.33) V (8.33) smo poenostavili zapise linearizacije hitrosti in pospeška čvg = D (vg[čug] čag = D (ag)u [čug] . Malo vec dela nam nalagata linearizaciji kotne hitrosti in pospeška. Posebej si oglejmo linearizacijo po osnovnih rotacijskih spremenljivkah. V enacbi (8.20) je ^ izrazena z osnovnimi rotacijskimi spremenljivkami, toda v tej izrazavi nastopa inverz eksponentne preslikave kvaternionskega argumenta. Preslikavo nacelno sicer znamo linearizirati (glej dodatek B), kadar je argument inverzne eksponentne preslikave posamicen rotacijski kvaternion. V izrazu (8.20) pa se v inverzu eksponentne preslikave pojavi produkt treh rotacijskih kvaternionov, kar zahteva posebno skrb. Izraz (8.20) preoblikujemo v obliko I —[n]* r-~-[n] exp gg ] o dg1 J i+1,G. [n + 1] +1 ggn]* o Ag„ o g [n+1] g in po definiciji smernega odvoda (dodatek A, definicija 3) lineariziramo obe strani enacbe: d_ de exp ( ggn] * o čg [n] i+1,G [n+1] +1 (e) £=0 d_ de ?gn]* o Agg (e) o ging+1] £=0 Na vsaki strani enacaja se pojavi odvod eksponentne preslikave kvaternionskega argumenta po skalarju e, saj je Agg (e) = exp ^e čgg o gj™+1]* j. Odvod eksponentne preslikave je izpeljan v dodatku B2; tako dobimo A g[n]* qg o čg [n] i+1,G. [n+1] +1 (e) ggn]* o A (ečgg o g [n+1] i,g g[n]* qg o čg[n] [n+1] (e) de (8.34) [n] gg ] o čg( —[n+1]* o gi,g ] o g J £=0 [n+1]' " i,g d Zaradi vecje preglednosti ločšeno obravnavamo posamezne dele enacbe (8.34). Vstavimo e = 0 in za A f gn]* o 5g[n] i+1,G, [n+1] +1 (eM dobimo A ( ?gn]* o 5g[n] i+1,G, [n+1] +1 (e) )1 = A (ggn]* o ggn] o exp-1 (gJn]* o Ag (0) o gig^)) / J £=0 = A (exp-1 (ggn]* o exp (o 5qg o g]j+1]*) o g[n+1] A ( exp-1 (gM* o 1 o g n+1] i,g [n] A( exP ( Agi G[n + 1] A ggn]* o 5g'nG [n+1]; . [n] Rezultat je smiselen, saj pove, da se v primeru, ko je trenutni popravek 5 gg enak nic, kolicini 5g[ [ +1] g i+1,Gi+1 ] in 5g[n^[n+1] ujemata. V naslednjem delu enacbe (8.34) "< ggn|* o 5gi+1,G;++.] (e) de nastopa iskana linearizacija 5g 4"'* ovt .gH 9 ^ ] £=0 [n] i+1,G A (e5?g o ?!ng+1]* [n + 1] +1 po 5gg. Za naslednji izsek enacbe (8.34), z zapisom £=0 , po definiciji (B.6) matrike A in za e = 0 dobimo A 0 = I. Torej celotnemu izrazu (8.34) pripada zapis a (ggn]* o 5gj[n+1 O (V o d (5g[n] i+1,Gi'++17 q[n+1] [n] [5gg [n+1] [n+1] g i+1,9 [n+1]* g[n+1] g ~ o gi,g , ggn]*o 5g0 o o g iz katerega izrazimo iskano linearizacijo 5g[ [n+1] ;+1>Gi+1 D ( 5g[n]i „[n+1] i+1,Gi'++17 V [n+1] "+1 7 Vi+1,g [5gg] = (ggn]) A-^ggn]* o 5^+1]) (ggn]*) 5gg. (8.35) Sedaj je tudi linearizacija popravkov kotne hitrosti in pospeška cisto preprosta: V 5r [n+1] i+1,G'" + 1] / V[n+1] "+1 7 yi+1,g V 5a [n+1] i+1,Gin+1] / V [n+1] "+1 7 yi+1,g [5g [5g (g[n]*) D( 5g[n] ,[n+1] [n+1] L" ai + 1,9 gJ PAr lV* / v *t+1,Gi,+"+17s[n+1] gJ _ PAt^L \ ^ 1 ~ \ ,'-1-1 ^[n + 1] L "^g r0L (g[n]*) D (W [n] »+1,G[ +1 / qi++1] [5gg M. (8.36) (8.37) Se preprostejši izraz pa dobimo, ce poenostavimo zapisa linearizacije kotne hitrosti in pospeška v trenutni bazi >+1] ^ [^ ] - xn„ -n I x^[n+1] D (5r [n+1] [n+1]^ [5gg] = 5rg D f 5a [n+1] ^) V ;+1,GI^1]y a [n+1] gJ v ;+1,GI^1]y ;+1'G!n+^1]7 a [n+1]r ''+1 7 ai+1,g ;+1,G^]; a [n+1] [5?g] =5aG + 7 ai+1,g ter upoštevamo, da se po vkljucitvi (8.35) v enacbi (8.36)-(8.37) faktorja (jgn]J in (g[n]* j po-krajsšata: čJG = A-1 (qgn]* o čgH [n+1]) (qgn]*) (gg0]) čk, (8.38) čaG = ^A-1 (qgn]* o čJ^) (ggn]*) (gg0]) čkg. (8.39) Koncna izraza za linearizacijo kotne hitrosti in pospeška sta zapisana z uporabo (5.61) v odvisnosti od variacije osnovne spremenljivke čJg = čqg o J0g. V enacbi (6.71) pa poleg j g = Jg* o gg o Jg in a G = Jg* o j g o Jg nastopa tudi gg, ki jo je prav tako potrebno linearizirati. Ker velja gg = Jg o j G o J*, dobimo čiJg = (gg0] o jG o J*) čkg + 0L (qg) (qg*) ČjG + 0L (Jg o jG o Jg0]*) Ačk 8.4 Linearizacija posameznih Členov enaCb Pripravili smo linearizacije posamicnih dinamicnih kolicin, zato lahko nadaljujemo z linearizacijo posa-micnih clenov glavnih enacb f 1, J2 in robnih enacb h1 do h4, ki vsebujejo dinamicne kolicine. Tem clenom dodelimo naslednje oznake dyf = -PAr ag f = -J o (J (J* o a g o J)) oJ* = -g o J a g) oJ* f = -Hg (J o (J (g* o g g o g)) o g* ) = S (g g) (g o (J j g) o J in L/2 L dyh = - / PAr ag dC h = - / PAr ag dC 0 L/2 L/2 ' L , ' dAh = - / g o( j a g) og* de dA h = - J j o( j a^ og* d£ 0L L/2- ( ' ' ( \ dBh = - J J (gg) (J o( J jg) o J* W dBh = - J SS (g^. 0L Antisimetricno matriko Hg oziroma njeno razširitev Hg smo zapisali kot splošno antisimetricno matriko (glej definicijo 3.21 in tocko 1 opombe 10) Hg = S (gg). Linearizacijo vseh dinamicnih kolicin izvedemo direktno. Vektor vseh diskretnih neznank oznacimo z y, njegovo variacijo pa z Sy: y0= r uo i f Su°| 1 50 5g Syg= S5° Sug S5g +1 5gV +1 . 5g . S<+1 Linearizacije dinamicnih clenov oznacimo z operatorjem S, na primer D (\y f J [Syg] = Sgy f. Linens arizacije trikomponentnih vektorskih clenov so preproste Sdyf = -PAr Sag L/2 L Sgyh = - J pAr Safl Sdyh3 = - J pAr Sag 0 L/2 Zaradi podobnosti prikazujemo postopek linearizacije le nekaterih štirikomponentnih clenov. Linearizacija clena ^f je SdAf = -S5g ◦ J a G J o 5* - 5 o J Sa G J o 5* - 5 o (J" a G) ◦ S (5 *) = -0R (5 a G o 5*) S5g - 0L (5) (5 *) 5 Sa G - 0l (5 o5 a g) a S5g. Matriko A smo spoznali ze v poglavju 5.3.3, enacba (5.53), prav tako pa ze poznamo linearizacijo celotnega rotacijskega kvaterniona ter kotne hitrosti in pospeška v smeri S°g. Ker je S linearen operator svojega argumenta, je njegova linearizacija enaka zacetnemu operatorju, DS^ [Saj = S (Sa), kar je po-kazano v primeru 1, glej dodatek A. Linearizacija dBf je tako SdB f = -S (Sig) o (j i G J o 5* j - S (Wg ) (^jfg o (J5 i5 G J o 5* - s (ig) (5 o (J1 Sig) o 5*j - s (ig) (5 o j £g) o S5g* = s (5 o (J5 £ g) o 5*) Sig - s (ig) (j £ g o 5*) S5„ - S (ig) 0L (J) (5 *) 5 Si G - S (ig) 0L 5 o 5 i G A S5 Pri preureditvi prvega clena iz kvaternionskega v matricni zapis smo uporabili lastnost (5.69) antisime-tricnih matrik: 5 (a) b = -5 (b J 5, kar velja za poljubna cista kvaterniona a, b. Linearizacija preostalih štirikomponentnih dinamicnih clenov ima obliko: L/2 5dAh J5aG o J* 5g + 0L (J) (J*) J5ag J o J 5ag A 5Jg d{ L/2 52B h S (J o ( J n g) o J* 5ng — S (ng) ( j n g o J*) 5J„ —S (ng) 0l (J) (J*) J 5LjG — S (ng) 0L j o J nG A 5J„ d{ 5dy h L = L =- L/2 +0l( J L c = S J L/2 ( J5aG o J*) 5Jg + 0L (J) (J*) J5ag 8% h —S (ng) 0l (J) (J*) J 5LjG — S (ng) 0L [J o J nG j A 5JgJ Z upoštevanjem interpolacije osnovnih neznanih funkcij (6.81) oziroma njihovih linearizacij N+1 N+1 8ug (x, t) = E 5u (t) L^ (x) 5J (x, t) = E 5 J (t) Pi (x) i=0 i=0 lahko variacije osnovnih kolicin v izrazih 5dyh1 — 5dyh3 in v 5dyf 2a — 5dyf 4b hitro izpostavimo iz integrala. 8.5 Numerična integracija po kraju Tako kot pri enacbah za staticno analizo nosilcev clene v integralski obliki integriramo z uporabo Gau-ssovih kvadraturnih formul (5.76). Izbira stopnja integracije je poljubna. Navajamo jo pri posameznih numericnih primerih. V primeru linearnih elementov brez internih kolokacijskih tock izberemo dvotoc-kovno integracijo. 8.6 Lokalna napaka in prilagajanje časovnega koraka metode V numericno implementacijo vgradimo preprosto, a ucinkovito dvokoracno preverjanje lokalne napake. Rezultate y[n+1'At] v casu tn+1, ki so bili izracunani s korakom At = tn+1 — tn iz rezultatov y[n] pri casu tn, primerjamo z rezultati, izracunanimi z dvema zaporednima casovnima korakoma polovicne dolzine At 2 y[n] —^ y[n+At/2] —^ y[n+1]. (8.40) g * i- o o o g g g Primerjavo izvedemo s socasnim preverjanjem absolutne in relativne lokalne napake s pogojem min y [n+1, At] _ y [n+1] |y[n+1,At] — y[n+1] II ] |y y ^ 0 in nasprotna smeri tangente za s < 0. Enacba (9.11) bi imela v primeru, daje s naravni parameter v poljubni legi (oznacimo ga s s krajevnega vektorja po naravnem parametru velja preprostejšo obliko: Vd = s et = Vd et, saj za odvod = 1 in s = Vd, glej na primer [Saje, 1994, str. 5]. Tudi splošnejšo enacbo za hitrost (9.11), kjer s ni naravni parameter v splošni legi, lahko ob izbiri oznake Vd = s jjr'j (9.12) zapisšemo enako preprosto: Vd = Vd et, (9.13) kjer je | Vd | velikost hitrosti, sign(Vd) smer gibanja delca glede na tangento in sicer: sign(V) > 0, ce se delec giblje v smeri tangente in sign(V) < 0, ce se delec giblje v nasprotni smeri. Odvod vektorja hitrosti po casu je pospešek: _d = _d = sllf'11 _ + s (lir'ID' _ + S ||r'y _ s ||r || + s- r' • r' |r = [s||V'|| + sS2 (_ • et)] et + S2 ||_'|| et'. Vektor et' smo z enacbo (9.7) ze izrazili z baznimi vektorji Frenetovega ogrodja. Ko upoštevamo še enacbo (9.10), lahko zapišemo _ . 11_y u d__ ds et S "r " ds dt _d = [s'|_+ s2 (_ • _t)] _t + s = [s'||_'|| + s2 (_' • et)] _t + s 2 | _r '| l_' x (r" x _') 2 | _ | 2 1 _ s y_ || — en. Pd (9.14) Enacba pove, da sta pri gibanju delca po krivulji nenicelni samo dve komponenti pospeška; prva je v tangentni smeri in druga v smeri glavne normale; komponenta pospeška v smeri binormale je nic. V kolikor odvajamo preprostejšo obliko enacbe hitrosti (9.13), je formalna oblika enacbe za pospešek preprostejsša _d = v d _t + Vd _'s _ ii_' x (r" x _')ii _ = Vd et + Vd s--—o-en 1 v det + vd—en Pd (9.15) in se ujema s formalnim zapisom za primer naravne parametrizacije osi, glej [Andzelic, Stojanovic, 1965, str. 28, en. 2] ali [Saje, 1994, str. 9, en. 23]. Pri zapisu pospeška v obliki (9.15) je jasnejši pomen posameznih clenov. Na pospešek v tangentni smeri vpliva spreminjanje velikosti hitrosti, na pospešek v glavni normalni smeri pa velikost hitrosti skupaj z glavno ukrivljenostjo osi. n 9.3 Gibalne enaCbe delca Gibalna enacba (6.1) doloca, daje vsota vseh sil enaka odvodu gibalne kolicine. Gibalna kolicina delca je enaka produktu mase in hitrosti, njen odvod pa je ob predpostavki o nespremenljivi masi delca enak produktu mase in pospeška. Rezultanto sil Rd, s katero opišemo delovanje nosilca v tocki T na delec, zapišemo kot linearno kombinacijo baznih vektorjev Bd (slika 9.1 (b)) Rd = Td_t + Nd_n + Bd_b = Td + Nd + Bd. (9.16) Silo Td imenujemo sila trenja, Nd in Bd pa normalna in binormalna sila podlage. Za reakcijske sile, zapisane po smereh naravnega triedra, po Coulumbovem zakonu [Andzelic, Stojanovic, 1965, str. 281] velja, daje velikost sile trenja sorazmerna velikosti rezultante normalnih sil Nd + Bd, in daje njena smer nasprotna smeri hitrosti delca Td = -ju Nd + Bd vd (9.17) (a) Sile, ki delujejo na nosilec: (b) Sile, ki delujejo na delec: x=0--- \Td T ms2 d B N -i ^ ms x=L mge„ Slika 9.1: (a) Sistem sil na nosilcu; (b) sistem sil na delcu Figure 9.1: (a) Forces acting on a beam; (b) forces on a particle Sorazmernostni faktor p imenujemo koeficient trenja. Koeficient trenja je odvisen od podlage, hitrosti podrsavanja in velikosti sticne površine teles ter se doloca eksperimentalno. Tukaj predpostavimo, daje koeficient trenja konstanten. Enacbo (9.17) lahko zapišemo tudi drugace et. Td = -sign (s) p Nd + Bd Poleg reakcijskih sil upoštevamo še gravitacijsko silo Fd, ki deluje na delec Fd = mgeg, kjer g oznacuje velikost teznostnega pospeška, eg pa enotski vektor smeri delovanja sile teznosti. Tudi Fd zapišemo z vsoto sil po smereh naravnega triedra Fd = F + F„ + F, = Ftet + F„e„ + F6e6 za Ft = Fd ■ et Fn = Fd ■ Fn F, = Fd ■ F,. (9.18) Ostale zunanje vplive zanemarimo. Gibalne enacbe delca dolocajo, daje vsota zunanjih sil (9.16) in (9.18) enaka produktu mase in pospeška Rd + Fd = mad, kar lahko ob upoštevanju (9.15) razbijemo na komponentne enacšbe /dl : Td + Ft = mv)d (9.19) fd2 2 1 : Nd + Fra = mv2 — (9.20) Pd /d3 : Bd + F, = 0. (9.21) Nove neznanke so poleg opravljene poti s , glede na nedeformirano os nosilca, tudi komponente reak- cijske sile Td, Nd in Bd. Dodatno, cetrto enacbo, predstavlja Coulumbov zakon (9.17). Neznanko Bd znamo takoj izracšunati in sicer Bd = -F,, zato enacbe (9.21) ne dodajamo sistemu enacb. Enacbi (9.19)-(9.20) pa lahko z uporabo enacbe (9.17) zdruzimo v skalarno diferencialno enacbo za neznano funkcijo s (t) fd : -sign (s) u 2 1 _ _ _ mvd—en - Fn - Fb Pd + Ft - mVd = 0. Ce upoštevamo, da lahko komponentam sile teznosti izpostavimo velikost mase Ft = mg (_g ■ _t) lahko maso pokrajšamo iz enacbe dobimo Fn = mg (_g ■ _n) Fb = mg (_g ■ _b), poleg tega oznako Vd nadomestimo z izrazom iz enacbe (9.12) in fd : -sign (s) u . 2 ii_' x (r" x _') "_n - g (_g ' _n) _n - g (_g ' _b) _b (9.22) + g (_g ■ _t) - s^'^ - s2 (r" ■ _t) = o. Teznostni pospešek g = 9.81 m/s je konstanta in za vektor _g na površini Zemlje predpostavimo, daje konstanten vektor. Ostali skalarji in vektorji so bili ze predhodno izrazeni z osnovnimi spremenljivkami. Enacbo (9.22) dodamo sistemu enacb prostorskega nosilca (8.26)-(8.31). 9.4 Upoštevanje vpliva gibajoCega se delca na nosilec V prejšnjem poglavju smo pri zapisu gibalnih enacb upoštevali silo Rd, s katero nosilec v tocki stika deluje na delec. Zato moramo tudi na nosilcu v tocki stika T upoštevati obratno silo -Rd, s katero delec deluje na nosilec. Vsaka tocka tezišcne osi nosilca zadošca lokalni ravnotezni enacbi (6.70). Lokalno ravnotezje sil velja tudi za tocko T s krajevnim vektorjem _(s (t) ,t), v kateri se trenutno nahaja delec. Robni enacbi (6.86) in (6.88) smo izpeljali iz enacb (6.45) in (6.47) z integracijo enacbe lokalnega ravnotezja (6.70) na obmocjih [0, L) in [L, L]. Notranja sila N (x, t) ima pri _(s (t) , t) tocko nezveznosti (skok), zato locimo levo in desno limito notranje rezultantne sile pri T: (s, t) = limN (x - M) £—>0 N_d (s, t) = lim N (x + M). £0 V skladu s sliko 9.1(a) je N_D - N_L = - Nd - Bd. Sile na desni izracunamo iz komponentnih gibalnih enacb (9.19)-(9.21): Td = Nd = ms ||_'|| + ms2 (r" • _t) - ( _d • _t • 2 u_/112 _ _ ms ||_ || Kd - Fd • e, et Bd = - ( Fd • e^ eb. (9.23) (9.24) (9.25) (9.26) n Ker so enacbe nosilca zapisane v referencni bazi Bg, tudi reakcijske sile med delcem in nosilcem Td, N in Bd zapišemo v referencni bazi Bg: Nd,g, Td,g in Bd,g. Locimo primera, ko se delec nahaja v prvi ali v drugi polovici elementa: 1. Naj velja s G [0, D. Potem integral enacbe lokalnega ravnotezja sil (6.86) po tem obmocju izracunamo kot: s L/2 L/2 J Ngdx + J Ngdx + J (ng - pArrg) dx = 0. 0 s 0 Z integracijo dobimo L/2 NL (s, t) - Ng (0, t) + N g (L - nD (s,t)+/ (ng - pAr rg) dx = 0. 0 Po urejanju sledi L/2 nL (s, t) - nD (s, t) + N g (L, - Ng (0, t) + J (ng - pAr rg) dx = 0. Nd,g +Td,g +Bd,g 0 Razliko leve in desne vrednosti notranje rezultantne sile pri x = s po enacbi (9.23) nadomestimo s tockovnimi silami, ki predstavljajo vpliv gibajocega se delca na nosilec. Izrazimo sedaj rezultantno notranjo silo pri x = 0 L/2 N0 = Td,g + Nd,g + Bd,g + Ng + (ng - pAr rg) dx in jo uporabimo v robnem pogoju (6.45) L/2 'g + T d,g + N d,g + Bd,g + N g + (ng - PAr 'g) S + T d,g + N d,g + Bd,g + N g + J (ng - pAr 'g) dx = 0. (9.27) Kadar se potujoci delec nahaja na obmocju [0, D , moramo torej robni pogoj (6.86) nadomestiti s pogojem (9.27). 2. Podobno postopamo v primeru s G [L, L]. Enacbo lokalnega ravnotezja sil (6.88) razbijemo na dva dela in dobimo L N L (s, t) - N D (s, t) + Ng (L, t) - N g (L ,t) + J (ng - pAr 'g) dx = 0. Td,9+Ard,9+Bd,9 L/2 Izrazimo rezultantno silo pri x = L L i r N L = -T d,g - N d,g - Bd,g + Ng2 - J (ng - pAr 'g) dx L/2 in jo uporabimo v robnem pogoju (6.47) SL + + + Bd>g _ Ng2 + J (n _ pAr rg) dx = 0 L/2 (9.28) Robni pogoj (6.88) nadomestimo s pogojem (9.28), kadar je potujoci delec na obmocju [L, L] 9.5 Numerično reševanje Vpliv gibajocega se delca dodamo implementaciji dinamicne analize linijskih konstrukcij iz poglavja 9, ki uporablja Newmarkovo integracijsko shemo. Zato tudi dodatno enacšbo resšujemo z uporabo Newmar-kove sheme. V skladu z enacbami iz preglednice (8.1) torej velja čsln] = Ats[n] + At2 s[n+1] = s[n] + 0s[n] 1 ^ _ 3 ) s[n] + /3 s • [n+1] s[n+1] = s[n] + At (1 _ y) s[n] + ys • [n+1] (9.29) (9.30) (9.31) za 3 G [0, 2] in y G [0,1]. Po vpeljavi teh zvez enacba (9.22) postane algebrajska nelinearna enacba in jo skupaj z enacbami nosilca rešujemo iterativno z Newtonovo iteracijo. Iterativna aproksimacija hitrosti in pospeška delca v smeri tangente je po enacbah iz preglednice 8.2 oblike s [n+1 = s [n+1] + si+1 = si + 3 At s[n+1] = >+1] + si+1 = s + 3 At2 kjer je 5s trenutni popravek Newtonove iteracije. Nov priblizek opravljene poti v novi iteraciji je sin+1]= sjn+1]+ (9.32) (9.33) (9.34) Prediktor pospeška v novem casovnem koraku dolocimo enako kot v poglavju 8.1.1, in sicer tako, da predpostavimo nicelno spremembo poti. Tako iz enacbe (9.29) sledi obrazec prediktorja za spreminjanja drugega odvoda poti delca s[n+1] =__L s[n] _ i _ 3\ s[n] s = /Ats A 2 3)s , enacba (9.31) pa sluzi za prediktor prvega odvoda poti delca. Prediktorji prvega casovnega koraka pri t = Ato imajo tako obliko 1 [1] =__s[o] _ ± L _ 11 3At 3 2 3s [o] s[1] = s[0] + At [(1 _ y) s[0] + 7s'[1] s[1] = s[o] = 0. L 9.6 Linearizacija Zaradi iterativega reševanja enacb z Newtonovo metodo moramo enacbo (9.22) linearizirati. Prav tako moramo linearizirati nove clene (9.24)-(9.26) v robnih enacbah (9.27)-(9.28). Linearizacijo moramo izraziti v odvisnosti od osnovnih diskretnih neznank, ki smo jim dodali še parameter poti s. Zaradi preglednosti najprej linearizirajmo posamezne kolicine pri argumentu (s (t) , t). 9.6.1 Linearizacija krajevnega vektorja in njegovih odvodov Krajevni vektor rg = rg0] + ug do tocke T je funkcija spremenljivk ug in s, kar poudarimo z zapisom rg (ug ,s) = rg0] (s) + ug (ug, s) . Posebej si oglejmo linearizacijo krajevnega vektorja do zacetne lege rg0] (x (s)): {fM(x (,)) = g«,. bodisi x lN m(t) skupnih dolzšin Zveza med x in s je bodisi oblike x = s, kadar smo na prvem ali edinem elementu, ki ga delec pretecše, N m(t) s — YI L™, kadar je delec ze pretekel elemente el™, el™, ...,elm i=1 N m(t) L™. V obeh primerih je df = 1 in linearizacijo krajevnega vektorja rg ] lahko zapišemo krajše i=1 [0] 5rg0] = ^5s = rg0]/ (s) te. Podobno postopamo tudi pri linearizaciji pomika N+1 N+1 čun dL, dx £L.K + £ ug^^ i=0 i=0 N+1 N+1 E KL, + E ugL,5s. i=0 i=0 Linearizacija krajevnega vektorja tocške T je tako enaka N +1 N+1 i=0 5rg = rg0]/čs + ^ 5ugL, + £ ug i=0 N +1 = rg*s + Li^ug. (9.35) i=0 Linearizacija njegovega prvega odvoda je čr' = čr' + 5 N+1 Eug L . i=0 N +1 rg'čs + E L,5ug i=0 drugega odvoda pa rw+i 5rg' = 5rg' + 5 £ug L i=0 + 1 rg''5s + £ Li'5u; i=0 9.6.2 Linearizacija baznih vektorjev Najprej se lotimo linearizacije tangente (9.1). Norma vektorja je kvadratni koren skalarnega produkta vektorja s seboj: I 'g I = ('g "g)2 in izracšunamo linearizacijo skalarnega faktorja I-'II/ ii'' j-3 v g 'g rgll/ ||'g| (5-g "g) . (9.36) Variacijo et g izrazimo z 5rg: 5et,g = -1 ('g "g)- 2 K "g +-g • 5-g) 'g + ^ 5-g g (5-g "g) 'g + jjrn 5-g. Il'g N ll'gll V levem clenu 5rg nastopa v obliki, ki ni primerna za uporabo Newtonove iterativne metode. Zato produkt vektorjev (5rg-'g) 'g zapišemo po pravilu (6.12): (5-g-g) 'g = 'g X ('g X 5'g) + jj-g II2 5'g. Rezultat uporabimo pri izrazavi 5et g in po preurejanju dobimo: 1 II'' II2 1 5et,g = -^^'g X ('g X 5rg) - ^^5rg + m-Tm 5rg llrg - jlrg - "'g'! S (et,g) S (et,g) 5-' I '' I g 'g Vektorski produkt smo nadomestili z operatorskim zapisom (5.7) antisimetricnega operatorja S. Pri linearizaciji glavne normale izhajamo iz zapisa (9.5), na katerem uporabimo enacbo (9.36), in dobimo: 5 = 5rg X (rg' X rg) + rg X (5rg' X rg) + rg X (rg' X 5rg) 5en,g = ^ rm r" x r' (-g X (-g' X -g)) l-g X lrg' X -g 7T-3 (5 (-g X (-g' X rg)) ■ (rg X (rg' X rg))) . 1 1 5 Prvi clen, oznacimo ga z A, z uporabo antisimetricnega operatorja in menjave faktorjev preuredimo A = [S (rg) S (rg') - S (rg' x rg)] črg - S (rg) S (rg) črg' ||rg x (rg' x rg)| , medtem ko drugi clen zahteva vec dela, zato vpeljemo oznake za posamicne dele r' x (r" h — v g x ^'g °g = - r' x (r',; x rg)) llrg x lrg' x rg)ll Bi = « x (rg' x rg)) ■ (rg x (rg' x rg)) B2 = (rg x«' x rg)) -(rg x (rg' x rg)) B3 = (rg x (rg' x črg)) ■ (rg x (rg' x rg)) . V B1 po pravilu mešanega produkta ciklicno zamenjamo faktorje Bi = [(rg' x rg) x (rg x (rg' x rg))] ■ črg. (9.37) Taka oblika zapisa ni primerna za numericno implementacijo (enacba (8.32)), ker variacija osnovne (vektorske) kolicine nastopa v skalarnem produktu. Zato je potrebno zapis splošne oblike bgB1 = bg [cg • ], ki se pojavi v linearizaciji en g po upoštevanju enacbe (9.37), nadomestiti z vsoto produktov bg (cg • čdg) = Cg x [bg x ] + [bg-Cg] čdg (9.38) in dobimo bgBi = S ([rg' x rg] x [rg x (rg' x rg)]) S(bg) črg + (bg • [(rg' x rg) x (rg x (rg' x rg))]) črg. Podobno preuredimo še B2 in B3 z dvojno zaporedno ciklicno zamenjavo B2 = [(rg x (rg' x rg)) x rg) ■ [< x rg] = [rg x ((rg x (rg'x rg)) x rg)] ■ < (9.39) B3 = [(rg x (rg' x rg)) x rg] ■ [rg' x črg] = [((rg x (rg' x rg)) x rg) x rg'] ■ črg, (9.40) s cemer zavzamata enako splošno obliko kot B1, zato lahko ponovno uporabimo (9.38) in dobimo bgB2 = S (rg x [(rg x (rg' x rg)) x rg]) S (bg) < + (bg • [rg x ((rg x(rg' x rg)) x rg)]) *rg' bgB3 = S ([(rg x (rg' x rg)) x rg] x rg') S(bg) črg + (bg • [((rg x (rg' x rg)) x rg) x rg']) črg. Zadnji bazni vektor ima preprosto linearizacijo g = ^et,g x en,g + et,g x ^en,g = -S (en,g) 5et,g + S (et,g) čen,g. 3 Opomba 18 V linearizaciji krajevnega vektorja (9.35) nastopa odvod interpolacijske funkcije po kraju. V linearizaciji baznih vektorjev krivočrtne baze nastopajo drugi odvodi linearne spremembe krajevnih vektorjev po času. Da zagotovimo vplive med delcem in nosilcem moramo zagotoviti, da so tretji odvodi interpolacijskih funkcij neničelni; za polinomsko interpolacijo je zato potrebno izbrati vsaj polinome stopnje tri kar ustreza izbiri dveh notranjih kolokacijskih točk. 9.6.3 Linearizacija ukrivljenosti 22 Ker v enacbah nastopa krivinski polmer v obliki M rg 11 /pd, pripravimo še linearizacijo izraza M rg 11 Kd, ki ga zapišemo iir/ m2 k = ||rgx (rg/ x rg llrgll Kd = FTm2 • r' , rg Linearizacijo števca ulomka, ki doloca ukrivljenost, se skriva ze v linearizaciji normale, saj je rgx «x rg) [(rg x (rg' x rg)) ■ (rg x (rg' x <))]-2 [S (rg x (rg' x rg))] ■ (rg x (rg' x rg)) 1 [S (rg x (rg' x rg))] ■ (rg x (rg' x rg)) , r' r'' r' kar bomo upoštevali direktno v enacbi, ki vsebuje faktor Kd. Posebej pripravimo linearizacijo imenovalca: S ^^ =- fTF rg • Srg • v|rgM) KM Linearizacija celega izraza ima torej obliko S (M2 Kd) = M B + b2 + B3 m2 -2 ||rg *(rg'x rg)M rg ■ Srg v ' |rgx wx rg) IIKII rgM za Bi, B2 in B3 definirane z enacbami (9.37), (9.39) in (9.40). (9.41) 9.6.4 Linearizacija Kontaktnih sil. Linearizacija reakcijske sile v smeri tangente je oblike ST, d = Sr' r' mSs M|r'M + ms g g - (Fd • Set,g) + mS {S2) (rg'-et,g) - ms2 (Srg'-e^g) - mS2 (rg' • Set,g)] et,g + [ms M|r'M + ms2 (rg-e^g) - (Fd-et,g)] Set,g. Ponovno nekajkrat uporabimo enacbo (9.38) in po urejanju dobimo ST d = m M|r'M et,g Ss + 2ms (rg'-et,g) et,g Ss + ms [S (et,g) S (et,g) + 1] Srg + ms2 [S (et,g) S (et,g) + 1] Srg - [S (Fd) S (et,g) + (Fd-et,g)] Set,g + ms2 [S (rg') S (et,g) + (et,g-rg')] Set,g + [ms M|r'M + ms2 (rg'-et,g) - (Fd-et,g)] Set,g. Pri linearizaciji reakcijske sile v smeri glavne normale za linearizacijo skalarnega faktorja ||r'y Kd upoštevamo ze izpeljano formulo (9.41). Tako je čNd,9 = 2msčs lir'II2 Kde„„ + ms • 2 B1 + B2 + B3 r.// rg rg 2 en,g — ms ^ vg"' g/ MM' g 2-1 lrg X (rg' X rgM" (rg • Srg) e n,g (Fd • čen,g ) en,g + r • 2 M^ M 2 /p, \ ms "r II Kd - (Fd^e„,g) če, n,g • Vec izrazov moramo preurediti v obliko, primerno za Newtonovo iteracijsko shemo (8.32). To storimo z uporabo enacbe (9.38); po preurejanju sledi čNd,g = 2ms "rg "2 Kde„,g čs + ms2 (bgBi + bgB2 + bg£3) 0 "rg X (rg' X rg)" - 2ms2 g " g (S (et,g) S (e„,g)) črg "rg" - (S (Fd) S (era,g) +2Fd •en,g -ms2 "r'" 2KdM če«,g, pri cemer so ustrezni zapisi produktov bgBj ze pripravljeni v poglavju 9.6.2, le daje vektor bn definiran drugacše in sicer en bg -ra,g r' r'' r' Lotimo se še linearizacije zadnje reakcijske sile Bd (9.26): čBd = - (Fd • če6,g) e6,g - (Fd^efc,g) čeS,s = -S (Fd) S (eft,g) če6,g - 2 [Fd-e6,g] če6,g. 9.6.5 Linearizacija gibalne enacbe delca Gibalna enacba delca (9.22) sestoji iz vec clenov; prvi med njimi zaradi kompleksnosti zahteva loceno obravnavo, zato vpeljemo oznako f d = -siSn (s) ^ = -sign (s) ^ s "r " Kd g (eg,g • en,g) en,g g (eg,g • eb,g) eb N d,g + B d,. m Ker je funkcija sign konstantna preslikava na obmocjih (-to, 0) in (0, to), v tocki 0 pa ni odvedljiva, velja č (sign (s)) = 0, za s = 0. 2 r Upoštevamo še izpeljano linearizacijo produkta ||r'|| Kd podano z enacbo (9.41) in dobimo 1 _ sign (s) p Nd,g + Bd,g 5Nd,g + d,g f d = — N d,9 m m sign(s) P f • 2 nr/ 112 N d,g { s2 ||rg II2 Kd — g (eg,g ■ era,g) e„,g — g (eg,g ■ e6,g) e6,g| 11 rg||2 Kde„,g čs + s2 (bgB1 + bgB2 + bgBa) -2s 2 Ilrg x (r£ x rg. LS (et,g) S (e„,g) 5rg — (gS (eg,g) S ( ) + 2ge "g,g °n,g _ • 2 || /|| 2 s r Kd) če, n,g — (gS (eg,g) S (eb,g) + geg,g ' eb,g)5eb,g}. (a) (b1) (b2) (b3) (b4) (b5) N Loceno si oglejmo posamezne clene brez skalarnega faktorja — (sign (s) p) / Nm^- . V clenu C1, ki pripada produktu izrazov (a) in (b1), upoštevamo, daje skalarni produkt baznih vektorjev eb,g in en,g enak nic in dobimo C = N d,g m ^2s 11rg||2 Kde„,g = 2smd |rg ||2 (Nd,g ' en,g) ^. V clenu C1, kije produkt izrazov (a) in prvega clena izraza (b2), upoštevamo, daje produkt eb,g ■ bg, za bg = — (rg x (rg/ x rg)) r/ r// r/ enak nic in dobimo C2 = % ' (s2bg B1) =s m •2 N? ■ (S (c1,g) S (bg) 5rg + (bg ■ c1,g) 5rg) za c1g = [rg/ x rg] x [rg x (rg/ x rg)]. Dobljen izraz preoblikujemo tako, da za mešani produkt vektorjev dvakrat zaporedoma uporabimo ciklicno permutacijo faktorjev: 1 _ • 2 Nd,g C21 = s m (S (c2,g) S (bg) 5rg + (bg ■ c2,g) 5rg) = s2- (—d,g ■ (c2,g x (bg x črg)) + (bg ■ c1,g) (—d,g ■ *rg)) m 1 = s2m ((—d,g x c2,J ■ (bg x črg) + (bg ■ (—d,g ■ *rg)) = s2m (((—d,g x c1,g) x bg) ■ črg + (bg ■ c2,g) (—d,g ■ *rg)) = s2m (((—d,g x c2,g) x bg) + (bg ■ c1,g) —d,g) ■ črg. m m a Enako postopamo še pri preostalih dveh clenih C| in C|, kijih dobimo s produktom izrazov (a) in (b2); tako je C2 = ^ ■ (s2bgB2) 2m = s2m (((Nd,g x c2,g) x bg) + (bg.c2,g) Nd,g) ■ *rg' za c2,g = rg x ((rg x (rg' x rg)) x rg) in C23 = ^ ■ (s2bgB3) 2m = s2m (((Nd,g x c2,g) x bg) + (bg-c2,g) Nd,g) ■ črg m 4g = ((rg x (rg' x rg)) x rg) x rg'. Tudi clen C3 = Nd,g + Bd,g ■ ( -2s2 |rg x„(rg'„x rg;i' S (et,g) S (en,g) črg m \ ||r' | rg ki je enak produktu izrazov (a) in (b3) preoblikujemo enako C3 = -s22 |rg x |rg|x rg) ^ (Nd,g + Bd,g) ■ (et,g x (en,g x črg)) m 11 rg 11 „2 ||rg x (rg' x rg )|| . „ = -s2 " g | g113 ((Nd,g + Bd,g) x et,g) ■ (en,g x črg) m 11 rg11 s • 2 2 Kx (rg'x rg) m 11 rg | (((Nd,g + Bd,g) x et,g) x en,g) ■ črg. Za poenostavitev zapisa upoštevamo še, daje produkt ((Nd,g x et g) x en g) vektor v smeri et g velikosti iN d,g i, medtem ko je ((Bd,g x etg) x eng) enak nic in dobimo „2 ||rg x (rg'x rg )|| . „ C3 = -s2 11 g | g| . 3 g;" iNd,gi (et,g ■ črg) . m | r 'g | Naslednji clen C4 = -Nd,gm Bdg ' (gS (eg,g) S (en,g) + 2geg,g ■ en,g - s2 ^r'^ 2Kd) čen,g, je enak produktu izrazov (a) in (b4). Mešani produkt ponovno preoblikujemo tako, da ciklicno permuti-ramo faktorje: g C4 = - m ((((Nd,g + Bd,g) x eg,g) x en,g)) ■ čen,g g - m2 (eg,g'en,g) (Nd,g + Bd,g) ■ ^en,g + mm s2 yr'y 2Kd (Nd,g + b d,g) ■ ^en,g. Podobno postopamo še pri C5, ki ga dobimo kot produkt izrazov (a) in (b5); dobimo g C5 = -m ((((Nd,g + Bd,g) X ) x e6;g)) ■ če6,g g - m (eg,g'eb,g) (Nd,g + Bd,g) ■ . Linearizacija cele enacbe (9.22) je s f d = f d + g (eg,g • ) - ||rgI Ss - (rg • ^rg) 11 rg II - 2s (rg'-et,g) čs - s2 (et,fl ■ <) - s2 (rg' ■ čet,g) . 9.7 Določanje lege delca v konstrukciji Vzemimo, daje potujoci delec do casa t pretekel pot cez vec elementov em konstrukcije dolzin Lm, za i = 1, 2, ...Nm (t) ter se nahaja nekje na elementu emm(t)+r Potem je zveza med parametrom dolzine poti s in krajevnim parametrom x elementa emm(t)+1 enaka N m(t) s (t)= £ Lm+x (t), i=1 kjer je Nm (t) število ze pretecenih elementov konstrukcije. Glede na vpetost spremenljivke s v robne enacbe prostorskega nosilca, moramo posebej dolociti lego delca glede na sredinsko tocko x = L elementa emm(t)+1, kar dosezemo s pogojem N m(t) Lm s (t) - £ Lm 0 za vsak t > 0, zato pri izracunu zacetnega pospeška izberemo sign (s) = 1. Tako dobimo zacetni vrednosti drugega odvoda poti delca po casu za oba testna primera: s = 6 za gladko podlago in s = 2 za hrapavo podlago. V primeru gladke podlage se numericni rezultati za s (t) izvrstno ujemajo z ana-liticnimi; napaka je enakega reda kot pomiki zaradi povesa nosilca (w 10-10), ki predstavljajo zacetno (neodstranljivo) napako glede na analiticno rešitev. V primeru hrapave podage je razhajanja med nu-mericnimi in analiticnimi rešitvami opaznejše, vendar so relativne napake še vedno majhne (w 10-6) in se s casom ne vecajo. 9.8.2 Delec na prostoležeCem nosilcu Z odzivom prostolezecega nosilca pri prehodu delca so se ukvarjali ze številni avtorji, med njimi Mofid in Akin (1996), Xu in sodelavci (1997), Yavari in sodelavci (2002) ter Lee (1996a), ki so obravnavali prehod delca s konstantno hitrostjo, medtem ko sta Lee (1996b) in Wang (1998) analizirala odziv nosilca pri pospesšenem prehodu delca. 93 m 9i X Slika 9.3: Zacetna lega za primer delca na prostolezecem nosilcu Figure 9.3: Simply supported beam: the initial configuration Posebno poenostavitev predstavljenega problema, kjer ne upoštevamo dinamicnih efektov mase delca temvecš zgolj obtezšbo nosilca s silo, ki jo povzrocši tezša delca, predstavlja prehod tocškovne sile preko nosilca (prehod sile). Za primer enakomernega gibanja sile obstaja tudi analiticna rešitev po linearni teoriji. Analiticne pomike v precni smeri (v smeri obtezbe) izracunamo po obrazcu (glej [Muršic, 1972]) 2 Fm ^ 1 / ^ Qra \ nnx y (x, t) = —— > —2-—2 sin Q„t--sin w„t sin —— (9.42) APL n=1 - Qn\ J L — L n2n2 /EJi kjer v[0] oznacuje konstantno hitrost prehoda sile in F njeno velikost. Namesto tocne rešitve y (x, t) za primerjavo pomikov vzamemo priblizek, pri katerem uporabimo prvih deset clenov vrste (9.42). V pricujocem poglavju izpeljano teorijo vpliva delca na nosilec je preprosto poenostaviti za primer prehoda sile s konstantno hitrostjo. V sistem enacb gibalne enacbe delca (9.22) ne dodamo, upoštevamo pa robni enacbi (9.27) oziroma (9.28); v njih prilagodimo sile, ki dolocajo vpliv delca na konstrukcijo tako, da produkt mase in teznostnega pospeška nadomestimo s potujoco silo F, preostale clene, ki vsebujejo maso, pa izlocimo in dobimo Fd = - (Fd • et) et = - (F • et) F NVd = - (Fd • e„) e„ = - (F • e„) e„ Bd = - (Fd • e6) eb = - (F • F,) e6. Trenutno pozicijo s potujoce sile dolocimo z obrazcem za enakomerno gibanje s = v t in v = v[0]. Primer 1. Vzemimo prostolezece podprt nosilec dolzine L = 4.352, prek katerega s konstantno hitrostjo v = 27.49 potuje delec z maso m = 21.83, slika 9.3. Ostale karakteristike nosilca so E = 2.02 ■ 1011 G = 7.7 • 1010 p = 1.5267 ■ 104 Ai = 1.31 ■ 10-3 A2 = A3 = 9.16 ■ 10-4 Ji = 1.142 ■ 10-6 J2 = J3 = 5.71 ■ 10-7. Podatki so privzeti po literaturi [Yavari et al., 2002]. Gravitacija deluje v nasprotni smeri tretjega baznega vektorja: eg = -g3, teznostni pospešek pa je velikosti g = 9.81. Numericnih podatkov racuna Yavari in sodelavci (2002) ne navajajo, medtem ko za racun istega primera uporabita Mofid in Akin (1996) mrezo 10 linearnih elementov (ostali numericni podatki prav tako niso navedeni). Za vse numericne izracune, razen za referencno rešitev, kjer izberemo mrezo stotih elementov, izberemo mrezo šestnajstih elementov z dvema internima tockama, numericno integracijo izvedemo s štirimi integracijskimi tockami. Izracuni potekajo z nespremenljivim casovnim korakom velikosti At = 0.001. Ob predpostavki o konstantni hitrosti delec (oziroma sila) pretece nosilec do casa 0.15, oziroma v 150 casovnih korakih. Potrebno je opozoriti, da pri numericšnem izracšunu prehoda sile predpostavimo konstantno hitrost, medtem ko pri prehodu delca podajamo le zacetno hitrost v[0] enako izbrani hitrosti v; nato se hitrost s casom spreminja v skladu z gibalno enacbo delca, vendar so spremembe ob izbiri podlage brez trenja tako majhne (odstopanja opravljene poti pri posameznih casih so velikostnega reda 10-4), da ne vplivajo na graficno natancnost. Na sliki 9.4, kjer so zbrani rezultati vertikalnih pomikov sredine nosilca, normirani z dolzino nosilca, lahko opazimo, da smo poleg numericne analize prehoda sile opravili še tri numericne analize za prehod delca, katerih rezultati se bistveno razlikujejo. Referencno rešitev predstavlja graf gosto pikcaste vijolicne krivulje; ti rezultati so skladni s predstavljeno teorijo in upoštevajo zacetno deformirano obliko nosilca - upogib zaradi lastne teze, vendar je na sliki 9.4 zaradi lazje primerjave z ostalimi rezultati prikazan le pomik zaradi dinamicne obtezbe z delcem. Rezultat, predstavljen s crtkano rumeno-rjavo crto, zacetne ukrivljene lege zaradi lastne teze nosilca ne upošteva. Oba rezultata bistveno odstopata od rezultatov iz literature [Mofid, Akin, 1996] in [Yavari et al., 2002], zato smo pri naslednjem numericnem izracunu dodatno upoštevali poenostavljen pristop avtorjev omenjenih clankov, ki ne zajema pospeška v normalni smeri. Tako redko pikcasta zelena krivulja predstavlja graf pomika sredine nosilca, kjer poleg upogiba zaradi lastne teze ne upoštevamo tudi pospeškov v normalni smeri, ki nastopajo v enacbah (9.22) in (9.25); ob upoštevanju teh poenostavitev in ob predpostavki, daje koeficient trenja enak nic, dobita enacbi (9.22) in (9.25) obliko fd : g (eg ■ et) - S'||r"|| - s2 (r" ■ et) =0 NVd = - (Fd ■ e„) e„. t Slika 9.4: Primerjava precnih pomikov sredine prostoleZecega nosilca s potujoco silo in potujocim delcem; vijolicna pikcasta crta prikazuje numericno rešitev ob upoštevanju zacetno deformiranega nosilca, v ostalih primerih je privzet zacetno raven nosilec Figure 9.4: Lateral displacements at the mid point of simply supported beam when moving particle/load is applied; the purple dotted-line represents the numerical solution considering initially deformed shape due to gravity, other numerical results consider initially straight beam Neupoštevanje dinamicnega efekta mase v normalni smeri vodi do takega precnega pomika na sredini nosilca, da se do graficne natancnosti ujema s pomiki zaradi prehoda sile; po tako poenostavljeni analizi se tudi najvecji vertikalni pomik sredine nosilca ujame z najvecjim vertikalnim pomikom sredine nosilca iz literature (glej [Mofid, Akin, 1996] in [Yavari et al., 2002]). Tudi pri numericni analizi prehoda sile (pikcasto-crtana modra krivulja) ne upoštevamo zacetnega deformiranega stanja zaradi vpliva lastne teze. S tem dosezemo ujemanje numericne rešitve z analiticno (9.42). Primerjava pomikov na sliki 9.4 kaze na precejšnjo podcenjenost pomikov na sredini nosilca v primeru prevec poenostavljenih analiz. V obravnavanem primeru dobimo po upoštevanju normalnega pospeška priblizno za 20% vecje maksimalne precne pomike sredine nosilca, kot ce pospeškov ne upoštevamo, medtem ko upoštevanje zacetno deformiranega stanja konstrukcije zaradi vpliva lastne teze pomik poveca za dodatnih 30%. Primer 2. Prostolezeci nosilec v ravnini, obtezen s premicnim delcem z maso, obravnava tudi Lee, in sicer v [Lee, 1996a] se delec giblje z enakomerno hitrostjo, v [Lee, 1996b] pa s konstantnim pospeškom. Podatki o nosilcu so E = 207 ■ 109 G = 77.6 ■ 109 L = 1 p = 7700. Lee (1996a) analizira celo paleto primerov enakomernega gibanja delca; izmed njih izberemo primer okroglega prereza s plošcino A = 0.00179 in z vztrajnostnim momentom J\ = J2 = 2.55115 ■ 10-7, za delec pa predpostavimo tri razlicne zacetne hitrosti (primeri A, B in C): va = 21.3876 VB = 97.217 vc = 213.877. Nosilec modeliramo z 32 elementi za primer A oziroma s 16 elementi za primera B in C; v vseh treh primerih za interpolacijske funkcije izberemo polinome stopnje šest, numericno integracijo izvršimo s sedmimi integracijskimi tockami. Zacetni korak je velikosti Ato = 2.5 ■ 10-4 in se spreminja v skladu z zahtevano natancnostjo etoi = 10-5. Zaustavitev Newtonove iteracije postavimo omejitev ^New = 10-8. Hitrost delca se spreminja v skladu z gibalno enacbo delca, vendar spremembe minimalno odstopajo od zacetne hitrosti delca in ne vplivajo na graficno natancnost rezultatov. 0.2 ■ ■ 0.5 0 u -0.2 0 u -0.5 -0.4 -1 -0.6 -1.5 -0.8 -2 -1 + -2.5 0.01 0.02 0.03 0.04 0.004 0.006 0.008 0.01 t — sila analitično (9.48) ■ ■ ■ ■ sila numerično — ■ — delec - - delec, lastna teža + značilna vrednost iz [Lee, 1996a] 0.001 0.002 0.003 0.004 t Slika 9.5: Primerjava rezultatov za drugi primer prostolezecega nosilca; zelena crtkana crta prikazuje rezultate v primeru, ko upoštevamo zacetno deformirano lego nosilca, pri ostalih rezultatih predpostavljamo zacetno raven nosilec Figure 9.5: The lateral displacements at the mid point of simply supported beam (case 2) when moving particle/load is applied; the green dashed-line represents the numerical solution considering initially deformed shape due to gravity, other numerical results consider initially straight beam Na grafu 9.5 so prikazani rezultati normiranih precnih pomikov u = U^L/''), kjer je za normiranje uporabljen staticni pomik us (L/2) = mgEj. Analiticni rezultati (9.42) pri enakomernem prehodu sile oznacuje polna rdeca crta, numericni izracun pri prehodu sile pa pikcasta modra crta. Pri obeh je predpostavljena ravna zacetna lega. Ostala dva grafa prikazujeta numericne rezultate pri prehodu delca: mešana vijolicna crta predpostavlja ravno zacetno lego, zelena crtkana crta pa prikazuje rezultate ob upoštevanju zacetno deformirane lege zaradi upogiba nosilca pod lastno tezo (ponovno prikazujemo in primerjamo le pomik zaradi dinamicne obtezbe z delcem). Na grafe je dodana še znacilna vrednost rezultatov iz [Lee, 1996a], kijih oznacuje crn krizec. V vseh treh primerih se numericni rezultati pri prehodu sile do graficne natancnosti ujemajo z analiticno rešitvijo. Odziv pri prehodu delca je praviloma vecji kot pri prehodu sile; ta vpliv se povecuje z vecanjem hitrosti prehoda. Veliko vecji, v primeru C celo še enkrat vecji odziv dobimo, ce upoštevamo zacetno deformirano lego nosilca zaradi lastne teze. Znacilni vrednosti iz [Lee, 1996a] se za primera B in C do graficne natancnosti ujameta z odzivom konstrukcije zaradi prehoda sile in sta v primerjavi z našimi izracuni podcenjena. 9.8.3 Delec na previsnem nosilcu Numericno rešitev prehoda delca prek enostransko togo vpetega nosilca smo zasledili le v delih Mofida in Akina (1996) ter Yavarija s sodelavci (2002), kjer je model Timoshenkovega nosilca poenostavljen v sistem togih palic in dvojnih vzmeti. Po [Yavari et al., 2002] povzemamo naslednje podatke: E = 2.068 ■ 1011 N/m2 L = 7.62 m Jt = Ji + J' Ji = J' = 4.5785 ■ 10-5 m4 G = 7.929 ■ 101U N/m p = 795.73 kg/m3 A = 5.89 ■ 10-3 m2 A1 = A2 = 4.91 ■ 10 -3 m2. »10 j(L,t) -2 -3 -4 naša rešitev [Yavari et al., 2002] _i_ 0.05 0.1 2 0 0 Slika 9.6: Primerjava rezultatov za primer prehoda delca prek previsnega nosilca Figure 9.6: Moving particle on cantilever beam: comparison of approximative solutions Yavari in sodelavci navajajo podatke v anglosaškem merilnem sistemu. Pri preracunu v standardne merske enote smo uporabili naslednje pretvorne faktorje 1 in. = 0.02540 m 1 psi = 6.894757 ■ 103 N/m2 1 lb = 4.44823 N. Delec z maso m = 53.57 kg pošljemo prek nosilca z zacetno hitrostjo v[0] = 50.8 m/s2 tako, da potuje od podprtega proti prostemu krajišcu. Podatkov o numericnem racunu v literaturi ne navajajo. Nosilec modeliramo z 32 elementi stopnje 3 s po štirimi integracijskimi tockami. Za zacetni korak izberemo Ato = 5 ■ 10-4, ki se spreminja v skladu z zahtevano lokalno natancnostjo etoi = 10 5, medtem ko za ustavitev Newtonove iteracije zahtevamo natancnost ^New = 10-8. Relativne pomike prostega krajišca prikazujemo na sliki 9.6. Primerjava z rezultatom iz literature kaze razlicen potek odziva prostega krajišca, vendar se velikostni razred pomikov ujema. Maksimalen pomik prostega krajišca konzole po poenostavljeni teoriji je nekoliko precenjen. 9.8.4 Vodni tobogani Na koncu obravnavajmo še dinamicno analizo dveh vodnih toboganov. Raven vodni tobogan. Podatke za prvi vodni tobogan crpamo iz internega porocila o staticni analizi vodnih toboganov [Zupan et al., 2005]. Tobogan je zgrajen iz polimernega materiala, kakršnega v svoji proizvodnji uporablja podjetje Veplas [Vedenik, 2005]. Deformacija pri mejni napetosti znaša emejni = 0.023. Za majhne deformacije je odvisnost med deformacijami in napetostmi linearna; raztros rezultatov enoosnih preizkusov kaze na elasticni modul E med 8 ■ 109 N/m in 1.2 ■ 1010 N/m2. Za racšun izberemo enega manj ugodnih rezultatov E = 8.2 ■ 109 N/m2. Upoštevamo strizno tog material, zato vzamemo G = 10E. Prerez toboganov je kolobar, ki ga dolocata kroga s polmeroma p1 = 0.605 m in p2 = 0.6 m (slika 9.7, levo). Za dolocitev deformacijskega stanja v precnem prerezu predpostavimo v skladu s predpostavkami teorije nosilcev linearen potek deformacij: e (y,z)= yi + + ZK2. (9.43) Raven tobogan z naklonom 15% glede na vodoravno ravnino je sestavljen iz treh enako dolgih ravnih segmentov, ki so med seboj togo povezani. Je obojestransko nepomicno clenkasto podprt tako, daje preprecen torzijski zasuk. Podpori sta v horizontalni smeri na razdalji 6 m in v vertikalni smeri na razdalji 0.9 m. Poleg obtezbe z lastno tezo (p = 2692.74 kg/m , gravitacija deluje v smeri vektorja —gg3) in tockovnih obtezb na mestu spojev F = -800 N v globalni vertikalni smeri upoštevamo še druge obtezbe, povzete po [Zupan et al., 2005]: 1. voda na toboganu: porazdeljena obtezba v smeri gravitacije velikosti 200 N/m po vsej dolzini; 2. obtežba z vetrom: porazdeljena obtezba v smeri lokalne osi, dolocene z vektorjem G2, velikosti 6000 N/m in v smeri lokalne osi, dolocene z vektorjem G3, velikosti —6000 N/m, obe vzdolz celotne tezišcne osi tobogana; I Slika 9.7: Precni prerez visokih vodnih toboganov; geometrijski podatki ravnega tobogana Figure 9.7: Cross-section of high water slides; geometrical data for straight water slide 3. obtežba (teža) uporabnikov tobogana: porazdeljena obtezba v smeri gravitacije v smeri, doloceni z vektorjem gg3, velikosti —1500N/m po vsej dolzini; 4. obtežba s snegom: porazdeljena obtezba v smeri gravitacije velikosti 2200 N/m po vsej dolzini. Najprej opravimo staticno analizo z enakimi obtezbami kot v [Zupan et al., 2005]. Obtezba uporabnikov tobogana Zupan in sodelavci upoštevajo konzervativno; kot sami navajajo, bi predpisom zadostili z upoštevanjem enake porazdeljene obtezbe zgolj na dolzini 1 m; vendar bi v tem primeru morali poiskati najmanj ugodno lego uporabnika. Za dinamicno analizo vpliv uporabnikov namesto s porazdeljeno obtezbo modeliramo s pomicnim delcem z maso 1500/g kg (g = 9.81 m/s2) in zacšetno hitrostjo v[0] = 10 m/s. Preglednica 9.1: Najvišje vrednosti deformacij in pomikov za primer ravnega tobogana Table 9.1: Straight water slide: maximum values of deformations and displacements deformacije [%] pomik [cm] statika 0.115 0.664 dinamika 0.106 0.568 [Zupan et al., 2005] 0.11 0.66 relativna razlika [%] 7.83 14.4 Model za staticno analizo tobogana izberemo enak kot Zupan in sodelavci (2005); modeliramo ga s štirimi elementi stopnje 6, za integracijo po kraju vzamemo sedem integracijskih tock; krajna segmenta modeliramo z enim, sredinskega pa z dvema elementoma (slika 9.7, desno). Dinamicno analizo opravimo z mrezo 15 elementov stopnje 6 s konstantnim korakom velikosti At = 0.04 s. Za ustavitev Newtonove iteracije postavimo zahtevo ^New = 10-8. Na sliki 9.8 prikazujemo potek normalnih deformacij prerezov in deformirano obliko tobogana za primer staticne analize, slika 9.9 pa prikazuje ovojnico normalnih deformacij prerezov po casu in koncno deformirano obliko tobogana za primer dinamicšne analize. Dinamicšno analizo smo opravili na zacšetno ukrivljeni legi, ki jo poleg teze povzrocijo še ostale stalne porazdeljene obtezbe (voda, sneg, veter). Potek deformacij je za oba nacina analize primerljiv. Najvišje vrednosti deformacij in pomikov prikazujemo Slika 9.8: Potek osnih deformacij po obodu prereza na deformirani legi za primer staticne analize ravnega tobogana; pomiki konstrukcije so pomnozeni s faktorjem 30, barvna skala oznacuje deformacije v procentih Figure 9.8: Static analysis of straight water slide: axial deformations on the deformed configuration; displacements are multiplied by factor 30, deformations are multiplied by factor 100 (in per cent) 0.1 0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1 Slika 9.9: Ogrinjaca poteka osnih deformacij na obodu prereza na deformirani legi za primer dinamicne analize ravnega tobogana; pomiki konstrukcije so pomnozeni s faktorjem 30, barvna skala oznacuje deformacije v procentih Figure 9.9: Dynamic analysis of straight water slide: evelope of axial deformations on the deformed configuration; displacements are multiplied by factor 30, deformations are multiplied by factor 100 (in per cent) v preglednici 9.1. Pri dinamicni analizi sta obe vrednosti nekoliko manjši kot pri staticni analizi, kar je pricakovano glede na to, daje teza uporabnika pri staticni analizi upoštevana kot porazdeljena obtezba vzdolz celotne dolzine tobogana. Najvišji vrednosti po staticni analizi se ujemata z vrednostima iz literature [Zupan et al., 2005]. Ukrivljen vodni tobogan. Tezišcna os ukrivljenega vodnega tobogana ima obliko osmine kroznice v navpicni ravnini s polmerom 10 m, slika 9.10. Konstrukcijo podpremo zgolj na krajišcih in sicer v levem krajišcu (vrh tobogana) poleg pomikov preprecimo torzijski zasuk (q° = 0), v desnem krajišcu (dno tobogana) pa poleg pomikov še torzijski zasuk in zasuk iz ravnine (qL = 0 in qL = 0). Gravitacija deluje v smeri, ki je nasprotna tretjemu baznemu vektorju referencnega koordinatnega sistema. Ostali geometrijski in materialni podatki so enaki kot v primeru ravnega tobogana, razen striznega modula, za katerega izberemo ustreznejšo vrednost: G = EE .Od ravnega vodnega tobogana privzamemo tudi obtezbe vkljucno z maso in zacetno hitrostjo delca, s katerim modeliramo kopalca. Konstrukcijo modeliramo s 30 ukrivljenimi elementi šeste stopnje, numericno integracijo izvedemo s sedmimi integracijskimi tockami. Ponovno locimo staticno in dinamicno analizo. Pri obeh analizah za zakljucek Newtonove iteracije postavimo zahtevo &New = 10-6. Dinamicno analizo opravimo z nespremenljivim korakom At = 0.03 s. Zacetno obliko tobogana za dinamicno analizo dolocimo s staticno analizo kot posledico obremenitve s staticnimi obtezbami (lastna teza, teza vode na toboganu in obtezba z vetrom). Bocna obtezba z vetrom povzroci prostorsko ukrivljeno zacetno lego nosilca. Nato tobogan dinamicno obremenimo še s pomicnim delcem. Potek normalnih deformacij prerezov vzdolz nosilca po staticni analizi je zelo podoben ovojnici teh deformacij po dinamicni analizi, zato potekov deformacij vzdolz nosilca graficno ne prikazujemo. V preglednici 9.2 primerjamo le največje deformacije in pomike po obeh analizah. Maksimalne deformacije in maksimalni pomiki po obeh analizah so skoraj enaki, kar je posledica dokaj nizke hitrosti kopalca, ki povzroci majhen dinamicen odziv konstrukcije, po drugi strani pa pri staticni analizi vpliv kopalca upoštevamo kar vzdolz celotne dolzine tezišcne osi. Gibanje delca primerjamo še z rezultati primera, ko za tobogan predpostavimo, daje tog nosilec. Diferencialna enacba, ki opisuje gibanje delca po kroznici, skupaj s pripadajocima zacetnima pogojema, je: St (t)= g sin ^ st (0)=0 s t (0) = 10. (9.44) R A Z Slika 9.10: Zacetna geometrija ukrivljenega tobogana Figure 9.10: Initially curved water slide: the initial configuration Zacetno nalogo (9.44) rešimo numericno z uporabo programskega paketa Mathematica [Wolfram, 2003]. Primerjavo poti, hitrosti in pospeška delca po deformabilnem in po togem toboganu (9.44) prikazujemo v preglednici 9.3. Po pricakovanjih majhne deformacije deformabilnega tobogana (glej preglednico 9.2) ne vplivajo bistveno na potek gibanja delca po toboganu. Preglednica 9.2: Najvišje vrednosti deformacij in pomikov za primer ukrivljenega tobogana Table 9.2: Curved water slide: maximum values of deformations and displacements deformacije [%] pomik [cm] statika 0.164 0.00933 dinamika 0.161 0.00931 relativna razlika [%] 1.83 0.21 Preglednica 9.3: Ukrivljen tobogan: pomiki, hitrosti in pospeški delca za primera deformabilnega in togega tobogana Table 9.3: Displacements, velocities and accelerations for a deformable and rigid curved water slide pomik hitrost pospešek cas s St v vt = St a at = St 0.9 0.91 0.90 10.12 10.04 0.73 0.88 1.8 1.82 1.81 10.23 10.16 1.68 1.77 2.7 2.75 2.73 10.42 10.36 2.62 2.65 3.6 3.70 3.68 10.70 10.64 3.54 3.53 4.5 4.68 4.65 11.06 10.99 4.51 4.40 5.4 5.70 5.66 11.51 11.43 5.43 5.26 6.3 6.76 6.71 12.04 11.94 6.31 6.10 10 ZakljuCek V disertaciji predstavljamo novo formulacijo in nove numericne metode za reševanje prostorskih linijskih konstrukcij. Poglavitne novosti predstavljene formulacije in postopkov so: • Zapis enacb za analizo prostorskih nosilcev po geometrijsko tocni Reissner-Simovi teoriji prostorskih nosilcev z elementi in operacijami kvaternionske algebre, pri cemer so rotacije para-metrizirane z rotacijskim kvaternionom. • Linearizacija enacb nosilca v kvaternionskem zapisu in algoritem za konsistentno upoštevanje linearnih popravkov rotacijskih kvaternionov. • Nova druzina koncnih elementov za staticno analizo prostorskih nosilcev po kolokacijski metodi diskretizacije, ki temelji na diskretizaciji konsistencnih enacbah v šibki obliki, kjer za osnovne neznanke izberemo pomike in rotacijske kvaternione. • Nova druzina koncnih elementov za dinamicno analizo prostorskih nosilcev po kolokacijski metodi, z izborom sšibkih konsistentnih enacšb za vodilne enacšbe problema ter pomikov in rotacijskih kvaternionov za osnovne neznanke, prirejeno za casovne integratorje druzine Runge-Kutta. • Druzina koncnih elementov za dinamicno analizo prostorskih nosilcev z enako diskretizacijo neznank in enacb po kraju ter s posplošeno casovno integracijo Newmark, prirejeno za rotacijske kvaternione. • Poseben pristop preoblikovanja robnih pogojev za reševanje zacetnega problema (po casu), ki algebrajske robne pogoje spremeni v diferencialne enacšbe po cšasu ter socšasno ohranja simetrijo koncnega elementa po kraju. • Razširitev zgornjih formulacij s povezanim reševanjem odziva konstrukcije ob hkratnem vplivu potujocega delca z maso. Poglavitne prednosti predstavljene kvaternionske parametrizacije rotacij so: • Kvaternionska parametrizacija rotacij odpravlja dvojnost rotacij, kiju v klasicnem pristopu z vektorsko parametrizacijo rotacij predstavljajo trije rotacijski parametri in rotacijska matrika, saj nam rotacijski kvaternioni sluzijo kot edina z rotacijami povezana kolicina. Rotacijski kvaternioni nam sluzijo za parametrizacijo rotacij, za rotacijo precnih prerezov in za koordinatno transformacijo. • Kvaternionski zapis enacb, skupaj z naravo rotacijskih kvaternionov, omogoca preprostejšo linearizacijo enacb v primerjavi z vektorsko parametrizacijo rotacij. Znacšilnosti in prednosti predstavljenih koncšnih elementov so: • Koncni elementi omogocajo poljubno ukrivljeno in zvito zacetno lego prostorskega nosilca. • Koncni elementi so objektivni in nimajo tezav s striznim blokiranjem ne glede na izbiro integracije. • Poleg nelinearne geometrije omogocajo vgradnjo razlicnih nelinearno-elasticnih materialov. • Dinamicno analizo lahko opravimo z dvema bistveno razlicnima postopkoma integracije po casu: - z uporabo ze razvitih integratorjev Runge-Kutta, kar omogoca izjemno enostavno implementacijo in vse prednosti teh metod; - z Newmarkovo integracijo za prostorsko rotacijo, posebej prirejeno za potrebe parametriza-cije rotacij z rotacijskimi kvaternioni. • Nov numericni postopek za povezano analizo dinamicnega odziva nosilca pri prehodu delca z maso ali pri prehodu sile ob hkratni nelinearni geometriji nosilca, kjer - omogocamo poljubno zacetno lego, poljubne zacetne pogoje in poljubno podpiranje nosilca; - dopušcamo moznost zacetne staticne analize pred dinamicno analizo za izracun zacetno deformirane lege zaradi lastne tezše in drugih stalnih obtezšb; - upoštevamo prepletenost enacb prek lege delca in medsebojnih sil med delcem in konstrukcijo. • Predstavljeni koncni elementi predstavljajo svez pristop na podrocju prostorskih nosilcev in odpirajo nove moznosti za razvoj podrocja numericnega modeliranja konstrukcij. 11 Povzetek V disertaciji smo prikazali teoreticno izpeljavo in numericno implementacijo nove druzine koncnih elementov za staticno in dinamicno analizo prostorskih linijskih konstrukcij, osnovano na kinematicno tocni teoriji prostorskih nosilcev. Linijsko konstrukcijo smo modelirali s tezišcno osjo in druzino precnih prerezov. Enacbe modela in izpeljani numericni algoritem omogocajo poljubno zacetno ukrivljeno in zvito obliko in nepravokoten kot med tezišcno osjo in prerezi nosilca. Za osnovne neznanke problema smo izbrali kinematicne kolicine, pomike tezišcne osi v obicajni vektorski parametrizaciji in zasuke precšnih prerezov v kvaternionski parametrizaciji. Zaradi izbrane kvaternionske parametrizacije rotacij smo enacšbe prostorskega nosilca lahko v celoti zapisali z operacijami in elementi kvaternionske algebre in se popolnoma izognili uporabi rotacijske matrike in rotacijskega vektorja. Tako smo se izognili dvojnosti rotacijska matrika-rotacijski vektor, ki je znacšilna za klasicšne elemente, in s tem povezanih pretvorb, kar je povecšalo racšunsko ucšinkovitost in preglednost posameznih izrazov. V delu smo podrobno predstavili algebro kvaternionov in izpeljali povezave med kvaternionsko in v teorijah nosilcev bolj uveljavljeno parametrizacijo rotacij z rotacijskim vektorjem. Enacbe smo reševali numericno po metodi koncnih elementov. Enacbe smo diskretizirali po metodi kolokacije, osnovne neznanke pa inter-polirali skozi diskretne vrednosti s polinomi poljubnih stopenj. Z izbiro izoparametricne interpolacije za vse štiri komponente rotacijskega kvaterniona in ob uporabi normalizacije v vmesnih, interpoliranih tockah, smo zagotovili ohranjanje objektivnosti deformacij. Tocke diskretizacije in kolokacije smo poenotili in s tem dodatno izboljšali natancnost elementa. V notranjih kolokacijskih tockah smo zadostili šibkim konsistencnim enacbam, kar pomeni, da smo zahtevali enakost odvodov rezultantnih ravnoteznih in materialnih sil in momentov v prerezu; na robu pa smo zadošcali ravnoteznim enacbam. Sistem diskretnih enacb za staticno analizo linijskih konstrukcij smo reševali z uporabo Newtonove metode. Izpeljali smo postopek za konsistentno upoštevanje linearnih iterativnih popravkov rotacijskih kvaternionov in ostalih rotacijskih kolicin. Popravek, kije dobljen z Newtonovo iteracijo in je element tangentnega prostora, je v splošnem neenotski kvaternion in neposredno ne doloca popravka rotacije. S pravilnim upoštevanjem lastnosti kolicin, ki so povezane z zasuki, izpeljemo postopek za multiplikativni popravek, ki je rotacijski kvaternion in ga lahko dodamo predhodnemu rotacijskemu kvaternionu z ustreznim kvaternionskim produktom. Konsistentno popravljanje odvodov rotacijskih kvaternionov preko aditivnosti upogibnih in torzijskih deformacij v ustreznih bazah prispeva k višji natancnosti elementa, omogoca objektivnost deformacij in neodvisnost numericne metode od poti obtezevanja. Linearizacija enacšb v kvaternionskem zapisu je bistveno manj zahtevna kot pri uporabi vektorske parametrizacije rotacij, ki vodi v Liejevo grupo in zahteva uporabo smernih odvodov v nelinearnih prostorih. Ustreznost racunskega algoritma za staticno analizo konstrukcij smo prikazali na znanih testnih primerih iz literature. S primeri smo pokazali splosšnost formulacije (poljubna zacšetna geometrija, poljubna stopnja interpolacije) ter natancšnost in ucšinkovitost predstavljenega algoritma. Numericšne simulacije so pokazale tudi neobcutljivost postopka na strizno blokiranje in neodvisnost rezultatov od postopka nalaganja konservativne obtezšbe. Za potrebe dinamicne analize smo robne ravnotezne enacbe preoblikovali v diferencialno obliko tako, da v njih nastopajo osnovne spremenljivke odvajane po casu. S tem smo se izognili numericno obcutljivemu reševanju diferencialno-algeebrajskega sistema enacb. Postopek diskretizacije po kraju smo povzeli po staticni analizi, za diskretizacijo in integracijo po casu pa smo uporabili dva zelo razlicna postopka. Za prvi postopek integracije po casu smo izbrali dve uveljavljeni metodi druzine Runge-Kutta, kot sta implementirani v programskem okolju Matlab. (Čeprav so klasicne metode Runge-Kutta namenjene reševanju diferencialnih enacb v aditivnih konfiguracijskih prostorih, se za izpeljano numericno implementacijo izkazeta ti dve metodi kot dovolj natancni tudi pri zahtevnejših prostorskih primerih. Vzrok temu je ucinkovitost in robustnost zapisa enacb v kvaternionski algebri. Prednost tega pristopa je v njegovi preprostosti, saj diskretizacija kolicin in linearizacija enacb po casu nista bili potrebni, poleg tega nismo potrebovali posebnega postopka za popravljanje rotacijskih kolicin. S tem smo se popolnoma izognili zahtevnim matematicnim strukturam in prostorom, ki se navezujejo na rotacijske operatorje in njihove parametrizacije. Metodi imata vgrajeno avtomatsko preverjanje lokalne napake, ki neposredno vpliva na dolzino casovnega koraka integracije. V literaturi tako preprostega algoritma za dinamicno analizo kinematicno tocnih prostorskih nosilcev nismo zasledili. Metoda uspešno rešuje tako probleme, ki vkljucujejo velike pomike in rotacije, kot tiste, ki vkljucujejo velike deformacije. Tocnost rezultatov smo preverili z rezultati drugih avtorjev. Drugi nacin integracije enacb predstavlja povezavo predstavljene metode koncnih elementov s casovnimi integratorji, ki so bili razviti posebej za dinamiko prostorskih nosilcev. To je izjemnega pomena, saj so razlicni avtorji na podrocju prostorskih nosilcev za metode, ki temeljijo na parametrizaciji rotacij z rotacijskim vektorjem, ze razvili številne casovne integratorje. Izbrali smo casovni integrator tipa Newmark, ki sta ga predstavila Simo in Vu-Quoc (1988), ter ga priredili za uporabo na rotacijskih kvaternionih. Ta pristop sicer iznici nekatere pozitivne lastnosti rotacijskih kvaternionov, na primer ponovno se pojavi dvojnost rotacijski operator-parametrizacija. Kljub temu so numericne študije in primerjave z drugimi avtorji pokazale, da lahko casovne integratorje, razvite za vektorsko parametrizacijo rotacij, brez tezav priredimo za uporabo na rotacijskih kvaternionih. Z obema metodama casovne integracije smo izracunali referencne primere iz literature. Rezultati po obeh metodah casovne integracije so medsebojno primerljivi, prav tako se ujemajo z rezultati drugih avtorjev. Ugotovili smo, da izbira integratorja mocno zaznamuje casovno in prostorsko zahtevnost numericne implementacije in natancnost rezultatov. Regulacija casovnega koraka prek preverjanja lokalne napake (Runge-Kutta) prispeva k lokalno natancnejšim rezultatom, ki v rešitev zajamejo tudi višje nihajne oblike, vendar se lahko socšasno bistveno povecša prostorska in cšasovna zahtevnost algoritma. Z uporabo integratorja, ki ne preverja lokalne napake in ob socasni izbiri fiksnega casovnega koraka, ki ne zajame višjih nihajnih oblik (Newmark), pa lahko dosezemo stabilnejši racun na daljšem casovnem obmocju, dobljeni rezultati pa dobro dolocajo globalno obnašanje konstrukcije. V zadnjem delu doktorske naloge smo obravnavali potovanje delca po prostorskem nosilcu. Predpostavili smo, da delec drsi po tezišcni osi nosilca skladno z gibalno enacbo delca in ne s predpisano hitrostjo ali pospeškom. Gibalno enacbo delca dodamo diskretnemu sistemu enacb za dinamicno analizo nosilca ob izbiri posplošene Newmarkove integracijske sheme, izbranim osnovnim kinematicnim neznankam pa dodamo novo neznanko - dolzino opravljene poti delca. Neznanke so v razširjenem sistemu enacb mocno prepletene in reševanje razširjenega sistema enacb predstavlja povezan problem. Pri formulaciji gibalne enacbe smo upoštevali, da parameter tezišcne osi v deformirani legi ni enak naravnemu parametru nedeformirane osi. Enacbe smo zato zapisali glede na krivocrtni Frenetov koordinatni sistem. Za modeliranje vpliva trenja smo uporabili Coulombov zakon. Racunski model nosilca omogoca razlicna podpiranja in poljubno zacetno obliko, delec pa ima poljubno maso in zacetne hitrosti. Racunski postopek zaradi primerjave z analiticnimi rezultati in z drugimi avtorji priredimo še za primer prehoda sile s konstantno hitrostjo. Numericni rezultati se dobro ujemajo z dosegljivimi analiticnimi rešitvami. V primerjavi z drugimi avtorji, ki uporabljajo preprostejše modele, v nekaterih primerih opazimo bistveno vecji dinamicni odziv konstrukcije. 12 Summary In the thesis new formulation of the kinematically exact three-dimensional beam theory for the static and dynamic analysis of beam structures has been theoretically derived and numerically implemented. Three-dimensional beam has been modeled by the line of centroids and by the family of cross-sections. The derived theoretical model along with the new algorithm enable us to consider an arbitrary initial shape of the beam. The kinematic quantities, i.e. displacements in a classical vector parametrization and rotations in the quaternion parametrization, have been chosen as the primary unknowns of the problem. For this reason we fully abandoned the rotational vector concept, and rewrote the system of the governing equations of the beam in the quaternion algebra description. In such a description of rotations, we completely avoided the rotation matrix-rotational vector duality. Consequently, the numerical efficiency and transparency of equations have been improved considerably. The quaternion algebra is presented in detail and the relations between quaternion and the classical rotational vector parametriza-tion of rotation have been derived. New finite element procedure for solving beam equations have been proposed. The collocation type of the method for the discretization of the equations has been used. We interpolated primary unknowns with polynomials of an arbitrary degree. The isoparametric interpolation of all four rotational quaternion components together with the normalization procedure leads to the objective strains. The locations of the discretization and the collocation points have been made equal which additionally improved the accuracy of the method. At the internal collocation points, the weak consistency conditions (the equality of the derivatives of the equilibrium stress resultants and the constitutive stress resultants) have been satisfied. At the two beam boundaries, the equilibrium equations have been satisfied. The system of discrete nonlinear equations of the statics problem has been solved by Newton's method. New consistent update procedure for linear corrections of rotational quaternions and for other rotational quantities has been derived. The update of the rotational quaternion, as obtained from Newton's procedure, is an element of the tangent space and is thus not a rotational (unit) quaternion. Therefore a multiplicative update procedure for rotational quaternions has been derived on the basis of the properties of the spatial rotations and their linearization. A consistent update of derivatives of rotational quaternions is based on the additivity of the rotational strain, yet with respect to appropriate bases only. The proposed update of the rotational quaternions and their derivatives leads to a higher accuracy of results and it seems to be crucial for the conservation of the objectivity of strain measures and path independency of conservative systems. Moreover, the linearization of equations using quaternion algebra has been shown to be considerably simpler compared to classical approaches, where the Lie algebra is used. The numerical algorithm has been validated by numerical examples given in literature. Numerical tests demonstrate the ability of the formulation to consider properly initially curved and twisted configurations of the beam and an arbitrary order of the finite element. We have presented the high accuracy and efficiency of the proposed method as well as the resistance to shear locking and path dependency. For the dynamics analysis, the boundary equilibrium equations have been transformed into the weak differential form. By such an approach we have completely avoided the mixed differential-algebraic system of equations which can be computationally instable. Discretization of the equations with respect to the centroidal axis is taken as in the static analysis. For the time discretization and the numerical integration of equations with respect to time, two independent approaches were applied. The first one employs two different methods among the Runge-Kutta family, as implemented in the Matlab environment. Despite the fact that the standard Runge-Kutta methods were developed for solving problems in additive configuration spaces, the two Runge-Kutta methods have proven to be capable of solving demanding three-dimensional beam problems when applied in the proposed quaternion-based finite-element discretization. The main advantage of this approach is in its simplicity, as the special linearization of the equations and the update procedure are not needed. A local error verification is automatically incorporated in the Runge-Kutta methods which results in a variable length of the time step. We have successfully solved demanding beam problems of finite rotations, displacements and strains. The results compare well with the solutions in literature. In another approach used here, we have modified the solver, developed previously specially for the three-dimensional beam dynamics. We have demonstrated that various solvers developed for the rotational-vector beam formulations can be suitably rearranged for the use with rotational quaternions. The Ne-wmark solver, presented by Simo and Vu-Quoc (1988), has been here modified to deal with the quaternion formulation. Despite the fact that some of the advantages of the quaternion formulation, such as a unique description of rotations without needing matrices, are lost, we have shown that it is possible to employ, with only tehnical modifications, a wide range of solvers developed for rotation-vector based beam formulations. The numerical results showed that the two approaches are comparable. The present results are in good agreement with the literature. The choice of the numerical solver, however, affects considerably the computational efficiency of the overall algorithm. By controlling the local error, the Runge-Kutta solvers can incorporate the influence of higher-order mode shapes. This can lead to high computational times; however the results include local and global behaviour of the structure. Solvers using fixed time steps such as Newmark's can behave numerically stable over longer time interval and well describe the global behaviour of the structure. The last chapter of the thesis deals with the moving mass problem. A mass particle is moving along the axis of the three-dimensional beam. Its behaviour is described by the equation of motion resulting in the velocity and the acceleration being dependent on the shape and the oscillation of the supporting beam and on the direction of the gravity. The equation of motion of the particle is added to the system of the beam dynamic equations and a new unknown is added, describing the arc-coordinate of the curvilinear path. As the contact forces are dependent on the deformation of the beam and as the length of the centroidal axis varies with beam deformation, the equations are strongly coupled. The equation of motion of the particle was expressed with respect to the Frenet frame and Coulomb's law was used to model the friction. The present numerical model is capable of considering arbitrary initial geometry and various types of supports. For comparison reasons, the model was further simplified to consider the moving force with a constant speed. Numerical results are in an excellent agreement with the analytical solutions. Comparisons with other authors and simplified models have shown that the consideration of initial displacements and non-linear behaviour of the beam can result in a significantly intensive dynamic response of the structure. Literatura Andelic, T., Stojanovic, R. 1965. Racionalna mehanika. Beograd, Zavod za izdavanje udzbenika SR Srbije. Argyris, J.H., Dunne,P.C., Malejannakis, G., Scharpf, D. W. 1978. On large displacements - small strain analysis of structures with rotational degrees of freedom. Comput. Methods Appl. Mech. Engrg. 14: 99-135. Argyris, J.H. 1982. An excursion into large rotations. Comput. Methods Appl. Mech. Engrg. 32: 85-155. Argyris J., Poterasu, V.F. 1993. Large rotations revised application of Lie algebra. Comput. Methods in Appl. Mech. Engrg. 103: 11-42. Ascione, L., Feo, L., Mancusi, G. 2000. On the statical behaviour of fibre-reinforced polymer thin-walled beams. Compos. Part B-Eng. 31: 643-654. Atluri, S. N., Cazzani, A. 1995. Rotations in computational solid mechanics. Arch. Comput. Methods Eng. 2: 49-138. Bathe, K.J., Bolourchi, S. 1979. Large displacement analysis of three-dimensional beam structures. Int. J. Numer. Meth. Eng. 14: 961-986. Bathe, K.J. 1996. Finite Element Procedures. Prentice-Hall International, Inc. Battini, J., Pacoste, C. 2002. Co-rotational beam elements with warping effects in instability problems. Comp. Methods Appl. Mech. Engrg. 191: 1755-1789. Bauchau, O.A.,. Damilano, G, Theron, N. J. 1995. Numerical integration of non-linear elastic multi-body systems. Int. J. Numer. Meth. Eng. 38: 2727-2751. Bauchau, O.A., Theron, N. J. 1996. Energy decaying scheme for non-linear beam models. Comput. Methods Appl. Mech. Engrg. 134: 37-56. Belytschko, T., Liu, W.K., Moran, B. 2000. Nonlinear Finite Elements for Continua and Structures. Chichester-New York-Weinheim-Brisbane-Singapore-Toronto, John Wiley & Sons. Betsch, P., Steinmann, P. 2002. Frame-indifferent beam finite elements based upon the geometrically exact beam theory. Int. J. Numer. Meth. Eng. 54: 1775-1788. Bohte, Z. 1991. Numerične metode. Ljubljana, Društvo matematikov, fizikov in Slovenije. Bottasso, C.L. 1997. A new look at finite elements in time: a variational interpretation of Runge-Kutta methodes. Appl. Numer. Math. 25: 355-368. Bottasso, C.L., Borri, M. 1997. Energy preserving/decaying schemes for non-linear beam dynamics using the helicoidal approximation. Comput. Methods Appl. Mech. Engrg. 143: 393-415. Bottasso, C.L., Borri, M. 1998. Integrating finite rotations. Comput. Methods Appl. Mech. Engrg. 164: 307-331. Bronštejn, J.N., Semendjajev, K.A. 1988. Matematižni prirožnik. prevedel A. Zabkar, Ljubljana, Tehnisška zalozšba Slovenije. Butcher, J.C. 1987. The Numerical Analysis of Ordinary Differential Equations. Runge-Kutta and General Linear Methods. Chichester-New York-Brisbane-Toronto-Singapore, John Wiley & Sons. Cardona, A., Geradin, M. 1988. A beam finite element non-linear theory with finite rotations. Int. J. Numer. Meth. Eng. 26: 2403-2438. Crisfield, M.A. 1990. A consistent co-rotational formulation for non-linear, three-dimensional beam elements. Comp. Methods Appl. Mech. Engrg. 81: 131-150. Crisfield, M.A. 1997. Non-linear Finite Element Analysis of Solids and Structures, Volume 2, Advanced Topics. Chichester-New York-Weinheim-Brisbane-Singapore-Toronto, John Wiley & Sons. Crisfield, M.A., Jelenic, G. 1999. Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation. Proc. R. Soc. Lond. 455: 1125-1147. Crivelli, L.A., Felippa, C.A. 1993. The three-dimensional non-linear Timoshenko beam based on the core-congruential formulation. Int. J. Numer. Meth. Eng. 36: 3647-3673. Dormand, J.R., Prince, P.J. 1980. A family of embedded Runge-Kutta formulae. J. Comp. Appl. Math. 6: 19-26. Esmailzadeh, E., Ghorashi, M. 1995. Vibration analysis of beams traversed by uniform partially distributed moving masses. J. Sound Vib. 184: 9-17. Esmailzadeh, E., Ghorashi, M. 1997. Vibration analysis of a Timoshenko beam subjected to a travelling mass. J. Sound Vib. 199: 615-628. Evans, G.A. 1995. Practical Numerical Analysis. Chichester - New York - Brisbane - Toronto -Singapore, John Wiley&Sons. Gams, M., Planinc, I., Saje M. 2007. Energy conserving time integration scheme for geometrically exact beam. Comput. Methods Appl. Mech. Engrg 196: 2117-2129. Geradin, M., Rixen, D. 1995. Parametrization of finite rotations in computational dynamics: a review. Revue europeenne des elements finis 4: 497-553. Ghosh, S., Roy, D. 2008. Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam. Comput. Methods Appl. Mech. Engrg. 198: 555-571. Gerald, C.F., Wheatley, P.O. 1994. Applied Numerical Analysis: fifth edition. Addison-Wesley publishing company. Hairer, E., Wanner, G. 1991. Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Berlin-Heidelrberg-New York, Springer-Verlag. Hosea, M.E., Shampine, L.F. 1996. Analysis and implementation of TR-BDF2. Appl. Num. Math. 20: 21-37. Ibrahimbegovic, A., Frey, F. 1993. Finite element analysis of linear and nonlinear deformations of elastic initially curved beams. Int. J. Numer. Methods Eng. 36: 3239-3258. Ibrahimbegovic, A. 1995. On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements. Comput. Methods Appl. Mech. Engrg. 122: 11-26. Ibrahimbegovic, A., Frey, F., Kozar, I. 1995. Computational aspects of vector-like parametrization of three-dimensional finite rotations. Int. J. Numer. Methods Eng. 38: 3653-3673. Ibrahimbegovic, A. 1997. On the choice of finite rotation parameters. Comput. Methods Appl. Mech. Engrg. 149: 49-71. Ibrahimbegovic, A., al Mikdad, M. 1998. Finite rotations in dynamics of beams and implicit time-stepping schemes. Int. J. Numer. Meth. Eng. 41: 781-814. Ibrahimbegovic, A., Mamouri, S. 1999. Nonlinear dynamics offleksible beams in planar motion: formulation and time-stepping scheme for stiff problems. Comput. & Struct. 70: 1-22. Iura, M., Atluri, S.N. 1989. On a consistent theory and variational formulation of finitely stretched and rotated 3-D space-curved beams. Comput. Mech. 4: 73-88. 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. Engrg. 120: 131-161. Jelenic, G., Crisfield, M.A. 1999a. Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics. Comput. Methods Appl. Mech. Engrg. 171: 141171. Jelenic, G., Crisfield, M.A. 1999b. Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation. Proc. R. Soc. Lond. A 455: 11251147. Johnson, S. M., Williams, J. R., Cook, B. K. 2007. Quaternion-based rigid body rotation integration algorithms for use in particle methods. Int. J. Numer. Meth. Eng. doi: 10.1002/nme.2210. Krizanic, F. 1990. Temelji realne matematične analize. Ljubljana, Drzavna zalozba Slovenije. Krizanic, F. 1993. Linearna algebra in linearna analiza. Ljubljana, Drzavna zalozba Slovenije. Lee, H.P. 1996a. The dynamic response of a Timoshenko beam subjected to a moving mass. J. Sound Vib. 198: 249-256. Lee, H.P., 1996b. Transverse vibration of a Timoshenko beam acted on by an accelerating mass. Appl. Acoust. 47: 319-330. Leyendecker, S., Betsch, P., Steinmann, P. 2006. Objective energy-momentum conserving integration for the constrained dynamics of geometrically exact beams. Comput. Methods Appl. Mech. Engrg. 195: 2313-2333. MacNeal, R.H., Harder, R.L. 1985. A proposed standard set of problems to test finite element accuracy. Finite Elem. Anal. Design 1: 3-20. Martys, N.S., Mountain, R. D. 1999. Velocity verlet algorithm for dissipative-particle-dynamics. Physical Review E 59: 3733-3736. Mata, P., Oller, S., Barbat, A.H. 2007. Static analysis of beam structures under nonlinear geometric and constitutive behavior Comput. Methods Appl. Mech. Engrg. 196: 4458-4478. Mata, P., Oller, S., Barbat, A.H. 2008 Dynamic analysis of beam structures considering geometric and constitutive nonlinearity. Comput. Methods Appl. Mech. Engrg. 197: 857-878. McRobie, F.A., Lasenby, J. 1999. Simo-Vu Quoc rods using Clifford algebra. Int. J. Numer. Methods Eng. 45: 377-398. Michaltsos, G., Sophianopoulos, D., Kounadis, A.N. 1996. The effect of a moving mass and other parameters on the dynamic response of a simply supported beam. J. Sound Vib. 19: 357-362. Mofid, M., Akin, J.E. 1996. Discrete element response of beams with traveling mass. Adv. Eng. Softw. 25: 321-331. Muršic, M. 1972. Uvod vkinetiko konstrukcij. Ljubljana, Fakulteta za arhitekturo, gradbeništvo in geodezijo, Univerza v Ljubljani. Nour-Omid, B., Rankin, C.C. 1991. Finite rotation analysis and consistent linearization using projectors. Comput. Methods Appl. Mech. Engrg. 93: 353-384. Park, S., Youm, Y. 2001. Motion of a moving elastic beam carrying a moving mass-analysis and experimental verification. J. Sound Vib. 240: 131-157. Phillips, W.F., Hailey, C.E., Gebert, G.A. 2001. Review of attitude representations used for aircraft kinematics. J. Aircraft 38: 718-737. Poreous, I.R. 1995. Clifford Algebras and the Classical Groups. Cambridge, Cambridge University Press. Reissner, E. 1981. On finite deformation of space-curved beams. J. Appl. Math. Phys. 32: 734-744. Ritto-Correa, M., Camotim, D. 2002. On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner-Simo beam theory. Int. J. Numer. Methods Eng. 55: 1005-1032. Romero, I., Armero, F. 2002. An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy-momentum conserving scheme in dynamics. Int. J. Numer. Methods Eng. 54: 1683-1716. Romero, I. 2004. The interpolation of rotations and its application to finite element models of geometrically exact rods. Comput. Mech. 34: 121-133. Saje, M. 1994. Kinematika in dinamika. UL FGG, Ljubljana. Saje, M., Srpcic, S. 1986. Large deformations of thin curved plane beam of constant initial curvature. Int. J. Mech. Sci. 28: 275-287. Shoemake, K. 1985. Animating rotation with quaternion curves. Computer Graphics (ACM) 19: 245254. Siddiqui, S.A.Q., Golnaraghi, M. F., Heppler, G. R. 2000. Dynamics of a flexible beam carrying a moving mass using perturbation, numerical and time-frequency analysis techniques. J. Sound Vib. 229: 1023-1055. Siddiqui, S.A.Q., Golnaraghi, M. F., Heppler, G. R. 2003. Large free vibrations of a beam carrying a moving mass. Int. J. Nonlinear Mech. 38: 1481-1493. Simo, J.C. 1985. A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Comput. Methods Appl. Mech. Engrg. 49: 55-70. Simo, J.C., Vu-Quoc, L. 1986. A three-dimensional finite-strain rod model. Part II: Computational aspects. Comput. Methods Appl. Mech. Engrg. 58: 79-116. Simo, J.C., Vu-Quoc, L. 1988. On the dynamics in space of rods undergoing large motions - A geometrically exact approach. Comput. Methods Appl. Mech. Engrg. 66: 125-161. Simo, J.C., Wong, K.K. 1991. Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. Int. J. Numer. Methods Eng. 31: 19-52. 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. Spurrier, R.A. 1978. Comment on "Singularity-free extraction of a quaternion from a direction-cosine matrix". J. Spacecraft 15: 255. Stanek, M., Turk, G. 1998. Osnove mehanike trdnih teles. Ljubljana, Fakulteta za gradbeništvo in geodezijo, Univerza v Ljubljani. The MathWorks, Inc. 1999. MATLAB, Using MATLAB, Natick, http://www.mathworks.com. Timoshenko, S.P., Gere, J.M. 1961. Theory of Elastic Stability New York, McGraw-Hill. Vedenik, U. 2005. Optimizacija proizvodnje elementov vodnih toboganov v družbi Veplas d.d. Diplomsko delo, Maribor, Fakulteta za gradbenisštvo, Unoverza v Mariboru. Vidav, I. 1989. Diferencialna geometrija. Ljubljana, Društvo matematikov, fizikov in astronomov Slovenije. Wang, Y.-M. 1998. The dynamical analysis of a finite inextensible beam with an attached accelerating mass. Int. J. Solids Struct. 35: 831-854. Ward, J. P. 1997. Quaternions and Cayley Numbers. Dordrecht-Boston-London, Kluwer Academic Publishers. Wolfram, S. 2003. Mathematica. Addison - Wesley Publishing Company. Wu, J.J., Whittaker, A.R., Cartmell, M.P. 2000. The use of finite element techniques for calculating the dynamic response of structures to moving loads. Comput. Struct. 78: 789-799. Wu, J.J., Whittaker, A.R., Cartmell, M.P. 2001. Dynamic responses of structures to moving bodies using combined finite element and analytical methods. Int. J. Mech. Sci. 43: 2555-2579. Xu, X., Xu, W., Genin, J. 1997. A non-linear moving mass problem. J. Sound Vib. 204: 495-504. Yavari, A., Nouri, M., Mofid, M. 2002. Discrete element analysis of dynamic response of Timoshenko beams under moving mass. Adv. Eng. Softw. 33: 143-153. Yun, X., Bachmann, E. R. 2006. Design, implementation, and experimental results of a quaternion-based Kalman filter for human body motion tracking. IEEE Trans. on Robotics 22: 1216-1227. Zupan, D. 2003. Rotacijsko invariantne deformacijske koližine v geometrijsko tožni teoriji prostorskih nosilcev. Doktorska disertacija, Ljubljana, Fakulteta za gradbeništvo in geodezijo, Univerza v Ljubljani. Zupan, D., Saje, M. 2003. Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures. Comput. Methods Appl. Mech. Engrg. 192: 52095248. Zupan, D., Saje, M. 2004. On A proposed standard set of problems to test finite element accuracy: The twisted beam. Finite Elem. Anal. Design 40: 1445-1451. Zupan, D. Cas, B., Planinc, I., Saje, M. 2005. Statična (geometrijsko in materialno) nelinearna analiza vodnih toboganov. interno porocilo, Ljubljana, Fakulteta za gradbeništvo in geodezijo, Univerza v Ljubljani. Zupan, D., Saje, M. 2006. The linearized three-dimensional beam theory of naturally curved and twisted beams: The strain vectors formulation. Comput. Methods Appl. Mech. Engrg. 195: 4557-4578. Zupan, E., Saje, M., Zupan, D. 2009. The quaternion-based three-dimensional beam theory Comput. Methods Appl. Mech. Engrg. 198: 3944-3956. Priloga Zupan, E., Saje, M., Zupan, D. 2009. The quaternion-based three-dimensional beam theory. Comput. Methods Appl. Mech. Engrg. 198: 3944-3956. Dodatki A Odvod, integral in variacija Naj oznaka V pomeni tri ali štirirazsezen vektorski prostor nad obsegom realnih števil, V V V izberimo referencno bazo (—i (x) , i = (0), 1, 2, 3) z izhodišcem v tocki O in druzino pomicnih baz Gi (x) , i = (0), 1, 2, 3 ), odvisnih od parametra x G A C IR. A1 Odvod in variacija operatorja Rotacija R je preslikava na trirazseznem, rotaciji in pa sta preslikavi na štirirazseznem vektorskem prostoru. Poleg rotacij imamo še druge operatorje, na primer antisimetricni operator S in njegovo štirirazsezno razširitev S ter eksponentno preslikavo kvaternionskega argumenta exp. Za odvajanje operatorjev veljajo drugacna pravila kot za odvajanje navadne funkcije. Definicija 2 Operator F : V ^ V, definiran na odprti mnočici A C V,je odvedljiv v točki a G A, če obstaja tak e > 0, da za vse b G V, | b | < e velja F (a + - = F (a) + F- b + O (a,b) , (A.1) kjerje F^ zvezen linearni operator in za O (a,b) velja lim O a, b b| 0. Potem F-t imenujemo krepki ali Frechetov odvod operatorja F v točki a. Vektor F-t variacija operatorja F v točki— v smeri b in označimo z 5F. imenujemo Ker definicija krepkega odvoda ni prakticna za konkreten izracun odvoda operatorja, vpeljemo še smerni odvod. Definicija 3 Naj za vse b G IR3 obstaja limita lim a0 F (- + aM — F (a) ^ a da F (a + ab a=0 b Vrednost limite imenujem o smerni ali Gateauxov odvod operatorja F v točki " v smeri b, kar kratko zapišemo z oznako DF— b VT- d_ da F (a + ab a=0 (A.2) Izrek 1 Ce je operator F odvedljiv v točki " G A, v tej točki obstaja smerni odvod in je enak variaciji operatorja v tej točki. Ker se tu omejujemo zgolj na odvedljive operatorje, ki imajo krepki in smerni odvod enak, bomo lahko variacije operatorjev racunali po formuli (A.2). Primer 1 Oglejmo si odvod zveznega linearnega operatorja F, kije vselej odvedljiv. Zato je smerni odvod enak krepkemu. Za linearen operator velja Po primerjavi z enačbo (A.1) sledi J— F " + b) = F(") + F b) . = T b in O ",b =0 torej je krepki odvod linearnega operatorja kar operator sam, variacija linearnega operatorja v toCki " v smeri b pa kar njegova vrednost v točki b bT = T (b) in J— = T. Primer 2 Zanimiv rezultat dobimo tudi pri krepkem odvodu identičnega operatorja F (") = "a, za katerega velja F (a + b j = "a + b, po definiciji smernega odvoda pa F (a + ab J — F (a) DJ- lim lim a — + ab — "a a = Ker je identiteta odvedljiv operator, po izreku 1 sledi, da se smerni odvod identitete ujema z variacijo bT = DF" = b, zato vektor b imenujemo kar variacija vektorja b" = b. Podobno kot variacijo operatorjev racunamo tudi variacijo vektorskih funkcij vektorskega argumenta x ^ "(s) . Njihov smerni odvod izracunamo kot b" = lim "(" + ab-) — "(s) = d ["(" + ab")]a=0 . a^0 a d^ a=0 b b b A2 Odvod in variacija vektorske funkcije skalarnega argumenta Odvod poljubne vektorske funkcije r (x) skalarnega argumenta po x v tocki x = x0 je definiran s predpisom dr , r (xo + Ax) — r (x0) — = lim ---, dx Ax^o Ax kadar limita obstaja. Vektorsko funkcijo r (x) lahko zapišemo v referencni ali v pomicni bazi r (x) = E rgi (x) gi = E rGi (x) Gi (x), i i kjer i vselej pretece števila od 0 ali 1 do 3, i = (0), 1, 2, 3. Odvod vektorske funkcije po x je neodvisen od izbire baze, zato je dr (x) sr^ drgi (x) — dx dx gi i ^ drGi (x) - dGi (x) = £ (x) + E rGi ii Odvod vektorske funkcije po x v zapisu s pomicno bazo je sestavljen iz dveh delov: iz spremembe komponent, brez upoštevanja spremembe baze; ta del imenujemo relativni odvod vektorja (Dre, = S ^ (x) • 'A 3) drugi del pa upošteva samo spremembo pomicne baze in ga zato imenujemo sistemski odvod vektorja (IL = E"Gi (x) **(A.4) Tako lahko odvod vektorske funkcije x zapišemo kot vsoto relativnega in sistemskega odvoda dr / dr \ ^ / dr \ (A 5) dx Vdx) rei Vdx) sist. . Variacijo vektorja smo izrazili kot variacijo identicnega operatorja, glej primer 2. Kadar vektorje razvijemo po pomicni bazi, se variacije vektorskih funkcij izrazajo popolnoma enako kot odvodi in sicer kot vsota relativne in sistemske variacije, ki skupaj tvorijo celotno ali absolutno variacijo Sr = (čf)rei + (č-) sist (A.6) (čf)rei = SrGiGi + Srg2G2 + 5rg3G3 (A.7) (Sr)sist = rGi sGi + rg2 SG2 + rg3 §G3. (A.8) B Eksponentna preslikava kvaternionskega argumenta Ker imamo v tem poglavju opravka le z izrazavo v referencni bazi, njen indeks (g) kar izpustimo. B1 Eksponentna oblika rotacijskega kvaterniona Vzemimo rotacijski kvaternion 1 = cos | + ■ sin | v polarnem zapisu. V njem ■ predstavlja vektor na osi rotacije velikosti $ (rotacijski vektor). Pokazali bomo, da lahko 1 zapišemo z neskoncno potencno vrsto ■ 1 ■ ■ 1 ■ ■ ■ ■ 1 (■) = 1 +—---1—-— o--1—-— o — o--h ... = exp — 1! 2 2! 2 2 3! 2 2 2 \ 2 (B.1) v kateri je ■ cisti kvaternion rotacijskega vektorja, ■ = [0 $]T. Izraz (B.1) obicajno imenujemo eksponentna oblika rotacijskega kvaterniona, operator exp (x) pa eksponentna preslikava kvaternionskega argumenta. Za potrditev pravilnosti izraza (B.1) uporabimo lastnosti kvaternionskega potenciranja cistih kvaternionov (4.11) na kvaternionu ■, na primer ■ o ■ = -$2 ■ o (d o = -$2$. (B.2) (B.3) Lastnosti (B.2)-(B.3) uporabimo v vrsti (B.1) in dobimo 1 1 1 /$\2 1 $21 1 /$ 1 = 1 + TT2 - 2!V2/ - 3! 231 + 41V2 = 1 - A 2! V 2 1 /$ + 4! 2 + 4 1/1 $ 1 /$x3 + ... + $ \ 1T2 - 3! \2 + $ 1 $ = COS2+ $ Sin2. Faktor s™//2, ki nastopa v polarnem zapisu rotacijskega kvaterniona (4.22), je za $ = 0 singularen, a koncen v okolici tocke 0. Ker je numericno popolnoma stabilno in tocno izracunljiv, lahko brez tezav za preracunavanje iz rotacijskih vektorjev v rotacijske kvaternione uporabimo kar polarno obliko (4.22). Izraz (B.1) je primeren tudi za racun multiplikativnega popravka rotacijskega kvaterniona Aq iz iterativnega popravka rotacijskega vektorja 5$: _____ . ~ 151 151 51 151 51 51 /51 A1 (5$) = 1 + —--h—— o--h—— o — o--h ... = exp — ' 1! 2 2! 2 2 3! 2 2 2 \ 2 (B.4) i T kjer je 51 = [ 0 5$ JT . Analogen izraz za multiplikativni popravek rotacijskega kvaterniona Aq, vendar tokrat v odvisnosti od iterativnega popravka rotacijskega kvaterniona 51, dobimo v obliki A1 (5q) = exp (5q o 1*), v kateri smo upoštevali zvezo med 51 in 51 po enacbi (4.57). Zapis multiplikativnega popravka s potencno vrsto (B.4) v odvisnosti od 51 pa je A1 (51) = 1 + 151 o 1 * + 2 51 o 1 * o 51 o 1 * + 3_ 51 o 1 * o 51 o 1 * o 51 o 1 * + 4 2 B2 Odvod in variacija rotacijskega kvaterniona {Prvi odvod rotacijskega kvaterniona g po x dobimo z odvajanjem vrste (B.1). Pri tem moramo dosledno uposštevati nekomutativnost kvaternionskega produkta: w ( « „n 1 0 1 0 ti 0 \ g' (0,0') = -7 — + ^ — o - + - o — 1ir 2M 2 2 2 2 J 1 /0' 0 0 0 0' 0 0 0 0'\ +--- — o — o--1--o — o--1--o — o — + .. 3! \ 2 2 2 2 2 2 2 2 2 (B.5) Ti1 + 2! M 0 + M 00 + UR(0o 0 +