Modelling contact with isotropic and anisotropic friction by the bipotential approach

Based on an extension of Fenchel's inequality, the bipotential approach is a non smooth mechanics tool used to model various non associative multivalued constitutive laws of dissipative materials (friction contact, soils, cyclic plasticity of metals, damage). Generally, such constitutive laws are given by a graph $M$. We propose a simple necessary and sufficient condition for the existence of a bipotential $b$ for which $M$ is the set of couples $(x,y)$ of dual variables such that $b(x,y) = \langle x,y \rangle$, and a method to construct such a bipotential by covering $M$ with cyclically monotone graphs which are not necessarily maximal (bipotential convex cover). As application, we show how to obtain the bipotential of the law of unilateral contact with Coulomb's friction by a bipotential convex cover. Introduced to extend the classical calculus of variation, the bipotential concept is also useful to construct numerical schemes for friction contact laws. In recents works, we extended the bipotential approach to a certain class of orthotropic frictional contact with a non-associated sliding rule proposed by Michalowski and Mroz. The bipotential suggests a predictor-corrector numerical scheme.

1. Basic tools. X and Y are topological, locally convex, real vector spaces of dual variables x ∈ X and y ∈ Y , with the duality product •, • : X × Y → R. We shall suppose that X, Y have topologies compatible with the duality product. We use the notation:R = R ∪ {+∞}. For any convex and closed set A ⊂ X, its indicator function, χ A , is defined by The subgradient of a function φ : X →R at a point x ∈ X is the (possibly empty) set: The conjugate is always convex and lower semi continuous (lsc).
We denote by Γ(X) the class of convex and lower semicontinuous functions φ : X →R. The class of convex and lower semicontinuous functions φ : X → R different of the constant +∞ is denoted by Γ 0 (X).

GÉRY DE SAXCÉ
A graph M is cyclically monotone if for all integer m > 0 and any finite family of couples (x j , y j ) ∈ M, j = 0, 1, . . . , m, A cyclically monotone graph M is maximal if it does not admit a strict prolongation which is cyclically monotone. By reindexing the couples, we easily recast the previous inequality as what shows that the graphs of a law and its dual law are simultaneously cyclically monotone. Rockafellar [25] Theorem 24.8 (see also Moreau [23] Proposition 12.2) proved a theorem that can be stated as: Theorem 1.1. Given a graph M , there exist potentials φ ∈ Γ 0 (X) such that M ⊂ Graph(∂φ) if and only if M is cyclically monotone. They are defined by where x 0 and φ(x 0 ) are arbitrarily fixed and the 'sup' is extended to any m > 0 and to any couples (x k , y k ) ∈ M, k = 1, 2, . . . , m.
Because the dual law is also cyclically monotone, we can apply once again the construction of the previous Theorem, giving the function ψ(y) = sup x m , y − y m + m k=1 x k−1 , y k − y k−1 + ψ(y 0 ), (4) such that M ⊂ M (∂ψ * ). Excepted when M is maximal, φ and ψ * are in general distinct function, as it will be seen further in the application.

2.
What is a bipotential? The constitutive laws of the materials can be represented, as in Elasticity, by a univalued mapping T : X → Y or, as in Plasticity, can be generalized in the form of a multivalued mapping T : X → 2 Y , but this representation is not necessarily convenient. Its graph M = Graph(T ) is a non empty part of X × Y . When the graph is maximal cyclically monotone, we can modelize it thanks to a convex and lower semi-continuous (l.s.c.) function φ : X →R, called a superpotential (or pseudo-potential), such that the graph M is the one of its subdifferential Graph(∂φ). The dissipative materials admitting a superpotential of dissipation are often qualified as standard [17] and the law is said to be a normality law, a subnormality law or an associated law. However, many experimental laws proposed these last decades, particularly in Plasticity, are non associated. For such laws, we proposed in [12] a suitable modelization thanks to a function called bipotential.
Definition 2.1. A bipotential is a function b : X × Y →R, with the properties: (a) b is convex and lower semicontinuous in each argument; (b) for any x ∈ X, y ∈ Y we have b(x, y) ≥ x, y ; (c) for any (x, y) ∈ X × Y we have the equivalences: MODELLING CONTACT WITH ISOTROPIC AND ANISOTROPIC

411
The graph of b is If the graph M of a law is the graph of a bipotential b, we say that the law (the graph) admits a bipotential. Particularly, to each superpotential φ we can associate the separable bipotential: where φ * is the Fenchel conjugate of φ. Hence, the cornerstone inequality of the bipotential (definition 2.1 (b)) is reduced to well known Fenchel's inequality [13]. We introduce also in [5], [6] the notion of a strong bipotential. Conditions (B1S) and (B2S) appear as relations (51), (52) in Laborde and Renard [21].
(a) b is convex and lower semicontinuous in each argument; Proposition 2.1. Any strong bipotential is a bipotential.
The notion of a strong bipotential (introduced in relations (51), (52) [21]) is motivated also by the fact that all bipotentials considered in applications to mechanics are in fact strong bipotentials.
The introduction of non separable bipotentials allows modelizing in a more general way the non associated constitutive laws. The laws admitting a bipotential are called laws of implicit standard materials because the relation y ∈ ∂b(•, y)(x) is a subnormality law but the relation between x and y is implicit. Linked to the structural mechanics and in particular with the Calculus of Variation, the bipotential theory offers an elegant framework to modelize a broad spectrum of non associated laws. Examples of such non associated constitutive laws are: non-associated Drücker-Prager [10] and Cam-Clay models [8] in soil mechanics, cyclic Plasticity ( [7], [2]) and Viscoplasticity [18] of metals with non linear kinematical hardening rule, Lemaitre's damage law [1], the coaxial laws ( [9], [26]), the Coulomb's friction law [12], [7], [3], [14], [16], [20], [10], [12], [21]. A complete survey can be found in [9]. In the previous works, robust numerical algorithms were proposed to solve structural mechanics problems.
3. Existence and non uniqueness of the bipotential. For all these particular constitutive laws, the bipotentials were heuristically constructed, without knowing beforehand the conditions under which ones the law admits a bipotential, nor a systematic algorithm to construct this bipotential. This question was answered latter. In [4] we solved two key problems: (a) when the graph of a given multivalued operator can be expressed as the set of critical points of a bipotentials, and (b) a method of construction of a bipotential associated (in the sense of point (a)) to a multivalued, typically non monotone, operator. The main tool was the notion of convex lagrangian cover of the graph of the multivalued operator, and a related notion of implicit convexity of this cover. The results of [4] apply only to bi-convex, bi-closed graphs (for short BB-graphs) admitting at least one convex lagrangian cover by maximal cyclically monotone graphs. This is a rather large class of graph of multivalued operators but important applications to the mechanics, such as the bipotential associated to contact with friction [12], are not in this class.
In more recent papers [5], [6], we proposed an extension of the method presented in [4] to a more general class of BB-graphs. This is done in two steps. In the first step we proved that the intersection of two maximal cyclically monotone graphs is the critical set of a bipotential if and only if a condition formulated in terms of the inf convolution of a family of convex lsc functions is true [6] . In the second step we extended the main result of [4] by replacing the notion of convex lagrangian cover with the one of bipotential convex cover (definition 4.2). In this way we were able to apply our results to the bipotential for the Coulomb's friction law.
The domain of M is by definition Hence the law T assigns at each x ∈ X the section M (x) and the inverse law assigns to each y ∈ Y the section M * (y). Let a constitutive law be given by a graph M . Does it admit a bipotential? The existence problem is easily settled by the following result. The proof can be found in [4]. Then we say that M is bi-convex and bi-closed, or in short that M is a BB-graph. This criterion, simple to verify, allows straightaway moving laws without bipotential aside. If the law is represented by a BB-graph, a closely related topics, is to know whether the bipotential is unique. The answer is no. The proof of the previous result is based on the introduction of the bipotential Therefore, here is a counterexample. If M is cyclically monotone maximal, it admits at least two distinct bipotentials, the separable bipotential defined by (7) and b ∞ . Therefore the graph of the law alone is not sufficient to uniquely define the bipotential. 4. Bipotential convex cover. For a given multivalued constitutive law, Theorem 3.1 does not give a satisfactory bipotential because the bipotential b ∞ is somehow degenerate. We would like to be able to find a bipotential b which is not everywhere infinite outside the graph M . We saw that the graph alone is not sufficient to construct interesting bipotentials. We need more information to start from. This is provided by the notion of bipotential convex cover.
Let Bp(X, Y ) be the set of all bipotentials b : X × Y →R. We shall need the following definitions.
Definition 4.1. Let Λ be an arbitrary non empty set and V a real vector space. The function f : Λ × V →R is implicitly convex if for any two elements (λ 1 , z 1 ), (λ 2 , z 2 ) ∈ Λ × V and for any two numbers α, β ∈ [0, 1] with α + β = 1 there exists λ ∈ Λ such that Definition 4.2. A bipotential convex cover of the non empty set M is a function λ ∈ Λ → b λ from Λ with values in the set Bp(X, Y ), with the properties: (a) The set Λ is a non empty compact topological space, Then for any x ∈ X and for any y ∈ Y the functions f (•, x, •) : Λ × Y →R and f (•, •, y) : Λ × X →R are lower semi continuous on the product spaces Λ × Y and respectively Λ × X endowed with the standard topology, A bipotential convex cover is in some sense described by the collection {b λ : λ ∈ Λ}. This is this point of view that we will adopt in the sequel. The next result defines under which conditions the notion of bipotential convex cover is independent of the choice of the parametrization [5].
be a bipotential convex cover and g : Λ → Λ be a continuous, invertible, with continuous inverse, function. Then The next theorem, proved in [5], is the key result needed further.
Then b is a bipotential and M = M (b).
The result is rather surprising because an inferior envelop of functions, even convex, is not generally a convex function. The property (d) of the Definition 4.2 is essential to ensure the convexity properties of b.

5.
Bipotential for cyclically monotone graphs. Maximal cyclically monotone graphs are critical sets of separable bipotentials. For a non maximal cyclically monotone graph M , Rockafellar's theorem [25] claims only that there exists a superpotential φ such that M ⊂ Graph(∂φ). Hence φ is not sufficient to define unambiguously M . In this section, we show that M can be characterized unequivocally by a bipotential b = max(b 1 , b 2 ) where b 1 and b 2 are separable bipotentials.
is the intersection of two maximal cyclically monotone graphs.
The demonstration of the Theorem is given in [5]. 6. Example. For any BB-graph M , let us show how to construct b ∞ . For any u ∈ dom(M ), the graph M u = {u} × M (u) is cyclically monotone. Indeed, for any finite family of couples (x j , y j ) ∈ M, j = 0, 1, . . . , m, we have x j = u for all j and The sets M u cover obviously the graph M when u runs in dom(M ). As x j = u in (3) and putting v = y m , the function defined by Theorem 1.1 is reduced to and taking into account x j = u for all j, the function defined by (4) is If x does not belongs to dom(M ), the function is equal to +∞. Otherwise, we choose u = x to minimise, that gives Of course, as discussed before, this bipotential is rather trivial but the method allows to get more interesting ones by choosing suitable bipotential convex covers as in the next Section. 7. Application to unilateral contact with Coulomb's dry friction. To be short, the space X = R 3 is the one of relative velocities between points of two bodies, and the space Y , identified also to R 3 , is the one of the contact reaction stresses. The duality product is the usual scalar product. We put whereu n is the gap velocity,u t is the sliding velocity, r n is the contact pressure and r t is the friction stress. The friction coefficient is µ > 0. The graph of the law of unilateral contact with Coulomb's dry friction is defined as the union of three sets, respectively corresponding to the 'body separation', the 'sticking' and the 'sliding'.
It is well known that this graph is not monotone, then not cyclically monotone. As usual, we introduce Coulomb's cone and its conjugate cone . Now, we define some sets useful in the sequel. Let us consider p > 0 and the closed convex disc obtained by cutting Coulomb's cone at the level r n = p D(p) = r t ∈ R 2 | r t ≤ µp .
Therefore, for each value of p > 0, we define a set of 'sticking couples' and a set of 'sliding couples' So, we can cover the graph M by the set of following subgraphs parameterized by All these subgraphs are cyclically monotone but none of them is maximal. by applying twice Rockafellar's Theorem [25] to M p , let us construct the corresponding superpotentials φ p : X →R such that M p ⊂ Graph(∂φ p ) and ψ p : Y →R such that M p ⊂ Graph(∂ψ p ). It is worthwhile to observe that ψ p is not the Fenchel conjugate of φ p because M p is not maximal. To fix the arbitrary constant in the definition of the superpotentials, we suppose that φ p (0) = ψ p (0) = 0. For p ∈ (0, +∞), the computations give φ p (−u) = −pu n + µp −u t , ψ p (r) = χ D(p) (r t ) .

Their Fenchel conjugates are
For p = 0, we obtain φ 0 (−u) = 0 , ψ 0 (r) = χ K0 (r) . Their Fenchel conjugates are . As an application of Theorem 5.1 we obtain that b p = max {b 1,p , b 2,p } is a bipotential. Indeed, we shall check only the point (ii') from Theorem (5.1) (the point (ii") is true by a similar computation). For λ ∈ [0, 1) and p = 0 we have: Also, by computation we obtain: < +∞ then in particular r n = p and we obtain (10) as an equality 0 = 0. All other cases, involving λ = 1 or p = 0 are solved in the same way.
8. Unilateral contact with orthotropic friction. For many industrial applications, the assumption of isotropic friction is unrealistic because of directional surface machining and finishing operations. For orthotropic frictional contact, Micha lowski and Mróz pointed out the non-associated nature of the sliding rule [22] (see also Mróz and Stupkiewicz [24]). In this model, the convex friction cone is defined by where • µ denotes the elliptic norm with the Euclidean norm • and The classical isotropic Coulomb's friction condition is recovered by setting µ x = µ y = µ. The dual norm • * µ associated to (15) is given by Micha lowski-Mróz model is twice non associated because, as in isotropic friction, the normal velocity is absent when sliding occurs, but also because the tangential velocity is not normal to the elliptical level curves r t µ = constant (Figure 1). This additional lack of associativity can be described introducing a slip potential defined by g(r tx , r ty ) = r t p (18) in which r t p is given by It is convenient to introduce the sliding non-associativity matrix x y r t y r t x Figure 1. Friction condition and sliding rule.
The complete form of the frictional contact law involves three possible states, which are separating, contact with sticking, and contact with sliding. It can be described using two overlapped "if...then...else" statements if r n = 0 theṅ u n ≥ 0 ! separating elseif r ∈ int K µ theṅ u n = 0 andu t = 0 ! sticking else (r ∈ bd K µ and r n > 0) u n = 0 ,λ ≥ 0 and −u t =λ ∂g ∂r t =λ P −2 r t r t p ! sliding endif (22) where "int K µ " and "bd K µ " denote the interior and the boundary of K µ , respectively.
It was proved in [19] that the graph of the previous law is the critical set of the following bipotential: For the isotropic case, we recover the bipotential (13) as particular case. 9. Implicit predictor-corrector scheme. The time interval [0, T ], within which the loading history is defined, is partitioned into N sub-intervals of size t, not necessarily equal. For the sake of clarity, we focus our attention on the first time increment [t 0 , t 1 ]. The value of any quantity a at the begining (resp. at the end) of the step is denoted by a 0 (resp. a 1 ). The corresponding increment is a. In order to ensure convergence and stability requirements, the implicit scheme is considered. Taking into account the fact that b is positively homogeneous of order one with respect to the first argument, the complete contact law is satisfied at the end of each time step: (25) takes into account the gap −h 0 at the beginning of the time-step.
Unfortunately, unlike elastoplastic problems, the incremental bipotential is not differentiable. For this reason, the incremental bipotential is not directly used and a regularization technique is applied in order to replace the unpleasant inequation by an equivalent equation, simpler to solve. Thus the implicit scheme (24) leads to The variational inequality (26) is rewritten in an obvious manner as: where ρ is a positive number that need to be chosen in a suitable range to ensure convergence. The variational inequality (27) means that r 1 is the proximal point of the augmented forcer = r 1 − ρ u, with respect to the function r 1 → ρb c (− u, r 1 ) (see [23]) The solution of (28) can be obtained using the Usawa's algorithm which involves two steps: By substituting b c (− u, r 1 ) by its expression (25) in (27), the inequality can be rearranged to obtain: are the components of the modified augmented surface traction. The last inequality means that the reaction r at the end of the time step is the projection of the augmented surface traction onto the convex friction cone K µ . With respect to algorithm with separated treatment of the friction and the unilateral contact, the main advantage of the method is that only one predictor-corrector step is required for the discrete frictional contact problem. At each step, the iterative algorithm is the following: (1) A value u of the displacement increment being known, the new value of the contact stress r 1 at the end of the time step is obtained by computing the augmented stress τ i+1 and the proximal point r i+1 1 : (2) Next, the displacement increment is updated and the procedure is iteratively repeated. In the solution of the projection problem, three different situations must be considered according to the position of the prediction τ in R 3 . The first case corresponds to a prediction τ located in the cone K µ . Its projection is the prediction itself, i.e. r 1 = τ . The second one relates to a prediction situated in the cone K * µ , where its projection turns out to be the origin (r 1 = 0). In the last case, the prediction is neither in K µ nor in K * µ and the corrector step requires computing the projection of the prediction. The projection of a point onto a convex set is equivalent to the minimization of the distance between this point and the convex set. Next, the problem is reformulated as an unconstrained minimization problem by means of the Lagrange multipliers technique and leads to find the intersection of a quartic and a straight line. Details are given in [19]. 10. Numerical applications. In previous papers [11,12], several applications involving frictional contact problems with isotropic friction condition and associated sliding rule have been carried out using an algorithm based on the bipotential approach. The examples treated have shown that the algorithm is very competitive as the augmentation phase involves only one prediction-correction step.
However, although a lot of works were devoted to the isotropic friction, this hypothesis is not often realistic. In fact, most of frictional contact are anisotropic.The source of the roughness anisotropy is technological; the industrial process used to fabricate the bodies can create striations along preferential directions. In fact, most machining, finishing and superfinishing operations are directional, and machined surfaces have particular striation patterns unique to type of machining. Also specific techniques of manufacture produce a surface with anisotropic frictional properties. For a large number of machining processes, the striation directions are mutually orthogonal. For such surfaces, an orthotropic friction condition will provide a better description of the frictional behavior.
In [20], we proposed a benchmark test to validate the algorithm for a class of nonassociated anisotropic friction laws. The test of such frictional contact laws requires a 3D finite element model. The problem under consideration is a deformable elastic cylinder in contact with a rigid surface (Figure 2). The radius and the height of the cylinder are both equal to 10 mm. The Young modulus E of the cylinder is taken equal to 210000 MPa and the Poisson ratio is 0.3. On the surface contact, the friction condition is assumed to be anisotropic. A vertical rigid motion is imposed on the upper surface of the cylinder by an amount of 0.1 mm. The displacement is applied in one step. The base of the cylinder is in contact with the rigid plate whose normal vector is (0, 0, 1). The cylinder is subdivided into 1280 eight-node brick-like elements as shown in Figure 2. Each element has 27 integration points.
: Y Z Figure 2. Compression of a cylinder in contact with a rigid plate.
Five different sets of frictional parameters, shown in Table 1, are considered.  Table 1. Frictional properties The first case corresponds to a classical isotropic friction condition, considered here for comparison with anisotropic cases. Case 2 and case 3 represent an anisotropic frictional model with an associated sliding rule. The anisotropy is mild in case 2. The last two cases consider a non-associated sliding rule. Figure 3 shows the contour plots of the slip and the relative displacements between the lower surface of the cylinder and the rigid plate in the x-direction. The slip correspond to the Euclidean norm of the tangential displacement. The frictional model being isotropic, the iso-values of the slips are circular. As expected, the stick zone is located around the base center and sliding increases monotonically as we get closer the periphery. The magnitude of the tangential displacement in the x-direction is the largest on the cylinder edge.   In all figures, the x-axis is horizontal and the y-axis is vertical. Since the model is isotropic, a simple rotation about the z-axis of 90 0 gives the the tangential component of the contact reaction in the y-direction. The normal reaction is higher in the periphery of the cylinder and decreases as we approach the center of the cylinder basis. However, the decrease is not monotonic along every radius. Indeed, the normal reaction attains its minimum in four tiny areas located at approximately the third of the radius from the base center. For all cases considered, contour plots of the normal component of the contact reaction have similar pattern. The tangential reaction is higher in the periphery of the cylinder as sliding is more important there. Now, the effect of anisotropy is considered. Slip contour plots for H n 4 9 5 . 8 9 3 3 8 5 . 6 9 5 2 7 5 . 4 9 6 1 6 5 . 2 9 8 5 5 . 0 9 9 3 5 5 . 0 9 9 3 1 6 5 .   criterion. The slips increases gradually from the stick area and are maximum on the periphery in the y-direction since for both cases µ x is greather then µ y . If the disparity between the friction coefficients µ x and µ y is significant (as it is for case 3) the distribution of r ty changes notably, but the iso-values pattern of r tx remains similar to the isotropic case 1 (see Figure 6). The algorithm is still convergent if the sliding rule is non-associated but requires a few more iteration. The number of iterations depends strongly on the degree of non-associativity. A ratio µ x µ y much larger or much smaller than the ratio p x p y gives a strong non-associated sliding rule.
In Case 4, an isotropic sliding potential is considered which corresponds here to a mild non-associativity. In practice, high non-associativity (case 5) is often observed ( [22,24]). Once a non-associative sliding rule is considered, the slip distribution can change drastically according to the degree of non-associativity. Indeed, if the slip rule is strongly non-associated, the iso-values of the slip become non-convex as shown Figure 7. With an isotropic sliding potential, the stick zone is elliptical and the maximum displacement is lower than the one obtained with an associated sliding rule. Although, the principal friction coefficients are the same for both cases 4 and 5, the maximum slip is larger for case 5. The distribution pattern of r ty (Figure 8) is similar to cases 2 and 3 (see Figure 5.a) but the iso-values distribution of r tx is totally different if the non-associativity is strong. Moreover, it is worth to observe that, unlike the isotropic case, there is a strong hysteretic behavior occurring when friction is anisotropic [15]. This effect is particularly important because the wear rate is strongly coupled to the friction dissipation.   Figure 9 shows the relationship between the equivalent applied force (resulting from imposed displacement) and the slip of point P at the intersection of the the boundary of the contact zone and the bisectrix of the quadrant Oxy. It can be seen that the relationship is linear during loading. The unloading stage starts with a sharp drop in the applied force followed by a non-linear decrease. The difference between loading path and unloading path appears to be evident. Further, there is a slight difference between the trajectories during the first cycle and the second one. Figure 9. Total force versus slip at point P

11.
Conclusions. In the theoretical part, we showed that the bipotential of Coulomb's friction law is related to a specific bipotential convex cover with the property that any graph of the cover is non maximal cyclically monotone. On this ground, we proposed a general algorithm to construct explicitly the bipotential for the modelling of a given constitutive law.
In the numerical part, a robust algorithm, developed by de Saxcé and Feng [12] for the frictional contact law with associated sliding rule, has been adapted to handle non-associated sliding rules occuring when the friction is orthotropic. This algorithm has been successfully tested on examples including cyclic loading for which ones a strong hysteretic behavior has been demonstrated..