# On the coupling of local multilevel mesh refinement and ZZ methods for unilateral frictional contact problems in elastostatics

Hao Liu, Isabelle Ramière, Frédéric Lebon

► To cite this version:

Hao Liu, Isabelle Ramière, Frédéric Lebon. On the coupling of local multilevel mesh refinement and ZZ methods for unilateral frictional contact problems in elastostatics. Computer Methods in Applied Mechanics and Engineering, 2017, 323, pp.1-26. 10.1016/j.cma.2017.04.011. hal-01694040

HAL Id: hal-01694040

https://hal.science/hal-01694040

Submitted on 29 Aug 2018

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
# On the coupling of local multilevel mesh refinement and ZZ methods for unilateral frictional contact problems in elastostatics

H. Liu<sup>a,b</sup>, I. Ramière<sup>a</sup>, F. Lebon<sup>b</sup>

<sup>a</sup> CEA, DEN, DEC, SESC, LSC, F-13108 Saint-Paul Lez Durance, France,
isabelle.ramiere@cea.fr

<sup>b</sup> Aix-Marseille University, CNRS, Centrale Marseille, LMA, Marseille, France,
lebon@lma.cnrs-mrs.fr

## Abstract

This paper addresses an efficient and robust automatic adaptive local multilevel mesh refinement strategy for unilateral frictional contact problems in elastostatics. The proposed strategy couples the Local Defect Correction multigrid method LDC (Hackbusch, 1984) with the ZZ (Zienkiewicz and Zhu, 1987) a posteriori error estimator. An extension of LDC method to frictional contact problems is introduced. An interesting feature of this extended LDC algorithm is that it still only lies on interpolations of displacement fields. Neither forces conservation nor exchange of contact status is required between the refinement levels. The ZZ a posteriori error estimator is exploited to automatically build the sub-grids of the LDC method. A criterion linked to the relative error in stress is used. The efficiency of the proposed strategy is analyzed on examples derived from nuclear engineering. Practical numerical choices are proposed and justified. The refinement process automatically varies and stops with respect to a given tolerance. Post-treatments show that the sub-grids focus around the contact areas and that the converged LDC solution always respects the prescribed tolerance.

**Keywords:** Automatic Adaptive Mesh Refinement, Multilevel Local Defect Correction method, Zienkiewicz and Zhu a posteriori error estimator, Frictional contact problems, Signorini-Coulomb, User tolerance.

## 1. Introduction

In most physical problems, local effects or singularities as stress concentration, crack, hot spot or vortex are encountered. The numerical resolution requires fine local mesh steps to precisely catch these effects. As in practice, it is not computationally possible to deal with fine uniform meshes, local mesh refinement methods have been extensively developed and used since the 80's, see for example [BS82, ZM83, BO84, Hac84, McC84, DDO85, Fis92, SH92, BT93, GM93, Bra94, KACM96, DR01, YR05, BR05, RAB07b, HR10, KL11, BRL14] and the references therein. At least two main families of refinement methods can be distinguished: the h-, p-, r-, s-refinement methods and the local multigrid methods. The first family consists in generating a locally refined global
mesh by adding elements (h-method [BS82, DDO85, SH92, BT93]), adding basis functions (p-method [BS82, ZM83, DR01]), moving nodes (r-method [GM93, HR10]), using a composite mesh (s-method [Fis92, YR05, BR05, JD13]) or combining the previous cited methods (hp-methods [BS87, DORH89, DNR07], hr-methods [SZ92, EOD93, ARF01], etc.). Note that the multiscale finite element approaches [HW97, DK08, EH09] can be seen as an extension of the p-refinement method, where the basis functions are defined from local fine problems. In counterpart, the second family [BO84, Hac84, McC84, Bra94, KACM96, RAB07b, KL11, BRL14] does not modify the initial global mesh but enriches the solution by a multigrid iterative process [Bra77] based on resolutions on local fine sub-meshes. These two families of adaptive refinement methods have shown their large efficiency. However in local multigrid methods the solver is a black-box making them so-called non-intrusive methods. Their implementation in any commercial computational software is then possible. Additionally, the resolutions on each grid are done separately which implies low numbers of degrees of freedom per resolution and often widely compensates for the iterative process cost. Finally, the possibility to use regular, often uniform, local meshes improves the conditioning of the problems. Hence, in this paper, we will focus on the extension of such approaches to frictional contact problems which has never been done to the best of our knowledge.

Frictional contact problems are widely encountered in mechanical engineering in particular in structure assembly or multibody interaction. These problems usually modeled by Signorini and Coulomb laws [Sig59, Cou85] are well known as non smooth mechanical ones [Mor63] and still nowadays non trivial to solve (solver robustness or efficiency, effective contact zones a priori unknown, etc.). Local multigrid methods will be used in order to build robust and efficient solving strategies. In our case, this mechanical problem will be discretized by the finite element method and it is chosen not to deal with composite problems, thus the Local Defect Correction (LDC) method introduced by Hackbusch [Hac84] seems the best candidate to solve the problem compared to Flux Interface Correction (FIC) [AL94] or Fast Adaptive Composite (FAC) [McC84] methods. The local Full MultiGrid (FMG) [CGLC05] method introduced recently returns to the LDC method where each new local level is added only after convergence of previous multigrid cycles. This strategy accelerates the speed of convergence of the global iterative process when the refinement is not localized (see [Bra94]) but seems less useful in the case of localized zones of interest [Liu16] and has not be retained in our purpose. In opposition to the literature where the potential contact zones are usually a priori refined [RCL88, MSTC08], the coupling of an adaptive mesh refinement method with a posteriori error estimator enables to find automatically the optimal mesh steps around the contact zones [WS98] and will improve the ratio accuracy over cost of the numerical method. In this paper, we propose an efficient and robust automatic adaptive multilevel mesh refinement strategy for unilateral frictional contact problems in elastostatics. The paper is decomposed as follow: in Section 2, the frictional contact problem in elastostatics as well as the classical formulations and associated numerical solvers are recalled. Section 3 is devoted to the mathematical formulation of the extension of LDC method to frictional contact problems and its coupling with a posteriori error estimator. The efficiency of the proposed automatic refinement strategy is deeply analyzed in Section 4 on examples derived from nuclear engineering. Practical numerical choices are also detailed in this section.
## 2. Frictional contact problem in elastostatics

### 2.1. Position of the problem

We consider two linearly elastic solids $\Omega_1$ and $\Omega_2$ in contact (see figure 1), and $\Omega = \Omega_1 \cup \Omega_2$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Model problem of contact between two deformable solids</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Decomposition on the contact zone</figcaption></figure>

We take $\partial\Omega_1$ and $\partial\Omega_2$ to denote the sufficiently smooth boundaries of these two solids, respectively and $\partial\Omega = \partial\Omega_1 \cup \partial\Omega_2$. The boundary $\partial\Omega$ is divided in three parts $\partial\Omega_D$, $\partial\Omega_N$, $\partial\Omega_C$ such that $\partial\Omega = \partial\Omega_D \cup \partial\Omega_N \cup \partial\Omega_C$ with $\partial\Omega_D \cap \partial\Omega_N = \emptyset$ and $\partial\Omega_D \cap \partial\Omega_C = \emptyset$:

* on $\partial\Omega_D$ displacement conditions are imposed (Dirichlet boundary conditions).

* on $\partial\Omega_N$ surface forces are applied (Neumann boundary conditions).

* $\partial\Omega_C$ is the possible and a priori unknown contact zone (contact boundary conditions). We have $\partial\Omega_C = \partial\Omega_C^1 \cup \partial\Omega_C^2$, with $\partial\Omega_C^i$, $i = 1, 2$ the potential contact boundary on $\Omega_i$, $i = 1, 2$.

Let $n^1$ (resp. $n^2$) denote the external unit normal to $\Omega_1$ (resp. $\Omega_2$) on $\partial\Omega_C$. We take $u^1$ (resp. $u^2$) to denote the displacement field in $\Omega_1$ (resp. $\Omega_2$). We take $\sigma^1$ (resp. $\sigma^2$) to denote the stress field in $\Omega_1$ (resp. $\Omega_2$). The displacement field $u^i$ and the unknown contact force $F^i$ on $\partial\Omega_C^i$ ($i=1,2$) are decomposed in normal and tangential parts:

$$ \begin{aligned} u^1 &= u_N^1 \cdot n^1 + u_T^1, \quad u^2 = u_N^2 \cdot n^2 + u_T^2 \\ F^1 &= \sigma^1 \cdot n^1 = F_N^1 \cdot n^1 + F_T^1, \quad F^2 = \sigma^2 \cdot n^2 = F_N^2 \cdot n^2 + F_T^2 \end{aligned} \quad (1) $$

Note that under the small perturbations hypothesis, $n^2 = -n^1$. The Signorini's law of unilateral contact [Sig59] is given on $\partial\Omega_C$ by

$$ \begin{aligned} F_N &= F_N^1 = F_N^2 \le 0 \quad (\text{compression condition, no adhesion}) \\ u_N &= u_N^1 + u_N^2 \le d \quad (\text{non penetration condition}) \\ (u_N - d)F_N &= 0 \quad (\text{complementarity condition}) \end{aligned} \quad (2) $$
where $d$ is a given function (initial gap). The quasi-static Coulomb's law of dry friction [Cou85, Amo99] is given on $\partial\Omega_C$ by

$$ F_T = F_T^1 = F_T^2, |F_T| \le -\mu F_N \quad \text{and} \quad \begin{cases} |F_T| < -\mu F_N \Rightarrow u_T = 0 \quad (\text{adherence}) \\ |F_T| = -\mu F_N \Rightarrow \exists \lambda > 0, u_T = -\lambda F_T \quad (\text{sliding}) \end{cases} \qquad (3) $$

where $\mu$ is the given friction coefficient and $|\cdot|$ the euclidian norm. Note that this formulation is valid for monotonous loadings. In conclusion, three status are possible on the contact zone: no contact, contact with adherence and contact with sliding. The strong formulation of the elastostatics frictional contact problem can be expressed as

$$ \left\{ \begin{array}{ll} -\operatorname{div} \sigma = f & \text{in } \Omega \quad (4a) \\ \sigma = C\varepsilon & \text{in } \Omega \quad (4b) \\ \varepsilon = \frac{1}{2} (\operatorname{grad} u + \operatorname{grad}^T u) & \text{in } \Omega \quad (4c) \\ u = u_0 & \text{on } \partial\Omega_D \quad (4d) \\ \sigma n = g & \text{on } \partial\Omega_N \quad (4e) \\ \sigma n = F_N n + F_T & \text{on } \partial\Omega_C \quad (4f) \\ u_N \le d, F_N \le 0, (u_N - d)F_N = 0 & \text{on } \partial\Omega_C \quad (4g) \\ |F_T| \le \mu |F_N| \text{ with } \begin{cases} |F_T| < \mu |F_N| \Rightarrow u_T = 0 \\ |F_T| = \mu |F_N| \Rightarrow \exists \lambda > 0, u_T = -\lambda F_T \end{cases} & \text{on } \partial\Omega_C \quad (4h) \end{array} \right. $$

where $C$ is the fourth order elasticity tensor (with usual properties of symmetry and positivity), $f$ are body forces, $g$ are surface forces and $u_0$ are given displacements. Formally, problem (4) can be written on the following form:

$$ \left\{ \begin{array}{l} \text{Find } u \in K_d \text{ such that} \\ \mathcal{L}(u) = \hat{f} + F(u) \end{array} \right. \qquad (5) $$

where $K_d = \{v \in (H^1(\Omega))^3, v_N \le d \text{ on } \partial\Omega_C\}$, $\mathcal{L}$ is the differential operator obtained by substituting (4b) and (4c) in Eqs. (4a) (4d) (4e) (4f), $\hat{f}$ is the right-hand side associated to Eqs. (4a) (4d) (4e), and $F(u)$ is the right-hand side of Eq. (4f) verifying conditions (4g) and (4h). For more details, the reader can referred for example to [AC91].

## 2.2. Weak formulation

The weak formulation of problem (4), a quasi-variational inequality, is given by [DL76],

$$ \left\{ \begin{array}{l} \text{Find } u \in K_d \text{ such that} \\ a(u, v-u) + j(u, v) - j(u, u) \ge l(v-u) \quad \forall v \in K_d \end{array} \right. \qquad (6) $$
where

$$
\begin{align}
a(u, v) &= \int_{\Omega} C \varepsilon(u) \cdot \varepsilon(v) d\Omega \nonumber \\
l(v) &= \int_{\Omega} f \cdot v d\Omega + \int_{\partial\Omega_N} g \cdot v dS \tag{7} \\
j(u, v) &= - \int_{\partial\Omega_C} \mu F_N(u) |v_T| dS \nonumber
\end{align}
$$

Note that the existence and uniqueness result in a general case is an open problem. It can be observed (see [AC91, Leb03, Wri06] for example) that the last equation (primal/displacement formulation) can also be written following the mixed formulation

$$
\left\{
\begin{array}{ll}
\text{Find } v = (u, \lambda) \in K_d \times M \text{ such that} & \\
\tilde{\mathcal{L}}(v) = \tilde{f} & (8)
\end{array}
\right.
$$

where $\lambda = -F$ and $M$ is taken to denote a subset of $(H^{-1/2}(\partial\Omega_C))^3$.

2.3. Numerical resolution

In this section, we briefly recall the main numerical methods to solve fric-
tional contact problems. Usually, relaxation procedures (see for example [RCL88])
are associated to the displacements formulation via a fixed point on the non dif-
ferentiable term $j(\cdot, \cdot)$ in Equation (6). These methods are well known to be
robust. Some authors, see [WB98, Leb95, CHP00, LRR07] as examples, have
proposed methods of $h$-refinement type to improve this kind of method. The
other main class of methods concerns the methods based on mixed formulations
[AC91, Wri06]. These methods are well known to be also efficient. Some com-
parisons between these two classes of methods can be found in [CDR98, Leb03]
and an improvement with $h$-techniques in [AL95]. Active set methods (see
[ABD16] and references therein) can be considered as methods of second kind.
These methods are based on a sequence of linear problem resolution with given
Dirichlet (contact displacement) or Neumann (contact force) boundary con-
ditions. Without any loss of generality, the latter will be used afterward as
frictional contact solver.

3. On the coupling of Local Defect Correction method and Zienkiewicz and Zhu a posteriori estimator for frictional contact problems

3.1. Local Defect Correction method for frictional contact problems

The Local Defect Correction (LDC) method has been introduced by Hack-
busch [Hac84]. This method is part of the multigrid process family. It consists
in relying several level of grids through prolongation and restriction operators.
As the Full MultiGrid process [Bra77], the initial mesh is a coarse triangulation
of the whole domain. The main characteristic of the LDC method is to generate
the additional fine meshes only in some local zones of interest (typically where
the error is large), see Figure 3 and to solve sequentially local problems on these
meshes. The multilevel LDC algorithm is usually represented by ∧-cycles, see
Figure 4.
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of sub-grids</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Local multigrid process representation for 4 levels of grids</figcaption></figure>

This adaptive refinement method has the great advantages of being non-intrusive, very flexible (the solver and/or the refinement ratio and/or the mathematical model can change between levels of refinement), and powerful (low degrees of freedom by level, possibility to use uniform meshes,...). It has been developed and efficiently performed in the field of fluid mechanics [BG03, AMB03, MAM06, KCMM09, RAB07a]. Some recent works extended successfully this method to linear elasticity [BRL14] and nonlinear material mechanics [BRL15] which open the way of the current work.

In this section, an extension of the standard algorithm to frictional contact problems is introduced. Let ($\mathcal{P}$) denote the partial differential problem under consideration defined on the domain $\Omega$ with boundary $\partial\Omega$:

$$ (\mathcal{P}) : \mathcal{L}(u) = f^{\text{ext}} \qquad (9) $$
with $\mathcal{L}$ a differential operator in space, $f^{ext}$ the right hand side and $u$ the displacement solution.

In frictional contact problems, we chose $f^{ext}$ as the sum between the given external forces and the unknown contact forces, see problem (5).

We then define embedded local domains $\Omega_l$, $0 \le l \le l^*$, $\Omega_l \subset \Omega_{l-1}$ with $\Omega_0 = \Omega$. On each local domain $\Omega_l$ with boundary $\partial\Omega_l$, we have to solve the following problem:

$$ (\mathcal{P}|_{\Omega_l}) : \mathcal{L}|_{\Omega_l}(u) = f^{ext}|_{\Omega_l} \qquad (10) $$

The boundary conditions on $\partial\Omega_l$ will be precised by the LDC algorithm detailed hereafter.

Each domain $\Omega_l$ is discretized by a triangulation $\mathcal{T}_l$ of mesh step $h_l$ with $h_l < h_{l-1}$. The boundary of the underlying grid $G_l$ is denoted $\Gamma_l$. Then the approximated problem $(\mathcal{P}_l)$ on $G_l$ holds:

$$ (\mathcal{P}_l) : \mathcal{L}_l(u_l) = f_l \qquad (11) $$

where $\mathcal{L}_l$ is the discrete form of the operator $\mathcal{L}|_{\Omega_l}$ and $f_l$ the discretization of the right-hand side $f^{ext}|_{\Omega_l}$.

The LDC method consists in recursively solving problem $(\mathcal{P}_l)$ with modified right-hand sides. Hence at the $k^{th}$ iteration, $(\mathcal{P}_l)$ writes:

$$ (\mathcal{P}_l^k) : \mathcal{L}_l(u_l^k) = f_l^k \qquad (12) $$

Generally speaking, the prolongation operator consists in defining Dirichlet boundary conditions (BCs) on the nodes of the fictive surface $\partial\Omega_l^f = \partial\Omega_l \setminus \partial\Omega$ whereas the restriction operator works on the correction of the right-hand side. In practice these operators are based on interpolation techniques. It is worth noting that contrary to the classical multigrid method [Bra77], the restriction operator has not to be the transpose of the prolongation operator, see [Hac84]. Contact problems require a change in the definition of fictive surface usually defined as $\partial\Omega_l^f$. Indeed, the gap can be evaluated only between paired boundaries which induces that the non penetration condition can also only be written between paired boundaries. We then propose the following definition of the potential contact boundary $\partial\Omega_{C,l}$ at level $l$:

$$ \left\{ \begin{array}{l} \partial\Omega_{C,l}^{i=1,2} \text{ is the restriction of } \partial\Omega_C^i \text{ to } \partial\Omega_l \text{ such that each point of } \partial\Omega_{C,l}^1 \text{ can be} \\ \text{paired to a point of } \partial\Omega_{C,l}^2 \\ \partial\Omega_{C,l} = \partial\Omega_{C,l}^1 \cup \partial\Omega_{C,l}^2 \end{array} \right. $$

and to extend the definition of the fictive boundary to avoid this issue:

$$ \begin{aligned} \partial\Omega_{C,l}^f &= (\partial\Omega_C \cap \partial\Omega_l) \setminus \partial\Omega_{C,l} \\ \Gamma_l^f &= (\partial\Omega_l^f \cup \partial\Omega_{C,l}^f) \cap \Gamma_l \end{aligned} \qquad (13) $$

The localization of the fine grids requires to define a set of restriction nodes $A_l = \{x \in (G_{l+1} \setminus \Gamma_{l+1}^f) \cap (G_l \setminus \Gamma_l^f)\}$ and a set of correction nodes $\tilde{A}_l = \{x \in A_l; \mathcal{L}_l(u_l(x)) \text{ involves only } y \in A_l\}$, see Figure 5. The subset $\tilde{A}_l$ enables to correct the right-hand side without any perturbation from the fictive BCs. It is worth noting that the subset $A_l$ proposed here (and consequently $\tilde{A}_l$) is different from the standard Hackbusch's LDC subset [Hac84]: the nodes on $\Gamma_l \cap \partial\Omega$ are
included into $A_l$. The main idea here is to benefit at level $l$ from the better discretization of the problem's boundary conditions on the finer grids. This modification is essential to accurately deal with contact problems.

<figure><img src="image_6.png" /><figcaption>Figure 5: Example of restriction and correction sets</figcaption></figure>

Once the subsets $A_l$ and $\hat{A}_l$ are defined as previously, the proposed LDC algorithm for frictional contact problem only differs from the usual one [Hac84, ABS05, BRL14] by the definition of the local residual on $\hat{A}_l$. With the following notation (see problem (5))

$$f^{ext} = \hat{f} + F(u) \qquad (14)$$

the right-hand side of problem $(\mathcal{P}_l^k)$ writes

$$f_l^k = \hat{f}_l^k + F_l(u_l^k) \qquad (15)$$

where $\hat{f}_l^k$ can be viewed as an update of the given forces as in the standard LDC algorithm.

The restriction step on $G_l$ consists in

* only restrict the displacement field with a given restriction operator $R_{l+1}^l$

$$\tilde{u}_l^k(x) = R_{l+1}^l u_{l+1}^k(x) \quad \forall x \in A_l$$

* compute the local residual

$$r_l^k(u)(x) = (\mathcal{L}_l(\tilde{u}_l^k) - \hat{f}_l^0 - (F(u))^k_l)(x) \quad \forall x \in \hat{A}_l$$

where $\hat{f}_l^0$ are the initial problem given forces and $(F(u))^k_l$ denotes the effective contact forces on $G_l$, derived from the prolongation step. Note that we have chosen to evaluate the contact forces as a given right-hand side at level $l$ (as $\hat{f}_l^0$) in order to be independent on the formulation and on the resolution of the frictional contact problem (see Section 2.2), i.e. to only restrict the displacement field. Otherwise, in case of a mixed formulation it would be necessary to restrict the contact forces, which may be avoided in the current approach (in particular as the force conservation would have to be ensured).

* update the right-hand side

$$\hat{f}_l^k(x) = \hat{f}_l^0(x) + \chi_{\hat{A}_l}(x) r_l^k(u)(x)$$
hence

$$f_l^k(x) = \hat{f}_l^0(x) + \chi_{\hat{A}_l}(x)r_l^k(u)(x) + F_l(u_l^k)(x)$$

where $\chi_{\hat{A}_l}$ denotes the characteristic function of $\hat{A}_l$: $\chi_{\hat{A}_l}(x) = \begin{cases} 1 & \text{if } x \in \hat{A}_l \\ 0 & \text{elsewhere} \end{cases}$

The whole LDC algorithm for frictional contact problems is detailed in Algorithm 1.

It is worth noting that contrary to classical multigrid methods for contact problems [Leb95], the proposed algorithm only exchanges the displacement field between levels, which is easier to implement and seems more robust.

Algorithm 1 is written in a sequential manner between the levels of grids similarly to a Gauss-Seidel process. This algorithm can easily be modified to exploit parallel computation using a Jacobi-like process from the previous iteration grids solutions. However, the speed of convergence will be lower [LQ92].
<caption>**Algorithm 1:** Multilevel LDC algorithm for frictional contact problems</caption>

**Input:** $k_{max}$: maximal number of iterations, $l^*$: number of levels of refinement, $G_l$: grid at level $l$, tol: given tolerance for convergence
**Output:** $u_l$ for $l = 0, ..., l^*$
**Initialization:** computation of the solution $u_0^0$ of initial problem ($\mathcal{P}_0^0$)
//Iteration: actualization of $u_l^k$//
**for** $k=1$ **to** $k_{max}$ **do**
  $u_0^k = u_0^{k-1}$
  //Prolongation step//
  **for** $l=1$ **to** $l^*$ **do**
    **if** $k=1$ **then**
      Definition of natural BCs on $\Gamma_l \setminus \Gamma_l^f$
      Definition of right-hand side $\hat{f}_l^0$ on $G_l$ from $\hat{f}$
    **end**
    Definition of Dirichlet BCs on $\Gamma_l^f$ thanks to prolongation operator $P_{l-1}^l$
    $u_{\Gamma_l^f}^k = (P_{l-1}^l u_{l-1}^k)_{|\Gamma_l^f}$
    Computation of $u_l^k$ by solving ($\mathcal{P}_l^k$) with $f_l^k = f_l^{k-1}$
  **end**
  //Restriction step//
  **for** $l=l^*-1$ **to** $0$ **do**
    Restriction of $u_{l+1}^k$ on $A_l$ thanks to restriction operator $R_{l+1}^l$
    $\tilde{u}_l^k(x) = R_{l+1}^l u_{l+1}^k(x) \quad \forall x \in A_l$
    Computation of local residual on $\tilde{A}_l$
    $r_l^k(u)(x) = (\mathcal{L}_l(\tilde{u}_l^k(x)) - \hat{f}_l^0 - (F(u))_l^k)(x) \quad \forall x \in \tilde{A}_l$
    Right-hand side update
    $\hat{f}_l^k(x) = \hat{f}_l^0(x) + \chi_{\tilde{A}_l}(x) r_l^k(u)(x)$
    $f_l^k(x) = \hat{f}_l^k(x) + F_l(u_l^k)(x)$
    Computation of $u_l^k$ solution of ($\mathcal{P}_l^k$)
  **end**
  **if** $\frac{\|u_0^k - u_0^{k-1}\|}{\|u_0^{k-1}\|} < \text{tol}$ **then**
    end the algorithm
  **end**
**end**
### 3.2. Coupling LDC and a local error estimator

The LDC method improves the accuracy of the coarsest grid solution if and only if the subgrids are localized around the zones where the solution errors are maximum. Thus, LDC efficiency strictly depends on well detected refinement zones.

The LDC method has then to use a potent error estimator. The genericity of the LDC algorithm enables us to easily couple this method with any error estimator that gives a local indication of the error on each element of the mesh. We denote by $e_K$ the local error indicator an element $K$ of the triangulation $\mathcal{T}_l$. The following criterion will then be used to detect the maximal error zones

$$ \|e_K\| > \alpha \qquad (16) $$

where $\|\cdot\|$ is a chosen norm and $\alpha$ is the user defined tolerance.

At the first prolongation step, the following embedded domains $\Omega_l$, $1 \le l \le l^*$ are thus constructed:

$$ \Omega_l = \{ \cup \bar{K}; K \in \mathcal{T}_{l-1} \text{ and } \|e_K\| > \alpha \} \qquad (17) $$

Criterion (16) may be self-sufficient to stop the level generation except in the presence of singularities. To deal with this case, an additional stopping criterion proposed in [LRL17] is used here to face this issue. This extra criterion is strictly related to the LDC algorithm. In order to remain independent of the considered problem, this criterion is no more based on a minimal number of elements [BRL14] but on geometrical considerations. It relies on the fact that if $\Omega_{l^*}$ is too small compared to $\Omega_0$, the correction on this level has only a really low influence on the global solution. At least two choices are possible to stop the refinement process at level $l-1$, see Figure 6:

* The surface $S_l$ of $\Omega_l$ is less than the surface $S_e$ of one element of $\mathcal{T}_0$

* The surface $S_l$ of $\Omega_l$ is less than a ratio of initial surface $S_0$ of $\Omega$

<figure><img src="image_7.png" /><figcaption>Figure 6: Illustration of the geometrical stopping criteria.</figcaption></figure>

The second choice is preferred as the first one depends on the initial triangulation $\mathcal{T}_0$ and thus will not stop at the same local mesh step, see [LRL17]. The multilevel LDC algorithm for frictional contact problems coupled with a local error estimator and a geometrical stopping criterion is reported in Algorithm 2.
<caption>**Algorithm 2:** Multilevel LDC algorithm for frictional contact problems coupled with a local error estimator</caption>

**Input:** $k_{max}$: maximal number of iterations, $l_{max}$: maximal number of refinement levels, $\alpha$: maximal user-desired error, $\beta$: surface ratio, $r_l$: refinement ratio at level $l$, tol: given tolerance for convergence

**Output:** $l^*$: number of levels of refinement, $u_l$ for $l = 0, ..., l^*$

**Initialization:** computation of the solution $u_0^0$ of initial problem ($\mathcal{P}_0^0$),

$$l^* = l_{max}$$

//Iteration : actualization of $u_l^k$//

for $k=1$ to $k_{max}$ do

$$u_0^k = u_0^{k-1}$$

//Prolongation step//

for $l=1$ to $l^*$ do

if $k=1$ then

$$\Omega_l = \{\cup \bar{K}; K \in \mathcal{T}_{l-1} \text{ and } \|e_K\| > \alpha\}$$

if $S_l = \text{meas}(\Omega_l) > \beta S_0$ then

Discretization of $\Omega_l$ with grid $G_l$ of mesh step $h_l = h_{l-1}/r_l$

else $l^* = l - 1$

Definition of natural BCs on $\Gamma_l \setminus \Gamma_l^f$

Definition of right-hand side $\hat{f}_l^0$ on $G_l$ from $\hat{f}$

end

Definition of Dirichlet BCs on $\Gamma_l^f$ thanks to prolongation operator

$P_{l-1}^l$

$$u_{\Gamma_l^f}^k = (P_{l-1}^l u_{l-1}^k)_{|\Gamma_l^f}$$

Computation of $u_l^k$ by solving ($\mathcal{P}_l^k$) with $f_l^k = f_l^{k-1}$

end

//Restriction step//

for $l = l^* - 1$ to $0$ do

Restriction of $u_{l+1}^k$ on $A_l$ thanks to restriction operator $R_{l+1}^l$

$$\tilde{u}_l^k(x) = R_{l+1}^l u_{l+1}^k(x) \quad \forall x \in A_l$$

Computation of local residual on $\tilde{A}_l$

$$r_l^k(u)(x) = (\mathcal{L}_l(\tilde{u}_l^k(x)) - \hat{f}_l^0 - (F(u))_l^k)(x) \quad \forall x \in \tilde{A}_l$$

Right-hand side update

$$f_l^k(x) = \hat{f}_l^0(x) + \chi_{\tilde{A}_l}(x) r_l^k(u)(x) + F_l(u_l^k)(x)$$

Computation of $u_l^k$ solution of ($\mathcal{P}_l^k$)

end

if $\frac{\|u_0^k - u_0^{k-1}\|}{\|u_0^{k-1}\|} < \text{tol}$ then
end the algorithm

end

end
As Algorithm 1, Algorithm 2 could be parallelized. However the first pro- longation step should remain sequential due to the sub-grids detection by an a posteriori error estimator.

### 3.3. Zienkiewicz and Zhu a posteriori error estimators

Since last decades several techniques have been proposed to estimate the discretization errors i.e. the difference, in a chosen norm, between the exact solution of the problem ($P$) and its discrete approximation given by the solution of problem ($P_l$). Contrary to a priori error estimators which provide only information on the asymptotic behavior of the discretization errors (typically orders of convergence), a posteriori error estimators employ the approximated solution itself to derive local estimates of the computed solution error. This latter is representative of the discretization error while the resolution error can be neglected.

Nowadays, various a posteriori error estimates have been developed for a large class of problems: Stokes [OWA94], Maxwell [RZ05], Cracks [PDH06, LFPD15], Contact [WS98], Coulomb friction [HL09, KL11], Plasticity [GLP00], etc. These estimates are mainly based on one of the three family of error estimates: explicit and implicit element or subdomain residual estimates [BR78b, BR78a, BW85], recovery-based error estimators [SB86, ZZ87, ZZ92a, ZZ92b], estimators based on constitutive equations [LL83]. As presented in [GB05], recovery-based error estimators lie on a comparison between the gradient of the approximated solution and a smoothed version of the gradient. This methodology relies on the idea that usually, using standard discretization methods (e.g. finite element method), the gradient is discontinuous along the boundary of elements. The family of Zienkiewicz and Zhu (ZZ) a posteriori error estimators [ZZ87, ZZ92a, ZZ92b] can be seen as particular recovery-based error estimators, where the smoothed version of the gradient $\nabla u_l^{smooth}$ is decomposed on the base $N$ related the interpolation functions of the discretization method:

$$ \nabla u_l^{smooth} = \sum_i N_i \nabla u_l^{smooth,i} \qquad (18) $$

where $\nabla u_l^{smooth,i}$ are the nodal values of $\nabla u_l^{smooth}$. The coefficients $\nabla u_l^{smooth,i}$ are usually determined either by a projection [ZZ87] or by a non nodal polynomial approximation [ZZ92a, ZZ92b]. Then, an error indicator on each element $K$ is given by the difference, in a chosen norm, between the gradient of the approximation and the smoothed gradient on this element.

**Remark.** The classical ZZ error estimator is initially written in terms of Cauchy stress tensor $\sigma = C\nabla_s u$, where $C$ is the fourth order elasticity tensor and $\varepsilon(u) = \nabla_s u$, the strain tensor, is the symmetric part of the gradient of displacement $u$. This estimator can be seen as a weighted version of the estimator on the gradient.

In most industrial computational softwares (Abaqus, Aster, etc.), ZZ error estimator had been chosen for its implementation simplicity and its good ratio precision over cost. In section 4, it will be coupled with the frictional contact LDC method. To avoid modifying the ZZ estimator to take into account contact
contributions as in [WS98], the errors will be estimated in each body separately as proposed by [NL92] for multimaterial problems. Then

$$ \Omega_l = \{ \cup \bar{K}; K \in \mathcal{T}_{l-1}; \{K \subset \Omega_1 \text{ and } \|e_K^1\| > \alpha\} \text{ or } \{K \subset \Omega_2 \text{ and } \|e_K^2\| > \alpha\} \} \quad (19) $$

where $e^1$ (respectively $e^2$) denotes the error estimator performed on $\Omega_1$ (respectively $\Omega_2$).

# 4. Numerical experiments

## 4.1. Pellet-Cladding interaction

The so-called Pellet-Cladding Interaction (PCI) phenomenon describes a mechanical interaction that occurs during the irradiation in Pressurized Water Reactor (PWR). As represented in Figure 7, a PWR fuel rod is essentially constituted by a cylindrical cladding tube containing a piling up of around 300 cylindrical pellets of approximately 8 mm in diameter by 13 mm high. The fuel pellet is made of uranium dioxide ($UO_2$) or mixed oxide of uranium and plutonium (MOX) and is a fragile ceramic with a failure stress value included between 100 MPa and 150 MPa. An initial gap about 85 $\mu m$ takes place between the fuel pellets and the cladding in order to facilitate the loading of the pellets.

<figure><img src="image_8.png" /><figcaption>Figure 7: Example of a PWR fuel rod, AFA3G® 900 MWe, AREVA</figcaption></figure>

Due to the low failure stress value, the fuel pellets crack at the early stages of irradiation. The fragmented pellets deform in a hourglass shape due to high temperature gradients. On the other hand, the accumulation of gaseous fission products induces fuel pellets swelling while the water's pressure around the fuel rods leads to the cladding's creep. The initial gap between the pellets and the cladding decreases and discontinuous contacts due pellets hourglass shape and cracks appear: it is the Pellet-Cladding interaction, see Figure 8.

The PCI phenomenon induce localized high stress concentrations on the cladding in front of both inter-pellet planes and cracks. These stress levels might lead to the cladding failure under certain circumstances, which want to be avoided as the cladding is the first confinement barrier of the fission products. Studies (experiments and simulations) to further understand and accurately simulate the failure technology of the cladding are still worldwide carried out.

## 4.2. Test cases

We propose here to apply the previously introduced adaptive refinement method, based on an automatic local multigrid process combining an extension
<figure><img src="image_9.png" /><figcaption>Figure 8: Pellet-Cladding interaction</figcaption></figure>

of the LDC method with the ZZ a posteriori error estimator, in order to accu- rately evaluate the stress field around the PCI contact zone.

In order to validate the proposed approach, we focus on two simplified PCI models representing respectively the two main mechanisms inducing contact: the hourglass shape deformation and the cracking of the fuel pellet. The mate- rials are assumed to be elastic. The cladding has a Young’s modulus of 190 GPa and a Poisson’s ratio of 0.3 while the fuel pellet has a Young’s modulus of 78 GPa and a Poisson’s ratio of 0.34.

* Hourglass shape deformation. This bidimensional axisymmetric test case, called the 2D(r,z) test case thereafter, enables to study the effect on the cladding of the hourglass shape deformation of the pellets. Due to sym- metry reasons, only a quarter of a pellet and the corresponding part of the cladding are represented, see Figure 9. In order to simply modelize the dis- continuous pellet-cladding contact, a discontinuous pressure is applied at the outer radius of the cladding, see Figure 9. We impose $P_1 = 0$ MPa and $P_2 = 200$ MPa. The pressure in the gap is also imposed at $P_0 = 0$ MPa. The geometry and boundary conditions of the problem are given in Fig- ure 9. An initial gap of $d = 2\mu m$ is set to simulate the first pellet-cladding contact.

* Pellet cracking. This bidimensional test case takes place in the (r, θ) plane. The plane strain conditions are applied. Due to symmetry conditions, only an eighth of the pellet and the corresponding opposite part of the cladding are represented, see Figure 10. The geometry and boundary conditions of this so-called 2D (r,θ) test case are also given in Figure 10. Here again, an initial gap of $d = 2\mu m$ is set.

In this test case, the contact between the pellet and the cladding is induced by a differential thermal strain in the pellet. A radial thermal field is imposed in the pellet while a constant thermal field is imposed in the cladding, see Figure 11:

$$ \begin{aligned} \Delta T &= 7.0 \times 10^7 \times (R_{\text{pellet}}^2 - r^2) + 561 - 293 \text{ K} & \text{in the pellet} \\ \Delta T &= 561 - 293 \text{ K} & \text{in the cladding} \end{aligned} $$
<figure><img src="image_10.png" /><figcaption>Figure 9: Geometry and boundary conditions for the 2D (r,z) test case.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Geometry and boundary conditions for the 2D (r,θ) test case.</figcaption></figure>

where $R_{pellet} = 4.1$ mm is the external radius of the pellet, see Figure 9.

<figure><img src="image_12.png" /><figcaption>Figure 11: Thermal field imposed in the 2D(r,θ) test case</figcaption></figure>

The dilatation coefficient is $\alpha_T = 1.0 \times 10^{-5} K^{-1}$ in the pellet and $\alpha_T = 6.0 \times 10^{-6} K^{-1}$ in the cladding. The corresponding forces are applied as source term $f$.

### 4.3. Numerical considerations

The LDC-ZZ strategy for frictional contact problems is performed on the two test cases detailed previously thanks to the sequential Algorithm 2 implemented in the Cast3M [CAS] software. A constant refinement ratio $r_l = r = 2$ is
applied between each level of grid. All the meshes are discretized using standard Lagrange $Q_1$ finite elements. The sub-grids are simply obtained by regularly dividing the next coarser detected elements in order to obtain hierarchical meshes, see Figure 12 for example. Indeed in this case, the prolongation and restriction operators are easier to define. Let us notice that the genericity of the LDC method enables the use of non-hierarchical meshes, see [NM02] for example.

<figure><img src="image_13.png" /><figcaption>Figure 12: Example of hierarchical element refinement with $r = 2$</figcaption></figure>

As we want to use structured meshes, as regular as possible to benefit from good solver performances, we chose them non-boundary fitted to the pressure discontinuity. The expected order of convergence of the method with respect to the mesh step is then linear for the $L^2$-norm, see [Ram08, BRL14]. Then, a bilinear interpolation (which reduces to a linear interpolation in 2D) based on the discretization nodes is used as prolongation operator while the canonical restriction is applied as restriction operator. The canonical restriction is possible if and only if all the coarse nodes of $A_l$ are also fine nodes, which is the case by construction for hierarchical meshes, see Figure 13.

<figure><img src="image_14.png" /><figcaption>Figure 13: Canonical restriction for hierarchical meshes</figcaption></figure>

Concerning the ZZ error estimator, the classical projection method [ZZ87] well-known as “ZZ1” is used. The following relative energy error semi-norm (which is equivalent to a norm for frictional contact problem (4) under consideration) is considered to detect the elements to be refined:

$$ \|e_K\| = \left( \frac{\int_K (\sigma_l^{\text{smooth}} - \sigma_l) : (C^{-1}(\sigma_l^{\text{smooth}} - \sigma_l)) \, dK}{\int_K \sigma_l : (C^{-1}\sigma_l) \, dK + \int_K (\sigma_l^{\text{smooth}} - \sigma_l) : (C^{-1}(\sigma_l^{\text{smooth}} - \sigma_l)) \, dK} \right)^{1/2} \quad (20) $$

This norm has been introduced in [ZZ87] in order to guarantee a non-zero denominator. The detected zone to be refined is first structured in order to keep
structured subgrids and to avoid artifacts in ZZ1 error estimations (e.g. entrant corners). They are then enlarged to guarantee the solution at all the detected nodes to be corrected by the restriction step (due to the definition of the sets $A_l$ and $\tilde{A}_l$), see Figure 14 for an example. Moreover for contact problems, the re-

<figure><img src="image_15.png" /><figcaption>Figure 14: Determination of the refinement zone. Example for $\alpha = 1\%$.</figcaption></figure>

finement zone is finally again enlarged in order to ensure potential contact zones to face each other at the next sub-level, see Figure 15 for example. Hence, $\forall l \partial\Omega_{C,l} = \partial\Omega_C \cap \partial\Omega_{l,0}$ and $\partial\Omega_{C,l}^f = \emptyset$.

<figure><img src="image_16.png" /><figcaption>Figure 15: Determination of the refinement zone for contact problems. Example for $\alpha = 5\%$.</figcaption></figure>

As recommended in [LRL17], the stopping criterion for subgrids generation is
set to $\beta = 0.5\%$.

Without loss of generality, an active set method is used at each level to solve the obtained frictional contact problems.

The LDC algorithm converges rapidly. For example, only 3 $\Lambda$-cycles are usually necessary to reach the chosen convergence tolerance $tol = 10^{-4}$. In order to assess the efficiency of the LDC-ZZ refinement process, the LDC solution after convergence will be compared to a reference solution obtained with a uniform very fine mesh (more than 1 billion of nodes). As the stress field may be singular around the contact zone, the maximal local error in this zone is not a suitable post-treatment as for the fracture mechanics field [Tay07]. We will then focus on the relative energy error norm $e_h$ evaluated on the whole LDC composite grid

$$e_h^2 = \frac{\sum_{l=1}^{l^*+1} \sum_{K \subset T_{l-1} \setminus T_l K} \int (\sigma_h - \sigma_{ref}) : (\varepsilon_h - \varepsilon_{ref}) dK}{\sum_{l=1}^{l^*+1} \sum_{K \subset T_{l-1} \setminus T_l K} \int \sigma_{ref} : \varepsilon_{ref} dK} \quad (21)$$

with $T_{l^*+1} = \emptyset$.

We will see in the following examples that this global relative norm may lead to a too important smoothing of the local errors (especially with imposed strain problems) and that the best post-treatment seems to be an evaluation of the relative error on a local zone located around the potential contact zone.

## 4.4. Numerical results

### 4.4.1. 2D(r,z) frictionless test case

A frictionless 2D(r,z) test (see Fig. 9) is first considered in order to validate the proposed methodology.

The robustness of the adaptive refinement method dedicated to contact problems will be appreciated thanks to

* various initial meshes

* various user defined tolerances $\alpha$

* relative energy norm errors with respect to a reference solution

* contact forces

An example of automatically detected levels of mesh for $h_0 = 200 \ \mu m$ and $\alpha = 5\%$ is given in figure 16. Let us note that in all following figures, the gap had been scaled by a factor 50 for a better visualization. We can see that the sublevels are more and more localized and than the submesh generation automatically stops. Looking at the deformed and the contact forces on the coarsest mesh, see Figs. 17 and 18, we can conclude that the refinement automatically localizes around the contact zone, especially around the zone of contact status change (gap, contact).

Table 1 summarizes the different numerical experiments done on this test case.

The first conclusion to be drawn from Table 1 is that the sub-levels generation always automatically stops. The number of nodes adding by level may be
<figure><img src="image_17.png" /><figcaption>Figure 16: Frictionless 2D(r,z) test case: detected sub-meshes for $h_0 = 200 \ \mu m$ and $\alpha = 5\%$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Frictionless 2D(r,z) test case: deformed for $h_0 = 200 \ \mu m$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: Frictionless 2D(r,z) test case: contact forces for $h_0 = 200 \ \mu m$</figcaption></figure>

<table><caption>Table 1: Frictionless 2D(r,z) test case: various coarse mesh steps $h_0$ and user tolerances $\alpha$.</caption><thead><tr><th rowspan="2">h<sub>0</sub></th><th rowspan="2">Initial</th><th colspan="3">user tolerance α</th></tr><tr><th>5%</th><th>2%</th><th>1%</th></tr></thead><tbody><tr><td rowspan="5">200 µm</td><td>l<sup>*</sup></td><td>0</td><td>3</td><td>4</td><td>5</td></tr><tr><td>e<sub>h</sub></td><td>19.7%</td><td>3.49%</td><td>1.27%</td><td>0.86%</td></tr><tr><td>N<sub>l</sub></td><td>1836</td><td>1836/2484/5290/11764</td><td>1836/5670/13152/22654/52500</td><td>1836/7070/21942/50050/87542/201008</td></tr><tr><td>NN</td><td>1836</td><td>21374</td><td>95812</td><td>369448</td></tr><tr><td>CPU</td><td>0.3s</td><td>19s</td><td>131s</td><td>773s</td></tr><tr><td rowspan="5">100 µm</td><td>l<sup>*</sup></td><td>0</td><td>2</td><td>3</td><td>4</td></tr><tr><td>e<sub>h</sub></td><td>7.73%</td><td>3.42%</td><td>1.27%</td><td>0.86%</td></tr><tr><td>N<sub>i</sub></td><td>7070</td><td>7070/5290/11764</td><td>7070/13066/22842/52500</td><td>7070/21942/50050/87542/201008</td></tr><tr><td>NN</td><td>7070</td><td>24124</td><td>95478</td><td>367612</td></tr><tr><td>CPU</td><td>1.2s</td><td>17s</td><td>133s</td><td>750s</td></tr><tr><td rowspan="5">50 µm</td><td>l<sup>*</sup></td><td>0</td><td>1</td><td>2</td><td>3</td></tr><tr><td>e<sub>h</sub></td><td>6.61%</td><td>3.18%</td><td>1.25%</td><td>0.86%</td></tr><tr><td>N<sub>l</sub></td><td>27738</td><td>27738/11764</td><td>27738/22842/52500</td><td>27738/50050/87542/201008</td></tr><tr><td>NN</td><td>27738</td><td>39502</td><td>103080</td><td>366338</td></tr><tr><td>CPU</td><td>8s</td><td>53s</td><td>173s</td><td>1044s</td></tr><tr><td rowspan="5">25 µm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>1</td><td>2</td></tr><tr><td>e<sub>h</sub></td><td>2.99%</td><td>2.99%</td><td>1.22%</td><td>0.85%</td></tr><tr><td>N<sub>i</sub></td><td>109874</td><td>109874</td><td>109874/51474</td><td>109874/87542/201008</td></tr><tr><td>NN</td><td>109874</td><td>109874</td><td>161348</td><td>398424</td></tr><tr><td>CPU</td><td>68s</td><td>68s</td><td>512s</td><td>1558s</td></tr></tbody></table>

$N_l$: number of nodes for $l = 0, ..., l^*$; $NN = \sum_{l=0}^{l^*} N_l$: total number of nodes.
really lower than the equivalent uniform local fine mesh (see 11764 compared to 109874 in the column 5 % for example). These sub-levels are moreover localized around the same zones whatever the initial mesh step: for a given $\alpha$ similar local meshes, in particular in the contact zones, are obtained for equivalent local mesh steps, see the number of sub-levels $l^*$ and nodes per level $N_l$ in each column of Table 1 as well as their localization in Figs. 19 and 20.

<figure><img src="image_20.png" /><figcaption>Figure 19: Frictionless 2D(r,z) test case: grids superposition for $h_0 = 200 \ \mu m$ and $\alpha = 5 \%$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20: Frictionless 2D(r,z) test case: grids superposition for $h_0 = 100 \ \mu m$ and $\alpha = 5 \%$</figcaption></figure>

The second conclusion is that the relative error rapidly decreases by adding local sub-levels. The LDC algorithm is hence really performing. Moreover the obtained errors are always lower than the prescribed user tolerance. Each sub-level seems necessary to respect the tolerance (see rows of Table 1). The coupling LDC-ZZ is then really efficient and robust.

Concerning the CPU time, it seems strongly linked to the localization of the sub-levels. The most interesting in terms of CPU remains to use a very coarse mesh with many sub-levels except if the user tolerance requires this coarse mesh to be totally or quite totally refined. In this case, one has to use a more refined initial mesh to start the LDC algorithm. For example, for a relative error around 3% there is a ratio in CPU time of 3.5 between a LDC multi-grid adaptive refinement from an initial coarse mesh of $h_0 = 200 \ \mu m$ and a uniform fine refinement with a mesh step $h_0 = 25 \ \mu m$.

### 4.4.2. 2D(r,z) frictional test case

The 2D(r,z) test case is now considered with a friction coefficient $\mu = 0.2$. We can observe that the generated sub-levels are similar than without friction, see Fig. 21 with respect to Fig. 19. On Fig. 22 we can observe that contact forces have, as expected, tangential components and a sliding zone is obtained: three status are observed (gap, adherence, sliding). Table 2 sums up the differences made for this test case.

The proposed refinement strategy seems to work very satisfactory for frictional contact problems. The main conclusions to be drawn from Table 2 are similar to those for Table 1:
<figure><img src="image_22.png" /><figcaption>Figure 21: Frictional 2D(r,z) test case: grids superposition for $h_0 = 200 \ \mu m$ and $\alpha = 5 \%$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: Frictional 2D(r,z) test case: contact forces for $h_0 = \mu m$ with a zoom in the red box.</figcaption></figure>

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

<table><caption>Table 2: Frictional 2D(r,z) test case ($\mu = 0.2$): various coarse mesh steps $h_0$ and user tolerances $\alpha$. $N_l$: number of nodes for $l = 0, ..., l^*$; $NN = \sum_{l=0}^{l^*} N_l$: total number of nodes.</caption><thead><tr><th rowspan="2">h<sub>0</sub></th><th rowspan="2">Initial</th><th colspan="3">user tolerance α</th></tr><tr><th>5%</th><th>2%</th><th>1%</th></tr></thead><tbody><tr><td rowspan="4">200 μm</td><td>l<sup>*</sup></td><td>0</td><td>3</td><td>4</td><td>5</td></tr><tr><td>e<sub>h</sub></td><td>19.6%</td><td>3.51%</td><td>1.27%</td><td>0.85%</td></tr><tr><td>N<sub>l</sub></td><td>1836</td><td>1836/2484/4972/12036</td><td>1836/5670/12604/21988/51474</td><td>1836/7070/21666/48594/84550/196574</td></tr><tr><td>NN</td><td>1836</td><td>21328</td><td>93572</td><td>360290</td></tr><tr><td></td><td>CPU</td><td>0.6s</td><td>31s</td><td>252s</td><td>1663s</td></tr><tr><td rowspan="4">100 μm</td><td>l<sup>*</sup></td><td>0</td><td>2</td><td>3</td><td>4</td></tr><tr><td>e<sub>h</sub></td><td>7.72%</td><td>3.44%</td><td>1.27%</td><td>0.85%</td></tr><tr><td>N<sub>l</sub></td><td>7070</td><td>7070/4884/11682</td><td>7070/12604/21988/50728</td><td>7070/21666/48594/84550/196574</td></tr><tr><td>NN</td><td>7070</td><td>23636</td><td>92390</td><td>358454</td></tr><tr><td></td><td>CPU</td><td>2s</td><td>26s</td><td>250s</td><td>1641s</td></tr><tr><td rowspan="4">50 μm</td><td>l<sup>*</sup></td><td>0</td><td>1</td><td>2</td><td>3</td></tr><tr><td>e<sub>h</sub></td><td>6.58%</td><td>3.14%</td><td>1.25%</td><td>0.84%</td></tr><tr><td>N<sub>l</sub></td><td>27738</td><td>27738/11682</td><td>27738/21988/50728</td><td>27738/48594/84550/196574</td></tr><tr><td>NN</td><td>27738</td><td>39420</td><td>100454</td><td>357456</td></tr><tr><td></td><td>CPU</td><td>11s</td><td>68s</td><td>305s</td><td>1651s</td></tr><tr><td rowspan="4">25 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>1</td><td>2</td></tr><tr><td>e<sub>h</sub></td><td>2.96%</td><td>2.96%</td><td>1.20%</td><td>0.83%</td></tr><tr><td>N<sub>l</sub></td><td>109874</td><td>109874</td><td>109874/50456</td><td>109874/84550/196574</td></tr><tr><td>NN</td><td>109874</td><td>109874</td><td>160330</td><td>390998</td></tr><tr><td></td><td>CPU</td><td>88s</td><td>88s</td><td>598s</td><td>1317s</td></tr></tbody></table>

* the LDC process considerably reduces the initial errors adding localized sub-levels;

* the sub-levels generation stops automatically and leads to respect the prescribed user tolerance;

* for a given $\alpha$, similar local meshes are obtained for equivalent local mesh steps;

* there is still an optimal initial mesh step with regards to the CPU time depending on the localization of the sub-levels and the speed of convergence of the $\Lambda$-cycles (related to the number of refinement levels $l^*$).
The coupling LDC-ZZ is now tested for different friction coefficients, see Table 3 for an user tolerance $\alpha = 5\%$. The obtained results are really satisfactory. These results reinforce the robustness of the proposed automatic refinement

<table><caption>Table 3: Frictional 2D(r,z) test case: various coarse mesh steps $h_0$ and friction coefficients $\mu$. User defined tolerance $\alpha = 5\%$; $N_l$: number of nodes for $l = 0, ..., l^*$; $NN = \sum_{l=0}^{l^*} N_l$: total number of nodes.</caption><thead><tr><th rowspan="2">h<sub>0</sub></th><th rowspan="2"></th><th colspan="3">friction coefficient μ</th></tr><tr><th>0.2</th><th>0.6</th><th>1.0</th></tr></thead><tbody><tr><td rowspan="5">200 μm</td><td>l<sup>*</sup></td><td>3</td><td>3</td><td>3</td></tr><tr><td>e<sub>h</sub></td><td>3.51%</td><td>3.43%</td><td>3.59%</td></tr><tr><td>N<sub>l</sub></td><td>1836/2484/4972/12036</td><td>1836/2278/4998/12096</td><td>1836/2208/4522/11426</td></tr><tr><td>NN</td><td>21328</td><td>21208</td><td>19992</td></tr><tr><td>CPU</td><td>31s</td><td>56s</td><td>38s</td></tr><tr><td rowspan="5">100 μm</td><td>l<sup>*</sup></td><td>2</td><td>2</td><td>2</td></tr><tr><td>e<sub>h</sub></td><td>3.44%</td><td>3.34%</td><td>3.49%</td></tr><tr><td>N<sub>l</sub></td><td>7070/4884/11682</td><td>7070/4998/12096</td><td>7070/4760/11542</td></tr><tr><td>NN</td><td>23636</td><td>24164</td><td>23372</td></tr><tr><td>CPU</td><td>26s</td><td>30s</td><td>56s</td></tr><tr><td rowspan="5">50 μm</td><td>l<sup>*</sup></td><td>1</td><td>1</td><td>1</td></tr><tr><td>e<sub>h</sub></td><td>3.14%</td><td>3.26%</td><td>3.32%</td></tr><tr><td>N<sub>l</sub></td><td>27738/11682</td><td>27738/12096</td><td>27738/11542</td></tr><tr><td>NN</td><td>39420</td><td>39834</td><td>39280</td></tr><tr><td>CPU</td><td>68s</td><td>72s</td><td>105s</td></tr><tr><td rowspan="5">25 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>0</td></tr><tr><td>e<sub>h</sub></td><td>2.96%</td><td>3.06%</td><td>3.14%</td></tr><tr><td>N<sub>l</sub></td><td>109874</td><td>109874</td><td>109874</td></tr><tr><td>NN</td><td>109874</td><td>109874</td><td>109874</td></tr><tr><td>CPU</td><td>88s</td><td>95s</td><td>104s</td></tr></tbody></table>

strategy for frictional contact problems as its efficiency does not depend on the friction coefficient. Note that the CPU time globally increases with the friction coefficient. This is due to the frictional contact solver and not to the refinement methodology. We can observe that this intrinsic solver resolution time has an influence on the optimum initial mesh versus number of sub-levels in terms of CPU time.

### 4.4.3. 2D(r,θ) thermal unilateral contact test case

The aim of this test case is to validate our automatic adaptive mesh refinement algorithm on a different configuration (geometry and thermal strain imposed). For simplicity, only a frictionless contact problem is focused on.

**Hierarchical mesh refinement.** We first apply the LDC-contact algorithm with a hierarchical refinement as detailed in section 4.3. For this test case, this strategy leads not to improve the coarsest polygonal approximation of the arc contact boundary. If we look at the contact forces obtained through the LDC algorithm, we obtained non coherent results from the second sub-level, see Figs. 23 to 26 for example. Some contact discontinuities appear on the added fine nodes which are located on the polygon (string of the arc) and no more on the arc. This is not due to the LDC algorithm as the same feature appears on an equivalent whole domain mesh build from a hierarchical mesh refinement procedure,
<figure><img src="image_25.png" /><figcaption>Figure 23: 2D(r,θ) test case: Contact forces for $h_0 = 200$ µm</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 24: 2D(r,θ) test case: Contact forces on the first sub-level ($l = 1$) for $h_0 = 200$ µm and $\alpha = 1\%$. Hierarchical refinement.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 25: 2D(r,θ) test case: Contact forces on the second sub-level ($l = 2$) for $h_0 = 200$ µm and $\alpha = 1\%$. Hierarchical refinement.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 26: 2D(r,θ) test case: Contact forces on the third sub-level ($l = 3$) for $h_0 = 200$ µm and $\alpha = 1\%$. Hierarchical refinement.</figcaption></figure>

see Figure 27. As mentioned in [WS98], discontinuous contacts are in this case only a result of the discretization of both arcs by polygons. However, the contact forces seem consistent if the finite element mesh of the arc contact boundaries has all discretization nodes located on the arc, see Figure 28 for example. As a

<figure><img src="image_29.png" /><figcaption>Figure 27: 2D(r,θ) test case: Contact forces for a whole domain mesh obtained with two hierarchical mesh refinements of $r = 2$ from $h_0 = 200$ µm.</figcaption></figure>

consequence, a hierarchical mesh refinement can not be applied on arc contact boundaries to obtain reliable results.

*Quasi-hierarchical mesh refinement.* To ensure the potential contact boundary nodes to be located on the arc during the refinement process, we propose a quasi-hierarchical mesh refinement procedure:

* perform a hierarchical mesh refinement with refinement ratio $r$, see Figure 12.
<figure><img src="image_30.png" /><figcaption>Figure 28: 2D(r,θ) test case: Contact forces for $h_0 = 50 \ \mu m$.</figcaption></figure>

• move the finite element discretization nodes of the potential contact boundary to the arc, see Figure 29.

<figure><img src="image_31.png" /><figcaption>Figure 29: Example of quasi-hierarchical contact boundary element refinement with $r = 2$</figcaption></figure>

This strategy is slightly different from a hierarchical mesh refinement procedure and hence easy to implement. Moreover this quasi-hierarchical mesh refinement procedure implies a better discretization of the arc contact boundary at each level. The quasi-hierarchical procedure is only applied on contact boundaries in order to conserve unchanged the bilinear interpolation prolongation operator (as $\partial\Omega_C^f = \emptyset$). Moreover, it had been shown in [BRL14] that on natural boundaries ($\partial\Omega_D$ and $\partial\Omega_N$), hierarchical mesh refinement does not induce first order errors. This strategy is also interesting as the canonical restriction operator is also unchanged since the shift nodes are by construction always additional fine nodes. Hence, Algorithm 2 remains generic, the quasi-hierarchical mesh refinement occurring only in the “Discretization of $\Omega_l$” step.

The obtained results are then really satisfactory, see Figure 30 to 32.

<figure><img src="image_32.png" /><figcaption>Figure 30: 2D(r,θ) test case: Contact forces on the first sub-level ($l = 1$) for $h_0 = 200 \ \mu m$ and $\alpha = 1\%$. Quasi-hierarchical refinement.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 31: 2D(r,θ) test case: Contact forces on the second sub-level ($l = 2$) for $h_0 = 200 \ \mu m$ and $\alpha = 1\%$. Quasi-hierarchical refinement.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 32: 2D(r,θ) test case: Contact forces on the third sub-level ($l = 3$) for $h_0 = 200 \ \mu m$ and $\alpha = 1\%$. Quasi-hierarchical refinement.</figcaption></figure>

Table 4 summarizes some post-treatments made for various numerical experiments.

The first conclusion to be drawn from Table 4 is that the automatic refinement process is very consistent as the same local mesh steps and local grids are obtained for a same user tolerance whatever the initial coarse mesh (see columns
<table><caption>Table 4: Thermal unilateral contact 2D(r,θ) test case: various coarse mesh steps $h_0$ and user tolerances $α$. $N_l$: number of nodes for $l = 0, ..., l^*$; $NN = \sum_{l=0}^{l^*} N_l$: total number of nodes. Quasi-hierarchical refinement process.</caption><thead><tr><th rowspan="2">h<sub>0</sub></th><th rowspan="2"></th><th colspan="3">User tolerance α</th></tr><tr><th>5%</th><th>2%</th><th>1%</th></tr></thead><tbody><tr><td rowspan="4">200 μm</td><td>l<sup>*</sup></td><td>0</td><td>2</td><td>3</td></tr><tr><td>e<sub>h</sub></td><td>1.21%</td><td>0.69%</td><td>0.47%</td></tr><tr><td>N<sub>l</sub></td><td>468</td><td>468/650/1274</td><td>468/900/2254/4462</td></tr><tr><td>NN</td><td>468</td><td>2392</td><td>8084</td></tr><tr><td rowspan="4">100 μm</td><td>l<sup>*</sup></td><td>0</td><td>1</td><td>2</td></tr><tr><td>e<sub>h</sub></td><td>0.61%</td><td>0.56%</td><td>0.37%</td></tr><tr><td>N<sub>l</sub></td><td>1750</td><td>1750/1274</td><td>1750/2254/4462</td></tr><tr><td>NN</td><td>1750</td><td>3024</td><td>8466</td></tr><tr><td rowspan="4">50 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>1</td></tr><tr><td>e<sub>h</sub></td><td>0.31%</td><td>0.31%</td><td>0.25%</td></tr><tr><td>N<sub>l</sub></td><td>6762</td><td>6762</td><td>6762/4462</td></tr><tr><td>NN</td><td>6762</td><td>6762</td><td>11224</td></tr><tr><td rowspan="4">25 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>0</td></tr><tr><td>e<sub>h</sub></td><td>0.16%</td><td>0.16%</td><td>0.16%</td></tr><tr><td>N<sub>l</sub></td><td>26578</td><td>26578</td><td>26578</td></tr><tr><td>NN</td><td>26578</td><td>26578</td><td>26578</td></tr></tbody></table>

of Table 4). The second conclusion is that the obtained global errors are always lower than the prescribed tolerance. However, in view of the relative error it seems that some refinement levels were not necessary. Actually, this is due to the global error post-treatment that is not relevant for this thermo-mechanical example. Indeed, the stress field due to the imposed thermal strain is very high, especially at the center of the pellet, but well approximated by the Finite Element Method, which explained why the global relative error is so small. In order to catch the effect of the automatic adaptive mesh refinement on the stress field around the contact zone and following the main idea the Theory of the Critical Distances in fracture mechanics [Tay07] and especially the Line Method, the LDC composite energy error is now evaluated over a zone ZC surrounding the potential contact zone, see Figure 33. The length of this zone is chosen to be equivalent to the length of the coarsest grid cell (mesh step $h_0 = 200$ µm) in each body.

Table 5 reports the LDC relative errors integrated over the domain ZC, called $e_{ZC}$, for this thermal unilateral contact 2D(r,θ) test case:

$$e_{ZC}^2 = \frac{\sum_{l=1}^{l^*+1} \sum_{K \subset T_{l-1} \setminus T_{l \subset ZC}} \int (\sigma_h - \sigma_{ref}) : (\varepsilon_h - \varepsilon_{ref}) dK}{\sum_{l=1}^{l^*+1} \sum_{K \subset T_{l-1} \setminus T_{l \subset ZC}} \int \sigma_{ref} : \varepsilon_{ref} dK} \quad (22)$$

with $T_{l^*+1} = \emptyset$.

The difference between Table 4 and Table 5 only concerns the relative errors.
<figure><img src="image_35.png" /><figcaption>Figure 33: Definition of the zone ZC around the potential contact boundaries for evaluating of the relative LDC energy error.</figcaption></figure>

<table><caption>Table 5: Thermal unilateral contact 2D(r,θ) test case: various coarse mesh steps $h_0$ and user tolerance $α$. $N_l$: number of nodes for $l = 0, ..., l^*$; $NN = \sum_{l=0}^{l^*} N_l$: total number of nodes. Quasi-hierarchical refinement process, LDC composite energy error $e_{ZC}$ on the zone ZC surrounding the potential contact boundaries.</caption><thead><tr><th rowspan="2">h<sub>0</sub></th><th rowspan="2"></th><th colspan="3">User tolerance α</th></tr><tr><th>5%</th><th>2%</th><th>1%</th></tr></thead><tbody><tr><td rowspan="4">200 μm</td><td>l<sup>*</sup></td><td>0</td><td>2</td><td>3</td></tr><tr><td>e<sub>ZC</sub></td><td>4.10%</td><td>1.12%</td><td>0.62%</td></tr><tr><td>N<sub>l</sub></td><td>468</td><td>468/650/1274</td><td>468/900/2254/4462</td></tr><tr><td>NN</td><td>468</td><td>2392</td><td>8084</td></tr><tr><td rowspan="4">100 μm</td><td>l<sup>*</sup></td><td>0</td><td>1</td><td>2</td></tr><tr><td>e<sub>ZC</sub></td><td>2.10%</td><td>1.13%</td><td>0.61%</td></tr><tr><td>N<sub>l</sub></td><td>1750</td><td>1750/1274</td><td>1750/2254/4462</td></tr><tr><td>NN</td><td>1750</td><td>3024</td><td>8466</td></tr><tr><td rowspan="4">50 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>1</td></tr><tr><td>e<sub>ZC</sub></td><td>1.14%</td><td>1.14%</td><td>0.61%</td></tr><tr><td>N<sub>l</sub></td><td>6762</td><td>6762</td><td>6762/4462</td></tr><tr><td>NN</td><td>6762</td><td>6762</td><td>11224</td></tr><tr><td rowspan="4">25 μm</td><td>l<sup>*</sup></td><td>0</td><td>0</td><td>0</td></tr><tr><td>e<sub>ZC</sub></td><td>0.63%</td><td>0.63%</td><td>0.63%</td></tr><tr><td>N<sub>l</sub></td><td>26578</td><td>26578</td><td>26578</td></tr><tr><td>NN</td><td>26578</td><td>26578</td><td>26578</td></tr></tbody></table>

We can see that $e_{ZC}$ enables us to much appreciate the effect of the refinement process around the contact zone. The prescribed user tolerance is still always respected but the necessity of the refinement levels appears clearly with this post-treatment. Table 5 shows moreover the gain in precision and confidence brought by the proposed automatic adaptive mesh refinement method for the prediction of the cladding failure (defined by a critical stress value).

*Remark.* For the previous 2D(r,z) test cases, the main contribution in the global relative error comes from the contact zone. Hence the error $e_{ZC}$ is in this case equivalent to $e_h$.
## 5. Conclusion and perspectives

In this paper, we introduced an innovative adaptive mesh refinement method for frictional contact problems. This method is based on the extension of the Local Defect Correction multigrid method for frictional contact problems. This extension is generic in the sense where it is independent to the solver used (black-box solver), non-intrusive and it only focuses on displacement field interpolations between levels of grids. The LDC algorithm for frictional contact problems is then coupled with the Zienkiewicz and Zhu a posteriori error estimator to automatically detect the high level error zones. This relative stress estimator is applied on each body separately in order to avoid taking into account contact contributions in the ZZ estimator. As the stress field may be singular, a geometric stopping criterion is also used in order to stop the levels generation.

The proposed strategy has been tested with success on various configurations (unilateral frictional or frictionless contact, thermal strain imposed, contact arc boundaries,...). A hierarchical mesh refinement is used between each level except for contact arc boundaries where a local quasi-hierarchical refinement is employed. The obtained subgrids are consistent whatever the initial coarse mesh and are localized around the effective contact zone. The obtained errors rapidly decrease with the addition of local subgrids. The prescribed user tolerance is always respected after the LDC-ZZ algorithm. It had moreover been shown that similarly to fracture mechanics post-treatments, the errors have to be evaluated on a zone surrounding the potential contact boundaries in order to appreciate the effect of the refinement process without being pollute by singular local stresses or imposed strains. The proposed automatic strategy hence guarantees a confidence interval on the calculated stress field around the contact zone, which may be very useful in safety studies for example (cladding failure,...) Concerning the CPU time, there is an optimal strategy in terms of initial coarse mesh step and number of sublevels depending on the localization of the subgrids. When the subgrids are localized, the most interesting in terms of CPU remains to use a coarse initial level and add many sublevels.

Future works concern the extension of the frictional contact LDC algorithm to time dependent problems. Indeed, in this case the zones of interest and hence the automatically generated sublevels may vary during the loading. A special treatment have then to be done for history conservation when nonlinear material properties with internal variables are under study. These questions return to the problematic encountered by all local mesh refinement strategies for nonlinear mechanics time dependent problems.

## Acknowledgement

This work has been achieved in the framework of the collaboration protocol between the CEA (Commissariat à l'Énergie Atomique et aux Energies Alternatives) and the LMA (Laboratoire de Mécanique et d'Acoustique, CNRS, Marseille). The authors are grateful to the PLEIADES project, financially supported by CEA, EDF (Électricité de France) and AREVA, that funded this research work.
References

[ABD16] S. Abide, M. Barboteu, and D. Danan. Analysis of two active set type methods to solve unilateral contact problems. *Applied Mathematics and Computation*, 284:286-307, 2016.

[ABS05] M.J.H. Anthonissen, B.A.V. Bennet, and M.D. Smooke. An adaptive multilevel local defect correction technique with application to combustion. *Combustion Theory and Modelling*, 9:273-299, 2005.

[AC91] P. Alart and A. Curnier. A mixed formulation for frictional contact problems prone to newton like solution methods. *Computer Methods in Applied Mechanics and Engineering*, 92(3):353-375, 1991.

[AL94] P. Angot and M. Laugier. The FIC method of conservative connection between nested subdomains for an ocean circulation model. *Comptes Rendus de l'Académie des Sciences de Paris - Série II - Part 1*, 319(9):993-1000, 1994.

[AL95] P. Alart and F. Lebon. Solution of Frictional Contact Problems Using ILU and Coarse/Fine Preconditioner. *Computational Mechanics*, 16:98-105, 1995.

[AMB03] M.J.H. Anthonissen, R.M.M. Mattheij, and J.H.M.T. Boonkkamp. Convergence analysis of the local defect correction method for diffusion equations. *Numerische Mathematik*, 95(3):401-425, 2003.

[Amo99] G. Amontons. Sur l'origine de la résistance dans les machines. *Mémoire de de l'Académie Royale*, page 206-222, 1699.

[ARF01] H. Askes and A. Rodriguez-Ferran. A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples. *International Journal for Numerical Methods in Engineering*, 51(3):253-273, 2001.

[BG03] M. Belliard and M. Grandotto. Local zoom computation of two-phase flows in steam generators using a local defect correction method. *Numerical Heat Transfer - Part A Applications*, 43(2):111-135, 2003.

[BO84] M.J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. *Journal of Computational Physics*, 53(3):484-512, 1984.

[BR78a] I. Babuška and W. Rheinboldt. Error estimates for adaptive finite element computations. *SIAM Journal on Numerical Analysis*, 15(4):736-754, 1978.

[BR78b] I. Babuška and W.C. Rheinboldt. A-posteriori error estimates for the finite element method. *International Journal for Numerical Methods in Engineering*, 12:1597-1615, 1978.

[BR05] H. Ben Dhia and G. Rateau. The Arlequin method as a flexible engineering design tool. *International Journal for Numerical Methods in Engineering*, 62(11):1142-1162, 2005.
[Bra77] A. Brandt. Multi-level adaptive solutions to boundary-value problems. *Mathematics of Computation*, 31:333–390, 1977.

[Bra94] A. Brandt. Rigorous quantitative-analysis of multigrid .1. constant-coefficients 2-level cycle with l2-norm. *SIAM Journal on Numerical Analysis*, 31(6):1695–1730, 1994.

[BRL14] L. Barbié, I. Ramière, and F. Lebon. Strategies around the local defect correction multi-level refinement method for three-dimensional linear elastic problems. *Computers and Structures*, 130:73–90, 2014.

[BRL15] L. Barbié, I. Ramière, and F. Lebon. An automatic multilevel refinement technique based on nested local meshes for nonlinear mechanics. *Computers & Structures*, 147:14–25, 2015.

[BS82] I. Babuska and B. Szabo. On the rates of convergence of the finite element method. *International Journal for Numerical Methods in Engineering*, 18(3):323–341, 1982.

[BS87] Ivo Babuška and Manil Suri. The h-p version of the finite element method with quasiuniform meshes. *RAIRO-Modélisation mathématique et analyse numérique*, 21(2):199–238, 1987.

[BT93] T. Belytschko and M. Tabbara. H-adaptive finite-element methods for dynamic problems, with emphasis on localization. *International Journal for Numerical Methods in Engineering*, 36(24):4245–4265, 1993.

[BW85] R.E. Bank and A. Weiser. Some a posteriori error estimators for elliptic partial differential equations. *Mathematics of Computation*, 44(170):283–301, 1985.

[CAS] CAST3M. www.cast3m.cea.fr.

[CDR98] P. Chabrand, F. Dubois, and M. Raous. Various numerical methods for solving unilateral contact problems with friction. *Mathematical and Computer Modelling*, 28(4):97–108, 1998.

[CGLC05] P. Cavin, A. Gravouil, A.A. Lubrecht, and A. Combescure. Efficient FEM calculation with predefined precision through automatic grid refinement. *Finite Elements in Analysis and Design*, 41:1043–1055, 2005.

[CHP00] P. Coorevits, P. Hild, and J.P. Pelle. A posteriori error estimation for unilateral contact with matching and non-matching meshes. *Computer Methods in Applied Mechanics and Engineering*, 186(1):65–83, 2000.

[Cou85] C.A. Coulomb. Théorie de machines simples. *Mémoire de Mathématiques et Physique de l'Académie Royale*, 10:161–342, 1785.
[DDO85] L. Demkowicz, P. Devloo, and J.T. Oden. On a h-type mesh-refinement strategy based on minimization of interpolation errors. *Computer Methods in Applied Mechanics and Engineering*, 53(1):67–89, 1985.

[DK08] C. A. Duarte and D.-J. Kim. Analysis and applications of a generalized finite element method with global–local enrichment functions. *Computer Methods in Applied Mechanics and Engineering*, 197(6-8):487–504, 2008.

[DL76] G. Duvaut and J.L. Lions. *Inequalities in mechanics and physics*, volume 219. Springer-Verlag, 1976.

[DNR07] A. Düster, A. Niggl, and E. Rank. Applying the hp-d version of the FEM to locally enhance dimensionally reduced models. *Computer Methods in Applied Mechanics and Engineering*, 196:3524–3533, 2007.

[DORH89] L. Demkowicz, J.T. Oden, W. Rachowicz, and O. Hardy. Toward a universal h-p adaptive finite-element strategy - I. Constrained approximation and data structure. *Computer Methods in Applied Mechanics and Engineering*, 77(1-2):79–112, 1989.

[DR01] A. Düster and E. Rank. The p-version of the finite element method compared to an adaptive h-version for the deformation theory of plasticity. *Computer Methods in Applied Mechanics and Engineering*, 190:1925–1935, 2001.

[EH09] Y. Efendiev and T. Y. Hou. *Multiscale Finite Element Methods: Theory and Applications*. Springer, 2009.

[EOD93] M.G. Edwards, J.T. Oden, and L. Demkowicz. An hr-adaptive approximate Riemann solver for the Euler equations in 2 dimensions. *SIAM Journal on Scientific Computing*, 14(1):185–217, 1993.

[Fis92] J. Fish. The s-version of the finite-element method. *Computers and Structures*, 43(3):539–547, 1992.

[GB05] T. Gratsch and K.J. Bathe. Review - A posteriori error estimation techniques in practical finite element analysis. *Computers and Structures*, 83:235–265, 2005.

[Git72] J.H. Gittus. Theoretical analysis of the strains produced in nuclear fuel cladding tubes by the expansion of cracked cylindrical fuel pellets. *Nuclear Engineering and Design*, 18(1):69–82, 1972.

[GLP00] L. Gallimard, P. Ladevèze, and J.P. Pelle. An enhanced error estimator on the constitutive relation for plasticity problems. *Computers and Structures*, 78:801–810, 2000.

[GM93] S. Ghosh and S.K. Manna. r-adapted arbitrary Lagrangian-Eulerian finite-element method in metal-forming simulation. *Journal of Materials Engineering and Performance*, 2(2):271–282, 1993.
[Hac84] W. Hackbusch. Local Defect Correction Method and Domain De-
composition Techniques. *Computing Suppl. Springer-Verlag*, 5:89-
113, 1984.

[HL09] P. Hild and V. Lleras. Residual error estimates for coulomb friction. *SIAM Journal on Numerical Analysis*, 47(5):3550–3583, 2009.

[HR10] W. Huang and R.D. Russell. *Adaptive moving mesh methods*, volume 174. Springer Science & Business Media, 2010.

[HW97] T. Y. Hou and X.-H. Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. *Journal of Computational Physics*, 134(1):169–189, 1997.

[JD13] M. Joulaian and A. Düster. Local enrichment of the finite cell method for problems with material interfaces. *Computational Mechanics*, 52(4):741–762, 2013.

[KACM96] K. Khadra, Ph. Angot, J.P. Caltagirone, and P. Morel. Concept de zoom adaptatif en architecture multigrille locale ; étude comparative des méthodes L.D.C., F.A.C. et F.I.C. *RAIRO - Modélisation Mathématique et Analyse Numérique*, 30(1):39-82, 1996.

[KCMM09] W. Kramer, H.J.H. Clercx, R.M.M. Mattheij, and R. Minero. A finite volume local defect correction method for solving the transport equation. *Computers and Fluids*, 38(3):533–543, 2009.

[KL11] F. Kuss and F. Lebon. Error estimation and mesh adaptation for Signorini-Coulomb problems using E-FEM. *Computers and Structures*, 89(11-12):1148-1154, 2011.

[KR69] G. Kjaerheim and E. Rolstad. In-Core Study of the Mechanical Interaction between Fuel and Cladding. *Nuclear Technology*, 7(4):347–360, 1969.

[Leb95] F. Lebon. Two-grid method for regularized frictional elastostatics problems. *Engineering Computations*, 12:657–664, 1995.

[Leb03] F. Lebon. Contact problems with friction: models and simulations. *Simulation Modelling Practice and Theory*, 11:449–463, 2003.

[LFPD15] R.M. Lins, M.D.C. Ferreira, S.P.B. Proença, and C.A. Duarte. An a-posteriori error estimator for linear elastic fracture mechanics using the stable generalized/extended finite element method. *Computational Mechanics*, 56(6):947–965, 2015.

[Liu16] H. Liu. *Automatic mesh refinement and local multigrid methods for contact problems : application to the Pellet-Cladding mechanical Interaction*. PhD Thesis, in French, Aix-Marseille Université, 2016.

[LL83] P. Ladevèze and D. Leguillon. Error estimate procedure in the finite-element method and applications. *SIAM Journal on Numerical Analysis*, 20:485–509, 1983.
[LQ92] M. Lemke and D. Quinlan. Fast adaptive composite grid methods on distributed parallel architectures. *Communications in Applied Numerical Methods*, 8(9):609-619, 1992.

[LRL17] H. Liu, I. Ramière, and F. Lebon. An original adaptive mesh refinement strategy in case of local singularities. *to be submitted*, 2017.

[LRR07] F. Lebon, M. Raous, and I. Rosu. Multigrid methods for unilateral contact problems with friction. In P. Wriggers and U. Nackenhorst, editors, *IUTAM Symposium on Computational Contact Mechanics*, pages 1–16. Springer, 2007. 5-8 November, Hannover, Germany.

[MAM06] R. Minero, M.J.H. Anthonissen, and R.M.M. Mattheij. A local defect correction technique for time-dependent problems. *Numerical Methods for Partial Differential Equations*, 22:128–144, 2006.

[McC84] S.F. McCormick. Fast Adaptive Composite Grid (F.A.C.) Methods: theory for the variational case. *Computing Suppl. Springer-Verlag*, 5:115–121, 1984.

[Mor63] J.J. Moreau. Fonctionnelles sous-différentiables. *Comptes-Rendus Académie des Sciences I*, 257:4117-4119, 1963.

[MSTC08] B. Michel, J. Sercombe, G. Thouvenin, and R. Chatelet. 3D fuel cracking modelling in pellet cladding mechanical interaction. *Engineering Fracture Mechanics*, 75:3581–3598, 2008.

[NL92] R.V. Nambiar and K.L. Lawrence. The zienkiewicz-zhu error estimator for multiple material problems. *Communications in Applied Numerical Methods*, 8(4):273-277, 1992.

[NM02] V. Nefedov and R.M.M. Mattheij. Local defect correction with different grid types. *Numerical Methods for Partial Differential Equations*, 18(4):545–468, 2002.

[OWA94] J.T. Oden, W. Wu, and M. Ainsworth. An a posteriori error estimate for finite element approximations of the navier-stokes equations. *Computer Methods in Applied Mechanics and Engineering*, 111(1-2):185-202, 1994.

[PDH06] N. Pares, P. Díez, and A. Huerta. Subdomain-based flux-free a posteriori error estimators. *Computer Methods in Applied Mechanics and Engineering*, 195:297-323, 2006.

[RAB07a] I. Ramière, P. Angot, and M. Belliard. A fictitious domain approach with spread interface for elliptic problems with general boundary conditions. *Computer Methods in Applied Mechanics and Engineering*, 196:766-781, 2007.

[RAB07b] I. Ramière, P. Angot, and M. Belliard. A general fictitious domain method with immersed jumps and multilevel nested structured meshes. *Journal of Computational Physics*, 225:1347-1387, 2007.
[Ram08] I. Ramière. Convergence analysis of the $Q_1$-finite element method for elliptic problems with non-boundary-fitted meshes. *International Journal for Numerical Methods in Engineering*, 75(9):1007-1052, 2008.

[RCL88] M. Raous, P. Chabrand, and F. Lebon. Numerical methods for frictional contact problems and applications. *Journal of Theoretical and Applied Mechanics*, 7:111-128, 1988.

[RZ05] W. Rachowicz and A. Zdunek. An hp-adaptive finite element method for scattering problems in computational electromagnetics. *International Journal for Numerical Methods in Engineering*, 62(9):1226-1249, 2005.

[SB86] T. Sussman and KJ Bathe. Studies of finite element procedures—stress band plots and the evaluation of finite element meshes. *Engineering Computation*, 3:178-191, 1986.

[SH92] T. Strouboulis and K.A. Haque. Recent experiences with error estimation and adaptivity. In Error estimation for h-adaptive approximations on grids of triangles and quadrilaterals. *Computer Methods in Applied Mechanics and Engineering*, 100(3):359-430, 1992.

[Sig59] A. Signorini. Questioni di elasticità non linearizzata e semilinearizzata. *Rendiconti di Matematica e delle sue applicazioni*, 5:95-139, 1959.

[SZ92] W. Sun and N.G. Zamani. An adaptive h-r boundary element algorithm for the Laplace equation. *International Journal for Numerical Methods in Engineering*, 33(3):537-552, 1992.

[Tay07] D. Taylor. *The Theory of Critical Distances : A New Perspectives in Fracture Mechanics*. Elsevier, Oxford, 2007.

[WB98] P. Wriggers and A. Boersma. A parallel algebraic multigrid solver for problems in solid mechanics discretized by finite elements. *Computers & Structures*, 69(1):129-137, 1998.

[Wri06] P. Wriggers. *Computational Contact Mechanics*. Springer Berlin Heidelberg, 2006.

[WS98] P. Wriggers and O. Scherf. Different a posteriori error estimators and indicators for contact problems. *Mathematical and Computer Modelling*, 28(4-8):437-447, 1998.

[YR05] Z. Yue and D.H. Robbins. Adaptive superposition of finite element meshes in elastodynamic problems. *International Journal for Numerical Methods in Engineering*, 63(11):1604-1635, 2005.

[ZM83] O.C. Zienkiewicz and K. Morgan. *Finite elements and approximation*. Wiley, 1983.

[ZZ87] O.C. Zienkiewicz and J.Z. Zhu. A simple error estimator and adaptive procedure for practical engineering analysis. *International Journal for Numerical Methods in Engineering*, 24:337-357, 1987.
[ZZ92a] O.C. Zienkiewicz and J.Z. Zhu. The superconvergent patch recovery and a posteriori error estimation. Part I: The recovery technique. *International Journal for Numerical Methods in Engineering*, 33:1331–1364, 1992.

[ZZ92b] O.C. Zienkiewicz and J.Z. Zhu. The superconvergent patch recovery and a posteriori error estimation. Part II: Error estimates and adaptivity. *International Journal for Numerical Methods in Engineering*, 33:1365–1382, 1992.
Highlights:

* A multilevel mesh refinement strategy for frictional contact problems is introduced

* The automatic adaptive mesh refinement algorithm is non-intrusive

* Zones of interest are detected automatically with respect to a given tolerance

* Relative errors after refinement are lower than the prescribed threshold in stress

* The efficiency of the proposed strategy is analyzed on nuclear engineering examples