Philippe Bussetta
email: p.bussetta@ulg.ac.be

Daniel Marceau
email: daniel_marceau@uqac.ca

Jean-Philippe Ponthot
email: jp.ponthot@ulg.ac.be

J.-P Ponthot

Ltas Mn

The adapted augmented Lagrangian method: a new method for the resolution of the mechanical frictional contact problem

Keywords:

may come

The adapted augmented Lagrangian method: a new method for the resolution of the mechanical frictional contact problem

good resolution of the frictional contact problem is disrupted by the nonlinearity and the non differentiability of the equations of contact (stick/unstick and the initiation of sliding) [START_REF] Wriggers | Computational contact mechanics[END_REF]. This paper deals with solving the mechanical contact problem, using the finite element method, in the general case of thermo-electro-mechanical problem (in two or three dimensions). Until now, no method can be used to solve all contact problems. The solving method can be divided in two parts: the formulation of the contact equations and the solving of these equations using the finite element method. The methods most frequently used to obtain equations system are the penalty method (PM) [START_REF] Sh | Contact friction simulation in powder compaction process based on the penalty approach[END_REF][START_REF] Bayram | Enriched finite element-penalty function method for modeling interface cracks with contact[END_REF][START_REF] Zavarise | Real contact mechanisms and finite element formulation-a coupled thermomechanical approach[END_REF][START_REF] Krstulović-Opara | A c1-continuous formulation for 3d finite deformation frictional contact[END_REF][START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF] or augmented Lagrangian method (ALM) [START_REF] Elguedj | A mixed augmented Lagrangian-extended finite element method for modelling elasticplastic fatigue crack growth with unilateral contact[END_REF][START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF][START_REF] Wriggers | Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces[END_REF][START_REF] Baba | An augmented Lagrangian treatment of the metal forming process[END_REF], or hybrid methods [START_REF] Wriggers | Finite element algorithms for contact problems[END_REF][START_REF] Chen | An effective two-dimensional frictional contact model for arbitrary curved geometry[END_REF][START_REF] Walter | Modélisation 3D par éléments finis du contact avec frottement et de l'endommagement du béton: application à l'étude de fixation ancrées dans une structure en béton[END_REF][START_REF] Carpenter | Lagrange constraints for transient finite elements surface contact[END_REF][START_REF] Oldenburg | The position code algorithm for contact searching[END_REF]. But, alternative methods can be used, like the mathematical programming [START_REF] Klarbring | Shakedown in elastic contact problems with coulomb[END_REF][START_REF] Ireman | Finite element algorithms for thermoelastic wear problems[END_REF]. These methods are simple but sometimes very difficult to use in practice, because of the choice of the value of penalty parameters (normal and tangential). The more the value of the penalty coefficients is low, the more the error on contact stress or the computational time is important. Nevertheless, the more the value of the penalty coefficient is large, the more the risks of numerical oscillations and divergence of the algorithm are important. To overcome this lack of performances, an adapted penalty method (APM) was proposed by Chamoret [START_REF] Chamoret | New smoothing procedures in mechanics[END_REF] or [START_REF] Chamoret | Modélisation du contact: nouvelles approches numériques[END_REF]. But this method deals only with normal penalty coefficient. The value of the tangential penalty coefficient is still user chosen. In addition, this method is not appropriate to elastoplastic problems [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF] (it is proved in the Sect. [START_REF] Zavarise | Real contact mechanisms and finite element formulation-a coupled thermomechanical approach[END_REF]

by one example).

This paper presents a new method for the resolution of the mechanical contact problem, called "adapted augmented Lagrangian method" (AALM). This new method is based on the augmented Lagrangian method coupled with an adaptation of both penalty coefficients (normal and tangential). This one has the speed of the PM and the reliability of the ALM. In opposition to classical methods (PM, ALM), the user does not choose the value of the penalty parameters. With this new method, like with the classical methods (PM, ALM), the contact equations can be solved thanks to a large a variety of techniques: like a node-to-segment algorithm [START_REF] Zavarise | Real contact mechanisms and finite element formulation-a coupled thermomechanical approach[END_REF][START_REF] Chamoret | New smoothing procedures in mechanics[END_REF][START_REF] Chamoret | Modélisation du contact: nouvelles approches numériques[END_REF][START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF][START_REF] Zavarise | The node-to-segment algorithm for 2d frictionless contact: classical formulation and special cases[END_REF], a segment-to-segment algorithm [START_REF] Zavarise | A segment-to-segment contact strategy[END_REF][START_REF] Wriggers | A formulation for frictionless contact problems using a weak form introduced by Nitsche[END_REF], the mortar element method [START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Puso | A mortar segment-to-segment frictional contact method for large deformations[END_REF], the boundary element method [START_REF] Oysu | Coupled fem-bem for elastoplastic contact problems using Lagrange multipliers[END_REF][START_REF] Chernov | A priori error estimates for hp penalty bem for contact problems in elasticity[END_REF][START_REF] Yamazaki | Penalty method for three-dimensional elastic contact problems by element method[END_REF][START_REF] Rodríguez-Tembleque | 3D frictional contact of anisotropic solids using bem[END_REF], or the X-FEM [START_REF] Bayram | Enriched finite element-penalty function method for modeling interface cracks with contact[END_REF][START_REF] Elguedj | A mixed augmented Lagrangian-extended finite element method for modelling elasticplastic fatigue crack growth with unilateral contact[END_REF][32][START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF][START_REF] Legrain | Stress analysis around crack tips in finite strain problems using the extended finite element method[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. But, in this paper only a segmentto-segment algorithm is used (see Sect. 3). The proposed method is first tested on academic problems involving elastic as well as elastoplastic behaviour. The adapted augmented Lagrangian method is then used for an industrial problem with mechanical frictional contact: the electrical preheating of a cathode block used for aluminium production.

Equations to be solved

This paper deals with the problem of mechanical contact between two solids in three dimensions in the quasi-static regime. However, the methods used can easily be applied in the case of multiple bodies and multiple interfaces of contact or/and in the dynamic regime. The two solids are potentially in contact with the boundary Γ 1 c (solid 1) and Γ 2 c (solid 2). The equations that represent the mechanical equilibrium can be written using the principle of virtual work such that:

W int (u, δu) + W ext (u, δu) + W c (u, δu) = 0. (1)
Where u and δu are respectively the real and the virtual displacement field. Moreover, W int (u, δu) is the virtual work associated to the internal forces and W ext (u, δu) is the virtual work related to the external forces except the contact ones.

The virtual work connected to the contact is noted W c (u, δu). This last one can be written as:

W c (u, δu) = Γ 1 c δu • t c dΓ 1 c + Γ 2 c δu • t c dΓ 2 c ( 2
)
where t c is the surface traction on the contact boundary (Γ 1 c or Γ 2 c ). This surface traction is called contact stress. Using the law of reciprocal actions (Newton's third law), the virtual work of the contact forces (W c ) can be written as one integral over the contact surface Γ 1 c , called the slave surface, as [START_REF] Wriggers | Computational contact mechanics[END_REF] or [START_REF] Marceau | Modélisation du contact tridimensionnel avec frottement en grande transformation et son application á l'étude des dispositifs d'ancrage multitorons[END_REF]:

W c (u, δu) = Γ 1 c δξ β t t β + δg t n dΓ 1 c (3)
In this paper the Greek letters (α and β) can take the value 1 or 2. On the contact boundary Γ 1 c , the field of virtual relative displacement between the two interfaces is split in the normal direction, δg and on the tangential plane, δξ β (for more information see [START_REF] Wriggers | Computational contact mechanics[END_REF] or [START_REF] Marceau | Modélisation du contact tridimensionnel avec frottement en grande transformation et son application á l'étude des dispositifs d'ancrage multitorons[END_REF]). In the same way, the contact stress is projected on the normal direction (normal contact stress, t n ) and on the tangential plane (tangential contact stress, t t = t t 1 + t t 2 ) [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF]. The contact's laws are the unilateral contact (relations (4)) and the Coulomb's friction law (relations ( 5)), such that:

g ≤ 0 t n ≥ 0 gt n = 0 ( 4 ) Φ = t t -μt n ≤ 0 ζ ≥ 0 Φζ = 0 ( 5
)
where ζ is the absolute rate of sliding, g is the distance between the two boundaries and μ is the Coulomb's friction coefficient. The relation between the rate of sliding and the speed of sliding (v t ) is given by:

v t -ζ t t t t = 0. ( 6
)

Solution techniques

This problem is solved using the finite element method. The computation of the integral of the virtual work of the contact forces over the contact surface Γ 1 c is done by numerical integration. The contact stress is computed thanks to a segment-to-segment algorithm [START_REF] Marceau | Modélisation du contact tridimensionnel avec frottement en grande transformation et son application á l'étude des dispositifs d'ancrage multitorons[END_REF]. In opposition to the node-to-segment algorithm (see in [START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF][START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF][START_REF] Zavarise | The node-to-segment algorithm for 2d frictionless contact: classical formulation and special cases[END_REF][START_REF] Wriggers | Finite element formulation of large deformation impact contact problem with friction[END_REF]), the contact condition is not verified for the nodes of the slave boundary. The equations of contact are solved for every integration point of each element on the contact surface Γ 1 c . In this algorithm, each integration point is similar to a virtual node in the node-to-segment algorithm. So, the contact stresses are known through their values at integration points. The link between one point of the slave surface (Γ 1 c ) and the master surface (Γ 2 c ) is done by orthogonal projection on the closer element of the boundary Γ 2 c (see Fig. 1). The search of the link between a point of the slave surface and an element of the master surface, can be divided into three cases. In the first one, the point can only be projected on one element of the slave boundary (see Fig. 1). In the second case, the point of the slave surface can be projected on two elements of the master surface. The matching element of the master surface is the one corresponding to the smaller gap (see Fig. 2). In the last case, the point of the slave surface can not be projected on any element of the master surface. The tolerance on the boundary of the element can by-pass this case (see Fig. 3). So, in all cases, the contact stress is computed with the matching element of the master surface. This projection is used for all integration points of each element of the contact surface Γ 1 c to compute the virtual work of the contact forces related to this element.

This computation is available for any number of integration points and any integration scheme (Gauss, Newton-Cotes, etc.; see Fig. 1 in a two dimensional problem or Fig. 4 in a three dimensional problem). More details on the formulation of the mechanical contact are available on [START_REF] Marceau | Modélisation du contact tridimensionnel avec frottement en grande transformation et son application á l'étude des dispositifs d'ancrage multitorons[END_REF] or on [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problem[END_REF].

Some works deal with the numerical integration scheme of the virtual work of the contact forces. Investigations about integration rules for the PM are presented in [START_REF] Qiu | Stiffness matrix integration rules for contact-friction finite elements[END_REF]. Another study [START_REF] Lei | Contact friction analysis with a simple interface element[END_REF], deals with the comparison of the Gauss integration scheme and the Newton-Cotes integration scheme. This research shows that the Gauss scheme is a possible source for numerical oscillation of the contact stress, and the use of the Newton-Cotes scheme improves the reliability of the contact stress. On the other hand, a composite integration scheme is proposed in [START_REF] El-Abbasi | Stability and patch test performance of contact discretizations and a new solution algorithm[END_REF] and used in [START_REF] Harnau | Algorithmic aspects in large deformation contact analysis using 'solid-shell' elements[END_REF]. With this g t n t n t t t n Fig. 5 Contact stress versus the penetration and the tangential displacement scheme, each element of the slave boundary is divided into sub-elements (intersection of this element with an element of the master boundary). In [START_REF] El-Abbasi | Stability and patch test performance of contact discretizations and a new solution algorithm[END_REF][START_REF] Harnau | Algorithmic aspects in large deformation contact analysis using 'solid-shell' elements[END_REF], the numerical integration can be performed on each sub-element using either the Gauss integration scheme or the Newton-Cotes integration scheme. However, in the present paper only the Newton-Cotes integration scheme has actually been employed on each element of the contact surface Γ 1 c .

Penalty method (PM)

The PM is the most used method to solve this kind of problems. This method ensures a continuous relation between the contact stress and the displacements. Therefore, two user defined coefficients are inserted, the normal penalty coefficient (ε n ) and the tangential penalty coefficient (ε t ).

The relations used to define the contact stress are given by:

t n = ε n < g > v t -ζ t t t t = 1 ε t L v t t (7)
with the Coulomb's friction law (relations (5)). The operator L v means the Lie's derivative and < > the positive part of the operand (Macauley brackets). The value of the contact stress in function of the penetration and the tangential displacement is shown in the Fig. 5. The displacement field (u) is solved using the Newton-Raphson solution technique in an implicit context. So, this method introduces at the contact interface, a residual penetration and a tangential displacement without sliding (called reversible tangential displacement). The choice of the penalty coefficients is of the utmost importance in order to get an effective solution in an efficient way. The more the value of these coefficients is small, the less the solution has physical meaning because of the penetration and the reversible tangential displacement at the contact interface. On the other hand, a large value of these coefficients may produce numerical oscillations and generally prevents the convergence of the algorithm. To overcome the problem of the choice of the values of penalty coefficients, one APM was proposed by Chamoret [START_REF] Chamoret | New smoothing procedures in mechanics[END_REF] or [START_REF] Chamoret | Modélisation du contact: nouvelles approches numériques[END_REF]. This method introduces two new parameters: the minimal and the maximal value of the gap (g min and g max ). These new parameters have more physical meaning (they can be connected to the roughness of the two surfaces), but, this adaptation is only about the normal penalty coefficient, and the initial value of ε n is still user chosen.

The value of the normal penalty coefficient is updated as a function of the value of the penetration between the two solids and the gap tolerance such that:

ε n i+1 = F (|g i |, g min , g max )ε n i (8)
where the function F is defined by the relation (9) see also Fig. 6.

F (|g i |, g min , g max ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ |g i | g max if |g i | > g max |g i | g min if |g i | < g min 1 else (9)
This method was used to solve mechanical contact problems without friction [START_REF] Chamoret | New smoothing procedures in mechanics[END_REF] or [START_REF] Chamoret | Modélisation du contact: nouvelles approches numériques[END_REF]. Actually, APM is not adapted for elastoplastical problems because, the important augmentation of the value of the normal penalty coefficient produce numerical oscillation and generally prevent the convergence of the algorithm [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF].

Augmented Lagrangian method (ALM)

The ALM introduces additional parameters called Lagrange multipliers (λ). These parameters are used to reduce the impact of the penalty coefficients. Furthermore, the value of tolerance on the gap (g max ) and the reversible tangential displacement (ξ max ) are introduced to reduce the error on the contact stress.

The technique used to solve the mechanical problem is almost the same than with the PM. But the relations that define the contact stress are given by:

t n =< λ n + ε n g > v t -ζ t t t t = 1 ε t L v t t -L v λ t (10)
with the Coulomb's friction law (relations ( 5)).

The Lagrange multipliers (λ) are split into two parts, a normal part (λ n ) and a tangential part (λ t ). The initial value of the Lagrange multipliers is the value of the contact stress at the end of the lastest converged time step.

After the convergence of the Newton-Raphson solution technique, if the value of the gap or the reversible tangential displacement is higher than the prescribed tolerance, the value of the Lagrange multipliers is updated using the relations [START_REF] Wriggers | Finite element algorithms for contact problems[END_REF]. This updating is called an augmentation. The equations are then solved with this new value of the Lagrange multipliers (see the Algorithm 1).

λ k+1 n = t k+1 n λ k+1 t = t k+1 t (11)

Algorithm 1 Computation of Lagrange multipliers

Require: ε n , g i and ξ i .

1. The problem is solved by the Newton-Raphson algorithm if for all integration points g i < g max and ξ i < ξ max then The problem is solved else Augmentation of the Lagrange multipliers (with the relations ( 11)) Go back to 1 end if

The augmentation of the Lagrange multipliers allows the use of smaller value of the penalty coefficient. Nevertheless, the smaller the value of the penalty coefficients, the lower the speed of convergence is and the computational time is thus increased. Like for the PM, large values of the penalty coefficients prevent convergence of the algorithm. In addition, the quality of the solution can be controlled by the tolerance on the gap, the tolerance on the reversible tangential displacement and the number of augmentation of the Lagrange multipliers [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF]. So, the larger value of the penalty coefficients, the smaller the number of augmentation of the Lagrange multipliers will be. But, in the same time, the error on the contact stress increases. Because the incremental value of the contact stress on the last augmentation of the Lagrange multipliers is not negligible in comparison to the value of the Lagrange multiplier (the value of the contact stress is not converged; it is proved in the Sect. 4 by one example). So, the choice of the value of the penalty coefficients remains the key part of the success of both, the ALM as well as the PM.

Adapted augmented Lagrangian method (AALM)

The AALM is based on ALM enhanced with an adaptation of the penalty parameters (normal and tangential). The resolution of the mechanical problem is the same as for the ALM (based on the Newton-Raphson solution technique). Like the ALM, the augmentation of the Lagrange multipliers is done after the convergence of the Newton-Raphson algorithm (see the Algorithm 1). However, the value of the penalty coefficients is updated during the computation. A specific adaptation is used with the penalty coefficients (normal and tangential).

The normal penalty coefficient is adapted at each iteration and at each contact point. In the same way, the computation of the contact stress is different from ALM, as described in Algorithm 2.

Algorithm 2 Computation of normal contact stress

Require: ε n , g i , g i-1 , λ n and contact status at the iteration i -1.

if no contact at the iteration i -

1 then Evalt n = 0 ε n = initial value of ε n else if g i × g i-1 < 0 then if the sign of g changes two times in a row then Evalt n = 0 else Evalt n = ε n g i-1 end if end if ε n = adaptation of ε n (see Algorithm 3) end if t n =< Evalt n + ε n g i + λ n >
The same relations that with the ALM are used to calculate the value of normal contact stress when the contact takes place. In the other cases, when the point is in contact at the previous iteration, the value of the normal contact stress is a function of contact history (see Algorithm 2).

This procedure is based on the following assumption: if the sign of the gap changes (g i × g i-1 < 0), this means that the last correction (g i-1 ×ε n ) was too large, or the new correction (g i × ε n ) has a different sign. In the general case, the value of the normal stress lies between the value of Lagrange multiplier (λ n ) and the value of the last iteration (ε n g i-1 +λ n ). So, the last correction is added to the normal contact stress value (Evalt n = g i-1 ε n ) to reduce the numerical oscillations (see Algorithm 2).

On the other hand, if the sign of the gap changes two times in a row, no history correction is added (Evalt n = 0), because the oscillation of the gap sign is touched off by a too important value of the normal penalty coefficient. So the decrease of the normal penalty coefficient stabilizes the gap sign (see Algorithms 2 and 3).

In addition, the value of the normal penalty coefficient is adapted according to the contact history (excepted when the contact takes place, see Algorithms 2 and 3).

Algorithm 3 Adaptation of normal penalty coefficient

Require: ε n , g i and g i- [START_REF] Wriggers | Computational contact mechanics[END_REF] .

if g i × g i-1 < 0 then if g i-1 > g max then ε n = |(ε n g i-1 )/g i × (|g i | + g max )/(g i -g i-1 )| else ε n = |ε n g i-1 /(10g i )| end if else if g i > g max then if |g i -g i-1 | > max(g i /10; g i-1 /10; 5g max ) then ε n = 2ε n else if |g i | = |g i-1 | ± 1% < 10g max then ε n = ε n ( √ (|g i |/g max -1) + 1) 2 else if |g i | > |g i-1 | + 1% then ε n = 2ε n (g i-1 /g i ) else ε n = ε n (( √ (|g i |/g max -1) + 1)) end if else ε n = ε n end if
The adaptation of the normal penalty is separated in three cases: either the sign of the gap changes (g i × g i-1 < 0), or the absolute value of the gap is more important than the prescribed limit (|g i | > g max ), or the absolute value of the gap is less than the maximal limit (|g i | < g max ).

For the last case, the normal penalty coefficient remains unchanged.

For the first case, when g i × g i-1 is negative, the normal penalty coefficient is adapted so that the absolute value of the gap decreases without changing the sign of the gap.

For the second case, when |g i | is greater than g max , the adaptation of ε n will increase the absolute value of the penalty contribution of the normal stress without changing the sign of the gap. In the context of optimizing the computational time, several cases are considered.

• In the case that the difference with the gap of the last iteration is too large, the normal penalty coefficient is multiplied by two to limit the increase of the absolute value of the penalty contribution (and thus reduce the numerical oscillations; see the Fig. 7). • In the case that the gap is quasi-constant, the normal penalty coefficient is adapted twice (the multiplier of the coefficient is squared) to increase the rate of convergence (see the Fig. 8). • In the case that the absolute value of the gap increases, the adaptation of ε n multiplies by two the penalty contribution in the normal stress (see the Fig. 9). 8) and ( 9)), which reduces the numerical oscillations.

This adaptation let the algorithm converge with an wide range of the initial value of the normal penalty coefficient.

The value of the tangential penalty coefficient is based on an empirical study that shows that the overvaluation of the tangential stress has a bad influence on the resolution of the problem [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF]. When the contact takes place, the value of the tangential penalty coefficient leads to a tangential stress equal to the half of the stress value corresponding to sliding threshold. After that, the value of the tangential penalty coefficient is adapted during the Lagrange multipliers augmentation (see Algorithm 4). The case of the reversal sliding way means that the tangential penalty coefficient was too important so the new value is equal to the half of the old one.

Algorithm 4 Evaluation of the value of the tangential penalty coefficient

Require: ε t , GlRe = Absolute value of reversible sliding.

if (no contact at the iteration i -1) or (first iteration after Lagrangian augmentation and GlRe < Tol) then ε t = μt n /(2 GlRe ) else if reversal sliding way then

ε t = ε t /2 end if
The condition of unchanged penalty coefficients is added to the convergence of the Newton-Raphson solution technique. Like the ALM, the augmentation of the Lagrange multipliers is performed after the convergence of the Newton-Raphson solution technique. At the time of the Lagrange multipliers augmentation, the tolerance on the gap in the function of normal penalty adaptation is divided by ten to increase the convergence rate.

Thereby, the user can prescribe the maximal value of the gap in the solution. So, the algorithm is determined by two of these three parameters: the number of Lagrange multipliers augmentation, the gap tolerance for the solution or the gap tolerance for the first Lagrange multipliers augmentation (the tolerance on the gap in the function of normal penalty adaptation before the first augmentation of the Lagrange multipliers).

Numerical examples

In this section, numerical examples are presented to demonstrate the efficiency of the proposed adapted augmented Lagrangian method (AALM) and the advantages over the most frequently used methods (PM and ALM). At first, some academic problems are presented. The two first academic examples have already been presented in [START_REF] Bussetta | Résolution du problème de contact mécanique frottant: méthode du lagrangien augmenté adapté[END_REF] (the Hertz problem and an elastoplastic problem). In addition, one academic frictional contact problem is shown in this section. Then, an original industrial problem, a thermoelectro-mechanical problem with frictional contact is treated. In all problems, the linear elements are used and the computation of the integral over the contact surface is done with two Newton-Cotes integration points in each direction (two integration points for the line and four for the quadrangle). With this integration scheme, the integration points correspond to the nodes of the finite element mesh.

Hertz problem

The Hertz problem is frequently used to check the mechanical contact algorithm. It is a classical problem of mechanical contact between two parallel, infinitely long cylinders compressed one against the other using concentrated forces (see Fig. 10).

The Hertz contact problem without friction between two identical cylinders is solved. The analytical solution can be found in [START_REF] Timoshenko | Theory of elasticity[END_REF]. This problem is solved in two dimensions with the assumption of plane strain state.

The radius of the cylinders is 0.25 m. The value of the Young's modulus and the Poisson's ratio of the two cylinders are respectively 200 GPa and 0.3. In addition, this problem is frictionless. Considering the symmetry and the hypothesis of the small deformations, the simulation is done with only two quarters of cylinder (Fig. 11). This problem could be solved with only one quarter of cylinder and one rigid plan, but with this configuration, the computation of the contact is too easy (only one solid in contact is deformable and nodes are simply projected into an infinite plane). So, in the aim of testing the reliability of the AALM and its robustness, this simplification is not used (the two solids in contact are deformable; see [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF]). This problem is solved in two equal steps, the final value of δ is 0.11854 mm (corresponding to the analytical solution with a maximum value of the normal contact stress of 1,000 GPa).

The two quarters of cylinder are identically meshed with triangular elements (see Fig. 12). Three meshes are used (noted mesh n where n means the number of elements in the contact zone, n = 2, 4 or 6).

Table 1 shows that the PM is faster but the value of normal penalty coefficient has a strong influence on the final solution (see Fig. 13a). A too small value of this coefficient (10 5 ) results a bad evaluation of the contact stress (20% error, see Fig. 13a). With this value of the penalty coefficient, the maximal value of the gap is about 8×10 -3 mm. On the other hand, a too large value of this coefficient (10 7 ) can cause numerical oscillations and the divergence of the algorithm (e.g. mesh 2, see Table 1). In addition, the approximation of the contact stress is not better (with mesh 4, see Fig. 13a, ε n = 10 7 ) and the maximal value of the gap is about 10 -4 mm. With the APM the computational time is more important than with the PM (see Table 1). The dependency between the computational time and the normal penalty coefficient is not clear. In addition, the value of the normal stress is not closer to the analytical solution than with the PM (see Fig. 13a,b), but the error is independent of the value of the normal penalty coefficient (2 ≤ ε 0 n ≤ 2 × 10 4 ). With the ALM the computational time is equivalent to the PM but it increases quickly with the decreases of the normal penalty coefficient (see Table 1). Moreover, when the normal penalty coefficient increases, the quality of the solution decreases (see Fig. 13c, ε n = 10 6 ). This fact is explained by the hardening of the interface. Because, the more the normal penalty coefficient is important, the more a little difference on the gap has a great influence on the normal contact stress (ε n g). So to keep the error on the normal contact stress, the more the normal penalty coefficient is important, the more the tolerance on the gap is small.

For this example, the AALM needs a little more of iteration than the PM (see Table 1), but with the AALM the initial value of the penalty coefficient (2 ≤ ε 0 n ≤ 2 × 10 4 ) does not have any influence on the solution accuracy while the influence of the computational time is quite limited (see Fig. 13d). The convergence of the contact normal stress (with the AALM and ε 0 n = 2 × 10 2 ) as a function of the size of the elements is proved in Fig. 14.

Figure 15 shows that the value of the convergence rate is the same with the classical methods (PM, ALM) and the AALM.

Penalty Method

Adapted Penalty Method Augmented Lagrangian Method (d) Adapted Augmented Lagrangian Method 13 Contact normal stress on Hertz problem (mesh 4) So, this example shows that with the PM and the ALM, the choice of the normal penalty coefficient is of the utmost importance in order to get a reliable solution. On the other hand, thanks to the penalty adaptation, the APM and the AALM are able to solve this problem with a ited influence of the initial value of the normal penalty coefficient.

Elastoplastic problem

In this problem, the same two dimensional geometry, initial conditions, boundary conditions and, meshes as in the Hertzian problem are used. But, it is a problem of mechanical contact between two infinitely thin cylinders compressed one against the other using concentrated forces (see Fig. 10). So, the hypothesis of plane stress state is done. In addition, the triangular element is not the best way to compute an elastoplastic problem, but, the aim of this example is to test the contact computational methods when the solids exhibit elastoplastic behaviour. The mechanical properties of the material are the following:

-yield stress: σ e = 472 MPa Y Z Fig. 16 von Mises stress distribution (MPa, mesh 6, AALM with

ε 0 n = 2 × 10 4 )
-isotropic hardening law: σ 0 (MPa) = σ e + K ε p with K = 640 MPa and ε p is the equivalent plastic strain.

Like the Hertz problem, this problem is solved in two equal steps and the final value of δ is the same (0.11854 mm).

The stresses are concentred at the contact interface (see Fig. 16). So, the displacements of the contact boundaries and the variation of the gap between two iterations are larger. Consequently, the contact computation is more challenging.

Just as for the Hertz problem, with the PM, the value of the normal penalty coefficient is very important. But the optimal value of the normal penalty coefficient used in the Hertz problem is not appropriate in this problem (10 6 see Table 2). With this value, this problem cannot be solved with mesh 6 because numerical oscillations of contact status appear.

This example shows the limit of the APM, this method is not adapted to solve this problem in two steps (see Table 2 and Fig. 17b), because the problem cannot be solved with APM for ε 0 n equal to 2×10 4 and for ε 0 n equal to 2×10 2 the solution is obtained only with meshes 2 and 4. So, the suitable range of the initial value for the normal penalty coefficient cannot be determined. This underachievement is induced by the very important increase of the value of the normal penalty coefficient and the perturbation caused by the change of mechanical contact status (contact/no contact).

Like in the Hertz problem, with the ALM the value of the normal penalty coefficient is very important. Like with the PM, the optimal value of the normal penalty coefficient (10 4 ) with the ALM, is not the same as the optimal value of the Hertz problem (10 5 ). In addition, the solution is very reliable but the computational time is important.

Also, with AALM, the solution accuracy is independent of the initial value of normal penalty coefficient (see Fig. 17d). In addition, the influence of this initial value on the calculation time is weak (2 ≤ ε 0 n ≤ 2 × 10 4 , see Table 2). The convergence of the contact normal stress (with the AALM and ε 0 n = 2 × 10 2 ) as a function of the size of the elements is shown in Fig. 18.

Finally, as observed for the Hertz problem, this example shows that with the PM and the ALM the choice of the normal penalty coefficient is of the utmost importance. In addition, this example gives an evidence of the limitations of the APM. The adaptation of the normal penalty coefficient for this method is not able to solve this kind of elastoplastic problem. So, using a better penalty adaptation, the AALM is the only method that can solve this problem quickly for any mesh (mesh 2, 4 and 6) and any initial value of the normal penalty coefficient (2 ≤ ε 0 n ≤ 2 × 10 4 ).

Frictional contact problem

This academic problem treats of a infinitely long nylon corner moving on a steel slide with mechanical frictional contact (see Fig. 19). This example is solved in two dimensions with the assumption of plane strain state. Considering the symmetry, only the half of the problem is solved. The mesh of the slide is composed of 5 × 20 elements and the mesh of the corner is made of 20 × 20 elements (see Fig. 20). The dimensions and the boundary conditions are defined on the Fig. 20. The downward displacement on the top of the corner surface is imposed to 4.9 mm. This problem is solved in [START_REF] Mirchi | Comparative characterisation of graphitized and graphitic cathode blocks[END_REF] In addition, the Coulomb's friction coefficient is set equal to 0.15. Table 3 shows that the fastest method is the PM. On the other hand, with the PM the value of the penalty coefficients is very important. A too small value of the normal penalty coefficient (10 3 or 10 4 ) results in a bad quality of the solution (see Figs. 21 and22). With these value of the normal penalty coefficient, the difference between the solution and the reference solution (ε n = 10 5 see Fig. 23) is about 4 or 35%. In addition, a too large value of this coefficient (10 6 ) introduces numerical oscillation of the contact status that prevents the convergence of the algorithm. Like the value of the normal penalty coefficient, the value of the Fig. [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF] Nylon corner (A) moving along a steel slide (B) and boundary conditions penalty coefficient is very important. With a too small or a too large value of the tangential penalty coefficient, this method cannot solve this problem (ε t = 10 or 10 4 with ε n = 10 5 ). Nevertheless, with the appropriate value of the penalty coefficients (ε n = 10 5 and ε t = 10 2 or 10 3 ) the quality of the solution is good (see Fig. 23). In addition, the value of the normal penalty coefficient used in this problem is too small to be used in the previous problems. Like in the previous problems, with the ALM, the value of the penalty coefficient is very important. The computational time increases when the value of the normal penalty coefficient decreases (see Table 3). Moreover, a too small value of the tangential penalty coefficient (ε t = 1) causes numerical oscillations preventing the convergence of the algorithm. Furthermore, when the ALM can find an equilibrated solution, the accuracy of the solution is equal to the one of the reference solution with the PM (see Fig. 23).

To solve this problem with the AALM, just an initial value of the normal penalty coefficient is required. This method can solve this problem with a large value of the initial normal penalty coefficient (1 ≤ ε 0 n ≤ 10 2 ). In addition, the accuracy of the solution is equal to the one of the reference solution with the PM (see Fig. 23) or the one of the solution with the ALM. Table 3 shows that the AALM requires more iterations that the other ones (PM and ALM). At last, this example shows that with the PM and the ALM the choice of the value of the tangential penalty coefficient has the same importance that the one of the normal penalty coefficient. Besides, to solve a frictional contact problem with the AALM, only the initial value of the normal penalty coefficient is needed. However, this kind of problem can be calculated with the AALM with a large range of the initial value of the normal penalty coefficient (1 ≤ ε 0 n ≤ 10 2 ). To put in a nutshell, these academic problems show that the most frequently used methods (PM, ALM, and APM) are not adapted to solve a large variety of problems without changing the value of the penalties coefficient (the normal and the tangential one). With these classical methods, the best value of the penalty coefficient depends of the problem. Nonetheless, these academic problems are solved in very simple way with the proposed method (AALM) (the choice of the initial value of ε n is not sensitive: 1 ≤ ε 0 n ≤ 10 2 ).

Industrial application

Generally, aluminium is obtained in electrolytic cells. The reduction of the aluminium oxides is a high temperature process using electrical current. To avoid damages in the components of the cell produced by thermal shocks as well as to ensure an adequate level of containment of the cathode panel before the addition of electrolyte and molten metal, the electrolytic cell is initially preheated at a temperature around Considering these concerns, aluminium producers strive to improve their position regarding energy consumption and production costs. To do so, mathematical modelling offers a good way to study the behaviour of the electrolytic cell during its life. Considering the phenomenon occurring in the cell, the solution of this problem corresponds to a fully coupled thermo-electro-mechanical model including material non linearities and multi-physical behaviour at interfaces allowing accurate evaluation of the stress distribution in the cathode blocks and surrounding components. Therefore, the proper evaluation of mechanical contact conditions at the interface is a key issue in this process because it is directly associated to the thermo-electrical contact resistance. Several numerical models have been proposed to simulate the behaviour of the electrolytic cell (see e.g. [START_REF] Marceau | Advanced numerical simulation of the thermo-electromechanical behaviour of hall-héroult cells under electrical preheating[END_REF] or [START_REF] Dupuis | Thermal study of the coke pre-heating technique for a hall-heroult cell[END_REF]). One schematic view is presented in the Fig. 24.

In this specific example, the thermo-electro-mechanical behaviour of the cathode block is studied during a 16 h preheating phase of the cell followed by the electrolyte addition (2 h). As shown in Figs. 24 and25, the cathode block is composed of a solid carbon block sealed to a steel collector bar using cast iron. The linear increase of temperature on the top of the cathode block is such as described in Fig. 26.

Considering the symmetry of the problem, only one half of the cathode block is considered. The mesh has been taken from the work of Goulet (see [START_REF] Goulet | Modélisation du comportement thermo-électromécanique des interfaces de contact d'une cuve de Hall-Héroult[END_REF]), and is presented in Fig. 27 (The bigger length of the cathode block is about 2 m). All the material and contact properties as well as the electrical, [START_REF] Goulet | Modélisation du comportement thermo-électromécanique des interfaces de contact d'une cuve de Hall-Héroult[END_REF] thermal and mechanical boundary conditions are presented in [START_REF] Goulet | Modélisation du comportement thermo-électromécanique des interfaces de contact d'une cuve de Hall-Héroult[END_REF].

All dimensions are in metre and the penalty coefficients are in Pa/m. The steel and the cast iron have thermo-elastic behaviour. The carbon has a quasi-brittle thermo-elastoplastic behaviour, as proposed by D'Amours [START_REF] Amours | Développement de lois constitutives thermomécaniques pour les matériaux á base de carbone lors du préchauffage d'une cuve d'électrolyse[END_REF] (for more information on the values of the carbon properties see [START_REF] Mirchi | Comparative characterisation of graphitized and graphitic cathode blocks[END_REF][START_REF] Zolochevsky | Creep and sodium expansion in a semigraphitic cathode carbon[END_REF]). Regarding the contact interfaces, this problem includes six elastic foundations (multi-point constraints) that ensure the continuity of the displacements between the components under investigation and the surrounding materials as shown in Fig. 28a. The last one, located between the carbon and the cast iron (as shown in Fig. 28b), is a frictional contact interface. This is motivated by the fact that thermal expansion coefficient of these two materials are so different that sliding resulting from the thermal strains will occur. This interface is the only link between the carbon block and the collector bar ensuring the exit of the electrical current. Since the heat generation by Joule effect is related to the current intensity, it is therefore important to ensure a proper assessment of contact stresses at this interface to allow a proper computation of the thermal strains. Regarding the electrical problem, a zero voltage is applied on the top of the cathode block and a representative current intensity is applied on the exit surface of the collector bar. Convection is applied on each face of the cathode block. The preheating is initiated by a convection flux on the top of the cathode bloc using large convection coefficient. The self weight is considered as well as the weight of the anode, electrolyte and molten metal using uniform pressure on the top of the cathode block.

The solution is obtained in 18 time steps of 1 h each using the PM and the AALM for the solution of the contact problem. The comparison is performed only at the carbon/cast iron interface (Fig. 28b) considering a constant stiffness (PM) for all the other elastic foundations (Fig. 28a). For the PM, the normal and tangential penalty numbers are set to 10 10 and 10 8 respectively. For the AALM, the initial normal penalty number is set to 10 8 , the maximal value of the gap for the penalty adaptation and the tolerance on the gap are set to 10 -3 m. In addition, the tolerance on the reversible tangential displacement is equal to 10 -3 m.

As shown in Fig. 29, the AALM is faster than the classical PM for both mechanical (M) and the thermo-electrical (TE) sub-problems. In particular, the PM is unable to solve the problem for the last 2 h corresponding to the electrolytic bath addition. The AALM allows a gain of computational time (number of iterations) of about 50% for the problem. The AALM allows to control the gap and to obtain a good solution for a large value of the contact stress. Consequently, with the AALM, the solution of this mechanical problem is better. That explains the faster convergence of this kind of problem with the AALM. Even if one resolution of the mechanical sub-problem needs more iterations, the total number of iteration to solve the coupled problem is smaller. Considering the maximum gap values at the carbon/cast iron interface, Fig. 30 shows that the AALM allows very small gaps much lower than the prescribed value. The largest values appear at the end of the preheating phase, during the bath addition. At this time, the computational effort is more important due to the rapid change of the thermal strains at the top surface of the cathode block allowing important changes in the contact conditions at the interface. Also, Fig. 31 shows that the PM is unable to adequately predict the normal contact stress and more specifically approaching the bath addition where the normal contact force reaches very large value from 12 h at interface 1 (see Fig. 28). The AALM presents a smoother increase of the normal contact stress during the heat-up. This can also be observed in Fig. 32 where the normal contact stress distribution is shown at the interfaces 1 and 2. The prediction obtained with the PM (see Fig. 32a) implies a significant reduction in the section where the electric current can pass from the carbon to the collector bar. From a numerical point of view, this erroneous contact stress distribution creates an increase of the current density and therefore, a fictitious local heat-up generated by Joule's effect. This thermal shock allows significant increase of the stress in the carbon block and therefore, possible crack initiation due to excessive thermal gradient. Figure 32b shows that the AALM provides a more efficient numerical tool to predict the response of complex thermo-electro-mechanical problems.

Conclusion

This paper presents a new method for the solution of the mechanical contact problem with friction, the adapted augmented Lagrangian It is that this method has the advantages of the most frequently used methods (i.e. penalty method and augmented Lagrangian method). It allows the gap control, has the robustness of the augmented Lagrangian method and the efficiency of the penalty method. Also, the value of the convergence rate is the same that with the classical methods (PM and ALM). In addition, the user does not choose the value of any penalty coefficient. The academic examples as well as the industrial problem solved show that this method enables to get more easily a better evaluation of the contact stresses. In opposition to the more used methods (PM and ALM), the user does not modify any parameter when the geometry or the mechanicals properties of the problem change. The adapted augmented Lagrangian method has been used in the multi-physical context like industrial thermo-electro-mechanical problem [START_REF] Bussetta | Modélisation et résolution du probléme de contact mécanique et son application dans un contexte multiphysiques[END_REF]. Moreover, this method can be used to consider the deformation of the ruggedness of the contact surface.

Fig. 1 2 cFig. 2 2 cFig. 3 2 cFig. 4

 1222324 Fig. 1 Link between the contact surfaces Γ 1 c and Γ 2 c

Fig. 6

 6 Fig. 6 Value of F versus the penetration

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Example of the evolution of the gap between two iterations

Fig. 10

 10 Fig. 10 The Hertz contact problem

Fig. 11 Fig. 12

 1112 Fig. 11 Modelling of the Hertzian problem

Fig. 14 Fig. 15

 1415 Fig. 14 Contact normal stress on Hertz problem with adapted augmented Lagrangian method (ε 0 n = 2 × 10 2 )

Fig. 17 Fig. 18 2 × 10 2 )Fig. 19

 17182219 Fig. 17 Contact normal stress (mesh 4; elastoplastic problem)

Figures 21 ,

 21 22 and 23 show the distribution of the von Mises stress at the last time step with the PM in function of the values of the penalty coefficients. The von Mises stress distribution is a good indication of the value of the contact stress. The value of the von Mises stress decreases with the value of the penalty coefficients (i.e. the value of the contact stress). The solution presented on the Fig. 23 (ε n = 10 5 and ε t = 10 2 ) is considered like the reference solution with the PM.

Fig. 21 vonFig. 23 von

 2123 Fig. 21 von Mises stress distribution (MPa, PM with ε n = 10 3 and ε t = 1)

Fig. 24

 24 Fig.[START_REF] Wriggers | A formulation for frictionless contact problems using a weak form introduced by Nitsche[END_REF] Schematic view of one quarter of an electrolytic cell[START_REF] Marceau | Advanced numerical simulation of the thermo-electromechanical behaviour of hall-héroult cells under electrical preheating[END_REF]

Fig. 25 Fig. 26 Fig. 27

 252627 Fig.[START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF] Steel collector bar sealed in the bottom of the solid carbon block with cast iron (storage position)[START_REF] Goulet | Modélisation du comportement thermo-électromécanique des interfaces de contact d'une cuve de Hall-Héroult[END_REF]

Fig. 28

 28 Fig.[START_REF] Oysu | Coupled fem-bem for elastoplastic contact problems using Lagrange multipliers[END_REF] Mesh of the contact interfaces[START_REF] Goulet | Modélisation du comportement thermo-électromécanique des interfaces de contact d'une cuve de Hall-Héroult[END_REF]

Fig. 29

 29 Fig. 29 Cumulated of iterations required in function of the time step

Fig. 30 Fig. 31

 3031 Fig.[START_REF] Yamazaki | Penalty method for three-dimensional elastic contact problems by element method[END_REF] Maximum gap value at carbon/cast iron interface for each time step (adapted augmented Lagrangian method)

Fig. 32

 32 Fig. 32 Normal contact stress distribution at iron (t = 10 h)

Table 1

 1 Number of iterations to reach final stage

	ε 0 n = initial value of ε n
	"osc cnt" means oscillations preventing numerical convergence

Table 2

 2 Number of iterations to reach final stage ε 0 n = initial value of ε n "osc cnt" means oscillations preventing numerical convergence

Table 3

 3 Number of iterations to reach the final stage

	ε 0 n = initial value of ε n
	"osc cnt" means oscillations preventing numerical convergence

Acknowledgments

The authors gratefully acknowledge the financial support for this research from "Fonds québécois de la recherche sur la nature et les technologies" (FQRNT). All the computer programming has been done in the numerical tools fesh++ and contacta, respectively used for the solution of the multi-physics problems (thermo-electromechanical problems) and the contact one.


R I Leine

U Aeberhard

The Euler-Maupertuis principle of least action as variational inequality

Starting from Hamilton's Principle, the current paper discusses how we can derive the Euler-Maupertuis Principle of Least Action in the context of non-smooth dynamics. This variational principle allows us to directly obtain the space curve y(x) of a point-mass in a potential field V (x, y) without referring to the temporal dynamics. This paper generalises the Euler-Maupertuis Principle of Least Action to systems with impact by formulating the principle as a variational inequality.

Introduction

In 1744 Euler stated that the motion of a point-mass m in potential field V is described by the variational principle

√ 2m x1 x0 E -V (x, y) 1 + y (x) 2 dx = stationary, E = T + V = const.,
where x and y are the coordinates of the point mass in an inertial frame. This principle is called the (Euler-Maupertuis) Principle of Least Action and is in this form sometimes known as the Principle of Jacobi. The aim of this paper is to put the Principle of Least Action in the context of non-smooth dynamics, giving the principle a novel application field.

Principle of d'Alembert-Lagrange as variational inequality

Let us consider a point-mass m which we can push with an external force F against a frictionless wall, or even into the corner of a wall. The wall reacts with a contact force λ. Newton's second law

mr = F + λ (1)
gives the equilibrium between the inertia force, the external force F and the contact force λ. The contact force λ for a frictionless wall is an element of the normal cone,

-λ ∈ N K (r), (2)
where the normal cone N K (r) can be expressed as the subderivative ∂Ψ K (r) of the indicator function on the domain K which is enclosed by the wall. For simplicity we will assume K to be convex. If we substitute the equation of motion (1) in the set-valued force law (2), then we obtain a differential inclusion

-(mr -F ) ∈ ∂Ψ K (r) (3)
for the non-impulsive dynamics of the point-mass with unilateral constraint. Using some concepts of convex analysis we obtain the following result:

Ψ K (r * ) ≥ Ψ K (r) -λ T (r * -r) ∀ r * 0 ≥ -λ T (r * -r) ∀ r * ∈ K, δr = r * -r λ T δr = (mr -F ) T δr ≥ 0 ∀ δr ∈ T K (r) (4)
The virtual work λ T δr of the contact force λ is nonnegative for all virtual displacements which are taken from the tangent cone T K (r), which is the cone of all kinematically admissible virtual displacements. The expression (mr -F ) T δr ≥ 0, ∀ δr ∈ T K (r) is called a variational inequality. Consider the principle of d'Alembert-Lagrange for a perfect bilateral constraint: the virtual work of a perfect bilateral constraint force vanishes for all kinematically admissible virtual displacements. The variational inequality which we just obtained is nothing else than the principle of d'Alembert-Lagrange for a perfect unilateral constraint.

We now integrate the principle of d'Alembert-Lagrange as variational inequality over a compact time-interval I = [t 0 , t 1 ]. This time-interval can contain collision times which we gather in the set D = {t i }. For the moment, we leave the impact times aside and we integrate over the interval I\D:

I\D (mr -F ) T δr dt ≥ 0 ∀ δr ∈ T K (r). (5)
As is customary in this kind of analysis we perform a partial integration.

m ṙδr t1 t0 + i m ṙδr t + i t - i - I\D m ṙT δ ṙ dt - I\D F T δr dt ≥ 0
The boundary terms for t 0 and t 1 vanish for fixed boundary conditions r(t 0 ) = r 0 , r(t 1 ) = r 1 . The boundary terms at the collision times D vanish if the so-called Weierstrass-Erdmann conditions are fulfilled. The first Weierstrass-Erdmann condition requires that the impulsive forces lie within the normal cone, -Λ ∈ N K (r). The second Weierstrass-Erdmann condition demands that the collisions are completely elastic such that there is no energy loss during the impact. The second condition is superfluous if we consider a special class of test functions which lead to Gâteaux derivatives (see [START_REF] Leine | Hamilton's principle as variational inequality for mechanical systems with impact[END_REF]). The remaining terms can be recognised as the variation of the kinetic energy δT = m ṙT δ ṙ and the potential energy δV = F T δr, i.e.

-

I\D δT dt + I\D δV dt ≥ 0. (6)
We are still integrating over the time-interval without collision times. However, the missing integration parts are Lebesgue integrals and the collision times are Lebesgue negligible:

D δT dt = 0 D δV dt = 0. (7)
We obtain the principle of Hamilton as variational inequality,

-δ I L dt ≥ 0 L := T -V ∀ δr ∈ T K (r), (8)
which also accounts for perfect unilateral constraints (see [START_REF] Leine | Hamilton's principle as variational inequality for mechanical systems with impact[END_REF]).

Principle of Least Action in inequality form

In the same way we are also able to generalise the principle of least action to unilateral constraints. We can do so by starting from the principle of Hamilton (8) in inequality form and by taking the conservation of energy E = T + V as auxiliary condition. This relates the variation of the potential energy δV = -δT to the variation of the kinetic energy. Substitution of the latter in the principle of Hamilton gives a Lagrange function which is 2 times the kinetic energy -2 δT dt ≥ 0 ∀ δq ∈ T K (q), in which we used generalised coordinates q. A reformulation of the kinetic energy T = 1 2 qT M (q) q = 1 2 q T M (q)q dθ dt 2 gives the principle of least action in inequality form:

-√ 2 δ E 0 -V (q) q T M (q)q dθ ≥ 0 ∀ δq ∈ T K (q).

Once again, we obtain a variational inequality.