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

HAL
open science

# A fictitious domain method for unilateral contact problems in non-destructive testing

Eliane Bécache, Patrick Joly, Gilles Scarella

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

## To cite this version:

Eliane Bécache, Patrick Joly, Gilles Scarella. A fictitious domain method for unilateral contact problems in non-destructive testing. First MIT Conference on Computational Fluid and Solid Mechanics, Jun 2001, Cambridge, MA, United States. hal-01529185

HAL Id: hal-01529185

https://hal.science/hal-01529185

Submitted on 30 May 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A fictitious domain method for unilateral contact
problems in non-destructive testing.

E. Bécache†¶ P. Joly‡¶ G. Scarella§¶

**Abstract**

In this work, we present a numerical method for solving the diffraction of transient elastic
waves by cracks of arbitrary shapes in complex media, with Signorini's boundary conditions
on the crack. We use a fictitious domain method based on a mixed displacement-stress
formulation for elastodynamics. We propose an off-centered time discretisation scheme for
enforcing the stability.

*Keywords:* elastodynamics, unilateral contact, fictitious domain method, non-destructive testing, cracks

# 1. Introduction

In this paper, we are interested in solving the diffraction of transient elastic waves by cracks
of arbitrary shapes in complex media, with Signorini's boundary conditions on the crack.
This is the continuation of a previous work [1] done on the linear problem, that is when
the boundary condition on the crack is a free surface boundary condition. To get an effi-
cient method, we want to use regular meshes and at the same time respect the geometry
of the crack. This is possible thanks to the fictitious domain method, which takes into
account the boundary condition via a Lagrange multiplier defined on the crack, which can

<footnote>†(eliane.becache@inria.fr)</footnote>

<footnote>‡(patrick.joly@inria.fr)</footnote>

<footnote>§(gilles.scarella@inria.fr)</footnote>

<footnote>¶INRIA, Domaine de Voluceau-Rocquencourt, BP 105, F-78153 Le Chesnay Cédex</footnote>
be interpreted as the jump of the displacement through the crack. This allows to work with a uniform mesh in the whole domain and an independent mesh on the crack. In order to consider the unilateral contact boundary condition as a constraint, we are led to use the mixed displacement-stress formulation for elastodynamics. We will present a fictitious domain formulation of this problem in which the boundary conditions are taken into account by a variational inequality for the Lagrange multiplier. For the space discretisation of this problem, we propose to use the mixed finite element using spaces of symmetric tensors for the stress[1]. This choice was shown to allow the obtention of an explicit time discretization scheme (mass-lumping) in the linear case. In the non-linear case, we cannot use a centered difference scheme for the time discretisation which would lead to an unconditionally unstable scheme. That is why we propose an off-centered scheme that we show to be stable. This scheme is explicit in the volume unknowns (displacement and stress) but implicit in the Lagrange multiplier: one has to solve an optimisation problem with bound constraints at each time step.

## 2. Presentation of the dynamic unilateral contact problem

We want to solve the displacement-stress formulation of elastodynamics in a domain $\Omega = C\setminus\Gamma \in \mathbb{R}^2$, where $C$ is a rectangle (see Fig. 1):

$$ \left\{ \begin{array}{ll} \rho \frac{\partial^2 u}{\partial t^2} - \operatorname{div} \sigma & = f \quad \text{in } \Omega \times ]0, T[, \\ A\sigma & = \varepsilon(u) \quad \text{in } \Omega \times ]0, T[, \end{array} \right. \qquad (1) $$

where $u$ is the displacement field, $\sigma$ the stress tensor, and $\varepsilon(u)$ the strain tensor defined as $\varepsilon_{ij}(u) = (\partial_i u_j + \partial_j u_i)/2$. We add to (1) Dirichlet boundary conditions on the exterior boundary: $u = 0$ on $\partial C \times ]0, T[$, Signorini's boundary conditions without friction on
the crack Γ, [3]:

$$
\left\{
\begin{array}{ll}
[u_N]_{\Gamma} & \ge 0 \text{ on } \Gamma, \\
\\
\sigma_N & \le 0 \text{ on } \Gamma, \\
\\
\sigma_N [u_N]_{\Gamma} & = 0 \text{ on } \Gamma, \\
\\
[\sigma n]_{\Gamma} & = 0 \text{ on } \Gamma, \\
\\
\sigma_T & = 0 \text{ on } \Gamma.
\end{array}
\right.
\qquad (2)
$$

where, if $n$ denotes the unit normal to $\Gamma$, we set: $u_N = u.n$, $\sigma_N = \sigma n.n$, $\sigma_T = \sigma n - \sigma_N n$
and prescribed initial data that we will systematically omit.

3. A fictitious domain formulation

The fictitious domain method consists in extending the two unknowns (u, σ) in the whole
domain C and in introducing Lagrange multipliers to take into account the Signorini's
boundary conditions (2). Setting G = H<sup>1/2</sup><sub>00</sub>(Γ), we introduce

$$
\begin{align*}
M &= \{u \in [L^2(C)]^2\}, \\
X &= \{\tau \in [L^2(C)]^4, \operatorname{div} \tau \in [L^2(C)]^2 / \tau_{ij} = \tau_{ji}\}, \\
L_N &= H_{00+}^{1/2}(\Gamma) = \{\mu_N \in \mathcal{G} / \mu_N \ge 0 \text{ a.e. on } \Gamma\}, \\
L_T &= [H_{00T}^{1/2}(\Gamma)]^2 = \{\mu_T \in \mathcal{G}^2 / \mu_T \cdot n = 0\}
\end{align*}
$$

The fictitious domain formulation consists in finding $(\sigma, u, \lambda_N, \lambda_T) : ]0, T[ \to X \times M \times L_N \times L_T$

$$
\left\{
\begin{array}{ll}
a(\sigma, \tau) + d(\tau, u) + b_T(\tau, \lambda_T) + b_N(\tau, \lambda_N) & = 0 \quad \forall \tau \in X, \\
(\rho \frac{\partial^2 u}{\partial t^2}, v) - d(\sigma, v) & = (f, v) \quad \forall v \in M, \\
b_T(\sigma, \mu_T) & = 0 \quad \forall \mu_T \in L_T, \\
b_N(\sigma, \mu_N - \lambda_N) & \le 0 \quad \forall \mu_N \in L_N.
\end{array}
\right.
\qquad (3)
$$
with

$$
\begin{align*}
a(\sigma, \tau) &= \int_C A\sigma : \tau \,dx, \quad d(\tau, w) = \int_C w \cdot \operatorname{div} \tau \,dx, \\
b_T(\tau, \mu_T) &= \langle \tau_T, \mu_T \rangle_{L'_T, L_T}, \\
b_N(\tau, \mu_N) &= \langle \tau_N, \mu_N \rangle_{G', G}.
\end{align*}
$$

The Lagrange multipliers can be interpreted as $\lambda_N = [u_N]$ and $\lambda_T = [u_T]$, with $u_T = u - u_N n$.

4. Discretisation

4.1. Semi-discretisation in space

We now introduce some finite dimensional spaces $X_h \subset X$, $M_h \subset M$, $G_H \subset G$ and $(L_T)_H \subset L_T$ of dimensions respectively $N_X, N_M, N_G, N_L^T$. The semi-discretisation in space can then be written in a matricial form as: finding $(U, \Sigma, \Lambda_N, \Lambda_T) \in \mathbb{R}^{N_M} \times \mathbb{R}^{N_X} \times \mathbb{R}_+^{N_G} \times \mathbb{R}_-^{N_L^T}$ such that

$$
\left\{
\begin{array}{l}
M_{\sigma} \Sigma + D^{*} U + B_{T}^{*} \Lambda_{T} + B_{N}^{*} \Lambda_{N} = 0, & (i) \\
M_{v} \frac{d^{2} U}{dt^{2}} - D \Sigma = F, & (ii) \\
B_{T} \Sigma = 0, & (iii) \\
(B_{N} \Sigma, \mu_{N} - \Lambda_{N}) \leq 0, \quad \forall \mu_{N} \in \mathbb{R}_{+}^{N_{G}}. & (iv)
\end{array}
\right.
\qquad (4)
$$

In practice, and this is the interesting point in the fictitious domain method, we introduce
two meshes: the volumic unknowns $U$ and $\Sigma$ are defined on a regular grid, $\mathcal{T}_h$ made of
squares $K_j$ of size $h$ while the surfacic unknowns $\Lambda_N$ and $\Lambda_T$ are computed on a nonuniform
mesh on $\Gamma$, $\mathcal{T}_H$ made of segments $S_j$ of size $H_j$, $H = \sup_j H_j$, see Fig.2.

**Remark 1** The inequality (4)-(iv) can be reinterpreted as

$$
\Lambda_N = \Pi(\Lambda_N + B_N \Sigma) \tag{5}
$$

where $\Pi$ is the orthogonal projection on $\mathbb{R}_+^{N_G}$.
**Choice of the finite elements.** We intend to use the same discretisation than for the linear problem (see [1]). For the lowest order element, this choice corresponds to :

$$X_h = \{\sigma_h \in X/\forall K \in \mathcal{T}_h, \sigma_h|_K \in (Q_1(K))^4\}$$

$$M_h = \{v_h \in M/\forall K \in \mathcal{T}_h, v_h|_K \in (Q_0(K))^2\}$$

$$\mathcal{G}_H = \{\mu_H \in \mathcal{G}/\forall S \in \mathcal{T}_H, \mu_H|_S \in P_1(S)\}$$

Its main interest is that it leads to block diagonal mass matrices (even diagonal for $M_v$) so that $M_v$ and $M_\sigma$ are very easy to invert.

## 4.2. The fully discretised scheme.

It would be tempting to discretize (4) using centered finite difference operators, for instance:

$$\frac{\Lambda_N^{n+1} + \Lambda_N^{n-1}}{2} = \Pi(\Lambda_N^n + B_N \Sigma^n)$$

which would give an explicit scheme. However, one can show that this choice leads to an unconditionally unstable scheme! That is why we propose the following off-centered scheme:

$$\left\{
\begin{array}{ll}
M_{\sigma} \Sigma^n + D^* U^n + B_T^* \Lambda_T^n + B_N^* \Lambda_N^n = 0, & (i) \\
M_v \frac{U^{n+1} - 2U^n + U^{n-1}}{\Delta t^2} - D \Sigma^n = F^n, & (ii) \\
B_T \Sigma^n = 0, & (iii) \\
\Lambda_N^{n+1} = \Pi(\Lambda_N^{n+1} + B_N \frac{\Sigma^n + \Sigma^{n+1}}{2}) & (iv)
\end{array}
\right.
\qquad (6)$$

Note that if the mass matrices are block diagonal, this scheme is only implicit in $\Lambda_N$ and is explicit in the other unknowns. We can show a stability result:

**Theorem 1** The scheme (6) is stable under the usual CFL stability condition

$$\frac{\Delta t^2}{4} \|D^*D\| \le 1, \quad \text{with} \quad \|D^*D\| = \sup_{\Sigma} \frac{(D\Sigma, D\Sigma)}{(M_{\sigma}\Sigma, \Sigma)} \qquad (7)$$

For proving this result we show the decay of an energy. The precise result is the following.
We set

$$ 
\begin{aligned}
V^{n+1/2} &= \frac{U^{n+1} - U^n}{\Delta t}, \\
E_1^{n+1} &= \frac{1}{2} \left( (M_\sigma \Sigma^{n+1}, \Sigma^{n+1}) + (M_v V^{n+3/2}, V^{n+1/2}) \right).
\end{aligned} 
\qquad (8) $$

This quantity is an energy under the CFL condition (7), and one has the identity

$$ E_1^{n+1} - E_1^n = \Delta t \left( \frac{\Lambda_N^n - \Lambda_N^{n+1}}{\Delta t}, B_N \frac{\Sigma^{n+1} + \Sigma^n}{2} \right), $$

which shows that $E_1^{n+1} \le E_1^n$ thanks to 6-(iv).

## 5. Numerical aspects

The implementation of the method amounts to combine an explicit scheme for the unknowns $U$ and $\Sigma$ with an optimisation problem (quadratic functional with bounds constraints) to be solved at each time step for the unknown $\Lambda$. The algorithm has been tested in 1D (comparison with analytic solutions). We are currently developing a 2D code, the optimisation procedure being handled by an algorithm combining the active set method with gradient projection method, [2]. Numerical results will be presented at the conference.

## References

[1] E. Bécache, P. Joly, and C. Tsogka. Fictitious domains, mixed finite elements and perfectly matched layers for 2d elastic wave propagation. to appear in J. Comp. Acous. (Tech. Report INRIA 3889, 2000).

[2] J. Nocedal and S. J. Wright. Numerical Optimization. Springer, 1999.

[3] J.R. Willis and V.P. Smyshlyaev. Effective relations for nonlinear dynamics of cracked solids. J.Mech.Phys.Solids, 44(1):49-75, 1996.
<figure><img src="image_4.png" /><figcaption>Figure 1: The geometry of the problem</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: The two meshes</figcaption></figure>