SHM via topological derivative

We consider a rectangular metallic plate of size $(2{L}_{1})\,\times (2{L}_{2})\times d$. As anticipated in section 1, we shall only consider symmetric Lamb waves [45]. Assuming that the waves wavelength is large compared to the plate thickness d and that this thickness is small compared to both 2L₁ and 2L₂, the elastodynamic state admits a plane stress assumption, which allows for a two-dimensional approximation. Using a Cartesian coordinate system with the origin at the center of the plate and the x and y axes parallel to the plate sides, the plate is the rectangle $D:-{L}_{1}\lt x\lt {L}_{1},-{L}_{2}\lt y\lt {L}_{2}$. With the usual notation, the in-plane displacement vector $\overrightarrow{u}=(u,v)$, assumed to be conveniently small, obeys the following linear elastodynamics equation

$$\begin{eqnarray}&&\rho \displaystyle \frac{{\partial }^{2}\overrightarrow{u}}{\partial {t}^{2}}-{\rm{\nabla }}(\lambda {\rm{\nabla }}\cdot \overrightarrow{u})-{\rm{\nabla }}\cdot [\mu ({\rm{\nabla }}\overrightarrow{u}+{\rm{\nabla }}{\overrightarrow{u}}^{\top })]=\overrightarrow{f}(x,y,t)\end{eqnarray} \tag{ 1 }$$

in the domain D, with boundary conditions

$$\begin{eqnarray}&&\lambda ({\rm{\nabla }}\cdot \overrightarrow{u})\overrightarrow{n}+\mu ({\rm{\nabla }}\overrightarrow{u}+{\rm{\nabla }}{\overrightarrow{u}}^{\top })\cdot \overrightarrow{n}=\overrightarrow{0}\quad \mathrm{at}\ \\ &&\quad x=\pm {L}_{1},| y| \lt {L}_{2}\ \mathrm{and}\ y=\pm {L}_{2},| x| \lt {L}_{1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\overrightarrow{u}=\overrightarrow{0}\quad \mathrm{at}\ (x,y)=(\pm {L}_{1},-{L}_{2}).\end{eqnarray} \tag{ 3 }$$

These impose pure reflection at the lateral walls, in (2), and eliminate rigid-body displacements, in (3). In this formulation, ρ is the density, t denotes time, λ and μ are the Lamé coefficients, ${}^{\top }$ denotes the transpose, the forcing term $\overrightarrow{f}$ depends on the actuators configuration, and $\overrightarrow{n}$ is the outward unit normal to the boundary of D. Note that equations (1) and (2) can also be written as $\rho {\partial }^{2}\overrightarrow{u}/\partial {t}^{2}-{\rm{\nabla }}\cdot \bar{\bar{\sigma }}=\overrightarrow{f}$ and $\bar{\bar{\sigma }}\cdot \overrightarrow{n}=\overrightarrow{0}$, respectively, where the stress tensor $\bar{\bar{\sigma }}$ is defined as

$$\begin{eqnarray}&&\bar{\bar{\sigma }}=\lambda ({\rm{\nabla }}\cdot \overrightarrow{u}){\bar{\bar{I}}}_{2}\,+\,2\mu \bar{\bar{\varepsilon }}.\end{eqnarray} \tag{ 4 }$$

Here, ${\bar{\bar{I}}}_{2}$ is the 2 × 2 unit tensor and $\bar{\bar{\varepsilon }}\equiv ({\rm{\nabla }}\overrightarrow{u}+{\rm{\nabla }}{\overrightarrow{u}}^{\top })/2$ is the strain tensor, whose components are given by

$$\begin{eqnarray}{\varepsilon }_{11} & = & \displaystyle \frac{\partial u}{\partial x},\quad {\varepsilon }_{12}={\varepsilon }_{21}=\displaystyle \frac{1}{2}\left[\displaystyle \frac{\partial u}{\partial y}+\displaystyle \frac{\partial v}{\partial x}\right],\\ {\varepsilon }_{22} & = & \displaystyle \frac{\partial v}{\partial y}.\end{eqnarray} \tag{ 5 }$$

Since the solution to the problem (1)–(3) is linear on the forcing term $\overrightarrow{f}$ (though not on the material properties ρ, λ, and μ!), for multi-frequency excitation produced by M forcing devices, of the form

$$\begin{eqnarray}&&\overrightarrow{f}(x,y,t)=\displaystyle \sum _{k=1}^{K}\displaystyle \sum _{m=1}^{M}{\overrightarrow{F}}_{{km}}(x,y)\cos ({\omega }_{k}t+{\phi }_{k}),\end{eqnarray} \tag{ 6 }$$

the dynamics of the displacement vector can be set into the form

$$\begin{eqnarray}&&\overrightarrow{u}(x,y,t)=\displaystyle \sum _{k=1}^{K}\displaystyle \sum _{m=1}^{M}{\overrightarrow{U}}_{{km}}(x,y)\cos ({\omega }_{k}t+{\phi }_{k}).\end{eqnarray} \tag{ 7 }$$

Substituting (6)–(7) into (1)–(3) and setting to zero the various coefficients of $\cos ({\omega }_{k}t+{\phi }_{k})$ yield the following purely spatial problems (for k = 1, ..., K)

$$\begin{eqnarray}&&{\rm{\nabla }}(\lambda {\rm{\nabla }}\cdot {\overrightarrow{U}}_{{km}})+{\rm{\nabla }}\cdot [\mu ({\rm{\nabla }}{\overrightarrow{U}}_{{km}}+{\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{\top })]\\ &&\quad +\rho {\omega }_{k}^{2}{\overrightarrow{U}}_{{km}}=-{\overrightarrow{F}}_{{km}}(x,y),\end{eqnarray} \tag{ 8 }$$

in D, with boundary conditions

$$\begin{eqnarray}&&\lambda ({\rm{\nabla }}\cdot {\overrightarrow{U}}_{{km}})\overrightarrow{n}+\mu ({\rm{\nabla }}{\overrightarrow{U}}_{{km}}+{\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{\top })\cdot \overrightarrow{n}=\overrightarrow{0}\ \mathrm{at}\ \\ &&\quad x=\pm {L}_{1},| y| \lt {L}_{2}\ \mathrm{and}\ y=\pm {L}_{2},| x| \lt {L}_{1},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\overrightarrow{U}}_{{km}}=\overrightarrow{0}\quad \mathrm{at}\ (x,y)=(\pm {L}_{1},-{L}_{2}).\end{eqnarray} \tag{ 10 }$$

Now, if the forcing terms ${\overrightarrow{F}}_{{km}}$ appearing in (6) are all real, equations (8)–(10) show that the displacement fields ${\overrightarrow{U}}_{{km}}$ are also real, which means that each monochromatic component in the expansion (7) represents a standing wave. However, depending on the frequencies ω_k and the phases ϕ_k, the whole spatio-temporal pattern (7) may also be a standing wave or a propagative wave. For instance, by applying Fourier transform, the propagative pulses considered in conventional methods can be expanded in the form (6)–(7). In the sequel, all computations will be based on the purely spatial problems (8)–(10), which means that the phases ϕ_k appearing in (6)–(7) will play no role. This makes a significant difference with conventional methods, in which the phases do play an important role. This strategy based on (8)–(10) will be followed for the sake of simplicity in the presentation, though a more general approach that takes the phases into account could improve the obtained results.

As anticipated in section 1, the aim of the present paper is to identify defects in the plate by appropriate post-process of the data obtained in some sensors when the plate is subject to appropriate forcing. Defects will manifest themselves in the values of the material properties of the plate, namely the density and the Lamé coefficients, which are assumed to be of the form

$$\begin{eqnarray}&&\rho ={\rho }_{0}-\widetilde{\rho },\quad \lambda ={\lambda }_{0}-\widetilde{\lambda },\quad \mu ={\mu }_{0}-\widetilde{\mu },\end{eqnarray} \tag{ 11 }$$

where the subscript ₀ denotes the values of the material properties for the healthy plate and the perturbations, denoted with tildes, are due to the presence of localized defects. These are assumed to be appropriately small in the L₂ norm, namely to be such that

$$\begin{eqnarray}&&{\displaystyle \int }_{D}| \tilde{\rho }{| }^{2}\,{dxdy}\ll {\displaystyle \int }_{D}{\rho }_{0}^{2}\,{dxdy},\quad {\displaystyle \int }_{D}| \tilde{\lambda }{| }^{2}\,{dxdy}\\ &&\quad \ll {\displaystyle \int }_{D}{\lambda }_{0}^{2}\,{dxdy},{\displaystyle \int }_{D}| \tilde{\mu }{| }^{2}\,{dxdy}\ll {\displaystyle \int }_{D}{\mu }_{0}^{2}\,{dxdy}.\end{eqnarray} \tag{ 12 }$$

In the present context, these conditions will hold because of the localization of the defects, in spite of the fact that the perturbations in the coefficient values need not be small in the defects themselves. In this case, equation (12) just means that the area occupied by the defects must be small compared to the area of the plate. For instance, circular defects with a diameter equal to one tenth of the plate sides, exhibit an area that is one hundred of that of the plate.

Concerning the in-plane forcing terms ${\overrightarrow{F}}_{{km}}$ appearing in (8), we shall consider a set of M simple forcing configurations, each (for m = 1, ..., M) of one of the two types:

• Point localized (PL) excitation, in which a concentrated force at a given point (x_m, y_m) is imposed. Thus, the forcing term appearing in (8) is given by
  $$\begin{eqnarray}&&{\overrightarrow{F}}_{{km}}^{\mathrm{PL}}(x,y)={\overrightarrow{a}}_{{km}}\delta (x-{x}_{m},y-{y}_{m}),\end{eqnarray} \tag{ 13 }$$
  where the vector ${\overrightarrow{a}}_{{km}}$ is the applied force and δ is the Dirac delta function. This excitation mimics the non-contact optical laser actuators.
• Rotationally symmetric radial (RSR) excitation, in which a RSR force is applied in a circle with center at a given point $({x}_{m},{y}_{m})$ and radius r. Thus, the forcing term appearing in (8) is given by
  $$\begin{eqnarray}&&{{\overrightarrow{{F}_{}^{}}}_{km}^{{\rm{R}}{\rm{S}}{\rm{R}}}}^{}{}_{}(x,y)={a}_{{km}}\overrightarrow{g}(x-{x}_{m},y-{y}_{m},r),\end{eqnarray} \tag{ 14 }$$
  where the scalar a_km is the overall intensity of the applied force and $\overrightarrow{g}(x-{x}_{m},y-{y}_{m},r)$ denotes a forcing term consisting in applying a unit force along the outward normals to the circle of radius r, centered at (x_m, y_m). This mimics piezoelectric actuators.

Likewise, the resulting synthetic data, mimicking experimental data, will be collected in N sensors (for n = 1, ..., N), each of one of the following two types (which are the counterparts of the two types of excitation considered above):

• Point localized (PL) sensing, in which the displacement vector at a given point (x_n, y_n) is measured. This sensing device mimics laser Doppler sensors and produces two scalar measurements (namely, the two components of the displacement vector) at the sensor.
• Rotationally symmetric radial (RSR) sensing in a circle with center at a given point $({x}_{n},{y}_{n})$ and radius r, measuring the spatial average of the radial displacement vector along the circle. This mimics piezoelectric sensors of radius r centered at (x_n, y_n) and produces just one scalar measurement, the averaged radial displacement vector.

The forcing and sensing types (either PL or RSR, as described above) will affect the computation of the topological derivative in the next subsection.

According to our comments above, these excitation/sensing types will be combined in section 3 in three ways, namely PL-forcing/PL sensing, RSR-forcing/RSR sensing, and RSR-forcing/PL sensing, considering M actuator configurations, which will be common to the K considered frequencies. For each of these configurations, appropriate data will be collected at N sensors, all of them of the same type (either PL or RSR). Thus, the available scalar data is K × M × (2N) for PL sensing and K × M × N for RSR sensing. However, the unknowns in the inverse problem are the spatial distributions of the perturbations $\tilde{\rho }$, $\tilde{\lambda }$, and $\tilde{\mu }$, which involve an infinite number of scalar unknowns from the purely mathematical point of view. Upon discretization via a spatial mesh, with P nodes, the number of unknowns is $3P$, which is usually much larger than the number of available uncorrelated data. This is one of the difficulties involved in inverse problems, which are generally ill-posed. An additional, even more important, difficulty is that the unknown material properties do not depend continuously on the given data.

As already explained in section 1, our approach will consist in somehow minimizing a scalar objective function, which accounts for the least squares difference between the measured data at the sensors and their computed counterpart for a given varying location and shape of the defects. The defects manifest themselves in the present context through the distributions of the perturbations $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$, as defined in (11). In principle, minimizing the objective function defined that way should provide an approximation of $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$, and thus the location and shape of the defects when these are localized.

The above mentioned objective function is defined as the sum of the contributions from the K forcing frequencies and the M actuators/sensors configurations, namely

$$\begin{eqnarray}{ \mathcal H }(\widetilde{\rho },\widetilde{\lambda },\widetilde{\mu }) & = & \displaystyle \sum _{k=1}^{K}\displaystyle \sum _{m=1}^{M}{\alpha }_{{km}}{{ \mathcal H }}_{{km}}^{\mathrm{PL}}\quad \mathrm{and}\\ { \mathcal H }(\widetilde{\rho },\widetilde{\lambda },\widetilde{\mu }) & = & \displaystyle \sum _{k=1}^{K}\displaystyle \sum _{m=1}^{M}{\alpha }_{{km}}{{ \mathcal H }}_{{km}}^{\mathrm{RSR}},\end{eqnarray} \tag{ 15 }$$

for PL and RSR sensing, respectively, which must be treated separately. Here, as above, $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$ denote the perturbations (from the healthy state, according to (11)) in these properties appearing in the damaged plate, the subscripts k and m label the frequencies and actuators/sensors configurations, respectively, the weights α_km > 0 (to be selected below) measure the individual contribution of each frequency and each emitting configuration, and ${{ \mathcal H }}_{{km}}^{\mathrm{PL}}$ and ${{ \mathcal H }}_{{km}}^{\mathrm{RSR}}$ depend on the sensing mode.

For PL sensing,

$$\begin{eqnarray}{{ \mathcal H }}_{{km}}^{\mathrm{PL}}(\widetilde{\rho },\widetilde{\lambda },\widetilde{\mu }) & = & \displaystyle \frac{1}{2}\displaystyle \sum _{n=1}^{N}\parallel {\overrightarrow{U}}_{{km}}^{\mathrm{comp}.}({x}_{{mn}},{y}_{{mn}})\\ & & -{{\overrightarrow{U}}_{{km}}^{\mathrm{meas}.}({x}_{{mn}},{y}_{{mn}})\parallel }^{2},\end{eqnarray} \tag{ 16 }$$

where $({x}_{{mn}},{y}_{{mn}})$ denote the sensor positions (which may depend on the actuators set configuration), ${\overrightarrow{U}}_{{km}}$ is the kth Fourier component of the displacement vector field produced by the mth excitation, as defined in (7), and the superscripts ${}^{\mathrm{comp}.}$ and ^meas. denote computed (by solving (8)–(10)) and measured data. In the present case, we are using synthetic data, mimicking actual measured data, which means that U^meas. will also be computed by solving (8)–(10) for the actual damaged plate. Since (either actual or synthetic) measurements are known beforehand, dependence of ${{ \mathcal H }}_{{km}}^{\mathrm{PL}}$ on $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$ comes from ${\overrightarrow{U}}_{{km}}^{\mathrm{comp}.}$, whose calculation requires solving (8)–(10), which in turn depend on $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$. In other words, U^meas. will be computed by solving (8)–(10) for the known distributions of $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$, U^comp. will be computed by solving (8)–(10) for varying distributions of while $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$, which should be selected to minimize the objective function.

Similarly, for RSR sensing,

$$\begin{eqnarray}&&{{ \mathcal H }}_{{km}}^{\mathrm{RSR}}(\widetilde{\rho },\widetilde{\lambda },\widetilde{\mu })=\displaystyle \frac{1}{2}\displaystyle \sum _{n=1}^{N}{\langle {\overrightarrow{U}}_{{km}}^{\mathrm{comp}.}-{\overrightarrow{U}}_{{km}}^{\mathrm{meas}.}\rangle }_{{x}_{{mn}},{y}_{{mn}},r}^{2},\end{eqnarray} \tag{ 17 }$$

where ${\overrightarrow{U}}_{{km}}^{\mathrm{comp}.}$ and ${\overrightarrow{U}}_{{km}}^{\mathrm{meas}.}$ are as above and for a vector field $\overrightarrow{U}$, ${\langle \overrightarrow{U}\rangle }_{{x}_{{mn}},{y}_{{mn}},r}$ denotes the average of the outward normal component of $\overrightarrow{U}$ on the circle of radius r centered at the point (x_mn, y_mn). Note that ${\langle \overrightarrow{U}\rangle }_{{x}_{{mn}},{y}_{{mn}},r}$ is a scalar quantity in spite of the fact that $\overrightarrow{U}$ is a vector field. Again, dependence of ${{ \mathcal H }}_{{km}}^{\mathrm{RSR}}$ on $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$ comes through ${\overrightarrow{U}}_{{km}}^{\mathrm{comp}.}$.

As anticipated in section 1, minimizing the objective function (15) should give a good guess of the unknowns $\widetilde{\rho }$, $\widetilde{\lambda }$, and $\widetilde{\mu }$. However, the objective function depends nonlinearly on the unknowns and this approach is problematic. Thus, it is substituted in this paper by a guess based on a more direct, non-iterative process that uses the topological derivative of the objective function ${ \mathcal H }$, which is a scalar vector field ${ \mathcal T }(x,y)$ defined at each $(x,y)\in D$ as follows. We consider a small circular defect, with center at (x, y) and radius r. The topological derivative is the limit

$$\begin{eqnarray}&&{ \mathcal T }(x,y)=\mathop{\mathrm{lim}}\limits_{r\to 0}\displaystyle \frac{{{ \mathcal H }}^{r}-{{ \mathcal H }}^{0}}{\pi {r}^{2}},\end{eqnarray} \tag{ 18 }$$

where ${{ \mathcal H }}^{r}$ and ${{ \mathcal H }}^{0}$ denote the values of the objective function for the damaged and healthy plates, respectively. In other words, the topological derivative at (x, y) measures the sensitivity of the objective function to an infinitesimal defect located at this point. The idea is that damage is most likely to occur at those points of the plate where the topological derivative exhibits negative peaks. Note that the defects considered in (18) are infinitesimal, which means that considering them as circular is only a slight restriction and the actual defects that are to be localized may have an arbitrary shape. They will be just located at those regions where the function ${ \mathcal T }$ exhibits the most pronounced negative values, which will be easily identified from convenient color maps of the topological ${ \mathcal T }$.

The topological derivative is a linear operator, which means invoking (15) that

$$\begin{eqnarray}&&{ \mathcal T }(x,y)=\displaystyle \sum _{k=1}^{K}\displaystyle \sum _{m=1}^{M}{\alpha }_{{km}}{{ \mathcal T }}_{{km}}(x,y),\end{eqnarray} \tag{ 19 }$$

where ${{ \mathcal T }}_{{km}}$ is the topological derivative of ${{ \mathcal H }}_{{km}}^{\mathrm{PL}}$ or ${{ \mathcal H }}_{{km}}^{\mathrm{RSR}}$ for PL sensing or RSR sensing, respectively, defined in (16) and (17), respectively. Using well-known computations in the literature [47], it follows that

$$\begin{eqnarray}&&{{ \mathcal T }}_{{km}}=\rho {{ \mathcal S }}_{{km}}^{\rho }+\displaystyle \frac{1}{\mu (1+\nu )}\left[\displaystyle \frac{\nu -1-2{\nu }^{2}}{2(1-{\nu }^{2})}{{ \mathcal S }}_{{km}}^{\lambda }+{{ \mathcal S }}_{{km}}^{\mu }\right],\end{eqnarray} \tag{ 20 }$$

where $\nu =\lambda /[2(\lambda +\mu )]$ is the Poisson ratio and

$$\begin{eqnarray}{{ \mathcal S }}_{{km}}^{\rho } & = & -{\omega }_{k}^{2}{\overrightarrow{U}}_{{km}}^{d}\cdot {\overrightarrow{U}}_{{km}}^{a},\\ {{ \mathcal S }}_{{km}}^{\lambda } & = & ({\rm{\nabla }}\cdot {\overrightarrow{U}}_{{km}}^{d})({\rm{\nabla }}\cdot {\overrightarrow{U}}_{{km}}^{a}),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{ \mathcal S }}_{{km}}^{\mu }=\displaystyle \frac{1}{2}({\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{d}+{\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{d\top }):\\ &&\quad ({\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{a}+{\rm{\nabla }}{\overrightarrow{U}}_{{km}}^{a\top }).\end{eqnarray} \tag{ 22 }$$

Here, ${\overrightarrow{U}}_{{km}}^{d}$ and ${\overrightarrow{U}}_{{km}}^{a}$ are the solutions of the unperturbed direct and adjoint problems, which are both given by (8)–(10) with ρ = ρ₀, λ = λ₀, μ = μ₀ (the unperturbed values of these properties for the healthy plate) and the following forcing term in the right hand side of (8), ${\overrightarrow{F}}_{{km}}$:

• For the unperturbed direct problem, ${\overrightarrow{F}}_{{km}}$ depends on the forcing type and is given by (13) and (14) for PL and RSR-forcing, respectively.
• For the unperturbed adjoint problem, ${\overrightarrow{F}}_{{km}}$ depends on the sensing type, either PL or RSR. For PL sensing at the points (x_mn, y_mn), for n = 1, ..., N, we have
  $$\begin{eqnarray}{\overrightarrow{F}}_{{km}} & = & \displaystyle \sum _{n=1}^{N}[{\overrightarrow{U}}_{{km}}^{d}({x}_{{mn}},{y}_{{mn}})-{\overrightarrow{U}}_{{km}}^{\mathrm{meas}.}({x}_{{mn}},{y}_{{mn}})]\\ & \times & \delta (x-{x}_{{mn}},y-{y}_{{mn}}),\end{eqnarray} \tag{ 23 }$$
  where δ is the Dirac delta function, as above. For RSR sensing in devices of radius r, centered at the points (x_mn, y_mn), for n = 1, ..., N, ${\overrightarrow{F}}_{{km}}$ is given by
  $$\begin{eqnarray}{\overrightarrow{F}}_{{km}} & = & \displaystyle \sum _{n=1}^{N}{\langle {\overrightarrow{U}}_{{km}}^{d}-{\overrightarrow{U}}_{{km}}^{\mathrm{meas}.}\rangle }_{{x}_{{mn}},{y}_{{mn}},r}\\ & \times & \overrightarrow{g}(x-{x}_{{mn}},y-{y}_{{mn}},r),\end{eqnarray} \tag{ 24 }$$
  where the rotationally symmetric vector field $\overrightarrow{g}$ and the mean value ${\langle \cdot \rangle }_{{x}_{{mn}},{y}_{{mn}},r}$ are as defined after (14) and (17), respectively. Recall that the azimuthal average ${\langle \cdot \rangle }_{{x}_{{mn}},{y}_{{mn}},r}$ is a scalar quantity.

Once the unperturbed direct and adjoint problems have been solved, equations (19)–(22) allow to compute the topological derivative field ${ \mathcal T }(x,y)$, which still depends on the weights α_km. Following ideas in [39], which gave very good results for the two-dimensional wave equation, we select

$$\begin{eqnarray}&&{\alpha }_{{km}}=\displaystyle \frac{1}{| \mathop{\min }\limits_{(x,y)\in {D}_{I}}{{ \mathcal T }}_{{km}}(x,y)| },\end{eqnarray} \tag{ 25 }$$

where ${D}_{I}\subset D$ is an interrogation window, which is the subset of the plate where the topological derivative will be calculated. Such interrogation window must exclude a vicinity of the actuators and sensors because the topological derivative takes very large (positive or negative) values in near these devices. This is because the solutions to the unperturbed direct and adjoint problems are very steep near the actuators and sensors, respectively, and the topological derivative includes terms (namely, the second and third terms in the right hand side of equation (20)) that are proportional to the product of the gradients of the unperturbed direct and adjoint problems. The large values of the topological derivative near the actuators/sensors could give negative peaks that would be spurious in the present context and must be avoided.

It is interesting to note that the fields ${{ \mathcal S }}_{{km}}^{\rho }$, ${{ \mathcal S }}_{{km}}^{\lambda }$, and ${{ \mathcal S }}_{{km}}^{\mu }$, defined in (21)–(22) admit another interpretation, as sensitivities of ${{ \mathcal H }}_{{km}}$ to perturbations in ρ, λ, and μ, respectively. In other words, if ρ, λ, and μ are perturbed as in (11), with conveniently small perturbations, according to (12), the individual components of the objective function defined in (16) can be linearized around $\tilde{\rho }=\tilde{\lambda }=\tilde{\mu }=0$, which (using a standard continuous adjoint formulation [48]) gives

$$\begin{eqnarray}&&{{ \mathcal H }}_{{km}}(\widetilde{\rho },\widetilde{\lambda },\widetilde{\mu })-{{ \mathcal H }}_{{km}}^{0}\simeq {\displaystyle \int }_{D}[{{ \mathcal S }}_{{km}}^{\rho }(x,y)\tilde{\rho }(x,y)\\ &&\quad +{{ \mathcal S }}_{{km}}^{\lambda }(x,y)\tilde{\lambda }(x,y)+{{ \mathcal S }}_{{km}}^{\mu }(x,y)\tilde{\mu }(x,y)]\,{dxdy},\end{eqnarray} \tag{ 26 }$$

where ${{ \mathcal H }}_{{km}}^{0}$ is the unperturbed value of ${{ \mathcal H }}_{{km}}$, quadratic and higher order terms in the small quantities $\tilde{\rho }$, $\tilde{\lambda }$, and $\tilde{\mu }$ have been neglected, and the sensitivities ${{ \mathcal S }}_{{km}}^{\rho }$, ${{ \mathcal S }}_{{km}}^{\lambda }$, and ${{ \mathcal S }}_{{km}}^{\mu }$ are precisely those given by (21)–(22), with the unperturbed direct and adjoint problems as defined above. Thus, equation (20) represents a particular linear combination of these sensitivities. Other linear combinations could give better results, as could do a selection of the weights α_km different from that in (25). However, these possible improvements are beyond the scope of the present paper, where for the sake of simplicity and brevity, the somewhat standard theory described above will be used.

To begin with, we consider a number M of PL dual emitting/receiving PL devices, which will be used to obtain (for each of the K driving frequencies) M emitting/receiving configurations. In each of these, signals are emitted from one of the devices (which acts as actuator and promotes a displacement parallel to the y axis) and received at the remaining N = M − 1 devices (where the whole displacement vector is measured), which act as sensors. Since two scalar data are produced at each PL sensor, the total number of available scalar data is 2K × M × (M − 1), which will be only moderately large (not larger than a few tens of thousands) in the applications below. In addition, these data are not completely uncorrelated, meaning that the effective number of data is even smaller. The number of unknowns, instead, are the values of ρ, λ, and μ at the nodes of the computational mesh, which is ∼3 × 0.96 = 2.88 Million unknowns, namely much larger than the number of data. This illustrates the ill-posed nature of the inverse problem. Using the topological derivative, described in section 2.2, somehow circumvents this difficulty.

3.1.1. Homogeneous plate

The first considered emitting/receiving configuration consists of 32 evenly located (as indicated with black × symbols in figure 3) PL emitting/receiving devices, which are fairly close to the boundary, namely 2 cm apart. The devices are equispaced along the lines parallel to the sides of the boundary. Concerning the driving frequencies, we use K = 15 equispaced values in the interval (28). For this emitting/receiving configuration, we consider three sets of through-hole defects, namely (a) a circular defect of radius 2.5 mm, with its center at (x, y) = (−0.1 m, 0.2 m), (b) an elliptical defect with center at (x, y) = (−0.15 m, −0.2 m) and major and minor axes 4.84 mm and 1.29 mm, parallel to the x and y axes, respectively, and (c) the circular and elliptical defects simultaneously located in the plate. Note that the two defects exhibit the same area. In figure 3 (as in all plots of the below), the actual position of the defects is indicated with white + symbols; the shapes of the defects cannot be seen in this plot due to their very small sizes. The topological derivative for these three sets of defects is as plotted in figure 3, where only negative values of the topological derivative are represented in the color map; points with positive values and points outside the interrogation window are indicated in gray. As can be seen, the negative peaks of the topological derivative locate fairly well the defects, especially in cases (a) and (b), in which only one defect is present; case (c), with two defects, is more demanding and thus the contrast of the negative peaks of the topological derivatives at the defects is less prominent. Also, the circular and elliptical defects promote comparable negative peaks in the topological derivative, which is due to the fact that both defects have been selected to exhibit the same area. If, instead, one of the through-holes exhibited a smaller area, then the (negative) peak of the topological derivative at this defect would be less prominent.

The number of emitting/receiving devices, M = 32, was chosen in the previous test problem to ensure obtaining good results, but most likely they produce highly correlated data and can be decreased. Thus, we drastically decrease this number from 32 to 4, maintaining the same defects sizes and locations as in the previous test problem. The four considered emitting/receiving devices are located close to each other and in the upper side of the plate, as indicated with black × symbols in figure 4; in fact, they are a subset of those considered in figure 3. Concerning the number K of forcing frequencies (which again are equispaced in the interval (28)), this is K = 15, as in the previous test problem, when only one defect is present, but this number is increased to K = 35 in the more demanding case in which both defects are present. As can be seen in figure 4, the topological derivative locates the defects fairly well using only four compactly located emitting/receiving devices. Thus, in the remaining of this subsection, we shall use only this reduced set of devices.

Zoom In Zoom Out Reset image size

3.1.2. Inhomogeneous plate

The test problems considered here are much more demanding than the previous ones. Thus, they will require increasing the number of frequencies if the number of emitting/receiving devices is maintained, as anticipated above.

As a first test problem, we consider the plate with an elongated through-slit, considering two cases (see figure 5):

• A centered parallel (to the x axis) slit with length 0.8 m and width 1.6 cm located along the line y = 0.17 m, and a circular defect of radius 2.5 mm and center at (x, y) = (0.1 m, −0.2 m).
• An oblique slit with length 0.7 m and width 1.4 cm, located along a line with slope 40.17° and centered at (x, y) = (0.060 9 m, 0.129 m). The defect is again circular, with radius 2.5 mm and center at (x, y) = (0.15 m, 0.1 m).

For this test problem, we consider K = 30 frequencies, equispaced in the interval (28). The topological derivatives for the parallel and oblique slits, plotted in figure 5, show that the method locates the defects reasonably well in spite of the presence of the slits. In both cases, especially for the oblique slit, the defect is not reachable by the primary emitted waves, but only by waves reflected at the boundary, which would make the application of standard SHM methods extremely problematic. In fact, through-slits are quite demanding because they totally reflect the waves, as the lateral walls do, due to the infinite jump in the acoustic impedance. Here, the slits leave a space at the borders, which allows the guided waves to pass to the other side. If, instead, the through-slits were longer, as to completely cover a section of the plate, from side to side, it would be absolutely impossible to detect any defect beyond the slits. However, this does not happen with elongated inclusions of a different material, which are considered in the next test problem, because these inclusions allow for a part of the waves to trespass the inclusion (though they are subject to diffraction).

Zoom In Zoom Out Reset image size

As a second test problem, we consider a plate with an elongated inclusion of a different material. The inclusion is parallel to the x axis and transversally covers the whole plate at y = 0.16 m from side to side. As for the thickness and material properties of the inclusion and the nature of the material, we consider two cases:

• Titanium inclusion: titanium is a fairly dense and fairly rigid material, with density ρ = 4500 kg m⁻³ and Lamé coefficients λ = 74.07 GPa and μ = 41.67 GPa (Young modulus E = 110 GPa and Poisson ratio ν = 0.32). The thickness of the inclusion is 15 mm. In this case, we consider K = 25 forcing frequencies, equispaced in the interval (28).
• High-density polyethylene (HDPE) inclusion: now, both the density (ρ = 960 kg m⁻³) and especially the rigidity (λ = 1.777 GPa and μ = 0.390 1 GPa, which correspond to Young modulus E = 1.1 GPa and Poisson ratio ν = 0.41) are much smaller than those for aluminum. The thickness of the inclusion is 2.5 mm, much smaller than in the previous case since, as explained below, this is a more demanding case. Because of this, the number of forcing frequencies (again equispaced in the interval (28)) is also larger than in the previous case, namely K = 60.

In both cases, the defect that is to be detected is circular, with radius 2.5 mm and centered at (x, y) = (0.1 m, −0.2 m).

The topological derivatives for these two cases are as plotted in figure 6. As can be seen in this figure, the most demanding case is that of HDPE, which involves a smaller thickness of the inclusion and requires a larger number of frequencies. The reason for being so demanding is that both the density and the Lamé coefficients are much smaller than that of the aluminum plate. This makes this case closer to the elongated through-slit, which as anticipated is the most demanding inhomogeneity.

In any event, it is remarkable how the method is able to capture defects that are beyond the elongated through-slits covering a wide portion of the transversal section of the plate and elongated inclusions of different material covering the whole section of the plate, from side to side.

Zoom In Zoom Out Reset image size

Let us now consider dual emitting/receiving devices of the RSR type. Since RSR sensors only capture one scalar quantity (instead of two), a larger (than in the PL case) number of devices and a larger number of frequencies are needed in the present case. Thus, we consider eight RSR devices of radius 5 mm and centers located at 2 cm from the y = 0.5 m side of the plate, equispaced in the x coordinate (along the horizontal line y = 0.48 m).

For the homogeneous plate, we consider the same defects configurations as in figure 4 and an appropriate number of forcing frequencies. The obtained topological derivatives are given in figure 7. As can be seen in this figure, the topological derivative identifies the defects reasonably well, especially when only one defect is present, though the number of required frequencies is larger than in the case considered in figure 4, as expected. Note that two additional blue regions appear in the case of two defects. However, the topological derivative peaks at these additional regions are less intense than near the actual position of the defects (which means that the topological derivative locates well the defects, as stated).

Zoom In Zoom Out Reset image size

Concerning the case of the inhomogeneous plates considered in section 3.1.2, with the same elongated slits or inclusions and defects as in figures 5 and 6, the topological derivatives for the present actuators/sensors configuration are given in figure 8. As can be seen in this figure, the topological derivative locates the defects quite well using an appropriate number of forcing frequencies (again, equispaced in the interval (28)), namely 60 and 40 frequencies for the parallel and oblique through-slits, respectively, and 40 and 60 frequencies for the titanium and HDPE inclusions, respectively. Note that the number of frequencies is larger for the HDPE inclusion than for the titanium inclusion because, as already mentioned, this material is much more demanding in connection with inclusions.

Zoom In Zoom Out Reset image size

Summarizing, the results obtained for this new type of dual RSR actuators/RSR sensors are consistent with those obtained in section 3.1.2 for dual PL actuators/PL sensors. In particular, as anticipated, since RSR sensors collect a smaller amount of data than PL sensors, the numbers of required devices and frequencies are larger for the present case.

Let us now consider the case of RSR actuators and PL sensors. Since, as anticipated, PL sensors collect a larger amount of data (namely, twice as much) than RSR sensors, the present case is less demanding than that considered in the last subsection. This means that we may decrease either the number of actuators or the number of sensors, which must be located in different positions (but again 2 cm apart from the boundary) because they are of a different type among each other. We select the same four RSR actuators as in section 3.2 and four PL sensors aligned with the actuators in the line y = 0.48 m and equispaced along the x coordinate in the right part of this line. Actuators and sensors are indicated with black asterisks and black × symbols, respectively, in figures 9 and 10, which show the topological derivatives for the homogeneous and inhomogeneous plates, respectively, considering the same defects (and slits or inclusions) as in figures 7 and 8, respectively. As can be seen, the topological derivative locates the defects quite well in all considered cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Summarizing the results obtained in this subsection, RSR emision and PL sensing gives results similar to PL emision and PL sensing because, as anticipated, these two methods provide comparable number of data (which mainly depends on sensing).