Fakulteta za gradbeništvo in geodezijo
A-L4
, ,TT| 1 1 i|. iTTl i 1111111I1 m|
I TTTTltTTlTTTTi
PETER ČEŠAREK, univ. dipl. inž. grad.
DINAMIKA PROSTORSKIH LINIJSKIH ELEMENTOV Z INTERPOLACIJO DEFORMACIJSKIH KOLIČIN
DOKTORSKA DISERTACIJA
Ljubljana, april 2013
Univerza	
v Ljubljani	
	Fakulteta za
	gradbeništvo in
	geodezijo
PODIPLOMSKI ŠTUDIJ GRADBENIŠTVA
DOKTORSKI ŠTUDIJ
Kandidat:
PETER ČEŠAREK, univ. dipl. inž. grad.
DINAMIKA PROSTORSKIH LINIJSKIH ELEMENTOV Z INTERPOLACIJO DEFORMACIJSKIH KOLIČIN
Doktorska disertacija štev.: 231
DYNAMICS OF SPATIAL BEAMS WITH INTERPOLATION OF STRAIN MEASURES
Doctoral thesis No.: 231
Temo doktorske disertacije je odobrila Komisija za doktorski študij na 20. redni seji, dne 21. septembra 2011. Za mentorja je bil imenovan izr. prof. dr. Dejan Zupan, za somentorja prof. dr. Miran Saje.
Ljubljana, 16. april 2013
Univerza	
v Ljubljani	
	Fakulteta za
	gradbeništvo in
	geodezijo
i i ■ ri i i i i i i1 nI 11 IIi
::TiliTO:
Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi:
-	izr. prof. dr. Dejan Zupan,
-	prof. dr. Miran Saje
-	prof. dr. Jože Korelc,
-	izr. prof. dr. Gordan Jelenic, Gradbena fakulteta Univerze na Reki,
je imenoval Senat Fakultete za gradbeništvo in geodezijo na 20. redni seji, dne 20. aprila 2011.
Poročevalce za oceno doktorske disertacije v sestavi:
-	izr. prof. dr. Gordan Jelenic, Gradbena fakulteta Univerze na Reki,
-	prof. dr. Jože Korelc,
-	prof. dr. Igor Planinc,
je imenoval Senat Fakultete za gradbeništvo in geodezijo na 37. redni seji, dne 30. januarja 2013.
Komisijo za zagovor doktorske disertacije v sestavi:
-	prof. dr. Matjaž Mikoš, dekan UL FGG, predsednik,
-	izr. prof. dr. Dejan Zupan, mentor,
-	prof. dr. Miran Saje, somentor,
-	izr. prof. dr. Gordan Jelenic, Gradbena fakulteta Univerze na Reki,
-	prof. dr. Jože Korelc,
-	prof. dr. Igor Planinc,
je imenoval Senat Fakultete za gradbeništvo in geodezijo na 39. redni seji, dne 27. marca 2013.
Univerza	
v Ljubljani	
	Fakulteta za
	gradbeništvo in
	geodezijo
I I ■ ri I i I I I ni il 11 IIi
::TiliTO:
IZJAVA O AVTORSTVU
Podpisani Peter ČEŠAREK, univ. dipl. inž. grad., izjavljam, da sem avtor doktorske disertacije z naslovom:
»DINAMIKA PROSTORSKIH LINIJSKIH ELEMENTOV Z INTERPOLACIJO DEFORMACIJSKIH KOLIČIN«.
Izjavljam, da je elektronska različica v vsem enaka tiskani različici.
Izjavljam, da dovoljujem objavo elektronske različice v repozitoriju UL FGG.
Ljubljana, 16. april 2013
(podpis)
POPRAVKI
Stran z napako
Vrstica z napako	Namesto
Naj bo
BIBLIOGRAFSKO-DOKUMENTACIJSKA STRAN IN IZVLEČEK
UDK
Avtor:
Mentor:
Somentor:
Naslov:
Tip dokumenta: Obseg in oprema: Ključne besede:
624.074(043.3)
Peter Češarek
izr. prof. dr. Dejan Zupan
prof. dr. Miran Saje
Dinamika prostorskih linijskih elementov z interpolacijo deformacijskih količšin
Doktorska disertacija
114 str., 48 sl., 7 pregl., 260 en.
prostorski nosilec, točna kinematika, dinamika, časovna integracija, masni delec, trk
Izvlecek
V doktorski nalogi predstavimo nove metode za dinamicno analizo prostorskih linijskih konstrukcij. Obravnavamo casovno integracijske sheme in formulacije koncnih elementov, zasnovane na geometrijsko tocni teoriji prostorskih nosilcev. Bistvena novost predstavljenih algoritmov je izbira odvodov kinematicnih kolicin - pomikov in rotacij, za osnovne neznanke pri reševanju problema v odvisnosti od parametrov prostora in casa. Za casovno integracijo dinamicnih enacb predstavimo Newmarkovo shemo, prirejeno za nelinearen prostor rotacij, in predlagamo novo - modificirano Newmarkovo shemo. Osnovne neznanke nove sheme so odvodi pomikov in rotacij po casu - hitrosti in kotne hitrosti. Predstavimo dve formulaciji linijskih koncnih elementov. Za osnovne neznanke v prvi formulaciji izberemo krajevne odvode pomikov in rotacij - deformacije, v drugi formulaciji pa krajevne odvode hitrosti in kotnih hitrosti. Osnovne neznanke interpoliramo, za diskretizacijo zveznih enacb po kraju pa izberemo kolokacijsko metodo. Pri deformacijskem elementu kinematicne zveze po casu diskretiziramo z Newmarkovo shemo, za zagotovitev stabilnosti pa uporabimo posplošeno-a metodo numericnega dušenja. Pri elementu z interpolacijo krajevnih odvodov hitrosti, kinematicne enacbe po casu diskretiziramo z modificirano Newmarkovo shemo. Enacbe rešujemo z Newtonovo iteracijsko metodo. Za obe formulaciji podrobno predstavimo linearizacijo in dodajanje popravkov. Nove elemente vgradimo v racunalniški program in testiramo na številnih primerih. Na koncu model nosilca povezemo z gibalnimi enacbami masnega delca. Podrobneje obravnavamo trk masnega delca ob nosilec.
BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION AND ABSTRACT
UDC
Author:
Supervisor:
Co-supervisor:
Title:
Document type:
Notes:
Keywords:
624.074(043.3) Peter Češarek
assoc. prof. Dejan Zupan, Ph.D. prof. Miran Saje, Ph.D.
Dynamics of spatial beams with interpolation of strain measures Doctoral Dissertation
114 p., 48 fig., 7 tab., 260 eq.
spatial beam, exact kinematics, dynamics, time integration, mass particle, impact
Abstract
In the present thesis we develop new numerical methods for a dynamic analysis of beam-like structures in three-dimensional space. We consider time integration schemes and finite-element formulations, ba-sed on geometrically exact beam theory. The main novelty of the present approach is the choice of the derivatives of displacements and rotations, as the primary unknowns when considering the problem with respect to space and time. For time integration of dynamic equations we present the Newmark scheme extended to non-linear space of rotations and propose a new - modified Newmark scheme. The modified Newmark scheme employs time derivatives of displacements and rotations - velocities and angular velo-cities, as the only unknowns of the scheme. We introduce two finite-element formulations for analysis of spatial beams. The first formulation employs spatial derivatives of displacements and rotations - strains, as the primary unknowns, while the second formulation employs spatial derivatives of velocities and angular velocities as the only unknown functions. The primary unknowns are interpolated and a collocation method is chosen for discretization of the continuous equations. The Newmark scheme is used for time discretization of kinematic quantities and Generalized-a method is used to assure numerical stability of the strain-based formulation. The velocity-based formulation employs modified Newmark scheme for time discretization of kinematic equations. Equations are solved by the iterative Newton's method, therefore the linearization of equations and the update procedure are presented. The computer code is generated and the performance of the formulations is tested on several numerical examples. Finally, the beam model is coupled with equations of motion of a moving particle. The impact of mass particle against the beam is studied in detail.
ZAHVALA
Zahvaljujem se: Javni agenciji za raziskovalno dejavnost Republike Slovenije za dodeljena financna sredstva, ki so mi omogocila nadaljevanje študija. Mentorju, somentorju in vsem ostalim sodelavcem na katedri za mehaniko za pomoc pri študiju, nastajanju tega dela in za vsakdanje stvari zaradi katerih sem rad opravljal svoje delo. Posebej pa mojim bliznjim za razumevanje in podporo.
Hvala vsem, naj vam bo z dobrim povrnjeno.
KAZALO VSEBINE
BIBLIOGRAFSKO-DOKUMENTACIJSKA STRAN IN IZVLEČEK	VIII
BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION AND ABSTRACT	IX
ZAHVALA	X
1	UVOD	1
1.1	Pregled stanja na podrocju dinamike prostorskih nosilcev................................1
1.2	Vsebina dela................................................................................3
2	OSNOVNE ENACBE DINAMIKE	5
2.1	Osnovni izreki..............................................................................5
2.2	Kinematika delca..........................................................................6
3	PROSTORSKE ROTACIJE	8
3.1	Rotacijska matrika........................................................................8
3.2	Kotna hitrost in kotni pospešek............................................................9
3.2.1 Kompozitum rotacij in seštevanje kotnih hitrosti ................................10
3.3	Odvod in integral vektorske funkcije v pomicšni bazi ......................................12
3.4	Parametrizacija rotacij ....................................................................14
3.4.1	Rotacija s konstantno kotno hitrostjo..............................................15
3.4.2	Zveza med kotno hitrostjo in odvodom rotacijskega vektorja po casu............16
3.4.3	Variacija rotacijske matrike, parametrizirane z rotacijskim vektorjem............17
3.4.4	Variacija rotacijske matrike, parametrizirane z aditivnim rotacijskim vektorjem	18
4	Casovno integracijske sheme	20
4.1 Newmarkova shema........................................................................21
4.1.1	Linearizacija......................................................................21
4.1.2	Dodajanje popravkov..............................................................22
4.1.3 Izbira parametrov Newmarkove sheme............................................23
4.2	Modificirana Newmarkova shema........................................................23
4.2.1	Linearizacija......................................................................25
4.2.2	Racun novih priblizkov............................................................26
4.2.3	Izbira parametrov modificirane Newmarkove sheme ............................26
4.3	Numericni primeri gibanja togega telesa..................................................28
4.3.1	Vrtenje telesa pod vplivom predpisanega poteka momentov ....................28
4.3.2	Tezka vrtavka v gravitacijskem polju ............................................31
4.3.3	Dvojno nihalo ......................................................................34
5	OSNOVNE ENACBE DINAMIKE PROSTORSKIH NOSILCEV	40
5.1	Matematicni model prostorskega nosilca..................................................40
5.2	Kinematicne enacbe........................................................................40
5.3	Gibalne enacbe............................................................................42
5.4	Konstitutivne enacbe......................................................................44
5.5	Enacbe konsistence........................................................................44
5.6	Povzetek enacb............................................................................44
5.7	Zveze med deformacijami in hitrostmi....................................................45
6	NUMERICNO REŠEVANJE ENACB	48
6.1	Kinematicne enacbe........................................................................48
6.1.1	Newmarkova shema, osnovne neznanke so deformacije..........................49
6.1.2	Modificirana Newmarkova shema, osnovne neznanke so hitrosti................50
6.2	Ravnotezne enacbe........................................................................52
6.3	Vodilne enacbe dinamike prostorskega linijskega elementa..............................53
6.3.1	Element z interpolacijo deformacij ................................................54
6.3.2	Element z interpolacijo krajevnih odvodov hitrosti..............................54
6.4	Reševanje diskretnih enacb................................................................55
6.4.1	Element z interpolacijo deformacij................................................56
6.4.2	Element z interpolacijo krajevnih odvodov hitrosti..............................62
6.4.3	Numericna integracija............................................................68
6.5	Numericno dušenje........................................................................69
6.6	Numericni testi............................................................................71
6.6.1	Vsiljeno nihanje prostoleZecega nosilca..........................................73
6.6.2	Konzola upognjena v spiralo......................................................75
6.6.3	Kolenasta konzola ................................................................76
6.6.4	Nepodprti podajni nosilec v prostem letu ........................................79
6.6.5	Problem štiriclenskega mehanizma................................................83
6.6.6	Enostavno nihalo..................................................................89
7	GIBANJE MASNEGA DELCA IN NOSILCA	93
7.1	Osnovne enacbe............................................................................93
7.2	Racunski primeri..........................................................................94
7.2.1	Interpolacija osnovnih neznank nosilca s konstantnimi funkcijami..............94
7.2.2	Nihajocša gugalnica ................................................................95
7.2.3	Trk masnega delca ob nosilec ....................................................99
8	ZAKLJUČEK	104 POVZETEK	106 SUMMARY	108 VIRI	110
KAZALO SLIK
2.1 Lega poljubnega delca trdnega telesa........................... 6
3.1 Kompozitum rotacij.................................... 11
4.1	Rotirajoce telo pod vplivom linearnih rotacij. Primerjava absolutnih napak komponent rotacijskega vektorja dobljenih z originalnim in modificiranim Newmarkovim algoritmom............................................ 29
4.2	Rotirajoce telo pod vplivom linearnih rotacij. Modificirani Newmarkov algoritem. Absolutna napaka komponent rotacijskega vektorja..................... 29
4.3	Rotirajoce telo pod vplivom kvadratnih rotacij. Primerjava absolutnih napak komponent rotacijskega vektorja dobljenih z originalnim in modificiranim Newmarkovim algoritmom............................................ 30
4.4	Tezkavrtavka v gravitacijskem polju......................................................32
4.5	Tezka vrtavka. Absolutna napaka celotne energije vrtavke................................32
4.6	Tezka vrtavka. Konvergenca druge norme vrtilne kolicine in rotacijske matrike..........33
4.7	Dvojno nihalo..............................................................................34
4.8	Dvojno nihalo. Rezultati dobljeni z Newmarkovo in modificirano Newmarkovo shemo,
z relativno velikim casovnim korakom h = 10-3..................... 36
4.9	Dvojno nihalo. Rezultati dobljeni z Matlabovo rutino ode15s.............. 37
4.10	Dvojno nihalo. Rezultati dobljeni z modificirano Newmarkovo shemo s kratkim casovnim korakom h = 10-6..................................... 38
5.1	Matematicni model nosilca................................ 41
5.2	Obtezba nosilca...................................... 43
6.1	Vsiljeno nihanje prostolezecega nosilca.......................... 73
6.2	Vsiljeno nihanje prostolezecega nosilca. Absolutne napaka pomika na sredini nosilca na intervalu t G [0,8] (h = casovni korak, ne = število elementov, N = število interpolacijskih tock)....................................... 74
6.3	Konzola upognjena v spiralo............................... 75
6.4	Konzola upognjena v spiralo. Rezultati analize s 25 elementi z 8 tockovno interpolacijo krajevnih odvodov hitrosti................................. 76
6.5	Kolenasta konzola..........................................................................77
6.6	Kolenasta konzola. Pomik kolena in prostega konca konzole iz ravnine pri analizi s 4 elementi....................................................................................78
6.7	Kolenasta konzola. Pomik kolena in prostega konca konzole iz ravnine pri analizi z 20 elementi. ..................................................................................78
6.8	Kolenasta konzola. Pomiki kolena in prostega konca konzole iz ravnine, pri analizi s korakom h = 0.1............................................................................79
6.9	Kolenasta konzola. Mehanska energija....................................................79
6.10	Nepodprti podajni nosilec v prostem letu..................................................80
6.11	Podajni nosilec v prostem letu. Projekcije deformiranih leg..............................80
6.12	Podajni nosilec v prostem letu. Pomiki obtezenega vozlišca..............................81
6.13	Podajni nosilec v prostem letu. Mehanska energija........................................82
6.14	Stiriclenski mehanizem....................................................................83
6.15	Enoosni clenek..............................................................................84
6.16	Stiriclenski mehanizem. Pomiki vozlišca C..............................................84
6.17	Stiriclenski mehanizem. Notranje sile in momenti na sredini prvega nosilca..............85
6.18	Stiriclenski mehanizem. Hitrosti in kotne hitrosti na sredini prvega nosilca..............85
6.19	Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti
-	pomiki vozlišca C........................................................................86
6.20	Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti
-	notranje sile in momenti na sredini prvega nosilca......................................87
6.21	Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti
-	hitrosti in kotne hitrosti na sredini prvega nosilca........................................88
6.22	Enostavno nihalo..........................................................................89
6.23	Enostavno nihalo. Čelotna energija nihala v odvisnosti od casa..........................90
6.24	Enostavno nihalo. Absolutne napake pomikov in hitrosti..................................91
7.1	Nihajoca gugalnica........................................................................95
7.2	Nihajoca gugalnica. Zaporedje deformiranih leg..........................................96
7.3	Nihajoca gugalnica. Pomiki v tocki B, izracunani z elementom s konstantnimi deformacijami in posplošeno-a shemo..........................................................97
7.4	Nihajoca gugalnica. Čelotna energija masnega delca in gugalnice, diskretizirane z elementi s konstantnimi deformacijami......................................................97
7.5 Nihajoca gugalnica. Pomiki tocke B in celotna energija sistema pri analizi z elementi s
konstantnimi krajevnimi odvodi hitrosti in korakom h = 5•10-5........................98
7.6	Nihajoca gugalnica. Notranje sile v nosilcu, levo od tocke A............................98
7.7	Matematicni model trka masnega delca ob nosilec........................................99
7.8	Trk masnega delca ob nosilec..............................................................100
7.9	Rezultati analize trka masnega delca ob prostolezeci nosilec..............................101
7.10	Rezultati analize trka masnega delca ob nepomicno podprti nosilec......................102
KAZALO PREGLEDNIC
4.1	Parametri casovno integracijskih shem uporabljeni v numericnih primerih..............28
4.2	Rotirajoce telo pod vplivom linearnih rotacij. Konvergenca Newtonove iteracije pri
t = 12......................................................................................30
4.3	Rotirajoce telo pod vplivom kvadratnih rotacij. Konvergenca Newtonove iteracije pri
t = 12......................................................................................31
4.4	Tezka vrtavka. Racunski casi..............................................................33
5.1 Pregled enacb za dinamiko prostorskega nosilca v skladu z geometrijsko tocno teorijo
nosilcev. ..................................................................................45
6.1	Parametri casovno integracijskih shem uporabljeni v numericnih primerih..............72
6.2	Pomik prostega krajišca konzole v smeri obtezbe pri casu t = 10........................76
LIST OF FIGURES
2.1 The position of an arbitrary particle of a solid body.......................................6
3.1 Čomposite of rotations......................................................................11
4.1	Rotating body under linear rotations. Comparison of absolute error of components of
the rotational vector obtained with original and modified Newmark algorithm............29
4.2	Rotating body under linear rotations. Modified Newmark algorithm. Absolute error of components of the rotational vector........................................................29
4.3	Rotating body under quadratic rotations. Comparison of absolute error of components
of the rotational vector obtained with original and modified Newmark algorithm. ...	30
4.4	The heavy symmetrical top in gravity field................................................32
4.5	The heavy symmetrical top. Absolute error of the total energy od the top................32
4.6	The heavy symmetrical top. Convergence of the second norm of the angular momentum
and rotation matrix..........................................................................33
4.7	Double pendulum............................................................................34
4.8	Double pendulum. Results obtained by Newmark scheme and modified Newmark scheme, with relatively large time step h = 10-3..........................................36
4.9	Double pendulum. Results obtained by the Matlab ode15s routine........................37
4.10	Double pendulum. Results obtained with modified Newmark scheme with small time
step h = 10-6................................................................................38
5.1	Mathematical model of the beam..........................................................41
5.2	Loading of the beam........................................................................43
6.1	Forced vibration of simply supported beam................................................73
6.2	Forced vibration of simply supported beam. Absolute error of midspan deflection on the interval t G [0,8] (h = time step, ne = number of elements, N = number of interpolation points)......................................................................................74
6.3	Čantilever bent to a helical form............................................................75
6.4	Čantilever bent to a helical form. Results of the analysis with 25 elements with 8 point interpolation of spatial derivatives of vlocities..............................................76
6.5	Right-angle cantilever......................................................................77
6.6	Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the cantilever obtained with analysis with 4 elements..........................................78
6.7	Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the cantilever obtained with analysis with 20 elements........................................78
6.8	Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the cantilever, obtained with the time step h = 0.1............................................79
6.9	Right-angle cantilever. Mechanical energy..................................................79
6.10	Flexible beam in free flight..................................................................80
6.11	Flexible beam in free flight. Side views of deformed shapes..............................80
6.12	Flexible beam in free flight. Displacements of the loaded end of the beam................81
6.13	Flexible beam in free flight. Mechanical energy............................................82
6.14	Four bar mechanism........................................................................83
6.15	Revolute joint................................................................................84
6.16	Four bar mechanism. Displacements at point C............................................84
6.17	Four bar mechanism. Internal forces and moments at the midpoint of the first beam. . .	85
6.18	Four bar mechanism. Velocities and angular velocities at the midpoint of the first beam.	85
6.19	Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of velocities - displacements at point C......................................................86
6.20	Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of velocities - internal forces and moments at the midpoint of the first beam................87
6.21	Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of velocities - velocities and angular velocities at the midpoint of the first beam............88
6.22	Simple Pendulum............................................................................89
6.23	Simple pendulum. Total energy of the pendulum versus time..............................90
6.24	Simple pendulum. Absolute error of displacements and velocities........................91
7.1	The swing pendulum........................................................................95
7.2	The swing pendulum. Sequence of deformed shapes of the swing........................96
7.3	The swing pendulum. Displacements at point B, obtained by constant-strain elements
and Generalized-a scheme..................................................................97
7.4	The swing pendulum. Total energy of the mass and swing, discretized with constant strain elements..............................................................................97
7.5	The swing pendulum. Displacements of point B and total energy of the system obtained
with elements with constant spatial derivatives of velocities and time step h = 10-5. . .	98
7.6	The swing pendulum. Internal forces on the left-hand side of point A....................98
7.7	Mathematical model of an impact of a mass particle against beam........................99
7.8	Impact of a mass particle against beam....................................................100
7.9	Results of the analysis of the impact of a mass particle against simply supported beam.	101
7.10	Results of the analysis of the impact of a mass particle against fully supported beam. .	102
LIST OF TABLES
4.1	Parameters of the time integration schemes used in numerical examples..................28
4.2	Rotating body under linear rotations. Čonvergence of Newton iteration at t = 12. ...	30
4.3	Rotating body under quadratic rotations. Čonvergence of Newton iteration at t = 12. .	31
4.4	The heavy symmetrical top. Computational time..........................................33
5.1	Summary of equations for dynamics of three-dimensional beam according to geometri-
cally exact beam theory.....................................................................45
6.1	Parameters of the time integration schemes used in numerical examples..................72
6.2	Displacement of the free-end of the cantilever in the direction of load at time t = 10. .	76
1 UVOD
Razumevanje in napoved obnašanja konstrukcij zaradi razlicnih vplivov je motiv za nenehno raziskovalno delo na podrocju gradbenih konstrukcij. Pri iskanju odgovorov se lahko zatecemo k eksperimentom, ceneje in z veliko vecjo ponovljivostjo pa lahko odgovore išcemo z racunskim, matematicnim modeliranjem konstrukcij. Taka analiza konstrukcij je primerna samo, ce ima raziskovalec na voljo dovolj natancen matematicen model in orodja za reševanje. Prakticna aplikacija razumevanja obnašanja konstrukcij se odraza v zahtevah za projektiranje in analizo ze obstojecih zgradb. Sodobni predpisi na podrocju projektiranja uveljavljajo zahteve, ki jim lahko zadostimo samo z dovolj natancno analizo. Obe perspektivi obravnave konstrukcij tako narekujeta nenehen razvoj sodobnih metod za analizo konstrukcij. Prva zahteva pri razvoju je natancnost modelov konstrukcij, ta pa s sabo pogosto prinese racunsko obcutljivost postopka. Zato je robustnost analize kljucna za splošno uporabnost modela.
Mnoga podrocja raziskovanja in projektiranja konstrukcij zahtevajo natancno dinamicno analizo. Primere najdemo na podrocju gradnje mostov, pri študiju vplivov vozil in ostalih vibracij na konstrukcije, pri analizi konstrukcijskih elementov v razlicnih strojih, pri rušenju objektov. Pri projektiranju konstrukcij je še posebej pomembna natancna analiza izrednih dogodkov, kot so udarci objektov v stavbe in potresna obremenitev konstrukcij. V gradbenih konstrukcijah pogosto nastopajo elementi, pri katerih ena dimenzija prevladuje nad ostalimi. Takim elementom pravimo nosilci in jih pogosto srecamo tudi v strojni, letalski in vesoljski industriji. Matematicni modeli linijskega nosilca vse kolicine reducirajo na tezišcno os nosilca; s tem se znatno zmanjša število neznank problema in poveca ucinkovitost numericnih metod. Zato so linijski elementi pomembni gradniki racunskega modeliranja konstrukcij in primerni za razvoj novih numericnih metod. Mnoge metode, razvite za linijske elemente, lahko razširimo tudi na podrocje analize lupin in tridimenzionalnih teles.
V	doktorskem delu prvi korak k natancni dinamicni analizi opravimo z upoštevanjem prostorske narave obnašanja konstrukcij in pripravimo orodja za tridimenzionalno modeliranje linijskih konstrukcij. V primerjavi s poenostavljenimi ravninskimi modeli konstrukcij tu nastopijo prvi zapleti, ki jih povzroa zahteven opis prostorskih rotacij. Nadalje se oddaljimo od priljubljenega razmišljanja, da so pomiki konstrukcij majhni (linearna teorija) in tezimo k upoštevanju cim bolj tocnih (nelinearnih) geometrijskih zvez. V sklopu izpeljav numericnih algoritmov za dinamicno analizo konstrukcij dopušcamo vkljucitev poljubnih (realnih) materialnih modelov, algoritme pa kasneje preizkusimo na primerih konstrukcij iz linearno elasticnih materialov. Pri dinamicni analizi konstrukcij rešujemo problem v odvisnosti od parametrov prostora in casa. Temeljno vprašanje vsake take analize je izbira osnovnih neznank problema. S tem vprašanjem se sprva ukvarjamo pri problemu casovne integracije, za tem pa pri analizi prostorskih nosilcev.
1.1 Pregled stanja na področju dinamike prostorskih nosilcev
V	strokovni literaturi lahko najdemo razlicne teorije nosilcev. Temeljno delo na podrocju tocne (nelinearne) kinematike nosilcev je predstavil Reissner [1], ki zveze med deformacijskimi in kinematicnimi
kolicinami zapiše tako, da je zadošceno principu virtualnega dela, pri tem pa ne omeji velikosti pomikov in deformacij. Prvo uspešno parametrizacijo in numericno implementacijo Reissnerjeve teorije najdemo v delih Sima, samostojno [2] in skupaj z Vu-Quocom [3,4]. Ker predstavljeni pristop nima poenostavitev geometrije, ga imenujemo 'geometrijsko tocna teorija' linijskih nosilcev. Teorijo so za osnovo numericnih formulacij linijskih elementov privzeli številni avtorji na podrocju statike: Cardona in Geradin [5], Ibrahimbegovic [6], Jelenic in Saje [7], Jelenic in Crisfield [8], Zupan in Saje [9] in dinamike: Simo s sodelavci [10], Bauchau in Theron [11], Bottasso in Borri [12], Ibrahimbegovic in Al Mikdad [13], Romero in Armero [14].
Pri numericnih formulacijah prostorskih linijskih nosilcev (pa tudi lupin in tridimenzionalnih teles) je še posebej pomembna izbira parametrizacije prostorskih rotacij. Prostorske rotacije lahko enolicno opišemo s tremi neodvisnimi parametri, primer take parametrizacije pa je zapis zasukov z rotacijskim vektorjem [15]. Kljub temu da tak opis ni brez singularnosti [16], se je izkazal za zelo primernega in ga pogosto srecšamo v formulacijah koncšnih elementov.
V	povezavi z izbiro parametrizacije rotacij je v formulacijah po metodi koncnih elementov pomembna izbira osnovnih neznank problema. Omenjena dela s podrocja linijskih elementov, z izjemo prispevkov, ki so jih objavili Jelenic in Saje [7] ter Zupan in Saje [9,17], za osnovne neznanke izberejo pomike in zasuke. Taka izbira osnovnih neznank v formulacije koncšnih elementov vnese nekaj tezšav. Prva tezava je pojav 'blokiranja' koncnih elementov, ki jo avtorji rešujejo na razlicne nacine. Primer uspešne formulacije, v kateri ne pride do blokiranja koncnih elementov, sta predstavila Jelenic in Saje [7], kjer za edine neznanke problema izbereta zasuke. Druga tezava je multiplikativna narava rotacij, zaradi katere rotacijskih parametrov ne moremo interpolirati z aditivnimi nastavki. Zaradi tega namesto interpolacije celotnih zasukov avtorji v svoje formulacije vpeljejo interpolacijo infinitezimalnih popravkov [3] ali interpolacijo inkrementalnih popravkov rotacij [5]. Tretji problem formulacij, osnovanih na pomikih in zasukih, na katerega sta opozorila Crisfield in Jelenic [18], je neobjektivnost v smislu ohranjanja deformacijskih kolicin pri togih premikih. Problem sta rešila z vpeljavo posebnih interpolacijskih funkcij za zasuke.
Teoreticno povsem drugacen pristop je izbira deformacijskih kolicin za osnovne neznanke formulacije po metodi koncnih elementov. To izbiro motivira dejstvo, da so deformacije invariantne kolicine pri togih premikih in da z njimi opisujemo deformacijsko energijo konstrukcije. S tako izbiro osnovnih neznank so odpravljene tudi vse zgoraj omenjene slabosti formulacij osnovanih na pomikih in zasukih. Koncni element za ravninske nosilce, osnovan na izbiri deformacijskih kolicin, je predstavil Planinc [19]. Uspešno, tocno in numericno ucinkovito formulacijo prostorskih linijskih koncnih elementov, osnovano na deformacijah, sta predstavila Zupan in Saje [9,17]. Njuno delo nas motivira, da omenjeni pristop uporabimo sše za dinamicšno analizo prostorskih nosilcev.
V	trirazseznem prostoru je krajevni odvod rotacije enak produktu transformacijske matrike, odvisne od rotacij, in vektorja ukrivljenosti. Formulacije, osnovane na deformacijskih kolicinah, se med sabo razlikujejo ravno v nacinu reševanja te zahtevne kinematicne vezi. Čhoi in Lim [20] ter Friedman in Kosmatka [21] rotacije dolocajo iz predpostavljenega poteka ukrivljenosti po elementu, Schulz in Fili-ppou [22] pa reševanje poenostavita z interpolacijo infinitezimalnih ukrivljenosti in zasukov. V predlaganem delu sledimo pristopu, ki sta ga uveljavila Zupan in Saje [9,17], in upoštevamo tocne kinematicne zveze, ki jih v nekaterih primerih zadostimo analiticno, v splošnem pa rešujemo z orodji za numericno
reševanje diferencialnih enacb.
Sistem enacb gibanja v nelinearni mehaniki deformabilnih teles je navadno tog, zato ze majhne oscilacije osnovnih neznank povzrocijo buren odziv in nihanja z visokimi frekvencami. Stabilne numericne metode za reševanje dinamicnih problemov morajo zato vsebovati sposobnost kontroliranja visokih frekvenc. Osnovna pristopa sta t.i. numericno dušenje Hilber et al. [23], Čhung in Hulbert [24] in ohranjanje mehanske energije [25] ter gibalne in vrtilne kolicine. V dinamiki prostorskih nosilcev avtorji vecinoma uporabljajo implicitne metode casovne integracije. Poišcemo lahko številne formulacije na osnovi Eu-lerjeve diskretizacije po casu, ki zadošcajo enacbam nosilca v vmesni tocki: npr. Bottasso in Borri [12], Ibrahimbegovic in Al Mikdad [13], Jelenic in Črisfield [8], Romero in Armero [14], Simo in Wong [26], Simo et al. [10], Ibrahimbegovic in Mamouri [27]. Avtorji radi posezejo tudi po osnovnih ali modificiranih Newmarkovih metodah [28]: npr. Simo in Vu-Quoc [4], Simo et al. [10]. Izjema med formulacijami za prostorske nosilce je delo Bottass-a in Borri-ja [29], ki za reševanje dinamicnih enacb uporabita metode druzine Runge-Kutta. Pri študiju omenjenih formulacij lahko izpostavimo, da v daljših casovnih intervalih analize razlicnih avtorjev dajejo obcutno razlicne rezultate. Eden od razlogov je zagotovo nelinearnost prostorskih rotacij.
Multiplikativna narava rotacij v trirazseznem prostoru in z rotacijami povezanih kolicin se namrec odraza tudi na zvezi med odvodi rotacij po casu in kotnimi hitrostmi. Zveza med obema kolicinama je dolocena z enako transformacijsko matriko kot zveza med krajevnim odvodom rotacij in ukrivljenostjo. Sodobne metode casovne integracije bi zato morale temeljiti na cim bolj tocnem upoštevanju neaditivne narave z rotacijami povezanih kolicin, v nasprotju z uveljavljenimi metodami, ki lokalno uporabljajo aditivno diskretizacijo kolicin po casu.
Stevilne aplikacije, ki uporabljajo linijske elemente, segajo na področje analize vpliva drugih teles na konstrukcije. Poenostavljen model takega vpliva je masni delec. Reševanje povezanega problema gibanja nosilca in masnega delca ima dolgo zgodovino, predvsem v povezavi z analizo odziva mostov ob prehodu vozil. Velika vecina del s tega podrocja uporablja poenostavljene modele linijskih nosilcev [30]. Uspešno implementacijo geometrijsko nelinearne teorije pri analizi prehoda masnega delca v doktorski disertaciji predstavi E. Zupan [31]. Primer uporabe geometrijsko nelinearnega modela nosilca s predpostavljenimi majhnimi deformacijami v dinamicni analizi povezanega sistema podajnih teles je predstavil Bauchau
[32].	Predstavljeno casovno integracijsko shemo je nadgradil še za primer kontakta s trenjem in drsenjem
[33]	in za primer trka teles [34]. Lens in Čardona [35] sta predstavila zanimivo formulacijo geometrijsko nelinearnega modela nosilca in casovno integracijsko shemo primerno za analizo dinamicnih problemov, ki vsebujejo trk teles.
1.2 Vsebina dela
V doktorskem delu izpeljemo nova matematicna orodja za analizo prostorskih linijskih konstrukcij, ki so podvrzene dinamicnim vplivom. Izpeljave gradimo na osnovnih izrekih mehanike trdnih teles: izrek o gibalni kolicini in izrek o vrtilni kolicini. Drugi temelj izpeljav je upoštevanje tocnih kinematicnih enacb, ki podajajo zveze med pomiki in rotacijami ter njihovimi odvodi po parametru prostora in parametru casa.
Analiza prostorskih problemov zahteva skrbno obravnavo rotacij. Vse kolicine in izraze, povezane z rotacijami, ki jih potrebujemo v nalogi, izpeljemo in zberemo v posebnem poglavju o prostorskih rota-
cijah. Rotacije obravnavamo kot linearne preslikave v trirazseznem vektorskem prostoru, ki pripadajo multiplikativni grupi SO (3). Za parametrizacijo rotacijskega operatorja izberemo rotacijski vektor.
V	tematskem delu naloge se najprej ukvarjamo s casovno integracijskimi shemami. Izberemo algoritme, ki delujejo zgolj na kinematicnih zvezah med pomiki, rotacijami, hitrostmi in kotnimi hitrostmi ter pospeški in kotnimi pospeški. Najprej predstavimo Newmarkovo shemo [28] in njeno razlicico, prirejeno za prostorske rotacije [4]. V nadaljevanju razvijemo povsem novo casovno integracijsko shemo, ki temelji na drugacšni izbiri osnovnih neznank. Za osnovni neznanki cšasovne diskretizacije kinematicšnih zvez izberemo odvod vektorja pomikov po casu - vektor hitrosti - in vektorsko kolicino, ki je najbolj naravno povezana z odvodom rotacij po casu - vektor kotne hitrosti. Zaradi podobnosti z Newmarkovim algoritmom novo shemo poimenujemo modificirana Newmarkova shema. Delovanje obeh algoritmov predstavimo na primerih analize gibanja togih teles.
V	naslednjem poglavju predstavimo geometrijsko tocno teorijo za nelinearno analizo prostorskih nosilcev in podamo osnovne kinematicne zveze med pomiki in rotacijami ter njihovimi odvodi po kraju -deformacijami. Iz izrekov o gibalni in vrtilni kolicini izpeljemo osnovne enacbe dinamike trirazseznih linijskih elementov. Enacbe so zapisane tako, da omogocajo obravnavo nosilcev s poljubno zacetno geometrijo in zacetnimi deformacijami. Pri izbiri materialnega modela se ne omejimo na specificen material, predpostavimo le, da so konstitucijske zveze znane funkcije deformacij. Enakost ravnoteznih in konstitucijskih notranjih sil in momentov dosezemo z vpeljavo konsistencnih enacb.
V	nadaljevanju predstavimo numericno reševanje dinamicnih enacb geometrijsko tocne teorije prostorskih nosilcev. Izpeljemo dve formulaciji koncnih elementov, ki temeljita na izbiri deformacijskih kolicin kot osnovnih neznank problema. V prvem pristopu izhajamo iz objektivne in numericno ucinkovite formulacije za statiko prostorskih linijskih konstrukcij, ki je osnovana na izbiri deformacijskih vektorjev za osnovne neodvisne neznanke problema [9, 36]. Za casovno diskretizacijo izberemo Newmarkovo shemo, numericno stabilnost pa dosezemo z implementacijo posplošene-a sheme. V drugem pristopu za osnovne neznanke izberemo krajevne odvode hitrosti in kotnih hitrosti, ki so ekvivalentne casovnim odvodom deformacij. V novem pristopu za casovno integracijo uporabimo modificirano Newmarkovo shemo. Diskretizacija enacb po kraju v obeh formulacijah koncnih elementov temelji na interpolaciji osnovnih neznank: v prvi deformacij in v drugi krajevnih odvodov hitrosti. Za diskretizacijo zveznih enacb izberemo kolokacijsko metodo in jim v krepki obliki zadostimo le v izbranih diskretnih tockah. Diskretne enacbe rešujemo z Newtonovo metodo, zato jih lineariziramo, zaradi neaditivnosti prostorskih rotacij pa posebej predstavimo tudi dodajanje popravkov. Predstavljeni formulaciji koncni elementov vgradimo v racunalniški program zapisan v programskem okolju Matlab. Delovanje algoritmov testiramo na primerih znanih iz literature. Rezultate primerjamo z rezultati drugih avtorjev in predstavimo iz&pne numericne študije.
Studij dinamike prostorskega linijskega nosilca razširimo z reševanjem povezanega problema gibanja masnega delca glede na nosilec. Model nosilca razširimo z enacbami gibanja masnega delca na nosilcu. Posebej obravnavamo odziv sistema pri trku delca ob nosilec in predstavimo rešitve numericne študije.
2 OSNOVNE ENAČBE DINAMIKE 2.1 Osnovni izreki
Gibanje nosilcev, tako kot gibanje drugih trdnih teles, ki so podvrzšena dinamicšnim vplivom, vodijo osnovni izreki dinamike. Osnovni izreki so izhodišcna tocka naših izpeljav, zato jih v najsplošnejši obliki, ki velja za poljubno trdno telo, povzamemo v tem odseku.
Prvi pomemben izrek je izrek o ohranitvi mase, ki pove, da se masa telesa ali njegovega poljubnega dela med gibanjem ohranja, ce mase ne dodajamo ali odvzemamo:
pdV = p[0]dV[0].	(2.1)
Z p je oznacena gostota materiala z dV pa elementarni delec volumna v poljubni legi telesa. p[0] in dV[0] sta isti kolicini v zacetni legi telesa pri casu t = t0.
Naslednji od osnovnih izrekov je izrek o gibalni količini. Izrek o gibalni kolicini pravi, da je odvod gibalne kolicine telesa ali njegovega poljubnega dela po casu enak rezultanti vseh sil, ki delujejo na telo:
K = Frez.	(2.2)
Rezultanto sil, ki delujejo na telo ali na njegov poljubni del, smo oznacili z Frez. K oznacuje gibalno kolicšino, ki je definirana kot
K = j pr dV.	(2.3)
V
r je lega poljubnega delca telesa glede na opazovališce prostora O. S piko (') je oznacen odvod po casu.
Z opazovališcem prostora je povezan izrek o vrtilni količini, ki pravi, da je odvod vrtilne kolicine telesa ali njegovega poljubnega dela glede na opazovališce O, po casu enak rezultanti vseh momentov, ki delujejo na telo oziroma njegov poljubni del:
LO = M Oz.	(2.4)
Z Mrez smo oznacili rezultanto momentov glede na opazovališce O. Vrtilna kolicina telesa glede na opazovalisšcše O je definirana kot
LO = y r x pr dV.	(2.5)
V
Za uporabo osnovnih izrekov moramo primerno opisati poljubno lego telesa oziroma njegovega poljubnega delca. Osnovne izreke smo zapisali dokaj splošno. Za prakticno uporabo pa moramo izreke in kolicine, ki nastopajo v njih, zapisati v izbranih bazah prostora.
2.2 Kinematika delca
Opazujemo poljubni delec P telesa B (slika 2.1). V zacetni legi, pri casu t = t0, delec zavzema tocko, ki ima v izbrani parametrizaciji prostora koordinate (x,y, z). Lego delca P pri nekem casu t >t0 lahko natancno opišemo s krajevnim vektorjem delca r (x, y, z, t) glede na nepomicno opazovališce O. Za opis vektorjev glede na nepomicno izhodišce vpeljemo nepomicni sistem ortonormiranih baznih vektorjev {91,92,g3}, s katerimi je definiran globalni koordinatni sistem (X,Y, Z).
Slika 2.1: Lega poljubnega delca trdnega telesa. Figure 2.1: The position of an arbitrary particle of a solid body.
Med gibanjem se spreminja lega delca glede na opazovališce O. Hitrost delca v (t) je dolocena z odvodom krajevnega vektorja delca po casu. Krajevne vektorje delcev obicajno zapisujemo v nepomicni bazi, hitrost delca je potem enaka:
Vg (t) = rs (t).	(2.6)
Pospešek delca glede na izhodišcne nepomicne baze je dolocen z odvodom vektorja hitrosti po casu:
aä (t)= V g (t).	(2.7)
Telo se med gibanjem lahko deformira in vrti, pri tem pa se relativna lega delca P glede na ostale delce telesa spreminja. Te spremembe najlazje opišemo, ce na delec pripnemo pomicno ortonormirano bazo {G1 (x,y, z, t), G2 (x,y,z,t), G3 (x,y, z, t)}, ki napenja lokalni koordinatni sistem (£, n,Z). Deformiranje telesa opisujejo kinematicne enacbe, odvisne od parametrov trirazseznega prostora. Zapis teh enacb v splošni obliki presega okvire tega dela. Kinematicne enacbe bomo podali kasneje, ko bomo predstavili racšunski model nosilca.
Standardni opis gibanja teles tako vpelje vsaj dve bazi: nepomicno ali globalno bazo in pomicno ali lokalno bazo. Slednja je povezana z materialnim delcem telesa, zato jo pogosto imenujemo tudi materialna baza, vendar avtorji pri poimenovanjih niso enotni. Preslikava med obema ortonormirana bazama je rotacija R:
R : {91,92,93}-^{G1, G2,G3}.
Rotacije so pri obravnavi prostorskih problemov pomembne in, kot se izkaze, tudi zelo zahtevne mate-maticne strukture. V naslednjem poglavju pregledno predstavimo rotacije in z njimi povezane operatorje.
3 PROSTORSKE ROTACIJE
V tem poglavju predstavimo rotacije in operacije, povezane z rotacijami, ki jih potrebujemo v nadaljevanju dela. Poglavje o rotacijah je del mnogih ucbenikov, ki obravnavajo prostorske probleme v mehaniki, na primer Črisfield [37], v literaturi pa lahko poišcemo tudi številne samostojne publikacije, ki obravnavajo prostorske rotacije, omenimo samo deli avtorjev: Argyris [15] in Atluri in Čazzani [16]. Vse potrebne izpeljave apamo iz obširnega in natancnega pregleda podrocja prostorskih rotacij, ki ga je v doktorski disertaciji predstavil Zupan [36].
3.1 Rotacijska matrika
Rotacija je preslikava, ki poljubno ortonormirano bazo vektorskega prostora preslika v drugo ortonormi-rano bazo istega vektorskega prostora. Rotaciji R, izrazeni v bazi prostora, lahko priredimo rotacijsko matriko R. Zaradi narave rotacije ima rotacijska matrika lepe lastnosti:
(i)	rotacijska matrika R je ortogonalna matrika:
RRt = Rt R = I,	(3.1)
R-1 = Rt ,	(3.2)
kjer je I enotska matrika, ki pripada identicni preslikavi I: a —> a;
(ii)	ohranja dolzino vektorja, na katerega deluje, in
(iii)	ohranja kote med vektorji.
V nadaljevanju se osredotocimo na rotacijo, ki bazne vektorje nepomicne prostorske baze {g1,g2,g3} slika v bazne vektorje pomicne baze {G1, G2, G3}. Poljuben vektor pomicne baze je zavrten vektor nepomicne baze:
Gi,g = Rfl gt.	(3.3)
Z indeksom g smo poudarili, da so rotacijska matrika R in vektorji G j izrazeni v nepomicni bazi. (Če vektorje lokalne baze, izrazšene v globalni bazi, zapisšemo po komponentah:
G
1,9
Ga1,ig1 + Gg2,ig2 + Gg3,ig3, i = 1,2,3,
lahko preberemo komponente rotacijske matrike Rg:
Rg
Gg1,1 Gg1,2 Gg1,3 G g2,1 G g2,2 Gg2,3 Gg3,1 Gg3,2 Gg3,3
(3.4)
Stolpci rotacijske matrike, ki nepomicno bazo slika v pomicno, so torej kar komponente vektorjev lokalne baze, izrazene glede na globalno bazo. Zaradi ze navedenih lastnosti rotacije, (i)-(iii), so komponente rotacijske matrike medsebojno odvisne. Povezuje jih šest skalarnih enacb:
G1,g ■ G2,g = Rg91 ■ Rg#2 = 0 G2g ■ G3,g = Rg#2 ■ Rg93 = 0 ^3,g ■ ^1,g = Rgg3 ■ Rgg1 = 0 ||G1|| = ||Rg gxH = 1 IG2II = ||Rg 92^ = 1 ||G3N = l|Rg 93II = 1.
(3.5)
S piko (■) je oznacen skalarni produkt vektorjev, z (|| ||) pa dolzina vektorja. Vezne enacbe povedo, da so stolpci rotacijske matrike enotski in paroma ortogonalni.
Poljuben vektor u lahko zapišemo glede na lokalno ali glede na globalno bazo:
u = Ug = uG, Ug = Ug1g1 + ug2g 2 + u g3 g 3, UG = UG1G1 + UG2G2 + UG3G3.
Transformacijo komponent vektorja u med obema bazama predstavlja ravno rotacijska matrika R:
Ug = Rg ug.	(3.6)
Podobno velja za poljuben linearni operator A, ki mu v nepomicni bazi pripada matrika Ag, v pomicni pa matrika AG. Transformacijo med obema zapisoma lahko zapišemo kot
Ag = Rg agrT .	(3.7)
Pravilo (3.7) uporabimo za zapis rotacijske matrike v pomicni bazi:
RG = RJ Rg Rg = Rg = R.
Rotacijska matrika ima torej enak komponentni zapis v obeh bazah, med katerima slika, zato indeks baze pri zapisu rotacijske matrike v nadaljevanju opušcamo.
3.2 Kotna hitrost in kotni pospešek
Spreminjanje orientacije lokalne baze s casom dolocajo odvodi baznih vektorjev po casu. Za poljubni bazni vektor Gi velja:
G i,g = Rgi.
Obicajno zelimo spreminjanje lokalne baze izraziti glede na isto bazo.	Upoštevamo obrnjen izraz (3.3) in dobimo
Gi,g = RRTGi,g.	(3.8)
Produkt RRT je antisimetricna matrika, saj zaradi lastnosti (3.1) velja
4- (RRT) = RRT + RRT = d (I) = 0 dt dt
rrt = - rrt ^ (irrt )t .
Antisimetricni operator RRt oznacimo z Q in imenujemo operator kotne hitrosti. Ker na levi in desni strani enacbe (3.8) nastopajo vektorji lokalne baze, zapisani v globalni bazi, tudi produkt RRt doloca operator kotne hitrosti glede na globalno bazo:
Ofl = RRt .	(3.9)
Za vsak antisimetricni operator lahko poišcemo tak pripadajoci osni vektor, da velja:
n« = u x u.
Odvodi baznih vektorjev po casu so potem doloceni s predpisom
Gi,g = Ug x Gi,g.	(3.10)
Ug doloca hitrost vrtenja vektorjev pomicne baze, zato ga imenujemo vektor kotne hitrosti. Zapis operatorja kotne hitrosti v lokalni bazi dobimo s koordinatno transformacijo
= Rt &g R = Rt RRt R = RtR ,	(3.11)
pripadajoci osni vektor pa je enak: ug = RTUg. Z vektorjem kotne hitrosti ug lahko zapišemo odvode baznih vektorjev pomicne baze, izrazene kar v pomicni bazi:
Gi,G = RTGi,g = Rt (Ug x Gi,g) = RtUg x RtGi,g = UG x Gi,G.	(3.12)
Zanimajo nas tudi drugi odvodi baznih vektorjev po casu. Z odvajanjem izraza Gi,g = Og Gi g dobimo
G i,g = & g Gi,g + &g G i,g = & g Gi,g + &g &g Gi,g = g + Gi,g. (3.13)
Osni vektor antisimetricnega operatorja Čg imenujemo kotni pospešek in oznacimo z ag. Kotni pospešek ag je torej kar casovni odvod vektorja kotne hitrosti Ug, enako pa velja tudi za vektorja zapisana v pomicni bazi:
U g = dt (RT Ug) = RT Ug + Rt U g = RT Rüg + Rt U g = (Rt R )T UG + Rt U g = -&GUG + RT U g
= Rt U g = «g.	(3.14)
3.2.1 Kompozitum rotacij in seštevanje kotnih hitrosti
Vzemimo tri baze prostora {g1,g2,g3}, jg11], g21], g31] j in j g12], g22, g32 j, predstavljene na sliki 3.1. Preslikava med bazama {g1, g2, g3} in | g1 ], (2 ], (3 ] 1 je dolocena z rotacijo R[1]:
R[1]: {g1, g2, g3}-^ {g11], g21], g31]} ,
ki ji pripada rotacijska matrika R[1]. Rotacija R[1-2], dolocena s predpisom:
Rt1-2! : {(11],g21],g31]} {g12],g22],g32]} ,
Slika 3.1: Kompozitum rotacij. Figure 3.1: Čomposite of rotations.
slika med bazama -j g[1] , g21], G1] j in j g[2] , g22] , g32] 1-, pripada pa ji rotacijska matrika R[1 2]. Z
, ro-
tacijsko matriko R[1 2] lahko poljuben vektor baze | G1], g2], G3 ] j izrazimo v bazij G11], G21], G31] Gf1 = R[1-2]Gi1]. Bazni vektor G^1] pa lahko zapišemo kotzavrten vektor gi baze {g1,g2,g3}: G^1] = R[1]gj. (Če zdruzimo oba zapisa, dobimo predpis, ki daje zvezo med baznimi vektorji baze {g1,g2,g3}
in{ g12], g22], g32]
Gi2]= R[1-2]R[1]gi. Rotaciji R[2], ki slika med bazama {g1, g2, g3} in j g12], g22], g32]
r[2]: {91, g2, g3}-^{ g12], g22], g!2]
1 , G2 , G3
torej pripada rotacijska matrika R[2]
R[2] = r[1-2]r[1] .	(3.15)
Operacija, ki zdruzi dve rotaciji v eno samo, torej ni seštevanje temvec mnozenje rotacijskih operatorjev, kar poudarimo z izrazom kompozitum rotacij. Ker ima rotacijska matrika enak komponentni zapis v bazah, med katerima slika, tudi kompozitum rotacij v izrazu (3.15) lahko zapišemo glede na eno od obeh baz:
r[2] = R[1-2]R[1] = r[2] = R[1-2]R[1] Rg = Rg Rg = rg!2] = rg!2] rg!2] .
Pogosto poznamo rotacijsko matriko R! 2] izrazeno v bazi |g11] , G21], G31]}. Z uporabo pravila (3.7) kompozitum rotacij (3.15) dobi obliko
R[2] = R[1]R[1-2]R[1]T R[1] = R[1]R[1-2]	(3 16)
Rg = Rg rgw Rg Rg = Rg RGw .	(3.i6)
Zanima nas kotna hitrost kompozituma rotacij. Naj bo {g1,g2,g3}, nepomicna baza prostora bazi |g11], g21], g31]} in |g12], g22], g32]} pa pomicni bazi. Kotno hitrost kompozituma rotacij R[2] = R[1-2]R[1] dolocimo po definiciji (3.9). Za odvajanje po casu izberemo zapis kompozituma v nepomicni bazi, s cimer se izognemo problemu odvajanja baze, v kateri je rotacijski operator zapisan. Kotna hitrost
kompozituma je enaka
d (u[1-2W1A (D[1-2]D[1AT
Ol2] = i I2] Rp] T = -d ri1-2] ri1] ri1-2] rw ^g g g dt v g g y v g g
R[1-2]r[1] + R[1-2]li[1]> R[1]TR[1-2]T g g ^ g g J g g
= R [1-2] R[1-2]T + R[1-2] i [1]R[1]T R[1-2]T g g ^ g g g g
= oi1-2]+ rmOwri1-2^
g ^ g g g
Kotna hitrost kompozituma rotacij je torej vsota kotne hitrosti druge rotacije in transformirane prve rotacije. Vektor kotne hitrosti kompozituma rotacij je osni vektor gornje vsote, dolocen z
J2] = J1-2^ ri1-2] J1]
g g ' g g
Lepši in bolj zgovoren zapis vsote kotnih hitrosti dobimo, ce kotno hitrost kopmozituma zapišemo v novi bazi:
u[2l, = ri2] t J2]= ri2] t ^-2]+ ri1] t ri1-2] t ri1-2]
^G[2] g Ug g Ug + AXg AXg AXg Ug
= ^[1-2]+ ^[1] (317) = wg[2] + WG[1]-	(317)
Izraz (3.17) pove, da so kotne hitrosti aditivne kolicine samo v posebnem primeru, ko jih zapišemo v pomicnih bazah. Opozorimo, da sta sumanda v izrazu (3.17) izrazena v razlicnih lokalnih bazah. Rezultat je zelo pomemben pri numericnem reševanju enacb, v katerih kot osnovne neznanke nastopajo kotne hitrosti.
3.3 Odvod in integral vektorske funkcije v pomični bazi
Vzemimo poljuben vektor u (t) in ga zapišimo v pomicni bazi, odvisni od casa t. Absolutni ali totalni odvod vektorja v pomicni bazi pravilno dolocimo tako, da odvajamo komponente vektorja in bazne vektorje:
UG,abs (t) = UG1G1 + UG2G2 + UG3G3 + Ug1<G 1 + Ug2<^2 + UG3G3. Vektorju odvodov komponent pravimo relativni odvod
U G,re1 (t) = U G1G1 + U G2G2 + U G3G3,	(3.18)
vektorju, ki ga sestavljajo odvodi baznih vektorjev, pa sistemski odvod. Odvode baznih vektorjev izrazimo s kotno hitrostjo, kot to doloca pravilo (3.12). Sistemski odvod vektorja je potem enak
UG,sis (t) = UG1G1 + UG2G2 + UG3 G3
= UG1^G X G1 + UG2^G X G2 + UG3^G X G3 = Ug X (UG1G1 + UG2 X G2 + UG3 X G3)
= ug X uG.	(3.19)
Kar smo napisali, je pomembno za nadaljevanje našega dela. V izrazih nastopajo komponentni zapisi vektorjev v pomicni in nepomicni bazi. Poljuben vektor u zapišimo s komponentami v pomicni bazi
ug (t) = R (t) uG (t) in izrazimo njegov odvod po parametru t
d
Ug (t) = - (R (t) UG (t))
= R (t) UG (t)+ R (t) Ug (t) = R (t) Og (t) ug (t) + R (t) Ug (t) = R (t) (^g (t) x ug (t)) + R (t) Ug (t) = R (t) (^g (t) x ug (t) + Ug (t)).
V oklepaju prepoznamo absolutni odvod vektorja zapisanega v pomicni bazi. uG x uG je sistemski odvod, UG pa relativni odvod. Tako smo opredelili pomen oznake UG. S tako oznako bomo oznacevali odvod komponent vektorja, izrazenega v pomicni bazi, torej njegov relativni odvod
Ug = UG,rel.	(3.20)
Podobno, kot velja za odvajanje vektorja, izrazenega v pomicni bazi, velja tudi za njegovo integracijo. Vzemimo zopet vektor uG (t). Pri integraciji vektorja uG po parametru t moramo upoštevati, da so tudi bazni vektorji odvisni od t. Integriramo tako komponente, kot tudi bazne vektorje:
t t
J Ug (t) df = J [UG1 (t) G1 (t) + UG2 (t) G2 (t) + UG3 (f G3 (t)] dt 00
t t t = JUG1 (t) G1 (t) df + JUG2 (t) G2 (t) df + JUG3 (t) G3 (t) dt. 000 Taka integracija je lahko zelo zahtevna, ker obicajno poznamo le spreminjanje komponent vektorja v odvisnosti od t. Spreminjanje baznih vektorjev v odvisnosti od t pa obicajno poznamo le posredno, preko rotacijske matrike R (t). Vektorje zato integriramo v nepomicni bazi t t
Ug (t) dt = J R (t) Ug (t) df.
g
00
V posebnih primerih pa lahko integracijo vektorjev opravimo kar v pomicni bazi. Prvi tak primer integracije v pomicni bazi je integral absolutnega odvoda vektorja. Naj bo UG,abs (t) = dtuG (t); potem je integral /q UG,abs (t) df enak
t
J UG,abs (t) df = UG1 (t) G1 (t) + UG2 (t) G2 (t) + UG3 (t) G3 (t) = Ug (t) - Ug (0) , 0
kar preprosto dokazšemo z odvajanjem.
Drug poseben primer integracije v pomicni bazi je integral relativnega odvoda, pri katerem integriramo samo komponente vektorja glede na pomicno bazo prostora. Naj bo UG,rel (t) relativni odvod vektorja uG (t), definiran z izrazom (3.18). Če integriramo zgolj komponente relativnega odvoda dobimo:
t t t t J UG,rel (t) df = J UG1 (t) dt G1 + J UG2 (t) dt G2 + J UG3 (t) df G3 0 0 0 0
= UG1 (t) G1 + UG2 (t) G2 + UG3 (t) G3 - UG1 (0) G1 - UG2 (0) G2 - UG3 (0) G3 = Ug (t) - Ug (0) .	(3.21)
Integral komponent relativnega odvoda vektorja uG je torej do konstante natancno enak vektorju uG.
	0	-U3	U2
S (u) =	U2	0	-U1
	-U2	U1	0
3.4 Parametrizacija rotacij
Zaradi lastnosti (3.5) so komponente rotacijske matrike medsebojno odvisne. Pri reševanju enacb, v katerih nastopajo rotacije kot neznanke, moramo hkrati zadostiti še vezne enacbe. Temu se izognemo, ce za vsako rotacijo dolocimo tri neodvisne parametre, s katerimi je rotacija enolicno dolocena. Iz literature so znane razlicne parametrizacije, med pogosteje uporabljanimi pa je rotacijski vektor. Rotacijski vektor — je dolocen s smerjo n in velikostjo rotacije fi:
- = fin.	(3.22)
Z rotacijsko matriko sta povezana preko matricnega eksponenta:
1 2 1 3
R (-) = exp (S (-)) = I + S (-) + ^ S ($)2 + 3 S ($)3 +...	(3.23)
V izrazu (3.23) smo vpeljali oznako S (u) za antisimetricno matriko, ki pripada poljubnemu vektorju u:
-S (u)T.	(3.24)
Rotacijsko matriko, parametrizirano z rotacijskim vektorjem, lahko zapišemo v zakljuceni obliki z dobro znano Rodriguesovo formulo [15,36]
R ($) = I + S ($) + i-^ S2 ($).	(3.25)
Rotacijski vektor $ je lastni vektor matrike R, zato velja —g = R—G = —G.
V nadaljevanju bomo pogosto srecali izraze, v katerih nastopajo antisimetricne matrike. Izpeljave nekaterih enacb precej poenostavimo, ce vnaprej pripravimo nekaj zvez med antisimetricnimi matrikami in njihovimi osnimi vektorji. Naj bosta u in v poljubna vektorja, potem velja
(u ■ v) S (u) = -S (u) S (v) S (u)	(3.26)
S (u) S (v) - S (v) S (u) = S (S (u) v)	(3.27)
(u ■ v) S (u) - (u ■ u) S (v) = S (S2 (u) v).	(3.28)
Uporaba antisimetricšnega operatorja je prikladna tudi za zapis vektorskega produkta poljubnih vektorjev u in v, ki ga nadomestimo z matricnim produktom:
u x v = S (u) v = -S (v) u.	(3.29)
Lastnosti (3.26)-(3.29) preprosto dokazemo, ce vektorja u in v zapišemo v komponentni obliki, opravimo ustrezne operacije in primerjamo leve in desne strani dobljenih izrazov.
Za nekatere izpeljave je prirocno, ce definiramo še operator A, ki antisimetricni matriki, priredi njen osni vektor:
A (S (u))= u.	(3.30)
3.4.1 Rotacija s konstantno kotno hitrostjo
Oglejmo si primer pomicne baze, ki se glede na nepomicno vrti s konstantno kotno hitrostjo:
Og = Rt R = konst.
Za tak primer lahko preprosto dolocimo spreminjanje rotacijske matrike s casom. Rešujemo zacetni problem:
R (t)= R (t) Og, R (0)= R[0],	(3.31)
kjer je R[0] zacetna vrednost rotacije. Ker je kotna hitrost konstantna, ima enak komponentni zapis v vsaki legi pomicne baze {G1 (t), G2 (t),G3 (t)}, zato oznako specificne lege baze opušcamo: OG[0] = O g [t] = OG. Rešitev zacetnega problema (3.31) je enaka:
R (t)= R[0] exp (tOG).	(3.32)
Dokaz je preprost. Z odvajanjem rešitve (3.32) po casu dobimo zacetno diferencialno enacbo: d /
R (t) = dt (R[0]exp (tOc)
d	1 4n2 i 1 +3r»3 i , 1 +non
R[0] 4 (I + tOc +1 t2OG +113 OG + ■ ■ ■ +1 tn OG +.
dtv — 2f- 'G+si1OG+■■■+nr1 OG
= R[0] (Og + tOG + 12OG + ■ ■ ■ + tn-1OG +...
= R[0] (I + tOG + 1t2OG + ■ ■ ■ + (n—jrtn-1OG-1 +...) Og
= R[0] exp (tOc) Og = R (t) O g.
Vrednost predpostavljene rešitve v zacetni legi R (0) je kar enaka zacetni vrednosti R[0]:
R (0) = R[0] exp (0 OG) = R[0]I = R[0], s cimer je dokaz zakljucen.
Če za zapis eksponentne mape v rešitvi (3.32) uporabimo Rodriguesovo formulo (3.25) in upoštevamo, daje osni vektor antisimetricnega operatorja OG enak uG, lahko zapišemo:
R (t)= R[0]R (twc).	(3.33)
Rotacijska matrika R (t^G) slika med zacetno lego pomicnega koordinatnega sistema in poljubno lego istega koordinatnega sistema pri casu t. Ze iz zapisa matrike sledi, da jo lahko parametriziramo z rotacijskim vektorjem:
6g (t) = t^G	(3.34)
Za rotacijsko matriko, ki slika med pomicnima bazama, v nadaljevanju pogosto uporabimo izraz in-krementna rotacijska matrika, za rotacijski vektor, ki jo parametrizira, pa inkrementni rotacijski vektor. Podobno kot za kotno hitrost tudi za inkrementni rotacijski vektor velja, da ima enak komponentni zapis v vsaki legi pomicne baze: 0G[0] = 0G\t] = 0G.
3.4.2 Zveza med kotno hitrostjo in odvodom rotacijskega vektorja po Casu
Naj bo tf g rotacijski vektor, odvisen od casa t. Zanima nas zveza med odvodom rotacijskega vektorja tf g in vektorjem kotne hitrosti ug, ki pripada rotaciji, doloceni z tfg. Kotna hitrost je osni vektor produkta RRT, zato najprej zapišemo odvod rotacijske matrike po casu. Rotacijsko matriko, parametrizirano z rotacijskim vektorjem tfg, zapišemo z Rodriguesovo formulo (3.25), odvajamo in dobimo
R (tfg ) = S (tfg ) +	S ( tf g ) S (tfg ) + S (tfg ) S (tfg
+	t S (tfg) + t sin t - - cos ^ t S2 (tfg).	(3.35)
Dolzino rotacijskega vektorja t zapišemo s pomocjo skalarnega produkta: t = tf g ■ tf g. Odvod dolzine rotacijskega vektorja je potem enak
* = tfg ■tfg + tfg ■ *g = tfgA..	^
Izraz (3.36) vstavimo v (3.35), uposštevamo lastnost (3.26) in dobimo
RR (tfg) = sitt S (tfg) +	[S (tfg) S (tfg) + S (tfg) S (tfg
t cos t - sin t	( • )
"S (tfg) S h? g S (tfg)


t3
t sin t - 2 (1 - cos t)	( • ^ 2
--S (tfg) ^tfgJ S (tfg) .	(3.37)
Odvod rotacijske matrike po casu (3.37) z desne pomnozimo z RT (tfg). Po daljšem urejanju produkt izrazimo z vsoto antisimetricšnih matrik
R (tfg) RT (tfg) = S (tfg) + ^to^ S (S (tfg) tfg) + ^^ S (S2 (tfg) tfg
Kotno hitrost potem lahko zapišemo kot produkt matrike T (tfg), neodvisne od odvoda rotacijskega vektorja, in tf g:
^g = A (R (tfg ) RT (tfg )) = T (tfg ) tfg ,	(3.38)
kjer je matrika T (tfg) enaka
T(tf) = I +	S (tf) + ^S2 (tf).	(3.39)
Matriko T v literaturi pogosto imenujejo transformacijska matrika (npr. [36, 38]) in ji, podobno kot rotacijski matriki, pripada enak zapis tako v fiksni kot v ustrezni lokalni bazi: T (tfg) = T (tfG). Transformacijska matrika ima še nekaj lepih lastnosti:
Rt (tf) T (tf) = T (tf) Rt (tf)= Tt (tf)	(3.40)
Rt (tf) T (tf) R (tf)= T (tf),	(3.41)
njen inverz pa je enak:
T-1 (tf) =1 - 1S (tf) +-1 (1 - S2 (tf) ■
t2 \ 2 tan f
Če v enacbi (3.38) upoštevamo ug = R^g in lastnost (3.40) transformacijske matrike, dobimo zvezo med kotno hitrostjo, zapisano v pomicni bazi, in odvodom rotacijskega vektorja v nepomicni bazi:
= TT (d g ) dg .	(3.42)
Zgornje izpeljave lahko posplošimo na odvode rotacijskih matrik po poljubnem parametru, od katerega je matrika odvisna. Vzemimo na primer, da je inkrementni rotacijski vektor iz rešitve (3.34) zapisan v pomicni bazi prostora, odvisni od prostorskega parametra x. Zanima nas, kako je odvod rotacijske matrike R (6G (x)) po parametru x povezan z odvodom vektorja 6G (x) po x.
Odvod rotacijske matrike, izrazene z rotacijskim vektorjem v pomicni bazi, je odvod njenih komponent in ga formalno lahko zapišemo na enak nacin, kot smo zapisali odvod rotacijske matrike po casu (3.37):
d R (6 g (x)) = ssf S (6G) +	[s (6G) S (6 g) + S (6 g) S (6')]
ecose - sineS (6g) S (6G) S (6g)

e3
e sin e - 2 (1 - cos t

S (6g) s (6'g) s (6g)2
e4
Potem je produkt R' (6G (x)) RT (6G (x)) enak
R' (6g) RT (6g) = S (6G) +	s (S (6g) 6'g) +	S (S2 (6g) 6'g) ,
osni vektor tega produkta pa je enak
A (R' (6g) RT (6g)) = T (6g) 6'.	(3.43)
Na koncu še enkrat poudarimo, da smo odvajali samo komponente rotacijske matrike, zapisane v pomicni bazi. Zato je odvod inkrementnega rotacijskega vektorja 6G po x relativni odvod vektorja, ki ga v skladu z dogovorom (3.20) pišemo kot 6'G.
3.4.3 Variacija rotacijske matrike, parametrizirane z rotacijskim vektorjem
Vzemimo rotacijski vektor d in njegovo poljubno spremembo eöd. Smerni odvod rotacijskega operatorja R (d) je limita spremembe operatorja, povzrocene z eöd, ko gre e proti nic. Rotacija, ki pripada spremembi rotacijskega vektorja eöd, je enaka R (eöd), skupna rotacija pa je kompozitum R (eöd) R (d). Smerni odvod potem izracunamo kot:
ÖR w = i
R (eöd) R (d).
€=0
R (eö'd) zapišemo z Rodriguesovo formulo (3.25), odvajamo po e, izvrednotimo pri e = 0 in dobimo:
ÖR (d) = S (öd) R (d).	(3.44)
Podobno je sprememba rotacije lahko dolocena tudi kot komponiranje z desne R (dg) R (eödG), od koder sledi:
ÖR (d)= R (dg ) S (ÖdG) .
Variaciji obeh izrazov morata biti enaki, potem velja
S (50g) = RTS (Čtfg) R ^ 5-&G = RTčtfg.	(3.45)
Variacijo rotacijske matrike, ki deluje na poljuben vektor u, lahko izrazimo z variacijo rotacijskega vektorja:
ČR (tf) u = S (Čtf) R (tf) u = Čtf x R (tf) u = -R (tf) u x Čtf = -S (R (tf) u)	(3.46)
3.4.4 Variacija rotacijske matrike, parametrizirane z aditivnim rotacijskim vektorjem
Pri reševanju kinematicnih enacb s casovno diskretnimi shemami imamo rotacije pogosto izrazene s kompozitumom znane rotacije in inkrementne rotacije, ki je obicajno zapisana v pomicni bazi prostora. Primera takih rotacij smo zapisali ze v izrazih (3.16) in (3.33).
Vzemimo inkrementno rotacijsko matriko, parametrizirano z rotacijskim vektorjem, zapisanim v pomicni bazi: R (og). Variacijo rotacijske matrike lahko poišcemo multiplikativno
ČR (Og) = d
ae
R (Og) R (eč«g)
e=0
Spremembo rotacijske matrike zapišemo z Rodriguesovo formulo, odvajamo po e in izvrednotimo pri e = 0. Dobimo
ČR (og)= R («g) S (5og)	(3.47)
Spremenjeno rotacijsko matriko lahko opišemo tudi tako, da rotacijskemu vektorju spremembo kar ?G). Vektor Č«G
da velja
prištejemo: R (Og + eč«G). Vektor 50G imenujemo "aditivni" rotacijski vektor. DoloCmo ga tako,
ČR (°g) = d
e=0
R (Og + eȫG)
R (Og + eč«G) zapišemo z Rodriguesovo formulo (3.25), odvajamo po e in izvrednotimo pri e = 0. Dobimo:
ČR («G) = ^S (Č«G) +(S (Č»G) S (»g) + S (»g) S (Č»G))
+ 0cos00-sin0 («g • Č»G) S («G) + 0sin0 - 20(1 - cos0) («G. Č«G) S2 (»G) . (3.48)
Podobno kot za odvod rotacijske matrike lahko tudi za njeno variacijo v smislu aditivnega rotacijskega vektorja pokazemo, daje produkt ČRRT antisimetricna matrika
RRt = I ^ ČRRt + RČRt = 0 ^ ČRRt = -RČRt = - (ČRRt)t .
Izraz (3.48) z desne pomnozimo z RT (»g) in po koncanem urejanju zapišemo
ČR («g) RT (»g) = S (Č«G +	S («g) Č«G + ^^^S2 (Og) ČOG) .	(3.49)
V oklepaju na desni strani izraza (3.49) prepoznamo produkt transformacijske matrike T (0G) in aditivne variacije ödG:
ÖR (6g) RT (0g)= S (t (6g) SOG).	(3.50)
Če izraz (3.47) z desne mnozimo z RT (OG) in primerjamo z (3.50), pokazemo da je povezava med aditivno variacijo SOG in multiplikativno variacijo 56G enaka
ÖdG = RT (6g ) t (Og) SOG = TT (Og ) SOG.	(3.51)
Tudi v primeru variacije rotacijske matrike v smislu aditivnega rotacijskega vektorja je zaradi reševanja enacb koristno variacijo SR (OG) u, kjer je u poljuben vektor, izraziti z variacijo aditivnega rotacijskega vektorja. Izraz (3.47) najprej pomnozimo z u, nato pa upoštevamo (3.51) in dobimo:
SR (Og) u = R (Og) S (ÖOg) u = -R (Og) S (u) SO g
= -R (Og) S (u) TT (Og) SOG-	(3.52)
4 ČASOVNO INTEGRACIJSKE SHEME
Matematicni modeli konstrukcij in mehanizmov so pogosto opisani s togimi zacetnimi problemi oblike y' = f (t,y), y (t0) = y0. Eno najstarejših definicij togih diferencialnih enacb sta v precej pragmaticni obliki ze leta 1952 podala Čurtiss in Hirschfelder [39], ki pravita: "toge diferencialne enacbe so tiste diferencialne enacbe, pri katerih implicitne metode reševanja delujejo bolje, obicajno mnogo bolje kot eksplicitne metode". Najpreprostejša implicitna metoda je Eulerjeva metoda
yn+1 = yn + hf (tn+1,yn+1) , h = tn+1 - tn.	(4.1)
(Če shemo (4.1) uporabimo za diskretizacijo zvez med pomiki in hitrostmi (2.6) ter med hitrostmi in pospeški (2.7), dobimo posebno skupino metod za reševanje diferencialnih enacb drugega reda. Eno od takih metod je leta 1959 predlagal Newmark [28]. Vse od takrat so klasicna Newmarkova metoda in njene modifikacije široko uporabljene za diskretizacijo gibalnih enacb v dinamiki togih in podajnih teles. Lepa lastnost takih metod je, da ne posegajo neposredno v formulacijo enacšb dinamicšnega sistema, uporabimo jih le za diskretizacijo zvez med kinematicnimi kolicinami, ki nastopajo v teh enacbah. Tak pristop izberemo tudi za casovno diskretizacijo kinematicnih kolicin v formulaciji enacb dinamike prostorskih nosilcev, ki jo predstavimo v naslednjem poglavju.
Poleg togosti diferencialnih enacb so numericne tezave pri casovni integraciji povezane tudi z multiplika-tivno naravo prostorskih rotacij in razlicnimi parametrizacijami rotacij. Razširitev Newmarkove metode na prostorske rotacije sta prva predstavila Simo in Vu-Quoc [4], ki za osnovno neznanko problema poleg pomikov izbereta tudi rotacijski vektor.
V tem poglavju najprej predstavimo Newmarkovo shemo in njeno razširitev na prostorske rotacije, v kateri kot osnovne neznanke nastopajo pomiki in zasuki. V nadaljevanju pa predstavimo povsem novo implicitno casovno integracijsko shemo, v kateri za osnovne neznanke izberemo odvode pomikov in zasukov po cšasu.
Diskretizirali bomo kinematicne enacbe (2.6)-(2.7) in (3.11)-(3.14):
r g (t)= v g (t)	(4.2)
v g (t)= ag (t)	(4.3)
R (t)= R (t) S (UG (t))	(4.4)
UG (t)= ÖG (t) .	(4.5)
Privzemimo, da pomike, zasuke, hitrosti in pospeške poznamo pri casu tn. Za posamezno kolicino (•), izvrednoteno pri poljubnem casu tm, uporabljamo oznako: (• )[m]. Oznaka velja tudi za bazne vektorje
-,..-,,-,,	.	[m]	[m]^-,[m] . [m]yy[m] . [m] /-»[m]
casovno spremenljive lokalne baze, na primer: «Gm = ug1 G1 + UG2 G2 + UG3 G3 .
4.1 Newmarkova shema
Zveze med translacijskimi kolicinami lahko diskretiziramo po vzoru standardne aditivne sheme [28], ki predpostavlja zveze med kinematicnimi kolicinami pri casih tn in tn+1:
>+1]
rgn] + U
u,
= h vgn] + h2
ß )agn] + ßagn+1]
(4.6)
vgn+1] = vgn] + h
(1 - y ) agn] + y ajn+1]
h = tn+1 - tn je casovni korak sheme, ß G [0,1] in y G [0,1] sta parametra sheme.
Pri diskretizaciji zvez med rotacijskimi kolicšinami moramo uposštevati naravo multiplikativne grupe prostorskih rotacij, njenega tangentnega prostora in njune povezanosti preko matricnega eksponenta. Newmarkova shema za rotacijske kolicine je tako ustrezno prirejena [4]
R[n+1] = R[n] exp (S (0g)) = R[n]R(og)
d g = h wJ^L + h
G[n]
[n+1]	[n]
WG[n+1 = WG[n] + h
!_ ^ a[n] + ßa[n+1] 2 ßl aG[n] + ßaG[n+1
(4.7)
(1 - Y) «|jw+ YaG++11]
Inkrementni vektor Og, ki slika med bazama {g^, G2n], G3n]} in {G1n+1], G2n+1], G3n+1]}, ima enak komponentni zapis v obeh bazah, med katerima slika: Ogh = 0G[n+1].
Dinamicne enacbe, v katerih so kinematicne kolicine diskretizirane z implicitno shemo, rešujemo ite-racijsko, zato jih lineariziramo. Linearizacijo Newmarkove sheme (4.6)-(4.7) in dodajanje popravkov predstavimo v nadaljevanju.
,
4.1.1 Linearizacija
Rotacije parametriziramo z rotacijskim vektorjem. Za osnovni neznanki implicitne sheme (4.6)-(4.7) izberemo krajevni vektor rgn+1] in rotacijski vektor tfgn+1]. Shemo najprej prepišemo v obliki, kjer so preostale neznanke diskretizacije izrazene z osnovnimi neznankami in kolicinami, znanimi pri tn, nato pa enacbe variiramo v smereh črgn+1], čtfgn+1]. Newmarkova shema z izbranima osnovnima neznankama rgn+1] in tf gn+1] ima obliko:
u
= r[n+1] _ r[n]
Jn+1] = _7
ßh U, ^ 1 --
vgn]+ ^ 1 - 2/^ a[n]
v
[n]
a[n+1] = 1 Tv[n+1] ag = 7h lvö
R (Og) = RT ftfgn]) R(tfgn+1]
+ ( 1 - agn]

[n+1] G[n+1]
ßh Og + (1 - ß
[n] , u / 1 Y \ [n] ^GN+ %1 - Yß) "GH
a [n+1] = 1 fl,[n+1] _ . ,[n] aG[n+1] = 7h ^G[n+1] WG[n]
H 1 -1
Y
a
[n] G[n]
(x)
(4.8)
(4.9)
(4.10)
(4.11)
(4.12)
(4.13)
,
Variacije translacijskih kolicin sledijo neposredno iz sheme:
5ua = Črf+11	(4.14) =ß-ö%=(4.15)
«ag"+1| = = ^	(4.16)
Pri variacijah rotacijskih kolicin imamo nekaj vec dela. Variacijo öR yßgn+11 j ze poznamo iz enacbe (3.44):
ÖR (ßgra+11) = S () R (ßgra+11) ,	(4.17)
vendar v izrazih za kotno hitrost in kotni pospešek ne nastopa rotacijski vektor celotne rotacije ßg, temvec inkrementni rotacijski vektor 6G. Variacijo inkrementnega rotacijskega vektorja je zaradi zveze s kotno hitrostjo (4.12) bolj naravno izraziti v aditivnem smislu v skladu s poglavjem 3.4.4. Vzrok bo jasno viden iz dodajanja popravkov, osnovno idejo pa pojasnimo ze sedaj. Če <G++!| izrazimo v dveh zaporednih iteracijah in izracunamo razliko obeh vrednosti, sledi
"LGfc+n- wi,Gri = -ßh (0i+1'G - 0i'G).
Išcemo torej "aditivni popravek" inkrementnega rotacijskega vektorja, oznacimo ga z ödG, za katerega velja
6i+1,G = 0*,G + Ö&G.
Aditivni nastavek upoštevamo v linearizaciji enacbe (4.11) na levi strani, za desni del pa izraz (4.17)
R {6G + eÖdG) = RT (ßgn| ) S (öß[n+1] ) R (ß[n+1]
= R (6g) RT (-g^1 ) S (ößg^ ) R (ßg^1
d_ de
= R (6g) S (ößlG+j^ . Z upoštevanjem zveze med aditivno in multiplikativno variacijo (3.51) dobimo
Ö6G = T-T (6g) öß|ni++!|.	(4.18)
Variacijo kotne hitrosti in pospeška potem izrazimo kot
ö^l = ß-hö6o = ßh T-T (6g) öߣ++||	(4.19)
öaGn++1| = Y-Ö<Gn|++1| = 1 T-T (6g ) ÖßG++1|.	(4.20)
4.1.2 Dodajanje popravkov
Z rešitvijo lineariziranega sistema vodilnih enacb v trenutni, i—ti iteraciji dobimo popravke osnovnih neznank örgn+1 in ößg111. Vrednosti osnovnih neznank v novi iteraciji i + 1 dobimo tako, da popravke
dodamo starim vrednostim. Pomike seštejemo, rotacije pa mnozimo:
rin+g=rin+1]+s^	(4.21)
R (*£+.!!) = R (Ötfgra+10 R (tfT])	(4.22)
i+1,g J \ g J \ i,g
Popravljene vrednosti hitrosti, pospeškov, kotnih hitrosti in kotnih pospeškov lahko izrazimo tako, da enacbe (4.9)-(4.10) in (4.12)-(4.13) zapišemo v zaporednih iteracijskih korakih in jih odštejemo. Dobimo
vi+1,g = vi,g + ßh ig
a[n+1] = J™+1]+_^ Ör[n+1] ai+1,g = ai,g + ßh2 Öri,g
^i+1,G[n+1] = ^i,G[n+1] + ^T (Oi+1,G - Oi,G)
ßh
ai+1,Gl"+^1] = ai,Gin+1] + ßh2 (Oi+1,G - °i,G)
Vrednost totalnega in inkrementnega rotacijskega vektorja v novi iteraciji izracunamo iz popravljenih vrednosti rotacij s pomocjo Spurrierjevega algoritma [40]:
<+;; = Spurr (Rin+^1] Oi+1,G = Spurr (R[n]T Ri++1]
4.1.3 Izbira parametrov Newmarkove sheme
Izbira prostih parametrov sheme odlocilno vpliva na stabilnost in natancnost casovne integracije z Newmarkovo shemo. Analizo vpliva parametrov ß in 7 pri obravnavi linearnih dinamicnih sistemov podajajo številni avtorji, na primer [25,28,41]. Skupne ugotovitve avtorjev so, daje metoda brezpogojno stabilna pri izbiri parametrov 2ß > 7 > 1 in vkljucuje numericno dušenje pri vrednostih 7 > Nu-mericnemu dušenju sheme se izognemo pri vrednosti parametrov ß = 1 in 7 = 2, pri katerih je metoda drugega reda natancnosti - dobro znano trapezno pravilo. Pri analizi nelinearnih dinamicnih sistemov ni dokazov o brezpogojni stabilnosti sheme, ugotovitve o natancnosti, ki veljajo za linearne sisteme, pa lahko posplošimo tudi na nelinearne sisteme [42]. Pri razširitvi Newmarkove metode na nelinearen prostor rotacij [4] sta avtorja dokazala, daje pri izbiri parametrov ß = 4 in 7 = 2 metoda konvergentna in drugega reda natancnosti.
4.2 Modificirana Newmarkova shema
Reševanje enacb v nelinearni dinamiki konstrukcij ostaja odprt problem, kar se odraza v velikem številu objav novih metod in pristopov k casovni integraciji dinamicnih enacb. Temeljna zahteva oziroma usmeritev je razvoj stabilnih in robustnih shem, pri razvoju takih shem pa imamo vec moznosti. Eden od pristopov je razvoj integratorjev, ki stabilnost zagotovijo z algoritmicnim ohranjanjem ali disipacijo energije. Povsem drugacen pristop k povecanju stabilnosti in robustnosti algoritmov pa je drugacna izbira osnovnih neznank problema. Hosea in Shampine [43] pokazeta, da je izbira odvodov osnovnih
spremenljivk za neznanke iteracijskega procesa kljucna za uspešnost implicitnih metod za reševanje diferencialnih enacb. Primer takega pristopa sta pri obravnavi prostorskih nosilcev predstavila Zupan in Saje [9,17], ki namesto klasicne izbire osnovnih neznank - pomikov in zasukov - za osnovne neznanke formulacije koncnega elementa izbereta njihove krajevne odvode - deformacijske kolicine. S tako izbiro je avtomaticno rešen problem objektivnosti deformacijskih kolicin, blokiranja koncnih elementov, itd., ki so znacilni za koncne elemente, zasnovane na pomikih in rotacijah. Pri dinamiki vlogo krajevnih odvodov pomikov in zasukov prevzamejo njihovi casovni odvodi - hitrosti in kotne hitrosti. Zato novo shemo, ki jo predstavimo v nadaljevanju, osnujemo na izbiri casovnih odvodov pomikov in rotacij za osnovne neznanke problema.
Diskretizacijo kinematicnih enacb (4.2)-(4.5) zapišemo tako, da so hitrosti in kotne hitrosti edine neznanke diskretiziranega sistema. Nova diskretizacijska shema je podobna Newmarkovi, zato jo imenujemo kar modificirana Newmarkova shema.
Po Lagrangejevem izreku, glej Vidav [44, str. 183], za vsako odvedljivo funkcijo f (t) na koncnem intervalu [tn ,tn+i ] obstaja na tem intervalu vsaj ena tocka t, kjer je
f (t) = f (tn+l) — f (tn).	(4.23)
tn+1 tn
Z uporabo Lagrangejevega izreka (4.23) na poljubnem intervalu [tn,tn+1] lahko kinematicne enacbe (4.2) in (4.3) zapišemo v diskretni obliki:
,[n+1] _r[n]
g 'g
(4.24)
gh
[n+1] _ [n]
ag = vg h vg .	(4.25)
S h smo oznacili dolzino intervala: h = tn+1 — tn, vg je povprecna vrednost hitrosti, ag pa povprecna vrednost pospeška na [tn, tn+1]. Aproksimiramo ju z vrednostmi na robu intervala [tn, tn+1]:
i> g = (1 — ß) vgn]+ ß vgn+1]	(4.26)
ag = (1 — y ) agn]+ y agn+1],	(4.27)
kjer sta ß in y poljubna parametra z intervala [0,1].
Za osnovno neznanko diskretiziranega sistema (4.24)-(4.27) izberemo hitrost vg. Preostale neznane kinematicne kolicine na koncu intervala izrazimo z Vg:
rgn+1] = rgn]+ hV g.
vgn+1] = 1 vgn] + I Vg	(4.28)
a[n+1] = 7 — 1 a[n] __±_v[n] + _±_v
ag = y ag ß7h vg + ß7hvg.
Dobili smo novo implicitno diskretizacijsko shemo za translacijske kinematicne kolicine. Za razliko od klasicne diskretizacije (4.8)-sheme povprecna hitrost -vg.
[ n+ 1 ]
klasicne diskretizacije (4.8)-(4.10), pri kateri je osnovna neznanka rg , je osnovna neznanka nove
Pri diskretizaciji kinematicnih enacb (4.4) in (4.5) ravnamo podobno kot pri diskretizaciji translacijskih zvez. Povprecne vrednosti kotnih hitrosti in pospeškov na intervalu [tn, tn+1] izracunamo kot
& g = (i - ß) +ß &G+]	(4.29)
a G = (1 - y ) ag„]+ Y a[n,++1].	(4.30)
Smiselnost seštevanja kotnih hitrosti v pomicnih bazah smo utemeljili v poglavju 3.2.1 z enacbo (3.17).
Pri kinematicšno konsistentni aproksimaciji rotacij se zaradi multiplikativne narave rotacij ne moremo sklicevati na Lagrangejev izrek. Uposštevali bomo tocšno kinematicšno zvezo (4.4) ob predpostavki, da je kotna hitrost &G konstantna na intervalu [tn,tn+1]. Za tak primer smo rešitev kinematicne enacbe ze poiskali v poglavju 3.4.1. Ko rešitev (3.32) izvrednotimo pri casu tn+1, dobimo aproksimacijo rotacij izrazšeno kot
R[n+1] = R[n] exp (hS (c_g)) = R[n]R(hc_g) .	(4.31)
Zveza med kotnimi hitrostmi in pospeški je preprosta, kolicini pa sta aditivni, ce ju zapišemo v pomicni bazi, zato za aproksimacijo zveze (4.5) zopet uporabimo Lagrangejev izrek:
_ &G - &G
a G =-T-.
h
Za osnovno neznanko diskretizacije rotacijskih kolicin izberemo povprecno kotno hitrost c_G. Preostale neznane kolicine pri casu tn+1 pa izrazimo z znanimi vrednostmi pri casu tn in povprecno kotno hitrostjo c_ g:
R[n+1] = R[n]R (hc_G)
,[n+1] _ ß - 1 &[n] , 1
—&GW + ß
" [n+1] _ 2-1 a[n] ]--L&[n]+J_
Y GW hß y G[W] + hß y
&G++1 = &GW + * g	(4.32)
[n+1] Y- 1 [n]	[n] .
aG++1 = aGW - hß2&gW + — <_G.
Dobili smo novo implicitno diskretizacijsko shemo za rotacijske kolicine, ki temelji na povprecni vrednosti kotne hitrosti c_G.
Diskretizirane enacbe rešujemo z iteracijskim postopkom, zato jih lineariziramo. Linearizacija sheme in dodajanje popravkov sta opisana v naslednjem razdelku.
4.2.1 Linearizacija
Diskretizirane kinematicne zveze (4.28) in (4.32) lineariziramo pri casu tn+1 v smereh Čvg in to g. Linearizacija aditivnih zvez je preprosta:
čr!n+1] = höVg
ČvIn+1] = 1Sv g tojnW11] = 1 to g	(4.33)
g ß g ^G[n+1] ß äoT] = ßY,aG:W,l = G.
Linearizacijo diskretizirane rotacijske matrike (4.31) dobimo s smernim odvodom v smeri Scg. Linearizacijo pri casu tn+1 lahko zapišemo v aditivnem smislu kot
d
SR[n+1] = R[n]SR (hc G) = R[n] d
de
R (hcG + ehčc g) .
€=0
Tako linearizacijo smo pokazali ze v poglavju 3.4.4, zaradi aditivnosti kotnih hitrosti pa oznako A za aditivno variacijo opušcamo. V lineariziranih enacbah variacija 5R[n+1] vedno deluje na nekem vektorju u. Po enacbi (3.52) produkt SR (hcG) in poljubnega vektorja u zapišemo kot SR (hcG) u = —R (hcg) S (u) Tt (hcg) höüjG. Variacija celotne rotacije SR[n+1] na vektorju u pa je potem enaka
SR[n+1]u = — R[n]R (hcg) S (u) Tt (hüg) hčcg
= —R[n+1]S (u) Tt (hcg) hčcg.	(4.34)
4.2.2 Racun novih približkov
Z rešitvijo lineariziranega sistema vodilnih enacb dobimo iterativne popravke osnovnih neznank Svg in Scg. Nove iterativne rešitve osnovnih neznank dobimo tako, da popravke prištejemo trenutnim iterativnim vrednostim
V i+1,g = V i,g + SVg c i+1,G = c i,G + Sc G.
Nove iterativne vrednosti preostalih kinematicnih kolicin lahko izracunamo neposredno iz shem (4.28) in (4.32). Druga moznost pa je, da trenutnim iterativnim vrednostim prištejemo njihove popravke, ki jih dolocimo iz zvez (4.33):
r[n+1] ri+1,g	= r[n+1] i,g	+ hčv g				
[n+1]	[n+1] =v i,g	1 + ßßSv g	. .[n+1] i+1,G[I+++l]	,,[n+1] = "i&W	+	ß S * G
ai+1,g	= a [n+1] = ai,g	1 + a i,Svg ß^h	a[a+1]	a[:a+1] = aiGn+1]	+	ß YhS W G
Izjema so rotacije, ki jih ne moremo popravljati iterativno. Nove iterativne vrednosti rotacijske matrike v iteraciji i + 1 dobimo iz vrednosti povprecne kotne hitrosti v iteraciji i + 1:
Rn+1]= rmr (hci+1G) .
4.2.3 Izbira parametrov modificirane Newmarkove sheme
Preveriti zelimo vpliv izbire parametrov ß in y na odziv dinamicnega sistema, ki ga diskretiziramo z modificirano Newmarkovo shemo. Obravnavamo nedušeno nihanje sistema vzmet-utez. Gibalno enacbo sistema zapišemo v odvisnosti od mase utezi m in togosti vzmeti k
k
ag (t) + mrg (t)= 0.	(4.35)
Kineticna energija sistema je enaka Wkin (t) = 0mvg (t), potencialna energija pa WPot (t) = 1 kr2g (t). Opazujemo vpliv diskretizacijske sheme na ohranjanje celotne mehanske energije sistema.
Diskretizirana enacba gibanja (4.35) pri casu tn+1 ima obliko
hI aM__l_ vn + — v + a rN + ahv = 0.
Y g ßYh g ßYh g m g m g Če upoštevamo, da so pospeški pri casu enaki agn] = - m rgn], se diskretizirana gibalna enacba glasi
1 1 k k
vN + — Vfl + — rgn] + -hvfl = 0.	(4.36)
ßYh g ß7h mY g m
Z rešitvijo enacbe (4.36) dobimo rešitev za povprecno vrednost hitrosti vg na intervalu [tn,tn+1]. Ko poznamo resšitev, lahko zapisšemo spremembo celotne energije v dveh zaporednih cšasovnih korakih in opazujemo spremembo, ki nastane zaradi cšasovne diskretizacije. Sprememba celotne energije je enaka
aw = awfcin+awpot = w£+1] - +WN+11 - WplnJ
= 1 m (vf+1] 2 - vgn] 2) + 1 k )rgn+1]2 - rgn]2
2)+1K •
= 1 m )vgn+1] - vgn] ) (vgn+1] + vgn] ) + 1 k )rgn+1] - rgn] ) )rgn+1] + rgn] ) .	(4.37)
Spremembo energije zapišemo v odvisnosti od rešitve vg v skladu s casovno diskretizacijo (4.28)
AW = 1 m ( - ß vgn] + 1 v g) (-vgn] + 1 v^ + 2 k (hvg) )2rgn] + hv g ) .	(4.38)
rg J izrazimo iz (4.36) in vstavimo v enacbo (4.38). Po krajšem urejanju lahko spremembo energije zapisšemo kot
AW = m i-2* (-1 vgn] + 1 + k ^ (hv g )2
= m)v[™+1] - vg-])2 + k^ )rlr+1] -r^
2ß g g 2 g g 1 2ß
= —aWfcm + (1 - 2y ) aWpot.	(4.39)
V izrazu (4.39) je jasno izrazen vpliv parametrov ß in 7 na diskretizacijo enacb linearnega dinamicnega sistema z modificirano Newmarkovo metodo. Parameter ß vpliva na kineticno energijo sistema. Pri vrednostih ß < ^ bo kineticna energija rasla, pri ß > ^ pa se bo zmanjševala. Parameter 7 vpliva na potencialno energijo sistema. Če je 7 < 2, s casovno diskretizacijo pridobimo potencialno energijo, ce pa je 7 > 2, jo dušimo. Algoritmicni rasti ali disipaciji energije se lahko izognemo samo, ce za vrednosti parametrov ß in 7 vzamemo vrednost 1.
Ugotovitev o ohranjanju energije, ki veljajo za linearne dinamicne sisteme, ne moremo posplošiti na po-drocje nelinearne dinamike, ki ga obravnavamo v disertaciji. Za nelinearne primere nam ostane še dokaz konvergence in reda natancnosti metode. Modificirano Newmarkovo metodo (4.28), (4.32) zapišemo
tako, daizlocimo osnovni neznanki vg in üg, dobimo:
yag
vgn+1] = vgn] + h ((1 - 7) agn] + 7agn+1]) gn+1] = rgn] + hvgn] + h2ß )(1 - 7) agn] + 7agn+1])
r rg
üjn+1] = + h (1 - 7) ag] + h7ajn+1]
r[n+1] = R[n]R )hüj-] + h2ß (1 - 7) «g + h2ß7ajn+1]) .
Če za vrednosti parametrov modificirane sheme privzamemo ß = y = 2 in upoštevamo v zgornjem zapisu, ugotovimo, da se diskretizacija pomikov, rotacij, hitrosti in kotnih hitrosti z modificirano shemo ujema z diskretizacijo po osnovni Newmarkovi shemi z vrednostma parametrov ß = 1 in 7 = Zato lahko tudi za modificirano Newmarkovo shemo s parametroma ß = 7 = 1 trdimo, da ohranja pomembno lastnost shem osnovanih na Newmarkovi, to je daje konvergentna in drugega reda natancnosti.
4.3 Numericni primeri gibanja togega telesa
Gibanje togega telesa je natancno doloceno z lego tezišca in orientacijo poljubne baze, ki se giblje skupaj s telesom. Naj bo rg krajevni vektor tezišca togega telesa glede na nepomicno bazo {g1,g2,g3} in naj pomicnabaza {G1, G2, G3} sovpada z glavnimi vztrajnostnimi osmi togega telesa. Izreka o gibalni (2.2) in vrtilni kolicini (2.4) za togo telo lahko zapišemo v obliki
mag = F g	(4.40)
JaG + UG X JUG = RTMg.	(4.41)
m je masa togega telesa, J je vztrajnostni tenzor zapisan glede na glavne vztrajnostne osi telesa. Fg je rezultanta sil, Mg pa rezultanta momentov na tezišce telesa.
V nadaljevanju prikazemo reševanje gibalnih enacb togih teles z integracijskima shemama prikazanima v poglavjih 4.1 in 4.2. Pri izbiri prostih parametrov shem sledimo ugotovitvam iz tock 4.1.3 in 4.2.3 tega poglavja in v vseh numericnih primerih uporabimo parametre podane v preglednici 4.1.
Preglednica 4.1: Parametri casovno integracijskih shem uporabljeni v numericnih primerih. Table 4.1: Parameters of the time integration schemes used in numerical examples.
Newmarkova shema [4,28] Modificirana Newmarkova shema
ß=4	ß=2
11
Y = -	Y = ö
4.3.1 Vrtenje telesa pod vplivom predpisanega poteka momentov
Delovanje algoritmov zelimo najprej preizkusiti za primere, kjer poznamo analiticne rešitve. Obravnavamo primer, v katerem predpostavimo rešitev za rotacijski vektor $G (t), analiticno rešitev za kotno hitrost dobimo z odvajanjem predpostavljenega rotacijskega vektorja po casu, s pomocjo kinematicne enacbe (4.4), kotni pospešek pa iz enacbe (4.5). Analiticen potek momenta potem sledi kar iz gibalne enacbe (4.41): M g (t) = R (t) (JaG (t) + S (ug (t)) Jug (t)).
Obravnavamo dva primera, v katerih se rotacijski vektor s casom spreminja: (i) linearno: $G (t) = [100t, 20, — 2t]T in (ii) kvadraticno: $G (t) = [4t2,0,t]T. Vztrajnostni tenzor glede na glavne vztrajno-
stne osi telesa je enak:
10 0 0 0 5 0 0 0 1
Simulacijo gibanja izvršimo z novo shemo (4.32) in Newmarkovim algoritmom (4.7) na casovnem intervalu t G [0,20]. Numericne rešitve za rotacijske vektorje dobimo iz izvrednotene rotacijske matrike s pomočjo Spurrierovega algoritma [40]. Rešitve primerjamo s predpostavljenim potekom rotacijskega vektorja.
10-
|io-12
J io-14
ZA
V16
10-
,-10
|io-12
J io-14
ZA
V16
0 4
(a) Modificirani Newmarkov algoritem, h = 0.025
8 12 16 20
t
0 4
(b) Newmarkov algoritem, h = 0.025
8 12 16 20
t
Slika 4.1: Rotirajoce telo pod vplivom linearnih rotacij. Primerjava absolutnih napak komponent rotacijskega vektorja dobljenih z originalnim in modificiranim Newmarkovim algoritmom. Figure 4.1: Rotating body under linear rotations. Comparison of absolute error of components of the rotational vector obtained with original and modified Newmark algorithm.
10-
iS
^ 12
£10 1=
Jl0-14
ZA
V16
10-
10
0 4 8 12 16 20
t
(a) Modificirani Newmarkov algoritem, h = 0.1
		.....#G1
s M &10-	-12	#G2 —?G3
		
a |10-	-14	,yi
'o m	16	
0 4 8 12 16 20
t
(b) Modificirani Newmarkov algoritem, h = 0.05
Slika 4.2: Rotirajoce telo pod vplivom linearnih rotacij. Modificirani Newmarkov algoritem. Absolutna napaka komponent rotacijskega vektorja. Figure 4.2: Rotating body under linear rotations. Modified Newmark algorithm. Absolute error of
components of the rotational vector.
V prvem primeru s predpostavljenim linearnim potekom rotacij dajeta analizi z Newmarkovo in modificirano Newmarkovo shemo za casovni korak h = 0.025 podobne rezultate. Primerjava absolutnih
napak komponent rotacijskega vektorja je prikazana na sliki 4.1. Za casovne korake vecje od h = 0.025 je bil racun z integratorjem (4.7), osnovanim na rotacijah, neuspešen ze v prvem koraku, nova shema, osnovana na kotnih hitrostih, pa daje dobre rezultate tudi za vecje casovne korake. Absolutne napake komponent rotacijskega vektorja, izracunanega z modificirano Newmarkovo shemo, za casovna koraka h = 0.1 in h = 0.05 sta predstavljena na sliki 4.2.
Preglednica 4.2: Rotirajoce telo pod vplivom linearnih rotacij. Konvergenca Newtonove iteracije pri
t = 12.
Table 4.2: Rotating body under linear rotations. Čonvergence of Newton iteration at t = 12.
	Newmarkova shema		Modificirana Newmarkova shema	
i 1	p*|| 2.49976e+000	| /| 8.00074e+005	i 1	||/|| 5.11662e-013 7.42577e-010
2	8.67602e-002	2.78657e+004	2	5.21053e-016 9.17057e-013
3	9.90810e-004	3.17823e+002		
4	2.73014e-007	9.89206e-002		
5	8.77458e-015	2.88030e-009		
i - Newtonova iteracija, ||fe|| - norma vektorja prirastkov, ||/1| - norma vektorja desnih strani
Zanimiva je tudi primerjava hitrosti konvergence obeh metod in pripadajoci racunski casi. V preglednici 4.2 je prikazana konvergenca Newtonove iteracije za tipicen casovni korak, pri casu t = 12. Konvergenca nove sheme je obcutno hitrejša od klasicne Newmarkove sheme, temu primerno pa je tudi celoten racunski cas obcutno manjši. Pri racunu z Modificirano Newmarkovo shemo je znašal Tcelotni = 0.598 s, pri racunu s klasicno shemo pa Tcelotni = 1.482 s.
100
iS
^ 2
ft 10 1=
J 10-4
ZA
^ -6 10
100
s
^ 2
ft 10-2 t? 1=
I 10-4
^ -6 10
12
16 20
(a) Modificirani Newmarkov algoritem, h = 0.025
0 4 8 12 16
t
(b) Newmarkov algoritem, h = 0. 025
20
Slika 4.3: Rotirajoce telo pod vplivom kvadratnih rotacij. Primerjava absolutnih napak komponent rotacijskega vektorja dobljenih z originalnim in modificiranim Newmarkovim algoritmom. Figure 4.3: Rotating body under quadratic rotations. Comparison of absolute error of components of the rotational vector obtained with original and modified Newmark algorithm.
Absolutne napake komponent rotacijskega vektorja za drugi preucevani primer s kvadraticnim potekom
0
4
8
momentov so predstavljene na sliki 4.3. Rezultati, dobljeni z obema integratorjema, so prikazani za h = 0.025. Pri reševanju enacb z Newmarkovim integratorjem pri casu t = 15.75 Newtonova metoda ne konvergira, zato jeracun casovnega koraka neuspešen. Če za racun uporabimo casovne korake h < 0.01, pa metoda stabilno deluje do konca opazovanega intervala. Integracija z novo shemo je stabilna celo pri racunu z velikim korakom h = 0.05, enako pa velja ce uporabimo krajše casovne korake. Tudi za analizo s kvadraticnim potekom rotacij veljajo podobne ugotovitve kot pri analizi z linearnim potekom rotacij. Konvergenca nove sheme je obcutno hitrejša od klasicne. Newtonova iteracija za tipicen korak sheme, pri casu t = 12, je prikazana v preglednici 4.3.
Preglednica 4.3: Rotirajoce telo pod vplivom kvadratnih rotacij. Konvergenca Newtonove iteracije pri
t = 12.
Table 4.3: Rotating body under quadratic rotations. Čonvergence of Newton iteration at t = 12.
Newmarkova shema	Modificirana Newmarkova shema
i 1	Px|| 2.39656e+000	| /| 7.67051e+005	i 1	Px|| 1.33530e-001	| /| 2.14741e+002
2	8.89002e-002	2.67986e+004	2	1.38087e-005	2.64714e-002
3	9.90810e-004	3.17823e+002	3	2.06281e-013	3.25974e-010
4	2.80708e-007	2.95233e+002			
5	8.85612e-015	2.90506e-009			
i - Newtonova iteracija, ||čx|| - norma vektorja prirastkov, ||/1| - norma vektorja desnih strani
S primerom smo pokazali prednost nove casovno integracijske sheme v primerjavi s klasicnim pristopom pri obravnavi dinamicnih enacb v nelinearnem prostoru rotacij. V obeh obravnavanih primerih smo analizo z modificirano Newmarkovo shemo uspešno zakljucili z daljšimi casovnimi koraki kot analizo s klasicno Newmarkovo shemo za rotacije. Prednost nove sheme je v iteracijskem postopku, pri katerem so osnovne neznanke - kotne hitrosti - in njihovi iteracijski popravki elementi istega aditivnega Evklidskega prostora, medtem ko so osnovne neznanke klasicne sheme, osnovane na rotacijskem vektorju, elementi multiplikativne grupe SO (3), njihovi iteracijski popravki pa elementi linearnega tangentnega prostora.
4.3.2 TeZka vrtavka v gravitacijskem polju
V tem primeru preucujemo gibanje tezke vrtavke, ki je na enem koncu vrtljivo podprta in postavljena v gravitacijsko polje. To je standarden primer za testiranje casovno integracijskih shem za rotacije. Lokalni in globalni koordinatni sistem sta definirana, kot je prikazano na sliki 4.4. Lokalni koordinatni sistem je vzporeden z glavnimi vztrajnostnimi osmi vrtavke, pritrjen pa je v podprti tocki vrtavke. Lastnosti vrtavke in zacetne pogoje smo povzeli po referenci [26]. Teza vrtavke je mg = 20, kjer je g = 9.81 teznostni pospešek, vztrajnostni tenzor pa je podan glede na lokalni koordinatni sistem in je enak

5 0 0 050 0 0 1
Slika 4.4: Tezka vrtavka v gravitacijskem polju. Figure 4.4: The heavy symmetrical top in gravity field.
[0]	T '
Vrtavko zavrtimo z zacetno kotno hitrostjo uG] = [0,0,50] , nanjo pa ves cas gibanja deluje gravitacijska
T
T
sila Fg = [0,0, —mg] . Krajevni vektor tezišca vrtavke se s casom ne spreminja, ce ga zapišemo v
lokalni bazi, in je enak rG = [0,0,1]T. Vrtavka je v zacetni legi zavrtena okrog osi X za kot ^ = —0.3. Zacetna vrednost rotacijske matrike je tako enaka
R[0]
1
0
0
0 cos (0.3) — sin (0.3) 0 sin (0.3) cos (0.3)
Če gibalne enacbe vrtavke zapišemo glede na lokalni koordinatni sistem, pripet v podprti tocki O, je izreku o gibalni kolicšini identicšno zadosšcšeno, izrek o vrtilni kolicšini pa je oblike
JOaG + UG X JO UG = RT M O,
O,
»T 7
rO
O
kjer je M „ rezultirajoci zunanji moment na tocko O: M
O
= RrG X Fg.
Gibanje vrtavke opazujemo na intervalu t G [0,10]. Numericne simulacije smo izvedli z Newmarko-vim algoritmom za rotacije (4.7) in modificiranim Newmarkovim algoritmom (4.32). Oba algoritma pri enakih casovnih korakih dajeta zelo podobne rešitve. Slika 4.5 prikazuje rezultate za absolutno napako
10-
J3 10
a
§ 10-7 ca
JI 10-9
o
ZA
ü 10-11
0 2 4 6 8
t
(a) Modificirani Newmarkov algoritem
10
10
J3 10
13 a
10-7
ca
JI 10-9
o
ü 10-11
10-

0 2 4 6 8
t
(b) Newmarkov algoritem
Slika 4.5: Tezka vrtavka. Absolutna napaka celotne energije vrtavke. Figure 4.5: The heavy symmetrical top. Absolute error of the total energy od the top.
celotne energije vrtavke za oba algoritma. Po pricakovanju se z manjšanjem casovnega koraka zmanjšuje tudi absolutna napaka. V preglednici 4.4 podajamo še racunske case posameznih analiz, ki spet potrdijo dominantnost nove sheme. Ta je posledica, podobno kot v prejšnjem primeru, hitrejše konvergence Newtonove iteracije, dodaten prihranek pa je posledica dejstva, da so osnovne neznanke modificirane sheme aditivne kolicine, kar omogoca ucinkovito implementacijo metode. Pri Newmarkovi shemi za rotacije, kjer so osnovne neznanke neaditivne, se del racunskega casa porabi za racunanje inkrementnega in totalnega rotacijskega vektorja iz inkrementnih in totalnih rotacij.
Preglednica 4.4: Tezka vrtavka. Racunski casi. Table 4.4: The heavy symmetrical top. Computational time.
algoritem	- =1■10-2	- =5■10-3	- =1■10-3	- = 5■10-4
Modificirani Newmarkov alg.	1.32	2.63	10.48	20.18
Newmarkov algoritem	1.98	3.56	15.99	29.99
- = casovni korak
10'
10
G 10-4
10
10-
10-
10
-4
10-
h
10-
10-1
(a) Konvergenca druge norme vrtilne kolicine - primerjava Newmarkovega in modificiranega Newmarkovega algoritma z LIEMID algoritmi [45]
10°
10-
eRiO"
-4
10-
a
10-
10
10
10
h
10
10-1
(b) Konvergenca druge norme rotacijske matrike - primerjava Newmarkovega in modificiranega Newmarkovega algoritma z LIEMID algoritmi [45]
Slika 4.6: Tezka vrtavka. Konvergenca druge norme vrtilne kolicine in rotacijske matrike. Figure 4.6: The heavy symmetrical top. Čonvergence of the second norm of the angular momentum and
rotation matrix.
0
Rezultate analiz z Newmarkovim in Modificiranim Newmarkovim algoritmom smo primerjali še z rezultati, dobljenimi z druzino algoritmov LIEMID, ki jih je predstavil Krysl [45]. Algoritmi so osnovani na aproksimaciji rotacijskih gibalnih enacb togega telesa v sredinski tocki casovnega koraka ('midpo-int'algoritmi). Druzino algoritmov sestavlja implicitna metoda LIEMID[I], eksplicitna metoda LIE-MID[E], dve izboljšani eksplicitni metodi drugega reda LIEMID[E1] in LIEMID[E2] ter sestavljena eksplicitna metoda LIEMID[EA], ki pol casovnega koraka izracuna s prvo izboljšano eksplicitno metodo, pol pa z drugo. Avtor v clanku [45] prikaze konvergencne krivulje za velikost vrtilne kolicine in drugo normo rotacijske matrike pri casu t = 10. Ker analiticnih rešitev za ta problem ne poznamo, za referenco izberemo rešitev, dobljeno z zelo majhnim casovnim korakom h = 5 ■ 10-6, ki gaje za racun referencne rešitve uporabil tudi avtor v [45]. Na sliki 4.6(a) je prikazana primerjava krivulj konvergence za velikost vrtilne kolicine n, na sliki 4.6(b) pa je prikazana primerjava druge norme rotacijske matrike pri casu t = 10. Modificirana Newmarkova shema daje boljše rezultate za vrtilno kolicino, rezultati za rotacijsko matriko pa so malo slabši od tistih dobljenih z LIEMID algoritmi.
4.3.3 Dvojno nihalo
Obravnavamo problem dvojnega nihala, sestavljenega iz dveh nihal, dolzine 12 in mase m2 ter 11 in m1 (Slika 4.7). Dinamicen odziv dvojnega nihala je zelo obcutljiv na razmerja mas in dolzin nihal ter na zacetne pogoje. Obravnavamo samo ravninsko gibanje dvojnega nihala.
Z
Slika 4.7: Dvojno nihalo. Figure 4.7: Double pendulum.
Če za osnovne spremenljivke izberemo krajevne vektorje obeh mas r1,g in r2,g, lahko gibalne enacbe dvojnega nihala zapišemo iz izreka o ohranitvi gibalne kolicine (4.40) skupaj s kinematicnimi vezmi, ki zagotavljajo ohranjanje dozine nihala:
$
	0
	0
r1,g -l2 (r2,g - r1,g)T (r2,g - r1,g) - 12 Kinematicne vezi vpeljemo v sistem vodilnih enacb z uporabo metode Lagrangejevih mnoziteljev:
(4.42)
m1a1,g - F 1,g - $T10 A
T1,g' 'S
T
m2a2,g - F2,g - $T2,g Aö = 0
(4.43)
(4.44)
0
kjer sta F 1g, F2,g gravitacijski sili, Ag so Lagrangejevi mnoZitelji, $r1 g in $r2 g sta komponenti Jaco-bijeve matrike kinematicnih vezi:
0T
2 (r2,g - r 1,g)T _ '
Enacbe (4.42)-(4.44) predstavljajo sistem diferencialno algebrajskih enacb (DAE). Stabilnost numericnega reševanja DAE lahko izboljšamo [46], ce vezne enacbe rešujemo v šibki obliki:
$ = [	v1,g
. 2 (r2,g - r1,g)T (v2,g - v1,g)
Z drugacno izbiro osnovnih spremenljivk nihala lahko gibalne enacbe zapišemo tudi brez kinematicnih vezi (npr. [47, poglavje 5]):
(m1 + m2) + m212cos ($1 — $2) $2 + m212sin(^1 — $2)$2 + (m1 + m2) g sin $1 = 0
m211cos ($1 — $2) $1 + m212$2 — m211sin(^1 — $2)$1 + m2g sin $2 = 0. (4.46)
$1 in $2 sta kota, ki ga oklepata nihali z navpicnico. Vse kinematicne neznanke takega problema so aditivne kolicine, zato jih lahko rešujemo z uporabo algoritmov za reševanje zacetnih problemov v Ev-klidskih vektorskih prostorih, kot so metode druzine Runge-Kutta, BDF metode viskokega reda in metode z numericnimi diferencami [48]. V našem primeru uporabimo kar metode druzine Runge-Kutta za reševanje togih problemov, ki so vgrajene v programski paket Matlab. Te rešitve bomo uporabili za referencšne resšitve problema dvojnega nihala.
Obravnavamo tog primer dvojnega nihala z zelo razlicnima masama (slika 4.7), ki sta ga za demonstracijo delovanja energijsko konservativne sheme predstavila Lens in Cardona [49]. Lastnosti nihala in zacetne pogoje izberemo enake, kot so podani v clanku [49]: m1 = 0.005, m2 = 1, g = 9.81. Zacetni pogoji
TTTT
so enaki: r1,g = [1,0,0] , r2,g = [2,0,0] , v1,g = [0,0,0] , v2,g = [0,0,0] . Simulacijo dinamicnega odziva nihala racunamo z Newmarkovo in modificirano Newmarkovo metodo na intervalu t G [0,10].
Za izbrane parametre obeh shem (preglednica 4.1) se diskretizacija translacijskih kinematicnih kolicin z modificirano Newmarkovo shemo ujema z Newmarkovo shemo. Zato ni presenetljivo, da analize z obema shemama dajo enake rezultate. Trajektorije in hitrosti obeh mas ter energija sistema za casovni korak h = 10—3 so prikazani na sliki 4.8 in se do graficne natancnosti ujemajo z rezultati, predstavljenimi v [49].
Oscilacije hitrosti na sliki 4.8(b) kazejo, daje odziv nihala s tako razlicnima masama impulzivne narave. Za tocnejšo rešitev problema moramo izvršiti analizo z manjšim casovnim korakom. V ta namen problem, zapisan v obliki (4.46), rešimo z metodami, vgrajenimi v Matlabu, ki omogocajo kontrolo lokalne napake in prilagajajo casovni korak. Med metodami v Matlabu se je za reševanje ob strogih zahtevah za relativno in absolutno napako (tol = 10—10) najbolje obnesla metoda ode15s. Rezultati analize so prikazani na sliki 4.9. Trajektorije mas (slika 4.9(a)) in potek hitrosti (slika 4.9(b)) kazejo, da nihalo kmalu preide v izrazito impulzivno gibanje, za katero je znacilno vecanje amplitud hitrosti. Celotna energija sistema se vseeno ohranja (slika 4.9(c)).
Najmanjši casovni korak metode ode15s, potreben za izpolnitev predpisanih toleranc za napako, je bil h = 4.5272 ■ 10 6. Na podlagi tega simulacijo za povezan sistem enacb (4.43)-(4.45) izvedemo z modificiranim Newmarkovim algoritmom in Newmarkovim algoritmom z dovolj majhnim cšasovnim korakom
$
n,
2rL
_ — 2 (r2,g — r 1,g )
t
$
r2,,
	0
	0
(4.45)
Z
0.5 0
-0.5 -1 -1.5
50 40 30 20 ! 10 I 0 -10 -20 -30 -40
"-2 -1.5	-1 -0.5 0 0.5 1 1.5 2
	X
	(a) Trajektorijemas
- —	
- -	
- "-^2	
- 2	
	
0 1 2 3 4 5 6
t
(b) Hitrosti mas
8 9 10
20 15 10 5
^ 0
<D v
u -5 -10 -15 -20
/ — kinetična energija \ / \ / \
/ — potencialna energija \ / ■ / \
- / — totalna energija \ / \ / \ -\ i \ / \ / \ / \ / \
~7T~
~rr
! \
! \
! \
i
i
i
i
\ i \ i \ / \ /
\ /
i i i i \ i \ i \ i
\ I \ i
0 123456789 10
t
(c) Energija sistema
Slika 4.8: Dvojno nihalo. Rezultati dobljeni z Newmarkovo in modificirano Newmarkovo shemo, z
relativno velikim casovnim korakom h = 10-3. Figure 4.8: Double pendulum. Results obtained by Newmark scheme and modified Newmark scheme,
with relatively large time step h = 10-3.
7
h = 10-6. Pri racunu z Newmarkovo shemo ze v prvem koraku Newtonove iteracije opazimo slabo pogojenost Jacobijeve matrike lineariziranega sistema vodilnih enacb (4.43)-(4.45). Posledicno so dobljeni rezultati napacni, metoda pa ne konvergira. Racun z novim modificiranim algoritmom je stabilen ves cas analize. Rezultati pa so predstavljeni na sliki 4.10 in se odlicno ujemajo z referencnimi rezultati (slika 4.9).
Z
0.5 0
-0.5 -1 -1.5
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
X
(a) Trajektorije mas
50 40 30 20 ! 10 I 0 -10 -20 -30 -40
0123456789 10
t
(b) Hitrosti mas
20 15 10 5
^ 0
<D v
u -5 -10 -15 -20
' — kinetična energija
r\ i \ i \
potencialna energija \ / \
- ■ — totalna energija
n i \ i i
/ \ i
\ i
\ /
\.
/	\	/ \
/	\	' 1
. / \	\	! > / »
\ ' \ !	\	' 1 / 1
/ \ i / \ > i i
i

i i
'' / \ / \ /
\J
/ \ / \
i \ i \ 1 1
i i '
\ I \ ' 1 ' i ' i . s /
\ /
\ i
\ J
1
0123456789 10
t
(c) Energija sistema
Slika 4.9: Dvojno nihalo. Rezultati dobljeni z Matlabovo rutino ode15s. Figure 4.9: Double pendulum. Results obtained by the Matlab ode15s routine.
Numericne tezave pri reševanju slabo pogojenih sistemov diferencialno algebrajskih enacb v nekaterih primerih lahko odpravimo z razlicnimi numericnimi ukrepi s katerimi posegamo v sistem enacb. V literaturi, npr. [50,51], avtorji kot ucinkovite ukrepe predstavljajo metode predpogojevanja diskretiziranega sistema enacb in skaliranja Lagrangejevih mnoziteljev. Tak pristop uporabimo pri reševanju numericnih tezšav pri analizi dvojnega nihala z Newmarkovo shemo pri majhnih cšasovnih korakih. Lagrangejeve mnozitelje v enacbah (4.43)-(4.44) mnozimo s faktorjem s: Afl = sAfl. Dobimo sistem enacb
m1«1,0 - F 1,g - s$T10 Äö = 0
0.5 0
-0.5
Z
-1 -
-1.5 -
-2
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
X
(a) Trajektorijemas
50 40 30 20 ! 10 I 0 -10 -20 -30 -40
0 1 2 3 4 5 6
t
(b) Hitrosti mas
7 8 9 10
20 15 10 5
^ 0
<D v
u -5 -10 -15 -20

/\ /\ A	r\
■ — kinetična energija \ / \ , , ,■ — potencialna energija \ j \	\
- ■ — totalna energija
/
\ i \ / \
/N
' v / \ I \ I \
1	'	>	/	I	'	1	I	I	\	'
" 1 '	\ . \ I 1 \	\ I
\ /	\ i \ \ '	I '
- \	« \ ' ' > /	\ /'
V/	, , w , v , i/,	V
0123456789 10
t
(c) Energija sistema
Slika 4.10: Dvojno nihalo. Rezultati dobljeni z modificirano Newmarkovo shemo s kratkim casovnim
korakom h = 10 6.
Figure 4.10: Double pendulum. Results obtained with modified Newmark scheme with small time step
h = 10—6.
m2a2,g — F 2,g — S$T2,9 Xg = 0
2rlg v1,g = 0 2 (r2,g — r1,g)T (v2,g — v1,g) = 0,
ki ga rešujemo v odvisnosti od neznanih pomikov r1g, r2,g in skaliranih Lagrangejevih mnoziteljev Xg. Če za skalirni faktor izberemo vrednost s = h, je analiza z Newmarkovo shemo s korakom h = 10—6
stabilna ves opazovani interval, rezultati analize pa so enaki rezultatom, ki smo jih dobili z modificirano Newmarkovo shemo (slika 4.10).
Odziv nihala z majhnim korakom (h = 10-6) se precej razlikuje od rezultatov dobljenih z večjim korakom h = 10-3 . S tem primerom smo torej pokazali, kako pomembna je primernost implicitne cšasovno integracijske sheme za primere, ko so za zadostno natancnost potrebni majhni casovni koraki in z ukrepi za stabilizacijo numericne rutine ne zelimo ali pa ne moremo posegati v sistem enacb. Pomembna lastnost novo predstavljene sheme, osnovane na hitrostih, je torej povecana stabilnost numericne integracije pri majhnih cšasovnih korakih brez dodatnih numericšnih ukrepov.
5 OSNOVNE ENACBE DINAMIKE PROSTORSKIH NOSILCEV
Linijski nosilec je idealizirano trirazsezno deformabilno telo, pri katerem je ena dimenzija prevladujoca. Geometrijo nosilca vzdolz prevladujoce smeri lahko opišemo s prostorsko krivuljo, za opis preostalih dimenzij pa vpeljemo predpostavke. Matematicni model nosilca predstavimo v zacetku tega poglavja, pri cemer izhajamo iz standardnega modela v skladu z geometrijsko tocno teorijo nosilcev [1,2]. Kinematicne kolicine tako postanejo odvisne od dveh parametrov kraja in casa. Enacbe dinamicnega rav-notezja nosilca izpeljemo iz izrekov o gibalni in vrtilni kolicini, ob koncu pa jim dodamo še konstitucijske enacbe, s katerimi opišemo lastnosti materiala.
5.1 Matematični model prostorskega nosilca
Opazujemo nosilec s poljubno zacetno lego in obliko v trirazseznem Evklidskem prostoru. Nosilec lahko opišemo z druzino togih precnih prerezov [1,2]. Tezišca precnih prerezov so povezana s krivuljo, ki ji pravimo tezišcna os nosilca. Poljubna lega nosilca je natancno dolocena s krajevnim vektorjem tezišcne osi nosilca in z orientacijo precšnih prerezov. Lega nosilca v zacšetni in poljubni deformirani legi je prikazana na sliki 5.1. Za opis vektorjev in rotacij v prostoru sledimo oznakam, ki smo jih vpeljali v poglavju 2.2:
(i) Krajevni vektor tezišcne osi nosilca oznacimo z r (x,t). Merimo ga glede na nepomicno refe-rencno tocko O v prostoru. Pri tem je x G [0,L] naravni parameter tezišcne osi v zacetni konfiguraciji nosilca. Za opis gibanja vpeljemo nepomicni sistem ortonormiranih baznih vektorjev
(ii) Orientacija poljubnega precnega prereza je dolocena s smermi baznih vektorjev ortonormirane pomicnebaze {G1 (x,t), G2 (x,t), G3 (x,t)}, ki jo pripnemo v tezišce precnega prereza. Vektorja G2 in G3 sta orientirana vzdolz glavnih vztranostnih osi prerezov, G1 je normala: G1 = G2 x G3. Preslikavo med obema bazama doloca rotacijska matrika R (x, t), tako daje: Gj (x, t) = R (x, t) g^
5.2 Kinematične enačbe
Gibanje nosilca je odvisno od casovnega parametra t in prostorskega parametra x. Kinematicne enacbe podajajo zveze med pomiki in rotacijami ter njihovimi odvodi po parametrih. Odvode osnovnih kine-maticnih spremenljivk po casu smo zapisali ze v kinematicnih enacbah, ki veljajo za delec in poljubno pomicno bazo pripeto na delec (2.6)-(2.7), (3.11) in (3.14). Te enacbe veljajo tudi za vsak delec na
Slika 5.1: Matematicni model nosilca. Figure 5.1: Mathematical model of the beam.
tezšisšcšni osi nosilca
Vg (X,t)= Vg (X,t)	(5.1)
V g (x,t) = a,g (x,t)	(5.2)
R (x,t) = R (x,t) S (wG (x,t))	(5.3)
üg (x,t) = aG (x,t).	(5.4)
Ponovimo: vg (x,t) in ag (x,t) sta hitrost in pospešek delca na tezišcni osi, dolocenega s parametrom x. üG (x,t) in aG (x,t) sta kotna hitrost in kotni pospešek pomicne baze, pripete na prerez, dolocen s tocko x na tezišcni osi.
Kinematicnim enacbam (5.1)-(5.4) pripadajo zacetni pogoji:
Vg (x,0)= rg0] (x)	(5.5)
Vg (x,0)= vg0 (x)	(5.6)
R(x,0)= R[0] (x)	(5.7)
WG (x,0)= J§ (x).	(5.8)
Geometrijsko tocna teorija nosilcev za opis deformiranja tezišcne osi nosilca vpelje deformacijski vektor y (x,t). Spreminjanje orientacije lokalne baze vzdolz osi pa opisuje z deformacijskim vektorjem k (x,t). Zveze med kinematicimi in deformacijskimi kolicinami so izpeljane tako, da zadošcajo izreku
o virtualnem delu. Splošna oblika teh enacb se glasi
rg (x, t) = R (x, t) (yg (x, t) - cg (x))	(5.9)
R' (x, t) = R (x, t) S (kg (x, t) - dG (x)).	(5.10)
S crtico (') oznacimo odvod po parametru prostora x. Enacbe (5.9)-(5.10) so tocne kinematicne enacbe prostorskih nosilcev. Utemeljitelja takega opisa kinematike nosilcev sta Reissner [1] in Simo [2]. Konstanti cg in dG sta odvisni od zacetne oblike nosilca.
cg (x) = -Rt (x, 0) rg (x, 0) S (dG (x)) = -Rt (x, 0) R' (x, 0).
Nosilec obravnavamo kot element konstrukcije, zato je njegovo deformiranje odvisno od umestitve v konstrukcijo. Element se na robovih lahko stika z vec drugimi elementi, robovi nosilca so lahko podprti ali prosti, lahko pa jim pomike in zasuke tudi predpišemo. Kinematicnim enacbam (5.9)-(5.10) tako pripadajo še bistveni robni pogoji, ki jih formalno zapišemo kot
rg (0,t) = rg (t)	(5.11)
R (0,t)= R0 (t)	(5.12)
rg (L,t) = rL (t)	(5.13)
R (L, t) = Rl (t).	(5.14)
Za parametrizacijo rotacij izberemo rotacijski vektor (3.22). (Časovno kinematicno enacbo za rotacije, parametrizirane z rotacijskim vektorjem, smo izpeljali ze v poglavju 3.4.2. Odvod rotacijskega vektorja po casu izrazimo iz enacbe (3.42):
tfg (x, t) = T-T (tfg (x, t)) ^G (x, t).	(5.15)
Povsem analogno lahko izpeljemo izraz za odvod rotacijskega vektorja po prostorskem parametru. Ki-nematicna enacba za rotacije, parametrizirane z rotacijskim vektorjem, ima obliko:
tfg (x, t) = T-T (tfg (x, t)) (kg (x, t) - dG (x)),	(5.16)
pripadajoci robni pogoji pa so enaki
tfg (0,t) = < (t)	(5.17)
tfg (L, t) = < (t).	(5.18)
5.3 Gibalne enačbe
Z uposštevanjem predpostavk matematicšnega modela nosilca lahko gibalno (2.3) in vrtilno (2.5) kolicšino nosilca ali njegovega poljubnega dela 0 < x1 < x2 < L glede na nepomicno izhodišce O zapišemo kot:
X2
K g = j pA (x) vg (x, t) dx	(5.19)
X1
X2
Lg = / [rg (x, t) x pA (x) Vg (x,t)+ R (x,t) J (x) RT (x,t) <^g (x,t)] dx.	(5.20)
X1
Produkt gostote p in površine precnega prereza A oznacuje maso na enoto dolzine nosilca. J je mehanska vztrajnostna matrika, zapisana glede na glavne vztrajnostne osi prereza {G1 (x, t), G2 (x, t), G3 (x, t)}. Za zapis osnovnih izrekov potrebujemo sše odvoda gibalne in vrtilne kolicšine po cšasu:
X2
Kg = j pA (x) ag (x,t) dx
X1
X2
(5.21)
* O C
Lg = [rg (x,t) x pA (x) ag (x,t) + R (x,t) J (x) RT (x,t) ag (x,t)
X1
+ (x, t) x R (x, t) J (x) Rt (x, t) (x, t) dx.
(5.22)
Nosilec je obtezen z zunanjimi linijskimi silami n (x,t) in momenti m (x, t). Na robovih nosilca delujejo tockovne sile S0 (t), SL (t) in momenti P0 (t), PL (t). Na elementaren del nosilca poleg zunanjih linijskih sil delujeta še rezultanti notranjih sil N (x, t) in momentov M (x, t) (slika 5.2). Gibalne enacbe
p 0
Slika 5.2: Obtezšba nosilca. Figure 5.2: Loading of the beam.
nosilca izpeljemo iz izrekov o gibalni (2.2) in vrtilni (2.4) kolicini, ki ju zapišemo za infinitezimalen del nosilca dolzine dx, glede na izhodišce O. Izrek o gibalni kolicini dobi obliko
x+dx
x+dx
J pA (X) ag (X,t) dX = — Ng (x,t) + J ng (X,t) dX + Ng (x + dx,t).
Izrek o vrtilni kolicini, zapisan za elementarni del nosilca, pa dobi obliko
x+dx
J [rg (X, t) x pA (X) ag (X, t) + R (X, t) J (X) RT (X, t) ag (X, t)
x
+ (X, t) x R (X, t) J (X) Rt (X, t) (X, t)] dX = —Mg (x, t)
x+dx	x+dx
— rg (x,t) x Ng (x,t)+ J mg (X,t) dX + J rg (X,t) x ng (X,t) dX
xx
+ Mg (x + dx, t) + rg (x + dx, t) x Ng (x + dx, t).
V limiti, ko gre dx proti nic, dobimo enacbe dinamicnega ravnotezja linijskega nosilca:
Ng (x,t) = pAag (x,t) — ng (x,t)	(5.23) Mg (x, t) = R (x, t) JRt (x, t) ag (x, t) + Ug (x, t) x R (x, t) JRT (x, t) Ug (x, t)
— rg (x, t) X Ng (x, t) — mg (x,t).	(5.24)
Ravnotezni pogoji morajo biti izpolnjeni tudi na robovih nosilca. Pripadajoci naravni robni pogoji so tako enaki:
Ng (0,t) = —S0 (t)	(5.25)
Mg (0,t) = —P0 (t)	(5.26)
Ng (L,t) = SL (t)	(5.27)
Mg (L,t) = PL (t).	(5.28)
5.4 Konstitutivne enacbe
Deformacijsko-napetostno stanje v nosilcu je odvisno od lastnosti snovi, iz katere je zgrajen nosilec. Zveze med rezultnantnimi silami in momenti ter deformacijami v poljubnem prerezu nosilca so podane s konstitutivnimi enacbami:
N G (x, t) = CN (y g (x, t), kg (x, t))	(5.29)
M G (x, t) = Cm (y g (x, t), kg (x, t))	(5.30)
V izpeljavi enacb se ne omejimo na specificen material. Dovolj je, da privzamemo, da sta CN in CM znana in vsaj enkrat odvedljiva operatorja.
5.5 Enacbe konsistence
Napetosti, dolocene iz ravnoteznih enacb, so teoreticno enake napetostim, ki ustrezajo konstitutivnim enacbam. V zapisanih enacbah tega dejstva še nismo upoštevali. Sledimo konceptu, ki sta ga predstavila Vratanar in Saje [52], in eksplicitno zahtevamo enakost ravnoteznih in konstitucijskih sil in momentov v prerezih nosilca
Ng (x, t) = R (x, t) NG (x, t)	(5.31)
Mg (x, t) = R (x, t) M G (x, t)	(5.32)
Velika vecina formulacij koncnih elementov te enakosti ne zagotavlja.
5.6 Povzetek enacb
Zapisali smo vse enacšbe, s katerimi je problem dinamike prostorskega nosilca natancšno definiran. Povzetek enacb je skupaj s pripadajocimi zacetnimi in robnimi pogoji zapisan v preglednici 5.1.
Preglednica 5.1: Pregled enacb za dinamiko prostorskega nosilca v skladu z geometrijsko tocno teorijo
nosilcev.
Table 5.1: Summary of equations for dynamics of three-dimensional beam according to geometrically
exact beam theory.
kinematicne enacbe po casu
r g (x,t) = v g (x,t) V g (x,t) = ag (x, t)
R (x, t) = R (x, t) S (wG (x, t)) üG (x, t) = aG (x, t)
zacetni pogoji
rg (x, 0)= rg0](x)
vg (x, 0)= vg0] (x)
R (x, 0)= R[0](x) (x, 0) = (x)
kinematicne enacbe po kraju
rg (x, t) = R (x, t) (yg (x, t) - CG (x)) R' (x, t) = R (x, t) S (kg (x, t) - dG (x))
robni pogoji
rg (0,t) = r0 (t)
R (0, t) = R0 (t) R (L, t) = Rl (t)
rg (L, t) = rL (t)
ravnotezne enacbe
N g (x, t) = pAag (x, t) — ng (x, t) Mg (x, t) = R (x, t) JRt (x, t) ag (x, t)
(x, t) x R (x, t) JRt (x, t) ug (x, t) —r' (x, t) x Ng (x, t) — mg (x, t)
robni pogoji
Ng (0,t) = —S0 (t) Ng (L, t) = SL (t) Mg (0,t) = —P0 (t) Mg (L, t) = PL (t)
konsistencne enacbe
Ng (x, t) — R (x, t) N G (x, t) = 0 Mg (x, t) — R (x, t) M G (x, t) = 0
konstitutivne enacbe
N G (x, t) = CN (y g (x, t), kg (x, t)) M G (x, t) = Cm (y g (x, t), kg (x, t))
V	skupno 12 vektorskih enacbah (6 kinematicnih, 2 ravnoteznih, 2 konstitucijskih in 2 konsistencnih) nastopa 12 neznanih vektorskih kolicin: rg (x,t), R (x,t), 7G (x,t), kg (x,t), vg (x,t), ag (x,t), (x,t), aG (x,t), Ng (x,t), Mg (x,t), NG (x,t) in M G (x,t). Na tem mestu se moramo odlociti, katere od neznank bomo izbrali za osnovne neznanke in katere bodo tiste, ki jih bomo izrazili z osnovnimi neznankami. Izbira je naceloma poljubna, zavedamo pa se, da so vse neznanke zvezne kolicine in jih bomo v postopku numericnega reševanja enacb iskali v diskretni obliki. To pa pomeni, da jih bomo interpolirali. Z interpolacijo nekaterih kolicin so lahko povezane tudi precejšnje numericne tezave pri reševanju enacb.
V	delih Zupana in Sajeta [9,17] so za staticno analizo prostorskih nosilcev za osnovne neznanke izbrane deformacijske kolicšine. Taka izbira se je izkazala za zelo uspesšno, zato ji sledimo tudi pri resševanju dinamicnih enacb prostorskih nosilcev.
5.7 Zveze med deformacijami in hitrostmi
Deformacije so krajevni odvodi pomikov in rotacij, hitrosti pa so casovni odvodi istih kolicin. Zanima nas zveza med deformacijami in hitrostmi.
Ponovno si oglejmo osnovne kinematicšne enacšbe.
rg (x, t) = 7g (x, t)	(5.33)
7^g (x, t) = Vg (x, t)	(5.34)
R' (x, t) = R (x, t) S (Kg (x,t))	(5.35)
R (x, t) = R (x, t) S (wg (x, t)).	(5.36)
Zaradi preglednejšega zapisa smo v enacbi (5.33) deformacijski vektor 7 zapisali v nepomicni bazi skupaj s konstanto c: 7g = R(7G — cg), v enacbi (5.35) pa smo vpeljali oznako: KG = kg — dG.
Oglejmo si najprej zvezo med hitrostjo vg in deformacijskim vektorjem Kg. Enacbo (5.33) odvajamo po casu, enacbo (5.34) pa po kraju:
rg (x,t) = 7g (x,t) rg (x, t) = vg (x, t).
Primerjava obeh rezultatov pove, da je hitrost deformacije 7g kar enaka krajevnemu odvodu vektorja hitrosti vg:
7 g (x, t) = vg (x, t).	(5.37)
Na podoben nacin izpeljemo zvezo med kotno hitrostjo ^g in deformacijskim vektorjem kg. Enacbo (5.35) odvajamo po casu, enacbo (5.36) pa po kraju in dobimo
R' (x, t) = R (x, t) S (Kg (x, t)) + R (x, t) S Kg (x, t)
= R (x, t) S (^G (x, t)) S (Kg (x, t)) + R (x, t) S ^Kg (x, t)
R' (x,t) = R' (x, t) S (wG (x,t)) + R(x,t) S (w'G (x, t))
= R (x, t) S (Kg (x, t)) S (wG (x, t)) + R (x, t) S (w'G (x, t)).
S primerjavo obeh izrazov dobimo
S (Wg) S (Kg) + S (Kg) = S (Kg) S (Wg) + S (<^G) S (Wg) S (Kg) - S (Kg) S (Wg) = S (<^G) - S (Kg
ČČe upoštevamo še lastnosti antisimetricnih matrik, lahko zapišemo
(x,t) x Kg (x,t) = JG (x,t) - Kg (x,t).	(5.38)
Krajevni odvod kotne hitrosti je torej v tesni povezavi s hitrostjo deformacije KG. Izraz lahko izpeljemo tudi za nepomicno bazo. V fiksni bazi velja
R' (x, t) = S (Kg (x, t)) R (x, t) R (x, t) = S (wg (x, t)) R (x, t),
mešana odvoda pa sta enaka
R' (x, t) = S (kg (x, t)) R (x, t)+ S (Kg (x, t)) R (x, t)
= S ^kg (x,t)j R (x,t) + S (kg (x,t)) S (üg (x,t)) R (x,t) R' (x, t) = S (wg (x, t)) R (x, t) + S (üg (x, t)) R' (x, t)
= S (üg (x, t)) R (x, t) + S (üg (x, t)) S (kg (x, t)) R (x, t).
S primerjavo obeh rezultatov ugotovimo
S ^ /kg (x, t)^ + S (kg (x, t)) S (üg (x, t)) = S (üg (x, t) ) + S (üg (x,t)) S (kg (x, t)) S (kg (x, t)) S (üg (x, t)) - S (üg (x, t)) S (kg (x, t)) = S (üg (x, t)) - S ^kg (x, t) ^ .
kg (x, t) X üg (x, t) = ü'g (x, t) — kg (x, t) .
Rezultata (5.37) in (5.38) nas napeljeta k razmišljanju, da pojem deformacijske kolicine v dinamiki nekoliko razširimo. K deformacijskim kolicinam štejemo tudi kolicine, ki merijo hitrost deformiranja: Y G, k G, vg in ü'g.
6 NUMERICNO REŠEVANJE ENACB
V preglednici 5.1 smo zapisali sistem enacb, ki vodi gibanje prostorskega nosilca. Rešitev enacb zelimo poiskati v krepki obliki brez dodatnih predpostavk in veznih enacb, ki so znacilne za nekatere formulacije prostorskih nosilcev (npr. šibka oblika kinematicnih enacb [32], šibka oblika ravnoteznih enacb [10, 27,53], vezne enacbe za ortogonalnost rotacij [54]). Časovno diskretizacijski shemi, predstavljeni v poglavjih 4.1 in 4.2, sta primerni za tako reševanje, saj ne posegata v strukturo enacb. Uporabimo ju za diskretizacijo kinematicnih kolicin pri neznanem casu tn+1 z znanimi vrednostmi pri casu tn. Čeloten problem se potem prevede na iskanje rešitev pri nekem diskretnem casu tn+1. Potek reševanja enacšb je odvisen od izbire osnovnih neznank problema. Odlocšimo se, da za osnovne neznanke izberemo deformacijske kolicine. Na&t reševanja enacb lahko strnemo v nekaj tock:
(i)	Ker smo za osnovne neznanke izbrali deformacijske kolicine, najprej upoštevamo kinematicne enacbe in izrazimo preostale kolicine v odvisnosti od izbranih deformacijskih kolicin.
(ii)	Ob znani kinematiki in zunanji obtezšbi lahko iz ravnotezšnih enacšb dolocšimo ravnotezšne sile in momente.
(iii)	Ostanejo še konsistencne enacbe, v katerih upoštevamo konstitutivne enacbe in izpeljane izraze za ravnotezne sile. Konsistencne enacbe skupaj s kinematicnimi in ravnoteznimi robnimi pogoji tvorijo sistem vodilnih enacb pri nekem diskretnem casu tn+1, iz katerih izracunamo osnovne neznanke.
(iv)	Iskanje zveznih rešitev nadomestimo z iskanjem diskretnih rešitev, zato zvezne enacbe diskreti-ziramo, osnovne neznanke pa interpoliramo. Izberemo kolokacijsko metodo diskretizacije, ki je bila uspešno implementirana ze v formulaciji deformacijskih linijskih elementov za staticno analizo [9,17,36]. Diskretizirani sistem enacb rešujemo iterativno, zato ga lineariziramo.
6.1 Kinematicne enacbe
Pri obravnavi kinematicnih enacb (5.1)-(5.4) in (5.9)-(5.10) izberemo dva razlicna pristopa, ki se razlikujeta v izbiri metode diskretizacije kinematicšnih enacšb po cšasu in v izbiri osnovnih neznank.
V	prvem pristopu za diskretizacijo kinematicnih kolicin po casu izberemo Newmarkovo shemo. Za osnovne neznanke diskretizacije kinematicnih kolicin po kraju pa izberemo deformacijska vektorja. S tako izbiro osnovnih neznank zadostimo enemu glavnih ciljev te naloge, saj deformacijske elemente [9, 36], ki so bili pripravljeni za staticno analizo konstrukcij, vpeljemo tudi napodrocje dinamicne analize.
V	drugem pristopu za casovno diskretizacijo izberemo modificirano Newmarkovo shemo, ki je zasnovana na hitrostih in kotnih hitrostih. V poglavju 5.7 smo pokazali, da so krajevni odvodi hitrosti in kotnih hitrosti tudi mere za hitrost deformacije. To nas motivira, da jih izberemo za osnovne neznanke diskretizacije kinematicnih kolicin po kraju.
6.1.1 Newmarkova shema, osnovne neznanke so deformacije
Diskretizacija kinematicnih kolicin z Newmarkovo shemo je zasnovana na izbiri pomikov in rotacij za osnovne neznanke. Pomike in rotacije lahko izrazimo z integracijo krajevnih kinematicnih enacb ob predpostavki, da sta deformacijska vektorja 7^++!] (x) in «^++1] (x) znani funkciji parametra X pri nekem diskretnem casu ti+1. Krajevni vektor poljubne tocke nosilca izrazimo z integracijo enacbe (5.9)
x
rgra+1] (x) = rgn+1] (0) + JR[i+1] (X) (7G++1] (X) — cG[„+1] (X)) dX.	(6.1)
Diferencialne enacbe (5.16) ne moremo neposredno integrirati, ker je zveza med rotacijskim vektorjem in vektorjem deformacij izrazito nelinearna. Resšitev lahko formalno zapisšemo kot
x
4ra+1] (X) = ^gra+1] (0) + /T—T (^gra+1] (XX)) (k^ (X) — dG[„+1, (XX)) dX,	(6.2)
0
dejansko pa jo poišcemo numericno z orodji za reševanje diferencialnih enacb prvega reda.
Enacbi (6.1)-(6.2) veljata za poljuben x, torej tudi za x = L. Zaradi zadošcanja robnim pogojem moramo enacbi, izrazeni pri x = L, upoštevati v sklopu robnih enacb problema
rg1'" iJ = rg"1 J (L) = rg1'"1 iJ -1- ' "L'"1 iJ' -'[ra+1]
L[n+1] = rgn+1] (L) = rg[n+1] + j R[i+1] (x) (7G+1 (x) — cG[„+1, (x)) dx	(6.3)
0
L
tfL[i+1] = ^gn+1] (l) = tfg[n+1] + j T—T (tf^1 (x)) (^+1, (x) — dG[„+1] (x)) dx.	(6.4)
Ko poznamo pomike in rotacije, lahko izrazimo preostale kinematicne kolicine: hitrosti in pospeške, kotne hitrosti in kotne pospeške. Uporabimo Newmarkovo shemo (4.6)-(4.7) in dobimo
Ug (x) = rgi+1] (x) — rg1 (x)
<+1] (X) = ßh Ug + (1 — vi (x)+h (1—agn] (x)
agn+1] (x) = ^ (<+1] (x)—<] (x)) + (1—1) agn] (x) 0G (x) = Spurr (r[i]t (x) R[i+1] (x)'
'GS (X) = ßh0G (X) + (1 — ß J ^Gt] (X) + h ( 1 — ^ ~[1
«GS (X) = ("GS] (X) — -Gin] (X)) + (1 — 1) «Gin] (X).
Če v zgornji shemi upoštevamo enacbe (6.1)-(6.2), dobimo Newmarkovo casovno integracijsko shemo, prirejeno za deformacijske kolicine. V novi shemi namesto pomikov rg (x) in zasukov tfg (x) kot osnovne neznanke nastopajo robne vrednosti pomikov in zasukov rg, tf in deformacijska vektorja 7g (x), kg (x). Izrazi so precej obsezni in neprimerni za numericno implementacijo metode, zato koncnih formul ne podajamo. V racunalniškem algoritmu pomike in rotacije racunamo posebej, dobljene vrednosti pa uporabimo za racun hitrosti, kotnih hitrosti, pospeškov in kotnih pospeškov.
Kolicine, ki nastopajo v kinematicnih enacbah, smo izrazili kot funkcije zveznih vrednosti deformacijskih vektorjev 7'++1] (x) in	(x) ter diskretnih vrednosti pomikov in rotacij na robovih vg[n+1] , $g[n+1] , rL[n+1] in tfL[n+1], ki jih prav tako proglasimo za osnovne neznanke. Zraven pa smo zaradi specificne izbire neznank pridobili nova pogoja (6.3) in (6.4), ki ju dodamo sistemu vodilnih enacb in rešujemo socasno z ostalimi enacbami.
6.1.2 Modificirana Newmarkova shema, osnovne neznanke so hitrosti
Za reševanje casovnih kinematicnih enacb tokrat izberemo modificirano Newmarkovo shemo (4.28), (4.32). Osnovni neznanki nove sheme sta povprecna hitrost vg (x) in povprecna kotna hitrost üG (x). Naša naloga je, da preostale kinematicne kolicine izrazimo z osnovnima neznankama. Za poljubno tocko osi nosilca lahko zapišemo casovne kinematicne enacbe v diskretni obliki v skladu z novo shemo (4.28), (4.32):
V[n+1
Vg
vgn+1 agn+1
R[n+1
[n+ 1
G[n+1
[n+ 1
ü
a
G[n+1
(x) = rgn] (x)+ -vg (x)
(x) =
^ vgn](x) + ß v g (x)
(x) =1agn](x) — ßY- vgn](x) + ßY- v g(x)
(x) = R[n] (x) R(-üG (x)).
(x) = ß^ßßüGm (x) +ß üg (x)
(x) = ^«Gw (x) — -ßYüG>"] (x) + -ßYG (x) .
(6.5)
(6.6)
(6.7)
(6.8) (6.9)
(6.10)
Potrebujemo še zvezo med deformacijskima vektorjema in vektorjema povprecne hitrosti. Ravnamo podobno kot pri izpeljavi izrazov (5.37) in (5.38). Zahtevamo, da rešitvi (6.5) in (6.8) ustrezata kine-maticnim enacbam po kraju:
g[n+1] (x) = R[n+1] (x) ( YGn+i!] (x) — CG[n+1] (x)
R'[n+1] (x) = R[n+1] (x) S [k^ (x) — d'[n+1] (x) Enacbi (6.5) in (6.8) odvajamo po parametru x
rg[n+1] (x) = rg[n](x) + -vg (x)
R'[n+1] (x) = R'[n] (x) R (-üG (x)) + R[n] (x) R' (-üüG (x)) in primerjamo z enacbama (6.11)-(6.12). Dobimo
YG;|++|] (x) — CG[n+1] (x) = R[n+1]T (x) rg[n] (x) + R[n+1]T (x) -vg (x)
(6.11) (6.12)
(6.13)
(6.14)
S K
.[n+1]
"G["+1]
(x) — dG[„+1] (x)j = R[n+1]T (x) R'[n] (x) R (-üG (x))
+ R[n+1]T (x) R[n] (x) R' (-üg (x)).
Odvode pomikov in rotacij pri casu tn izrazimo iz kinematicnih enacb (6.11)-(6.12), ki jih zapišemo pri casu tn. Dobimo izraze, ki deformacije povezujejo s hitrostmi, oziroma njihovimi krajevnimi odvodi.
V
Deformacijski vektor 7G++1] (x) izrazimo kot
YK (x) - eG[n+i] (x) = R (hčg (x))T (7Gnn (x) - cGW (x)) + R[n+1]T (x) hvg (x). (6.15)
Deformacijski vektor kG++1] (x) je izrazen z antisimetricno matriko
S (KG++1] (x) - dG[n+1] (x)) = R (hčG (x))T S («On (x) - dG[„] (x)) R (hčg (x))
+ R (hč g (x))T R' (hč g (x)).
Produkt matrik R (hčG (x))TR' (hčG (x)) je antisimetricna matrika, njen osni vektor pa izrazimo s pomocjo zveze (3.43) in lastnosti (3.40). kG++1] (x) je tako enak
[n+1 ]	T f [n]	\	T /
«G;[n+J1](x) - dG[n+1](x) = R (hčg (x)) (^Gm(x) - dGw(x)J + T (hčG (x)) hčG (x). (6.16)
Enacbi (6.15) in (6.16) predstavljata diskretizacijo deformacijskih vektorjev po casu v skladu z modificirano Newmarkovo shemo. Diskretizacija je kinematicno konsistentna, ker smo jo izpeljali s pomocjo tocnih kinematicnih zvez (6.11)-(6.12).
Pokazali smo, da sta krajevna odvoda povprecne hitrosti v g (x) in povprecne kotne hitrosti č G (x) kolicini, ki (podobno kot deformacijska vektorja) opisujeta deformacijsko stanje v poljubnem prerezu nosilca. Izberemo ju za osnovni neznanki pri reševanju kinematicnih enacb z modificirano Newmarkovo shemo. Hitrosti in kotne hitrosti vzdolz nosilca dobimo z integracijo:
x
v g (x) = v g (0) + J v g (x) dx	(6.17)
0
x
č G (x) = č G (0) + J č G (x) dx.	(6.18)
0
V izrazu (6.18) integriramo vektor zapisan v pomicni bazi, ki je odvisna od parametra x, taka integracija pa je v splošnem lahko zahtevna. V skladu z dogovorom (3.20) je z č G oznacen relativni odvod č G po parametru x. Z enacbo (3.21) smo pokazali, daje integral relativnega odvoda vektorja v pomicni bazi do konstante natancno enak vektorju samemu, zato je v enacbi (6.18) dovolj integrirati komponente relativnega odvoda č'G.
Nosilec je v konstrukcijo umešcen z robnimi vrednostmi pomikov in rotacij: rg (t), r^ (t), R0 (t) in Rl (t). Robne vrednosti pri nekem diskretnem casu tn+1 lahko izrazimo v skladu s casovno diskretiza-cijsko shemo (4.28),(4.32):
r°[n+1] = rg[n]+ h v0	(6.19)
rL[n+1]= rL[n]+ hvL	(6.20)
R0[n+1] = R0[n]R (hčG)	(6.21)
RL[n+1] = RL[n]R (hčG).	(6.22)
rg[n], r^, R0[n] in RL[n] so znane vrednosti robnih pomikov in rotacij pri casu tn. Če predpostavimo, da so pri casu t = 0 robni pogoji zadošceni in so robne vrednosti pomikov in rotacij znane funkcije casa,
jih lahko nadomestimo z robnimi vrednostmi hitrosti vg, vL in kotnih hitrosti &G, & G. Robni pogoji, ki so primernejši za neposredno uporabo skupaj z modificirano Newmarkovo shemo, so tako oblike:
v g (0) = vg	(6.23)
v g (L) = vL	(6.24)
c G (0) = ccG	(6.25)
C g (L)= c G.	(6.26)
Robni pogoji na levem robu nosilca (6.23)-(6.26) so povezani prek enacb (6.17)-(6.18). Tako morajo relativni odvodi hitrosti in kotnih hitrosti zadostiti tudi pogoje na robovih, kar zapišemo s sledecima enacšbama
v g (L) = v0 +1 vg (x) dx = vLa (6.27) 0
L
C G (L) = cG + J & G (x) dx = C G,	(6.28)
0
ki ju dodamo naboru vodilnih enacšb problema.
Predstavili smo resšitve kinematicšnih enacšb nosilca, v katerih kot osnovne neznanke nastopajo zvezne funkcije krajevnih odvodov vektorjev hitrosti in kotne hitrosti v'g (x), CG (x) in robne vrednosti hitrosti in kotnih hitrosti vg, &G, v L in & G.
6.2 Ravnotežne enacbe
Ob predpostavki, da poznamo vse kinematicne kolicine in da so vektorji zunanjih linijskih sil in momentov znane funkcije parametra x, lahko rešitve ravnoteznih enacb (5.23)-(5.24) dobimo z integracijo. Ravnotezne sile pri diskretnem casu tn+1 so enake
x	x
N gn+1l (x) = N gn+1l (0) -J ngn+1l (x) dx + j pA (x) agn+1] (xx) dx,	(6.29)
00
ravnotezšni momenti pa
x	x
M gn+1l (x) = M gn+1l (0) -J rg[n+1] (x) x N gn+1l (x) dx - J mgn+1l (x) dx
00
x
+ J R[n+1] (x) J (xx) Rttn+1l (xx) agn+1l (x) dx 0
x
+ J &gn+1] (xx) x R[n+1] (xx) J (xx) Rt[n+1] (xx) &gn+1] (xx) dx.	(6.30)
0
V enacbah (6.29) in (6.30) nastopata neznani robni vrednosti ravnoteznih sil in momentov Ng (0, t) in Mg (0,t). Obe neznani kolicini dodamo k naboru neznank problema.
Iz strukture integriranih ravnoteZnih enacb lahko identificiramo komponente, ki so posledica notranjih vplivov (l), zunanjih vplivov (e) in vztrajnostnih vplivov (m). Ravnotezno silo in moment lahko potem zapišemo kot vsoto
Ngn+1l (x) = Ng[n+1] (x) + N(x) + N^1 (x)	(6.31)
Mgn+1l (x) = Mg[n+1] (x) + Mg[n+1] (x) + Mm[n+1] (x),	(6.32)
kjer so komponente ravnotezne sile enake
Ng[n+1] (x) = Ngn+1l (0)	(6.33)
x
N eIn+1] (x) = -J ngra+1] (xx)	(6.34) 0
x
Nm[n+1] (x) = j pA (£) agn+1] (xx) ,	(6.35) 0
komponente ravnotezšnega momenta pa
x
Mg[n+1] (x) = M gn+1] (0) -J S (rf+1] (X )) Ng[n+1] (xx) dXX	(6.36)
0
x	x
Mg[n+1] (x) = -J S (rf+1] (XX )) Ng[n+1] (X) -J mgn+1] (X)	(6.37)
00
x	x
M m[n+1] (x) = -J S (Vf+1] (XX )) n m[n+1] (X) + J R[n+1] (xx) J (xx) aj^ (X) dX
00
x
+ J S (R[n+1] (XX) ^1] (XX)) R[n+1] (XX) J (X) ^++1] (XX) dX.	(6.38)
V gornjih enacbah smo vektorske produkte zapisali z uporabo pravila (3.29). Kotne hitrosti in pospeške pa smo zapisali v pomicni bazi: RT[ra+1]^gra+1] = <^jn++1] in RT[ra+1]agn+1] = ajGj++1]. Komponentni zapis ravnoteznih sil in momentov omogoca preglednejšo linearizacijo, pomemben pa je tudi za izvedbo numericnega dušenja, ki ga prikazemo v poglavju 6.5.
6.3 Vodilne enacbe dinamike prostorskega linijskega elementa
Opravili smo delo, ki smo ga zapisali v tockah (i) in (ii) na zacetku tega poglavja. Izmed celotnega nabora enacb dinamike prostorskega nosilca smo del enacb formalno zadostili, ostale so enacbe konsi-stence, konstitutivni zakoni ter naravni robni pogoji ravnotezšja sil na robovih. V obeh pristopih resševanja kinematicnih enacb smo dobili dodatni enacbi, kiju dodamo preostalim enacbam. Skupaj tvorijo vodilni sistem enacb dinamike prostorskega nosilca, ki jih rešujemo s pomočjo pomoznih enacb, predstavljenih v tockah 6.1 in 6.2 tega poglavja.
Zaradi razlicnih metod za diskretizacijo kinematicnih enacb smo dobili tudi dva razlicna nabora osnovnih neznank. V nadaljevanju zapišemo vodilna sistema enacb za vsako izbiro osnovnih neznank posebej.
6.3.1 Element z interpolacijo deformacij
Pri resševanju kinematicšnih enacšb z Newmarkovo shemo smo za osnovni neznanki izbrali deformacijska vektorja 7G (x) in kg (x), neznane pa so tudi robne vrednosti pomikov in rotacijskih vektorjev: rg, tf, rL, tf. Z integracijo ravnoteznih sil in momentov smo pridobili še neznana vektorja ravnoteznih sil in momentov na levem robu elementa N g in M 0
Vodilni sistem enacb sestavljajo konsistencni pogoji (5.31)-(5.32), robni enacbi za pomike in zasuke krajišc nosilca (6.3)-(6.3) ter naravni robni pogoji (5.25)-(5.28):
f 1 = R[i+1] (x) N G[i+1] (x) — N g1+1] (x) = 0	(6.39)
f 2 = R[i+1] (x) MG[i+1] (x) — Mgi+1] (x) = 0	(6.40)
f 3 = rL[i+1]— rg1+1](L)= 0	(6.41)
f 4 = tf [i+1] — tf1+1] (L) = 0	(6.42)
f 5 = s0[i+1] + Ng1+1] (0) = 0	(6.43)
f 6 = P0[i+1] + Ml1+1] (0) = 0	(6.44)
f 7 = S L[i+1] — N g1 +1] (L) = 0	(6.45)
f 8 = PL[1+1] — M[1+1] (L) = 0.	(6.46)
Enacbe rešujemo numericno, zato vpeljemo interpolacijo komponent deformacijskih vektorjev. Funkciji 7g++1] (x) in kG++1] (x) nadomestimo z diskretnimi vrednostmi v poljubno izbranih interpolacijskih tockah xp, p = 1,..., N in interpoliramo z interpolacijskimi funkcijami Ip (x):
N
7G+1](x) = I Ip (x) 7PGr„^1]]	(6.47)
p=1
N
kG+1](x) = I Ip (X) kG!1K;].	(6.48)
p=1
Diskretizacija osnovnih neznank omogoca diskretizacijo zveznih enacb (6.39)-(6.40) glede na parameter x po metodi koncnih elementov. Izberemo kolokacijsko metodo in enacbam (6.39)-(6.40) zadostimo le v N izbranih diskretnih tockah xq, q = 1,..., N:
f 1,q = R[1+1] (Xq) NG[1+1] (Xq) — NI^1 (Xq) = 0	(6.49)
f 2,q = R[1+1] (Xq) MG[1+1] (Xq) — M[1+1] (Xq) = 0.	(6.50)
p[1+ 1 ]
Dobimo sistem 2N + 6 diskretnih nelinearnih vektorskih enacšb za 2N + 6 osnovnih neznank 7G[n+1] ,
kP[1+1] r0[™+1] tf0[1+1] rL[1+1] tfL[1+1] N0[1+1] in M0[1+1] p = 1 N KG[n+1], 'g , tfg ,'g , tfg , Ng in Mg ,p = 1,...,N.
6.3.2 Element z interpolacijo krajevnih odvodov hitrosti
V drugem pristopu smo za cšasovno diskretizacijo kinematicšnih enacšb uporabili modificirano Newmar-kovo shemo, za osnovne neznanke pa smo izbrali krajevne odvode povprecne vrednosti hitrosti v g (x) in kotne hitrosti — G (x) ter robne vrednosti povprecnih hitrosti in kotnih hitrosti: v -G, v L, L. Z
integracijo ravnoteZnih sil in momentov smo pridobili še neznana vektorja ravnoteZnih sil in momentov na levem robu elementa N g in M 0.
Vodilni sistem enacb tudi v tem primeru sestavimo iz konsistencnih pogojev (5.31)-(5.32), kinematicnih robnih enacb (6.27)-(6.28), ki smo jih zapisali v odvisnosti od hitrosti in kotnih hitrosti, in naravnih robnih pogojev (5.25)-(5.28):
hi = R[n+1] (x) NGK1 (x) - Ngn+11 (x) = 0	(6.51)
h2 = R[n+11 (x) MGK1 (x) - Mgn+11 (x) = 0	(6.52)
h3 = vL[n+11 - vgn+11 (L) = 0	(6.53)
h4 = c£[n+11- ^!n+11(L)= 0	(6.54)
h5 = s0[n+11 + n gn+11 (0) = 0	(6.55)
h6 = P0[n+11 + M gn+11 (0) = 0	(6.56)
h7 = SL[n+11 - Ngn+11 (L) = 0	(6.57)
h8 = PL[n+11 - Mgn+11 (L) = 0.	(6.58)
Diskretizacija po metodi koncnih elementov temelji na interpolaciji osnovnih zveznih neznank, ki jih nadomestimo z diskretnimi vrednostmi in interpolacijo:
N
v g (x) = £ Ip (x) <	(6.59)
p=1 N
'G (x) = £ Ip (x) č£'.	(6.60)
p=1
■VP' in u>G so vrednosti osnovnih neznank v poljubno izbranih interpolacijskih tockah xp, p = 1,... ,N, Ip (x) pa so interpolacijske funkcije. Zveznim enacbam (6.51)-(6.52) tudi v tem primeru zadostimo v poljubno izbranih tockah xq, q = 1,..., N:
h1,q = R[n+11 (Xq) NG[n+11 (Xq) - Ng*11 (Xq) = 0	(6.61)
G
h2,q = R[n+11 (Xq) MG[n+11 (Xq) - Mg^11 (Xq) = 0.	(6.62)
Sistem enacb (6.61)-(6.62), (6.53)-(6.58) je sistem 2N + 6 diskretnih nelinearnih enacb, v katerem kot osnovne neznanke nastopajo: vp', 'G', vg, 'G, vL, 'L, Ng[n+11 in Mg[n+11, p = 1,..., N.
6.4 Reševanje diskretnih enacb
Diskretna sistema enacb (6.49)-(6.50), (6.41)-(6.46) in (6.61)-(6.62), (6.53)-(6.58) rešujemo iterativno z Newtonovo metodo. Diskretne enacbe lineariziramo in pojasnimo dodajanje iterativnih popravkov, ki jih dobimo z rešitvijo lineariziranega sistema. Linearizacijo zapišemo za obe formulaciji vodilnih enacb elementa v pregledni obliki, ki je primerna za programiranje. Zaradi boljše preglednosti izrazov opušcamo oznako za diskreten cas tn+1, pri katerem rešujemo enacbe.
6.4.1 Element z interpolacijo deformacij Linearizacija
Enacbe (6.49)-(6.50) in (6.41)-(6.46) variiramo okrog trenutne iterativne konfiguracije v smereh osnovnih neznank: čyG, čkG, čr0, č^, črL, čtfj, čNg in čM0. Preden zapišemo variacije celotnega sistema enacb, pripravimo variacije kolicin, ki nastopajo v enacbah.
Variacije zveznih vrednosti deformacijskih vektorjev dobimo z upoštevanjem interpolacije (6.47)-(6.48):
N
ČYG (x) = £ Ip (x) čyG	(6.63)
P=i N
čkg (x) = £ /p (x) čkG .	(6.64)
p=i
Nadaljujemo z variiranjem krajevnih kinematicnih enacb. Variacija odvoda krajevnega vektorja po kraju (5.9) je enaka
črg (x) = -S (R (x) (7g (x) - cg (x))) čtfg (x) + R (x) č7g (x)
= -S (rg (x)) čtfg (x) + R (x) č7g (x).	(6.65)
Variirajmo še kinematicno enacbo za rotacije (5.10):
čR' (x) = čR (x) S (kg (x) - dG (x)) + R (x) S (čkg (x))
= S (čtfg (x)) R (x) S (KG (x) - dG (x)) + R (x) S (ČKG (x)) = S (čtfg (x)) R' (x) + R (x) S (čkg (x)).
Izraziti zelimo variacijo rotacijskega vektorja č$g (x), zato po kraju odvajamo izraz za variacijo rotacijske matrike (3.44):
čR' (x) = S (čtfg (x)) R (x) + S (čtfg (x)) R' (x).
S primerjavo obeh zapisov variacij odvoda rotacijske matrike dobimo
S (čtfg (x)) = R (x) S (čkg (x)) Rt (x)
čtfg (x) = R (x) čkg (x).	(6.66)
Izraza (6.65)-(6.66) predstavljata tocne kinematicne enacbe Reisnerjevega nosilca v šibki obliki. Variacije pomikov in rotacij dobimo kar z integracijo šibkih zvez (6.65)-(6.66):
x	x
črg (x) = črg -J S (rg (X)) čtfg (x) dx + J R (x) čy g (x) dx	(6.67)
00
x
č^g (x) = čtfg + J R (x) čkg (x) dx.	(6.68)
Integrirane izraze poenostavimo tako, da upoštevamo šibka interpolacijska nastavka (6.63)-(6.64) in vpeljemo nekaj pomozšnih kolicšin:
0
x
N f
R (X) G (X) dX = £ / R (X) Ip (X) dX ^yg p=10 N x
/ R (X) 5kg (XX) dX = £ / R (X) Ip (X) dX
0 p=10
x	x	x
/s (rg (X)) čtfg (XX)dX = JS (rg (X)) drččtfj +/S (rg (XX)) Wp(X) drč 0 0 0
Dobimo kompakten zapis variacij pomikov in rotacij:
N	N
5rg (X) = 5rg - U (x) - £ Vp (x) + £ Wp (x) 5yg	(6.69)
p=1 p=1
N
čtfg (x)= ^ + £ Wp (x) <jk£,	(6.70)
p=1
kjer so
U (x) = J S (rg (X)) dX 0
x
Wp (x) = J R (X) Ip (X) dX 0
x
Vp (x) = J S (rg (X)) Wp (X) dx.
'p
0
Iz znanih variacij pomikov in rotacij lahko izrazimo variacije hitrosti in pospeškov, ki smo jih v zapisali v enacbah (4.15)-(4.16) in (4.19)-(4.20):
y (	n	n	\
(X) = ßh r r0 - U (x) - £ Vp (x) + £ Wp (x) J	(6.71)
i /	N	N	\
1 0 0 p p
p=1 p=1
5ag (x) = ßh2 ( 5rg - U (x) - £ Vp (x) + £ Wp (x) )	(6.72)
^ (x) = ßht(0G (x))-TRt(x) f ^ + £ Wp (x) j	(6.73)
1N
1	,	, £ Wp(x).
5aG (x) = ßh2 T (Ög (x))-T RT (x) ^+ £ Wp (x) .	(6.74)
Pripravili smo vse potrebno za variiranje ravnoteznih sil in momentov. Zaradi preglednosti izkoristimo aditivni razcep ravnoteznih sil in momentov po enacbah (6.31)-(6.32). Variaciji sumandov, ki sta posledica notranjih (6.33) in zunanjih sil (6.34), sta trivialni:
5Ng (x) = 5Ng	(6.75)
5Ng (x) = 0.	(6.76)
x
Variacijo vztrajnostnega dela ravnotezne sile (6.35) formalno zapišemo kot
x
SNm (xx) = J pASag (X)
0
N	N
: NN m (x) Sr0+nn m (x) s^g+£ n m (x skg+£ n m (x) s7g,	(6.77)
p=i p=i
kjer je
x
N^ (xx) = J pAIdxX	(6.78)
0
x
1
NN m (x) = - ßh2 I pAU (XX)	(6.79)
0
x
NmmP (x) = -ßh-2 J pAVp (XX)	(6.80)
0
x
N(x) = ßh2 J PAWp ) .	(6.81)
0
Z variiranjem komponent ravnoteznega momenta (6.32) imamo vec dela. Najprej pripravimo osnovne izraze:
x	x
SMg (xx) = SM g (0) -J S (rg (xx)) SNg (XX) + J S (Ng (XX)) Srg (xx) dxx	(6.82)
00
x	x
SMg (xx) = -J S (rg (xX)) SNg (x) dxX + J S (Ng (xX)) Srg (xX) dxXX	(6.83)
00
x	x
SMm (xx) = -J S (rg (x)) SNm (x) dxx + J S (Nm (x)) Srg (x) dxx
00
x	x
+ / SR (xx) J (x) aG (xx) dxx + J R (xx) J (x) SaG (xx) dxx
00
x
-	J S (R (xx) J (x) ^g (xx)) SR (xx) ^g (xx) dxx 0
x
-	J S (R (xx) J (x) (xx)) R (xx) S^g (xx) dxx 0
x
+ / S (R (xx) (xx)) SR (xx) J (x) (xx) dxx
0
x
+ J S (R (xx) ^G (xx)) R (xx) J (xx) S^G (xx) dxx,	(6.84)
0
nato pa upoštevamo ze izpeljane izraze za variacije kinematicnih kolicin. Variacijo SMg lahko zapišemo
kot
N	N
ČMg (x) = ČMg + MNČNg + M^tfg + £ Mk,p (x) + £ m^ (x) Č7G,
p=1 p=1
pomozšne matrike pa so enake
x
Mi N (x) = -J S (rg (X)) dX
0
x
M (x) = -J S (Ng (X)) S (rg (X)) dX 0
x
MK,p (x) = -J s (Ng (X)) S (rg (X)) Wp (X) dX
g
0
x
(x) = J S (Ng (X)) /p (X) R (X) dX. 0
Podoben je tudi formalen zapis variacije clena, kije posledica zunanjih sil
NN
čMg (x) =	+ £ mK,p (x) + £ m^ (x) Č7G,
p= 1 p= 1
kjer je
x
Mi: (x) = -J S (Ng (X)) S (rg (X)) dX 0
x
MeK,p (x) = -J S (N g (X)) S (rg (X)) Wp (X) dx 0
x
MY,p (x) = J S (Ng (XX)) /p (XX) R (XX) dX. 0
Ostane še variacija vztrajnostnega clena ravnoteznega momenta, ki jo zapišemo kot
NN
(x)=M^g+M?wg+£ m™ (x) ^g+£ m™ (x) ^g .
p= 1 p= 1
Pomozšne matrike pri variacijah so enake
x
= -J S (rg (X)) Nm (XX) dX 0
x	x
Mm (x) = -J S (Nm (XX)) S (rg (XX)) dX -J S (rg (XX)) Nm (X) dX
00
x	x
- S (R (X) J(X) aG (XX)) dX + / S (R (XX) J (XX) (XX)) S (R (XX) ^g (XX)) dX
00
x
- J S (R (XX) (XX)) S (R (XX) J(XX) (XX)) dX 0
x
+ ßh2 /R (x) J (x) T (Ög (X))-TRT (X) dx
-	f S (R (xx) J (X) čg (X)) R (X) T (Ög (X))-T RT (X) dx ßh J
0
x
+ ßh J S (R (XX) čG (X)) R (X) J (X) T (ÖG (X))-T Rt (XX) dx ßh 0
x	x
M^ (x) = -J S (Nm (X)) S (rg (X)) Wp (X) dx -J S (rg (X)) N^ (X) Wp (X) dx
00
x
-/ S (R (X) J (X) aG (X)) Wp (X) dx 0
x
+ J S (R (X) J (X) čG (X)) S (R (X) čG (X)) Wp (X) dx 0
x
-	J S (R (X) čG (X)) S (R (X) J (X) čG (X)) Wp (X) dx 0
x
+ ßh2 / R (XX) J (X) T (ÖG (XX))-T RT (X) Wp (xx) dx
0
x
-	f S (R (X) J (X) čG (X)) R (X) T (ÖG (X))-T Rt (X) Wp (X) dx ßh
0
x
+ ßh J S (R (XX) čG (xx)) R (xx) J (xx) T (Ög (X))-T RT (X) Wp (xx) dx ßh 0
x	x
(x) = / S (N^ (X)) /p (X) R (X) dx -J S (rg (X)) N :p (X) dx.
Za linearizacijo konstitucijskih zvez (5.29)-(5.30) zadošca dolociti parcialne odvode operatorjev in v odvisnosti od komponent deformacijskih vektorjev 7G in kg. Parcialne odvode zaradi preglednosti oznacimo kot
C
77
dCf 97.
J J
C

d
j j
C

97.
J J
C
dKJ J
i, j = 1,2,3. (6.85)
Variacije konstitucijskih notranjih sil in momentov so potem enake
čN G (x) = C77 (x) ^7g (x) + C7k (x) čkg (x) ČMG (X) = Ck7 (X) g (X) + Ckk (x) čkg (x) .
(6.86) (6.87)
V konsistencnih enacbah nastopajo variacije konstitutivnih sil v nepomicni bazi. Z uporabo izpeljanih izrazov za variacijo rotacijske matrike in prostorske diskretizacije variacij deformacijskih vektorjev
dobimo
N
S (R (x) N G (x)) = —S (R (x) N G (x)) Sd°g — S (R (x) N G (x)) £ Wp (x) S kG
p=1
NN
+ £ Ip (x) R (x) C77 (x) SjPg + £ Ip (x) R (x) C7K (x) SkG	(6.88)
p=1 p=1
N
S (R (x) M G (x)) = —S (R (x) M G (x)) Sd°g — S (R (x) M G (x)) £ Wp (x) SkGg
p=1
NN
+ £ Ip (x) R (x) CK1 (x) SjPg + £ Ip (x) R (x) Ckk (x) S KG.	(6.89)
p=1 p=1
Koncno lahko zapišemo linearizacijo celotnega sistema enacb (6.49)-(6.50) in (6.41)-(6.46) kot linearno kombinacijo variacij osnovnih neznank sistema:
Sf 1,q = — N m (xq) Sr°g — [S (R (xq) N G (xq)) + NN Z (xq)] Sd°g — SN 0
N
+ £ [Ip (xq) R (xq) Cyk (xq) — S (R (xq) N G (xq)) Wp (xq) — NN ™p (xq)] SkG
p=1
N [ ]
+ £ [Ip (xq) R (xq) Cyy (xq) — NN™ (xq)] S^G p=1
Sf 2,q = —Mm (xq) Sr°g — [S (R (xq) M G (xq)) + M (xq) + M (xq) + Mg (xq)] Sd0g
N
—	M N (xq) SN g — SM°g + £ [Ip (xq) R (xq) CKK (xq)
p=1
—	S (R (xq) M G (xq)) Wp (xq) — M^ (xq) — M^ (xq) — M™p (xq)] SKG
N [ ]
+ £ [Ip (xq) R (xq) CK1 (xq) — M^p (xq) — M^ (xq) — M^ (xq)] SjG
p=1
NN
Sf 3 = —Sr°g + SrLg + U (L) S^ + £ Vp (L) SkG — £ Wp (L) SYG
p= 1 p= 1
N
Sf 4 = —+ S&L — £ Wp (L) SkG
p=1
Sf 5 = sn g Sf 6 = sm g
NN
Sf 7 = — N m (L) Sr°g — N m (L) S< — SNg — £ N Z (L) SkG — £ N ™ (L) SyG
p=1 p=1
Sf 8 = — Mg (L) Sr0 — [M (L) + M (L) + Mg (L)] Sd°g — M N (L) SNg N [ ]
—	SM g + £ [—MK,p (L) — M %,p (L) — MZ (L)] S kG
p=1
N [ ]
+ £ [—Ml(,p (L) — MY,p (L) — MZp (L)] Syg.
p=1
Dodajanje popravkov
Vzemimo, da poznamo vrednosti vseh kolicin pri i-ti iteraciji. Z rešitvijo lineariziranega sistema vodilnih enacb dobimo iterativne popravke osnovnih neznank: Č7P G., čk'g. , Čr0 g, Čtf0g, črLg, čtfLg, ČN0g in 5M0g .	. • . • • ' , '	'
Pomiki, ravnotezne sile in deformacijski vektor 7g so aditivne kolicine, zato lahko njihove popravke kar prištejemo trenutnim iterativnim vrednostim:
r0 = r0 + Čr0
i+1,g i,g ' J i,g
rL = rL + črL r i+1,g r i,g + Čr i,g
N 0+1,g = N 0g + ČN 0g M 0+1g = M 0g + ČM 0g
7 p+1,Gi+1 = 7 p,Gi + Č7 'iGi
Zaradi nelinearnosti prostorskih rotacij moramo vecš pozornosti nameniti kolicšinam, povezanim z rotacijami. (Če zelimo zadostiti robnemu pogoju za rotacije na desnem robu elementa, se morajo iterativni popravki rotacijskega vektorja ujemati z njihovimi variacijami. Iz kinematicne enacbe (6.68) potem sledi:
x
Čtfi,g (x) = Čtf0g + J Ri (X) ČKi)Gi (X) dX. 0
Aditivni popravek deformacijskega vektorja AaKi,Gi+1 potem izracunamo iz krepke oblike kinematicne enacbe (5.16), v kateri upoštevamo, da se variacijska konstanta dG ne spreminja. Podrobnosti izpeljave najdemo v clanku [9], rezultat pa je enak
ČAKi,Gi+1 = TT (RTČtfj;g) ČKijGi.
Tako dobljen popravek deformacijskega vektorja kg preprosto prištejemo trenutni iterativni vrednosti:
Ki+1,Gi+1 = Ki,Gi + ČAKi,Gi+1.
Ker smo izbrali, da se popravki rotacij ujemajo z njihovimi variacijami, popravke rotacij kar dodamo trenutnim iterativnim vrednostim v skladu z multiplikativno naravo rotacij:
R (tf0+1g) = R (Čtf0g ) R (tf0,g )
R (tfL+1,g ) = R (ČtfLg ) R (tfLg ) .
Popravljene vrednosti rotacijskih vektorjev tfi+1,g pa izracunamo iz rotacijskih matrik s pomocjo Spurri-erjevega algoritma [40]. Racun popravljenih vrednosti hitrosti in pospeškov pri diskretizaciji z Newmarkovo shemo smo predstavili v poglavju 4.1.2.
6.4.2 Element z interpolacijo krajevnih odvodov hitrosti Linearizacija
Sistem enacb (6.61)-(6.62), (6.53)-(6.58) variiramo v smereh osnovnih neznank Čv'', č—g, Čv^, č—g,
ČvL, Č— L, ČNg in ČMj.
Najprej zapišimo variacije zveznih vrednosti osnovnih neznank, ki jih dobimo z variiranjem interpolacijskih nastavkov (6.59)-(6.60):
N
ČVg (x) = £ /p (X) čvp' p=1 N
Čč G (x) = £ /p (x) ČčG'. p=1
Variacije povprecnih hitrosti in kotnih hitrosti dobimo z integracijo gornjih izrazov:
ČVg (x) = ČVg (0)+ £ f /p (X) dXČVp' = ČVg (0) + £ Qp (x) ČVp'
p=1 0 p=1
N x	N
ČčG (x) = ČčG (0) + £ / /p (X) dxČčg = ČčG (0) + £ Qp (x) Čč
p=1 p=1
Variacije preostalih kinematicnih kolicin smo predstavili ze v poglavju 4.2.1. Ponovno jih zapišimo z uposštevanjem diskretizacije po kraju:
N
črg (x) = hČVg (x) = h £ Qp (x) ČVp'	(6.90)
p=1 g
11
N
ČVg (x) = -ČVg (x) = - £ Qp (x) ČVp'	(6.91)
P	P p=1
11
N
Čag (X) = ^YhČVg (X) = ßYh £ qp (X) ČVp'	(6.92)
ČR (x) = R[n1 (x) ČR (hč g (x))	(6.93)
N
ČR (x) u = R (x) S (u)T T (hčg (x))t £ Qp (x) hČčg	(6.94)
p=1	G
1	1 n
ččg (x) = ^Čccg (x) = - £ Qp (x) ČccG'	(6.95)
ß ^ ß p=1
ČčG (x) = "J1, £ ßYh	ßYh p=1
p=
N
ČaG (x) =	G (x) = ß-h £ qp (x) ČčG'.	(6.96)
Deformacije izrazimo iz kinematicnih zvez (5.9)-(5.10)
YG (x) = Rt (x) rg (x) + CG (x) S (kg (x)) = R (x)T R' (x)+ S (dG (x)),
izraza variiramo, pri tem pa upoštevamo, da sta cG in dG variacijski konstanti:
ČYG (x) = ČRT (x) rg (x) + Rt (x) Črg (x) S (ČKG (x)) = ČRT (x) R' (x) + Rt (x) ČR' (x).
Variacijo transponirane rotacijske matrike izrazimo z variiranjem lastnosti (3.1)
rti	rTI	rTI	rr~T	rrT
ČRRT + RČRT = 0 ^ ČRT = -Rt ČRRT.
Zadnji clen variacije deformacijskega vektorja kg vsebuje variacijo krajevnega odvoda rotacije. Operaciji variiranja in odvajanja po kraju lahko zamenjamo, zato najprej zapišemo krajevni odvod produkta matrik RT (x) čR (x)
^ (rt (x) čR (x)) = Rt (x) čR' (x) + Rt' (x) čR (x), od koder sledi, daje produkt matrik RT (x) čR' (x) enak
RT (x) čR' (x) = d (RT (x) čR (x)) - RT' (x) čR (x) dx V	)
dx
= — f RT (x) čR (x)) + RT (x) R' (x) RT (x) čR (x).
dx v
Variacije deformacij potem izrazimo kot
d
čY G (x) = -RT (x) čR (x) RT (x) rg (x) + RT (x) — črg (x)
S (čkg (x)) = -RT (x) čR (x) RT (x) R' (x) + RT (x) R' (x) RT (x) čR (x) + | (rt (.)čR (.)).
Iz kinematicnih enacb (5.9)-(5.10) sledi, daje produkt RT rg enak deformaciji yg - cg, produkt RT R' pa S (kg - dG). Variacije pomikov in rotacij smo ze zapisali z enacbami (6.90) in (6.93). Poenostavimo še produkt RTčR s pomočjo pravila (3.50) za variacijo rotacijske matrike, parametrizirane z aditivnim rotacijskim vektorjem
RTčR = RT (hčg) R[n]TR[n]čR (hčg) = RT (hčg) čR (hčg) = S (Tt (hčg) hččg) .
Variacijo deformacijskega vektorja yg potem zapišemo kot
( ) d
čYg (x) = -S (Tt (hčg (x)) hččg (x)) (yg (x) - cg (x)) + RT (x) — (hč«g (x))
= S (yg (x) - cg (x)) T (hčg (x)) hččg (x) + RT (x) hčvg (x).	(6.97)
Antisimetricno matriko, sestavljeno iz komponent variacije deformacijskega vektorja kg, pa izrazimo kot
S (čkg (x)) = - S ^Tt (hč g (x)) hčč G (x) j S ^ kg (x) - dG (x) + S ^kg (x) - dG (x) jS ^Tt (hčg (x)) hččG (x) + S (TT (hčg (x)) hččG (x)) Če upoštevamo lastnost antisimetricnih matrik (3.27), se zgornji izraz poenostavi: S (čkg (x)) = S ^S (kg (x) - dG (x)) Tt (hčG (x)) hččG (x) + S ((TT (hčg (x)) hččG (x)) ) . Od tod pa sledi, daje variacija deformacijskega vektorja čkg enaka
čkg (x) = S (kg (x) - dG (x)) Tt (hčG (x)) hččG (x) + Tt' (hčG (x)) hččG (x)
+ Tt (hč g (x)) hčč G (x).	(6.98)
Ko poznamo variacije kinematicnih kolicin, je variiranje preostalih enacb preprosto. Ravnotezne enacbe zopet variiramo glede na adicijski razcep. Variacije sumandov v izrazu za ravnotezno silo so enake:
ČNg (x) = ČNg ČNg (x) = 0
x	N
ČN? (x) = f pAčag (XX) dX = IN? (x) Čvg + £ N(x) Čvp', 0 p=1
kjer je
x
Nm (x) = ß1^/pAIdX 0
x
NN m,p (x) = ß^f pAQp (x) Idx.
Variacije ravnoteznega momenta lahko ponovno zapišemo v splošni obliki (6.82)-(6.84), nato pa upoštevamo variacije kinematicnih kolicin. Variacijo clena ravnoteznega momenta, ki pripada notranjim silam, izrazimo z
N
ČMg (x) = ČMg + MN (x) ČNg + £ M^,p (x) Čvp',
p=1
pri tem sta pripadajoci matriki enaki:
x
Mi N (x) = -J S (rg (XX)) dX
0
x
MiV',p (x) = h J S (Ng (X)) /p (X) dx. 0
Mocno se poenostavi tudi variacija clena ravnoteznega momenta zaradi zunanjih sil
N
ČMg (x)= £ MlV',p (x) čvp',
p=1
kjer je
x
MV',p (x) = h J S (Ng (XX)) /p (XX) dX. 0
Razcep variacije vztrajnostnega cšlena ravnotezšnega momenta zapisšemo z vsoto
NN
čMm (x)=m? (x) čvg+£ Mim?,p (x) č«p+m? (x) č*g+£ Mm,p (x) č*g ,
p= 1 p= 1
pripadajocše matrike pa so enake
x
M $ (x) = -J S (rg (X)) N: dx 0
x	x
M m?>p (x) = -J S (rg (X)) N ^dx + h J S (N$ (X)) /p (X) dx
x	x
M^i (x) = -h J R (X) S (J (x) «G (xXX)) Tt (hč G (X)) dx + 7^ J R (X) J (x) dx 0	0
x
+ hj R (X) S (S (J (x) čG (X)) čG (X)) Tt (hč g (X)) dx 0
x
+1J R (X) (S (čG (X)) J (x) - S (J (x) čG (X))) dx
ß0
x	x
M:m,p (x) = -h JR (X) S (J (x) «G (X)) Tt (hčG (X)) Qp (X) dx + — JR (X) J (x) Qp (X) dx 0	0
x
+ h R (X) S (S (J (x) čG (X)) čG (X)) Tt (hčG (X)) Qp (X) dx
+1 / R (X) (S (čg (X)) J (x) - S (J (x) čg (X))) Qp (X) dx.
0
x
0
Variacije konstitucijskih sil tako kot v izrazih (6.86) in (6.87) zapišemo s pomocjo parcialnih odvodov konstitucijskih operatorjev CN in po komponentah deformacijskih vektorjev. Ko upoštevamo še izraze za variacije deformacijskih vektorjev (6.97) in (6.98), dobimo
čNG (x) = C77 (x) Rt (x) hčvg (x) + C77 (x) S (7g (x) - cg (x)) T (hčg (x))t
TT
+C7K (x) S (kg (x) - dG (x)) T (hčg (x)) + C7K (x) T (hčg (x)) hččg (x) + C7K (x) T (hč g (x))T hčč G (x) ČM G (x) = Cky (x) Rt (x) hčv g (x) + [ck7 (x) S (y g (x) - cg (x)) T (hč g (x))t
+Ckk (x) S (kg (x) - dG (x)) T (hčg (x)) + CKK (x) T (hčg (x)) hččg (x) + Ckk (x) T (hčč g (x))T hčč G (x)
Variacije konstitucijskih sil potrebujemo v konsistencnih enacbah zapisane v nepomicni bazi. Zaradi boljšega pregleda pa jih zopet zapišemo z vsotami
N	N
č (R (x) N G (x)) = £ N ^ (x) čvp' + NN C0 (x) ČččGG + £ N£,p (x) ČččG	(6.99)
p=1 p=1 N	N
č (R (x) MG (x)) = £ M^ (x) čvp' + MC0 (x) čč>G + £ MC',p (x) ččg.	(6.100)
p= 1 p= 1
Pripadajoce matrike, neodvisne od variacij, so enake: N (x) = R (x) C77 (x) Rt Ip (x)
INC0 (x) = hR (x) [C77 (x) S (7g (x) - cg (x)) + C7K (x) S (kg (x) - dg (x))
-S (NG (x))] Tt (hčg (x)) + hR(x) C7K (x) Tt' (hčg (x)) NNC')P (x) = hR (x) [C77 (x) S (7g (x) - cg (x) ) + C7K (x) S (kg (x) - dg (x))
-S (NG (x))] Tt (hčg (x)) Qp (x) + hR (x) C7k (x) Tt' (hčg (x)) Qp (x) + hR (x) C7K (x) Tt (hč g (x)) Ip (x) MM (x) = R (x) Ck7 (x) Rt Ip (x)
MlC0 (x) = hR (x) [Ck7 (x) S (7g (x) - cg (x)) + CKK (x) S (kg (x) - dG (x))
-S (MG (x))] Tt (hčg (x)) + hR (x) Ckk (x) Tt' (hčg (x)) M^C',p (x) = hR (x) [Ck7 (x) S (7g (x) - cg (x)) + CKK (x) S (kg (x) - dG (x))
-S (MG (x))] Tt (hčg (x)) Qp (x) + hR (x) CKK (x) Tt' (hčg (x)) Qp (x) + hR (x) Ckk (x) Tt (hč g (x)) Ip (x).
Pripravili smo vse potrebne kolicine za zapis variacij diskretnih enacb (6.61)-(6.62), (6.53)-(6.58). Vari-irane enacbe tvorijo linearen sistem enacb elementa z interpoliranimi krajevnimi odvodi hitrosti in kotnih hitrosti, ki ga v odvisnosti od variacij osnovnih neznank Süp', Sčg, S^j, Sčg, S1«L, SčL, SNg in SM^ zapisšemo kot
N ( )
Sh1 (xq) = -NNm? (xq) Sü0 + NC0 (xq) Sc^g - SN^ + £ (iNC?,p (x,) - NN^p (x,)) S<
p=1
N
+ £ iNC',p (xq) ScčG
p=1
Sh2 (xq) = -mmms (xq) S«0+(NC0 (xq) - »m (%)) Scč0 - MMN (%) SNg - SMg N ( )
+ £ (M^p (xq) - MV',p (xq) - M^p (^g) - »/.p (x,)) Süg
p=1
N ( )
+ £ (MlC',p (x,) - Mmm,p (x,)) Scčg
p=1
N
Sh3 = -Sü0+SvL - £ Qp (L) Süg
p=1 N
Sh4 = -Scč0 + SčL - £ Qp (L) Sčg
p=1
Sh5 = SNg Sh6 = SMg
N
Sh7 = -iNmm (L) Sv0 - SNg - £ iNm?,p (L) S<
p=1
Sh8 = -»mm (L) Svg - »m (l) Scg - min (l) SNg - SMg N ( ) N
+ £ (-M;,p (L) - MV,p (L) - M^p (L)) S«G - £ jMm?;p (L) Sčg. p=1 p=1
Dodajanje popravkov
Vse osnovne neznanke elementa z interpolacijo krajevnih odvodov hitrosti so aditivne kolicine. Racun novih priblizkov osnovnih neznank se zato zelo poenostavi. Denimo, da poznamo vse kolicine v i-ti iteraciji reševanja vodilnih enacb. Z rešitvijo lineariziranega sistema enacb dobimo popravke osnovnih neznank čvPg, 5c-P'G, 5-V0g, 5"0G, čv^, , 5N0,fl in 5M0,g, ki jih prištejemo trenutnim iterativnim vrednostim:
v 1,g = v tfl + 5v &
" p+ 1,g = " p,G+5" pG v0 _ v0 + 5.-0
"0+1,G _ "°,G + 5"0,G -T	-T c - T
vm,ö _ vt,ö + 5vt,ö "t+1,g _ "ttg + 5"ttg Ni+1,g (0) _ Ni,g (0) + 5Ni,g (0) M i+1,g (0) _ M i,g (0)+ 5M i,g (0) .
Nove iterativne vrednosti pomikov, hitrosti in pospeškov, dolocimo v skladu z modificirano Newmarkovo shemo, kot smo pokazali v poglavju 4.2.2. Nove iteracijske vrednosti deformacijskih vektorjev dolocimo iz absolutnih zvez (6.15)-(6.16).
6.4.3 NumeriCna integracija
V vodilnih in pomoznih enacbah ter koeficientih tangentne matrike obeh formulacij koncnih elementov, ki smo jih predstavili, nastopajo integrali, ki so v splošnem prezahtevni, da bi jih lahko izracunali analiticno. Zato za racun integralov uporabimo numericno integracijo.
Integral poljubne funkcije h(x) vzdolz celotne dolzine elementa nadomestimo s kvadraturnim pravilom
T	N
h (x) dx _ £ wr h (xr),
0	r_1
kjer so wr utezi, xr pa integracijske tocke, ki jih izberemo glede na dolzino elementa L in izbrano metodo numericne integracije. Za integracijo vzdolz celotne dolzine elementa uporabimo Gaussova kvadraturna pravila.
Obe formulaciji koncnih elementov sta zasnovani tako, da potrebujemo vrednosti integralov tudi v dis-kretizacijskih tockah. Integral funkcije h (x) na intervalu [0,xq] izracunamo rekurzivno
xq	xq-1	xq
J h (X) dX _ j h (X) dX + J h (X) dX.
0	0	xq-1
Lokalno integracijo med dvema diskretizacijskima tockama zopet nadomestimo s kvadraturnim pravilom
xq	Nq
/ h (X) dX _ £ Wqs h (Xqs )
xq 1	qs_1
kjer so Xqs tocke lokalne integracije med dvema diskretizacijskima tockama, wqs pa ustrezne uteZi, odvisne od števila izbranih lokalnih tock Nq in dolzine intervala [Xq—1,Xq]. Tudi v tem primeru lokalne integracije za dolocitev tock in utezi izberemo Gaussova kvadraturna pravila. Integral Hq = /0Xq h (X) dX (H0 = 0) je potem enak
Nq
Hq = Hq-1 + £ Wqs h (Xqs ) .
qs=l
Pri racunu vgnezdenih integralov moramo integral funkcije h (x) izvrednotiti tudi v lokalni integracijski tocki Xq. Tudi v tem primeru uporabimo rekurzivno formulo
xqs	xq-1	xqs
J h (X) dX = J h (X) dX + J h (X) dX.
0	0	xq-1
Integral Hqs = h (X) dX med tocko globalne in lokalne integracije lahko zapišemo kot
xqs—1	xqs
Hqs = J h (X) dX + J h (X) dX.
xq—1	xqs— 1
Integral med dvema tocškama lokalne integracije tokrat izracšunamo s trapeznim pravilom
xqs
t	h(Xqs)+ h(Xqs — 1) ,	,
h (XX) dxx =	qs ^ qs	(Xqs — Xqs —1) .
xqs —1
Integral funkcije h(X) na intervalu [0, Xqs] je tako enak
/h (X)dX = H q—1 + H 1 + h(Xqs )+2h(Xqs—1 '(X,s — x„—1).
0
Predstavljena numericna integracija se je izkazala za natancno in racunsko ucinkovito.
6.5 Numericno dušenje
Dinamicni odziv konstrukcij je sestavljen iz neskoncnega števila nihajnih oblik, ki razlicno prispevajo k celotnemu odzivu. Praviloma je vpliv nihajnih oblik z višjimi frekvencami manjši od tistih z nizjimi frekvencami, ne pa nujno. Konstrukcije modeliramo z diskretnimi matematicnimi modeli, ki imajo koncno število prostostnih stopenj, s tem pa tudi koncno število nihajnih oblik. Mreza koncnih elementov, s katero diskretiziramo konstrukcijo, pa praviloma ne zadošca za opis kompleksnih višjih nihajnih oblik. Pri casovni diskretizaciji dinamicnega odziva uporabljamo koncno velike casovne korake, s katerimi zanemarimo vpliv nihajnih oblik, ki imajo manjši nihajni cas od izbranega koraka.
Tako casovna kot tudi krajevna diskretizacija konstrukcij vplivata na pomanjkljiv opis nihajnih oblik konstrukcije, ki ga v numericni analizi zaznamo kot numericno motnjo ali šum, ki vpliva na natancnost in stabilnost reševanja. Numericne motnje, ki so posledica vpliva višjih nihajnih oblik, izlocimo z umetno povzrocšenim numericšnim dusšenjem.
V literaturi lahko najdemo veliko število metod z numericnim dušenjem. V primeru analize linearnih dinamicnih sistemov sta taki tudi Newmarkova shema in modificirana Newmarkova shema, ce izberemo vrednost parametra 7 > 2 (tocki 4.1.3, 4.2.3). Taka izbira ni optimalna, saj vpliva na vse frekvence, z numericnim dušenjem pa zelimo vplivati le na višje frekvence. Za obravnavo linearnih dinamicnih sistemov so bile zato razvite metode, ki poskušajo doseci optimalno razmerje med disipacijo visokih in nizkih frekvenc, taki metodi sta na primer HHT metoda Hilbera, Hughesa in Taylorja [23] in WBZ-a metoda Wooda, Bossaka in Zienkiewicza [55]. Čhung in Hulbert [24] sta predstavila posplošeno a-metodo, v kateri kot posebna primera nastopata HHT in WBZ-a metodi. Posplošena-a metoda in njena posebna primera se pogosto uporabljajo tudi v nelinearni dinamiki kljub dejstvu, da metode niso brezpogojno stabilne in v splošnem ne zagotavljajo disipacije energije [56].
Posplošeno a-metodo uporabimo za numericno dušenje v formulaciji koncnega elementa z interpolacijo deformacij v primerih, ko rešitev ni mogoce dobiti z Newmarkovo shemo. Posplošeno a-metodo izberemo zaradi njene razširjenosti in dokaj preproste implementacije, ki temelji na Newmarkovi shemi. Osnovna ideja metode je, da dinamicno ravnotezje pri casu tn+1 nadomestimo z ravnotezjem pri nekem vmesnem casu tn < tn + {h < tn+1, tako da notranje in zunanje ravnotezne sile in momente izvrednotimo pri casu tn+1-af, vztrajnostne sile in momente pa pri casu tn+1-am. a/ in am sta parametra metode, ki dolocata vmesne case za izracun sil in momentov:
W1-af = (1 - a/) W1 + a/
^n+1-am — (1 am) ^n+1 + amtn.
Ravnotezne sile in momente pri vmesnih casih izracunamo kot linearno kombinacijo vrednosti pri casih
tn in tra+1:
N g
i[n+1- a f
e[n+1-af N g
n m[ra+1-am N g
i[n+1- a f
M g
M g[n+1-af
M m[n+1-am
—	(1 - a/) Ng[n+11+ a/Ng[n1
—	(1 - a/) Ng[n+11+ a/Ng[n1
—	(1 - am) Nm[n+11+ amNm[n1
—	(1 - a/) Mg[n+11+ a/Mg[n1
—	(1 - a/) Mg[n+11+ a/Mg[n1 ^ —(1 - am) Mm[n+11+ amM;[n1
Tako kombiniranje sil sta za implementacijo poplošene a-metode uporabila tudi Kuhl in Crisfield [57] pri implementaciji posplošene a-sheme na dinamicne enacbe nosilcev, povprecenje sil pa je tudi sicer obicajen pristop energijskih metod.
Diskretizirane vodilne enacbe elementa z interpolacijo deformacij (6.49)-(6.50) in (6.41)-(6.46), na katerih implementiramo posplosšeno a-metodo, imajo potem obliko
f 1, — R[n+11 ) NGln+111 ) - Ng[n+1-af1 ) - Ng^^1 (xg)
-	n^1-^1 ) — 0
f 2, — R[n+11 ) MGK11 ) - Mg[n+1-af 1 ) - Mg^^1 (xg)
-	Mm[n+1-am1 )— 0
f 3 =	- rgra+1] (L) = 0
f 4=^L[n+1]- 4n+1] (l) = o
f 5 = s0[n+1-af] + N[n+1-af] (0)= o
f 6 = p0[n+1-«f ] + m[n+1-a/] (0)= o
f 7 = sL[n+1-«/] - Ng[n+1-af] (L) - Ng[n+1-a/] (L) - N^1-^ (L) = 0
f 8 = PL[n+1-«/] - Mg[n+1-af] (L) - Mg[n+1-a/] (L) - M^1-^ (L) = 0.
Optimalno delovanje numericnega dušenja, pri katerem doseZemo dušenje visokih frekvenc in minimi-ziramo vpliv dušenja na nizke frekvence, je pogojeno z izbiro parametrov integracijske sheme ß, 7 in parametrov metode numericnega dušenja a/ in am. Formule za dolocitev parametrov posplošene-a sheme pri analizi linearnih dinamicnih sistemov sta podala Chung in Hulbert [24]
2P~ - 1	P~ n 1 U	\ 2	1	,, 1r>1.
am =-—, a/ =-—, ß = - (1 - am + a A , 7 = - - am + a/.	(6.101)
p» +1	p» +1	4 v	y	2
Parametri so zapisani v odvisnosti od vrednosti spektralnega radija p» G [0,1] tangentne matrike dis-kretiziranega sistema enacb v limiti, ko gre največja realna lastna vrednost matrike proti neskoncno. V linearni dinamiki tak izbor parametrov zagotavlja brezpogojno stabilnost metode, v nelinearni dinamiki pa ne. Formule (6.101) kljub temu optimirajo delovanje numericnega dušenja, zato jih uporabljamo za dolocitev parametrov sheme tudi v nelinearni dinamiki. Parameter p» je dobra mera za numericno dis-ipacijo, saj neposredno vpliva na njeno velikost. Manjša vrednost parametra pomeni vecjo disipacijo in obratno, vecja vrednost parametra pomeni manjše numericno dušenje. Zato v kontekstu posplošene-a sheme parameter p» vecškrat imenujemo tudi parameter numericšnega dusšenja.
Povsem analogno bi posplošeno a-metodo lahko implementirali v formulacijo koncnega elementa z interpolacijo krajevnih odvodov hitrosti. Poglaviten namen implementacije dušenja po a- metodi je izboljšanje numericne stabilnosti racunskega postopka. V numericnih primerih bomo pokazali, da z uporabo formulacije z interpolacijo krajevnih odvodov hitrosti in hkratno uporabo modificirane Newmar-kove sheme povecamo numericno stabilnost algoritma, zato se v tem primeru z implementacijo dušenja ne ukvarjamo.
6.6 Numericni testi
Predstavili smo dve formulaciji koncnih elementov za dinamicno analizo prostorskih nosilcev. V prvi formulaciji smo za osnovne neznanke diskretizacije po kraju izbrali deformacijska vektorja 7 in k, za cšasovno integracijo enacšb pa smo prilagodili Newmarkovo shemo. Formulacijo smo dodatno opremili s posplošeno-a metodo za numericno dušenje. Druga formulacija je osnovana na izbiri krajevnih odvodov hitrosti v' in kotnih hitrosti d za osnovne neznanke, za casovno integracijo pa je izbrana modificirana Newmarkova shema.
Obe formulaciji koncnih elementov smo vgradili v racunalniški program, zapisan v programskem paketu Matlab. Program omogoca analizo prostorskih linijskih konstrukcij, ki so podvrzene dinamicnim vplivom, s poljubno, tudi ukrivljeno in zvito, zacetno geometrijo in upoštevanje razlicnih konstitucijskih
zakonov. Za interpolacijo osnovnih neznank v obeh formulacijah koncnih elementov uporabljamo La-grangeve polinome. Koncni elementi so zasnovani tako, da so interpolacijske in diskretizacijske tocke racuna poenotene
{xp,p = 1,...,N} = {xq ,q = 1,...,N} ,
s cimer pridobimo na ucinkovitosti zapisanih algoritmov. Sistem nelinearnih algebrajskih enacb v vsakem casovnem koraku rešujemo z Newtonovo iteracijsko metodo. Za oceno natancnosti rezultatov v vsakem koraku Newtonove iteracije uporabimo Evklidsko normo vektorja neznank ||%|| in vektorja desnih strani ||f ||. Newtonovo iteracijo prekinemo, ko sta obe normi manjši od predpisane tolerance: ||%|| < £y in ||f || < £f. V numericnih študijah, ki jih predstavljamo v nadaljevanju, smo za toleranci uporabili vrednosti: ey = 10—10 in £f = 10—8.
Z numericnimi študijami zelimo pokazati lastnosti in delovanje predstavljenih formulacij koncnih elementov in casovno integracijskih shem. Rezultate v prvem in zadnjem primeru primerjamo z analiticnimi rešitvami, v preostalih primerih pa z rešitvami drugih avtorjev.
V vseh numericnih primerih privzamemo linearno elasticne konstitucijske zveze
N G
M G
EAi 0 0 0 GA2 0 0 0 GA3
GIi 0 0 0 EI2 0 0 0 EI3
Y G
kg.
E in G sta elasticšni in strizšni modul materiala, A1 je povrsšina precšnega prereza, A2 in A3 pa sta povrsšini striznih prerezov v smereh glavnih vztrajnostnih osi prereza G2 in G3. Ii je torzijski vztrajnostni moment prereza, I2 in I3 pa glavna vztrajnostna momenta prereza okrog osi G2 in G3.
Za parametre Newmarkove in modificirane Newmarkove sheme izberemo vrednosti, priporocene za ana-lizno nelinearnih dinamicnih sistemov, o katerih smo pisali v tockah 4.1.3 in 4.2.3. V primerih, v katerih uporabimo numericno dušenje, parametre posplošene-a metode dolocimo glede na izbrano vrednost spektralnega radija p«, kot to dolocajo enacbe (6.101). Pregled parametrov casovno integracijskih shem, ki jih uporabljamo v numericnih primerih, je podan v preglednici 6.1.
Preglednica 6.1: Parametri casovno integracijskih shem uporabljeni v numericnih primerih.
Table 6.1: Parameters of the time integration schemes used in numerical examples.
	Element z interpolacijo deformacij	Element z interpolacijo krajevnih odvodov hitrosti
Newmarkova shema Posplošena-a shema		Modificirana Newmarkova shema
ß=4	2p« — 11 , n 2 am = ß = (1 am + af) p«, +1 4 '	ß=2
1 7 = 2	p« 1 af = p« + 1 7 = 2 ^ + af	1 7 = 2
6.6.1 Vsiljeno nihanje prostoleZeCega nosilca
Obravnavamo primer prostolezecega nosilca, prikazanega na sliki 6.1. Nosilec je izpostavljen vplivu harmonicne sile FZ = F0sin w0t, ki deluje na sredini nosilca.
Z;
V literaturi [58] lahko poišcemo analiticne rešitve za tak problem ob predpostavki, da so pomiki in zasuki dovolj majhni, da lahko upoštevamo linearizirano teorijo nosilcev. Rešitev za vertikalne pomike nosilca je podana z neskoncšno vrsto
Formulaciji koncnih elementov, predstavljeni v tockah 6.3.1 in 6.3.2, sta osnovani na geometrijsko tocni teoriji nosilcev. Kljub temu morata pri amplitudah sile vzbujanja, pri katerih nastanejo majhni pomiki, dati rezultate, primerljive z analiticnimi rešitvami, ki veljajo za linearizirano teorijo nosilcev. Primer izkoristimo za sštudij konvergence numericšne resšitve v odvisnosti od:
•	števila elementov, ne, izbranih za modeliranje nosilca,
•	števila interpolacijskih tock, N, v posameznem elementu in
•	izbranega casovnega koraka h.
Numericšne resšitve primerjamo med sabo in z analiticšnimi resšitvami.
Obravnavamo nosilec z dolzino L = 10 in maso pA = 1. Analiticna rešitev ne upošteva vpliva vztrajno-stnih karakteristik precnega prereza nosilca, zato jih v racunu zanemarimo: J1 = J2 = J3 = 0. Preostale materialne in geometrijske lastnosti so enake:
Slika 6.1: Vsiljeno nihanje prostolezecega nosilca. Figure 6.1: Forced vibration of simply supported beam.
u
kjer je wn lastna frekvenca n-te nihajne oblike
EA1 = GA2 = GA3 = 105, G/1 = E/2 = E/3 = 103.
Za amplitudo sile vzbujanja izberemo vrednost F0 = 0.3, za frekvenco vzbujanja pa w0 = 2. Maksimalen pomik, ki ga pri teh vrednostih obtezbe predvideva analiticna rešitev, nastopi pri casu t = 1.34 na sredini nosilca in je enak u (2) = 0.01076. Velikost maksimalnega pomika je priblizno tisocinka dolzine
nosilca, zato je primerjava pomikov, dobljenih z linearizirano in geometrijsko tocno teorijo nosilcev, smiselna.
-to
10
10
10-
A 10-
10
10-
h = 0.1 h = 0.05 h = 0.025
0 2 4 6 8
t
(a) Element z interpolacijo deformacij, N = 6, ne = 20 10-
I
-to
10
10-
-4
^ 10-
10-
0 2 4 6 8
t
(c) Element z interpolacijo deformacij, h = 0.05, ne = 20 10-
10
f 10-jC io
,-4
-6
10-
ne = 10 n = 20 = 40
0 2 4 6 8
t
(e) Element z interpolacijo deformacij, h = 0.05, N = 6
-to
10
10-
-4
A 10-
10-
h = 0.1 h = 0.05 h = 0.025
4
t
6
(b) Element z interpolacijo krajevnih odvodov hitrosti, N = 6, ne = 20
10-
I
-to
10
10-
-4
^ 10-
10-
4
t
(d) Element z interpolacijo krajevnih odvodov hitrosti, h = 0.05, ne = 20
10-
10
S
f 10-jC io
,-4
-6
10-
4
t
(f) Element z interpolacijo krajevnih odvodov hitrosti, h = 0.05, N = 6
Slika 6.2: Vsiljeno nihanje prostolezecega nosilca. Absolutne napaka pomika na sredini nosilca na intervalu t G [0,8] (h = casovni korak, ne = število elementov, N = število interpolacijskih tock). Figure 6.2: Forced vibration of simply supported beam. Absolute error of midspan deflection on the interval t G [0,8] (h = time step, ne = number of elements, N = number of interpolation points).
0
2
8
0
2
6
8
0
2
6
8
Analize opravimo s tremi casovnimi koraki, h _ 0.1,0.05,0.025, s tremi razlicnimi mrezami koncnih elementov, ne _ 10,20,40, in tremi stopnjami interpolacije, N _ 2,6,10. Numericne rešitve primerjamo z analiticnimi na casovnem intervalu t G [0,8]. Na sliki 6.2 prikazujemo absolutne napake vertikalnih pomikov na sredini nosilca. Kratkorocni odziv nosilca je dobro opisan ze z redkejšo mrezo koncnih elementov in elementi nizjega reda. Napake se s casom seštevajo, največji vpliv na zmanjšanje napak pa ima izbira manjših casovnih korakov. Del napak je posledica dejstva, da primerjamo rešitve dobljene z razlicnima teorijama nosilce. Zato povecanje števila elementov in števila interpolacijskih tock ne vpliva bistveno na velikost napak, ki so enakega velikostnega reda za obe formulaciji koncnih elementov.
6.6.2 Konzola upognjena v spiralo
Deformacijski elementi, razviti za staticno analizo [9], odlicno delujejo pri velikih pomikih in, kar je še bolj pomembno, pri poljubno velikih prostorskih rotacijah. Tako obnašanje pricakujemo tudi od novega elementa z interpoliranimi krajevnimi odvodi hitrosti in kotnih hitrosti.
Slika 6.3: Konzola upognjena v spiralo. Figure 6.3: Cantilever bent to a helical form.
Element preizkusimo na primeru konzole, ki je na prostem koncu obremenjena z upogibnim momentom My in silo FY, ki konzolo potiska iz ravnine . Zacetna lega in obtezba konzole sta prikazani na sliki 6.3.
Dolzina konzole je L _ 10, preostali geometrijski in materialni podatki pa so:
EA1 _ GA2 _ GA3 _ 104 G11 _ E12 _ E13 _ 102.
Primer konzole na sliki 6.3 s staticno obtezbo je prvi predstavil Ibrahimbegovic [38]. Rezultate staticne analize pa sta predstavila tudi Zupan in Saje [9]. V naši analizi dinamicne vplive zanemarimo tako, da porazdeljeni masi elementov in komponentam vztrajnostnega tenzorja pripisšemo majhne vrednosti: pA _ 10-15, J1 _ J2 _ J3 _ 10-15.
Konzolo modeliramo z razlicnimi mrezami koncnih elementov in razlicnimi stopnjami interpolacije osnovnih spremenljivk. Za vse analize uporabimo casovni korak h _ 0.01. V preglednici 6.2 primerjamo pomike konca konzole v smeri obtezbe pri casu t _ 10, dobljene z obema formulacijama koncnih elementov. Opazimo, daje konvergenca pomikov dobljenih z elementi z interpolacijo krajevnih odvodov hitrosti celo nekoliko hitrejša od konvergence deformacijskih elementov.
Preglednica 6.2: Pomik prostega krajišca konzole v smeri obteZbe pri casu t = 10. Table 6.2: Displacement of the free-end of the cantilever in the direction of load at time t = 10.
Element z interpolacijo krajevnih
Element z interpolacijo deformacij
odvodov hitrosti
		N=3	N = 5	N = 8	N=3	N = 5	N=8
ne	= 10				-0.04625	-0.07844	-0.07640
ne	= 25				-0.07711	-0.07634	-0.07642
ne	= 50	-0.03611	-0.04406	-0.05281	-0.07645	-0.07644	-0.07644
ne	= 100	-0.07042	-0.07159	-0.07298	-0.07644	-0.07644	-0.07644
ne	= 200	-0.07565	-0.07580	-0.07598			
ne = število elementov, N = stopnja interpolacije
Pri casu t = 10 komponenti obtezbe dosezeta vrednosti = 50 in MY = 200n, konzola je zvita v spiralo z desetimi obroci, pomik v smeri sile FY pa konvergira k vrednosti = -0.07644. Deformirana lega konzole, modelirane s 25 elementi z 8 tockovno interpolacijo hitrosti pri casu t = 10, je prikazana na sliki 6.4(a). Slika 6.4(b) prikazuje spreminjanje komponente pomika v smeri obtezbe s casom. Prikazani rezultati se zelo dobro ujemajo z rezultati predstavljenimi v literaturi [9,38]. S tem smo pokazali, da tudi nova formulacija z interpolacijo krajevnih odvodov hitrost zmore natancno opisati deformiranje konstrukcij, z izrazito velikimi rotacijami.
z
0 -0.1 -0.2
-0.3 -0.15
-0.08 0.15~0 Y
(a) Deformirana lega konzole zvite v spiralo pri casu t = 10
4 3 2
*
6 1
o
a
0
10
(b) Pomik prostega krajišca konzole v smeri obtezbe
Slika 6.4: Konzola upognjena v spiralo. Rezultati analize s 25 elementi z 8 tockovno interpolacijo
krajevnih odvodov hitrosti. Figure 6.4: Čantilever bent to a helical form. Results of the analysis with 25 elements with 8 point
interpolation of spatial derivatives of vlocities.
0
2
4
6
8
6.6.3 Kolenasta konzola
Obravnavamo kolenasto konzolo, predstavljeno na sliki 6.5. Primer sta prva predstavila Simo in Vu-Quoc [4], za primerjavo z novo razvitimi formulacijami linijskih elementov pa so ga povzeli tudi drugi
avtorji npr. Ibrahimbegovic in Al Mikdad [13]. Konzolaje med gibanjem izpostavljena velikim pomikom in rotacijam in je zato zelo primerna za študij tako prostorske kot tudi casovne diskretizacije.
Slika 6.5: Kolenasta konzola. Figure 6.5: Right-angle cantilever.
Geometrijske in materialne lastnosti konzole so:
L — 10, EA1 — GA2 — GA3 — 106, G/1 — E12 — E13 — 103, pA — 1, J1 — 20, J2 — J3 — 10.
Obtezba deluje v kolenu konzole in je usmerjena izven ravnine konzole. S casom se spreminja v obliki trikotne sheme, predstavljene na sliki 6.5. Po koncu delovanja obtezbe konzola prosto niha. Nihanje opazujemo do casa t — 30. Avtorja v [4] predstavita rezultate, dobljene za konzolo, modelirano s štirimi elementi - dva elementa za vsak nosilec in dvajsetimi elementi - deset elementov za vsak nosilec. Pomiki in inkrementalni zasuki so interpolirani s kvadratnimi polinomi. Simulacija je bila opravljena z Newmarkovim integratorjem s konstantnim casovnim korakom h — 0.25.
Eksperiment ponovimo z novima formulacijama koncnih elementov. Deformacije in hitrosti interpoli-ramo s kvadratnimi polinomi. Pri racunu z deformacijskim koncnim elementom izberemo Newmarkov integrator, pri racšunu s koncšnim elementom z interpolacijo hitrosti pa modificirano Newmarkovo shemo.
V obeh primerih izberemo konstanten casovni korak h — 0.25, simulacijo pa ponovimo tudi z manjšim korakom h — 0.1. Konzolo modeliramo s štirimi in dvajsetimi koncnimi elementi. Pomiki kolena in prostega konca konzole pri analizi s štirimi koncnimi elementi so prikazani na sliki 6.6. Racun z deformacijskim elementom z Newmarkovo shemo pri casu t — 17.75 divergira, racun z elementom z interpolacijo hitrosti pa je stabilen ves cas integracije. Na sliki 6.7 so predstavljeni pomiki kolena in prostega konca konzole iz ravnine, dobljeni z analizo z 20 koncšnimi elementi. V tem primeru je racšun z obema formulacijama stabilen ves opazovani interval. Rezultati, dobljeni z gosto mrezo, se dobro ujemajo za obe formulaciji koncnih elementov in z rezultati predstavljenimi v [4]. Pomiki, ki jih izracunamo z redkejšo mrezo koncnih elementov z interpolacijo hitrosti, se razlikujejo od pomikov, ki sta jih za mrezo štirih elementov pokazala avtorja v [4], zelo dobro pa se ujemajo z rezultati, dobljenimi z mrezo dvajsetih koncšnih elementov.
Pomiki kolena in prostega konca konzole pri natancnejši analizi z mrezo 20 koncnih elementov in casovnim korakom h — 0.1 so prikazani na sliki 6.8. Z obema formulacijama dobimo povsem enake
koleno: Simo in Vu-Quoc, [4]
—	koleno: interpolacija deformacij prosti konec: Simo in Vu-Quoc, [4]
—	prosti konec: interpolacija deformacij
10 i-■-■-■-■-■-
* 2
0 5 10 15 20 25 30
t
(a) Newmarkova shema
koleno: Simo in Vu-Quoc, [4]
—	koleno: interpolacija hitrosti prosti konec: Simo in Vu-Quoc, [4]
—	prosti konec: interpolacija hitrosti 10
0 5 10 15 20 25 30
t
(b) modificirana Newmarkova shema
Slika 6.6: Kolenasta konzola. Pomik kolena in prostega konca konzole iz ravnine pri analizi s 4
elementi.
Figure 6.6: Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the cantilever obtained with analysis with 4 elements.
koleno: Simo in Vu-Quoc, [4] koleno: interpolacija deformacij prosti konec: Simo in Vu-Quoc, [4] ■ prosti konec: interpolacija deformacij
10 6
M 2
10 15 20 25 30
t
05
(a) Newmarkova shema
koleno: Simo in Vu-Quoc, [4] koleno: interpolacija hitrosti prosti konec: Simo in Vu-Quoc, [4] ■ prosti konec: interpolacija hitrosti
0 5 10 15 20 25 30
t
(b) modificirana Newmarkova shema
Slika 6.7: Kolenasta konzola. Pomik kolena in prostega konca konzole iz ravnine pri analizi z 20
elementi.
Figure 6.7: Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the cantilever obtained with analysis with 20 elements.
rezultate. Opazimo tudi, da se pomiki, dobljeni z manjšim casovnim korakom, v zadnji tretjini opazovanega casa razlikujejo od pomikov, dobljenih z vecjim korakom.
Na sliki 6.9 prikazujemo še casovni potek mehanske energije kolenaste konzole, modelirane z 20 elementi. V fazi prostega nihanja je energijsko stanje konzole stabilno pri obeh formulacijah koncnih ele-
10 6
* 2 6
a-2 -6 -10
	— prosti konec
/V f X	' \—koleno \ / t V \ \ / ' M ■
/ / s 1 / / / J J	\ I 1 \ \ 1 \ N \ - i X \ \/\ i N
	\\ r-'
	\ fc = 0.1
	\ ne = 20
0 5 10 15 20 25 30
t
	10
	6
*	2
S	
a	-2
	-6
	10
	— prosti konec
''v / X ' / * 1	' \—koleno \ / ' N \ \ 1 M \ / ' vi \ 1 1 \
t //	\ V 1 \\ r-'
	\ fc = 0.1
	\ ne = 20
5 10 15 20 25 30
t
(a) Element z interpolacijo deformacij, Newmarkova (b) Element z interpolacijo krajevnih odvodov hitrosti,
shema
modificirana Newmarkova shema
Slika 6.8: Kolenasta konzola. Pomiki kolena in prostega konca konzole iz ravnine, pri analizi s
korakom h = 0.1.
Figure 6.8: Right-angle cantilever. Out-of-plane displacements of the elbow and the free-end of the
cantilever, obtained with the time step h = 0.1.
0
mentov. Pri deformacijski formulaciji z Newmarkovo shemo v drugi polovici analize opazimo vecje oscilacije energije okrog konstantne vrednosti, medtem ko so oscilacije energije pri formulaciji z modificiranim Newmarkovim integratorjem skoraj neopazne tudi za večji casovni korak. V tem smislu lahko govorimo o povecani stabilnosti casovne integracije s formulacije koncnih elementov z modificirano shemo v primerjavi s formulacijo z Newmarkovo shemo.
120 « 100
80
u Ö u
ja 60
ZA
Ü 40
Jh
u
S 20
00
h = 0.25 — h = 0.1
ne = 20
5 10 15 20 25 30
t
(a) Newmarkova shema
120 « 100
s?
c
u tf M
m Š
80 60 40
S 20
00
h = 0.25 — h = 0.1
10
15
t
ne = 20
20 25 30
(b) Modificirana Newmarkova shema
Slika 6.9: Kolenasta konzola. Mehanska energija. Figure 6.9: Right-angle cantilever. Mechanical energy.
5
6.6.4 Nepodprti podajni nosilec v prostem letu
Opazujemo gibanje nepodprtega podajnega nosilca, ki se med gibanjem mocno deformira in doseze izjemno velike pomike in rotacije. Primer sta prva obravnavala Simo in Vu-Quoc [4], kasneje pa tudi Simo s sodelavci [10] in mnogi drugi. Primer je dober test za obnašanje linijskih elementov pri velikih
pomikih in rotacijah, predvsem pa je dober test za stabilnost casovnih integratorjev.
Slika 6.10: Nepodprti podajni nosilec v prostem letu. Figure 6.10: Flexible beam in free flight.
10 8 6
N 4 2 0 -2
10 8 6 4 2 0 -2
	. t = 3		5. 8
V^=2		t = 5\	t = 6.1
t=0\\	J t = 3.8-V		/-
-	/ t = 4.4......'		^ t = 7"".....-"' _
- / >			
10
15
20
25
10 8 6 4 2 0 -2
X
6 4 2 0 -2 -4 -6
r
(a) Element z interpolacijo deformacij, projekcija na ravnini XZ in YZ
10
	. t = 3		5. 8
V^=2		t = 5\	t = 6.1 "
t=0\\			-
/ >	/ t = 4.4......'		. t = 7"......_
10	15
X
20
25
8 6 4 2 0 -2
. t = 0 t = 3.5/ "tA.	^t = 2.5
t = 4.5[	
6 4 2 0 -2 -4 -6 Y
(b) Element z interpolacijo krajevnih odvodov hitrosti, projekcija na ravnini XZ in YZ
Slika 6.11: Podajni nosilec v prostem letu. Projekcije deformiranih leg. Figure 6.11: Flexible beam in free flight. Side views of deformed shapes.
0
5
0
5
Raven nosilec dolzine L = 10 je v zacetni legi postavljen v ravnino , koordinati robnih tock tezišcne osi nosilca sta enaki (6,0,0) in (0,0,8). Na spodnji rob deluje tockovna sila in par momentov MY in M z . Zacetna lega nosilca in spreminjanje obtezbe s casom sta predstavljena na sliki 6.10.
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(a) Element z interpolacijo deformacij, posplošena-a shema: _ 0.5, h _ 0.1
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(c) Element z interpolacijo deformacij, posplošena-a shema: _ 0.9, h _ 0.1
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(e) Element z interpolacijo krajevnih odvodov hitrosti, modificirana Newmarkova shema: h _ 0.1
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(b) Element z interpolacijo deformacij, posplošena-a shema: _ 0.5, h _ 0.01
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(d) Element z interpolacijo deformacij, posplošena-a shema: px _ 0.9, h _ 0.01
6 4 2
M
6 o
o a
-2
20
40 60
t
80 100
(f) Element z interpolacijo krajevnih odvodov hitrosti, modificirana Newmarkova shema: h _ 0.01
Slika 6.12: Podajni nosilec v prostem letu. Pomiki obtezenega vozlišca. Figure 6.12: Flexible beam in free flight. Displacements of the loaded end of the beam.
0
0
0
0
0
0
Nosilec modeliramo z 10 koncšnimi elementi, prvicš z linearno interpoliranimi deformacijami in drugicš z linearno interpoliranimi krajevnimi odvodi hitrosti. Izberemo casovni korak, kije priporocen v litereturi: h _ 0.1. Za analizo gibanja z deformacijskim elementom izberemo Newmarkovo shemo, analizo z elementom z interpolacijo hitrosti pa izvajamo z modificirano Newmarkovo shemo.
V zacetku gibanja obe analizi dajeta enake rezultate. Na sliki 6.11 so prikazane projekcije deformiranih leg nosilca na ravnini XZ in YZ. Projekcije, dobljene z razlicnima elementoma, so identicne, dobro pa se ujemajo tudi s projekcijami predstavljenimi v literaturi [4, Fig. 17 in Fig. 19].
Analizo z Newmarkovo shemo priblizno ob casu t = 16 zaznamuje nenaden porast mehanske energije, ki vodi k odpovedi numericne analize, kot kaze slika 6.13(a). Porast energije se zgodi zaradi napak, ki izvirajo iz casovne diskretizacije z Newmarkovo shemo, posledicno pa pride do prevelikega vpliva visokih frekvenc na odziv konstrukcije. Pri analizi z modificirano Newmarkovo shemo teh tezav ne zaznamo, casovna integracija pa gladko tece do casa t = 100, pri katerem prenehamo z opazovanjem gibanja.
Analizo gibanja z deformacijskim elementom ponovimo z uporabo posplošene a-metode. Opazujemo vpliv izbire parametra p„ in velikost izbranega koraka h na odziv nosilca. Na sliki 6.12 primerjamo pomike obtezenega vozlišca nosilca dobljene s casovnima korakoma h = 0.1, h = 0.01 in vrednostma parametra p«, = 0.5 in p«, = 0.9. Predstavljena je tudi referencna rešitev, dobljena z modificirano Newmarkovo shemo. Vpliv numericnega dušenja na odziv konstrukcije je velik, uporaba prevelikega dušenja in prevelikega casovnega koraka pa se odraza v napacnem odzivu konstrukcije. Pomiki, dobljeni z modificirano Newmarkovo shemo, so enaki za oba izbrana cšasovna koraka.
900
:^890
jjf
§ 880 iS
Ü 870
s
| 860
850

0
20
40 60
t
80 100
(a) Element z interpolacijo deformacij, Newmarkova shema
9001-■-■-■-■-
890
jjf
§ 880
s
Ü 870
s
| 860
850
— fc = 0.01
20
40 60
t
80 100
(c) Element z interpolacijo deformacij, posplošena a-shema: p^ = 0.9
900
890
Is
§ 880
s
Ü 870
s
| 860
850
fc = 0.1 — fc = 0.01
20
40 60
t
80 100
(b) Element z interpolacijo deformacij, posplošena a-shema: px = 0.5
900
:^890
g
Ü 880
s
Ü 870
s
| 860
850
— fc = 0.01

20
40 60
t
80 100
(d) Element z interpolacijo krajevnih odvodov hitrosti, modificirana Newmarkova shema
Slika 6.13: Podajni nosilec v prostem letu. Mehanska energija. Figure 6.13: Flexible beam in free flight. Mechanical energy.
0
0
Na sliki 6.13 prikaZemo spreminjanje mehanske energije nosilca s casom za vse opravljene analize. Potrdimo ugotovitev [57], da numericna disipacija za poljubne vrednosti parametra p„ ne zagotavlja nujno tudi disipacijo mehanske energije, slika 6.13(c), to pa lahko popravimo z uporabo manjših casovnih korakov. Ponovno pa se pokaze prednost analize z modificirano Newmarkovo shemo, pri kateri se mehanska energija pri konservativnih problemih obnaša izrazito stabilno. Za vecje casovne korake lahko razberemo manjše oscilacije okrog konstantne vrednosti, ki pa ne vplivajo na dolgorocni casovni odziv.
6.6.5 Problem štiričlenskega mehanizma
Obravnavamo primer, ki gaje prvi predstavil Bauchau [32]. Mehanizem, predstavljen na sliki 6.14, je sestavljen iz treh nosilcev, ki v zacetni legi lezijo v ravnini . Nosilci so v tockah B in C med sabo povezani z enoosnima clenkoma ter podprti v tockah A in D. Obe podpori sta nepomicni in dovoljujeta rotacijo mehanizma okrog normale na ravnino . Os vrtenja enoosnega clenka v tocki B je v zacetni legi pravokotna na ravnino v kateri lezi mehanizem. Posebnost primera je, daje os vrtenja clenka v tocki C v zacetni legi zavrtena tako, da oklepa z normalo na ravnino mehanizma kot <p = 5° in lezi v ravnini vzporedni z ravnino . Z rotacijo clenka je modelirana zacetna nepopolnost mehanizma. Mehanizem se vrti s predpisano kotno hitrostjo okrog osi clenka v podpori A. Posledica nepopolnosti je prostorski odziv mehanizma, pri katerem nastopajo velike notranje sile.
C
V G2
*
Gi D
L
X
prečni prerez nosilca 1 in 2:
^Gv
b L
b
prečni prerez nosilca 3: aG?,
bG2
Slika 6.14: Stiriclenski mehanizem. Figure 6.14: Four bar mechanism.
Enoosne clenke modeliramo s kinematicnimi vezmi, ki jih v sistem enacb konstrukcije vpeljemo z metodo Lagrangejevih mnoziteljev. Vzemimo enoosni clenek, predstavljen na sliki 6.15, v katerem se stikata nosilca, oznacena z 1 in 2. Naj bo {Gr (t), Gn (t), Gb (t)} desnosucna ortonormirana baza pripeta na clenek, tako daje bazni vektor Gr usmerjen v smeri osi vrtenja clenka. V robno tocko nosilca 1 pripnemo bazo { Gj (t), (t), G^ (t)}, v robno tocko nosilca 2 pa bazo { G2 (t), (t), G2 (t)}, tako da bazni vektorji obeh baz sovpadajo z baznimi vektorji clenka:
Gr (t) = G1 (t) = G2 (t), Gn (t) = Gn (t) = Gn (t), Gb (t) = G1 (t) = G2 (t). Algebrajski vezni enacbi, s katerima je enoosni clenek enolicno opisan, zahtevata da sta normali Gn in
Slika 6.15: Enoosni clenek. Figure 6.15: Revolute joint.
G2 pravokotni na os vrtenja, Gj, v poljubni deformirani legi:
O = Gj (t) ■ Gn (t) = 0 O = Gj (t) ■ g2 (t) = 0.
Primer štiriclenskega mehanizma je obseZno dokumentiran na spletnem naslovu [59], kjer so prikazani rezultati, dobljeni s programskim paketom Dymore [60]. Materialne in geometrijske karakteristike mehanizma zato povzamemo iz reference [59]:
E = 207 ■ 109, v = 0.3, p = 7800, L = 0.12, 61 = 0.016, 62 = 0.008.
Isti vir navaja predpisano velikost kotne hitrosti = 0.6. Vsak od nosilcev je modeliran s štirimi koncnimi elementi s kubicno interpolacijo osnovnih neznank. Analizo opravimo z obema formulacijama, odziv mehanizma pa opazujemo do cšasa t = 12.
(a) Element z interpolacijo deformacij, posplošena-a	(b) Dymore [59]
shema: = 0
Slika 6.16: Stiriclenski mehanizem. Pomiki vozlišca C. Figure 6.16: Four bar mechanism. Displacements at point C.
1000 0
-1000 rt -2000 " -3000 -4000 -5000 -6000
		
-O	f -">>	r~
		— NG1 -
		— NG2 ■
		NG3
0
2
4
6 t
8 10
(a) Element z interp. def., posplošena-a shema, p»
140
iS 60 6
g 20
12
0
(c) Element z interp. def., posplošena-a shema, p»
12
0
1000 0
-1000 rt -2000 " -3000 -4000 -5000 -6000
140 100
iS 60 6
g 20
0
246
t
(b) Dymore [59]
8 10 12
12
(d) Dymore [59]
Slika 6.17: Stiriclenski mehanizem. Notranje sile in momenti na sredini prvega nosilca. Figure 6.17: Four bar mechanism. Internal forces and moments at the midpoint of the first beam.
o
a £
0.04 0.02 0
-0.02 -0.04
0
2
4
6 t
8 10 12
(a) Element z interp. def., posplošena-a shema: p» = 0 0.1
^ -0.1
m O
a
3 -0.3
a o
^ -0.5
-0.7
6 t
8 10
(c) Element z interp. def., posplosšena-a shema: p»
12
0
o
a £
0.04 0.02 0
-0.02 -0.04
0
o
^ -0.5
-0.7
246
t
(b) Dymore [59]
8 10 12
2 4 6
t
(d) Dymore [59]
8 10 12
Slika 6.18: Stiriclenski mehanizem. Hitrosti in kotne hitrosti na sredini prvega nosilca. Figure 6.18: Four bar mechanism. Velocities and angular velocities at the midpoint of the first beam.
0
2
4
0
Rezultati predstavljeni v viru [59] so pridobljeni s casovnim korakom h = 0.004 in z asimptoticnim numericnim dušenjem, pri katerem je vrednost parametra p«, = 0. Potreba po numericnem dušenju se pokaze tudi pri analizi s koncnimi elementi z interpolacijo deformacij, pri katerih uporabljamo Newmarkovo shemo. Zaradi motenj, ki so posledica zavrtenega clenka v vozlišcu C, pride do velikih oscilacij deformacijskih vektorjev in posledicne slabe konvergence numericne rešitve ne glede na izbran casovni korak. Pri velikosti koraka h = 0.004 racun preneha konvergirati pri casu t = 4.7. Analizo ponovimo s posplošeno-a shemo z vrednostjo parametra p«, = 0. Z uporabo dušenja postane racun stabilen ves opazovani casovni interval, ce uporabimo casovni korak manjši ali enak h = 0.004. Primerjavo rezultatov naše analize z rezultati, predstavljenimi v [59], prikazujemo na slikah 6.16, 6.17 in 6.18. Na sliki 6.16 so prikazani pomiki vozlišca Č v odvisnosti od casa. Primerjava notranjih sil in momentov na sredini nosilca z oznako 1 je prikazana na sliki 6.17, primerjava hitrosti in kotnih hitrosti v isti tocki pa je prikazana na sliki 6.18. Vse rešitve se zelo dobro ujemajo z referencnimi rešitvami [59].
Analizo opravimo še z elementi z interpolacijo krajevnih odvodov hitrosti in modificirano Newmarkovo shemo s casovnimi koraki h = 0.04, h = 0.004 in h = 0.0004. Kljub odsotnosti numericnega dušenja je analiza v vseh treh primerih stabilna ves opazovani cas. Rezultate analiz primerjamo na slikah 6.19, 6.20 in 6.21.
(a) h = 0.04	(b) h = 0.004
Slika 6.19: Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti -
pomiki vozlisšcša C.
Figure 6.19: Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of
velocities - displacements at point C.
Razlika med pomiki vozlišca C, ki jih dobimo pri racunu s tremi razlicnimi casovnimi koraki, je neopazna - na sliki 6.19 prikazani pomiki vozlišca C, dobljeni pri racunu s korakoma h = 0.04 in h = 0.004. Rešitve se tudi zelo dobro ujemajo z referencnimi rezultati prikazanimi na sliki 6.16(b). Vecje razlike v odzivu konstrukcije glede na izbran racunski korak opazimo na sliki 6.20, na kateri prikazujemo notranje sile in momente na sredini prvega elementa, in na sliki 6.21, na kateri so prikazane hitrosti in kotne hitrosti na sredini prvega elementa. V odzivu so lepo vidne oscilacije notranjih sil in predvsem oscilacije hitrosti, ki so posledica modeliranja zacetne nepopolnosti sistema z algebrajskimi veznimi enacbami. Z zmanjševanjem casovnega koraka se povecuje natancnosti racuna, zato se amplitude oscilacij znatno zmanjšajo. Pomembno je, da oscilacije ne vplivajo na stabilnost racuna, kar zopet prica o robustnosti elementov z interpolacijo krajevni odvodov hitrosti.
1000 0
-1000 « -2000 " -3000 -4000 -5000 -6000
0
1000 0
-1000 « -2000 " -3000 -4000 -5000 -6000
0
1000 0
-1000 « -2000 " -3000 -4000 -5000 -6000
J		■"'■• .HllfiMWi P.............
	............... 1	
		— NG1 '
		NG2 ■
		Ng3
246 t
(a) h = 0.04
246 t
(c) h = 0.004
8 10 12
■nGi
■nG2
nG3
8 10 12
		
--\		r"
		— Ni '
		— NG2 ■
		Ng3
0 2 4 6 8
t
(e) h = 0.0004
10 12
140 100
iS 60 6
I 20
140 100
iS 60 6
I 20
140 100
iS 60 6
I 20
246 t
(b) h = 0.04
246
t
(d) h = 0.004
8 10 12
8 10 12
2 4 6 8 10 12
t
(f) h = 0.0004
Slika 6.20: Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti -
notranje sile in momenti na sredini prvega nosilca. Figure 6.20: Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of velocities - internal forces and moments at the midpoint of the first beam.
o £
'JE
0.05 0.025
i 0
i
-0.025 -0.05
0.05 0.025 0
-0.025 -0.05
0.05 0.025
i 0 -0.025
0 2 4 6
(a) h= 0.04
8 10 12
0 2 4 6 8 10 12
t
(c) h = 0.004
-0.05
0 2 4 6
t
(e) h = 0.0004
8 10 12
0.4 0.2 0
m
o
£ -0.2 | -0.4 -0.6
-0.8 -1
A4
0
246 t
(b) h = 0.04
-0.
|-0 J3 -0 -0.
(d) h= 0.004
0 -0.2
13 -0.4
o M
-0.6 -0.8 -1
0
246 t
(f) h = 0.0004
8 10 12
10 12
10 12
Slika 6.21: Stiriclenski mehanizem. Analiza z elementi z interpolacijo krajevnih odvodov hitrosti -
hitrosti in kotne hitrosti na sredini prvega nosilca. Figure 6.21: Four bar mechanism. Analysis with elemens with interpolation of spatial derivatives of velocities - velocities and angular velocities at the midpoint of the first beam.
8
6.6.6 Enostavno nihalo
Številni avtorji za preizkus casovno integracijskih shem in formulacij linijskih koncnih elementov obravnavajo enostavna nihala. Posebno obcutljivi so primeri zelo togih nihal. Primer analize togega nihala z Newmarkovo casovno integracijsko shemo podrobno obravnava Bathe [41], za njim pa številni drugi avtorji. Za našo analizo izberemo togo nihalo, predstavljeno na sliki 6.22, ki so ga obravnavali Crisfield in Shi [61] ter Kuhl in Crisfield [57]. Ugotovitve avtorjev so, daje za casovno stabilno integracijo z Newmarkovo shemo potreben zelo majhen casovni korak. Pri analizi z večjimi casovnimi koraki pride do velikih oscilacij osnih deformacij nosilca, ki zaradi velike togosti nihala privedejo do porasta celotne energije nihala in napacšnih resšitev. Primer zato izkoristimo za oceno delovanja implementacije po-splošene-a metode za numericno dušenje, ki smo jo predstavili v tocki 6.5 tega poglavja. Obravnavamo
L = 3.0443 EA = 1010 pA = 6.57 v-0 =7.72
Slika 6.22: Enostavno nihalo. Figure 6.22: Simple Pendulum.
ravninsko gibanje nihala s porazdeljeno maso pA = 6.57, ki ga iz mirovanja pozenemo v gibanje z zacetno hitrostjo v tangentni smeri vf = 7.72. Nihalo modeliramo s štirimi elementi s po dvema interpolacijskima tockama. Analizo gibanja racunamo s tremi casovnimi koraki h = 0.1, h = 0.05 in h = 0.025, gibanje nihala pa opazujemo do casa t = 30. Racun izvedemo z Newmarkovo shemo, s posplošeno-a shemo z vrednostima parametra p„ = 0.8 in p„ = 0.6 in z modificirano Newmarkovo shemo. Ker je nihalo zelo togo, lahko numericšne resšitve primerjamo s tocšnimi resšitvami za popolnoma togo nihalo. Tocne rešitve za komponente krajevnega vektorja in vektorja hitrosti prostega konca nihala so enake
v0 t	v0 t
rx (L, t) = L sin^f-, rz (L, t) = —L cos^-,	(6.102)
LL
v0 t v0 t
vx (L, t) = vfcos vX_, vz (L, t) = vfsin -f-.	(6.103)
LL
Na sliki 6.23 prikazujemo celotno energijo nihala za vse opravljene analize. Racun s koncnimi elementi z interpoliranimi deformacijskimi vektorji je za vse tri izbrane casovne korake neuspešen, kar je v skladu s pricakovanji, saj formulacija uporablja Newmarkovo shemo. Po nekaj sekundah gibanja pride do velikega porasta celotne energije nihala, zato racun prekinemo. Za ta zahteven racunski primer tudi racun s formulacijo koncnih elementov z interpolacijo krajevnih odvodov hitrosti, ki uporablja modificirano Newmarkovo shemo, ni uspesšen. Tudi resševanje z numericšnim dusšenjem s posplosšeno-a shemo ne zagotavlja brezpogojne stabilnosti racuna. Pri analizi s korakom h = 0.1 in parametrom p„ = 0.8, slika 6.23(a), kljub numericnemu dušenju pride do porasta celotne energije nihala. Stabilnost racuna pri izbranem cšasovnem koraku lahko dosezšemo le s primerno izbranim parametrom numericšne disipacije.
(a) Čelotna energija nihala, h = 0.1
(b) Čelotna energija nihala, h = 0.05
(c) Čelotna energija nihala, h = 0.025
Slika 6.23: Enostavno nihalo. Čelotna energija nihala v odvisnosti od casa. Figure 6.23: Simple pendulum. Total energy of the pendulum versus time.
101
> o
& 1
g 10-1 o
& 10 &
§ 10-5
ü 10
10
20
30
101 o i
£ 10-1
si 10-3 j| 10-5 ^ 10-7,
10
20
30
(a) Element z interpolacijo deformacij, posplošena-a shema, px = 0.6, h = 0.1
101
> o
& 1
g 10-1 o
&10-3 &
§ 10-5
ü 10-
10
20
30
101
o
£ 10 si 10-3 j| 10-5 ü 10-7
10
20
30
(b) Element z interpolacijo deformacij, posplošena-a shema, px = 0.6, h = 0.05
101
>
o M
'Š 10 o
&10-3 &
§ 10-5
ü 10
10
20
30
101
10-i10-j 10-5 ü 10-7
10
20
30
(c) Element z interpolacijo deformacij, posplošena-a shema, px = 0.6, h = 0.025
101
> o
£ 1
g 10-1 o
&10-3 a
§ 10-5
Ü 10
-7
— uX
uz
10
20
30
101
ä 10 si 10
j| 10-5 ü 10-7
(d) Element z interpolacijo deformacij, Newmarkova shema, h =
t
0.001
30
10'
> o
g 10-1 o
&10-3 &
§ 10-5
Ü 10-
ux
uz
20
30
101
ä 10-
si 10
ca
10-5 ca 10-7
■Vx VZ
w^i^^^nnnnnnnnnnnnnnr
10
20
(e) Element z interpolacijo odvodov hitrosti, modificirana Newmarkova shema, h = 0.001
0
0
0
0
0
0
0
Slika 6.24: Enostavno nihalo. Absolutne napake pomikov in hitrosti. Figure 6.24: Simple pendulum. Absolute error of displacements and velocities.
O podobnem obnašanju posplošene-a sheme, ko jo uporabimo za reševanje problemov v nelinearni dinamiki, so porocali številni avtorji, na primer [56,62]. Naši rezultati za togo nihalo pa se ujemajo z rezultati predstavljenimi v clanku [57, slika 6]. Od tod sklepamo, daje implementacija posplošene-a sheme, ki jo predstavimo v tem delu, kvalitativno enakovredna formulacijam s posplošeno-a shemo, ki jih srecamo v literaturi.
Uporaba numericne disipacije ima lahko velik vpliv na natancnost dobljenih rezultatov. Na sliki 6.24 prikazujemo absolutne napake pomikov in hitrosti pri analizi s posplošeno-a shemo z p„ = 0.6 v primerjavi s tocnimi rešitvami (6.102)-(6.103). Z manjšanjem casovnih korakov se natancnost racuna s posplosšeno-a shemo povecšuje. Primerjavi dodamo tudi rezultate analiz z Newmarkovo in modificirano Newmarkovo shemo, dobljene s casovnim korakom h = 0.001, s katerim je bila casovna integracija stabilna skozi celoten opazovani interval. Zanimiva je primerjava natancnosti rezultatov, dobljenih z analizama brez numericnega dušenja. Opazimo, da je natancnost analize z modificirano Newmarkovo shemo za velikostni red vecja od natancnosti analize z Newmarkovo shemo.
7 GIBANJE MASNEGA DELCA IN NOSILCA
Dinamicen odziv konstrukcije je lahko posledica casovno spremenljivih obtezb, ki delujejo na konstrukcijo, lahko je posledica predpisanih hitrosti in pospeškov posameznih delov ali pa kar celotne konstrukcije. Primere takih vplivov na konstrukcije smo obravnavali v prejšnjem poglavju. K casovno spremenljivim vplivom pa prištevamo tudi vpliv drugih teles na konstrukcije, pa naj bo to trk dveh konstrukcij, gibanje vozil na konstrukciji, letea delci, ki udarijo ob konstrukcijo, in mnogi drugi.
Analiza in numericšno modeliranje tako splosšnih vplivov je mnogo preobsezšna tema za namen tega dela. Zato se omejimo na zunanje vplive, ki jih lahko zajamemo s poenostavljenim modelom masnega delca. Z masnim delcem imamo v mislih togo telo s konstantno maso, ki je tako majhno, da lahko zanemarimo njegovo vrtilno kolicino. Zanima pa nas predvsem dogajanje pri trku masnega delca ob nosilec.
V literaturi lahko poišcemo razlicne pristope k modeliranju trkov, zdruzimo pa jih lahko v dve vecji skupini glede na predpostavke o trajanju trka:
•	V prvo skupino prištevamo pristope, ki trk obravnavajo kot nezvezen dogodek. Trajanje trka je tako kratko, da lahko predpostavimo, da se konfiguracija teles v kontaktu med trkom ne spremeni, stiskanje in razmikanje teles ob trku pa prevzame fiktiven deformabilen delec na stiku obeh teles. Tak pristop je bil sprva uporabljen pri obravnavi trkov med togimi telesi [63], kasneje pa je bil prirejen tudi za trk proznih teles [64].
•	V drugo skupino prištevamo pristope, ki trk obravnavajo kot zvezen dogodek, ki ima koncno dolgo trajanje. Telesi se v obmocju kontakta lahko deformirata, deformacije pa so prek konstitutivnega zakona povezane s kontaktnimi silami.
Za modeliranje trka masnega delca in nosilca izberemo pristop s koncnim trajanjem trka, predvsem zato, ker omogoca izracun kontaktnih sil. Kontakt masnega delca in nosilca v sistem enacb vpeljemo prek kinematicne vezi, kot to predlaga Bauchau [34].
7.1 Osnovne enacbe
Obravnavamo povezan problem gibanja masnega delca in nosilca. Gibanje masnega delca vodi izrek o gibalni kolicini (4.40), iz katerega izpeljemo vodilno gibalno enacbo masnega delca:
f
delca
mX — £ F

0.
(7.1)
md je masa delca, ad je pospešek delca, Fd pa poljubna sila, ki deluje na delec. Vodilni sistem enacb povezanega sistema delec-nosilec dobimo tako, da gibalno enacbo delca (7.1) dodamo vodilnemu sistemu enacb konstrukcije f konst
f
sist
f fcorast fdelca
0.
(7.2)
Delec in nosilec sta lahko v kontaktu ze na zacetku gibanja, delec na primer drsi po nosilcu ali pa je trdno pritrjen na nosilec, do kontakta pa lahko pride med samim gibanjem - masni delec trci ob nosilec. Pogoje kontakta (stik, drsenje, trk) opišemo s kinematicnimi veznimi enacbami, ki jih formalno zapišemo kot
$ ^ rnosilec f nosilec fd f d ^ = q	(7 3)
Vezne enacbe vpeljemo v vodilni sistem enacb z uporabo metode Lagrangejevih mnoziteljev
f sist - 5$TX = 0	(7.4)
$ = 0.	(7.5)
Z smo oznacili Jacobijevo matriko, ki pripada kinematicnim vezem, X pa so Lagrangejevi mnozitelji. Za primeren opis vezi prepoznamo v Lagrangejevih mnoziteljih ravno kontaktne sile, tem pa v vecini primerov kontakta pripadajo konstitutivni zakoni, ki opisujejo lastnosti kontakta, na primer Čoulumbov zakon trenja za drsec kontakt ali pa konstitutivni zakon trka.
7.2 RaCunski primeri
7.2.1 Interpolacija osnovnih neznank nosilca s konstantnimi funkcijami
Trku masnega delca ob konstrukcijo obicajno sledi buren odziv konstrukcije, pri katerem zaradi mase delca pomembno vlogo igrajo tudi višje nihajne oblike. Za analizo takih dogodkov je zato nujna robustna diskretizacija konstrukcije po metodi koncnih elementov in uporaba kratkih casovnih korakov pri integraciji dinamicšnih enacšb po cšasu.
Pri formulaciji koncnih elementov z interpolacijo deformacij (6.3.1) in tudi pri formulaciji z interpolacijo krajevnih odvodov hitrosti (6.3.2) smo za interpolacijske funkcije izbrali Lagrangeve polinome. Pri analizah konstrukcij v stiku z masnimi delci smo pri obeh formulacijah koncnih elementov opazili velike oscilacije interpoliranih kolicin, predvsem v koncnih elementih, na katere vpliva masni delec. Posledicno je bila konvergenca Newtonovih iteracij pri reševanju enacb pocasnejša, za uspešen racun (tudi z numericnim dušenjem) pa smo morali uporabiti gostejše mreze koncnih elementov z nizko stopnjo interpolacije. Podoben fenomen je znan iz staticne analize konstrukcij, pri katerih pride do lokalizacije deformacij, o cemer porocajo avtorji v [65,66]. Ena od rešitev problema je uporaba koncnih elementov s predpostavljenim konstantnim potekom interpoliranih kolicin. Tak pristop se je izkazal kot zelo uspešen, saj omogoca racunsko zelo ucinkovito formulacijo. Linijski koncni element s predpostavljenim konstantnim potekom deformacij za staticno analizo konstrukcij je podrobno predstavljen v clanku [66]. Podoben pristop izberemo tudi za resševanje povezanega problema gibanja masnega delca in nosilca in pripravimo koncne elemente za dinamicno analizo s predpostavljenim konstantnim potekom osnovnih neznank.
Pri formulaciji koncnega elementa z interpolacijo deformacij interpolacijske nastavke deformacij (6.47)-(6.48) nadomestimo s konstantnim potekom deformacijskih vektorjev:
Y [nt%{x)= Y G+l	(7.6)
^[++1 (X) = KG++1 .	(7.7)
Diskretiziranim vodilnim enacbam (6.49)-(6.50) pa zaradi konstantnosti deformacij zadostimo na sredini elementa
f 1 = R[n+1] (2) NG[n+1] (2) - Nin+1] (2) = 0 f 2 = R[n+1] (2) MG[n+1] (2) - Mf+1] (2) = 0.
(7.8)
(7.9)
Enako storimo tudi pri elementu z interpolacijo krajevnih odvodov hitrosti. Interpolacijske nastavke za krajevne odvode hitrosti (6.59)-(6.60) nadomestimo z
vf+1] (x) = vf+1]
* Gn+1] (x) = * Gn+1],
vodilnim enacbam (6.61)-(6.61) pa tudi v tem primeru zadostimo na sredini elementa
0 0.
h1 = R[n+1] (2) NG[n+1] (2) - Nf+1] (2) h2 = R[n+1] (2) MG[n+1] (2) - Min+1] (2)
(7.10)
(7.11)
(7.12)
(7.13)
7.2.2 Nihajoča gugalnica
Obravnavamo gugalnico, prikazano na sliki 7.1. Gugalnica je sestavljena iz togih veznih palic in nosilca, na katerem je pripet masni delec. Racunski primer je prvi predstavil Bauchau s sodelavci [11,67], kasneje pa so ga obravnavali številni avtorji, na primer [27,68]. Posebnost primera so toge vezi in masni delec, ki ima izredno veliko maso v primerjavi z maso nosilca.
z
	r
L	
-1-	
FxO 2
0
0 0.128 0.256 t
prerez:
b
G2
bi
Slika 7.1: Nihajoca gugalnica. Figure 7.1: The swing pendulum.
Na nosilec, ki je clenkasto vpet med togi palici, je v tocki A na sredini nosilca togo pripet delec z maso m = 0.5. Materialne in geometrijske lastnosti nosilca so enake
E = 73 ■ 109, G = 28.08 ■ 109, p = 2700, 61 = 0.005, 62 = 0.001, 2 = 0.72.
PreCni prerez togih palic je enak precnemu prerezu nosilca. Maso palic zanemarimo, togost pa modeliramo tako, da elasticnemu modulu palic pripišemo stokrat vecjo vrednost od elasticnega modula nosilca.
Nihalo v gibanje pozene tockovna sila FX s prijemališcem v tocki A, ki se s casom spreminja v obliki trikotnega pulza z maksimalno vrednostjo = 2 pri casu t = 0.128. Pri casu t = 0.256 je vrednost sile enaka nic, sistem pa od tega trenutka dalje prosto niha. Nihalo se v zacetku gibanja vrti v nasprotni smeri urinega kazalca. Priblizno pri casu t = 0.64 sta leva palica in nosilec povsem poravnana, desna palica pa se zacne gibati v smeri urinega kazalca. To povzroci nenadno spremembo v smeri hitrosti gibanja masnega delca, kar zaradi velike mase delca na gibanje sistema ucšinkuje podobno kot trk masnega delca. Primer je tako dober test za robustnost krajevne in cšasovne diskretizacije vodilnih enacšb.
Ker je masni delec pritrjen na konstrukcijo, njegov vpliv upoštevamo kot negativno obtezbo v tocki A: Fm = ma^. Stik delca in nosilc masnega delca in tocke A enaka
Fm = ma^. Stik delca in nosilca pa bi lahko upoštevali tudi z vezno enacbo, ki zahteva, da sta pomika
$=rA (t) - rm (t)=o.
Za prostorsko diskretizacijo nihala izberemo koncne elemente s konstantnimi deformacijami in konstantnimi krajevnimi odvodi hitrosti. Vsako togo palico modeliramo s štirimi koncnimi elementi, podajni nosilec pa z 52 elementi. Gibanje nosilca opazujemo do casa t = 2.5.
0.8

0.6
0.4
0.2
-0.2
-0.2 0	0.2 0.4 0.6 0.8
X
1.2
Slika 7.2: Nihajoca gugalnica. Zaporedje deformiranih leg. Figure 7.2: The swing pendulum. Sequence of deformed shapes of the swing.
0
1
Zacetne analize opravimo s casovnim korakom, ki ga priporoca literatura h = 0.0005. Analiza z elementi, ki uporabljajo Newmarkovo shemo se predcasno zakljuci zaradi pojava divergence pri casu t = 0.76, analiza z modificirano Newmarkovo shemo pa pri casu t = 0.86. V obeh primerih je analiza neuspešna zaradi velikega porasta celotne energije sistema. Če uporabimo manjši casovni korak h = 0.0001, opa-
zimo podobno obnašanje. Analizo z Newmarkovo shemo je uspešna do casa t = 0.93, analiza z modificirano Newmarkovo shemo pa do casa t = 1.56. Stabilnost racuna z Newmarkovo shemo se tudi pri racunih s krajšimi casovnimi koraki ne izboljša, zato racun ponovimo še s posplošeno-a shemo z dvema razlicnima vrednostima parametra numericnega dušenja p«, = 0.9 in p«, = 0.5. Uporaba numericnega dušenja stabilizira racun, ki ga uspešno zakljucimo pri casu t = 2.5. Pri analizi z elementi s konstantnimi krajevnimi odvodi hitrosti, ki uporabljajo modificirano Newmarkovo shemo, pa opazimo, da se z manjšanjem racunskega koraka stabilnost povecuje. Analizo pa uspešno zakljucimo s casovnim korakom h = 5 ■ 10-5. Rezultate analiz predstavljamo na slikah 7.2-7.6.
(a) h = 0.0005	(b) h = 0.0001
Slika 7.3: Nihajoca gugalnica. Pomiki v tocki B, izracunani z elementom s konstantnimi
deformacijami in posplošeno-a shemo. Figure 7.3: The swing pendulum. Displacements at point B, obtained by constant-strain elements and
Generalized-a scheme.
0.01 0
0 0.5 1 1.5 2 2.5
t
(a) h= 0.0005
0.01 0
0 0.5 1 1.5 2 2.5
t
(b) h= 0.0001
Slika 7.4: Nihajoca gugalnica. Čelotna energija masnega delca in gugalnice, diskretizirane z elementi s
konstantnimi deformacijami. Figure 7.4: The swing pendulum. Total energy of the mass and swing, discretized with constant strain
elements.
Slika 7.2 prikazuje zaporedje deformiranih leg gugalnice na intervalu t G [0,2.5], ki jih dobimo pri analizi z modificirano Newmarkovo shemo s korakom h = 5 ■ 10-5. Na sliki 7.3 primerjamo pomike v tocki B,
M
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2
0 0.5 1 1.5 2 2.5
t
(a) Pomiki tocke B
0.01 0
(b) Čelotna energija sistema
Slika 7.5: Nihajoca gugalnica. Pomiki tocke B in celotna energija sistema pri analizi z elementi s
konstantnimi krajevnimi odvodi hitrosti in korakom h = 5 ■ 10-5. Figure 7.5: The swing pendulum. Displacements of point B and total energy of the system obtained with elements with constant spatial derivatives of velocities and time step h = 10-5.
150
100
ca 50 ö
-50
0 0.5 1 1.5
t
2.5
(a) Element s konstantnimi deformacijami, posplošena-a shema, px = 0.9, h = 0.0001
150
100
ca 50 ö
-50
0 0.5 1 1.5
t
(c) Element s konstantnimi krajevnimi odvodi hitrosti, modificirana Newmarkova shema h = 0. 00005
0 0.5 1 1.5
t
(b) Element s konstantnimi deformacijami, posplosšena-a shema, p^ = 0.9, h = 0.0001
5
0 0.5 1 1.5 2 2.5
t
(d) Element s konstantnimi krajevnimi odvodi hitrosti, modificirana Newmarkova shema h = 0. 00005
0
0
Slika 7.6: Nihajoca gugalnica. Notranje sile v nosilcu, levo od tocke A. Figure 7.6: The swing pendulum. Internal forces on the left-hand side of point A.
dobljene z analizami z elementi s konstantnimi deformacijami in posplošeno-a shemo, na sliki 7.4 pa pripadajoco celotno energijo nihala. Pri casu t = 0.64, ko se spremeni smer gibanja desne palice, kar povzrocši visokofrekvencšne oscilacije v odzivu nosilca, opazimo velik padec celotne energije sistema. Disipacija energije se nato umiri in ostane majhna do konca analize. Pri pomikih opazimo drugacen vpliv numericnega dušenja. V zacetnem delu gibanja, priblizno do casa t = 1, je ujemanje pomikov, dobljenih z razlicnimi vrednostmi parametra p„, kar dobro, prav tako se pomiki dobro ujemajo z rezultati predstavljenimi v [11,27]. V nadaljevanju gibanja pa razlike postanejo vecje.
Pomike tocke B in celotno energijo sistema pri analizi z elementi s konstantnim potekom krajevnih odvodov hitrosti prikazujemo na sliki 7.5. Pomiki na sliki 7.5(a) se do graficne natancnosti ujemajo s pomiki, ki jih v doktorski disertaciji pokaze Gams [69]. Čelotna energija je med prostim nihanjem gugalnice konstantna, slika 7.5(b), kar je posledica stabilnega racuna te metode. Na sliki 7.6 so prikazane osne in precšne sile v prvi diskretizacijski tocški v elementu levo od masnega delca. Opazimo velik vpliv numericnega dušenja na velikost sil, kljub racunu s krajšim casovnim korakom in nizko stopnjo dušenja. Sile izracunane pri analizi z elementi s konstantnim potekom krajevnih odvodov hitrosti se zelo dobro ujemajo z rezultati, ki jih za nedušeno nihanje predstavljajo avtorji v [69,70].
Elementi s predpostavljenim konstantnim potekom osnovnih neznank so se v obravnavanem primeru izkazali kot zelo robustni. Z dovolj gosto mrezo koncnih elementov smo lahko zajeli tudi vpliv višjih nihajnih oblik, ki so zaradi velike mase delca pomembno vplivale na odziv konstrukcije. Visokofrekvencni odziv je sicer kriv za neuspešne analize z Newmarkovo shemo, medtem ko se je analiza z modificirano Newmarkovo shemo izkazala za zelo stabilno, ce smo izbrali dovolj majhen casovni korak.
7.2.3 Trk masnega delca ob nosilec
Zanima nas dogajanje pri trku masnega delca ob nosilec. Trk nastopi v trenutku, ko postane razdalja med delcem in nosilcem, q, enaka nic. Izberemo elasticen model trka, predstavljen na sliki 7.7, ki privzame lokalno deformiranje teles v tocki kontakta: q < 0.
F
a
prosti let
trk
(a) Kontaktni model, q-razdalja med delcem in nosilcem, a-penetracija delca v nosilec
(b) Hertzov konstitutivni zakon trka [71]
Slika 7.7: Matematicni model trka masnega delca ob nosilec. Figure 7.7: Mathematical model of an impact of a mass particle against beam.
Velikost deformacije opišemo s spremenljivko a, ki meri penetracijo teles v kontaktu. Stanje kontakta
potem opišemo z algebrajsko vezno enacbo
O = q + a = 0.
Velikost penetracije je odvisna od materialnih in geometrijskih lastnosti teles, opisuje pa jo konstitutivni zakon trka. Konstitutivni zakon podaja zvezo med deformacijo teles v kontaktu, a, in kontaktno silo, Fc. V primeru trka masnega delca in nosilca, ki ga obravnavamo v nadaljevanju, izberemo konstitutivni zakon, ki gaje predlagal Hertz [71]
Fc = cxa3,	(7.14)
kjer je ci konstanta, odvisna od materialnih lastnosti teles v trku in geometrije kontaktne površine.
Obravnavamo primer nosilca, ki gaje predstavil Bauchau [34]. Nosilec ima dolzino L = 2 in kvadratni precni prerez s površino A = 10-4. Elasticni modul nosilca je E = 73 ■ 109, Poissonov kolicnik v = 0.3, gostota pa je enaka p = 2700. Nosilec na sredini razpona zadene krogla z maso m = 0.27 in zacetno hitrostjo vZ = -4. Opazujemo razlike v odzivu sistema glede na robne pogoje nosilca. Obravnavamo prostolezeci in nepomicno podprti nosilec (slika 7.8). Kroglo modeliramo kot masni delec, v analizi pa zanemarimo vpliv sile tezše in sile trenja med delcem in nosilcem in uposštevamo le kontaktno silo v smeri normale na nosilec. Vrednost konstante c1 v konstitutivnem zakonu trka (7.14) je c1 = 9.07 ■ 109.
(a) Prostolezšecš nosilec	(b) Nepomicšno podprt nosilec
Slika 7.8: Trk masnega delca ob nosilec. Figure 7.8: Impact of a mass particle against beam.
Za modeliranje nosilca izberemo gosto mrezo 80 elementov s konstantnim potekom krajevnih odvodov hitrosti. Analizo opravimo z modificirano Newmarkovo shemo s kratkim casovnim korakom h = 10-6. Rezultate analiz za oba nacina podpiranja nosilca prikazujemo na slikah 7.9 in 7.10.
Odziv masnega delca in nosilca je podoben v obeh primerih podpiranja. Bistveno razliko opazimo v dolzini casovnega intervala, ki ga merimo od zacetnega kontakta do trenutka, v katerem se delec popolnoma loci od nosilca in v nadaljevanju gibanja ne pride vec v stik z nosilcem. Dolzino tega intervala proglasimo za trajanje trka. Pri trku delca in prostolezecega nosilca je ta dolzina priblizno enaka t = 0.12, pri trku delca ob nepomicno podprt nosilec pa je ta dolzina priblizno štirikrat krajša t = 0.03. Na zacetku trka sta delec in nosilec v tesnem stiku, hitrost delca se zmanjšuje, v odzivu nosilca pa opazimo narašcajoce oscilacije. Oscilacije povzrocijo, da se delec in nosilec locita, v nadaljevanju pa veckrat trcita.
Delec se pri trku s prostolezecim nosilcem giblje navpicno navzdol priblizno do casa t = 0.075, kar je priblizno 62% trajanja celotnega trka, potem pa spremeni smer gibanja. To fazo trka imenujemo
0.08 0.06
0.02 0.04 0.06 0.08 0.1 0.12 0.14
t
(a) Vertikalni pomik masnega delca in tocke na sredini nosilca: element s konstantnimi krajevnimi odvodi hitrosti (-), Bauchau [34] ( )
3
-4
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
t
(b) Hitrost masnega delca: element s konstantnimi krajevnimi odvodi hitrosti (-), Bauchau [34] (-)
5000 4000 S 3000
JI
32000 Ö
J3 1000 0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
t
(c) Kontaktna sila med masnim delcem in nosilcem: element s konstantnimi krajevnimi odvodi hitrosti (-), Bauchau [34] (-)
Slika 7.9: Rezultati analize trka masnega delca ob prostoleZeci nosilec. Figure 7.9: Results of the analysis of the impact of a mass particle against simply supported beam.
1				
				
kompresijska faza, preostal del trka, ko se delec giblje navzgor, pa restitucijska faza. Tudi pri trku delca ob nepomicno podprti nosilec znaša trajanje kompresijske faze, ki se konca priblizno pri casu t = 0.0185, 62% dolzine celotnega trka. V obeh obravnavanih primerih ob prehodu kompresijske v restitucijsko fazo trka nastopijo največje kontaktne sile, ki so podobnih velikosti. Pri trku delca in prostolezecega nosilca je najvecja kontaktna sila priblizno enaka Fc = 5000, pri trku delca in nepomicno podprtega nosilca pa
Fc = 4500.
Ob koncu trka je zanimiva primerjava hitrosti delca, ki sta podobni v obeh primerih podpiranja nosilca. Ob koncu trka s prostolezecim nosilcem ima delec hitrost priblizno vz = 2.6, ob koncu trka z nepomicno podprtim nosilcem pa vZ = 2.2.
M
O
a
J3
M
tr >
0.015 0.01 0.005 0
-0.005 -0.01 0.015 -0.02 -0.025 -0.0
0 0.005 0.01
0.015 0.02
t
0.025 0.03 0.035
(a) Vertikalni pomik masnega delca in tocke na sredini nosilca: element s konstantnimi krajevnimi odvodi hitrosti (-), Bauchau [34] ( )
3
0.005 0.01 0.015 0.02 0.025 0.03 0.035
t
(b) Hitrost masnega delca
0
5000 ^ 4000 S 3000
JI
32000 Ö
J3 1000 0
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035
t
(c) Kontaktna sila med masnim delcem in nosilcem: element s konstantnimi krajevnimi odvodi hitrosti (-), Bauchau [34] ( )
Slika 7.10: Rezultati analize trka masnega delca ob nepomicno podprti nosilec. Figure 7.10: Results of the analysis of the impact of a mass particle against fully supported beam.
dilll I
Na slikah 7.9 in 7.10 so prikazane tudi rešitve, ki jih je predstavil Bauchau [34] za simulacijo polovicne
dolzine nosilca z upoštevanjem simetrije. Model nosilca je bil sestavljen iz 8 elementov s kubicno interpolacijo pomikov in rotacij. Simulacija je bila izvršena z energijsko konservativno shemo, podatek o uporabljenem casovnem koraku pa v clanku ni podan. Rezultati analiz za prostolezeci nosilec se bistveno razlikujejo v izracunanih kontaktnih silah, posledicno pa tudi v odzivu mase in nosilca v restitucijski fazi trka. Boljše je ujemanje rezultatov obeh simulacij za nepomicno podprt nosilec.
8 ZAKLJUČEK
V	doktorski nalogi smo obravnavali numericne metode za analizo konstrukcij, ki so podvrzene di-namicnim vplivom. Posebej smo obravnavali metode casovne integracije in reševanje dinamicnih enacb nosilcev. Obe temi v idejnem smislu povezuje razprava o izbiri osnovnih neznank problema, bodisi pri reševanju enacb v odvisnosti od casa bodisi pri reševanju enacb v odvisnosti od parametra kraja.
V	prakticnem smislu obe temi zdruzimo in pripravimo nabor koncnih elementov za dinamicno analizo prostorskih linijskih konstrukcij, ki za diskretizacijo enacb po casu uporabljajo predstavljene casovno integracijske sheme.
V	sklopu analize casovno integracijskih shem smo predstavili znano in splošno uveljavljeno Newmarkovo shemo, prilagojeno za nelinearen prostor rotacij. Poleg tega smo izpeljali povsem novo diskretizacijo kinematicnih enacb po casu, ki smo jo poimenovali modificirana Newmarkova shema. Glavne znacšilnosti nove implicitne sheme so:
•	Osnovne neznanke nove sheme so odvodi pomikov in rotacij po casu - hitrosti in kotne hitrosti.
•	Rezultat take izbire osnovnih neznank je povecana numericna stabilnost casovne integracije dis-kretiziranega sistema enacšb.
•	Poenostavi se linearizacija sheme in dodajanje popravkov, posledicno pa tudi racunska ucinkovi-tost algoritma.
•	Splošnost sheme omogoca njeno neposredno vkljucitev v ze obstojeca racunska okolja.
V	sklopu reševanja dinamicnih enacb prostorskih nosilcev smo predstavili novi formulaciji koncnih elementov. Glavne znacšilnosti formulacij so:
•	Formulaciji sta skladni z geometrijsko tocno teorijo nosilcev, saj smo pri izpeljavi enacb upoštevali tocne zveze med pomiki, zasuki in deformacijami.
•	Izpeljani formulaciji omogocata upoštevanje zahtevne (ukrivljene, zvite) zacetne geometrije.
•	V obeh formulacijah lahko uporabimo poljuben nelinearen materialni model, za katerega poznamo tangentno materialno matriko. V racunskih primerih smo obravnavali zgolj primere z linearno elasticšnim materialnim modelom.
•	Vodilna sistema enacb za obe formulaciji temeljita na obliki, ki stajo za staticno analizo prostorskih linijskih konstrukcij uveljavila Zupan in Saje [9]. Vodilni sistem enacb sestavljajo enacbe konsistence v krepki obliki ter kinematicni in ravnotezni robni pogoji. Z upoštevanjem konsi-stencnih pogojev dosezemo ujemanje ravnoteznih in konstitucijskih notranjih sil in momentov, kar je redkost v formulacijah koncšnih elementov.
•	Osnovne interpolirane neznanke v obeh formulacijah so deformacijske kolicine. Zaradi aditivnosti deformacijskih kolicin enake interpolacijske nastavke uporabimo tudi za njihove variacije. V prvi
formulaciji interpoliramo krajevne odvode pomikov in zasukov - deformacije - v drugi formulaciji pa krajevne odvode hitrosti in kotne hitrosti, za katere pokazemo, da so hkrati tudi mere za hitrost deformacij.
•	Zvezne enacbe konsistence diskretiziramo s kolokacijsko metodo, tako da ujemanje ravnoteznih in konstitucijskih sil in momentov zagotovimo v vnaprej dolocenih diskretnih tockah na tezišcni osi nosilca. Pri numericnem reševanju poenotimo diskretizacijske in interpolacijske tocke, s cimer pridobimo na racunski ucinkovitosti in natancnosti.
•	Pri elementu z interpolacijo deformacij kinematicne enacbe po casu diskteriziramo z Newmarkovo shemo. S tako izbiro v formulacijo koncšnega elementa vpeljemo lastnosti sheme, ki je pri nelinearnih problemih zgolj pogojno stabilna. Stabilnost izboljšamo tako, da v formulacijo koncnega elementa vpeljemo numericno disipacijo s posplošeno-a metodo.
•	V formulaciji koncnega elementa z interpolacijo krajevnih odvodov hitrosti za diskretizacijo kine-maticnih zvez po casu uporabimo modificirano Newmarkovo shemo. Izboljšana stabilnost sheme omogoca racun zahtevnejših primerov brez uporabe numericnega dušenja. Element zelo dobro deluje pri analizi primerov, v katerih rešitve išcemo na daljšem casovnem intervalu, in v primerih, kjer so za zadostno natancnost potrebni kratki casovni koraki.
•	Z numericšnimi primeri pokazšemo, da obe formulaciji koncšnih elementov zelo dobro opisšeta obnašanje konstrukcij pri neomejeno velikih pomikih, zasukih in deformacijah.
V zadnjem delu doktorske naloge smo na kratko obravnavali povezan problem gibanja masnega delca in nosilca. Gibalne enacbe delca smo dodali sistemu enacb konstrukcije, kontaktne pogoje pa smo vpeljali prek kinematicšnih vezi z metodo Lagrangejevih mnozšiteljev. Pri aplikaciji numericšnega postopka smo prisšli do naslednjih ugotovitev:
•	Pri analizi primerov smo opazili slabši odziv elementov, pri katerih za interpolacijo kolicin uporabljamo Lagrangejeve polinome. Zato smo za povezano dinamicšno analizo prostorskih nosilcev ob upoštevanju razlicnih vplivov zunanjih delcev z maso, pripravili nabor koncnih elementov s predpostavljenim konstantnim potekom deformacijskih kolicin. Za posebej robustne so se izkazali elementi s predpostavljenim konstantnim potekom krajevnih odvodov hitrosti in kotnih hitrosti.
•	Podrobneje smo obravnavali trk masnega delca ob prostolezec in nepomicno podprt nosilec. Nacin podpiranja nosilca je bistveno vplival le na trajanje trka. Velikost kontaktnih sil in hitrost delca med trkom pa sta bila predvsem odvisna od uporabljenega kontaktnega modela, mase in zacetne hitrosti delca.
POVZETEK
V disertaciji se ukvarjamo z dinamicno analizo prostorskih linijskih konstrukcij. Glavni namen dela je priprava novih matematicnih orodij za natancno in numericno stabilno reševanje enacb, s katerimi opišemo odziv konstrukcij, ki s podvrzene dinamicnim vplivom.
Matematicni model linijskega nosilca postavimo skladno z geometrijsko tocno Reissner-Simovo [1,2] teorijo nosilcev. Vodilo izpeljav je pristop, ki sta ga za staticšno analizo prostorskih nosilcev uveljavila Zupan in Saje [9]. Glavne znacilnosti pristopa so: (i) omogoca ujemanje konstitucijskih in ravnoteznih sil in momentov v krepki obliki, (ii) zadošca tocnim kinematicnim zvezam med pomiki, rotacijami in deformacijami, (iii) osnovne neznanke so deformacijske kolicine in (iv) zvezne enacbe so diskretizirane s kolokacijsko metodo. Pri formulaciji koncnih elementov za dinamiko zelimo ohraniti enako strukturo enacb in izbiro osnovnih neznank, s cimer se pri izbiri casovno integracijskih shem omejimo na algoritme, ki diskretizirajo kinematicne zveze med pomiki, hitrostmi in pospeški ter rotacijami, kotnimi hitrostmi in kotnimi pospeški.
Problemi v dinamiki konstrukcij in mehanizmov so obicajno togi, zato za casovno integracijo dinamicnih enacb izberemo implicitne sheme. Najprej predstavimo Newmarkovo shemo [28] in njeno razširitev na nelinearen prostor rotacij [4], v kateri kot osnovne neznanke nastopajo pomiki in zasuki. Pri reševanju problemov v nelinearni dinamiki je Newmarkova shema zgolj pogojno stabilna, zato predlagamo novo -modificirano Newmarkovo shemo. Osnovne neznanke nove sheme so odvodi pomikov in rotacij po casu - hitrosti in kotne hitrosti. Tak izbor osnovnih neznank vpliva na povecano stabilnost casovne integracije pri analizah, v katerih rešitve išcemo na daljšem casovnem intervalu, in še posebej pri analizah, ki jih racšunamo s kratkim cšasovnim korakom. Osnovne neznanke nove sheme so aditivne kolicšine, kar mocšno poenostavi iteracijsko reševanje enacb in racunsko ucinkovitost algoritma.
Enacbe dinamicnega ravnotezja nosilca izpeljemo iz izrekov o gibalni in vrtilni kolicini. Za opis materialnih lastnosti nosilca privzamemo povsem splosšen materialni model, predpostavimo le, da poznamo tangentno materialno matriko poljubnega prereza nosilca. Pri numericnem reševanju enacb nosilca kon-stitucijska sila in moment ne zadošcata nujno ravnoteznim enacbam. Sledimo konceptu, ki sta ga uveljavila Vratanar in Saje [52], in z enacbami konsistence eksplicitno zahtevamo ujemanje ravnoteznih in konstitucijskih sil in momentov. Predstavimo dve formulaciji linijskih koncnih elementov, ki temeljita na izbiri deformacijskih kolicin za osnovne neznanke problema. Vodilne enacbe obeh formulacij so sestavljene iz konsistencnih enacb, kinematicnih in ravnoteznih robnih pogojev. V prvi formulaciji za osnovne neznanke izberemo krajevne odvode pomikov in rotacij - deformacije. Kinematicne zveze po casu dis-kretiziramo z Newmarkovo shemo, za zagotovitev stabilnosti pa vpeljemo posplošeno-a metodo [24] numericnega dušenja. V drugi formulaciji za diskretizacijo kinematicnih enacb po casu izberemo modificirano Newmarkovo shemo. Iz kinematicnih zvez med hitrostmi in deformacijami sledi, da so krajevni odvodi hitrosti in kotnih hitrosti tudi mere za hitrost deformacij. To nas motivira, da jih v formulaciji z modificirano Newmarkovo shemo izberemo za osnovne neznanke diskretizacije kinematicšnih kolicšin po kraju. Osnovne neznanke v obeh formulacijah so aditivne kolicine, zato jih izberemo za edine kolicine, ki jih interpoliramo. Zvezne enacbe konsistence po kraju dikretiziramo s kolokacijsko metodo, tako da
enakost ravnoteznih in konstitucijskih sil in momentov zahtevamo v naprej izbranih diskretnih tockah.
V	numericšni implementaciji koncšnih elementov tocške interpolacije in kolokacije poenotimo, s cšimer pridobimo na natancnosti in racunski ucinkovitosti algoritma.
Vodilni sistem enacb obeh formulacij koncnih elementov s krajevno in casovno diskretizacijo prevedemo na sistem nelinearnih algebrajskih enacb, ki ga rešujemo z Newtonovo iteracijsko metodo. Za obe formulaciji podrobno predstavimo linearizacijo in dodajanje popravkov. Nove elemente vgradimo v racunalniški program in testiramo na številnih primerih, rezultate pa primerjamo z rešitvami drugih avtorjev. Z numericnimi študijami pokazemo, da obe formulaciji dobro opišeta obnašanje konstrukcije pri poljubno velikih pomikih, zasukih in deformacijah. Izboljšana stabilnost modificirane Newmarkove sheme pride do izraza v formulaciji elementov z interpolacijo krajevnih odvodov hitrosti, saj omogoca racun zahtevnih racunskih primerov brez numericnega dušenja.
V	zadnjem delu doktorske naloge obravnavamo povezan problem gibanja masnega delca in nosilca. Gibalne enacbe delca dodamo sistemu enacb konstrukcije, kontaktne pogoje pa vpeljemo prek kinematicnih vezi [34] z metodo Lagrangejevih mnoziteljev. Za povezano dinamicno analizo prostorskih nosilcev ob upoštevanju razlicnih vplivov zunanjih delcev z maso pripravimo nabor koncnih elementov s predpostavljenim konstantnim potekom deformacijskih kolicin. Za posebej robustne se izkazejo elementi s predpostavljenim konstantnim potekom krajevnih odvodov hitrosti in kotnih hitrosti. Podrobneje analiziramo trk masnega delca ob nosilec. Trk obravnavamo kot zvezen dogodek s koncnim trajanjem, pri katerem se telesi v obmocšju kontakta deformirata.
SUMMARY
In the present dissertation we study dynamics of beam-like structures in three-dimensional space. The main goal of the present work is the development of new mathematical tools for accurate and numerically stable solution of equations which describe the response of the structures under dynamic influences.
We formulate a mathematical model of the beam based on the Reissner-Simo [1,2], geometrically exact beam theory. We follow the approach presented by Zupan and Saje [9] for static analysis of spatial beams. The main properties of their approach are: (i) the equality of the equilibrium and constitutive internal forces is satisfied in strong form, (ii) kinematic relations between displacements, rotations and strains are exactly satisfied, (iii) strain measures are the primary variables and (iv) continuous equations are discretized by the use of collocation algorithm. In the finite-element formulation for dynamic analysis we want to keep the same structure of equations and the selection of the primary unknowns. By doing so, we limit the choice of the possible time integration schemes to the algorithms which discretize kinematic relations between displacements, velocities and accelerations and the relations between rotations, angular velocities and angular accelerations.
Problems in dynamics of structures and mechanisms are often described by stiff equations, therefore we select implicit schemes for time integration of dynamic equations. First we introduce Newmark scheme [28] and its extension to the non-linear space of rotations [4], where displacements and rotations are chosen for the primary unknowns. The Newmark scheme is only conditionally stable in non-linear problems, therefore we develop a new - modified Newmark scheme, with improved stability properties. The primary unknowns of the new scheme are time derivatives of displacements and rotations - velocities and angular velocities. Such a selection of primary unknowns improves the stability of the long term analysis and the analysis where small time steps are required. The primary unknowns of the modified scheme are additive quantities which simplifies the iteration process and increases numerical efficiency of the algorithm.
The equations of the dynamic equilibrium of the beam are derived from linear and angular momentum conservation laws. The material properties of the beam are represented by a general form of constitutive equations, where we assume that the tangent constitutive matrix of an arbitrary cross-section of the beam is known. In the numerical solution of the equations of the beam the constitutive force and moment vectors do not automatically satisfy the equilibrium equations. We follow the concept introduced by Vratanar and Saje [52] and explicitly demand the equality of the constitutive and equilibrium force and moment vectors. We introduce two finite-element formulations based on strain measures. The governing system of equations of each formulation consist of consistency conditions, kinematic and equilibrium boundary conditions. The first formulation employs spatial derivatives of displacements and rotations - strains, as the primary unknowns. The Newmark scheme is used for time discretization of kinematic quantities and Generalized-a method is used to assure numerical stability of the strain-based formulation. The second formulation employs the modified Newmark scheme for time discretization of kinematic equations. It follows from kinematic relations between strains and velocities that spatial derivatives of velocities and angular velocities are also strain-rate measures. Thus, we select them as the
primary unknown functions of the spatial discretization of the formulation with the modified Newmark scheme. The primary unknowns of both formulations are additive quantities and are selected to be the only interpolated quantities of the formulations. The continuous consistency equations are discretized by the use of a collocation method, such that the equality of the equilibrium and constitutive forces is satisfied at the chosen discrete points. In numerical implementation of the finite-elements we make the interpolation points to coincide with collocation points, which improves the accuracy and computational efficiency of the algorithm.
When the governing equations of the two formulations are discretized with respect to space and time we obtain non-linear system of equations, which is solved by the iterative Newton's method. Therefore the linearization of equations and the update procedure are presented. The computer code is generated and the performance of the finite-elements is tested on several numerical examples. The results are compared to results presented by other authors. We show that both formulations give accurate results and are able to describe the behaviour of structures where finite displacements, rotations and strains occur. The improved stability of the modified Newmark scheme, implemented in velocity-based formulation, enables us to obtain results of more demanding problems without numerical dissipation.
In the last part of the thesis we consider coupled motion of a mass particle and a beam. The system of equations of the beam is supplemented with equations of motion of a particle; contact conditions are enforced as constraints [34], by the use of Lagrange multiplier technique. For coupled dynamic analysis of spatial beams and various influences of mass particles, we present a special set of finite-elements with assumed constant strain measures. The formulation with assumed constant spatial derivatives of velocities and angular velocities shows great robustness. The impact of a mass particle against the beam is studied in detail. The impact is modelled as an event of finite duration, where local deformation of the contacting bodies occur.
VIRI
[1]	Reissner, E. 1981. On finite deformation of space-curved beams. Z. Angew. Math. Phys. 32: 734-744.
[2]	Simo, J. Č. 1985. A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Čomput. Meth. Appl. Mech. Eng. 49: 55-70.
[3]	Simo, J. Č., Vu-Quoc, L. 1986. A three-dimensional finite-strain rod model. Part II: Čomputational aspects. Čomput. Meth. Appl. Mech. Eng. 58: 79-116.
[4]	Simo, J. Č., Vu-Quoc, L. 1988. On the dynamics in space of rods undergoing large motions - A geometrically exact approach. Čomput. Meth. Appl. Mech. Eng. 66: 125-161.
[5]	Čardona, A., Geradin, M. 1988. A beam finite element non-linear theory with finite rotations. Int. J. Numer. Meth. Eng. 26: 2403-2438.
[6]	Ibrahimbegovic, A. 1995. On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements. Čomput. Meth. Appl. Mech. Eng. 122: 11-26.
[7]	Jelenic, G., Saje, M. 1995. A kinematically exact space finite strain beam model - finite element formulation by generalized virtual work principle. Čomput. Meth. Appl. Mech. Eng. 120: 131-161.
[8]	Jelenic, G., Črisfield, M. A. 1999. Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics. Čomput. Meth. Appl. Mech. Eng. 171: 141-171.
[9]	Zupan, D., Saje, M. 2003. Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures. Čomput. Meth. Appl. Mech. Eng. 192: 5209-5248.
[10]	Simo, J. Č., Tarnow, N., Doblare, M. 1995. Non-linear dynamics of three-dimensional rods: exact energy and momentum conserving algorithms. Int. J. Numer. Meth. Engng. 38: 1431-1473.
[11]	Bauchau, O. A., Theron, N. J. 1996. Energy decaying scheme for nonlinear elastic multi-body systems. Čomput. Struct. 59: 317-331.
[12]	Bottasso, Č. L., Borri, M. 1997. Energy preserving/decaying schemes for non-linear beam dynamics using the helicoidal approximation. Čomput. Meth. Appl. Mech. Eng. 143: 393-415.
[13]	Ibrahimbegovic, A., al Mikdad, M. 1998. Finite rotations in dynamics of beams and implicit time-stepping schemes. Int. J. Numer. Meth. Engng. 41: 781-814.
[14]	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. Meth. Engng. 54: 1683-1716.
[15]	Argyris, J. H. 1982. An excursion into large rotations. Čomput. Meth. Appl. Mech. Eng. 32: 85-155.
[16]	Atluri, S., Čazzani, A. 1995. Rotations in computational solid mechanics. Arch. Čomput. Meth. Eng. 2: 49-138.
[17]	Zupan, D., Saje, M. 2003. The three-dimensional beam theory: Finite element formulation based oncurvature. Čomput. Struct. 81: 1875-1888.
[18]	Jelenic, G., Črisfield, M. A. 1999. Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation. Proc. Roy. Soc. Lond. Math. Phys. Sci. 455: 1125-1147.
[19]	Planinc, I. 1998. Racun kriticnih tock konstrukcij s kvadraticno konvergentnimi metodami. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo (samo-zalozba I. Planinc): 83 f.
[20]	Čhoi, J., Lim, J. 1995. General curved beam elements based on the assumed strain fields. Čomput. Struct. 55: 379-386.
[21]	Friedman, Z., Kosmatka, J. B. 1998. An accurate two-node finite element for shear deformable curved beams. Int. J. Numer. Meth. Engng. 41: 473-498.
[22]	Schulz, M., Filippou, F. Č. 2001. Non-linear spatial Timoshenko beam element with curvature interpolation. Int. J. Numer. Meth. Engng. 50: 7361-785.
[23]	Hilber, H. M., Hughes, T. J. R., Taylor, R. L. 1977. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake. Eng. Struct. Dynam. 5: 283-292.
[24]	Čhung, J., Hulbert, G. M. 1993. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-a method. J. Appl. Mech. 60: 371-375.
[25]	Hughes, T. J. R. 1983. Analysis of transient algorithms with particular reference to stability beha-vior. V: Belytscho, T. (ur.), Hughes, T. J. R. (ur.). Čomputational Methods for Transient Analysis. Amsterdam, North Holland: str. 67-155.
[26]	Simo, J. Č., Wong, K. K. 1991. Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. Int. J. Numer. Meth. Eng. 31: 19-52.
[27]	Ibrahimbegovic, A., Mamouri, S. 2002. Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations. Čomput. Meth. Appl. Mech. Eng. 191: 4241-4258.
[28]	Newmark, N. M. 1959. A method of computation for structural dynamics. J. Eng. Mech. Div. 85: 67-94.
[29]	Bottasso, Č. L., Borri, M. 1998. Integrating finite rotations. Čomput. Meth. Appl. Mech. Eng. 164: 307-331.
[30]	Lee, H. P. 1996. The dynamic response of a Timoshenko beam subjected to a moving mass. J. Sound. Vib. 198: 249-256.
[31]	Zupan, E. 2010. Dinamika prostorskih nosilcev s kvartenionsko parametrizacijo rotacij. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo (samozalozba E. Zupan): 163 str.
[32]	Bauchau, O. A. 1998. Čomputational Schemes for Flexible, Nonlinear Multi-Body systems. Mul-tibody System Dynamics 2: 169-225.
[33]	Bauchau, O. A. 1999. On the Modeling of Friction and Rolling in Flexible Multi-Body systems. Multibody System Dynamics 3: 209-239.
[34]	Bauchau, O. A. 2000. Analysis of Flexible Multibody systems with Intermittent Čontacts. Multibody System Dynamics 4: 23-54.
[35]	Lens, E. V., Čardona, A. 2008. A nonlinear beam element formulation in the framework of an energy preserving time integration scheme for constrained multibody system dynamics. Čomput. Struct. 86: 67-94.
[36]	Zupan, D. 2003. Rotacijsko invariantne deformacijske kolicine v geometrijsko tocni teoriji prostorskih nosilcev. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo (samozalozba D. Zupan): 210 str.
[37]	Črisfield, M. A. 1997. Non-linear finite Element Analysis of Solids and Structures, Volume 2, Advanced Topics. Čhichester - New York - Weinheim - Brisbane - Singapore - Toronto, John Wiley & Sons: 494 str.
[38]	Ibrahimbegovic, A. 1997. On the choice of finite rotation parameters. Čomput. Meth. Appl. Mech. Eng. 149: 49-71.
[39]	Čurtiss, Č. F., Hirscfelder, J. 1952. Integration of stiff equations. Proc. Natl Acad. Sci. USA 38: 235-243.
[40]	Spurrier, R. A. 1978. Čomment on 'singularity-free extraction of a quaternion from a direction-cosine matrix'. J. Spacecraft Rockets 15: 255.
[41]	Bathe, K.-J. 1982. Finite Element Procedures in Engineering Analysis. Englewood Čliffs, New Jersey, Prentice-Hall: 735 str.
[42]	Hughes, T. J. R. 1976. Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics. Čomput. Struct. 6: 313-324.
[43]	Hosea, M. E., Shampine, L. F. 1996. Analysis and implementation of TR-BDF2. Appl. Numer. Math. 20: 21-37.
[44]	Vidav, I. 1990. Višja matematika I. Ljubljana, DMFA: 477 str.
[45]	Krysl, P. 2005. Explicit momentum-conserving integrator for dynamics of rigid bodies approxima-ting the midpoint Lie algorithm. Int. J. Numer. Meth. Eng. 63: 2171-2193.
[46]	Hairer, E., Wanner, G. 1991. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Berlin, Springer: 594 str.
[47]	Arnold, V. I. 1989. Mathematical Methods of Člassical Mechanics, New York, Springer: 516 str.
[48]	Shampine, L. F., Reichelt, M. W. 1997. The MATLAB ODE Suite. SIAM J. Sci. Čomput. 18: 1-22.
[49]	Lens, E., Čardona, A. 2007. An energy preserving/decaying scheme for nonlinearly constrained multibody systems. Multibody System Dynamics 18: 435-470.
[50]	Petzold, L., Lotstedt, P. 1986. Numerical solution of nonlinear differential equations with algebraic constraints II: Practical implications. SIAM J. Sci. Statist. Čomput. 7: 720-733.
[51]	Bottasso, Č. L., Bauchau, O. A., Čardona, A. 2007. Time-step-size-independent conditioning and sensitivity to perturbations in the numerical solution of index three differential algebraic equations. SIAM Journal on Scientific Čomputing 29: 397-414.
[52]	Vratanar, B., Saje, M. 1999. A consistent equilibrium in a cross-section of an elastic-plastic beam. Int. J. Solid. Struct. 36: 311-337.
[53]	Bottasso, Č., Borri, M. 1997. Energy preserving/decaying schemes for non-linear beam dynamics using the helicoidal approximation. Čomput. Meth. Appl. Mech. Eng. 143: 393-415.
[54]	Betsch, P., Steinmann, P. 2002. A DAE Approach to Flexible Multibody Dynamics. Multibody System Dynamics 8: 367-391.
[55]	Wood, W. L., Bossak, M., Zienkiewicz, O. Č. 1981. An Alpha Modification of Newmark's Method. Int. J. Numer. Meth. Eng. 15: 1562-1566.
[56]	Erlicher, S., Bonaventura, L., Bursi, O. S. 2002. The analysis of the Generalized-a method for non-linear dynamic problems. Čomput. Mech. 28: 83-104.
[57]	Kuhl, D., Črisfield, M. A. 1999. Energy-conserving and decaying algorithms in non-linear structural dynamics. Int. J. Numer. Meth. Eng. 45: 569-599.
[58]	W. Weaver, J., Timoshenko, S. P., Young, D. H. 1990. Vibration Problems in Engineering. 5. izdaja. Čhichester - New York - Weinheim - Brisbane - Singapore - Toronto, John Wiley & Sons: 610 str.
[59]	Bauchau, O. A. Dymore Solutions. Simulation Tools for Flexible Multibody Systems. Benchmark test cases: four-bar mechanism. http://www.dymoresolutions.com/Benchmarks/FourBar/ HeadFourBar.html (Pridobljeno 1.10.2012).
[60]	Bauchau, O. A. Dymore Solutions. Simulation Tools for Flexible Multibody Systems. http://www. dymoresolutions.com (Pridobljeno 1.10.2012).
[61]	Črisfield, M. A., Shi, J. 1994. A co-rotational element/time-integration strategy for non-linear dynamics. Int. J. Numer. Meth. Eng. 37: 1897-1913.
[62]	Črisfield, M. A., Galvanetto, U., Jelenic, G. 1997. Dynamics of 3-D co-rotational beams. Čomput. Mech. 20: 507-519.
[63]	Wehage, R. A., Haug, E. J. 1981. Dynamic analysis of mechanical systems with intermittent motion. ASME J. Mech. Design 103: 560-570.
[64]	Khulief, Y. A., Shabana, A. A. 1986. Dynamic analysis of constrained systems of rigid and flexible bodies with intermittent motion. ASME J. Mech. Transmissions Automat. Design 108: 38-44.
[65]	Bratina, S. 2003. Odziv armiranobetonskih linijskih konstrukcij na pozarno obtezbo. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo (samozalozba S. Bratina): 159 str.
[66]	(Češarek, P., Saje, M., Zupan, D. 2012. Kinematically exact curved and twisted strain-based beam. Int. J. Solid. Struct. 49: 1802-1817.
[67]	Bauchau, O. A., Damilano, G., Theron, N. J. 1995. Numerical integration of non-linear elastic multi-body systems. Int. J. Numer. Meth. Engng. 38: 2727-2751.
[68]	Ibrahimbegovic, A., Mamouri, S. 1999. Nonlinear dynamics of flexible beams in planar motion: formulation and time-stepping scheme for stiff problems. Comput. Struct. 70: 1-22.
[69]	Gams, M. 2008. Geometrijsko tocna dinamicna analiza elasticnih in armiranobetonskih okvirnih konstrukcij. Doktorska disertacija. Ljubljana, Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo (samozalozba M. Gams): 107 str.
[70]	Bauchau, O. A., Theron, N. J. 1996. Energy decaying scheme for non-linear beam models. Čomput. Meth. Appl. Mech. Eng. 134: 37-56.
[71]	Timoshenko, S. P., Gere, J. M. 1961. Theory of Elastic Stability. New York - Toronto - London, McGraw-Hill: 541 str.