Indice 1
Risoluzione della cinematica
5
1.1 1.1 Quad Quadri rila late tero ro prin princi cippale ale - Sos Sostegn tegnoo ruot ruotaa . . . . . . . . 1.1.1 Equazione in posizione . . . . . . . . . . . . . . 1.1.2 Equazione in velocità . . . . . . . . . . . . . . . 1.1.3 Equazione in accelerazione . . . . . . . . . . . . 1.2 1.2 Seco Second ndoo quad quadri rila late tero ro - Coll Colleg egam amen ento to barr barraa di tors torsio ione ne 1.2.1 Equazione in posizione . . . . . . . . . . . . . . 1.2.2 Equazione in velocità . . . . . . . . . . . . . . . 1.2.3 Equazione in accelerazione . . . . . . . . . . . . 2
. . . . . . . .
. . . . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . 5 . . 5 . . 7 . . 7 . . 8 . . 8 . . 9 . . 10 11
Energia cinetica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Energia pot potenziale gravitazionale ale . . . . . . . . . . . . . . . . . . . . Energi Energiaa potenz potenzial ialee e dissip dissipati ativva del del siste sistema ma moll molla-s a-smor morzat zatore ore a terra terra Ener Energi giaa pote potenz nzia iale le elas elasti tica ca dell dellaa mo moll llaa tors torsio iona nale le . . . . . . . . . . . Energia dissipativa dello smorzatore . . . . . . . . . . . . . . . . . . Equazione di moto e forma di stato . . . . . . . . . . . . . . . . . . . Precar Precarico ico static staticoo ∆γ 20 20 . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . .. . . . . . . . .. . .. . . ..
Equazione di moto linearizzata
3.1 3.2 3.3 3.4 4
. . . . . . . .
Scrittura dell’equazione di moto
2.1 2.2 2.3 2.4 2.4 2.5 2.6 2.7 3
. . . . . . . .
Energia cinetica . . . Energia potenziale . Funzione dissipativa Equazione di moto .
. . . .
. . . .
. . . .
. . . .
. . . .
11 14 16 18 18 20 21 22
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
Analisi del comportamento del sistema
4.1 Frequenza propria . . . . . . . . . . 4.2 Moto libero . . . . . . . . . . . . . . 4.2.1 Piccole per perturbazioni . . . . . 4.2.2 Perturbazioni di media entità 4.2.3 Grandi per perturbazioni . . . . . 4.3 Moto forzato - zona quasi-statica . . 4.4 Moto forzato - zona di risonanza . . 1
. . . . . . .
22 22 24 24 25
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .. . . .. ..
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
25 26 26 26 27 28 28
2
INDICE
4.5 Moto forzato - zona sismografica . . . . . . . . . . . . . . . . . . . . . . . 29 5
Scripts del programma Matlab
5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9
main.m . . . . . . . . . . . cinematica.m . . . . . . . . precarico.m . . . . . . . . . energiacinetica.m . . . . . . energiapotenziale.m . . . . . smorzatore.m . . . . . . . . manovellismoNONlineare.m manovellismolineare.m . . . rk4.m . . . . . . . . . . . .
. . . . . . . . .
31
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
31 34 35 35 36 37 37 38 38
Elenco delle figure 1 2
Vettura di F1 Lotus72 del 1971 . . . . . . . . . . . . . . . . . . . . . . . . Disegno della sospen pensione in analisi . . . . . . . . . . . . . . . . . . . . . .
4 4
1.1 1.1 Chiu Chiusu sura ra del del quad quadri rila late tero ro arti artico cola lato to prin princi cipa pale le . . . . . . . . . . . . . . . . 1.2 1.2 Chiu Chiussura ura del del secon econdo do quadr uadril ilat ateero arti artico cola lato to . . . . . . . . . . . . . . . . .
6 9
2.1 2.1 2.2 2.2 2.3 2.3 2.4
Posi Posizi zion onee gene generi rica ca dei dei bari barice cenntri tri per i due due quad quadri rila late teri ri . Quot Qu otee dei dei bari barice cenntri tri rispe rispett ttoo ad ad un asse asse oriz orizzo zont ntal alee fisso fisso Rapp Rappre rese senntazi tazion onee semp sempli lific ficat ataa dell dellaa sosp sospen ensi sion onee . . . . . Metodi Metodi di calcol calcoloo di di ∆l2 . . . . . . . . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
12 15 16 19
4.1 4.2 4.3 4.4 4.4 4.5 4.6
deg . . . . . . Risposta Risposta alla alla condizione condizione iniziale iniziale (α − αs ) = +3 deg Risposta Risposta alla alla condizione condizione iniziale iniziale (α − αs ) = −26deg . . . . . Risposta Risposta alla alla condizione condizione iniziale iniziale (α − αs ) = −65deg . . . . . Ris Rispost postaa in condi ondizzioni ioni quasi uasi-s -sta tati ticche . . . . . . . . . . . . . . Rispost Rispostaa in condiz condizion ionii di risona risonanza nza (Ω = ω0 = 51.1465 rad/s) Rispos posta in condizioni sismografiche . . . . . . . . . . . . . .
. . . . . .
. . . .. . . .. . .. . .. . .. .. .. .. . . . . . . . .. .. .
26 27 27 28 29 30
3
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
4
ELENCO DELLE FIGURE
Figura 1: Vettura di F1 Lotus72 del 1971
(a) Terminologia relativa a generica sospensione a quadrilatero
(b) Sospensione in analisi (rappresentazione speculare)
Figura 2: Disegno della sospensione in analisi
Capitolo 1
Risoluzione della cinematica La sospensione in esame risulta caratterizzata, dal punto di vista cinematico, da un doppio quadrilatero articolato (Figura 2 nella pagina precedente). Più in particolare, un primo quadrilatero (considerato principale) permette il sostegno ed il collegamento della ruota alla cassa della vettura; su di esso quindi agirà direttamente il “forzamento” associato alle vibrazioni provenienti dal terreno. L’elemento smorzatore della sospensione è inoltre incernierato alla biella di tale cinematismo. Il secondo quadrilatero ha il compito invece di azionare la barra torsionale, l’elemento elastico della sospensione. Si nota che il braccetto superiore costituisce la manovella superiore di entrambi i quadrilat drilateri eri.. Si decide decide pertanto pertanto di assume assumere re la rotazi rotazione one di tale tale elemen elemento to (α) come unica coordin coordinata ata libera di tutto tutto il sistem sistemaa in analisi: analisi: le equazion equazionii di moto moto (non (non linear linearee e linearizzata) saranno ricavate in funzione di essa.
1.1
Quadr Quadrila ilater tero o princi principal pale e - Sosteg Sostegno no ruota ruota
Per la risoluzione della cinematica si utilizzerà il metodo delle equazioni di chiusura, scegliendo come asse reale quello passante per le due cerniere a terra del quadrilatero. 1.1. 1.1.1 1
Equa Eq uazio zione ne in in posi posizio zione ne
Chiusura del quadrilatero articolato: aeiα + beiβ = ceiγ + d
Proiettando sui due assi reale e immaginario si ha:
b cos β = −a cos α + c cos γ + d b sin β = −a sin α + c sin γ
L’equazio L’equazione ne è non lineare in α, β e γ e quindi si ricorrerà alla soluzione numerica implementata nei capitoli seguenti. Si propone però anche un metodo algebrico che consente 5
CAPITOLO 1. RISOLUZIONE RISOLUZIONE DELLA CINEMATICA CINEMATICA
6
Figura 1.1: Chiusura del quadrilatero articolato principale ugualmente di ricavare le grandezze di interesse, ma che non verrà utilizzato. Approccio analitico
Quadrando e sommando è possibile eliminare β : a2 − b2 + c2 + d2 − 2ac cos α cos γ − 2ac sin γ sin α − 2ad cos α + 2cd cos γ = 0
si pone: A = −2ac sin α B = 2cd − 2ac cos α C = a2 − b2 + c2 + d2 − 2ad cos α D=
A2 + B 2 − C 2
si nota che i quattro parametri ottenuti hanno dipendenza soltanto dal gdl α, quindi si cercherà cercherà di esprimere esprimere le grandezze grandezze cinematich cinematichee ricercate ricercate in funzione funzione di essi. essi. Attrav Attraverso erso passaggi algebrici si ottiene: AC + BD A2 + B 2 BC + BD cos γ = − 2 A + B2
sin γ = −
7
CAPITOLO 1. RISOLUZIONE RISOLUZIONE DELLA CINEMATICA CINEMATICA
quindi noto γ : c sin γ − a sin α b d + c cos γ − a cos α cos β = b
sin β =
1.1. 1.1.2 2
Equa Eq uazio zione ne in in veloc velocit ità à
Derivando rispetto al tempo i termini dell’equazione in posizione si ottiene: ˙ iβ = iγce ˙ iα + iβbe ˙ iγ iαae iβ
moltiplicando l’equazione trovata una volta per e
−
˙ i(α iαae ˙ i(α iαae
e una volta per e
iγ
−
si ottiene:
˙ − iγce + iβb ˙ i(γ β ) = 0 γ ) ˙ i(β γ ) − iγc + iβbe ˙ =0 β)
−
−
−
−
Proiettando sull’asse reale:
˙ sin(α − β ) − γc ˙ sin(γ − β ) = 0 αa ˙ sin(β − γ ) = 0 ˙ sin(α − γ ) + βb αa
da cui a sin(α − γ ) b sin(γ − β ) a sin(α − β ) ˙ = α˙ γ c sin(γ − β )
˙ = α˙ β
∂β α˙ Poichè β = β (α) e γ = γ (α), applicando la regola di derivazione a cascata si ha β ˙ = ∂α ∂γ ∂β ∂γ e γ ˙ = ∂α α˙ , quindi si possono esplicitare i coefficienti f b = ∂α e f c = ∂α che risultano costanti ad ogni data posizone α:
1.1.3 1.1.3
f b =
a sin(α − γ ) b sin(γ − β )
(1.1)
f c =
a sin(α − β ) c sin(γ − β )
(1.2)
Equazi Equ azione one in acceler accelerazio azione ne
Derivando una seconda volta l’equazione di chiusura rispetto al tempo si ha: ¨ iβ ¨ aeiα − α˙ 2 aeiα + iβbe iαae α
−
˙ 2 beiβ = iγce ¨ce iγ − γ ˙ 2 ceiγ β γ
come in precedenza si moltiplica una volta per e
¨ aei(α iαae α ¨ aei(α iαae α
β)
α˙ 2 aei(α γ ) ˙ 2 aei(α −α
−
−
−
iβ
−
e una volta per e
¨ − β ˙ 2 b − iγce + iβb ¨ce i(γ γ γ ) ¨ i(β γ ) − β ˙ 2 bei(β + iβbe β)
−
−
β)
−
−
iγ
−
:
+ γ ˙ 2 cei(γ β ) = 0 γ ) ¨c + γ ˙ 2c = 0 − iγc γ
−
−
CAPITOLO 1. RISOLUZIONE RISOLUZIONE DELLA CINEMATICA CINEMATICA
8
Proiettando sull’asse reale si determinano le accelerazioni angolari: ¨ a sin(α − γ ) + α˙ 2 a cos(α − γ ) + β ˙ 2 b cos(β − γ ) − γ ˙ 2c αa α ¨ β = b sin(γ − β ) ¨ a sin(α − β ) + α˙ 2 a cos(α − β ) + β ˙ 2 b − γ ˙ 2 c cos(γ − β ) αa α ¨ = γ γ c sin(γ − β )
Ricorrendo ancora alle regole di derivazione in cascata si ha ∂β ¨+ α ∂α ∂γ ¨ = ¨+ γ γ α ∂α
¨= β
∂ 2 β 2 α˙ ∂α 2 ∂ 2 γ 2 α˙ ∂α 2
E’ possibile quindi riscrivere le equazioni in accelerazione come a cos(α − γ ) − f c2 c + f b2 b cos(β − γ ) 2 a sin(α − γ ) ¨ ¨+ β = α α˙ b sin(γ − β ) b sin(γ − β ) a cos(α − β ) + f b2 b − f c2 c cos(γ − β ) 2 a sin(α − β ) ¨ = ¨+ γ γ α α˙ c sin(γ − β ) c sin(γ − β )
dove si riconoscono i coefficienti f b e f c determinati in precedenza (1.1 e 1.2 nella pagina ∂ 2 β ∂ 2 γ precedente) e si introducono i coefficienti f f b = ∂α 2 e f f c = ∂α 2 :
1.2
a cos(α − γ ) − f c2 c + f b2 b cos(β − γ ) f f b = b sin(γ − β )
(1.3)
a cos(α − β ) + f b2 b − f c2 c cos(γ − β ) f f c = c sin(γ − β )
(1.4)
Second Secondo o quadril quadrilate atero ro - Collegam Collegamen ento to barra barra di torsion torsione e
L’analisi L’analisi della cinematica cinematica per p er il secondo secondo quadrilatero quadrilatero risulta simile alla preceden precedente. te. Si mantiene l’angolo α come coordinata libera, tenendolo in evidenza grazie all’introduzione dell’angolo costante α che rappresenta l’inclinazione tra i braccetti superiori dei due quadrilateri. Pertanto si mantiene lo stesso sistema di riferimento Re-Imm . ∗
1.2. 1.2.1 1
Equa Eq uazio zione ne in in posi posizio zione ne
Anche in questo caso, si scrive l’equazione di chiusura del quadrilatero: a2 ei(α
α )
−
∗
+ b2 eiβ2 = c2 eiγ 2 + d2 eiδ2
Proiettando sui due assi reale e immaginario si ha:
b2 cos β 2 = −a2 cos(α − α ) + c2 cos γ 2 + d2 cos δ 2 b2 sin β 2 = −a2 sin(α − α ) + c2 sin γ 2 + d2 sin δ 2 ∗
∗
Per la determinazione degli angoli β 2 e γ 2 si sfrutta il calcolo numerico mediante il software Matlab.
9
CAPITOLO 1. RISOLUZIONE RISOLUZIONE DELLA CINEMATICA CINEMATICA
Figura 1.2: Chiusura del secondo quadrilatero articolato 1.2. 1.2.2 2
Equa Eq uazio zione ne in in veloc velocit ità à
Derivando rispetto al tempo i termini dell’equazione in posizione si ottiene: ˙ 2 ei(α iαa
α ) ∗
−
+ iβ ˙2 b2 eiβ2 = iγ ˙2 c2 eiγ 2
moltiplicando l’equazione trovata una volta per e
iβ2
−
˙ 2 ei(α iαa ˙ 2 ei(α iαa
α
∗
α
∗
−
−
e una volta per e
iγ 2
−
si ottiene:
+ iβ ˙2 b2 − iγ ˙2 c2 ei(γ 2 β2 ) = 0 γ 2 ) + iβ ˙2 b2 ei(β2 γ 2 ) − iγ ˙2 c2 = 0 β2 )
−
−
−
−
Proiettando sull’asse reale:
˙ 2 sin(α − α αa ˙ 2 sin(α − α αa
∗
∗
β 2 ) − γ ˙2 c2 sin(γ 2 − β 2 ) = 0 ˙2 b2 sin(β 2 − γ 2 ) = 0 − γ 2 ) + β −
da cui a2 sin(α − α − γ 2 ) β ˙2 = α˙ b2 sin(γ 2 − β 2 ) a2 sin(α − α − β 2 ) γ ˙2 = α˙ c2 sin(γ 2 − β 2 ) ∗
∗
Poiché β 2 = β 2 (α) e γ 2 = γ 2(α), applicando la regola di derivazione a cascata si ha 2 2 2 2 ˙ ˙ = ∂γ ˙ = ∂β = ∂γ β ˙2 = ∂β ∂α α e γ 2 ∂α α, quindi si possono esplicitare i coefficienti f b2 ∂α e f c2 ∂α
10
CAPITOLO 1. RISOLUZIONE RISOLUZIONE DELLA CINEMATICA CINEMATICA
che risultano costanti ad ogni data posizione α: a2 sin(α − α − γ 2 ) b2 sin(γ 2 − β 2 )
(1.5)
a2 sin(α − α − β 2 ) f c2 = c2 sin(γ 2 − β 2 )
(1.6)
∗
f b2 =
∗
1.2.3 1.2.3
Equazi Equ azione one in acceler accelerazio azione ne
Derivando una seconda volta l’equazione di chiusura rispetto al tempo si ha: ¨ a2 ei(α iαa α
α )
−
∗
−
α˙ 2 a2 ei(α
α ) ∗
−
+ iβ ¨2 b2 eiβ2
come in precedenza si moltiplica una volta per e
¨ a2 ei(α iαa α
α
−
¨ a2 ei(α iαa α
α
−
−
−
∗
∗
−
α˙ 2 a2 ei(α
α
−
γ 2 )
−
α˙ 2 a2 ei(α
α
−
∗
−
∗
iβ2
−
−
e una volta per e
2 + iβ ¨2 b2 − β ˙2 b2 − iγ ¨2 c2 ei(γ 2 2 γ 2 ) + iβ ¨2 b2 ei(β2 γ 2 ) − β ˙2 b2 ei(β2
β2 )
−
2
¨2 c2 eiγ 2 β ˙2 b2 eiβ2 = iγ
−
β2 )
−
iγ 2
−
:
+ γ ˙2 2 c2 ei(γ 2 β2 ) = 0 γ 2 ) ¨2 c2 + γ ˙2 2 c2 = 0 − iγ
β2 )
−
γ ˙2 2 c2 eiγ 2
−
−
Proiettando sull’asse reale si determinano le accelerazioni angolari: ¨2 = β
¨ 2 sin(α − α αa
−
¨ a2 sin(α − α αa α
−
∗
∗
¨2 = γ
2
γ 2 ) + α˙ 2 a2 cos(α − α − γ 2 ) + β ˙2 b2 cos(β 2 − γ 2 ) − γ ˙2 2 c2 b2 sin(γ 2 − β 2 ) ∗
2
β 2 ) + α˙ 2 a2 cos(α − α − β 2 ) + β ˙2 b2 − γ ˙2 2 c2 cos(γ 2 − β 2 ) c2 sin(γ 2 − β 2 ) ∗
Ricorrendo ancora alle regole di derivazione in cascata si ha ∂β 2 ∂ 2 β 2 2 ¨ ¨+ β 2 = α α˙ 2 ∂α ∂α ∂γ 2 ∂ 2 γ 2 2 ¨2 = ¨+ γ α α˙ ∂α ∂α 2
E possibile quindi riscrivere le equazioni in accelerazione come a2 cos(α − α a2 sin(α − α − γ 2 ) ¨+ α b2 sin(γ 2 − β 2 )
−
a2 cos(α − α a2 sin(α − α − β 2 ) ¨2 = ¨+ γ α c2 sin(γ 2 − β 2 )
−
¨2 = β
∗
∗
∗
∗
γ 2 ) − f c22 c2 + f b22 b2 cos(β 2 − γ 2 ) 2 α˙ b2 sin(γ 2 − β 2 ) β 2 ) + f b22 b2 − f c22 c2 cos(γ 2 − β 2 ) 2 α˙ c2 sin(γ 2 − β 2 )
dove si riconoscono i coefficienti f b2 f c2 determinati in precedenza (1.5 e 1.6) e si intro2 2 ducono i coefficienti f f b2 = ∂ ∂αβ22 e f f c2 = ∂ ∂αγ 22 : a2 cos(α − α
−
a2 cos(α − α
−
∗
f f b2 = f f c2 =
∗
γ 2 ) − f c22 c2 + f b22 b2 cos(β 2 − γ 2 ) b2 sin(γ 2 − β 2 )
(1.7)
β 2 ) + f b22 b2 − f c22 c2 cos(γ 2 − β 2 ) c2 sin(γ 2 − β 2 )
(1.8)
Capitolo 2
Scrittura dell’equazione di moto L’equazione di moto verrà ricavata mediante l’equazione di Lagrange, scritta rispetto alla coordinata libera scelta, ovvero l’angolo α indicato in figura: d dt
∂E c ∂ α˙
−
∂E c ∂V ∂D + + = Qα . ∂α ∂α ∂ α˙
Con E c si indica l’energia cinetica totale, con V e D rispettivamente l’energia potenziale totale e quella dissipativa totale; Qα è la componente lagrangiana delle forze esterne rispetto rispetto alla coordinata libera α. Si nota fin da subito che non sono applicate forze esterne non conservative: si può quindi scrivere: Qα = 0 Convenzioni di segno
Si considerano positivi gli spostamenti verso destra, verso l’alto e le rotazioni antiorarie. Gli allungamenti degli elementi elastici sono stati considerati positivi se di trazione, negativi se di compressione.
2.1 2.1
Ener Energi gia a cine cineti tica ca
L’energia cinetica del sistema è la somma delle energie cinetiche dei vari corpi: 1 2
1 2
1 2
1 2
2 2 2 2 + J A ωA + mB vG + J B ωB + E c = mA vG A B
1 1 2 2 + mC vG J C C ωC + C 2 2
1 1 1 1 2 2 2 2 + mB 2 vG + + + J ω m v J C B 2 C 2 C 2 ωC 2 B2 GC 2 B2 2 2 2 2
Si procede ricavando ricavando le velocità assolute assolute dei baricent baricentri. ri. Per il corpo A l’equazione l’equazione di chiusura è: dA ei(α+α0 ) = xGA + iyGA
11
12
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
(a) Quadrilatero principale
(b) Secondo quadrilatero
Figura Figura 2.1: Posizion Posizionee generica generica dei baricent baricentri ri per p er i due quadrilateri quadrilateri Proiettando sugli assi:
Derivando:
xGA = dA cos(α + α0 ) yGA = dA sin(α + α0 ) x˙ GA = −dA sin(α + α0 ) y˙GA = dA cos(α + α0 )
Elevando al quadrato e sommando: 2 2 2 2 2 = x˙ G + y˙ G = dA vG α˙ A A A
Equazione di chiusura per il corpo B: aeiα + dB ei(β +β0 ) = xGB + iyGB
Proiettando sugli assi:
xGB = a cos α + dB cos(β + β 0 ) yGB = a sin α + dB sin(β + β 0 )
Derivando rispetto al tempo:
Da cui:
˙ x˙ GB = −a sin αα˙ − dB sin(β + β 0 )β ˙ y˙ GB = a cos αα˙ + dB cos(β + β 0 )β
2 2 2 2 2 = x˙ G + y˙G = a2 + dB vG f b + 2adB f b cos(α − β − β 0 ) α˙ 2 . B B B
Per il corpo C, con ragionamento analogo a quello per vGA , e considerando che γ ˙ è la velocità angolare assoluta per l’asta considerata, si ottiene: ˙ vGC = dC γ
13
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
da cui:
2 2 2 2 = dC vG f C α˙ C
Equazione di chiusura per il corpo B2: a2 ei(α
α )
−
∗
+ dB2 ei(β2 +β02 ) = xGB2 + iyGB2
Proiettando sugli assi:
xGB2 = a2 cos(α − α ) + dB2 cos(β 2 + β 02 02 ) yGB2 = a2 sin(α − α ) + dB2 sin(β 2 + β 02 02 ) ∗
∗
Derivando rispetto al tempo:
Da cui: 2
2
˙ x˙ GB2 = −a2 sin(α − α )α˙ − dB 2 sin(β 2 + β 02 02 )β 2 ˙ y˙GB2 = a2 cos(α − α )α˙ + dB2 cos(β 2 + β 02 02 )β 2 ∗
∗
2
vGB2 = x˙ GB2 + y˙ GB2 =
a22
+
2
dB2 f b22
+ 2a2 dB2 f b2 cos((α − α ) − β 2 − β 02 ˙ 2. 02 ) α ∗
Per il corpo C2, considerando che γ ˙2 è la velocità angolare assoluta per l’asta considerata, si ottiene: vGC 2 = dC 2 γ ˙2
da cui:
2 2 2 = dC ˙2 vG 2 f C 2 α C 2
Passiamo infine al calcolo delle velocità angolari delle aste: 2 = α˙ 2 ωA 2 = β ˙ 2 = f b2 α˙ 2 ωB 2 = γ ˙ 2 = f c2 α˙ 2 ωC 2
2 = β ˙2 = f b22 α˙ 2 ωB 2 2 = γ ˙2 2 = f c22 α˙ 2 ωC 2
É possibile ora esplicitare l’espressione dell’energia cinetica:
1 1 2 + J A α˙ 2 = J A (α)α˙ 2 E cA = mA dA 2 2 J A (α)
1 1 2 2 E cB = mB a2 + dB f b + 2adB f b cos(α − β − β 0 ) +J B f b2 α˙ 2 = J B (α)α˙ 2 2 2
E cC
J B (α)
1 1 2 2 2 = ˙ 2 = J C ˙2 mC dC f c + J C C f c α C (α)α 2 2 J C C (α)
14
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
E cB 2
E cC 2
1 2 2 2 = ˙2 = mB 2 a22 + dB 02 ) +J B 2 f b2 α 2 f b2 + 2a2 dB 2 f b2 cos((α − α ) − β 2 − β 02 2 ∗
J B2 (α)
1 = J B2 (α)α˙ 2 2 1 1 2 2 2 = ˙ 2 = J C ˙2 mC 2 dC C 2 f c2 α C 2 (α)α 2 f c2 + J C 2 2
J C C 2 (α)
dove i termini J i (α) sono i momenti d’inerzia generalizzati dei corpi rispetto alla coordinata libera α. Si ha quindi J (α) = i J i (α) = J A (α) + J B (α) + J C C (α) + J B 2 (α) + J C C 2 (α). Il contributo dell’energia cinetica all’equazione di Lagrange è dato dai seguenti termini: d dt
∂E c ∂ α˙
−
1 ∂J (α) 2 ∂E c = J (α)α ¨+ α˙ 2 ∂α ∂α
(2.1)
occorre pertanto esplicitare la derivata parziale dei momenti d’inerzia generalizzati, fatta ∂ 2 β ∂ 2 γ rispetto ad α. Ricordando la definizione dei coefficienti f f b = ∂α 2 e f f c = ∂α 2 si ottiene: ∂J A (α) =0 ∂α ∂J B (α) 2 = mB 2dB f b f f b + 2adB f f b cos(α − β − β 0 ) − 2adB f b sin(α − β − β 0 )(1 − f b ) + ∂α
+ 2J B f B f f B ∂J C C (α) = 2mC d2c f c f f c + 2J C C f C C f f C C ∂α ∂J B2 (α) 2 = mB2 2dB 2 f b2 f f b2 + 2a2 dB 2 f f b2 cos(α − α ∂α
−
∗
2a2 dB 2 f b2 sin(α − α
∗
−
2.2
∂J (α) ∂α
=
i
∂J i (α) ∂α
β 2 − β 02 02 )+
β 2 − β 02 02 )(1 − f b2 ) +2J B 2 f B 2 f f B 2
∂J C C 2 (α) = 2mC 2 d2c2 f c2 f f c2 + 2J C C 2 f C C 2 f f C C 2 ∂α
Si può scrivere infine
−
.
Energi Energia a poten potenzia ziale le gra gravit vitazi aziona onale le
Si considereranno separatamente i due contributi (elastico e gravitazionale) all’energia potenziale del sistema (V = V k + V g ). L’equazio L’equazione ne che esprime esprime l’energia l’energia potenziale potenziale gravitazionale è: V g = mA gh GA + mB gh GB + mC gh GC + mB2 gh GB2 + mC 2 gh GC 2
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
(a) Quadrilatero principale
(b) Secondo quadrilatero
Figura 2.2: Quote dei baricentri rispetto ad un asse orizzontale fisso Le altezze dei baricentri, calcolate rispetto ad un asse orizzontale fisso, risultano: hGA = h1 − dA cos(θ + α + α0 ) hGB = h1 − a cos(θ + α) − dB cos(θ + β + β 0 ) hGC = h2 − dC cos(θ + γ + γ 0 ) hGB2 = h1 − a2 cos(θ + α − α ) − dB2 cos(θ + β 2 + β 02 02 ) ∗
hGC 2 = h3 − dC 2 cos(θ + γ 2 + γ 02 02 )
le cui derivate rispetto ad α sono: ∂h GA = dA sin(θ + α + α0 ) ∂α ∂h GB = a sin(θ + α) + dB sin(θ + β + β 0 )f b ∂α
15
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
(a) Schem Schema a comple completo to della della sospens sospension ionee
16
(b) Chiusu Chiusura ra usata usata per l’ana l’analis lisii del sistem sistema a mollamollasmorzatore a terra
Figura 2.3: Rappresentazione semplificata della sospensione ∂h GC = dC sin(θ + γ + γ 0 )f c ∂α ∂h GB2 = a2 sin(θ + α − α ) + dB2 sin(θ + β 2 + β 02 02 )f b2 ∂α ∂h GC 2 = dC 2 sin(θ + γ 2 + γ 02 02 )f c2 ∂α ∗
Derivando l’energia potenziale rispetto alla coordinata libera si ottiene il termine utile all’equazione di Lagrange: ∂V g ∂h Ga ∂h Gb ∂h Gc ∂h GB2 ∂h GC 2 = mA g + mB g + mC g + mB 2 g + mC 2 g ∂α ∂α ∂α ∂α ∂α ∂α
2.3
(2.2)
Energia Energia potenzi potenziale ale e dissipati dissipativ va del del sistema sistema molla-smor molla-smorzato zatore re a terra
Ai fini energetici, si rappresenta il pneumatico con un sistema di molla e smorzatore (figura (figura 2.3). Tali elementi elementi vengono vengono definiti con il pedice 1, mentre z è lo spostamento del vincolo a terra. Si osserva che il vincolamento a terra avviene mediante un pattino, quindi l’allungamento della molla è dato dalla sola traslazione verticale degli estremi. L’allungamento della molla 1 è così definito: ∆l1 = ∆l01 + ∆ld1 (α) = ∆l01 + lF (α) − lI (αs )
Dove il pedice s rappresenta il valore degli angoli calcolati nella posizione di equilibrio statico. Il precarico statico ∆l01 è noto e si ottiene considerando che la molla in condizione
17
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
statiche equilibria un quarto del peso della vettura, più il peso della sospensione stessa: ∆l01 = −
0.25P auto auto + P sosp sosp k1
Si determina determina quindi la lunghezza lunghezza iniziale iniziale e finale: finale: )] − z lF (α) = [h1 − a cos(θ + α) − b cos(θ + β + β )] ∗
∗
lI (αs ) = [h1 − a cos(θ + αs ) − b cos(θ + β + β s )] ∗
∗
Sostituendo si ottiene: ∆ld1 = −a cos(θ + α) − b cos(θ + β + β ) + a cos(θ + αs ) + b cos(θ + β + β s ) − z ∗
∗
∗
∗
∆l1 = ∆l01 − a cos(θ + α) − b cos(θ + β + β ) + a cos(θ + αs ) + b cos(θ + β + β s ) − z ∗
∗
∗
∗
Si ha quindi un contributo all’energia potenziale pari a: V k1 =
1 k1 ∆l12 2 ∂ ∆l1 ∂α
L’equazione di Lagrange richiede il calcolo del termine
:
∂ ∆l1 = a sin(θ + α) + b sin(θ + β + β )f b ∂α ∗
∗
Il termine utile all’equazione di moto risulta pertanto: ∂V k1 ∂ ∆l1 = k1 ∆l1 ∂α ∂α
Molla e smorzatore sono in parallelo, quindi al fine di determinare la quota parte di dissipazione dello smorzatore 1 sarà sufficiente derivare rispetto al tempo l’allungamento già trovato. Noto ∆l1 = ∆l1 (α(t), z (t)), si ottiene: d∆l1 ∂ ∆l1 ∂ ∆l1 ∆˙l1 = = α˙ + z˙ dt
Determinato 2 ∆˙l1 =
∂ ∆l1 ∂z
∂α
∂z
= −1, si ricava il termine
∂ ∆l1 ∂α
2
∂ ∆l1 ∂ ∆l1 α˙ + 2 α˙ z˙ + ∂α ∂z 2
∂ ∆l1 ∂z
2
z˙ 2
da cui la funzione dissipativa
1 1 ∂ ∆l1 2 D1 = r1 ∆˙l1 = r1 2 2 ∂α
2
ed il termine richiesto dall’equazione di Lagrange:
∂D 1 ∂ ∆l1 = r1 ∂ α˙ ∂α
2
α˙ + r1
1 ∂ ∆l1 ∂ ∆l1 ∂ ∆l1 α˙ 2 + r1 α˙ z˙ + r1 2 ∂α ∂z ∂z
∂ ∆l1 ∂ ∆l1 z˙ ∂α ∂z
2
z˙ 2
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
2.4
18
Energi Energia a potenzia potenziale le elasti elastica ca della della molla molla torsio torsional nale e
Si conclude l’analisi dell’energia potenziale inserendo il contributo della seconda molla, rappresentante la rigidezza k2 della barra torsionale. torsionale. Detto ∆γ l’allungamento angolare della molla: V k2 =
1 K 2 ∆γ 22 2
Si procede quindi esplicitando l’allungamento della molla come: ∆γ 2 = ∆γ 20 20 + γ 2 (α) − γ 2 (αs )
dove compare il precarico statico incognito ∆γ 20 20 , che potrà essere calcolato inponendo ∂V = 0, essendo noto il precarico relativo alla molla k 1. ∂α Per ottenere il termine relativo all’equazione di Lagrange occorre valutare la derivata parziale dell’allungamento rispetto alla coordinata libera: ∂ ∆γ 2 ∂γ 2 = = f c2 ∂α ∂α
Il termine utile all’equazione di moto risulta pertanto: ∂V k2 ∂ ∆γ 2 = k2 ∆γ 2 ∂α ∂α
2.5
Energi Energia a dissip dissipati ativ va dello dello smorz smorzato atore re
L’energia dissipativa associata al vero e proprio smorzatore della sospensione è espressa 2 come: D2 = 12 r2 ∆˙l2 . La variazione della lunghezza dello smorzatore nel tempo viene calcolata mediante il teorema di Pitagora: l2 (α) =
(xR2 − xR1 )2 + (yR2 − yR1 )2
in cui xR1 e yR1 sono noti e costanti; xR2 e yR2 si calcolano invece come:
xR2 = a sin(θ + α) + b00 sin(θ + β + β 00 00 ) yR2 = −a cos(θ + α) − b00 cos(θ + β + β 00 00 )
mentre le loro derivate parziali rispetto ad α valgono:
∂x R2 ∂α ∂y R2 ∂α
= a cos(θ + α) + b00 cos(θ + β + β 00 00 )f b = a sin(θ + α) + b00 sin(θ + β + β 00 00 )f b
Ora, ∆˙l2 può essere scritta come
∂l 2 ˙ ∂α α
, per cui:
1 ∂l 2 = (xR2 − xR1 )2 + (yR2 − yR1 )2 2 ∂α
1 2
−
∂x R2 ∂y R2 2(xR2 − xR1 ) + 2(yR2 − yR1 ) ∂α ∂α
19
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
(a) Risolu Risoluzio zione ne con equ equazi azioni oni di chius chiusura ura
(b) Risolu Risoluzio zione ne median mediante te Teorema eorema di Pitago Pitagora ra
Figura 2.4: Metodi di calcolo di ∆l2 ed infine:
1 1 ∂ ∆l2 D2 = r2 ∆˙l22 = r2 2 2 ∂α
2
2
α˙
⇒
∂D 2 ∂ ∆l2 = r2 ∂ α˙ ∂α
2
α˙
Osservazione
Si propone un’altro metodo per valutare ∆˙l2 mediante un’equazione di chiusura. π
yc eiπ + xc ei 2 + ∆l2 ei = aeiα + b00 ei(β +β00 )
Si proietta sugli assi ottenedo il seguente sistema:
yc + ∆l2 cos = a cos α + b00 cos(β + β 00 00 ) xc + ∆l2 sin = a sin α + b00 sin(β + β 00 00 )
−
20
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
Da questo questo sistema sistema grazie all’utilizzo all’utilizzo della funzione funzione fsolve di Matlab è possibile ricavare ∆l2 (α) e (α). Si passa quindi alla soluzione in velocità: ˙ 00 ei(β +β00 ) ∆˙l2 ei + i˙∆l2 ei = iαae ˙ iα + iβb
moltiplicando per e
i
−
e proiettand proiettandoo sull’asse sull’asse reale si ottiene ottiene ∆˙l2 :
∆˙l2 = [−a sin(α − ) − b00 f b sin(β + β 00 ˙ = f ∆l2 α˙ 00 − )]α
dove f ∆l2 = [−a sin(α − ) − b00 f b sin(β + β 00 00 − )] =
1 ∂ ∆l2 D2 = r2 2 ∂α
2
α˙ 2 =
∂ ∆l2 ∂α
. Si può quindi scrivere
1 r2 f 2 α˙ 2 2 ∆l2
da cui, derivando rispetto ad α˙ , si ottiene il termine richiesto dall’equazione di Lagrange: ∂D 2 2 = r2 f ∆ ˙ l2 α ∂ α˙
2.6 2.6
Equa Equazi zion one e di mot moto o e form forma a di sta stato to
Nell’applicare l’equazione di Lagrange, è ora sufficiente sommare i termini evidenziati nelle sezioni precedenti per ottenere l’equazione del moto in grande del sistema: d dt
∂E c ∂ α˙
−
∂E c ∂V g ∂V k1 ∂V k2 ∂D 1 ∂D 2 + + + + + = Qα ∂α ∂α ∂α ∂α ∂ α˙ ∂ α˙
(2.3)
Ricordando che Qα = 0 e sostituendo i singoli contributi si ottiene l’equazione di moto: ¨+ J (α)α
1 ∂J (α) 2 α˙ + 2 ∂α
∂h i ∂ ∆l1 ∂ ∆γ 2 + k1 ∆l1 + k2 ∆γ 2 + ∂α ∂α ∂α
mi g
i
+ r1
∂ ∆l1 ∂α
2
∂ ∆l1 ∂ ∆l1 ∂ ∆l2 α˙ + r1 z˙ + r2 ∂α ∂z ∂α
2
α˙ = 0
(2.4)
Isolando ora α¨ si può scrivere l’equazione nella forma di stato, come richiesto per la risoluzione numerica con il metodo di Runge-Kutta.
¨=− α α˙ = α˙
1 ∂J (α) 2 ˙ 2 ∂α α
+
∂V g ∂α
+
∂V k1 ∂α
+
∂V k2 ∂α
+
∂D 1 ∂ α ˙
+
∂D 2 ∂ α˙
1
J (α)
21
CAPITOL CAPITOLO O 2. SCRITTUR SCRITTURA A DELL’EQUAZ DELL’EQUAZIONE IONE DI MOTO
2.7 2.7
∆γ 20
Prec Pr ecar aric ico o stat static ico o
Il precarico statico della molla torsionale viene calcolato imponendo l’annullarsi della derivata dell’energia potenziale rispetto alla coordinata libera nella posizione di equilibrio statico. ∂V ∂α
αs
∂V g = ∂α
Da cui:
∂V k1 + ∂α αs
∂V k2 + ∂α αs
=−
∆l1 k1 ∆l01 ∂ ∂α −
γ 2 k2 ∂ ∆ ∂α
αs
+
αs
αs
αs
αs
i
∆l1 k1 ∆l01 ∂ ∂α
γ 2 k2 ∂ ∆ ∂α
mA g ∂h∂αGa
=
∂h i i mi g ∂α
∆γ 20 20 = −
αs
+ mB g ∂h∂αGb
αs
∂h i mi g ∂α
∂ ∆l1 +k1 ∆l01 ∂α αs
∂ ∆γ 2 +k2 ∆γ 20 20 ∂α αs
αs
αs
+ mC g ∂h∂αGc
+
αs
γ 2 k2 ∂ ∆ ∂α
αs
∂h G mB 2 g ∂αB2
αs
+
∂h G mC 2 g ∂αC 2
αs
Capitolo 3
Equazione di moto linearizzata Si procede alla linearizzazione dell’equazione di moto nell’intorno della posizione di equilibrio statico rappresentata dal valore αs della coordinata libera. Prima di procedere, è opportuno evidenziare la necessità di considerare come coordinata indipendente nel processo di linearizzazione, l’angolo α¯ = α − αs che permette di descrivere il sistema in funzione di spostamenti rispetto alla configurazione di equilibrio. Questo permette di fatto la linearizzazione delle forme di energie tramite l’espansione in serie di Taylor.
¯ = α − αs α ¯˙ = α˙ α ¨¯ = α ¨ α
3.1 3.1
Ener Energi gia a cine cineti tica ca
Dallo sviluppo in serie dell’energia cinetica si può mostrare che, trascurando i termini di ordine superiore che impedirebbero di ottenere una forma quadratica, risulta: Ec
≈
1 ∂ 2 Ec 2 ∂ α˙ 2
e quindi d dt
3.2 3.2
∂E c ∂ α˙
−
¯˙ 2 = α
αs
1 ¯˙ 2 J (αs )α 2
∂E c ¨¯ = J (αs )α ∂α
Ener Energi gia a poten potenzi zial ale e
Ricordando che V = V (α, z ), si procede a sviluppare fino al secondo ordine l’energia potenziale rispetto ad entrambe le variabili. ∂V V ≈ V (αS , zs ) + ∂α
1 ∂ 2 V 2 ∂V ¯+ ¯ + α α 2 2 ∂α ∂z αs ,zs αs ,zs
22
1 ∂ 2 V ∂ 2 V 2 + z+ z 2 ∂z 2 αs ,zs ∂α∂z αs ,zs
¯ zα αs ,zs
23
CAPITOLO 3. EQUAZIONE DI MOTO LINEARIZZA LINEARIZZAT TA
Eliminando i termini nulli e quelli che non dipendono da α:
1 ∂ 2 V ∂ 2 V 2 ¯ + V ≈ + α 2 ∂α 2 αs ,zs ∂α∂z
Sviluppando i singoli termini:
∂ 2 V ∂ ∆l1 = k1 2 ∂α ∂α
2
¯ zα αs ,zs
∂ 2 ∆l1 ∂ ∆γ 2 +k1 ∆l1 + k2 2 ∂α ∂α
2
∂ 2 ∆γ 2 +k2 ∆γ 2 + ∂α 2
i
∂ 2 hi mi g ∂α 2
Il secondo termine diventa: ∂ 2 V ∂ ∆l1 ∂ ∆l1 ∂ ∆l1 ∂ ∆l1 = k1 = k1 (−1) = −k1 ∂α∂z ∂α ∂z ∂α ∂α
Si analizzano quindi le singole derivate da inserire: ∂ ∆l1 = a sin(θ + α) + b sin(θ + β + β )f b ∂α ∂ ∆2 l1 = a cos(θ + α) + b cos(θ + β + β )f b2 + b sin(θ + β + β )f f b ∂α 2 ∂ ∆γ 2 a2 sin(α − α − β 2 ) = f c2 = ∂α c2 sin(γ 2 − β 2 ) ∗
∗
∗
∗
∗
∗
∗
a2 cos(α − α ∂ ∆2 γ 2 = = f f c 2 ∂α 2
∗
−
β 2 ) + f b22 b2 − f c22 c2 cos(γ 2 − β 2 ) c2 sin(γ 2 − β 2 )
Per la parte di energia potenziale gravitazionale: ∂ 2 hGA = dA cos(θ + α + α0 ) ∂α 2 ∂ 2 hGB = a cos(θ + α) + dB cos(θ + β + β 0 )f b2 + dB sin(θ + β + β 0 )f f b ∂α 2 ∂ 2 hGC = dC cos(θ + γ + γ 0 )f c2 + dC sin(θ + γ + γ 0 )f f c 2 ∂α ∂ 2 hGB2 2 = a2 cos(θ + α − α ) + dB2 cos(θ + β 2 + β 02 02 )f b2 + dB 2 cos(θ + β 2 + β 02 02 )f f b2 2 ∂α ∂ 2 hGC 2 2 = dC 2 cos(θ + γ 2 + γ 02 02 )f c2 + dC 2 sin(θ + γ 2 + γ 02 02 )f f c2 2 ∂α ∗
24
CAPITOLO 3. EQUAZIONE DI MOTO LINEARIZZA LINEARIZZAT TA
3.3
Funzion unzione e dissip dissipati ativ va
L’energia dissipativa risulta funzione di 4 grandezze: il suo sviluppo in serie completo al secondo ordine sarebbe
1 ∂ 2 D ∂D ∂D ∂D ∂D ¯+ ¯˙ + ¯2+ D ≈D(αs , zs , α˙s , z˙s ) + α z+ α z˙ + α 2 2 ∂α s ∂α s ∂z s ∂ α˙ s ∂ z˙ s
1 ∂ 2 D 2 1 ∂ 2 D 1 ∂ 2 D 2 ∂ 2 D ∂ 2 D ∂ 2 D 2 + ¯˙ + ¯+ ¯α ¯˙ + ¯z+ z + α z˙ + zα α αz α 2 ∂z 2 s 2 ∂ α˙ 2 s 2 ∂ z˙ 2 s ∂α∂z s ∂α∂ α˙ s ∂α∂ z˙ s ∂ 2 D ∂ 2 D ∂ 2 D ˙ + ¯+ ¯˙ zα z z˙ + z˙ α ˙ z˙ s ∂z∂ ∂z ∂ α˙ s ∂z ∂ z˙ s ∂ α∂
dove con il pedice s si è indicato che le espressioni sono calcolate delle espressioni nelle condizioni condizioni di equilibri equilibrioo statico. statico. Consideran Considerando do solo i termini termini non nul nulli li che dipendono dipendono da α˙ si ottiene: D
≈
e quindi:
1 ∂ 2 D 2 ∂ α˙ 2
1 ∂ 2 D 2 ˙α ¯ + 2 ∂ z˙ 2 α˙s ,z˙s
∂D ∂ 2 D = ∂ α˙ ∂ α˙ 2
dove
¯˙ + α α˙s ,z˙s
∂ 2 D ˙ z˙ ∂ α∂
∂ 2 D 2 z˙ + ˙ z˙ ∂ α∂ α˙s ,z˙S
2
∂ ∆l2 ∂α
+r2
¯˙ z˙ α
α˙s ,z˙S
z˙
α˙s ,z˙s
∂ 2 D ∂ ∆l1 = r 1 ∂ α˙ 2 ∂α
2
∂ 2 D ∂ ∆l1 ∂ ∆l1 ∂ ∆l1 = r1 = −r1 ˙ z˙ ∂ α∂ ∂α ∂z ∂α
3.4 3.4
Equa Equazi zion one e di moto moto
L’equazione del moto linearizzata risulta:
∂ ∆l1 ∂α
¨¯ + r1 J (αs )α + k2
∂ ∆γ 2 ∂α
s
2
2
+r2
s
∂ ∆l2 ∂α
∂ 2 ∆γ 2 +k2 ∆γ 20 20 ∂α 2
2
s
+
s
mi g
i
∂ ∆l1 ∂ ∆l1 ¯˙ − r1 α z˙ + k1 ∂α s ∂α ∂ 2 h
i
∂α 2
¯ − k1 α
s
s
2
∂ 2 ∆l1 +k1 ∆l01 + ∂α 2 s
∂ ∆l1 z=0 ∂α s
Esplicitando α¨ e mettendo a sistema l’identità α¯˙ = α¯˙ si può scrivere l’equazione nella forma di stato, come richiesto per la risoluzione numerica con il metodo di Runge-Kutta.
Capitolo 4
Analisi del comportamento del sistema La risposta del sistema in esame è valutata in varie condizioni di forzamento e perturbazione iniziale. Si confrontano in particolare il comportamento del sistema non lineare con quello linearizzato, osservando la variazione nel tempo della posizione angolare associata alla coordinata libera α utilizzata nel corso della trattazione. Nella caratterizzazione del sistema sono stati utilizzati i dati numerici forniti, inseriti in unità di misura SI. L’angolo α nelle condizioni di equilibrio risulta pari a 106 deg.
4.1
Frequen requenza za propri propria a
La frequenza propria del sistema in analisi viene calcolata a partire dall’equazione linearizzata nearizzata del sistema considerato considerato non smorzato smorzato e non forzato. Tale equazione equazione risulta essere:
∂ ∆l1 ∂α
¨ + k1 J (αs )α + k2
∂ ∆γ 2 ∂α
s
2
2
s
∂ 2 ∆l1 +k1 ∆l01 + ∂α 2 s
∂ 2 ∆γ 2 +k2 ∆γ 20 + 20 ∂α 2 s
i
∂ 2 hi mi g ∂α 2
α=0
s
Denominando keq (αs ) il termine racchiuso nelle parentesi quadre si ottiene le scrittura semplificata: ¨ + keq (αs )α = 0 J (αs )α
Per risolvere l’equazione differenziale omogenea imponiamo la soluzione α(t) = Aeiω0 t :
−
J (αs )ω0 + keq (αs ) Aeiω0 t = 0
25
26
CAPITOLO 4. 4. ANALISI DEL COMPOR COMPORTAMENTO DEL SISTEMA SISTEMA
Da cui si ottiene (tramite Matlab): ω0 =
f 0 =
ω0 = 8, 1402H z 2π
4.2 4.2
keq (αs ) rad = 51, 1465 J (αs ) s
Moto Moto libe libero ro
4.2.1 4.2.1
Piccole Piccole perturba perturbazio zioni ni
109 non lineare linearizzato ] 108 g e d [ α
a 107 l l e v o n a m106 a l l e d e 105 n o i z a t o r
104
103
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
4.5
5
deg Figura Figura 4.1: Risposta Risposta alla condizione condizione iniziale (α − αs ) = +3 deg
Si comincia sollecitando il sistema senza forzamento e con condizioni iniziali prossime deg). Si osserva correttamente (Figura alla condizione di equilibrio statico ( (α − αs ) = +3 deg 4.1) che in questo caso le risposte lineare e non lineare tendono a coincidere. 4.2.2 4.2.2
Pertu Perturba rbazion zionii di media media ent entità ità
Aumentando la sollecitazione iniziale ((α − αs ) = −26deg) si nota giustamente come la risposta non lineare si discosti in ampiezza da quella linearizzata (Figura 4.2). Inoltre gli effetti delle non linearità determinano una variazione della frequenza di oscillazione nel tempo, mentre correttamente la risposta linearizzata mantiene costante la pulsazione.
27
CAPITOLO 4. 4. ANALISI DEL COMPOR COMPORTAMENTO DEL SISTEMA SISTEMA
130 non lineare linearizzato
125 ] g 120 e d [ α
115
a l l e v 110 o n a m105 a l l e d 100 e n o 95 i z a t o 90 r
85 80
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
4.5
5
Figura 4.2: Risposta Risposta alla condizione condizione iniziale (α − αs ) = −26deg 180 non lineare linearizzato
160 ] g 140 e d [ α
a 120 l l e v o n 100 a m a l l 80 e d e n 60 o i z a t o 40 r
20 0
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
4.5
5
Figura 4.3: Risposta Risposta alla condizione condizione iniziale (α − αs ) = −65deg 4.2.3 4.2.3
Grandi Grandi perturb perturbazi azioni oni
Imponendo uno spostamento angolare iniziale di elevata entità, si raggiunge (Figura 4.3) una condizione condizione tale per cui la derivata derivata seconda seconda dell’energ dell’energia ia potenziale potenziale divent diventaa negativ negativaa
28
CAPITOLO 4. 4. ANALISI DEL COMPOR COMPORTAMENTO DEL SISTEMA SISTEMA
e quindi il sistema (non lineare) diverge verso una nuova condizione di equilibrio ( α = 18deg). La risposta linearizzata non risente di tale fenomeno e pertanto mantiene inalterato il suo andamento, andamento, portandosi portandosi sull’inizial sull’inizialee condizione condizione di equilibrio equilibrio.. Ciò mostra come per elevate perturbazioni la linearizzazione non è più accettabile. Si nota infine che l’andamento della risposta non lineare non è più sinusoidale, ed in particolare le frequenze di oscillazione non sono più confrontabili con la frequenza propria del sistema linearizzato.
4.3 4.3
Moto Moto forz forzat ato o - zona zona quasi quasi-s -sta tati tica ca
120 non lineare linearizzato ] g 115 e d [ α
a l l e v 110 o n a m a l l e d 105 e n o i z a t o 100 r
95
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
4.5
5
Figura 4.4: Risposta in condizioni quasi-statiche Applicando una forzante con pulsazione Ω = 10 rad/s e ampiezza 0.03 metri, si ottiene la risposta in figura 4.4. E’ stata imposta una perturbazione iniziale piccola: questo ha genera generato to rispost rispostee simili simili nei due casi (lineare (lineare e non lineare) lineare).. Si osserv osservaa che che nei primi primi istanti la risposta transitoria si somma all’integrale particolare, il quale poi permane a regime. Nel caso si applicassero forzanti di ampiezza maggiore, le due risposte diventerebbero sensibilmente sensibilmente differenti.
4.4 4.4
Moto Moto for forza zato to - zon zona a di riso risona nanz nza a
Le condizioni di risonanza si ottengono imponendo una pulsazione del vincolo pari alla pulsazione pulsazione propria del sistema linearizza linearizzato to (Ω = 51.1465 rad/s) e ampiezza 0.01 m. La
29
CAPITOLO 4. 4. ANALISI DEL COMPOR COMPORTAMENTO DEL SISTEMA SISTEMA
150 non lineare linearizzato
140 ] g e 130 d [ α
a 120 l l e v o n 110 a m a 100 l l e d e 90 n o i z a 80 t o r
70 60
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
4.5
5
Figura 4.5: Risposta in condizioni di risonanza (Ω = ω0 = 51.1465 rad/s) risposta linearizzata ottenuta (Figura 4.5) presenta un aumento dell’ampiezza di oscillazione lazione fino al valore massimo, massimo, corrisponde corrispondente nte a circa 40 gradi. Correttame Correttamente nte tale ampiezza risulta finita, grazie alla presenza di elementi smorzatori nel sistema. La risposta reale (non lineare) coincide con la precedente nei primi istanti, in cui le ampiezze di oscillazione sono piccole; successivamente gli effetti delle non linearità determinano uno scostamento dalla soluzione approssimata. L’oscillazione a regime (anch’essa alla stessa frequenza della forzante) è caratterizzata da ampiezze inferiori (circa 28 gradi).
4.5 4.5
Moto Moto for forza zato to - zon zona a sism sismog ogra rafic fica a
Si impone in questo caso un forzamento a Ω = 300 rad/s e ampiezza 0.03 m, con condizioni iniziali mediamente perturbate. Terminato il transitorio, in entrambi i casi (lineare e non lineare) le oscillazioni a regime permangono ma le ampiezze risultano correttamente molto minori del caso quasi-statico (Figura 4.6). Tale riduzione aumenta all’aumentare della frequenza della forzante.
30
CAPITOLO 4. 4. ANALISI DEL COMPOR COMPORTAMENTO DEL SISTEMA SISTEMA
140 non lineare linearizzato ] 130 g e d [ α
a l l 120 e v o n a m110 a l l e d e 100 n o i z a t o r
90
80
0
0.5
1
1.5
2
2.5 tempo [s]
3
3.5
4
Figura 4.6: Risposta in condizioni sismografiche
4.5
5
Capitolo 5
Scripts del programma Matlab Gli scripts originari forniti sono stati modificati inserendo tutti i risultati ottenuti nel corso della trattazione, realizzando quindi un programma in grado di simulare il comportamento della sospensione in esame qualora siano note tutte le condizioni geometriche, inerziali, elastiche e viscose. Ovviamente molti dei parametri e delle grandezze introdotte risulteranno nulli in questo caso, data la scarsità di informazioni riguardo le caratteristiche inerziali degli elementi costituenti. Per l’analisi dello smorzatore, è stato implementato il primo dei due metodi proposti, e sono state usate in tutto il programma unità di misura SI. Per la risoluzione numerica delle equazioni differenziali è stato preferito l’utilizzo della funzione ode45 di Matlab alla rk4 fornita: anch’essa utilizza un metodo di Runge-Kutta del quarto ordine (insieme ad uno del quinto ordine usato per stimare l’errore ad ogni passo), ma risulta più efficiente e meglio integrata nell’ambiente Matlab. Infine sono state apportate alcune modifiche alla funzione main.m volte a rendere più semplice il cambiamento delle condizioni di perturbazione ad ogni esecuzione, e la visualizzazione delle singole risposte nei casi lineare e non lineare.
5.1 1 2 3 4 5 6 7 8 9 10 11 12 13 14
main.m
c le l e ar ar a l l c lo l o se se a l l clc s e t (0 , ’ Defa ultF igur eWin dowS tyle ’ , ’ Docked Docked ’ ) ;
w a rn rn in in g o f f global global global global global global global global
g a b c d a2 b2 c2 d2 bst b00 %geometria d a db d c d b2 b2 d c2 c2 %barice ntr i (lunghezze ) al fa 0 beta0 gamma0 beta02 gam gamma02 % b a r i c e n t r i ( a n g o l i ) a l f a s t t et e t a be b e ta t a st s t be be t a 0 0 d el e l ta ta 2 %a n g o li li c o s t a n t i ma J a mb mb J b mc mc J c mb2 J b2 b2 mc2 J c 2 k 1 k 2 r 1 r 2 z0 wz psiz a l f a s b et e t a s f b s f f b s ga g ammas f c s f f c s
31
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
15 global b e t as a s 2 f b 2 s f f b 2 s gammas2 f c 2 s f f c 2 s 16 global bet a_o ld gamm gamma_ a_ol old d bet a2_ old gamma gamma2_ 2_ol old d 17 global DL01 Dgamma20 Dgamma20 18 global x 1 y 1 19 20 % d a t i 21 g = 9 . 8 1 ; 22 a = 0.2254855; 23 b = 0.2640701; 24 c = 0 . 3 2 2 0 0 6 2 ; 25 d=0.2931100; 26 a2=0.1259365; 27 b2=0.2585440; 28 c2=0.078; 29 d2=0.3310906; 30 b s t = 0 . 2 4 6 6 0 8 5 ; 31 b00=0.2141454; 32 da=1; 33 db=0.2466085; 34 dc=1; 35 db2=1; 36 dc2=1; 37 a l f a 0 = 0; 0; 38 beta0 =58 =58∗ pi / 1 8 0 ; 39 gamma gamma0=0; 0=0; 40 beta02=0; 41 gamma0 gamma02=0 2=0 ; 42 a l f a s t =54∗ p i / 1 8 0 ; 43 b e t a s t =b =b e t a 0 ; 4 4 d e l t a 2 = 1∗ pi / 1 8 0 ; 45 t e t a = 34 34 0∗ pi / 1 8 0 ; % −20 46 ma ma = 0 ; 47 J a = 0; 0; 48 mb mb =40; 49 J b = 5; 5; 50 mc mc =0; 51 Jc=0; 52 mb2=0; mb2=0; 53 Jb2=0; 54 mc2=0; mc2=0; 55 Jc2=0; 56 k 1 = 1 e 5 ; % rigidezza 57 k2=10; 58 r1=1e2 r1=1e2 ; 59 r2 = 10; % smorzamento 60 x1=0.0559407; 61 y1=0.0782934; 62 DL01= DL01= − 0.023544; 63 beta00 =6∗ pi / 1 8 0 ; 64 65 % f o r z a m e n t o 66 disp ( ’FORZAM ’FORZAMENTO ENTO’’ ) 6 7 z 0 =input ( ’ A mp mp i e z z a d i o s c i l l a z i o n e d e l v i n c o l o ( m e tr t r i )= ’ ) ;
32
% ampiezza
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
68 wz =input ( ’ P u l s a z i o n e d e l v i n c o l o ( r i s o n a n z a : 5 1 . 1 4 6 5 r ad a d / s )= ’ ) ; % pulsazione 69 p s i z = input ( ’ F as a s e d e l v i n c o l o ( r a d )= )= ’ ) ; % fase 70 71 % e q u i l i b r i o 72 a l f a s =10 6∗ p i / 1 8 0 ; % an a n go g o lo lo d i e q u i l i b r i o 7 3 b e t a _ o l d = 2 0 ∗ pi / 1 8 0 ; 74 gamma_ol gamma_old=110 d=110 ∗ pi / 1 8 0 ; 7 5 b e t a 2 _ o l d = 3 5 8 ∗ pi / 1 8 0 ; %% −2 76 gamma2_old=86 gamma2_old=86∗ pi / 1 8 0 ; 77 [ b e t a s , f b s , f f b s , g am ammas , f c s , f f c s , b e t a s 2 , f b 2 s , f f b 2 s , g am ammas2 , f c 2 s , f f c 2 s ] = cinematica ( al fa s ) ; 78 Dgamma20 = pr ec ar ic o () ; 79 80 % D a t i d i s i m u l a z i o n e % c o n d i z i o n i i n i z i a l i [ VE VELOCITA ; POSIZIONE ] 81 % a l f a d i e q u i l i b r i o =106 g r ad ad i 82 x0= x0=zeros ( 2 , 1 ) ; 83 disp ( ’CONDI ’CONDIZI ZION ONII INI ZIA LI ’ ) 8 4 x 0 ( 2 ) =input ( ’ A ng ng ol ol o i n i z i a l e ( e q u i l i b r i o : 1 0 6 g r a d i )= ’ ) ; 8 5 x 0 ( 2 ) = x 0 ( 2 ) ∗ pi / 1 8 0 ; 8 6 x 0 ( 1 ) =input ( ’ V e l o c i ta t a a n g o la l a r e i n i z i a l e ( r a d/ d / s )= ’ ) ; 87 disp ( ’ Ca lc ol o . . . ’ ) ; 88 %x 0 = [ + 1 . 5 ; 6 0∗ p i / 1 8 0 ] ; % d u r at at a 89 T = 5 ; % pa p a ss s s o d i i n t e g ra ra z i o n e 9 0 d t = 5 e − 3; 91 a z i o ne n e d e l s is i s t em e m a n on l i n e a r e 92 % % I n t e g r az 93 %[ t_nl , y_nl ] = rk4 ( ’ manove llismoN ONlinea re ’ ,T, dt , x0 ) ; 94 [ t_nl , y_nl]= ode45 ( ’manovellismoNONlineare ’ ,[0 T] , x0) ; 95 % In I n t e g ra r a z i on o n e d e l s is i s te te m a l i n e a r i z z a t o 96 %[ t_l , y_l ] = rk4 ( ’ man ove lli smo lin ear e ’ ,T, dt , x0) ; 97 [ t_l , y_l]= ode45 ( ’ m a n o v e l l i s m o l i n e a r e ’ , [ 0 T ] , x 0 ) ; 98 99 % C o n f r o n t o 100 posiz _non _lin=y_nl _lin=y_nl ( : ,2 ) . ∗ 1 8 0 / pi ; 101 veloc _non_ lin=y_nl lin=y_nl ( : ,1 ) . ∗ 1 8 0 / pi ; 102 po siz _li n=y_ n=y_ll ( : ,2) . ∗ 1 8 0 / pi ; 103 ve lo c_ lin=y_ lin=y_ll ( : ,1) . ∗ 1 8 0 / pi ; 104 f i g u r e 105 plot ( t_nl , posiz_ non_lin , ’ b ’ , ’LineWidth ’ ,2) ; 106 ylabel ( ’ r o t a z i o n e d e l l a m a no no v el e l la l a \ a l p ha h a [ d eg eg ] ’ ) 107 xlabel ( ’tempo ’tempo [ s ] ’ ) 108 legend ( ’ n on on l i n e a r e ’ ) 109 f i g u r e 110 plot ( t_l , posiz _lin , ’ g ’ , ’LineWidth ’LineWidth ’ ,2 ) ; 111 ylabel ( ’ r o t a z i o n e d e l l a m a no no v el e l la l a \ a l p ha h a [ d eg eg ] ’ ) 112 xlabel ( ’tempo ’tempo [ s ] ’ ) 113 legend ( ’ l i n e a r i z z a t o ’ ) 114 f i g u r e 115 plot ( t_nl , posiz_ non_lin , ’ b ’ ) ; 116 hold on 117 plot ( t_l , posiz _lin , ’ g ’ ) ; 118 ylabel ( ’ r o t a z i o n e d e l l a m a no no v el e l la l a \ a l p ha h a [ d eg eg ] ’ )
33
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
119 120 121 122
34
xlabel ( ’tempo ’tempo [ s ] ’ ) legend ( ’ n on on l i n e a r e ’ , ’ l i n e a r i z z a t o ’ )
%%%%%% f re qu en za pr op ri a [M M,, dM dM ] = e n e r g i a c i n e t i c a ( a l f a s , b e t a s , f b s , f f b s , f c s , f f c s , b e t a s 2 , f b 2 s , f f b 2 s , fc2s , ff c2 s ) ;
123 124 124 [ dhga , ddhga , dhgb , ddhgb , dhgc , ddhgc , dhgb2 , ddhgb2 , dhgc2 , ddhgc 2 , DLdin1 , dDL1, dDL1, ddD ddDL1 ,Dga ,Dgamma2, dDgamma2, ddDgamma2] = e ne rg ia po t en zi al e ( al fa s , bet as , betas2 , gam gammas ,gamm ,gammas2 as2 , fbs , ffb s , fb2s , ffb 2s , fc s , ff cs , fc2 s , ff c 2 s ) ; 125 126 dDL2 = smorz ator e ( al fa s , betas , fb s ) ; 127 128 128 keq=(k1 ∗dDL1^2+k1 ∗ DL01 ∗ddDL1+k2 ∗dDgamma2^2+k2 ∗Dgamma20∗ddDgamma2+ma∗ g ∗ ddhga +mb∗ g ∗ ddhgb+mc∗ g ∗ ddhgc+mb2∗ g ∗ ddhgb2+mc2 ∗ g ∗ ddhgc2) ; 129 130 wpropria=sqrt ( keq /M) 131 132 f p r o p r i a =w =w p r o p r i a / ( 2 ∗ pi )
5.2 5.2
cine cinema mati tica ca.m .m
1 function [ beta , fb , ffb , gamma, fc , ffc , beta2 , fb2 , ffb 2 , gamma2, fc2 , f fc 2 ] = cinematica ( alfa ) 2 3 global a b c a 2 b 2 c 2 4 global a lf lf a_ a_ tm tm p a l f a s t 5 global bet a_o ld gamm gamma_ a_ol old d bet a2_ old gamma gamma2_ 2_ol old d 6 7 a lf lf a_ a_ tm tm p = a l f a ; 8 %%%%%%% 1 q u a d r i l a t e r o %%%%%%%%%%%%%%%%%% 9 o p t i o n s = o p t i m s e t ( ’ T o lF lF un un ’ , 1 e − 6, ’ Dis pla y ’ , ’ of f ’ ) ; 10 so l= fs ol ve ( @posiz ione , [ beta_old gam gamma_ol a_old d ] , opt ion s ) ; 11 beta= s o l ( 1 ) ; 12 gamma= s o l ( 2 ) ; 1 3 f b = a ∗ s i n ( a l f a −gamma) /(b ∗ s i n (gamma−beta ) ) ; 14 fc=a fc=a ∗ s i n ( a l f a −beta ) / ( c ∗ s i n (gamma−beta ) ) ; 1 5 f f b = ( a ∗ co s ( a l f a −gamma)− f c ^ 2∗ c+fb^2 ∗ b ∗ co s ( beta−gamma) ) /( b ∗ s i n (gamma−beta ) ) ; 16 f f c = (a (a ∗ co s ( a l f a −beta )+fb^2∗ b− f c ^ 2 ∗ c ∗ co s (gamma−beta ) ) /( c ∗ s i n (gamma−beta ) ) ; 1 7 b e t a _ o l d = beta ; 18 gamma_old= gamma_old=gamma ; 19 20 %%%%%%% 2 q u a d r i l a t e r o %%%%%%%%%%%%%%%%%% 21 o p t i o n s = o p t i m s e t ( ’ T o lF lF un un ’ , 1 e − 6, ’ Dis pla y ’ , ’ of f ’ ) ; 22 s o l 2 = f s o l v e ( @ p o s i z i o n e 2 , [ be be t a 2 _ o ld ld ga mm mma2_old ] , o p t i o n s ) ; 23 b e t a 2= 2= s o l 2 ( 1 ) ; 24 gam gamma2 ma2=so l 2 (2 ) ; 2 5 f b 2 = a 2 ∗ s i n ( a l f a − a l f a s t −gamma2) /(b2 ∗ s i n (gamma2−bet a2 ) ) ; 26 fc 2=a2 2=a2 ∗ s i n ( a l f a − a l f a s t −beta2 ) /( c2 ∗ s i n (gamma2−bet a2 ) ) ; 2 7 f f b 2 = ( a 2 ∗ co s ( a l f a − a l f a s t −gamma2)− f c 2 ^ 2∗ c2+fb2^2 ∗ b2 ∗ co s ( b e t a 2 −gamma2) ) /(b2 ∗ s i n (gamma2−bet a2 ) ) ;
35
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
28 f f c 2 =( a2 a2 ∗ co s ( a l f a − a l f a s t −beta2 )+fb2^2∗ b2− f c 2 ^ 2∗ c 2 ∗ co s (gamma2−beta2 ) ) /( c2 ∗ s i n (gamma2−bet a2 ) ) ; 29 b e t a 2 _ ol ol d = b e t a 2 ; 30 gamma2_old=gam gamma2_old=gamma2; ma2; 31 32 function f = p o s i z i o n e ( a ng ng ) 33 34 global a b c d 35 global alfa_tmp 36 37 f = [ a ∗ co s (alfa_tmp)+b ∗ co s ( a n g ( 1 ) )−c ∗ co s ( a n g ( 2 ) )−d ; 38 b ∗ s i n ( ang( 1) )+a )+a ∗ s i n (alfa_tmp) −c ∗ s i n (ang (2 ) ) ] ; 39 40 41 function f = p o s i z i o n e 2 ( a ng ng 2 ) 42 43 global a 2 b 2 c 2 d 2 d e lt lt a 2 44 global a lf lf a_ a_ tm tm p a l f a s t 45 46 f = [ a 2 ∗ co s (alfa_tmp − a l f a s t ) +b +b 2 ∗ co s ( a n g 2 ( 1 ) )−c 2 ∗ co s ( a n g 2 ( 2 ) )−d2 ∗ co s ( d e l t a 2 ) ; 47 b 2 ∗ s i n ( ang2( 1) )+a2 )+a2 ∗ s i n (alfa_tmp − a l f a s t )−c 2 ∗ s i n ( a n g 2 ( 2 ) )−d2 ∗ s i n ( d e l t a 2 ) ];
5.3 5.3 1 2 3 4 5 6
prec precar aric ico. o.m m
function Dgamma20 = pr ec ar ic o ()
g DL01 DL01 ma mb mc mb2 mc2 mc2 k1 k2 a l f a s b et e t a s b et e t as a s 2 ga g ammas ga gammas2 f b s f f b s f b 2s 2s f f b 2 s f c s f f c s ffc2s
global global global global
fc2s
7 8 [ dhga , ddhga , dhgb , ddhgb , dhgc , ddhgc , dhgb2 , ddhgb2 , dhgc2 , ddhgc 2 , DLdin1 , dDL1, dDL1, ddD ddDL1 ,Dga ,Dgamma2, dDgamma2, ddDgamma2] = e ne rg ia po t en zi al e ( al fa s , bet as , betas2 , gam gammas ,gamm ,gammas2 as2 , fbs , ffb s , fb2s , ffb 2s , fc s , ff cs , fc2 s , ff c 2 s ) ; 9 %[ dya , ddya ddya , dyb , ddyb ,DL,d ,DL,dD DL,ddD L,ddDL] = en er gi a_ po te nz ia le ( alf a0 , beta0 , dbeta0 , ddbe ta0 ) ; 10 11 Dgamm Dgamma20 a20 = −1 ∗( g ∗ (ma∗ dhga+mb∗ dhgb+mc∗ dhgc+mb2∗ dhgb2+mc2∗ dhgc2)+k1 ∗ DL01 ∗ dDL1 ) /(k2 ∗dDgamma2) ;
5.4
energi energiaci acinet netica ica.m .m
1 function [ M, M, dM dM ] = e n e r g i a c i n e t i c a ( a l f a , beta , fb , ffb , fc , ffc , beta2 , fb2 fb2 , ffb2 , fc2 , ffc 2 ) 2
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
3 4 5 6 7 8 9 10 11
36
global a d a d b d c d b 2 d c 2 a 2 b 2 c 2 b c global ma J a mb mb J b mc J c mb2 J b2 b2 mc2 J c 2 global b et e t a 0 b e ta t a 02 02 a l f a s t
Ma = ma ∗ da^2/4+Ja ; dM dMa = 0 ; Mb = mb ∗ (a^2+db^2 ∗ fb^2+2∗ a ∗ db ∗ f b ∗ co s ( a l f a −beta−be ta 0 ) )+Jb )+Jb ∗ fb ^2; dMb = mb ∗ ( 2 ∗ db^2 ∗ f b ∗ f f b + 2∗ a ∗ db ∗ f f b ∗ co s ( a l f a −beta−beta0 ) −2∗ a ∗ db ∗ f b ∗ s i n ( a l f a − beta−beta0 ) ∗ (1 − fb ) )+2∗ Jb ∗ f b ∗ f f b ;
12 13 Mc = mc ∗ dc^2∗ fc^2+Jc ∗ f c ^ 2 ; 14 dMc = 2 ∗ mc∗ dc^2 ∗ f c ∗ f f c + 2 ∗ Jc ∗ f c ∗ f f c ; 15 16 Mb2 = mb2 ∗ (a2^2+db2^2 ∗ fb2^2+2 ∗ a2 ∗ db 2 ∗ f b 2 ∗ co s ( ( a l f a − a l f a s t )−beta2 −bet a02 ) )+ )+ Jb 2 ∗ fb2 ^2; 17 dMb2 = mb2 ∗ ( 2 ∗ db2^2∗ f b 2 ∗ f f b 2 + 2∗ a2 ∗ db2 db 2 ∗ f f b 2 ∗ co s ( a l f a − a l f a s t −beta2 −beta0 2 ) −2∗ a2 ∗ db2 db 2 ∗ f b 2 ∗ s i n ( a l f a − a l f a s t −beta2 −beta02 ) ∗ (1 − fb 2 ) )+2∗ Jb 2 ∗ f b 2 ∗ f f b 2 ; 18 19 Mc2 = mc2 ∗ dc2^2∗ fc2^2+Jc2 ∗ f c 2 ^ 2 ; 20 dMc2 = 2 ∗ mc2 ∗ dc2^2 ∗ f c 2 ∗ f f c 2 + 2∗ J c 2 ∗ f c 2 ∗ f f c 2 ; 21 22 M = Ma+ Ma+Mb+Mc+ Mc+Mb2+Mc2 Mb2+Mc2 ; 23 dM = dMa+ dMa+dMb+dMc+dMb2+ dMc+dMb2+dMc dMc2 2;
5.5
energi energiapot apoten enzia ziale. le.m m
1 function [ dhga , ddhga , dhgb , ddhgb , dhgc , ddhgc , dhgb2 , ddhgb2 , dhgc2 , ddhgc2 , DLdin1 ,dDL ,dDL1 ,ddDL1,Dg ,ddDL1,Dga amma2, dDgamma2, ddDgamma2] = e ne rg ia po te nz ia le ( al fa , beta , beta2 , gamma,gam ,gamma2, fb , ff b , fb2 , ff b2 , fc , ff c , fc2 , f f c 2 ) 2 3 global a b s t d a d b d c a 2 d b2 b2 d c2 c2 4 global a l f a s b e t a s gammas2 t e t a a l f a 0 b e t a 0 b e t a 0 2 ga gamma0 ga gamma02 a l f a s t betast 5 6 dhga=da dhga=da ∗ s i n ( t e t a +a +a l f a + a l f a 0 ) ; 7 ddhga=da ddhga=da ∗ co s ( t e t a +a +a l f a + a l f a 0 ) ; 8 9 dhgb= dhgb=a ∗ s i n ( te ta+al fa )+db )+db∗ s i n ( t e t a +beta+beta0 ) ∗ f b ; 10 ddhgb= ddhgb=a ∗ co s ( te ta+al fa )+db )+db∗ co s ( t e t a +beta+beta0 ) ∗ ( fb ^2)+db ^2)+db ∗ s i n ( t e t a +beta+ beta0 ) ∗ f f b ; 11 12 dhgc=dc dhgc=dc ∗ s i n ( t e t a +gamma+gamma0) ∗ f c ; 13 ddhgc=dc ddhgc=dc ∗ co s ( t e t a +gamma+gamma0) ∗ ( f c ^2)+dc ^2)+dc ∗ s i n ( t e t a +gamma+gamma0) ∗ f f c ; 14 15 dhgb2=a2 dhgb2=a2 ∗ s i n ( t e t a +a +a l f a − a l f a s t ) +d +d b2 b2 ∗ s i n ( te ta+beta2+ ta+beta2+beta02 ) ∗ fb2 ; 16 ddhgb2=a2 ddhgb2=a2 ∗ co s ( t e t a +a +a l f a − a l f a s t )+ )+d b2 b2 ∗ co s ( te ta+beta2+ ta+beta2+beta02 ) ∗ ( fb2 ^2)+db2 ^2)+db2 ∗ s i n ( te ta+beta2+ ta+beta2+beta02 ) ∗ f f b 2 ; 17 18 dh gc2=dc2 gc2=dc2 ∗ s i n ( t e t a+gam a+gamma2 ma2+ +gamma02 gamma02 ) ∗ f c 2 ;
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
37
19 d dhgc2=dc2 dhgc2=dc2 ∗ co s ( t e t a+gam a+gamma2 ma2+ +gamma02 gamma02 ) ∗ ( fc 2 ^2)+ ^2)+dc2 ∗ s i n ( t e t a+gam a+gamma2 ma2+ +gamma02 gamma02 ) ∗ ffc2 ; 20 21 DLdin1= DLdin1=−a ∗ co s ( t e t a +a +a l f a )− b s t ∗ co s ( t e t a +b +b e t a s t +beta )+a ∗ co s ( t e t a + a l f a s )+ )+b s t ∗ co s ( te ta+bet as t+bet t+bet as ) ; 22 dDL1=a dDL1=a ∗ s i n ( t e t a +a +a l f a ) +b +b s t ∗ s i n ( t e t a +b +b e t a s t +beta ) ∗ f b ; 23 ddDL1=a ddDL1=a ∗ co s ( t e t a +a +a l f a ) +b +b s t ∗ co s ( t e t a +b +b e t a s t +beta ) ∗ ( fb ^2)+bst ^2)+bst ∗ s i n ( t e t a +b +b e t a s t +beta ) ∗ f f b ; 24 % d e r iv iv a t e r i s p e t t o ad a l f a 25 Dgamma2= Dgamma2=gamma2 gamma 2−gammas gammas2 2; 26 dDgam dDgamma2 ma2=f =f c 2 ; 27 ddDg ddDgam amma ma2 2=f f c 2 ;
5.6 5.6 1 2 3 4 5 6 7 8 9 10
function dDL2 = smorz ator e ( al fa , beta , fb ) e t a0 a0 0 t e t a x 1 y 1 global a b 0 0 b et
x2=a x2=a ∗ s i n ( tet a+al fa )+b00 )+b00∗ s i n ( t e t a +beta+beta0 0 ) ; y2=−a ∗ co s ( t e t a +a +a l f a )−b0 0 ∗ co s ( t e t a +beta+beta0 0 ) ; dx2=a dx2=a ∗ co s ( tet a+al fa )+b00 )+b00∗ co s ( t e t a +beta+beta0 0 ) ∗ f b ; dy2=a dy2=a ∗ s i n ( tet a+al fa )+b00 )+b00∗ s i n ( t e t a +beta+beta0 0 ) ∗ f b ; dDL2=0.5 ∗ ( ( x 2−x1)^2+(y2−y 1 ) ^ 2 ) ^ ( − 0.5) ∗ ( 2 ∗ ( x2 −x1 ) ∗ dx2+2 ∗ ( y2 −y1 ) ∗ dy2) ;
5.7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
smor smorza zato tore re.m .m
manove manovellism llismoNON oNONlinea lineare.m re.m
function d y = m a n o v e l l i s m o N O N l i n e a r e ( t , y ) global global global global global
g ma mb mc mb2 mc2 mc2 k1 k 2 r 1 r 2 z0 wz psiz DL01 Dgamma20 Dgamma20
alfap = y(1) ; alfa = y(2) ; sum ( z 0 . ∗ s i n ( wz ∗ t+psi z ) ) ; z = sum sum ( z 0 . ∗ wz . ∗ co s ( wz ∗ t+ps z p = sum t+ps iz ) ) ;
[ beta , fb , ffb , gamma, f c , f f c , b e t a 2 , f b 2 , f f b 2 , g am amma2 , f c 2 , f f c 2 ] = c i n e m a t i c a ( a l f a );
16 17 [ M, M, dM dM ] = e n e r g i a c i n e t i c a ( a l f a , beta , fb , ffb , fc , ffc , beta2 , fb2 , ffb 2 , fc2 , f fc 2 ) ; 18
38
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
19 [ dhga , ddhga , dhgb , ddhgb , dhgc , ddhgc , dhgb2 , ddhgb2 , dhgc2 , ddhgc 2 , DLdin1 , dDL1, dDL1, ddD ddDL1 ,Dga ,Dgamma2, dDgamma2, ddDgamma2] = e ne rg ia po t en zi al e ( al fa , beta , beta2 , gamma,gam ,gamma2, fb , ff b , fb2 , ffb 2 , fc , ff c , fc2 , ff c 2 ) ; 20 21 dDL dDL2=smorz at ore ( al fa , beta , fb ) ; 22 23 dV= g ∗ (ma∗ dhga+mb∗ dhgb+mc ∗ dhgc+mb2∗ dhgb2+mc2∗ dhgc2)+k1 ∗ (DLdin1−z+DL01) ∗ dDL1 +k2 ∗ (Dgamma2+Dgamma20) ∗dDgamm dDgamma2 a2 ; 24 25 dD = r1 ∗dDL1^2∗ a l f a p −r 1 ∗ dDL1∗ zp+r2 ∗ dDL2^2∗ a l f a p ; 26 27 dy=[ −1 ∗ (dM/2 ∗ a l fa p^2+ p^2+dD+dV) /M; a l fa p ] ;
5.8 1 2 3 4 5 6 7 8 9 10 11 12 13 14
mano manovellism ellismoli olinea neare. re.m m
re ( t , y ) function d y = m a n o v e l l i s m o l i n e a re g DL01 DL01 Dgam Dgamma ma20 20 ma mb mb mc mc m mb b2 mc2 k 1 k 2 r 1 r 2 z0 wz psiz a l f a s b et e t a s b et e t as a s 2 ga g ammas ga gammas2 f b s f f b s f b 2s 2s f f b 2 s f c s f f c s ffc2s
global global global global
fc2s
alfap = y(1) ; alfa = y(2) ; sum ( z 0 . ∗ s i n ( wz ∗ t+psi z ) ) ; z = sum sum ( z 0 . ∗ wz . ∗ co s ( wz ∗ t+ps z p = sum t+ps iz ) ) ;
[ M, M, dM dM ] = e n e r g i a c i n e t i c a ( a l f a s , b e ta ta s , f b s , f f b s , f c s , f f c s , b e t a s 2 , f b 2 s , f f b 2 s , fc2s , ff c2 s ) ;
15 16 [ dhga , ddhga , dhgb , ddhgb , dhgc , ddhgc , dhgb2 , ddhgb2 , dhgc2 , ddhgc 2 , DLdin1 , dDL1, dDL1, ddD ddDL1 ,Dga ,Dgamma2, dDgamma2, ddDgamma2] = e ne rg ia po t en zi al e ( al fa s , bet as , betas2 , gam gammas ,gamm ,gammas2 as2 , fbs , ffb s , fb2s , ffb 2s , fc s , ff cs , fc2 s , ff c 2 s ) ; 17 18 dDL2 = smorz ator e ( al fa s , betas , fbs ) ; 19 2 0 d V = ( k 1 ∗dDL1^2+k1 ∗ DL01 ∗ddDL1+k2 ∗dDgamma2^2+k2 ∗Dgamma20∗ddDgamma2+ma∗ g ∗ ddhga+mb∗ g ∗ ddhgb+mc∗ g ∗ ddhgc+mb2 ∗ g ∗ ddhgb2+mc2∗ g ∗ ddhgc2) ∗ ( a l f a − a l f a s )−k1 ∗ dDL1 ∗ z ; 21 dD dD = ( r1 ∗dDL1^2+r2 ∗ dDL2^2) ∗ a l f a p −r 1 ∗ dDL1 ∗ zp ; 22 2 3 d y = [ − (dD+dV (dD+dV)) /M; 24 alfap ] ;
5.9
rk4.m
CAPITOL CAPITOLO O 5. SCRIPTS SCRIPTS DEL PROGRAMM PROGRAMMA A MATLAB MATLAB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
function [ t ,x ,xp]=rk4 ( sta tof un , tf , dt , x0) % i n t e g r a z i o n e c on on i l m et et od od o d i Run ge− K u tt tt a d e l 4 . o rd r d in i n e a p as a s so so f i s s o
t =[0: dt : tf ] ; np= np=length ( t ) ;
kut=[dt/6;2∗ d t / 6 ; 2∗ dt /6; dt /6] ; x ( : , 1 ) =x =x 0 ; x p ( : , 1 ) = f e v a l ( sta tof un ,0 , x0) ; a t t= t = w a it i t b a r ( 0 , ’ A t t e nd nd e re re , r i s o l u z i o n e i n c o r s o . . . ’ ) ; f o r i =2:np
t e mp x ( : t e mp x ( : t e mp x ( : t e mp x ( :
, 1 )= )=f e v a l ( stato fun , 2 )= )=f e v a l ( stato fun , 3 )= )=f e v a l ( stato fun , 4 )= )=f e v a l ( stato fun
, t ( i −1) ,x ( : , i −1)) ; , t ( i −1)+dt 1)+dt /2 ,x (: , i −1)+dt/2 ∗ tempx tempx ( : , 1 ) ) ; , t ( i −1)+dt 1)+dt /2 ,x (: , i −1)+dt/2 ∗ tempx tempx ( : , 2 ) ) ; , t ( i −1)+ 1)+dt , x (: , i −1)+dt ∗ tempx tempx ( : , 3 ) ) ;
x ( : , i )=x ( : , i −1)+tempx ∗ kut ; x p ( : , i )=f e v a l ( stat of un , t ( i ) ,x ( : , i ) ) ; w ai a i t ba b a r ( i / np ) end c l o s e ( att ) ;
39
