# The inequality level-set approach to handle contact: membrane case

Matthieu Graveleau, Nicolas Chevaugeon, Nicolas Moes

► To cite this version:

Matthieu Graveleau, Nicolas Chevaugeon, Nicolas Moes. The inequality level-set approach to handle contact: membrane case. Advanced Modeling and Simulation in Engineering Sciences, 2015, 2 (1), 10.1186/s40323-015-0034-8. hal-01398280

HAL Id: hal-01398280

https://hal.science/hal-01398280

Submitted on 17 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

## RESEARCH ARTICLE

## Open Access

<figure><img src="image_3.png" /></figure>

CrossMark

# The inequality level-set approach to handle contact: membrane case

Matthieu Graveleau, Nicolas Chevaugeon* and Nicolas Moës

### Abstract

**Background:** Contact mechanics involves models governed by inequality constraints. Even for the simplest contact problem, inequalities arise from the lack of information on the contact zone position. In addition to increasing the difficulty to solve such problems, an unknown contact zone makes it difficult to use an appropriate mesh and to represent efficiently phenomena on the contact zone boundary. Nevertheless these phenomena are often crucial to have an accurate representation of the problem such as weak discontinuity of the displacement.

**Methods:** In this paper, we propose a method specifically designed to solve inequality constraint problems linked to an unknown domain without remeshing. In order to do so, level sets coupled with X-FEM is used to define the unknown domain and take into account the specific behavior at the contact zone boundary. The key idea of the method is to split the problem involving inequality constraints into two problems. In the first problem, the unknown domain is set and therefore it only involves equalities. Nevertheless, the constraints might be violated, meaning the set domain has to be changed. Then, the other problem is a shape optimization of this domain and leads to an updated set domain. These two problems are iterated up to convergence of the algorithm. Moreover, the addition of adhesion to the problem will be considered.

**Results:** The studied case in this paper is a membrane in the context of small deformations. First, a 1D example will be given to illustrate the method with and without adhesion. Then 2D cases will be studied. Finally an example with an evolving load will be given. Comparison will be made with a classical active-set method.

**Conclusions:** The ILS is proved to be an efficient method giving a convincing accuracy for the contact boundary without need of re-meshing. It is also able to naturally handle adhesion.

**Keywords:** X-FEM, Level-set, Contact, ILS, Adhesion, Membrane

### Background

Being able to predict how bodies in contact are going to behave is important. Indeed, such situations are omnipresent in daily life. Thus, contact mechanics has been intensively studied. Nowadays several methods are available to tackle contact problems. In this paper we will focus on the methods based on the finite element method [1] even if other options are available such as iso-geometric methods, for which a review can be found in [2]. Furthermore, only contact with a rigid body will be studied, therefore avoiding, for now, the problem of non-matching meshes. Methods have been designed to overcome this challenge.

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence: nicolas.chevaugeon@ec-nantes.fr
GeM, UMR CNRS 6183, Ecole Centrale de Nantes, 1 rue de la Noë, 44321 Nantes, France</footnote>
Among them, the mortar methods [3, 4] have shown great effectiveness. For our focus of study, the most common methods are the Lagrange multipliers method [5, 6], the penalty method [5, 6] and the augmented Lagrangian method [6, 7]. One of the main difficulty of contact is that the contact area is a *priori* unknown. This leads to the introduction of an inequality and consequently, computational challenges even in the simplest contact problem [5]. They usually require iterative algorithms and can be conditionally convergent. These difficulties are common to any inequality problems. Another issue to be addressed is the representation of the contact boundary area and the phenomena on it. In fact, the finite element method is strongly dependent on a predefined mesh. However, as the contact area is unknown, it is impossible to design a suitable mesh at first. In addition, low regularity of the solution are often present at the boundary of the contact area. The X-FEM [8] is a definite asset in representing these discontinuities, leading to a higher order convergence rate with respect to element size. We propose in this paper a new approach for contact, the ILS, which is using X-FEM and is designed to handle these problems.

In the first section of this paper we will state a general framework for the class of problem we are going to study. Then we will give a brief review of the most used methods. In the next section, we will set up the ILS framework applied to contact. We will explain how the variational inequality problem is transformed into an equality problem coupled with a shape optimization. We will also emphasize the particular tools required such as the notion of level-set and the X-FEM. The addition of adhesion behavior to the method will be discussed. Finally we will give examples of our method applied to membrane problems. First, a 1D case will be presented to allow an easier understanding of the method. Then the 2D case will be studied in order to use the full extent of the method.

## Description of the problem and classical methods

### The problem

Overall, the goal of contact mechanics is to be able to handle the physics governing multiple bodies coming into contact with each other. In this paper we want to focus on the understanding of the ILS method. Therefore, we will base our explanations on a simple problem even if the ILS could be generalized. The studied problem is a deformable membrane in contact with a rigid surface considered at least $C^1$. We will assume the contact to be friction-less and we will focus only on the static problem with small deformation hypothesis. The deformable body occupies the volume $\Omega$, its frontier is $\partial\Omega$ and $\boldsymbol{x} = (x, y, z)$ is the coordinate of a point. We will denote $u(\boldsymbol{x})$ the displacement of the point $\boldsymbol{x}$ and $T$ the tension of the membrane. The distance between the solid and deformable body without any load applied is represented by a function $d(\boldsymbol{x}) > 0$. The non penetration of the two bodies can be written $u + d \ge 0$. Without loss of generality zero displacements are imposed on $\partial\Omega$ and only continuous body forces $f_d$ will be applied. Finally, the part of the boundary in contact with the rigid surface will be called $\Omega^-$ and its boundary $\Gamma$. We will denote everything defined only on the contact side of $\Gamma$ by a “-” and on the other side a “+”. If we denote the closure of a domain $\Omega$ by $\bar{\Omega}$ then the domains must satisfy $\bar{\Omega}^+ \cup \bar{\Omega}^- = \bar{\Omega}$ and $\Omega^+ \cap \Omega^- = 0$. We will define on every point of the boundary $\Gamma$ the normal $n$ to it contained in the membrane surface. It is chosen to be directed towards $\Omega^+$. Figure 1 illustrates this description.
<figure><img src="image_5.png" /><figcaption>**Figure 1** Membrane into contact with a rigid surface.</figcaption></figure>

Our problem being parametrized, we can build the mathematical framework needed to find the equilibrium solution of this problem. Of course, this framework has already been settled and the reader can refer to [5] for more details. Indeed, we will only outline the main results needed here.

The governing equations of the equilibrium of a membrane are well known. Without contact, every point of the membrane has to fulfill the strong equations:

$$ \begin{cases} T \Delta u + f_d = 0 & \text{on } \Omega, \\ u = 0 & \text{on } \partial \Omega. \end{cases} \qquad (1) $$

With contact, the set of equations is more complex and involves inequalities. If there is contact the classical equivalent of (1) can be written with an addition of the Signorini equations (3), also called Karush Kuhn Tucker (KKT) conditions:

**Problem 1** Find $u$ and $\tilde{p}$ such as:

$$ \begin{cases} T \Delta u + f_d + \tilde{p} = 0 & \text{on } \Omega, \\ u = 0 & \text{on } \partial \Omega. \end{cases} \qquad (2) $$

$$ \left\{ \begin{array}{l} u + d \geq 0 \text{ on } \Omega, \\ \tilde{p} \geq 0 \text{ on } \Omega, \\ \tilde{p}(u+d) = 0 \text{ on } \Omega. \end{array} \right. \qquad (3) $$

where $\tilde{p}$ is the normal reaction exerted by the rigid surface to the membrane.

In Problem 1, all the equations have been generalized to the whole domain. However, if we split the domain into $\Omega^+$ and $\Omega^-$ (being the actual contact zone), as described previously, Problem 1 is equivalent to:
**Problem 2** Find $u, \tilde{p}$ and $\Omega^-$ such that:

<table><thead><tr><th></th><th>ln Ω<sup>+</sup></th><th>ln Ω<sup>-</sup></th><th>On Γ</th></tr></thead><tbody><tr><td>Equality</td><td>T Δu + f<sub>d</sub> = 0<br>u = 0 on ∂Ω</td><td>T Δu + f<sub>d</sub> + p̃ = 0<br>- u = d</td><td>u<sup>+</sup> = u<sup>-</sup></td></tr><tr><td>Inequality</td><td>-u ≤ d</td><td>p̃ ≥ 0</td><td></td></tr></tbody></table>

The separation of the two domains will be a key point for the method exposed in "Method: the ILS applied to contact problems".

Another problem, which is closer to what is solved in practice, is to find the equilibrium solution assuming a contact zone.

**Problem 3** Find $u, \tilde{p}$ and $\tilde{\lambda}$ for a given $\Omega^-$ such that:

<table><thead><tr><th></th><th>ln Ω<sup>+</sup></th><th>ln Ω<sup>-</sup></th><th>On Γ</th></tr></thead><tbody><tr><td>Equality</td><td>T Δu + f<sub>d</sub> = 0<br>u = 0 on ∂Ω</td><td>T Δu + f<sub>d</sub> + p̃ = 0<br>- u = d</td><td>T[∇u] · n̄ = p̃<br>u<sup>+</sup> = u<sup>-</sup></td></tr></tbody></table>

with $[\nabla u(s)] = \nabla u(s)^+ - \nabla u(s)^-$ the jump of $\nabla u$ at a point of $\Gamma$ located by its curvilinear abscissa s and $\tilde{\lambda}$ the reaction on the boundary of the rigid solid (= 0 if $\Omega^-$ is exact). To be equivalent to Problems 2, 3 needs to be solved for the right $\Omega^-$.

Let us show that if the contact conditions are fulfilled, i.e. the contact zone is the right one, we must satisfy $[\nabla u] \cdot n = 0$. First, let us look at the case $[\nabla u] \cdot n > 0$, thus $\tilde{\lambda} > 0$. This means that the membrane is attracted by the rigid body. Without adhesion this is not a physical solution. Now, if $[\nabla u] \cdot n < 0$ and therefore $\nabla u^+(s) \cdot n < \nabla u^-(s) \cdot n$. As by definition $\nabla d(s)$ is continuous and $\nabla u^-(s) = \nabla d(s)$ we must have $\nabla u^+(s) \cdot n < \nabla d(s) \cdot n$ and therefore there is penetration of the rigid surface.

Hence,

$$[\nabla u] \cdot n = 0 \quad \forall s \in \Gamma \qquad (4)$$

is a necessary condition to have the right contact zone. This property will reveal itself quite useful in the next section.

In any case, a strong formulation is often not suited for computational purpose. The classical approach is to derive the weak form of this formulation. Let us first define the admissible displacements space that will be used in this article:

$$V_d = \{u \in H^1(\Omega^+ \cup \Omega^-) | u = 0 \text{ on } \partial\Omega\} \qquad (5)$$

Therefore if $u \in V_d$ it can admit a strong discontinuity on $\Gamma$. Problem 3 can be written as the weak formulation:

**Problem 4** Find $u \in V_d, p \in P$ and $\lambda \in H^{-1/2}(\Gamma)$ such that:

$$\left\{
\begin{array}{ll}
a(u, u^*) + \langle p, u^* \rangle_P + \langle \lambda, [u^*] \rangle_{H^{-1/2}} = f(u^*) & \forall u^* \in V_d, \\
\langle p^*, u + d \rangle_P = 0 & \forall p^* \in P, \\
\langle \lambda^*, [u] \rangle_{H^{-1/2}} = 0 & \forall \lambda^* \in H^{-1/2}(\Gamma).
\end{array}
\right.
\qquad (6)$$
with:

$$ 
\left| 
\begin{array}{l}
a(u, u^*) = \int_{\Omega} T \nabla u \cdot \nabla u^* \, d\Omega \qquad f(u^*) = \int_{\Omega} f_d u^* \, d\Omega \\
\langle \lambda^*, \llbracket u \rrbracket \rangle_{H^{-1/2}} = c(\lambda^*, u) = \int_{\Gamma} \lambda^* \llbracket u \rrbracket \, d\Gamma
\end{array}
\right. 
\qquad (7) 
$$

Two options for P are available. The most obvious one is $L^2(\Omega^-)$ meaning that the duality pairing reads therefore $\langle p, u \rangle_{L^2} = \int_{\Omega^-} pu \, d\Omega$. Even if this formulation appears as the straightforward way to impose a Dirichlet condition it is not well suited to satisfy the Babuška–Brezzi (BB) condition. The other choice is $H^1(\Omega^-)$ and the duality pairing becomes $\langle p, u \rangle_{H^1} = \int_{\Omega^-} l_m \nabla p \cdot \nabla u \, d\Omega + \int_{\Omega^-} pu \, d\Omega$. The addition of $l_m$ is only there for dimensional homogenization and will be set to $l_m = 1m^2$. This last formulation fulfills naturally the Babuška–Brezzi (BB) condition as the operator is the natural scalar product associated to the space where both u and p live. Consequently, it is chosen for this paper and Problem 4 writes:

**Problem 5** Find $u \in V_d, p \in H^1(\Omega^-)$ and $\lambda \in H^{-1/2}(\Gamma)$ such that:

$$ 
\left\{ 
\begin{array}{ll}
a(u, u^*) + b(p, u^*) + c(\lambda, u^*) = f(u^*) & \forall u^* \in V_d, \\
b(p^*, u) = l_p(p^*) & \forall p^* \in H^1(\Omega^-), \\
c(\lambda^*, u) = 0 & \forall \lambda^* \in H^{-1/2}(\Gamma).
\end{array}
\right. 
\qquad (8) 
$$

with:

$$ 
\left| 
\begin{array}{l}
b(p, u^*) = \int_{\Omega^-} l_m \nabla p \cdot \nabla u^* + pu^* \, d\Omega \qquad l_p(p^*) = - \int_{\Omega^-} l_m \nabla p \cdot \nabla d + pd \, d\Omega
\end{array}
\right. 
\qquad (9) 
$$

A drawback of this formulation is the loss of equivalence between both $\tilde{p}$ and $p$ just as $\tilde{\lambda}$ and $\lambda$. Therefore, post-processing will be needed to compute the reaction of the ground. In this section, several ways of expressing a contact problem have been shown. To summarize, it is convenient to break down the full problem given by Problem 1 into two problems: solve Problem 5 for a given $\Omega^-$ on one hand and find an admissible $\Omega^-$ on the other hand. We will now discuss a brief state of the art of the methods used to solve these problems.

### Some existing approaches

As we said, this problem is numerically difficult due to the inequality constraint and iterative algorithms are needed. The most common approach is to iterate on the set of degrees of freedom for which the constraint is active. Such methods are called active-set methods. Nonetheless, the contribution of the contact constraint have to be taken into account. The most spread methods in the industry are the penalty method and the Lagrange multipliers methods. An overview of these methods can be found in [5, 6]. Some methods combine both of them. This is the case of the perturbed Lagrange formulation [9] or the augmented Lagrange formulation [7]. Another approach is to use tools from mathematical programming like in [10].

The penalty method being the easier to implement is widely used. It also does not need to compute extra degrees of freedom as in the Lagrange multipliers method.
Nevertheless, it relies on the choice of a penalty parameter. If the latter is too small, penetration is allowed whereas if it is too big, the problem might become ill-conditioned. Another spread method is the Lagrange multipliers method which leads to an exact fulfillment of the contact condition. However it does require additional degrees of freedom and a special care is to be taken for rigid body motion between two contacting bodies. To overcome these difficulties, an augmented Lagrange formulation is often coupled to an Uzawa algorithm. Combining both advantages of the penalty and Lagrange multiplier methods is more reliable but also more difficult to implement. Again, a comprehensive survey of these methods can be found in [6] or [11] for instance.

Even with the latter one, other difficulties arise. One of the most obvious one is that the representation of the contact zone is strongly linked with the discretization. As the contact zone is a *priori* unknown it is impossible to use an appropriate mesh at first. A solution to this problem is to move the nodes of the mesh to the approximated boundary of the contact zone [12]. A comparison of this method against more classical ones can be found in [13]. Obviously this method needs to update the mesh at each iteration. In the next section, we propose a new method to tackle contact especially designed to solve inequality constrained problems on a fixed mesh.

## Method: the ILS applied to contact problems

The ILS method has been designed to be able to treat problems involving inequality constraints arising from an unknown constrained region. This means that if we can predict the constrained region there is no inequality anymore. Several problems fall into this category. Between them, a differentiation can be made. In the case of 3D solids the constrained domain can be a volume or a boundary of the domain (thus a surface). An example of the first case was treated in [14] where the ILS was first applied. It dealt with volumetric kinematic constraints inside a 3D domain. Contact problems are clearly part of the other category. Indeed if we know the contact area on the boundary we can prescribe the distance of the body to the rigid surface to be strictly zero and no more greater than zero. Nevertheless, even if the problem we are dealing with is also a contact problem we neglected the thickness of the membrane. Therefore the dimension of the contact zone is the same as the dimension of the membrane. This does not change the main idea of the ILS method but only local details. The Dirichlet boundary condition on $Ω^−$ can be enforced using Lagrange multipliers for instance. Therefore the challenge is to find this contact zone. Starting from an initial contact zone, a shape optimization is set to make the contact zone evolve toward one which allows all the contact conditions to be fulfilled. In short, the main idea of the method is to transform the inequality problem to a sequence of equality problems involving a shape optimization at each step.

## The equilibrium problem

To be able to solve the equilibrium problem we need to differentiate the domains and thus to define the boundary of the contact zone $Γ$. In order to do so, we are going to use a level-set [15, 16] and particularly its iso-zero set to fit $Γ$ (cf. Figure 2b). A particular case of level-set is a signed distance function $\phi(x)$ (cf. Figure 2c). There are several advantages to using a level-set. The first one is that it is easy to define if a point is inside or outside the contact zone by looking at the sign of the level-set. By convention we choose $\phi$ to be negative inside $Ω^−$. Furthermore, it is easy to compute set theory operations such as the union
<figure><img src="image_6.png" /><figcaption>**Figure 2** The level-set has to fulfill the properties illustrated in (b). a Shows one of the possible choice and (c) the particular case of a signed distance function (dashed line represents the iso-zero).</figcaption></figure>

or the difference of sets represented by a level-set. Usually it is only a question of finding the maximum or the minimum of the related level-sets. Nevertheless, in order for $\phi$ to stay a signed distance function it needs to be reset. To reset $\phi$ the iso-zero is kept and the distance to it is re-computed. Last but not least, a level-set can evolve on a fixed mesh so we do not have to re-mesh.

To solve the problem with a given contact zone, a numerical strategy needs to be set. A first step is to find a good discretization of the fields. We will start with classical shape functions called $N_i$. For instance it could be the well known hat functions. With such functions we can approximate the value of the level-set at each point of the domain by:

$$ \phi(x) = \sum_{i \in I} \phi_i N_i(x) \qquad (10) $$

where $I$ represents the number of shape functions. However in order to represent the displacement field, this classical approximation is not enough. Indeed, the solution might present weak discontinuities on the boundary. Hence, we are using the X-FEM method [8] to enrich the displacement. If we choose to use an enrichment function $F$ suited to our problem an approximation of the field $u$ is:

$$ u_h(x) = \sum_{i \in I} u_i N_i(x) + \sum_{j \in J} u_j^e N_j(x) F(x) \qquad (11) $$

with $J$ the subset of $I$ which is enriched.

A straightforward choice for $F$ would be a ridge function allowing a weak discontinuity of the displacement field. Then Lagrange multipliers have to be defined on $\Gamma$ in order to represent the linear force needed to physically allow a slope discontinuity in the membrane. Nonetheless, an equivalent strategy is to set $F$ to be an Heaviside function and to use the above mentioned Lagrange multipliers to cancel the jump of the displacement $[u]$. A better numerical behavior of the Heaviside function, especially for high order finite element, motivates this choice.

Lagrange multipliers also need to be discretized. In order to have good convergence properties we need to use the right approximation for $p$ and $\lambda$ such as the BB condition is satisfied. This problem was handled in [17] for Lagrange multiplier on a boundary (i.e. for $\lambda$) coupled with X-FEM. Rather than degree of freedoms linked directly to the front, a combination of degree of freedoms attached to the elements cut by the level-set is used. Hence this discretization is used here. Furthermore, the choice of the formulation guarantees a good behavior of $p$. We can then approximate $p$ and $\lambda$:
$$p_h(\underline{x}) = \sum_{i \in I^p} p_i N_i^p(\underline{x}) \qquad (12)$$

$$\lambda_h(\underline{x}) = \sum_{i \in I^{\lambda}} \lambda_i N_i^{\lambda}(\underline{x}) \qquad (13)$$

with $I^p$ and $I^{\lambda}$ the set of degrees of freedom for each field. Finally, with such a discretization we can apply a Galerkin method to approximate Problem 5 and write it in a matrix form:

**Problem 6** Find $\mathbf{u}$, $\mathbf{u}^e$, $\mathbf{p}$ and $\lambda$ such that:

$$[\mathbf{K}] \begin{Bmatrix} \mathbf{u} \\ \mathbf{u}^e \\ \mathbf{p} \\ \lambda \end{Bmatrix} = \begin{Bmatrix} \mathbf{f} \\ \mathbf{f}^e \\ l_p \\ 0 \end{Bmatrix} \quad \text{where } [\mathbf{K}] = \begin{bmatrix} \mathbf{A}^{cc} & \mathbf{A}^{ce} & \mathbf{B} & \mathbf{C} \\ \mathbf{A}^{ceT} & \mathbf{A}^{ee} & \mathbf{B}^e & \mathbf{C}^e \\ \mathbf{B}^T & \mathbf{B}^{eT} & 0 & 0 \\ \mathbf{C}^T & \mathbf{C}^{eT} & 0 & 0 \end{bmatrix} \qquad (14)$$

with:

$$\left| \begin{array}{lll} \mathbf{u} = \{u_i\} & A_{ij}^{cc} = \int_{\Omega} \nabla N_i T \nabla N_j d\Omega & \mathbf{f} = \left\{ \int_{\Omega} N_i f d\Omega \right\} \\ \mathbf{u}^e = \{u_i^e\} & A_{ij}^{ee} = \int_{\Omega} \nabla N_i^e T \nabla N_j^e d\Omega & \mathbf{f}^e = \left\{ \int_{\Omega} N_i^e f d\Omega \right\} \\ \mathbf{p} = \{p_i\} & B_{ij} = \int_{\Omega^{-}} N_i N_j^p d\Omega & l_p = \left\{ \int_{\Omega^{-}} N_i^p d\Omega \right. \\ & \left. + \int_{\Omega} l_m \nabla N_i \nabla N_j^p d\Omega \right. & \left. + \int_{\Omega} l_m \nabla N_i^e \nabla N_j^p d\Omega \right. \\ & C_{ij} = \int_{\Gamma} N_i N_j^{\lambda} d\Gamma & C_{ij}^e = \int_{\Gamma} N_i^e N_j^{\lambda} d\Gamma \end{array} \right. \qquad (15)$$

The size of the matrix $\mathbf{A}^e$ is usually very small compare to the size of $\mathbf{A}^c$. This last matrix is not changing through the iteration, which is numerically convenient. Furthermore, from a broader point of view, when it will come to deal with plain solid both the size of $\mathbf{B}$ and $\mathbf{C}$ will be modest compare to $\mathbf{A}^c$. Therefore a great deal of computational effort can be spared by working on $\mathbf{A}^c$ only once.

## Shape optimization

At this point, we can solve Problem 2 for a given contact zone. However the solution might not fulfill the contact conditions (i.e. the contact zone is not the right one). We are going to search for the exact contact zone with an iterative process. As with an active-set method, we are using the fields computed at the current iteration to determine the evolution of the contact zone. The main difference is that we are not using the same criterion as in active-set method to make the contact zone evolve. An active-set method uses directly the value of the fields. For instance a criterion would be if the Lagrange multipliers are positive then remove the constraint at the node. The ILS is mainly using another criterion $\rho$ which depend on the problem it is dealing with. This criterion is chosen strictly related to the boundary of the contact zone rather than whole of the solid boundary. Let us denote $\Gamma^k$
the boundary of the contact zone at the kth iteration. To lighten the notations we will usu-
ally only write $\Gamma$. The abovementioned criterion must satisfy:

$$
\Gamma = \Gamma^{\text{ex}} \Rightarrow \varrho = 0 \text{ on } \Gamma \qquad (16)
$$

where $\Gamma^{\text{ex}}$ is the exact boundary of the contact zone. Thus if $\varrho \neq 0$, the position of $\Gamma$ is not the right one. This leads to try to find at each iteration the evolution of the contact zone which cancel $\varrho$. Nevertheless, having $\varrho = 0$ is not always sufficient to insure that $\Gamma = \Gamma^{\text{ex}}$. Indeed the KKT conditions might not be fulfilled everywhere. Therefore, time to time, these conditions are checked on the full domain in order to take into account contact zone nucleation and disappearance.

Here, the obvious choice for $\rho$ is $[\nabla u] \cdot \mathbf{n}$, see Eq. (4).

An efficient way to find the zero of a functional is to use a Newton method. Therefore, if we can compute the derivative of $\rho$ with respect to a variation of the contact zone shape we can apply a Newton algorithm to update $\Gamma$. The derivative which is going to be used is the directional derivative $D[\cdot]$. Physically, this derivative is the variation of the quantity considered, following the moving interface. Without loss of generality, as we are only interested in the change of shape of the boundary, we can only consider movement of the interface in the normal direction. This special case is called normal directional derivative. To define this derivative, let us introduce a normal displacement of the interface $w_T = w_T - \mathbf{n}$. Then we denote $\Gamma(0)$ the original position of the boundary and $\Gamma(\zeta w_T)$ the new position after a displacement $w_T$ (Figure 3). The normal directional derivative of $\rho$ is then defined as:

$$
D\rho[\underline{w}_T] = \lim_{\zeta \to 0} \frac{\rho(\Gamma(\zeta \underline{w}_T)) - \rho(\Gamma(0))}{\zeta} \quad (17)
$$

Then the zero of $\rho$ can be found by solving:

**Problem 7** Find $\underline{w}_T \in L^2(\Gamma)$ such that:

$$
\rho + D\rho[\underline{w}_T] = 0 \text{ on } \Gamma \tag{18}
$$

$$
\Leftrightarrow \int_{\Gamma} \rho \underline{w}_{T}^{*} \, d\Gamma + \int_{\Gamma} D\rho[\underline{w}_{T}] \underline{w}_{T}^{*} \, d\Gamma = 0 \quad \forall \underline{w}_{T}^{*} \in L^{2}(\Gamma) \quad (19)
$$

<figure><img src="image_7.png" /><figcaption>**Figure 3** Evolution of the boundary following a path $\underline{w}_T$.</figcaption></figure>
We will denote the normal directional derivative of a quantity a: $D(a)[w_T] = \dot{a}$. In this particular case, $\dot{\rho} = \overline{\overline{\nabla u}} \cdot \underline{n} = \overline{\overline{\nabla \dot{u}}} \cdot \underline{n} - \overline{\overline{\nabla u}} \cdot \nabla w_T \cdot \underline{n} + \overline{\overline{\nabla u}} \cdot \nabla \underline{n} \cdot w_T$. To compute it, we are using the same trick as in [18] or [19]. We assimilate the change of shape to a body motion which is a classical problem with plenty of tools at our disposal.

With these new tools we are going to take the directional derivative of Problem 5 to have a system with $(\dot{u}, \dot{p}, \dot{\lambda})$ as unknowns. In order to compute $\dot{u}, \dot{p}$ and $\dot{\lambda}$ we need to extend the definition of $D[\cdot]$ to the whole domain. The boundary displacement $w_T$ is then extended onto this domain. This extended displacement is noted $w$. But, as we are only, in practice, interested by the evolution of $\rho$ on the contact boundary, we are going to introduce a function $r(\phi)$. This means that this function depends on the distance from $\Gamma$ in the normal direction. This function allows us to set the displacement value to zero after a certain distance $\delta$ and thus to reduce the computational cost (cf. Figure 4). Finally, in order to keep the extended displacement normal to the boundary of the contact area we are going to use the gradient of the level-set $\nabla\phi(\underline{x})$ as its norm is one and it is normal to the iso-level-set.

$$ w(\underline{x}) = w_{\Gamma}(s_c) \nabla\phi(\underline{x}) r(\zeta) \qquad (20) $$

where $s_c$ represent the curvilinear abscissa of the closest point of $\underline{x}$ on $\Gamma$.

Now with this extended displacement, the normal directional derivative can be defined for any quantity $f$ in the domain in the same fashion as before.

$$ Df[w] = \lim_{\zeta \to 0} \frac{f(\Gamma(\zeta w)) - f(\Gamma(0))}{\zeta} \qquad (21) $$

It is shown in Appendix A that the derivative of Problem 5 is:

**Problem 8** Find $\dot{u} \in V_d, \dot{p} \in H^1(\Omega^-)$ and $\dot{\lambda} \in H^{-1/2}(\Gamma)$ such that:

$$ \left\{ \begin{array}{ll} a(\dot{u}, u^*) + b(\dot{p}, u^*) + c(\dot{\lambda}, u^*) = f^s(u^*, w) & \forall u^* \in V_d, \\ b(\dot{p}, \dot{u}) = l_p^s(\dot{p}^*, w) & \forall \dot{p}^* \in L^2(\Omega^-), \\ c(\dot{\lambda}^*, \dot{u}) = l_\lambda^s(\dot{\lambda}^*, w) & \forall \dot{\lambda}^* \in H^{-1/2}(\Gamma). \end{array} \right. \qquad (22) $$

<figure><img src="image_8.png" /><figcaption>**Figure 4** Function $r$ in the case of a circular level-set.</figcaption></figure>
where:

$$
\begin{align}
f^s(u^*, w) &= - \int_{\Omega} T (\nabla u^T \nabla u^* \operatorname{div} \underline{w} - (\nabla u^T (\nabla \underline{w}^T + \nabla \underline{w}) \nabla u^*)) \, d\Omega \nonumber \\
&\quad - \int_{\Omega^-} l_m (\nabla p^T \nabla u^* \operatorname{div} \underline{w} - (\nabla p^T (\nabla \underline{w}^T + \nabla \underline{w}) \nabla u^*)) \, d\Omega \nonumber \\
&\quad - \int_{\Omega^-} p u^* \operatorname{div} \underline{w} \, d\Omega - \int_{\Gamma} \lambda \llbracket u^* \rrbracket (\operatorname{div} \underline{w} - \underline{n} \cdot \nabla \underline{w} \cdot \underline{n}) \, d\Gamma + \int_{\Omega} f_d u^* \operatorname{div} \underline{w} \, d\Omega \nonumber \\
l_p^s(p^*, w) &= - \int_{\Omega^-} p^* \, d\Omega - \int_{\Omega^-} l_m \nabla p^* \nabla \, d\Omega \nonumber \\
l_\lambda^s(\lambda^*, w) &= - \int_{\Gamma} \lambda^* \llbracket u \rrbracket (\operatorname{div} \underline{w} - \underline{n} \cdot \nabla \underline{w} \cdot \underline{n}) \, d\Gamma \tag{23}
\end{align}
$$

It is important to note that Problem 8 has the same left-hand side as Problem 5. Using
the discretization given in (11, 12, 13) Problem 8 allows the following discrete form:

**Problem 9** Find $\hat{\mathbf{u}}$, $\hat{\mathbf{p}}$ and $\hat{\mathbf{\lambda}}$ such that:

$$
[K] \begin{Bmatrix} \hat{\mathbf{u}} \\ \hat{\mathbf{p}} \\ \hat{\mathbf{\lambda}} \end{Bmatrix} = \begin{Bmatrix} \mathbf{f}^s \\ \mathbf{l}_p^s \\ \mathbf{l}_\lambda^s \end{Bmatrix} \qquad (24)
$$

and thus the left-side matrice is the same as the one of Problem 8. Thereby we can use
the matrix decomposition already computed to speed up the process.

However, as $w_T$ is unknown, we cannot directly compute Problem 9. Therefore, $w_T$ is
going to be approximated by $M$ given modes $w_m$. For each mode, a $\tau_m$ is associated and
represent the $m$th mode contribution to $w_T$.

$$
w_T = \sum_m \tau_m w_m. \tag{25}
$$

There is plenty of possibilities to construct the modes. In this paper we start for instance
with the Fourier modes $w_m = \mathcal{F}_m(\alpha)$ ($\alpha$ being a parametrization angle for the contact
area boundary as in Figure 4). To conclude one step of the algorithm we have to solve
$M$-times Problem 9 in order to get all the $\hat{\varrho}_m$ and then solve Problem 10.

Problem 7 can thereby be approximated using a Galerkin approximation:

**Problem 10** Find $\tau = \{\tau_m\} \in \mathbb{R}^M$ such that:

$$
\sum_n \int_{\Gamma} \rho w_n \, d\Gamma + \sum_m \sum_n \int_{\Gamma} \tau_m D\rho [w_m \underline{n}] w_n \, d\Gamma = 0 \quad (26)
$$

$$
\Leftrightarrow \sum_n \int_{\Gamma} \rho w_n \, d\Gamma + \sum_m \sum_n \int_{\Gamma} \tau_m \partial_m w_n \, d\Gamma = 0 \quad (27)
$$

Again, Problem 10 can be written as:

**Problem 11** Find $\tau \in \mathbb{R}^m$ such that:

$$
[K_{\Gamma}] \tau = \rho^{\pi} \tag{28}
$$
<figure><img src="image_9.png" /><figcaption>Figure 5 Iterative algorithm to find Γ for a given load and a fixed number of mode.</figcaption></figure>

where:

$$K_{nm} = \sum_{i} \int_{\Gamma} \rho_{i} N_{i}^{j} [w_{m}] w_{n} d\Gamma \qquad (29)$$

$$\rho_{n}^{\pi} = \sum_{i} \int_{\Gamma} \rho_{i} N_{i}^{j} w_{n} d\Gamma \qquad (30)$$

These steps (Problem 6, M-times Problem 9 and Problem 11), illustrated in Figure 5, are iterated up to convergence. To be able to compute the new Γ we need to update the level-set using the displacement w. In order to do so, an Hamilton–Jacobi equation $(\partial\phi/\partial t = w \cdot \nabla\phi)$ is solved [15]. In the case where a fixed number of modes is used, the algorithm is stopped when the norm of projection of ϱ noted π(ϱm) onto the mode space is small enough (i.e. inferior to a scalar πmin), which means that we could do no better with the modes used. It is also interesting to have an increasing number of modes through the iterations. Indeed, at the beginning, having only a few modes allow a fast computation to get closer to the solution. Then a greater number of mode allows a better precision. The projection of ρ̂m is checked at every iteration, if it is too small—meaning that the modes are too poor to represent the variation of ϱ—we increase the number of modes. Thus the criterion for convergence is to look directly at the norm of ϱ or to stop when a maximal number of modes Mmax is reached.

However, it was already pointed out that, by itself, the algorithm might not converge to a solution as ϱ could be zero for another configuration. That is why, at the end of the convergence of the algorithm (or if there is no convergence at all) extra steps are added. For instance, let us say that the converged level-set is the one of Figure 6a. First, as in classical active set methods, we check if there is penetration and also that the reaction of the support is not attracting the membrane. This defines areas where the contact conditions are violated and they should be added in case there is penetration or removed for attraction. Therefore creating level-sets with the values of the gap φg Figure 6b or the reaction of the support φp Figure 6c allows us to do easy set operations. We only compute these level-sets where it is necessary, namely outside the contact zone for the gap and inside for the reaction. A special care as to be taken to avoid the creation of unnecessary iso-zero (dashed lines in Figure 6d). An easy way to do it is to shift slightly up (resp.
<figure><img src="image_10.png" /><figcaption>**Figure 6** Level-set values of the starting level-set (a), penetration (b), reaction (c), the newly computed level-set (d) and the associated signed distance function (e). Solid lines represent the iso-zero of each one.</figcaption></figure>

down) the reaction (resp. gap) in order to avoid a strict zero value on this lines. As it was said, taking the minima (resp. maxima) of two level-sets allows the creation of a new level-set which represents the union (resp. intersection) of the former level-sets. Nevertheless it was also highlighted that the property of signed distance function was lost in process. Thus we reset the level-set after such operation. This means keeping the iso-zero of the level-set then computing the distance to this front on the whole domain (cf. Figure 6e). It is to noticed that these operations are weakly dependent on the mesh. The full method is summed up in the algorithm beneath.

**Algorithm: ILS algorithm**
Initialize: level-set $\phi \leftarrow \phi_0$; number of modes $M \leftarrow M_0$;
**while** $M < M_{max}$ **do**
Solve the direct problem $[K]U = F$;
**if** criterium not small enough : $\|\rho\|_{L^2} > \rho_{min}$ **then**
**for** $m=1$ **to** $M$ **do**
Compute modes $w_m$ on $\Gamma$;
Extend $w_m$ on $\Omega$;
Solve the $m^{th}$ sensibility problem $[K]U_m^s = F_m^s$;
Solve the shape optimization problem $[K_\Gamma]\tau = \rho^\pi$;
**if** $|\tau| < \tau_{min}$ **then**
$M \leftarrow M + 1$;
update $\phi$ with $\tau$ and $w_m$;
**else criterium small enough :** $\|\rho\|_{L^2} < \rho_{min}$
**if** check for penetration and traction $\rightarrow$ OK **then**
return $1$;
**else**
compute $\phi_g$ and $\phi_l$;
update $\phi$ with $\phi_g$ and $\phi_l$;
## An easy adaptation to adhesion problems

Let us now consider the influence of adhesion in the contact zone and how it affects the Inequality Level Set approach. For instance, consider that the floor in Figure 1 is covered with a fluid with surface tension $\gamma$. At equilibrium, the configurational force on the boundary of the contact zone is no longer zero but needs to match the value of surface tension. Indeed, with adhesion, it takes energy to create new surfaces as the contact zone is reduced in size. As a consequence, at equilibrium the pressure on the boundary of the contact zone is no longer zero but we have

$$ \frac{T}{2} \rho^2 = \gamma \rightarrow \rho = \rho_a = \sqrt{\frac{2\gamma}{T}} \qquad (31) $$

The ILS algorithm described in "Method: the ILS applied to contact problems" is still the same except that the pressure is driven to $\rho_a$ and not zero.

It is interesting to note that the introduction of adhesive energy changes the solution from an inflection point to a minimum point.

Regarding the literature on adhesion, contact under adhesion has been studied in the papers [20–22] as well as the book [23] and the more recent book [24] for the mathematical aspects. On the boundary of the contact zone, pressure is no longer zero. Few papers exist in the literature on the numerical simulation of contact with adhesion. Lately, micro/nano-scale adhesion has been investigated, motivated for instance by microelectromechanical systems [25], nanoindentation [26] or biological [27, 28] purpose. Numerical improvements of classical methods have been proposed for example by Sauer [29]. Here, the macro-scale adhesion is investigated. Therefore we do not aim to represent precisely the micro-scale behavior. Note also the paper from the computational graphics community [30] which adapts classical contact algorithms to adhesion. An example of this application can be found in the next section (see Figure 7).

## Results

In this section we are going to study the case of a membrane with a homogeneous loading. Then, the full 2D problem will be dealt with.

<figure><img src="image_11.png" /><figcaption>Figure 7 Solution given by the algorithm with an without adhesion.</figcaption></figure>
## The axi-symmetric problem

We recall that the problem we are looking at is a circular membrane (Figure 8) with radius $R$ (set to 1 m for numerical purpose) and homogeneous tension $T$ (set to $1 \text{ Nm}^{-1}$). This membrane is loaded with an homogeneous loading $f_d$. The quantity of interest is the vertical displacement $u$ which is prescribed to be zero along the border. This displacement is limited by an undeformable horizontal plane of altitude $d$.

The problem can be summarized into the equations:

$$ \begin{cases} T \Delta u + f_d = 0 \quad \text{in } \Omega \\ -u \le d \quad \text{in } \Omega \\ u = 0 \quad \text{for } r = R \end{cases} \qquad (32) $$

The analytical solution of the equilibrium for a given $r_c$ is:

$$ u = -\frac{f_d r^2}{4T} + A \ln\left(\frac{r}{R}\right) + B \qquad (33) $$

where:

$$ \begin{aligned} A &= \frac{f_d(R^2 - r_c^2) - 4dT}{4T \ln(R/r_c)} \\ B &= \frac{f_d R^2}{4T} \end{aligned} \qquad (34) $$

$r_c$ is the radius of the contact zone

However, the only solution which fulfills the contact condition imposes that the exact value of $r_c$ is:

$$ r_c^{\text{ex}} = R \exp \left( \frac{1}{2} \left( 1 + W_{-1} \left( \frac{4dT - f_d}{ef_d} \right) \right) \right) \qquad (35) $$

With $W_{-1}(x)$ the branch of the Lambert function defined for $-1/e \le x < 0$ and $e = \exp(1)$. In this section, for every numerical simulation $f_d = 1 \text{ N}$ and $h = 0.1 \text{ m}$.

### With axi-symmetric hypothesis

As said before, if we take into account the axi-symmetry the problem is simplified in a 1D problem (Figure 9). At a given contact area, we need to find $u$ such that:

$$ \begin{cases} T \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r} \right) + f_d = 0 & \text{on } [0, r_c] \\ -u = d & \text{on } [r_c, R] \\ u = 0 & \text{for } r = R \end{cases} \qquad (36) $$

Let us illustrate what was described in "Method: the ILS applied to contact problems" with this simplified problem. In order to impose $u = d$ on $[0, rc]$ we use Lagrange multipliers $p$. Even if $p$ is going to be defined on $[0, R]$ we will be able to take into account their influence only on the contact zone. Indeed, taking advantage of the level-set framework we integrate only on the contact zone.
**Problem 12** Find $u \in V_d, p \in H^1(\Omega^-)$ and $\lambda \in \mathbb{R}$ such that:

$$ \left\{ \begin{array}{ll} a(u, u^*) + b(p, u^*) + c(\lambda, u^*) = f(u^*) & \forall u^* \in V_d, \\ b(p^*, u) = l_p(p^*) & \forall p^* \in H^1(\Omega^-), \\ c(\lambda^*, u) = 0 & \forall \lambda^* \in \mathbb{R}. \end{array} \right. \qquad (37) $$

with:

$$ \left| \begin{array}{ll} a(u, u^*) = \int_0^R T r \frac{\partial u}{\partial r} \frac{\partial u^*}{\partial r} d\Gamma & f(u^*) = \int_0^R f_d r u^* d\Gamma \\ b(u, p) = \int_0^{r_c} \left( l_m \frac{\partial p}{\partial r} \frac{\partial u}{\partial r} + p u \right) r d\Gamma & l_p(p^*) = - \int_0^{r_c} \left( l_m \frac{\partial p}{\partial r} \frac{\partial d}{\partial r} + p d \right) r d\Gamma \\ c(u, \lambda) = \lambda r_c \llbracket u \rrbracket & \end{array} \right. \qquad (38) $$

As we are dealing with a 1D problem we will use only one mode for $w$ which corresponds to a unit displacement of $r_c$ in $e_r$ direction. If we take the normal directional derivative of these equations we obtain the following sensitivity problem.

**Problem 13** Find $\hat{u} \in V_d, \hat{p} \in H^1(\Omega^-)$ and $\hat{\lambda} \in \mathbb{R}$ such that:

$$ \left\{ \begin{array}{ll} a(\hat{u}, \hat{u}^*) + b(\hat{p}, \hat{u}^*) + c(\hat{\lambda}, \hat{u}^*) = f^s(\hat{u}^*) & \forall \hat{u}^* \in V_d, \\ b(\hat{p}^*, \hat{u}) = l_p^s(\hat{p}^*) & \forall \hat{p}^* \in H^1(\Omega^-), \\ c(\hat{\lambda}^*, \hat{u}) = l_\lambda^s(\hat{\lambda}^*) & \forall \hat{\lambda}^* \in \mathbb{R}. \end{array} \right. \qquad (39) $$

with:

$$ \left| \begin{array}{ll} f^s(u^*) = -(r_c f_d + \lambda)u^{s+}(r_c) & \\ l_p^s(p^*) = 0 & (40) \\ l_\lambda^s(\lambda^*) = -\lambda^* r_c \frac{\lambda}{T} & \end{array} \right. $$

First we will use a linear approximation. The {$N_i$} family is therefore the family of the hat functions. The displacement field is enriched at $r_c$ with an Heaviside function. Thus the enriched nodes are the ones surrounding $r_c$ or the node itself if $r_c$ is on a node. One shall keep in mind that the integration of $p$ only on the contact zone is enabled by the

<figure><img src="image_12.png" /><figcaption>**Figure 8** Axi-symmetric membrane with a homogeneous load.</figcaption></figure>
knowledge of the position of $r_c$ between two nodes. This approximation fields are illustrated in Figure 10.

Once these two problems are solved we have access to $\rho$ and $\dot{\rho}$ and we can solve the Problem 10 associated to our case.

$$ \rho + \dot{\rho}\tau = 0 \quad \Leftrightarrow \quad \tau = -\frac{\rho}{\dot{\rho}} \qquad (41) $$

Thus if we were at the $k^{th}$ iteration the next $r_c^{k+1}$ would be:

$$ r_c^{k+1} = r_c^k + \tau = r_c^k - \frac{\rho}{\dot{\rho}} \qquad (42) $$

An illustration of the results of the ILS method can be seen in Figure 7. This same problem was also treated adding adhesion as described in "An easy adaptation to adhesion problems". We need to apply the Newton algorithm on the function $\rho - \rho_a$. The weak

<figure><img src="image_13.png" /><figcaption>**Figure 9** Axi-symmetric simplification of the membrane (right side) and illustration of the solution.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10** Approximation of the fields.</figcaption></figure>
discontinuity of the displacement is well represented by the X-FEM Heaviside enrichment combined with Lagrange multipliers.

We compared this method to an active set method with nodal Lagrange multipliers (which will be referred as “classical”) in the case without adhesion. The set where the constraint is active (namely where the contact is imposed) is a set of $N_c$ nodes. The coordinates of the $k^{th}$ node in this set is noted $\underline{x}_c^k$. In order to make this contact zone evolve, the penetration and the reaction of the support are checked as explained in the last section. Then, every node which violates the penetration is added in the active set whereas the ones violating the support reaction are removed. The same discretization of $u$ is used without the enrichment. The problem for a given $\Omega^-$ is then:

**Problem 14** Find $u \in \mathcal{V}_d, p \in \mathbb{R}^{N_c}$ such that:

$$ \left\{ \begin{array}{l} a(u, u^*) + \sum_{k=1}^{N_c} p_k u^*(\underline{x}_c^k) = f(u^*) \quad \forall u^* \in \mathcal{V}_d, \\ \\ \sum_{k=1}^{N_c} p_k^* u(\underline{x}_c^k) = - \sum_{k=1}^{N_c} p_k^* d(\underline{x}_c^k) \quad \forall p^* \in \mathbb{R}^{N_c}. \end{array} \right. \qquad (43) $$

The initial $r_c$ is chosen as $r_{c0} = 0.3m$. In Figure 11a we are interested in the evolution of the relative error on the position of $r_c$ with the size of the elements $h$. With ILS method we observe a better convergence rate than with the other one. The stages of the other method are also avoided. They arise because as long as the refinement of the mesh does not add a node closer to $r_c^{ex}$ there is no better approximation available. This problem is avoided by the enrichment of the ILS method. Figure 11b on the other hand shows the power of having an evolution driven by the sensitivity. Indeed, by adding more degree of freedom we increase the number of iterations needed with the classical active set as there is more sets possibilities. Whereas it is not increasing in our case as it is sensitivity driven.

### Without axi-symmetric hypothesis

Without this hypothesis, we recall that the direct problem is Problem 5. Then, if we use the modal base $\{w_m\}$ to approximate the displacement field and note $\underline{w}_m = w_m \underline{n}$, the sensitivity problem becomes similar to Problem 13 with the above notation and:

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11** Evolution of **a** the error on $r_c$, **b** the number of iterations, with the size of the mesh.</figcaption></figure>
$$
\begin{align}
f^s(u^*, w) &= -\int_{\Omega} T (\nabla u^T \nabla u^* \operatorname{div} \underline{w} - (\nabla u^T (\nabla \underline{w}^T + \nabla \underline{w}) \nabla u^*)) \, d\Omega \nonumber \\
&\quad -\int_{\Omega-} l_m (\nabla p^T \nabla u^* \operatorname{div} \underline{w} - (\nabla p^T (\nabla \underline{w}^T + \nabla \underline{w}) \nabla u^*)) \, d\Omega \nonumber \\
&\quad -\int_{\Omega-} p u^* \operatorname{div} \underline{w} \, d\Omega - \int_{\Gamma} \lambda \llbracket u^* \rrbracket (\operatorname{div} \underline{w} - \underline{n} \cdot \nabla \underline{w} \cdot \underline{n}) \, d\Gamma + \int_{\Omega} f_d u^* \operatorname{div} \underline{w} \, d\Omega \nonumber \\
l_p^s(p^*, w) &= 0 \nonumber \\
l_\lambda^s(\lambda^*, w) &= -\int_{\Gamma} \lambda^* \llbracket u \rrbracket (\operatorname{div} \underline{w} - \underline{n} \cdot \nabla \underline{w} \cdot \underline{n}) \, d\Gamma \tag{44}
\end{align}
$$

In order to improve the efficiency of the algorithm a limitation is added to the basic algo-
rithm described Figure 5. First, as the starting position of the interface can be far from the
converged one, the evolution of Γ is limited such that it does not evolve further than the
limit set by the value of δ (which defined the restrained computational area specified in
Figure 4).

The primary task is to set the initial contact area boundary to be the solution of the problem without any contact conditions. Therefore, the first guessed contact area is a centered circle of radius 0.77 m. Only the first mode should be activated and we should observe a proportional growth of the contact zone. Indeed, the first mode represent an homogeneous growth of the front. The combination of the second and the third allows a translation of the contact zone. Latter modes become more and more localize. A first simulation is done with a non axi-symmetric mesh. The result can be seen Figure 12. The jump in the membrane slope, clearly visible in the initial solution, is canceled at the end of the iterations. Here, a small δ was taken on purpose. This explained the constant advance of the front for the first iterations.

In Figure 13a it can be seen that other modes than the first one are activated. This is
due to the numerical error during the level-set evolution which breaks the axi-symme-
try. Therefore the algorithm needs extra modes to compensate these errors. To ensure
that this is due to these numerical aspects an axi-symmetric mesh was used in another
simulation. Figure 13b confirmed that no other modes but the first one are activated at
first. The activation of the other modes at the 8th iteration is just due to the fact that the
convergence stage as been reached as Figure 14 shows.

A stage in the convergence can be observe after some iterations (cf. Figure 14). It is due
to the geometrical approximation of the level set on the mesh. Therefore, it is decreasing
with the size of the mesh as we get a better representation of the level-set.

<figure><img src="image_17.png" /><figcaption>**Figure 12** Evolution of the contact interface from initial to final configuration (b) with the displacement associated to these to extreme position (a, c).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 13** Evolution of the norm of $\tau$ and percentage of the contribution of each mode.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14** Error on the position of the contact interface compare to a reference solution for each problem.</figcaption></figure>

Let us compare again the ILS method to the active set method. The behavior of these methods is quite similar to the axi-symmetric hypothesis case, retrieving the same properties of convergence for the location of $\Gamma$ (cf. Figure 15a).

In relevance to the efficiency of the algorithm, Table 1 gives the total computational cost for each method up to convergence as well as the average time per iteration. First

<figure><img src="image_20.png" /><figcaption>**Figure 15** Evolution of **a** the error on $\Gamma$, **b** the number of iterations, with the size of the mesh.</figcaption></figure>
<table><caption>Table 1 Time comparison between ILS and active set method. First on fixed meshes then for a targeted error</caption><thead><tr><th rowspan="2">Mesh<br>1/h</th><th colspan="2">ILS</th><th colspan="2">Classical</th></tr><tr><th>t<sup>ILS</sup><sub>tot</sub> (s)</th><th>t<sup>ILS</sup><sub>iter</sub> (s)</th><th>t<sup>class</sup><sub>tot</sub> (s)</th><th>t<sup>class</sup><sub>iter</sub> (s)</th></tr></thead><tbody><tr><td>2</td><td>0.4</td><td>0.08</td><td>0.06</td><td>0.02</td></tr><tr><td>4</td><td>1.06</td><td>0.212</td><td>0.39</td><td>0.08</td></tr><tr><td>8</td><td>2.53</td><td>0.506</td><td>1.06</td><td>0.151</td></tr><tr><td>16</td><td>7.59</td><td>1.52</td><td>5.24</td><td>0.476</td></tr><tr><td>32</td><td>36.4</td><td>7.28</td><td>48.7</td><td>2.43</td></tr><tr><td>64</td><td>243.9</td><td>48.78</td><td>617.8</td><td>19.92</td></tr><tr><td>128</td><td>2,062</td><td>412.4</td><td>20,971</td><td>349.5</td></tr></tbody></table><table><thead><tr><th rowspan="2">Error on Γ</th><th colspan="2">ILS</th><th colspan="2">Classical</th></tr><tr><th>t<sup>ILS</sup><sub>tot</sub> (s)</th><th>t<sup>ILS</sup><sub>iter</sub> (s)</th><th>t<sup>class</sup><sub>tot</sub> (s)</th><th>t<sup>class</sup><sub>iter</sub> (s)</th></tr></thead><tbody><tr><td>5 × 10<sup>-2</sup></td><td>0.4</td><td></td><td>1.06</td><td></td></tr><tr><td>5 × 10<sup>-3</sup></td><td>2.53</td><td></td><td>617.8</td><td></td></tr><tr><td>5 × 10<sup>-4</sup></td><td>36.4</td><td></td><td>×</td><td></td></tr><tr><td>10<sup>-2</sup></td><td>2,062</td><td></td><td>×</td><td></td></tr></tbody></table>

All computations have been done on an AMD Opteron 6328 3.2 GHz octocore with 4Go of RAM.

the comparison is made for several fixed meshes and then for a targeted error on the location of $\Gamma$. First, let us discuss the computational times on fixed meshes. The ILS uses more time on coarse mesh than the classical method. Nonetheless for regular sized or fine meshes this trend is reversed. This is due to both the increase of iterations needed by the classical method and the reduction of the difference of time needed by each method for one iteration. This last point is explained by the change in the dominant time consuming function. For a coarse mesh, the assembling takes the greatest computational effort. Whereas, for a fine mesh, the solving is the main expense. Also, the ILS algorithm is more elaborate. Thus, it takes more efforts to be optimized and we think there is still room for improvement. But the comparison should also be made with a precision criterion on $\Gamma$. Indeed, we have already emphasized the importance of the contact boundary. For this particular purpose, the ILS is far more efficient and does not require an overkill mesh to obtain sufficient accuracy for the contact boundary position.

## Non-axi-symmetric problems

In the previous section, we have dealt with axi-symmetric problems. The ILS has been shown to be consistent when keeping the axi-symmetric features. However in these cases, only the first mode was needed to converge. We will now seek examples which need extra modes. In this part, we are going to see to the resolution of a membrane coming into contact with a non-constant floor following the mathematical expression $d(x,y) = -1/1,000(15 + 5 \cos(2\pi x) \cos(2\pi y))$. A representation of the floor and the used mesh can be found in Figure 16.

The loading will be increased and then decreased such that contact zones appear, then merge and split ($|f_d| = 0.1 \rightarrow 0.2 \rightarrow 0.3 \rightarrow 0.4 \rightarrow 0.3 \rightarrow 0.2 \rightarrow 0.1$). Results of this process can be found in Figure 17.
<figure><img src="image_21.png" /><figcaption>**Figure 16** Floor shape (2D and 3D representation) and associated mesh.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

a $f = -0.1, \text{RT, 3 ILS}$

<figure><img src="image_23.png" /></figure>

b $f_d = -0.2, \text{RT, AP, 5 ILS}$

<figure><img src="image_24.png" /></figure>

c $f_d = -0.3, \text{AP, 2 ILS}$

<figure><img src="image_25.png" /></figure>

d $f_d = -0.4, \text{AP, 2 ILS}$

<figure><img src="image_26.png" /></figure>

e $f_d = -0.3, \text{RT, 2 ILS}$

<figure><img src="image_27.png" /></figure>

f $f_d = -0.2, \text{RT, AP, 5 ILS}$

<figure><img src="image_28.png" /></figure>

g $f_d = -0.1, \text{RT, 4 ILS}$

<figure><img src="image_29.png" /><figcaption>**Figure 17** Evolution of the solution and the contact boundary (green line) with the loading (AP meaning addition of the area with penetration to the contact zone, RT meaning removal of the area where there is traction on the membrane and ILS is short for ILS iteration).</figcaption></figure>
The method seems to be consistent as the results are similar for the loading and the unloading of the membrane. The differences come from the different initial level-set between the loading and the unloading. Indeed, using Fourier modes, their number is limited by the number of elements cut by $\Gamma$ and therefore local errors (mainly arising when penetration or release checks are done) are difficult to overcome. In a future work localized modes will be developed to deal with this issue. Nevertheless, it is notable that few iterations are needed between each steps.

## Conclusion

The use of level-sets coupled with X-FEM seems to be a great asset to handle contact problems. The two key advantages are the capacity to represent:

- precisely the contact zone without mesh constraints,

- the low regularity of the solution on the boundary of the contact zone.

This leads to a good quality solution without use of an excessively refined mesh. On top of it, the ILS takes advantage of configurational mechanics to make the contact zone evolve. It has been shown in this paper to be an efficient way to find the contact zone with and without adhesion. Numerical experiments did show that the convergence is fast.

In a future work, we will investigate the case of plain solids, particularly punch problems. Also the ILS framework offers a natural way to numerically analyze the bifurcation of systems involving contact and adhesion. An example of such system is given in [31].

## Abbreviations

ILS: inequality Level-Set; X-FEM: eXtended Finite Element Method; BB condition: Babuška-Brezzi condition.

## Authors' contribution

MG worked on the algorithms, performed the computations and drafted the manuscript. NC worked on the algorithms, performed the computations and carried out detailed revision. NM worked on the algorithms and carried out detailed revision. All authors read and approved the final manuscript.

## Acknowledgements

The support of the ERC Advanced Grant XLS no 291102 is greatly acknowledged. We would also like to thank Peter Wriggers for his advices and Anthony Nouy for his help on the mathematical formulation.

## Compliance with ethical guidelines

## Competing interests

The authors declare that they have no competing interests.

## Appendix A: calculus of the sensibility

Lets first recall some basic formulas:

$$ \bar{\nabla} \dot{u} = \nabla \dot{u} + \nabla \nabla u \cdot w = \nabla(\dot{u} + \nabla u \cdot w) - \nabla u \cdot \nabla w = \nabla \dot{u} - \nabla u \nabla w \quad (45) $$

$$ D\left(\int_{\Omega} f \, d\Omega\right)[\underline{w}] = \int_{\Omega} D(f)[\underline{w}] + f \operatorname{div}\underline{w} \, d\Omega \quad (46) $$

$$ D\left(\int_{\Gamma} f \, d\Gamma\right)[\underline{w}] = \int_{\Gamma} D(f)[\underline{w}] + f(\operatorname{div}\underline{w} - \underline{n} \cdot \nabla \underline{w} \cdot \underline{n}) \, d\Gamma \quad (47) $$
We have

$$ 
\begin{aligned}
D(a(u, u^*))[\underline{w}] &= \int_{\Omega} T(\dot{\bar{\nabla}}u \cdot \nabla u^* + \nabla u \cdot \dot{\bar{\nabla}}u^* + \nabla u \cdot \nabla u^* \operatorname{div}_{\underline{w}}) \, d\Omega \\
&= \int_{\Omega} T((\nabla \dot{u} - \nabla u \cdot \nabla_{\underline{w}}) \cdot \nabla u^* + \nabla u \cdot (\nabla \dot{u}^* - \nabla_{\underline{w}}\nabla u^*) + \nabla u \cdot \nabla u^* \operatorname{div}_{\underline{w}}) \, d\Omega
\end{aligned}
\quad (48) $$

$$ 
\begin{aligned}
D(b(p, u^*))[\underline{w}] &= \int_{\Omega^{-}} l_m(\nabla \dot{p} \nabla u^* + \nabla p \nabla \dot{u}^* - \nabla p(\nabla \underline{w}^T + \nabla \underline{w})\nabla u^* + \nabla p \nabla u^* \operatorname{div}_{\underline{w}}) \, d\Omega \\
&\quad + \int_{\Omega^{-}} \dot{p}u^* + pu^* + pu^* \operatorname{div}_{\underline{w}} \, d\Omega
\end{aligned}
\quad (49) $$

$$ D(c(\lambda, u^*))[\underline{w}] = \int_{\Gamma} \dot{\lambda}[\underline{u}^*] + \lambda[\underline{u}^*] + \lambda[\underline{u}^*](\operatorname{div}_{\underline{w}} - \underline{n} \cdot \nabla_{\underline{w}} \cdot \underline{n}) \, d\Gamma \quad (50) $$

$$ D(f(u^*))[\underline{w}] = \int_{\Omega} \dot{f}u^* + fu^* + fu^* \operatorname{div}_{\underline{w}} \, d\Omega \quad (51) $$

$$ 
\begin{aligned}
D(l_p(p^*))[\underline{w}] &= - \int_{\Omega^{-}} l_m(\nabla \dot{p}^* \nabla d + \nabla p^* \nabla \dot{d} - \nabla p^*(\nabla \underline{w}^T + \nabla \underline{w})\nabla d + \nabla p^* \nabla d \operatorname{div}_{\underline{w}}) \, d\Omega \\
&\quad - \int_{\Omega^{-}} \dot{p}^* d + p^* \dot{d} + p^* d \operatorname{div}_{\underline{w}} \, d\Omega
\end{aligned}
\quad (52) $$

The test fields $u^*$, $p^*$ and $\dot{u}^*$ are chosen such that for all $\underline{x}$ where the fields are defined, $\dot{u}^*(\underline{x}) = 0$, $\dot{p}^*(\underline{x}) = 0$ and $\dot{\lambda}^*(\underline{x}) = 0$. This means that the test fields follow the movement of the interface. Therefore if the imposed force field $f_d$ to also follow the movement of the interface (i.e. $\dot{f}_d=0$) we find the formulation of Problem 8.

Received: 28 November 2014 Accepted: 2 June 2015

Published online: 15 July 2015

References

1. Zienkiewicz OC, Taylor RL (2005) The finite element method for solid and structural mechanics. Butterworth-Heinemann

2. Temizer I, Wriggers P, Hughes TJR (2011) Contact treatment in isogeometric analysis with NURBS. Comput Methods Appl Mech Eng 200:1100–1112

3. Belgacem FB, Hild P, Laborde P (1998) The mortar finite element method for contact problems. Math Comput Model 28(4):263–271

4. Popp A, Wohlmuth B, Gee M, Wall W (2012) Dual quadratic mortar finite element methods for 3d finite deformation contact. SIAM J Sci Comput 34:B421–B446

5. Kikuchi N, Oden JT (1988) Contact problems in elasticity: a study of variational inequalities and finite element methods, vol 8. Siam

6. Wriggers P (2006) Computational contact mechanics, vol 30167. Springer, New York

7. Alart P, Curnier A (1991) A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput Methods Appl Mech Eng 92(3):353–375

8. Belytschko T, Parimi C, Moës N, Sukumar N, Usui S (2003) Structured extended finite element methods for solids defined by implicit surfaces. Int J Numer Methods Eng 56:609–635

9. Simo JC, Wriggers P, Taylor RL (1985) A perturbed Lagrangian formulation for the finite element solution of contact problems. Comput Methods Appl Mech Eng 50(2):163–180
10. Conry TF, Seireg A (1971) A mathematical programming method for design of elastic bodies in contact. J Appl Mech 38(2):387-392

11. Laursen TA (2002) Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis. Springer, New York

12. Franke D, Rank E, Düster A (2011) Computational contact mechanics based on the rp-version of the finite element method. Int J Comput Methods 08:493-512

13. Franke D, Düster A, Nübel V, Rank E (2010) A comparison of the h-, p-, hp-, and rp-version of the FEM for the solution of the 2d Hertzian contact problem. Comput Mech 45:513-522

14. Bonfils N, Chevaugeon N, Moës N (2012) Treating volumetric inequality constraint in a continuum media with a coupled X-FEM/level-set strategy. Comput Methods Appl Mech Eng 205-208:16-28

15. Osher S, Fedkiw R (2003) Level set methods and dynamic implicit surfaces. Springer, New York

16. Sethian JA (1996) A fast marching level set method for monotonically advancing fronts. Proc Natl Acad Sci 93(4):1591-1595

17. Moës N, Béchet E, Tourbier M (2006) Imposing Dirichlet boundary conditions in the extended finite element method. Int J Numer Methods Eng 67(12):1641-1669

18. Taroco E (2000) Shape sensitivity analysis in linear elastic fracture mechanics. Comput Methods Appl Mech Eng 188:697-712

19. Pradeilles-Duval R-M, Stolz C (1995) Mechanical transformations and discontinuities along a moving surface. J Mech Phys Solids 43:91-121

20. Johnson KL, Kendall K, Roberts AD (1971) Surface energy and the contact of elastic solids. Proc R Soc Lond A Math Phys Sci 324(1558):301-313

21. Derjaguin BV, Muller VM, Toporov YP (1975) Effect of contact deformations on the adhesion of particles. J Colloid Interface Sci 53:314-326

22. Fremond M (1982) Adhérence des solides. J De mécanique Théorique et Appliquée 6(3):383-407

23. Johnson KL (1987) Contact mechanics. Cambridge University Press, Cambridge

24. Sofonea M, Han W, Shillor M (2010) Analysis and approximation of contact problems with adhesion or damage. CRC Press, London

25. Zhao YP, Wang LS, Yu TX (2003) Mechanics of adhesion in MEMS-a review. J Adhes Sci Technol 17:519-546

26. Fischer-Cripps AC (2011) Nanoindentation. Springer Science & Business Media

27. Chu Y-S, Dufour S, Thiery JP, Perez E, Pincet F (2005) Johnson-Kendall-Roberts theory applied to living cells. Phys Rev Lett 94:028102

28. Autumn K, Sitti M, Liang YA, Peattie AM, Hansen WR, Sponberg S et al (2002) Evidence for van der Waals adhesion in gecko setae. Proc Natl Acad Sci 99(19):12252-12256

29. Sauer RA (2011) Enriched contact finite elements for stable peeling computations. Int J Numer Methods Eng 87(6):593-616

30. Gascón J, Zurdo JS, Otaduy MA (2010) Constraint-based Simulation of Adhesive Contact. In: Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA '10. Eurographics Association, (Aire-la-Ville, Switzerland, Switzerland), pp 39-44

31. Neukirch S, Roman B, de Gaudemaris B, Bico J (2007) Piercing a liquid surface with an elastic rod: Buckling under capillary forces. J Mech Phys Solids 55(6):1212-1235

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Immediate publication on acceptance

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com