´ Uvod do simulace dynamiky sypk´ych hmot Luk´ aˇs Posp´ıˇsil Katedra aplikovan´ e matematiky, ˇ - Technick´ VSB a Univerzita Ostrava
Semin´ aˇr aplikovan´e matematiky, 1. 4. 2014
Osnova
formulace probl´emu transformace lok´ aln´ıch souˇradnic Eulerovy u ´hly jednotkov´ e kvaterninony rotace od pozice k rychlosti
kontaktn´ı u ´lohy bez tˇren´ı, reformulace na QP kr´ atce o Coulombovsk´ em tˇren´ı
Pˇredn´ aˇska obsahuje 4 videa, jeden pˇr´ıklad a jeden d˚ ukaz.
Formulace probl´emu
2
1
5
0
−1
4
−2 −2
0
2 −2
−1
0
1
2
3 2
2
1
1
0
0
−1
−2 −2
−1 0
2 −2
−1
0
1
2
−2 −2
2
0
1
0
−1
2
−2 −2
4 0
2 −2
−1
0
1
2
2
4
0 −2
Formulace probl´emu
Vztah mezi lok´ aln´ım a glob´ aln´ım syst´emem T := [Tx , Ty , Tz ]T - poloha tˇeˇziˇstˇe v glob´ aln´ım souˇradnicov´em syst´emu, φ, θ, ψ - Eulerovy u ´hly rotace. rglob = T + Rrlok rlok ∈ R3 souˇradnice bodu v lok´ aln´ım syst´emu, rglob ∈ R3 souˇradnice bodu v glob´ aln´ım syst´emu. R ∈ R3,3 matice rotace (z φ, θ, ψ).
Eulerovy u ´hly rotace
2 1 0 −1
rlok ∈ R3
−2 −2 0 2 −2
−1
0
1
2
Eulerovy u ´hly rotace
2 1 0
cos φ R1 = sin φ 0
−1
− sin φ cos φ 0
−2
R1 rlok −2 0 2 −2
−1
0
1
2
0 0 1
Eulerovy u ´hly rotace
2 1 0
1 R2 = 0 0
−1
0 cos θ sin θ
−2
R2 R1 rlok −2 0 2 −2
−1
0
1
2
0 − sin θ cos θ
Eulerovy u ´hly rotace
2 1 0
cos ψ R1 = sin ψ 0
−1
− sin ψ cos ψ 0
−2
R3 R2 R1 rlok −2 0 2 −2
−1
0
1
2
0 0 1
Eulerovy u ´hly rotace
rglob = T + R3 R2 R1 rlok | {z } =:R
cos ψ cos φ − cos θ sin φ sin ψ R = cos ψ cos φ + cos θ cos φ sin ψ sin θ sin ψ
− sin ψ cos φ − cos θ sin φ cos ψ − sin ψ sin φ + cos θ cos φ cos ψ sin θ cos φ
sin θ sin φ − sin θ sin φ cos θ
´ Uloha se mˇen´ı v ˇcase (Tx (t), Ty (t), Tz (t), φ(t), θ(t), ψ(t)), tj. okamˇzit´ a poloha bodu v ˇcase t ∈ R+ rglob (t) = T (t) + R(t)rlok
Rychlost
rglob (t) = T (t) + R(t)rlok Okamˇzit´ a rychlost Okamˇzitou rychlost z´ısk´ ame derivac´ı dle t ˙ lok r˙glob = T˙ + Rr
R˙ =?
Rychlost
R je ortogon´ aln´ı, tj. RT R = I , derivac´ı ˙ T + R R˙ T = 0. RR ´ Upravou ˙ T RR ˙ RR T ˙ T je antisymetrick´ ⇒ matice RR a. Je to matice vektorov´eho souˇcinu.
= =
−R R˙ T ˙ T) −(RR
Matice vektorov´eho souˇcinu
Mˇejme u = [u1 , u2 , u3 ]T ∈ R3 , v = [v1 , v2 , v3 ]T ∈ R3 . Pak i u×v = u1 v1
j u2 v2
k u3 v3
u2 v3 − v2 u3 = (u2 v3 −v2 u3 )i+(u3 v1 −u1 v3 )j+(u1 v2 −u2 v1 )k = u3 v1 − u1 v3 u1 v2 − u2 v1
Lze vˇsak tak´e maticovˇe
u2 v3 − v2 u3 u×v =u ˜v = u3 v1 − u1 v3 u1 v2 − u2 v1 kde u ˜ ∈ R3,3 je matice vektorov´eho souˇcinu 0 −u3 0 u ˜ := u3 −u2 u1 Tato matice je antisymetrick´ a, tj. u ˜ = −˜ uT
u2 −u1 . 0
Rychlost
Oznaˇcme ˙ ˜˙ := R T R. ω Pak u ´pravou ˜˙ ω ˜˙ Rω ˜˙ Rω Dosazen´ım r˙glob r˙glob r˙glob
= = =
= = =
R T R˙ RR T R˙ R˙
˙ lok T˙ + Rr ˙ ˜˙ lok T + R ωr ˙ T + R(ω˙ × rlok )
˙ θ, ˙ ψ] ˙ T ∈ R3 je u kde ω˙ = [φ, ´hlov´ a rychlost.
Motivace rglob r˙glob
= =
T + Rrlok T˙ + R(ω˙ × rlok )
Vyˇc´ıslit R znamen´ a vyˇc´ıslit goniometrick´e funkce. Jednotkov´y kvaternion rotace R3 3 [φ, θ, ψ] := ω
↔
e := [e0 , e1 , e2 , e3 ] ∈ R4 , kek = 1
D´ıky t´eto reprezentaci nebude nutno vyˇc´ıslovat goniometrick´e funkce.
Kvaternion rotace - motivace
[φ]
↔
[e0 , e1 ]
Pol´ arn´ı souˇradnice e1
e0
e0 e1
= =
cos φ sin φ
Kvaternion rotace - motivace
[φ, θ]
↔
[e0 , e1 , e2 ]
Sf´erick´e souˇradnice e2
e0
e1
e0 e1 e2
= = =
sin φ cos θ sin φ sin θ cos φ
Kvaternion rotace - motivace
[φ, θ, ψ]
↔
[e0 , e1 , e2 , e3 ]
Kvaternion rotace e0
=
e1
=
e2
=
e3
=
+ θ) cos 12 ψ cos 21 (φ − θ) sin 21 ψ sin 12 (φ − θ) sin 12 ψ sin 12 (φ + θ) cos 21 ψ cos
1 (φ 2
Kvateriony m´ısto Eulerov´ych u ´hl˚ u
Zaved’me matice
−e1 e0 −e3 e2 e3 e0 −e1 = [−e123 , e˜123 + e0 I ] E := −e2 −e3 −e2 e1 e0 −e1 e0 e3 −e2 e0 e1 = [−e123 , −˜ e123 + e0 I ] G := −e2 −e3 −e3 e2 −e1 e0 kde T T e123 = [e1 , e2 , e3 ]T , e = [e0 , e1 , e2 , e3 ]T = [e0 , e123 ]
Kvateriony m´ısto Eulerov´ych u ´hl˚ u
Plat´ı
−e1 Ee = −e2 −e3
e0 e3 −e2
−e3 e0 e1
e0 −e3 e2
e3 e0 −e1
−e1 Ge = −e2 −e3
e0 e2 e1 −e1 e2 e0 e3 e0 −e2 e1 e1 e2 e0 e3
=0 =0
Kvateriony m´ısto Eulerov´ych u ´hl˚ u
Plat´ı
EE T
−e1 = −e2 −e3
e0 e3 −e2
−e3 e0 e1
−e1 e2 e0 −e1 −e3 e0 e2
(Matice E m´ a ortonorm´ aln´ı ˇr´ adky) −e1 −e2 −e3 −e1 e0 e3 −e2 −e2 ETE = −e3 e0 e1 −e3 e2 −e1 e0
e0 e3 −e2
−e2 e3 e0 −e1
−e3 e0 e1
∈ R3
e2 −e1 = I − ee T e0
Obdobnˇe GG T = I ,
−e3 −e2 =I e1 e0
G T G = I − ee T
∈ R4
Kvateriony m´ısto Eulerov´ych u ´hl˚ u
Plat´ı (!) −e1 −e2 −e3 −e1 e0 −e3 e2 e0 −e3 e2 −e2 e3 e0 −e1 e3 e0 −e1 −e3 −e2 e1 e0 −e2 e1 e0 2 e0 − e12 − e22 − e32 2e1 e2 − 2e0 e3 2e1 e3 + 2e0 e2 2e1 e2 + 2e0 e3 e02 − e12 + e22 − e32 2e2 e3 − 2e1 e0 2e1 e3 − 2e0 e2 2e2 e3 + 2e1 e0 e02 − e12 − e22 − e32 ··· = R
EG T
=
= =
Obdobnˇe se d´ a uk´ azat E G˙ T = E˙ G T Pak R˙ = E˙ G T + E G˙ T = 2E G˙ T
Kvateriony m´ısto Eulerov´ych u ´hl˚ u
Pro derivaci kvaternin˚ u plat´ı (d´ a se uk´ azat z vlastnost´ı vektorov´eho souˇcinu) e˙ =
1 T E ω˙ 2
Proto d´ ale budeme popisovat polohu lok´ aln´ıch souˇranic v glob´ aln´ım syst´emu pomoc´ı zobecnˇen´eho vektoru pozice a zobecnˇen´eho vektoru rychlost´ı Tx T˙ x Ty T˙ y Tz ˙ Tz 7 vi := ˙ ∈ R6 qi := e0 ∈ R , φ e1 ˙ θ e2 ψ˙ e3
Od rychlosti k poloze
ˇ Casov´ e sch´ema Aproximace derivace diferenc´ı q(t) ˙ := lim
h→0
q(t + h) − q(t) h
≈ q(t) ˙ =
q(t + h) − q(t) . h
ˇ Casov´ e sch´ema m´ a pak tvar q(t + h) = q(t) + hq(t) ˙ , kde h ∈ R+ je dostateˇcnˇe mal´y ˇcasov´y krok. S kvaternionem rotace pak q(t + h) = q(t) + h |
I 0
0 v (t) . 1 T E {z2 }
=:Q
Algoritmus
D´ ano q(0) ∈ R7nb , v (0) ∈ R6nb . Zvol ˇcasov´y krok h > 0 t := 0 dokud t < tmax q(t + h) := q(t) + hQv (t) spoˇcti v (t + h) t := t + h konec dokud
Prvn´ı test
pust’ animaci.
V´ypoˇcet rychlosti ´ Ukolem je tedy nal´ezt vektor aktu´ aln´ı rychlosti v (t), jehoˇz pˇr´ır˚ ustek oproti v (t − h) je z´ avisl´y na hmotnosti jednotliv´ych tˇeles, p˚ usob´ıc´ıch vnˇejˇs´ıch sil´ ach F (t, q, v ), kontaktech a omezuj´ıc´ıch podm´ınk´ ach.
M v˙ = FC + Fext
M(v (t + h) − v (t)) = h(Fext + FC ) ,
kde M je zobecnˇen´ a matice hmotnosti (diagon´ aln´ı), FC je vektor sil vyvolan´ych omezen´ım, Fext je vektor vnˇejˇs´ıch sil. ´ Upravou v (t + h) = v (t) + hM −1 (Fext + FC )
Vektor vnˇejˇs´ıch sil - gravitace
Fg = m.g ,
m = ρ.V
Fext = [0, 0, −Fg , 0, 0, 0]T D´ ano q(0) ∈ R7nb , v (0) ∈ R6nb . Zvol ˇcasov´y krok h > 0 t := 0 dokud t < tmax q(t + h) := q(t) + hQv (t) v (t + h) := v (t) + M −1 (Fext + FC ) t := t + h konec dokud
Druh´y test
pust’ dalˇs´ı skvˇelou animaci.
Kontakt a omezen´ı
Uvaˇzujme kontakt mezi tˇelesy TA a TB .
TA FA
C
FB
TB
Kontakt a omezen´ı Oznaˇcme nA (C ) jako vnˇejˇs´ı jednotkovou norm´ alu tˇelesa TA v bodˇe kontaktu C v glob´ aln´ım souˇradnicov´em syst´emu. Zmˇena pozice tˇeˇziˇstˇe Tˇeleso TA p˚ usob´ı na tˇeleso TB silou FB := γ¯ nA (C ) , kde γ¯ ≥ 0 je nezn´ am´ a velikost p˚ usob´ıc´ı s´ıly. Tato s´ıla ovlyvn´ı zmˇenu polohy tˇelesa TB (sloˇzky rychlosti tˇelesa odpov´ıdaj´ıc´ı poloze tˇeˇziˇstˇe).
Zmˇena rotace tˇelesa Zmˇenu rotace tˇelesa TB ovlyvˇ nuje krout´ıc´ı moment MB := C B × FB , kde C B jsou souˇradnice bodu C v lok´ aln´ım souˇradnicov´em syst´emu tˇelesa TB .
Kontakt a omezen´ı
´ Upravou e B nA (C ) . MB = C B × FB = γ¯ (C B × nA (C )) = γ¯ C Naopak na tˇeleso TA p˚ usob´ı opaˇcn´ a s´ıla FA := −¯ γ nA (C ) , a krout´ıc´ı moment e A nA (C ) . MA := C A × FA = −¯ γC
Kontakt a omezen´ı
Pˇripomeˇ nme −¯ γ nA (C ) FA e A nA (C ) MA −¯ γ C = FC := FB γ¯ nA (C ) e B nA (C ) MB γ¯ C
kde −nA (C ) e A nA (C ) −C . D := nA (C ) B e C nA (C )
= γ¯ D ,
Kontakt a omezen´ı
Pokud se n´ am podaˇr´ı urˇcit velikost p˚ usob´ıc´ı s´ıly γ˜ , pak v´yslednou zmˇenu rychlosti lze spoˇc´ıtat z rovnice M(v (t + h) − v (t)) = hFext + Dγ , kde γ := h¯ γ.
Kontakt a omezen´ı
Pˇr´ıklad Uvaˇzujme syst´em tˇr´ı tˇeles T(1) , T(2) , T(3) .
T(1)
C(1,2)
C(1,3) C(2,3) T(3)
T(2)
Kontakt a omezen´ı
T T T T vektor zobecnˇen´e polohy q = [q(1) , q(2) , q(3) ] ∈ R3∗7 T T T T vektor zobecnˇen´e rychlosti v = [v(1) , v(2) , v(3) ] ∈ R3∗6 .
Hled´ ame velikosti ˇsesti nezn´ am´ych sil: s´ıla, kterou p˚ usob´ı tˇeleso T(1) na tˇeleso T(2) , s´ıla, kterou p˚ usob´ı tˇeleso T(2) na tˇeleso T(1) , s´ıla, kterou p˚ usob´ı tˇeleso T(2) na tˇeleso T(3) , s´ıla, kterou p˚ usob´ı tˇeleso T(3) na tˇeleso T(2) , s´ıla, kterou p˚ usob´ı tˇeleso T(1) na tˇeleso T(3) , s´ıla, kterou p˚ usob´ı tˇeleso T(3) na tˇeleso T(1) , ˇcili hled´ ame γ = [γ(1,2) , γ(2,1) , γ(2,3) , γ(3,2) , γ(1,3) , γ(3,1) ]T ∈ R6 .
Kontakt a omezen´ı
Matice D ∈ R3∗6,6 m´ a blokovou strukturu D =
−n(1) (C(1,2) ) e (1) n (C −C ) (1,2) (1) (1,2) n(1) (C(1,2) ) (2) e C n (C ) (1,2) (1) (1,2) 0
n(2) (C(1,2) ) e (1) n (C C ) (1,2) (2) (1,2) −n(2) (C(1,2) ) (2) e −C n (C ) (1,2) (2) (1,2) 0
0
0
0 0 n(3) (C(2,3) ) e (2) n (C C ) (2,3) (3) (2,3) −n(3) (C(2,3) ) e (3) n (C −C ) (2,3) (3) (2,3)
0 0 −n(2) (C(2,3) ) e (2) n (C −C ) (2,3) (2) (2,3) n(2) (C(2,3) ) e (3) n (C C ) (2,3) (2) (2,3)
−n(1) (C(1,3) ) e (1) n (C −C ) (1,3) (1) (1,3) 0
n(3) (C(1,3) ) e (1) n (C C ) (1,3) (3) (1,3) 0
0
0
n(1) (C(1,3) ) e (3) n (C C ) (1,3) (1) (1,3)
−n(3) (C(1,3) ) e (3) n (C −C ) (1,3) (3) (1,3)
Kontakt a omezen´ı Podm´ınku nepronik´ an´ı tˇeles TA a TB popisuje tzv. gap funkce Φ : R7+7 → R ≈ ”vzd´ alenost tˇeles”(nejmenˇs´ı vzd´ alenost bod˚ u jednoho a druh´eho tˇelesa). Plat´ı Φ([qA , qB ]) = 0 pokud jsou tˇelesa v kontaktu, Φ([qA , qB ]) > 0 pokud tˇelesa nejsou v kontaktu, Φ([qA , qB ]) < 0 pokud tˇelesa ”pronikla do sebe”. Pro norm´ alov´e napˇet´ı v bodˇe kontaktu plat´ı pokud jsou tˇelesa v kontaktu, norm´ alov´e napˇet´ı existuje, pokud tˇelesa nejsou v kontaktu, norm´ alov´e napˇet´ı neexistuje. Z´ısk´ ame podm´ınku komplementarity a tedy uˇzit´ım gap funkce lze definovat podm´ınky pro v´ypoˇcet velikosti p˚ usob´ıc´ıch sil 0 ≤ Φ(q) ⊥ γ ≥ 0 . ⇔ (numericky stabilnˇejˇs´ı) 0≤
1 Φ(q) + D T v (t + h) ⊥ γ ≥ 0 . h
Linear complementarity problem ´ Ukolem je tedy ˇreˇsit n´ asleduj´ıc´ı probl´em. Linear complementarity problem (LCP)
M(v (t + h) − v (t)) = hFext + Dγ 1 Φ(q) + D T v (t + h) ≥ 0 h 1 Φ(q) + D T v (t + h) ⊥ γ h γ≥0
Ne.
(1a) (1b) (1c) (1d)
QP s line´arn´ım omezen´ım
Theorem ˇ sen´ı optimalizaˇcn´ı u Reˇ ´lohy 1 min γ T Nγ + r T γ , γ≥0 2
(2)
kde N := D T M −1 D 1 r := Φ + D T M −1 k h k := Mv (t) + h.Fext je ekvivalentn´ı s ˇreˇsen´ım p˚ uvodn´ı u ´lohy (LCP).
(3a) (3b) (3c)
QP s line´arn´ım omezen´ım D˚ ukaz Lagrangeova funkce m´ a tvar L(γ, λ) =
1 T γ Nγ + r T γ − λT γ 2
a odpov´ıdaj´ıc´ı KKT podm´ınky Nγ + r − λ = 0
(4a)
γ≥0
(4b)
λ≥0
(4c)
γ λ=0
(4d)
T
Podm´ınky (4b) a (1d) jsou stejn´e. Z podm´ınky (4a) pˇr´ımo plyne λ = Nγ + r a spolu s (3a),(3b),(3c) a dosazen´ım do (4d) z´ısk´ ame 0
= = = =
γ T λ = γ T (Nγ + r ) γ T (D T M −1 Dγ + h1 Φ + D T M −1 k) γ T (D T M −1 Dγ + h1 Φ + D T M −1 (Mv (t) + h.Fext )) γ T ( h1 Φ + D T (v (t) + M −1 Dγ + h.M −1 Fext ))
(5)
QP s line´arn´ım omezen´ım
N´ aslednˇe uˇzit´ım (1a) z´ısk´ ame 1 0 = γ T ( Φ + D T (v (t) + M −1 Dγ + h.M −1 Fext )) h coˇz je podm´ınka (1c). Nakonec dosazen´ım podm´ınky (5) do (4c) lze z´ıskat 0 ≤ λ = Nγ + r = coˇz je podm´ınka (1b).
1 Φ + D T v (t + h) , h
Algoritmus
D´ ano q(0) ∈ R7nb , v (0) ∈ R6nb . Zvol ˇcasov´y krok h > 0 t := 0 dokud t < tmax q(t + h) := q(t) + hQv (t) Nalezeni kontakty. Pokud neexistuj´ı kontakty, poloˇz γ := 0 Pokud existuj´ı kontakty, sestav N, r a vyˇreˇs optimalizaˇcn´ı probl´em v (t + h) := v (t) + M −1 (hFext + Dγ) t := t + h konec dokud
Tˇret´ı test
Nefl´ akej se a pust’ dalˇs´ı animaci.
Tˇren´ı
Uvaˇzujme (opˇet) kontakt mezi tˇelesy TA a TB . Oznaˇcme n ∈ R3 jednotkov´y norm´ alov´y vektor kontaktu v glob´ aln´ıch souˇradnic´ıch, u, w ∈ R3 jednotkov´e smˇerov´e vektory kontaktn´ı roviny v glob´ aln´ıch souˇradnic´ıch. Oˇcividnˇe {n, u, w } je mnoˇzina ortonorm´ aln´ıch vektor˚ u (a tvoˇr´ı b´ azi teˇcn´eho prostoru). Tˇrec´ı s´ılu p˚ usob´ıc´ı v bodˇe kontaktu na tˇeleso A pak m˚ uˇzeme vyj´ adˇrit jako F = Fn + FT = γ n n + γ u u + γ w w , kde Fn = γn n ∈ R3 je norm´ alov´ a sloˇzka tˇrec´ı s´ıly, FT = γu u + γw w ∈ R3 je teˇcn´ a sloˇzka tˇrec´ı s´ıly, γn > 0 je velikost norm´ alov´e sloˇzky tˇrec´ı s´ıly F , γu , γw ∈ R jsou velikosti teˇcn´ych sloˇzek tˇrec´ı s´ıly F .
Tˇren´ı
Vztah mezi jednotliv´ymi sloˇzkami tˇrec´ı s´ıly popisuje napˇr´ıklad Coulomb˚ uv model tˇren´ı. Coulomb˚ uv model tˇren´ı
γn ≥ 0,
Φ(q) ≥ 0, Φ(q)γn = 0 , p γu2 + γw2 ≤ µγn , p kvT k µγn − γu2 + γw2 = 0 ,
(6b)
hFT , vT i = −kFT k.kvT k .
(6d)
Objasnˇeme v´yznam jednotliv´ych rovnic.
(6a)
(6c)
Tˇren´ı
γn ≥ 0,
Φ(q) ≥ 0,
Φ(q)γn = 0
podm´ınka komplementarity gap funkce Φ a velikosti norm´ alov´e sloˇzky γn norm´ alov´ a sloˇzka tˇrec´ı s´ıly je vˇzdy nez´ aporn´ a (pokud uvaˇzujeme kontakt, pak tˇeleso A nem˚ uˇze p˚ usobit na tˇeleso B silou, kter´ a smˇeˇruje ”od tˇelesa” B), hodnota gap funkce je vˇzdy nez´ aporn´ a (vzd´ alenost tˇeles A a B je nez´ aporn´ a), pokud je norm´ alov´ a sloˇzka nenulov´ a, pak je nulov´ a vzd´ alenost tˇeles, pokud je vzd´ alenost tˇeles nenulov´ a, pak je norm´ alov´ a sloˇzka tˇrec´ı s´ıly nulov´ a, m˚ uˇze nastat situace, ˇze vzd´ alenost tˇeles je nulov´ a i norm´ alov´ a sloˇzka tˇrec´ı s´ıly je nulov´ a (tˇelesa jsou ”v klidu” v kontaktu).
Tˇren´ı
p γu2 + γw2 ≤ µγn rovnice vztahu mezi velikost´ı norm´ alov´e sloˇzky tˇrec´ı s´ıly a teˇcn´e sloˇzky tˇrec´ı s´ıly, tj. kFT k ≤ µkFN k , kde µ ∈ h0, 1i je koeficient tˇren´ı takov´y, ˇze pokud µ = 1, pak velikost teˇcn´e sloˇzky tˇrec´ı s´ıly m˚ uˇze b´yt maxim´ alnˇe rovna velikosti norm´ alov´e s´ıly, pokud µ = 0, pak velikost teˇcn´e sloˇzky tˇrec´ı s´ıly je vˇzdy nulov´ a, ˇcili s´ıla p˚ usob´ı pouze ve smˇeru norm´ aly kontaktu a tedy nedoch´ az´ı ke tˇren´ı.
Tˇren´ı
p kvT k µγn − γu2 + γw2 = 0 podm´ınka komplementarity mezi velikost´ı rychlosti pohybu tˇelesa v teˇcn´em smˇeru a dosaˇzen´ı maxim´ aln´ı moˇzn´e velikosti tˇrec´ı s´ıly vzhledem k omezen´ı pˇredchoz´ı rovnic´ı pokud nedojde k maxim´ aln´ı moˇzn´e velikosti tˇrec´ı s´ıly vzhledem k µ, tj. kFT k < µkFN k, pak velikost teˇcn´e rychlosti tˇelesa je nulov´ a. pokud naopak dojde k maxim´ aln´ı moˇzn´e velikosti teˇcn´e s´ıly, pak doch´ az´ı k pˇresmyku a velikost teˇcn´e rychlosti m˚ uˇze b´yt jak´ akoliv.
Tˇren´ı
hFT , vT i = −kFT k.kvT k podm´ınka smˇeru tˇrec´ı s´ıly a teˇcn´e rychlosti, konkr´etnˇe pro u ´hel mezi tˇemito vektory plat´ı hFT , vT i cos ϕ = = −1, tj. ϕ = −π . kFT k.kvT k ˇ teˇcn´ Cili a sloˇzka tˇrec´ı s´ıly m´ a opaˇcn´y smˇer jako vektor rychlosti pˇresmyku (kluzn´ a rychlost).
Tˇren´ı
Obdobnˇe jako v u ´loze bez tˇren´ı FA FB
:= :=
−γn n − γu u − γw w γn n + γu u + γw w
= =
kde γ := [γn , γu , γw ]T ∈ R3 S := [n, u, w ] ∈ R3,3
−Sγ , Sγ ,
Tˇren´ı
MA MB
:= :=
e A Sγ , C A × F¯A = −C A × (Sγ) = −C B B B e ¯ C × FB = C × (Sγ) = C Sγ ,
Pak −Sγ FA e A Sγ MA −C FC := FB = Sγ e B Sγ MB C
kde −S e AS −C . D := S eBS C
= Dγ ,
Quadratic Cone Complementarity problem
Theorem ˇ sen´ı optimalizaˇcn´ı u Reˇ ´lohy 1 min γ T Nγ + r T γ , 2
γ≥Ω
kde N := D T M −1 D 1 r := [ Φ, 0, 0]T + D T M −1 k h k := Mv (t) + h.Fext
(7a) (7b) (7c) (7d)
Ω := Ω1 × · · · × Ωnc
(7e)
p Ωi := {[x, y , z] ∈ R : y 2 + z 2 ≤ µi x} T
3
je ekvivalentn´ı s ˇreˇsen´ım p˚ uvodn´ı u ´lohy.
(7f)
Quadratic Cone Complementarity problem
D˚ ukaz Pˇr´ıˇstˇe.
Otevˇren´e probl´emy
h =? efektivn´ı implementace (GPU) detekce kontakt˚ u domain decomposition (?) QP ˇreˇsiˇc
Dˇ ekuji za pozornost Dost´ al, Z.: Optimal Quadratic Programming Algorithms, with Applications to Variational Inequalities, SOIA, first ed., vol. 23, Springer, US, New York, 2009. Heyn T.: On the Modeling, Simulation, and Visualization of Many-Body Dynamics Problems with Friction and Contact. Ph.D. Thesis, 2013. online Heyn T., Anitescu M., Tasora A., Negrut D.: Using Krylov Subspace and Spectral Methods for Solving Complementarity Problems in Many-Body Contact Dynamics Simulation. Int. J. Numer. Meth. Engng, 2012. Shabana, A.: Computational Dynamics, John Wiley & Sons, Inc., New York, Third Edition, 2010. Haug, E. J.: Computer-Aided Kinematics and Dynamics of Mechanical Systems, Volume-I. Prentice-Hall, Englewood Cliffs, New Jersey, 1989.
