# Mode Selection Model for Rail Crack Detection Based on Ultrasonic Guided Waves

Authors: Bo Xing; Zujun Yu; Xining Xu; Liqiang Zhu; Hongmei Shi
Journal: Shock and Vibration
 (2020)

---

## Abstract

The cross-sectional geometry of a rail is complex, and numerous guided wave modes can be propagated in rails. In order to select a mode which is the most suitable for detecting a specific crack on a rail, a mathematical model of guided wave mode selection is constructed. The model is composed of a modal vibration factor and a modal orthogonal factor. By setting a reasonable vibration coefficient and orthogonal coefficient, the mode with the highest sensitivity to cracks is selected for crack detection. Taking a vertical crack on the rail bottom as an example, mode 1 at a frequency of 60 kHz is selected as the most suitable detection mode. At the same time, mode 7 and mode 11 are selected as comparative modes, and these three modes are simulated to detect rail cracks. Among them, mode 1 is the best, which verifies the correctness of the mode selection model. In addition, vertical cracks are manufactured artificially on the side of the rail bottom. The cracks are successfully detected by mode 1, and the positioning error is 0.07 m. After correction, the error is reduced to 0.02 m. The model can effectively select guided wave modes suitable for detecting arbitrary cracks on rails, which provides a theoretical solution for rail crack detection.


---

## Body

## 1. Introduction

Ultrasonic guided waves are widely used in rail integrity detection because of the advantages of long transmission distance and wide coverage [1–3]. The rail cross-sectional geometry is complex. At the same time, the higher the frequency of the wave is, the shorter its wavelength is, the stronger the particle character is, and the stronger the echo produced when the wave meets with the crack. Compared with large cracks, the reflective wave intensity produced by small cracks is relatively low, and the guided wave of low frequency and long wavelength may not produce detectable reflective echo when encountering small cracks. Therefore, rail defect detection needs to be analyzed under high-frequency conditions. Traditional analytical methods, such as the method of establishing multiple complex equations, cannot analyze the propagation regularity of ultrasonic guided waves in rails. For this reason, Gavric [4] introduced the numerical method and proposed the semianalytical finite element (SAFE) method for the first time. SAFE method was originally used for numerical modeling of guided wave propagation to replace the precise solution method based on the superposition of partial bulk waves (SPBW). When using the SAFE method to extract frequency dispersion, finite element discretization is only needed on the cross section of waveguide, while the displacement along the propagation direction of waveguide is expressed by simple harmonic vibration. Since then, SAFE has been widely used in many studies [5–9]. In this paper, the SAFE method is used to calculate the mode shape data of guided wave modes.Lee et al. [10] selected several modal points on phase velocity dispersion curves and obtained the energy distribution of four modes in a rail at different frequencies (30–315 kHz) through a large number of three-dimensional model simulation calculations. The energy distribution regularity of each mode propagating in the rail was obtained, and the appropriate frequencies and modes were selected to detect the transverse cracks of the rail head, which guided the design of sensors. Long and Loveday [3, 11] applied the SAFE method to solve the mode shapes of guided waves propagating in rails in the frequency range of 25–45 kHz. According to the vibration situation of different regions in the modal diagram, four modes that can be propagated over a long distance were selected in the mode shapes at 30 kHz frequency. They also use a hybrid SAFE-FE method to model the defects in the rail, which combines 2D SAFE models of the semi-infinite waveguide and a 3D solid FE model of defect only. The modes suitable for detecting defects at the rail head and rail web were selected from four modes. Coccia et al. [12] calculated the strain energy of different modes of ultrasonic guided waves in rails. The distribution of strain energy of some specific modes at 200 kHz frequency was solved, and suitable modes for detecting defects in different areas of the rail head were obtained. From different points of view, all these studies obtained suitable modes for crack detection of rail heads. The results are different because of the different frequency ranges. However, there is no detailed study on the selection of detection modes for different types of cracks at all locations. Therefore, the focus of this paper is to establish a guided mode selection model for rail crack detection that can be applied to cracks of any angle and size in all regions of the rail.In this paper, a mathematical model of guided wave mode selection is established to select the modes suitable for detecting cracks in all regions of rails. First, the SAFE method is applied to analyze the propagation characteristics of ultrasonic guided waves in rails, and the relationship between guided waves and crack detection is discussed in detail. Then, a mathematical model is constructed based on the orthogonal relationship between directions of modal vibration at the crack location and the shape of the crack. In addition, taking a vertical crack on the rail bottom as an example, appropriate and comparative modes are selected to detect the crack with the established mathematical model. The validity of the mode selection model is verified by simulations and experiments.In the future, the research of mechanical vibration monitoring system [13] will be added, and wireless sensor network will be established to realize the real-time monitoring of rail defects.


## 2. Mode Shapes

Because of the low accuracy of low-frequency ultrasonic guided waves for rail crack detection, conventional crack detection needs to be carried out under high-frequency conditions. In some waveguides with simple cross-sectional geometry, many modes can propagate at a high frequency. However, due to the complexity of the rail cross-sectional geometry, much more guided wave modes can be propagated in the rail under high-frequency conditions. Therefore, it is impossible to analyze the guided wave propagation in the rail by simple analytical calculation. The numerical methods, such as semianalytical finite element (SAFE) method, can be used to describe the dispersion curves of ultrasonic guided waves in rails, obtain the relationship between guided wave modes and frequencies, and provide data for the construction of a modal selection model [4, 7]. Taking CHN60 rail as an example, its coordinate system is established, as shown in Figure 1. The x-axis is the longitudinal direction of rail, the y-axis is the lateral direction of rail, and the z-axis is the vertical direction of rail.Figure 1
CHN60 rail coordinate system.The relationship between displacementu and strain ε at any point in the rail can be seen as follows [14]:(1)ε=Lx∂u∂x+Ly∂u∂y+Lz∂u∂z,among which(2)Lx=100000000000001010,Ly=000010000001000100,Lz=000000001010100000.When the SAFE method is used to calculate the dispersion curves of ultrasonic guided waves in rails, it is assumed that they propagate along the rails in the form of simple harmonics, which can be expressed instead ofeiξx−ωt [4]. The displacement of an elastic waveguide extending to infinity in the x-direction is expressed as(3)ux,y,z,t=uxx,y,z,tuyx,y,z,tuzx,y,z,t=Uxy,zUyy,zUzy,zeiξx−ωt.The triangular element was used to discretize the rail cross section of CHN60. There are totally 177 nodes and 255 elements generated. The discrete cross section is shown in Figure2.Figure 2
Discretization of cross section of CHN60 rail.According to the shape function of the elements, the displacement, strain, and stress of any part can be obtained. The equation of guided wave propagation in rails is that of Gavric, as proposed in [4]:(4)K1+iξK2+ξ2K3−ω2MMU=0,where ξ is the wavenumber of the guided wave; ω is the angular frequency; K1,K2,andK3 are the node stiffness matrix; M is the node mass matrix; and U is the node displacement matrix.Given the angular frequencyω, which is in the range of 2π× [0, 100] kHz, the eigenvalue equation is solved by equation (4). The eigenvalue wavenumber ξ and the eigenvector U^ contain the rail mode shapes. The frequency wavenumber curves of guided wave modes in rails are plotted as shown in Figure 3.Figure 3
Frequency wavenumber curves of ultrasonic guided waves in rails.As can be seen in Figure3, with increased frequency, the corresponding wavenumber points and the number of modes increase at the same frequency. At the frequency critical point where a new mode is generated, the wavenumber of the mode is the smallest of all the modes at this frequency. For example, at the frequency of 60 kHz in the figure, the selected one is a new mode, which appears from the frequency point and exists in all frequencies whose value is larger than that frequency.The phase velocityCp and group velocity Cg of modes can be calculated by the following equations, respectively:(5)Cp=ωξ,(6)Cg=dωdξ.Thus, the dispersion curves of phase velocity and group velocity of ultrasonic guided waves in rails can be drawn, as shown in Figure4.Figure 4
Dispersion curves of ultrasonic guided waves in rails: (a) phase velocity; (b) group velocity.
(a)(b)With the increase in frequency, the number of modes increases greatly. Multiple modes at the same frequency are numbered in the order of phase velocity from small to large. According to equation (5), at the same frequency, the wavenumber is inversely proportional to the phase velocity. The smaller the wavenumber, the larger the phase velocity. Therefore, moving along the transverse axis in the positive direction in Figure 3, a new mode has the smallest wavenumber and the largest phase velocity when it is generated. There is no direct proportional relationship between the group velocity and the wavenumber or phase velocity.At the same time, according to the eigenvector obtained by solving equation (4), the modal data of each mode can be extracted. Taking guided waves with a frequency of 60 kHz as an example, the mode shapes of 35 modes are shown in Figure 5.Figure 5
Mode shapes diagram (60 kHz).Among the first 10 modes, the vibration energy of modes 1, 2, 4, 5, and 9 is concentrated on the bottom of the rail, but the vibration positions on the bottom and the torsion forms are different. Modes 1, 2, and 9 are more suitable for the detection of corner cracks at the bottom. Modes 4 and 5 also have some vibration in the middle of the rail bottom. At the same time, the mode with the torsional direction perpendicular to the shape direction of the crack has a larger reflection coefficient, which is easier to detect in practical application. Similarly, the vibration energy of modes 3 and 6 is concentrated on the rail web, which is more suitable for detecting cracks in the rail web. The vibration energy of mode 7 is concentrated on the rail head, which is more suitable for detecting cracks in the rail head. However, the vibration energy of modes 8 and 10 is dispersed, which is not the preferred mode for crack detection.


## 3. Establishing a Mode Selection Model

In order to detect rail cracks, it is necessary to select suitable modes for detecting corresponding cracks. The number of modes that can be propagated in rails is related to the frequency. The greater the frequency is, the more modes there are, and the more difficult it is to study. When the frequency is more than 60 kHz, the guided wave mode propagating in the rail tends to be a surface wave, which is not conducive to the detection of internal cracks [10]. When the frequency is lower than 20 kHz, the accuracy of crack detection is lower. In order for the mathematical model to be applied to the selection of modes for crack detection of the whole rail, taking into account both the difficulty in research and the accuracy of detection, the frequency range considered in this paper is limited to [20, 60] kHz.In order to locate cracks in the rail by using ultrasonic guided waves, it is necessary to analyze the reflection signals generated by the interaction between guided waves and cracks. Therefore, the modes sensitive to specific cracks should be selected.Crack sensitivity W is defined to represent the sensitivity of different modes to a specific crack. The sensitivity of a mode to a crack is affected not only by the vibration of the mode at the crack but also by the orthogonality of the mode vibration and the crack trend. They have different effects on crack sensitivity, so it is necessary to allocate their respective weight coefficient. Crack sensitivityW is expressed by the following equation:(7)W=α×W1+β×W2,where W1 is the modal vibration factor; W2 is the modal orthogonal factor; and α and β are the vibration coefficient and orthogonal coefficient, respectively, and α+β=1.The crack sensitivityW of all modes is calculated, and the most sensitive mode is the most suitable to detect the crack.


### 3.1. Modal Vibration FactorW1

In order to eliminate the influence of defects in other parts of a rail, it is hoped that the mode with obvious vibration at the crack location and no obvious vibration at other locations will be selected when using ultrasonic guided waves to detect specific cracks. The modal vibration factor is defined to represent the vibration of each mode near the crack location.First, the corresponding nodes near the crack location should be selected to calculate their average vibration displacement in different modes. According to the coordinate system established in Figure1, it is assumed that the crack ends are AxA,yA,zAandBxB,yB,zB. According to the meshing method in Figure 2, the nearest nodes A′xA′,yA′,zA′ and B′xB′,yB′,zB′ on the cross section of points A and B are found according to (8). Na is the number of nodes of the rail cross section:(8)miny1−yA2+z1−zA2,y2−yA2+z2−zA2,⋯,yNa−yA2+zNa−zA2,miny1−yB2+z1−zB2,y2−yB2+z2−zB2,⋯,yNa−yB2+zNa−zB2.TakeC′yC′,zC′ as the center of a circle, which is the midpoint of A'B', and make the radius of the circle r. Making a circle as shown in Figure 6, the mode shapes of Nc nodes covered by the circle are selected to represent the vibration magnitude of the crack, with the node number Ns_1,Ns_2,…,Ns_Nc.Figure 6
Selected nodes schematic.Among them,(9)yC′=yA′+yB′2,zC′=zA′+zB′2,r=yA′−yB′2+zA′−zB′2.According to equation (4), the mode shapes of all nodes in three directions of each mode can be solved, which are expressed by MSxfi,N,m, MSyfi,N,m, and MSzfi,N,m, where fi is the frequency, m is the mode number at frequency fi, and N is the node number. The average vibration displacements of the joints near the crack V1m and the whole section V2m are, respectively, as follows:(10)V1mfi=∑N=Ns_1Ns_NcMSxfi,N,m2+MSyfi,N,m2+MSzfi,N,m2Nc,where,N=Ns_1,Ns_2,⋯,Ns_Nc.V2mfi=∑N=1NaMSxfi,N,m2+MSyfi,N,m2+MSzfi,N,m2Na.The corresponding modal vibration factor is expressed by equation (11), which is used to characterize the vibration trend of each mode in the region near the crack. When the average vibration displacement of the nodes near the crack is less than that of the whole section, that is, V1mfi<V2mfi, W1mfi<0, and when the average vibration displacement of the nodes near the crack increases to much greater than that of the whole section, that is, V1mfi>V2mfi, W1mfi approaches 1:(11)W1mfi=V1mfi−V2mfiV1mfi.The greater the modal vibration factorW1m(fi), the greater the modal vibration near the crack, and the greater the sensitivity to the crack. When W1m(fi) approaches 1, the mode vibration near the crack is much greater than that at other locations. Therefore, the mode whose W1m(fi) is larger should be selected to detect rail cracks.


### 3.2. Modal Orthogonal FactorW2

In the process of receiving signals of rail crack detection, the results will be affected by the acquisition equipment and field noise. If the amplitude of the received crack reflection signal is low, there may be missed detection. Therefore, the mode used to detect rail cracks should produce relatively large crack reflections as much as possible. In the extreme case, the mode whose vibration direction is orthogonal to the crack trend has a larger reflection amplitude than the mode whose vibration direction is parallel to the crack trend. If the vibration direction and crack trend are neither orthogonal nor parallel, the amplitude will be between them. According to the vibration regularity of modes, the mode orthogonal factorW2 is defined to represent the degree of orthogonality between the vibration direction of modes and the trend of cracks. When the frequency is fi, the orthogonal factor of mode m is(12)W2mfi=1−Q×P,where Q is the direction vector of vibration and P is the direction vector of the crack trend. When the direction of the crack is orthogonal to the direction of the mode vibration, Q×P=0, W2mfi=1. When the direction of the crack is parallel to the direction of the mode vibration, Q×P is the largest, and W2mfi is the smallest. Therefore, the more W2mfi approaches 1, the better the orthogonality of the two.In order to prevent a sudden change of results, some additional nodes should be added to the vector to make the results smoother. According to the finite element discretization method, the mode shapesI of any point in the discretized element can be expressed as follows:(13)Iex,y,z,t=∑k=1nNky,zIxk∑k=1nNky,zIyk∑k=1nNky,zIzkeeiξx−ωt=Ny,zqeeiξx−ωt,where Ny,z is the shape function matrix, qe is the nodal mode shape vector, the parameters with superscript are element parameters, and the degree of freedom n = 3 when triangular elements are used.(14)Ny,z=N1N2⋅⋅⋅NnN1N2⋅⋅⋅NnN1N2⋅⋅⋅Nn,(15)qe=Ix1Iy1Iz1Ix2Iy2Iz2⋯IxnIynIzn.The shape function of the finite elementNi is presented, and the triangular element composed of nodes i,j,k is shown in Figure 7.Figure 7
Triangle element.The coordinates of each node areyi,zi, yj,zj, and yk,zk, and the shape functions are as follows:(16)Ni=12Sαi+βiy+δiz,Nj=12Sαj+βjy+δjz,Nk=12Sαk+βky+δkz,where S is the area of the triangle and αi,βi,δi are the intermediate coefficients. (17)S=12yizj−zk+yjzk−zi+ykzi−zj,αi=yjzk−ykzj,βi=zj−zk,δi=yk−yj,αj=ykzi−yizk,βj=zk−zi,δj=yi−yk,αk=yizj−yjzi,βk=zi−zj,δk=yj−yi.When the frequency isfi, the mode shapes of mode m at points A′ and B′ are recorded as IA′ and IB′, respectively. The line segment between A′B′ is divided into four parts. The middle three parts are A1′xA1′,yA1′,zA1′, A2′xA2′,yA2′,zA2′, and A3′xA3′,yA3′,zA3′. The known vibration mode shapes of the three vertices of the triangular unit where the three parts are can be applied to equations (13)–(15). Then, the mode shapes of these three parts expressed by IA1′, IA2′, and IA3′ are solved, together with IA′ and IB′, forming a vibration vector, which is recorded as Q:(18)Q=IA′,IA1′,IA2′,IA3′,IB′.The crack line segment AB is divided into four parts. The middle three points areA1xA1,yA1,zA1, A2xA2,yA2,zA2, and A3xA3,yA3,zA3, together with A and B, and a coordinate vector is composed, which is recorded as P:(19)P=A,A1,A2,A3,BT.If the direction of mode vibration is orthogonal to the crack trend |Q×P| = 0, W2 = 1. If the direction of mode vibration is not orthogonal to the crack trend |Q×P| > 0, W2 < 1. Therefore, the more the orthogonal factor W2 tends to be 1, the more orthogonal the mode is. When choosing guided wave modes for rail crack detection, the mode with a larger orthogonal factor W2 should be selected.When the frequency isfi, given the mode vibration coefficient α and the mode orthogonal coefficient β, the crack sensitivity of modes can be solved by using equation (7). The detection of rail cracks can be realized by selecting modes with high sensitivity. It should be noted that the value of α here should be greater than that of β, which will be proved in Section 5.


## 3.1. Modal Vibration FactorW1

In order to eliminate the influence of defects in other parts of a rail, it is hoped that the mode with obvious vibration at the crack location and no obvious vibration at other locations will be selected when using ultrasonic guided waves to detect specific cracks. The modal vibration factor is defined to represent the vibration of each mode near the crack location.First, the corresponding nodes near the crack location should be selected to calculate their average vibration displacement in different modes. According to the coordinate system established in Figure1, it is assumed that the crack ends are AxA,yA,zAandBxB,yB,zB. According to the meshing method in Figure 2, the nearest nodes A′xA′,yA′,zA′ and B′xB′,yB′,zB′ on the cross section of points A and B are found according to (8). Na is the number of nodes of the rail cross section:(8)miny1−yA2+z1−zA2,y2−yA2+z2−zA2,⋯,yNa−yA2+zNa−zA2,miny1−yB2+z1−zB2,y2−yB2+z2−zB2,⋯,yNa−yB2+zNa−zB2.TakeC′yC′,zC′ as the center of a circle, which is the midpoint of A'B', and make the radius of the circle r. Making a circle as shown in Figure 6, the mode shapes of Nc nodes covered by the circle are selected to represent the vibration magnitude of the crack, with the node number Ns_1,Ns_2,…,Ns_Nc.Figure 6
Selected nodes schematic.Among them,(9)yC′=yA′+yB′2,zC′=zA′+zB′2,r=yA′−yB′2+zA′−zB′2.According to equation (4), the mode shapes of all nodes in three directions of each mode can be solved, which are expressed by MSxfi,N,m, MSyfi,N,m, and MSzfi,N,m, where fi is the frequency, m is the mode number at frequency fi, and N is the node number. The average vibration displacements of the joints near the crack V1m and the whole section V2m are, respectively, as follows:(10)V1mfi=∑N=Ns_1Ns_NcMSxfi,N,m2+MSyfi,N,m2+MSzfi,N,m2Nc,where,N=Ns_1,Ns_2,⋯,Ns_Nc.V2mfi=∑N=1NaMSxfi,N,m2+MSyfi,N,m2+MSzfi,N,m2Na.The corresponding modal vibration factor is expressed by equation (11), which is used to characterize the vibration trend of each mode in the region near the crack. When the average vibration displacement of the nodes near the crack is less than that of the whole section, that is, V1mfi<V2mfi, W1mfi<0, and when the average vibration displacement of the nodes near the crack increases to much greater than that of the whole section, that is, V1mfi>V2mfi, W1mfi approaches 1:(11)W1mfi=V1mfi−V2mfiV1mfi.The greater the modal vibration factorW1m(fi), the greater the modal vibration near the crack, and the greater the sensitivity to the crack. When W1m(fi) approaches 1, the mode vibration near the crack is much greater than that at other locations. Therefore, the mode whose W1m(fi) is larger should be selected to detect rail cracks.


## 3.2. Modal Orthogonal FactorW2

In the process of receiving signals of rail crack detection, the results will be affected by the acquisition equipment and field noise. If the amplitude of the received crack reflection signal is low, there may be missed detection. Therefore, the mode used to detect rail cracks should produce relatively large crack reflections as much as possible. In the extreme case, the mode whose vibration direction is orthogonal to the crack trend has a larger reflection amplitude than the mode whose vibration direction is parallel to the crack trend. If the vibration direction and crack trend are neither orthogonal nor parallel, the amplitude will be between them. According to the vibration regularity of modes, the mode orthogonal factorW2 is defined to represent the degree of orthogonality between the vibration direction of modes and the trend of cracks. When the frequency is fi, the orthogonal factor of mode m is(12)W2mfi=1−Q×P,where Q is the direction vector of vibration and P is the direction vector of the crack trend. When the direction of the crack is orthogonal to the direction of the mode vibration, Q×P=0, W2mfi=1. When the direction of the crack is parallel to the direction of the mode vibration, Q×P is the largest, and W2mfi is the smallest. Therefore, the more W2mfi approaches 1, the better the orthogonality of the two.In order to prevent a sudden change of results, some additional nodes should be added to the vector to make the results smoother. According to the finite element discretization method, the mode shapesI of any point in the discretized element can be expressed as follows:(13)Iex,y,z,t=∑k=1nNky,zIxk∑k=1nNky,zIyk∑k=1nNky,zIzkeeiξx−ωt=Ny,zqeeiξx−ωt,where Ny,z is the shape function matrix, qe is the nodal mode shape vector, the parameters with superscript are element parameters, and the degree of freedom n = 3 when triangular elements are used.(14)Ny,z=N1N2⋅⋅⋅NnN1N2⋅⋅⋅NnN1N2⋅⋅⋅Nn,(15)qe=Ix1Iy1Iz1Ix2Iy2Iz2⋯IxnIynIzn.The shape function of the finite elementNi is presented, and the triangular element composed of nodes i,j,k is shown in Figure 7.Figure 7
Triangle element.The coordinates of each node areyi,zi, yj,zj, and yk,zk, and the shape functions are as follows:(16)Ni=12Sαi+βiy+δiz,Nj=12Sαj+βjy+δjz,Nk=12Sαk+βky+δkz,where S is the area of the triangle and αi,βi,δi are the intermediate coefficients. (17)S=12yizj−zk+yjzk−zi+ykzi−zj,αi=yjzk−ykzj,βi=zj−zk,δi=yk−yj,αj=ykzi−yizk,βj=zk−zi,δj=yi−yk,αk=yizj−yjzi,βk=zi−zj,δk=yj−yi.When the frequency isfi, the mode shapes of mode m at points A′ and B′ are recorded as IA′ and IB′, respectively. The line segment between A′B′ is divided into four parts. The middle three parts are A1′xA1′,yA1′,zA1′, A2′xA2′,yA2′,zA2′, and A3′xA3′,yA3′,zA3′. The known vibration mode shapes of the three vertices of the triangular unit where the three parts are can be applied to equations (13)–(15). Then, the mode shapes of these three parts expressed by IA1′, IA2′, and IA3′ are solved, together with IA′ and IB′, forming a vibration vector, which is recorded as Q:(18)Q=IA′,IA1′,IA2′,IA3′,IB′.The crack line segment AB is divided into four parts. The middle three points areA1xA1,yA1,zA1, A2xA2,yA2,zA2, and A3xA3,yA3,zA3, together with A and B, and a coordinate vector is composed, which is recorded as P:(19)P=A,A1,A2,A3,BT.If the direction of mode vibration is orthogonal to the crack trend |Q×P| = 0, W2 = 1. If the direction of mode vibration is not orthogonal to the crack trend |Q×P| > 0, W2 < 1. Therefore, the more the orthogonal factor W2 tends to be 1, the more orthogonal the mode is. When choosing guided wave modes for rail crack detection, the mode with a larger orthogonal factor W2 should be selected.When the frequency isfi, given the mode vibration coefficient α and the mode orthogonal coefficient β, the crack sensitivity of modes can be solved by using equation (7). The detection of rail cracks can be realized by selecting modes with high sensitivity. It should be noted that the value of α here should be greater than that of β, which will be proved in Section 5.


## 4. Simulation Verification

### 4.1. Selection of Detection Modes

Taking the vertical crack on the rail bottom as an example, the feasibility of the modal selection model is verified. The crack’s location is the same as the line shown in Figure6. The two points A and B coincide with A′ and B′, numbered 29 and 32, and the longitudinal coordinates of the two ends of the crack are zero.The frequency range is [20, 60] kHz. Five frequencies were selected: 20 kHz, 30 kHz, 40 kHz, 50 kHz, and 60 kHz. According to equation (4), all modes at these five frequencies are solved, totaling 122 modes.Equations (9)–(11) are solved to obtain the vibration factors W1 of all modes, which are shown in Figure 8 after normalization.Figure 8
The vibration factorsW1.The factors are sorted from large to small, as shown in Table1. Only the first 20 modes are listed.Table 1
Modal vibration factorW1.
Frequency (Hz)Mode no.W160,0002160,0001150,00020.9950,00010.9940,00020.9940,00010.9830,00020.9830,00010.9720,00020.9720,00010.9660,000110.9550,000110.9320,000120.9240,000100.9030,000100.8660,000310.8650,00050.8530,00040.8530,00070.8440,00040.84From Table1, it can be seen that the vibrations of modes 1 and 2 at [30, 60] kHz are very obviously near 29 and 32 nodes on the rail bottom, and the results are consistent with those in Figure 5. In the modes at a frequency of 20 kHz, there is no obvious vibration mode near 29 and 32 nodes.According to equations (12)–(19), the modal orthogonal factors W2 of all modes are calculated, as shown in Figure 9.Figure 9
Modal orthogonal factorsW2.These factors are sorted from large to small, as shown in Table2. Only the first 20 modes are listed.Table 2
Modal orthogonal factorW2.
Frequency (Hz)Mode no.W240,00070.9999850,00070.9999660,00070.9999660,00030.9995850,00030.9987130,00050.9985340,000150.9978440,000230.9970940,00030.9963560,000140.9959850,000280.9901860,00060.9899550,000140.9892540,000240.9890330,00030.9883850,00060.9831220,00040.9804230,000160.9775360,000340.9756260,000190.97243Table2 shows that the vibration directions of mode 7 at [40, 60] kHz have the best orthogonality with the vertical crack on the rail bottom, but the corresponding modal vibration factors of these three modes rank 120–122, respectively, and the vibration magnitude is the worst. According to Figure 5, mode 7 at a frequency of 60 kHz has only rail head vibration, which is consistent with the calculation results, while modes 1 and 2 at a frequency of 60 kHz have the largest vibration factors and the most obvious vibration. The orthogonal factors of these two modes are 0.9327 and 0.9326, respectively, which rank 43 and 44. The orthogonality of these two modes is good. There is little difference among them and the orthogonality of the mode 7 at a frequency of 40 kHz.Letα = 0.7 and β = 0.3. The crack sensitivity W of each mode is solved according to equation (7). The results are shown in Figure 10.Figure 10
Crack sensitivityW.After sorting from large to small, the first 20 modes are listed in Table3.Table 3
Crack sensitivityW.
Frequency (Hz)Mode no.W60,00020.835460,00010.835250,00020.830350,00010.830040,00020.823640,00010.823130,00020.813930,00010.813220,00020.797520,00010.796150,00050.766530,00070.761930,00040.760940,00040.758360,00050.758140,00050.757950,00040.754860,00040.751760,000120.744140,000210.7363From Table3, it can be seen that modes 1 and 2 at a frequency of 60 kHz have the best sensitivity to the vertical crack on the bottom of the rail, and there is little difference between them. Since mode 1 is symmetrical and mode 2 is antisymmetrical, the excitation method of mode 1 is relatively simple in the process of single-mode excitation, so mode 1 is selected as the detection mode of the vertical crack on the rail bottom side.In order to verify the correctness of the mode selection model and the influence of the two factors on the detection results, a mode with a larger vibration factor but a smaller orthogonal factor (mode 11 at a frequency of 60 kHz) is selected in Table1, and a mode with a larger orthogonal factor but a smaller vibration factor (mode 7 at a frequency of 60 kHz) is selected in Table 2. Vibration factors, orthogonal factors, and sensitivity data of these three modes are listed in Table 4. The detection results of the three modes on the vertical crack of rail bottom are simulated and analyzed.Table 4
Vibration factors, orthogonal factors, and sensitivity data of modes 1, 7, and 11 at 60 kHz.
Mode no.W1W2W110.960.997010.30110.950.550.83


### 4.2. Simulation and Analysis

Referring to the single-mode excitation algorithm proposed by Xu et al. [15], the excitation direction and location of the selected mode are determined. The rail model with vertical cracks on the rail bottom is established, and the corresponding single mode is excited on the rail according to the solved excitation method. The received reflection waves of the crack and the model end surface are collected, and the crack detection results are analyzed.


#### 4.2.1. Mode Excitation

(1) Determining the Excitation Direction. The excitable area on the rail is limited to the periphery. There are 97 periphery joints of the rail section, shown in Figure 2. The Euclidean distance is used to represent the real distance between modes. The Euclidean distance of any two modes m and n in the longitudinal direction can be expressed as follows: (20)Xρmn=∑i=1pxim−xin2,p=97.There are 35 propagating modes at a frequency of 60 kHz, and the longitudinal Euclidean distance of all modes is(21)Xρ=Xρ11Xρ12⋯⋯Xρ1kXρ21Xρ22⋯⋯Xρ2k⋱⋱Xρk1Xρk2⋯⋯Xρkk,k=35.Similarly, the Euclidean distances of all modes at transverse and vertical directions of the rail can be calculated. The direction with the maximum Euclidean distance of the selected mode is the excitation direction of the mode.The calculated Euclidean distances of modes 1, 7, and 11 inx, y, and z directions are shown in Table 5.Table 5
Euclidean distances of modes at a frequency of 60 kHz.
Mode no.XYZ10.00480.00470.006970.00500.00450.0071110.00510.00790.0046From Table5, it can be seen that modes 1 and 7 have the largest vertical Euclidean distance, and vertical excitation is chosen. Mode 11 has the largest transverse Euclidean distance, and transverse excitation is chosen.(2) Determining the Excitation Point. The torsion degree of rail modes at any position can be expressed by the covariance corresponding to their modes. The covariance of the two nodes t and e of mode m in the transverse, longitudinal, and vertical directions of the rail bottom can be calculated as follows:(22)GD_Covmt,e=∑s=13ts−t¯es−e¯3−1,where ts and es are mode shapes of three directions at nodes t and e, respectively, and t¯ and e¯ are the mean values of the mode shapes in three directions at nodes t and e. Then, the covariance of the q nodes of mode m on the bottom can be expressed by the following equation:(23)GD_Covm=GD_Cov1,1GD_Cov1,2⋯⋯GD_Cov1,qGD_Cov2,1GD_Cov2,2⋯⋯GD_Cov2,q⋱⋱GD_Covq,1GD_Covq,2⋯⋯GD_Covq,q.Similarly, the covariance matrix of the rail head and rail web can be solved separately.If the values of all items of the covariance matrix are positive, the node with the largest autocovariance is selected as the positive excitation point; if there is any negative number in the covariance matrix, the point with the smallest covariance is selected as the reverse excitation point according to the modal vibration characteristics.The excitation points of modes 1, 7, and 11 are solved by equations (22) and (23), and the results are shown in Table 6. The node location is shown in Figure 11.Table 6
Mode excitation points.
Mode no.Excitation node no.129760, 69 (same direction)1130, 3 (reverse direction)Figure 11
Selected node locations.Therefore, at a frequency of 60 kHz, a single mode 1 can be excited by vertically exciting node 29 along the rail, a single mode 7 can be excited by vertically exciting nodes 60 and 69 in the same direction along the rail, and a single mode 11 can be excited by reversely exciting nodes 3 and 30 along the transverse direction of the rail.


#### 4.2.2. Crack Detection

Because the SAFE method is only applicable to the excitation response of an intact waveguide medium, it is impossible to solve the echo signal caused by the crack. Therefore, ANSYS three-dimensional finite element commercial software is used to simulate the cracked rail. First, the rail model is built by using three-dimensional modeling software. The rail length is 3.5 m. The section meshing method is consistent with Figure2, and the width of the meshing element in the length direction is 3 mm. The crack is located at 1.73 m, as shown in Figure 12.Figure 12
Rail model.The excitation signal is the sinusoidal wave modulated by the Hanning window with a frequency of 60 kHz and five cycles. It is received in the same direction at 60 mm, and the receiving point is the node at this section with the same point as excitation. Mode 1 is vertical excitation and receiving node is 29; mode 7 is vertical excitation and nodes 60 and 69 and receiving node 69; and mode 11 is lateral excitation and nodes 3 and 30 and receiving node 30. The total simulation time is 3 ms, with a total of 3600 calculation steps. The simulation results are shown in Figure13.Figure 13
Received signals of (a) mode 1, (b) mode 7, and (c) mode 11 for crack detection.
(a)(b)(c)In Figure13(a), the first wave packet of mode 1 is the direct wave, the second wave packet is the reflection of the rail bottom crack, and the third wave packet is the reflection of the rail end. The distance between the end face and the receiving position is d = 3.44 m, the peak time of direct wave packet t0 is 6.167e − 5 s, the peak time of rail crack reflection packet t1 is 1.157e − 3 s, the amplitude is 0.018 mm, and the peak time of end face reflection packet t2 is 2.308e − 3 s. From equation (24), the group velocity Vg1 = 3062.8 m/s can be calculated, which is only 13.2 m/s different from the theoretical group velocity of 3076.0 m/s. The error is within the allowable range. Therefore, the excited mode is mode 1. According to the received crack reflection signal, the crack position d0 is calculated by equation (25), where Vg0 is the theoretical group velocity. The calculated result is 1.74 m, which is only 0.01 m different from the actual crack location of 1.73 m. Therefore, mode 1 can detect the vertical crack on the rail bottom:(24)Vg1=2×dt2−t0,(25)d0=0.06+t1−t0×Vg02.Similarly, the simulated group velocities of modes 7 and 11 and the detected crack locations are calculated. The results are shown in Table7.Table 7
Simulated group velocity and crack location.
Mode no.Vg0 (m/s)Vg1 (m/s)Group velocity error (m/s)Crack location (m)Location error (m)Amplitude (mm)13076.03062.813.21.740.010.017873148.23119.528.7——0112866.72848.518.21.700.030.0099Note: “—“ means undetected.As shown in Table7, single modes 7 and 11 can be successfully stimulated. In Figure 13(b), the received signal of mode 7 has no crack reflection, which indicates that mode 7 cannot detect the vertical crack on the rail bottom side. In Figure 13(c), the received signal of mode 11 has crack reflection, which can detect the vertical crack of rail bottom with an error of only 0.03 m. However, the amplitude of crack reflection is nearly half of mode 1. Therefore, the detection effect of mode 1 is the best among the three modes, followed by mode 11, with a larger vibration factor, while mode 7, with a larger orthogonal factor but the smallest vibration factor, cannot detect the vertical crack at the bottom of the rail, which is consistent with the results of the model selection. In addition, the simulation results show that the vibration factor has greater influence than the orthogonal factor, so the weight setting of the vibration coefficient should be greater than that of the orthogonal coefficient, that is, α > β.


## 4.1. Selection of Detection Modes

Taking the vertical crack on the rail bottom as an example, the feasibility of the modal selection model is verified. The crack’s location is the same as the line shown in Figure6. The two points A and B coincide with A′ and B′, numbered 29 and 32, and the longitudinal coordinates of the two ends of the crack are zero.The frequency range is [20, 60] kHz. Five frequencies were selected: 20 kHz, 30 kHz, 40 kHz, 50 kHz, and 60 kHz. According to equation (4), all modes at these five frequencies are solved, totaling 122 modes.Equations (9)–(11) are solved to obtain the vibration factors W1 of all modes, which are shown in Figure 8 after normalization.Figure 8
The vibration factorsW1.The factors are sorted from large to small, as shown in Table1. Only the first 20 modes are listed.Table 1
Modal vibration factorW1.
Frequency (Hz)Mode no.W160,0002160,0001150,00020.9950,00010.9940,00020.9940,00010.9830,00020.9830,00010.9720,00020.9720,00010.9660,000110.9550,000110.9320,000120.9240,000100.9030,000100.8660,000310.8650,00050.8530,00040.8530,00070.8440,00040.84From Table1, it can be seen that the vibrations of modes 1 and 2 at [30, 60] kHz are very obviously near 29 and 32 nodes on the rail bottom, and the results are consistent with those in Figure 5. In the modes at a frequency of 20 kHz, there is no obvious vibration mode near 29 and 32 nodes.According to equations (12)–(19), the modal orthogonal factors W2 of all modes are calculated, as shown in Figure 9.Figure 9
Modal orthogonal factorsW2.These factors are sorted from large to small, as shown in Table2. Only the first 20 modes are listed.Table 2
Modal orthogonal factorW2.
Frequency (Hz)Mode no.W240,00070.9999850,00070.9999660,00070.9999660,00030.9995850,00030.9987130,00050.9985340,000150.9978440,000230.9970940,00030.9963560,000140.9959850,000280.9901860,00060.9899550,000140.9892540,000240.9890330,00030.9883850,00060.9831220,00040.9804230,000160.9775360,000340.9756260,000190.97243Table2 shows that the vibration directions of mode 7 at [40, 60] kHz have the best orthogonality with the vertical crack on the rail bottom, but the corresponding modal vibration factors of these three modes rank 120–122, respectively, and the vibration magnitude is the worst. According to Figure 5, mode 7 at a frequency of 60 kHz has only rail head vibration, which is consistent with the calculation results, while modes 1 and 2 at a frequency of 60 kHz have the largest vibration factors and the most obvious vibration. The orthogonal factors of these two modes are 0.9327 and 0.9326, respectively, which rank 43 and 44. The orthogonality of these two modes is good. There is little difference among them and the orthogonality of the mode 7 at a frequency of 40 kHz.Letα = 0.7 and β = 0.3. The crack sensitivity W of each mode is solved according to equation (7). The results are shown in Figure 10.Figure 10
Crack sensitivityW.After sorting from large to small, the first 20 modes are listed in Table3.Table 3
Crack sensitivityW.
Frequency (Hz)Mode no.W60,00020.835460,00010.835250,00020.830350,00010.830040,00020.823640,00010.823130,00020.813930,00010.813220,00020.797520,00010.796150,00050.766530,00070.761930,00040.760940,00040.758360,00050.758140,00050.757950,00040.754860,00040.751760,000120.744140,000210.7363From Table3, it can be seen that modes 1 and 2 at a frequency of 60 kHz have the best sensitivity to the vertical crack on the bottom of the rail, and there is little difference between them. Since mode 1 is symmetrical and mode 2 is antisymmetrical, the excitation method of mode 1 is relatively simple in the process of single-mode excitation, so mode 1 is selected as the detection mode of the vertical crack on the rail bottom side.In order to verify the correctness of the mode selection model and the influence of the two factors on the detection results, a mode with a larger vibration factor but a smaller orthogonal factor (mode 11 at a frequency of 60 kHz) is selected in Table1, and a mode with a larger orthogonal factor but a smaller vibration factor (mode 7 at a frequency of 60 kHz) is selected in Table 2. Vibration factors, orthogonal factors, and sensitivity data of these three modes are listed in Table 4. The detection results of the three modes on the vertical crack of rail bottom are simulated and analyzed.Table 4
Vibration factors, orthogonal factors, and sensitivity data of modes 1, 7, and 11 at 60 kHz.
Mode no.W1W2W110.960.997010.30110.950.550.83


## 4.2. Simulation and Analysis

Referring to the single-mode excitation algorithm proposed by Xu et al. [15], the excitation direction and location of the selected mode are determined. The rail model with vertical cracks on the rail bottom is established, and the corresponding single mode is excited on the rail according to the solved excitation method. The received reflection waves of the crack and the model end surface are collected, and the crack detection results are analyzed.


### 4.2.1. Mode Excitation

(1) Determining the Excitation Direction. The excitable area on the rail is limited to the periphery. There are 97 periphery joints of the rail section, shown in Figure 2. The Euclidean distance is used to represent the real distance between modes. The Euclidean distance of any two modes m and n in the longitudinal direction can be expressed as follows: (20)Xρmn=∑i=1pxim−xin2,p=97.There are 35 propagating modes at a frequency of 60 kHz, and the longitudinal Euclidean distance of all modes is(21)Xρ=Xρ11Xρ12⋯⋯Xρ1kXρ21Xρ22⋯⋯Xρ2k⋱⋱Xρk1Xρk2⋯⋯Xρkk,k=35.Similarly, the Euclidean distances of all modes at transverse and vertical directions of the rail can be calculated. The direction with the maximum Euclidean distance of the selected mode is the excitation direction of the mode.The calculated Euclidean distances of modes 1, 7, and 11 inx, y, and z directions are shown in Table 5.Table 5
Euclidean distances of modes at a frequency of 60 kHz.
Mode no.XYZ10.00480.00470.006970.00500.00450.0071110.00510.00790.0046From Table5, it can be seen that modes 1 and 7 have the largest vertical Euclidean distance, and vertical excitation is chosen. Mode 11 has the largest transverse Euclidean distance, and transverse excitation is chosen.(2) Determining the Excitation Point. The torsion degree of rail modes at any position can be expressed by the covariance corresponding to their modes. The covariance of the two nodes t and e of mode m in the transverse, longitudinal, and vertical directions of the rail bottom can be calculated as follows:(22)GD_Covmt,e=∑s=13ts−t¯es−e¯3−1,where ts and es are mode shapes of three directions at nodes t and e, respectively, and t¯ and e¯ are the mean values of the mode shapes in three directions at nodes t and e. Then, the covariance of the q nodes of mode m on the bottom can be expressed by the following equation:(23)GD_Covm=GD_Cov1,1GD_Cov1,2⋯⋯GD_Cov1,qGD_Cov2,1GD_Cov2,2⋯⋯GD_Cov2,q⋱⋱GD_Covq,1GD_Covq,2⋯⋯GD_Covq,q.Similarly, the covariance matrix of the rail head and rail web can be solved separately.If the values of all items of the covariance matrix are positive, the node with the largest autocovariance is selected as the positive excitation point; if there is any negative number in the covariance matrix, the point with the smallest covariance is selected as the reverse excitation point according to the modal vibration characteristics.The excitation points of modes 1, 7, and 11 are solved by equations (22) and (23), and the results are shown in Table 6. The node location is shown in Figure 11.Table 6
Mode excitation points.
Mode no.Excitation node no.129760, 69 (same direction)1130, 3 (reverse direction)Figure 11
Selected node locations.Therefore, at a frequency of 60 kHz, a single mode 1 can be excited by vertically exciting node 29 along the rail, a single mode 7 can be excited by vertically exciting nodes 60 and 69 in the same direction along the rail, and a single mode 11 can be excited by reversely exciting nodes 3 and 30 along the transverse direction of the rail.


### 4.2.2. Crack Detection

Because the SAFE method is only applicable to the excitation response of an intact waveguide medium, it is impossible to solve the echo signal caused by the crack. Therefore, ANSYS three-dimensional finite element commercial software is used to simulate the cracked rail. First, the rail model is built by using three-dimensional modeling software. The rail length is 3.5 m. The section meshing method is consistent with Figure2, and the width of the meshing element in the length direction is 3 mm. The crack is located at 1.73 m, as shown in Figure 12.Figure 12
Rail model.The excitation signal is the sinusoidal wave modulated by the Hanning window with a frequency of 60 kHz and five cycles. It is received in the same direction at 60 mm, and the receiving point is the node at this section with the same point as excitation. Mode 1 is vertical excitation and receiving node is 29; mode 7 is vertical excitation and nodes 60 and 69 and receiving node 69; and mode 11 is lateral excitation and nodes 3 and 30 and receiving node 30. The total simulation time is 3 ms, with a total of 3600 calculation steps. The simulation results are shown in Figure13.Figure 13
Received signals of (a) mode 1, (b) mode 7, and (c) mode 11 for crack detection.
(a)(b)(c)In Figure13(a), the first wave packet of mode 1 is the direct wave, the second wave packet is the reflection of the rail bottom crack, and the third wave packet is the reflection of the rail end. The distance between the end face and the receiving position is d = 3.44 m, the peak time of direct wave packet t0 is 6.167e − 5 s, the peak time of rail crack reflection packet t1 is 1.157e − 3 s, the amplitude is 0.018 mm, and the peak time of end face reflection packet t2 is 2.308e − 3 s. From equation (24), the group velocity Vg1 = 3062.8 m/s can be calculated, which is only 13.2 m/s different from the theoretical group velocity of 3076.0 m/s. The error is within the allowable range. Therefore, the excited mode is mode 1. According to the received crack reflection signal, the crack position d0 is calculated by equation (25), where Vg0 is the theoretical group velocity. The calculated result is 1.74 m, which is only 0.01 m different from the actual crack location of 1.73 m. Therefore, mode 1 can detect the vertical crack on the rail bottom:(24)Vg1=2×dt2−t0,(25)d0=0.06+t1−t0×Vg02.Similarly, the simulated group velocities of modes 7 and 11 and the detected crack locations are calculated. The results are shown in Table7.Table 7
Simulated group velocity and crack location.
Mode no.Vg0 (m/s)Vg1 (m/s)Group velocity error (m/s)Crack location (m)Location error (m)Amplitude (mm)13076.03062.813.21.740.010.017873148.23119.528.7——0112866.72848.518.21.700.030.0099Note: “—“ means undetected.As shown in Table7, single modes 7 and 11 can be successfully stimulated. In Figure 13(b), the received signal of mode 7 has no crack reflection, which indicates that mode 7 cannot detect the vertical crack on the rail bottom side. In Figure 13(c), the received signal of mode 11 has crack reflection, which can detect the vertical crack of rail bottom with an error of only 0.03 m. However, the amplitude of crack reflection is nearly half of mode 1. Therefore, the detection effect of mode 1 is the best among the three modes, followed by mode 11, with a larger vibration factor, while mode 7, with a larger orthogonal factor but the smallest vibration factor, cannot detect the vertical crack at the bottom of the rail, which is consistent with the results of the model selection. In addition, the simulation results show that the vibration factor has greater influence than the orthogonal factor, so the weight setting of the vibration coefficient should be greater than that of the orthogonal coefficient, that is, α > β.


## 4.2.1. Mode Excitation

(1) Determining the Excitation Direction. The excitable area on the rail is limited to the periphery. There are 97 periphery joints of the rail section, shown in Figure 2. The Euclidean distance is used to represent the real distance between modes. The Euclidean distance of any two modes m and n in the longitudinal direction can be expressed as follows: (20)Xρmn=∑i=1pxim−xin2,p=97.There are 35 propagating modes at a frequency of 60 kHz, and the longitudinal Euclidean distance of all modes is(21)Xρ=Xρ11Xρ12⋯⋯Xρ1kXρ21Xρ22⋯⋯Xρ2k⋱⋱Xρk1Xρk2⋯⋯Xρkk,k=35.Similarly, the Euclidean distances of all modes at transverse and vertical directions of the rail can be calculated. The direction with the maximum Euclidean distance of the selected mode is the excitation direction of the mode.The calculated Euclidean distances of modes 1, 7, and 11 inx, y, and z directions are shown in Table 5.Table 5
Euclidean distances of modes at a frequency of 60 kHz.
Mode no.XYZ10.00480.00470.006970.00500.00450.0071110.00510.00790.0046From Table5, it can be seen that modes 1 and 7 have the largest vertical Euclidean distance, and vertical excitation is chosen. Mode 11 has the largest transverse Euclidean distance, and transverse excitation is chosen.(2) Determining the Excitation Point. The torsion degree of rail modes at any position can be expressed by the covariance corresponding to their modes. The covariance of the two nodes t and e of mode m in the transverse, longitudinal, and vertical directions of the rail bottom can be calculated as follows:(22)GD_Covmt,e=∑s=13ts−t¯es−e¯3−1,where ts and es are mode shapes of three directions at nodes t and e, respectively, and t¯ and e¯ are the mean values of the mode shapes in three directions at nodes t and e. Then, the covariance of the q nodes of mode m on the bottom can be expressed by the following equation:(23)GD_Covm=GD_Cov1,1GD_Cov1,2⋯⋯GD_Cov1,qGD_Cov2,1GD_Cov2,2⋯⋯GD_Cov2,q⋱⋱GD_Covq,1GD_Covq,2⋯⋯GD_Covq,q.Similarly, the covariance matrix of the rail head and rail web can be solved separately.If the values of all items of the covariance matrix are positive, the node with the largest autocovariance is selected as the positive excitation point; if there is any negative number in the covariance matrix, the point with the smallest covariance is selected as the reverse excitation point according to the modal vibration characteristics.The excitation points of modes 1, 7, and 11 are solved by equations (22) and (23), and the results are shown in Table 6. The node location is shown in Figure 11.Table 6
Mode excitation points.
Mode no.Excitation node no.129760, 69 (same direction)1130, 3 (reverse direction)Figure 11
Selected node locations.Therefore, at a frequency of 60 kHz, a single mode 1 can be excited by vertically exciting node 29 along the rail, a single mode 7 can be excited by vertically exciting nodes 60 and 69 in the same direction along the rail, and a single mode 11 can be excited by reversely exciting nodes 3 and 30 along the transverse direction of the rail.


## 4.2.2. Crack Detection

Because the SAFE method is only applicable to the excitation response of an intact waveguide medium, it is impossible to solve the echo signal caused by the crack. Therefore, ANSYS three-dimensional finite element commercial software is used to simulate the cracked rail. First, the rail model is built by using three-dimensional modeling software. The rail length is 3.5 m. The section meshing method is consistent with Figure2, and the width of the meshing element in the length direction is 3 mm. The crack is located at 1.73 m, as shown in Figure 12.Figure 12
Rail model.The excitation signal is the sinusoidal wave modulated by the Hanning window with a frequency of 60 kHz and five cycles. It is received in the same direction at 60 mm, and the receiving point is the node at this section with the same point as excitation. Mode 1 is vertical excitation and receiving node is 29; mode 7 is vertical excitation and nodes 60 and 69 and receiving node 69; and mode 11 is lateral excitation and nodes 3 and 30 and receiving node 30. The total simulation time is 3 ms, with a total of 3600 calculation steps. The simulation results are shown in Figure13.Figure 13
Received signals of (a) mode 1, (b) mode 7, and (c) mode 11 for crack detection.
(a)(b)(c)In Figure13(a), the first wave packet of mode 1 is the direct wave, the second wave packet is the reflection of the rail bottom crack, and the third wave packet is the reflection of the rail end. The distance between the end face and the receiving position is d = 3.44 m, the peak time of direct wave packet t0 is 6.167e − 5 s, the peak time of rail crack reflection packet t1 is 1.157e − 3 s, the amplitude is 0.018 mm, and the peak time of end face reflection packet t2 is 2.308e − 3 s. From equation (24), the group velocity Vg1 = 3062.8 m/s can be calculated, which is only 13.2 m/s different from the theoretical group velocity of 3076.0 m/s. The error is within the allowable range. Therefore, the excited mode is mode 1. According to the received crack reflection signal, the crack position d0 is calculated by equation (25), where Vg0 is the theoretical group velocity. The calculated result is 1.74 m, which is only 0.01 m different from the actual crack location of 1.73 m. Therefore, mode 1 can detect the vertical crack on the rail bottom:(24)Vg1=2×dt2−t0,(25)d0=0.06+t1−t0×Vg02.Similarly, the simulated group velocities of modes 7 and 11 and the detected crack locations are calculated. The results are shown in Table7.Table 7
Simulated group velocity and crack location.
Mode no.Vg0 (m/s)Vg1 (m/s)Group velocity error (m/s)Crack location (m)Location error (m)Amplitude (mm)13076.03062.813.21.740.010.017873148.23119.528.7——0112866.72848.518.21.700.030.0099Note: “—“ means undetected.As shown in Table7, single modes 7 and 11 can be successfully stimulated. In Figure 13(b), the received signal of mode 7 has no crack reflection, which indicates that mode 7 cannot detect the vertical crack on the rail bottom side. In Figure 13(c), the received signal of mode 11 has crack reflection, which can detect the vertical crack of rail bottom with an error of only 0.03 m. However, the amplitude of crack reflection is nearly half of mode 1. Therefore, the detection effect of mode 1 is the best among the three modes, followed by mode 11, with a larger vibration factor, while mode 7, with a larger orthogonal factor but the smallest vibration factor, cannot detect the vertical crack at the bottom of the rail, which is consistent with the results of the model selection. In addition, the simulation results show that the vibration factor has greater influence than the orthogonal factor, so the weight setting of the vibration coefficient should be greater than that of the orthogonal coefficient, that is, α > β.


## 5. Experimental Verification

A vertical crack was artificially cut on the side of the rail bottom, as shown in Figure14. The rail length was 3.5 m, and the crack was located at 1.74 m. Modes 1 and 7 were excited in the rail, respectively, to verify the detection ability of the modes to the crack.Figure 14
Vertical crack at rail bottom.In order to excite mode 1, two transducers with a central frequency of 60 kHz were vertically fastened at node 29 of the rail bottom at one end, and the receiving transducer was fastened at 0.06 m. The signal generator was used to send out the excitation signal, which acts on the rail by an amplifier. The excitation signal was at a frequency of 60 kHz, and the sinusoidal wave was modulated by the Hanning window in five cycles. The device is shown in Figures15(a) and 15(b).Figure 15
Experimental installation: (a, b) mode 1; (c) mode 7.
(a)(b)(c)Ten groups of experiments were conducted. The received signals were collected by an oscilloscope. The results after taking the average value are shown in Figure16(a).Figure 16
Signal acquisition results: (a) mode 1; (b) mode 7.
(a)(b)In Figure16(a), the crack reflection and end reflection signals are obvious. We calculated the time difference between the crack reflection signal and the peak value of the direct wave and extracted the group velocity of mode 1 in Table 7. According to equation (25), the distance between the crack and the excitation transducer is calculated to be 1.80 m, which is only 0.06 m different from the actual distance of 1.74 m. Therefore, the selected mode 1 at a frequency of 60 kHz can detect the vertical crack on the rail bottom side, which is consistent with the simulation results.In order to excite mode 7, two transducers with a center frequency of 60 kHz were vertically pasted at nodes 60 and 69 at one end of the rail, and the receiving transducer was pasted at node 69 at 0.07 m. The signal generator was used to send out the excitation signal, which acts on the rail by an amplifier. The excitation signal was at a frequency of 60 kHz, and the sinusoidal wave was modulated by the Hanning window in five cycles. The device is shown in Figure15(c). Similarly, 10 groups of experimental data were collected and averaged. The results are shown in Figure 16(b). In the waveform, there is only an end reflection wave packet and no crack reflection wave packets. Therefore, mode 7 is not suitable for detecting this crack, which is consistent with the simulation results.


## 6. Discussion

### 6.1. Error Correction

As for the causes of errors, five factors were analyzed. First, the stated center frequency of the transducer is 60 kHz, but the actual center frequency is not exactly the same, which is fastened tightly on the rail by the coupling agent and clamp. The frequency shift is large, and the group velocity of modes at this frequency is a little different from that at a frequency of 60 kHz, which results in a calculation error. Second, changes in temperature and aging of the rail will cause a change in elastic modulus, thus affecting the propagation velocity of guided wave modes, resulting in systematic errors. Third, multimodal aliasing will cause time error, which will lead to detection deviation. Forth, the contact surface between the transducer and the rail is a circle with a diameter of 20 mm, and the measuring position is calculated as its center point, but the real force acting on the rail at every point may not be uniform, which will result in a measurement error. Fifth, the contact between the transducer and the rail is connected by the bonding of the coupling agent. Each time the experiment is carried out, the amount and uniformity of the coupling agent cannot be controlled artificially, which will cause a random error.Through the frequency component analysis of the received signal after mean processing, the actual center frequency of the signal can be obtained, and the second calculation error can be avoided. The system error can be corrected according to many experiments. The magnitude of the third error is related to the mode excitation method and signal processing method. The simpler the mode is, the smaller the error is. The last two uncontrollable errors cannot be corrected.


#### 6.1.1. The Correction of Center Frequency

For the received signal obtained in the experiment, the frequency component distribution is calculated, as shown in Figure17.Figure 17
The frequency component distribution of received signals.According to Figure17, the center frequency of the received signal is 56.5 kHz. Therefore, the theoretical group velocity corresponding to the frequency is used to locate the crack again, and the calculation positioning error is reduced by 0.005 m. Since the group velocity of mode 1 changes little around 60 kHz, the resulting error is relatively small. However, for the mode whose group velocity changes greatly in the detection frequency band, the positioning error caused by the center frequency offset is relatively large.


#### 6.1.2. The Correction of Elastic Modulus

First, the SAFE method was used to solve the group velocity of mode 1 at 56.5 kHz at a different elastic modulus. The calibration curve is shown in Figure18. Then, receiving probe 2 was installed 0.78 m from the experimental receiving point described in Section 5, and 10 groups of experiments were carried out again. The received waveform was processed by signal processing, as shown in Figure 19.Figure 18
Elastic modulus calibration curve (52.73 kHz).Figure 19
Received signals of contrast experiment.According to the time difference of the direct wave between the two signals, the group velocity of guided wave propagation in the rail is 2991.9 m/s, and the corresponding elastic modulus is 204 GPa, that is, the real elastic modulus of the rail. Using this group velocity, the crack location is 1.75 m, and the corrected error is 0.02 m, which is 67% less than the previous one.When mode 1 with a frequency of 60 kHz is used to detect rail crack, the error caused by the change in elastic modulus is much larger than that caused by the deviation of signal center frequency. However, for the modes with serious frequency dispersion, the error caused by the deviation of signal center frequency cannot be ignored.


### 6.2. Optimization and Improvement

In the same frequency range, the group velocity trends of different modes are different. Taking mode 7 and mode 14 with a frequency of 60 kHz as examples, it can be seen from Figure5 that both are modes with obvious rail head vibration. At the same time, in the factors calculation results of transverse crack of rail head, the crack region vibration factors of the two modes are ranked first and sixth, respectively, but the group velocity of the two modes has totally a different trend around 60 kHz. Extract group velocity curves of two modes are shown in Figure 20.Figure 20
Group velocity curves of mode 7 and mode 14.It can be seen from Figure20 that the group velocity curve of mode 7 tends to be horizontal between [50, 70] kHz, and the group velocity value is concentrated around 3227 m/s, which is called as having good nondispersion characteristics. The group velocity curve of mode 14 rises monotonously in the [50, 70] kHz range, and the group velocity value changes from 2415 m/s to 2875 m/s, which is called as poor nondispersion characteristics.The detection signal has a wide frequency band; the frequency band of the received signal is similar to the excitation signal, but the frequency component may change. In the waveform of the received time-domain signal, the frequency component of the main mode will shift. The displacement of the modal frequency components used to detect cracks will lead to the change in group velocity and affect the final positioning accuracy. Therefore, for the crack detection method which depends on time-domain signal or nonsingle-point frequency component, the mode with better nondispersion characteristic should be selected as far as possible.The nondispersion characteristic is considered as the evaluation index in the model. For the five frequency points of 20 kHz, 30 kHz, 40 kHz, 50 kHz, and 60 kHz, the frequency range of 6 kHz near each frequency point is taken to calculate the nondispersion characteristic evaluation index coefficient and frequency stepΔf=10Hz, and the sum of group velocity variation values corresponding to each frequency step is used to describe the nondispersion characteristics of the mode:(26)Fgm=∑i=1NfVgf−3000+Δf×i,m−Vgf−3000+Δf×i−1,mNf.Then, the evaluation index coefficient of modal nondispersion characteristics can be expressed as(27)W3m=lgFgmlgminFg1,Fg2,⋯,Fgm,⋯,Fg122.The evaluation indexes of nondispersion characteristics of 122 modes are calculated, and the results after normalization are shown in Figure21. The nondispersion characteristic coefficients of the first 20 modes sorted from large to small are shown in Table 8.Figure 21
Nondispersion characteristics.Table 8
Nondispersion characteristicsW3.
Frequency (Hz)Mode no.W360000716000080.935000070.725000080.706000090.614000070.5960000100.5760000110.576000030.545000090.5360,00040.5350,000110.5260,00010.5260,00020.5160,000220.5140,00080.5060,00060.5040,00090.5050,00030.4950,00040.47In the first 20 modes of nondispersion, the frequencies of 60 kHz, 50 kHz, and 40 kHz account for 55%, 30%, and 15%, respectively. It can be seen that the higher the frequency is, the better the probability of the nondispersion characteristics is, and the modes 7 and 8 with frequencies of [50, 60] kHz have the best nondispersion characteristics.By introducing the nondispersion evaluation index into the crack sensitivity evaluation model, not only the mode suitable for crack detection can be selected but also the accuracy of crack location detection can be guaranteed.


### 6.3. Other Crack Patterns

Firstly, the rail section is divided into seven areas according to the areas of the crack in the rail, as shown in Figure22. Because the rail is symmetrical, the symmetrical areas are represented by the same label.Figure 22
Rail section division areas.Then, according to the trend of cracks in the rail, the cracks in each area are divided into four categories: transverse cracks, vertical cracks, longitudinal cracks, and diagonal cracks, totally 28 kinds of cracks. The improved crack action mechanism analysis model is applied to select the detection modes of each crack, and the selection results are shown in Table9.Table 9
Detection mode of all cracks.
Crack areaCrack tendFrequency (kHz)Mode numberArea ITransverse607VerticalLongitudinalObliqueArea IITransverseVerticalLongitudinalObliqueArea IIITransverse603VerticalLongitudinalObliqueArea IVTransverseVerticalLongitudinalObliqueArea VTransverse606VerticalLongitudinalObliqueArea VITransverse304VerticalLongitudinal604ObliqueArea VIITransverse601VerticalLongitudinalObliqueIt can be seen from Table9 that all kinds of cracks in the rail head area are suitable for the application of mode 7 with a frequency of 60 kHz; cracks in the rail waist area III and IV are suitable for the application of mode 3 with a frequency of 60 kHz; cracks in the lower part of the rail waist area V are suitable for the application of mode 6 with a frequency of 60 kHz; transverse and vertical cracks in the middle part of the rail bottom area VI are suitable for the application of mode 4 with a frequency of 30 kHz, mode 4 of 60 kHz is suitable for detection of longitudinal and oblique cracks, and mode 1 of 60 kHz is suitable for detection of cracks in foot area VII of rail. It can be seen that the crack detection of the whole rail involves six modes at two frequencies. Therefore, it is the best way to detect cracks in different regions by simultaneous or time-sharing excitation of corresponding modes.


## 6.1. Error Correction

As for the causes of errors, five factors were analyzed. First, the stated center frequency of the transducer is 60 kHz, but the actual center frequency is not exactly the same, which is fastened tightly on the rail by the coupling agent and clamp. The frequency shift is large, and the group velocity of modes at this frequency is a little different from that at a frequency of 60 kHz, which results in a calculation error. Second, changes in temperature and aging of the rail will cause a change in elastic modulus, thus affecting the propagation velocity of guided wave modes, resulting in systematic errors. Third, multimodal aliasing will cause time error, which will lead to detection deviation. Forth, the contact surface between the transducer and the rail is a circle with a diameter of 20 mm, and the measuring position is calculated as its center point, but the real force acting on the rail at every point may not be uniform, which will result in a measurement error. Fifth, the contact between the transducer and the rail is connected by the bonding of the coupling agent. Each time the experiment is carried out, the amount and uniformity of the coupling agent cannot be controlled artificially, which will cause a random error.Through the frequency component analysis of the received signal after mean processing, the actual center frequency of the signal can be obtained, and the second calculation error can be avoided. The system error can be corrected according to many experiments. The magnitude of the third error is related to the mode excitation method and signal processing method. The simpler the mode is, the smaller the error is. The last two uncontrollable errors cannot be corrected.


### 6.1.1. The Correction of Center Frequency

For the received signal obtained in the experiment, the frequency component distribution is calculated, as shown in Figure17.Figure 17
The frequency component distribution of received signals.According to Figure17, the center frequency of the received signal is 56.5 kHz. Therefore, the theoretical group velocity corresponding to the frequency is used to locate the crack again, and the calculation positioning error is reduced by 0.005 m. Since the group velocity of mode 1 changes little around 60 kHz, the resulting error is relatively small. However, for the mode whose group velocity changes greatly in the detection frequency band, the positioning error caused by the center frequency offset is relatively large.


### 6.1.2. The Correction of Elastic Modulus

First, the SAFE method was used to solve the group velocity of mode 1 at 56.5 kHz at a different elastic modulus. The calibration curve is shown in Figure18. Then, receiving probe 2 was installed 0.78 m from the experimental receiving point described in Section 5, and 10 groups of experiments were carried out again. The received waveform was processed by signal processing, as shown in Figure 19.Figure 18
Elastic modulus calibration curve (52.73 kHz).Figure 19
Received signals of contrast experiment.According to the time difference of the direct wave between the two signals, the group velocity of guided wave propagation in the rail is 2991.9 m/s, and the corresponding elastic modulus is 204 GPa, that is, the real elastic modulus of the rail. Using this group velocity, the crack location is 1.75 m, and the corrected error is 0.02 m, which is 67% less than the previous one.When mode 1 with a frequency of 60 kHz is used to detect rail crack, the error caused by the change in elastic modulus is much larger than that caused by the deviation of signal center frequency. However, for the modes with serious frequency dispersion, the error caused by the deviation of signal center frequency cannot be ignored.


## 6.1.1. The Correction of Center Frequency

For the received signal obtained in the experiment, the frequency component distribution is calculated, as shown in Figure17.Figure 17
The frequency component distribution of received signals.According to Figure17, the center frequency of the received signal is 56.5 kHz. Therefore, the theoretical group velocity corresponding to the frequency is used to locate the crack again, and the calculation positioning error is reduced by 0.005 m. Since the group velocity of mode 1 changes little around 60 kHz, the resulting error is relatively small. However, for the mode whose group velocity changes greatly in the detection frequency band, the positioning error caused by the center frequency offset is relatively large.


## 6.1.2. The Correction of Elastic Modulus

First, the SAFE method was used to solve the group velocity of mode 1 at 56.5 kHz at a different elastic modulus. The calibration curve is shown in Figure18. Then, receiving probe 2 was installed 0.78 m from the experimental receiving point described in Section 5, and 10 groups of experiments were carried out again. The received waveform was processed by signal processing, as shown in Figure 19.Figure 18
Elastic modulus calibration curve (52.73 kHz).Figure 19
Received signals of contrast experiment.According to the time difference of the direct wave between the two signals, the group velocity of guided wave propagation in the rail is 2991.9 m/s, and the corresponding elastic modulus is 204 GPa, that is, the real elastic modulus of the rail. Using this group velocity, the crack location is 1.75 m, and the corrected error is 0.02 m, which is 67% less than the previous one.When mode 1 with a frequency of 60 kHz is used to detect rail crack, the error caused by the change in elastic modulus is much larger than that caused by the deviation of signal center frequency. However, for the modes with serious frequency dispersion, the error caused by the deviation of signal center frequency cannot be ignored.


## 6.2. Optimization and Improvement

In the same frequency range, the group velocity trends of different modes are different. Taking mode 7 and mode 14 with a frequency of 60 kHz as examples, it can be seen from Figure5 that both are modes with obvious rail head vibration. At the same time, in the factors calculation results of transverse crack of rail head, the crack region vibration factors of the two modes are ranked first and sixth, respectively, but the group velocity of the two modes has totally a different trend around 60 kHz. Extract group velocity curves of two modes are shown in Figure 20.Figure 20
Group velocity curves of mode 7 and mode 14.It can be seen from Figure20 that the group velocity curve of mode 7 tends to be horizontal between [50, 70] kHz, and the group velocity value is concentrated around 3227 m/s, which is called as having good nondispersion characteristics. The group velocity curve of mode 14 rises monotonously in the [50, 70] kHz range, and the group velocity value changes from 2415 m/s to 2875 m/s, which is called as poor nondispersion characteristics.The detection signal has a wide frequency band; the frequency band of the received signal is similar to the excitation signal, but the frequency component may change. In the waveform of the received time-domain signal, the frequency component of the main mode will shift. The displacement of the modal frequency components used to detect cracks will lead to the change in group velocity and affect the final positioning accuracy. Therefore, for the crack detection method which depends on time-domain signal or nonsingle-point frequency component, the mode with better nondispersion characteristic should be selected as far as possible.The nondispersion characteristic is considered as the evaluation index in the model. For the five frequency points of 20 kHz, 30 kHz, 40 kHz, 50 kHz, and 60 kHz, the frequency range of 6 kHz near each frequency point is taken to calculate the nondispersion characteristic evaluation index coefficient and frequency stepΔf=10Hz, and the sum of group velocity variation values corresponding to each frequency step is used to describe the nondispersion characteristics of the mode:(26)Fgm=∑i=1NfVgf−3000+Δf×i,m−Vgf−3000+Δf×i−1,mNf.Then, the evaluation index coefficient of modal nondispersion characteristics can be expressed as(27)W3m=lgFgmlgminFg1,Fg2,⋯,Fgm,⋯,Fg122.The evaluation indexes of nondispersion characteristics of 122 modes are calculated, and the results after normalization are shown in Figure21. The nondispersion characteristic coefficients of the first 20 modes sorted from large to small are shown in Table 8.Figure 21
Nondispersion characteristics.Table 8
Nondispersion characteristicsW3.
Frequency (Hz)Mode no.W360000716000080.935000070.725000080.706000090.614000070.5960000100.5760000110.576000030.545000090.5360,00040.5350,000110.5260,00010.5260,00020.5160,000220.5140,00080.5060,00060.5040,00090.5050,00030.4950,00040.47In the first 20 modes of nondispersion, the frequencies of 60 kHz, 50 kHz, and 40 kHz account for 55%, 30%, and 15%, respectively. It can be seen that the higher the frequency is, the better the probability of the nondispersion characteristics is, and the modes 7 and 8 with frequencies of [50, 60] kHz have the best nondispersion characteristics.By introducing the nondispersion evaluation index into the crack sensitivity evaluation model, not only the mode suitable for crack detection can be selected but also the accuracy of crack location detection can be guaranteed.


## 6.3. Other Crack Patterns

Firstly, the rail section is divided into seven areas according to the areas of the crack in the rail, as shown in Figure22. Because the rail is symmetrical, the symmetrical areas are represented by the same label.Figure 22
Rail section division areas.Then, according to the trend of cracks in the rail, the cracks in each area are divided into four categories: transverse cracks, vertical cracks, longitudinal cracks, and diagonal cracks, totally 28 kinds of cracks. The improved crack action mechanism analysis model is applied to select the detection modes of each crack, and the selection results are shown in Table9.Table 9
Detection mode of all cracks.
Crack areaCrack tendFrequency (kHz)Mode numberArea ITransverse607VerticalLongitudinalObliqueArea IITransverseVerticalLongitudinalObliqueArea IIITransverse603VerticalLongitudinalObliqueArea IVTransverseVerticalLongitudinalObliqueArea VTransverse606VerticalLongitudinalObliqueArea VITransverse304VerticalLongitudinal604ObliqueArea VIITransverse601VerticalLongitudinalObliqueIt can be seen from Table9 that all kinds of cracks in the rail head area are suitable for the application of mode 7 with a frequency of 60 kHz; cracks in the rail waist area III and IV are suitable for the application of mode 3 with a frequency of 60 kHz; cracks in the lower part of the rail waist area V are suitable for the application of mode 6 with a frequency of 60 kHz; transverse and vertical cracks in the middle part of the rail bottom area VI are suitable for the application of mode 4 with a frequency of 30 kHz, mode 4 of 60 kHz is suitable for detection of longitudinal and oblique cracks, and mode 1 of 60 kHz is suitable for detection of cracks in foot area VII of rail. It can be seen that the crack detection of the whole rail involves six modes at two frequencies. Therefore, it is the best way to detect cracks in different regions by simultaneous or time-sharing excitation of corresponding modes.


## 7. Conclusions

With increased frequency of ultrasonic guided waves, the detection accuracy will also increase, but too high frequency will lead to an increased number of modes, which will also affect the detection range. Moreover, the cross-sectional geometry of rail is complex, so there are many guided wave modes that can be propagated along the rail. Therefore, for cracks in different locations and directions, the first prerequisite is to select a suitable mode and frequency to detect the corresponding cracks. At present, crack detection mostly focuses on cracks in the upper part of the rail and usually chooses the mode by observing the mode shapes at the same frequency and analyzing the guided wave dispersion curves qualitatively. There is no further distinction between the trends of the crack.A guided wave mode selection model is constructed based on the theory that the higher the sensitivity of the mode to cracks, the more efficient the crack detection will be. It can realize the quantitative calculation, sequencing, and selection of optimal modes of detection of cracks in different locations of rails and different directions, and the selection process is scientific and rapid. In the model, there are modal vibration and modal orthogonal factors. Taking the vertical crack on the rail bottom as an example, mode 1 is selected as the detection mode of the crack. The correctness of the model is verified by three-dimensional simulation analysis, and the crack location is directly, accurately, and successfully detected by experiments. At the same time, the reasons of the location errors of crack detection are analyzed and the method of avoiding the error is proposed. For the algorithm of crack location which depends on multifrequency points, the evaluation index of nondispersion characteristic is introduced and the model of crack sensitivity is improved. Finally, according to the location areas and trends, the cracks are divided into 28 types, and the optimal detection mode of each type of crack is calculated, which lays a foundation for further study of mode excitation and crack detection.



---

*Source: 8045626-2020-09-08.xml*