A Kriging surrogate model of an electromagnetic acoustic transducer (EMAT) generating omnidirectional Lamb waves

Abstract

Electromagnetic acoustic transducers (EMATs) are non-contact ultrasonic transducers used to generate ultrasonic waves directly in the tested metal samples, suitable for special situations like testing moving or hot metal solids. EMATs rely on electromagnetic coupling, and the transduction efficiency of EMATs is relatively low, so it’s crucial to improve their performance through building accurate models and implementing structural optimisations. Modern evolutionary optimization algorithms involve a number of evaluations of the objective function of the optimization problem, so it’s necessary to explore surrogate models to reduce the time consumption during evaluation of the computationally expensive objective function, especially when the function corresponds to complex numerical models, like the EMAT model. We focused on the application of Kriging surrogate model (DACE) in EMAT modeling in this work. The formulations of the DACE model were reviewed, then the model of an omnidirectional EMAT was introduced. Three methods to calculate the amplitudes of the displacement components were discussed and compared, including the time-domain method, the multifrequency method (with FFT/IFFT processing), and the single frequency method applying directly the phasors. The DACE model of the EMAT was then built to predict outputs on a discrete grid of the design variables. Finally the obtained DACE model was used as the surrogate model in the optimization of the EMAT driven by the simulated annealing algorithm. With the DACE model, the total time consumption of optimization of the EMAT was reduced greatly.

Keywords

Kriging surrogate model DACE omni-directional electromagnetic acoustic transducer EMAT Lamb waves optimization simulated annealing

1. Introduction

Conventionally, piezoelectric transducers are used in ultrasonic testing to generate ultrasonic waves in the solid sample under investigation. These transducers rely on liquid coupling to transfer the generated waves into the sample, and uncertainty exists because of the coupling. Other limitations include the required preparation of the surface of the sample when it’s rough, and difficulty in testing moving and hot objects.

As a promising alternative to the piezoelectric transducer, electromagnetic acoustic transducer (EMAT) is a kind of non-contact transducer, usually composed of a bias magnet, a coil fed with high frequency alternating current, and the tested metal sample itself. EMATs rely on electromagnetic mechanism to transfer energy directly into the metal sample, without requirement of liquid coupling, and they are becoming widely used in metal material testing, including some special situations like testing moving or hot samples. Another advantage of EMATs is that with different configurations of the bias magnetic field and the coil, various types of ultrasonic waves could be generated. In spite of these benefits, the energy transduction efficiency of EMATs is relatively low compared with that of the piezoelectric transducers. So it’s always imperative to study the mechanism of EMATs further and build accurate and efficient models, and strive to improve their performance.

The model of an EMAT is multiphysics in nature, involving coupling of the electromagnetic and elastodynamic fields. This fact renders modeling of EMATs relatively difficult, and it has been an attractive topic in the previous years. Here we mainly focus on the numerical models. Ludwig and Dai conducted transient analysis of a meander coil EMAT placed on isotropic non-ferromagnetic half-space, assuming uniform static magnetic field [11]. Jafari-Shapoorabadi et al. studied in detail the controlling eddy current equations and argued that the previous work using the total current divided by the cross section area of the conductor as the source current density is equivalently applying the incomplete equation, and this means ignoring the skin effect and the proximity effect [7], while we proved the opposite in [25] via customising the underlying integro-differential or normal differential equations. Shi compared three different definitions of the current density with weak formulations [19]. Garcia-Rodriguez et al. implemented both un-coupled and coupled FEM codes for EMAT simulations in 2D geometries [5]. Dhayalan and Balasubramaniam used the FEM package COMSOL to build the electromagnetic model of a meander EMAT, and the simulated Lorentz force was exported to another package Abaqus as the driving force to excite Lamb waves [2]. The method was later extended to spiral coil EMAT generating bulk waves [3]. The above-mentioned modeling work only involves non-magnetic materials. There is also some initial work on modeling EMATs used to test magnetic material, while we will not discuss further here.

One possible way to improve the performance of EMATs is through optimal design of the structural parameters of the EMATs. The work on optimizations of EMATs is still rare. Mirkhani et al. conducted a parametric study of an EMAT composed of a racetrack coil, by varying the ratio of the width of the magnet to the width of the coil, and found that if this ratio was set at 1.2, the amplitude of the ultrasonic beam would be improved [12]. One design variable and one objective function were used in this optimization, accomplished only through observation of a set of curves corresponding to different design variables, instead of using a real optimization algorithm. Seher et al. optimized a spiral coil EMAT using genetic algorithm optimization procedure in the global optimization toolbox of Matlab [17, 18]. The ratio of the amplitude of the A0 mode Lamb waves to that of the S0 mode waves was selected as the objective function to be maximised, i.e. preferably generating the A0 mode. We extended the problem to multi-objective optimization in [23].

Modern evolutionary optimization algorithms involve a large number of evaluations of the objective function of the optimization problem. Quite often the evaluation of the objective function is time-consuming, for example, when the function corresponds to an underlying numerical model, like the EMAT model. Instead of sticking to the original computationally expensive model, an alternative plan is to build an approximate model from some known input/output data sets first, and then substitute it for the original model in optimizations. The approximate model, called surrogate model, metamodel, etc., should mimic the behavior of the simulation model as closely as possible while being computationally cheaper to evaluate. The surrogate models are already used in electromagnetics problems [22], but still very rare in ultrasonic applications, not to mention our hybrid and multiphysics modeling of EMATs. This situation motivates us to study and apply the surrogate models to the optimization of EMATs.

One particular surrogate model that we’re interested in is Kriging by Danie G. Krige, originating from geostatistics. In this method mathematical statistics was applied to the spatial evaluation of orebodies, which helped improve ore evaluation techniques and reduce the financial risk of investing in mining projects. Kriging surrogate model is flexible because of the wide range of spatial correlation functions available, and it can approximate linear and non-linear functions equally well [20]. The DACE (Design and analysis of computer experiments) is one implementation of Kriging in computer experiments with deterministic outputs, proposed by Sacks et al. [15]. Here ’deterministic’ means that for one set of inputs, no matter how many times we repeat the computer experiment (here the numerical model of the EMAT), it will always generate the same set of outputs. The DACE is data-driven, since some known input/output data pairs are required as prior knowledge.

In this paper, we review the theory of the DACE model first, then summarize the underlying equations of an axisymmetric EMAT working on the Lorentz force mechanism and introduce the specifications of an omni-directional EMAT model (implemented in COMSOL) and possible approaches to calculate the amplitudes of the displacement components, to be used to compute the output of the model. The DACE model is then applied to implement prediction of the output (the objective function). Finally the DACE model is used in optimization of the proposed EMAT driven by the simulated annealing algorithm. The formation of the optimization problem and the parameters of the EMAT are similar as [17, 18], but this is only for convenience of comparison with these previous work. The methodology could be easily applied to other EMATs with different structures and parameters, if some conditions are satisfied.

2. DACE – The Kriging surrogate model

For convenience, we restrict the number of output variables to one, i.e. the model has a scalar output. Supposing $\bm{x}$ , an $n\times 1$ column vector, is one point of the input data. Then in the DACE model, the response corresponding to $\bm{x}$ is,

$\displaystyle Y(\bm{x})=\sum^{k}_{j=0}\beta_{j}f_{j}(\bm{x})+Z(\bm{x})=\bm{f(x% )}^{T}\bm{\beta}+Z(\bm{x})$ (1)

in which $f_{j}(\bm{x}),j=0,\ldots,k$ are chosen functions, $\beta_{j},j=0,\ldots,k$ are the coefficients. $\bm{f}(\bm{x})=(f_{0}(\bm{x}),\ldots,f_{k}(\bm{x}))^{T}$ is the column vector of the functions, while $\bm{\beta}=(\beta_{0},\ldots,\beta_{k})^{T}$ is the column vector of the coefficients. The chosen functions could be one constant function as $f_{0}(\bm{x})=1$ ( $k=0$ ), or a set of linear or quadratic functions. This selection of the chosen functions is not important in the DACE model so in fact we will only use the constant function in this work. In Eq. (1), $\bm{f}(\bm{x})^{T}\bm{\beta}$ is the regression part of the model. $Z(\bm{x})$ is the random process representing the deviation from the regression part. $Y(\bm{x})$ is a transformation of $Z(\bm{x})$ , so it is also a random process. We will use $y(\bm{x})$ to represent the real black box function, i.e. the computer model, and it’s what we want to predict, for any untried input point $\bm{x}$ . $y(\bm{x})$ is one realization of the random process $Y(\bm{x})$ , and this idea is critical in the DACE model.

The random process $Z(\bm{x})$ has zero mean ( $E[Z(\bm{x})]=0$ ) and constant variance $\sigma^{2}$ ( $\textit{Var}[Z(\bm{x})]=\sigma^{2}$ ). The covariance of two samples of the $Z(\bm{x})$ random process, $Z(\bm{w})$ and $Z(\bm{v})$ , is,

$\displaystyle\textit{Cov}(Z(\bm{w}),Z(\bm{v}))=\sigma^{2}R(\bm{w},\bm{v})$ (2)

in which $R(\bm{w},\bm{v})$ is the (spatial) correlation function of $Z(\bm{w})$ and $Z(\bm{v})$ . One popular definition of $R(\bm{w},\bm{v})$ is,

$\displaystyle R(\bm{w},\bm{v})=\prod^{n}_{i=1}\mathrm{e}^{-\theta_{i}\left|w_{% i}-v_{i}\right|^{p_{i}}}$ (3)

in which $n$ is the number of the input variables (dimension of $\bm{w}$ and $\bm{v}$ ), $w_{i}$ and $v_{i}$ are components of $\bm{w}$ and $\bm{v}$ , respectively, and $\theta_{i}$ and $p_{i}$ are the extra parameters of the correlation function. When $p_{i}=1$ , the correlation function is an exponential function. When $p_{i}=2$ , the correlation function, called Gauss function, is relatively smooth and used widely in the literature.

Before prediction of the output at untried input data, suppose the experiment is repeated $m$ times. The input data set is $\bm{S}=(\bm{x}_{1},\ldots,\bm{x}_{m})$ , and the corresponding known output data is $\bm{y}_{S}=(y_{1},\ldots,y_{m})^{T}=(y(\bm{x}_{1}),\ldots,y(\bm{x}_{m}))^{T}$ . $\bm{y}_{S}$ is a realization of the random vector $Y_{\bm{S}}=(Y(\bm{x}_{1}),\ldots,Y(\bm{x}_{m}))^{T}$ , which itself is composed of samples of the $Y(\bm{x})$ process. Corresponding to $Y_{\bm{S}}$ , the samples of the random process $Z(\bm{x})$ at the input data set form the random vector $Z_{\bm{S}}=(Z(\bm{x}_{1}),\ldots,Z(\bm{x}_{m}))^{T}$ .

Evaluated at the input data set $\bm{S}$ , the vector of chosen functions $\bm{f}(\bm{x})$ forms the design matrix $\mathbf{F}$ as,

$\displaystyle\mathbf{F}=\begin{pmatrix}\bm{f}(\bm{x}_{1})^{T}\\ \vdots\\ \bm{f}(\bm{x}_{m})^{T}\end{pmatrix}=\begin{pmatrix}f_{0}(\bm{x}_{1})&\cdots&f_% {k}(\bm{x}_{1})\\ \vdots&\vdots&\vdots\\ f_{0}(\bm{x}_{m})&\cdots&f_{k}(\bm{x}_{m})\\ \end{pmatrix}$ (4)

Then from Eq. (1), we have $Y_{\bm{S}}=\mathbf{F}\bm{\beta}+Z_{\bm{S}}$ .

Suppose the predictor is linear,

$\displaystyle\hat{y}(\bm{x})=\bm{c}(\bm{x})^{T}\bm{y_{S}}$ (5)

in which $\bm{c}(\bm{x})$ is $m\times 1$ column vector function of $\bm{x}$ .

With the random vectors ( $Y_{\bm{S}}$ and $Z_{\bm{S}}$ ) and processes ( $Y(\bm{x})$ and $Z(\bm{x})$ ), the error of the prediction is,

$\displaystyle\bm{c}(\bm{x})^{T}Y_{\bm{S}}-Y(\bm{x})=\left[Z_{\bm{S}}^{T}\bm{c}% (\bm{x})-Z(\bm{x})\right]+\left(\mathbf{F}^{T}\bm{c}(\bm{x})-\bm{f(x)}\right)^% {T}\bm{\beta}$ (6)

The mean squared error (MSE) is,

$\displaystyle\mbox{MSE}=E\left[\left(\bm{c}(\bm{x})^{T}Y_{\bm{S}}-Y(\bm{x})% \right)^{2}\right]$ (7)

The unbiasedness condition of the predictor requires that $E\left[\bm{c}(\bm{x})^{T}Y_{\bm{S}}-Y(\bm{x})\right]=0$ [16], also, we have $E[Z(\bm{x})]=0$ , so the unbiasedness condition is,

$\displaystyle\mathbf{F}^{T}\bm{c}(\bm{x})-\bm{f(x)}=\bm{0}$ (8)

With this condition, the MSE is simplified as,

$\displaystyle\text{MSE}=E\left[\left(Z_{\bm{S}}^{T}\bm{c}(\bm{x})-Z(\bm{x})% \right)^{2}\right]=\bm{c}(\bm{x})^{T}E\left[Z_{\bm{S}}Z_{\bm{S}}^{T}\right]\bm% {c}(\bm{x})+E\left[Z^{2}(\bm{x})\right]-2E\left[Z(\bm{x})Z_{\bm{S}}^{T}\right]% \bm{c}(\bm{x})$ (9)

If $\mathbf{R}$ is the correlation matrix of the known input data set with the element $R_{ij}$ as,

$\displaystyle R_{ij}=R(\bm{x}_{i},\bm{x}_{j}),\ i,j=1,\ldots m$ (10)

and the correlation vector between the input data set and $\bm{x}$ ,

$\displaystyle\bm{r}(\bm{x})=\left(R(\bm{x}_{1},\bm{x}),\ldots,R(\bm{x}_{m},\bm% {x})\right)^{T}$ (11)

then the MSE is,

$\displaystyle\mbox{MSE}=\sigma^{2}\left(1+\bm{c}(\bm{x})^{T}\mathbf{R}\bm{c}(% \bm{x})-2\bm{r}(\bm{x})^{T}\bm{c}(\bm{x})\right)$ (12)

To minimize the MSE with respect to $\bm{c}(\bm{x})$ , with equality constraint Eq. (8), the method of Lagrange multipliers is used. Suppose $\bm{\lambda}$ is the column vector of Lagrange multipliers, the Lagrange function is,

$\displaystyle L(\bm{c},\bm{\lambda})=\sigma^{2}\left(1+\bm{c}(\bm{x})^{T}% \mathbf{R}\bm{c}(\bm{x})-2\bm{r}(\bm{x})^{T}\bm{c}(\bm{x})\right)-\bm{\lambda}% ^{T}\left(\mathbf{F}^{T}\bm{c}(\bm{x})-\bm{f(x)}\right)$ (13)

The gradient vector with respect to $\bm{c}(\bm{x})$ is,

$\displaystyle\nabla L(\bm{c},\bm{\lambda})=2\sigma^{2}\left(\mathbf{R}\bm{c}(% \bm{x})-\bm{r}(\bm{x})\right)-\mathbf{F}\bm{\lambda}$ (14)

The necessary condition for optimality is that this gradient is equal to zero vector, together with the unbiasedness constraint, we have,

$\displaystyle\begin{pmatrix}\mathbf{R}&\mathbf{F}\\ \mathbf{F}^{T}&\bm{0}\end{pmatrix}\begin{pmatrix}\bm{c}(\bm{x})\\ \hat{\bm{\lambda}}\end{pmatrix}=\begin{pmatrix}\bm{r}(\bm{x})\\ \bm{f(x)}\end{pmatrix}$ (15)

in which $\hat{\bm{\lambda}}=-\frac{\bm{\lambda}}{2\sigma^{2}}$ . The solution to this system of equations is,

$\displaystyle\hat{\bm{\lambda}}=\left(\mathbf{F}^{T}\mathbf{R}^{-1}\mathbf{F}% \right)^{-1}\left(\mathbf{F}^{T}\mathbf{R}^{-1}\bm{r}(\bm{x})-\bm{f(x)}\right)% ,\bm{c}(\bm{x})=\mathbf{R}^{-1}\left(\bm{r}(\bm{x})-\mathbf{F}\hat{\bm{\lambda% }}\right)$ (16)

Then the predictor in Eq. (5) is,

$\displaystyle\hat{y}(\bm{x})=\bm{c}(\bm{x})^{T}\bm{y_{S}}=\bm{f(x)}^{T}\hat{% \bm{\beta}}+\bm{r}(\bm{x})^{T}\mathbf{R}^{-1}\left(\bm{y_{S}}-\mathbf{F}\hat{% \bm{\beta}}\right)$ (17)

in which,

$\displaystyle\hat{\bm{\beta}}=\left(\mathbf{F}^{T}\mathbf{R}^{-1}\mathbf{F}% \right)^{-1}\mathbf{F}^{T}\mathbf{R}^{-1}\bm{y_{S}}$ (18)

Now the MSE is,

$\displaystyle\mbox{MSE}=\sigma^{2}\left[1+\bm{u(x)}^{T}\left(\mathbf{F}^{T}% \mathbf{R}^{-1}\mathbf{F}\right)^{-1}\bm{u(x)}-\bm{r}(\bm{x})^{T}\mathbf{R}^{-% 1}\bm{r}(\bm{x})\right]$ (19)

in which $\bm{u(x)}=\mathbf{F}^{T}\mathbf{R}^{-1}\bm{r}(\bm{x})-\bm{f(x)}$ .

Assuming $Z(\bm{x})$ is a Gaussian random process, the joint PDF of the random vector $Z_{\bm{S}}$ is a multivariate Gaussian PDF [9],

$\displaystyle p_{Z_{\bm{S}}}(\bm{z})=\frac{1}{(2\pi\sigma^{2})^{m/2}\mbox{det}% ^{1/2}(\mathbf{R})}\mathrm{e}^{-\frac{1}{2\sigma^{2}}\bm{z}^{T}\mathbf{R}^{-1}% \bm{z}}$ (20)

in which $\mbox{det}(\mathbf{R})$ is the determinant of $\mathbf{R}$ . In short hand, this PDF could be written as $Z_{\bm{S}}\sim\mathcal{N}(\bm{\mu},\mathbf{C})$ where the mean vector is $\bm{\mu}=\bm{0}$ and the covariance matrix is $\mathbf{C}=\sigma^{2}\mathbf{R}$ . $Y_{\bm{S}}$ is a transformation of $Z_{\bm{S}}$ , so the joint PDF of $Y_{\bm{S}}$ is,

$\displaystyle p_{Y_{\bm{S}}}(\bm{y})=\frac{1}{(2\pi\sigma^{2})^{m/2}\mbox{det}% ^{1/2}(\mathbf{R})}\mathrm{e}^{-\frac{1}{2\sigma^{2}}(\bm{y}-\mathbf{F}\bm{% \beta})^{T}\mathbf{R}^{-1}(\bm{y}-\mathbf{F}\bm{\beta})}$ (21)

If the PDF $p_{Y_{\bm{S}}}(\bm{y})$ is viewed as a function of the unknown parameters $\bm{\beta}$ and $\sigma^{2}$ , it’s the likelihood function [8]. Then maximum likehood estimations (MLEs) of the parameters could be solved via differentiating the log-likelihood function $\ln p_{Y_{\bm{S}}}(\bm{y})$ with respect to $\bm{\beta}$ and $\sigma^{2}$ and setting the derivatives to zero.

$\displaystyle\frac{\partial\ln p_{Y_{\bm{S}}}(\bm{y})}{\partial\bm{\beta}}=\bm% {0}$ (22)

gives the same result as Eq. (18). While

$\displaystyle\frac{\partial\ln p_{Y_{\bm{S}}}(\bm{y})}{\partial\sigma^{2}}=0$ (23)

leads to

$\displaystyle\hat{\sigma^{2}}=\frac{1}{m}(\bm{y}_{\bm{S}}-\mathbf{F}\hat{\bm{% \beta}})^{T}\mathbf{R}^{-1}(\bm{y}_{\bm{S}}-\mathbf{F}\hat{\bm{\beta}})$ (24)

For implementation of the DACE model, we choose a Matlab toolbox bearing the same name of ‘DACE’, developed by Lophaven et al. [14].

3. Model of an omni-directional EMAT

3.1 Formulations of the EMATs

The basic equations for the electromagnetic field simulation in the study of EMATs are Maxwell’s equations [6],

$\displaystyle\nabla\times\mathbf{E}=-\frac{\partial{\mathbf{B}}}{\partial{t}}$ (25) $\displaystyle\nabla\times\mathbf{H}=\mathbf{J}_{f}+\frac{\partial{\mathbf{D}}}% {\partial{t}}$ (26) $\displaystyle\nabla\cdot\mathbf{D}=\rho_{f}$ (27) $\displaystyle\nabla\cdot\mathbf{B}=0$ (28)

These equations are Faraday’s law, Ampere-Maxwell law, Gauss’s law for electric fields and Gauss’s law for magnetic fields, respectively. $\mathbf{E}$ is the electric field, $\mathbf{B}$ is the magnetic flux density, $t$ is the time variable, $\mathbf{H}$ is the magnetic field strength, $\mathbf{J}_{f}$ is the free current density, $\mathbf{D}$ is the electric flux density, $\rho_{f}$ is the free charge density. $\nabla\times$ is curl of a vector, and $\nabla\cdot$ is divergence of a vector. Since the frequency in EMAT operation is no higher than several MHz, the term $\frac{\partial{\mathbf{D}}}{\partial{t}}$ (the displacement current density) in Maxwell’s equations could be neglected.

To solve Maxwell’s equations, another set of equations called the constitutive equations is needed,

$\displaystyle\mathbf{B}=\mu\mathbf{H}$ (29) $\displaystyle\mathbf{D}=\epsilon\mathbf{E}$ (30)

in which $\mu$ is the magnetic permeability, and $\epsilon$ is the dielectric constant.

With the magnetic vector potential (MVP) $\mathbf{A}$ defined as $\mathbf{B}=\nabla\times\mathbf{A}$ , the equation describing the eddy current phenomenon is,

$\displaystyle-\frac{1}{\mu}\nabla^{2}\mathbf{A}+\sigma\frac{\partial\mathbf{A}% }{\partial t}=\mathbf{J}_{s}$ (31)

in which $\nabla^{2}$ is the vector Laplacian operator, $\sigma$ is the conductivity, and $\mathbf{J}_{s}$ is the source current density. The eddy current density, not written explicitly in Eq. (31), is $\mathbf{J}_{e}=-\sigma\frac{\partial\mathbf{A}}{\partial t}$ . This equation holds where there exists a conductor and a source current flows inside this conductor. At the region without the source current, the $\mathbf{J}_{s}$ term is dropped. Where there is no conductor, like in the air, the negative eddy current density term $\sigma\frac{\partial\mathbf{A}}{\partial t}$ (or $-\mathbf{J}_{e}$ ) is also dropped.

For 2D axisymmetric field simplification assuming $\mathbf{J}_{s}$ is perpendicular to the $r$ - $z$ plane, we have $\mathbf{J}_{s}=J_{s\phi}\mathbf{e}_{\phi}$ , with $J_{s\phi}$ as the $\phi$ component of vector $\mathbf{J}_{s}$ , and $\mathbf{e}_{\phi}$ as the unit vector along the $\phi$ axis. We also have $\mathbf{A}=A_{\phi}\mathbf{e}_{\phi}$ . From now on, $A_{\phi}$ will be written as $A$ and $J_{s\phi}$ be written as $J_{s}$ for simplicity. With this notation, the vector Eq. (31) is transformed to the $\phi$ component scalar equation,

$\displaystyle-\frac{1}{\mu}\left(\frac{\partial^{2}A}{\partial r^{2}}+\frac{% \partial^{2}A}{\partial z^{2}}+\frac{1}{r}\frac{\partial A}{\partial r}-\frac{% A}{r^{2}}\right)+\sigma\frac{\partial A}{\partial t}=J_{s}$ (32)

This equation is a diffusion equation describing the eddy current phenomenon in cylindrical coordinates. Now the eddy current density term is $J_{e}=-\sigma\frac{\partial A}{\partial t}$ .

With total current $i$ , the original eddy current Eq. (32) in cylindrical coordinates could now be written as [24],

in which $S$ is the cross section of the source conductor.

For steady state or frequency-domain analysis, the phasor notation is adopted,

$\displaystyle-\frac{1}{\mu}\left(\frac{\partial^{2}\dot{A}}{\partial r^{2}}+% \frac{\partial^{2}\dot{A}}{\partial z^{2}}+\frac{1}{r}\frac{\partial\dot{A}}{% \partial r}-\frac{\dot{A}}{r^{2}}\right)+\mbox{j}\omega\sigma\dot{A}=\frac{% \dot{i}+\mbox{j}\omega\sigma\iint_{S}\dot{A}dS}{r\iint_{S}\frac{1}{r}dS}$ (34)

in which the dots on $A$ and $i$ indicate they are complex phasors. $\omega$ is the angular frequency. j is the imaginary unit. This steady-state equation was proposed by Preis [13], following the work of Konrad on the similar steady state integro-differential equation in a 2D planar model [10]. These equations are solved in this work with a mid-level approach as introduced in [24] so the underlying physics could be controlled completely.

The equations describing the generation and propagation of the ultrasonic waves in an elastic solid are [1],

$\displaystyle\nabla\cdot\mathbf{T}=\rho\frac{\partial^{2}{\mathbf{u}}}{% \partial{t^{2}}}-\mathbf{F}$ (35) $\displaystyle\mathbf{T}=\mathbf{c}:\mathbf{S}$ (36) $\displaystyle\mathbf{S}=\nabla_{s}\mathbf{u}$ (37)

These equations are equation of motion, Hooke’s law and strain-displacement relation respectively. $\mathbf{T}$ is the stress tensor, $\nabla\cdot\mathbf{T}$ is its divergence, $\rho$ is the density, $\mathbf{F}$ is the body force density, $\mathbf{c}$ is the stiffness tensor, $\mathbf{S}$ is the strain tensor, and $\mathbf{u}$ is the displacement vector. $:$ is the double dot product of a fourth rank tensor $\mathbf{c}$ and a second rank tensor $\mathbf{S}$ . $\nabla_{s}\mathbf{u}$ is the symmetric part of the gradient of the vector $\mathbf{u}$ . It should be noted that the same symbol $\rho$ was also used to represent charge density in Maxwell’s equations, and this shouldn’t cause misunderstanding in this paper.

For homogenous and isotropic media, we have Navier’s equation,

$\displaystyle\mu\nabla^{2}\mathbf{u}+\left(\lambda+\mu\right)\nabla\left(% \nabla\cdot\mathbf{u}\right)=\rho\frac{\partial^{2}{\mathbf{u}}}{\partial{t^{2% }}}-\mathbf{F}$ (38)

Here $\lambda$ and $\mu$ are Lamé constants. Like the symbol $\rho$ , the symbol $\mu$ was also used to represent the magnetic permeability.

The link between the electromagnetic field and the elastodynamic field is the volume Lorentz force density defined as,

$\displaystyle\mathbf{F}_{L}=\mathbf{J}\times\mathbf{B}=\mathbf{J}\times\left(% \mathbf{B}_{0}+\mathbf{B}_{d}\right)$ (39)

in which $\mathbf{B}$ is the total magnetic flux density composed of the static flux density $\mathbf{B}_{0}$ of the bias magnet and the dynamic flux density $\mathbf{B}_{d}$ generated by the excitation coil. Usually for the excitation current with moderate magnitude, the $\mathbf{B}_{d}$ term is very small compared with $\mathbf{B}_{0}$ , so it is ignored in this work.

Further details of the formulations given here could be found in [24].

3.2 Specifications of the model

In this work we consider an omnidirectional EMAT composed of a spiral coil and a cylindrical permanent magnet, to generate Lamb waves in a plate. The coil is composed of tightly wound copper wires, so both S0 mode and A0 mode Lamb waves will be generated. In this work, the aim is to generate A0 mode Lamb waves, so we build the model bearing this preference in mind.

The complete EMAT model is composed of one magnetostatic sub-model describing the magnetic field of the permanent magnet, one eddy current sub-model analysing the eddy current phenomenon accompanied by the skin and proximity effects, and one elastodynamic sub-model for the simulation of the generation and propagation of Lamb waves in the plate. The two electromagnetic sub-models share one geometry containing the air, the inner section of the plate, the copper wires, and the permanent magnet, as in Fig. 1. Note that in this geometry only the inner section of the full plate is modeled. The elastodynamic sub-model has its own geometry, only containing the full plate. The Lorentz force calculated from the two electromagnetic sub-models is transferred to the elastodynamic sub-model as the driving force of Lamb waves.1 With these two geometries, the structure of the model is very clear. Additionally, we can use different meshing rules for these two geometries, according to the respective physics. This two-geometry treatment is valid because the Lorentz force is local in the region of the plate just under the transducer.

Figure 1.

The geometry of the electromagnetic sub-models. This geometry is used for magnetostatic analysis and eddy current analysis.

In Fig. 1, it’s only necessary to consider the region where $r>0$ , since this is an axisymmetric model. The testing frequency is 50 kHz. As in [17], the relative magnetic permeability of the steel plate is 160, i.e. it’s a simplified linear material.2 The conductivity of the plate is 4.032 MS/m. The thickness of the plate is 10 mm. The remanent magnetic flux density of the magnet is set to 1.3 T along the positive direction of the $z$ axis. $R_{M}$ is the radius of the magnet. $l_{M}$ is the liftoff distance of the magnet from its bottom to the top of the coil. The two parameters of $R_{M}$ and $l_{M}$ will be used as the design variables in the optimizations based on the surrogate model, while all the other parameters are fixed.

$R_{C}$ is the average radius of the coil decided as,

$\displaystyle R_{C}=(2n-1)\frac{\lambda}{4},n=1,2,\ldots$ (40)

in which $\lambda$ is the wavelength of the desired Lamb wave mode.3 For the EMAT on the steel plate, $n$ is chosen to be 1, i.e. $R_{C}=\frac{\lambda}{4}$ , similar as in [18]. As stated previously, we want to selectively generate A0 mode Lamb waves. For A0 mode Lamb waves in a steel plate of 10 mm thickness at 50 kHz, from the dispersion curves generated with a program we developed, the phase velocity is 1867.78 m/s, then the wavelength is $\lambda_{A0}=$ 37.36 mm.

$W_{C}$ is the radial width of the coil (difference between the outer and inner radii of the coil). The coil is composed of two layers of copper wires with conductivity as 5.998 $\times$ 10 ${}^{7}$ S/m. The wires form an array of 23 columns and 2 rows, as in [26]. The wires have rectangular cross sections. The radial width of each wire is 0.3 mm, and the radial gap between adjacent wires in the same layer is 0.1 mm. The axial height of the wire and the gap between the two layers are both 0.1 mm. We have chosen to model each wire individually for a reason to be discussed later.

The geometry of the elastodynamic sub-model simply contains a full plate with a radius of 1.2 m. Young’s modulus is 200 $\times$ 10 ${}^{9}$ Pa, Poisson’s ratio is 0.33, and the density is 7850 kg/m ${}^{3}$ . The observation point to record the displacement components in the simulations is located at 60 cm from the $z$ axis, in the middle plane of the plate. From the displacement wave structures of Lamb waves with the specified frequency and plate thickness, at the middle plane of the plate, the displacement component $u=u_{r}$ only corresponds to the S0 mode, while the other component $w=u_{z}$ only corresponds to the A0 mode.

The boundaries of the sub-models must be handled with care. In the geometry for the electromagnetic sub-models, there is a layer of infinite elements at the air boundary simulating an air region extending to infinitely far away. In the geometry of the elastodynamic sub-model, the top and bottom boundaries of the full plate are free boundaries without constraints or loads. For a transient analysis, the outer end edge of the plate (at $r=$ 1.2 m) is also a free boundary. If the full plate is long enough in the radial direction and the total time of simulation is limited to a proper value, the reflections from the end of the plate can be avoided. While for a frequency domain analysis, an extra perfectly matched layer (PML) must be added to the end of the plate so that the energy in the plate can dissipate.

To further increase the accuracy of the model, fillets are added to the sharp corners of the magnet and the wires, so that the singularities are removed, while at the same time the number of elements and hence the complexity of the model is also increased.

4. Approaches to evaluate the amplitudes of the displacement components

The amplitudes of the displacement components at the observation point will be used to calculate the objective function of optimization. So we must explore how to calculate the amplitudes.

The first approach is evaluating the amplitudes via the time-domain model. The bias magnetic field comes from the magnetostatic simulation. The eddy current distribution, and the generation and propagation of the Lamb waves are from the time-dependent simulations. An implicit time-stepping scheme (generalized- $\alpha$ ) is used for this simulation. To ensure the convergence of the time-dependent solver in COMSOL, a very small time step must be used, which means the simulation will be time-consuming. In this work, the number of time steps is usually set as 6000, for the tone-burst excitation signal $x(t)$ composed of 5 sinusoidal periods modulated with a Hanning window function. Once the time waveforms $u(t)$ and $w(t)$ at the observation point are simulated, the amplitudes/peaks of the envelops of these waveforms will be solved.

The second approach is to transform the input time-continuous excitation signal $x(t)$ to its frequency components via the Fourier transform (implemented with the FFT on a computer), feed them into a frequency domain model, transform the output back into the time-response with the inverse Fourier transform (implemented with the IFFT), and finally solve the peaks of the envelopes of the time waveforms. The time waveforms obtained in this way can be expressed as,

$\displaystyle u(t)=\mathcal{F}^{-1}\left\{\mathcal{F}[x(t)]H_{u}(\omega,R_{M},% l_{M})\right\}$ (41) $\displaystyle w(t)=\mathcal{F}^{-1}\left\{\mathcal{F}[x(t)]H_{w}(\omega,R_{M},% l_{M})\right\}$ (42)

in which $\mathcal{F}$ represents the Fourier transform, $\mathcal{F}^{-1}$ is the inverse Fourier transform, $x(t)$ is the input tone burst signal, $H_{u}(\omega,R_{M},l_{M})$ is the system function for the displacement component $u$ , and $H_{w}(\omega,R_{M},l_{M})$ is the system function for the displacement component $w$ . $R_{M}$ and $l_{M}$ are included to stress that these system functions change with the design variables, while the input signal $x(t)$ is fixed. From the time waveforms the amplitudes/peaks could be solved just like in the first approach. Because the spectrum of the input burst signal is concentrated around the center frequency, we can select only the frequency components bigger than some threshold value as an acceptable approximation. Normally tens of (or even fewer) frequency components are enough, so this approach will be less time-consuming than the first approach applying time-domain simulation. We build a frequency-domain model of the EMAT, in which the bias magnetic field is again from the magnetostatic simulation, but the eddy current sub-model and the elastodynamic sub-model are completely in the frequency-domain. Then we implement this approach by connecting the frequency-domain model in COMSOL with the Matlab environment, in which the FFT/IFFT processing are accomplished. The time waveforms from this approach are carefully compared with the waveforms from the previous time domain simulations, which serve as a reference. From many test simulations, it was found that not every type of current or current density excitation in COMSOL could satisfy our requirement that the time waveforms from the second approach be the same as those from the first approach. For this reason, we choose to model every wire of the coil individually by specifying the total current in the wire, which is more complex than specifying other types of excitations. This is necessary because it satisfies our requirement. As an example, $u$ and $w$ waveforms from the time-dependent simulation and the frequency domain model with FFT/IFFT processing are compared in Fig. 2. The design variables are selected as $R_{M}=$ 8 mm and $l_{M}=$ 1 mm. The threshold to select the frequency components is 10%, that is, only the frequency components higher than 10% of the peak value of the spectrum are used, and others are discarded. With this threshold, 11 components around the center frequency are kept.

One important requisite to validate the frequency domain model is that the whole model must be linear. With the frequency domain model we cannot handle the dynamic part of the Lorentz force easily, because two complex phasors cannot be multiplied to obtain another phasor, so we have to specify a small value of the input current so that the Lorentz force component from the dynamic magnetic field could be ignored.

Since we are mostly concerned with the amplitudes of the $u$ and $w$ waveforms, yet another approach exists, where only one single frequency is used in the frequency model. That is, we only consider the center frequency of the burst signal (50 kHz for the omni-directional EMAT), and use the absolute values of the complex phasors to approximate the amplitudes of the waveforms. This process could be formulated as,

$\displaystyle|\dot{u}|=|H_{u}(\omega_{c},R_{M},l_{M})|$ (43) $\displaystyle|\dot{w}|=|H_{w}(\omega_{c},R_{M},l_{M})|$ (44)

in which $\omega_{c}$ is the center frequency in radian, $\dot{u}$ is the complex phasor of $u$ , and $\dot{w}$ is the complex phasor of $w$ .

We have to further prove that the amplitudes from the single frequency approach (Approach 3) are acceptable approximations of those from the multifrequency approach (Approach 2). From Eq. (43), we can see that $|\dot{u}|$ is the system function evaluated at the center frequency. While in Eq. (41), the spectrum of the tone burst signal $\mathcal{F}[x(t)]$ is bell-shaped, i.e. narrow-banded, and doesn’t change with the design variables, so the result of $\mathcal{F}^{-1}[\cdot]$ operation is roughly proportional to the value of the system function $H_{u}(\omega,R_{M},l_{M})$ at the center frequency $\omega_{c}$ , if the system function is smooth (slowly changing with frequency) with respect to the spectrum of the input burst signal. Then a bigger $H_{u}(\omega_{c},R_{M},l_{M})$ means bigger magnitude of the time waveform, and thus bigger peak value of its envelope. So the phasors could be used to approximately compute the objective functions. Corresponding to Fig. 2, the $u$ and $w$ system functions are shown in Fig. 3. These functions are indeed slowly changing with frequency compared with the spectrum of the input signal which is concentrated around 50 kHz.

The amplitudes of $u$ , $w$ and their ratio from Approaches 2 and 3 are solved numerically with fixed $l_{M}=1$ mm and different $R_{M}$ values, to further validate Approach 3. The results are shown in Fig. 4. $p_{u}$ is the peak value of the envelope of $u$ waveform solved with the frequency domain model and FFT/IFFT processing, while $|\dot{u}|$ is the magnitude of $u$ phasor solved with one single frequency in the frequency domain model. It could be observed that the amplitudes from these two approaches are similar. In fact, the curves from the two approaches are not required to be the same. What’s important is that they have similar shapes and reach respective optimal values at the same set of design variables.

Figure 2.

$u$ and $w$ waveforms from time-dependent simulation and frequency domain model with FFT/IFFT processing.

The third approach is the fastest, since only one frequency is used.

5. Prediction of the objective function

In the previous sections, we reviewed the DACE model as an implementation of the Kriging surrogate model for deterministic computer experiments, built an omnidirectional EMAT model, discussed the possible approaches to compute the amplitudes of the displacement components. Next, we can apply the DACE model to predict the objective function, as a preparation for the optimizations.

Firstly, to test the performance of the prediction, we use the single frequency approach to calculate the amplitudes, because it’s the least time-consuming. We solved the objective function ( $|\dot{u}|/|\dot{w}|$ ) on a 60 $\times$ 60 discrete grid of the design variables as a reference. The contour plot is shown at left in Fig. 5. The surface plot is also shown at right in the figure so that its shape is clearer. There is a valley in the plots approximately along the $l_{M}$ axis. The optimal point should be in this valley.

Figure 3.

$u$ and $w$ system functions.

Figure 4.

$u$ , $w$ amplitudes and their ratio at different $R_{M}$ values, from Approaches 2 and 3. $l_{M}=1$ mm is fixed.

Next, we generate different number of known input/output data pairs according to Latin Hypercube Sampling [15], build the DACE model with these known data, and predict the objective function using the DACE model on a 60 $\times$ 60 discrete grid of the design variables. The exponential correlation function is chosen ( $p_{i}=1$ in Eq. (3)), instead of the commonly used Gauss function, because with the exponential correlation function, the predicted data are more close to the true data. In Fig. 6 are the contour plots of the predicted objective function, using the single frequency model, with different values of $m$ , the number of the known input/output data pairs. It’s obvious that with increasing $m$ , the performance of prediction is improving, since more data are used as known information for prediction.

Figure 5.

Contour plots of the predicted objective function, single frequency model, on a 60 $\times$ 60 discrete grid of the design variables with different $m$ .

Figure 6.

Contour plots of the predicted objective function, multi-frequency model, on a 60 $\times$ 60 discrete grid of the design variables with different $m$ .

Figure 7.

Contour plots of the predicted objective function, multi-frequency model, on a 60 $\times$ 60 discrete grid of the design variables with different $m$ .

Then we explore prediction with the multifrequency model (Approach 2, frequency-domain model with FFT/IFFT processing). This method is more time-consuming, so we only build the DACE model and ignore the true model data like in the contour and surface plots in Fig. 5 which require 3600 evaluations of the objective function. The results are shown in Fig. 7, with different values of $m$ , the number of known input/output data pairs. Clearly the shape of the predicted data on the grid of design variables is similar to that in Fig. 6, showing a valley. From the figure we can see that when $m=$ 800, the prediction should be accurate enough to be used in the optimization. For the case of $m=$ 800, the total time to generate the known data is about 109927 s (30.5 hours), on a computer installed with Intel Core i7-3770 CPU @ 3.40 GHz, a RAM of 32 GB, running x64 version of Windows 10 operating system. Note that with each pair of input/output data, 11 frequency components are used, as in Section 4. With these known data, the time to build the DACE model is 1.249 s, and the time to predict the data on the grid of design variables is 1.208 s, on the same computer.

Figure 8.

MSE for prediction, multi-frequency model, $m$ =800.

With the DACE Kriging model, we can also evaluate the MSE directly. The MSE for prediction using the multi-frequency model with $m=$ 800 is shown in Fig. 8. From various tests, we found that with increasing $m$ value, the magnitude of the MSE is decreasing, meaning the prediction is more accurate with bigger $m$ values.

6. Optimizition of EMAT based on the Kriging surrogate model

Finally we consider the problem of optimization of the EMAT. Formally, the problem is stated as,

$\displaystyle\text{minimize}\ \ f(R_{M},l_{M})=\frac{A_{u}}{A_{w}}$ (45)

in which $A_{u}$ is the amplitude of $u$ and $A_{w}$ is the amplitude of $w$ . They could be regarded as the amplitudes of the S0 and A0 modes respectively. So the form of the optimization problem indicates that the A0 mode is preferred with respect to the S0 mode, representing a tendency in the study of Lamb waves [4]. The design variables $R_{M}$ and $l_{M}$ have upper and lower bounds as, $R_{M}\in[2.5,15]$ mm, $l_{M}\in[1,3]$ mm.

In this section, we apply the surrogate model obtained in the previous section for optimization of the omni-directional EMAT. The specific optimization algorithm selected is simulated annealing (SA) [21]. The idea of SA came from the process of cooling molten materials to solid state, where the system tends to reach the minimum energy state. The energy of a particle jumps randomly, as governed by the temperature of the system composed of many particles. The probability $P$ of transition from energy state $e_{i}$ to energy state $e_{j}$ with temperature $T$ is,

$\displaystyle P\left(e_{i},e_{j},T\right)=\mathrm{e}^{\left(e_{i}-e_{j}\right)% /\left(k_{B}T\right)}$ (46)

in which $k_{B}$ is Boltzmann’s constant. If energy $e_{i}$ is higher than $e_{j}$ ( $e_{i}>e_{j}$ ), we have $P>1$ , i.e. the transition is adopted whatever since the energy gets lower. While if $e_{i}$ is lower than $e_{j}$ ( $e_{i}<e_{j}$ ), we have $0<P<1$ , although this transition is not welcome, it’s still possible that the transition will be granted, as decided by the probability $P$ . In this way, hopefully we will not be trapped at local optimal solutions, so that a global optimum could be found. Also note that with higher temperature $T$ , the transition to a higher energy state is more likely to occur, so that in the early stage of the optimization, a more global search is encouraged, i.e. we want to explore the search space. The temperature $T$ is updated according to the relation $T_{k+1}=\alpha T_{k}$ in which $\alpha$ is the cooling fraction.

We developed a SA program in Matlab for optimization of the EMAT with the Kriging surrogate model (the DACE model), applying the multifrequency model (Approach 2). The initial temperature is 1, Boltzmann’s constant is $k_{B}=$ 0.01, and the cooling fraction is $\alpha=$ 0.97. The number of cooling steps is 100, the number of steps of each temperature is 100, so the number of evaluations of the objective function is 10000 for each run of the program, for which the total time consumption is 1.7 s, on the computer described in the previous section. Without the surrogate model, the total time would be around 388 hours.

We run the SA programs multiple times. The reason is that the SA algorithm is stochastic so generally with every run the obtained result is slightly different. While for 5 runs of the program solving our EMAT optimization problem, we always get the same results. The minimum of the objective function is 0.0057, the design variables are $R_{M}=$ 10.74 mm and $l_{M}=$ 1.14 mm. In [18], the optimization result was $2R_{M}=$ 21.05 mm and $l_{M}=$ 1.47 mm, which is close to our result obtained with the Kriging surrogate model (DACE). This proves the effectiveness of the surrogate model used in the optimization of EMATs. With the Kriging surrogate model, the total time consumpition is greatly decreased, while at the same time acceptable accuracy is achieved for the optimization problem.

Finally we can summarize the strategy to optimize an EMAT, possibly with a surrogate model,

The time domain model is too time-consuming. One run of the model normally takes several hours. Even with the surrogate model, the time to generate enough known data pairs is too much, so generally we do not want to use it in optimizations.

The multifrequency method is much faster. Without a proper surrogate model, the total time of optimization will be too long. With the surrogate model, although some time is needed to generate the known data, the total time is greatly decreased and acceptable. One requisite is that the model must be linear or a frequency model cannot be applied.

The single frequency method, using the phasors diretly, is the fastest. Even without a surrogate model, the total time of optimization is not too long and acceptable. The downside is that the solved amplitudes of the displacement components are approxmations of the peaks of the envelopes of the time waveforms solved in the other two methods.

7. Conclusion

In this paper we introduced the Kriging surrogate model, used as a substitute for the black box function of an omnidirectional EMAT, in the optimization of this EMAT. The specific implementation of the Kriging surrogate model adopted is the DACE, suitable to model the deterministic output of a computer experiment.

Firstly the formulations of the DACE were reviewed, in which the output is modeled as sum of a regression part and a deviation part. Both the deviation and the total output are treated as random processes. The essential idea behind DACE is that the real response is one realisation of the output random process, and correlation between two samples of the deviation random process decreases with increasing distance between the two samples.

Then the formulations of an omnidirectional EMAT working on the Lorentz force mechanism were summarized. The underlying integro-differential equations describing the eddy current phenomenon were solved with a mid-level method proposed previously by the authors, with the benefit that the physics could be controlled completely. After this, the specifications of the model of the omnidirectional EMAT for generation of Lamb waves on a linearised steel plate were introduced. The axisymmetric model was divided into 2 geometries and 3 sub-models. Radial and axial displacement components, corresponding to the S0 mode and the A0 mode Lamb waves respectively, of an observation point at the middle plane of the plate were recorded to solve the objective function, defined as the ratio of the amplitude of the S0 mode to that of the A0 mode.

Next three different approaches to calculate the amplitudes of the displacement components were discussed. The first one is solving the peaks of the envelopes of the time waveforms from time-domain simulation, which involves small time steps and is very time-consuming. The second one is transforming the input excitation signal into its spectrum and feeding the spectrum to the frequency model, then transforming the modified spectrum back to the time domain and solving the peaks of the envelopes of the time waveforms. The third one is using the magnitudes of the phasors of the displacement components directly, i.e. using a single frequency in the frequency model. The waveforms from the first and second approaches were compared. The reason why the third approach could be applied as an approximation when evaluating the amplitudes was also explained.

Then different numbers ( $m$ ) of the known input/output data pairs chosen with the Latin Hypercube Sampling rule were used to construct the DACE model. The DACE model was then used to predict the output of the EMAT model on a 60 $\times$ 60 grid of the design variables. The single frequency model (Approach 3) was used first, with true outputs drawn on a contour plot and a surface plot as the reference. Then the multi-frequency model (Approach 2) was used. When $m=$ 800, by observation, the prediction was accurate enough. The MSE surface was also solved for the case of the multi-frequency model.

Finally the obtained DACE model of the EMAT with the multifrequency approach was applied for optimization of the EMAT. As already discussed, the objective function to minimize is the ratio of the amplitude of the S0 mode to that of the A0 mode, i.e. preferably generating the A0 mode Lamb waves. The chosen global optimization algorithm was the simulated annealing. The program was run 5 times and always gave the same results. The design variables at the optimal point were compared with previous work to validate the Kriging surrogate model (DACE). With the surrogate model, the total time of the optimization was decreased greatly, while at the same time, the accuracy of optimization was acceptable.

Footnotes

This crucial step of force transfer is implemented with Extrusion coupling operator in COMSOL, which could be used to map an expression from a source to a destination.

Obviously the material is ferromagnetic, nevertheless, the magnetostriction effect is not considered here, as in [].

$\lambda$ was also used to represent one of the Lamé constants.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (grant No. 51677093 and 51777100), National Key Scientific Instrument and Equipment Development Project (grant No. 2013YQ140505), and China Scholarship Council (grant No. 201506215055).

References

Auld

B.A.

, Acoustic Fields and Waves in Solids. Krieger Publishing Company, 1990.

Dhayalan

and Balasubramaniam

, A hybrid finite element model for simulation of electromagnetic acoustic transducer (EMAT) based plate waves, NDT&E International 43(6) (SEP 2010), 519–526.

Dhayalan

Kumar

Rao

B.P.C.

and Jayakumar

, A hybrid finite element model for spiral coil electromagnetic acoustic transducer (EMAT), International Journal of Applied Electromagnetics and Mechanics 46(3) (2014), 491–500.

Fromme

, Noncontact Excitation of Guided Waves (A(0) Mode) Using an Electromagnetic Acoustic Transducer (EMAT), In Chimenti

D.E.

and Bond

L.J.

, editor, 42Nd Annual Review Of Progress In Quantitative Nondestructive Evaluation: Incorporating The 6th European-American Workshop On Reliability Of Nde, volume 1706 of AIP Conference Proceedings. Ctr Nondestruct Evaluat; QNDE Programs, 2016. 42nd Annual Review of Progress in Quantitative Nondestructive Evaluation (QNDE), Minneapolis, MN, JUL 26–31, 2015.

Garcia-Rodriguez

Mihalache

and Ueda

, EMAT simulations based on a two-dimensional FEM coupled electro-mechanical formulation, International Journal of Applied Electromagnetics and Mechanics 45(1-4) (2014), 543–549. 16th International Symposium on Applied Electromagnetics and Mechanics (ISEM), Univ Laval, Pavill Alphonse Desjardins, Quebec City, CANADA, Jul 31–Aug 02, 2013.

Ida

, Engineering Electromagnetics, Springer, 2007.

Jafari-Shapoorabadi

Konrad

and Sinclair

, Improved finite element method for EMAT analysis and design. IEEE Transactions On Magnetics 37(4, 1) (Jul 2001) 2821–2823. 8th Joint Magnetism and Magnetic Materials International Magnetics Conference (MMM-INTERMAG), SAN ANTONIO, TEXAS, Jan 7–11, 2001.

Kay

S.M.

, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory. Prentice Hall, 1993.

Kay

S.M.

, Intuitive Probability and Random Processes using MATLAB. Springer, 2005.

10.

Konrad

, The numerical-solution of steady-state skin effect problems - an integrodifferential approach, IEEE Transactions On Magnetics 17(1) (1981), 1148–1152.

11.

Ludwig

and Dai

, Numerical-simulation of electromagnetic acoustic transducer in the time domain, Journal of applied physics 69(1) (1 Jan 1991), 89–98.

12.

Mirkhani

Chaggares

Masterson

Jastrzebski

Dusatko

Sinclair

Shapoorabadi

Konrad

and Papini

, Optimal design of EMAT transmitters, NDT&E International 37(3) (Apr 2004), 181–193.

13.

Preis

, A contribution to eddy-current calculations in plane and axisymmetric multiconductor systems, IEEE Transactions on Magnetics 19(6) (1983), 2397–2400.

14.

Lophaven

S.N.

Nielsen

H.B.

and Søndergaard

, DACE A MATLAB Kriging Toolbox, Technical report, Technical University of Denmark, 2002.

15.

Sacks

Welch

W.J.

Mitchell

T.J.

and Wynn

H.P.

, Design and analysis of computer experiments, Statist. Sci. 4(4) (1989), 409–423.

16.

Santner

T.J.

Williams

B.J.

and Notz

W.I.

, The Design and Analysis of Computer Experiments, Springer, 2003.

17.

Seher

Huthwaite

Lowe

Nagy

and Cawley

, Numerical design optimization of an EMAT for A0 Lamb wave generation in steel plates, in: Chimenti

Bond

and Thompson

, eds, 40TH Annual Review Of Progress In Quantitative Nondestructive Evaluation: Incorporating The 10Th International Conference On Barkhausen Noise And Micromagnetic Testing, Vols 33A& 3B, volume 1581 of AIP Conference Proceedings, pages 340–347. Quantitat Nondestruct Evaluat Programs; Amer Soc Nondestruct Testing; NDE Centers, World Federat; Stresstech Grp & Int Conf Barkhausen Noise & Micromagnet Testing Org Comm; Int Comm NDT; Natl Sci Fdn Ind Univ Cooperat Res Ctr; Iowa State Univ, Ctr Nondestruct Evaluat, 2014.

18.

Seher

Huthwaite

Lowe

M.J.S.

and Nagy

P.B.

, Model-based design of low frequency Lamb wave EMATs for mode selectivity, Journal Of Nondestructive Evaluation 34(3) (Sep 2015).

19.

Shi

Gong

Zhao

Fan

and Tan

, Comparison of three formulations for eddy-current and skin and proximity problems in a meander coil electromagnetic acoustic transducer, International Journal of Applied Electromagnetics and Mechanics 49(2) (2015), 239–249.

20.

Simpson

T.W.

Mauery

T.M.

Korte

J.J.

and Mistree

, Kriging models for global approximation in simulation-based multidisciplinary design optimization, AIAA Journal 39(12) (Dec 2001), 2233–2241.

21.

Skiena

S.S.

, The Algorithm Design Manual. Springer, 2nd edition, 2008.

22.

Sykulski

, New trends in optimization in electromagnetics, In Institution of Engineering and Technology, 7th International Conference on Computation in Electromagnetics, CEM 2008, pages 44–9, 2008.

23.

Wang

Huang

Velichko

Wilcox

and Zhao

, A multi-objective structural optimization of an omnidirectional electromagnetic acoustic transducer, Ultrasonics 81 (2017), 23–31.

24.

Wang

Huang

Zhang

and Zhao

, Modeling of an omni-directional electromagnetic acoustic transducer driven by the Lorentz force mechanism, Smart Materials and Structures 25(12) (2016), 25029.

25.

Wang

Huang

Zhang

and Zhao

, Multiphysics modeling of a Lorentz force-based meander coil electromagnetic acoustic transducer via steady-state and transient analyses, IEEE Sensors Journal 16(17) (1 Sep 2016), 6641–6651.

26.

Wilcox

Lowe

and Cawley

, The excitation and detection of Lamb waves with planar coil electromagnetic acoustic transducers, IEEE Transactions on Ultrasonics Ferroelectrics And Frequency Control 52(12) (Dec 2005), 2370–2383.