# Elastodynamic shape derivative: focus on the sampling of a circular distribution

Axel Thomas, Tom Druet, Arnaud Recoquillay, Bastien Chapuis

► To cite this version:

Axel Thomas, Tom Druet, Arnaud Recoquillay, Bastien Chapuis. Elastodynamic shape derivative: focus on the sampling of a circular distribution. 11th European Workshop on Structural Health Monitoring (EWSHM 2024), Jun 2024, Postdam, Germany. pp.ISSN 1435-4934, 10.58286/29844. cea-04746274

HAL Id: cea-04746274

https://cea.hal.science/cea-04746274v1

Submitted on 21 Oct 2024

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
# Elastodynamic shape derivative: focus on the sampling of a circular distribution

Axel THOMAS¹, Tom DRUET², Arnaud RECOQUILLAY³, Bastien CHAPUIS⁴
¹, ², ³, ⁴ CEA LIST, Paris-Saclay, France, [axel.thomas@cea.fr](mailto:axel.thomas@cea.fr)

**Abstract.** In the field of Structural Health Monitoring (SHM), the use of tomographic methods employing guided elastic waves has been extensively investigated. Despite their performance, these existing methods are limited by issues such as resolution and their applicability to only specific structures, caused by the approximations on which they are based. This paper presents an adaptation of the shape derivative method for elastodynamic case, enabling the high-precision reconstruction of complex structural thicknesses. Moreover, as of the number of sensor is of utmost importance in the context of SHM, the paper will focus on reconstructions with under-sampled data.

**Keywords:** Shape derivative, Tomography, Elastic Guided Waves, SHM

## Introduction

The health of structures is crucial for safety in sensitive sectors, such as monitoring the residual thickness of structural components in nuclear power plants. Structural Health Monitoring (SHM) aims to continuously oversee these critical infrastructures. A promising physical mean to do so is the use of ultrasonic guided waves, which propagate over large distances, enabling a sparse array of sensors.

Tomography methods are employed for structural imaging. This implies the physical acquisition of data within the domain, which is then processed to generate an image. The processing is frequently presented in the form of an inverse problem. In order to resolve, it is necessary to model the considered physics to a greater or lesser degree of detail. To use elastic guided waves, existing literature predominantly explores methods based on acoustic modelling [1], where simplifying assumptions are made to facilitate rapid thickness estimations. These include single-mode propagation and the absence of mode conversion at interfaces. Although these approximations are effective under specific conditions, they can limit both the resolution and applicability of the methods. It is notable that these assumptions are generally applicable when the wave packet of interest can be clearly distinguished from the remainder of the signal and the structure can be considered as a perfect waveguide.

In this paper, we introduce an application of the shape derivative method, designed to eliminate the reliance on prior assumptions. The shape derivative method is an iterative
method, which it resolves an inverse problem. The objective is to minimize a cost function established between observation data and synthetic data. The particularity of this method lies in the fact that the optimization process is directly based on the form of the domain. The optimal shape is such that it minimises the associated cost function through the resolution of a partial differential equation on the domain $\Omega$. We propose an adaptation to the elastodynamic case. 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 the use of a spectral finite elements solver, reducing the cost of each iteration.

The SHM context implies to pay attention to the intrusiveness. Indeed, the number of sensors necessary to monitoring the structures must not compromise its integrity. A consequence of relying on waves is that in order to reconstruct properly a defect inside a sensor distribution, it is necessary to respect a spacing of less than half the wavelength used [2]. In practice, for tomography configurations this represents a large number of sensors. It is therefore common to work with under-sampled data [3] in the guided wave tomography method based on an acoustic assumption. This paper investigates this for the shape derivative. The first part of the paper outlines the theoretical basis of shape derivative and spectral finite element method, while the second part evaluates the performance of shape derivative with under-sampled data. The behaviour of the algorithm will be analysed when the number of sensors falls below the minimum sampling criterion and compared to the literature.

# 1. Method

The concept of the shape derivative method is introduced here. The method presented in this paper is based on the shape derivative existing in the literature for the static case [4] [5]. The so-called shape optimisation problem is defined as $\min_{\Omega} J(\Omega)$, for a given cost function $J$. The objective is to minimize the function $J(\Omega)$ over the domain $\Omega$. Frequently, the cost function $J(\Omega)$ depends on $\Omega$ through the resolution of a partial differential equation, namely the state equation.

In the case of guided elastic waves, the state equation system under consideration is such that $u(x, t)$ satisfies

$$ egin{cases} 
\rho(x) \partial_{tt} u(x, t) - \nabla \cdot \sigma(u(x, t)) = 0, \quad \text{in } \Omega \times [0, T] \\ 
\sigma(u(x, t)) n = g(x, t), \quad \text{on } \partial\Omega \times [0, T] \\ 
u(x, 0) = \partial_t u(x, 0) = 0, \quad \text{on } \Omega, 
\end{cases} 
\quad (1) $$

where $u(x, t)$ is the elastic displacement, $\sigma(x, t)$ is the stress tensor which satisfies the Hooke's law such that $\sigma(x, t) = C: \varepsilon(x, t)$ where $\varepsilon(x, t) = (\nabla u(x, t) + \nabla^T u(x, t))$ is the strain tensor and $C$ is the Hooke's tensor. The term $g(x, t)$ is the surface loads applied on the boundary $\partial\Omega \in \Omega$, $n$ is the unit outward normal to $\partial\Omega$ and $T$ is the maximum observation time.

The objective is to modify the surface $\partial\Omega \in \Omega$ to match the observation and hence recover the domain's thickness. The information obtained by the $N_s$ sensors is the elastic displacement $u_{\text{obs}}$ and the domain $\Omega$ gives the synthetic solution $u_{\Omega}$ through the resolution of (1). Therefore, the least square minimization is a suitable choice to fit them. The associated $L_2$ cost function is defined as

$$ J(\Omega) = \frac{1}{2} \sum_{n=0}^{N_s} \int_{0}^{T} \left\| d\mathbf{u}_{\Omega}^{n} - \mathbf{d}_{\text{obs}}^{n} \right\|^{2} dt, $$
where $n$ designates the $n$-th illumination of the domain, $d$ is the observation function reporting, through the state, the part corresponding to the measures and $d_{\text{obs}}$ is the physical measures. The function $d$ hence represents the measurement by the sensor. For example, for a simple sensor, it can be an averaging of the normal displacement over the surface of the sensor. In this paper, we consider only simulated data, which means that the same observation function $d$ is used to generate $d_{\text{obs}}$.

## 1.1 Boundary Variation Method

The boundary variation method, introduced by Hadamard [6] in 1908, lies in the variations of a given domain $\Omega \in \mathbb{R}^3$ such as

$$ \Omega_{\theta} := (1 + \theta)(\Omega), $$

where $\theta: \mathbb{R}^3 \to \mathbb{R}^3$ is a mathematical object that represents a vector field used to update the domain $\Omega$. It gives the whole deformation of the domain in order to reduce the cost function under the given constraints.

The shape derivative $\theta \mapsto \nabla J(\Omega)(\theta)$ of a domain's function $J$ in $\Omega$ is defined as the differential at $\theta = 0$ of the application $\theta \mapsto J(\Omega_{\theta})$, such that its Taylor's expansion is

$$ J(\Omega_{\theta}) = J(\Omega) + J'(\Omega)(\theta) + o(\theta), \quad \text{where} \quad \frac{|o(\theta)|}{||\theta||} \xrightarrow{\theta \to 0} 0. $$

The shape derivative $J'(\Omega)(\theta)$ is necessary to develop an optimization algorithm based on the domain's shape. If we choose a vector field $\theta$ such as $J'(\Omega)(\theta) < 0$, we have, for an $\alpha > 0$ sufficiently small

$$ J(\Omega_{\theta}^{\alpha}) \approx J(\Omega) + \alpha J'(\Omega)(\theta) < J(\Omega). \quad (2) $$

According to the equation (2), information about the gradient of $J$ is essential for minimisation. As demonstrated, the relationship between the domain's shape and its cost function is not straightforward. To derive the gradient, a numerical approximation method, such as the finite difference method, may be employed. Nevertheless, this method results in a computational cost that is proportional to the problem's dimension, making it a significant burden for high-dimensional problems. To address this issue, the adjoint state method is a suitable choice.

## 1.2 Adjoint Method to Shape Gradient Derivation

The adjoint state method is a general variational method to compute the gradient of a functional that depends on a set of state variables, which is the solution of the state equations. The adjoint state variable $v_{\Omega}$ is the solution of an adjoint state system. A didactic introduction to this method is given in Ref. [7].

The literature has shown that in most cases, with the adjoint state method, the shape derivative of the function $J(\Omega)$ defined on the domain $\Omega$ has the following form:

$$ \nabla J(\Omega)(\theta) = \int_{\partial \Omega} \boldsymbol{p}_{\Omega} \boldsymbol{\theta} \cdot \boldsymbol{n} \, \mathrm{d}\boldsymbol{s}, \quad (3) $$
where $p_\Omega: \partial\Omega \to R$ is a function that depends on the cost function $J$, the state variable $u_\Omega$ and the adjoint state variable $v_\Omega$.

The adjoint state formulation of (1) is such that $v(x, t)$ satisfies

$$ egin{cases} 
\rho(x)\partial_{tt}v(x,t) - \nabla \cdot \sigma(v(x,t)) = 0, & x, t \in \Omega \times [0, T] \\ 
\sigma(v(x,t)) \cdot n = \partial_u J, & x, t \in \partial\Omega \times [0, T] \\ 
v(x,T) = \partial_t v(x,T) = 0, & x \in \Omega. 
\end{cases} 
\quad (4) $$

We can then show that

$$ p_{\Omega} = \int_{0}^{T} \left( \rho(x)v_{\Omega}(x, t)\partial_{tt}u_{\Omega}(x, t) + \varepsilon(u_{\Omega}(x, t)) : C : \varepsilon(v_{\Omega}(x, t)) \right) dt. $$

Given that the elastodynamic system is self-adjoint, the adjoint formulation (4) is the same as the state equation (1) but with a different source term. Indeed, the derivative of $J$, with respect to the state variable $u$, is the residual between the observation $u_{obs}$ and the simulated signals $u_\Omega$. Moreover, the adjoint system has a final condition in time but, the elastodynamic problem being revertible in time, it is possible to solve it as usually, considering the final condition as an initial condition and reverting time afterwards. The resulted field is the so-called residual backpropagated field.

Finally, the analytic form (3) allows to identify a descent direction by choosing a vector field $\theta$ such as

$$ \theta = -p_{\Omega}n \quad \text{on} \quad \partial\Omega. $$

The equation (2) becomes

$$ J(\Omega_{\theta}^{\alpha_0}) \approx J(\Omega) - \alpha_0 \int_{\partial\Omega} p_{\Omega}^2 ds < J(\Omega_0) $$

and the updated domain satisfies

$$ \Omega_{\theta}^{\alpha} = (1 + \alpha_{0}\theta) \Omega_{0}, \text{ with } \alpha_{0} > 0. $$

Using the descent direction $\theta$ as the inverse of the gradient is the same as the steepest descent method, which suffers from a very low convergence rate, all the more with a fixed descent step. To address this issue, we have chosen to use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method with the first weak Wolfe condition (Armijo condition) [8]. In this context, the BFGS algorithm is particularly interesting and combined with the Armijo condition, it ensures a robust and fast convergence.

The convergence criterion is defined as

$$ |\nabla J_{k-1} - \nabla J_k| \le \varepsilon \nabla J_0, $$

where $\varepsilon = 10^{-3}$, the subscript $k$ represents the current iteration and $0$ represents the first iteration.
## 1.2 Adapted Shape Derivative Based on Spectral Element Method

As it has been shown previously, the method still needs to solve $2N_s$ problems for each iteration, where $N_s$ the number of sensors, hence the need of an efficient solver. To optimize this step, the resolution is performed with the spectral finite elements in space and an explicit leapfrog scheme in time. The solver is presented in the Ref. [9] and used in the SHM module of CIVA [10].

The particularity of the spectral element method is its capacity to work with Legendre polynomial to perform the numerical integration inside the mesh element. Increase the order of the polynomial allows for a better consideration of the small and complex variations of the field. Thus, compared to adding elements, this method is more efficient in terms of computational cost and memory management because it limits the increase in the number of global degrees of freedom. Compared to the traditional finite element method, the linear system to be solved in this case involves a diagonal mass matrix. This reduces the cost of solving the system, remove the need to assemble the full system and enables parallel implementation of the calculations. In addition, the method demonstrates high accuracy in the context of guided wave propagation. The performance of this solver for wave propagation has been evaluated in this article [11] and shows remarkable performance compared to other finite element solvers. This makes it a wise choice for our application.

## 2. Numerical results

In this paper, the influence of the number of sensors was investigated. The literature demonstrates [2], in the case of guided wave tomography with an acoustic assumption, that in order to capture all spatial frequency information of a defect with an effective radius $r_0$ placed at the centre of a circular distribution, the distance between two adjacent sensors must satisfy

$$d_{\max} < \frac{\lambda R}{2r_0},$$

where $\lambda$ is the working wavelength, $R$ is the radius of the sensor distribution and $r_0$ is the radius of a defect placed at the centre of the distribution.

In our case, it can be seen as the largest imaging area without reconstruction artefacts. This ensures that the circle centred on the distribution with a radius $r_0$ is free from grating lobes. Therefore, to obtain the full view imaging, the distance between sensors must not fall under $\lambda/2$, i.e. the minimum number of sensor must satisfy

$$N_s > \frac{4\pi R}{\lambda}, \qquad (5)$$

where $N_s$ is the number of sensors.

The same methodology will be employed in order to investigate this further in the context of the shape derivative.
## 2.1 Simulation Setup

### 2.1.1 Configuration

The domain $\Omega$ is a bounded aluminium plate of size $600 \times 580 \times 2.95 \text{ mm}^3$. A defect, represented by a thickness loss, is implemented at the coordinates [360, 280] mm. The defect radius is $r = \frac{3}{2}\lambda$ composed by a flat base of radius $r_{\text{base}} = \frac{1}{2}r$ at the thickness $e_{\text{def}} = 2.55 \text{ mm}$ and a smooth transition zone (described by a Hann half-window) of radius $r_{\text{Hann}} = \frac{1}{2}r$, which varies from $e_{\text{def}}$ to $e_{\text{healthy}}$, see Fig. 1. The used source consists in a 5 cycles tone burst of centre frequency 40 kHz and a Gaussian-shaped normal load at the surface of the plate of standard deviation 2.5 mm. The wavelength of the $A_0$ mode at the center frequency is $\lambda = 25.97 \text{ mm}$. The simulated sequence is a round robin capture of punctual emissions placed on a circle of radius $r_{\text{circle}} = 150 \text{ mm}$. According to the equation (5), the minimal number of sensor to the full view is $N_s \approx 73$ sensors. In a SHM context, this is too many sensors to get a system with an acceptable intrusiveness, so we want to evaluate the ability of the shape derivative to work with under-sampled data.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Thickness loss reference map.</figcaption></figure>

### 2.1.2 Mesh and SEM Configuration

To perform the finite element resolution on the domain, we need to discretize it. We have used a full structured hexahedral mesh, see Fig. 2, with a mesh size such that there are five nodes per wavelength. The length of an edge is $= \frac{\lambda_{\min}}{5\sqrt{2}}$, where $\lambda_{\min} = 21.5 \text{ mm}$ represents the minimal wavelength in the emission signal. The order of the spectral elements is 2 in the two directions of the plane and 3 in the thickness direction. Therefore, there are ten DoFs in the plane and three in the thickness, which allows a good convergence of the finite element.

<figure><img src="image_3.png" /><figcaption>Fig 2. Defect mesh.</figcaption></figure>

Using our formulation, it would be possible to use higher order finite elements to further reduce the number of DoFs. However, to properly reconstruct the domain, the geometrical order of the elements should be also chosen accordingly. In this first study, we restricted ourselves to first order geometrical elements that is straight polyhedral, hence limiting the order of the used finite elements. Future studies will use higher order geometrical elements for a more precise description of the geometry and faster computations.

### 2.1.3 Shape Derivative Results

The shape derivative is performed for several configurations. We started with the optimal configuration, i.e. with the optimal number of sensor $N_s = 73$. Then, we have reduced the number of sensors to evaluate the robustness of the shape derivative to under sampled data.

<footnote n="1">¹ The square root $\sqrt{2}$ is used to ensure that the diagonal of the quadrangular elements is sufficiently sampled in relation to the wavelength.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Shape Derivative reconstruction with 73, 50, 20, 15, 10, 5 sensors. The red circles represent the sensors position and the white circles represents the $r_0$ circle.</figcaption></figure>

In the optimal configuration, see Fig. 3a, we see that the imaging area is free from artefacts. The localisation of the defect is well positioned until ten sensors Fig. 3e, but the quality of the reconstruction is greatly degraded below twenty sensors Fig. 3c.

The white circles on Fig. 1 represent the $r_0$ circle with respect to the number of sensor. For the optimum configuration $r_0 = R$, therefore the white circle is confounded with the sensor circle. We can see that the artefacts appear when the number of sensors is less than the optimum resulting on a noisier background and a less visible defect signature. Below twenty sensors, the defect is outside of the $r_0$ circle. Those results seem to indicate that the relation established under the acoustic assumption also holds for the shape derivative.

The thickness loss for each configuration is plotted in Fig. 4, where we can clearly see that the diminution of the number of sensor deteriorates the result. The relative error on the thickness loss is below 15% up to ten sensors that is very encouraging, because it is one of the most crucial information for the SHM of corrosion. Up to the configuration of fifteen sensors, the thickness loss is slightly overestimated. This is a preferred scenario in an SHM context. Actually, it is preferable to overestimate the residual thickness slightly and carry out maintenance, for example, rather than to underestimate it and risk a failure. However, we also need information that is more global. While the residual thickness provides insight into specific points on the reconstruction map, we also require information on the global image but also on the overall shape of the defect.

In order to enable a quantitative and global comparison of the reconstructions to the reference, the following error calculation is introduced.

$$ \text{err} = \frac{|\Delta_{\text{image}}|}{|e_{\text{ref}}^{\text{healthy}} - \min e_{\text{ref}}|} \times 100, \qquad (6) $$

where $\Delta_{\text{image}} = e - e_{\text{ref}}$ is the difference between the reconstruction map e on Fig. 3 and the reference map $e_{\text{ref}}$ on Fig. 1 and $e_{\text{ref}}^{\text{healthy}}$ is the structure's reference thickness.

The figure 5 illustrates the evolution of the mean of the relative error computed by the equation (6). As expected, the relative error increase as the number of sensors decreases.
It seems to reach an asymptote when the defect is in the $r_0$ area. This can be justified by the fact that at this level, the algorithm has all the data it needs to reconstruct the structure. Therefore, the addition of a new sensor does not result in any further information being produced. When the defect is located outside of the $r_0$ area, it can be observed a "linear" increase of the relative error.

If we now concentrate the calculation of the relative error on the corrosion zone, see Fig. 6, the evolution seems to have a steeper slope than on the entire map, and the degradation of the reconstruction becomes significant from fifteen sensors onwards. Finally, the study on the relative error shows that up to fifteen sensors, the defect is globally well reconstructed.

In a SHM context, some infrastructures need less-intrusive instrumentation and thus keep the number of sensors relatively low. This study reinforces the interest in this method, since with a relatively small number of sensors, twenty in this case, the algorithm continues to converge and the error in terms of residual thickness and global thickness remains low enough. This result must be confirmed in other cases to be sure of its performance.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Plot of the minimum residual thickness with associated relative error to residual thickness of reference.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Plot of the mean of the relative error with the one sigma error bar.</figcaption></figure>

Furthermore, in the configuration with twenty sensors, the defect is located outside of the $r_0$ area. However, the algorithm is still able to reconstruct the defect, demonstrating its robustness against under-sampled data.

### 3. Conclusion

In this paper, the formalism of the shape derivative method adapted to guided elastic waves has been shown. In order to measure the performance of the method, we have chosen to make a study on the

minimal number of sensors to maintain a good shape reconstruction and a precise residual thickness. The results show that the method seems to be robust to under-sampled data. Indeed, even with a relatively small number of sensors compared to the ideal number, the relative error in the defect is kept below 15% for up to ten sensors. Finally, it has been shown that the appearance of grating lobes for this method seems to follow the established law observed for guided wave tomography under an acoustic assumption.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Plot of the mean of the relative error in the defect area with the one sigma error bar.</figcaption></figure>
A recent work established the performance of the method with experimental data [12]. Moreover, the result of the tomography under the acoustic assumption is employed as an initialisation of the shape derivative in order to reduce the number of iterations and the sensitivity of the cost function to local minima. Based on these preliminary results, future works will handle the reconstruction of more complex structural thicknesses.

## Acknowledgments

The author would like to thank Alexandre IMPERIALE (CEA LIST, Paris-Saclay, France) for his help with numerical resolution considerations, particularly with the spectral finite element solver [9].

## References

[1] P. Huthwaite and F. Simonetti, "High-resolution guided wave tomography," *Wave Motion*, vol. 50, no. 5, pp. 79-993, 2013.

[2] F. Simonetti, L. Huang and N. Duric, "On the spatial sampling of wave fields with circular ring apertures," *Journal of Applied Physics*, vol. 101, April 2007.

[3] 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, 2019.

[4] G. Allaire, F. Jouve and A.-M. Toader, "Structural optimization using sensitivity analysis and a level-set method," *Journal of Computational Physics*, vol. 194, p. 363-393, February 2004.

[5] C. Dapogny, "Optimisation de formes," *Gazette des Mathématiciens*, vol. 176, 2023.

[6] J. Hadamard, Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques encastrées., Paris: Imprimerie nationale, 1908.

[7] R.-E. Plessix, "A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.," *Geophysical Journal International*, vol. 167, no. 2, pp. 495-503, 2006.

[8] S. Wright, J. Nocedal and others, Numerical optimization, vol. 35, 1999, p. 7.

[9] 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, p. 964-990, May 2019.

[10] Extende, "STRUCTURAL HEALTH MONITORING WITH CIVA," [Online]. Available: https://www.extende.com/structural-health-monitoring-with-civa. [Accessed 25 04 2024].

[11] O. Mesnil, A. Recoquillay, T. Druet, V. Serey, H. T. Hoang, A. Imperiale and E. Demaldent, "Experimental Validation of Transient Spectral Finite Element Simulation Tools Dedicated to Guided Wave-Based Structural Health Monitoring," *Journal of Nondestructive Evaluation, Diagnostics and Prognostics of Engineering Systems*, vol. 4, April 2021.

[12] A. Thomas, T. Druet, A. Recoquillay and B. Chapuis, "Simulation-Assisted Guided Waves Imaging for Shm : Tomography and Shape Derivative," in *51st Annual Review of Progress in Quantitative Nondestructive Evaluation*, Colorado, 2024.