A dual Lagrange method for contact problems with regularized frictional contact conditions: Modelling micro slip

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487 www.elsevier.com/locate/cma

A dual Lagrange method for contact problems with regularized frictional contact conditions: Modelling micro slip Saskia Sitzmann a,∗ , Kai Willner b , Barbara I. Wohlmuth c a Central Institute for Scientific Computing, University Erlangen-Nuremberg, Martensstraße 5a, 91058 Erlangen, Germany b Chair of Applied Mechanics, University Erlangen-Nuremberg, Egerlandstraße 5, 91058 Erlangen, Germany c Institute for Numerical Mathematics, Technical University Munich, Boltzmannstraße 3, 85748 Garching, Germany

Received 15 July 2014; received in revised form 2 November 2014; accepted 13 November 2014 Available online 26 November 2014

Abstract This paper presents an algorithm for solving quasi-static, non-linear elasticity contact problems with friction in the context of rough surfaces. Here, we want to model the transition from sticking to slipping also called micro slip in a physically correct way in order to reproduce measured frictional damping. The popular dual Mortar method is used to enforce the contact constraints in a variationally consistent way without increasing the algebraic system size. The algorithm is deduced from a perturbed Lagrange formulation and combined with a serial–parallel Iwan model. This leads to a regularized saddle point problem, for which a non-linear complementary function and thus a semi-smooth Newton method can be derived. Numerical examples demonstrate the applicability to industrial problems and show good agreement to experimentally obtained results. c 2014 Elsevier B.V. All rights reserved. ⃝

Keywords: Contact with friction; Dual Mortar methods; Micro slip; Rough surfaces; Constitutive contact equations; Iwan model

1. Introduction Since solving non-linear contact problems within the FEM framework is still a challenging task from both the mathematical and the engineering point of view, a lot of research has been done in this field in the past decades. For an overview to contact problems in general, we refer to the monographs by Johnson [1] or Bowden/Tabor [2] for contact between rough surfaces, the monographs by Kikuchi/Oden [3] and Eck [4] for existence and uniqueness results and to the monographs by Willner [5], Laursen [6] and Wriggers [7] for computational aspects. The most popular approaches to enforce the contact condition are the penalty method and the Lagrange multiplier method. Within the penalty framework, a spring between a slave location and a master surface is generated modelling the contact traction tcs = G n (u) = an · u n as a function of the displacements with the spring constant an , if the non-penetration condition is violated. This approach is often combined with a node-to-surface formulation. Though easy to implement, since the contact forces do not enter as an extra variable in the system, it has the drawback of an ∗ Corresponding author.

E-mail address: [email protected] (S. Sitzmann). http://dx.doi.org/10.1016/j.cma.2014.11.022 c 2014 Elsevier B.V. All rights reserved. 0045-7825/⃝

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

469

ill-conditioned matrix system for stiff contact laws, i.e., an → ∞ and does not satisfy the patch test. The Lagrange multiplier method introduces the contact traction as an extra variable with the benefit of being able to model arbitrary stiff contact laws. Combined with a surface-to-surface discretization, the Lagrange multiplier formulation leads to more robust algorithms, since the contact condition is not enforced point-wise in the strong form, but weakly in an integral formulation over the surrounding surface of a point. The drawback, however, is the need for solving a saddle point system with respect to the displacements and the contact traction in every Newton–Raphson iteration. If contact problems between two deformable bodies and non-matching meshes are considered, the Lagrange multiplier setting is strongly related to mortar methods introduced originally in the context of domain decomposition techniques for non-matching meshes by Bernardi [8] and Belgacem [9]. The Mortar methods are also characterized by the introduction of an additional variable to model the traction across the interface and a weak formulation of the boundary constraints leading to a surface-to-surface approach. It was adapted to small deformation contact problems by Belgacem [10]. The idea of the weak formulation of the contact condition was later expanded to large deformation contact [11,12], for curved interfaces and large sliding [13–15] and dynamical problems [16,17]. Among the Mortar methods, the dual Mortar method [18,19] has become of great interest for solving contact problems in recent years [20–23] since it enforces the non-penetration condition in a variationally consistent way without increasing the system size. Due to a smart choice of the basis functions for the contact traction, the Lagrange multiplier can be condensed from the global system before solving without losing the optimality of the solution. Combined with a semi-smooth Newton method [24,25], which interprets the contact conditions as a semi-smooth non-linear function, one obtains a powerful tool for contact problems since all non-linearities of the system (material and contact conditions) can be handled within the same iteration loop. So far a perfect Coulomb law has been considered in most publications using the Mortar formulation. More complicated friction models like thermomechanically coupled friction were explored by H¨uber [26] and Temizer [27]. Yang and Laursen [28] have considered the full Reynolds approximation for lubricated contact. But for use in structural dynamics it is important to model the transition between sticking and slipping, also called micro slip, in a physical correct way [29,30]. With this smooth friction model one can reproduce measured frictional damping [31], like for example in the bolted connection in the casing of a jet turbine. Thus there is a need for more general contact laws, such as constitutive contact laws for rough surfaces, in combination with suitable numerical algorithms. Traditionally constitutive contact laws are realized within a penalty framework with all resulting drawbacks mentioned above. Thus we propose a new generalized dual Mortar formulation to solve numerically problems with regularized frictional contact conditions. In the Lagrange multiplier framework, the smooth transition between sticking and slipping can be easily realized and due to the dual formulation the Lagrange multipliers can be condensed from the global system. Therefore, we extend the ideas presented in our previous paper [32] to regularized frictional contact laws. The outline of the present work is as follows: In Section 2, the fundamental equations in the strong form and their weak formulation are summarized. The contact behaviour of rough surfaces and their realization within the dual Mortar framework is presented in Section 3. In Section 3.1 an elastic stick part is added to the Coulomb friction law and a semi-smooth Newton for the active set iteration is derived. Finally, the Lagrange multipliers are locally eliminated from the matrix representation and the displacement based reduced algebraic system is given. In Section 3.2 the ideas from Section 3.1 are combined with a serial–parallel Iwan model [33] in order to get a smooth transition from sticking to slipping. Again a semi-smooth Newton is derived. Section 4 illustrates the robustness and applicability to real life contact problems of the newly derived dual Mortar formulation. First, we consider a simple mathematical benchmark for which analytical solutions exist, and then we present a resonator example with a bolted connection, for which measurements exist. Some conclusions are given in Section 5. 2. Problem definition 2.1. Strong form We start with a brief overview of the 3D frictionless two-body contact problem in a non-linear elasticity setting [3,5–7]. As one can see in Fig. 1, two contacting bodies are considered, a slave body, denoted by s, and a master body, denoted by m, with domains Ω α ∈ R3 , α = s, m in the reference configuration. The bodies are undergoing motion

470

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

Fig. 1. Two-body contact problem.

during the time interval [0, T ] described by the mapping ϕ α : Ω α × [0, T ] →R3 X α → x α = ϕ α (X α , t)

α = s, m,

Xα

which maps material points of the reference configuration to the current configuration x α . Let us define the displacements as u(X, t) = ϕ(X, t)− X with u = (u s , u m ) defined on Ω = Ω s ×Ω m . The boundary of the bodies ∂Ω α is divided into three disjoint sets, the Dirichlet boundary Γ Dα with given displacements, the Neumann boundary Γ Nα with given traction and the potential contact boundary Γcα . Their spatial counterparts on ∂Ωtα are denoted by γ Dα , γ Nα and γcα . It holds Γ¯ Dα ∪ Γ¯ Nα ∪ Γ¯ cα = ∂Ω α . In terms of the displacements u, the deformation gradient F = ∇[X +u(X, t)], the right Cauchy–Green tensor C = F T F and the Green–Lagrange strain tensor E = 21 (C − 1) = 12 (∇u + ∇u T + ∇u ∇u T ) can be defined. Definition 1 (The Strong Form of the Two-Body Contact Problem). Find u α : Ω α × [0, T ] → R3 satisfying  dynamic equation ∇ · P α + ρ0α b¯ α = ρ0α v˙ α , in (0, T ] × Ω α ,  α in (0, T ] × Γ Dα , u = u¯ α , α α ¯ boundary conditions P n = t , in (0, T ] × Γ Nα ,  α P n = tcα , in (0, T ] × Γcα ,  α u (0, ·) = u α0 , in Ω α , initial conditions α v (0, ·) = v0α ,

(1a) (1b)

(1c)

where P is the first Piola–Kirchhof stress tensor. The dynamic equation describes the balance of momentum with ∇ · P α being the internal forces, ρ0α b¯ α being the body forces and ρ0α v˙ α being the inertia forces acting on Ω α . On the Dirichlet boundary, we assume u¯ α = 0 and on the Neumann boundary P α n = t¯α , where t¯α is a given surface force. The contact traction tcα and the actual contact zone are not known a priori and have to be determined numerically. Remark 1. Material non-linearities are considered assuming the existence of a strain–energy W with ∂∂W E = S being the second Piola–Kirchhoff stress tensor. Since we will focus on quasi-static calculations, the term modellingthe iner tia forces will be set to zero, and the time interval [0, T ] will be split into a sequence of small load increments ti , ti+1 . For the contact constraints, we define for each point on the slave surface the orthogonal system (n, τ1 , τ2 ), with n = n(X, t) being the current normal on the slave surface in the point X at time t and τ = (τ1 , τ2 ) the corresponding tangents, see Fig. 1. In terms of this nomenclature, we define the projection operator Pt : Γcs → Γcm , X −→ Pt (X ) which takes a point X of the slave surface in the current configuration and projects it onto the master surface in direction of the actual normal n(X, t) in X of the slave surface, see Fig. 1. With the Projection operator Pt the gap function in the normal direction can be defined for the current configuration: G n (X s , t) := [X s + u(X s , t)]n = n T · (ϕ s (X s , t) − ϕ m (Pt (X s ), t)). Now the normal contact conditions for hard contact on Γcs can be stated as follows: s s G n ≤ 0 ∧ − tcn ≥ 0 ∧ − tcn · Gn = 0

471

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

s = n T · t s . The hard contact conditions are given here for simplicity. For a regularized version of the normal with tcn c contact, we refer to our previous work [32] or to [1] for more general aspects. For the tangential contact, we use the incremental form of Coulombs law [6] to resolve the path dependency: s s s s Φ(tcn , tcτ ) := ∥tcτ ∥ + µtcn ≤0

∧

˙ τ = −γ˙ [u]

γ˙ ≥ 0

∧

γ˙ · Φ = 0,

s tcτ , s ∥ ∥tcτ

(2)

s = (t , t ), t = τ T · t s , the friction coefficient µ, 0 < µ ≤ 1 and [u] with tcτ ˙ τ = τ T · (u˙ s (X s , t) − u˙ m (Pt (X s ), t)). i 1 2 c i The perfect Coulomb law is given here for the sake of completeness. A regularized version will be used in Section 3.

2.2. From weak formulation to discretization Here we recall briefly the weak formulation of the dynamic equation, the contact conditions will be handled separately in Section 3. We start with the balance of momentum equation (1a), multiply it with test functions δu α and integrate over the domain Ω α . Performing partial integration and inserting the boundary conditions yields                 α α α α α α α ¯α α α tc δu d S = 0. P · Gradδu d V − ρ0 (b − v˙ ) · δu d V − t¯ δu d S −   Ωα Γcα Ωα Γ Nα  α∈{s,m}           α α δ Πint,ext

δ ΠC

It follows a short functional analytical classification of the used function spaces, for further considerations see [18]. The displacements u α and the test functions for the displacements δu α are lying in the Sobolev-space    d V0α := δu α ∈ H 1 (Ω α ) : u α |Γ Dα = 0 .  d  d On the boundary ∂Ω s the displacements u s are elements of the trace space of H 1 (Ω s ) , namely W := H 1/2 (ΓCs ) . The contact traction, later identified with the negative Lagrange multiplier −λ, is set to be in the dual space of W , M := W −1 , in order to use the duality pairing  ⟨η, ν⟩ = η νd S, η ∈ W ν ∈ M ΓCs

in the discrete weak formulation. Term for internal and external forces. The term s m (u m , δu m )) δΠint/ext (u, δu) := (δΠint/ext (u s , δu s ), δΠint/ext

for the virtual internal and external work can now be reformulated and discretized with standard FE basis functions. We recall that for this term there is no coupling between the master and the slave side except if one considers self contact. To solve the non-linear discrete equation due to material non-linearities or non-linear strain calculation, a Newton– Raphson algorithm is used, resulting in a linear system DδΠint/ext · 1u k = f k , u k+1 = u k + 1u k in every Newton iteration with DδΠint/ext = DδΠint/ext (u k ) also known as the tangent matrix K . For details, we refer to the monographs [5–7]. Term for contact forces. The virtual contact work δΠc := (δΠcs (tcs , δu s ), δΠcs (tcm , δu m )) can be reformulated with the relation tcs (t, X ) = −tcm (t, Pt (X )) and the introduction of the dual Lagrange multiplier λ as the negative contact traction on the slave side λ = −tcs :  δΠC = λ [δu] d S = ⟨[δu] , λ⟩. Γcs

472

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

Hence the virtual contact work couples the slave body with the master body. The space for the dual Lagrange multiplier λ considering friction is a closed convex cone of M, namely M(λ) := {ν ∈ M : ⟨η, νn ⟩ ≤ ⟨∥ητ ∥, µλn ⟩η ∈ W with ηn ≤ 0}, see [21]. To get to the discretized term for the contact forces, the test functions for the displacements δu  α are discretized with standard basis functions δu αh = np=1 δu αp · φ p and the Lagrange multipliers λ with dual basis  s functions λh = np=1 λ p · ψ p . Separating the nodes into slave nodes S, master nodes M and all other nodes N and using the definition of the well-known Mortar matrices   Dd [i, j] = ψi · φ j dS Idd = δi j φi dS Idd i, j ∈ S, Γcs

 Bd [i, j] = −

Γcs

Γcs

ψi · (φ j ◦ Pt )dS Idd

(3) i ∈ S, j ∈ M,

the discrete version of the virtual contact work can be stated as δΠc = (δu S )T DdT λ + (δu M )T BdT λ. Remark 2. As in our previous paper [32], we assume small deformations on the contact zone in every load increment and use the linearized version of the normal contact law. This assumption also leads to Dd = Dd (u 0 ) and Bd = Bd (u 0 ) with u 0 being the solution of the previous load increment. The same holds for the calculation of the normals and tangents at the slave nodes n = n(u 0 ) and τ = τ (u 0 ). This formulation has the advantage of having to calculate the coupling of the master and slave surface and evaluating the discrete gap only once per increment. For a fully deformation dependent formulation of the contact law we refer to [34,22]. Nevertheless the large deformation formulation for the balance of momentum equation is used to capture the effects of material non-linearities and geometrical non-linearities beyond the contact zone in the numerical calculations. Remark 3. For the calculation of the coupling matrices, a triangulation of the slave surface is needed. Details on the clipping algorithm can be found in [35]. Algebraic representation. In total we solve DδΠint/ext ·1u k + δΠc ·λk+1 = δ f     K

[Bd Dd ]T

in every Newton–Raphson iteration. We omit the subscripts here for better readability. 3. Regularized Coulomb law The standard way of realizing friction in the Mortar framework is to model a perfect Coulomb law. So, each point on the slave surface is either sticking and no relative displacement in tangential direction between the contacting bodies is allowed or it is slipping and the norm of the shear stress equals the friction coefficient µ multiplied by the normal pressure. When one has a closer look at the tangential contact of rough surfaces on the micro scale a different behaviour can be observed, see Fig. 2. Applying first normal pressure, the asperities on the micro scale are pressed against each other, deform elastically and get caught-up with each other. For small shear stresses, these caught up asperities deform elastically in tangential direction and lead to an elastic stick region. With increasing shear stress, parts of the micro contacts are starting to slip while other micro contacts are still sticking. This effect is called micro slip [2,31]. Increasing the shear stress further leads to macro slip, where all micro contacts are slipping. This smooth transition from sticking to slipping can be observed locally on the micro scale but also globally on the macro scale when we consider the contact zone as a whole. In order to model this behaviour numerically, we start in Section 3.1 with a perturbed Lagrange formulation in tangential direction introduced by Simo, Wriggers and Taylor [36] in order to model the elastic stick. In Section 3.2 this approach will be adapted to model the region of micro slip. The latter approach can also be used to model a scale-dependent friction model when one wants to catch the global contact behaviour without modelling the contact zone with fine meshes.

473

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Bolted connection.

(b) Constitutive contact law.

Fig. 2. Example of bolted connection with surface roughness on the micro scale and the resulting frictional constitutive contact law taking the surface roughness into account.

3.1. Coulomb friction with elastic stick We start with the perfect Coulomb law (2) and add a regularization term depending on λτ to the stick branch, see Fig. 3: Φ(λn , λτ ) := ∥λτ ∥ − b ≤ 0 ∧ [u] ˙ τ−

1 λτ λ˙ τ = γ˙ ∧ γ˙ ≥ 0 ∧ γ˙ · Φ = 0, aτ ∥λτ ∥

with aτ being the tangential stiffness also called stick slope and the friction bound b = µ · λn . u¯ ˙ Using a backward Euler method for the displacements u˙ = u− dt and the Lagrange multiplier λ = being the solution of the last load-increment and dt being the load increment, one gets Φ(λn , λτ ) := ∥λτ ∥ − b ≤ 0

∧

(γ − γ¯ ) ≥ 0

∧

λ−λ¯ dt

¯ with (u, ¯ λ)

1 λτ , (λτ − λ¯ τ ) = (γ − γ¯ ) aτ ∥λτ ∥ (γ − γ¯ ) · Φ = 0. ([u]τ − [u] ¯ τ) −

It has to be noted, that the term 1/dt can be neglected here due to the linear relation in the stick branch. 3.1.1. Weak formulation and discretization The weak formulation for the regularized contact law is now derived from the quasi-variational inequality formulated by Klarbring [37]: Find λ ∈ M such that ⟨[u]n − G cn (λn ), νn − λn ⟩ + ⟨([u]τ − [u] ¯ τ ) − G cτ (λτ − λ¯ τ ), ντ − λτ ⟩ ≤ ⟨gn , νn − λn ⟩,

ν ∈ M(λ).

Here the linearized and regularized version of the normal contact condition is used and gn defines the gap in normal direction evaluated at the start of the increment. As before, the displacements are discretized with standard basis functions and the Lagrange multiplier with dual basis functions. The regularization term G cτ (λτ ) = 1/aτ · λτ is discretized again with standard basis functions following the same argument as in [32]: Find λh ∈ M h such that ⟨[u]nh − G cn (λnh ), νnh − λnh ⟩ + ⟨([u]τh − [u] ¯ τh ) − G cτ (λτh − λ¯ τh ), ντh − λτh ⟩ ≤ ⟨gn , νnh − λnh ⟩

ν h ∈ M h (λh ).

This leads to point-wise decoupled contact constraints in tangential direction for each slave node p ∈ S:  ∥λ pτ ∥ − b p ≤ 0       1  ∥λ ∥ − b < 0 ⇒    ˜ · λ u ˜ − pτ p pτ  = 0  pτ a τ       1    ˜ ∥λ pτ ∥ − b p = 0 ⇒ u˜ pτ − · λ pτ  >0 aτ

(4)

474

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

Fig. 3. Coulomb with linear regularization.

with the friction bound b p = µλ pn , defining u˜ := u− ˆ u, ¯ λ˜ := λ−λ¯ and using the well known basis transformation u → −1 S s m β β {M, }. u, ˆ uˆ = u + Dd Bd u , uˆ = u β ∈ N Eq. (4) can be reformulated as a projection operator similar to [38,39]:   λ pτ + c · u˜ pτ − a1τ λ˜ pτ     . λ pτ = b p (5)   max b p , λ pτ + c · u˜ pτ − a1τ λ˜ pτ  3.1.2. Semi-smooth Newton for regularized Coulomb Starting with the projection operator (5), one can easily construct a non-linear complementary function and derive a semi-smooth Newton similar to that in [40]. We obtain for the Coulomb law with elastic stick:          1 1 PL   ˜ ˜ Cτ (λ p , uˆ p ) = max b p , λ pτ + c · u˜ pτ − λ pτ  · λ pτ − b p · λ pτ + c · u˜ pτ − λ pτ aτ aτ with c

b p = µ · max{0, λ pn + c(uˆ pn − g˜ p − gn pl/mpl (λ pn ))}. Remark 4. Defining b p in the semi-smooth Newton as above, the non-linear complementary function in normal direcc tion is reflected and a better numerical convergence is achieved. With gn pl/mpl we denote the regularization functions for the normal contact defined in [32]. Here we restrict our considerations to a linear or a piece-wise linear regularization of the normal contact, since this covers most of the relevant cases. The derived Newton for an exponential regularization of the normal contact would be slightly different, since we use the linearity of gnc in the derivation. Given the previous iterate (u k , λk ), one obtains the new iterate k+1 k k k k (uˆ k+1 p , λ p ) = (uˆ p , λ p ) + (1uˆ p , 1λ p )

by solving the Newton-equation DCτP L (uˆ kp , λkp )(1uˆ kp , 1λkp ) = −CτP L (uˆ kp , λkp ).

(6)

The two max functions now decompose the slave nodes in each semi-smooth Newton iterations into three disjoint sets:   • Inactive nodes In = p ∈ S : bp ≤ 0 ,      • Sticking nodes Iτ = p ∈ S : λ pτ + c · u˜ pτ − a1τ λ˜ pτ  < b p ,        • Slipping nodes Anτ = p ∈ S : λ pτ + c · u˜ pτ − a1τ λ˜ pτ  ≥ b p ≥ 0 . A straightforward computation of (6) shows that the equality constraints for (1u, ˆ λk+1 ) are as follows:

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

475

• Inactive nodes: k+1 λk+1 pn = 0 ∧ λ pτ = 0.

The choice for λ pτ is not a rigorous consequence of the Newton iteration, but directly derived from CτP L (λ p , uˆ p ) =    λ pτ + c · u˜ pτ − a1τ λ˜ pτ  λ pτ = 0. Another possibility to derive the equality constraint for the inactive nodes is to take the projection operator from [38],        bp 1   , 1 · λ pτ + c · u˜ pτ − λ˜ pτ λ pτ = min  ,   aτ λ pτ + c · u˜ pτ − a1τ λ˜ pτ   and construct an alternative non-linear complementary function. • Sticking nodes 1uˆ kpn − 1gnc (λkpn ) · λk+1 = rk,  pn    1 k+1 1 k 1 k k k ˜ 1uˆ pτ − λ pτ + µ · u˜ pτ − λ pτ λ . bkp · λk+1 pn = − aτ aτ aτ pτ    gr c | p∈S

For details how to derive the normal equality constraint, we refer to [32]. • Slipping nodes: k 1uˆ kpn − 1gnc (λkpn ) · λk+1 pn = r , k k+1 −λk+1 pτ + L p · 1uˆ pτ + µv p · λ pn = h p ,

with   bkp 1 d p = λkpτ + c · u˜ kpτ − λ˜ kpτ , ep =   , d p  aτ   T λkpτ d p Fp = k   , M p = e p · (I d − F p ), b p d p    c Hp = I d − 1 − Mp, L p = cH p−1 M p , a τ   dp 1 1 − L p · λkpτ , v p = H p−1   . hp =  c aτ dp Here d p , v p and h p are two-dimensional vectors, e p is a scalar and F p , M p , H p and L p are 2 × 2 matrices. Since all of these quantities can be built point-wise, the computational cost can be neglected. Remark 5. Within the semi-smooth Newton it has to be assured that the matrix H p is invertible for each (uˆ k , λk ). Therefore we make the same modifications to the Robin system as in [40]. Remark 6. For a low tangent stiffness aτ , c should be set to aτ . In this case it holds H p = I d and h p = 0. Thus the non-linear complementary function in tangential direction looses its non-linearity in λ, since d p reduces to d p = λ¯ pτ + aτ u˜ pτ being dependent only on the previous solution λ¯ pτ . The same is achieved by choosing the simple radial return mapping λ¯ pτ + aτ u˜ pτ  λ pτ = b p  ¯λ pτ + aτ u˜ pτ 

(7)

as the projection operator for the slip case. 3.1.3. Condensation of the Lagrange multipliers Now we exploit the full advantage of the dual Mortar formulation and condense the Lagrange multiplier before solving the linear system using λk+1 = D −T ( f k − K k 1u k ).

(8)

476

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

In total we get: • Sticking nodes: k T k −T k Rst D −T p K p · 1u p + τ p · 1uˆ p = gr c + Rst D p f p

with   u˜ kpτ − a1τ λ˜ kpτ 1 T Rst = τ p − µ · n Tp . aτ bp • Slipping nodes: k k k −T k Rsl D −T f p + h kp p K p · 1u p + L t · 1uˆ p = Rsl D

with L t = L kp τ pT ,

Rsl = (τ pT − µv kp n Tp ).

The condensed part of the normal contact and the full matrix representation can be found in [32] with modifications in the tangential part as noted above. 3.2. Iwan regularization Now the concepts from Section 3.1 are expanded in order to construct a Coulomb law with elastic stick and micro slip (Fig. 4(a)). The main idea is that the Coulomb law with a linear regularization can be regarded point-wise as an elastic-slip element consisting of a spring with spring constant k = aτ and a slider element with sliding bound b p = µ · λn in series. This elastic–slip element is also called Iwan element in literature, see e.g. [33]. The smooth transition from sticking to slipping is then realized by parallel arrangement of n iwan of these elastic–slip elements, see Fig. 4(b), resulting in a serial–parallel Iwan model. In the dual Mortar framework, we therefore introduce a separate local Lagrange multiplier λiτ for  every Iwan element i. The same aτ is used for every spring, but different friction bounds are chosen, satisfying b p = bip . The global Lagrange multiplier is finally assembled by   n iwan i λτ · τ. λ = λn · n + i=1

Performing the weak formulation and the discretization analogue to Section 3.1, but using the sum over all Iwan elements for the tangential part instead of one Iwan element, one ends up with point-wise decoupled contact constraints in tangential direction for every Iwan element i:   ∥λipτ ∥ − bip ≤ 0    i n iwan i ∥λ pτ ∥ − bip < 0 ⇒ ∥u˜ pτ − λ˜ ∥ = 0 aτ τ  n  i iwan ˜ i  ∥λ pτ ∥ − bip = 0 ⇒ ∥u˜ pτ − λ ∥>0 aτ τ  with bip = α i · b p and α i = 1. This is equivalent to using the projection operator λipτ = α i b p

aτ λ¯ ipτ + n iwan · u˜ pτ     aτ max α i b p , λ¯ ipτ + n iwan · u˜ pτ 



(9)

for every Iwan element i. We deliberately use the simpler projection operator depending only linearly on λ for the Iwan elements since the case of aτ = ∞ makes no sense when modelling micro slip. 3.2.1. Semi-smooth Newton for Iwan model Starting with Eq. (9), a non-linear complementary function is constructed for every Iwan Element i:        CτP L ,i (λip , uˆ p ) = max α i b p , λ¯ ipτ + aτ /n iwan · u˜ pτ  · λipτ − α i b p λ¯ ipτ + aτ /n iwan · u˜ pτ c

with b p = µ · max{0, λ pn + c(uˆ pn − g˜ p − gn pl/mpl (λ pn ))}.

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Micro slip model.

477

(b) Iwan model.

Fig. 4. Scale dependent friction model realized by a serial–parallel Iwan model.

As one can easily see, the choice whether a node is active depends on the normal component of the global Lagrange multiplier and the choice whether an active node is sticking or slipping depends on the local tangential part of the Lagrange multiplier for the Iwan element i. The nodes are thus decomposed into three sets,   • Inactive nodes Ini = p ∈ S : b p ≤ 0       • Stick nodes Iτi = p ∈ S : λ¯ ipτ + aτ /n iwan · u˜ pτ  < α i b p       • Slip nodes Ainτ = p ∈ S : λ¯ ipτ + aτ /n iwan · u˜ pτ  ≥ α i bip ≥ 0 for each Iwan element i. Now, the computation of the equality constraints for the Iwan element i is straightforward: • Inactive: λi,k+1 = 0. pτ • Sticking: − λi,k+1 −µ pτ

  k λ˜ i,k pτ − aτ /n iwan · u˜ pτ bp

λk+1 pn +

aτ 1uˆ kpτ = −λi,k pτ . n iwan

(10)

• Slipping: − λi,k+1 + pτ

di i k µα i · λk+1 pn + L p · 1uˆ pτ = 0 ∥d i ∥

(11)

with d i = λ¯ ipτ + F pi =

aτ n iwan

i,T λi,k pτ d , i α b p ∥d i ∥

u˜ pτ ,

eip =

αi b p , ∥d i ∥

M ip = eip (I d − F pi ),

L ip = M ip

aτ . n iwan

3.2.2. Condensation of the Lagrange multiplier Before one can condense the global Lagrange multiplier from the algebraic system, one has to calculate the equality constraints for the global Lagrange multiplier. Having the equality constraints for every λipτ , the constraint for the global Lagrange multiplier is calculated as the sum over all Iwan elements. In total we get:

478

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

• Inactive nodes: k+1 λk+1 pn = 0 ∧ λ pτ = 0.

• Active nodes (sum over all Iwan elements): n T · 1uˆ kp − 1gnc (λkpn )n T · λk+1 = rk, p −R pτ λk+1 + L pτ · 1uˆ kp = h pτ p with L pτ =

n slip  i=1

i=n slip

R pτ = τ T −

n slip  i=1

N 

h pτ = −

n iwan

L ip τ T +

aτ T τ n iwan +1

n iwan λ˜ ipτ − aτ /n iwan · u˜ kpτ T di i T µα n + n µ ∥d i ∥ bp i=n +1 slip

λi,k pτ

i=n slip +1

with n slip being the number of Iwan elements per node which are slipping. Remark 7. For the global Lagrange multiplier, the active nodes are no longer separated into sticking and slipping, because each node can be in n iwan + 1 different states, l · 100/n iwan % slipping with 0 ≤ l ≤ n iwan . Condensating the global Lagrange multiplier, one ends up with k k k −T k R pτ D −T f p + h pτ p K p · 1u p + L pτ · 1uˆ p = R pτ D

for all active nodes. The condensed normal part can again be found in [32]. Remark 8. It has to be noted, that the implementation of the Iwan model requires more memory than the approach ¯ ipτ have to be stored for each Iwan element i apart from the global with only one elastic–slip element since λi,k pτ and λ Lagrange multiplier λkp for each node p ∈ S. λkp is updated using Eq. (8) while λi,k p is updated via Eqs. (10) and (11) for each Iwan element i and each active slave node p. But since the Iwan model is typically used to catch the global contact behaviour in case of only a few elements in the contact zone in order to model a scale dependent friction model, these extra costs can be neglected. 4. Numerical examples All numerical results are carried out with the open-source FEM package CalculiX (www.calculix.de). A comparison between the different contact models is made. Moreover, we compare experimental data with simulation results to verify our model and algorithms. 4.1. Menq example 4.1.1. Set-up To test the impact of the regularized tangential contact law, we first have a look at a quite simple numerical benchmark example which was studied analytically by Menq [41] in 2D. A steel block of 10.0 × 1.0 × 0.5 mm (L × B × H ) with the material characteristics of steel E = 2 · 105 MPa and ν = 0.0 is pressed against a rigid surface with p = 4 MPa, see Fig. 5. Here the Poisson’s ratio is set to zero to be able to compare the FEM calculations with the analytically obtained results. We assume a friction coefficient of µ = 0.15 and the bar is meshed with 40 linear HEX8 elements over the length of the bar. Then we apply the following steps in a linear-geometrical quasi-static calculation: In the first step the bar is pressed against the surface. Afterwards we start to pull at one side of the bar denoted by A (initial loading), see Fig. 5. In the third step we push at side A in tangential direction (load reversal) and finally we

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

479

Fig. 5. Menq example, set-up.

pull again (load reversal 2). The shear stress over the length of the bar for initial loading can be calculated analytically for aτ < ∞ by  cosh(κ x L ) t 0 ≤ xl < 1 − a max cosh(κ(1 − a)) tτ = (12)  tmax 1 − a ≤ xl ≤ 1, with x L = x/L being the normalized length of the bar and tmax the friction bound. The ratio a describes the portion of the contact zone, which is already slipping and the parameter κ is set to be  aτ · L 2 · b κ := . E·A Here we compare a perfect Coulomb law (aτ = ∞) with a linear regularized Coulomb law (aτ = 105 ) and an Iwan model with 10 elements per node with the same tangential stiffness. That means we use the tangential stiffness aτ /n iwan for every Iwan element and the friction bounds bip = b p · 2i/(n iwan · (n iwan + 1)). In normal direction we assume hard contact. 4.1.2. Results In Fig. 6, the normalized shear stress over the normalized length of the bar is plotted for different load increments of the initial loading step (from t = 1 to t = 2) and different contact laws: the perfect Coulomb law in Fig. 6(a), the linear regularization in Fig. 6(b) and the Iwan regularization in Fig. 6(c). One can easily see the differences in the transition between sticking and slipping for the different tangential contact laws and how the ratio a of the slipping region is increasing over time. At the end of step 2 (t = 2.0) around 80% of the contact zone is slipping for the perfect Coulomb law, for the linear regularization around 70% and for the Iwan model around 65%. After side A was pulled to a maximal displacement of 2 · 10−4 mm, macro slip is not reached in the whole contact zone. When the tangential displacement in now reversed, one can observe the transition from slipping back to sticking and to slipping again but in the inverse direction, see Fig. 7. Also the doubling of the shear stress function due to the load reversal can be seen nicely at x L = 0.0. For the Iwan model, the smooth transition between sticking and slipping can be seen in Fig. 7(c). Now we have a closer look at the local hysteresis for this simple example. Fig. 8 shows the resulting local hysteresis in one node in the middle of the contact zone. Here, the tangential displacement is plotted versus the normalized shear stress. For all tangential contact laws, the contact conditions are fulfilled and the transition from sticking to slipping is clearly visible. Especially the smooth transition from sticking to slipping in the case of the Iwan regularization can be observed. This behaviour changes if one looks at the global behaviour. In Fig. 9, the total reaction force on side A is plotted versus the displacements. One sees that with fine divisions in the contact zone one can reproduce globally a smooth transition from sticking to slipping even in the case of a locally linear regularized Coulomb law or even a perfect Coulomb law. Differences in the maximum value of the force can be observed because one is still in the region of partial slip at the points of load reversal. The hysteresis for the Iwan regularization exactly matches the solution of a linear regularized Coulomb law, because they are both based on the same tangential stiffness, and thus the transmitted tangential contact force is the same for both regularizations in this simple example. When one reduces the division over the length of the bar to only two elements, the numerical results change significantly. Fig. 10 shows, that there is no smooth transition from sticking to slipping any more in case of a perfect Coulomb law and a linear regularization. Only the Iwan regularization can still reproduce globally a smooth behaviour. But it has to be noted that the tangential stiffness had to be reduced in this case in order to approximate the solution for the fine mesh (aτ = 2 · 104 ). In

480

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Perfect Coulomb.

(b) Linear regularization.

(c) Iwan model with 10 elements. Fig. 6. Menq example, initial tangential loading, normalized shear stress over normalized length of the bar for different frictional contact laws.

summary one can observe that in the case of large scales and consequently coarse meshes in the contact zone a more complex friction model like the Iwan regularization is needed in order to reproduce a physically correct behaviour. 4.2. Resonator 4.2.1. Set-up The second example is an industrial scenario for which measurements exist. We consider a resonator consisting of two discs with a bolted connection under pretension with a M10 bolt in it, see Fig. 11(a). The resonator is excited at one end with an excitation force of 100 N, and the accelerations of the two discs are measured at points A and B, see Fig. 11(b). The relative displacements and the transmitted tangential force are subsequently calculated with a reduced 1D model leading to the measured global hysteresis seen in Fig. 11(c). The tests have been performed by Dominik S¨uß at the Chair of Applied Mechanics at the Friedrich-Alexander University Erlangen-Nuremberg [42]. In order to reproduce the measured behaviour (Fig. 11(c)), the linear H ex8 mesh from Fig. 11(b) is used and the measured relative displacement between the points A and B is applied in a quasi-static calculation as a first approximation. Therefore one has to add quasi-static mounting in order to avoid singularities during the calculation. All parts of the resonator are modelled with the material properties of steel E = 2.1 · 105 MPa and ν = 0.3. A friction coefficient of µ = 0.8 is deduced from the measurements. Then the following steps are applied: • A pretension force Fn = 1188 N is applied in the bolt (Fig. 12(a)). • A tangential displacement of 4.8 · 10−3 mm between A and B is applied in 20 load increments (Fig. 12(b)).

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Perfect Coulomb.

481

(b) Linear regularization.

(c) Iwan model with 10 elements. Fig. 7. Menq example, tangential load reversal, normalized shear stress over normalized length of the bar for different frictional contact laws.

Fig. 8. Menq example, local hysteresis at point in the middle of the contact zone, normalized shear stress over relative tangential displacement.

• The inverse tangential displacement is applied in 40 load increments (Fig. 12(c)). • The tangential displacement is reapplied in 40 load increments (Fig. 12(d)).

482

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

Fig. 9. Menq example, global hysteresis, total tangential contact force over applied displacement at side A, fine mesh.

Fig. 10. Menq example, global hysteresis, total tangential contact force over applied displacement at side A, coarse mesh.

As a naive approach, we model a perfect Coulomb law in tangential direction. Then we compare the solution with our new dual regularized contact formulation, the linear regularization with tangential stiffness aτ = 550 and a Iwan model with 5 elements per node and the same tangential stiffness. Again hard contact in normal direction is assumed. 4.2.2. Results Fig. 13(a) indicates clearly that the perfect Coulomb law does not catch the hysteresis effects properly and fails to reproduce the measurements. From this fact it can be concluded that the set-up is dominated by the micro slip part of the physical tangential contact law (Fig. 2(b)). Therefore we identify a suitable tangential contact stiffness aτ and apply our newly derived dual Mortar formulation. Fig. 13(b) shows, that with a tangential stiffness of aτ = 5.5 · 102 , the hysteresis for the linear regularized Coulomb law and the Iwan model with five elements can approximate the measured behaviour. Of course there are still some differences in the solutions since we perform a quasi-static calculation with quasi-static mounting leading to some unsymmetrical behaviour not observed in the real model (see Fig. 14 which shows the distribution of the normal contact pressure for different steps). Another observation made during the FEM calculation was that the model with linear regularization and the Iwan model needed one more cycle (tangential displacement and inverse tangential displacement) in order to reach the limit hysteresis shown in Fig. 13(b) even in the case of a quasi-static calculation. However, from the numerical point of view the regularization in the tangential contact law leads to a faster convergence as one can see in Table 1. Averaged over all increments of a cycle, the linear regularization needed one iteration

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Experimental set-up.

483

(b) FEM set-up.

(c) Measured global hysteresis. Fig. 11. Resonator example: experimental set-up, FEM mesh and measured global hysteresis, relative displacement between measured points vs. total tangential contact force.

less per increment and the Iwan model needed even 1.4 increments per iteration less than the perfect Coulomb law in order to reach convergence. The maximal number of iterations per increment also reduced significantly using the regularized versions of the tangential contact laws. This better numerical performance can be deduced from the decreased non-linearity in the semi-smooth Newton iteration due to the simpler projection operator (7). In summary this example illustrates that there is a need for this newly derived algorithms being able to deal with constitutive contact laws in tangential direction. 5. Conclusion In this paper, a dual Mortar formulation with a regularization for the tangential contact in the context of rough surfaces or a scale dependent friction model was presented. In terms of these regularizations, the perfect Coulomb law is enriched in order to model elastic stick. The deformation of the contacting asperities on the micro scale is considered in a constitutive contact law formulated in terms of Lagrange multipliers on the contact interface resulting in regularized contact conditions for the tangential behaviour. In order to exploit the full advantages of the duality pairing, a mass lumping in the regularization term is proposed resulting in point-wise decoupled contact constraints. In a second step, the linear regularization used for the elastic stick is combined with an Iwan model in order to model micro slip or alternatively a scale-dependent friction model. With the new dual mortar formulation one can approximate a wide range of constitutive contact laws in tangential direction and, in particular, model a smooth transition

484

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Step 1.

(b) Step 2.

(c) Step 3.

(d) Step 4.

Fig. 12. Resonator example, steps in FEM calculation, step 1 pretension of bolt, step 2 tangential movement, step 3 inverse tangential movement, step 4 tangential movement (displacements (mm) multiplied by 103 for better visualization).

Table 1 Total number of iterations for one cycle, average number of iterations and minimum and maximum number of iterations per increment for resonator example. Reg.

# it per cycle

∅ # it per inc

min # it

max # it

Perfect Coulomb Linear reg. Iwan with 5 el.

363 284 253

4.53 3.55 3.16

2 2 2

13 6 5

from sticking to slipping in a variationally consistent way without increasing the system size due to the duality of the Lagrange multiplier. Combining this formulation with a semi-smooth Newton method having one iteration loop for all system non-linearities, one obtains a robust and powerful tool for industrial contact problems. This was shown by selected numerical examples. Example one, the Menq bar, demonstrates the need of a scale-dependent friction model in case of coarse meshes in the contact zone. The second example, the resonator, demonstrates the need of regularized tangential contact laws in order to approximate measured frictional damping.

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

(a) Perfect Coulomb.

485

(b) Regularization.

Fig. 13. Resonator example: global hysteresis, relative displacement between measured points vs. total tangential contact force.

(a) Step 2.

(b) Step 3. Fig. 14. Resonator example, contact pressure for steps 2 and 3.

Acknowledgements The work in this paper was partly funded by the Federal Republic of Germany, Federal Ministry of Economics and Technology on the basis of a decision by the German Bundestag through the project Lufo 4-4 R&E Turb (grant number 20T1104A) in cooperation with MTU Aero Engines AG. The authors gratefully acknowledge this support. Special thanks go to Dominik S¨uß from the Chair of Applied Mechanics at the Friedrich-Alexander University ErlangenNuremberg for the provided measurements for the resonator example and his support during the FEM-calculations. Special thanks also go to Guido Dhondt from MTU Aero Engines AG for his support concerning the implementation of the presented numerical algorithms within the program package CalculiX. References [1] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1987. [2] F.P. Bowden, D. Tabor, The Friction and Lubrication of Solids, Clarendon Press, Oxford, 2001. [3] N. Kikuchi, J. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, in: SIAM Studies in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1988. [4] C. Eck, J. Januˇsek, M. Krbec, Unilateral Contact Problems: Variational Methods and Existence Theorems, in: Chapman and Hall/CRC Pure and Applied Mathematics Series, Chapman & Hall/CRC Press, London, 2005.

486

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

[5] K. Willner, Kontinuums- und Kontaktmechanik: Synthetische und analytische Darstellung, in: Engineering Online Library, Springer, Berlin, Heidelberg, 2003. [6] T.A. Laursen, Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer, Berlin, Heidelberg, 2002. [7] P. Wriggers, Computational Contact Mechanics, Wiley, New York, 2002. [8] C. Bernardi, Y. Maday, A. Patera, A new nonconforming approach to domain decomposition: the Mortar element method, in: H. Brezis, J. Lions (Eds.), Nonlinear Partial Differential Equations and their Applications, Pitman/Wiley, London, New York, 1994, pp. 13–51. [9] F. Belgacem, The Mortar finite element method with Lagrange multipliers, Numer. Math. 84 (2) (1999) 173–197. [10] F. Belgacem, P. Hild, P. Laborde, The Mortar finite element method for contact problems, Math. Comput. Modelling 28 (4–8) (1998) 263–271. [11] K. Fischer, P. Wriggers, Frictionless 2D contact formulations for finite deformations based on the Mortar method, Comput. Mech. 36 (3) (2005) 226–244. [12] M. Tur, F. Fuenmayor, P. Wriggers, A Mortar-based frictional contact formulation for large deformations using Lagrange multipliers, Comput. Methods Appl. Mech. Engrg. 198 (37–40) (2009) 2860–2873. [13] M. Puso, T. Laursen, J. Solberg, A segment-to-segment Mortar contact method for quadratic elements and large deformations, Comput. Methods Appl. Mech. Engrg. 197 (2008) 555–566. [14] M. Puso, T. Laursen, A Mortar segment-to-segment frictional contact method for large deformations, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4891–4913. [15] T.A. Laursen, J.C. Simo, A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems, Internat. J. Numer. Methods Engrg. 36 (20) (1993) 3451–3485. [16] C. Hesch, P. Betsch, Transient three-dimensional domain decomposition problems: frame-indifferent Mortar constraints and conserving integration, Internat. J. Numer. Methods Engrg. 82 (3) (2010) 329–358. [17] P. Deuflhard, R. Krause, S. Ertel, A contact-stabilized Newmark method for dynamical contact problems, Internat. J. Numer. Methods Engrg. 73 (9) (2008) 1274–1290. [18] B. Wohlmuth, A Mortar finite element method using dual spaces for the Lagrange multiplier, SIAM J. Numer. Anal. 38 (1998) 989–1012. [19] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numer. 20 (2011) 569–734. [20] S. H¨ueber, B. Wohlmuth, A primal–dual active set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg. 194 (2005) 3147–3166. [21] C. Hager, S. H¨ueber, B. Wohlmuth, A stable energy conserving approach for frictional contact problems based on quadrature formulas, Internat. J. Numer. Methods Engrg. 73 (2008) 205–225. [22] M. Gitterle, A. Popp, M. Gee, W. Wall, Finite deformation frictional Mortar contact using a semi-smooth Newton method with consistent linearization, Internat. J. Numer. Methods Engrg. 84 (2010) 543–571. [23] T. Dickopf, R. Krause, Efficient simulation of multi-body contact problems on complex geometries: a flexible decomposition approach using constrained minimization, Internat. J. Numer. Methods Engrg. 77 (2009) 1834–1862. [24] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92 (3) (1991) 353–375. [25] M. Hinterm¨uller, K. Ito, K. Kunisch, The primal–dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13 (3) (2002) 865–888. [26] S. H¨ueber, B. Wohlmuth, Thermomechanical contact problems on non-matching meshes, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1338–1350. [27] I. Temizer, Thermomechanical contact homogenization with random rough surfaces and microscopic contact resistance, Tribol. Int. 44 (2) (2011) 114–124. [28] B. Yang, T.A. Laursen, A Mortar-finite element approach to lubricated contact problems, Comput. Methods Appl. Mech. Engrg. 198 (47–48) (2009) 3656–3669. [29] P. Rabier, J. Martins, J. Oden, L. Campos, Existence and local uniqueness of solutions to contact problems in elasticity with nonlinear friction laws, Internat. J. Engrg. Sci. 24 (11) (1986) 1755–1768. [30] W. Tworzydlo, W. Cecot, J. Oden, C. Yew, Computational micro- and macroscopic models of contact and friction: formulation, approach and applications, Wear 220 (2) (1998) 113–140. [31] K. Willner, Constitutive contact laws in structural dynamics, CMES Comput. Model. Eng. Sci. 48 (2009) 303–336. [32] S. Sitzmann, K. Willner, B.I. Wohlmuth, A dual Lagrange method for contact problems with regularized contact conditions, Internat. J. Numer. Methods Engrg. 99 (3) (2014) 221–238. [33] W.D. Iwan, A distributed-element model for hysteresis and its steady-state dynamic response, J. Appl. Mech. 33 (4) (1966) 893–900. [34] A. Popp, M. Gitterle, M. Gee, W. Wall, A dual Mortar approach for 3D finite deformation contact with consistent linearization, Internat. J. Numer. Methods Engrg. 83 (2010) 1428–1465. [35] M. Puso, A Mortar segment-to-segment contact method for large deformation solid mechanics, Comput. Methods Appl. Mech. Engrg. 193 (6–8) (2004) 601–629. [36] J.C. Simo, P. Wriggers, R.L. Taylor, A perturbed Lagrangian formulation for the finite element solution of contact problems, Comput. Methods Appl. Mech. Engrg. 50 (2) (1985) 163–180. [37] A. Klarbring, J.-S. Pang, Existence of solutions to discrete semicoercive frictional contact problems, SIAM J. Optim. 8 (2) (1998) 414–442. [38] P.W. Christensen, A. Klarbring, J.S. Pang, N. Str¨omberg, Formulation and comparison of algorithms for frictional contact problems, Internat. J. Numer. Methods Engrg. 42 (1) (1998) 145–173. [39] P.W. Christensen, A semi-smooth Newton method for elasto-plastic contact problems, Internat J. Solids Struct. 39 (8) (2002) 2323–2341. [40] S. H¨ueber, G. Stadler, B. Wohlmuth, A primal–dual active set algorithm for three-dimensional contact problems with Coulomb friction, SIAM J. Sci. Comput. 30 (2008) 572–596.

S. Sitzmann et al. / Comput. Methods Appl. Mech. Engrg. 285 (2015) 468–487

487

[41] C.-H. Menq, J. Bielak, J. Griffin, The influence of microslip on vibratory response, part I: a new micro slip model, J. Sound Vib. 107 (2) (1986) 279–293. [42] D. S¨uß, M. Jerschl, K. Willner, Investigation of jointed structures using the multiharmonic balance method, in: Conference Procedings of the Scociety of Experimental Mechanics Series, in: Topics in Nonlinear Dynamics, vol. 1, 2013, pp. 223–228.