2-D Elastodynamic Time-Reversal Analysis for Surface Defects on Thin Plate Using Topological Sensitivity

Abstract

In recent years, there has been increasing attention on the development of non-destructive evaluation (NDE) methods using guided waves for long-length materials such as thin plates and pipes. The guided waves are capable of long-distance propagation in thin plates and pipes, and they exhibit properties such as multimodality and dispersion. These characteristics of the guided waves make inspection using guided waves challenging. In this study, we apply a 2-D elastodynamic time-reversal method to detect surface breaking cracks of a thin plate where guided waves are present. The finite element method (FEM) is used to calculate the scattered waves from surface breaking cracks and their corresponding time-reversal waves. We also employ topological sensitivity as an assessment index for detecting surface breaking cracks using the time-reversal method. As numerical examples, we demonstrate guided wave propagation, scattering, and the time-reversal wave propagation obtained by using the FEM. Finally, we present the results of surface breaking crack detection in a thin plate and discuss the validity and effectiveness of the proposed method.

1. Introduction

In recent years, there has been increasing attention on the development of non-destructive evaluation (NDE) methods using guided waves [1] for long-length materials such as thin plates and pipes. The guided waves are capable of long-distance propagation in thin plates and pipes, and they exhibit properties such as multimodality and dispersion [2,3]. These characteristics contribute to the complexity of the guided wave propagation, making it challenging to estimate the position of defects in thin plates and pipes. Numerical simulations, including the finite difference method (FDM) [4], finite element method (FEM) [5], and boundary element method (BEM) [6], play a crucial role in understanding complex ultrasonic guided wave propagation phenomena. In fact, numerous simulations of guided waves using such numerical analysis methods have been conducted [7,8,9,10,11,12,13,14].

On the other hand, the objective of ultrasonic NDE is to detect defects in the material being inspected and to estimate their size, position, and other characteristics. Therefore, several studies have been conducted to reconstruct defects in the material being inspected using guided waves. Gunawan et al. developed a method called the mode exciting method [15] and have determined the relationship between reflection and transmission coefficients, and defects for each guided wave mode. Saitoh et al. [16] have been working on estimating a defect in lap joints using a fast boundary element method called the fast multipole boundary element method (FMM) [17]. Cho et al. [18] investigated the interactions of elastic guided waves with various defects to explore defect characterization possibilities. They employed a hybrid boundary element method (BEM) combined with the elastodynamic boundary integral equation and the Lamb wave normal mode expansion technique. Liu et al. [19] have proposed a two-stage reconstruction algorithm using PWL (Lamb wave tomography), aiming for rapid inspection. In addition, Wang et al. [20] attempted to detect defects within plates using the SH wave by applying the inverse scattering technique based on the Born approximation [21,22]. Additionally, Han et al. [23] used neural networks to estimate defects within plates.

Thus, though many studies have been conducted on inverse analysis using guided waves, the recent proposal of NDE4.0 [24] has drawn attention to NDE utilizing data science and information technology. The time-reversal method can be pointed out as an inverse analysis technique that is relatively close to the concept of NDE4.0. The time-reversal method was first developed by Fink [25]. The time-reversal method is an inverse analysis technique that fully exploits the reversibility and reciprocity of wave propagation. The scattered waves from a defect observed at the measurement points are time-reversed and then virtually re-introduced into the test material [26,27,28]. During this process, the propagation of the time-reversed waves is replicated using numerical simulations such as FEM and BEM. These operations are similar to the concept of digital twins, making them well-matched with the NDE4.0 concept. In general, the time-reversal method determines the position of defects as the convergence point of the time-reversed wave. However, when applying this time-reversal method to defect detection in thin plates, the estimation of the convergence position of the time-reversed waves becomes extremely challenging due to their repeated reflections by the top and bottom surfaces of the plate. This problem has been addressed using topological sensitivity, which is widely used in the field of structural optimization and topology optimization [29,30]. In fact, successful inverse analysis of surface cracks in a thin plate using SH waves has been achieved [31].

Thus, in this paper, the previous work [31] by the author using the SH wave is extended to a 2-D elastic wave problem considering in-plane wave propagation. In the following, we first describe the problem to be solved. However, instead of using scattered waves obtained from actual measurement experiments, we use scattered waves calculated by the FEM. Consequently, the problem to be addressed involves FEM simulations for both scattered waves and time-reversal waves. Next, we briefly explain the FEM formulation for 2-D elastodynamics and the concept of topological sensitivity. Then, we will present the results of estimating the position of defects using the time-reversal method with the topological sensitivity. By comparing and examining these results with the actual positions of defects, we can evaluate the applicability and performance of the proposed method.

2. Problem Statement

In this section, we discuss the analysis model, conditions, and the FEM for Lamb wave propagation and scattering in the plate.

2.1. Analysis Model and Conditions

Let us consider 2-D guided wave propagation and scattering by surface breaking cracks in an infinite and homogeneous plate with a thickness of h, as illustrated in Figure 1. Incident waves are excited from source points and propagate in the infinite plate, guided waves are excited, and some of them are scattered by surface braking cracks and the rest are transmitted in a thin plate. Source points correspond to ultrasonic transducers in experiments and the incident waves are represented by providing an appropriate function to the displacement at source points. The scattered waves from surface breaking cracks are received at an observation area located to the left of the surface breaking cracks, as shown in Figure 1. The problem is determining the position and length of surface breaking cracks using scattered waves obtained at observation points and their time-reversed waves. However, surface breaking cracks are assumed to be straight and perpendicular to the surface of the plate.

As mentioned in the previous section, for the analysis model of an infinitely long plate, as shown in Figure 1, we perform FEM calculations for wave propagation inside the plate, scattering due to surface breaking cracks, and the computation of time-reversal waves with the final goal of defect detection.

2.2. FEM Formulation for Lamb Wave Propagation and Scattering

In this section, we briefly discuss the FEM used for calculating the waves propagating inside the plate. However, there are many references available regarding the formulation of the FEM. Therefore, in this section, we will only describe the essential aspects of the formulation.

2.2.1. Governing Equation and Selection of Numerical Method

We modelize the infinite plate as an isotropic homogeneous elastic solid with Lamé constants $\lambda$ and $\mu$, and the mass density $\rho$. The in-plane displacement, $u_i(\mathit{x},t)(i=1,2)$, in a thin plate, $\mathrm{\Omega }$, satisfies the equation of motion at time t and the position $\mathit{x}$ as follows:

$$\begin{array}{cc}\hfill \rho {\ddot{u}}_i& =\mu u_{i,jj}+(\lambda +\mu )u_{j,ij}(\mathit{x}\in \mathrm{\Omega },0<t)\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & u_i(\mathit{x},0)={\dot{u}}_i(\mathit{x},0)=0(\mathit{x}\in \mathrm{\Omega })\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & t_i(\mathit{x},t)=0(\mathit{x}\in \partial \mathrm{\Omega },0<t)\hfill \end{array}$$

(3)

where $\dot{\left(\right)}$ is the partial derivative with respect to time t and $t_i(\mathit{x},t)$ is the traction component that corresponds to the displacement $u_i(\mathit{x},t)$. In addition, ${\left(\right)}_{,i}$ indicates the partial derivative with respect to the $x_i$, $\partial /\partial x_i$, and $\partial \mathrm{\Omega }$ shows the boundary of the domain $\mathrm{\Omega }$. Equation (3) represents the traction-free condition on the top and bottom surfaces of the thin plate.

The objective of this research is to numerically reconstruct surface breaking cracks using the received scattered waves $u_i^{\mathrm{sc}}(\mathit{x},t)$ at observation points. In general, the FEM and the BEM are often used for performing time-domain wave analysis. The classical time-domain BEM can be time-consuming and may present numerical instability for small time increments [6]. As a method to improve the numerical stability of the classical time-domain BEM, the convolution quadrature time-domain BEM [32,33] has been developed. However, it requires the implementation of somewhat complex techniques such as the FMM [17] for the acceleration [34,35]. Therefore, in this study, the time-domain FEM is utilized to simulate the wave propagation and scattering by surface breaking cracks and the time-reversal wave propagation. In particular, the pixel-based modeling is used for the space discretization of the FEM. This calculation step for obtaining scattered waves at observation points is called “forward analysis” in this paper.

2.2.2. Space and Time Discretization Using Pixel-Based Modeling

The pixel-based FEM [36] in 2-D elastodynamics is utilized to obtain the wave field data inside the plate and the scattered wave data at the observation area, as shown in Figure 1. Typically, FEM meshing tools such as Gmsh [37] are used for finite element modeling. However, if image data for the analysis model can be prepared, considering that images typically consist of a large number of pixels, it is possible to easily create a finite element model corresponding to the analysis model by associating one pixel with one finite element. This method, known as pixel-based modeling, is employed to simplify the creation of FEM models for analysis.

The Galerkin method is applied to the spatial discretization of Equation (1). By multiplying the shape function, $N_{\alpha }(\alpha =1,\cdots ,4)$ [5], as the weight function by Equation (1), and integrating the resulting equation over the analysis domain $\mathrm{\Omega }$ to be solved, we can obtain the following equation:

$$\begin{array}{c}\hfill \sum _{e=1}^m{\int }_{v^e}{N}_{\alpha }(\mu u_{i,jj}+(\lambda +\mu )u_{j,ij}-\rho {\ddot{u}}_i)d{v}^e=0\end{array}$$

(4)

where m is the number of finite elements and $v_e(e=1,\dots ,m)$ is each finite element region. The shape function, $N_{\alpha }$, is given by $N_{\alpha }=1/4(1-{\xi }_1)(1-{\xi }_2)$ with the local coordinate system of ${\xi }_1$ and ${\xi }_2$ for the case of considering the use of pixel elements in this study [5]. Here, Green’s theorem is written as follows:

$$\begin{array}{c}\hfill {\int }_{v^e}{N}_{\alpha }{u}_{i,jj}d{v}^e={\int }_{S^e}{N}_{\alpha }{u}_{i,j}{n}_{j}d{S}^e-{\int }_{v^e}{N}_{\alpha ,j}{u}_{i,j}d{v}^e\end{array}$$

(5)

$$\begin{array}{c}\hfill {\int }_{v^e}{N}_{\alpha }{u}_{j,ij}d{v}^e={\int }_{S^e}{N}_{\alpha }{u}_{j,j}{n}_{i}d{S}^e-{\int }_{v^e}{N}_{\alpha ,i}{u}_{j,j}d{v}^e\end{array}$$

(6)

where $n_k$ is the $x_k$ direction component of the outward unit normal vector $\mathit{n}$ on a boundary of a finite element. In addition, $S^e$ represents the boundary of the finite element domain $v^e$. By substituting Equations (5) and (6) into Equation (4), we have:

$$\begin{array}{cc}\hfill & \sum _{e=1}^m\sum _{\alpha ,\beta =1}^4[\mu {\int }_{v^e}{N}_{\alpha ,j}{N}_{\beta ,j}d{v}^e{u}_{i\beta }^e+\lambda {\int }_{v^e}{N}_{\alpha ,i}{N}_{\beta ,j}d{v}^e{u}_{j\beta }^e+\mu {\int }_{v^e}{N}_{\alpha ,j}{N}_{\beta ,i}d{v}^e{u}_{j\beta }^e]\hfill \\ \hfill & +\sum _{e=1}^m\sum _{\beta =1}^4\left[\rho {\int }_{v^e}{N}_{\alpha }{N}_{\beta }d{v}^e{\ddot{u}}_{i\beta }^e\right]-\sum _{e=1}^m\sum _{\beta =1}^4\left[{\int }_{S^e}{N}_{\alpha }{N}_{\beta }d{S}^e{t}_{i\beta }^e\right]=0\hfill \end{array}$$

(7)

where $u_{i\beta }^e$ denotes the displacement $u_i(\mathit{x},t)$ at the $\beta$-th node of the finite element $v^e$. The material constants, $\rho$, $\lambda$, and $\mu$, are assumed to be constant within each finite element. Therefore, these material constants have been taken outside of the integral. Generally, in the FEM, the first term of Equation (7) represents the stiffness matrix, the second term represents the mass matrix, and the third term shows the traction vector components. Therefore, Equation (7) can be reduced to the following well known matrix representation of the finite element equation:

$$\begin{array}{c}\hfill \left[K\right]\left\{u_i\right\}+\left[M\right]\left\{{\ddot{u}}_i\right\}-\left\{T_i\right\}=0\end{array}$$

(8)

where $\left[K\right]$ is the global stiffness matrix, $\left\{u_i\right\}$ is the node displacement vector, $\left[M\right]$ is the global mass matrix, and $\left\{T_i\right\}$ is the traction vector. The second term in Equation (8) has the partial derivative with respect to time t. In this research, this term with respect to time t is discretized by the central difference scheme as follows:

$$\begin{array}{c}\hfill \left\{{\ddot{u}}_i\right\}\simeq \frac{{\left\{u_i\right\}}_{n+1}-2{\left\{u_i\right\}}_n+{\left\{u_i\right\}}_{n-1}}{{\left(\Delta t\right)}^2}\end{array}$$

(9)

where $\Delta t$ is the time increment. Assuming the matrix $\left[M\right]$ is a diagonal matrix, Equation (8) is explicitly approximated using the following equation:

$$\begin{array}{c}\hfill {\left\{u_i\right\}}_{n+1}=-\left[{(\Delta t)}^2{\left[M\right]}^{-1}\left[K\right]+2\left[E\right]\right]{\left\{u_i\right\}}_n-{\left\{u_i\right\}}_{n-1}+{(\Delta t)}^2{\left[M\right]}^{-1}\left\{T_i\right\}\end{array}$$

(10)

where $\left[E\right]$ is the identity matrix. Equation (10) at the $(n+1)$-th time step can be solved explicitly with the solutions at the previous n and $(n-1)$-th time steps. Thus, Equation (10) can be solved step-by-step starting from the initial and boundary conditions on the plate surface at $n=0$. The computational order is $O\left(mn\right)$.

By using the FEM described in this section, it is possible to determine the displacement field $u_i(\mathit{x},t)$ inside the plate, as shown in Figure 1, as well as the scattered wave $u_i^{\mathrm{sc}}(\mathit{x},t)$ by surface breaking cracks in the observation area.

As previously mentioned, in this study, pixel elements are used as finite elements. However, in general, it is difficult to represent precise cracks using pixel elements. Therefore, in this study, surface breaking cracks are treated as slits with a finite thickness.

3. Time-Reversal and Topological Sensitivity

In this section, the time-reversal method and the topological sensitivity for 2-D elastodynamics are briefly described for the surface breaking crack detection in a plate.

3.1. Time-Reversal Method

We consider detecting the surface breaking cracks in the domain $\mathrm{\Omega }$ using the time-reversal method described in this section. The waves scattered by surface breaking cracks, observed at the receiver points in the observation area, as illustrated in Figure 1, are time-reversed. These time-reversed waves are then sent back into the inside of the plate $\mathrm{\Omega }$ as incident waves from the receiver points and propagate to both the left and right sides of the plate. The propagation of this time-reversed wave is analyzed using the FEM in the same manner as the forward analysis discussed in Section 2.2.

The boundary value problem for the time-reversal wave fields, corresponding to the forward analysis defined in Equations (1)–(3), is formulated as follows:

$$\begin{array}{cc}\hfill \rho {\ddot{u}}_i^{\mathrm{tr}}& =\mu u_{i,jj}^{\mathrm{tr}}+(\lambda +\mu )u_{j,ij}^{\mathrm{tr}}(\mathit{x}\in \mathrm{\Omega },0<t<T)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill u_i^{\mathrm{tr}}(\mathit{x},0)& ={\dot{u}}_i^{\mathrm{tr}}(\mathit{x},0)=0(\mathit{x}\in \mathrm{\Omega })\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill u_i^{\mathrm{tr}}({\mathit{x}}^m,t)& =u_i^{\mathrm{sc}}({\mathit{x}}^m,T-t)(0<t<T)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill t_i^{\mathrm{tr}}(\mathit{x},t)& =0(\mathit{x}\backslash {\mathit{x}}^m\in \partial \mathrm{\Omega },0<t)\hfill \end{array}$$

(14)

where $u_i^{\mathrm{tr}}(\mathit{x},t)$ is the time-reversal wave field and $t_i^{\mathrm{tr}}(\mathit{x},t)$ represents the traction fields corresponding to $u_i^{\mathrm{tr}}(\mathit{x},t)$. In addition, T is the measurement time for the forward FEM analysis, and ${\mathit{x}}^m(m=1,\dots ,M)$ represents the observation points. M denotes the total number of observation points. The boundary condition, given in Equation (13), is derived from the time-reversal of the scattered wave $u_i^{\mathrm{sc}}({\mathit{x}}^m,t)$ at the receiver points ${\mathit{x}}^m(m=1,\dots ,M)$, which can be obtained by solving the forward analysis discussed in Section 2.2. $u_i^{\mathrm{tr}}({\mathit{x}}^m,t)$ in Equation (13) plays the role of the incident wave for this time-reversal analysis. This time-reversal problem defined in Equations (11)–(14) is also solved by the FEM described in Section 2.2. The time-reversal waves transmitted from the receiver points, ${\mathit{x}}^m$, travel inside the plate, as shown in Figure 2, by repeating the reflection on the upper and lower surface of the plate. Thus, the time-reversal wave analysis can be simplified to the problem described by Equations (11)–(14).

If the wave propagation model is an infinite space instead of a thin plate, we would be able to visually recognize the time-reversed waves converging back to the original scattering sources, which are the surface breaking cracks. At that time, the convergence location of the time-reversed waves can be considered as the position of the surface breaking cracks. However, unfortunately, in this problem involving guided waves, the convergence location cannot be visually determined because the time-reversed waves propagate as guided waves and their multiple modes are excited, resulting in complex wave fields in the time domain. To address this challenge, topological sensitivity is utilized to find the focusing location of the time-reversed waves.

3.2. Topological Sensitivity

We will now provide a brief overview of topological sensitivity. The topological sensitivity, denoted as $\mathcal{T}\left(\mathit{x}\right)$, is a metric employed in topology optimization [30] to measure design sensitivity. The topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ is defined in the context of two scenarios. We firstly consider a defect-free domain ${\mathrm{\Omega }}^0$, as depicted in Figure 3a. The introduction of an infinitesimally small circular void with radius $\epsilon$ at a point $\mathit{x}$ in this domain ${\mathrm{\Omega }}^0$ results in a variation in the objective function. This newly introduced circular void is illustrated in Figure 3b, where ${\Gamma }^{\epsilon }$ represents its boundary. Indeed, when an infinitesimal circular void is introduced, we can define the region as ${\mathrm{\Omega }}^{\epsilon }$, and the corresponding objective function as $J\left({\mathrm{\Omega }}^{\epsilon }\right)$. Given this scenario, the topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ can be expressed by normalizing with the voids area, given by $\pi {\epsilon }^2$, as follows:

$$\begin{array}{c}\hfill \mathcal{T}\left(\mathit{x}\right)=\underset{\epsilon \to 0}{\lim }\frac{J\left({\mathrm{\Omega }}^{\epsilon }\right)-J\left({\mathrm{\Omega }}^0\right)}{\pi {\epsilon }^2}\end{array}$$

(15)

where ${\mathrm{\Omega }}^0$ and ${\mathrm{\Omega }}^{ϵ}$ show the reference domain without any defects and with an infinitesimal circular void, respectively. Here, the objective function is defined as the difference between two received wave data at observation point ${\mathit{x}}^m(m=1,\dots ,M)$ on the observation area $S^{\mathrm{obs}}$. The first received wave data, $u_i^0({\mathit{x}}^m,t)$, correspond to the scenario without surface breaking cracks, whereas the other wave data, $u_i^{\mathrm{obs}}({\mathit{x}}^m,t)$, reflect the presence of surface breaking cracks at the actual location. The latter case corresponds to the actual measured waveform data in NDE measurements. Thus, the objective function can be expressed by the following equation:

$$\begin{array}{cc}\hfill & J\left({\mathrm{\Omega }}^0\right)={\int }_0^T{\int }_{S^{\mathrm{obs}}}^{}\phi (u_i^0({\mathit{x}}^m,t),{\mathit{x}}^m,t)d{S}_{x}dt\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \phi (u_i^0({\mathit{x}}^m,t),{\mathit{x}}^m,t)=\frac{1}{2}\mid u_i^0({\mathit{x}}^m,t)-u_i^{\mathrm{obs}}({\mathit{x}}^m,t)|\hfill \end{array}$$

(17)

To evaluate the objective function $J\left({\mathrm{\Omega }}^{\epsilon }\right)$, the solution $u_i^{\epsilon }(\mathit{x},t)$ for the case in which the domain $\mathrm{\Omega }$ has an infinitesimal circular void with radius $ϵ$ is solved by the following initial boundary value problems:

$$\begin{array}{cc}\hfill \rho {\ddot{u}}_i^{\epsilon }(\mathit{x},t)& =\mu u_{i,jj}^{\epsilon }(\mathit{x},t)+(\lambda +\mu )u_{j,ji}^{\epsilon }(\mathit{x},t)(\mathit{x}\in {\mathrm{\Omega }}^{\epsilon },0<t)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill u_i^{\epsilon }(\mathit{x},0)& ={\dot{u}}_i^{\epsilon }(\mathit{x},0)=0(\mathit{x}\in {\mathrm{\Omega }}^{\epsilon })\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill t_i^{\epsilon }(\mathit{x},t)& =0(\mathit{x}\in {\Gamma }^{\epsilon },0<t)\hfill \end{array}$$

(20)

where $t_i^{\epsilon }(\mathit{x},t)$ is the traction component corresponds to the displacement $u_i^{\epsilon }(\mathit{x},t)$. Solution $u_i^{\epsilon }(\mathit{x},t)$ for the initial boundary value problems defined in Equations (18)–(20) can be divided into $u_i^0(\mathit{x},t)$ and the other solutions, ${\tilde{u}}_i^{\epsilon }(\mathit{x},t)$, as follows:

$$\begin{array}{c}\hfill u_i^{\epsilon }(\mathit{x},t)=u_i^0(\mathit{x},t)+{\tilde{u}}_i^{\epsilon }(\mathit{x},t)\end{array}$$

(21)

Note that the physical meaning of $u_i^0(\mathit{x},t)$ is that $u_i^0(\mathit{x},t)$ is the free field without any surface breaking cracks, as mentioned in below Equation (15). The integrand $\phi$ of Equation (16) can be expanded as follows:

$$\begin{array}{c}\hfill \phi (u_i^{\epsilon }(\mathit{x},t),\mathit{x},t)=\phi (u_i^0(\mathit{x},t),\mathit{x},t)+\frac{\partial \phi }{\partial u_i^0}(u_i^0(\mathit{x},t),\mathit{x},t){\tilde{u}}_i^{ϵ}(\mathit{x},t)+o\left(\right|{\tilde{u}}_i^{ϵ}\left|\right)\end{array}$$

(22)

By taking the integration with respect to time t and observation area $S^{\mathrm{obs}}$, Equation (22) can be written as follows:

$$\begin{array}{c}\hfill J\left({\mathrm{\Omega }}^{\epsilon }\right)=J\left({\mathrm{\Omega }}^0\right)+{\int }_0^T{\int }_{S^{\mathrm{obs}}}\frac{\partial \phi }{\partial u_i^0}(u_i^0(\mathit{x},t),\mathit{x},t){\tilde{u}}_i^{\epsilon }(\mathit{x},t)d{S}_{x}dt\end{array}$$

(23)

By substituting Equation (23) into Equation (15), the topological sensitivity, $\mathcal{T}\left(\mathit{x}\right)$, at point $\mathit{x}$ is derived as:

$$\begin{array}{c}\hfill \mathcal{T}\left(\mathit{x}\right)=\underset{\epsilon \to 0}{lim}\frac{1}{\pi {ϵ}^2}{\int }_0^T{\int }_{S^{\mathrm{obs}}}\frac{\partial \phi }{\partial u_i^0}(u_i^0(\mathit{x},t),\mathit{x},t){\tilde{u}}_i^{\epsilon }(\mathit{x},t)d{S}_{x}dt\end{array}$$

(24)

Considering the time-reversal analysis defined in Equations (11)–(14), the reciprocity theorem, and Gauss’s divergence theorem, the topological sensitivity, $\mathcal{T}\left(\mathit{x}\right)$, can be obtained as follows [38]:

$$\begin{array}{cc}\hfill & \mathcal{T}\left(\mathit{x}\right)=({\widehat{\sigma }}_{ij}^0\ast \left(A_{ijkl}{\sigma }_{kl}^0\right)+\rho {\dot{u}}_i^0\ast {\dot{\widehat{u}}}_i^0)(\mathit{x},T)\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & A_{ijkl}=\frac{1-\nu }{\mu }\left[({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})-\frac{1}{2(1+\nu )}{\delta }_{ij}{\delta }_{kl}\right]\hfill \end{array}$$

(26)

where ${\delta }_{ik}$ is the Kronecker delta and ${\sigma }_{ij}^0$ is the stress field corresponding to the wave field $u_i^0$. Moreover, $\nu$ is the Poisson’s ratio and ∗ represents the convolution integral.

In Equation (25), the displacement, $u_i^0$, and its corresponding particle velocity, ${\dot{u}}_i^0$, can be evaluated by the forward FEM analysis for the domain ${\mathrm{\Omega }}^0$ without any surface breaking cracks. However, the displacement ${\widehat{u}}_i^0$, its corresponding particle velocity ${\dot{\widehat{u}}}_i^0$, and the stress field ${\widehat{\sigma }}_{ij}^0$ can be calculated using the time-reversal analysis introduced in Section 4.1. The displacement field ${\widehat{u}}_i^0$ is the solution to the adjoint problem defined for the forward analysis described in Section 2.2. Hence, the topological sensitivity, $\mathcal{T}\left(\mathit{x}\right)$, in domain $\mathrm{\Omega }$ can be determined by executing both forward and time-reversal analyses, as defined by Equations (1)–(3) and Equations (11)–(14), respectively. Topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$, defined in Equation (15), can be considered as an index that shows how much the received waveform is affected when a small defect exists at a certain place inside the area. In practical application to the ultrasonic NDE, an actual surface breaking crack can be considered as a superposition of the small infinitesimal voids defined in Figure 3. The topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ has a large negative value at the defect position, as shown in Equation (15).

Thus, if the topological sensitivity, $\mathcal{T}\left(\mathit{x}\right)$, inside the inspection area of a thin plate can be obtained, there is a possibility that the position and shape of the defect can be reconstructed.

4. Numerical Example

In this section, some numerical examples are presented. The thickness of the plate is $h=4$ [mm]. Furthermore, the horizontal analysis region of the plate is 160 [mm], and beyond that, absorption boundaries are set on both sides to assume the infinite plate. The time increment $\Delta t$ and the finite element width $\Delta x$ for the FEM analysis are set as $\Delta t=10.0$ [ns] and $\Delta x=0.1$ [mm], respectively. The density $\rho$ of the plate is $\rho =2700$ [kg/m${}^3$]. The velocities of the longitudinal wave (P-wave) $c_L$ and transverse wave (S-wave) $c_T$ are given by $c_L=6350$ [m/s] and $c_T=3130$ [m/s], respectively. The incident wave is excited by applying the following Ricker wave at each incident point ${\mathit{x}}^i(i=1,\dots ,M^{\mathrm{in}})$:

$$\begin{array}{cc}\hfill u_1^{\mathrm{in}}({\mathit{x}}^i,t)& =\frac{\sqrt{\pi }}{2}(\alpha -0.5)exp(-\alpha )\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \alpha & ={\left(\frac{\pi (t-t_s)}{t_p}\right)}^2\hfill \end{array}$$

(28)

In the following numerical examples, the incident wave period $t_p$ and the time parameter $t_s$ for adjusting the peak of the incident wave are given by $t_p=2.0\times {10}^{-6}$ [s] and $t_s=2.0\times {10}^{-6}$ [s], respectively. As with the BEM, the FEM cannot accurately represent cracks with zero thickness. Therefore, in this study, each crack is modeled as a 2-[mm]-long slit. In addition, the boundary condition of the plate surface is considered as the traction free, $t_i(\mathit{x},t)=0$.

4.1. Reconstruction of Single Surface Breaking Crack from Upper Side of a Plate

As the initial numerical example, the reconstitution of a single surface breaking crack on the top of the plate is conducted as shown in Figure 4. The incident points where the Ricker wave in Equation (27) is applied are set at 7 points ($M^{\mathrm{in}}=7$) with a 1.0 [mm] interval, as shown in Figure 4. On the other hand, the observation points where the scattered waves are received are set at 41 points ($M=41$) with a 1.0 [mm] interval. The surface breaking crack is located at $x_1=130$ [mm], as shown in Figure 4.

4.1.1. Forward Analysis for Single Surface Breaking Crack

First, the scattering analysis of the incident wave due to a single surface crack, as shown in Figure 4, is conducted using the FEM. Figure 5 shows the results for elastic wave scattering by a single surface breaking crack in Figure 4 obtained by the FEM analysis at several time-steps. In each of Figure 5, the actual single surface breaking crack is indicated by a white line.

The elastic waves excited at the incident points ${\mathit{x}}^i(i=1,\dots ,7)$ propagate by undergoing reflections on the top and bottom surfaces. Due to this fact, the symmetric and antisymmetric plate waves are generated, as shown in Figure 5a. The first wave of symmetric mode arrives at the surface breaking crack in Figure 5b. Moreover, due to the surface breaking crack on the top surface of the plate, it can be observed that some of the excited waves are transmitted, whereas others are scattered, and propagate within the interior of the plate, as shown in Figure 5c,d. The scattered waves generated by the surface breaking crack are received at 41 observation points, as shown in Figure 4, and they are utilized for time-reversal analysis discussed in the next section.

4.1.2. Time-Reversal Analysis for Single Surface Breaking Crack

Time-reversal analysis is performed using the results obtained by the forward analysis in the previous section. The received waves at the observation points ${\mathit{x}}^m(m=1,\dots ,41)$ are time-reversed and sent back into the plate again. In the case of time-reversal analysis, note that the crack existence is not considered. Therefore, we solely focus on propagating the time-reversal wave within the interior of the plate.

Figure 6 shows the results for the time-reversal FEM analysis at several time-steps. For reference, the approximate location of the actual crack is also shown in Figure 6. As shown in Figure 6, we can observe the propagation of the time-reversal wave within the interior of the plate. However, in the calculation of the time-reversal wave, it is difficult to clearly distinguish between the symmetric and antisymmetric modes and confirm the convergence point of time-reversal waves. The scattered waves are propagated in the reverse direction, making it difficult to distinguish between each mode. In general, if the time-reversal method is applied to a defect in an infinite domain, it is easy to visually confirm that the time-reversed waves converge to the scattering source (defect). Therefore, in this study, we attempt the surface breaking crack reconstruction by conducting topological sensitivity calculations for the scenario involving plate waves that undergo repeated reflections on the upper and lower surfaces.

4.1.3. Reconstruction of Single Surface Breaking Crack

Using the results from the previous forward analysis and time-reversal analysis, we perform topological sensitivity calculations to reconstruct the single surface breaking crack from the top surface of the plate.

Figure 7 shows the obtained topological sensitivity, $\mathcal{T}\left(\mathit{x}\right)$, using the forward and time-reversal analysis results, as derived from Equation (25). The visualization range for the topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ is within the visualization area of Figure 4 (i.e., $x_1$ = 120 [mm] to $x_1$ = 140 [mm]). The actual single surface breaking crack is indicated by the blue line in Figure 7.

From Figure 7, we can see that the topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ exhibits significant negative values near the position of the surface breaking crack. Indeed, it can be observed that the topological sensitivity values are higher on the incident side than at the actual position of the crack. The reconstruction result shows that the single surface breaking crack is well reconstructed, though it is slightly shorter in length compared to the actual crack length.

4.2. Reconstruction of Two Surface Breaking Cracks from Upper and Lower Side of a Plate

Next, the reconstruction of two surface breaking cracks from the upper and lower surface of the plate are considered. Two surface breaking cracks are located at $x_1=110$ [mm] (lower surface) and $x_1=130$ [mm] (upper surface), as shown in Figure 8. The positions of the source points and observation points remain the same as those in the Section 4.1. Under these conditions, the following forward analysis, time-reversal analysis, and the topological sensitivity analysis are conducted.

4.2.1. Forward Analysis for Two Surface Breaking Cracks

Figure 9 shows the guided wave propagation and the scattering by two surface breaking cracks in Figure 8. As seen in Figure 9a, the guided waves of the symmetric and antisymmetric modes can be confirmed. As shown in Figure 9b, the excited guided waves are first scattered by the surface breaking crack from the bottom surface, and the displacement indicates large values just before the crack plane. Then, scattered waves are generated by the interaction between the incident waves and two surface breaking cracks, as shown in Figure 9c. After that, the scattered waves generated by two crack planes propagate towards the left side of the plate, as shown in Figure 9d. Similar to Section 4.1.1 for the single surface breaking crack, these scattered waves are received at 41 observation points in Figure 8.

4.2.2. Time-Reversal Analysis for Two Surface Breaking Cracks

Next, the time-reversal analysis corresponding to the forward analysis presented in Section 4.1.2 is carried out. Figure 10 shows the results for the time-reversal analysis in which no two surface breaking cracks are considered at several time-steps. As seen in Figure 10, the time-reversal wave propagates through the infinite plate as time progresses. However, as mentioned in the previous time-reversal results shown in Figure 6, the convergence points of the time-reversal waves cannot be determined because of the complicated wave fields inside the thin plate.

4.2.3. Reconstruction of Two Surface Breaking Cracks

Here, the reconstruction result for two surface breaking cracks on the top and bottom surfaces are shown. Figure 11 shows the topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$, calculated using Equation (25), around both surface breaking cracks. The topological sensitivity $\mathcal{T}\left(\mathit{x}\right)$ is plotted for 100 [mm] $\le x_1\le$ 140 [mm] in Figure 8. The actual two surface breaking cracks are indicated by the blue line in Figure 11.

As well as the results in Figure 7, the reconstruction in this case also shows that the incident side with respect to the crack plane can be successfully reconstructed, as shown in Figure 11. Furthermore, it can be observed that the upper surface breaking crack, which is located on the same plane as the incident and observation points, shows a more accurate indication of the defect size (crack depth) due to the influence of surface waves compared to the lower surface crack.

Thus, we demonstrated that our proposed method has a potential to effectively detect surface breaking cracks even when considering guided wave propagation.

5. Discussion

A comprehensive discussion on the proposed method is provided in this section based on the reconstruction results of surface breaking cracks shown in Figure 7 and Figure 11. In the proposed method, as mentioned in Section 3, a time-reversal analysis must be conducted to detect defects. The time-reversal analysis is generally performed using numerical techniques such as FEM or BEM. Therefore, the computational time for defect detection depends on their respective computational order. From this point of view, the proposed method typically takes more time to complete the inverse analysis compared to other inverse analysis methods often used in NDE, such as SAFT (synthetic aperture focusing technique) [39] and inverse scattering techniques [40]. However, the time-reversal analysis, as indicated by Equations (11)–(14), takes into account the boundary conditions of the area under inspection. Therefore, if an analysis model can be prepared, it is believed that defects inside a domain with complex boundaries can be detected with high precision. This method can be applied not only to plates, but also to various complex shapes, such as the curved sections of pipes and layered objects. The computational time of direct calculations, such as FEM and BEM, is expected to decrease with the advancement of computers in the future. Therefore, it is believed that the benefit of being able to handle complex problems will be greater than the computational time concerns.

6. Conclusions

We performed defect detection of surface breaking cracks in a thin plate using the FEM. By employing topological sensitivity as a defect detection index in the framework of the 2D elastodynamics time-reversal method, we confirmed the capability of defect detection for surface breaking cracks in a thin plate. The numerical results demonstrated the effectiveness of the proposed method in cases with up to two surface breaking cracks. In the future, we plan to further improve the defect detection accuracy, which means determining the length of surface breaking cracks and pinpointing their positions, and to investigate the applicability of this method to anisotropic materials such as CFRP (carbon fiber reinforced plastics). Furthermore, we will extend this approach for 2-D to 3-D cases.