A simple weak form for contact problems with coulomb friction

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 4 (2011), pp. 259 – 282 A SIMPLE WEAK FORM FOR CONTACT PROBLEMS WITH COULOMB FRICTION Nguyen Huynh Tan Tai1 and Le Van Anh2 1Faculty of Civil Engineering, Thu Dau Mot University, Vietnam 2Faculty of Science, University of Nantes, France Abstract. This paper proposes a weak form for the contact problem with Coulomb fric- tion, written as extension of the standard virtual work principle and involving both the displacements and the multipliers defined on the reference contact surface. The mixed relationship is shown to be equivalent to the strong form of the initial/boundary value contact problem, and it can be discretized by means of the finite element method in a simple way. Typical numerical examples are given to assess the efficiency of the formu- lation in statics and dynamics. Key words: Contact, friction, finite element method. 1. INTRODUCTION Contact problems are of practical significance in mechanics and their solution is complex because of the nonlinear and non-smooth laws governing the interface response between the contacting bodies. When large deformations and large slips are involved, the problem becomes highly nonlinear and robust formulations must be designed with a view to efficiently solve it. Over more than three decades, a great number of formulations and solution meth- ods have been proposed in the literature in order to deal with different kinds of contact problems encountered, from frictionless contact in linear elastostatics up to frictional con- tact in large deformation with complex material constitutive and interfacial contact laws, either in statics or dynamics. A complete review can be found in the books of Wriggers [1] and Laursen [2], where one can find a great deal of references related to the existing methods for the treatment of the contact constraints, in particular the classical Lagrange multiplier, the penalty and the augmented Lagrangian methods. The augmented Lagrangian methods are the most interesting to date to cope with the contact inequality constraints inasmuch as it combines the regularizing effect of the penalty method and the exact satisfaction of the constraints ensured by the Lagrange multiplier method, without having the ill-conditioning problem inherent to the penalty method. Kinematically linear problems with frictional contact were solved using the aug- mented Lagrangian approach by Alart and Curnier [3], Zavarise, Wriggers and Schrefler [4]. Large deformation contact without friction was also successfully addressed by the same 260 Nguyen Huynh Tan Tai, Le Van Anh method in the early work of Glowinski and Le Tallec [5] considering a plane rigid obstacle, in Heegard and Curnier’s paper [6] considering contact between two elastic bodies, and later in Zavarise and Wriggers [7]. The solution of large deformation frictional contact problems between deformable bodies was accomplished by Simo and Laursen [8, 9] in- cluding the algorithmic symmetrization of the contact tangent stiffness matrix, Wriggers [10], Pietrzak and Curnier [11, 12] and more recently by Feng, Peyraut and Labed [13] for Blatz-Ko hyperelastic bodies, and Mijar and Arora [14]. All previous works considered the Coulomb friction model or an extension of it to incorporate micro-mechanical effects for instance. This paper proposes a simple formulation for the contact problem with Coulomb friction, which is close to Curnier et al’s augmented Lagrangian formulations [11, 12]. First, the basic equations and inequalities governing the strong form of the initial/boundary value contact problem are formulated in Section 2. Then, a new weak form - similar to the virtual work principle - is proposed in Section 3 and it is shown to be equivalent to the strong form of the contact problem. The main characteristics of the weak form is that it contains the inequality constraints specific to contact laws, so that the contact laws are automatically satisfied regardless of the constitutive laws of the contacting bodies and no local integration is needed at the contact surface level. The semi-discretization in space of the weak form by means of the finite element method is conducted in Section 4, which leads one to solve a coupled nonlinear semi - discrete equation system in two unknowns: the nodal displacements and the nodal multipliers. Eventually, some numerical examples in quasi - statics and dynamics assessing the efficiency of the proposed formulation are presented in Section 5. 2. STRONG FORM OF THE CONTACT PROBLEM Consider two bodies whose positions at any time are represented by some regions in a three - dimensional Euclidean affine space. The position at any time of any particle belonging to the bodies is represented by a point in this space. As usually done, a fixed origin being chosen, any point as well as the associated position vector will be identified to a three-component vector in R3. We shall describe the mechanical problem using a Lagrangian description. The reference configuration of the two bodies are represented by the regions Ω (1) o and Ω (2) o in R3. Let us assume that in the time interval [O, T ] of interest, the two bodies may come into contact with each other through the motions denoted φ(1) and φ(2) ∀i ∈ {1, 2}, ∀t ∈ [O, T ], φ(i)(., t) : Ω(i)o 3 X (i) 7→ x(i) = φ(i)(X(i), t) (1) where X(i) and x(i) are the reference and current positions of any particle of body i, respectively. The displacement field in body i is defined as U(i) = x(i) −X(i). The strong form of the problem is expressed by the local equations summarized below, which must be satisfied at any time t ∈ [O, T ] and at any point X(i) ∈ Ω (i) o , i ∈ {1, 2}. The linear momentum balance equation can be written for each body i via divΠ(i) + ρ(i)o f (i) = ρ(i)o U¨ (i) (2) A simple weak form for contact problems with Coulomb friction 261 where Π(i) is the first Piola-Kirchhoff stress tensor, ρ (i) o is the reference density of body i and f(i) the prescribed body force per unit mass. The constitutive relationships relating stress to strain are left unspecified until the numerical examples are considered. This means that the weighted residual relationship proposed hereafter holds for all constitutive materials, as is the case with the standard virtual work principle. The reference boundary S (i) o of body i is partitioned into three non-overlapping and pairwise disjoint parts denoted S (i) oU , S (i) oT and S (i) oc , where S (i) oU and S (i) oT are the parts where the displacement and the tractions are prescribed, respectively, and S (i) oc is the part where contact potentially takes place between the two bodies during some portion of the time interval [O, T ]. The boundary conditions on S (i) oU and S (i) oT on each body i can be summarized as U(i) = U¯ (i) in S (i) oU , Π (i)N(i) = T¯ (i) in S (i) oT (3) where U¯ (i) , T¯ (i) are the prescribed displacement and nominal traction, respectively, N(i) the unit outward normal at point X(i) to surface S (i) o . The boundary conditions on the contact surface will be presented in detail in Section 2.2. The initial conditions for body i are ∀X(i) ∈ Ω(i)o , U (i)(X(i), 0) = 0, U˙ (i) (X(i), 0) = V¯ (i) 0 (X (i)) (4) where V¯ (i) 0 is the prescribed initial velocity field for body i. 2.1. Contact kinematics Let us denote the spatial counterparts of Ω (i) o , S (i) o , S (i) oU , S (i) oT and S (i) oc by Ω(i), S(i), S (i) U , S (i) T and S (i) c , respectively, e.g., S (i) c = φ (i)(S (i) oc , t). In order to describe the relative motion of the two bodies, we arbitrarily choose one body (the target), say body 2, as the reference and to evaluate the proximity of the other (the contactor), thus body 1, with respect to it. It is assumed that the reference surface S (i) oc , i = 1 or 2, is parametrized by the mapping θ(i) and the current surface S (i) c at time t by the mapping ψ (i)(., t)=φ(i)(., t)◦θ(i): A region of R2 → S(i)oc ⊂ R 3 → S(i)c ⊂ R 3 ξ(i) = (ξ(i)1, ξ(i)2) 7→ X(i) = θ(i)(ξ(i)) 7→ x(i) = φ(i)(X(i), t) = ψ(i)(ξ(i), t) (5) Consider a particular point x ∈ S (1) c and define the opposite point y ∈ S (2) c as the closest point to x in the Euclidean sense (Fig. 1): y = arg min x(2)∈S (2) c ‖x− x(2)‖ (6) As may be seen, the notation x(i), i ∈ {1, 2}, in (5) designates any point on S (i) c , whereas the notation x designates a particular point of interest on S (1) c and y the opposite point of x on S (2) c defined by (6). Considering the point X ∈ S (1) oc related to the point x considered by x = φ(1)(X, t), then the point Y ∈ S (2) oc related to point y by y = φ (2)(Y, t) 262 Nguyen Huynh Tan Tai, Le Van Anh




              
 




 
 
 
 
Fig. 1. Contact kinematics satisfies Y = arg min X (2) ∈S (2) oc ‖φ(1)(X, t)−φ(2)(X(2), t)‖ (7) A point Y thus being associated to point X, the proximity g is defined by g = −ν(x− y) = −ν(φ(1)(X, t)− φ(2)(Y, t)) (8) where ν is the normal at point y. The sign convention chosen in (8) implies that the prox- imity g is positive when there is interpenetration and negative when there is a gap between the two bodies. Moreover, the inverse image of X (resp. Y) under the parametrization θ(1) (resp. θ(2)) will be denoted ξ (resp. η). By making x(2) = ψ(2)(ξ(2), t) in (6), one finds that η is characterized by η = argmin ξ(2) ‖x− ψ(2)(ξ(2), t)‖ (9) For notational conveniences, the dependencies on φ(1), φ(2) and ψ(2) will be omitted in the material kinematics variables and only the dependency on (X, t) will be shown. Thus, for instance, one writes η = η(X, t) and g = g(X, t). Assuming that all the mappings introduced are smooth enough, let us define the natural basis (a1, a2) at point y ∈ S (2) c as a1 = ∂ψ(2) ∂ξ(2)1 (η(X, t), t), a2 = ∂ψ(2) ∂ξ(2)2 (η(X, t), t) (10) The surface S (2) c is parametrized conveniently so that the normal n = a1×a2/‖a1× a2‖ at y to S (2) c is directed toward the contactor. Use will also be made of the dual basis (a1, a2) lying in the tangent plane of the current target surface as follows ∀α ∈ {1, 2}, aα = aαβaβ (11) where implicit summation from 1 to 2 is assumed over repeated Greek indices and the 2 × 2 matrix [aαβ] is the inverse of matrix [aαβ] = [aα.aβ]. Eventually, the contact laws A simple weak form for contact problems with Coulomb friction 263 with friction involve the following objective relative velocity at point X and time t [15, 16] V = V(1)(X, t)−V(2)(Y(X, t), t) + gν˙ = −g˙ν + η˙αaα (12) where V(1)(X, t) and V(2)(Y, t) are the velocities of the particles X and Y(X, t), respec- tively. The rates g˙ and η˙α are given by [1, 2]: g˙ = − [ V(1)(X, t)−V(2) (Y(X, t), t) ] ν η˙α = a¯αβaβ [V (1)(X, t)−V(2)(Y(X, t), t)]− ga¯αβνV (2) ,β (η, t) (13) where matrix [a¯αβ] is the inverse of matrix [a¯αβ], with a¯αβ = aαβ + gκαβ, καβ = νψ, (2) αβ is the curvature of the surface ψ(2) at η. On the basis of (12) one can define the slip velocity or the relative tangential velocity VT as the projection of V into the tangent plane at y to S (2) oc : VT = (I− ν ⊗ ν).V = η˙ αaα (14) It should be noted that the slip rate VT associated to the material point X ∈ S (1) oc is resolved in terms of the spatial basis (a1, a2) at point y ∈ S (2) c . As done with other kinematical quantities, use will be made of the abbreviated notations V = V(X, t) and VT = VT (X, t), omitting the dependencies on motions φ (1), φ(2). 2.2. Contact laws Let T(X, t) = Π(1).N be the nominal Piola-Kirchhoff traction vector at point X ∈ S (1) oc (N is the normal vector at X to S (1) oc ). It is resolved in terms of the spatial basis at y ∈ S (2) c as follows T(X, t) = TN(X, t)ν −TT (X, t) = TN(X, t)ν − TTαa α (15) The normal contact law reads ∀t ∈ [O, T ], ∀X ∈ S (1) oc , g(X, t) ≤ 0 { if g(X, t) < 0, then TN(X, t) = 0 if g(X, t) = 0, then TN(X, t) ≥ 0 (16) The tangential contact response is modeled by the Coulomb friction law: ∀t ∈ [O, T ], ∀X ∈ S(1)oc , (a) g(X, t) < 0⇒ TT (X, t) = 0 (no relationship for VT (X, t)) (b) g(X, t) = 0⇒ ‖TT (X, t)‖ ≤ µTN(X, t) where { if ‖TT ‖ < µTN , then VT = 0 (stick) if ‖TT ‖ = µTN , then VT ×TT = 0,VT .TT ≥ 0 (slip) (17) where µ > 0 is the coefficient of friction. It can be verified that the contact laws (16)-(17) are objective. 264 Nguyen Huynh Tan Tai, Le Van Anh 3. WEAK FORM OF THE CONTACT PROBLEM In this section, the strong form of the contact problem - described by equations (1)-(2), the constitutive law, the boundary conditions (3), the contact laws (16)-(17) and the initial conditions (4) - will be converted into an equivalent weak form which is more convenient for the finite element implementation. Our approach is inspired by the augmented Lagrangian formulation of Pietrzak and Curnier [12] and results in a weighted residual relationship which involves the displacement fields U(i), i ∈ {1, 2}, defined in Ω (i) o , and the multiplier field λN(X) and λT (X) = λTαa α defined on S (1) oc . Accordingly, the weighting functions are the virtual displacements U (1)∗, U(2)∗, and the virtual multipliers λ∗N , λ ∗ T = λ ∗ Tαa α. One needs the following preliminary result: Proposition 1. Two positive constants N , T being chosen, the contact laws (16)- (17) are equivalent to the following equalities on S (1) oc : λN = 〈λN + N g〉 λT = ( 1− 〈1− µ〈λN + Ng〉 ‖λT + TVT‖ 〉 ) (λT + TVT ) (18) where 〈.〉 is the Macauley bracket: 〈a〉 = a if a ≥ 0, = 0 if a < 0. The proof is given in [17, 18]. Using the equalities (18) then enables one to formulate the weak form for the contact problem as follows: Proposition 2. (i) The solution fields of the strong problem – namely U(1), U(2) in Ω (1) o ,Ω (2) o and TN , TT , VT on S (1) oc - satisfy the following relationship, provided one makes λN = TN and λT = TT : ∀t ∈ [O, T ], ∀U (1)∗, ∀U(2)∗, ∀λ∗N , ∀λ ∗ T 1, ∀λ ∗ T 2, 2∑ i=1 { − ∫ Ω (i) o Π (i)T : ∇ X (i)U (i)∗dΩo + ∫ Ω (i) o ρ(i)o f (i) U (i)∗dΩo + ∫ S (i) oU ∪S (i) oT T (i) U (i)∗dSo } + ∫ S (1) oc [ 〈λˆN〉ν − ( 1− 〈1 − µ〈λˆN 〉 ‖λˆT ‖ 〉 ) (λˆT ) ]( U (1)∗(X)−U(2)∗(Y(X)) ) dSo + ∫ S (1) oc { (λN − 〈λˆN 〉) λ∗ N N + [ λT − ( 1− 〈1− µ〈λˆN 〉 ‖λˆT ‖ 〉 ) λˆT ] λ∗ T T } dSo = 2∑ i=1 ∫ Ω (i) o ρ(i)o U¨ (i) U (i)∗dΩo (19) where λˆN ≡ λN + Ng, λˆT ≡ λT + TVT , ∇X(i)U (i)∗ is the gradient tensor of U(i)∗ with respect to variables X(i), and 〈1 − µ〈λˆN〉 ‖λˆT ‖ 〉 must be replaced by 0 at any point on S (1) oc where λˆT = 0. (ii) Conversely, the contact surfaces S (1) c and S (2) c being parametrized by ξ 7→ x ∈ S (1) c and η 7→ y ∈ S (2) c , it is assumed that there is a bijective (or piecewise bijective) A simple weak form for contact problems with Coulomb friction 265 mapping S (1) c 3 x 7→ y ∈ S (2) c between them. Then, relation (19) implies at any time t ∈ [O, T ] the following local equations: (1) The momentum balance equation (2) for the two bodies 1 and 2. (2) Relation (3b) on the boundary portion S (i) o \S (i) oc = S (i) oT ∪ S (i) oU : Π (i)N(i) = T(i). (3) The following equalities between the components of the nominal traction vectors and the multipliers on the contactor surface S (1) oc ∀X ∈ S(1)oc , T (1) = Π(1)N(1) = λNν − λT ⇔ TN = λN and TT = λT (20) (4) The normal and tangential contact laws (16) and (17). (5) The following relationship which expresses the equilibrium of the traction vectors at the contact interface ∀X ∈ S(1)oc ,T (2)(Y)dS(2)oc =−(λNν − λT ) ‖X,ξ1 ×X,ξ2 ‖ ‖Y,η1 ×Y,η2 ‖ Dξ Dη dS(2)oc =−T (1)(X)dS(1)oc (21) where point Y = Y(X, t) ∈ S (2) oc corresponds to point X through the above-mentioned bijection, dS (1) oc is a differential reference area in S (1) oc and dS (2) oc its counterpart in S (2) oc related to dS (1) oc . 4. FINITE ELEMENT DISCRETIZATION AND SEMI-DISCRETE EQUATION SYSTEM The spatial discretization of the weighted residual relationship (19) by means of the finite element method will be detailed for the terms pertaining to contact in equality (19) only, since the other terms can be discretized in a standard way. Let the contactor surface S (1) oc (resp. the target surface S (2) oc ) be subdivided into isoparametric finite elements e(1) (resp. e(2)) (the superscripts (1) and (2) are used throughout to remind of surfaces S (1) oc and S (2) oc ). One obtains the following semi-discrete coupled equation system [M]{U¨}+ {Ψ(U)} = {Φ}+ {Φcontact(U,Λ)} {RΛ(U,Λ)} = {0} (22) The unknowns in (22) are the nodal displacement vector {U} = { {U(1)} {U(2)} } in the two bodies 1 and 2, and the multiplier vector {Λ} on the reference contact surface S (1) oc . [M] is the mass matrix, {Ψ(U)} the internal force vector resulting from the strain work in the bodies, {Φ} the external force vector excluding contact force. The contact force vector {Φcontact(U,Λ)} and the residual vector {RΛ(U,Λ)} for the multipliers are the assemblies of the element vectors {Φcontact}ej(ξ) and {RΛ}ej(ξ), where j(ξ) is the area jacobian times the weight corresponding to the parent coordinates ξ = (ξ1, ξ2) ∈ R 2 of quadrature point X ∈ e(1), see Fig. 1 j(ξ) = ‖X,ξ1 (ξ)×X,ξ2 (ξ)‖.(weight at location ξ) (23) 266 Nguyen Huynh Tan Tai, Le Van Anh and {Φcontact} e =   [N (1)(ξ)]e (1) T ( 〈λˆN〉{ν} − ( 1− 〈1− µ〈λˆN〉 ‖λˆT ‖ 〉 ) {λˆT } ) −[N (2)(η)]e (2) T ( 〈λˆN 〉{ν} − ( 1− 〈1− µ〈λˆN 〉 ‖λˆT ‖ 〉 ) {λˆT } )   (24) {RΛ} e = [N (1)(ξ)]e (1) T ( 1 N ( 〈λˆN 〉 − λN ) {ν}+ 1 T { λT − ( 1− 〈1− µ〈λˆN 〉 ‖λˆT ‖ 〉 ) λˆT }) (25) Note that 1− 〈 1− µ〈λˆN〉 ‖λˆT ‖ 〉 =   0 if λˆN < 0 (algorithmic gap) µλˆN ‖λˆT‖ if λˆN ≥ 0 and ‖λˆT ‖ ≥ µλˆN (algorithmic contact with slip) 1 if λˆN ≥ 0 and ‖λˆT ‖ < µλˆN (algorithmic contact with stick) (26) In (26), use has been made of the adjective "algorithmic" in order to emphasize the difference with the "real" situation. For instance, the "algorithmic gap" case means λˆN = λN + Ng < 0, in contrast with the real "gap" case which merely means g < 0; and the "algorithmic slip" case means ‖λˆT ‖ ≥ µλˆN , as opposed to the ‘slip’ case which means ‖λT ‖ = µλN . The nonlinear semi-discrete equations (22) will be solved as a whole using the stan- dard full Newton-Raphson method to simultaneously obtain the nodal displacements and the multipliers, as was done in [6, 12, 7]. • The slip velocity VT = η˙αaα defined in (14) and involved in (24)-(25) has to be computed here by means of an adequate temporally discrete expression. For this purpose, let us first note that Propositions 1 and 2 remain valid if VT is replaced with VT∆t, where ∆t is an arbitrary finite time interval. The quantity VT∆t having the dimension of a length will be referred to as the slip, its use instead of the slip velocity is more consistent in the numerical context because the penalty parameters N and T involved in the weak form then have the same dimension of a force by unit volume, i.e. the same SI units of N/m3. Let us now look at a typical iteration of the current time step n and introduce some new notations as summarized in Fig. 2. For the sake of brevity, any quantity computed at the current step n will not be indexed by n, whereas any quantity related to the previous step will be distinguished by subscript n − 1. Consider two times t and tn−1, which are respectively the current time corresponding to the current iteration of the current step n and the time corresponding to the previous step (n−1). Given a particle on the contactor surface S (1) c with reference position X, its position at time tn−1 is xn−1 = φ (1)(X, tn−1) and its current position is x = φ(1)(X, t). We denote the following quantities related to time tn−1: yn−1 the projection of xn−1 on the target surface S (2) c (tn−1) at time tn−1, Yn−1 and ηn−1 the inverse images of yn−1 under the chain of parametrizations (5). A simple weak form for contact problems with Coulomb friction 267




  
  

 
 
 

     
         





    
  
 

Fig. 2. Illustration of the slip VT∆t It should be emphasized that (i) y and yn−1 are not the respective positions at times t and tn−1 of a same particle in body 2 and (ii) the current position of the particle Yn−1 is φ (2)(Yn−1, t), in general distinct from y = φ (2)(Y, t). With these notations in hand, using the chain of parametrizations (5), which gives y = ψ(2)(η, t) and ψ(2)(ηn−1, t) = Φ (2)(Yn−1, t), and applying the backward Euler in- tegration scheme with time step ∆t = t − tn−1 leads to the following expression for the slip VT∆t = y−Φ (2)(Yn−1, t) (27) • The tangent stiffness matrix related to contact, required for the Newton-Raphson solution of the coupled equation system (22), is the square matrix: [Kcontact] = [ [KUU ] [KUΛ] [KΛU ] [KΛΛ] ] (28) obtained by assembling the element contact stiffnesses [Kcontact] e: [Kcontact] e = [ [KUU ] e [KUΛ] e [KΛU ] e [KΛΛ] e ] (29) These matrices are rectangular and given in [19]. 5. NUMERICAL EXAMPLES This section presents the numerical examples which are carried out on the basis of the previous finite element formulation in order to assess the efficiency of the proposed formulation. All numerical results presented here have been obtained by a home made finite element code. 268 Nguyen Huynh Tan Tai, Le Van Anh 5.1. Rolling disk This example involves quasi-static finite displacements and rotations of a disk rolling between two rigid planes. One of the interesting features of this example is that it allows to check the accuracy of the discrete expression (27) for the slip rate VT (or the slip VT∆t), which is the most essential ingredients in the contact kinematics with friction. Consider an elastic cylinder with radius R = 1m compressed between two rigid parallel plates. The problem is solved assuming that the cylinder is made of an isotropic hyperelastic Saint-Venant Kirchhoff material characterized by the strain energy function in terms of the Green strain E w(E) = 1 2 λ(trE)2 + µ¯tr(E2) (30) where the Lamé constants λ, µ¯ are related to Young’s modulus E and Poisson’s ratio ν by the same relations as in linear elasticity λ = Eν/(1+ ν)/(1− 2ν) µ¯ = E/2/(1+ ν) (31) Here, one takes E = 2.1011N/m2 et ν = 0.3. Friction takes place in the contact between the bodies with the friction coefficient µ = 0.1, a value high enough to ensure the rolling without slip. The lower plate is at rest, while the upper is submitted to a prescribed horizontal translation parallel to the x-axis which makes the disk roll, Fig. 3-a. x y -2 -1 0 1 2 3 4 -1 0 1 2 3 4 5 Elastic Cylinder Fixed plate Rigid plate Prescribed y-displacement x-displacement Prescribed (a) R x y -2 -1 0 1 2 3 4 -1 0 1 2 3 4 5 (b) Fig. 3. Reference configuration of the rolling disk problem • Analytical solution. The problem is solved in plane strain. Before prescribing the horizontal movement, the upper plate is slightly lowered in the y-direction in such a way that it comes into contact with the disk without crushing it much. Thus, the deformation of the disk is small and the following analytical solution from rigid bodies mechanics is acceptable: if there is no slip, the displacement in x of the disk center is equal to half that of the upper plate. A simple weak form for contact problems with Coulomb friction 269 For instance, when the upper plate is translated by piR, the disk rotates by a quarter and the disk center moves by piR/2. • Numerical results. x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 Fig. 4. Initial position and positions at increments 20, 40, 60, 80, 100, respectively. The plates are discretized using one Q8 element 270 Nguyen Huynh Tan Tai, Le Van Anh x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 x y -2 0 2 4 6 0 2 4 6 Fig. 5. Initial position and positions at increments 20, 40, 60, 80, 100. The plates are discretized in 12 Q8 elements The disk is chosen as the contactor and the plates as the targets. The numerical computations are carried out using two meshes: (i) in the first one, each rigid plate is discretized in one single eight-node (Q8) element, and the disk in 148 six-node triangle elements and 333 nodes. (ii) in the second mesh, the plates are discretized in more elements, A simple weak form for contact problems with Coulomb friction 271 12 Q8 elements, whereas the mesh of the disk remains the same as before, Fig. 3-b. In order to prevent the rigid body movement of the disk in the y-direction, a small initial penetration is prescribed at the beginning of the computations. The penalty coefficients are N = 1.10 10N/m3 and T = 1.10 9N/m3. The computations are performed in 100 increments, in the last increment the pre- scribed displacement of the upper plate is Ux = 100× 0.0314 = 3.14m. The convergence is achieved in 4 iterations in average. The results obtained with the two meshes are identical to the 8-th digit. This agree- ment is not so obvious since the projection points y and Yn−1 in (27), in the case of the finer mesh, often belong to two distinct elements of the target surface, and yet the discrete expression (27) gives the same result as with the coarse mesh. Fig. 4 and Fig. 5 show the disk positions at increments 20, 40, 60, 80 and 100. As can be predicted, the displacement of the upper plate is found to be twice that of the disk center, as shown in Fig. 6. Table 1 gives the disk rotation values as a function of the increments, which virtually coincide with the theoretical values since the error is about 0.02% after a quarter of a revolution. These results fully validate the discrete expression (27). Increments D is p la ce m en t (m ) 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 Fig. 6. Displacements as a function of increments : •: horizontal displacement of the upper plate; ◦: horizontal displacement of the disk center Table 1. Disk rotation θ Increment Computed θ(rad) Theoretical θ(rad) Error (%) 20 3.141664E− 01 3.141593E− 01 0.002 40 6.281763E− 01 6.283185E− 01 -0.023 60 9.422779E− 01 9.424778E− 01 -0.021 80 1.256446E+ 00 1.256637E+ 00 -0.015 100 1.570471E+ 00 1.570796E+ 00 -0.021 272 Nguyen Huynh Tan Tai, Le Van Anh 5.2. Interference fit problem This example deals with a cylinder which is forced into a tapered hole in a deformable block through application of a prescribed displacement on its top, Fig. 7. This example is interesting in that both bodies are deformable and may undergo significant plastic strains. The geometry is the same as in [2]. In the reference configuration, the cylinder has a radius of 4cm and a length of 30cm. The outer surface of the block is cylindrical, 30cm in radius and 50cm in length. The tapered hole is the union of a truncated cone (large radius = 4cm, small radius = 3.5cm, length = 15cm) and a cylindrical part (radius = 3.5cm), so that the wall angle of the conical part is 1.9◦. 50 cm 30 cm 30 cm 3.5 cm 4 cm 15 cm Fig. 7. Interference fit problem Both bodies are made of an elastoplastic material with a linear isotropic hardening. The elastoplastic formulation and the material data are the same as in [11] : Young’s modulus E = 100000MPa, Poisson’s ratio ν = 0.3, yield stress σ0 = 300MPa and hard- ening modulus H = 3560MPa. In the numerical computations, the cylinder is taken as the contactor and the block as the target. The problem is modelled in axisymmetry. The whole mesh contains 120 8-node quadrilaterals (Q8) and 430 nodes. The 10 elements on the contact surface which is the outer surface of the cylinder are 3-node line elements. Roller boundary conditions are assumed on the outer and lower surfaces of the block. The numerical quadrature is made using 3×3 Gauss points in Q8 elements, 7 Gauss points in the line elements belonging to the contact surface. Computations are performed without or with friction (µ = 0.2). The penalty parameters are N = 10 11, T = 5.10 10. The axial prescribed displacement is controlled via the arc length method, with the arc length regularly incremented by factor 1.05, until the top of the cylinder is moved through a distance equal to the original cylinder length. The convergence tolerances are A simple weak form for contact problems with Coulomb friction 273 10−6 for the force residual {R}, the arc-length residual as well as the multiplier residual {RΛ}. Without friction, the computation is achieved in 58 steps and the convergence re- quires 4-5 iterations in average. With friction, the computation is achieved in 301 steps and the convergence requires 5-6 iterations in average. Fig. 8 shows the deformed configu- (a) No friction (b) Friction, µ=0.2 Fig. 8. Deformed configurations for the interference fit problem: (a) No friction, (b) Friction µ = 0.2 rations at the last step without or with friction, Fig. 9 the contours of the effective plastic strain and Fig. 10 the contact stresses exerted on the cylinder. 0.21 0.19 0.17 0.15 0.13 0.11 0.08 0.06 0.04 0.02 Effective plastic strain (a) 0.42 0.38 0.34 0.30 0.25 0.21 0.17 0.13 0.08 0.04 (b) Effective plastic strain Fig. 9. Effective plastic strain in the interference fit problem: (a) No friction, (b) Friction µ = 0.2 As expected, without friction the cylinder penetrates further into the block than with friction. Fig. 8b shows that the friction induces larger strains at the upper base of the cylinder as well as an overflow of the block material surrounding the hole. In the 274 Nguyen Huynh Tan Tai, Le Van Anh (b) (b) Fig. 10. Contact stresses on the cylinder for the interference fit problem: (a) No friction, (b) Friction µ = 0.2 frictionless case, the effective plastic strain in the cylinder reaches its highest value (21%) at the top, whereas in the frictional case its maximal value is twofold (42%) near the end of the cylinder and the spread of the plastic region in the block is more significant. The effective plastic strain in the block is smaller than in the cylinder and concentrated around the contact interface. In accordance with the extent of the plasticity observed, Fig. 10 shows that the contact stresses are more significant in the presence of friction. It is found that the cylinder volume remains almost unchanged (maximal change at the last step less than 0.7%, either with or without friction). This may be due to the fact that the wall angle is small (1.9◦) and the block itself is deformable, with the result that the contact stresses are relieved. Fig. 11 plots the thrust force F versus the displacementW of the top of the cylinder, |W| (m) F (N ) 0 0.1 0.2 0.3 0 5E+06 1E+07 1.5E+07 2E+07 2.5E+07 Friction, µ=0.2 No friction Fig. 11. Axial force F versus top displacement W for the interference fit problem. Not all the steps are shown for clarity A simple weak form for contact problems with Coulomb friction 275 in frictional and frictionless cases. Given a displacement of the cylinder top, the thrust force is always higher in frictional than in frictionless case. When the cylinder is completely pushed into the hole, force F in the frictional case is more than 20 times as big as in the frictionless case, this ratio agrees well with [2]. Without friction, force F is practically constant after the prescribed displacement reaches about half the reference length of the cylinder (|W | ≈ 15cm), this result is in good agreement with those reported in [11, 2]. 5.3. Disk projected into a wedge Consider the dynamical interaction between an elastic disk and a wedge bounded by two rigid walls, Fig. 12. This example is based on Wriggers et al’s study [20] related to a 2D ring and on Magnain’s [21] related to a 2D disk. Here, the disk is modeled in 3D, with radius 0.1m and thickness 0.02m. The rigid obstacle is defined in plane xz by the four vertex co-ordinates (in m): A(0.05,-0.305), B(0.15,-0.305), C(0.15,0.095), D(0.13,0.095). With these dimensions, the half-angle of the wedge is 90◦ − 78.32◦ = 11.68◦. Y X Z V 0 (a) Rigid wedge Y X Z 78.32° (b) 0.1m R=0 .1m V 0 A B D C Fig. 12. Impact of a disk on a wedge: a) 3D view, b) front view The disk is made of a Saint-Venant Kirchhoff material with E = 5.107N/m2, ν = 0.4 and ρ0 = 700kg/m 3. Before the impact, it is projected in uniform translation along the opposite direction to z at the velocity V0 = 30m/s. The coefficient of friction between the disk and the walls is either µ = 0 or µ = 0.4. The disk is chosen as the contactor and the walls as the target. The penalty parameters are N = 10 11N/m3 and T = 10 9N/m3. There is no gravity. In order to check the precision of the numerical computations, the disk and both walls are discretized without taking into account the symmetry with respect to the central plane, and the finite element program should provide itself symmetrical results. The disk contains 88 twenty-node brick elements and 8 fifteen-node prismatic elements. Each wall is discretized with only one brick element. The time stepping is carried out using the Moreau scheme with parameters θ = ξ = 1 2 , time step ∆t = 10−4s, and the whole computation takes 130 steps. The convergence necessitates 3-5 iterations in average without friction and more iterations - 4-7 in average - in the case of friction. 276 Nguyen Huynh Tan Tai, Le Van Anh Y X Z t=0.000s Y X Z t=0.001s Y X Z t=0.002s Y X Z t=0.003s Y X Z t=0.004s Y X Z t=0.005s Y X Z t=0.006s Y X Z t=0.007s Y X Z t=0.008s Y X Z t=0.009s Y X Z t=0.010s Y X Z t=0.011s Y X Z t=0.012s Y X Z t=0.013s Fig. 13. Impact of a disk on a wedge. Deformed configurations through the process - Frictionless case A simple weak form for contact problems with Coulomb friction 277 • Verification of the inequality g ≤ 0 The numerical results show that the inequality g ≤ 0 is satisfied with high precision at nodes in contact and at any time: in the increments where contact occurs, the maximum proximity g over all the contacting nodes is found to be about 10−4m, i.e. 0.1% of the disk radius. • Movement of the disk Fig. 13 displays the disk configurations in the frictionless case. The contact takes place during about 0.09s, between times t = 0.0016s and t = 0.0106s. The disk plunges into the wedge and rebounds back, leaving the wall for good after time t > 0.016. Fig. 14 and Fig. 15 representing the position of the disk center as a function of time t (s) U z (m ) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 -0.1 -0.05 0 0.05 no friction friction, µ=0.4 Fig. 14. Impact of a disk on a wedge. The disk center displace- ment versus time t (s) V z (m /s ) 0.002 0.004 0.006 0.008 0.01 0.012 0.014 -30 -20 -10 0 10 20 30 40 Fig. 15. The disk center velocity versus time in the frictionless case shows that the disk center stops at t = 0.006s approximately (the middle of the contact interval), and then goes back at the average velocity of 30m/s, which is precisely opposite to the incident velocity. Comparing Fig. 14 with Fig. 19 hereafter giving the kinetic energy, one finds that the whole disk comes to a standstill at time t = 0.006s when the disk center reaches the lowest level. Fig. 16 displays the disk configurations in the frictional case, µ = 0.4. This time, the disk does not go so far as in the frictionless case and at about time t = 0.0045s when the disk center is at lowest level, the whole disk comes to a standstill. Afterwards, the disk rebounds a little before its vibration is stopped by the strong friction and the disk is finally stuck in the wedge: Fig. 14 and Fig. 20 show that the disk is almost at rest after tc = 0.008s. All these results are in good accordance with those reported in [20] and [21]. A closer examination of the numerical results shows that once the disk is stuck in the wedge, it keeps vibrating slightly. However, the vibrations diminish rapidly due to 278 Nguyen Huynh Tan Tai, Le Van Anh Y X Z t=0.000s Y X Z t=0.001s Y X Z t=0.002s Y X Z t=0.003s Y X Z t=0.004s Y X Z t=0.005s Y X Z t=0.006s Y X Z t=0.007s Y X Z t=0.008s Y X Z t=0.009s Y X Z t=0.010s Y X Z t=0.011s Y X Z t=0.012s Y X Z t=0.013s Fig. 16. Impact of a disk on a wedge. Deformed configurations through the process - Frictional case with µ = 0.4 A simple weak form for contact problems with Coulomb friction 279 small frictional slips at the contact interface. Eventually, one may say that the disk is in equilibrium between the walls. • Verification of the inequalities TN ≥ 0 and ‖TT‖ ≤ µTN The inequality TN ≥ 0 in (16) is fully satisfied, except at some increments and some nodes belonging to those elements overlapping the contact and the non-contact zones. The Coulomb inequality ‖TT ‖ ≤ µTN in (17) is met with high precision (wherever TN ≥ 0). • Contact stresses The contact stress distribution in space and time is rather complex, as showed in Fig. 17 and Fig. 18 for some typical increments. Y X Z Y X Z Y X Z Fig. 17. Impact of a disk on a wedge. Frictionless case. Contact stresses With friction µ = 0.4, the contact stresses present three profiles: (i) In the first stage of the movement when the disk sinks into the wedge, Fig. 18a, there is slip at every contact nodes. The contact stresses are upward, in the opposite direction of the movement. (ii) In the second stage when the disk slightly rebounds backwards, Fig. 18b, there is still slip everywhere but this time the contact stresses are downward. (iii) In the last stage when the disk is stuck in the wedge, Fig. 18c, there is stick at the contacting nodes, the contact stress profile is complex with their sum being zero in y- and z-directions and thus satisfying the global equilibrium of the disk. Y X Z Y X Z Y X Z Fig. 18. Impact of a disk on a wedge. Frictional case with µ = 0.4. Contact stresses • Energies The evolution in time of the system energies (kinetic energy, elastic strain energy and total energy) is given in Fig. 19 for the frictionless case. 280 Nguyen Huynh Tan Tai, Le Van Anh When the contact begins, the kinetic energy gradually drops to zero as the strain energy increases, since the disk slows down and is more and more deformed. In the next stage when the disk rebounds backwards, the kinetic energy is recovered whereas the strain energy decreases. After the contact is over, the disk becomes circular again and the strain energy stored during the impact is released in the form of kinetic energy. During the whole process, the total energy is constant. It can be noted that the strain energy is small yet not zero, which means that there are small oscillations of the disk after the rebound. t (s) E (J ) 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 50 100 150 200 Deformation energy Total energy Kinetic energy Fig. 19. Impact of a disk on a wedge. Energy evolution in the frictionless case t (s) E (J ) 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 50 100 150 200 Deformation energy Total energy Kinetic energy Fig. 20. Impact of a disk on a wedge. Energy evolution with friction µ = 0.4 In the frictional case, the energy evolution is more intricate than in the frictionless case, as can be seen in Fig. 20. The curves present the same trends as in [21], although the hypotheses adopted are different (the problem was solved in [21] in plain strain using the Blatz-Ko hyperelastic constitutive law). The total energy decreases steadily and markedly due to friction. At the end of the computation, there is a strong total energy loss and the remaining energy is stored in the elastic strain energy. 6. CONCLUSIONS The weak form (19) has been designed as an extension of the standard virtual work principle to deal with the contact problem with Coulomb friction. Since it is a genuine weighted residual relationship involving arbitrary (regular) virtual displacements and multipliers, it can be discretized by means of the finite element method in a simple way. The finite element discretization process holds for any isoparametric finite element on the contact surface. As the weak form already incorporates the contact laws, no local integration of the frictional equations is required in the numerical solution procedure. The A simple weak form for contact problems with Coulomb friction 281 numerical examples presented show the ability of the proposed formulation to deal with quasistatic and dynamic contact problems, and many other examples can be found in [19]. ACKNOWLEDGMENTS The authors would like to thank Prof. Nguyen Quoc Son, of Laboratoire de Mecanique des Solides, Ecole Polytechnique in France, for his instructive remarks which have enabled them to refine Proposition 2. REFERENCES [1] P. Wriggers, Computational Contact Mechanics, Wiley, (2002). [2] T. A. Laursen, Computational Contact and Impact Mechanics - Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer, (2003). [3] P. Alart and A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Computer Methods in Applied Mechanics and Engineering, 92, (1991), 353-375. [4] G. Zaravise, P Wriggers, and B. A. Schrefler, On augmented Lagrangian algorithms for ther- momechanical contact problems with friction, International Journal for Numerical Methods in Engineering, 38, (1995), 2929-2949. [5] R. Glowinski and P. Le Tallec, Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics, SIAM Publications, Philadelphia, (1989). [6] J.-H. Heegaard and A. Curnier, An augmented Lagrangian method for discrete large-slip contact problems, International Journal for Numerical Methods in Engineering, 36, (1993), 569-593. [7] G. Zavarise and P. Wriggers, A superlinear convergent augmented Lagrangian procedure for contact problems, Engineering Computations, 16, (1999), 88-119. [8] J. C. Simo and T. A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, Computers and Structures, 42, (1992), 97-116. [9] T. A. Laursen and J. C. Simo, Algorithmic symmetrization of Coulomb frictional problems using augmented Lagrangians, Computer Methods in Applied Mechanics and Engineering, 108, (1993), 133-146. [10] P. Wriggers, Finite element methods for contact problems with friction, Tribology Interna- tional, 29, (1996), 651-658. [11] G. Pietrzak, Continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems, PhD thesis, Ecole Polytechnique Fédérale de Lau- sanne, (1997). [12] G. Pietrzak and A. Curnier, Large deformation frictional contact mechanics: continuum for- mulation and augmented Lagrangian treatment, Comput. Methods Appl. Mech. Engrg., 177, (1999), 351-381. [13] Z. -Q. Feng, F. Peyraut and N. Labed, Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies, International Journal of Engineering Science, 41, (2003), 2213-2225. [14] A. R. Mijar and J. S. Arora, An augmented Lagrangian optimization method for contact anal- ysis problems, 1: formulation and algorithm, Structural and Multidisciplinary Optimization, 28, (2004), 99-112. 282 Nguyen Huynh Tan Tai, Le Van Anh [15] A. Curnier, Qi.-Ch. HE and A. Klarbring, Continuum mechanics modelling of large deforma- tion contact with friction, In M. Raous, M. Jean and J. J Moreau, eds, Contact Mechanics (Plenum, New York), (1995), 145-158. [16] A. Klarbring, Large displacement frictional contact: a continuum framework for finite element discretisation, European Jounal of Mechanics - A/Solids, 14, (1995), 237-253. [17] A. Le Van and Tai H. T. Nguyen, Une formulation variationnelle du problème de contact avec frottement de Coulomb, Comptes Rendus Mécanique, 336, (2008), 606-611. [18] A. Le Van and Tai H. T. Nguyen, A weighted residual relationship for the contact problem with Coulomb friction, Computers and Structures, 87, (2009), 1580-1601. [19] H. T. T. Nguyen, Contribution à l’étude du contact avec frottement de Coulomb en statique et dynamique des structures, PhD thesis, Université de Nantes, France, (2009). [20] P. Wriggers, T. Vu Van and E. Stein, Finite element formulation of large deformation impact- contact problems with friction, Computers and Structures, 37, (1990), 319-331. [21] B. Magnain, Développement d’algorithmes et d’un code de calcul pour l’étude des problèmes de l’impact et du choc, PhD thesis, Université Evry-Val d’Essonne, (2006). Received July 1, 2011