**Applications on Ultrasonic Wave**

Editor

**Jaesun Lee**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Jaesun Lee Changwon National University Korea

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/Applications Ultrasonic Wave).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1120-7 (Hbk) ISBN 978-3-0365-1121-4 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

**Jaesun Lee** is a professor in the School of Mechanical Engineering and currently serves as the Director of Smart Manufacturing Engineering on the regional innovation platform and Vice Dean of Planning, Changwon National University, Korea. He obtained his PhD from Pusan National University, Korea, in the year 2015. He worked at Korea Railroad Research Institute (KRRI), Korea, and he joined the school of mechanical engineering of Changwon National University (CWNU), South Korea, as a tenure track faculty. His research includes ultrasonic NDE and structural health monitoring, long range guided wave inspection technology, theoretical wave scattering analysis, tomographic defect imaging, structural analysis and nonlinear ultrasonics material characterization. He had a global PhD fellowship from the Korean government in 2011 for 3 years. Moreover, he has won more than 11 best paper and presentation awards from the Korean Society for Precision Engineering, the Korean Society of Mechanical Engineers, the Korean Society for Nondestructive Testing and the Korean Society for Railway. He is currently leading the NDE research group of CWNU and several national funding projects, including the National Regional Leading Research Center.

## **Preface to "Applications on Ultrasonic Wave"**

This Special Issue presents the latest advances in the field of applications on ultrasonic waves, including nondestructive methodology to evaluate the mechanical properties, damage states, and material condition of engineering structures. This Special Issue is composed of five papers covering new insights into inspection applications on road, pipeline and material characterization.

The advent of material characterization and non-destructive evaluation by ultrasonic wave applications is detailed in this book. Features related to propagation and the scattering of ultrasonic waves, such as wave velocities, dispersion, scattered amplitudes, and attenuation, show high sensitivity to material condition. Ultrasonic waves are also widely used in medicine to obtain images of internal body structures, such as muscles, tendons, blood vessels, joints, and organs. The mechanism is that ultrasound pulses are sent into the tissue using a probe, and the reflected signals are recorded and analyzed in order to build the desired images of internal structures.

> **Jaesun Lee** *Editor*

*Article*

## **Dynamic Non-Destructive Evaluation of Piezoelectric Materials to Verify on Accuracy of Transversely Isotropic Material Property Measured by Resonance Method**

#### **Yu-Hsi Huang 1,\*, Chien-Yu Yen <sup>1</sup> and Tai-Rong Huang <sup>2</sup>**


Received: 23 June 2020; Accepted: 21 July 2020; Published: 23 July 2020

**Featured Application: This paper provides the measurement methods and corresponding calculating formulas for piezoelectric material constants. The completely orthotropic material property, i.e., dielectric, piezoelectric, and elastic constants, is obtained and verified by piezoelectric dynamic characteristics. By means of use in anisotropic material property, the mode shape and corresponding natural frequency were determined by finite element numerical calculation. The di**ff**erence is rare in comparison with the results obtained by optically dynamic vibration measurement. The influence on self-heating phenomenon of soft and hard piezoelectric ceramics is also discussed in the study.**

**Abstract:** In many engineering applications of piezoelectric materials, the design and prediction of the dynamic characteristics depends on the anisotropic electromechanical material property. Through collecting the complete formula in literature and listing all the prepared specimens, transversely isotropic material constants were obtained and verified by dynamic non-destructive evaluation in the paper. The IEEE (Institute of Electrical and Electronics Engineers) resonance method was applied to measure and calculate the orthotropic material constants for piezoelectric ceramics. Five specimens need to be prepared for the measurements using an impedance analyzer, in order to obtain the resonant and anti-resonant frequencies from the modes of thickness extension, length-extension, thickness-shear extension, length-thickness extension, and radial extension. The frequencies were substituted into the formulas guided on the IEEE standard to determine the elastic, dielectric, and piezoelectric constants. The dynamic characteristics of soft and hard piezoelectric ceramics in the results from the finite element method (FEM), which is analyzed from the anisotropic material constants of the resonance method, were verified with the mode shapes and natural frequencies found by experimental measurements. In self-heating, considered as operating on resonant frequencies of piezoelectric material, the resonant frequency and corresponding mode shape calculated by the material constants from resonance method in FEM are more accurate than the material property provided by the manufacturer and literature. When the wide-bandwidth frequency is needed to design the application of piezoelectric ceramics, this study completely provided the measurement method and dynamic verification for the anisotropic electromechanically material property.

**Keywords:** piezoelectric materials; transversely isotropic material; IEEE resonance method; resonant and anti-resonant frequencies; dynamic non-destructive evaluation

#### **1. Introduction**

Piezoelectric materials are smart materials which serve as crucial components in sensors or actuators. When the need arises for vibration control in lightweight flexible structures with complex shapes, piezoelectric sensors and actuators can be used to prevent adding to the overall mass of the structure and to achieve precision control. Piezoelectric materials have superior electromechanical properties, are lightweight and easy and inexpensive to produce, and have low power consumption [1–4]. For instance, the puncture systems used in cellular biology research have precise and fine vibration control [5]. In addition to applying the vibration characteristics of inverse piezoelectricity, the puncture systems also apply energy extraction systems with direct piezoelectric effects. Recently, vibration energy harvesting devices using piezoelectric materials have garnered the attention of researchers. These devices convert the vibration energy of mechanical structures into usable electrical energy, which can detect the inhibiting effects of vibration signals [6]. Energy harvesting systems based on lead zirconate titanate (PZT) can be embedded into vibrating objects to generate tens to hundreds of microwatts of energy for low-power integrated circuits (ICs). Thus, piezoelectric materials can be combined with low-power ICs to create continuously-working miniature wireless sensors, such as implantable biomedical sensors and wearable devices [7], wireless sensor networks [8], and smart buildings and special structures [9–11]. The vibration characteristics of piezoelectric materials can be analyzed in advance using finite element software for greater convenience in the optimization of piezoelectric effects. To obtain accurate analysis results, the boundary conditions of the model must be clearly defined, and the right piezoelectric material parameters must be used. The parameters provided with commercially-available piezoelectric materials are often inadequate for numerical calculations. Complete dynamic characteristics cannot be derived from the mechanical parameters of anisotropic materials, and thus, measuring effective and complete material parameters is key to analysis accuracy.

In this study, the material constants of the piezoelectric specimens were measured using the resonance method (RM). This method has been widely applied. In 1954, Mason et al. [12] used the resonance characteristics of an equivalent circuit simulating a piezoelectric crystal and employed static methods, quasi-static methods, dynamic methods, and hydrostatic methods to measure the elastic, piezoelectric, and dielectric constants of 45◦ Rochelle salt. Based on these measurement methods, they then developed RM to measure the material properties of piezoelectric crystals. In this method, the elastic and piezoelectric constants of the piezoelectric crystal are measured and the dielectric constants are inversely calculated from measured capacitance. Although the compliance *SE* <sup>11</sup> and *SD* <sup>11</sup>, electromechanical coupling coefficients *<sup>k</sup>*<sup>33</sup> and *<sup>k</sup>*, and permittivity <sup>ε</sup>*<sup>T</sup>* <sup>33</sup> were obtained, the measurement was lack of guideline for prepared the standard specimens. In 1961, the IRE Standards on Piezoelectric Crystals [13] presented five specimens with different geometric shapes, aspect ratios, polarization orientations, and electrode surface orientations. The mode shapes corresponding to these five specimens were the thickness extension (TE) mode, length extension (LE) mode, thickness shear extension (TS) mode, length thickness extension (LTE) mode, and radial extension (RAD) mode. As polarized PZT has a hexagonal lattice (6 mm), it has 10 material constants, including five elastic constants, three piezoelectric constants, and two dielectric constants. With the TE, LE, TS, LTE, and RAD modes, the resonant and anti-resonant frequencies can be measured, and then the material constants of its various forms can be inversely calculated.

In 1987, the IEEE (Institute of Electrical and Electronics Engineers) Standards [14] referred to the RM in the IRE Standards and established a complete resonance measurement method with complete assumptions regarding constitutive equations (d-form, e-form, g-form, and h-form) and simplified specimen dimensions for piezoelectric materials and five specimens with different geometric shapes, aspect ratios, polarization orientations, and electrode surface orientations. Some formula listed in Refs [13,14] are more complicated, therefore the simplified form was proposed by subsequent researcher. In 1989, Wang et al. [15] applied the thickness vibration theory in RM to the vibrations in a piezoelectric plate in which the length and width are assumed to be much greater than the thickness that they approach infinity. They used α-quartz with two different cuts to measure the elastic, piezoelectric, and

dielectric constants of α-quartz. The material constants of Y-cut quartz were studied in the paper; hence, the method is not provided for obtained the complete transversely piezoelectric isotropic material constants. In 1998, Cao et al. [16] conducted a detailed study on specimens with two types of shear mode shapes: thickness shear and length shear. They examined specimens with different aspect ratios and compared their measurements. The experiment results revealed that only the length-to-thickness ratio has a direct relationship and influences the values of the resonant and anti-resonant frequencies. A pure shear mode shape can be generated at l/20t, thereby reducing the errors in the resonant and anti-resonant frequency measurements.

We used an impedance analyzer to measure the resonant and anti-resonant frequencies of piezoelectric materials and substituted them into the equations listed in RM to obtain 10 independent material constants. Although the material parameters of piezoceramics are practically independent from frequency in the order of several MegaHertz [17–19], in the frequencies of the order of kilohertz—several tens of kilohertz—the dependence on frequency is very strong. It performs on the main resonances in longitudinal, transverse, shear, or radial modes. We then substituted the measured material constants into finite element numerical calculations to verify the in-plane and out-of-plane resonant frequencies and mode shapes of the piezoelectric materials. The experimental measurements of the out-of-plane resonant frequencies and mode shapes were obtained using electronic speckle pattern interferometry (ESPI), whereas the in-plane vibrations were measured using ESPI and an impedance analyzer. A comparison of the experiment results and the numerical calculations revealed fairly good accuracy in the material constants obtained using RM.

#### **2. Theory Based on Resonance Method**

From the linearly electro-thermo-elastic model [20], the constitutive equations are represented as

$$T\_{ij} = \varepsilon\_{ijkl} S\_{kl} - \varepsilon\_{ijk} E\_k - \alpha\_{ij} \Theta \tag{1}$$

$$D\_i = \varepsilon\_{i\bar{j}k} S\_{kl} + \varepsilon\_{i\bar{j}}^S E\_{\bar{j}} - \lambda\_i \Theta \tag{2}$$

$$
\eta = \alpha\_{kl} S\_{kl} + \lambda\_k E\_k + \mathcal{C}\_\vartheta \theta \tag{3}
$$

$$q\_{\dot{i}} = -\mathbb{K}\_{\dot{i}\dot{j}}\theta\_{\dot{,}\dot{j}}\tag{4}$$

In these equations, *Tij*, *Skl*, *Ek*, *Di*, *qi* and θ are the components of stress, strain, electric field, electric displacement, heat flux and temperature difference with initial condition, respectively. The *cijkl*, *eijk*, α*ij*, ε*<sup>S</sup> ij*, λ*i*, *C*<sup>θ</sup> and *Kij* are the components of elastic, piezoelectric, thermal expansion, permittivity, pyroelectric coefficients, and the thermal conductivity tensor, respectively. If neglecting the thermo-elastic and thermo-electrical effects, the piezoelectric *e*-form is

$$\begin{cases} \begin{array}{l} T\_{ij} = \mathfrak{c}\_{ijkl} S\_{kl} - \mathfrak{e}\_{ijk} E\_k \\ D\_i = \mathfrak{e}\_{ijk} S\_{kl} + \mathfrak{e}\_{ij}^S E\_j \end{array} \end{cases} \tag{5}$$

and the *d*-form constitutive equations are, respectively,

$$\begin{cases} S\_{ij} = s\_{ijkl} T\_{kl} + d\_{kij} E\_k \\\ D\_i = d\_{ikl} T\_{kl} + \varepsilon\_{ik}^T E\_k \end{cases} \tag{6}$$

The linear piezoelectric constitutive equations for a piezoceramic material with crystal symmetry class *C6mm* are presented in the components of elastic (*c<sup>E</sup> pq*), piezoelectric (*epq*), and dielectric coefficients (ε*<sup>S</sup> pq*) respectively.

The IEEE Standards [14] indicate that measuring the complete material constants of piezoelectric ceramic material requires piezoelectric specimens with five different mode shape characteristics (TE, LE, TS, LTE, and RAD), each with its own geometric shape, aspect ratio, polarization orientation, and electrode surface orientation. We used an impedance analyzer (model: HP-4194A) to measure their resonant and anti-resonant frequencies and inversely calculate the material constants of the piezoelectric materials. Based on Equations (7), (8), (12), (13), (19), (20), (22), and (23), the characteristics of the impedance measurements for TE, LE, TS, and LTE are as shown in Figure 1. Only one set of resonant and anti-resonant frequency values (*fr*, *fa*) needs to be substituted into Equations (25) and (26). In Equations (25) and (26) for RAD, the impedance measurement characteristics, as shown in Figure 1, have two sets of resonant and anti-resonant frequency values (*f<sup>r</sup>* <sup>1</sup> , *<sup>f</sup> <sup>a</sup>* <sup>1</sup> , *<sup>f</sup><sup>r</sup>* <sup>2</sup> , *<sup>f</sup> <sup>a</sup>* <sup>2</sup> ). Below, we obtain the material constants from the specimens of the five vibration modes according to the explanations in Ref [14].

**Figure 1.** Schematics in resonant and anti-resonant frequencies obtained by impedance analyzer.

#### *2.1. Thickness Extension (TE)*

Table 1 presents the geometric shape and material constants of a thin square plate in which the directions of polarization and vibration are oriented along the depth (thickness). When the aspect ratios of the component satisfy *l* > 10*t* and *w* > 10*t*, we use five specimens with dimensions *l* = 46 ± 0.1 mm, *w* = 46 ± 0.1 mm, and *t* = 1.2 mm, and the following equations can be used to calculate their material constants:

$$k\_t^2 = \frac{\pi}{2} \frac{f\_r}{f\_a} \cot \left(\frac{\pi}{2} \frac{f\_r}{f\_a}\right) \tag{7}$$

$$c\_{33}^D = 4\rho \left(tf\_a\right)^2\tag{8}$$

$$c\_{33}^E = c\_{33}^D \left(1 - k\_t^2\right) \tag{9}$$

where *kt* is the electromechanical coupling coefficient by thickness extension mode, *fr*is resonant frequency, *fa* is anti-resonant frequency, *t* is the thickness of specimen, and ρ is the density of specimen.

For the dielectric constant, the capacitance measured using an impedance analyzer at a specific frequency (1 kHz) can be substituted into the following equations, which will produce ε*<sup>T</sup>* <sup>33</sup>, the dielectric constant in constant stress, and ε*<sup>S</sup>* <sup>33</sup>, the dielectric constant in constant strain:

$$
\varepsilon\_{33}^T = \frac{C\_{l'}t}{A} \tag{10}
$$

$$
\varepsilon\_{33}^{\rm S} = \frac{\mathcal{C}\_{\rm h} \cdot \mathbf{t}}{A} \tag{11}
$$

where *A* is the area of the electrode, *Cl* is the capacitance under 1kHz, and *Ch* is the capacitance at twice the anti-resonant frequency (2 *fa*).

**Table 1.** Specification on vibration modes of standard piezoelectric elements and their correspondent measured material constants [14].


**E**: Electric-field direction; **P**: Polarization direction.

#### *2.2. Thickness Shear Extension (TS)*

In this mode, we have a thin square plate in which the polarization is oriented along the length, and shear vibrations occur along the depth. In other words, the orientation of polarization is perpendicular to the direction of the electric field, causing the piezoelectric component to display shear vibrations. Table 1 presents the geometric shape and material constants. When the aspect ratios of the component

satisfy *l* > 20*t* and *w* > 10*t*, we use five specimens with dimensions *l* = 9 ± 0.05 mm, *w* = 8 ± 0.05 mm, and *t* = 0.52 mm, and the following equations can be used to calculate their material constants:

$$k\_{15}^2 = \frac{\pi}{2} \frac{f\_r}{f\_a} \cot \left(\frac{\pi}{2} \frac{f\_r}{f\_a}\right) \tag{12}$$

$$c\_{44}^D = 4\rho (tf\_a)^2\tag{13}$$

$$
\sigma\_{44}^E = \sigma\_{44}^D \left(1 - k\_{15}^2\right) \tag{14}
$$

For PZT material, the following relations can be used to derive *S<sup>D</sup>* <sup>44</sup>, the elastic compliance constant with constant electric displacement, and *SE* <sup>44</sup>, the elastic compliance constant in a constant electric field:

$$S\_{44}^D = 1/c\_{44}^D\tag{15}$$

$$S\_{44}^{E} = 1/c\_{44}^{E} \tag{16}$$

where *k*<sup>15</sup> is the electromechanical coupling coefficient of the thickness shear mode shape.

For the dielectric constant, the capacitance measured using an impedance analyzer at a specific frequency (1 kHz) can be substituted into the following equations, which will produce ε*<sup>T</sup>* <sup>11</sup>, the dielectric constant in constant stress, and ε*<sup>S</sup>* <sup>11</sup>, the dielectric constant in constant strain:

$$
\varepsilon\_{11}^T = \frac{\mathbf{C}\_L \cdot \mathbf{t}}{A} \tag{17}
$$

$$
\varepsilon\_{11}^{\mathcal{S}} = \frac{\mathcal{C}\_{H} \cdot t}{A} \tag{18}
$$

*CL* is the capacitance measured at 1 kHz, and *CH* is the capacitance measured at twice the anti-resonant frequency (2 *fa*).

#### *2.3. Length Extension (LE)*

In this mode, we have a thin square column in which the polarization and vibrations are oriented along the direction of the length. Table 1 presents the geometric shape and material constants. When the aspect ratios of the component satisfy *l* > 5*w*<sup>1</sup> and *w* > 5*w*2, we use five specimens with dimensions *l* = 10 ± 0.1 mm, *w*<sup>1</sup> = 1.4 ± 0.1 mm, and *w*<sup>2</sup> = 0.51 mm, and the following equations can be used to calculate their material constants, where *k*<sup>33</sup> denotes the electromechanical coupling coefficient of the longitudinal length mode shape, and *l* is the length of specimen.

$$k\_{33}^2 = \frac{\pi}{2} \frac{f\_r}{f\_a} \cot \left(\frac{\pi}{2} \frac{f\_r}{f\_a}\right) \tag{19}$$

$$S\_{33}^D = \frac{1}{4\rho \left(lf\_a\right)^2} \tag{20}$$

$$S\_{33}^{E} = \frac{S\_{33}^{D}}{1 - k\_{33}^{2}}\tag{21}$$

#### *2.4. Length Thickness Extension (LTE)*

In this mode, we have a thin square plate in which the polarization is oriented along the depth and lateral vibrations occur along the length. In other words, the vibrations occur along the length and perpendicular to the polarization orientation. Table 1 presents the geometric shape and material constants. When the aspect ratios of the component satisfy *l* > 10*t*, *l* > 3*w*, and 3*t* > *w*, we use five specimens with dimensions *l* = 31.8 ± 0.1 mm, *w* = 3 ± 0.1 mm, and *t* = 1.21 mm, and the following equations can be used to calculate their material constants:

$$\frac{k\_{31}^2}{k\_{31}^2 - 1} = \frac{\pi}{2} \frac{f\_a}{f\_r} \cot \left(\frac{\pi}{2} \frac{f\_a}{f\_r}\right). \tag{22}$$

$$S\_{11}^E = \frac{1}{4\rho \left(lf\_r\right)^2} \tag{23}$$

$$S\_{11}^D = S\_{11}^E \left(1 - k\_{31}^2\right) \tag{24}$$

where *k*<sup>31</sup> denotes the electromechanical coupling coefficient of the thickness shear mode shape.

#### *2.5. Radial Extension (RAD)*

In this mode, we have a thin round plate in which the polarization and vibrations are both oriented along the depth. Table 1 presents the geometric shape and material constants. The aspect ratios of the component must satisfy <sup>∅</sup> <sup>&</sup>gt; <sup>20</sup>*t*, and we use five specimens with dimensions <sup>∅</sup> <sup>=</sup> 19.9 <sup>±</sup> 0.05 mm and *t* = 0.81 mm. In 1973, Meitzler, O'Bryan, and Tiersten [21] derived the constants of radial vibrations in round plates. The formula for *kp*, the electromechanical coupling coefficient of the radial vibrations in a piezoelectric round plate is

$$\frac{k\_P^2}{1 - k\_P^2} = \frac{\left(1 - \sigma^E\right) l\_1 \left[\eta\_1 \left(1 + Hf/f\_I\right)\right] - \eta\_1 \left(1 + \Delta f/f\_I\right) l\_0 \left[\eta\_1 \left(1 + \Delta f/f\_I\right)\right]}{\left(1 - \sigma^E\right) l\_1 \left[\eta\_1 \left(1 + \Delta f/f\_I\right)\right]}\tag{25}$$

where *kp* denotes the electromechanical coupling coefficient of the radial mode shape; Δ*f* = *fa* − *fr* is the difference between the anti-resonant frequency and the resonant frequency; *J*<sup>0</sup> is Bessel Function of First Kind and Zero Order; *J*<sup>1</sup> is Bessel Function of First Kind and First Order; σ*<sup>E</sup>* is Poisson's Ratio; and η<sup>1</sup> is the first positive root from the function 1 − σ*<sup>E</sup> J*1(η) = η*J*0(η). However, in 2005, Zhang, Alberta, and Eitel [22] modified the electromechanical coupling relation above using approximation into

$$k\_P^2 = \frac{1}{P} \cdot \frac{f\_a^2 - f\_r^2}{f\_a^2} \tag{26}$$

where

$$P = \frac{2\left(1 + \sigma^E\right)}{\left\{\eta\_1^2 - \left[1 - \left(\sigma^E\right)^2\right]\right\}}\tag{27}$$

Based on RM by IEEE standard [14,15], several formula were simplified also by Zhang et al. [22] and we used them to calculate some parameters, i.e., ε*<sup>T</sup>* <sup>33</sup> <sup>ε</sup>*<sup>S</sup>* <sup>33</sup>, <sup>ε</sup>*<sup>T</sup>* 11.

RM uses the vibration characteristics corresponding to the five types of piezoelectric components above to derive the elastic, piezoelectric, and dielectric constants. Furthermore, Poisson's ratio is defined as the ratio of lateral strain to longitudinal strain:

$$
\sigma^E = -\frac{S\_{12}^E}{S\_{11}^E} \tag{28}
$$

where the superscript *E* indicates values measured in a fixed electric field. The negative indicates that the two vibrations are opposite in direction, which means that when the longitudinal vibration is extending, the lateral vibration is contracting, and when longitudinal vibration is contracting, the lateral vibration is extending. Obtaining the Poisson's ratio generally requires the measurement of the fundamental frequency and one overtone frequency of the piezoelectric round plate. The formula used is 1 − σ*<sup>E</sup> J*1(η) = η*J*0(η), the source and usage of which are explained in Ref [23,24]. Then, η<sup>1</sup> and η<sup>2</sup> are calculated, and we can look up the fundamental frequency *fr* = *f<sup>r</sup>* <sup>1</sup> and the overtone frequency *<sup>f</sup><sup>r</sup>* 2 in the table provided in the Appendix to determine the Poisson's ratio σ*E*, which is substituted into Equation (25) to derive *kp*, the electromechanical coupling coefficient of the radial mode shape.

Although RM can be used to obtain the material constants of piezoelectric materials, there are still some constants that cannot be measured via experiment, and some conversions between constants are required to derive the constants that cannot be directly measured before all of the elastic, piezoelectric, and dielectric constants are obtained. Below is the solved equation including the piezoelectric strain constant:

$$d\_{ij}^2 = k\_{ij}^2 \varepsilon\_{ii}^T s\_{jj}^E \tag{29}$$

where *ij* = 31, 33, 15. Based on the equation above, if *kij*, ε*<sup>T</sup> ij*, and *<sup>s</sup><sup>E</sup> ij* (i.e., the electromechanical coupling coefficient, the dielectric constant in constant stress, and the elastic compliance constant in a constant electric field) are known, then *d*31, *d*33, and *d*15, (i.e., the piezoelectric strain constants) can be obtained. Next, we can use the following equations to obtain the other elastic and compliance constants:

$$\mathbf{s}\_{12}^{E} = -\boldsymbol{\sigma}^{E} \cdot \mathbf{s}\_{11}^{E} \tag{30}$$

$$\mathbf{s}^{D}\_{12} = \mathbf{s}^{E}\_{12} - k\_{31}^{2} \cdot \mathbf{s}^{E}\_{11} \tag{31}$$

$$s\_{\theta\theta}^{E} = s\_{\theta\theta}^{D} = \frac{1}{c\_{\theta\theta}^{E}} = \frac{1}{c\_{\theta\theta}^{D}} = \mathcal{Z} \{s\_{11}^{E} - s\_{12}^{E}\} \tag{32}$$

$$s\_{13}^D = -\left\{\frac{s\_{11}^D + s\_{12}^D}{2} \mathbf{\dot{s}} \begin{pmatrix} s\_{33}^D - \frac{1}{c\_{33}^D} \end{pmatrix} \right\}^{0.5} \tag{33}$$

$$s\_{13}^E = s\_{13}^D + k\_{31} \cdot k\_{33} \left(s\_{33}^E \cdot s\_{11}^E\right)^{0.5} \tag{34}$$

$$\mathbf{c}\_{11}^{E} = \frac{\mathbf{s}\_{11}^{E}\mathbf{s}\_{33}^{E} - \left(\mathbf{s}\_{13}^{E}\right)^{2}}{\left(\mathbf{s}\_{11}^{E} - \mathbf{s}\_{12}^{E}\right)\left[\mathbf{s}\_{33}^{E}\left(\mathbf{s}\_{11}^{E} + \mathbf{s}\_{12}^{E}\right) - 2\left(\mathbf{s}\_{13}^{E}\right)^{2}\right]} \tag{35}$$

$$\mathbf{c}\_{12}^{E} = \frac{-\mathbf{s}\_{12}^{E}\mathbf{s}\_{33}^{E} + \left(\mathbf{s}\_{13}^{E}\right)^{2}}{\left(\mathbf{s}\_{11}^{E} - \mathbf{s}\_{12}^{E}\right)\left[\mathbf{s}\_{33}^{E}\left(\mathbf{s}\_{11}^{E} + \mathbf{s}\_{12}^{E}\right) - 2\left(\mathbf{s}\_{13}^{E}\right)^{2}\right]} \tag{36}$$

$$\mathbf{c}\_{13}^{E} = \frac{-\mathbf{s}\_{13}^{E}}{\mathbf{s}\_{33}^{E} \left(\mathbf{s}\_{11}^{E} + \mathbf{s}\_{12}^{E}\right) - 2 \left(\mathbf{s}\_{13}^{E}\right)^{2}} \tag{37}$$

Using Equations (30) through (37), we can obtain all of the elastic and compliance constants. We can then use ε*<sup>T</sup> ik* <sup>=</sup> <sup>ε</sup>*<sup>S</sup> ik* + *dipekp* to obtain the piezoelectric stress constant *eij*. The primary calculation formulas are as follows:

$$d\mathbf{c}\_{31} = d\_{31} \left( c\_{11}^{E} + c\_{12}^{E} \right) + d\_{33} c\_{13}^{E} \tag{38}$$

$$\mathbf{e}\_{33} = 2d\_{31}\mathbf{c}\_{12}^{E} + d\_{33}\mathbf{c}\_{33}^{E} \tag{39}$$

$$c\_{15} = d\_{15} c\_{44}^E \tag{40}$$

#### **3. Specimens and Experimental Techniques**

This study performed the experimental measurements and material constant calculations of two piezoelectric ceramic materials using RM. The in-plane and out-of-plane vibration characteristics were measured using instruments to verify the accuracy of the material constants obtained using finite element analysis. Below, we introduce the piezoelectric materials and instruments used in this study.

#### *3.1. Prepation on Piezoelectric Ceramics*

To establish complete piezoelectric ceramic material constants, we used RM to perform measurements of five specimens with specific geometric shapes, aspect ratios, polarization orientations, and electrode surface orientations. As shown in Table 1, the specific mode shapes include the TE, LE, TS, LTE, and RAD modes.

We procured piezoelectric ceramic specimens that fit the specifications of the TE and RAD modes in RM from American Piezo Ceramics (APC) International, Ltd. and the Fuji Ceramics Corporation. The piezoelectric specimens for the LE and LTE modes were cut to the required aspect ratios using ultrasonic cutting. Finally, for the TS mode, we ground our own specimen, removed the electrodes, and redistributed the electrodes in the direction normal to that of polarization so that the aspect ratio, polarization orientation, and the direction of the electric field created by the electrodes would meet specifications for TS mode testing.

#### *3.2. Impedance Analysis*

We derived material constants using the theoretical equations in RM. The resonant and anti-resonant frequencies substituted into the equations were obtained by measuring the mode shapes with specific structures, namely the TE, LE, TS, LTE, and RAD modes, using an impedance analyzer.

An impedance analyzer uses fixed voltages to excite vibrations in the target object and then measures the continuous changes in electrical impedance at various frequencies, as shown in Figure 1. The electrical impedance can be regarded as the total impedance applied by the resistors, inductors, and capacitors in the target object on the alternating current. Coupling characteristics are evident in the conversion between electrical and mechanical energy in piezoelectric materials. Thus, the electrical impedance of piezoelectric materials differs from that of components without piezoelectric properties. Generally, the results will display a resonant frequency (*fresonance*, *fr*) or anti-resonant frequency (*fanti*−*resonance*, *fa*) at certain resonant frequencies in frequency-impedance characteristic diagrams.

#### *3.3. Electronic Speckle Pattern Interferometry (ESPI)*

We input the material constants derived from measurements and calculations using RM and the material parameters provided by manufacturers or literature into finite element software to analyze vibration characteristics and obtain resonant frequency values and mode shapes. Next, we used the vibration characteristics obtained using optical non-destructive testing to verify the accuracy of RM.

Derived from holography, ESPI involves the casting of two coherent light beams on the target object and then the differences and changes in the optical path lengths are used to obtain the deformation data of the target object. The measurement method commonly used when ESPI is applied to vibrations is the time-averaging method. Within the exposure time of the CCD (charge-coupled device), images of the vibrating object are captured at different times, and a zero-order Bessel function is used to modulate the interference fringes. The brightest region in the interference images is the nodal line, where the displacement is zero, and the light and dark fringes show the contours of displacement. The resolution of the measurement system depends on the light source wavelength. We set up a sub-micron level ESPI system so as to obtain the full-field distributions of micro-vibration displacements in the interference patterns [25,26]. We then compared the ESPI measurement results of the low-frequency out-of-plane vibrations and the high-frequency in-plane vibrations in the piezoelectric specimens with the FEM analysis results derived from the piezoelectric material constant inputs.

#### 3.3.1. Out-Of-Plane Vibration

As shown in Figure 2, a beam splitter is used to create two coherent laser beams. The one applied to the object surface is the object beam, whereas the other, known as the reference beam, is projected on a reference plate so that it scatters to create a speckled image. The two coaxial beams reflect and focus on the photosensitive plane of a CCD camera and create interference images. Using differencing calculations to process the images then swiftly produces interference fringe patterns on the computer screen, thereby providing real-time measurements during experiments.

**Figure 2.** Out-of-plane setup in ESPI.

#### 3.3.2. In-Plane Vibration

As shown in Figure 3, a beam splitter creates two laser beams which reach the target object from different sides but at identical angles of incidence with equal optical path lengths. Spatial filters then scatter the laser beams, and a CCD camera records the optimal interference images which undergo computer processing and produce vibration interference data.

**Figure 3.** In-plane setup in electronic speckle pattern interferometry (ESPI).

#### **4. Finite Element Method**

We used the commercial software ABAQUS to perform finite element method (FEM) analysis calculations. For the material constants, we used the calculation results obtained using RM and the parameters provided by manufacturers and literature. It is noted that errors of about 6%-30% are produced on the natural frequencies of the first three modes, when the element used only one-layer thickness or the linear element is applied on the perturbation. To enhance the accuracy of the numerical calculation results and to confirm the convergence in calculation, we used high-order three-dimensional 20-node piezoelectric coupling elements with reduced integration (C3D20RE). The grid size was 1 mm in the length and width directions and was divided into five equal parts along the depth, resulting in 10,580 hexahedral elements comprising 50,243 nodes. We then compared the out-of-plane and in-plane mode shapes of the piezoelectric materials using FEM with experiment measurements to verify the accuracy of the material constants.

#### **5. Results in Material Constants**

The values required for RM were measured using an impedance analyzer and then converted. With the APC-855 material as an example, we took the average of the results from the five specimens. Table 2 presents the resonant and anti-resonant frequencies measured by the impedance analyzer in the TE mode, which were substituted into Equations (7) through (10) to derive *Kt* , *C<sup>E</sup>* <sup>33</sup>, and *<sup>C</sup><sup>D</sup>* <sup>33</sup>. Table 3 presents the resonant and anti-resonant frequencies measured by the impedance analyzer in the TS mode, which were substituted into Equations (12) through (16) to derive *<sup>K</sup>*15, *<sup>C</sup><sup>E</sup>* <sup>44</sup>, and *CD* <sup>44</sup>. Table 4 shows the capacitance values measured at 1 kHz in the TE and TS modes, which were substituted into Equations (10), (11), (17), and (18) to derive ε*<sup>T</sup>* <sup>33</sup> and <sup>ε</sup>*<sup>T</sup>* <sup>11</sup>. Table 5 displays the resonant and anti-resonant frequencies measured by the impedance analyzer in the LE mode, which were substituted into Equations (19) through (21) to derive *<sup>K</sup>*33, *SE* <sup>33</sup>, and *SD* <sup>33</sup>. Table 6 contains the resonant and anti-resonant frequencies measured by the impedance analyzer in the LTE mode, which were substituted into Equations (22) through (24) to derive *<sup>K</sup>*31, *<sup>S</sup><sup>E</sup>* <sup>11</sup>, and *<sup>S</sup><sup>D</sup>* <sup>11</sup>. Table 7 shows the resonant and anti-resonant frequencies measured by the impedance analyzer in the RAD mode, which were substituted into Equations (25) through (27) to derive *KP* and <sup>σ</sup>*E*.

**Table 2.** Resonant and anti-resonant frequencies and the converted material constants in thickness extension (TE) mode.


**Table 3.** Resonant and anti-resonant frequencies and the converted material constants in thickness shear extension (TS) mode.



**Table 4.** Capacitance and dielectric constant measured at 1 kHz in TE and TS mode.

**Table 5.** Resonant and anti-resonant frequencies and the converted material constants in length extension (LE) mode.


**Table 6.** Resonant and anti-resonant frequencies and the converted material constants in length thickness extension (LTE) mode.


**Table 7.** Resonant and anti-resonant frequencies and the converted material constants in radial extension (RAD) mode.


Table 8 compares the material constants obtained using RM and those provided by the manufacturer for piezoelectric material APC-855 from APC International, Ltd., and Table 9 compares the material constants for piezoelectric material C-2 from the Fuji Ceramics Corporation. Some differences exist between the material constants obtained using RM in this study and those provided by the manufacturer or in literature for the two piezoelectric ceramic materials. Although the difference is not very large in piezoelectric *d*-form parameters, the finite element analysis produces calculating inaccuracy or it cannot obtain the results because of the missing terms in the manufacturer provided one. Due to transformation in calculation, the *e*-form parameters in the manufacturer provided one have a larger difference with those obtained by the resonance method in this study. We therefore also substituted these material constants into FEM to calculate the high-frequency in-plane low-frequency out-of-plane vibration characteristics and performed experimental measurements to determine the accuracy of applying these material constants to dynamic analysis.


**Table 8.** Material constant of APC-855 by resonance method, manufacturer, and Ref [27].

**Table 9.** Material constant of FUJI C-2 by resonance method and manufacturer.


\* Combined calculation results from RM experimental measurements and material parameters provided by manufacturer.

#### **6. Discussion and Verification**

We used ESPI to measure the out-of-plane and in-plane vibration characteristics of piezoelectric materials. The white regions in the images were the nodal lines, and the first lines next to them show the sub-micron vibration displacements [25,26]. The resonant frequencies measured by the impedance analyzer were also listed. Next, we applied the material constants obtained using the mathematical calculations in RM to finite element analysis software (ABAQUS) to derive the out-of-plane and in-plane resonant frequencies and mode shapes. The bold black lines in the mode shape results were the nodal lines, and we compared these numerical calculation results with experimental measurements to verify the reliability of the material constants.

#### *6.1. APC-855*

The vibration verification specimen used in this study was the thin square plate in the TE mode, in which the polarization and the vibrations are oriented along the depth. For the geometric dimensions, we used *l* = 46 mm, *w* = 46 mm, and *t* = 1.2 mm in the FEM calculations and conducted experimental verification.

#### 6.1.1. Out-Of-Plane Vibration

As shown in Figure 4, nine mode shapes were measured in the out-of-plane vibrations, reaching around 10 kHz. "RM" indicates the results of inputting the parameters derived using RM into FEM calculations, and "Ref [27]" indicates the analysis results obtained in literature for APC-855. The simulation results of these two sets of material constants were both consistent with the mode shapes measured in the piezoelectric specimen experiments. However, the RM calculation results for the resonant frequencies were precise, the errors remaining within 6.472%. The errors in Ref [27], on the other hand, were all greater than 10% and increased as the frequency increased. This demonstrates that the material constants obtained using RM in this study have a certain degree of reliability.


**Figure 4.** *Cont.*

**Figure 4.** Out-of-plane vibration modes for APC-855.

#### 6.1.2. In-Plane Vibration

As shown in Figure 5, five mode shapes could be obtained and the in-plane resonant frequency was in the high-frequency range from 20 kHz to 100 kHz. Mode shapes could be obtained via mode shape analysis of the material parameters derived in RM and Ref [27]. A comparison of the resonant frequency measurement results of impedance analysis and the ESPI frequency measurements revealed that they were almost identical. However, we found that the error between the RM analysis results and ESPI measurements increased as the frequency increased and exceeded 10%. In contrast, the analysis results in Ref [27] were relatively accurate. We speculate that this is because the impedance is smallest at resonant frequency vibrations; however, ESPI measurement of in-plane vibrations requires high voltages for excitation, which generates a self-heating effect and material property changes during resonance. We therefore speculate that the material constants derived in Ref [27] are better suited to the prediction of the in-plane vibrations in APC-855 following the self-heating effect. We explain this phenomenon in further detail in the following section.

**Figure 5.** In-plane vibration modes for APC-855, after self heating in resonance.

#### 6.1.3. Influence on Self-Heating in Resonance

For the out-of-plane modes, the piezoceramics was operated at lower frequencies and excited on almost one-hundred voltages. The piezoelectric material remained at almost the same temperature, which was measured by thermography at room temperature. However, we speculate that the in-plane parameters in Ref [27] are more accurate than the material constants obtained using RM in this study because APC-855 is a soft piezoelectric ceramic that is more prone to self-heating at high in-plane resonant frequencies. Thus, we used new specimens to conduct the in-plane vibration measurements. A comparison of Figures 5 and 6 show that the main difference lies in the duration of excitation. In Figure 6, the in-plane vibrations in the experiment were controlled so as to prevent self-heating caused by over-excitation [28,29]. The lower impedance occurs at the resonant frequency. The experimental measurement of in-plane vibrations requires higher excitation voltage, since ESPI

measurement is sub-micrometer resolution. Under high-voltage, high-frequency, and low-impedance excitation conditions, high temperatures are often generated in piezoelectric components, and the material properties of piezoelectric components change at high temperatures, which may alter the capacitance and material rigidity of piezoelectric ceramic. The results in Tables 5 and 6 show that after self-heating, the errors between the resonant frequency measurements of APC-855 and the results of FEM calculations using the material constants obtained with RM become significant. When the voltage and excitation time are controlled under only several seconds to confirm slight temperature rising, the resonant frequencies from the experimental ESPI measurements and impedance analysis compared to FEM are almost identical. The RM results also become accurate, with errors less than 5.314%. This verifies that the parameters measured using RM can also be used to calculate in-plane vibration values with accuracy. The errors of the Ref [27] parameters input into FEM instead become higher than they were before the self-heating of the specimen, ranging from 7% to 11%. From Equations (1) through (4), the electro-thermal-elastic model illustrated that the physical parameters and material constants were influenced by a thermal effect. In the experimental results, it is concluded that material constants and dynamic characteristics of the piezoelectric material depend on temperature. Zhang et al. [22] also show the influence on dielectric loss, dielectric constant, shifting resonant frequency of piezoelectric material through the thermal effect in their results. Hence, the RM used in the measurement of piezoelectric material constants has to be considerate of the temperature dependency in the actual application. This study therefore demonstrated that the material constants obtained using RM can be used to accurately calculate the out-of-plane and in-plane vibration characteristics of piezoelectric materials. However, the impact of self-heating must be taken into account.

**Figure 6.** In-plane vibration mode for APC-855 before self heating in resonance.

#### *6.2. FUJI C-2*

In this paper, we discuss the accuracy of using RM to measure material constants of piezoelectric material APC-855, which is a soft ceramic in which self-heating is more severe during in-plane vibrations. We also examined the accuracy of using RM to measure material constants of piezoelectric material FUJI C-2, which is a hard ceramic. The dimensions of the test specimen were identical to those of the APC-855 specimen.

#### 6.2.1. Out-of-Plane Vibration

As shown in Figure 7, ten mode shapes were obtained in the out-of-plane vibrations up to 12 kHz. We compared the experiment results and the FEM results. "Manufacturer provided" indicates the results of FEM calculations using the material constants provided by the manufacturer. However, not all of the parameters were available from the manufacturer, so we had to use parameters measured using RM to make up for the missing *C<sup>E</sup>* <sup>13</sup>, *<sup>e</sup>*31, *<sup>e</sup>*33, and <sup>ε</sup>*<sup>S</sup>* <sup>33</sup> in the FEM calculations of the vibration characteristics. Comparisons of the out-of-plane mode shape revealed high accuracy in all modes except Mode 5, which could not be derived using the parameters provided by the manufacturer. The RM errors were all less than 3.765%, thereby indicating that the RM results were all highly accurate. In contrast, only the errors of the results from the parameters provided by the manufacturer in Modes 1, 2, 3 and 6 were less than 10%. The errors in the remaining modes were greater than 10%. The manufacturers provided their piezoelectric material constants by measured in standard specimens; however, the quality control in turns of mass production influenced the difference of material properties. We can thus confirm that this study successfully used RM to obtain reliable material constants for piezoelectric ceramic materials.

**Figure 7.** *Cont.*

**Figure 7.** Out-of-plane vibration modes for FUJI C-2.

#### 6.2.2. In-Plane Vibration

As shown in Figure 8, five mode shapes were obtained in the in-plane vibration experiments. The FEM results from RM material constants were all relatively accurate, whereas the parameters provided by the manufacturer could not produce any results at all. Because finite element analysis inputs the *e*-form material constants to calculate the results, it causes the transferred *e*-form material constants to have a larger difference with those that are manufacturer provided, as listed in Table 9. The impedance analysis results were also relatively consistent with the ESPI measurements, with only slight errors in Mode 5. The errors of the FEM calculations from the RM parameters were all within 4%, thereby indicating that the piezoelectric ceramic material constants obtained in this study were relatively reliable.

**Figure 8.** In-plane vibration modes for FUJI C-2.

#### *6.3. Discussion on APC-855 and FUJI C-2*

APC-855 and FUJI C-2 are soft and hard piezoelectric ceramics, respectively, and the applications of the piezoelectric effects of soft and hard piezoelectric ceramics differ. Soft piezoelectric ceramics present stronger direct piezoelectric effects and are thus more suitable for sensors, whereas hard piezoelectric ceramics show more distinct inverse piezoelectric effects and are thus more suitable for actuators. We examined the vibration measurements of these two materials via the drive voltage, and found that lower voltages can be applied to FUJI C-2 than to APC-855 while still producing ESPI measurement results at the sub-micro level. This verifies that actuators comprising hard ceramic piezoelectric materials may have stronger inverse piezoelectric effects. In contrast, soft ceramic piezoelectric materials require greater drive voltages to make the vibration fringes clear enough to observe, but the self-heating issue severely affects the resonant frequency measurement results. Once the temperature reaches half of the Curie Temperature and remains there, thermal deterioration takes place in the material and is accompanied by depolarization [30–33]. The vibration characteristics of the piezoelectric material are also affected, thereby affecting the accuracy of material constants derived using RM. Thus, this study successfully used dynamic characteristic measurement to verify the applicability of material constants obtained using RM.

#### **7. Conclusions**

Due to the anisotropic mechanical properties of piezoelectric materials, their material constants cannot be obtained using general material testing methods. Literature indicates that RM, an approach developed over years of research using structural and polarization characteristics, can effectively measure the anisotropic parameters of piezoelectric materials. This study contains complete details on the specimens and all of the equations that should be prepared, and we used dynamic experiments to verify the accuracy of the material constants that were obtained.

RM measurements require TE, LE, TS, LTE, and RAD mode specimens. Using impedance analysis, their resonant and anti-resonant frequencies were obtained and then substituted into the RM equations. We also used the conversional relationships in constitutive equations to obtain the complete elastic, piezoelectric, and dielectric coefficients of piezoelectric ceramic materials.

We used a non-destructive testing method for dynamic characteristic measurement, verifying FEM analysis results from anisotropic material constants obtained using RM from soft and hard piezoelectric ceramic materials. The errors between the out-of-plane and in-plane mode shapes, which is frequency-dependent on material constants, derived from the RM material constants and those measured in experiments, were very small. We also discussed the impact of the temperature changes caused by self-heating in soft and hard piezoelectric ceramics on the properties of the piezoelectric materials. Whether the piezoelectric ceramic is prone to self-heating affects the accuracy of RM material constants. The dynamic analysis of non-destructive testing revealed that soft ceramic is more prone to self-heating, which shifts the original resonant frequency. This indicates that the question of whether the properties of piezoelectric materials change due to self-heating must be noted, which will in turn affect the original operating frequency and voltage designs used in transducers.

**Author Contributions:** This study was initiated and designed by Y.-H.H., set-up experimental platform and conducted the experiments by C.-Y.Y. and T.-R.H. T.-R.H. were also responsible in numerical simulations. C.-Y.Y. and Y.-H.H. wrote the paper. All authors contributed in analyzing the simulation and experimental results. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Ministry of Science and Technology (Republic of China) under Grant MOST 107-2221-E-002 -193-MY3.

**Acknowledgments:** The authors gratefully acknowledge the financial support of this research by the Ministry of Science and Technology (Republic of China) under Grant MOST 107-2221-E-002 -193-MY3.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Table A1 is listed in the material property provided by APC International, Ltd. Table A2 is listed in the material property provided by Fuji Ceramics Corporation. Table A3 is listed in the Poisson's ratio, σE, correspondent to the ratio of the first two resonant frequencies.



**Table A2.** Material constants of APC-855 piezoceramic provided from Fuji Ceramics Corporation.


**Table A3.** Poisson's ratio corresponds to the ratio of the first two resonant frequencies.



**Table A3.** *Cont.*

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Application Study on Road Surface Monitoring Using DTW Based Image Processing and Ultrasonic Sensors**

#### **Sunil Kumar Sharma 1, Haidang Phan 2,3 and Jaesun Lee 4,\***


Received: 15 May 2020; Accepted: 28 June 2020; Published: 29 June 2020

#### **Featured Application: Road surface monitoring for severe conditions.**

**Abstract:** Road surface monitoring is an essential problem in providing smooth road infrastructure to commuters. This paper proposed an efficient road surface monitoring using an ultrasonic sensor and image processing technique. A novel cost-effective system, which includes ultrasonic sensors sensing with GPS for the detection of the road surface conditions, was designed and proposed. Dynamic time warping (DTW) technique was incorporated with ultrasonic sensors to improve the classification and accuracy of road surface detecting conditions. A new algorithm, HANUMAN, was proposed for automatic recognition and calculation of pothole and speed bumps. Manual inspection was performed and comparison was undertaken to validate the results. The proposed system showed better efficiency than the previous systems with a 95.50% detection rate for various road surface irregularities. The novel framework will not only identify the road irregularities, but also help in decreasing the number of accidents by alerting drivers.

**Keywords:** ultrasonic sensors; GPS; surface monitoring; image processing; dynamic time warping

#### **1. Introduction**

Treacherous road surface conditions are a significant problem for safe and comfortable transportation. The importance of road infrastructure for society can be compared with the importance of blood vessels for humans. For better road surface quality, it should be monitored continuously and repaired as necessary. The optimal distribution of resources for road repairs is possible by providing the availability of comprehensive and objective real-time data about the state of the roads. Developing countries with an escalating economy have an enormously extensive network of roads. Events within the past have shown that countries like the USA have invested about \$68 billion a year on the maintenance of roads, bridges, and highways. These potholes also contribute to the degradation of suspension and fuel efficiency. In India, roads are the most common means of mobility; roads carry 90% of passengers and 65% of freight. Even with such a massive network, traveling is tiring in India as the roads are not sufficiently wide and the maintenance is inadequate. Regardless of where you are, driving is a stressful and potentially life-threatening affair. Over the last two decades, with an exponential rise in population, the number of vehicles has increased at a tremendous rate. The road infrastructure is not at par with the number of vehicles to support it, resulting in traffic congestion and

accidents. As the ECONOMIC TIME's report of January 2018 suggests, at the very least, 400 people lose their lives every day in the country to road accidents due to potholes [1].

According to the report, "Road Accidents in India, 2016 [1] ", by the Ministry of Road Transport and Highways, there were 480,652 road accidents in 2016 as opposed to 501,423 in 2015. However, fatalities resulting from these accidents have risen by about 3.2%, which means that nearly 150,785 people were killed in 2016 compared to 46,133 in 2015 [1]. In India, most of the vehicles (unusually small hatchbacks) sold are either a 0-star rating or a 1-star rating in correspondence to the European standards, which implies that there could be fatal injuries even at moderate speeds. A considerable amount of research has been conducted on reducing traffic congestion and the number of ever-growing accidents. Vast sums of money are spent on repairing these potholes and bumps, however, small developing countries cannot afford to invest these amounts, which are a significant cause of road degradation. Irregularities on the road such as potholes caused by activities like erosion, weather, traffic, and some other factors may lead to accidents. Moreover, human-induced bumps such as speed breakers, which are made in the form of speed bumps, are used to control the speed of the vehicles, but they are not distributed at even distances, and the heights of the speed bumps have not been scientifically justified, which could also cause accidents. While said obstacles ought to be signaled according to specific road regulations, they are not always correctly labeled. One of the primary goals of this research was to detect road irregularities and decrease the number of accidents and fatalities. This approach is cost-efficient and can be applied worldwide. Generically, in modern-day field practices, the evaluation of the road irregularity data is done entirely via the raw data assessment [2]. Several road irregularity stats are used for the same. The existing technique can be categorically classified into three main types: imaging-based, sensor, and manual techniques.

There has been a substantial amount of research done in the past on the detection of potholes [3–5] using various techniques like reconstruction using laser [6], stereo vision [7,8], Light Detection and Raanging (LIDAR) [9], ultrasonic [10–13], etc. Forrest et al. [10] proposed a method of detecting road surface disruptions based on ultrasonic sensors in which the pothole detection rate of 62% was achieved on a paved road. Singh et al. [11] used a global positioning system receiver and ultrasonic for the identification of geographical location coordinates of the detected potholes and speed breaker with 59.23% efficiency. Madli et al. [12] worked on the automatic detection and notification of potholes and humps on roads to aid drivers using ultrasonic and a Global Positioning System (GPS). It was found that the system was 74% efficient and the data could be stored in the cloud. Shivaleelavathi et al. [13] also proposed an intelligent system using an ultrasonic sensor based on Raspberry Pi, GPS, and ultrasonic to take corrective measures. Devapriya et al. [14] proposed a real-time system that detected speed bumps by the analysis of photos of the road that were further enhanced using a Gaussian Filter (computer vision technique). This system had a correct rate ranging from 30% to 92%; however, failure was witnessed for faded speed bumps, and for a higher price, they had to be painted and maintained; this was a significant drawback of the system. Eriksson et al. [8] established the patrol system that comprised of a GPS and a 3-axis accelerometer for the detection of potholes. Likewise, Chen et al. [15] used the very same system and evaluated the power spectral density (PSD) for the detection of pavement irregularities using Fourier transformations. Inaccurate location of the irregularities was a shortcoming of the system. Nericell was a system developed by Mohan et al. [7] that comprised of amici, global system for mobile communication (GSM), accelerometer, and GPS sensors for road surface monitoring. The method yielded a False Negative Rate (FNR) of 37% and 41% for low speed and high-speed bumps, respectively. The use of smartphones, an adaptive fuzzy classifier, GPS, and inertial sensors for the detection of sudden driving events was presented by Arroyo et al. [16] and a precision performance of 0.87 and 0.91 was recorded. A fuzzy inference system was proposed by Aljaafreh et al. [17] for the detection of speed bumps using accelerometers embedded in smartphones. Nevertheless, the actual results of the tests were missing and accuracy only ranged from 70–80%. A multi-mobile phone system was developed by Astarita et al. [18], who also used accelerometers

for the detection of bumps and potholes. A whopping 90% detection rate was obtained for bumps whereas the False Positive Rate (FPR) was about 35%.

Furthermore, the use of a smartphone accelerometer with six different algorithms to categorize road surface conditions was done by González et al. [19], who reported an Area Under Curve (AUC) ranging from 0.823 to 0.9445. Nature driven approaches such as genetic algorithm have also been explored: Yu and Yu [20] explored the use of genetic algorithms using images as a source of information to detect road irregularities. Moreover, abnormalities or distress data collection are increasingly being automated by using various imaging systems. However, analysis of the collected raw video clips for distress assessment is still predominantly being done manually, which is expensive, time-consuming, and slows down road maintenance management [21]. Pre-existing research suggests that three-dimensional image-based techniques are quite expensive; on the other hand, vibration sensor-based techniques are not as accurate and reliable [22] for road surface monitoring. The majority of the existing two-dimensional imaging methods have solely focused on either crack distress detection or pothole detection.

In this paper, an experimental analysis was conducted to detect potholes and speed bumps in quasi-real-time conditions using three different techniques, namely, an ultrasonic sensor with DWT, imaging-based systems, and manual inspection for validation from various roads of Noida city. In the first technique, which consists of an ultrasonic sensor, a GPS and Arduino connected to the Head Up Display (HUD) device are introduced. This system is a cost-effective method for a better, faster, and more reliable way to avoid accidents caused by potholes and raised humps. During the conduct of this experimental analysis, a collection of an exclusive data pattern of speed bumps and potholes was undertaken and directly fed to dynamic time warping (DTW) to find the likeness and similarity. A robust method for automated segmentation of frames with and without distress from the image processing technique is presented. These photo processing methods are responsible for the detection of potholes and speed bumps, which are supported by user-defined decision logic. The derivation of the decision logic is dependent on three characteristic visual properties of these irregularities (i.e., the standard deviation for calculating pixels intensities, circularity (CIRC) for calculating the shape, and average width for the dimension). An agile video separation algorithm called irregularity frame selection (IFS) was used to separate the visual frames with irregularities or not categorically. A new algorithm named the Hanuman algorithm took care of the second phase for automatically detecting and assessing potholes and speed bumps in one go. Finally, in the third technique, a manual visual inspection was undertaken to validate the above two methods.

#### **2. Materials and Methods**

The proposed threshold-based efficient road monitoring system detects the pothole and a speed break. Figure 1 shows the methodology overview of the information taken from the ultrasonic sensor. Collection and analysis of data were done using the primary server where the application of the DTW algorithm and noise filtering techniques was made. After the observation of the limitations of the existing systems, the present framework introduces the use of an ultrasonic sensor, DTW, IFS (algorithms), and filters for the detection of road irregularities. The ultrasonic sensor is fit in the moving direction of the car, as explained in Section 4. Information gathered with the help of sensors is sent to the central server for further analysis. Furthermore, noises, vibrations, and gravity components are eliminated using the filters. DTW helps to obtain the closeness score by analyzing the patterns for potholes and speed bumps. Ergo, road irregularity detection is undertaken with its GPS signature. Eventually, the final output is displayed on the mobile HUD for ease of the driver, making it even more ergonomic.

**Figure 1.** Flow chart of the proposed novel system includes ultrasonic sensor with dynamic time warping (DTW).

#### *2.1. Hardware Setup*

The main aim of the research was to create inexpensive hardware that runs on open software. Ergo, we used an Arduino Uno with an ultrasonic sensor (in Figures 2 and 3) and a GPS. Table 1 shows the hardware specification and Figure 4 shows the components used in the setup. The architecture of the proposed system consisted of four parts: an ultrasonic sensor, Arduino, a server module, and a heads-up display module. The ultrasonic sensor and Arduino are used to gather information about potholes and speed bumps and their geographical locations, and this information is sent to the server. The server module receives information from the microcontroller module, processes, and stores in the database.

**Figure 2.** Working of ultrasonic sensor, (**a**) Plain road, (**b**) Speed bump, and (**c**) pothole.

**Figure 3.** Timing response of the ultrasonic sensor.


**Table 1.** The hardware specification of the proposed framework.

**Figure 4.** Components used in the road surface monitoring system: (**a**) Arduino with HY-SRF05; (**b**) GPS Receiver; (**c**) Head-up display module.

Ultrasonic sensors are utilized to determine the separation between the vehicle frame and the street surface numerically, and this value is passed to the microcontroller. Threshold distance is defined as the distance between the body of the vehicle and the smooth road surface, which depends upon the ground clearance of vehicles. If the value obtained by the ultrasonic sensor is more than the threshold distance, it recognizes it as a pothole, and if it is less, it is recognized as an unimportant irregularity, otherwise the street is smooth.

The microcontroller module (Arduino) is the main controller of the entire system and is responsible for incorporating the ultrasonic sensors and sending and receiving signals between the GPS receiver and the sensors. The microcontroller module is merely the processing power of the whole system and is utilized to accumulate data about potholes and their topographical areas, and these data are passed onto the server. The server module gets this data from the microcontroller and uploads its respective location signatures to the cloud server. This location will appear on the GPS system whenever any vehicle passes through or near it. The Arduino contains several filters to smoothen out the output from the ultrasonic sensor. The basic variations from the ultrasonic sensor are bolstered to the filters, where the information is smoothed and standardized. The filters are to filter the vibration due to some irregularities from the road surface apart from speed bumps and potholes. These vibrations are removed from the ultrasonic sensor information to lessen the impact of other outside powers and for the system to be able to differentiate between a constant/normal variation and an actual pothole or a bump.

#### *2.2. Automated Assessment of Potholes and Speed Bumps Using Imaging Processing*

The inception of the motivation for this research paper was induced directly or indirectly by our personal first-hand experiences with the Indian road conditions. The validation was for the determination of the efficiency of the proposed system. Initially, the pothole and bumps were detected

with the help of using ultrasonic sensors paired with an Arduino attached with a heads-up display feeding the information to the GPS receiver. Irregularity frame selection is a recursive algorithm for searching all the vertices of a graph or tree data structure. This algorithm uses raw visual clips that are further put through several image processing methods; this is done using user-defined logics, which ensures the precise detection of irregularities. Segmentation of the irregular pixels in a visual frame is done by the algorithm using user defined decision logics that categorize those frames by the area covered by the irregular pixel. This segmentation is categorized into two parts: (1) Frames with irregularities, (2) frames without irregularities.

The Hanuman algorithm is used in the presented technique for computerized location and estimation of speed bumps and potholes in a single go from a grouping of video outlines. This algorithm is generated by the arrangement of two visual properties of bumps and potholes, which are distinguished from video photographs of Indian streets. This aggregate arrangement of recognized visual properties incorporates the accompanying, including the following:

**Picture Texture**: The picture surface of the interior of a pothole and speed bump is more different and distinguishing than the encompassing zone (i.e., irregularity less street surface region).

**Shape Factor**: The state of a pothole and speed bump is distinguishable. Their shape factor is estimated as far as circularity, which changes depending on their surface territory and border.

**Measurement**: The element of pothole and speed bump is expressed in terms of average width. This particular algorithm fundamentally includes five segments:


The calculation produced for the robotized identification, estimation, and arrangement of casings with potholes and speed bumps without basic irregularity from a grouping of street video outlines is recorded underneath.

	- (a) Pothole, if STD ≥ 12 & CIRC ≥ 0.14 & W ≥ 70 mm;
	- (b) Speed Bump, if Convexity ≥ 5;

The procedure applied to extract the pothole information using the proposed algorithm is shown in Figure 5. Figure 5a is the original image with potholes; (b) is the binary image after median filtering (i.e., after random noise reduction in column (a) image), and adaptive thresholding (i.e., after object segmentation in column (b) image); (c) is binary image after erosion (i.e., after joining disparate black pixels in column (b) image); column (d) image after the Hanuman algorithm and (e) is the report of the pothole. Similarly, Figure 6 shows the speed bumps.

**Figure 5.** Procedure to categorize the image using the proposed method: (**a**) original image; (**b**) processed image; (**c**) further processed image; (**d**) image after Hanuman algorithm; (**e**) report of the Hanuman algorithm.

**Figure 6.** Procedure to categorize the image using the proposed method: (**a**) original image; (**b**) processed image; (**c**) further processed image; (**d**) image after Hanuman algorithm; (**e**) report of Hanuman algorithm.

#### **3. Experimental Design and Location**

To investigate the road surface, the framework used for the testing is shown in Table 2. The research was performed in the city of Noida, which comes under the National Capital Region of New Delhi. The city has an elevation of 200 m above sea level and a flat topographical region. Like most of the country, there are extremely contrasting roads in Noida, ranging from excellent roads to terrible roads. An expressway just outside of Noida (Figure 7a) and a slum road in Noida (Figure 7b) were considered for the experiment. Several test runs were performed, yielding different results. The experiment was conducted for three months covering around 150–170 km. For the evaluation of the performance of the framework, the algorithm presented in this research was applied in a Windows environment (Dell XPS 15Z with Intel i7 2600, 8GB DDR3 RAM on Windows 10 Pro) with the use of Visual Studio 2017 (VC ++ Redistributable) and Open CV3.3 Library. Hanuman, IFS, and DTW algorithms were used across the entire framework. For test runs, a Honda Amaze was used and image processing was done by a cost effective 4k resolution action camera. Table 2 shows the software and hardware specifications used in the experiment

**Table 2.** Software and hardware specifications used in the experimental setup


**Figure 7.** Field test sites for road surface monitoring: (**a**) Yamuna expressway (from Noida to Agra, India); (**b**) potholes on the main roads (Khora, Noida)

#### **4. Results and Discussion**

The evaluation of the three different techniques is presented while traveling over the R1 (Yamuna expressway) and R2 (Khora Village) roads of Noida city and was also compared with other existing techniques.

#### *4.1. Proposed Novel Framework*

The ultrasonic sensor detects the pothole and speed bumps; these data are saved on the server and used to determine the separation amid the vehicle body and the street surface numerically, and this value is passed by the Arduino. Threshold distance is defined as the distance between the body of the vehicle and the smooth road surface and depends upon on the ground clearance of vehicles. If the value obtained by the ultrasonic sensor is more than the threshold distance, it will recognize it as a pothole, and if it is less, it is a mound, otherwise the street is smooth. The proto test was conducted using the design described in the Methodology Section above, and the results for R1 and R2 are shown in Tables 3 and 4.


**Table 3.** Information of pothole (P) and bumps (B) collected during real-time testing using an ultrasonic sensor with dynamic time warping (DTW) for the R1 test site.

**Table 4.** Information of pothole (P) and bumps (B) collected during real-time testing using ultrasonic sensor with DTW for R2 test site


As can be clearly seen, the uneven roads of the rural areas showed many potholes and bumps with some specific dimensions. On the other hand, the motorway road only showed irregularities up to 6.2 cm, which were filtered by filters such as Speed, Virtual Re-orientation, Filtering Z-pivot, SMA, and Band-Pass on the basis of the threshold limit being 8 cm. Figure 8 shows the experimental comparison in the z-axis displacement for R1 and R2.

**Figure 8.** Experimental comparison of R1 (Yamuna expressway, urban site) and R2 (potholes on the main road (Khora, Noida, rural site)

#### *4.2. Image Processing Result*

Identification of the unique visuals of irregularities like illumination, size, and shape are made using image processing techniques. These accurately identify and help classify those irregularities. The irregularity frame selection algorithm is used as a quick video separation of the real-life visual clips captured of the Indian National Highways at several locations, which are divided into two frame categories (frames with irregularities and frames without irregularities). The Hanuman algorithm is responsible for the evaluation of the database of frames with irregularities, which is an automatic technique for the recognition and analysis of potholes and speedbumps. The proto test was conducted using the design described in the image processing methodology and the results for R1 and R2 are shown in Tables 5 and 6.



**Table 6.** Information of pothole (P) and bumps (B) collected during real-time testing using the image processing method for the R2 test site.


#### *4.3. Manual Inspection*

Manual inspection was performed using a band meter, which is used to measure the depth of an object with precision and accuracy. Figure 9 shows the manual inspection at R1 and R2, respectively. The proto test was conducted using manual inspection, and the results for R1 and R2 are shown in Tables 7 and 8. Generally, manual inspection is performed using some technical equipment, but in our case, we performed a very traditional manual inspection. Different scales were used to measure the different potholes and speed bumps for the validation.

**Figure 9.** Manual inspection on a rough patch in the middle of R1 (**a**) and R2 (**b**).


**Table 7.** Information of pothole (P) and bumps (B) collected by manual inspection for the R1 test site.

**Table 8.** Information of pothole (P) and bumps (B) collected by manual inspection for the R2 test site.


#### *4.4. Validation of Proposed System with Other Techniques*

The same experiment was conducted using three different methods several times on the very same roads; the rural road of Khora village and the urban road of the Yamuna Expressway. Each method, ultrasonic sensor, image processing, and manual inspection had slight deviations, as shown below. The comparison between three different techniques for depth/height, latitude, and longitude with their error are given in Tables 9–14 for both roads.


**Table 9.** Comparison of depth/height between the proposed system (PS), image processing (IP), and manual inspection (MI) for the R1 test site.

**Table 10.** Comparison of latitude between the proposed system (PS), image processing (IP), and manual inspection (MI) for the R1 test site.


**Table 11.** Comparison of longitude between the proposed system (PS), image processing (IP), and manual inspection (MI) for the R1 test site.


**Table 12.** Comparison of depth/height between the proposed system (PS), image processing (IP), and manual inspection (MI) for the R2 test site.


The deviation of the location of the pothole is relevant in different methods. The ultrasonic sensor and image processing detection deviated slightly from the actual position. Ergo, to calculate the accuracy of the system, the mean of the output using different method was calculated. The maximum error was 4.60% and the efficiency of the proposed system was about 95.48%.

Table 15 shows a full comparison between the results obtained with the proposed methodology and those yielded by similar approaches. As can be seen, our approach clearly outperformed the others, where TPR is the True Positive Rate, FPR is the False Positive Rate, and FNR is the False Negative Rate. One of the strongest points presented in this research is the use of a multi-variate model for pothole and speed bump detection. The presented methodology builds a model by combining different methods for the detection of road surface irregularities which provides greater accuracy to the entire framework. The collaboration of different methods with our novel framework targets a much bigger audience, and since it is very cost effective, it makes this system ideal for road surface monitoring in all regions.

**Table 13.** Comparison of latitude between the proposed system (PS), Image processing (IP), and Manual inspection (MI) for the R2 test site.


**Table 14.** Comparison of longitude between the proposed system (PS), image processing (IP), and manual inspection (MI) for the R2 test site.


**Table 15.** Detection rate comparison of different approaches on previous research.


#### **5. Conclusions**

It has been highlighted that the current road anomaly detection frameworks used by the government or the existing authorities are very tedious and overwhelmingly time consuming, which is why this device eradicates all the problems and complications with the existing systems with much better efficiency. The analysis results demonstrated the use of ultrasonic sensors mated with an Arduino, GPS receivers, and HUD module, which runs on several filters. Dynamic time warping was also one of the algorithms used, which makes the validation process more manageable and is also one of the reasons for the increased efficiency. The presented Hanuman algorithm could detect the potholes and speed bumps accurately and efficiently. The calculated accuracy of the system was 95.48%, and the use of inexpensive components makes it the most cost effective framework among the existing systems. The model proposed in this paper serves two critical purposes (i.e., automatic detection of potholes, bumps, and alerting vehicle drivers to evade potential accidents). The HUD application provides real-time alerts about potholes and bumps right in front onto the windshield, making it more ergonomic. The framework works perfectly in all seasons, especially in monsoon as the signal/warning is generated from a database. The novel framework ha the highest efficiency when compared to all preexisting systems. Since it is a cost effective long time solution for pothole management, it can be easily provided to the maximum number of people; which can furthermore save many lives. Integration of GMaps and Satnav with the system can be done for ease of use. This framework can be applied to different types of road conditions and will be the future scope of this work.

**Author Contributions:** S.K.S. prepared the inspection system for the experiment and performed the experiments of this research. H.P. and J.L. conceived the original idea. J.L. designed the methodology and gave guidance and helped to improve the quality of the manuscript. S.K.S. wrote the original draft preparation in consultation with H.P. J.L reviewed and edited the final manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2019R1A5A808320111, No. NRF-2019R1G1A1004577).

**Acknowledgments:** The authors would like to thank Amrit Seth, Ansh Mittal, and Vansh Malik, students at Amity University, who helped in conducting this research.

**Conflicts of Interest:** The authors declare that they have no conflict of interest. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

#### **List of Abbreviations**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Defects Inspection in Wires by Nonlinear Ultrasonic-Guided Wave Generated by Electromagnetic Sensors**

#### **Junpil Park 1, Jaesun Lee 2, Junki Min <sup>3</sup> and Younho Cho 1,\***


Received: 26 May 2020; Accepted: 25 June 2020; Published: 28 June 2020

#### **Featured Application: The second harmonic measurement method using EMAT proposed in this paper has established a system that can detect microdefects in wires without destroying objects before wire processing.**

**Abstract:** Steel wires are widely used as raw materials for spring valves in engines. Considering the quality and safety issues of their structure, there is a demand to develop nondestructive inspection approaches to detect initial damages in steel. In this study, nonlinear ultrasonic-guided waves generated by an electromagnetic acoustic transducer (EMAT) were used to inspect the defects in steel wires. As one of the noncontact testing methods, the use of EMAT has significant advantages to decrease the nonlinearity induced by instruments and transducer contact condition. The principles of design and manufacturing of EMAT are first introduced. The fundamental theory of nonlinear guided waves is also briefly discussed in this investigation. Phase-matched guided wave modes were generated and measured by using EMAT. Variations of acoustic nonlinearity corresponding to existing defects in specimens were obtained. A scanning electron microscope (SEM) was used to check the existence of microdefects in specimen. The results indicate that the use of EMAT can be an effective means to generate and measure nonlinear ultrasonic-guided waves for inspection of microdefects.

**Keywords:** steel wire rod; nonlinear ultrasonic; EMAT; guided wave

#### **1. Introduction**

Steel wire rods are used as raw materials for various mechanical components. This research can be applied to the steel wires, raw material of valve springs and have an influence on keeping smooth cycle motion in a gasoline engine. Under periodic fatigue state, failure of spring valves may occur before expectation of fatigue failure. Therefore, we need a reliable method for steel wires that can evaluate microdefects or early stage damage detection. Up to now, for the research work related with this study not so many articles were published. Gao et al. developed nondestructive evaluation (NDE) techniques using magnetic flux for macro defect [1].

The geometric shape of steel wires before processing of valve springs is cylindrical such as rod and solid shaft. Considering shape characteristics, guided waves have been applied to efficiently inspect long specimens such as wire rods. Research on guided ultrasound has been going on for a long time, but it has not progressed actively due to the complexity of theory and the variety of wave modes. Gazis performed basic guided ultrasonic analysis of cylindrical structures [2,3]. Based on this, the test for guided ultrasound was started by Rose and Cawley's research works [4,5].

The nonlinear technique has the advantage that it can examine the microstructural changes in the nonlinear elastic section. It is based on information of distance, time and velocity of ultrasonic wave used in linear technique. In addition, the acquired ultrasonic signal is further analyzed in the frequency domain. The acquisition method of this technique is as follows: Information on the defect can be obtained by comparing the change at the second or third harmonic frequency. Through this, the internal microdefect can be diagnosed early. It is also possible to predict the life expectancy based on quantitative results. Buck et al. [6] conducted a study on the change of function of the stress and the second harmonic in which the transition occurs in the material. Subsequently, Jacobs and colleagues [7–9] studied the plasticization of materials using nonlinear guided ultrasound. Furthermore, according to using acoustic nonlinear guided wave technique, phase-matching mode is required to generate higher order harmonic amplitude. The basic concept of this research, concerning which factor is related with phase velocity and group velocity, has been improved by M. Deng [10,11].

EMAT is a noncontact technique. Research on EMAT has been conducted by M. Hiaro and H. Ogi of Osaka University [12–14]. Quantitative changes of physical nonlinearity should be measured by ultrasonic nonlinear method. Therefore, there is a need to continuously carry out experiments under the same conditions. That is, the same experimental apparatus and conditions are required to minimize the influence of the nonlinearity of the system. At this time depending on the contact conditions of the excitation and receiving transducers the reliability of this experiment may be greatly affected. This may change the nonlinearity. Therefore, to overcome the nonlinearity, a noncontact EMAT was applied.

Currently, research on microdefects in the material is proceeding smoothly by combining the contact method and the nonlinear method. However, the contact test method using piezoelectric transducer (PZT) requires a certain size of contact area, but it is impossible to test using PZT in a circular form in the case of the wire rod studied in this study. In addition since the PZT patch requires a contact medium such as an ultrasonic couplant, the signal is irregular depending on the amount or contact strength of the contact medium, so reproducibility is not constant when performing repeated experiments. Therefore, research using noncontact technology is necessary, but the result of the research is insufficient. In particular, there have been no reports on the application of the noncontact technique to the internal microdefects of steel wire rods. Therefore, this study applies the nonlinear method, guided wave and noncontact method.

Steel wire rod, which is the raw material of the engine piston valve spring, is typically manufactured from three materials. The nonlinear tendency according to the property of each material was analyzed to secure the suspected defective area. Flexible, torsional and longitudinal modes exist for guided waves in cylindrical structures such as piping and wire rods. Longitudinal mode is applied to ensure smooth secondary harmonics. In addition, we tried to obtain the reliability of the experimental results by setting the same geometric conditions. The nonlinear tendency of each specimen can be grasped through experiments. According to the result of the experiment, the suspected defect area was secured. This study provides a technical and academic understanding of nonlinear guided wave. By using noncontact ultrasonic technique rather than contact acoustic nonlinear guided waves, we can study for the purpose of precise diagnosis.

#### **2. Nonlinear Guided Waves**

#### *2.1. Nonlinear Ultrasonic Test*

The method for detecting nonlinearity in a material is to send and receive induced ultrasonic wave at a specific distance from the material. After the nonlinearity of the material propagates at a certain distance, the underlying frequency wave is distorted to produce a higher harmonic [15–21]. For guided wave, the acoustic nonlinearity can be written as shown in Equation (1).

$$
\beta = \frac{8}{K^2 a} \frac{A\_2}{A\_1^{-2}} f \tag{1}
$$

where *a* is the wave propagation distance and *f* is a frequency independent function. In addition, the ratio *A*2/*A*<sup>2</sup> <sup>1</sup> with the relative acoustic nonlinear parameter β applied in this research is considered as a measure of acoustic nonlinearity, where *A*<sup>1</sup> and *A*<sup>2</sup> are the measured amplitudes of the fundamental and second harmonics.

$$
\beta' = \frac{A\_2}{A\_1} \propto \beta a \tag{2}
$$

The relative acoustic nonlinear parameter increases linearly as a function of propagation distance. The cumulative effect of the second harmonic amplitude is a great advantage for detection in experiment. It is found that second harmonic amplitude linearly grows with propagation distance, when there are synchronism and nonzero power flux. If the wave mode chosen satisfies these two conditions, the second harmonic amplitude will be cumulative. Even a series of double frequency wave components will be generated by the driving sources of nonlinearity. In practice, interest is focused on the second harmonic generation with the cumulative effect, since the cumulative second harmonic wave plays a dominant role after second harmonics propagating some distance.

Figure 1a,b shows the typical waveform of the received signal and its spectrum. The received time domain signal is processed in the frequency domain with the fast Fourier transform (FFT) to obtain its spectrum. Figure 1b shows the fundamental frequency and the second harmonic.

**Figure 1.** Typical received signal. (**a**) Time-domain waveform; (**b**) Fourier spectra of the fundamental and second harmonic signal.

As indicated in Equations (1) and (2), the normalized nonlinear parameter can be calculated with Equation (2) by using the amplitudes of the fundamental and the second harmonic modes.

Figure 2 shows the nonlinear parameters according to the propagation distance. The cumulative effect means that the nonlinearity increases with distance under certain frequency conditions. The double frequency mode which satisfies a specific condition such as the expression shown in the diagram is selected based on the frequency.

**Figure 2.** Nonlinear parameters according to distance.

#### *2.2. EMAT*

EMAT is composed of a generating coil for generating an eddy current of the electromagnet and the transmission and a receiving detecting coil. When the probe is approached to the metal specimen, the inside of the specimen is affected by the electromagnet. As Shown in Figure 3, the eddy current is formed on the surface of the specimen by an alternating current flowing in the generating coil inside the probe. Lorentz force is generated between the eddy current and the low power line, and mechanical displacement occurs on the surface of the specimen. This mechanical displacement is called ultrasound.

**Figure 3.** Principle of electromagnetic acoustic transducer (EMAT) with magnetic flux.

The purpose of this study was to investigate the nonlinearity of materials by inducing an ultrasonic wave in the longitudinal mode in steel wire rod. The permanent magnet was magnetically biased in the axial direction on the specimen. Then, when the coil was wound in the circumferential direction and the current was applied, the longitudinal mode of the guided ultrasonic wave was generated in the cylindrical structure such as the pipe or the wire rod.

EMAT was generated by Lorentz force effect in meander coil. The alternative current induce in the coil and the magnetic field is built in coil by high-power permanent magnet. Here **J** is eddy currents field and B is the magnetic field. The Lorentz force is generated by the eddy currents and magnetic field

$$
\overrightarrow{F} = \overrightarrow{J} \times \overrightarrow{B} \tag{3}
$$

As shown in Figure 4, the geometry used for describing the reception of EMAT is considering various factors. Equation (3) explain the relationship between magnetic field and eddy current. The influence of magnetic field for inducing eddy current is indicated by Equations (4) and (5). Produce a voltage pick up by receiving coil is Equation (6). EMAT does not require a contact medium in flaw detection. As a noncontact method using electromagnetic, it shows higher performance than

PZT for measuring nonlinearity of materials. However, it should be considered that the receiving sensitivity is lower than PZT.

$$\begin{array}{ccc} E\_{\text{cd}} = \nu \times B & \xrightarrow{\cdot} & \xrightarrow{\cdot} & \xrightarrow{\cdot} & \xrightarrow{\cdot} \\ J\_{\text{ed}} = \sigma \times E\_{\text{ed}} \end{array} \tag{4}$$

$$\begin{aligned} \nabla^2 \mathbf{A}\_0 &= j\omega \sigma\_0 \mu\_0 \mathbf{A}\_0\\ \nabla^2 \mathbf{A}\_1 &= -\omega^2 \sigma\_1 \mu\_1 \mathbf{A}\_1\\ \nabla^2 \mathbf{A}\_2 &= j\omega \sigma\_2 \mu\_2 \mathbf{A}\_2\\ \nabla^2 \mathbf{A}\_3 &= j\omega \sigma\_3 \mu\_3 \mathbf{A}\_3 \end{aligned} \tag{5}$$

$$\begin{array}{ll} \mathbf{E}\_r = -\frac{\partial \mathbf{A}}{\partial t} & \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \\ \end{array} \\ V\_r = N \times l \times \mathbf{E}\_r \end{array} & \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \\ \end{array} \\ \end{array} \\ \end{array} = \begin{array}{c} \begin{array}{c} \begin{array}{c} \end{array} \\ \end{array} \\ \end{array} \tag{6}$$

where ν is the velocity of particle vibration, B is the static magnetic field, σ is the conductivity of the test piece, E*ed* and **J***ed* are the induced electric field and the induced eddy current sheet due to the particle vibration at the depth *zed* respectively.

**Figure 4.** Geometry of steel wire and EMAT.

#### **3. Experiments**

Three different kinds of specimens with total 9 samples were tested in this investigation. As shown in Figure 5, the size information of the three specimens was the same, but there were differences in material properties and tensile strength. Detailed information of the samples can be found in in Table 1.

**Figure 5.** Specimens of steel wire rods.


**Table 1.** Material properties of specimens.

The frequency and modes were selected by frequency multiply wire thickness (*f* ×*d*). Phase matching must be performed through the dispersion curve. The conditions for generation of second harmonic in a cylindrical structure such as wire rod and pipe are as follows [13].

The power flux at the first harmonic should not be 0 in the second harmonic frequency range. (*fn sur f* + *fn vol* -0)

Phase and group velocities of the first and second harmonics should be the same through phase matching (*kn* ∗ = 2*k*).

Figure 6 shows the dispersion diagram of steel wire rod. It can be confirmed that the phase velocity and the group velocity were the same at 4 MHz · mm, which was a double frequency of 2 MHz · mm. Through this phase-matching, we could use a nonlinear guided wave technique. L1 (f·d = 2 MHz·mm) and L2 (f·d = 4 MHz·mm) were selected because the resulted of the L1 and L2 modes were clearer for the second harmonic component than the other modes.

**Figure 6.** Dispersion curve of steel wire rod.

EMAT was designed and manufactured through impedance and LC-matching. Impedance matching means to eliminate the phase difference by minimizing reflection or loss of input power by buffering impedance difference at input and output. Selecting the appropriate inductor and capacitor element values through Equation (7) and choosing the resonance frequency is called LC-matching.

$$Z = R + jwL + \frac{1}{jwC}\overrightarrow{F} = \overrightarrow{J} \times \overrightarrow{B} \tag{7}$$

The frequency at which the imaginary part becomes 0 in the Equation (7) through matching between the inductor (L) and the capacitor (C) is called a resonance frequency. Figure 7 shows the resonance frequency due to LC-matching of the capacitor and the inductor.

**Figure 7.** Resonant frequency setting by LC-matching of capacitor and inductor.

In order to realize the inductance in the ultrasonic generating coil, a long meander coil was used in a limited space. A resonant circuit was constructed to maximize current flow in the meander coil. Figure 8 shows a meander coil with a wavelength of 0.38 mm.

**Figure 8.** Meander coils with a wavelength of 0.38 mm.

Figure 9 shows the impedance-matching EMAT transducer used in this study. The EMAT was manufactured at 1 MHz as shown in Figure 10 and excited at 0.625 MHz. The wavelength (λ = 5.82 mm) of 1 MHz can be derived by substituting the desired frequency and longitudinal wave velocity (*C<sup>L</sup> Steel* = 5.82 mm/μs) in the steel into Equation (8).

$$
\lambda = \frac{c}{f} \stackrel{\rightarrow}{F} = \stackrel{\rightarrow}{J} \times \stackrel{\rightarrow}{B} \tag{8}
$$

**Figure 9.** Impedance-matched EMAT transducer.

**Figure 10.** Impedance matching.

Practically, 5.82 mm cannot be applied to steel wire with a diameter of 3.2 mm. It can be applied with lower wavelength as λ/2, λ/4, λ/8. Since 0.38 mm is an approximate value of 1/16 of wavelength, EMAT was fabricated using 0.38-mm meander coil.

Experiments were conducted to investigate the behavior of nonlinear induced ultrasound. Figure 11 shows the experimental setup. The RITEC RAM-4000 was used as a tone-burst ultrasonic transmitter and receiver and it was processed by a pitch–catch method.

**Figure 11.** Experimental setup for using EMAT.

In order to obtain the reliability of the experimental results, the same geometric conditions were set, and two kinds of experiments were conducted. As shown in Figure 12, a meander line and permanent magnets were installed on the steel wire rod. As Shown in Figure 13a, the first experiment was divided into four sections of 60 mm each, and the relative nonlinearity of each wire was confirmed. The second experiment was conducted by increasing the propagation distance by 60 mm as shown in Figure 13b. The schematic diagram of each experiment is shown in Figure 13. Experiments were conducted to determine whether defects could be detected through the relative nonlinearity.

**Figure 12.** Setup of the meander line and permanent magnet on steel wire rod.

**Figure 13.** Schematic of experimental: (**a**) Experimental method with each detection area; (**b**) experimental method with increasing distance.

#### **4. Results and Discussion**

The relative nonlinearity of the fabricated EMAT sensor was measured for each specimen in the longitudinal direction of the steel wire. The relative nonlinearity tendency was studied. As mentioned above, the nonlinear guided wave technique is a technique to diagnose changes in material properties and microdefects by using harmonic components unlike the conventional technique. EMAT shows higher performance than PZT to measure nonlinearity of material because it does not need contact medium and uses electromagnetic in flaw detection. In case of using PZT, the error due to the individual difference of the experimenter cannot be ignored. Therefore, the method using EMAT is very appropriate. Figure 14 shows the signals for linear ultrasonic wave in three specimens. The time domain signal for each specimen can distinguish between the distance and the speed relative to the time, however, microdefects cannot be confirmed. In addition, the analysis of the modes showed that the L1 mode was well implemented for SWOSC-V and SWOCS-VHV, but SWOSC-VHS was able to confirm that the various modes were mixed. Therefore, the experiment was conducted using a nonlinear guided ultrasonic technique using second harmonic.

**Figure 14.** Time domain signal. (**a**) SWOSC-VHS #1; (**b**) SWOSC-VHV #2; (**c**) SWOSC-VHS #3.

The fundamental amplitude A1 and second harmonic amplitude A2 amplitude values can be derived by measuring the time domain signal of 60 mm intervals from three specimen for each material and then performing FFT. Substituting the derived A1 and A2 into Equation (2), nonlinear parameters can be obtained for each specimen. Figure 15 is the result of comparing the relative nonlinearity in each zone. In all three specimens of SWOSC-V, the relative nonlinearity increases with increasing propagation distance. The relative nonlinearity of all three specimens decreases at the propagation distance of 181–240 mm, which can be regarded as the material properties of the specimen. SWOSC-VHV also shows the same tendency of three specimens. However, SWOSC-VHS # 2 has a tendency to be different from the # 1 and # 3 specimens at 61–180 mm propagation distance. Three specimens with the same physical properties should show the same tendency toward in each section. However, since VWOSC-VHS #2 has a different tendency from the VWOSC-VHS #1 and #3 samples, it can be determined that there is some abnormality inside the material.

**Figure 15.** Relative nonlinearity for each detection area (**a**) Results of the three-specimen test in SWOSC-V (**b**) results of the three-specimen test in SWOSC-VHV (**c**) results of the three-specimen test in SWOSC-VHS.

Because of the tendency to suspect microscopic defects, additional experiments were conducted to ensure the reliability of the experimental data. Figure 16 is the relative nonlinearity with increasing distance, which was an additional experimental result. SWOSC-V and SWOSC-VHS could obtain a tendency for relative nonlinearity change as the distance increased. However, # 3 of SWOSC-VHV showed a large difference in relative nonlinear factor values from # 1 and # 2 at a propagation distance of 120 mm, but it was very difficult to be considered as a suspected defect region due to the agreement of the tendency for each specimen. At the propagation distance of 180 mm, SWOSC-VHH # 2 had the same tendency as the previous experimental results. A total of two experiments were used to consider suspected microdefects in SWOSC-VHV # 2 at a propagation distance of 180 mm.

**Figure 16.** Relative nonlinearity variation of considering propagation distance (**a**) Results of the three-specimen test in SWOSC-V (**b**) results of the three-specimen test in SWOSC-VHV (**c**) results of the three-specimen test in SWOSC-VHS.

Figure 17 is cross-sectional photographs of microdefect suspected areas and defective free area. Figure 17b shows a clear difference in the distribution of microdefects compared to Figure 17a. This is consistent with the experimental results and SEM photographs. Therefore, it can be seen that the reliability of the experimental result was secured. Figure 18 shows the result of enlarged the micro defect.

(**a**) (**b**)

**Figure 17.** Comparison of SEM result for steel wire. (**a**) SWOSC-VHS #1: non-defect area; (**b**) SWOSC-VHS #2: defect area.

**Figure 18.** Enlarged microdefect.

The experimental results were not quantitative because the distance between the coils and the specimen changed finely when the reinstallation was repeated after separating the meander coil and the specimen from each other during the experiment. In other words, the lift-off changed. In order to solve this problem, we had to build a proper system to fix the meander coil and the specimen firmly. In addition, since the nonlinearity was very sensitive, the ultrasound path difference may have occurred due to the minute scratches on the surface of the specimen, resulting in energy loss and affecting the damping coefficient. Therefore, it is necessary to acquire precision machined specimens.

#### **5. Conclusions**

In this study, the nonlinearity tendency of microdefects was investigated by applying guided wave and nonlinear method to steel wire using EMAT. The relative nonlinearity of each specimen is specific for each material. This is considered to be a characteristic of the content of elements constituting the material. Experimental results show that specimens with large errors and specimens without specimens exist. This is considered to be an error due to the distance between the meander coil and the specimen, that is, the lift-off, when the specimen and the coil are separated and remounted. SEM images were taken to verify the reliability of the experimental results. The results were very satisfactory. The possibility of diagnosing microdefects in steel wire production line was confirmed. In order to utilize the nonlinear technique, there is a need to acquire defect-free specimens. We have identified the need for a system that can acquire clear signals. However, there were experimental errors, there is no critical dispute in confirming the tendency for each specimen. It is possible to observe relatively the suspicion of microdefect. It is expected that quantitative evaluation of each specimen will be possible if the parameters of nonlinear factor are removed through a system which can fix the defect-free specimen and lift-off later. In addition, if this research technique is used, it is possible to quickly check noncontact for defects of steel wire rods passing through the inspection line, and to detect fine defects early before processing the steel wire rod to enhance the safety of the material. It was suitable for inspecting the initial fine defects in steel wires before spring processing by utilizing EMAT, a noncontact method without contact medium. The noncontact method is a method capable of minimizing the deformation of the signal that changes depending on the amount or compression state of the contact medium. If the technique is applied to the industrial site, it is judged that it will be possible to test for microdefects quickly and consistently.

**Author Contributions:** J.P. and J.M. designed and performed the experiments. J.L. and Y.C. conceived the original idea. J.P. developed the theory and performed the computations with support from J.L. and Y.C. verified the analytical methods. J.M. and J.L. wrote the draft study, J.P. completed the final study and Y.C. made the final review. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Korean Evaluation Institute of Industrial Technology (KEIT) grant funded by the Korean government (MOTIE) (NO.10085576).

**Acknowledgments:** This research was carried out with the support of materials from Daewon Kang Up Co., Ltd. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **High-Precision Noncontact Guided Wave Tomographic Imaging of Plate Structures Using a DHB Algorithm**

#### **Junpil Park 1, Jaesun Lee 2, Zong Le <sup>3</sup> and Younho Cho 1,\***


Received: 26 May 2020; Accepted: 24 June 2020; Published: 25 June 2020

**Abstract:** The safety diagnostic inspection of large plate structures, such as nuclear power plant containment liner plates and aircraft wings, is an important issue directly related to the safety of life. This research intends to present a more quantitative defect imaging in the structural health monitoring (SHM) technique by using a wide range of diagnostic techniques using guided ultrasound. A noncontact detection system was applied to compensate for such difficulties because direct access inspection is not possible for high-temperature and massive areas such as nuclear power plants and aircraft. Noncontact systems use unstable pulse laser and air-coupled transducers. Automatic detection systems were built to increase inspection speed and precision and the signal was measured. In addition, a new Difference Hilbert Back Projection (DHB) algorithm that can replace the reconstruction algorithm for the probabilistic inspection of damage (RAPID) algorithm used for imaging defects has been successfully applied to quantitative imaging of plate structure defects. Using an automatic detection system, the precision and detection efficiency of data collection has been greatly improved, and the same results can be obtained by reducing errors in experimental conditions that can occur in repeated experiments. Defects were made in two specimens, and comparative analysis was performed to see if each algorithm can quantitatively represent defects in multiple defects. The new DHB algorithm presented the possibility of observing and predicting the growth direction of defects through the continuous monitoring system.

**Keywords:** tomographic imaging; DHB algorithm; noncontact; air-coupled transducer

#### **1. Introduction**

It is well known that Lamb waves can propagate long distances with little energy loss. Lamb waves have been shown to be sensitive to most types of defects. They are also efficient for detecting tiny subsurface damages and inspecting large areas. Such noncontact nondestructive testing (NDT) methods have great advantages for online inspection because high temperature or other negative factors can restrict the application of contact methods. In earlier research, ultrasonic Rayleigh and Lamb waves were used to reconstruct images in metal plates using water immersion techniques [1]. Other studies have used electromagnetic acoustic transducer (EMAT)-based methods [2–4] and contact techniques [5–9]. Air-coupled transducers have been used for tomography imaging since the end of the last century [10]. The application of purely air-coupled transducers has always been limited for the detection of nonmetallic materials [11,12] due to the high acoustic impedance between air and solid materials. Other methods, such as disc sensors, have also been widely used for the detection of metallic materials [13], but those methods are not truly noncontact methods. The research team of

Rymantas J. K. developed PMN-PT synthesized with lead magnesium niobate–lead titanate crystals, and applied PMN-PT, which has a strong piezoelectric effect, to an air-coupled transducer to improve performance [14,15]. Therefore, it is believed that the hybrid application of a pulse laser and air-coupled transducer will have good performance as a fully noncontact detection system [16,17].

Most of the works mentioned above are based on the traditional filtered back projection (FBP) algorithm. The FBP algorithm is mainly used for radiographic inspection using X-ray [18,19]. Although high resolution is shown in the radiographic inspection, there is a disadvantage that image artifacts are generated when ultrasonic tests are applied. Other works have been based on the reconstruction algorithm for the probabilistic inspection of damage (RAPID) [20–23]. The RAPID algorithm is actively used in ultrasonic tomographic techniques as a technique for imaging defects using a signal difference coefficient. However, the RAPID algorithm has the advantage of being able to visualize faults with small transmission and reception points, but due to some negative influences such as shape factor (factor of influence), there is a limit to quantifying the defects. It is also not suitable for multidefect imaging.

In this study, in order to overcome the limitations of the resolution of the RAPID algorithm, which is typical of the existing ultrasonic tomographic algorithm, the FBP algorithm, an imaging technique used in the radiography technique, was applied to the ultrasonic technique. We confirmed the problem with the image artifacts generated in the FPB algorithm using ultrasonic waves. Difference of the projection data, Hilbert transform of the difference data, distance weight, and back projection were applied to the FBP algorithm to reconstruct the Difference Hilbert Back Projection (DHB) algorithm. The radiographic imaging techniques require many data points and reliable data for high-resolution imaging. However, manual operation is still being widely used for data acquisition. The work efficiency of manual operation is very low, and the precision of signal acquisition cannot be guaranteed. Although automatic control was used in some studies, many of the devices can only be used in some special circumstances. In this work, a laser switch has been adopted in an automatic detection system to improve the precision of signal acquisition. The automatic detection system was developed based on the hybrid application of a pulse laser and an air-coupled transducer. This is one of the best combinations among many different kinds of noncontact NDT methods. The detection system consists of two main parts. One part is an automatic scanning robot arm for the air-coupled transducer. Another part is an automatic clamping device for different kinds of specimens to be tested. The tomography imaging can be conducted efficiently by using this automatic detection system. High-precision signal acquisition can be guaranteed with the help of the laser switch. The iterative reconstruction method used in this work is the fan-beam method. The practicability of this detection system has been confirmed by experimental results. The experiment results matched well with both the theoretical simulation and the actual defects. The imaging quality has been improved to a great extent compared with other methods used in earlier studies.

#### **2. Tomography Theorem**

Tomography imaging was used to obtain the section view of the specimen under testing. Image reconstruction from the data projection is a mathematical problem. The original object can be copied from a large number of projection signals. The results of computed tomography are the distribution image of the linear attenuation coefficients.

The DHB algorithm used in this work has evolved from the traditional FBP algorithm. The experiment results show that the imaging quality obtained from the DHB algorithm is much better than that obtained from other analogous algorithms such as RAPID or FBP. These two frequently used algorithms are also discussed to emphasize the advantages of the new DHB algorithm. The theoretical foundation and the corresponding experimental results will be illustrated in the following part of this paper.

The RAPID is based on a damage index known as the signal difference coefficient (SDC). This SDC can be taken as a measurement of how statistically different the defect signals are from a baseline nondefect signal. Damage indexes for all transducer pairs are spatially distributed. Those damage indexes are summed to generate an image of the damage area. The first step of this algorithm is to calculate SDC between the current signal *xij*(*t*) and the reference *yij*(*t*). The index *i* is for the transmitter and *j* is for the receiver. Mathematically, the SDC can be stated as:

$$SDC\_{ij} = 1 - \left| \frac{\int\_{t\_0}^{t\_0 + \Delta T} [x\_{ij}(t) - \mu\_x][y\_{ij}(t) - \mu\_y]^2 dt}{\sqrt{\int\_{t\_0}^{t\_0 + \Delta T} [x\_{ij}(t) - \mu\_x]^2 dt \int\_{t\_0}^{t\_0 + \Delta T} [y\_{ij}(t) - \mu\_y]^2 dt}} \right|, \tag{1}$$

where *t*<sup>0</sup> is the direct arrival time for each transducer pair, μ*<sup>x</sup>* and μ*<sup>y</sup>* are the mean values of the corresponding current signal and baseline signal, and Δ*T* is the time window. The next step after the SDC values for all sensor pairs are calculated is image reconstruction. An image is generated by distributing each SDC value spatially on the image plane in an elliptical pattern. The two foci of the ellipses are the two transducers' locations. A parameter β controls the size of the ellipse. The amplitude linearly tapers from its maximum value along the line connecting the ellipse foci to zero on the periphery of the ellipse. The parameter β is a shape factor that controls the size of the elliptical distribution. Its value can be arbitrary, but greater than 1.0. Finally, the image amplitude at each pixel is a linear summation of the location probabilities from each transmitter–receiver pair (M). In Equation (2), the coordinate position of the image is displayed. *xRK*, *yRK* is the coordinate position of the sensor with *xRK*, *yRK* is the coordinate position of the receiving sensor. *PK*(*x*, *y*) is the image value of image obtained from a pair of sensors and sensors. *P*(*x*, *y*) is the sum of the image values obtained from the sensors of all pairs (M) and the receiving sensor. *R*(x, y, x*TK*, *yTK*, *xRK*, *yRK*) is a formula showing the relationship between the coordinate position of the image and the position of the sensor. The larger the *R*(x, y, x*TK*, *yTK*, *xRK*, *yRK*), the farther away from the shortest distance of the pair of sensors. β is a shape factor that limits the coverage of *R*(x, y, x*TK*, *yTK*, *xRK*, *yRK*).

$$P(\mathbf{x}, y) = \sum\_{i=-1}^{N-1} \mathbf{P}\_k(\mathbf{x}, y) \; \;= \sum\_{K=-1}^{M} \frac{\text{SDC}}{\beta - 1} \|\beta - R(\mathbf{x}, y, \mathbf{x} \mathbf{x}, \mathbf{y}\_{TK'}, \mathbf{x} \mathbf{x} \mathbf{x}, \mathbf{y} \mathbf{x} \mathbf{x})\|\tag{2}$$

The RAPID has been widely used in the field of nondestructive testing. However, some disadvantages of this algorithm are gradually being exposed with the increasing demands of imaging quality. Those disadvantages are listed as follows. In signal extraction, the time window Δ*T* in Equation (1) should cover at least two modes of Lamb waves because it is necessary to know the difference between the current signal and baseline signal. A smaller correlation coefficient value means a stronger indication of the defect's appearance. This will undoubtedly make the guided wave tomographic imaging become more complex. Particularly, this adverse effect will become more serious when the air-coupled transducer is used because air-coupled transducers are insensitive to the symmetric modes of a Lamb wave. On the other hand, only one mode of the Lamb wave is needed if the DHB algorithm is used. In defect quantitative analysis, the parameter β determines the shape of the elliptical path. Then, the speckle of the defect shape can also be influenced. There is one optimal value for specific defects with a determined size. Other defects of different sizes cannot be imaged successfully. The influence of β on the imaging results is shown in Figure 1.

**Figure 1.** Image reconstruction results by using different values of parameter β.

The FBP algorithm is also widely used in the NDT field. Its reconstruction formula can be written as:

$$f(\overrightarrow{\textbf{x}}) = \int\_{0}^{2\pi} \frac{D}{\left(D + \overrightarrow{\textbf{x} \cdot \textbf{e}}\_{w}\right)^{2}} \int\_{-\gamma\_{m}}^{\gamma\_{m}} (q(\lambda, \gamma) \cos \gamma) \otimes h\_{H}(\gamma) d\gamma d\lambda,\tag{3}$$

where <sup>→</sup> *x* is the vector of the reconstruction points, λ is the angle between the starting coordinate axis and the line that connects the trigger points and the origin of the starting coordinate system, λ*<sup>m</sup>* is the angle range of the fan beam, *q*(λ, γ) is the projection value, <sup>→</sup> *e <sup>w</sup>* is unit orthogonal vector, *hH* is kernel of the Hilbert transform, and D is the distance between the trigger points and the intermediate point of the receiving sensors' array.

Many different kinds of iterative reconstruction methods can be used for the data acquisition, as shown in Figure 2. Only the fan-beam method was used in this work to emphasize the comparison of the imaging quality of different tomographic imaging algorithms.

**Figure 2.** Three common iterative reconstruction methods: (**a**) Crosshole; (**b**) Double-crosshole; (**c**) Fan-beam.

As shown in Figure 3, it is based on the signal of a nondefective ultrasonic signals, and the nondefective signal becomes a reference signal. When the ultrasonic wave passes through the defect, the ultrasonic amplitude is attenuated, so the defect or object is imaged by utilizing the difference in amplitude between the two signals. The simulated results were made based on the attenuation law of a Lamb wave. As shown in Figure 4, results based on different detection points were made for comparison. The FBP algorithm is a technique used in radiographic inspection. It is a method of transmitting X-rays through an object and receiving it through a detector. It shoots 360 degrees around the object and displays the image. However, radiographic inspection is limited by shielding facilities or radiation regulations. In this study, ultrasonic technology using a Lamb wave was applied to compensate for such shortcomings. Figure 4 is the result of imaging by applying the simulation signal of Figure 3 to the FBP algorithm. The imaging quality can be improved by increasing the emitter and receive detection points. This method has the same characteristic as the new DHB algorithm in this respect. The number of the detection points adopted in the actual experiment was 90, 120, and 180.

**Figure 3.** Guided wave signal: (**a**) No defect signal; (**b**) Defect signal.

**Figure 4.** Simulated results from the filtered back projection (FBP) algorithm: (**a**) 180 Points; (**b**) 120 Points; (**c**) 90 Points.

The FBP algorithm has a big advantage in the aspect of imaging precision compared with RAPID. However, the problem of this algorithm is in its image artifacts. Defect imaging close to the scanning edge cannot be displayed clearly, as shown in Figure 4. These are called image artifacts. The image artifacts will undoubtedly affect the precise identification of the defects near the detection boundary.

The new DHB algorithm is proposed to solve the shortcomings of the algorithms discussed above. The new DHB algorithm has good performance in the detection of plate-like structures. The imaging quality can be significantly improved compared with other algorithms, such as RAPID or FBP. Particularly, it can give accurate size and shape of different kinds of defects.

The geometry of the corresponding fan-beam iterative reconstruction method is shown in Figure 5. Its reconstruction formula is Equation (4) [24]. The relevant parameters in Equation (4) are the same as in Equation (3). The unit orthogonal vector <sup>→</sup> *e <sup>w</sup>* is represented by <sup>→</sup> *a* (λ)/ → *a* (λ) and is substituted into Equation (3). The notation *gm*(λ, γ) is used to denote the measured projections. Three steps are included in using this method: the difference of the projection data, the Hilbert transform of the difference data, the distance weight, and back projection. The corresponding geometry schematic can be referenced in Figure 5.

$$f(\overrightarrow{\mathbf{x}}) = \frac{1}{4\pi} \int\_0^{2\pi} d\lambda \frac{D}{\|\overrightarrow{\mathbf{x}} + \overrightarrow{a}(\lambda)\|} \int\_{-\gamma\_m}^{\gamma\_m} d\gamma h \mu \sin(\gamma' - \gamma) \left(\frac{\partial}{\partial\lambda} + \frac{\partial}{\partial\gamma}\right) \mathbf{g}\_m(\lambda, \gamma). \tag{4}$$

**Figure 5.** Fan-beam geometry based on the equiangular rays.

The same simulated results have been made by using the DHB algorithm. As shown in Figure 6, the problem of image artifacts can be solved perfectly. The imaging quality is also improved to a large extent.

**Figure 6.** Simulated result from the Difference Hilbert Back Projection (DHB) algorithm: (**a**) 180 Point; (**b**) 120 Point; (**c**) 90 Point.

In addition, short-scan conditions can also be relaxed once the new DHB algorithm is used. It is generally believed that accurate and stable reconstruction of any region of interest (ROI) requires line integrals of the object density for all lines passing through the object. This condition leads to the notion of circular short scan. As shown in Figure 7, π + θ is referred to as a short scan.

**Figure 7.** Schematic of the short-scan problem: (**a**) A conventional short-scan of π + θ allows reconstruction of the whole object; (**b**) A continuous scan of less than π + θ allows reconstruction of all object points inside the convex hull of the scan; (**c**) A scan of three equally spaced segments of 80◦ each allows reconstruction of a triangular ROI in the center of the object.

All line integral measurements through the object are required for any ROI reconstruction using the conventional FBP algorithm. However, this problem can be relaxed using the new DHB algorithm. As shown in Figure 7c, small central ROIs can be reconstructed using several short arcs if only an ROI inside the object needs to be reconstructed. Reconstruction is possible for any ROI inside the central triangular region obtained by joining the beginning point of each arc with the end point of the next arc.

#### **3. Experimental Set-Up and Specimen**

The hybrid application of an air-coupled transducer and pulse laser is used in this work. Air-coupled sensors can have increased sensitivities compared to optical transducers. Well-designed capacitive air transducers can be expected to be as sensitive as laser interferometric sensors [25,26]. In addition, the air-coupled transducers can receive a specific mode of guided waves. This can be realized by changing the receiving angle to a leak direction of the selected modes. It is well-known that a static scanning process of the receiver is very important for high-precision acquisition of the signals. The application of air-coupled transducers is one of the most appropriate choices to realize high-precision acquisition. However, the mismatch of acoustic impedances between solid materials and air is a huge weakness of a purely air-coupled transducer detection system. Consider the scenario shown in Figure 8. The ultrasound passes through the interface between the solid specimen and air two times. Only a tiny fraction (10−8) of the original acoustic energy is transmitted through both solid–air interfaces to reach the receiver.

**Figure 8.** Transmission in air using two air-coupled transducers. Values shown are approximate fractions of the original energy.

Consider the second scenario shown in Figure 9. The ultrasound is generated directly in the solid using a pulsed laser. The solid–air interface on the generation part has effectively been removed. The energy reaching the receiver could increase by a factor of 104 [27]. As a receiver, the air-coupled transducer can also have good performance if the proper parameters are set [28]. The hybrid application of an air-coupled transducer and pulse laser can make material characterization more feasible. It can also be used where the air-coupled transducer and pulse laser are put on the opposite sides of the specimen [29,30]. The laser ultrasonic technique can also be used in the detection of a curved surface specimen [31,32].

A line laser source is much better than a point source for long-range detection [33,34]. Many researchers used spatial array illumination sources produced by several means. Those methods include slit masks, lenticular arrays, optical diffraction gratings, multiple lasers, and interference patterns. The linear slit array mask was used in this study, as shown in Figure 10. A Lamb wave of specific modes can be generated selectively by adjusting the line spacing. The distance between two adjacent line laser sources is defined as the line spacing or element gap marked as Δ*S*. The value of the line spacing is equal to the wavelength of the generated Lamb wave. The difference between it and the element width equals the width of the single line laser source. The width of the single line light source is also known as the line width marked as w. The value of *w*/Δ*s* is known as the duty ratio. It is found that an A0 mode Lamb wave has the optimum signal–noise ratio when the value of *w*/Δ*s* is 0.5 [35].

**Figure 9.** Improved transmission using a laser source. Values shown are approximate fractions of the original energy.

**Figure 10.** Schematic of the line laser sources of 1-D and Gaussian model: (**a**) One-dimensional theory; (**b**) Gaussian model.

As shown in Figure 10b, the generated wave should be analyzed by using a Gaussian model. Specifically, the spatial energy profile of each line beam is assumed as Gaussian. A Gaussian model is also suitable for the focused laser beams. However, the diffraction effect is negligible when the slit mask is very close to the target surface. Then, the spatial energy profile of each line beam is closer to a square pulse than Gaussian in the line-arrayed slit mask shown in Figure 10a. The energy of each line laser beam can be assumed to have a uniform distribution, and the temporal profile is a delta function. The spatial energy profile of each line beam can be taken as a simple one-dimensional square pulse. For an array of N equally spaced identical laser sources, the final displacement can be taken as a summation of N single-source surface waves displaced in time by the surface wave propagation interval between two adjacent sources, as shown in Equation (5).

$$g(t) \;=\sum\_{n=1}^{N} f(t - n\Delta t) \tag{5}$$

In the above equation, Δ*t* = Δ*s*/*c* and *f*(*t*) is the out-of-plane displacement of the guided wave generated by a single line laser source. It is expressed as a function of time in the thermoelastic regime. The displacement is proportional to the derivative of the absorbed laser energy function. The energy distribution of the line laser beams was regarded as a square pulse, and *f*(*t*) can be approximated by

the square pulse derivative. The relationship between the Lamb wave frequency dispersion curve and the corresponding parameters is given in Figure 11.

**Figure 11.** The frequency dispersion curve of a Lamb wave and the corresponding parameters: (**a**) Phased velocity; (**b**) Group velocity.

The signal is generated on the specimen surface when the illuminated laser beam transmits through the slit mask. The wavelength corresponds to the element gap of the slit array. The relation between the wavelength, phase velocity, and wave frequency is shown in Equation (6), where *d* is the specimen thickness, *f* is frequency, and λ is wavelength. The principle of generating the desired mode is shown in Figure 11.

$$\mathbb{C}\_{ph} = f \cdot \lambda \, = \, f \cdot \Delta \mathbf{s} = f \cdot d \cdot \left(\frac{\Delta s}{d}\right). \tag{6}$$

As shown in Figure 11, a diagonal line with a slope of Δs/d is described in the dispersion curves of the Lamb wave. The modes of the curves that have cross points with the diagonal line can be generated. Therefore, various modes with different frequencies can be generated simultaneously.

The experimental setup is shown in Figure 12. A Nd:YAG pulsed laser with a wavelength of 1064 nm and a pulse duration of 6 ns was used as the energy source. A Quantel Brilliant B pulsed YAG laser was used with principal radiation at 1064 nm, max energy/pulse 850 mJ, max average power 7.5 W, and pulse duration 6 ns. This laser beam was expanded and illuminated on a slit mask. The slit mask is located close to the specimen surface to suppress the diffraction phenomena.

**Figure 12.** Schematic diagram of the experimental setup.

An unfocused planar air-coupled transducer with a circular aperture was used as the receiver. It was placed on an automatic scanning robot arm. The diameter of its aperture area is 10 mm. The air-coupled transducer used in the experiment is a transducer capable of transmitting and receiving broadcast frequencies. It can detect an ultrasonic wave in air with a frequency range of 0.2 to 0.5 MHz. The lift-off distance of this air-coupled transducer from the specimen surface is 7 mm. The air-coupled transducer has a small amplitude because it receives ultrasonic waves through the air noncontacting the specimen. Therefore, an amplifier is used to amplify the signal. There are two key points that contribute to the imaging quality of the tomography. One is the algorithm, and another one is the high-accuracy acquisition of signals. By using the automatic scanning robot arm, the air-coupled transducer can be positioned automatically with the right receiving angle. The details of the automatic position will be illustrated. The wave propagation direction and the receiving direction of the air-coupled transducer should be kept in the same straight line in practical application. Although the line source is based on a one-dimensional model, the pulse laser has a Gaussian model [36]. This makes the accurate positioning of the air-coupled transducer very important. Thus, an automatic positioning robot arm was developed in this work.

Figure 13 shows the characteristics of the Gaussian beam. In the wave traveling toward the focus (*z* = 0), the beam width is continuously reduced, so that the beam width is minimized at the focus. Wave past focus increases beam width almost linearly. Suppose that the electromagnetic wave is traveling in the +z direction, and also assume that the Gaussian beam changes like a Gaussian function in the *x* and *y* axis directions perpendicular to the *z* axis. Then, Equation (7) can be obtained.

*<sup>E</sup>*(r) = *<sup>A</sup>*(*z*) exp! *i k*ρ<sup>2</sup> 2*q*(*z*) " ·*e ikz*, *where* ρ<sup>2</sup> = *x*<sup>2</sup> + *y*<sup>2</sup> (7)

**Figure 13.** Gaussian beam propagation factor.

Here, *K* is the wave number, and the gauss beam moves in the +z direction, so the beam width widens and the phase changes, so if *q*(*z*) in Equation (8) is defined as plural, *ZR* can be defined by considering Figure 13. Referring to Figure 13, *Wo* has the smallest beam width (*z* = 0), so measurement is possible. That is, since Wo is half the beam width, *ZR* is also determined by *W*(*z*). *ZR* is called Rayleigh length, and *<sup>W</sup>*(*z*) is increased by <sup>√</sup> 2 times at *Z* = *ZR*, and the position of *Z* = *ZR*, which is doubled by the area of the beam, becomes Rayleigh length.

$$\left|A(z)\right| = \sqrt{2\eta P\_0} \sqrt{\frac{2}{\pi}} \frac{1}{w(z)}\tag{8}$$

Even if the Gaussian beam travels in the +z direction, the total power must not change, so *A*(*z*) can be expressed as Equation (8). Substituting 2 into Equation (7), the Gaussian beam formula is

defined as Equation (9). By installing the lens using the formula, a plane wave is generated through the lens, and when passing through the slit, a desired wave can be created.

$$E(\mathbf{\tilde{r}}) = E\_0 \sqrt{\frac{2}{\pi}} \frac{w\_0}{w(z)} \exp\left[-\frac{\rho^2}{w^2(z)} + i\frac{k\rho^2}{2R(z)}\right] \mathbf{e}^{ikz} \tag{9}$$

As shown in Figure 14, the scanning process consisted of four steps, and the initial condition is given in step 1. The air-coupled transducer should firstly revolve at a certain angle α around the test center when it is necessary to change the detection points. Then, the slit mask is rotated to make the wave propagation direction stay in the same straight line with the new receiving angle. This process can be very complicated with manual computation and operation. Thus, automatic control was adopted in this work. The LED laser line sources and photoresistors were installed on the edge of the slit mask and the rotation center of the air-coupled transducer. They were marked as LED1/<sup>2</sup> and R1/2, respectively. All the automatic movement in this work is based on a step motor. To rotate the slit mask, LED1 will be turned on at the same time as the step motor M1. The illumination direction of LED1 is the same as that of the wave propagation. The step motor M1 will stop by itself once the light of LED1 points to the photoresistor R1. Then, LED2 will turn on together with M2, and the step motor M2 will also stop automatically once the light of LED2 points to the photoresistor R2. The new detection point can be fixed with the right receiving angle automatically based on the process discussed above. The complex calculations and manual operation can be omitted.

**Figure 14.** Schematic of the automatic scanning robot arm.

The design of the automatic scanning is shown in Figure 15. The key part here is the photoresistor. Its electrical resistance will decrease immediately if it is illuminated by the light of the laser LED. Then, the corresponding circuit that was used to control the step motor will close. Conversely, the circuit will turn off. In this work, a single chip computer was used to control the step motor. The step motor will be fixed by its magnetic pole once the circuit is controlled by the photoresistor. It can guarantee that the experimental equipment will never be influenced by potential factors from outside. In other words, the movement of the transducer can only be controlled by a program. Then, the detection points can be positioned with very high accuracy with the right receiving angle of the air-coupled transducer.

**Figure 15.** The working principle of the automatic scanning robot arm.

In a practical environment, the photoresistor can also be influenced by sunlight or the bulbs in the lab. Those negative factors can influence the normal operation of the photoresistor, as shown in Figure 16. Adjustable resistors R3 and R4 were installed in series with the photoresistors to avoid these influences. They can improve the adaptability of the equipment in different working environments. The photoresistor will only respond to the laser LED by using the method mentioned above. Both the air-coupled transducer and the pulse laser are on the same side of the specimen in the experimental setup. The automatic clamping device will be put on the other side of the specimen. The accurate collection of the signals cannot only rely on the automatic positioning of the air-coupled transducer. The precise movement and fixation of the specimen are equally as important as they are for the air-coupled transducer. Accurate collection can be achieved by using an automatic clamping robot. As shown in Figure 17, it can clamp different kinds of specimens automatically by using four robot arms. The specimens include plates, curved plates, or even pipes. There are four lead-screw mechanisms that are used to control the robot arms. Those robot arms can stop by themselves once they come into contact with the specimen boundaries. The movement of those robot arms will never be influenced due to the self-locking phenomenon of the lead-screw mechanisms.

**Figure 16.** The actual settings of the automatic scanning robot arm.

**Figure 17.** The 3-D model of the automatic clamping robot arms.

The key part of the automatic clamping robot arm is an automatic clutch. The robot arm can clamp the plate specimen with different shapes automatically by using this part, and the clamping force is adjustable. As shown in Figure 18, two transmission shafts are connected together based on elastic force. Those two shafts will always be connected if *Fp* is not bigger than a maximum value. Otherwise, no force can be transferred through this mechanism. The corresponding calculation formulas are shown in Equations (10) and (11).

$$F\_p := F\_\ell \prec F'\_\ell. \tag{10}$$

$$F\_p' = F\_x' = -k \cdot \mathbf{x} \tag{11}$$

**Figure 18.** The force analysis of the automatic clutch.

Each of the four robot arms can stop by itself if it comes into contact with the specimen's boundaries. At the same time, the other robot arms will continue to clamp the specimen until they also come into contact with the specimen's boundaries.

In this work, all of the linear motion is based on the lead screw mechanisms, and all of the circular motion is based on worm and gear mechanisms. A self-locking effect can be realized by using both mechanisms. The reliability of the necessary movement in actual operation can be guaranteed in a further step based on the self-locking effect. The actual experimental platform is shown in Figure 19. The air-coupled transducer was placed on the same side of the specimen together with the pulse laser. The clamping robot arms was fixed on the other side of the specimen. The clamping robot arms were installed on a scanning device through which high-accuracy positioning of the specimen can be achieved.

**Figure 19.** The actual experimental platform.

#### **4. Results and Discussion**

In this work, the test specimen is an aluminum plate with a thickness of 1 mm. The wavelength was designed as 6 mm. Both the A0 and S0 modes of a Lamb wave can be excited based on the laser ultrasonic technique. Because the air-coupled transducer is more sensitive to the antisymmetric modes than the symmetric modes, the detection range of the air-coupled transducer was chosen as 0.2 to 0.5 MHz to detect the A0 mode Lamb wave.

An experiment was done to see the appearance of a defect, as shown in Figure 20. There is a penetrating hole in the center of the specimen. The diameter of this through hole is 20 mm.

**Figure 20.** The schematic of the defect-detecting experiment.

The signal of both the defect-free and with the defect areas have been obtained in both time and frequency domains. As shown in Figure 21, the signal attenuates clearly due to the appearance of the defect. Only an A0 mode Lamb wave can be detected, as mentioned above. The generated S0 mode was filtered by the air-coupled transducer due to its detection range. Only one mode of the Lamb wave was needed for the tomography once the new DHB method was used. This is also another advantage of the new algorithm. For example, the probabilistic algorithm requires at least two modes of Lamb waves for comparison. The tomographic imaging technique used in this work is based on the wave attenuation. The attenuation law of a Lamb wave is subject to Equation (12), and the equation for data projection is given in Equation (13). The parameter α in those two equations is the attenuation coefficient, *A* is the signal amplitude, and *x* is the propagation distance.

$$A = A\_0 e^{-ax} \tag{12}$$

$$\int\_{-\infty}^{\infty} \alpha(\mathbf{x})d\mathbf{x} = \ln \frac{A\_0}{A} \tag{13}$$

**Figure 21.** The signals of the defect detecting experiment.

The attenuation curve and the corresponding scanning sketch based on the fan-beam iterative method are shown in Figure 22a,b. The signal attenuation becomes evident within the scanning range of 83◦ to 90◦ in the area where a defect starts to appear. Particularly, the illumination area can cause a "dead zone" due to the geometrical features (the diameter of the slit mask used in this work is 60 mm). In other words, the detection cannot be achieved effectively along a full circle arc when the line laser source is used. As shown in Figure 22, the signal within the scanning range of 0◦ to 13◦ was invalid data.

**Figure 22.** The attenuation curve of Lamb wave obtained from the experiment: (**a**) Attenuation curve; (**b**) Schematic of scanning direction.

It is unnecessary to test the "dead zone" when the DHB algorithm is used because the short-scan condition can be relaxed to a great extent. The real data shown in Figure 22 also supports the advantage of this new DHB algorithm compared with FBP or other algorithms. Figure 23a,b is a sample drawing and specimen of the initial single defect. The imaging of the penetrating hole is given in Figure 24. There are 36 emitter and receive points. If there are 19 receiving points from a single emitter point to a semicircle, and inspection is performed up to the 20th emitter point, 361 data points can be collected. Complex shapes or multiple defects should be inspected up to the 36th emitter point, but a single circular defect, as shown in Figure 24, was able to quantitatively represent the defect. The imaging of the defect is displayed clearly compared with the defect-free area.

**Figure 23.** Circular defect specimen #1: (**a**) Circular defect specimen; (**b**) Circular defect specimen drawing (Unit: mm).

**Figure 24.** Circular defect imaging.

Another experiment was conducted to further verify the effectiveness of the new DHB algorithm in the aspect of detection precision. Different types of defects were made in different positions of specimen #2. As shown in Figure 25, the defects include 10 micropores, one narrow through defect (width 2 mm, length 52 mm), and two nonpenetrating defects (circle defect: diameter 52 mm, depth 0.5 mm) (rectangular defect: width 52 mm, length 26 mm, depth 0.2 mm) with different depth, which can be considered as an imitation of corrosion defects. There are 106 emitter and receiver points. If there are 53 receiving points from a single emitter point to a semicircle, and inspection is performed up to the 106th emitter point, 5618 data points can be collected. Specimen #2 has multiple defects in various shapes, so the more imitators and reception points, the more quantitatively the defects can be represented. Although it has so many emitters and receiving points, it was able to receive uniform data through the robot arm. Figure 26 shows the imaging results of specimen #2. The imaging results agree very well with the actual defects.

**Figure 25.** Image reconstruction of specimen #2: (**a**) Multidefect specimen; (**b**) Multidefect specimen drawing (Unit: mm).

**Figure 26.** Multidefect imaging using the DHB algorithm.

#### **5. Conclusions**

The hybrid application of an air-coupled transducer and pulse laser has been adopted in this paper. The FBP algorithm, which is mainly used in radiographic inspection, was applied to ultrasonic testing and tomographic imaging, and the disadvantages of image artifacts were confirmed. The factors affecting image quality, such as image artifacts, were overcome with the DHB algorithm through three projection data processing steps; Hilbert transform of the difference data, distance weight, and back projection. A new DHB algorithm has been successfully applied to give accurate imaging of different kinds of defects in plate-like structures. The imaging results agree very well with the actual defects. Particularly, the detection efficiency has been improved to a large extent by using the automatic scanning system. Further research will focus on the detection of different kinds of specimens with different geometries, such as curved surfaces. The algorithm should be further improved to make high-quality images with fewer detection points. The presented techniques can be applied to structural health monitoring to predict faults by quantitatively imaging the growth status of the faults.

**Author Contributions:** J.P. and Z.L. designed and performed the experiments. Y.C. conceived the original idea. J.P. developed the theory and performed the computations with support from J.L. and Y.C. verified the analytical methods. J.L. and Z.L. wrote the draft paper, J.P. completed the final paper, and Y.C. made the final review. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by National Research Foundation of Korea (NRF) grant funded by Korea government (MSIT)(NRF-2020M2D2A1A02069933).

**Acknowledgments:** Authors appreciate the support from Wang Yu at University of Electronic Science and Technology of China.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Acoustic Scattering Models from Rough Surfaces: A Brief Review and Recent Advances**

**Michel Darmon 1,\*, Vincent Dorval <sup>1</sup> and François Baqué <sup>2</sup>**


Received: 22 September 2020; Accepted: 18 November 2020; Published: 23 November 2020

**Abstract:** This paper proposes a brief review of acoustic wave scattering models from rough surfaces. This review is intended to provide an up-to-date survey of the analytical approximate or semi-analytical methods that are encountered in acoustic scattering from random rough surfaces. Thus, this review focuses only on the scattering of acoustic waves and does not deal with the transmission through a rough interface of waves within a solid material. The main used approximations are classified here into two types: the two historical approximations (Kirchhoff approximation and the perturbation theory) and some sound propagation models more suitable for grazing observation angles on rough surfaces, such as the small slope approximation, the integral equation method and the parabolic equation. The use of the existing approximations in the scientific literature and their validity are highlighted. Rough surfaces with Gaussian height distribution are usually considered in the models hypotheses. Rather few comparisons between models and measurements have been found in the literature. Some new criteria have been recently determined for the validity of the Kirchhoff approximation, which is one of the most used models, owing to its implementation simplicity.

**Keywords:** acoustic wave; scattering; random rough surface; impedance

#### **1. Introduction**

The scattering of waves from rough surfaces is an important phenomenon in diverse areas of science, such as optics, geophysics, communications, acoustics, hydraulics [1], elastodynamics for non-destructive testing, etc. For instance, in the particular case of ultrasonic waves, scattering from rough surfaces leads to many applications in non-destructive testing [2], in underwater acoustics for high resolution ultrasonic sonars [3] and for medical sciences [4,5]. The modeling of such phenomena is a century-old and still active research topic.

Scattering models of acoustic waves from rough surfaces could be classified into analytical models, numerical ones and combinations of both analytical and numerical methods. However, such a classification is rather ambiguous, since some analytical models (usually called semi-analytical) employ numerical calculations and since some numerical methods start from analytical approximations. The current review focuses on analytical or semi-analytical methods and the validity of each method will be emphasized as well as some applications examples. Some analytical methods are presented or used in the following in more detail: the Kirchhoff approximation (KA), the perturbation theory (PT), the integral equation (IE) method, the parabolic equation (PE) method and the small slope approximation (SSA). The Kirchhoff and PT methods represent early classical approaches but are still very often employed in the literature, notably the Kirchhoff one. The other methods represent more modern approaches which have larger domains of validity but are more complex and therefore more difficult to implement or can lead to a higher calculation time. All the previously cited methods have

been found to be amongst the most common in the literature and many of the other methods of the literature are based on or have much in common with these approaches.

This review focuses thus only on the scattering of acoustic waves and does not deal with the transmission through a rough surface of waves within a solid material. Indeed the models used in acoustics are often also employed in elastodynamics. For instance, the elastic Kirchhoff model has been devised to deal with rough defects [6–8] and to take into account mode-conversions between compressional and shear waves. A variety of analytical models involving reflection and transmission waves (for instance through a rough fluid–solid interface), such as the phase screen approximation [9–12] and Kirchhoff theory, have been developed to understand the elastic fields in solids or scattered from different embedded scatterers.

Several reviews of scattering from rough surfaces exist in the literature, notably in general physics (diverse areas of sciences) for random rough surfaces [13–15] or in electromagnetism [16,17]. The main review describing works in acoustic scattering from rough surfaces has been written by Ogilvy in the 80s [18], although its review does not only deal with acoustics. This review has notably summarized the two classical approximations (Kirchhoff and the perturbation theory) but it is now out of date: since this review, important models have been devised for grazing incidence or observation, such as the small slope approximation, and new criteria have notably been established for defining the validity of the Kirchhoff approximation, which is probably the most used analytical model for roughness scattering. Some other textbooks also deal with acoustic scattering from rough surfaces notably in the view of ocean applications [19,20]. At Ogilvy's time, the accuracy of Kirchhoff approximation came from Thorsos' work [21] and was determined by the ratio of the height correlation length *L* of the random rough surface to the wavelength. It is also the case in Voronovich's more recent textbook (see Figure 1) [20], where approximations KA, PT and SSA are all said to require this ratio less than a constant in order to be valid; another different criterion has even been proposed for radar applications [22]. Recent studies—described in the present review and based on numerical calculations—have shown for elastodynamic problems that the accuracy of the Kirchhoff theory should depend on the parameter σ2/*L* (σ and *L* being the RMS surface height and the height correlation length). Kirchhoff is considered accurate when this parameter is smaller than an upper bound c depending on the wavelength: (σ2/*L* < c). Kirchhoff accuracy is also dependent on incidence and scattering angles. In the recent years new developments or comparisons of the previous methods have also been proposed in the literature and for different applications, such as non-destructive evaluation (NDE), underwater acoustics or even medicine.

Section 1 of the current paper is devoted to the parametrization of rough surfaces, Gaussian height distribution being usually assumed in modeling. The two historical approximations (Kirchhoff and the perturbation theory) are then described (Section 2) as well as their use in the scientific literature, their validity and an application example for non-destructive testing. Then, more modern approaches with larger domains of validity (notably suitable for grazing observation angles and generally having applications for underwater acoustics) are presented in Section 3.

#### **2. Rough Surfaces Parametrization**

This section is devoted to mathematical methods used to characterize randomly rough surfaces. A rough surface is described using its height deviation from a mean smooth surface. Hereafter are presented the statistical distributions (height and correlation functions) and associated parameters usually employed to properly define the rough surface properties. We will here take the example of Gaussian height distributions since most of the used models of rough surface scattering assume such a distribution. However, although many surfaces (both natural and engineered) exhibit Gaussian height distributions, there are also many others that do not [23,24]. Examples of generating non-Gaussian rough surfaces can be found in the literature [25,26].

#### *2.1. Statistical Properties*

Surface roughness is described here as the variation of height *z* = *h*(*x*, *y*) above or below a mean surface. <*h*>, the average of the height variation over the surface, is by definition 0. *h* is assumed to be random and described statistically using parameters such as the RMS height and spatial correlation lengths. In the theory presented hereafter and developed by Ogilvy [6,27], it is assumed that the height distribution is well represented by a Gaussian.

The probability of any measured height being between h and h + dh above or below the mean plane is then given by p(*h*) dh, where:

$$p(h) = \frac{1}{\sigma (2\pi)^{\frac{1}{2}}} e^{-h^2/2\nu^2} \tag{1}$$

where σ is the RMS height:

$$
\sigma = \left[ \left< h^2 \right> \right]^{\frac{1}{2}} \tag{2}
$$

Random rough surfaces show some spatial correlation, as the height at a location is not totally independent of surrounding heights. This may be described by the correlation function of a surface, *C*(**R**) where:

$$\mathcal{C}(\mathbf{R}) \, \, = \frac{1}{\sigma^2} \langle h(\mathbf{r})h(\mathbf{r} + \mathbf{R}) \rangle \, \, \, \tag{3}$$

which represents the extent to which the surface height at any one point **r** is, on average, related to the surface height at point **r** + **R** away. Gaussian functions are often used as an approximation of the correlation of real surfaces. Assuming the same Gaussian correlation in every direction, *C* can be written:

$$C(\mathbf{R}) \, := \varepsilon^{\aleph \downarrow \natural^z} \tag{4}$$

where *L* is called the correlation length of the surface. Anisotropic correlation functions can also be used to describe surface with a direction-dependent correlation. Such surfaces often occur in practice in industrial structures, as roughness can be generated in anisotropic manner during the manufacturing process.

#### *2.2. Generation Process of an Individual Gaussian Rough Surface*

The method of generating these surfaces is described in detail in [6,27]. It first generates uncorrelated random numbers with a Gaussian distribution of zero mean and specified RMS value. Then, it performs a moving average with weights that ensure that the correlation function is Gaussian with the desired correlation length. The following weights yield a surface with two independent lengths *Lx* and *Ly* in the two perpendicular directions:

$$\mathcal{W}\_{ij} = e^{-2\left[\left(i\Lambda x\right)^{2}/L\_{x}^{2} + \left(i\Lambda y\right)^{2}/L\_{y}^{2}\right]} \tag{5}$$

Δ*x* is the interval between adjacent surface points in the x direction and Δ*y* is the corresponding interval in the y direction. The surface correlation function is defined for a reference point (*x*0, *y*0) by:

$$C(\mathbf{x}, y) \;= \frac{\langle h(\mathbf{x}\_0, y\_0)h(\mathbf{x}\_0 + \mathbf{x}, y\_0 + y) \rangle}{\langle h^2 \rangle} \tag{6}$$

The weights given in Equation (5) yield the following surface correlation function:

$$\mathcal{L}(\mathbf{x}, \mathbf{y}) \;= \mathcal{e}^{-[\mathbf{x}^2/L\_x^2 + \mathbf{y}^2/L\_y^2]} \tag{7}$$

This method allows generating different realizations of random surfaces sharing the same statistical parameters (see Figure 7d for an application example).

#### *2.3. Implementation Example for Generation of an Individual Gaussian Rough Surface*

The parameters to be chosen as inputs for the parametrization of a Gaussian rough surface are:


The number of mesh points of the rough surface along the local axes x and y has then to be automatically fixed by the chosen mesh algorithm. The previously described Gaussian rough surface model is currently implemented for planar and cylindrical surfaces in 3D computations, and for their linear cross-sections in 2D, in the CIVA software [28] (platform for non-destructive evaluation simulation [29–33]). Examples of generated surfaces for small specimens with relatively high roughness are shown in Figure 1.

**Figure 1.** Exaggerated examples of rough surfaces generated for CIVA calculations with standard deviation of the height set to 1 mm for different kinds of rough surfaces: (**a**–**c**) linear profile for respective ratio *L*/σ = 1, *L*/σ = 1.5 (case of Figure 6) and *L*/σ = 4.3 (case of Figure 11); (**d**) planar for *L*/σ = 1; (**e**) cylindrical for *L*/σ = 1. Scales in mm.

#### **3. Historical Approximations**

This section deals with the perturbation and Kirchhoff approximations. Both have been extensively applied to rough surfaces and have their own strengths and weaknesses. Other approaches will be dealt with in later sections.

Figure 2 shows the geometry used in this section for rough surface scattering.

**Figure 2.** Variations of heights around an average plane.

We consider a rough surface whose roughness is defined by z = *h*(x,y) where *h* is the height variation above or below the mean plane, z = 0. The xz plane contains the incident wavevector.

#### *3.1. The Perturbation Theory*

The perturbation theory (PT) assumes some restrictive hypotheses on the local height *h* and gradient ∇*h* at points of a rough surface:

$$|k|!\eta \ll 1\tag{8}$$

$$|\nabla h| << 1\tag{9}$$

*k* being the incident wavevector modulus.

It is assumed that if previous conditions are satisfied, the roughness is weak compared to acoustic wavelength and slowly varies. It expresses the scattered field at an observation point *r* as an expansion series:

$$
\psi(r) = \psi^{inc}(r) + \sum\_{n=0}^{\infty} \psi\_n^{sc}(r) \tag{10}
$$

where ψ*inc* is the incident wave field and ψ*sc <sup>n</sup>* is the nth-order approximation to the scattered field. The *n* = 0 term is the scattered field from a smooth surface. Then the main principle of the perturbation Theory consists in assuming that the boundary conditions may be expanded in a Taylor series about the mean plane *z* = 0. Mathematically, the boundary conditions on the rough surface expressed as:

$$f(\mathbf{x}, y, z) = 0\tag{11}$$

where *f* can be, for instance, the acoustical pressure for Dirichlet conditions—which can be expanded in a Taylor series about the mean plane *z* = 0 as follows:

$$f(\mathbf{x}, y, z)|\_{z=h(\mathbf{x}, y)} = f(\mathbf{x}, y, z)|\_{z=0} + h \frac{\partial f}{\partial z}|\_{z=0} + \frac{h^2}{2} \frac{\partial^2 f}{\partial z^2}|\_{z=0} + \dots \tag{12}$$

Some applications of the perturbation theory consider first scattering from canonical rough surfaces [18]. For instance for the Lamb's problem (scattering from a surface excited by a load force applied at one particular surface point) [34], the solution is expressed as convolution integrals evaluated only for a triangular irregularity of small height compared to its width and to the wavelength. The application of PT to random rough surfaces leads to consider separately a mean coherent contribution and an incoherent one. For the coherent part, reflection and transmission coefficients are determined. For the incoherent part, the mean intensity of the scattered field is calculated [35]. The determined first-order PT term can be substituted back into the Helmholtz scattering formula to iteratively deduce higher order terms.

Blakemore [36] has then presented a Perturbation Theory at the first order for a fluid/solid interface. PT simulations in backscattering for rigid and elastic solids are very close except at the vicinity of the Rayleigh angle as shown in Figure 3.

**Figure 3.** Blakemore's Figure 8 [36]. Backscattered intensity with the perturbation theory on a steel plate immersed in water. The used parameters are respectively frequency = 2.25 MHz, angular beam width = 3◦, observation distance = 300 mm, RMS surface height σ = 10 μm and the two correlation lengths of the rough surface both equal to 100 μm. Reproduced with permission from [36], Elsevier, 1993.

In this region a small peak in the intensity of the incoherently scattered field is predicted, consistent with experiment. Nevertheless, in the case of the proposed experimental validation in pulse echo configuration [37] (Figure 4), Kirchhoff prediction (not shown here) are about 3 dB lower than PT one at backscattering angles theta = 40◦, and 5.5 dB lower at theta = 70◦.

**Figure 4.** Blakemore's Figures 2 and 11 [36]. Experimental validation of Blakemore's perturbation theory at first order: at left measurements from de Billy et al. [37] on a brass plate immersed in water; at right, Blakemore's results with frequency = 2.2 MHz, angular beam width = 3◦, observation distance = 300 mm, RMS surface height σ = 14.9 μm and the correlation lengths *Lx* = 90 μm and *Ly* = 105 μm. Reproduced with permission from [36], Elsevier, 1993.

Some applications and references of PT for scattering from fluid/fluid, fluid/solid interfaces and fluid/stratified layers [38–40] for inhomogeneous media [41,42] can be found in the literature, generally for ocean applications. An acoustic model in the time domain based on first-order perturbation theory for scattering from layered sediments has recently been proposed [43].

Only a few authors have been interested in the scattering from canonical non-planar rough targets. Notably Bjorno [44,45] recently investigated the scattering problem for an elastic rough sphere. He developed in that case a second-order perturbation model. He compared the scattering patterns obtained versus the adimensional wave-number ka and for different values of the ratio of the rms height to the sphere radius with its perturbation solution to the classical acoustic scattering solution for a smooth sphere [46,47]. Measurements of the backscattering response of a cast iron sphere versus ka led to a good experimental validation of the proposed model: the roughness yields a decrease in the backscattered amplitude notably for large ka values. The accuracy is limited to small values of the ratio rms height to wavelength.

Note also, for instance, a recent development in electromagnetism of the perturbation theory for a slightly deformed cylinder [48].

We also have to cite the phase perturbation method [49,50] in which the phase (rather than the field) is expanded in a series of different orders in kh. This theory exhibits links with both PT and Kirchhoff [51]. A validation [52] of the phase perturbation method has been proposed and has consisted in comparing its results to numerical ones in a region where PT and Kirchhoff are found ineffective. The phase perturbation method is found to work well in such a region except for low grazing observation angles. However, applications of this new theory seem not numerous and are more focused on radar cross sections calculations [53].

Although the perturbation approach can take into account of the scattering and re-radiation of Rayleigh waves, it appears possibly more restrictive than the Kirchhoff method in its range of applicability. The perturbation method is notably an approach valid for low roughness (when compared to acoustic wavelength). Several validations of the PT are provided in the next sections; they also show that PT is rather accurate far from the specular direction.

#### *3.2. The Kirchho*ff *Approximation*

#### 3.2.1. Theory and History

The main hypothesis of the Kirchhoff approximation (also referred to as the physical optics approximation) is to assume that each point of the surface behaves as if it belonged to the infinite tangent plane to this surface. Then it is practical to employ geometrical acoustics to approximate the exact fields on the local surface: the field on the illuminated part of surface is the sum of the incident field and the field reflected by the local tangential plane. This supposed surface field is then propagated far from the surface using a propagator (called Green function). Indeed the scattered far field is expressed by the way of an integral over the flaw surface (the Kirchhoff–Rayleigh integral) whose integrand is—for the example of Dirichlet conditions—the product of the Green function by the derivative of the geometrical field on the flaw surface. For complex incident fields or scatterer shapes, this integral is evaluated numerically.

The Kirchhoff approximation has been extensively used to model scattering of waves in optics (early introduced by Meecham [54]), electromagnetism [55,56], acoustics, etc. In both acoustics [57] and elastodynamics [58,59], the Kirchhoff approximation also has the advantage in dealing with anisotropic materials and impedant surfaces [57].

Note that when doing an asymptotic evaluation of the KA integral on a finite surface, the KA resulting leading terms describe the reflected waves, and the other terms describe diffracted waves by the surface edges which have the same nature (cylindrical or conical wave fronts) but not the same amplitude as those predicted by the geometrical theory of diffraction (GTD) [55,60]. The main advantage of the Kirchhoff approximation is to provide finite results for the scattered field in the whole space, contrary to GTD [61–65]. One of the main inconveniences of the Kirchhoff approximation (KA) is its incorrect quantitative prediction of edges diffraction, since the approximation of the surface field by the geometrical field is invalid near the edges. To handle this drawback, the physical theory of diffraction

(PTD) has been proposed in electromagnetism [66], acoustics [59,67] and elastodynamics [68–70] and compared to other uniform scattering models [67,70,71].

In the literature, the Kirchhoff approximation has been studied by different authors. A review has been carried out by Ogilvy [18]. The first one to introduce the Kirchhoff approximation is Eckart [72].

Ogilvy [27] has well described the theory for 1D rough profile (Figure 2) and Dirichlet conditions (soft interface) extending Eckart's works. De Billy et al. proposed a Kirchhoff-based theory extended to shadowing effects for the backscattering of acoustic waves by randomly rough surfaces [37]. De Billy et al. employed in fact the potential-method formulation developed by Welton which showed by iterating the scattering integral equation that the scattered pressure can be expressed as a series: the first term corresponds to the singly and the others to the multiply scattered field components. De Billy et al. considered only the first term. Their experimental validations of their model using Gaussian rough (liquid–solid interfaces) surfaces have shown that the experimental measurements of the intensity backscattered by randomly rough interfaces (liquid–solid) are in good agreement with the theory until 40◦ backscattering angles (with respect to the normal to the mean plane). Nevertheless, an anomaly has been observed in Kirchhoff's predictions for the Rayleigh's angle of incidence. In fact, to model the impedance surfaces used in their experimental validations, they simply add a reflection coefficient as a factor inside the Kirchhoff integral, and for backscattering the local reflection coefficient is the one at normal incidence for a plane surface: this formalism does not take into account any mode conversion in the elastic medium and is unable to predict any perturbation at the different critical angles. Following de Billy's works which studied the Kirchhoff approximation using monochromatic plane waves, Imbert [73] has proposed a Kirchhoff model in time domain for 2D rough profiles and Neumann conditions (rigid interface).

We can sum up the main steps of Kirchhoff approximation application as follows:

1. Use of the Helmholtz equation: the scattered field is:

$$
\psi(r) = \psi^{\rm inc}(r) + \int\_{S\_0} \left( \psi(r\_0) \frac{\partial G(r, r\_0)}{\partial n\_0} - G(r, r\_0) \frac{\partial \psi(r\_0)}{\partial n\_0} \right) dS\_0 \tag{13}
$$

The equation is also called Green theorem and *<sup>G</sup>*(*r*,*r*0) <sup>=</sup> *eik*|*r*−*r*0<sup>|</sup> <sup>4</sup>π|*r*−*r*0<sup>|</sup> is called the Green function, where *r*<sup>0</sup> is a current point on the scattering rough surface *S*<sup>0</sup> and *r* the observation point. The field in the fluid medium can be obtained thanks to the knowledge of its value on the surface *S*0.

2. Boundary Conditions on the rough surface: Ogilvy [27] considered for instance Dirichlet boundary conditions (soft case) so that ψ(*r*0) = 0. So one term in the previous Green integral disappears:

$$
\psi(r) = \psi^{inc}(r) + \int\_{S\_0} \left( -G(r, r\_0) \frac{\partial \psi(r\_0)}{\partial n\_0} \right) dS\_0 \tag{14}
$$


Note that the Kirchhoff approximation is mainly used to predict far-field scattering; a model has been proposed to evaluate the scattering pattern in near-field [75] after the interaction of an acoustic wave with a rough interface.

#### 3.2.2. Validity

Historically KA has been considered—in the point of view of its local tangential plane assumption—as accurate when the rough surface can be seen planar by the incident wave. In Ogilvy's review [7], KA is thus said valid for krc cos3 θi >> 1, with rc the surface local radius of curvature and θi the incident angle. This historical criterion forgets that multiple scattering and shadowing effects are not taken into account in KA leading to large errors for grazing incidences or higher roughness.

KA prediction has been compared to Integral Equation (IE) results by Thorsos [21] for 1D rough profiles. Thorsos proposed to use the surface correlation length *L* as main parameter for defining the Kirchhoff validity away from the low grazing angle region: KA is said to be accurate for a correlation length greater than the wavelength or even for lower correlation lengths if shadowing effects are not important. We notice this first criterion proposed by Thorsos at this step is too simple and not involve the rms height. In the low grazing angle region some examples of validation with comparison to another methods are given in next Section 4; at low grazing angles, additional validity dependence enters through the RMS slope which is for Gaussian rough surfaces: *s* = tan γ = √2 <sup>σ</sup> *<sup>L</sup>* i.e., the rms height is now involved. However, at shorter correlation lengths, the needed shadowing correction proposed by Thorsos becomes greater to improve Kirchhoff prediction at low grazing angles in the case of backscattering.

Note that in order to improve its prediction for low grazing angles, Chapman [76] provided an improved Kirchhoff formula (using a smoothing special function) notably by including the correlation length of the scattering surface into the reflection coefficient. The modification impact is not noticeable since shadowing and multiple scattering will dominate for low grazing angles. Chapman showed that the PT and KA results differ except in the limit of infinite correlation length (smooth surface).

Shi et al. [77] recently proposed for elastic materials (also for shear waves [78,79]) and for compressional waves a different surface parameter than the surface correlation length *L* used by Thorsos, a function σ2/*L* of both *L* and σ, to define a criterion when the KA use is valid. They found that, with a normal incidence angle θi = 0◦, KA is valid with a scattering angle −70◦ < θs < 70◦ and when <sup>σ</sup>2/*L* <sup>≤</sup> *c*, the constant upper bound *c* being in two dimensions 0.20 <sup>λ</sup>p and in three dimensions 0.14 λp (λp being the compressional wavelength). Figure 5 shows a comparison between KA and a finite element (FE) method for a 2D case when σ2/*L*~*c.* A small incidence angle of 30◦ can improve the KA accuracy when σ2/*L* exceeds *c* as in the case shown in Figure 6. We can however regret that the numerical evaluation of this proposed criterion was limited to small incidences angles: however, this criterion is likely to be suitable for many NDE configurations for which incidence is close to normal. A study of this criterion near grazing incidence should also be useful.

As said in the introduction, another different KA validity criterion has been proposed for radar applications [22] and is also based on numerical calculations: it is given by Equation (1) in [22] (see also their Figure 3, Kirchhoff being called the scalar approximation) and ensures a KA validity for *L*2/σ more than a constant times the wavelength and a condition *kL* > 6 is also associated. We find the condition on *L*2/σ very disbelieving and in contradiction with the previously cited studies and also that given by Voronovich (σ/*L* < *C* and large kL; see Figure 1.1 in [20]). The condition *kL* >6 is also found in other references as [52]. Notably, Figure 2 in [52] represents the stated validity of first-order perturbation theory and KA in *k*σ–*kL* space obtained for bistatic scattering cross sections simulation: PT is said valid for small *k*σ < 2 to 4 depending on *kL* value, KA valid for *kL* > 6 and both KA and first-order PT needs σ/*L* < *C* (*C* being a constant value).

In our opinion, this KA condition *kL*>6 is too restrictive for many applications, such as non-destructive evaluation, which often employs near specular observation inspection. First, we previously study or carry out numerical comparisons of KA simulation for cracks of finite size *L* inspected with longitudinal waves in solids with numerical methods: the method of optimal truncation (MOOT) [80] and a finite element method (FEM) [81,82]. These comparisons with MOOT [83] and FEM [69,84] both show that KA is valid for near specular observation inspection for *kL* > 2. Shi et al. [77] did not propose a validity criterion on *kL* (but only the condition on σ2/*L);* they nevertheless obtained an acceptable disagreement of 2dB

(see their Figure 10 [77]) for *kL* = π/*2* (corresponding to L = λp/4) and for an important roughness σ = λp/4 = *L* equal to the correlation length. It is a pity that that they did not provide numerical validations below *kL* = π/*2.*

**Figure 5.** Shi's Figure 11 [77]. Scattered amplitude versus observation angle using 2D finite element (FE) and Kirchhoff approximation (KA) simulations, for θi = 0◦, σ = λp/3, *L* = 0.6 λp. Reproduced with permission from [77], Royal Society, 2015.

**Figure 6.** Shi's Figure 12 [77]. (**a**) Shi's Figure 11 [77]. Scattered amplitude versus observation angle using 2D FE and KA simulations, for θi = 30◦, σ = λp/3, *L* = λp/2 < 0.6 λp. (**b**) Error between 2D FE and KA with respect to θi in the specular direction for different values of σ. Reproduced with permission from [77], Royal Society, 2015.

Our conclusive opinion is that for many applications using rather near specular observation and non-grazing incident angles, two criteria can be used for KA validity: <sup>σ</sup>2/*L* <sup>≤</sup> *C*λp and kL > 1 or 2. For more important observation angles to the normal the criterion on kL and the constant C value are likely to be a little more restrictive: Zhang [85] using 2D FEM calculations recently showed that KA can be considered accurate for kL > π and for incidence and scattering angles over the range from −80◦ to 80◦.

#### 3.2.3. Elastic Targets and Applications

The Kirchhoff approximation models bulk waves, and does not take into account surface acoustics waves (SAW) as Rayleigh waves neither head waves. Head waves led to recent researches [70,86–88] and a modification [89] of the Kirchhoff approximation is henceforth proposed for head waves account. Head waves on irregular surfaces can also be modelled [90].

The Kirchhoff approximation for elastic rough targets is proposed by Ogilvy [13] and then revisited by Bjorno [91] and Dacol [92]. The Kirchhoff approximation validity have also been specifically considered in the particular cases of rigid and soft targets in a modeling study [93] of the reflected field including multiple scattering. One example of KA application to impedance surfaces [57] for non-destructive evaluation (NDE) is the modeling of ultrasonic telemetry of immersed targets [59,94].

We report here an application of KA simulation to impedance rough surfaces for NDE. Ultrasonic inspection of a planar steel component immersed in water (water path 25 mm) at frequency = 2 MHz (wavelength λ = 0.75 mm) is simulated with the Kirchhoff approximation for elastic rough entry surfaces of varying roughness using CIVA software inspection simulations tools [57,95,96]. A phased-array (represented in yellow in Figure 7a,b) is used to emit an incident beam directed at 117◦ with respect to the surface thanks an emission delay law (applied in red in Figure 7a to only some left elements) and all elements receives the scattered waves without delay law in the angular domain (118◦, 58◦).

In Figure 7a,b, the true B-scans (true-to-geometry imaging) are surperimposed to the specimen in color codes; they represent the echographic signal versus time for all elements in reception (corresponding to scattering angles in the angular domain (118◦, 58◦)). Figure 7c shows that increasing the roughness deteriorates the specular reflection echo and causes appearance of echoes at non-specular angles. Figure 7d shows that the scattering directivity patterns obtained for three different generated rough surfaces of the same statistical properties are slightly different.

A recent study [97] has compared the accuracy of KA and PT using a finite element (FE) model for fluid and elastic ocean bottoms. For a fluid bottom, KA is close to FE with maximal gap of 2 dB except for small grazing incident angles for which PT is more effective. PT behaves globally well except near specular observation direction notably for small grazing incident angles. For the elastic case, the disagreement between KA and FE increases to 5 dB and KA is inaccurate for incident angles less than the longitudinal critical angle: the interferences between transversal, longitudinal and surface waves are not taken into account in KA. The PT prediction is generally good, but another limitation appears in the elastic case for incident angles near the longitudinal critical one. However, we regret the provided validations are established only for a small roughness k σ= 0.6 favorable to PT.

#### 3.2.4. Periodic Surfaces

As it will be discussed in the Section 5, Perspectives, it can also be useful to model scattering from periodic surfaces—such surfaces can notably appear in additive manufacturing parts or in oceans. For the case of periodic surfaces, one exact solution is the classical Rayleigh–Fourier method (RFM) [98]. Impinged by a plane wave, a periodic surface leads to Bragg scattering with a finite number of plane wave propagation angles in far field. A very recent study evaluates the validity of the Kirchhoff approximation (KA) for the scattering problem from a sinusoidal surface—mimicking the ocean one—with comparison to the RFM solution [99]. An asymptotic evaluation of the Kirchhoff integral which is consequently equivalent to geometrical acoustics [100] provides a ray theory which is also evaluated to better understand the Kirchhoff prediction. This notably shows that Kirchhoff has the advantage to exhibit a smooth behaviour near caustics generated by the sinusoidal surface contrary to geometrical acoustics. The KA validity is also studied versus the observation range (projected distance between the source and receiver points): for a short range, KA is very close to RFM. At moderate range, KA sometimes over-estimated the peak pressure due to interferences between several rays. For longer ranges, the KA disagreement increases strongly. This study [99] is very useful for better understanding

the KA limitations due to multiple paths arrivals and grazing incidence/observation. Nevertheless, KA remains a very useful tool even for periodic surfaces when source and observations are not too far.

**Figure 7.** Ultrasonic inspection of a planar steel component (orange box) immersed in water (water path 25 mm) at frequency = 2 MHz (wavelength λ = 0.75 mm): Bscans simulated with KA for a (**a**) smooth and (**b**) a rough entry surface (RMS surface height σ = 1 mm and correlation lengths *Lx* = *Ly* = 10 mm). Scattering patterns (maximum of the echographic signal versus observation angle) (**c**) for different roughnesses (smooth, 10 μm, 0.1, 0.5 and 1 mm) and (**d**) for three different generated rough surfaces of the same statistical properties (roughness σ = 0.2 mm and correlation lengths *Lx* = *Ly* = 10 mm).

#### **4. More Modern Approaches with Larger Domains of Validity**

We have presented the two classical approximations: the perturbation theory (PT) and the Kirchhoff approximation (KA). In the view of main application to scattering from ocean, it is important to cite the two-scale roughness theories [101–104]. They consist in assuming that scattering can be considered as the combination of two phenomena: that produced by large facets and well approximated by KA and that due to small roughness and approximated by PT. One major inconvenient of this technique is that it requires to fix a wavenumber parameter [105] to separate the two scattering domains. The small slope approximation, described in the next section, is a much more convenient technique to link the classical approximations and it reduces to their results in some limiting cases.

In this section, some sound propagation models which have been notably devised for underwater acoustic applications and providing a better validity notably for grazing observation angles are studied. Underwater acoustic propagation depends on many factors. When sound propagates underwater the sound intensity decreases when increasing the propagation distance. Propagation loss is a quantitative measure of the reduction in sound intensity between two points, normally the sound source and a distant receiver. Modeling the acoustic reflection loss [106] at the rough ocean surface needs simulating forward scattering for incident grazing angles and is of interest for sonar especially for shallow oceans or oceans with a surface duct (sound waves trapped in a duct layer below the sea surface) [106]. Several models described hereafter (small slope approximation—SSA, and parabolic equation—PE) were notably designed to handle grazing angles. The integral equation (IE) method is also used as reference solution for validating SSA and comparing it to other approximations as KA and geometrical optics (GO).

#### *4.1. Small Slope Approximation (SSA)*

The SSA technique was described first by Voronovich [107–109]. An alternative small slope approximation has originally been developed by Urusovskii and then studied by McDaniel (1994) [110]. All these methods attempt to exploit the assumed smallness of the surface slope without restricting the RMS surface height to wavelength ratio.

The SSA model from Voronovich has been investigated by Thorsos and Broschat [111,112]. SSA consists in introducing a roughness correction to the Helmholtz integral. This correction is expanded into a functional Taylor series in slope whose terms depend on the integral of the wave number spectrum of the surface height function. The approximation of each order of SSA is based on the approximations of previous orders. Owing to its principles SSA takes into account shadowing and diffraction.

Thorsos and Broschat showed that, at lowest order the small slope approximation of Urusovskii/McDaniel is limited to forward scattering [110] but that the lowest order term in Voronovich's SSA was incorrect in McDaniel's SSA.

SSA tends to the perturbation approximation for small k σ (adimensional wave numbers in terms of RMS surface height σ). The second order SSA (SSA (2)) reduces to the Kirchhoff approximation (KA) at high frequencies, hence does not take into consideration the shadowing effects. The predictions of SSA(2), SSA(3), and SSA(4) are compared [112] with Monte Carlo integral equation (IE) calculations devised earlier by Thorsos [21], and with other approximations: the perturbation approximation, the Kirchhoff approximation (KA), and geometrical optics (GO). The scattering geometry is described in Figure 8.

**Figure 8.** Figure 1. Thorsos' Figure 1 [112]. Scattering geometry. Reproduced with permission from [112], AIP Publishing, 1997.

SSA and IE scattering predictions are presented in Figure 9 for respective k σ and kL 0.5 and 4.0, for 45◦ incident angle and RMS slope s <sup>=</sup> 0.177 (*<sup>s</sup>* <sup>=</sup> <sup>√</sup> 2 σ *<sup>L</sup>* for Gaussian surface), and the RMS slope angle γ = tan−<sup>1</sup> *s* is 10◦. The contribution of specular reflection—leading to a peak in the IE prediction—has not been integrated in the SSA models.

**Figure 9.** Thorsos' Figure 2 [112]. Comparison of the second-order SSA (2), third-order SSA(3), and fourth-order SSA(4) small slope approximation and Monte Carlo integral (IE) scattering strengths (=10 log of the scattering cross section) for a modest rms surface slope angle. Here, k σ = 0.5, kL = 4.0, the RMS slope angle γ is 10◦, and the incident angle θ*<sup>i</sup>* is 45◦. Reproduced with permission from [112], AIP Publishing, 1997.

In Figure 10, k σ and kL have both been increased by a factor of three to 1.5 and 12.0, respectively. The increase in the surface roughness (resp. correlation length) causes disappearance of the specular peak (resp. decrease of the backscattered strength). A higher SSA order improves the prediction in backscattering.

**Figure 10.** Thorsos' Figure 3 [112]. Comparison of SSA(2), SSA(3), SSA(4), Kirchhoff approximation (KA) and IE scattering strengths. Here k σ = 1.5, kL = 12.0, γ is 10◦ and θ*<sup>i</sup>* is 45◦. Reproduced with permission from [112], AIP Publishing, 1997.

According to Thorsos [21], the Kirchhoff approximation KA is accurate away from low grazing angles for rough surfaces smooth on the scale of a wavelength, assuming small surface slopes. For a Gaussian roughness spectrum, the KA is said to be generally accurate for kL >6, that is when *<sup>L</sup>* <sup>λ</sup> > 1. We have seen in Section 3.2.2 the recent definition of a more precise criterion for the KA validity. In the following of this section the SSA and KA are thus considered for several examples with kL > 6.

In Thorsos' previous Figure 3 [112], KA is valid away from grazing angles, whereas PT (not shown) is not. The SSA(2) result reduces to the KA result in the specular region; the third and fourth-order results do as well. KA behaves better than SSA for backscattering, but it is the contrary in the forward direction at grazing angles.

Finally, several cases have been examined beyond the region of validity of the fourth-order Perturbation Theory PT(4) and of KA. For example, in their Figure 6 [112] since k σ = 1, PT(4) breaks down whereas KA is globally correct and SSA(4) the more realistic (Figure 11).

**Figure 11.** Thorsos' Figure 6 [112]. Comparison of the SSA(4), KA, PT(4) and IE scattering strengths for a case when neither PT(4) nor the KA is accurate. Here, k σ = 1.0, kL = 4.3, γ = 18.2◦, and θ*<sup>i</sup>* = 45◦. Reproduced with permission from [112], AIP Publishing, 1997.

Other validations showed that the fourth-order SSA (SSA(4)) mimics some shadowing phenomena and the SSA results are globally accurate for all observation angles for incident angles up to 80◦, but the validity decreases with the incident angle.

Another study [113] of the Bragg scattering from surface waves by a rough interface recently provided comparisons between the Kirchhoff approximation, the perturbation theory and the integral equation method; Kirchhoff leads to results relatively close to the integral equation method in the studied configurations (in forward and backscattering).

Voronovich (see Figure 1.1 in [20]) argued that SSA is valid for σ/*L* less than a constant value. According to him, KA and PT also require this condition but PT needs in addition a small kσ and KA a great kL.

SSA has also been developed for rough solid elastic interfaces [114]. Berman [115] showed that it provides an effective solution compared to Monte-Carlo integral equations simulations for grazing incident angles, whereas Kirchhoff leads to an overall important underestimation. SSA and first-order PT give rather good predictions and exhibit minima or peaks at the Rayleigh angle that can be exaggerated or discontinuities at critical angles which are not physical. The authors noted that Kirchhoff is rather smooth or does not have the valid behaviour for these particular angles: it is logical because the associated phenomena (Rayleigh or head waves) are not accounted on the simplified Kirchhoff theory (see also the perspectives section). We nevertheless notice that Kirchhoff yields an overall rather good prediction versus scattering angle despite some quantitative discrepancies. We also emphasise that one comparison shows an overestimation by first-order PT of the specular amplitude, which is a strong inconvenient for numerous applications as NDE.

SSA for layered fluid seafloors has only been studied recently. A study published this year [116] reviewed three different approximations. Jackson [117] has proposed an acoustic SSA method for rough water–sediment interfaces with a fluid layer of finite thickness overlying a semi-infinite fluid layer. A second method consists in using in acoustics the SSA developed for the general layered case in electromagnetics [118–120]. A third approach is developed in [116] by translating usual small-slope ansatz from k-space to coordinate space. It is then shown in [116] using different examples that the different approximations give different results in the case of strong internal reflections and also differ from Kirchhoff in specular direction. In an example with a mud layer overlying a harder basement, the first approximation differs strongly from the two others and consequently seems erroneous. Unfortunately, no numerical validation is provided to define the best ansatz.

#### *4.2. Parabolic Equation (PE)*

Acoustic losses at the rough ocean surface can be simulated utilizing a propagation code which can use as input an explicitly defined roughness profile. The PE (Parabolic equation [121–123]) model Ramsurf [124] is often used and enables to model shadowing from the surface roughness. Ramsurf is based on the split-step Padé method [121] (a pseudo-spectral numerical method for solving nonlinear partial differential equations).

Firstly, based on the rough ocean surface of Gaussian spectrum, the validity of Ramsurf acoustic propagation model has been studied [125] by comparison with measurements after an explosion at receivers located on the ocean bottom around 1000 m depth. The height for the sea surface is from 1.5 m to 2 m. Transmission loss curve under a smooth sea surface and that under the rough sea surface are simulated respectively by Ramsurf. The sea surface reflection loss caused by the surface roughness is about 4 to 6 dB. The transmission loss simulated by Ramsurf model is very close to the experimental one. Kirchhoff approximation (KA) and small slope approximation (SSA) abilities to simulate the reflection loss are also compared: SSA and Ramsurf match well for very small grazing angles (<6◦) for which KA breaks down. We nevertheless notice that the gap between KA and Ramsurf becomes less than 3 dB above 6◦: KA is thus inaccurate only in a very small region.

Another study [126] focuses on the modeling of underwater acoustical reflection from a wind-roughened ocean surface and compares the surface loss between KA, a small-slope model SSA(2) and a stochastic modeling using Ramsurf. SSA(2) matches well with Ramsurf for the entire range of frequency and wind speed combinations studied (1.5–9 kHz, 5–12.5 m/s) and for sound incidence grazing angles from 1 to 11◦. KA can be used effectively from an angle about 2.5 to 7◦.

An extension of the parabolic equation method to deal with forward scattering has recently been proposed for both Dirichlet and Neumann boundary conditions [127] but future works are needed to make the connection between forward and backward 3D scattering and treat other boundary conditions.

#### **5. Summary, Discussion and Perspectives**

A review of acoustic wave scattering from rough surfaces has been done. Several analytical models proposed in the literature have been described in more detail as well as their validity. We are also focused in specific acoustic applications as scattering from ocean, non-destruction evaluation, medical sciences.

Two classical approximations have been early devised: the perturbation theory and the Kirchhoff approximation.

The perturbation theory (PT) consists in expressing the scattered field in an expansion series using as small parameter the roughness height. This method leads to a major inconvenient: it is valid for low roughness when compared to acoustic wavelength (at first order for small *k*σ < 2 to 4 depending on *kL* value, σ, *L* and *k* being respectively the rms roughness height, the correlation height and the wave number). Existing validations also showed that PT is rather accurate far from the specular direction. Consequently, it appears possibly more restrictive than the Kirchhoff method for numerous applications notably non-destruction evaluation (NDE) which usually employs inspection in specular direction.

The Kirchhoff approximation (KA) leads to a simple analytical expression; it is a relatively powerful and inexpensive approximation in terms of calculations and easily extended to impedance rough surfaces. One important fact is that KA is the most used method in NDE for rough surfaces: its implementation of the rough structure echo consists in simply modifying that of the smooth surface; the "rough" KA algorithm can thus simply mesh the smooth mean surface (surface without roughness) in small surface elements and then modify the contribution of each mesh point by including some local geometrical parameters of the rough surface.

KA has the advantage of handling surfaces of roughness smaller than a constant depending on the wavelength and on the correlation length with a relatively satisfying validity except at the vicinity of grazing observation and critical angles. Indeed we have connected and merged recent studies of the KA validity and conclude that for many applications using rather near specular observation and nongrazing incident angles, two criteria can be used for KA validity: <sup>σ</sup>2/*L* <sup>≤</sup> *C*λp (C and <sup>λ</sup>p being a constant value and the wavelength) and kL >1 or 2. If we do not want to limit the KA validity to near specular observation, the previous criteria are likely to be a little more restrictive to reach the KA validity for incidence and scattering angles over the range from −80◦ to 80◦(kL > π is for instance needed)

It is well admitted that KA validity for application to rough surfaces decreases for grazing angles and the region of inaccuracy in terms of scattering angles is much more important for a fluid/solid interface. One has to pay attention at the KA validity in the vicinity of particular angles (Rayleigh and critical ones). In NDE inspections, these angles are often avoided by the inspector; nevertheless it may be possible to meet them in some geometrical configurations. Improving KA or analytical approximations to better take into account head waves occurring at critical angles is an important challenge already for planar surfaces [128]. Modeling head waves on irregular surfaces is even more complicated and seems intractable using analytical methods [90,129]. Finite elements methods [130] can be useful to model both the head waves generated at critical angles, the bulk waves radiated near the surface and their possible interferences [90,129].

For planar surfaces, some smoothing modifications of KA can be proposed [69] in the vicinity of critical angles for better validity, but they do not simulate head waves. A modification of KA using an alternative integral formula for the Kirchhoff integral has just been developed [89], and it could be extended to scattering from a rough surface. Nevertheless, this perspective has to be tested and validated, since such a model will generate head waves as if they have been created on successive local non-finite tangential plane.

PT and KA are usually seen as complementary, PT being valid for small roughness scales and Kirchhoff for large scales. Several other approximations have been developed in the view of handling the inconveniences of the classical methods (PT and KA) and enlarge the simulation validity. Some of them (phase perturbation theory, two-scales models) found some applications in electromagnetism. A recent study [131] has just began in order to try to improve the two-scale models. We are not convinced that improving these models is the best perspective because they inherently own a certain complexity and there are probably better candidates (evoked in the following) in simulation methods.

The small slope approximation (SSA) is a model based on an integral field equation; it is more efficient than Kirchhoff for grazing angles but is only valid for small surface slopes. Indeed, the SSA is of great importance in ocean applications for modeling scattering from ocean surfaces on which the incidence is generally near from grazing: in such a configuration the SSA accuracy is well established. One main future challenge is to extend this approximation to layered fluid seafloors. In a study published this year [116], different approximations are proposed but they give different results in the case of strong internal reflections and also differ from Kirchhoff in specular direction. Numerical validations

are clearly needed to try to find the best approximation among them; another problem is that only one extended from electromagnetism [118] can deal with multiple layers.

The parabolic equation describes all physical effects including shadowing, but this method needs a complex and heavy numerical technique to solve the parabolic equation. The parabolic equation method has recently been extended to forward scattering, but numerous efforts should be necessary in the future to build a more robust model applicable for both forward and backward 3D scattering. Such a strategy can be seen complex when compared to other methods as fully numerical ones which are generic for all scattering directions.

There is also recent progress in methods based on the integral equation; for instance the new development described in [132] leads to impressive results, notably for grazing angles; but the calculation time is still for the moment a limitation by comparison with simpler methods, if we are rather interested in near normal incidence.

Numerical methods are obviously the main rivals of the existing analytical models. For numerical techniques, obtaining 3D simulations with a reasonable computation time remains a challenge. In a very recent study, the boundary element method has been adapted [133] to reduce the 3D calculation time; the proposed adaptation is validated by using approximated models (KA and PT) showing that the classical approximations are also useful for validation purposes.

In the field of NDE, one major future challenge is to model scattering from roughness in components fabricated by additive manufacturing (AM). The fatigue life of AM specimens is very influenced by their surface topography [134]. AM processes for metallic parts are notably the selective laser melting (SLM), selective electron beam melting (SEBM), laser metal deposition (LMD) and also wire arc additive manufacturing (WAAM [135]). The corresponding roughness value of each process is discussed in [136,137]: it is 4–11 μm Ra for SLM, 15–35 μm Ra for SEBM, 10 μm Ra and 200 μm Ra for LMD, and can be very important for WAAM. Note that models are searched for the prediction of surface roughness (see a list of studies in the introduction of [138]), the optimization of process parameters as well as the roughness minimization (as in [138]) for manufacturing using WAAM and milling.

In AM processes there a lot of phenomena causing roughness [139,140], notably the sticking of nonmelted or partially melted particles on the surfaces and the formation of menisci. The latter phenomenon is more pronounced for the WAAM process (see Figure 12, below).

**Figure 12.** Surface topography of a specimen manufactured by wire arc additive manufacturing (WAAM). Reproduced from [141] (Open access, no copyright).

The balling effect (Figure 13) is a detrimental phenomenon that may result in breaking up the liquid scan track during SLM and produce particles in spherical shapes.

**Figure 13.** Surface topography of a specimen manufactured by selective laser melting (SLM). Reproduced from [137] (Open access, no copyright).

The rough surface of AM parts can thus exhibit a kind of periodicity and elementary egg-shaped irregularities [142]. Finite element models [143] are likely to be the best solution to both model the bulk macrostructure of such components and their surface roughness—such models have been recently used to model scattering in virtual simulated microstructures of polycrystalline materials either with equiaxed grains [144–146] or transverse isotropic as in austenitic welds [147]. The increasing power brought by the spectral finite elements (SEM) method is noticeable [148–150] in terms of calculation time and could be very useful for roughness with complex shape irregularities; SEM has recently been used for modeling ultrasonic propagation in granular materials, such as concrete [151].

In another domain, a FEM model has just been employed to model the ultrasonic response of bone–dental implant interfaces [5,152] whose irregularities can generate multiple scattering. Such a numerical method also presents the advantage of being easily adaptable in the future, for instance to treat non-linear effects at the contact interface. All the previous examples show the increase of applications for numerical methods.

**Author Contributions:** Software, V.D. and M.D.; validation, M.D. and V.D.; investigation, M.D.; resources, M.D.; writing—original draft preparation, M.D.; writing—review and editing, M.D. and V.D.; project administration, F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported through the "Generation IV/Sodium Fast Reactors/ASTRID" program with CEA, EDF, and FRAMATOME funded by the French government.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18