A weighted residual relationship for the contact problem with Coulomb friction

In this work, a weighted residual relationship is proposed as an extension of the standard virtual work principle to deal with the large deformation contact problem with Coulomb friction. This weak form is a mixed relationship involving the displacements and the multipliers defined on the reference contact surface of the contactor and is shown to be equivalent to the strong form of the initial/boundary value contact problem. The discretization in space by means of the finite element method is carried out on the mixed relationship in a simple way in order to obtain the semi-discrete equation system. The contact tangent stiffness is derived and numerical examples are presented to assess the efficiency of the formulation.


Introduction
Problems where contact takes place between two or more bodies are widespread and of practical significance in mechanics. Their solution is complex because of the nonlinear and nonsmooth laws governing the interface response between the contacting bodies. In the context of large deformations and large slips, the problem becomes even more 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 methods 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 contact in large deformation with complex material constitutive and interfacial contact laws, either in statics or dynamics. The interested reader is referred to the book of Kikuchi and Oden [1] and the more recent books of Wriggers [2] and Laursen [3] for a review of the numerous 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. As the formulation proposed in the present work is inspired from methods of the augmented Lagrangian type, we shall limit ourselves to a brief review of these methods in the literature.
The augmented Lagrangian method is often chosen 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. It was applied to different contact problems in many vari-ants. Kinematically linear problems with frictional contact were solved using the augmented Lagrangian approach by Alart and Curnier [4], Zavarise et al. [5]. Large deformation contact without friction was also successfully addressed by the same method in the early work of Glowinski and Le Tallec [6] considering a plane rigid obstacle, in Heegard and Curnier's paper [7] considering contact between two elastic bodies, and later in Zavarise and Wriggers [8]. The solution of large deformation frictional contact problems between deformable bodies was accomplished by Simo and Laursen [9,10] including the algorithmic symmetrization of the contact tangent stiffness matrix, Wriggers [11], Pietrzak and Curnier [12,13] and more recently by Feng et al. [14] for Blatz-Ko hyperelastic bodies, and Mijar and Arora [15,16]. All previous works considered the Coulomb friction model or an extension of it to incorporate micro-mechanical effects for instance.
In most of the formulations, the tangential contact stresses were computed by locally integrating the contact laws and the contact stiffness matrix is consistently derived from these stresses, following the spirit of the return mapping scheme in elastoplasticity as described by Giannakopoulos in [17]. In some formulations, e.g., [5], the tangential contact law adopted was not of the rate type and the above-mentioned local integration was not necessary. Curnier and co-workers [4,7,13] chose another approach by directly incorporating the contact laws into their augmented Lagrangians, which enabled them to obtain the correct tangential contact stresses without resorting to the local integration.
As regards the global solution of the discrete equation system, Uzawa's method was usually chosen because of its simplicity: the multipliers are updated in an outer augmentation loop and the displacements computed in a standard way inside an inner loop while keeping the multipliers fixed. Otherwise, the discrete equation system can be solved simultaneously for the displacements and the multipliers, as was done in [13,8], for instance. One can also consult several different global algorithms for solving augmented Lagrangian problems in [6, Chapters 3 and 7].

Motivation
The purpose of this paper is to propose a simple formulation for the contact problem with Coulomb friction which combines the benefits of the aforementioned formulations of the augmented Lagrangian type. The sought formulation is based on a new weak form -similar to the virtual work principle -which should fulfill the following conditions: (i) The weak form should incorporate the contact laws. More precisely, although it is an equality, it must imply the inequality constraints specific to contact laws. The consequence of this is 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. (ii) The weak form may look like an augmented Lagrangian principle in that it involves both the displacements and the multipliers, but, more importantly, it should be stated in the form of a weighted residual relationship, with arbitrary (regular) virtual displacements and multipliers (an overview of the weighted residual methods can be found in [18]). As such, the weak form can be discretized by means of the finite element method in a straightforward way.

Outline of the paper
The paper is organized as follows. The strong form of the initial/ boundary value contact problem is summarized in Section 2. In Section 3, a weighted residual relationship is proposed as an extension of the standard virtual work principle with the aim of encompassing the large deformation contact response with Coulomb friction. This mixed weak form which involves both the displacements and the multipliers defined on the reference contact surface of the contactor will be shown to be equivalent to the strong form of the contact problem. Section 4 deals with the semi-discretization in space of the weak form by means of the finite element method, which results in a coupled nonlinear semi-discrete equation system in two unknowns: the nodal displacements and the nodal multipliers. In Section 5, the contact tangent stiffness matrix which is necessary to the Newton-Raphson iterative solution is derived by the exact linearization of the weak form. Different numerical examples assessing the efficiency of the proposed formulation are presented in Section 6.

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 R 3 .
We shall describe the mechanical problem using a Lagrangian description. The reference configuration of the two bodies are rep- 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 2 ½O; T where V ðiÞ 0 is the prescribed initial velocity field for body i.

Contact kinematics
Let us now focus on the equations governing the response on the contact interface, which are the essential ingredients in the large deformation contact problem. Let the spatial counterparts of X ðiÞ Consider a particular point x 2 S ð1Þ c and define the opposite point y 2 S ð2Þ c as the closest point to x in the Euclidean sense ( Fig. 1): In writing (6), we assume that point y is unique without discussing the case of non-uniqueness, which is not of major consequence in numerical computations. The interested readers may consult comprehensive comments on this issue in the literature, e.g., [ is smooth at y, then y is the projection of x on S ð2Þ c . As may be seen, the notation x ðiÞ ; i 2 f1; 2g, 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).
The point x 2 S ð1Þ c being associated to (at least) one point y 2 S ð2Þ c , the proximity g is defined by where m is the normal at point y. The sign convention chosen in (7) implies that the proximity g is positive when there is interpenetration and negative when there is a gap between the two bodies. Thus, symbol g is the opposite of what is referred to as the gap in the literature. It can be checked that g is an objective quantity.
Let us now convert all the spatial quantities defined above into material ones defined on the reference configurations, by means of the motions / ð1Þ and / ð2Þ in (1). Considering the point X 2 S ð1Þ oc related to the point x considered by x ¼ / ð1Þ ðX; tÞ, then the point Y 2 S ð2Þ oc related to point y by y ¼ / ð2Þ ðY; tÞ satisfies: A point Y thus being associated to point X, the proximity g defined in (7) can be rewritten as Eventually, the inverse image of X (resp. Y) under the parametrization h ð1Þ (resp. h ð2Þ ) will be denoted n (resp. g). By making x ð2Þ ¼ w ð2Þ ðn ð2Þ ; tÞ in (6), one finds that g is characterized by For notational conveniences, the dependencies on / ð1Þ , / ð2Þ and w ð2Þ will be omitted in the material kinematics variables and only the dependency on ðX; tÞ will be shown. Thus, for instance, one writes g ¼ gðX; tÞ and g ¼ gðX; tÞ.
Assuming that all the mappings introduced are smooth enough, let us define the natural basis ða 1 ; a 2 Þ at point y 2 S ð2Þ c as a 1 ¼ ow ð2Þ on ð2Þ2 gðX; tÞ; t The surface S ð2Þ c is parametrized conveniently so that the normal n ¼ a 1 Â a 2 =ka 1 Â a 2 k at y to S ð2Þ c is directed toward the contactor. Use will also be made of the dual basis ða 1 ; a 2 Þ lying in the tangent plane of the current target surface as follows: where implicit summation from 1 to 2 is assumed over repeated Greek indices and the 2 Â 2 matrix ½a ab is the inverse of matrix The contact kinematics with friction is described here by using the following objective relative velocity at point X and time t [19,20]: where V ð1Þ ðX; tÞ and V ð2Þ ðY; tÞ are the velocities of the particles X and YðX; tÞ, respectively. The rates _ g; _ g a and _ m are given by [2,3]: where matrix ½ a ab is the inverse of matrix ½ a ab , with a ab ¼ a ab þ gj ab ; j ab ¼ mw ð2Þ ;ab is the curvature of the surface w ð2Þ at g, the symbol designates the tensor product defined as ða bÞ Á c ¼ aðb Á cÞ for all vectors a, b and c. Note that tensor j ac a cb a a a b ¼ a ac j cb a a a b in (14c) is symmetric.
On the basis of (13) one can define (i) the relative normal velocity V N as the projection of V on the normal vector m, and (ii) the slip velocity or the relative tangential velocity V T as the projection of V into the tangent plane at y to S ð2Þ oc : It should be noted that the slip rate V T associated to the material point X 2 S ð1Þ oc is resolved in terms of the spatial basis ða 1 ; a 2 Þ at point y 2 S ð2Þ c . As done with other kinematical quantities, use will be made of the abbreviated notations V ¼ VðX; tÞ and V T ¼ V T ðX; tÞ, omitting the dependencies on motions / ð1Þ ; / ð2Þ .

Contact tractions
Let the contact Cauchy traction vector tðx; tÞ at a typical point x 2 S ð1Þ c be resolved in terms of the local basis ða 1 ; a 2 ; mÞ at the projection point y 2 S ð2Þ c as follows: The normal component t N of traction vector t in the direction of normal m is the contact pressure, positive in compression. The tangential component t T ¼ t Ta a a is the opposite of the frictional traction. It can be verified that the pressure t N and the frictional traction t T are objective variables. Let TðX; tÞ ¼ P ð1Þ :N be the nominal Piola-Kirchhoff traction vector at point X 2 S ð1Þ oc (N the normal vector at X to S ð1Þ oc ). Like the contact Cauchy traction, the nominal traction T is resolved in terms of the spatial basis at y 2 S ð2Þ c as follows: From relation tdS ð1Þ ¼ TdS o is a differential reference area, dS ð1Þ its spatial counterpart), one derives the following relationships between the spatial and material traction components:

Contact laws
As a means of describing the interfacial response at current time t, the contact laws are naturally expressed in the spatial form. However, one may use (18) in order to recast them into the mixed form involving the relative material velocity and the nominal contact tractions, which are more suitable for the numerical implementation within the Lagrangian description framework. The normal contact law then reads: In this work, the tangential contact response is modeled by the Coulomb friction law: where l > 0 is the coefficient of friction.
It can be verified that the contact laws (19) and (20) are objective.

Weighted residual relationship for the contact problem
The strong form of the contact problem is described by Eqs. (1) and (2), the constitutive law, the boundary conditions (3), the contact laws (19) and (20) and the initial conditions (4). In this section, we convert the strong form into an equivalent weak one which is more convenient for the finite element implementation.
The approach chosen is inspired by the augmented Lagrangian formulation of Pietrzak and Curnier [13] and amounts to rewrite the augmented Lagrangian formulation of Alart and Curnier [4] in a generic weak form. More specifically, a weighted residual relationship is proposed which involves the displacement fields U ðiÞ ; i 2 f1; 2g, defined in X ðiÞ o , and the multiplier field k N ðXÞ and k T ðXÞ ¼ k Ta a a defined on S ð1Þ oc . Accordingly, the weighting functions are the virtual displacements U ð1ÞÃ , U ð2ÞÃ , and the virtual multipliers Ta a a . Two positive constants N ; T being chosen, the weighted residual relationship is stated in the following proposition [21].

Proposition 1
I The solution fields of the strong problem -namely U ð1Þ ; U ð2Þ in X ð1Þ o ; X ð2Þ o and T N ; T T ; V T on S ð1Þ oc -satisfy the following relationship, provided one makes k N ¼ T N and k T ¼ T T : where r X ðiÞ U ðiÞÃ is the gradient tensor of U ðiÞÃ with respect to variables X ðiÞ ; hÁi is the Macauley bracket: hai ¼ a if a P 0; ¼ 0 if a < 0, and h1 À lhk N þ N gi kk T þ T V T k i must be replaced by 0 at any point oc where k T þ T V T ¼ 0.
I Conversely, the contact surfaces S ð1Þ c and S ð2Þ c being parametrized by n # x 2 S ð1Þ c and g # y 2 S ð2Þ c , it is assumed that there is a bijective (or piecewise bijective) mapping S ð1Þ c 3 x#y 2 S ð2Þ c between them. Then, relation (21) implies at any time t 2 ½O; T the following local equations: (1) The momentum balance Eq. (2) for the two bodies 1 and 2. The following equalities between the components of the nominal traction vectors and the multipliers on the contactor surface S ð1Þ oc .
The normal and tangential contact laws (19) and (20). (5) The following relationship which expresses the equilibrium of the traction vectors at the contact interface.

Remarks
(i) In general, k T þ T V T differs from zero since k T and V T do not vanish simultaneously: if the slip rate vanishes, then the tangential traction does not, and vice versa. However, in some problems with a specific symmetry, there may be some points where k T þ T V T is zero. In practice, the set of these points is either isolated points or lines on the contact surface S ð1Þ oc . (ii) In dynamics, when shocks and impacts are expected, the acceleration € U ðiÞ should be understood as the derivative of the velocity in the sense of the distributions [22].
Proof. The first part of the above proposition, which states that the strong form implies the weak one, can be readily checked. Let us prove the reciprocal statement for an arbitrary fixed time t 2 ½O; T.

Since (21) holds for all virtual fields
The first consequence of (24) is k N P 0. Next, two mutually exclusive cases arise: -either k N þ N g < 0, then (24) entails: -or k N þ N g P 0, then (24) entails: The results for the two cases above can be gathered together as follows: which is the normal contact law (19) provided k N ¼ T N which will be shown later.

Consider (21) again and choose now
The localization theorem can then be applied assuming that the coefficient of k Ã T is a continuous function in space. At any point on S ð1Þ The first consequence of (28) is that the slip velocity V T is parallel to k T . Let us show a more precise result, specifically V T is oriented in the same direction as k T . To do this, let us temporarily As a P 0 and 1 À h1 À ai P 0; 8a, one has the following mutually exclusive cases depending on a's value: Thus, in any case V T is oriented in the direction of k T . Now, let us distinguish two cases: -If k N þ N g < 0, then one gets from (25) and (28) g < 0 and k T ¼ 0.
-If k N þ N g P 0, then one gets from (26) g ¼ 0. Two mutually exclusive sub-cases arise here: By taking the norm of both sides of the last equality, recalling here k N þ N g P 0 and g ¼ 0, one finds kk T k ¼ lk N .The last results can be summarized as follows: which is merely the tangential contact law (20) provided Eventually, let us make in (21) and U ð2ÞÃ arbitrary. According to the divergence theorem, relationship (21) leads to the following equality: where N ðiÞ is the unit outward normal at point X ðiÞ 2 S ðiÞ o . Note that the last integral of the left-hand side of (32) is taken over S ð1Þ oc . It can be recast into an integral over S ð2Þ oc by using the chain of bijections n#X#Y#g which allows one to change the integration variables: By substituting (34) into (32), one gets relations (2), (3b), (23) and Comparing (35) with the definition (17), T ð1Þ ¼ T N m À T T , and making use of (24) and (28) lead to equality (22) between the contact tractions and the multipliers. Next, the normal and tangential contact laws (19) and (20) are derived from (22), (27) and (31). It remains to consider the points on S ð1Þ The subsequent reasoning is analogous yet simpler than the case k T þ T V T -0. It can easily be checked that the tangential contact law (20) as well as (Relations 22) and (23) are satisfied at the points where k T þ T V T ¼ 0. The whole Proposition 1 has been proved. j

Semi-discrete equation system
In this section, the weighted residual relationship (21) will be discretized in space by means of the finite element method. The discretization process will be described in the 3D framework only but the results for 2D problems can be derived in a similar manner. The following shorthand notations will be used for brevity: We shall confine ourselves to the terms pertaining to contact in equality (21), denoted W Ã contact , since the other terms can be discretized in a standard way. Let the contactor surface S ð1Þ oc be subdivided into a number of isoparametric finite elements e ð1Þ (the superscript '(1)' is used throughout to remind of the contactor surface S ð1Þ oc ), then expression W Ã contact can be recast as the sum of the element contact terms ðW Ã contact Þ e ð1Þ related to elements e ð1Þ and evaluated using the Gauss quadrature rule: where all terms are computed at Gauss points and jðnÞ is the area jacobian times the weight corresponding to the parent coordinates n ¼ ðn 1 ; n 2 Þ 2 R 2 of quadrature point X 2 e ð1Þ , see Fig. 1: In the sequel, the following notational convention is adopted for the matrix representations: curly brackets {} designate a column vector, the angle brackets < > an arrow vector, and the square brackets [ ] a general matrix. In any element e ð1Þ & S ð1Þ oc , the reference coordinates X, the displacement U ð1Þ and the multipliers k are interpolated by In (39b) for instance, fU ð1Þ g is the column vector with the three components of the displacement U ð1Þ related to a fixed Cartesian orthonormal basis ðe 1 ; e 2 ; e 3 Þ. The column vector fU ð1Þ g e ð1Þ contains the 3nne ð1Þ nodal displacement components associated with element e ð1Þ , where nne ð1Þ is the number of nodes of element e ð1Þ . The matrix ½N ð1Þ ðnÞ e ð1Þ contains the shape functions of element e ð1Þ and is of dimension 3 Â 3nne ð1Þ .
The normal and tangential components of multiplier k can be derived from (39c) via: where fk T g is a 3-component column vector.
Like the contactor surface, the target surface S ð2Þ oc is subdivided into a number of isoparametric finite elements. Given a quadrature point X in element e ð1Þ , let e ð2Þ the (or an) element of S ð2Þ oc containing the reference point Y corresponding to the projection In element e ð2Þ , the reference coordinates Y and the displacement U ð2Þ are interpolated by where g ¼ ðg 1 ; g 2 Þ 2 R 2 is the parent coordinate of Y and the superscript (2) is used to remind of element e ð2Þ & S ð2Þ oc . The column vector fU ð2Þ g e ð2Þ contains the 3nne ð2Þ nodal displacement components associated with element e ð2Þ , where nne ð2Þ is the number of nodes of element e ð2Þ . The 3 Â 3nne ð2Þ matrix ½N ð2Þ ðgÞ e ð2Þ contains the shape functions of element e ð2Þ .
As in [3], given a quadrature point X in element e ð1Þ & S ð1Þ oc , there is at least one element e ð2Þ & S ð2Þ oc containing the reference point Y of Both element e ð2Þ and the number of nodes nne ð2Þ , which depend on the projection y, vary as a function of the quadrature point X. Moreover, the whole set of quadrature points of element e ð1Þ generally involves several elements e ð2Þ in S ð2Þ oc . All the virtual quantities are interpolated in a similar way. Eventually, the element contact expression ðW Ã contact Þ e ð1Þ in (37) on the element e ð1Þ level can be written in the discretized form: where fU contact g e is the element contact force vector and fR K g e the element residual vector for the multipliers, defined as Vector fU contact g e in (43) has nne ð1Þ þ nne ð2Þ components, where nne ð2Þ varies as a function of the considered quadrature point; vector fR K g e has nne ð1Þ components. Note that In (45), use has been made of the adjective 'algorithmic' in order to emphasize the difference with the 'real' situation. For instance, the 'algorithmic gap' case meansk N ¼ k N þ N g < 0, in contrast with the real 'gap' case which merely means g < 0; and the 'algorithmic slip' case means kk T k P lk N , as opposed to the 'slip' case which means kk T k ¼ lk N .
By adding the element contributions in (42), one obtains the following semi-discrete coupled equation system: where vectors fU contact ðU; KÞg and fR K ðU; KÞg are the assembly of the element vectors fU contact g e jðnÞ and fR K g e jðnÞ given in (43)  in the two bodies 1 and 2, and the multiplier vector fKg on the reference contact surface S ð1Þ oc . Except for the above-mentioned vectors fU contact ðU; KÞg and fR K ðU; KÞg which are specific to contact, the other vectors and matrices are standard and result from the discretization of terms other than W Ã contact in relationship (21): [M] is the mass matrix, fWðUÞg the internal force vector resulting from the strain work in the bodies, fUðUÞg the external force vector (excluding contact force) which may depend on displacement {U} if follower loads (e.g., a pressure) are applied on either body. The nonlinear semi-discrete Eq. (46) will be solved as a whole using the standard full Newton-Raphson method to simultaneously obtain the nodal displacements and the multipliers, as was done in [4,7,13,8] (Curnier and co-workers even used a socalled generalized Newton method which is more sophisticated).
The slip velocity V T ¼ _ g a a a defined in (15) and involved in (43) and (44) has to be computed here by means of an adequate temporally discrete expression. For this purpose, let us first note that Proposition 1 remains valid if V T is replaced with V T Dt, where Dt is an arbitrary finite time interval. The quantity V T Dt 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=m 3 . 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 t nÀ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 t nÀ1 is x nÀ1 ¼ / ð1Þ ðX; t nÀ1 Þ and its current position is x ¼ / ð1Þ ðX; tÞ. We denote the following quantities related to time t nÀ1 : -y nÀ1 the projection of x nÀ1 on the target surface S ð2Þ c ðt nÀ1 Þ at time t nÀ1 , -Y nÀ1 and g nÀ1 the inverse images of y nÀ1 under the chain of parametrizations (5). It should be emphasized that (i) y and y nÀ1 are not the respective positions at times t and t nÀ1 of a same particle in body 2, and (ii) the current position of the particle Y nÀ1 is / ð2Þ ðY nÀ1 ; tÞ, in general distinct from y ¼ / ð2Þ ðY; tÞ.With these notations in hand, using the chain of parametrizations (5), which gives y ¼ w ð2Þ ðg; tÞ and w ð2Þ ðg nÀ1 ; tÞ ¼ U ð2Þ ðY nÀ1 ; tÞ, and applying the backward Euler integration scheme with time step Dt ¼ t À t nÀ1 leads to the following expression for the slip: It should be noted that other expressions such as V T Dt ¼ w ð2Þ ðg; tÞ À w ð2Þ ðg nÀ1 ; tÞ would not be suitable as surface S ð2Þ c is divided into finite elements and the parametrization (5) is in fact defined piecewisely. It may happen that the projection points y nÀ1 and y do not belong to the same element on S ð2Þ c although they are close to each other; in this case their local parent coordinates g nÀ1 and g may be quite distant even though the time step Dt tends to zero. The algebraic expression (47) does not involve the parent coordinates and thus avoids the aforementioned discontinuity problem. It is adopted to compute the slip in the contact force vector (43) and the residual (44), as well as the tangent stiffness in the next section. The important issues of convergence and stability of a mixed finite element formulation were investigated in the literature for continuum problems and later for contact problems. They are related to the inf-sup condition which is expressed in contact problems as a compatibility condition between the displacement and multiplier spaces on the contact surfaces.
Following Chapelle and Bathe's developments [23] on incompressible elasticity, El-Abbasi and Bathe [24] examined the small deformation frictionless contact problem and proposed an infsup test procedure where a series of refined meshes is used in order to evaluate the evolution of the numerical inf-sup value. A further mathematical study was carried out by Bathe and Brezzi in [25], giving rise to precise relations between the degrees of the polynomial interpolations for the displacements and the multipliers. In general, the mathematical study of the inf-sup condition is confined to the small deformations framework since it is difficult when the structure undergoes large deformations and the contact surface is unknown a priori. No attempt has been made in this work in order to prove-analytically or numerically-that the proposed method passes the inf-sup condition. This issue should be investigated in the future. In the meanwhile, it can be observed through the numerical examples presented in Section 6 that the convergence is achieved when the mesh is refined and there is neither locking nor spurious contact pressure.

Contact tangent stiffness matrix
This section is devoted to the derivation of the tangent stiffness matrix required for the Newton-Raphson solution of the coupled equation system (46). We will confine ourselves to that part of the tangent stiffness related to contact, since the remaining part related to non contact terms can be linearized in a standard way. Let dU and dk be the independent variations of the displacement and the multiplier, respectively defined by where Latin indices take the values 1, 2, 3, Greek indices take the values 1, 2, repeated indices imply summation and the variations dU i ; i 2 f1; 2; 3g; dk N ; dk T1 ; dk T2 are arbitrary. One has to compute the exact linearization of W Ã contact corresponding to these variations and extract the contact stiffness ½K contact from the relation: The variation ðdW Ã contact Þ e ð1Þ over any element e ð1Þ & S ð1Þ oc is obtained by taking the directional derivatives of (37) in the directions dU ðiÞ ; i 2 f1; 2g; dk N and dk T . The variations of the basic kinematic variables are expressed in terms of the displacement variations via [2,3]: ;b g ð Þ; where dU ð2Þ ;a ðgÞ denotes the partial derivative with respect to g a of the composite function dU ð2Þ ðgÞ ¼ dU ð2Þ ðYðgÞÞ.
As for the variation of the slip, it is derived from (47): where W and Z a ; a 2 f1; 2g, are non symmetric second-order ten- In the above relationship, fdkg for instance is the 3-component column vector representing dk in the fixed basis ðe 1 ; e 2 ; e 3 Þ and the column vector fdkg e ð1Þ with 3nne ð1Þ components contains the triplets of the nodal components of dk in the same basis. One also needs to interpolate the variation fdU ð2Þ ðY nÀ1 Þg involved in (51): where the element contact stiffness ½K contact e is a rectangular matrix of the form: The explicit expression of the above matrix is given in the Appendix. It should be emphasized that (i) despite the superscript e which stands for element, ½K contact e is actually defined for each quadrature point of element e ð1Þ & S ð1Þ oc and (ii) although the element matrices ½K contact e jðnÞ are rectangular, their assembly eventually results in the square contact stiffness ½K contact of the form:

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. Although the semi-discrete equation system (46) holds in dynamics, from now on we will confine ourselves to those quasistatic cases that display enough specific features of the contact phenomenon to consider the formulation validated. The dynamic cases present additional physical and numerical difficulties and are deliberately disregarded at this stage. All numerical results presented here have been obtained by a home made finite element code.

Contact patch test
The first example we consider is related to the contact patch test which was originally conceived in [26,27], with the aim of checking whether the finite element contact formulation allows transmission of a constant pressure through an arbitrary nonconforming contact surface.
In this work, we used the patch test proposed in [28,24] and shown in Fig. 3a. An elastic body is pressed against another one by a uniform distributed load p = 1000 N/m 2 applied on its top surface. Assuming a frictionless interface ðl ¼ 0Þ between the two bodies, the exact solution of the patch test entails a constant contact pressure equal to p along the contact line. The patch test is passed if the numerical results yield such a uniform stress distribution, regardless of the mesh used. In practice, one has to verify the test using some typical meshes.
Here the problem is solved assuming that the two bodies are made of an isotropic hyperelastic Saint-Venant Kirchhoff material characterized by the strain energy function in terms of the Green strain E: where the Lamé constants k; l are related to Young's modulus E and Poisson's ratio m by the same relations as in linear elasticity: The problem is modeled in plane strain conditions, taking the dimensions of the bodies a properties E ¼ 10 6 N=m 2 and m ¼ 0. The structure is discretized using two meshes -the first one is regular and the second distorted -as shown in Fig. 3b-c, in both cases the mesh is nonconforming along the contact interface. The two bodies are discretized with 8-node quadrilaterals (Q8), so that the elements on the contact surface are 3-node line elements. The numerical quadrature is performed using 3 Â 3 Gauss quadrature points in the Q8 elements and 7 Gauss points in line elements belonging to the contact surface. The penalty parameter is N ¼ 10 12 N=m 3 . Fig. 4 shows the deformed configurations and contact stress distributions along the contact line when the upper, respectively the lower, body is selected to be the contactor. The arrows represent the contact stresses exerted by the target body upon the contactor. In all cases, the contact stresses are found to be constant and equal to p up to 5 or 6 significant digits.
All the nodes on the contact interface have the same vertical displacement U y ¼ À0:0010015 m. The stress field inside the two bodies is constant, with the only non-zero stress component R yy ¼ À1001 N=m 2 ðR is the second Piola-Kirchhoff stress tensor).
These values agree with the analytical solution of the nonlinear elasticity (nonlinear elasticity must be used for the comparison since the numerical results are obtained by a nonlinear finite element code).

Hertzian contact problem
In this example, we consider a Hertzian contact problem in order to investigate the accuracy of the proposed contact formulation. Referring to Fig. 5a, a 2D elastic cylinder with radius R ¼ 10 m is pressed against a rigid foundation by a force F which is perpendicular to the rigid plane and applied on the symmetry axis of the cylinder. The cylinder is made of the same isotropic hyperelastic Saint-Venant Kirchhoff material as in (58), here we take E ¼ 100 N=m 2 ; m ¼ 0:3. The computations are performed assuming plane strain conditions and either a frictionless ðl ¼ 0Þ or a frictional ðl ¼ 0:3Þ contact between the cylinder and the foundation.
In the frictionless case, the analytical solution assuming the linear small strain theory of elasticity is well known. The pressure distribution in the contact region is given by where b is the half contact region width: The computation is carried out with two values of the total load: F = 1N/m and F = 2 N/m. Due to symmetry, only half of the structure is considered. Two meshes are considered: a coarse one containing 198 elements and 455 nodes, and a finer one containing 345 elements and 773 nodes (the latter is shown in Fig. 5b). In both cases, the rigid foundation is chosen as the target and modeled by one single element with fixed nodes. Clearly, the meshes are nonconforming.
The numerical quadrature employs 7 points in the triangles and the elements belonging to the contact line. The penalty parameter is N ¼ 100 N=m 3 . The (relative) convergence tolerances are taken to be 10 À6 for the force residual fRðU; KÞg ¼ fWðUÞgÀ fUðUÞg À fU contact ðU; KÞg in (46a) and the multiplier residual fR K ðU; KÞg in (46b). Fig. 6 plotting the numerically computed contact pressure acting on the cylinder shows that the finite element solution converges as the mesh is refined.
It should be noted that for a given mesh and a given load, there is in general one element of the cylinder which overlaps the contact area and the free surface, as can be seen in Fig. 6. In such an element the contact stress vector-or the multiplier k-differs from zero in the contact zone and vanishes outside, and the interpolation of k from the nodal multipliers, see Relation (39c), is not precise. Despite this difficulty, the numerical results obtained here are in good agreement with the theory. With the finer mesh and F = 1 N/m, the maximal pressure at x ¼ 0 is 2% lower than the theoretical value. With the larger load F = 2 N/m, the numerical results deviate more from the theory, the maximal pressure at x ¼ 0 is 2:8% lower than the theoretical value.
In actual fact, the contact stresses obtained here are the nominal Piola-Kirchhoff ones (i.e. the multipliers according to (22)) whereas the analytical pressure in (60) comes from the linear elasticity theory. Although the deformation is small, there still is a slight difference between these stresses. An elaborate calculation shows that the Cauchy stresses corresponding to the numerical Piola-Kirchhoff ones are twice closer to the theoretical values. Having said this, we decide to only present the Piola-Kirchhoff stresses throughout this work.
In the frictional case, the computations are made with the fine mesh only. The analytical solution was obtained by Spence [29,30], assuming that the influence of friction on the normal stress is negligible in order to uncouple the equations governing the normal and tangential stresses. The contact region is subdivided into a stick zone jxj 6 c and a slip zone c 6 jxj 6 b. The extent b of the contact region is given by the same relation (61)  The normal pressure is given by the same relation (60) as in the frictionless case while the tangential contact stresses is computed via: where Kðk; wÞ is the incomplete elliptic integral of the first kind, The total load F = 1 N/m is applied in 50 equal increments, the penalty parameters are N ¼ T ¼ 100 N=m 3 . In each time step, the solution is obtained after 3 iterations in average, with a maximum of 7 iterations at the first time step. The shape of the stresses along the contact interface is displayed in Fig. 7a, which clearly shows the progressive transition from stick to slip over the contact line. The normal T y and tangential T x contact stresses are plotted in Fig. 7b, showing a good agreement with Spence's theoretical solution, except in the central stick region where the friction seems to have the effect of raising the normal contact pressure. The value at x ¼ 0 of the normal stress exceeds the analytical one by about 4%, showing the same trend as in [12] where a similar problem of the contact between a rigid punch and an elastic halfspace was investigated.
In order to check whether the finite element procedure enables one to satisfy the local inequality constraints, the ratio T x =lT y is computed at every node on the contact line for F ¼ 0:5 and 1 N/ m and it is verified that the Coulomb inequality jT x j 6 lT y holds everywhere, see Fig. 7c. At the sliding nodes, the equality jT x j ¼ lT y is met within 0.1%.
x(m) Eventually, two additional computations are carried out with only one and 10 increments instead of 50. As shown in Fig. 8, the contact width, the normal contact pressure as well as the frictional shear in the slip region can be recovered with good precision using one or 10 time steps, whereas the frictional shear distribution in the stick zone is rather sensitive to the load increment and requires up to 50 time steps to be accurate. This result is in full accord with Pietrzak's observation in [12].

Contact between two beams
In this example, we are concerned with the quasistatic contact between two beams, one of which being subjected to either a dead or a follower force, with a view to showing the ability of the proposed formulation to deal with large deformation contact. The two beams have identical reference geometries, with rectangular cross section, a length of L ¼ 1 m along x, a height of 0.05m along y and a thickness of 0.1m along z, see Fig. 9. They are both made of an isotropic hyperelastic Saint-Venant Kirchhoff material with E ¼ 2:10 6 N=m 2 and m ¼ 0:3.
The beams are clamped at end X ¼ 0. In the reference configuration, the lower face of the upper beam (contactor) coincide with the upper face of the lower beam (target). Each beam is modeled with 15 20-node brick (i.e. hexahedral) elements (H20) and the contact surface with 15 8-node quadrilaterals (Q8). The numerical quadrature employs 3 Â 3 Â 3 Gauss points in the H20's and 6 Â 6 points in the Q8 elements belonging to the contact surface.
In the following, two types of loading, a dead and a follower force, are applied on the contactor. For each loading, the computations are carried out assuming either frictionless or frictional contact, l ¼ 0:3. In any case, the penalty parameters are N ¼ T ¼ 10 5 N=m 3 .

Dead load
A dead surface load ÀFy is uniformly distributed over the cross section X ¼ 1 m of the upper beam. Force F is progressively in-creased from zero using an arc length method with the arc length regularly incremented by a factor 1.1. The (relative) convergence tolerances are taken to be 10 À6 for the force residual fRðU; KÞg in (46), 10 À6 for the arc length residual and 10 À10 for the multiplier residual fR K ðU; KÞg in (46b).
In the frictionless case, the computation is carried out in 29 steps, the resultant force at the last step is F ¼ 62:62 N and the convergence is achieved in 7-8 iterations in average. Fig. 10 shows the deformed configurations and the contact traction vectors acting on the upper beam. Only x-y in-plane views are shown since the complete 3D views are not legible. In the frictional case, the entire calculation requires 38 steps, the final resultant force is F ¼ 61:20 N and the convergence achieved in 11-12 iterations in average. The corresponding deformed configurations and the contact tractions are shown in Fig. 11.
The contact immediately takes place after force F is applied and persists over the entire lower surface of the upper beam until the last step. During the first steps, the contact tractions are higher near the end cross section where force F is applied, whereas far from this section, the contact tractions are smaller and barely visible in Figs. 10 and 11. As the deformations grow up, the contact tractions become significant in the neighborhood of the built-in section. It should be mentioned that there are possible spatial oscillations of the contact stresses due to the use of quadratic shape functions which take negative values at certain points and due to coarse meshes. Thus, during the bending process, there may be some nodes on the contact surface where the contact law is not satisfied, namely T N < 0 or kT T k > l Á T N in frictional case.
Fortunately, these violations only occur at the nodes where T N and kT T k are very small, specifically smaller than about 1% of the maximal T N and kT T k values on the whole contact surface. Moreover, the negative contact pressures disappear if the mesh is fine enough. In Figs. 10 and 11, the contact tractions violating the contact law are hardly visible. In presence of the friction, the contact state is rather complex. Roughly speaking, it is found that there  is slip at most nodes and every step, and stick at some nodes and some steps. Fig. 12 shows that the deflection at point A (central point of the cross section X ¼ 1 m of the upper beam, as shown in Fig. 9) versus force F is almost the same, with or without friction. If the friction prevented the slip to occur at some steps, force F necessary to reach the same deflection would be higher than in the frictionless case and the two plots would be distinct.

Follower force
The dead force considered above is now replaced with a uniform pressure distributed over the last Q8 element of the upper face of the upper beam, next to the cross section X ¼ 1 m. The pressure is a follower force in the sense that it remains continually perpendicular to the area on which it is applied, so that its direction changes as the beam deforms, contrary to the dead force which keeps a fixed direction. We choose the same automatic control of the pressure and convergence tolerances as in the dead load case.
In the frictionless case, the computation is carried out in 34 steps leading to the final pressure resultant P = 73.84 N, the convergence is achieved in 9-10 iterations in average and the deformed configurations together with the contact stresses are shown in Fig. 13. In the frictional case ðl ¼ 0:3Þ, the computation is carried out in 43 steps, the final pressure resultant is P = 74.83 N, the convergence is achieved in 10-11 iterations in average and the results are shown in Fig. 14.
As shown in Fig. 13, at the beginning the contact occurs on the whole lower surface of the upper beam. At an advanced deformed state, the absence of friction leads to a release over a large portion of the former contact surface, as shown in Fig. 11e and f. On the other hand, in the presence of sufficiently high friction, the contact is maintained all over the lower surface of the upper beam, as shown in Fig. 14.
As in the dead load case, the contact state in the presence of friction is rather complex. It is found that at every step there is slip at most nodes, except the last step where there is more stick than slip. Fig. 15 shows the deflection at point A of the upper beam versus pressure resultant P. Note that under the follower loading there is a limit point for the deflection, the maximum absolute value for the deflection being equal to about three-quarter of the beam length.

Contact between two membranes
This example relating to the contact between two pressurized membranes involves large displacement, large strain and large slip as well. The problem is modeled in plane strain with respect to the x-y plane. In the reference configuration, both membranes are 1 mm thick and represented by line segments with a width of 1 m along x. The membranes are Dy ¼ 1 m apart in the direction y and their centers are shifted by a distance Dx ¼ 1 m along x, as shown in Fig. 16. The edges of the membranes parallel to z-axis are fixed.  The membranes are made of an isotropic hyperelastic material of the neo-Hookean type, defined by the strain energy function in terms of the invariants I 1 ¼ trC and I 3 ¼ det C of the right Cauchy-Green tensor C [31]: We take E ¼ 2:10 6 N=m 2 ; m ¼ 0:3, where Young's modulus E and Poisson's ratio m are related to constants k and l by the same relations (59). The second Piola-Kirchhoff stress tensor R is The membranes are subjected to the same pressure value p which brings them into contact. The computations are carried out assuming either no or with friction, l ¼ 0:5, a rather significant coefficient of friction in order to simulate the contact between rubber membranes.
Each membrane is divided into 20 3-node line elements of the membrane type (i.e. with no bending stiffness), in plane strain with respect to x-y plane and plane stress with respect to the membrane thickness direction. Here, the mesh of the contactor (upper membrane) is also used as the mesh of the contact surface. In order to avoid the singularity of the stiffness matrix at the beginning of the computation-which is a common phenomenon in membrane problems-a small artificial prestress is applied at the first step and removed after.
The numerical quadrature is carried using 3 Gauss points in the line elements to compute the internal force vector, 6 Gauss points to compute the contact force vector. The penalty parameters are N ¼ T ¼ 2:10 4 N=m 3 . The pressure p is controlled by the arc length method with the arc length regularly incremented by 1.05. The convergence tolerances are 10 À6 for the force residual fRg; 10 À6 for the arc length residual and 10 À8 for the multiplier residual fR K g.
Without friction, the computation is made in 43 steps, the maximum pressure value is p ¼ 4023 N=m 2 and the convergence achieved in 4 iterations in average. With friction, the computation is made in 124 steps, the maximum pressure value is p ¼ 5420 N=m 2 and the convergence achieved in 4 iterations in average. Fig. 17 shows the deformed configurations and the contact stresses exerted on the upper membrane in the frictionless case. Let us mark out two particular particles belonging respectively to the two membranes, identified by letters A and B in Fig. 17. At the contact onset ðp ¼ 2049 N=m 2 Þ, their locations A 1 and B 1 are almost the same, see Fig. 17a. When the pressure is p ¼ 4023 N=m 2 , their locations A 2 and B 2 are quite distinct as indicated in Fig. 17b. The distance A 2 B 2 is a (poor) approximation of the slip V T Dt, computed at the particle marked by A 1 and between the two time instants corresponding to the two quoted pressures, see (47) and Fig. 2. Fig. 18 shows the deformed configurations and the contact stresses exerted on the upper membrane in the frictional case. Note that the contact stresses in Figs. 17 and 18 are not drawn at the same scale in order to optimize the drawings.
While the shape of final deformed configuration is similar to that in the frictionless case, the displacements of the particles are quite different. To see this, consider the above-mentioned pair of particles again, with almost identical locations A 1 and B 1 at the beginning of the contact (pressure p ¼ 2049 N=m 2 ). This time, due to significant friction, it is found that, up to pressure p ¼ 5420 N=m 2 corresponding to an advanced deformed shape, these two particles almost stick together, as shown in Fig. 18b. Otherwise stated, the slip velocity at the particle marked by A 1 is virtually zero. An in-depth analysis at different steps shows that there is stick in the neighborhood of particle A 1 , and there may be slip at contact nodes outside this neighborhood. Fig. 19 plots the deflection of the middle point of membrane 2 (point C in Fig. 16) versus pressure p. As predicted, the deflection is much smaller in the presence of friction. The resultant contact forces F x and F y between the membranes along x-and y-direction, respectively, together with the magnitude jFj ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi shown in Fig. 20 versus pressure p. It can be seen that for a given pressure the magnitude jFj is a little smaller in frictional than in frictionless case whereas the force direction, determined by F x and F y , differs significantly.

Conical extrusion of a cylinder
This section deals with the quasistatic extrusion of a cylinder into a rigid conical die, Fig. 21. The geometry data are taken from [1,3]: the wall angle of the die is 5°and in the reference configuration, the cylinder is 25.4 cm long and 5.08 cm in radius.
In this example, use should be made of a suitable constitutive model in order to overcome computational difficulties due to extremely high compressive stresses in the bodies at some advanced deformed stage. In the classical elastoplasticity model, the elastic response is represented by the Saint-Venant Kirchhoff law which has an unstable behavior under compression, so that in the present example the cylinder is found unable to move forward when it is compressed enough. Modification of the plastic part does not help to improve the situation and this difficulty only disappears when the Saint-Venant Kirchhoff law is replaced with an appropriate one, e.g., by using the hyperelastic stored energy given in [3] p37. Here, we choose the so-called Mooney rate-independent elastoplasticity model instead, which is described in [12] and proved to  behave correctly even for huge deformations both in compression and tension. This model is chosen merely because it is easier to be implemented in our finite element code. The question about which elastoplasticity model best predicts the extrusion response when compared to experimental results is crucial but out of the scope of this work. Here, the extrusion computations using the chosen elastoplasticity model only aim to show the ability of the formulation to successfully deal with contact in large deformation with a nonlinear constitutive law of the elastoplastic type. The cylinder is made of an elastoplastic material with a linear isotropic hardening. The elastoplastic formulation and the material data are the same as in [12]: Young's modulus E = 100000 MPa, Poisson's ratio m ¼ 0:3, yield stress r 0 = 300 MPa and hardening modulus H = 3560 MPa. Note that the geometry in [12] is not considered since it is not described with enough precision (the die is not conical and its geometry is rather complex).
The problem is modeled in axisymmetry. The cylinder (taken as the contactor) is discretized with 20 8-node quadrilaterals (Q8) interconnected by 74 nodes. The 10 elements on the contact surface are 3-node line elements. The rigid die (the target) is modeled with one single Q8 element. The numerical quadrature is performed using 3 Â 3 Gauss quadrature points in the Q8 elements, 7 Gauss points in line elements belonging to the contact surface. The computation are made assuming either frictionless or frictional contact between the cylinder and the die wall ðl ¼ 0:1Þ. The penalty parameters are N ¼ 10 11 N=m 3 ; T ¼ 5:10 10 N=m 3 .   The numerical computation is carried out by regularly incrementing the prescribed axial displacement on the left end of the cylinder by factor 1.1. The convergence tolerances are 10 À6 for the force residual {R} and the multiplier residual fR K g.
In the frictionless case, the computation is achieved in 31 steps, the axial displacement of the pushed end is 18.2 cm at the last step and the convergence requires 4-5 iterations in average. Fig. 22 displays the deformed configurations together with the contours of the effective plastic strain while Fig. 23 shows the contact stresses exerted on the extruded cylinder. The plastic strain is concentrated near the cylinder head (where the contact stresses are highest, too) and reaches about 130% at some points.
With friction, the computation is achieved in 109 steps, the axial displacement at the last step is 20.4 cm and the convergence requires 4 iterations in average. Fig. 24 displays the deformed configurations together with the contours of the effective plastic strain while Fig. 25 shows the contact stresses.
From the numerical results it is found that all the nodes on the contact line are sliding at every loading step and the equality jT T j ¼ l Á T N is accurately rendered to the last displayed digit (6th digit).
Note that at the final step the contact stresses near the cylinder head increase more than twice as compared with the frictionless case, as shown in Fig. 23d and Fig. 25d. The plastic strain is concentrated near the cylinder head where the contact stresses are strongest, with a maximal value of about 150%. Fig. 26 plots the extrusion force F versus the displacement U of the left end of the cylinder, in frictional and frictionless cases. As can be predicted, the extrusion force increases as the cylinder is moving forward and the contact surface between the cylinder and the die is larger. Also, given an amount prescribed displacement U, force F is always higher in frictional than frictionless case. When the end displacement U is about 18 cm, F is about 4 times as big as in the frictionless case. The extrusion force F in the frictional  case plotted in Fig. 26 shows the same trend as that obtained in [3] whereas the numerical values differ rather significantly since the material model and the material data are not identical. Fig. 27 shows that the volume of the cylinder diminishes gradually during the extrusion process and the relative volume change is maximal in absolute value at the last step, equal to À6.8% in the frictionless case and À15.9% in the frictional case. The volume decrease is basically the same in both cases as long as the prescribed displacement is less than about half the reference length of the cylinder. When the prescribed displacement exceeds this value, the volume decreases more rapidly in the presence of friction.

Conclusions
The main features of the present formulation can be summarized as follows: A weighted residual relationship is proposed in (21) as an extension of the standard virtual work principle to cope with the large deformation contact problem with Coulomb friction. It is inspired by the augmented Lagrangian formulation of Pietrzak and Curnier [13] and amounts to rewrite the augmented Lagrangian formulation of Alart and Curnier [4] in a generic weak form. Proposition 1 shows that the proposed weak form is equivalent to the strong form of the contact problem. It is noteworthy that equality (21) yields the inequality constraints of contact. Another interesting fact is that Relation (23) which expresses the equilibrium of the traction vectors at the contact interface is found to be a consequence of the weak form (21).Apart from the analogy with an augmented Lagrangian principle, the weak form (21) is stated as a genuine weighted residual relationship in that it involves arbitrary (regular) virtual displacements and multipliers. As a consequence, it can be discretized by means of the finite element method in a simple way. The finite element discretization process holds for any finite element type on the contact surface. It is only required that the finite elements are isoparametric, so that the current and reference coordinates can be interpolated in the same way. Any discrete contact quantity -namely, the element contact force vector fU contact g e in (43), the element residual vector for the multipliers fR K g e in (44) or the element contact stiffness matrice ½K contact e in (56) -is defined at a quadrature point X of an element e ð1Þ belonging to the contactor surface S ð1Þ oc . Its numerical computation involves the nodal quantities associated to two elements: element e ð1Þ and the (or any) element e ð2Þ in the target surface S ð2Þ oc containing the projection point Y of the quadrature point X in question. Element e ð2Þ , which depends on Y, varies as a function of X. This feature is common to the formulation given by Laursen [3]. The temporally discrete relation (47) proposed for the slip V T Dt does not involve the parent coordinate g of the projection point Y and thus avoids difficulties due to the discontinuity in variable g related to the finite discretization of S ð2Þ oc . Instead, it involves variables y and Y nÀ1 which have intrinsic kinematic meanings independent of the finite discretization of S ð2Þ oc , as clearly shown in Fig. 2. The element contact stiffness matrices related to relative tangential displacement, ½k 4 contact e and ½k 5 contact e in (A6) and (A7), are rectangular since they also involve the projection point Y nÀ1 determined at the previous time step. Nevertheless, their assembly eventually yields the global contact stiffness matrix which is square as usual. As regards the smoothness of the kinematic functions, it should be noted that the present formulation involves differentiations of order 2 at most (see, for instance, the definition of the curvature j ab in (14)), whereas other formulations may require higher order differentiations. As the weak form proposed in the present formulation already incorporates the contact laws, no integration of the frictional equations are required on the local level of the numerical solution procedure. This feature resembles, in some ways, an implicit property in Curnier and co-workers' formulations [4,13] and differs from most of other developments in the literature.
The numerical examples presented in quasistatics shows the ability of the proposed formulation to deal with large deformation contact. Investigations are in progress to deal with numerical examples in dynamic contact.

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 1.

Appendix. The element contact stiffness matrix
This Appendix gives the explicit expressions for the four blocks of the element contact stiffness matrix (56), namely ½K UU e ; ½K UK e ; ½K KU e and ½K KK e .
(1) The rectangular matrix ½K UU e in (56) pertains to the displacements fU ð1Þ g e ð1Þ on e ð1Þ ; fU ð2Þ g e ð2Þ on e ð2Þ and fU ð2Þ g e ð2Þ nÀ1 on e ð2Þ nÀ1 . It is of dimension ð3nne ð1Þ þ 3nne ð2Þ Þ Â ð3nne ð1Þ þ 3nne ð2Þ þ 3nne In (A4), the notation [T] for any second order tensor T designates the 3 Â 3 matrix representing T in the fixed basis ðe 1 ; e 2 ; e 3 Þ. For instance, if T ¼ T b a a a a b ¼ T ab a a a b , then the (p,q)-th element of matrix ½T is T pq ¼ e p Á T Á e q ¼ T b a ða a Þ p ða b Þ q ¼ T ab ða a Þ p ða b Þ q , where ða a Þ p is the p-th component of a a in the fixed basis. As j ac a cb a a a b ¼ a ac j cb a a a b is a symmetric tensor, so is matrix ½k : : ðA5Þ The contribution of point g nÀ1 is found in the rectangular matrices ½k contact e and ½k 5 contact e of dimension ð3nne ð1Þ þ 3nne ð2Þ Þ Â ð3nne ð1Þ þ 3nne ð2Þ þ 3nne ð2Þ nÀ1 Þ: In (A6) and (A7), ½W; ½Z a are 3 Â 3 matrices representing the tensors W; Z a defined by (52), in the fixed basis ðe 1 ; e 2 ; e 3 Þ : W ij ¼ e i :W:e j ¼ a ab ða a Þ i ða b Þ j ; Z a ij ¼ e i :Z a Á e j ¼ g a ab ða b Þ i m j . The last matrix in (A1), ½k 6 contact e , is square of the same dimension as ½k 1 contact e , yet not symmetric: (2) The matrices ½K UK e and ½K KU e in (56) representing the coupling between the displacements and the multipliers are rectangular, of dimension ð3nne ð1Þ þ 3nne ð2Þ Þ Â 3nne ð1Þ and 3nne ð1Þ Â ð3nne ð1Þ þ 3nne ð2Þ Þ, respectively: and ; ðA6Þ : ðA7Þ : ðA8Þ