Skip to main content

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

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.

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. 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 \(\underline{x}=(x,y,z)\) is the coordinate of a point. We will denote \(u(\underline{x})\) the displacement of the point \(\underline{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( \underline{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 \(\overline{\Omega }\) then the domains must satisfy \(\overline{\Omega }^+\cup \overline{\Omega }^-=\overline{\Omega }\) and \(\Omega ^+ \cap \Omega ^-=0\). We will define on every point of the boundary \(\Gamma\) the normal \(\underline{n}\) to it contained in the membrane surface. It is chosen to be directed towards \(\Omega ^+\). Figure 1 illustrates this description.

Figure 1
figure 1

Membrane into contact with a rigid surface.

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{aligned} \left\{ \begin{array}{l} T\Delta u +f_d=0 \quad \text {on } \Omega ,\\ u = 0\quad \text {on }\partial \Omega. \end{array}\right. \end{aligned}$$
(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{aligned} \left\{ \begin{array}{l} T\Delta u +f_d+\tilde{p} =0 \quad \it {on} \;\Omega ,\\ u = 0\quad {on } \; \partial \Omega. \end{array}\right. \end{aligned}$$
(2)
$$\begin{aligned} \left\{ \begin{array}{rcl} u+d\,\ge \,0 \; { on } \;\Omega \, ,\\ \tilde{p} \,\ge \,0 \; { on } \; \Omega \, ,\\ \tilde{p}(u+d) \,=\,0\; { on }\; \Omega. \end{array}\right. \end{aligned}$$
(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:

 

In \(\Omega ^+\)

In \(\Omega ^-\)

On \(\Gamma\)

Equality

\(\begin{aligned}T\Delta u +f_d =0 \\ u=0 \;{on} \;\partial \Omega \end{aligned}\)

\(\begin{aligned} &T\Delta u +f_d+\tilde{p} =0\\ &-u =d \end{aligned}\)

\(u^+=u^-\)

Inequality

\(-u\le d\)

\(\tilde{p}\ge 0\)

 

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:

 

In \(\Omega ^+\)

In \(\Omega ^-\)

On \(\Gamma\)

Equality

\(\begin{aligned} T\Delta u +f_d =0 \\ u=0 \;on\; \partial \Omega \end{aligned}\)

\(\begin{aligned} &T\Delta u +f_d+\tilde{p} =0 \\ & -u =d \end{aligned}\)

\(\begin{aligned} &T[\![\nabla u]\!]\cdot \underline{n}=\tilde{\lambda }\\& u^+=u^- \end{aligned}\)

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 \underline{n}=0\). First, let us look at the case \([\![\nabla u]\!]\cdot \underline{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 \underline{n}< 0\) and therefore \(\nabla u^+(s)\cdot \underline{n}<\nabla u^-(s)\cdot \underline{n}\). As by definition \(\nabla d(s)\) is continuous and \(\nabla u^-(s)=\nabla d(s)\) we must have \(\nabla u^+(s)\cdot \underline{n}<\nabla d(s)\cdot \underline{n}\) and therefore there is penetration of the rigid surface.

Hence,

$$\begin{aligned}{}[\![\nabla u]\!]\cdot \underline{n} =0 \quad \forall s\in \Gamma \end{aligned}$$
(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:

$$\begin{aligned} \mathcal {V}_d = \left\{ u\in H^1(\Omega ^+ \cup \Omega ^-) | u=0 \text { on } \partial \Omega \right\} \end{aligned}$$
(5)

Therefore if \(u \in \mathcal {V}_d\) it can admit a strong discontinuity on \(\Gamma\). Problem 3 can be written as the weak formulation:

Problem 4

Find \(u\in \mathcal {V}_d\), \(p \in P\) and \(\lambda \in H^{-1/2}(\Gamma )\) such that:

$$\begin{aligned} \left\{ \begin{array}{l} a(u,u^{\!\star })+\langle p,u^{\!\star }\rangle _P+\langle \lambda , [\![u^{\!\star }]\!]\rangle _{H^{ -1/2}}= f(u^{\!\star })\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ \langle p^{\!\star },u+d\rangle _P = 0 \quad \forall p^{\!\star }\in P,\\ \langle \lambda ^{\!\star }, [\![u ]\!]\rangle _{H^{ -1/2}} = 0 \quad \forall \lambda ^{\!\star }\in H^{-1/2}(\Gamma ). \end{array}\right. \end{aligned}$$
(6)

with:

(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}=\displaystyle \int _{\Omega ^-}\!\!\!pu\,\text {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}=\displaystyle \int _{\Omega ^-}\!\!\!l_{\text {m}}\nabla p \cdot \nabla u\,\text {d}\Omega +\int _{\Omega ^-}\!\!\!pu\,\text {d}\Omega\). The addition of \(l_{\text {m}}\) is only there for dimensional homogenization and will be set to \(l_{\text {m}}=1\text {m}^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 \mathcal {V}_d\), \(p \in H^1(\Omega ^-)\) and \(\lambda \in H^{-1/2}(\Gamma )\) such that:

$$\begin{aligned} \left\{ \begin{array}{l} a(u,u^{\!\star })+b(p,u^{\!\star })+c(\lambda ,u^{\!\star })= f(u^{\!\star })\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ b( p^{\!\star },u) = l_p( p^{\!\star }) \quad \forall p^{\!\star }\in H^1(\Omega ^-),\\ c( \lambda ^{\!\star },u) = 0 \quad \forall \lambda ^{\!\star }\in H^{-1/2}(\Gamma ). \end{array}\right. \end{aligned}$$
(8)

with:

(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 \(\Omega ^-\) 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 \(\Gamma\). In order to do so , we are going to use a level-set [15, 16] and particularly its iso-zero set to fit \(\Gamma\) (cf. Figure 2b). A particular case of level-set is a signed distance function \(\phi (\underline{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 \(\Omega ^-\). Furthermore, it is easy to compute set theory operations such as the union 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.

Figure 2
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).

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:

$$\begin{aligned} \phi (\underline{x})=\sum _{i\in I} \phi _i N_i(\underline{x}) \end{aligned}$$
(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:

$$\begin{aligned} u_h(\underline{x}) = \sum _{i\in I} u_iN_i(\underline{x})+\sum _{j\in J} u^e_jN_j(\underline{x})F(\underline{x}) \end{aligned}$$
(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\):

$$\begin{aligned} p_h(\underline{x}) = \sum _{i\in I^p} p_iN^p_i(\underline{x}) \end{aligned}$$
(12)
$$\begin{aligned} \lambda _h(\underline{x}) = \sum _{i\in I^{\lambda }} \lambda _iN^{\lambda }_i(\underline{x}) \end{aligned}$$
(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 matricial form:

Problem 6

Find \(\mathbf {u}\), \(\mathbf {u}^e\), \(\varvec{p}\) and \(\bf {\varvec{\lambda }}\) such that:

$$\begin{aligned} \left[ {\mathbf {K}} \right] \left\{ {\begin{array}{*{20}c} {\mathbf {u}} \\ {{\mathbf {u}}^{e} } \\ {\mathbf {p}} \\ {\varvec {\lambda }} \\ \end{array} } \right\} = \left\{ {\begin{array}{*{20}c} {\mathbf {f}} \\ {{\mathbf {f}}^{e} } \\ {{\mathbf {l}}_{p} } \\ {\mathbf {0}} \\ \end{array} } \right\} \quad {\text {where }}\left[ {\mathbf {K}} \right] = \left[ {\begin{array}{*{20}c} {{\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}} }&{} {\mathbf {0}} &{} {\mathbf {0}}\\ {{\mathbf {C}}^{T} } &{{\mathbf {C}}^{{eT}} }&{} {\mathbf {0}} &{} {\mathbf {0}} \\ \end{array} } \right] \end{aligned}$$
(14)

with:

(15)

The size of the matrix \(\mathbf {A}^{\!\!ee}\) is usually very small compare to the size of \(\mathbf {A}^{\!\!cc}\). 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}^{\!\!cc}\). Therefore a great deal of computational effort can be spared by working on \(\mathbf {A}^{\!\!cc}\) 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 \(\varrho\) 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 usually only write \(\Gamma\). The abovementioned criterion must satisfy:

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

where \(\Gamma ^{\text {ex}}\) is the exact boundary of the contact zone. Thus if \(\varrho \ne 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 \(\varrho\) is \([\![\nabla u ]\!]\cdot \underline{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 \(\varrho\) 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[\,\,]\). 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 \(\underline{w}_{\!\Gamma }=w_{\!\Gamma } \underline{n}\). Then we denote \(\Gamma (0)\) the original position of the boundary and \(\Gamma (\zeta\underline{w}_{\!\Gamma })\) the new position after a displacement \(\underline{w}_{\!\Gamma }\) (Figure 3). The normal directional derivative of \(\varrho\) is then defined as:

$$\begin{aligned} D\varrho [\underline{w}_{\!\Gamma }]=\lim _{\zeta \rightarrow 0}\frac{\varrho (\Gamma (\zeta \underline{w}_{\!\Gamma }))-\varrho (\Gamma (0))}{\zeta } \end{aligned}$$
(17)
Figure 3
figure 3

Evolution of the boundary following a path \(\underline{w}_{\!\Gamma }\).

Then the zero of \(\varrho\) can be found by solving:

Problem 7

Find \(\underline{w}_{\!\Gamma } \in L^2(\Gamma )\) such that:

$$\begin{aligned}&\varrho +D\varrho [\underline{w}_{\!\Gamma }]=0 \text { on } \Gamma \end{aligned}$$
(18)
$$\begin{aligned} \Leftrightarrow&\int _{\Gamma }\varrho \underline{ w}_{\!\Gamma }^{\star }\,\text {d}\Gamma +\int _{\Gamma }D\varrho [\underline{w}_{\!\Gamma }]\underline{w}_{\!\Gamma }^{\star }\,\text {d}\Gamma =0 \quad \forall \,\underline{w}_{\!\Gamma }^{\star }\in L^2(\Gamma ) \end{aligned}$$
(19)

We will denote the normal directional derivative of a quantity a: \(D(a)[\underline{w}_{\!\Gamma }]=\mathring{a}\). In this particular case, \(\mathring{\varrho }= \mathring{\overline{[\![\nabla u ]\!]\cdot \underline{n}}}=[\![\nabla \mathring{u} ]\!]\cdot \underline{n}-[\![\nabla u ]\!]\cdot \nabla \underline{w}_{\!\Gamma }\cdot \underline{n} +[\![\nabla u ]\!]\cdot \nabla \underline{n}\cdot \underline{w}_{\!\Gamma }\). 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 \((\mathring{u},\mathring{p},\mathring{\lambda })\) as unknowns. In order to compute \(\mathring{u}\), \(\mathring{p}\) and \(\mathring{\lambda }\) we need to extend the definition of \(D[\,\,]\) to the whole domain. The boundary displacement \(\underline{w}_{\!\Gamma }\) is then extended onto this domain. This extended displacement is noted \(\underline{w}\). But, as we are only, in practice, interested by the evolution of \(\varrho\) 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.

$$\begin{aligned} \underline{w}(\underline{x})=w_{\Gamma }(s_c)\nabla \phi (\underline{x})r(\varsigma ) \end{aligned}$$
(20)

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

Figure 4
figure 4

Function r in the case of a circular level-set.

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.

$$\begin{aligned} Df[\underline{w}]=\lim _{\zeta \rightarrow 0}\frac{f(\Gamma (\zeta \underline{w}))-f(\Gamma (0))}{\zeta } \end{aligned}$$
(21)

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

Problem 8

Find \(\mathring{u}\in \mathcal {V}_d\) , \(\mathring{p}\in H^1(\Omega ^-)\) and \(\mathring{\lambda }\in H^{-1/2}(\Gamma )\) such that:

$$\begin{aligned} \left\{ \begin{array}{l} a(\mathring{u},u^{\!\star })+b(\mathring{p},u^{\!\star })+c(\mathring{\lambda },u^{\!\star })= f^s(u^{\!\star },w)\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ b( p^{\!\star },\mathring{u}) = l^s_p( p^{\!\star },w) \quad \forall p^{\!\star }\in L^2(\Omega ^-),\\ c( \lambda ^{\!\star },\mathring{u}) = l^s_{\lambda }( \lambda ^{\!\star },w) \quad \forall \lambda ^{\!\star }\in H^{-1/2}(\Gamma ). \end{array}\right. \end{aligned}$$
(22)

where:

(23)

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 \(\mathring{\mathbf {u}}\), \(\mathring{\varvec{p}}\) and \(\mathring{\varvec{\lambda }}\) such that:

$$\begin{aligned} \left[ \mathbf {K}\right] \left\{ \begin{array}{l} \mathring{\mathbf {u}} \\ \mathring{\varvec{p}}\\ \mathring{\varvec{\lambda }} \end{array}\right\} =\left\{ \begin{array}{l} \mathbf {f}^s \\ \mathbf {l}_p^s\\ \mathbf {l}_{\lambda }^s \end{array}\right\} \end{aligned}$$
(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_{\!\Gamma }\) is unknown, we cannot directly compute Problem 9. Therefore, \(w_{\!\Gamma }\) is going to be approximated by M given modes \(w_m\). For each mode, a \(\tau _m\) is associated and represent the \(m{\text {th}}\) mode contribution to \(w_{\!\Gamma }\).

$$\begin{aligned} w_{\!\Gamma }=\sum _m\tau _mw_m. \end{aligned}$$
(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 \(\mathring{\varrho }_m\) and then solve Problem 10.

Problem 7 can thereby be approximated using a Galerkin approximation:

Problem 10

Find \(\varvec{\tau }=\{\tau _m\} \in \mathbb {R}^M\) such that:

$$\begin{aligned}&\sum _n \int _{\Gamma }\varrho w_n\,\text {d}\Gamma +\sum _m\sum _n\int _{\Gamma }\tau _m D\varrho [w_m\underline{n}]w_n\,\text {d}\Gamma =0 \end{aligned}$$
(26)
$$\begin{aligned} \Leftrightarrow&\sum _n \int _{\Gamma }\varrho w_n\,\text {d}\Gamma +\sum _m\sum _n\int _{\Gamma }\tau _m \mathring{\varrho }_mw_n\,\text {d}\Gamma =0 \end{aligned}$$
(27)

Again, Problem 10 can be written as:

Problem 11

Find \(\varvec{\tau } \in \mathbb {R}^m\) such that:

$$\begin{aligned} \left[ \mathbf {K}_{\Gamma }\right] \varvec{\tau }=\varvec{\varrho }^{\pi } \end{aligned}$$
(28)

where:

$$\begin{aligned} K_{nm}=\sum _{i}\int _{\Gamma }\varrho _iDN^{\lambda }_i[w_m]w_n\,\text {d}\Gamma \end{aligned}$$
(29)
$$\begin{aligned} \varrho ^{\pi }_{n}=\sum _{i}\int _{\Gamma }\varrho _i N^{\lambda }_iw_n\,\text {d}\Gamma \end{aligned}$$
(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 \(\Gamma\) we need to update the level-set using the displacement \(\underline{w}\). In order to do so, an Hamilton–Jacobi equation (\(\partial \phi /\partial t=\underline{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 \(\varrho\) noted \(\pi (\varrho _m)\) onto the mode space is small enough (i.e. inferior to a scalar \(\pi _{\text {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 \(\mathring{\varrho }_m\) is checked at every iteration, if it is too small—meaning that the modes are too poor to represent the variation of \(\varrho\)—we increase the number of modes. Thus the criterion for convergence is to look directly at the norm of \(\varrho\) or to stop when a maximal number of modes \(M^{\text {max}}\) is reached.

Figure 5
figure 5

Iterative algorithm to find \(\Gamma\) for a given load and a fixed number of mode.

However, it was already pointed out that, by itself, the algorithm might not converge to a solution as \(\varrho\) 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 \(\phi _g\) Figure 6b or the reaction of the support \(\phi _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. 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.

Figure 6
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.

figure a

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

$$\begin{aligned} \frac{T}{2} \varrho ^2 = \gamma \rightarrow \varrho = \varrho _a = \sqrt{\frac{2 \gamma }{T}} \end{aligned}$$
(31)

The ILS algorithm described in "Method: the ILS applied to contact problems" is still the same except that the pressure is driven to \(\varrho _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 [2022] 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).

Figure 7
figure 7

Solution given by the algorithm with an without adhesion.

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.

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\text { 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.

Figure 8
figure 8

Axi-symmetric membrane with an homogeneous load.

The problem can be summarized into the equations:

$$\begin{aligned} \left\{ \begin{array}{rll} T\displaystyle \Delta u + f_d\,=0 \; \text {in } \Omega \\ - u \,\le d \; \text {in } \Omega \\ u \, =0 \; \text {for } r=R\\ \end{array}\right. \end{aligned}$$
(32)

The analytical solution of the equilibrium for a given \(r_c\) is:

$$\begin{aligned} u=-\frac{f_dr^2}{4T}+A\ln \left( \frac{r}{R}\right) +B \end{aligned}$$
(33)

where:

(34)

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

$$\begin{aligned} r_c^{\text {ex}}=R\exp \left( \frac{1}{2}\left( 1+\text {W}_{-1}\left( \frac{4 d T-f_d}{e f_d}\right) \right) \right) \end{aligned}$$
(35)

With \(\text {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{aligned} \left\{ \begin{array}{rll} T\displaystyle \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 } [rc,R]\\ u =0 &{} \text {for } r=R\\ \end{array}\right. \end{aligned}$$
(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 \mathcal {V}_d\), \(p\in H^1(\Omega ^-)\) and \(\lambda \in \mathbb {R}\) such that:

$$\begin{aligned} \left\{ \begin{array}{l} a(u,u^{\!\star })+b(p,u^{\!\star })+c(\lambda ,u^{\!\star })= f(u^{\!\star })\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ b( p^{\!\star },u) = l_p( p^{\!\star }) \quad \forall p^{\!\star }\in H^1(\Omega ^-),\\ c( \lambda ^{\!\star },u) = 0 \quad \forall \lambda ^{\!\star }\in \mathbb {R}. \end{array}\right. \end{aligned}$$
(37)

with:

(38)

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

Problem 13

Find \(\mathring{u}\in \mathcal {V}_d\), \(\mathring{p}\in H^1(\Omega ^-)\) and \(\mathring{\lambda }\in \mathbb {R}\) such that:

$$\begin{aligned} \left\{ \begin{array}{l} a(\mathring{u},u^{\!\star })+b(\mathring{p},u^{\!\star })+c(\mathring{\lambda },u^{\!\star })= f^s(u^{\!\star })\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ b( p^{\!\star },\mathring{u}) = l^s_{p}( p^{\!\star }) \quad \forall p^{\!\star }\in H^1(\Omega ^-),\\ c( \lambda ^{\!\star },\mathring{u}) = l^s_{\lambda }( \lambda ^{\!\star }) \quad \forall \lambda ^{\!\star }\in \mathbb {R}. \end{array}\right. \end{aligned}$$
(39)

with:

(40)

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 enable by the knowledge of the position of \(r_c\) between two nodes. This approximation fields are illustrated in Figure 10.

Figure 9
figure 9

Axi-symmetric simplification of the membrane (right side) and illustration of the solution.

Once these two problems are solved we have access to \(\varrho\) and \(\mathring{\varrho }\) and we can solve the Problem 10 associated to our case.

$$\begin{aligned} \varrho +\mathring{\varrho }\tau =0 \quad \Leftrightarrow \quad \tau =-\frac{\varrho }{\mathring{\varrho }} \end{aligned}$$
(41)

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

$$\begin{aligned} r_c^{k+1}=r_c^{k}+\tau =r_c^{k}-\frac{\varrho }{\mathring{\varrho }} \end{aligned}$$
(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 \(\varrho -\varrho _a\). The weak discontinuity of the displacement is well represented by the X-FEM Heaviside enrichment combined with Lagrange multipliers.

Figure 10
figure 10

Approximation of the fields.

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^{\text {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:

$$\begin{aligned} \left\{ \begin{array}{l}\displaystyle a(u,u^{\!\star })+\sum _{k=1}^{N_c}p_ku^{\!\star }(\underline{x}_c^k)= f(u^{\!\star })\quad \forall u^{\!\star }\in \mathcal {V}_d , \\ \displaystyle \sum _{k=1}^{N_c} p^{\!\star }_k u(\underline{x}_c^k) = \displaystyle -\sum _{k=1}^{N_c} p^{\!\star }_k d(\underline{x}_c^k) \quad \forall p^{\!\star }\in \mathbb {R}^{N_c}. \end{array}\right. \end{aligned}$$
(43)

The initial \(r_c\) is chosen as \(r_{c_0} =0.3\text {m}\). 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^{\text {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.

Figure 11
figure 11

Evolution of a the error on \(r_c\), b the number of iterations, with the size of the mesh.

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:

(44)

In order to improve the efficiency of the algorithm a limitation is added to the basic algorithm described Figure 5. First, as the starting position of the interface can be far from the converged one, the evolution of \(\Gamma\) is limited such that it does not evolve further than the limit set by the value of \(\delta\) (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 \(\delta\) was taken on purpose. This explained the constant advance of the front for the first iterations.

Figure 12
figure 12

Evolution of the contact interface from initial to final configuration (b) with the displacement associated to these to extreme position (a, c).

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

Figure 13
figure 13

Evolution of the norm of \(\tau\) and percentage of the contribution of each mode.

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 14
figure 14

Error on the position of the contact interface compare to a reference solution for each problem.

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

Figure 15
figure 15

Evolution of a the error on \(\Gamma\), b the number of iterations, with the size of the mesh.

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

Table 1 Time comparison between ILS and active set method. First on fixed meshes then for a targeted error

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.

Figure 16
figure 16

Floor shape (2D and 3D representation) and associated mesh.

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 17
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).

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

References

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

    MATH  Google Scholar 

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

    Article  MathSciNet  Google Scholar 

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

    Article  MathSciNet  MATH  Google Scholar 

  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

    Article  MathSciNet  MATH  Google Scholar 

  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

    Book  Google Scholar 

  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

    Article  MathSciNet  Google Scholar 

  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

    Article  MATH  Google Scholar 

  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

    Article  MathSciNet  Google Scholar 

  10. Conry TF, Seireg A (1971) A mathematical programming method for design of elastic bodies in contact. J Appl Mech 38(2):387–392

    Article  MATH  Google Scholar 

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

    Google Scholar 

  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

    Article  Google Scholar 

  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

    Article  MATH  Google Scholar 

  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

    Article  Google Scholar 

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

    Book  Google Scholar 

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

    Article  MathSciNet  MATH  Google Scholar 

  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

    Article  MATH  Google Scholar 

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

    Article  MathSciNet  MATH  Google Scholar 

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

    Article  MathSciNet  Google Scholar 

  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

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    MathSciNet  Google Scholar 

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

    Google Scholar 

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

    Google Scholar 

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

    Article  MATH  Google Scholar 

  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

    Article  Google Scholar 

  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

    Article  Google Scholar 

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

    Article  MATH  Google Scholar 

  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

    Article  MathSciNet  MATH  Google Scholar 

Download references

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

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Nicolas Chevaugeon.

Appendix A: calculus of the sensibility

Appendix A: calculus of the sensibility

Lets first recall some basic formulas:

$$\begin{aligned} \mathring{\overline{\nabla u}}=\nabla \dot{u} +\nabla \nabla u \cdot w= \nabla (\dot{u} +\nabla u \cdot w)-\nabla u\cdot \nabla w=\nabla \mathring{u} -\nabla u\nabla w\end{aligned}$$
(45)
$$\begin{aligned} D\left( \int _{\Omega } f\,\text {d}\Omega \right) [\underline{w}]=\int _{\Omega } D(f)[\underline{w}]+f\text {div}\underline{w}\,\text {d}\Omega \end{aligned}$$
(46)
$$\begin{aligned} D\left( \int _{\Gamma } f\,\text {d}\Gamma \right) [\underline{w}]=\int _{\Gamma } D(f)[\underline{w}]+f(\text {div}\underline{w}-\underline{n}\cdot \nabla \underline{w}\cdot \underline{n})\,\text {d}\Gamma \end{aligned}$$
(47)

We have

$$\begin{aligned} D(a(u,u^{\!\star }))[\underline{w}]&=\int _{\Omega }T\left( \mathring{\overline{\nabla u}}\cdot \nabla u^{\!\star }+ \nabla u\cdot \mathring{\overline{\nabla u^{\!\star }}} +\nabla u\cdot \nabla u^{\!\star }\,\text {div}\underline{w}\right) \,\text {d}\Omega \nonumber \\&= \int _{\Omega }T\left( (\nabla \mathring{u}- \nabla u\cdot \nabla \underline{w})\cdot \nabla u^{\!\star }+ \nabla u\cdot (\nabla \mathring{u^{\!\star }}- \nabla \underline{w}\nabla u^{\!\star })+\nabla u\cdot \nabla u^{\!\star }\,\text {div}\underline{w}\right) \,\text {d}\Omega \end{aligned}$$
(48)
$$\begin{aligned} D(b(p,u^{\!\star }))[\underline{w}]&= \int _{\Omega ^-}\!\!\!l_{\text {m}}(\nabla \mathring{p}\nabla u^{\!\star }+ \nabla p \nabla \mathring{u^{\!\star }}-\nabla p(\nabla \underline{w}^T+\nabla \underline{w})\nabla u^{\!\star }+\nabla p \nabla u^{\!\star }\text {div}\underline{w})\,\text {d}\Omega \nonumber \\&\quad +\int _{\Omega ^-}\!\!\!\mathring{p}u^{\!\star }+p \mathring{u^{\!\star }} +p u^{\!\star }\text {div}\underline{w}\,\text {d}\Omega \end{aligned}$$
(49)
$$\begin{aligned} D(c(\lambda ,u^{\!\star }))[\underline{w}]= \int _{\Gamma } \mathring{\lambda }[\![u^{\!\star }]\!]+\lambda [\![\mathring{u^{\!\star }}]\!]+\lambda [\![u^{\!\star }]\!](\text {div}\underline{w}-\underline{n}\cdot \nabla \underline{w}\cdot \underline{n})\,\text {d}\Gamma \end{aligned}$$
(50)
$$\begin{aligned} D(f(u^{\!\star }))[\underline{w}]= \int _{\Omega }\mathring{f}u^{\!\star }+f\mathring{u^{\!\star }}+fu^{\!\star }\,\text {div}\underline{w}\,\text {d}\Omega \end{aligned}$$
(51)
$$\begin{aligned} D(l_p( p^{\!\star }))[\underline{w}]&= -\int _{\Omega ^-}\!\!\!l_{\text {m}}(\nabla \mathring{ p^{\!\star }}\nabla d+ \nabla p^{\!\star }\nabla \mathring{d} -\nabla p^{\!\star }(\nabla \underline{w}^T+\nabla \underline{w})\nabla d+\nabla p^{\!\star }\nabla d\,\text {div}\underline{w})\,\text {d}\Omega \nonumber \\&\quad -\int _{\Omega ^-}\!\!\!\mathring{ p^{\!\star }}d+ p^{\!\star }\mathring{d} + p^{\!\star }d\,\text {div}\underline{w}\,\text {d}\Omega \end{aligned}$$
(52)

The test fields \(u^{\!\star }\), \(p^{\!\star }\) and \(\lambda ^{\!\star }\) are chosen such that for all \(\underline{x}\) where the fields are defined, \(\mathring{u^{\!\star }}(\underline{x})=0\), \(\mathring{ p^{\!\star }}(\underline{x})=0\) and \(\mathring{ \lambda ^{\!\star }}(\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. \(\mathring{f}_d\)= 0) we find the formulation of Problem 8.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Graveleau, M., Chevaugeon, N. & Moës, N. The inequality level-set approach to handle contact: membrane case. Adv. Model. and Simul. in Eng. Sci. 2, 16 (2015). https://doi.org/10.1186/s40323-015-0034-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-015-0034-8

Keywords