# Simulation of an Electromagnetic Acoustic Transducer Array by Using Analytical Method and FDTD

Authors: Yuedong Xie; Sergio Rodriguez; Wenbo  Zhang; Zenghua Liu; Wuliang Yin
Journal: Journal of Sensors
 (2016)

---

## Abstract

Previously, we developed a method based on FEM and FDTD for the study of an Electromagnetic Acoustic Transducer Array (EMAT). This paper presents a new analytical solution to the eddy current problem for the meander coil used in an EMAT, which is adapted from the classic Deeds and Dodd solution originally intended for circular coils. The analytical solution resulting from this novel adaptation exploits the large radius extrapolation and shows several advantages over the finite element method (FEM), especially in the higher frequency regime. The calculated Lorentz force density from the analytical EM solver is then coupled to the ultrasonic simulations, which exploit the finite-difference time-domain (FDTD) method to describe the propagation of ultrasound waves, in particular for Rayleigh waves. Radiation pattern obtained with Hilbert transform on time-domain waveforms is proposed to characterise the sensor in terms of its beam directivity and field distribution along the steering angle, which can produce performance parameters for an EMAT array, facilitating the optimum design of such sensors.


---

## Body

## 1. Introduction

There are a variety of nondestructive testing (NDT) techniques employed in industries, such as magnetic particle inspection (MPI), electromagnetic methods (EM), eddy current methods, and ultrasonic methods [1–7]. Due to its advantages of having good penetration depth and mechanical flexibility, the piezoelectric ultrasonic method is widely used for thickness measurement, flaw evaluation, and material characterization [8–13]. The transducer frequently used is made of piezoelectric ceramics or crystals [12–15]. However, one primary disadvantage of the piezoelectric ultrasonic testing is the need to have good sonic contact with the test piece, typically by means of a couplant for acoustic impedance matching [16].Electromagnetic acoustic transducers (EMATs) are becoming increasingly popular due to their noncontact nature [17, 18]. An EMAT sensor typically consists of a permanent magnet providing a large static magnetic field and a coil carrying an alternating current which is placed next to the test piece [19–21]. There are two EMAT interactions which can produce ultrasound: magnetostriction for magnetic materials and the Lorentz force mechanism for conducting metallic materials [18, 20, 22]. Because an EMAT generates ultrasonic waves directly into the testing piece instead of coupling through the transducer, an EMAT has advantages in applications where surface contact is not possible or desirable [23, 24]. Another attractive feature of EMAT which is a variety of waves modes can be produced based on different combinations of coils and magnets [18, 25]. In this paper, only EMAT based on Lorentz force mechanism to generate surface waves is discussed.Considerable works have been reported on EMAT modelling [26–30]. All of these papers divide EMAT modelling into two parts, electromagnetic simulation to obtain Lorentz force density and ultrasonic simulation to model the ultrasonic wave propagation due to Lorentz force. Electromagnetic simulation can be achieved by the finite element method (FEM) and the analytical method; ultrasonic simulation can be performed with the finite element method (FEM), finite-difference time-domain (FDTD), and the analytical method. Some papers combine the finite element method (FEM) and the analytical method to model EMAT; that is, the finite element method (FEM) is used to carry out electromagnetic simulation and the analytical method is to achieve ultrasonic simulation [30–33]. Some papers model EMAT arrays with the finite element method (FEM) for both electromagnetic and ultrasonic simulations, that is, the implicit finite element software COMSOL for the electromagnetic simulation and the explicit finite element software Abaqus for the ultrasonic simulation [18, 22]. The summary of the state of the art methods used for EMAT modelling is shown in Table 1; the method combining finite element method (FEM) and finite-difference time-domain (FDTD), the method wholly using analytical solutions, and the method combining the analytical method and the finite-difference time-domain (FDTD) method to model EMAT arrays have not been studied before.Table 1
Summary of methods used for modelling EMAT.
People
Electromagnetic simulation
Ultrasonic simulation
FEM
Analytical method
FEM
FDTD
Analytical method
[30–33]
[18, 22]
Authors
Authors
AuthorsThe method using finite element method (FEM) and finite-difference time-domain (FDTD) to model EMATs has been reported recently by authors in [34]. This paper proposes a new method using an analytical method and the finite-difference time-domain (FDTD) method to model EMAT arrays. On one hand, the EM analytical approach is used to calculate Lorentz force density for a given coil and a particular DC biased magnet configuration, which then can be fed through to ultrasonic simulations. On the other hand, the finite-difference time-domain (FDTD) method is exploited to describe the ultrasonic wave propagation due to the EM Lorentz force density acting upon the metallic sample.


## 2. EMAT Modelling

An EMAT sensor consists basically of a coil carrying an alternating current, a permanent magnet providing a large static magnetic field, and the test piece, as shown in Figure1. The coil induces eddy currents J in the surface layers of the testing material, and the interaction between the static magnetic field B and eddy currents J produces a Lorentz force density F based on (1), which in turn generates ultrasound waves propagating within the testing sample:(1)F=J×B.In this work, the EMAT modelling consists of two simulations, electromagnetic simulation and ultrasonic simulation. Electromagnetic simulation is used to obtain the eddy current and further the Lorentz force density, which is the link between electromagnetic simulation and ultrasonic simulation.Figure 1
The configuration of a typical EMAT.The testing sample used is a stainless steel plate with a dimension of 1000 × 1000 × 80 mm3, and the permanent magnet used is NdFeB35, whose size is 80 × 80 × 30 mm3. The meander coil carries an alternating current with the peak of 50 A, the lift-off is 1 mm, the operation frequency is 500 kHz, and the skin depth calculated is 0.679 mm. The Rayleigh wave velocity is 3.033 mm/μs in the stainless steel plate used, so the centre-to-centre distance between two adjacent lines of the meander coil is 3.033 mm to form the constructive interference.


### 2.1. Electromagnetic Modelling

For electromagnetic calculation, an analytical solution is adapted from the Deeds and Dodd formula to obtain the magnetic vector potential and the eddy current density. Other analytical solutions are available as well [36, 37]; the analytical solution proposed by Deeds and Dodd is employed as it is of highly satisfactory accuracy and minimum model difference between [35] and authors built.


#### 2.1.1. The Governing Equations

Dodd and Deeds proposed analytical solutions to the circular coil over a layered conductor in [35]. In this work, the geometry used for EM calculation is shown in Figure 2, where a circular coil is placed above the test piece made of stainless steel plate. r1 and r2 are the inner and outer radius of the circular coil, h is the height of the coil, l is the lift-off distance, and R and d are the length and height of stainless steel plate, respectively.Figure 2
The geometry used in [35] for eddy current calculation.The governing equations for induced eddy current calculation are(2)∇2A=-μI+μσ∂A∂t+μϵ∂2A∂t2+μ∇1μ×∇×A,E=-jωA,J=σE,where A is the vector potential, μ, σ, and ϵ are the permeability, conductivity, and permittivity of the material, respectively, I is the applied current density, ω is the angular frequency of the applied alternating current, E is the electric field, and J is the induced eddy current.From (2), for eddy current calculation, the main thing is to calculate the vector potential. Dodd and Deeds provided the final analytical solutions to the vector potential calculation as follows [35]:(3)Ar,z=M∫0∞1a2∫ar1ar2xJ1xdxJ1are-al-e-al+hea1za+a1da,M=μ0NI,a1=a2+jωμ1σ,where N is the turns of the coil, a and x are the integration variables, and J1ar and J1x are the Bessel functions of first kind. μ0 and μ1 are the permeability of air and metal, respectively.We build a model to study the analytical solutions to the vector potential problem; the test piece used is stainless steel, and the parameters used are listed in Table2. By calculating with Matlab, the magnitude distribution of the vector potential A along the surface of stainless steel (z=0) is shown in Figure 3, where the vector potential is mainly concentrated under the circular coil; the unit of the vector potential is tesla ∗ meter. The blue square in Figure 3 means the maximum vector potential, which shows the distribution of the vector potential is not symmetrical with r=2.5mm because the wire of the circular coil is not straight.Table 2
Parameters used for studying the analytical solutions.
Description
Symbol
Value
The length of  the stainless steel
R
5 mm
The height of  the stainless steel
d
5 mm
Inside radius of the circular coil
r
1
2.45 mm
Mean radius of the circular coil
r
1
+
r
2
2
2.5 mm
The height of the coil
h
1 mm
Permeability of air
μ
0
1.2566 × 10−6 H/m
Lift-off
l
1 mm
Current density
I
1 A/m2
Outside radius of the circular coil
r
2
2.55 mm
Frequency
f
10 kHz
Conductivity of stainless steel
σ
1100000 siemens/m
Permeability of stainless steel
μ
1
1.26 × 10−6 H/mFigure 3
The magnitude distribution of the vector potentialA under a circular coil.


#### 2.1.2. Analytical Method for Meander Coil

In this work, the coil used in EMAT is a meander coil, so the analytical solutions to a straight wire are needed. Based on the analytical solutions proposed by Dodd and Deeds, we proposed an assumption; that is, when the radius of the circular coil is very large, the bent wire of the circular coil can be approximated to a straight wire, and the distribution of the vector potential would be symmetrical. To verify this assumption, we build a model with the same parameters used in Table2, except that the mean radius of the circular coil is 20.05 m, and the length of the stainless steel sample is 20.1 m. The magnitude distribution of the vector potential A along the surface (z=0) is shown in Figure 4. From this figure, the magnitude distribution is symmetrical with r=20.05m, where the wire of the coil is located. This verifies the assumption; that is, when the radius of the spiral coil is very large, the bend wire can be approximated to a straight wire and the solution of Deeds and Dodd can be adapted for this application.Figure 4
The magnitude distribution of the vector potentialA under a straight wire.


#### 2.1.3. Comparison with FEM at 10 kHz

For a straight wire solution, in order to investigate the accuracy of the adapted analytical solutions, the comparison between the analytical solution and the finite element method (FEM) is needed. Maxwell Ansoft is used to construct a model with the same parameters used in Section2.1.2. The vacuum region to be solved is four times as large as the stainless steel plate; the total number of mesh elements is 43985. The computation of the FEM solver is based on minimising the energy error, which is as low as 0.068% in the final iteration (the fifth). The vector distribution on the y-z section of the stainless steel plate with the analytical method and the numerical method is shown in Figures 5(a) and 5(b), respectively. In Figure 5, both the analytical method and the finite element method (FEM) show that the maximum vector potential occurs at the surface of stainless steel plate and decreases along the depth. The magnitude range of the analytical method is consistent with that of the finite element method (FEM), which is from 0 to 3.5×10-7 tesla ∗ meter.Figure 5
The magnitude distribution of the vector potential based on different methods.
(a)
Analytical method
(b)
Finite element method (FEM)On the surface of stainless steel (z=0), the distribution of the vector potential is shown in Figure 6, where the magnitude, the real part, and the imaginary part of the vector potential distribution are shown in Figures 6(a), 6(b), and 6(c), respectively. In Figures 6(a) and 6(b), for magnitude and real part distribution, the analytical method and the finite element method are consistent. However, for the imaginary part of the vector potential, as shown in Figure 6(c), the analytical method shows a more accurate result than FEM because the vector potential based on FEM is not approaching zero when A is away from the wire.Figure 6
With an operating frequency of 10 kHz, the distribution of the vector potential along the surface of stainless steel; (a) the magnitude distribution, (b) the real part distribution, and (c) the imaginary distribution.
(a)
(b)
(c)


#### 2.1.4. Comparison with FEM at 1 MHz

In this part, the analytical solutions at a high operating frequency are studied. This is because EMAT normally operates at high frequencies and the eddy current is typically limited near the surface. The model used is the same as that in Section2.1.3, except that the frequency used is 1 MHz.Along the surface of the stainless steel (z=0), the real part distribution of the vector potential based on the analytical method and the finite element method (FEM) is shown in Figure 7. The curve obtained by the analytical method is smoother than that of FEM; that is because of the numerical nature of the FEM; numerical approximation due to finite mesh density and element interpolation are inevitable. With comparison vector potential distribution at 10 kHz and 1 MHz, as shown in Figures 6 and 7, respectively, the analytical method is of satisfactory accuracy to describe the electromagnetic properties and is more precise than FEM, in particular for the higher frequency regime.Figure 7
With an operating frequency of 1 MHz, the real part distribution of the vector potential along the stainless steel surface. Real part ofA under a straight wire.


#### 2.1.5. Analytical EMAT-EM Modelling

The meander coil used in this study has a dimension of 56 × 34.163 × 0.036 mm3, which is very small compared to the size of the stainless steel plate. In order to improve modelling time, only the area (100 × 100 × 2 mm3) where the meander coil mainly has an effect is picked to study the Lorentz force distribution.As mentioned before, the distribution of the induced eddy current under a straight wire can be obtained by the analytical solutions. For a meander coil, the total induced eddy current is the sum of the induced eddy current caused by each wire segment; the distribution of the induced eddy current ony-z section is shown in Figure 8, where fields between two adjacent wires are opposite due to opposite directions of the alternating currents. In addition, the values of the eddy current under the outmost wires are the largest, because the outmost wires are only affected by the fields on one side.Figure 8
The distribution of the induced eddy current based on the analytical method. Real part of the induced eddy current.Along the surface of the stainless steel plate, eddy current distribution is shown in Figure9(a), which confirms the observation that the amplitude along the outmost wire is the largest. Because there are six pairs of adjacent wires with different current directions, the eddy current plot has 6 crests and 6 troughs. Figure 9(b) shows the distribution of Lorentz force density on the surface of the stainless steel plate. It can be seen that the Lorentz force density on the outmost lines is larger than that on the inner lines; that is because both the maximum magnetic field and the maximum eddy current occur at the places corresponding to the edges of the EMAT sensor.Figure 9
Fields distribution along the surface of the stainless steel plate; (a) the distribution of the induced eddy current and (b) the distribution of Lorentz force density.
(a)
(b)


### 2.2. Ultrasonic Modelling

#### 2.2.1. Governing Equations

Elastodynamic equations are a set of partial differential equations describing how material deforms and becomes internally stressed as shown in the following [38, 39]:(4)ρx∂vi∂tx,t=∑j=1d∂Tij∂xjx,t+fix,t,(5)∂Ti∂t=∑j=1d∑i=1dcijklx∂vk∂xl+θijx,t,where ρ is the mass density and cijkl is the 4th stiffness tensor of the testing sample; fi and θij are the force source and strain tensor rate source, respectively. The parameters to be calculated are velocity vi and stress tensor Ti. Equation (4) is Newton’s Second Law: when a force is applied to a testing sample, stress and deformation are generated, as well as particle displacement. Equation (5) is, based on Hooke’s Law, describing the relationship of stress tensor rate and strain tensor rate when deformation occurs.The finite-difference time-domain (FDTD) method is a numerical method to solve differential equations by discretion of the differential form to the finite-difference form [40]. In this paper, forward difference and centre difference methods are used to calculate the unknown parameters velocity vi and stress tensor Ti [39].


#### 2.2.2. Combination of Electromagnetic Simulation and Ultrasonic Simulation

In this work, Lorentz force density obtained from the electromagnetic model is imported to ultrasonic model to generate ultrasound waves. As shown in Figure10, the 12 alternating Lorentz force densities are added to the ultrasonic model along the surface. The original y-z cross-section is 1000∗80 mm2; in order to save modelling geometry and improve modelling time, only the area (400∗80 mm2) where the EMAT sensor mainly has an effect and surface waves are mainly propagating along is chosen for modelling. Two receivers, R1 and R2, are placed on (50,79) and (50,77), respectively, to inspect the arrival signals.Figure 10
Transformation from electromagnetic model to ultrasonic model.The propagation of ultrasound waves is shown in Figure11: at 23 μs after firing, the bulk waves, head waves, and surface waves can be identified. Bulk waves contain longitudinal waves and shear waves, which are obliquely propagating into the material; the velocity of the longitudinal waves is larger than that of the shear wave, so the longitudinal wave arrives earlier than the shear wave. Surface waves, which are Rayleigh waves in this work, are propagating along the surface and the subsurface of the material. The velocities of Rayleigh waves and shear waves are slightly different; in most of situations, the velocity of Rayleigh waves is 90 percent of that of the shear waves. So the propagation of Rayleigh waves is slightly delayed than that of shear waves, as shown in Figure 12, where the Rayleigh waves can be identified more clearly at 45 μs.Figure 11
Wave propagation at 23μs. Four types of waves are produced: longitudinal waves, shear waves, head waves, and Rayleigh waves.Figure 12
Wave propagation at 45μs. Rayleigh waves’ propagation is slightly delayed compared to the shear waves’ propagation.The receiving signals from R1 and R2 are shown in Figure13; the first arrival signal is the longitudinal wave arriving at about 25.5 μs, and the latter arrival signal is the Rayleigh wave at about 49.5 μs with a large amplitude compared to that of longitudinal waves. The Rayleigh waves can be validated by the distance-of-flight calculation: the distance from R1 and the EMAT sensor is a constant, so the flight distance of the longitudinal waves and the Rayleigh waves should be the same. The velocities of the longitudinal waves and the Rayleigh waves used are 5.9 mm/μs and 3.033 mm/μs, respectively, and the arrival times of these two waves are 25.5 μs and 49.5 μs, respectively, so the flight distance of the longitudinal waves and the Rayleigh waves calculated is 150.45 mm and 150.13 mm, which are almost the same. In addition, the amplitude of R1 is slightly larger than that of R2, which confirms that the energy of the Rayleigh waves is mainly distributed in the surface of the material.Figure 13
Receiving signals from receivers R1 and R2.The receiving coil used is the same as the transmitting coil; the induced voltage in the receiving coil can be calculated from the receiving velocity fields [41]. In this work, receiving signals calculated are in millivolts. In practical applications, a high power amplifier is essential for EMAT transduction; RITEC RPR4000 and RITEC RAM5000, which are typically used in an EMAT system, are capable of amplifying the receiving signals to volts.


### 2.3. The Radiation Pattern of EMAT-Rayleigh Waves

Most of the previous works calculate the radiation pattern based on analytical equations [42–44]. In this paper, we define the radiation pattern by combining the FDTD method and the Hilbert transform. For a geometry containing M∗N points, each point on the geometry has a time series signal. The max amplitude of the time series signal indicates the arrival time of the signal, and the Hilbert transform is used to calculate the envelope of the signal (6), resulting in identification of the signal arrival time more clearly, as shown in Figure 14:(6)ht=ft∗1πt,zt=ft+jht,et=ft2+ht2,where ft is the time series signal, ht is the signal after Hilbert transform, zt is the analytical signal, and e(t) is the envelope of zt.Figure 14
The amplitude and envelope of receiving signals from R1.The radiation pattern of the EMAT-Rayleigh waves is shown in Figure15(a), which shows that Rayleigh waves are mainly distributed along the surface and subsurface of the material; longitudinal and shear waves are propagating obliquely into the material, and shear waves have a stronger energy than longitudinal waves. In this work, only Rayleigh waves are of interest; the beam features of Rayleigh waves are studied by means of Figure 15(b).Figure 15
The radiation pattern of EMAT-Rayleigh waves for studying beam features.
(a)
The radiation pattern of EMAT-Rayleigh waves
(b)
Radiation pattern for beam featuresThere are two beam features to be analysed, beam directivity and field distribution along the steering angle. Beam directivity is, at a radial length from the centre of the array, the velocity or pressure distribution, as shown in the green circle arc in Figure15(b). Field distribution along the steering angle is the velocity or pressure distribution along the steering angle, as shown in the red line in Figure 15(b). r denotes the radial distance from the middle of the sensor and θ denotes the steering angle; define the steering angle parallel to the surface as 0°; in order to minimise the effects of the longitudinal waves and shear waves, keep the radial distance r 140 mm and the steering angle θ from 0° to 20°. The beam directivity calculated is shown in Figure 16, where the Rayleigh waves are mainly distributed along the steering angles from 0° to 2.5°. Because the radial distance r used is 140 mm, the depth d of the Rayleigh waves’ energy distribution can be calculated by(7)d=tan⁡θ∗r.The depth d calculated is 6.1 mm, and the wavelength of the Rayleigh waves used in this work is 6.066 mm, which validates that the energy of the Rayleigh waves is mainly distributed within one wavelength of the Rayleigh waves. In addition, when choosing another value of the radial distance r, the steering angle θ changes correspondingly, but the depth d of the Rayleigh waves’ energy distribution is a constant.Figure 16
The beam directivity and field distribution along the steering angle 0° of EMAT-Rayleigh waves.Another beam feature, field distribution along the steering angle, is studied as well. The radial distancer used is 200 mm from the centre of EMAT arrays, and the steering angle used is 0° (along the surface). The field distribution along the steering angle 0° is shown in Figure 16 and the length period of the maximum magnitudes in this figure is about 17.5 mm, which is consistent with the modelling geometry, where the distance from the centre to the end of the EMAT sensor is 17.1 mm. As a result, the field energy is maximum at the places where the sensor is placed and decreases along the radial distance due to the attenuation of Rayleigh waves.


## 2.1. Electromagnetic Modelling

For electromagnetic calculation, an analytical solution is adapted from the Deeds and Dodd formula to obtain the magnetic vector potential and the eddy current density. Other analytical solutions are available as well [36, 37]; the analytical solution proposed by Deeds and Dodd is employed as it is of highly satisfactory accuracy and minimum model difference between [35] and authors built.


### 2.1.1. The Governing Equations

Dodd and Deeds proposed analytical solutions to the circular coil over a layered conductor in [35]. In this work, the geometry used for EM calculation is shown in Figure 2, where a circular coil is placed above the test piece made of stainless steel plate. r1 and r2 are the inner and outer radius of the circular coil, h is the height of the coil, l is the lift-off distance, and R and d are the length and height of stainless steel plate, respectively.Figure 2
The geometry used in [35] for eddy current calculation.The governing equations for induced eddy current calculation are(2)∇2A=-μI+μσ∂A∂t+μϵ∂2A∂t2+μ∇1μ×∇×A,E=-jωA,J=σE,where A is the vector potential, μ, σ, and ϵ are the permeability, conductivity, and permittivity of the material, respectively, I is the applied current density, ω is the angular frequency of the applied alternating current, E is the electric field, and J is the induced eddy current.From (2), for eddy current calculation, the main thing is to calculate the vector potential. Dodd and Deeds provided the final analytical solutions to the vector potential calculation as follows [35]:(3)Ar,z=M∫0∞1a2∫ar1ar2xJ1xdxJ1are-al-e-al+hea1za+a1da,M=μ0NI,a1=a2+jωμ1σ,where N is the turns of the coil, a and x are the integration variables, and J1ar and J1x are the Bessel functions of first kind. μ0 and μ1 are the permeability of air and metal, respectively.We build a model to study the analytical solutions to the vector potential problem; the test piece used is stainless steel, and the parameters used are listed in Table2. By calculating with Matlab, the magnitude distribution of the vector potential A along the surface of stainless steel (z=0) is shown in Figure 3, where the vector potential is mainly concentrated under the circular coil; the unit of the vector potential is tesla ∗ meter. The blue square in Figure 3 means the maximum vector potential, which shows the distribution of the vector potential is not symmetrical with r=2.5mm because the wire of the circular coil is not straight.Table 2
Parameters used for studying the analytical solutions.
Description
Symbol
Value
The length of  the stainless steel
R
5 mm
The height of  the stainless steel
d
5 mm
Inside radius of the circular coil
r
1
2.45 mm
Mean radius of the circular coil
r
1
+
r
2
2
2.5 mm
The height of the coil
h
1 mm
Permeability of air
μ
0
1.2566 × 10−6 H/m
Lift-off
l
1 mm
Current density
I
1 A/m2
Outside radius of the circular coil
r
2
2.55 mm
Frequency
f
10 kHz
Conductivity of stainless steel
σ
1100000 siemens/m
Permeability of stainless steel
μ
1
1.26 × 10−6 H/mFigure 3
The magnitude distribution of the vector potentialA under a circular coil.


### 2.1.2. Analytical Method for Meander Coil

In this work, the coil used in EMAT is a meander coil, so the analytical solutions to a straight wire are needed. Based on the analytical solutions proposed by Dodd and Deeds, we proposed an assumption; that is, when the radius of the circular coil is very large, the bent wire of the circular coil can be approximated to a straight wire, and the distribution of the vector potential would be symmetrical. To verify this assumption, we build a model with the same parameters used in Table2, except that the mean radius of the circular coil is 20.05 m, and the length of the stainless steel sample is 20.1 m. The magnitude distribution of the vector potential A along the surface (z=0) is shown in Figure 4. From this figure, the magnitude distribution is symmetrical with r=20.05m, where the wire of the coil is located. This verifies the assumption; that is, when the radius of the spiral coil is very large, the bend wire can be approximated to a straight wire and the solution of Deeds and Dodd can be adapted for this application.Figure 4
The magnitude distribution of the vector potentialA under a straight wire.


### 2.1.3. Comparison with FEM at 10 kHz

For a straight wire solution, in order to investigate the accuracy of the adapted analytical solutions, the comparison between the analytical solution and the finite element method (FEM) is needed. Maxwell Ansoft is used to construct a model with the same parameters used in Section2.1.2. The vacuum region to be solved is four times as large as the stainless steel plate; the total number of mesh elements is 43985. The computation of the FEM solver is based on minimising the energy error, which is as low as 0.068% in the final iteration (the fifth). The vector distribution on the y-z section of the stainless steel plate with the analytical method and the numerical method is shown in Figures 5(a) and 5(b), respectively. In Figure 5, both the analytical method and the finite element method (FEM) show that the maximum vector potential occurs at the surface of stainless steel plate and decreases along the depth. The magnitude range of the analytical method is consistent with that of the finite element method (FEM), which is from 0 to 3.5×10-7 tesla ∗ meter.Figure 5
The magnitude distribution of the vector potential based on different methods.
(a)
Analytical method
(b)
Finite element method (FEM)On the surface of stainless steel (z=0), the distribution of the vector potential is shown in Figure 6, where the magnitude, the real part, and the imaginary part of the vector potential distribution are shown in Figures 6(a), 6(b), and 6(c), respectively. In Figures 6(a) and 6(b), for magnitude and real part distribution, the analytical method and the finite element method are consistent. However, for the imaginary part of the vector potential, as shown in Figure 6(c), the analytical method shows a more accurate result than FEM because the vector potential based on FEM is not approaching zero when A is away from the wire.Figure 6
With an operating frequency of 10 kHz, the distribution of the vector potential along the surface of stainless steel; (a) the magnitude distribution, (b) the real part distribution, and (c) the imaginary distribution.
(a)
(b)
(c)


### 2.1.4. Comparison with FEM at 1 MHz

In this part, the analytical solutions at a high operating frequency are studied. This is because EMAT normally operates at high frequencies and the eddy current is typically limited near the surface. The model used is the same as that in Section2.1.3, except that the frequency used is 1 MHz.Along the surface of the stainless steel (z=0), the real part distribution of the vector potential based on the analytical method and the finite element method (FEM) is shown in Figure 7. The curve obtained by the analytical method is smoother than that of FEM; that is because of the numerical nature of the FEM; numerical approximation due to finite mesh density and element interpolation are inevitable. With comparison vector potential distribution at 10 kHz and 1 MHz, as shown in Figures 6 and 7, respectively, the analytical method is of satisfactory accuracy to describe the electromagnetic properties and is more precise than FEM, in particular for the higher frequency regime.Figure 7
With an operating frequency of 1 MHz, the real part distribution of the vector potential along the stainless steel surface. Real part ofA under a straight wire.


### 2.1.5. Analytical EMAT-EM Modelling

The meander coil used in this study has a dimension of 56 × 34.163 × 0.036 mm3, which is very small compared to the size of the stainless steel plate. In order to improve modelling time, only the area (100 × 100 × 2 mm3) where the meander coil mainly has an effect is picked to study the Lorentz force distribution.As mentioned before, the distribution of the induced eddy current under a straight wire can be obtained by the analytical solutions. For a meander coil, the total induced eddy current is the sum of the induced eddy current caused by each wire segment; the distribution of the induced eddy current ony-z section is shown in Figure 8, where fields between two adjacent wires are opposite due to opposite directions of the alternating currents. In addition, the values of the eddy current under the outmost wires are the largest, because the outmost wires are only affected by the fields on one side.Figure 8
The distribution of the induced eddy current based on the analytical method. Real part of the induced eddy current.Along the surface of the stainless steel plate, eddy current distribution is shown in Figure9(a), which confirms the observation that the amplitude along the outmost wire is the largest. Because there are six pairs of adjacent wires with different current directions, the eddy current plot has 6 crests and 6 troughs. Figure 9(b) shows the distribution of Lorentz force density on the surface of the stainless steel plate. It can be seen that the Lorentz force density on the outmost lines is larger than that on the inner lines; that is because both the maximum magnetic field and the maximum eddy current occur at the places corresponding to the edges of the EMAT sensor.Figure 9
Fields distribution along the surface of the stainless steel plate; (a) the distribution of the induced eddy current and (b) the distribution of Lorentz force density.
(a)
(b)


## 2.1.1. The Governing Equations

Dodd and Deeds proposed analytical solutions to the circular coil over a layered conductor in [35]. In this work, the geometry used for EM calculation is shown in Figure 2, where a circular coil is placed above the test piece made of stainless steel plate. r1 and r2 are the inner and outer radius of the circular coil, h is the height of the coil, l is the lift-off distance, and R and d are the length and height of stainless steel plate, respectively.Figure 2
The geometry used in [35] for eddy current calculation.The governing equations for induced eddy current calculation are(2)∇2A=-μI+μσ∂A∂t+μϵ∂2A∂t2+μ∇1μ×∇×A,E=-jωA,J=σE,where A is the vector potential, μ, σ, and ϵ are the permeability, conductivity, and permittivity of the material, respectively, I is the applied current density, ω is the angular frequency of the applied alternating current, E is the electric field, and J is the induced eddy current.From (2), for eddy current calculation, the main thing is to calculate the vector potential. Dodd and Deeds provided the final analytical solutions to the vector potential calculation as follows [35]:(3)Ar,z=M∫0∞1a2∫ar1ar2xJ1xdxJ1are-al-e-al+hea1za+a1da,M=μ0NI,a1=a2+jωμ1σ,where N is the turns of the coil, a and x are the integration variables, and J1ar and J1x are the Bessel functions of first kind. μ0 and μ1 are the permeability of air and metal, respectively.We build a model to study the analytical solutions to the vector potential problem; the test piece used is stainless steel, and the parameters used are listed in Table2. By calculating with Matlab, the magnitude distribution of the vector potential A along the surface of stainless steel (z=0) is shown in Figure 3, where the vector potential is mainly concentrated under the circular coil; the unit of the vector potential is tesla ∗ meter. The blue square in Figure 3 means the maximum vector potential, which shows the distribution of the vector potential is not symmetrical with r=2.5mm because the wire of the circular coil is not straight.Table 2
Parameters used for studying the analytical solutions.
Description
Symbol
Value
The length of  the stainless steel
R
5 mm
The height of  the stainless steel
d
5 mm
Inside radius of the circular coil
r
1
2.45 mm
Mean radius of the circular coil
r
1
+
r
2
2
2.5 mm
The height of the coil
h
1 mm
Permeability of air
μ
0
1.2566 × 10−6 H/m
Lift-off
l
1 mm
Current density
I
1 A/m2
Outside radius of the circular coil
r
2
2.55 mm
Frequency
f
10 kHz
Conductivity of stainless steel
σ
1100000 siemens/m
Permeability of stainless steel
μ
1
1.26 × 10−6 H/mFigure 3
The magnitude distribution of the vector potentialA under a circular coil.


## 2.1.2. Analytical Method for Meander Coil

In this work, the coil used in EMAT is a meander coil, so the analytical solutions to a straight wire are needed. Based on the analytical solutions proposed by Dodd and Deeds, we proposed an assumption; that is, when the radius of the circular coil is very large, the bent wire of the circular coil can be approximated to a straight wire, and the distribution of the vector potential would be symmetrical. To verify this assumption, we build a model with the same parameters used in Table2, except that the mean radius of the circular coil is 20.05 m, and the length of the stainless steel sample is 20.1 m. The magnitude distribution of the vector potential A along the surface (z=0) is shown in Figure 4. From this figure, the magnitude distribution is symmetrical with r=20.05m, where the wire of the coil is located. This verifies the assumption; that is, when the radius of the spiral coil is very large, the bend wire can be approximated to a straight wire and the solution of Deeds and Dodd can be adapted for this application.Figure 4
The magnitude distribution of the vector potentialA under a straight wire.


## 2.1.3. Comparison with FEM at 10 kHz

For a straight wire solution, in order to investigate the accuracy of the adapted analytical solutions, the comparison between the analytical solution and the finite element method (FEM) is needed. Maxwell Ansoft is used to construct a model with the same parameters used in Section2.1.2. The vacuum region to be solved is four times as large as the stainless steel plate; the total number of mesh elements is 43985. The computation of the FEM solver is based on minimising the energy error, which is as low as 0.068% in the final iteration (the fifth). The vector distribution on the y-z section of the stainless steel plate with the analytical method and the numerical method is shown in Figures 5(a) and 5(b), respectively. In Figure 5, both the analytical method and the finite element method (FEM) show that the maximum vector potential occurs at the surface of stainless steel plate and decreases along the depth. The magnitude range of the analytical method is consistent with that of the finite element method (FEM), which is from 0 to 3.5×10-7 tesla ∗ meter.Figure 5
The magnitude distribution of the vector potential based on different methods.
(a)
Analytical method
(b)
Finite element method (FEM)On the surface of stainless steel (z=0), the distribution of the vector potential is shown in Figure 6, where the magnitude, the real part, and the imaginary part of the vector potential distribution are shown in Figures 6(a), 6(b), and 6(c), respectively. In Figures 6(a) and 6(b), for magnitude and real part distribution, the analytical method and the finite element method are consistent. However, for the imaginary part of the vector potential, as shown in Figure 6(c), the analytical method shows a more accurate result than FEM because the vector potential based on FEM is not approaching zero when A is away from the wire.Figure 6
With an operating frequency of 10 kHz, the distribution of the vector potential along the surface of stainless steel; (a) the magnitude distribution, (b) the real part distribution, and (c) the imaginary distribution.
(a)
(b)
(c)


## 2.1.4. Comparison with FEM at 1 MHz

In this part, the analytical solutions at a high operating frequency are studied. This is because EMAT normally operates at high frequencies and the eddy current is typically limited near the surface. The model used is the same as that in Section2.1.3, except that the frequency used is 1 MHz.Along the surface of the stainless steel (z=0), the real part distribution of the vector potential based on the analytical method and the finite element method (FEM) is shown in Figure 7. The curve obtained by the analytical method is smoother than that of FEM; that is because of the numerical nature of the FEM; numerical approximation due to finite mesh density and element interpolation are inevitable. With comparison vector potential distribution at 10 kHz and 1 MHz, as shown in Figures 6 and 7, respectively, the analytical method is of satisfactory accuracy to describe the electromagnetic properties and is more precise than FEM, in particular for the higher frequency regime.Figure 7
With an operating frequency of 1 MHz, the real part distribution of the vector potential along the stainless steel surface. Real part ofA under a straight wire.


## 2.1.5. Analytical EMAT-EM Modelling

The meander coil used in this study has a dimension of 56 × 34.163 × 0.036 mm3, which is very small compared to the size of the stainless steel plate. In order to improve modelling time, only the area (100 × 100 × 2 mm3) where the meander coil mainly has an effect is picked to study the Lorentz force distribution.As mentioned before, the distribution of the induced eddy current under a straight wire can be obtained by the analytical solutions. For a meander coil, the total induced eddy current is the sum of the induced eddy current caused by each wire segment; the distribution of the induced eddy current ony-z section is shown in Figure 8, where fields between two adjacent wires are opposite due to opposite directions of the alternating currents. In addition, the values of the eddy current under the outmost wires are the largest, because the outmost wires are only affected by the fields on one side.Figure 8
The distribution of the induced eddy current based on the analytical method. Real part of the induced eddy current.Along the surface of the stainless steel plate, eddy current distribution is shown in Figure9(a), which confirms the observation that the amplitude along the outmost wire is the largest. Because there are six pairs of adjacent wires with different current directions, the eddy current plot has 6 crests and 6 troughs. Figure 9(b) shows the distribution of Lorentz force density on the surface of the stainless steel plate. It can be seen that the Lorentz force density on the outmost lines is larger than that on the inner lines; that is because both the maximum magnetic field and the maximum eddy current occur at the places corresponding to the edges of the EMAT sensor.Figure 9
Fields distribution along the surface of the stainless steel plate; (a) the distribution of the induced eddy current and (b) the distribution of Lorentz force density.
(a)
(b)


## 2.2. Ultrasonic Modelling

### 2.2.1. Governing Equations

Elastodynamic equations are a set of partial differential equations describing how material deforms and becomes internally stressed as shown in the following [38, 39]:(4)ρx∂vi∂tx,t=∑j=1d∂Tij∂xjx,t+fix,t,(5)∂Ti∂t=∑j=1d∑i=1dcijklx∂vk∂xl+θijx,t,where ρ is the mass density and cijkl is the 4th stiffness tensor of the testing sample; fi and θij are the force source and strain tensor rate source, respectively. The parameters to be calculated are velocity vi and stress tensor Ti. Equation (4) is Newton’s Second Law: when a force is applied to a testing sample, stress and deformation are generated, as well as particle displacement. Equation (5) is, based on Hooke’s Law, describing the relationship of stress tensor rate and strain tensor rate when deformation occurs.The finite-difference time-domain (FDTD) method is a numerical method to solve differential equations by discretion of the differential form to the finite-difference form [40]. In this paper, forward difference and centre difference methods are used to calculate the unknown parameters velocity vi and stress tensor Ti [39].


### 2.2.2. Combination of Electromagnetic Simulation and Ultrasonic Simulation

In this work, Lorentz force density obtained from the electromagnetic model is imported to ultrasonic model to generate ultrasound waves. As shown in Figure10, the 12 alternating Lorentz force densities are added to the ultrasonic model along the surface. The original y-z cross-section is 1000∗80 mm2; in order to save modelling geometry and improve modelling time, only the area (400∗80 mm2) where the EMAT sensor mainly has an effect and surface waves are mainly propagating along is chosen for modelling. Two receivers, R1 and R2, are placed on (50,79) and (50,77), respectively, to inspect the arrival signals.Figure 10
Transformation from electromagnetic model to ultrasonic model.The propagation of ultrasound waves is shown in Figure11: at 23 μs after firing, the bulk waves, head waves, and surface waves can be identified. Bulk waves contain longitudinal waves and shear waves, which are obliquely propagating into the material; the velocity of the longitudinal waves is larger than that of the shear wave, so the longitudinal wave arrives earlier than the shear wave. Surface waves, which are Rayleigh waves in this work, are propagating along the surface and the subsurface of the material. The velocities of Rayleigh waves and shear waves are slightly different; in most of situations, the velocity of Rayleigh waves is 90 percent of that of the shear waves. So the propagation of Rayleigh waves is slightly delayed than that of shear waves, as shown in Figure 12, where the Rayleigh waves can be identified more clearly at 45 μs.Figure 11
Wave propagation at 23μs. Four types of waves are produced: longitudinal waves, shear waves, head waves, and Rayleigh waves.Figure 12
Wave propagation at 45μs. Rayleigh waves’ propagation is slightly delayed compared to the shear waves’ propagation.The receiving signals from R1 and R2 are shown in Figure13; the first arrival signal is the longitudinal wave arriving at about 25.5 μs, and the latter arrival signal is the Rayleigh wave at about 49.5 μs with a large amplitude compared to that of longitudinal waves. The Rayleigh waves can be validated by the distance-of-flight calculation: the distance from R1 and the EMAT sensor is a constant, so the flight distance of the longitudinal waves and the Rayleigh waves should be the same. The velocities of the longitudinal waves and the Rayleigh waves used are 5.9 mm/μs and 3.033 mm/μs, respectively, and the arrival times of these two waves are 25.5 μs and 49.5 μs, respectively, so the flight distance of the longitudinal waves and the Rayleigh waves calculated is 150.45 mm and 150.13 mm, which are almost the same. In addition, the amplitude of R1 is slightly larger than that of R2, which confirms that the energy of the Rayleigh waves is mainly distributed in the surface of the material.Figure 13
Receiving signals from receivers R1 and R2.The receiving coil used is the same as the transmitting coil; the induced voltage in the receiving coil can be calculated from the receiving velocity fields [41]. In this work, receiving signals calculated are in millivolts. In practical applications, a high power amplifier is essential for EMAT transduction; RITEC RPR4000 and RITEC RAM5000, which are typically used in an EMAT system, are capable of amplifying the receiving signals to volts.


## 2.2.1. Governing Equations

Elastodynamic equations are a set of partial differential equations describing how material deforms and becomes internally stressed as shown in the following [38, 39]:(4)ρx∂vi∂tx,t=∑j=1d∂Tij∂xjx,t+fix,t,(5)∂Ti∂t=∑j=1d∑i=1dcijklx∂vk∂xl+θijx,t,where ρ is the mass density and cijkl is the 4th stiffness tensor of the testing sample; fi and θij are the force source and strain tensor rate source, respectively. The parameters to be calculated are velocity vi and stress tensor Ti. Equation (4) is Newton’s Second Law: when a force is applied to a testing sample, stress and deformation are generated, as well as particle displacement. Equation (5) is, based on Hooke’s Law, describing the relationship of stress tensor rate and strain tensor rate when deformation occurs.The finite-difference time-domain (FDTD) method is a numerical method to solve differential equations by discretion of the differential form to the finite-difference form [40]. In this paper, forward difference and centre difference methods are used to calculate the unknown parameters velocity vi and stress tensor Ti [39].


## 2.2.2. Combination of Electromagnetic Simulation and Ultrasonic Simulation

In this work, Lorentz force density obtained from the electromagnetic model is imported to ultrasonic model to generate ultrasound waves. As shown in Figure10, the 12 alternating Lorentz force densities are added to the ultrasonic model along the surface. The original y-z cross-section is 1000∗80 mm2; in order to save modelling geometry and improve modelling time, only the area (400∗80 mm2) where the EMAT sensor mainly has an effect and surface waves are mainly propagating along is chosen for modelling. Two receivers, R1 and R2, are placed on (50,79) and (50,77), respectively, to inspect the arrival signals.Figure 10
Transformation from electromagnetic model to ultrasonic model.The propagation of ultrasound waves is shown in Figure11: at 23 μs after firing, the bulk waves, head waves, and surface waves can be identified. Bulk waves contain longitudinal waves and shear waves, which are obliquely propagating into the material; the velocity of the longitudinal waves is larger than that of the shear wave, so the longitudinal wave arrives earlier than the shear wave. Surface waves, which are Rayleigh waves in this work, are propagating along the surface and the subsurface of the material. The velocities of Rayleigh waves and shear waves are slightly different; in most of situations, the velocity of Rayleigh waves is 90 percent of that of the shear waves. So the propagation of Rayleigh waves is slightly delayed than that of shear waves, as shown in Figure 12, where the Rayleigh waves can be identified more clearly at 45 μs.Figure 11
Wave propagation at 23μs. Four types of waves are produced: longitudinal waves, shear waves, head waves, and Rayleigh waves.Figure 12
Wave propagation at 45μs. Rayleigh waves’ propagation is slightly delayed compared to the shear waves’ propagation.The receiving signals from R1 and R2 are shown in Figure13; the first arrival signal is the longitudinal wave arriving at about 25.5 μs, and the latter arrival signal is the Rayleigh wave at about 49.5 μs with a large amplitude compared to that of longitudinal waves. The Rayleigh waves can be validated by the distance-of-flight calculation: the distance from R1 and the EMAT sensor is a constant, so the flight distance of the longitudinal waves and the Rayleigh waves should be the same. The velocities of the longitudinal waves and the Rayleigh waves used are 5.9 mm/μs and 3.033 mm/μs, respectively, and the arrival times of these two waves are 25.5 μs and 49.5 μs, respectively, so the flight distance of the longitudinal waves and the Rayleigh waves calculated is 150.45 mm and 150.13 mm, which are almost the same. In addition, the amplitude of R1 is slightly larger than that of R2, which confirms that the energy of the Rayleigh waves is mainly distributed in the surface of the material.Figure 13
Receiving signals from receivers R1 and R2.The receiving coil used is the same as the transmitting coil; the induced voltage in the receiving coil can be calculated from the receiving velocity fields [41]. In this work, receiving signals calculated are in millivolts. In practical applications, a high power amplifier is essential for EMAT transduction; RITEC RPR4000 and RITEC RAM5000, which are typically used in an EMAT system, are capable of amplifying the receiving signals to volts.


## 2.3. The Radiation Pattern of EMAT-Rayleigh Waves

Most of the previous works calculate the radiation pattern based on analytical equations [42–44]. In this paper, we define the radiation pattern by combining the FDTD method and the Hilbert transform. For a geometry containing M∗N points, each point on the geometry has a time series signal. The max amplitude of the time series signal indicates the arrival time of the signal, and the Hilbert transform is used to calculate the envelope of the signal (6), resulting in identification of the signal arrival time more clearly, as shown in Figure 14:(6)ht=ft∗1πt,zt=ft+jht,et=ft2+ht2,where ft is the time series signal, ht is the signal after Hilbert transform, zt is the analytical signal, and e(t) is the envelope of zt.Figure 14
The amplitude and envelope of receiving signals from R1.The radiation pattern of the EMAT-Rayleigh waves is shown in Figure15(a), which shows that Rayleigh waves are mainly distributed along the surface and subsurface of the material; longitudinal and shear waves are propagating obliquely into the material, and shear waves have a stronger energy than longitudinal waves. In this work, only Rayleigh waves are of interest; the beam features of Rayleigh waves are studied by means of Figure 15(b).Figure 15
The radiation pattern of EMAT-Rayleigh waves for studying beam features.
(a)
The radiation pattern of EMAT-Rayleigh waves
(b)
Radiation pattern for beam featuresThere are two beam features to be analysed, beam directivity and field distribution along the steering angle. Beam directivity is, at a radial length from the centre of the array, the velocity or pressure distribution, as shown in the green circle arc in Figure15(b). Field distribution along the steering angle is the velocity or pressure distribution along the steering angle, as shown in the red line in Figure 15(b). r denotes the radial distance from the middle of the sensor and θ denotes the steering angle; define the steering angle parallel to the surface as 0°; in order to minimise the effects of the longitudinal waves and shear waves, keep the radial distance r 140 mm and the steering angle θ from 0° to 20°. The beam directivity calculated is shown in Figure 16, where the Rayleigh waves are mainly distributed along the steering angles from 0° to 2.5°. Because the radial distance r used is 140 mm, the depth d of the Rayleigh waves’ energy distribution can be calculated by(7)d=tan⁡θ∗r.The depth d calculated is 6.1 mm, and the wavelength of the Rayleigh waves used in this work is 6.066 mm, which validates that the energy of the Rayleigh waves is mainly distributed within one wavelength of the Rayleigh waves. In addition, when choosing another value of the radial distance r, the steering angle θ changes correspondingly, but the depth d of the Rayleigh waves’ energy distribution is a constant.Figure 16
The beam directivity and field distribution along the steering angle 0° of EMAT-Rayleigh waves.Another beam feature, field distribution along the steering angle, is studied as well. The radial distancer used is 200 mm from the centre of EMAT arrays, and the steering angle used is 0° (along the surface). The field distribution along the steering angle 0° is shown in Figure 16 and the length period of the maximum magnitudes in this figure is about 17.5 mm, which is consistent with the modelling geometry, where the distance from the centre to the end of the EMAT sensor is 17.1 mm. As a result, the field energy is maximum at the places where the sensor is placed and decreases along the radial distance due to the attenuation of Rayleigh waves.


## 3. Conclusion

A method combining the analytical method for EM simulation and the finite-difference time-domain (FDTD) method for UT simulation to model an EMAT system is proposed. For electromagnetic simulation, analytical solutions to a meander coil are proposed and verified with FEM. By comparing with FEM, the analytical method proposed provides several advantages, especially at high frequencies. The calculated Lorentz force density is used as the excitation source in the ultrasonic simulation, which exploits the finite-difference time-domain (FDTD) method to describe ultrasound wave propagation. The radiation pattern shows that the maximum energy of surface waves is concentrated at the steering angle 0° and at the position where the sensor is placed.



---

*Source: 5451821-2016-01-06.xml*