Printer Friendly

A simulation study of flaw detection for rail sections based on high frequency magnetic induction sensing using the boundary element method.


The speed and loads of trains have been increasing greatly in recent years, and these factors inevitably raise the risk of producing rail defects [1]. Damages on rails increasingly originate from the surface as a result of rolling contact fatigue (RCF) [2]. Such damage can became dangerous for the operation of rail traffic, unless it is carefully monitored and managed. It is challenging to develop effective and efficient rail flaw detection techniques, which are needed to improve the safety of the railway transportation systems. Ultrasonic sensors (UT), eddy current (EC) and magnetic particle inspection (MPI) are the conventional non-destructive evaluation (NDE) techniques which are currently being used in the rail industry. On line inspection systems tend to favour ultrasonic methods and a comprehensive review of these techniques is given in [3]. Ultrasonic inspection has the best performance for detecting internal cracks [4, 5]. However, its inspection speed is limited to 75 km/h due to the transit speed of ultrasonic waves [3]. Eddy current testing which is based on the law of electromagnetic induction identifies defects using magnetic field generated by eddy currents [6]. It has relatively high inspection speed and is able to detect surface defects, so it is widely combined with ultrasonics for rail inspection.

Magnetic induction sensing technique is a non-invasive and non-contact detection approach, suitable for industrial and biomedical applications [7-12]. Magnetic induction sensing system applies a magnetic field from the excitation coil to induce eddy currents in the material and the scattering field is then detected by the receiving coil [13, 14].

For typical magnetic induction sensing technique for industrial applications, the driving frequency is normally in the range of kHz to 100 kHz. Imaging samples with low conductivities such as biological tissues or ionized water is a difficult task for the low frequency magnetic induction sensing system because the eddy currents induced in the target is very weak and the receiving coil voltage is very small. High frequency magnetic induction sensing system (> MHz) emerged to solve this problem [7, 15].

When the driving frequency is significantly high, the conductivity of the rail in the internal field can be treated as infinite and the rail behaves essentially like a perfect electric conductor (PEC). The incident magnetic field is almost totally reflected, namely the magnetic field only exists on the surface of the rail. For a PEC, the FEM would require discretisation of the volume exterior to the object between the conductor and a suitably chosen outer boundary, where a far field boundary condition should be imposed. In contrast, Boundary Element Method (BEM) based on integral formulations becomes an effective way to analyse this kind of scattering problems since meshes are only required on the surface of the object [16].

In this paper, a high frequency flaw detection system is proposed, complemented by BEM, which is based on scalar potential. When the investigated system, comprising one excitation coil and one receiving coil, scans on the surface of the rail, the primary field from the excitation coil induces currents in the rail, which in turn radiates a scattered field. We can retain a simple integral equation formulation in scalar potential for the region outside the rail, where magnetic fields are irrotational. By calculating the gradient of the scalar potential, the distributions of the magnetic field outside the rail are derived. The voltage obtained will be changed when the system comes across a flaw. Then, the flaws with different shapes and positions on the surface can be detected. The simulation results show that the proposed system has possibility to test the surface flaw of rail. Based on the simulation and mathematical analysis, the hardware system can be built in the future to verify the proposed detection strategy.


The flaw detection system consists of sensor array and a rail model with surface flaws. The rail under test is a cuboid with its center set in the origin of the coordinate system. The intercepts of the rail along X, Y and Z axis are a = 0.25 m, b = 0.5 m, c = 0.45 m, respectively. So the length, width and height are 2b = 1m, 2a = 0.5 m, 2c = 0.9 m. The sensor array contains a pair of circular coils with the same size which are used as excitation coil and receiving coil respectively. The radius of the coils is 0.1 m and the sensor array moves along the positive direction of the Y axis with a measurement range Y [member of](-0.5 m, 0.5 m).

2.1. Sensor Array

The configurations of sensor array are of three different choices:

(1) Sensor_1: The two coils are symmetrical of the XY plane, as shown in Figure 1(a). This is the most common sensor distribution in two-coil magnetic induction sensing system. The centers of the coils are (0, 0, 0.8) and (0, 0, -0.8), respectively. respectively. Of course this is not physically realistic for a deployable sensor, but is used for illustration purposes.

(2) Sensor_2: The two coils are parallel to the metallic surface and the ligature of two centers is perpendicular to the Y axis, as shown in Figure 1(b). The centers of the coils are (0.15, 0, 0.8) and (-0.15, 0, 0.8), respectively.

(3) Sensor_3: The two coils are parallel to the metallic surface, as shown in Figure 1(c). The excitation coil and receiving coil are positioned along the Y axis in tandem and the centers are (0, -0.15, 0.8) and (0, 0.15, 0.8), respectively.

2.2. Different Surface Flaws

The parameters of the flaws contains: position p, shape s, area a, depth d.

Different kinds of flaws with different characteristics are listed as follows:

(1) Flaw_1: Flaw_1 is a cuboid with the center in the Z axis. The intercepts on the X, Y axis are [a.sub.0] = a/2 and [b.sub.0] = b/2. The depth is 0.1c or 0.2c.

(2) Flaw_2: Flaw_2 is a cuboid with the center in the Z axis, as shown in Figure 2(a). The intercepts on the X, Y axis are [a.sub.0] = a/5 and [b.sub.0] = b/5. The depth is 0.1c or 0.2c.

(3) Flaw_3: Flaw_3 is a cylinder with the center in the Z axis, as shown in Figure 2(b). The radius is [r.sub.0] = 0.04m and the depth is 0.1c or 0.2c. The base area is s = [pi][r.sup.2.sub.0] = 0.005 [m.sup.2].

(4) Flaw_4: Flaw_4 is an spheroid with the center in the Z axis. The semi minor and major axis is [a.sub.0] = 0.1/[pi]m and [b.sub.0] = 0.05 m, respectively, and the depth is 0.1c or 0.2c. The base area is s = [pi][a.sub.0][b.sub.0] = 0.005 [m.sup.2].

(5) Flaw_5: Flaw_5 is a regular triangle cylinder with the center in the Z axis. The length is l = 0.107456993 m and the depth is 0.1c or 0.2c. The base area is s = 0.5[l.sup.2] sin(60[degrees]) = 0.005 [m.sup.2].

Flaw_3, Flaw_4 and Flaw_5 have the same base area. As all flaws are cylinders, they have the same volume when their depths are the same. The signals with the same volume but different shapes can be obtained.

(6) Flaw_6: Flaw_6 is a cuboid which can be seen as Flaw_2 moved along X axis for [a.sub.0] and the depth is 0.1c or 0.2c.

(7) Flaw_7: Flaw_7 and Flaw_6 are symmetrical of YZ plane.

(8) Flaw_8: Flaw_8 is a cuboid which can be seen as Flaw_2 moved along Y axis for [b.sub.0] and the depth is 0.1c or 0.2c.

(9) Flaw_9: Flaw_9 and Flaw_8 are symmetrical of XZ plane.


For a highly conducting and simply connected object in high frequency magnetic induction sensing system, where it behaves as a PEC [17], the skin depth [delta] = [square root of 2/[sigma][mu][omega]] approaches to zero as [sigma] [right arrow] [infinity]. The reference [18] has given a simple model for PEC induction problem, which is available to our research. We begin with Ampere's law, applied to the external region of the target, and Gauss's magnetic divergence law

[nabla] x [H.sup.e] = [[sigma].sup.e][E.sup.e] + i[omega][epsilon][E.sup.e] (1)

[nabla] x H = [nabla] x [[mu].sup.e][H.sup.e] = 0 (2)

Here H is the magnetic field and E the electric field. [mu], [epsilon] and [sigma] are magnetic permeability, permittivity and electrical conductivity of the media, respectively. The exterior (air) environment is assumed to have the magnetic permeability of free space. The superscript e indicates the quantities in the external environment.

The total magnetic field consists of two parts: the primary magnetic field generated by the excitation coil and the scattered magnetic field radiated by the eddy current, that is [H.sup.e] = [] + [], where the superscripts pr and sc indicate the primary field and the scattered field, respectively.

We assume that the magnetic field is quasi-static in the external region and hence the displacement current can be negligible. Also, the space exterior to the target is non-conductive so that the electric field is almost zero and the region investigated doesn't include the given current. So the right side of (1) is eliminated which means that the exterior total magnetic field is irrotational and thus can be represented efficiently using a simple scalar potential. Together with the divergence theorem and the jump condition on the normal magnetic field across the surface, the magnetic scalar potential [[phi].sup.e] satisfies a boundary integral equation written as:

c[[phi].sup.e](r) + [[integral].sub.[GAMMA]][[[phi].sup.e](r')[[partial derivative]g/(r, r'/[partial derivative]n] d[GAMMA] = [[phi]](r) (3)

where g(r, r') = [1/4[pi]][1/[absolute value of r - r']] is the Green's function for the Laplace equation, and [[phi]] is the magnetic scalar potential of the primary magnetic field. Vectors r and r' represent the observation point and source point respectively. [GAMMA] denotes the surface of the object with the unit normal vector n, and the parameter c depends on the location of observation point. c is 1/2 on a smooth boundary of [GAMMA], 1 in [GAMMA] and 0 elsewhere [19].

3.1. Numerical Implementation of the BEM

Solving (3) involves evaluating several integrals which imposes an enormous computational burden if the number of mesh is too large. There are a number of methods emerged to facilitate the computation by building a local coordinate system on the basis of which the parameter transformation is accomplished to make the integral equations more easier solved [20-22].

For numerical purposes, the surface of the target could be subdivided into NE patches with the eth patch known as [[GAMMA].sub.e]. Then the integration could be transformed into equations in discrete collocation points and NP is the number of discrete points.


We used the similar parameters and functions introduced in [21], where a local coordinate system is introduced to facilitate the calculation of 3-D numerical integrations based on Green's function or its gradient on a plane triangle. The new coordinate system parallels with triangular element and has origin on the first vertex of the triangle. One edge of the triangle is set as the abscissa, as shown in Figure 3. Here a triangular patch was taken as an example.

[V.sub.1], [V.sub.2] and [V.sub.3] are three vertexes of the triangular patch [[GAMMA].sub.e] in terms of the global coordinates. [partial derivative][S.sub.i], which is opposite to [V.sub.i], is the ith edge of [[GAMMA].sub.e]. [l.sub.i], [s.sub.i] and [m.sub.i] (i = 1, 2, 3) are the length, unit tangent vector and unit normal vector of [partial derivative][S.sub.i] respectively. Definitions of the local coordinate system in terms of the global coordinates are [zeta] = [OV.sub.2] - [OV.sub.1]/l3, [eta] = [??] x [zeta], [xi] = [??], where n is the normal vector of [[GAMMA].sub.e].

In BEM, although the accuracy increases with higher order interpolation schemes, so does the number of the nodes used to divide the boundary and formulation complexity. Hence, the number of equations in the system of linear equations and the computational time increase correspondingly with high order interpolation schemes. The first order boundary elements is less complicated than the second order boundary elements, at the same time it has higher accuracy than zero order boundary elements. For the problem we proposed, we have found that the first order boundary elements can satisfy the requirements in terms of accuracy and computational speed.

Over the surface, we interpolated [[phi].sup.e](r') with simple linear interpolation


where N is the linear nodal function and [[phi].sup.k] is the magnetic scalar potential on the kth vertex of [[GAMMA].sub.e].

So the final integral equation was given by:


Here [[phi].sup.e.sub.NPx1] are the magnetic scalar potentials required on the boundary, and [[phi]] are the magnetic scalar potentials of discrete points generated by a given magnetic field. [[delta].sub.kj] is introduced to transform the local number k into the global number j.

Let R = [absolute value of R] = [absolute value of r - r'] be the distance between an arbitrarily located observation point r and a source point r' on [[GAMMA].sub.e]. The relationships of parameters are defined as follows:

[[zeta]'.sub.a] = [zeta]' - [[zeta].sub.0]; [[eta]'.sub.a] - [eta]' - [[eta].sub.0]; [lf.sub.i] = ln([R.sup.+.sub.i] + [s.sup.+.sub.i]/[R.sup.-.sub.i] + [s.sup.-.sub.i]) (7)



The integrals of the function are now easily obtained from [21,Eqs. (26), (27)]:



With the aid of (9)-(11), it is convenient to computer the scalar potential. After that, if we know a given additional magnetic field, then we can calculate (4) to obtain the magnetic potentials of discrete points on the target surface.

3.2. Computation of the Scattered Magnetic Field and Primary Magnetic Field

Through the similar transformation and then take the gradient of both sides of (4), the magnetic field outside the target can be derived:


Here [] is the additional magnetic field and [[phi].sup.k] the magnetic potentials of discrete points on the surface of target. With the similar procedure we can solve the integral equation.

If the current on the circular excitation coil is I, the magnetic field generated by the coil can be easily derived from the Biot-Savart law:

[] = [integral] d[] = []/[[mu].sub.0] = [I/4[pi]] [[??].sub.l] [dl x R/[[absolute value of R].sup.3]] = [I/4[pi]] [[integral].sup.2[pi].sub.0] [t x R/[[absolute value of R].sup.3]]rd[theta] (13)

Note that [] = -[nabla][[phi]], the potential of the primary magnetic field is written as [23]:

[[phi]] = [I/4[pi]] [integral][[integral].sub.s][R/[[absolute value of R].sup.3]] x [n.sub.0]ds (14)

Here B is magnetic flux density and s denotes the plane enclosed by the excitation coil l. R = r - r', t is the unit tangent vector of l and [n.sub.0] is the unit normal vector of s. [[mu].sub.0] is the permeability of the free space. Further, according to Maxwell's equations, we have


where U is the voltage on the receiving coil. s', the plane enclosed by the receiving coil l', is divided into small elements. [s'.sub.i] is the area of the ith element of s' and [n'.sub.0] is the unit normal vector of s'. The angular frequency is [omega] = 2[pi]f, where f is the driving frequency and N is the number of grids of s'.


4.1. Evaluation of Different Flaw and Sensor Configurations

Simulations were carried out with the aid of MATLAB. Combine different patterns of the flaws and the sensors, we can obtain the simulation results and figures. Assume a sinusoidal current I = 1 x sin([omega]t) is applied to the excitation coil, where [omega] = 2[pi]f. Here the frequency is f = 10 MHz and the unit normal vector of the excitation coil [n.sub.0] is directed along the z axis. The direction of the current and [n.sub.0] satisfy the right-handed corkscrew rule.

(1) Flaw_1 + Sensor_1; Flaw_2 + Sensor_1.

Change the depth d = 0 (there is no flaw), d = 0.1c or d = 0.2c, the simulation results in Figure 4 show that the curves go by the same trend and they have the same shape despite they are different on the magnitude for both Flaw_1 and Flaw_2.

Flaw_1 is larger than Flaw_2 when the depth is the same. Figure 4 shows that the flaw size has an effect on the signal variance. By decreasing the depth of the flaw, there has a decrease on the signal variance. When the flaw is larger, the signal variance is larger.

(2) Flaw_2 + Sensor_2; Flaw_2 + Sensor_3.

Figures 4 and 5(a) show that, because of the symmetry of the Sensor_1 and Sensor_2 of the direction of motion, the measured voltages are symmetrical. For the same flaw, the signal variances obtained by Sensor_2 and Sensor_3 are larger than that obtained by Sensor_1.

At the same time, Figure 5(b) shows that for Sensor_3, the flaw position can influence the peak point of the curve. The peak point will move if the flaw moves.

From figures above, we know that Sensor_3 not only can produce larger signal variance, but have a relationship between the peak point of the curve and flaw position, so the simulations following will use Sensor_3.

(3) Flaw_3 + Sensor_3; Flaw_4 + Sensor_3; Flaw_5 + Sensor_3. Figure 6 shows that, if the flaw volume and the position are the same, the shapes of the curves are the same on the whole with a small difference in value.

(4) Flaw_6 + Sensor_3; Flaw_7 + Sensor_3 (shown in Figure 7).

(5) Flaw_8 + Sensor_3; Flaw_9 + Sensor_3 (shown in Figure 8).

4.2. The Relationship of the Peak Point and the Flaw Position for Sensor_3

If other conditions such as flaw area, depth (here d = 0.2c) and shape are the same, the curves obtained by different positions with Sensor_3 are described as Figure 9.

(1) Curves of Flaw_2, 6 and 7

The curves in Figure 9(a) present that when the flaw moves only along axis X, the peak point of the curve doesn't change with a little difference in the magnitude of data.

(2) Curves of Flaw_2, 8 and 9.

4.3. Curves with Different Flaw Shape

The measurement curves of Flaw_3, 4 and 5 are obtained with Sensor_3, as shown in Figure 10. The result indicates that flaws with different shapes but same volume are likely to produce similar curves.


5.1. The Effect of Different Flaws on the Results

From the simulation results we can summarize how the flaw parameters affect the curves, and here we use the single variable analysis.

(1) Flaw volume: Figure 4 shows how the result changes with the flaw volume. When the flaw is larger, the signal variance is larger.

(2) Flaw position: For Sensor_1 and 2, the position has no influence in the curve shape. For Sensor_3, the peak point of curves changes when the flaw moves along Y axis. Figure 9 shows that when the flaw moves in the positive direction of axis Y, the peak point will move to the right.

(3) Flaw shape: Figure 10 shows that when material conditions except for the shape are same, the curves of different flaws are basically same.

5.2. The Effect of the Sensor Distribution on the Results

Figures 4(b) and 5 show that, signal variance for distinguishing different flaws obtained by Sensor_3 are larger than that obtained by Sensor_1 and Sensor_2. At the same time, different flaw positions will contribute to different peak points of curves for Sensor_3. So Sensor_3 is the best choice for flaw detection.


The paper presented a BEM model based on a simple high frequency magnetic induction sensing system to inspect the surface flaw of rail which can be treated as a PEC. The magnetic fields were obtained numerically and the voltages on the receiving coil were calculated for different kinds of flaws. In the simulation, the curves obtained by the BEM are plotted in MATLAB, then how the flaw features and the sensor array affect the results was discussed. The simulation results of flaw detection system, with the help of the BEM, suggest that magnetic induction sensing technique has the potential to inspect the surface flaw of the rail. Sensor_3 is the best configuration for the measurement of flaws because it not only can effectively characterize the flaw volume, but also is sensitive to other parameters such as the flaw position. The flaw volume has a great effect on the magnitude of the curve while the peak point of the curve is determined by the flaw position. Future work will involve building an experimental system to verify the observations obtained through the mathematical analysis in this paper.


The authors would like to express their gratitude to Professor Anthony Peyton for his encouragements and guidance in the field of MIT.


[1.] Li, Q. Y. and S. W. Ren, "A real-time visual inspection system for discrete surface defects of rail heads," IEEE Transactions on Instrumentation and Measurement, Vol. 61, 2189-2199, 2012.

[2.] Rowshandel, H., G. L. Nicholson, C. L. Davis, and C. Roberts, "A robotic system for non-destructive evaluation of RCF cracks in rails using an ACFM sensor," 5th IET, 29-30, 2011.

[3.] Papaelias, M., C. Roberts, and C. L. Davis, "A review on non-destructive evaluation of rails: State-of-the-art and future development," Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit, Vol. 222, 367-384, 2008.

[4.] Clark, R., S. Singh, and C. Haist, "Ultrasonic characterisation of defects in rails," Insight, Vol. 44, 341-347, 2002.

[5.] Edwards, R., S. Dixon, and X. Jian, "Characterisation of defects in the railhead using ultrasonic surface waves," NDT & E. Int., Vol. 39, 468-475, 2006.

[6.] Cacciola, M., F. C. Morabito, D. Polimeni, and M. Versaci, "Fuzzy characterization of flawed metallic plates with eddy current tests," Progress In Electromagnetics Research, Vol. 72, 241-252, 2007.

[7.] Watson, S., R. J. Williams, W. A. Gough, and H. Griffiths, "A magnetic induction tomography system for samples with conductivities less than 10S[m.sup.-1]," Measurement Science & Technology, Vol. 19, 045501, 2008.

[8.] Yin, W. L. and A. J. Peyton, "Simultaneous measurements of thickness and distance of a thin metal plate with an electromagnetic sensor using a simplified model," IEEE Transactions on Instrumentation and Measurement, Vol. 53, 1335-1338, 2004.

[9.] Yin, W. L., A. J. Peyton, G. Zysko, and R. Denno, "Simultaneous noncontact measurement of water-level and conductivity," IEEE Transactions on Instrumentation and Measurement, Vol. 57, 2665-2669, 2008.

[10.] Ma, L., H.-Y. Wei, and M. Soleimani, "Planar magnetic induction tomography for 3D near subsurface imaging," Progress In Electromagnetics Research, Vol. 138, 65-82, 2013.

[11.] Wei, H.-Y. and M. Soleimani, "Three-dimensional magnetic induction tomography imaging using a matrix free Krylov subspace inversion algorithm," Progress In Electromagnetics Research, Vol. 122, 29-45, 2012.

[12.] Wei, H.-Y. and M. Soleimani, "Four dimensional reconstruction using magnetic induction tomography: Experimental study," Progress In Electromagnetics Research, Vol. 129, 17-32, 2012.

[13.] Ma, X., A. J. Peyton, S. R. Higson, A. Lyons, and S. J. Dickinson, "Hardware and software design for an electromagnetic electromagnetic induction tomography (EMT) system for high contrast metal process applications," Measurement Science & Technology, Vol. 17, 111-118, 2006.

[14.] Griffiths, H., "Magnetic induction tomography," Measurement Science & Technology, Vol. 12, 1126-1131, 2001.

[15.] Wei, H.-Y. and M. Soleimani, "Two-phase low conductivity flow imaging using magnetic induction tomography," Progress In Electromagnetics Research, Vol. 131, 99-115, 2012.

[16.] Wu, K. L., G. Y. Delisle, D. G. Fang, and M. Lecours, "Coupled finite element and boundary element methods in electromagnetic scattering," Progress In Electromagnetics Research, Vol. 02, 113-132, 1990.

[17.] Liao, S. and R. J. Vernon, "On the image approximation for electromagnetic wave propagation and PEC scattering in cylindrical harmonics," Progress In Electromagnetics Research, Vol. 66, 65-88, 2006.

[18.] Sun, K. L., K. O'Neill, F. Shubitidze, S. A. Haider, and K. D. Paulsen, "Simulation of electromagnetic induction scattering from targets with negligible to moderate penetration by primary fields," IEEE Transactions on Geoscience and Remote Sensing, Vol. 40, 910-927, 2002.

[19.] Pham, M. H. and A. J. Peyton, "A model for the forward problem in magnetic induction tomography using boundary integral equations," IEEE Transactions on Magnetics, Vol. 44, 2262-2267, 2008.

[20.] Morrison, J. A., "Integral equations for electromagnetic scattering by perfect conductors with two-dimensional geometry," Bell Syst. Tech. J., Vol. 58, 409-425, 1979.

[21.] Graglia, R. D., "On the numerical integration of the linear shape functions times the 3-D Green's function or its gradient on a plane triangle," IEEE Transactions on Antennas and Propagation, Vol. 41, 1448-1455, 1993.

[22.] Graglia, R. D., "Static and dynamic potential integrals for linearly varying source distributions in two- and three-dimensional problems," IEEE Transactions on Antennas and Propagation, Vol. 35, 662-669, 1987.

[23.] Zhang, Z. M. and Y. R. Den, "A new method using Biot-Savart law to derive magnetic scalar potential notation," Journal of Chongqing Institute of Civil Engineering and Architecture, Vol. 4, 99-103, 1985.

Qian Zhao (1), *, Jian Na Hao (1), and Wu Liang Yin (2)

(1) School of Electrical Engineering and Automation, Tianjin University, Tianjin 300072, China

(2) School of Electrical and Electronic Engineering, University of Manchester, Manchester M60 1QD, UK

Received 27 April 2013, Accepted 31 May 2013, Scheduled 22 July 2013

* Corresponding author: Qian Zhao (
COPYRIGHT 2013 Electromagnetics Academy
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2013 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Zhao, Qian; Hao, Jian Na; Yin, Wu Liang
Publication:Progress In Electromagnetics Research
Article Type:Report
Geographic Code:9CHIN
Date:Sep 1, 2013
Previous Article:A shaped-beam series-fed aperture-coupled stacked patch array antenna.
Next Article:A near-field 3D circular SAR imaging technique based on spherical wave decomposition.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters