Energy frictional dissipating algorithm for rigid and ellastic body’s contact problems

135 ENERGY FRICTIONAL DISSIPATING ALGORITHM FOR RIGID AND ELLASTIC BODY’S CONTACT PROBLEMS R. B a o∗and J.L. P´e ez–Apa icio† ∗Depa men S uc u al Mechanics and Hyd aulic Enginee ing Uni e si y o G anada Campues de Fuen enue a, 18071 G anada, Spain e-mail: b a [email p o ec ed] †Depa men o Con inuum Mechanics and Theo y o S uc u es Uni e sidad Poli ´ecnica de Valencia 46022 Valencia, Spain e-mail: jopeap@up ne .up .es Key wo ds: Con ac , Time in eg a ion scheme, Consis en Dissipa ion, F ic ion Abs ac . An Ene gy F ic ional Dissipa ing Algo i hm (EFDA) o ime in eg a ion o Coulomb ic ional impac –con ac p oblems is p esen ed. Using he Penal y Me hod, and in he con ex o a conse ing amewo k, linea and angula momen a a e conse ed and ene gy is consis en ly dissipa ed. Published o mula ions we e s able, o cing he ene gy dissipa ion o be mono onic in o de o p e en uns able ene gy g ow h. The sho coming o many was ha hey we e no able o ep oduce he eal kinema ics and dissipa ion o physical p ocesses, p o ided by analy ical o mula ions and expe imen s. EFDA o mula es a conse ing amewo k based on a physical ene gy dissipa ion es ima o . This amewo k uses an enhanced Penal y con ac model based on a sp ing and a dashpo , en o cing physical ic ional ene gy dissipa ion, con olling gap ib a ions and modi ying he eloci ies and con ac o ces du ing each ime s ep. The esul is ha he dissipa ed ene gy, kinema ics and con ac o ces a e consis en wi h he expec ed physical beha io . 1 INTRODUCTION The nume ically accu a e analysis o ic ional dynamic con ac p oblems has been a challenge o he las 30 yea s. Complex p oblems do no ha e analy ical solu ion and due o hei high nonlinea i y, non–smoo h unila e al es ic ion and he p esence ic ion, hey a e ha d o model. The e o e, nume ical ime–s epping schemes a e de eloped o emula e he conse a i e p ope ies o he co esponding con inuous p oblem. P e ious au ho s ha e add essed ic ionless con ac p oblems, o ins ance [7] ocused on i e a i e bu no ime–s epping o mula ions, [1] and [3] o implici . These au ho s 1 XI In e na ional Con e ence on Compu a ional Plas ici y. Fundamen als and Applica ions COMPLAS XI E. Oña e, D.R.J. Owen, D. Pe ic and B. Suá ez (Eds) 136 R. B a o and J.L. P´e ez–Apa icio in ended o c ea e obus and s able algo i hms o he en o cemen o he con ac con- s ain s, while ecen o mula ions ha e ocused on ic ional o mula ion and p oposed uncondi ionally posi i e ene gy dissipa ion. Re . [3] de eloped a posi i e ene gy dissipa - ing algo i hm, s able o ic ion wi h he Penal y Me hod and showed an a i icial ene gy ans e be ween bodies–penal y sp ings in which he inal ene gy was always lowe han he ini ial o S ick and Slip cases. The e o e he beha io o he simula ion was no consis en wi h he physical con ac p oblem. Re . [1] minimized ha a i icial ene gy ans e be ween body and penal y sp ing. Fo he non–sliding si ua ion he ene gy a e con ac was equal o he ini ial and lowe du ing con ac while o he Slip case ob ained a igo ous posi i e ene gy dissipa ion. The dissipa ion in bo h e e ences was no based on a consis en conse ing amewo k: he dissipa ion, al hough dec easing mono onically, was no in acco dance o ha o he con inuous p oblem. Re . [5] de eloped a conse - a i e amewo k ha en o ced he impene abili y condi ion, elimina ed he a i icial ene gy ans e be ween body–penal y sp ing and ook in o accoun he ic ional dissipa- ion h ough an ene gy es ima o . This o mula ion used a con ac eloci y ha modi ied a p edic o –co ec o scheme, hen he con ac esponse ag eed in eloci ies bu no in o ces and posi ions, no being accu a e o pe sis en con ac . This a icle p esen s an Ene gy F ic ional Dissipa ing Algo i hm (EFDA) based on he ic ionless algo i hm o [2]. Fo Penal y con ac p oblems, he new o mula ion conse es momen a, simula es he kinema ics, con ac o ces and dissipa es ene gy consis en ly, ac- co ding o he physical p oblem in each ime s ep. The algo i hm key is a conse a i e amewo k based on upda ing con ac o ces and momen a o e e y con ac . The ame- wo k akes in o accoun dissipa ion by an ene gy es ima o based on ic ional Coulomb law, and is able o en o ce ene gy conse a ion o he S ick con ac and he igh dissi- pa ion o he Slip con ac . 2 DEFINITION OF THE PROBLEM AND GOVERNING EQUATIONS 2.1 Hamil onian desc ip ion o mo ion The Hamil onian Mechanics pe mi o ob ain he equa ions o mo ion o mul iple bod- ies ha in e ac by con ac . This subsec ion b ie ly desc ibes he Hamil onian equa ions o a con inuous p oblem. Conside a mani old Q ha desc ibes he con igu a ion o a mechanical sys em whose phase space is P=T⋆Q, he angen space o Q. This space is composed o each poin o body iby posi ions Qi(x,y, ) and linea momen a Pi(x,y, ) as unc ion o ime . The Hamil onian unc ion HQi(x,y, ),Pi(x,y, )de ines he o al ene gy o he sys em and is assumed o be sepa able in kine ic K(Pi(x,y, )) and po en ial V(Qi(x,y, )) ene gies, Eqs. 1. 2 137 R. B a o and J.L. P´e ez–Apa icio HQi(x,y, ),Pi(x,y, )= nbd  i=1 KPi(x,y, )+VQi(x,y, ) KPi(x,y, )=1 2Ωi Pi(x,y, )2 ρdΩ (1) whe e nbd is he o al numbe o bodies, Ωi he domain o body iand ρ he densi y. The kine ic ene gy is a eal unc ion K:P→Rand he po en ial V:Q→Ris an a bi a y unc ion. The mo ion is go e ned by he Hamil onian canonical equa ions, Eqs. 2. ˙ Qi(x,y, )= ∂H(x,y, ) ∂Pi=Ωi Pi(x,y, ) ρdΩ ˙ Pi(x,y, )=−∂H(x,y, ) ∂Qi=−∇VQi(x,y, ) (2) The con inuum a iables om Eqs. 2 may be disc e ized, gi ing Eqs. 3. Qi(x,y, )= nnod  A=1 NA(x, y)qA( ); Pi(x,y, )= nnod  A=1 NA(x, y)pA( ) (3) Fo he igid bodies used in he cu en pape , his disc e iza ion is based on i s o de shape unc ions NA(x, y); o he elas ic case hey may be ex ended o highe o de . Fo hese i s o de unc ions, he disc e iza ion is based on only one node, nnod = 1, usually de ined a he cen e o g a i y xi,yio each pa icle. The nodal displacemen s and linea momen a o all bodies i o ka e g ouped in he ec o s q( ), p( ). Fo each body i, he disc e iza ion Eqs. 3 applied o Eqs. 2 p oduce he sys em o equa ions: ˙ qi=M−1 ipi;˙ pi= i c+ i ex (4) whe e Miis a diagonal mass ma ix, wi h en ies: Mi=Ωiρ[NA(x, y)] NA(x, y) dΩ. Al hough con ac o ces i ca e applied in he con ac poin s, he disc e iza ion conside s an equi alen o ce applied o he nodes. The same hing can be said o he ex e nal o ces i ex . 3 NEW ALGORITHM FORMULATION AND ENERGY–MOMENTUM CONSERVATION The aim o his sec ion is he disc e iza ion in ime o Eqs. 4. The new equa ions will en o ce he impene abili y condi ion and disc e ely inhe i he conse a ion p ope ies 3 138 R. B a o and J.L. P´e ez–Apa icio h ough he conse ing amewo k o sec ion 4. The main h ee cha ac e is ic o his algo i hm a e: i) ene gy conse a ion o no mal con ac , ii) consis en dissipa ion o angen ial Slip and iii) conse a ion o angen ial S ick. 3.1 Time–disc e e o mula ion The ic ional de elopmen o EFDA is based on he Simo–Ta now’s algo i hm om [6], an ene gy–momen um conse ing ime in eg a ion scheme. This scheme is a disc e e app oxima ion o a Hamil onian sys em (Eqs. 5) in ime con igu a ion n+1/2. Conside ing he in e al [ n, n+1], he i s se o equa ions o his algo i hm ela es displacemen s qi n, qi n+1 and linea momen a pi n,pi n+1; he second, is he disc e e app oxima ion o second New on’s law a n+1/2: ˙qi=M−1 ipi→qi n+1 −qi n ∆ =M−1 ipi n+1 2 ˙pi= i cN + i cT →pi n+1 −pi n ∆ = i cN n+1 2+ i cT n+1 2 (5) whe e ∆ = n+1 − n,qi n≈qi( n), pi n≈pi( n), qi n+1 ≈qi( n+1), pi n+1 ≈pi( n+1) and pi n+1/2=(pi n+1+pi n)/2. The e ms i cN n+1/2and i cT n+1/2a e he disc e e app oxima ions o he esul ing no mal and angen ial con ac o ce ec o s. In o de o ob ain a conse ing and a igh kinema ic esponse o con ac be ween wo bodies i, k, in EFDA addi ional linea momen a pik cN n+1/2,pik cT n+1/2and con ac o ces ′ik cN n+1/2, ′ik cT n+1/2(upda ing a iables) a e added o Eqs. 5, gi ing Eqs. 6. Fo hese ou a iables, in he ollowing he subsc ip n+1/2 will be omi ed o simplici y. The ole o hese new a iables is o en o ce bodies’ ene gy conse a ion o no mal con ac and conse a ion–consis en dissipa ion o angen ial, espec i ely: qi n+1 −qi n ∆ =M−1 ipi n+1 2 + nbd  k=1 k�=ipik cN +pik cT =M−1 ipi n+1 2 +pi cN +pi cT  pi n+1 −pi n ∆ = nbd  k=1 k�=i ik cN + ′ik cN + ik cT + ′ik cT = i cN + ′i cN + i cT + ′i cT (6) The exp essions o he upda ing a iables a e now o mula ed in local con ac coo - dina es and ans o med o global by he uni no mal and angen ial ec o s Nik n+1/2, Tik n+1/2, bo h a he con ac poin . To ob ain om Eqs. 6 a conse a i e solu ion o S ick and dissipa i e o Slip, he upda ing a iables mus ul ill he disc e e conse ing equa ions de ined in sec ion 4. These a iables a e de ined o bo h con ac di ec ions as: 4 139 R. B a o and J.L. P´e ez–Apa icio NORMAL pik cN =ψik 2N 2Nik n+1 2 Nik n+1 2 (pi n+1 −pi n); ′ik cN =ψik 1N 2Nik n+1 2KN(gik Nn+1 −gik Nn) TANG. STICK pik cT =ψik 2T 2Tik n+1 2 Tik n+1 2 (pi n+1 −pi n); ′ik cT =ψik 1T 2Tik n+1 2 KT(gik Tn+1 −gik Tn) TANG. SLIP pik cT =0; ′ik cT = 0; ik cT =−µΨ ik cN + ′ik cN Tik n+1 2 (7) whe e KN,KTa e use –de ined penal ies o no mal and angen ial con ac , gik Nn+1,gik Nn, gik Tn+1,gik Tn no mal and angen ial gaps a nand n+ 1, Ψ = ±1 he Slip di ec ion, µ he ic ion coe icien and ψik 1N,ψik 1T,ψik 2Nand ψik 2Tp opo ionali y pa ame e s o he upda ing a iables. No ice ha pik cN , ′ik cN,pik cT , ′ik cT (a n+1/2) en o ce he conse a i e esponse o no mal and angen ial S ick con ac s. On he o he hand, o Slip pik cT = ′ik cT = 0 since no angen ial penal y sp ing is p esen ; he new Coulomb ic ion o ce ik cT is compu ed wi h he absolu e alue o he o al (con ac plus upda ing) no mal con ac o ces. The no mal and angen ial–S ick con ac o ces ik cN , ik cT a e de ined in Eqs. 8 using he [4] de i a i e, p o iding a disc e e exp ession ha conse es he a i icial penal y ene gy. i cN = nbd  k=1 k�=i ik cN = nbd  k=1 k�=i V(gik Nn+1)−V(gik Nn) gik Nn+1 −gik Nn Nik n+1 2= nbd  k=1 k�=i KNNik n+1 2(gik Nn+1 +gik Nn) i cT = nbd  k=1 k�=i ik cT = nbd  k=1 k�=i V(gik Tn+1)−V(gik Tn) gik Tn+1 −gik Tn Tik n+1 2= nbd  k=1 k�=i KTTik n+1 2(gik Tn+1 +gik Tn) whe e V(gik Nn+1)=KN(gik Nn+1)2/2,V(gik Nn)=KN(gik Nn)2/2 a e he no mal and V(gik Tn+1)= KT(gik Tn+1)2/2, V(gik Tn)=KT(gik Tn)2/2, he angen ial penal y po en ial con ac ene gies. 4 DISCRETE LINEAR, ANGULAR MOMENTUM CONSERVATION AND CONSISTENT BODY ENERGY DISSIPATION This sec ion de elops he disc e e conse ing amewo k o EFDA o ob ain he body ene gy conse a ion o no mal con ac and conse a ion–dissipa ion o angen ial. 4.1 Disc e e linea momen um balance The disc e e a ia ion o he linea momen um o EFDA is de ined h ough he second o Eqs. 6, he disc e e coun e pa o second New on’s law. The e o e, o body i he esul an o he no mal con ac o ces i cN , ′i cN plus he angen ial i cT , ′i cT is equal o he disc e e linea momen um balance be ween nand n+1. Also, he o al linea momen um 5 140 R. B a o and J.L. P´e ez–Apa icio balance (Eq. 8) is he summa ion o e ha o each body and equals he esul an o he con ac o ces on nbd. p o n+1 −p o n ∆ = nbd  i=1  i cN + ′i cN + i cT + ′i cT = nbd  i=1 nbd  k=1 k�=i ik cN + ′ik cN + ik cT + ′ik cT (8) Gi en wo bodies i, k in con ac , due o he AR p inciple: ik cN =− ki cN, ′ik cN =− ′ki cN , ik cT =− ki cT , ′ik cT =− ′ki cT o S ick and ik cT =− ki cT , ′ik cT = ′ki cT = 0 o Slip. Then, he igh e m o Eq. 8 is ze o in all si ua ions and p o n+1 =p o n. 4.2 Disc e e angula momen um balance We again ede ine he a iables qi n+1,qi nas he posi ions o he con ac poin . F om he second o Eqs. 6 and mul iplying by he c oss p oduc ×(qi n+1 −qi n), he disc e e angula momen um balance o a body iis: pi n+1 −pi n ∆ ×(qi n+1 −qi n)=Ji n+1 −Ji n ∆ = nbd  k=1 k�=i ik cN + ′ik cN + ik cT + ′ik cT ×(qi n+1 −qi n) (9) In oking he AR p inciple and exp essing he posi ion’s inc emen s as unc ion o he no mal gap (qi n+1 −qi n)−(qk n+1 −qk n)=gik Nn+1/2(Nki) , he o al angula momen um balance o nbd is: J o n+1 −J o n ∆ = nbd  i=1 nbd  k=i+1  ik cN + ′ik cN + ik cT + ′ik cT ×gik Nn+1 2(Nki) (10) Since ec o s ik cN + ′ik cN , and he no mal gap a e collinea , hei c oss p oduc is ze o. The p oduc be ween ik cT + ′ik cT and his gap is also ze o since he angen ial con ac o ces depend on he angen ial gap: a e some algeb a we a i e o he iple p oduc (Tik n+1/2) Tik n+1/2×gik Nn+1/2(Nki) ha is ze o since he i s ec o is o hogonal o he c oss p oduc . 4.3 Disc e e o al bodies’ ene gy balance This equa ion is ob ained by p emul iplying bo h Eqs. 6 by (pi n+1 −pi n) and −(qi n+1 − qi n) espec i ely, hen added o all con ac ing bodies nbd. A e some algeb a: 6 141 R. B a o and J.L. P´e ez–Apa icio ∆Ekin = nbd  i=1 (pi n+1 −pi n) M−1 ipi n+1 2   ∆Ei kin =− nbd  i=1 nbd  k=1 k�=i−(qi n+1 −qi n) ik cT   ∆E ik cT  + nbd  i=1 nbd  k=1 k�=i(qi n+1 −qi n) ik cN   ∆E ik cN −(pi n+1 −pi n) M−1 ipik cN   ∆Epik cN (ψik 2N) −(pi n+1 −pi n) M−1 ipik cT   ∆Epik cT (ψik 2T) −(qi n+1 −qi n) ′ik cN   ∆E ′ik cN (ψik 1N) −(qi n+1 −qi n) ′ik cT   ∆E ′ik cT (ψik 1T) (11) whe e ∆Ei kin =Ei n+1 −Ei nis he kine ic body ene gy balance be ween nand n+1 when bodies a e igid and ex e nal o ces a e no applied. Then, ∆E ik cN ,∆E ik cT a e he con ac o ces ene gy balance and ∆Epik cN (ψik 2N), ∆Epik cT (ψik 2T), ∆E ′ik cN (ψik 1N), ∆E ′ik cT (ψik 1T), all unc ions o he p opo ionali y pa ame e s, a e he upda ing a iables ene gy balance. This equa ion is he conse ing amewo k ha ela es he o al bodies’ ene gy balance wi h ha o he upda ing a iables. The ole o ene gy conse a ion o no mal con ac is included in he e ms ∆E ik cN ,∆E ′ik cN (ψik 1N), ∆Epik cN (ψik 2N), while dissipa ion–conse a ion o Slip and S ick is included in he e ms ∆E ik cT ,∆E ′ik cT (ψik 1T), ∆Epik cT (ψik 2T). The e o e, he ene gy loss is always consis en since he dissipa ion is included in he ene gy balance. No ice ha ∆Ekin is he o al ene gy o all bodies. Since he ene gy ela ed o no mal con ac is conse ed, ψik 1N,ψik 2Nmay be posi i e o nega i e, and ∆Epik cN ,∆E ′ik cN add o sub ac ene gy. The same can be said o ψik 1T,ψik 2T o he S ick and Slip cases: STICK: he ene gy is conse ed since he con ac o ce does no c ea e–dissipa e angen ial wo k. This condi ion is en o ced by ze oing he igh pa o Eq. 11: 0= nbd  i=1 nbd  k=1 k�=i∆E ik cT +∆E ik cN +∆Epik cN +∆Epik cT +∆E ′ik cN +∆E ′ik cT (12) This equa ion p o ides in ini e ela ionships ha sa is y he o al bodies’ ene gy con- se a ion o ψik 1N,ψik 2N,ψik 1T,ψik 2T. Using he AR p inciple and he ecip oci ies ψik 1N=ψki 1N, ψik 2N=ψki 2N,ψik 1T=ψki 1Tand ψik 2T=ψki 2T, Eq. 12 may be decoupled o he no mal (Eq. 13) and angen ial (Eq. 14) con ac s be ween bodies i, k as: (pi n+1 −pi n) M−1 ipik cN   ∆Epik cN (ψik 2N) +(pk n+1 −pk n) M−1 kpki cN   ∆Epki cN (ψik 2N) +(qi n+1 −qi n) −(qk n+1 −qk n)  ik cN   ∆E ik cN +(qi n+1 −qi n) −(qk n+1 −qk n)  ′ik cN   ∆E ′ik cN (ψik 1N) =0 (13) 7 142 R. B a o and J.L. P´e ez–Apa icio (pi n+1 −pi n) M−1 ipik cT   ∆Epik cT (ψik 2T) +(pk n+1 −pk n) M−1 kpki cT   ∆Epki cT (ψik 2T) +(qi n+1 −qi n) −(qk n+1 −qk n)  ik cT   ∆E ik cT +(qi n+1 −qi n) −(qk n+1 −qk n)  ′ik cT   ∆E ′ik cT (ψik 1T) =0 (14) Bo h imply ha he ene gy ans e ed o he no mal and angen ial penal y sp ings is eco e ed by he upda ing a iables. The ene gies ∆Epik cN ,∆Epik cT en o ce he o al bodies’ ene gy conse a ion, while ∆E ′ik cN ,∆E ′ik cT adjus he con ac o ces o he conse a i e solu ion. SLIP: om physical conside a ions, he o al ene gy dissipa ed by ic ion mus be equal o he inc emen o o al ene gy, En+1 −En=−nbd i=1 nbd k=1 k � =i ∆E ik cT . This equal- i y is en o ced by EFDA ze oing he las summa ion o Eq. 11. Also, ∆Epik cT (ψik 2T)= ∆E ′ik cT (ψik 1T) = 0 and ψik 1T=ψik 2T= 0 since he e is no penal y sp ing in he angen ial di- ec ion, see he Slip condi ion in Eq. 7. The e o e, he equa ion ha p o ides he in ini e ( he e o e unde e mined) ela ions be ween ψik 1N,ψik 2Nis: nbd  i=1 nbd  k=1 k�=i (pi n+1 −pi n) M−1 ipik cN   ∆Epik cN (ψik 2N) +(qi n+1 −qi n) ik cN   ∆E ik cN +(qi n+1 −qi n) ′ik cN   ∆E ′ik cN (ψik 1N)= 0 (15) ha ep esen s a no mal con ac ene gy balance o all bodies. As in he p e ious case, his balance is en o ced o e e y con ac ; using again AR, ecip oci ies ψik 1N=ψki 1N, ψik 2N=ψki 2Nand decoupling Eq. 15 o e e y con ac , one a i es o ela ionships be ween ψik 1N,ψik 2N: (pi n+1 −pi n) M−1 ipik cN   ∆Epik cN (ψik 2N) +(pk n+1 −pk n) M−1 kpki cN   ∆Epki cN (ψik 2N) + (qi n+1 −qi n) −(qk n+1 −qk n)  ik cN   ∆E ik cN +(qi n+1 −qi n) −(qk n+1 −qk n)  ′ik cN   ∆E ′ik cN (ψik 1N) =0 (16) The condi ion a he beginning o he SLIP i em and he las equa ion, en o ce bo h he ene gy conse a ion o he no mal con ac and he dissipa ion o angen ial con ac . 5 DYNAMIC CONTACT, ENHANCED PENALTY METHOD Fo e e y con ac , Eqs. 13, 14, 16 p o ide a non–unique ela ion be ween ψik 1N,ψik 2N and ψik 1T,ψik 2T espec i ely ha au oma icaly conse e o dissipa e consis en ly he o al 8 143 R. B a o and J.L. P´e ez–Apa icio ene gy o he angen ial S ick and Slip cases. Using modal analysis decomposi ion, [2], i is possible o ob ain he second o de dynamic equa ion associa ed wi h he desc ip i e algo i hm o Eqs. 6, de ining an enhanced penal y con ac model. The model desc ibed by Eqs. 17 consis s o a sp ing and dashpo ha con ol he gap and he pene a ion eloci y, espec i ely. 0=¨qNn+1/2+ 2ξNωN   ω2 N∆ ψ1N+ψ2N 2˙qNn+1/2+ω2 NqNn+1/2 0=¨qTn+1/2   Ine ia + 2ξTωT   ω2 T∆ ψ1T+ψ2T 2˙qTn+1/2   Dashpo +ω2 TqTn+1/2   Sp ing (17) The a iables qNn+1/2,qTn+1/2 ep esen he pa icle mo ion in no mal and angen ial di ec ion. The dashpo s a e con olled by he use –de ined pa ame e s ξN,ξT, a penaliza- ion o pene a ion eloci ies ha app oxima ely en o ce he consis ency Kuhn–Tucke condi ion. Eqs. 17 p o ide he ela ions ha can be easily gene alized o any con ac be ween igid bodies i, k: ψik 1N+ψik 2N=4ξN Ωik N ;ψik 1T+ψik 2T=4ξT Ωik T (18) whe e Ωik N=ωik N∆ ,Ω ik T=ωik T∆ ,ωik N=KN/m,ωik T=KT/m, and mis he la ges o he wo con ac ing masses. The combina ion o Eqs. 13, 14, 16 wi h Eq. 18 p o ide he unique explici exp essions o ψik 1N,ψik 2Nand ψik 1T,ψik 2T. The inse ion o hese exp essions in Eqs. 6 en o ces a conse a i e esponse o S ick and consis en dissipa i e o Slip. 6 NUMERICAL SIMULATIONS 6.1 Ellip ical pa icle Ca om p oblem In his subsec ion, he ajec o y o he successi e impac s o a igid ellipse inside a one–me e squa e is simula ed. The ellipse, o axes 15/6 cm is ini ially posi ioned a (0.45,0.1) m, inclina ion α= 50◦as seen in Fig. 1 op, and i is subjec ed o ini ial eloci y Vx=1,V y=−0.4 m/s in di ec ion θ=−22◦, wi hou spin. The ic ion angle is φ= 15◦and he es o he nume ical pa ame e s a e he same as hose in he p e ious simula ion. To isualize he o a ion o he ellipse, he o ien a ion is de ined by he la ges semiaxis. Figs. 1 depic he e olu ion o ajec o y ( op), linea eloci ies (le ), o a ional e- loci y and ene gy ( igh ) om EFDA and om he analy ical solu ion ep oduced in he 9