A Fast Finite Element Simulation Method of Phased Array Ultrasonic Testing and Its Application in Sleeve Fillet Weld Inspection

Abstract

Numerical simulation can provide quantitative information on ultrasonic beam propagation and plays an important role in analyzing its detection ability and in optimizing the corresponding parameters of the phased array ultrasonic testing (PAUT). In this paper, a fast finite element simulation method of PAUT is developed using an improved explicit integration algorithm and a non-zero element compressed storage method. The new method is applied for the simulation of PAUT of a type-B sleeve weld, and compared with the commercial finite element software and experiment results. We found that the computation time and memory consumption of the new method is only about a 15th and a 40th of the commercial finite element software, respectively.

1. Introduction

Phased array ultrasonic testing (PAUT) technology can flexibly deflect and focus the beam through electronic scanning, so as to realize the rapid imaging detection of the internal defects of materials [1,2]. Compared with the traditional ultrasonic testing technology, it has the advantages of a faster detection speed, higher resolution, and being more suitable for the detection of complex parts [3,4,5,6]. At present, it has been widely used in industrial nondestructive testing fields such as pipeline welds. However, during the implementation of type-B sleeve welding repair and subsequent pipeline operation, some welding defects often appear at the fillet weld [7], which affects the application of type-B sleeve repair technology in the defect repair project of large-diameter and high-strength pipelines. Therefore, the integrity of the type-B sleeve fillet weld structure needs to be tested and evaluated using the phased array ultrasonic testing method [8,9]. As a result of the special welding structure of the type-B sleeve, it is still difficult to detect and quantitatively evaluate the weld defects using phased array ultrasound testing. In order to further study the detection ability of the phased array ultrasonic detection method for internal defects of welding structures such as the type-B sleeve weld, using the numerical simulation method can simulate the propagation characteristics of the phased array ultrasonic in this structure, analyze its detection ability, and optimize its relevant parameters, which can greatly reduce the blindness and workload of the test [10,11].

In the field of numerical simulation of phased array ultrasonic testing, CIVA based on the wave line method is one of the most commonly used simulation softwares [12,13,14]. Based on geometric acoustics, this simplified approximation method has the advantages of having a shorter calculation time and easy implementation [15], while it ignores the fluctuation of sound, resulting in a low accuracy and reliability. Compared with the wave line method, the finite difference time domain method (FDTD) calculated using MATLAB has a higher accuracy [16]. This method obtains the approximate solution of the differential equation by configuring the discrete grid nodes in a continuous region in a certain way, and directly expanding the differential operator with truncated Taylor series on the grid nodes [17]. However, the calculation accuracy of this method is still not high enough, and it is difficult to deal with the problems of the free interface and complex structure [18]. Compared with the wave line method and FDTD, the finite element method (FEM) and the commercial software based on it, such as COMSOL [19] and Abaqus [20,21], have a higher calculation accuracy and adaptability [22,23]. However, in the process of ultrasonic finite element numerical simulation, the element size is usually less than a 10th of the wavelength, resulting in a large amount of calculations and storage, which limits its practical application of the present ultrasonic wave simulation method based on FEM.

In order to solve the above problem, a fast FEM-based simulation method of ultrasonic testing is developed in this paper. Using a new explicit integration algorithm and non-zero element one-dimensional compressed storage of the coefficient matrixes, the amount of calculation and storage of the new fast FEM-based simulation method can be significantly reduced. To verify the calculation efficiency and accuracy of the new method, it is applied for the simulation of PAUT of the type-B sleeve weld and compared with the commercial FEM software and experiment results. These studies are shown in detail in the following three sections: Fast Finite Element Simulation Method for PAUT, Simulation Results and Comparison, and Conclusion.

2. Fast Finite Element Simulation Method for PAUT

2.1. Fast Finite Element Simulation Method

In ultrasonic testing, the probe generates an impulse force on the surface of the specimen under the excitation of pulse voltage, so as to generate an ultrasonic wave in the specimen. According to the basic theory of elasticity, the basic control equation of an ultrasonic wave in solid materials can be obtained as follows [24]:

$$\mu {\nabla }^{2}u+(\lambda +\mu )\nabla (\nabla \cdot u)+f=\rho \ddot{u}+\gamma \dot{u}$$

(1)

where $\lambda$ and $\mu$ are the Lame elastic constant of the material, $\rho$ is the density of the material, $\gamma$ is the ultrasonic damping coefficient, $u$ is the particle displacement vector, and $f$ is the excitation force vector.

According to the basic principle of the finite element method, the above differential control equation can be transformed into the following discrete finite element control equation:

$$[M]\left\{\ddot{U}\right\}+[D]\left\{\dot{U}\right\}+[K]\left\{U\right\}=\left\{F\right\}$$

(2)

where $[M]$ is the mass matrix, $[D]$ is the damping matrix, $[K]$ is the stiffness matrix, $\left\{F\right\}$ is the node load vector, and $\left\{U\right\}$ is the displacement vector of the node to be solved. In the calculation of the wave problem, the consistent (coordinated) mass matrix of an element can be transformed into a diagonally agglomerated mass matrix to obtain more reliable calculation results. It is assumed that the mass of the element is agglomerated at the node of the element; that is, the element mass matrix becomes a diagonal matrix [25]:

$$M_{ij}^{\mathrm{e}}=\{\begin{array}{cc}{\displaystyle \sum _{k=1}^8{M}_{ik}^{\mathrm{e}}}\hfill & (i=j)\hfill \\ 0\hfill & (i\ne j)\hfill \end{array}(i,j=1,\dots ,8)$$

(3)

Equation (2) is actually a discrete dynamic initial value problem. For the general dynamic initial value problem, the time-domain step-by-step integration algorithm can be used to solve it. The current time-domain integration algorithms are roughly divided into two categories, namely the implicit integration method and explicit integration method [26,27]. The implicit integration method is characterized by a high computational stability and unconditional stability, but it needs to solve the coupled linear equations, which require a large amount of storage and calculations. The characteristic of the explicit integration method is that it does not need to solve the coupled linear equations, and can be combined with one-dimensional compressed storage technology. Its storage and calculation can be greatly reduced, but its calculation is conditionally stable. In order to further reduce the amount of storage and calculations in the process of ultrasonic numerical simulation and to meet the need for more than second-order calculation accuracy in engineering applications, a new and more concise explicit integration improved algorithm based on the central difference method [26] and Newmark average velocity method [27] was established as follows:

$$\{\begin{array}{l}{\left\{\dot{U}\right\}}_{t+\mathsf{\Delta }t}={\left\{\dot{U}\right\}}_{t-\mathsf{\Delta }t}-2\mathsf{\Delta }t\left[\tilde{\mathrm{C}}\right]{\left\{\dot{U}\right\}}_t-2\mathsf{\Delta }t\left[\tilde{K}\right]{\left\{U\right\}}_t+2\mathsf{\Delta }t{\left[M\right]}^{-1}{\left\{F\right\}}_t\\ {\left\{U\right\}}_{t+\mathsf{\Delta }t}={\left\{U\right\}}_t+\frac{{\left\{\dot{U}\right\}}_{t+\mathsf{\Delta }t}+{\left\{\dot{U}\right\}}_t}{2}\mathsf{\Delta }t\end{array}$$

(4)

where $\left[\tilde{\mathrm{K}}\right]={\left[M\right]}^{-1}\left[K\right]$, $\left[\tilde{\mathrm{C}}\right]={\left[M\right]}^{-1}\left[C\right]$. It can be seen that solving Equation (4) does not need to solve the coupled linear equations, and the amount of calculations are greatly reduced compared with the traditional implicit integration method. Compared with the existing explicit integration method (Lee’s method) [28], the storage and calculation amounts in the calculation process are also reduced by about half.

All numerical algorithms must have more than second-order accuracy in order to meet the needs of practical engineering applications [29]. Therefore, in order to verify that the calculation accuracy of the improved explicit integration algorithm can meet the actual needs, we used the single free system as an example to perform a preliminary analysis of the accuracy and stability of this method. Firstly, it is assumed that a single degree of freedom system satisfies the following dynamic equation:

$$\ddot{\mathrm{x}}\left(t\right)+2\xi \omega \dot{x}\left(t\right)+{\omega }^{2}x\left(t\right)=f\left(t\right)$$

(5)

According to Equation (4), the following relationship is obtained in the case of a single degree of freedom:

$${\dot{x}}_{t+\mathsf{\Delta }t}=-2\mathsf{\Delta }t{\omega }^2{x}_t-4\mathsf{\Delta }t\xi \omega {\dot{x}}_t+{\dot{x}}_{t-\mathsf{\Delta }t}+2\mathsf{\Delta }t{f}_t+o_1$$

(6)

$$x_{t+\mathsf{\Delta }t}=x_t+\frac{\mathsf{\Delta }t}{2}({\dot{x}}_{t+\mathsf{\Delta }t}+{\dot{x}}_t)+o_2$$

(7)

where $o_1$ and $o_2$ represent the calculation errors of velocity and displacement introduced by the difference schemes, respectively. According to the exact expression of the central difference method and Taylor expansion,

$$\frac{{\dot{x}}_{t+\mathsf{\Delta }t}-{\dot{x}}_{t-\mathsf{\Delta }t}}{2\mathsf{\Delta }t}={\ddot{x}}_t+O\left(\mathsf{\Delta }{t}^3\right)$$

(8)

$$x_{t+\mathsf{\Delta }t}=x_t+\mathsf{\Delta }t{\dot{x}}_t+\frac{\mathsf{\Delta }t}{2}{\ddot{x}}_t+O\left(\mathsf{\Delta }{t}^3\right)$$

(9)

Substitute Equation (8) into Equation (6) to obtain $o_1=\mathrm{O}\left(\mathsf{\Delta }{t}^4\right)$, and substitute Equation (9) into Equation (7) to obtain $o_2=\mathrm{O}\left(\mathsf{\Delta }{t}^3\right)$. Therefore, the calculation cut-off error of Equation (5) is between $\mathrm{O}\left(\mathsf{\Delta }{t}^3\right)$ and $\mathrm{O}\left(\mathsf{\Delta }{t}^4\right)$ in the case of a single degree of freedom. Through the above theoretical analysis, it can be seen that the calculation accuracy of the improved time-domain explicit integration algorithm is between the second-order and the third-order, and may even reach third-order calculation accuracy. It can be seen that, compared with the traditional explicit integration method (Lee’s method), this method not only needs less storage and calculations, but also improves the calculation accuracy, to a certain extent.

Moreover, this paper combines the non-zero element one-dimensional compressed storage method with the new explicit integration algorithm to minimize the amount of calculations and storage, because the coefficient matrixes $\left[\tilde{\mathrm{K}}\right]$ and $\left[\tilde{\mathrm{C}}\right]$ in Equation (4) are sparse matrixes. Most of the elements in the matrixes are zero elements. As shown in Figure 1, only the non-zero elements in the coefficient matrix are stored, and their corresponding row and column numbers are recorded. Compared with the traditional half bandwidth storage method, the storage capacity of this method is further reduced. At the same time, when solving Equation (4), only non-zero elements need to be calculated. Therefore, the number of calculations is also greatly reduced.

2.2. Numerical Model of PAUT of Type-B Sleeve Fillet Weld

In the process of welding repair and the application of the type-B sleeve, the fillet weld has welding defects, which are tested using the phased array ultrasonic testing method. During detection, the phased array ultrasonic probe is placed on the sleeve, and the internal defects in the triangular area of the weld can be detected through phased array ultrasonic sector scanning. According to the basic principle of phased array ultrasonic testing and the phased array ultrasonic testing method of the type-B sleeve fillet weld, a two-dimensional numerical calculation model of the phased array ultrasonic testing of the type-B sleeve fillet weld defects was established, as shown in Figure 2. The phased array ultrasonic excitation source can be approximated as multiple impulse P_i(t) loads arranged according to a certain delay sequence, and the expression of the ith excitation load can be written as follows:

$$P_i(t)=\{\begin{array}{cc}A\cdot \sin [2\pi f\cdot (\mathrm{t}-\mathsf{\Delta }{t}_i)]& 0\le \mathrm{t}-\mathsf{\Delta }{t}_i\le 1/f\\ 0& Others\end{array}$$

(10)

where A is the peak value of the excitation signal amplitude, f is the excitation signal frequency, and $\mathsf{\Delta }{t}_i$ is the delay time of the ith excitation signal. The size array element position x_i, phased array ultrasonic focusing depth d_f, deflection angle θ, and ultrasonic wave velocity v can be written as follows:

$$\mathsf{\Delta }{t}_i=\left(\sqrt{d_f{}^2+x_i{}^2-2\ast d_f\ast x_i\ast \cos (\theta +\pi /2)}-d_f\right)/v$$

(11)

The size of each array element of the phased array ultrasonic probe used in this paper was 0.6 mm, the spacing between each array element was 0.2 mm, the center frequency of the probe was 2 MHz, and the number of excited array elements was 16. The type-B sleeve and fillet weld were made of X70 pipeline steel, with an elastic modulus of 210 GPA, Poisson’s ratio of 0.3, and density of 7.85 g/mm³. The chemical compositions of the X70 pipeline steel are shown in the

Table 1. Chemical compositions of the X70 pipeline steel and the weld.

Element	C	Mn	P	S
Proportion of mass fraction (%)	0.06	1.12	0.01	0.004

. The thickness of the pipe and sleeve was 12.8 mm and 20 mm, respectively, and the gap between them was 2 mm.

3. Simulation Results and Comparison

3.1. Simulation of Phased Array Ultrasonic Sound Field and Detection Signal

In order to study the reliability and efficiency of the fast finite element numerical simulation method, a numerical code based on this method was developed in Fortran Language. This fast FEM code was used to simulate the phased array ultrasonic testing of the internal defects of the type-B sleeve fillet weld, as shown in Figure 3. The material of the type-B sleeve and its weld was X70 stainless steel. The parameters of the simulation model are shown in

Table 2. Parameters of the simulation model.

Parameters	Value
Material density	7.85 kg/m3
Material elastic modulus	$2.1\times {10}^{11}$ N/m2
Material Poisson’s ratio	0.3
Grid size	0.1 mm
Number of probe array elements	16
Array element width	0.6 mm
Space between array elements	0.2 mm
Thickness of pipe	12.8 mm
Thickness of sleeve	20.0 mm
Position of defect relative to probe	(22 mm, 15 mm)
Diameter of defect	1.0 mm

. To validate the reliability and efficiency of the new fast FEM code, commercial finite element software was also used with the same simulation model for comparison.

Figure 4 shows the simulated PA longitudinal wave focusing field with a focusing angle of 56° and a focusing depth of 27 mm at a different given time using the fast FEM code and commercial finite element software. It can be seen that the ultrasonic wave fields simulated by the fast FEM code and commercial finite element software were almost the same as each other. In the ultrasonic wave fields, both the longitudinal wave and the transverse wave could be observed. However, because the velocity of the longitudinal wave was used in the focusing rule of Equations (10) and (11), the longitudinal wave had a much larger amplitude and better focus. The actual focusing angle and depth were in good agreement with the set values. The PA longitudinal wave signal of the defect simulated using the fast FEM code and commercial finite element software is shown in Figure 5. It can also be seen that the waveforms of the two signals were very similar to each other.

Figure 6 shows the simulated PA transverse wave focusing field with a focusing angle of 56° and focusing depth of 27 mm at different given times using the fast FEM code and commercial finite element software. The velocity of the transverse wave was used in the focusing rule of Equations (10) and (11). It was found that the longitudinal wave almost disappeared in the ultrasonic field as the focusing angle of the transverse wave was larger than the critical refraction angle. There are only focused transverse waves in the tested specimen during the transverse wave phased array detection. Moreover, when the transverse wave propagated to the location of the defect, it mainly produced an obvious reflected transverse wave rather than an obvious mode converted longitudinal wave signal. These reduced the interference to the identification of the defect transverse wave echo signal. It can also be seen that the ultrasonic wave fields simulated using the fast FEM code and commercial finite element software were almost the same as each other. The PA transverse wave signal of the defect simulated using the fast FEM code and commercial finite element software is shown in Figure 7. It can also be seen that the waveforms of the two signals were very similar to each other.

In the finite element model, the distance between the defect and the center of the PAUT probe was 26.63 mm. Using the arrival time of the defect signals in Figure 6, the velocity of both the longitudinal wave and transverse wave simulated using the fast FEM method and commercial finite element software could be calculated and compared with the theoretical values, as shown in

Table 3. Comparison of the results of the ultrasonic wave velocity using the different methods.

	Theoretical Value	Commercial Finite Element Software		Fast FEM Code
		Value	Error	Value	Error
Longitudinal wave	6020.18 m/s	5964.17 m/s	0.93%	5875.34 m/s	2.4%
Transverse wave	3217.92 m/s	3196.88 m/s	0.65%	3166.47 m/s	1.6%

. It can be seen that ultrasonic wave velocity from the results of both the fast FEM code and commercial finite element software are in good agreement with the theoretical values.

After verifying the accuracy of the fast finite element numerical simulation method, the rapidity of this method was verified by comparing the storage and calculation time of the fast finite element simulation program and commercial finite element software. It can be seen from Figure 8 that, under the same grid size and time step, the program based on the fast finite element numerical simulation method occupied less CPU and computer memory and had a higher computational efficiency compared with the commercial finite element software.

3.2. Simulation of PAUT Imaging and Experiment Validation

In practical detection, the accuracy of phased array ultrasound imaging is also very important. Therefore, the type-B sleeve fillet weld specimen with defects (in Figure 3) was tested using 16-channel transverse wave phased array ultrasonic testing equipment to obtain the phased array fan scan image of the actual defects (as shown in Figure 9a). According to the fast finite element numerical simulation program established above, the echo signals of the focused sound field at different angles were simulated, and the simulated phased array fan scan imaging was finally generated (as shown in Figure 9b). In the numerical calculation model, the position coordinates of the defect center relative to the center of the rightmost array element were 22 mm and 15 mm, and in the fan scan simulation imaging results, the position coordinates of the defect imaging center were 22.44 mm and 15.13 mm. The absolute error between the defect position of the fan scan imaging and the defect position set by the numerical calculation model was within 0.5 mm in the horizontal and vertical directions, and the relative error was within 2.0% in these two directions. In addition, the fan scan images obtained through the simulation and experiment in Figure 9 had a high consistency.

4. Conclusions

• In this paper, through combining an improved explicit integration algorithm with the non-zero element one-dimensional compressed storage method, a fast finite element numerical simulation method for PAUT was established. Based on this new numerical method, a fast FEM simulation program for PAUT in Fortran Language was developed.
• To verify the reliability and efficiency of this new fast simulation method and program, a simulation of PAUT of the type-B sleeve weld was realized and compared with the commercial finite element software and experiment measurements. Through the simulation results, it was found that both the ultrasonic fields and measured defect signals from the fast FEM simulation method and the commercial finite element software were in good agreement with each other.
• Furthermore, the computation time and memory consumption of the new method was only about a 15th and a 40th of the commercial finite element software, respectively. Because of its high efficiency and low storage space, this method could not only greatly shorten the calculation time, but could also be used for the simulation of phased array ultrasonic testing of large structures, which is difficult to carry out through conventional numerical simulations.
• Finally, the PA fan scan imaging of the defect in the type-B sleeve fillet weld was simulated using this fast FEM simulation method and compared with the experiment results. The reliability of this fast numerical simulation method was verified again. At present, this method has only been compiled into a preliminary program and the operability of this program can be further optimized in the future. In brief, this method has broad application prospects in the field of phased array ultrasonic testing simulations.