A Nitsche nite element method for dynamic contact

(1)

contact

Franz Chouly

Laboratoire de Mathématiques de Besançon UMR CNRS 6623, France

Patrick Hild

Institut de Mathématiques de Toulouse UMR CNRS 5219, France

Yves Renard

Université de Lyon, CNRS, INSA-Lyon ICJ UMR5208, LaMCoS UMR5259, France

FBP 2017, Saint-Etienne, 20th of June, 2017

1/30

(2)

Analysis of the semi-discrete problem Time-marching schemes

Numerical experiments Conclusions

2/30

(3)

Setting: contact in linear elastodynamics

For the sake of simplicity, we consider the classical simple situation of the equilibrium of a linearly elastic body in contact with (static) Coulomb friction with a rigid foundation. Note that the main interest of this problem is to be similar to the one obtained by applying a time integration scheme on a quasi-static problem.

Let Ω⊂R^d, d∈ {2,3}be a bounded regular domain which represents the reference configu- ration of a linearly elastic body submitted to a Neumann condition on Γ_N, a Dirichlet condition on Γ_Dand a unilateral contact condition with friction on Γ_Cwith a rigid foundation, where Γ_N, Γ_Dand Γ_Care non-overlapping open parts of∂Ω, the boundary of Ω (see Fig. 1). The part Γ_D is supposed of non-zero measure in∂Ω.

.

Ω

Rigid foundation Γ_C ΓD

g Γ_N

Γ_N n0

n

.

Figure 1:Linearly elastic bodyΩin contact with a rigid foundation.

The displacementu(t, x) of the body obeys the following equations:

−divσ(u) =f, in Ω, σ(u) =Aε(u), in Ω, σ(u)n=k, on Γ_N, u= 0, on Γ_D, (1) whereσ(u) is the stress tensor,ε(u) = (∇u+∇u^T)/2 is the linearized strain tensor,Ais the fourth order elasticity tensor which satisfies usual conditions of symmetry, coercivity and boundedness, nis the outward unit normal to Ω on∂Ω andf,kare given force densities. On Γ_C, it is usual to decompose the displacement and the stress in normal and tangential components as follows, assuming the shape of the rigid foundation to have theC¹regularity:

u_N=−u.n0, u_T=u+u_Nn0, σ_N(u) =−(σ(u)n).n0, σ_T(u) =σ(u)n+σ_N(u)n0, wheren0is the unit outward normal to the obstacle (see Fig. 1). Denoting bygthe initial normal gap between the solid and the rigid obstacle, the unilateral contact condition is expressed by the following complementary condition:

u_N−g≤0, σ_N(u)≤0,(u_N−g)σ_N(u) = 0, (2) while the static Coulomb friction condition is expressed as follows forFthe friction coefficient:

|σ_T| ≤ −Fσ_N, ifu_T&= 0 then σ_T=−Fσ_N u_T

|u_T|. (3)

A classical weak formulation (see [14]) can be obtained introducing V=!

v∈H¹(Ω;R^d) : v= 0 on Γ_D"

, K={v∈V:v_N−g≤0 on Γ_C}, 3

Ω⊂R^d (d=2,3)

Find u: [0,T)×Ω→R^d s.t.:

ρ¨u−divσ(u) = f in(0,T)×Ω u= 0 on(0,T)×Γ_D σ(u)n= g on(0,T)×Γ_N u(0,·) = u0 inΩ

u˙(0,·) = u˙0 inΩ

Contact without friction on(0,T)×Γ_C:

( un≤0 σ_n(u)≤0 σ_n(u)un=0 (contact)

σ_t(u) =0 (no friction)

with u=u_nn+u_t ,σ(u)n=σ_n(u)n+σ_t(u).

3/30

(4)

Coulomb frictional case (σ_t(u)6=0) on (0,T)×Γ_C:

u˙_t=0 =⇒ |σ_t(u)| ≤ −µ σ_n(u) u˙t6=0 =⇒ σ_t(u) =µ σ_n(u) u˙t

|u˙_t| µ≥0 : friction coecient (µ=0↔ no friction)

Remark:

I No friction→ energy conservation.

I Friction→ energy dissipation when slip.

4/30

(5)

(Scalar) wave equation with unilateral constraints:

1. One dimensional case→ one contact point:

Existence and uniqueness.

(Schatzman, 1980 ; Dabaghi, Petrov, Pousin, Renard, 2014) 2. N ≥2 dimensional case:

Existence and uniqueness (half space).

(Lebeau, Schatzman, 1984)

Existence (bounded domain with smooth boundary).

(Kim, 1989)

5/30

(6)

Discretization: method of lines.

Finite elements for space variables + time-marching schemes.

Contact condition: Lagrange multipliers or penalty.

Diculties for the space semi-discrete problem:

1. illposedness with mixed FE discretization of contact conditions

(Khenous, Laborde, Renard, 2008)

→ modied mass method

2. wellposedness of discrete systems with springs (Ballard, 2000 ; Ballard, Charles, 2014)

6/30

(7)

Discretization: method of lines.

Finite elements for space variables + time-marching schemes.

Contact condition: Lagrange multipliers or penalty.

Diculties for the space-time discrete problem and for some conservative time-stepping (e.g., Crank-Nicolson).

1. Spurious oscillations on contact pressure, velocity and displacement.

2. No energy conservation and possible blow up.

(Khenous, Laborde, Renard, 2008 ; Doyen, Ern, Piperno, 2011 ; Krause, Walloth, 2012)

7/30

(8)

(Non-exhaustive) list of strategies:

1. Modify the model Impact law. (Paoli, 2001) 2. Modify the time-discretization

2.1 Implicit (dissipative) time-discretization of contact term.

(Carpenter-et-al, 1991 ; Kane-et-al, 1999 ; Dumont & Paoli, 2006 ; Deuhard-et-al, 2008)

2.2 Velocity update method / contact condition with velocity.

(Laursen & Love, 2002 ; Laursen & Chawla, 1997) 2.3 Penalty with energy conservative schemes.

(Armero & Petöcz, 1998 ; Hauret & Le Tallec, 2006) 3. Modify the space discretization

Modied mass method and extensions.

(Khenous, Laborde & Renard, 2008 ; Hager, Hüeber & Wohlmuth, 2008 ; Hauret, 2010 ; Renard, 2010)

8/30

(9)

8/30

(10)

8/30

(11)

Key idea: contact condition can be reformulated as σ_n(u) =−_γ¹[u_n−γσ_n(u)]₊ (Alart & Curnier, 1988)

NitscheFEM for (frictionless) dynamic contact:











Find u^h: [0,T]→V^h s.t.:

hρ¨u^h(t),v^hi+AΘγ_h(u^h(t),v^h) + Z

Γ_C

1

γ_h [Pγ_h(u^h(t))]₊PΘγ_h(v^h)dΓ

= L(t)(v^h) ∀v^h∈V^h u^h(0,·) =u^h₀ u˙^h(0,·) = ˙u^h₀

Notations:

I γ_h:=γ₀h

I PΘγ_h(v^h) :=v_n^h−Θγ_hσ_n(v^h)

I AΘγ_h(u^h,v^h) :=a(u^h,v^h)− Z

Γ_C

Θγ_hσ_n(u^h)σ_n(v^h)dΓ V^h: H¹-conformal FE space (Pkcont, k=1,2).

9/30

(12)











Γ_C

1

Notations:

I γ_h:=γ₀h

Γ_C

9/30

(13)











Γ_C

1

Notations:

I γ_h:=γ₀h

Γ_C

9/30

(14)

NitscheFEM for dynamic contact:











Γ_C

1

Two numerical parameters:

1. γ₀ (inγ_h:=γ₀h): Nitsche's parameter.

Should be small enough for well-posedness (in the fully discrete framework).

2. Θ: Select a variant.

Skew-symmetric (−1), non-symmetric (0) or symmetric (1).

(F.C., Hild & Renard, 2014)

Remark: consistency for allγ₀ >0 (6=penalty).

10/30

(15)











Γ_C

1

10/30

(16)











Γ_C

1

10/30

(17)

No additional diculty ...

since the frictional contact conditions can be rewritten:

σ_n(u) = −1

γ[u_n−γσ_n(u)]₊ σ_t(u) = −1

γ[ ˙u_t−γσ_t(u)]_−γµσ_n₍_u₎

= −1

γ[ ˙u_t−γσ_t(u)]_µ[_u_n_−γσ_n₍_u_)]₊

11/30

(18)

NFEM for dynamic contact can be recasted as:







M^h¨u^h(t) +B^hu^h(t) =L^h(t) u^h(0,·) =u^h₀ u˙^h(0,·) = ˙u^h₀ with

(M^hv^h,w^h)_γ_h =hρv^h,w^hi (B^hv^h,w^h)γ_h =AΘγ_h(v^h,w^h) +

Z

Γ_C

1

γ_h [Pγ_h(v^h)]+PΘγ_h(w^h)dΓ (L^h(t),w^h)_γ_h =L(t)(w^h)

for all v^h,w^h∈V^h, and

(v^h,w^h)_γ_h := (v^h,w^h)₁_,Ω+ (γ_h⁻¹²v_n^h, γ_h⁻¹²w_n^h)₀_,Γ_C.

12/30

(19)

(?)







Find u^h: [0,T]→V^h s.t. : M^h¨u^h(t) +B^hu^h(t) =L^h(t) u^h(0,·) =u^h₀ u˙^h(0,·) = ˙u^h₀ Theorem (F.C., Hild, Renard)

B^his Lipschitz-continuous:

kB^hv^h₁−B^hv₂^hkγ_h ≤C(1+γ₀)(1+|Θ|)kv^h₁−v^h₂kγ_h

for all v₁^h,v^h₂ ∈V^h; C >0 independent of h,Θandγ₀. So, for allΘ∈Randγ₀ >0, Problem (?) admits one unique solution u^h∈C²([0,T],V^h).

Remark: the same result applies in the frictional case.

13/30

(20)

(?)







Find u^h: [0,T]→V^h s.t. : M^h¨u^h(t) +B^hu^h(t) =L^h(t) u^h(0,·) =u^h₀ u˙^h(0,·) = ˙u^h₀ Theorem (F.C., Hild, Renard)

B^his Lipschitz-continuous:

kB^hv^h₁−B^hv₂^hkγ_h ≤C(1+γ₀)(1+|Θ|)kv^h₁−v^h₂kγ_h

for all v₁^h,v^h₂ ∈V^h; C >0 independent of h,Θandγ₀. So, for allΘ∈Randγ₀ >0, Problem (?) admits one unique solution u^h∈C²([0,T],V^h).

Remark: the same result applies in the frictional case.

13/30

(21)

2. Energy conservation

Proposition (F.C., Hild & Renard)

For a conservative system, i.e., L(t)≡0 for all t ∈[0,T], the solution u^h of NFEM semi-discretization satises:

d

dtE_Θ^h = (Θ−1) Z

Γ_C

1

γ_h [Pγ_h(u^h)]₊u˙^h_ndΓ where

E_Θ^h :=E^h−Θ 2

hkγ_h¹²σ_n(u^h)k²₀,Γ_C − kγ_h⁻¹²[Pγ_h(u^h)]+k²₀,Γ_C

i

E^h:= 1

2ρku˙^hk²₀,Ω+ 1

2a(u^h,u^h) WhenΘ =1, the augmented energy E_Θ^h is conserved.

14/30

(22)

2. Energy conservation

For a conservative system, i.e., L(t)≡0 for all t ∈[0,T], the solution u^h of NFEM semi-discretization satises:

d

dtE_Θ^h = (Θ−1) Z

Γ_C

1

γ_h [Pγ_h(u^h)]₊u˙^h_ndΓ where

E_Θ^h :=E^h−Θ 2

hkγ_h¹²σ_n(u^h)k²₀,Γ_C − kγ_h⁻¹²[Pγ_h(u^h)]+k²₀,Γ_C

i

E^h:= 1

2ρku˙^hk²₀,Ω+ 1

2a(u^h,u^h) WhenΘ =1, the augmented energy E_Θ^h is conserved.

14/30

(23)

Parameter isθ∈[0,1].

For n≥0, the fully discretized problem reads:











Find u^h^,ⁿ⁺¹,u˙^h^,ⁿ⁺¹,u¨^h^,ⁿ⁺¹ ∈V^hsuch that:

u^h^,ⁿ⁺¹ =u^h^,ⁿ+τu˙^h^,ⁿ^+θ, u˙^h,n⁺¹ = ˙u^h,n+τ¨u^h,n^+θ,

hρ¨u^h^,ⁿ⁺¹,v^hi+AΘγ_h(u^h^,ⁿ⁺¹,v^h) + Z

Γ_C

1

γ_h [Pγ_h(u^h^,ⁿ⁺¹)]₊PΘγ_h(v^h)dΓ

=Lⁿ⁺¹(v^h), ∀v^h∈V^h.

Well posedness whenθ=0 or when(1+ Θ)²γ₀ ≤C

1+ _τ^ρ₂^h_θ²₂ . Stability (augmented energy) whenΘ =1 (symmetric variant) and θ=1 (backward Euler).

15/30

(24)

Parameters areβ ∈[0,1/2],γ ∈[0,1]. For n≥0, the fully discretized problem reads:











Find u^h^,ⁿ⁺¹,u˙^h^,ⁿ⁺¹,u¨^h^,ⁿ⁺¹ ∈V^hsuch that:

u^h^,ⁿ⁺¹ =u^h^,ⁿ+τu˙^h^,ⁿ+τ²

2 ¨u^h^,ⁿ⁺²^β, u˙^h,n⁺¹ = ˙u^h,n+τ¨u^h,n^+γ,

hρ¨u^h,n+1,v^hi+AΘγ_h(u^h,n+1,v^h) + Z

Γ_C

1

γ_h [Pγ_h(u^h,n+1)]₊PΘγ_h(v^h)dΓ

=Lⁿ⁺¹(v^h), ∀v^h∈V^h.

Well posedness whenβ =0 or when(1+ Θ)²γ₀ ≤C

1+_τ^ρ^h₂_β² . Stability (augmented energy) whenΘ =1 (symmetric variant) and β=1/2,γ =1.

16/30

(25)

To summarize:

1. For the θ-scheme:

numerical stability for arbitrary time-steps τ only forθ=1.

2. Newmark:

numerical stability for arbitrary time-steps τ only forγ =1, β = ¹₂.

3. Especially for Crank-Nicolson: no theoretical proof and no numerical evidence of (unconditional) stability.

4. Can we design a time-marching scheme that is stable and (almost) conservative ?

17/30

(26)

To summarize:

2. Newmark:

17/30

(27)

To summarize:

2. Newmark:

17/30

(28)

To summarize:

2. Newmark:

17/30

(29)

Hybrid time-discretization of NFEM for dynamic contact:











Find u^h,n⁺¹,u˙^h,n⁺¹,¨u^h,n⁺¹∈V^h s.t.:

u^h^,ⁿ⁺¹=u^h^,ⁿ+τu˙^h^,ⁿ⁺¹² u˙^h,n⁺¹= ˙u^h,n+τ¨u^h,n⁺¹²

hρ¨u^h,n+¹²,v^hi+AΘγ_h(u^h,n+¹²,v^h) + Z

Γ_C

1

γ_h Φ(u^h,n,u^h,n+1)PΘγ_h(v^h)dΓ =Lⁿ⁺¹²(v^h)

∀v^h∈V^h

Initial conditions: u^h,0=u^h0,u˙^h,0= ˙u^h0,¨u^h,0= ¨u^h0.

Hybrid Crank-Nicolson/Midpoint discretization of contact term:

Φ(u^h,n,u^h,n+1) :=H(Pγ_h(u^h,n))[Pγ_h(u^h,n+¹²)]++H(−Pγ_h(u^h,n))[Pγ_h(u^h)]ⁿ⁺₊ ¹² H(·)is the Heaviside function.

x^h,n⁺¹² =¹₂x^h,n+¹₂x^h,n⁺¹

(Gonzalez, 2000, Hauret & Le Tallec, 2006, F.C., Hild & Renard, 2015)

18/30

(30)











u^h^,ⁿ⁺¹=u^h^,ⁿ+τu˙^h^,ⁿ⁺¹² u˙^h,n⁺¹= ˙u^h,n+τ¨u^h,n⁺¹²

Γ_C

1

∀v^h∈V^h

Hybrid Crank-Nicolson/Midpoint discretization of contact term:

Φ(u^h,n,u^h,n+1) :=H(Pγ_h(u^h,n))[Pγ_h(u^h,n+¹²)]++H(−Pγ_h(u^h,n))[Pγ_h(u^h)]ⁿ⁺₊ ¹² H(·)is the Heaviside function.

x^h,n⁺¹² =¹₂x^h,n+¹₂x^h,n⁺¹

(Gonzalez, 2000, Hauret & Le Tallec, 2006, F.C., Hild & Renard, 2015)

18/30

(31)

(HN)











u^h^,ⁿ⁺¹=u^h^,ⁿ+τu˙^h^,ⁿ⁺¹² u˙^h,n+1= ˙u^h,n+τ¨u^h,n+¹²

Γ_C

1

∀v^h ∈V^h

Proposition (F.C., Hild, Renard) If(1+ Θ)²γ₀ ≤C

1+^ρh_τ₂²

, then for all n, Problem (HN) admits one unique solution.

19/30

(32)

(HN)











u^h^,ⁿ⁺¹=u^h^,ⁿ+τu˙^h^,ⁿ⁺¹² u˙^h,n+1= ˙u^h,n+τ¨u^h,n+¹²

Γ_C

1

∀v^h ∈V^h

Proposition (F.C., Hild, Renard) If(1+ Θ)²γ₀ ≤C

1+^ρh_τ₂²

, then for all n, Problem (HN) admits one unique solution.

19/30

(33)

Suppose Lⁿ ≡0 for all n≥0 and that Problem (HN) is well-posed. Ifγ₀small (orΘ =−1), then, for all n≥0:

EΘ^h,n⁺¹−EΘ^h,n

=−Θ Z

Γ_C

1

2γ_h H(Pⁿ)H(Pⁿ+Pⁿ⁺¹)[Pⁿ⁺¹]²−

+H(Pⁿ)H(−Pⁿ−Pⁿ⁺¹)[Pⁿ]²₊+ [Pⁿ]−[Pⁿ⁺¹]+

dΓ + (Θ−1)Z

Γ_C

1

2γh H(Pⁿ)[Pⁿ+Pⁿ⁺¹]++H(−Pⁿ) [Pⁿ]++ [Pⁿ⁺¹]+

u_n^h,n⁺¹−u^h,n_n dΓ

with Pⁿ=Pγ_h(u^h,n), Pⁿ⁺¹=Pγ_h(u^h,n⁺¹)and EΘ^h^,ⁿ:=E^h^,ⁿ−Θ

2

hkγ_h¹²σ_n(u^h^,ⁿ)k²_0,Γ_C − kγ_h⁻¹²[Pγ_h(u^h^,ⁿ)]+k²_0,Γ_Ci E^h^,ⁿ:=1

2ρku˙^h^,ⁿk²_0,Ω+1

2a(u^h^,ⁿ,u^h^,ⁿ)

20/30

(34)

Corollary (F.C., Hild & Renard)

Suppose Lⁿ≡0 for all n≥0 and that Problem (HN) is well-posed.

Ifγ₀ small and Θ =1, then, for all n≥0:

E_Θ^h,n⁺¹−E_Θ^h,n

=−Θ Z

Γ_C

1

2γ_h H(Pⁿ)H(Pⁿ+Pⁿ⁺¹)[Pⁿ⁺¹]²₋ +H(Pⁿ)H(−Pⁿ−Pⁿ⁺¹)[Pⁿ]²₊+ [Pⁿ]−[Pⁿ⁺¹]₊

dΓ≤ 0 So the hybrid scheme is unconditionally stable (for all h, τ >0).

21/30

(35)

Multiple impacts of an elastic bar

x= 0 L

t= 0 t1= 1 t2= 2 t3= 3

I Parameters:

f =0, E =1,ρ=1, L=1, u₀(x) = ¹₂ −^x₂,u˙₀(x) =0

I Periodic analytical solution with multiple impacts.

(Dabaghi-et-al, 2014)

I Numerical experiments with GetFEM++.

http://download.gna.org/getfem/html/homepage/

I Generalized Newton to solve the non-linear problem.

22/30

(36)

Mixed(+Crank-Nicolson) vs. Nitsche(+Crank-Nicolson) vs.

Nitsche-Hybrid

0 2 4 6 8 10 12

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

t

u

exact solution mixed Nitsche Nitsche−Hybrid

9.8 10 10.2 10.4 10.6 10.8 11 0

0.005 0.01 0.015 0.02 0.025 0.03

t

u

exact solution mixed Nitsche Nitsche−Hybrid

Displacement

Symmetric variantΘ =1, with γ₀ =10⁻⁶, 100 nite elements (h=0.01), τ =0.015.

23/30

(37)

Nitsche-Hybrid

0 2 4 6 8 10 12

80 100 120 140 160 180 200 220

t

Eh

mixed Nitsche Nitsche−Hybrid

Energy

24/30

(38)

Nitsche-Hybrid

0 2 4 6 8 10 12

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2 0

t σC

exact solution mixed

0 2 4 6 8 10 12

−1

−0.8

−0.6

−0.4

−0.2 0

t σC

exact solution Nitsche

0 2 4 6 8 10 12

−1

−0.8

−0.6

−0.4

−0.2 0

t σC

exact solution Nitsche−Hybrid

Contact pressure

25/30

(39)

Convergence when h, τ →0

10⁻³ 10⁻² 10⁻¹ 10⁰

10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

τ kehτkl2(τ;L2(Ω))

mixed (slope=0.5) Nitsche (slope=0.5) Nitsche−Hybrid (slope=0.68) modified mass (slope=0.73)

10⁻³ 10⁻² 10⁻¹ 10⁰

10⁻¹ 10⁰ 10¹

τ kehτkl2(τ;H1(Ω))

mixed (slope=−0.2) Nitsche (slope=−0.2) Nitsche−Hybrid (slope=0.24) modified mass (slope=0.36)

Error curves for u: error e^τh in norm l²(τ;L²(Ω))and l²(τ;H¹(Ω)) Ratioτ /h is kept constant.

26/30

(40)

Convergence when h, τ →0

10⁻³ 10⁻² 10⁻¹ 10⁰

10⁻¹ 10⁰ 10¹

τ kehτ σkl2(τ)

mixed (slope=0.011) Nitsche (slope=0.021) Nitsche−Hybrid (slope=0.29) modified mass (slope=0.36)

10⁻³ 10⁻² 10⁻¹ 10⁰

10⁻² 10⁰ 10² 10⁴

τ kehτEkl2(τ)

mixed (slope=−0.42) Nitsche (slope=−0.42) Nitsche−Hybrid (slope=1.1) modified mass (slope=0.95)

Error curves for the contact pressure σ_C and for the energy E.

Ratioτ /h is kept constant.

27/30

(41)

Impact of a ball

0 10 20 30 40

0 200 400 600 800 1000

t

Eh

mixed Nitsche Nitsche−Hybrid modified mass

Deformed conguration and von Mises strain at t =18, for Nitsche-Hybrid. Energy E^h for dierent methods.

Parameters: D =40,λ=20,µ=20, ρ=1,kfk=0.1.

P2 Lagrange FEM, h'8 (400 elements), τ =0.1.

Symmetric variantΘ =1, with γ₀ =0.001.

28/30