Elastodynamic adapted shape derivative
for residual thickness reconstruction in SHM context

Axel Thomas¹,*, Tom Druet¹, Arnaud Recoquillay¹, Bastien Chapuis¹

¹University Paris-Saclay, CEA, List, F-91120, Palaiseau, France

*Email: axel.thomas@cea.fr

**Abstract**

In the context of Structural Health Monitoring (SHM), tomographic methods using guided elastic waves have been studied. However, the existing methods face certain limitations, such as resolution or limited potential application cases, due to the approximations they rely on.

Therefore, this paper introduces an adapted shape derivative method for the elastodynamic case, which allows for the reconstruction of complex structure thicknesses with high precision.

**Keywords:** Elastodynamic Tomography, Structural Health Monitoring (SHM), Guided Waves, Shape Derivative

**1 Introduction**

The remaining thickness of a structure is a critical information for many sensitive industries, such as the thickness of a cooling pipe in a nuclear plant. To this end, SHM aims to establish a regular monitoring of the critical infrastructure. We are particularly interested in tomographic reconstruction methods which used guided elastic waves.

In the literature, the most studied methods are based on an acoustic model. Therefore, the underlying physics relies on approximations to enable rapid thickness estimation. For instance, the propagation considers only one mode, or there is no mode conversion at the interfaces. While these assumptions are valid in many cases, they inevitably affect the resolution and the potential applications. Significantly, these assumptions are generally only valid if the wave packet of interest can be extracted from the rest of the signal.

In this approach, we present an adaptation of the shape derivative to eliminate the need of previous hypothesis. This method attempts to iteratively solve the inverse problem by searching for an optimal study domain shape. The main limitation of such a method is the need to solve the elastodynamic system for each iteration, inducing a large computation burden

in practice. We propose to alleviate this issue with two tools: first, we will use an acoustic tomography algorithm to obtain a good initial guess, thus limiting the number of iterations and the sensitivity to local minima. Second, the resolution of the elastodynamic problem is based on spectral finite elements, reducing the cost of each iteration.

**2 Formalism**

The mathematical problem is based on static shape derivative described in Ref. [3] and adapted to the dynamic case.

We consider a reference regular open bounded domain $\Omega_0$. The boundary variation method modifies the domain such as $\Omega_\theta := (1 + \theta) \Omega$, where $\Omega$ a given domain and $\theta$ a mathematical object that represents a vector field used to update the domain. In order to obtain the medium response to the propagation of guided waves, we need to solve the elastodynamic system defined such as $\mathbf{u}(\mathbf{x}, t)$ satisfies

$$ 
\left\{ 
\begin{array}{ll} 
\rho(\mathbf{x})\partial_{tt}\mathbf{u}(\mathbf{x}, t) - \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}(\mathbf{x}, t)) = 0 & \text{in } \Omega \times [0, T] \\ 
\boldsymbol{\sigma}(\mathbf{u}(\mathbf{x}, t)) = \mathbf{n} = \mathbf{g}(\mathbf{x}, t) & \text{on } \Gamma_N \times [0, T] \\ 
\mathbf{u}(\mathbf{x}, t) = \partial_t \mathbf{u}(\mathbf{x}, t) = 0 & \text{in } \Omega \times \{0\}. 
\end{array} 
\right. 
\quad (1) 
$$

where $\boldsymbol{\sigma}(\mathbf{u}) = \mathbf{C} : \boldsymbol{\varepsilon}(\mathbf{u})$, $\rho(\mathbf{x})$ is the density, $\mathbf{g}(\mathbf{x}, t)$ is the surface load, $\Gamma_N$ is the Neumann boundary.

The spectral finite element solver of CIVA SHM software is used. This solver is based on an explicit leapfrog scheme in time and spectral finite elements in space. Compared to the traditional finite element method, the spectral finite elements method leads to mass lumping. The linear system to be solved implies a diagonal mass matrix. Thus, the inversion at each time step as a very small computational cost and the system does not need to be assembled, making it even more efficient.

In practice, we used the least-square minimization to solve this inverse problem. The associated $L_2$ shape derivative cost function is formulated as
$$J(\Omega) = \frac{1}{2} \iint_{\Omega \times [0, T]} \|\mathbf{u}_{\Omega} - \mathbf{u}_{\text{obs}}\|^2 dx dt, \quad (2)$$

where $\mathbf{u}_{\Omega}$ and $\mathbf{u}_{\text{obs}}$ are, respectively, the elastic displacement resulting from the resolution of (1) and the measured displacement.

Finally, the adjoint state method is used to compute the associated gradient of (2) with respect to the shape. It is expressed as

$$\nabla J(\Omega)(\boldsymbol{\theta}) = - \iint_{\partial\Omega \times [0, T]} \boldsymbol{\theta} \cdot \mathbf{n}(\rho v \partial_{tt} \mathbf{u} + \boldsymbol{\varepsilon}(\mathbf{u}) : \mathbf{C} : \boldsymbol{\varepsilon}(\mathbf{v})) \, ds \, dt \quad (3)$$

where $\mathbf{v}$ is the adjoint variable that represents the elastic displacement obtained by solving the elastodynamic adjoint formulation.

As (1) is self-adjoint, the formulation of its adjoint state is identical, but with a different source term. However, to maintain the symmetry in the adjoint operator, the initial condition in (1) becomes a final condition in the adjoint formulation, making it unsolvable directly. Therefore, the equation is computed in inverse time, and the result is obtained in time due to the reversibility property of waves. According to Eq. (3), we see that the gradient only depends on the normal component of the boundary, this means that the surface is moved in the opposite direction to the gradient, leading to a reconstruction of the domain thickness. In practice, this step involves moving the vertices of the mesh along the $\boldsymbol{\theta} \cdot \mathbf{n}$ vector at each iteration. The next iteration is then computed on the updated domain $\Omega_{\theta}$ and repeated accordingly. The shape of the domain represents the optimization variable of the problem through its discretization.

The quasi-newton methods, namely BFGS, linked to the Armijo criterion are used to compute a descent direction from the value of the computed gradient (3) and an adapted descent step, enabling rapid and robust convergence.

An acoustic diffraction tomography was used [2] to perform the domain thickness initialization.

## 3 Numerical Results

The results of the adapted shape derivative algorithm on simulated forward data are being discussed. The domain is a bounded aluminium plate of size 600 × 580 × 3 mm³. To solve the direct and adjoint formulation, the domain is discretised in a form of structured hexahedral mesh

<figure><img src="image_1.png" /><figcaption>Figure 1: Thickness maps with a mesh size of $l = \lambda/(5\sqrt{2})$. The order of the spectral elements is 2 in the two direction of the plane and 3 in the thickness direction. A flaw of radius $r_{\text{flaws}} = 1.5\lambda$, with $\lambda = 25.97$mm, composed by a flat base of radius $r_{\text{base}} = \frac{3}{4}\lambda$ and a smooth transition zone (described by a Hann half-window) of radius $r_{\text{Hann}} = \frac{3}{4}\lambda$, is positioned at coordinates [370, 270] mm, see Fig. 1b. The emission signals are 5-cycles Tone bursts centred at the frequency of $f = 40$ kHz with a Gaussian shape in space polarized normal to the plane and having a standard deviation of $\sigma = 2.5$mm. The simulated sequence is a round robin capture of 30 punctual emissions placed on a circle of radius $r_{\text{circle}} = 150$mm.</figcaption></figure>

The shape derivative was computed with an initial guess being the acoustic algorithm results. Figure 1a presents the shape derivative results after 8 iterations. The results are promising. Indeed, the estimated position, size, and residual thickness are close to the real ones.

## References

[1] A. Imperiale and E. Demaldent, “A macro-element strategy based upon spectral finite elements and mortar elements for transient wave propagation modeling. Application to ultrasonic testing of laminate composite materials,” *International Journal for Numerical Methods in Engineering*, vol. 119, pp. 964–990, 2019.

[2] T. Druet, J.-L. Tastet, B. Chapuis, and E. Moulin, “Autocalibration method for guided wave tomography with undersampled data,” *Wave Motion*, vol. 89, pp. 265–283, Jun. 2019. DOI: 10.1016/j.wavemoti.2019.04.002.

[3] G. Allaire, F. Jouve, and A.-M. Toader, “Structural optimization using sensitivity analysis and a level-set method,” *Journal of Computational Physics*, vol. 194, no. 1, pp. 363–393, Feb. 2004. DOI: 10.1016/j.jcp.2003.09.032.