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