The geometry of the problem under investigation can be described as a hollow circular tube surrounded by a semi-infinite cladding with an inner radius of a and an outer radius of b. This is illustrated in Figure 1 using a cylindrical coordinate system (r, θ, z), where the z-direction extends infinitely, and the circumferential SH wave propagates in the matrix along the θ-direction.

The intrinsic structure relationship for orthotropic isotropic materials in the column coordinate system is: σ r r σ θ θ σ z z σ θ z σ r z σ r θ = C 11 C 12 C 13 0 0 0 C 22 C 23 0 0 0 C 33 0 0 0 C 44 0 0 C 55 0 s y m m e t r y C 66 ε r r ε θ θ ε z z 2 ε θ z 2 ε r z 2 ε r θ (1)

In the small deformation hypothesis, the geometric relationship of the structure in the column coordinate system is: ε r r = ∂ u r ∂ r   ε r θ = 1 2 1 r ⋅ ∂ u r ∂ θ + ∂ u θ ∂ r − u θ r ε θ θ = 1 r ⋅ ∂ u θ ∂ θ + u r r   ε r z = 1 2 ∂ u r ∂ z + ∂ u z ∂ r ε z z = ∂ u z ∂ z   ε θ z = 1 2 ∂ u θ ∂ z + ∂ u z r ⋅ ∂ θ (2)

In Eqs 1, 2, u _i , σ _ij and ε _ij represent the displacement, stress and strain, respectively; C _ij is the elastic constant of the medium.

The equation of motion in the column coordinate system under the condition without body force is: ∂ σ r r ∂ r + 1 r ∂ σ r θ ∂ θ + ∂ σ r z ∂ z + σ r r − σ θ θ r = ρ ∂ 2 u r ∂ t 2 ∂ σ θ r ∂ r + 1 r ∂ σ θ θ ∂ θ + ∂ σ θ z ∂ z + 2 σ r θ r = ρ ∂ 2 u θ ∂ t 2 ∂ σ r z ∂ r + 1 r ∂ σ z θ ∂ θ + ∂ σ z z ∂ z + σ r z r = ρ ∂ 2 u z ∂ t 2 (3)

To transform the basic equations for elastic waves into a state matrix, the displacement and stress field quantities mentioned above must be expressed in vector form. As circumferential SH waves propagate in the circumferential direction, the harmonic factor is defined as e ^i(kθ-ωt) , and the displacements and stresses are reformulated as state vectors in the following manner: u = u r u θ u z T e i k θ − ω t τ i = σ i r σ i θ σ i z T e i k θ − ω t (4)

In Eq. 4, k is the wave number and ω is the angular frequency.

By utilizing the state vector format, the equations of motion can be expressed in terms of the stress and displacement vectors, which are derived as follows: ∂ r τ r ∂ r = − r ρ ω 2 u − ∂ τ θ ∂ θ − r ∂ τ z ∂ z − M τ θ (5)

Where the expression for M is M = 0 − 1 0 1 0 0 0 0 0 (6)

The matrix form of the stress vector can be obtained by combining the intrinsic and geometric equations, expressed as follows: τ r = D 11 ∂ u ∂ r + 1 r D 12 ∂ u ∂ θ + M u + D 13 ∂ u ∂ z (7) τ θ = D 21 ∂ u ∂ r + 1 r D 22 ∂ u ∂ θ + M u + D 23 ∂ u ∂ z (8) τ z = D 31 ∂ u ∂ r + 1 r D 32 ∂ u ∂ θ + M u + D 33 ∂ u ∂ z (9) Where the [ D _ij ] is the matrix of elastic constants, and more information about it can be found in the Appendix.

By combining Eqs 5–9, the dispersion equation for the state vector can be obtained as follows: − k 2 D 22 u + i k r D 12 + D 21 ∂ u ∂ r + i k D 22 M + M D 22 u + r D 11 + D 12 M + M D 21 ∂ u ∂ r + r 2 D 11 ∂ 2 u ∂ r 2 + M D 22 M u + r 2 ρ ω 2 u = 0 (10)

The derivation of Eq. 10 is discussed in detail in the literature (Gao et al., 2019).

In order to solve the above dispersion characteristic equations, the displacement of the circumferential wave is expanded using the Legendre polynomial and Laguerre polynomial, Legendre series expansion for the displacement in the steel pipe and Laguerre series expansion for the displacement in the cement. As the circumferential SH wave only vibrates along the z-direction, only the z-directional displacement component needs to be expanded.

The displacement component in the steel matrix can be expressed as: u s r = ∑ n = 0 N − 1 Ψ n Q n χ (11)

Where Q _n (χ) is the Legendre polynomials of order n over χ ϵ [−1,1], ψ _n is the Legendre polynomials expansion factor, N is the cut-off order of the Legendre polynomials, and s represents the displacement in the steel matrix. Since the interval of the Legendre polynomials is [−1,1], the coordinates r (a ≤ r ≤ b) need to be transformed to the coordinates χ and further deduced as follows: χ = l r − R ,   l = 2 h ,   R = a + b / 2 ,   h = b − a (12) r 2 = h 2 χ + R 2 = h 2 4 χ 2 + R 2 + h R χ r 2 ∂ 2 u s ∂ r 2 = χ 2 ∂ 2 u s ∂ χ 2 + l 2 R 2 ∂ 2 u s ∂ χ 2 + l 2 h R χ ∂ 2 u s ∂ χ 2 r 2 ∂ u s ∂ r = h 2 χ 2 ∂ u s ∂ χ + l R 2 ∂ u s ∂ χ + 2 R χ ∂ u s ∂ χ r ∂ u s ∂ r = χ ∂ u s ∂ χ + l R ∂ u s ∂ χ (13)

Substituting Eqs 11, 13 into Eq. 10, multiplying both sides of the equation by Legendre polynomials of order M and then integrating χ from −1 to 1 yield. − k 2 D 22 s ∑ n = 0 N − 1 ∫ − 1 1 Ψ n Q n χ Q m χ d χ + i k D 12 s + D 21 s χ + l R ∑ n = 0 N − 1 ∫ − 1 1 Ψ n ∂ Q n χ ∂ χ Q m χ d χ + i k D 22 s M + M D 22 s ∑ n = 0 N − 1 ∫ − 1 1 Ψ n Q n χ Q m χ d χ + D 11 s + D 12 s M + M D 21 s χ + l R ∑ n = 0 N − 1 ∫ − 1 1 Ψ n ∂ Q n χ ∂ χ Q m χ d χ + D 11 s χ 2 + l 2 R 2 + l 2 h R χ ∑ n = 0 N − 1 ∫ − 1 1 Ψ n ∂ 2 Q n χ ∂ χ 2 Q m χ d χ + M D 22 s M ∑ n = 0 N − 1 ∫ − 1 1 Ψ n Q n χ Q m χ d χ + ρ ω 2 χ 2 l 2 + R 2 + h R χ ∑ n = 0 N − 1 ∫ − 1 1 Ψ n Q n χ Q m χ d χ = 0 (14)

Using the recurrence relation of the Legendre polynomials, χ ∂ Q n χ ∂ χ , ∂ Q n χ ∂ χ , ∂ 2 Q n χ ∂ χ 2 , χ ∂ 2 Q n χ ∂ χ 2 , χ 2 ∂ 2 Q n χ ∂ χ 2 , χ 2 Q n χ , χ Q n χ can be expressed in analytic form by a series of Legendre polynomials, as described in reference (Zheng et al., 2018). In order to simplify the integral calculation in Eq. 14, the orthogonality of Legendre polynomials can be used for derivation, which can be expressed as: ∫ − 1 1 Q n χ Q m χ d χ = 2 2 m + 1 n = m 0 n ≠ m (15)

The displacement component in the cement cladding can be expressed as: u c r = e − 1 2 χ ∑ n = 0 N − 1 F n L n χ (16)

Where L _n (χ) is the Laguerre polynomials of order n over χ ϵ [0, +∞], F _n is the Laguerre polynomials expansion factor, N is the cut-off order of the Laguerre polynomials and c represents the displacement in the cement cladding. Since the interval of Laguerre polynomials is [0, +∞], the coordinates r (r ≥ b) need to be transformed to the coordinates χ, and further calculations are performed as follows: χ = r − b (17) r 2 = χ + b 2 = χ 2 + b 2 + 2 χ b r 2 ∂ 2 u c ∂ r 2 = χ 2 ∂ 2 u c ∂ χ 2 + b 2 ∂ 2 u c ∂ χ 2 + 2 χ b ∂ 2 u c ∂ χ 2 r 2 ∂ u c ∂ r = χ 2 ∂ u c ∂ χ + b 2 ∂ u c ∂ χ + 2 χ b ∂ u c ∂ χ r ∂ u c ∂ r = χ ∂ u c ∂ χ + b ∂ u c ∂ χ (18)

Substituting Eqs. 16, 18 into equation Eq. 10, multiplying both sides of the equation by Laguerre polynomials of order m, and integrating over χ from 0 to +∞, the result obtained is: − k 2 D 22 c ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ + i k D 12 c + D 21 c χ + b − 1 2 ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ + ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ ∂ L n χ ∂ χ L m χ d χ + i k D 22 c M + M D 22 c ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ + D 11 c + D 12 c M + M D 21 c χ + b − 1 2 ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ + ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ ∂ L n χ ∂ χ L m χ d χ + D 11 c χ + b 2 ( 1 4 ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ − ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ ∂ L n χ ∂ χ L m χ d χ + ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ ∂ 2 L n χ ∂ χ 2 L m χ d χ + M D 22 c M ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ + ρ ω 2 χ + b 2 ∑ n = 0 N − 1 ∫ 0 + ∞ F n e − χ L n χ L m χ d χ = 0 (19)

Using the recurrence relation for Laguerre polynomials, ∂ L n χ ∂ χ , χ ∂ L n χ ∂ χ , χ 2 ∂ L n χ ∂ χ , ∂ 2 L n χ ∂ χ 2 , χ ∂ 2 L n χ ∂ χ 2 , χ 2 ∂ 2 L n χ ∂ χ 2 , χ 2 L n χ , χ L n χ can be expressed in analytic form by a series of Laguerre polynomials, details of which can be found in the Appendix. In order to simplify the integral calculation in Eq. 19, the orthogonality of the Laguerre polynomials is used, which can be expressed as: ∫ 0 + ∞ e − χ L n χ L m χ d χ = 1 n = m 0 n ≠ m (20)

In addition, in order to solve for circumferential SH waves in the cementing structures of cementing casing, the inner surface should satisfy stress-free boundary conditions, and the interface between steel and cement layers should satisfy stress and displacement continuity conditions, the boundary conditions under the state vector as well as the continuity conditions can be expressed as: D 11 s l ∑ n = 0 N − 1 ψ l n − 1 n n n + 1 2 + 1 a D 12 s M ∑ n = 0 N − 1 ψ l n − 1 n + 1 a D 12 s i k ∑ n = 0 N − 1 ψ l n − 1 n = 0 (21) D 11 s l ∑ n = 0 N − 1 ψ l n 1 n n n + 1 2 + 1 b D 12 s M ∑ n = 0 N − 1 ψ l n 1 n + 1 b D 12 s i k I ∑ n = 0 N − 1 ψ l n 1 n − D 11 c ∑ n = 0 N − 1 − n + 0.5 F n − 1 b D 12 c M ∑ n = 0 N − 1 F n − 1 b D 12 c i k ∑ n = 0 N − 1 F n = 0 (22) ∑ n = 0 N − 1 ψ n − ∑ n = 0 N − 1 F n = 0 (23)

Combining Eqs 14, 19, and Eqs 21–23, the linear characteristic equation can be obtained as: k 2 R Ω + i k S Ω + T Ω = 0 (24)

Where R, S, T is the corresponding coefficient matrix, and Ω is the Legendre and Laguerre polynomials coefficient matrix. To simplify the calculation, the unit matrix I and the auxiliary variable Λ = ikRΩ are introduced, transforming Eq. 24 into the linear characteristic problem: S − I − R 0 − i k T 0 0 I Ω Λ = 0 (25)

Ultimately by solving for the eigenvalues of Eq. 25, the wave number k can be obtained, and the dispersion curve of the circumferential SH wave can be calculated.

The steel pipe was calculated with an inner radius a of 80 mm, an outer radius b of 90 mm, and the mechanical property parameters of the steel pipe and cement are presented in Table 1. In order to improve the convergence of higher-order modes in the dispersion curve, the order N of the series cut-off term is set to 10.

Figure 2 shows a comparison between the phase velocity dispersion curves obtained in this research and the results reported in the literature (Paweł, 2017). The latter used a combination of semi-analytical finite elements and fully matched layer techniques to analyze the propagation characteristics of circumferential waves, and the mechanical property parameters and geometric parameters of the materials used in this research are consistent with this method. This demonstrates that the method employed in this research is suitable for evaluating the propagation characteristics of circumferential SH waves in the cementing casing, which verifies the accuracy of the method. In addition, the method uses the recursive and orthogonal properties of the Legendre and Laguerre polynomials to derive analytical solutions to integral expressions, avoiding large-scale integration operations and reducing matrix construction time.

The dispersion curve shows that the velocity of each mode tends to be the same as the frequency increases, with the circumferential SH0 mode being nearly a straight line over the whole frequency range, indicating that the dispersion of this mode is weak, so the circumferential SH0 mode is chosen for the numerical simulation of cement bonding defects in the cementing casing.

A three-dimensional simulation model of the cementing casing was established for finite element analysis with the help of the finite element software ABAQUS, as shown in Figure 3.

The simulation area consists of a steel pipe and cement cladding, which are rigidly connected. The thickness of the cement cladding is considered semi-infinite in the theoretical calculation. However, because of the actual situation and the calculation efficiency, the thickness of the cement cladding is set to 40 mm in the model. The loading area consists of two rectangles of 40 mm in length and 10 mm in width. The excitation signal of the same size and opposite direction is applied to both areas simultaneously along the z-direction only, as shown in Figure 4. The receiving point is set at the inner wall of the casing, in the same plane as the center of the loading area, and the arc length between the two points is 1/4 of the circumference of the inner diameter of the pipe.

This research investigates the effect of the cementitious defect’s radial depth and circumferential angle on the circumferential SH wave, setting the defect in the area between the loading position and the reception point using the part function module in the software. Both defects run the entire length of the cement cladding in the z-direction, so the axial depth is set at 100 mm. Firstly, the circumferential angle of the defect is fixed at 15°, and the radial depth is set in steps of 10 mm, increasing from 10 mm to 40 mm; secondly, the radial depth of the defect is fixed at 10 mm and the circumferential angle is set in steps of 15°, increasing from 15° to 60°. A schematic diagram of the two defects is shown in Figure 5.

Considering the relationship between casing thickness and wavelength, as well as between frequency and detection sensitivity, the modal wavelength of SH0 was selected as 20 mm, twice the wall thickness. Thus, combined with the theoretical derivation of the SH wave dispersion curve, the frequency of the waveguide can be determined as 159 kHz. The excitation signal is a sinusoidal signal modulated by a five-period Hanning window. f t = sin 2 π f t × 0.5 − 0.5 ⁡ cos 2 π f t n 0 , 0 < t ≤ τ , t > τ (26)

Where n is the pulse period; f is the center frequency; τ is the signal pulse width; and τ = n/f. The waveform is shown in Figure 6.

The simulation uses a dynamic display solver, and the solution time step should be less than 1/8 of the period, setting the time step to 0.1 μs and the total solution time to 200 µs In addition, the mesh should be as fine as possible when meshing to ensure that the wave propagation process can be accurately modeled. However, a perfect mesh size will dramatically increase the number of cells and affect computational efficiency. Therefore, it is essential to have a suitable mesh division for the model.

According to the wave transmission effect, the maximum mesh size must satisfy the following conditions. λ min max Δ x , Δ y , Δ z > 8 (27)

In Eq. 27: λ _min is the minimum wavelength corresponding to the excitation of the generated modal wave; Δx, Δy, and Δz are the spacing between two adjacent nodes in the x, y, and z directions, respectively. Taking into account the accuracy of the calculation results and the speed of the calculation, the mesh cell length is chosen as 1 mm and the cell type is the eight-node linear hexahedral cell C3D8R, as shown in Figure 7.

The model is set up in finite element simulation software ABAQUS as described above, and simulations are carried out. As the circumferential SH guide wave only vibrates in the z-direction and propagates in the θ-direction, only the displacement in the z-direction at the receiving point needs to be extracted. To verify the correctness of the proposed model and determine whether circumferential SH waves can be excited, the radial depth and circumferential angle of the bonding defect are first set to 0 to simulate an excellent gluing condition. Then the displacement signal at the receiving point is extracted, and the wave speed is calculated according to the crossing time and compared with the theoretical calculation.

When the casing is effectively cemented, the signal obtained is shown in Figure 8, where the first wave packet represents the initial excitation signal of the loading area, the second wave packet indicates the displacement signal at the receiving point, and the transition time between the peaks of the two wave packets is 38.7 µs The spacing between the two packets is 1/4 of the circumference of the inner diameter of the steel pipe, from which the guided wave speed can be calculated as 3247 m/s. The speed of the dispersion curve in the theoretical calculation is 3194 m/s coincides. This indicates that the model developed above produces a circumferential SH₀ guided wave, and also lays the foundation for subsequent investigation of the interaction between the bonding defects and the circumferential SH wave.

The displacement signal obtained by changing the radial depth when cement bonding defects occur is shown in Figure 9. The four wave packets in the figure basically coincide, and the transition time of the four wave packets does not change significantly compared with the wave packets in intact cementation, which indicates that the change of the radial depth of the defect does not affect the circumferential SH0 wave velocity. In order to explore the interaction between radial depth of the defect and circumferential SH wave amplitude, the wave packet amplitudes at radial depths of 10 mm, 20 mm, 30 mm, and 40 mm are extracted and fitted, respectively. The resulting variation pattern is shown in Figure 10.

Comparing the peak value of the wave packet at the reception point when the cementation is intact, it can be seen that the amplitude of the direct wave is significantly higher. This is because the circumferential SH₀ wave propagates in the circumferential direction of the casing and also diffuses energy into the cement cladding. However, when cement bonding defects exist between the cement and the casing, the circumferential SH₀ wave no longer diffuses energy into the defects. The direct wave energy increases, and the amplitude becomes larger. As the radial depth increases, the change in direct wave amplitude is no longer significant, suggesting that the attenuation of circumferential SH₀ waves by cement bonding defects is independent of the missing volume of cement.

The simulation results for different circumferential angles are shown in Figure 11. Consistent with the simulation results obtained by changing the radial depth, the transition time of the four wave packets does not change, indicating that the wave velocity of the circumferential SH0 wave is not affected by the circumferential angle of the defect. However, it is clear from the local magnification that the amplitude of the direct wave increases as the angle increases. In order to visualize the interaction between the circumferential angle of the defect and the circumferential SH wave, the direct wave amplitudes corresponding to the four angles were fitted, and the results are shown in Figure 12. The attenuation of the circumferential SH waves is related to the missing surface area of the cement.