Singular Integral Solutions of Analytical Surface Wave Model with Internal Crack

: In this study, singular integral solutions were studied to investigate scattering of Rayleigh waves by subsurface cracks. Deﬁning a wave scattering model by objects, such as cracks, still can be quite a challenge. The model’s analytical solution uses ﬁve di ﬀ erent numerical integration methods: (1) the Gauss–Legendre quadrature, (2) the Gauss–Chebyshev quadrature, (3) the Gauss–Jacobi quadrature, (4) the Gauss–Hermite quadrature and (5) the Gauss–Laguerre quadrature. The study also provides an e ﬃ cient dynamic ﬁnite element analysis to demonstrate the viability of the wave scattering model with an optimized model conﬁguration for wave separation. The obtained analytical solutions are veriﬁed with displacement variation curves from the computational simulation by deﬁning the correlation of the results. A novel, veriﬁed model, is proposed to provide variations in the backward and forward scattered surface wave displacements calculated by di ﬀ erent frequencies and geometrical crack parameters. The analytical model can be solved by the Gauss–Legendre quadrature method, which shows the signiﬁcantly correlated displacement variation with the FE simulation result. Ultimately, the reliable analytic model can provide an e ﬃ cient approach to solving the parametric relationship of wave scattering.


Introduction
The propagation of disturbances known as wave propagation in solids has been studied in many branches of physical science and engineering. Wave propagation carries energy (e.g., kinetic and potential energies) through a solid medium and energy transmission over a radiation pattern is through the motion of particles in the wave motion phenomenon. Although wave motion in elastic solids is well-studied, adequately defining a wave scattering model by objects such as cracks still can be quite a challenge.
Over the past decades, crack evaluations using a wave scattering approach in nondestructive evaluations (NDEs) have gained much attention [1,2]. Many researchers in the early 1980s offered analytical solutions to define and explain the scattered surface waves caused by cracks (e.g., subsurface cracks) [3][4][5]. A vertical subsurface crack is considered capable of generating an efficiently scattered field since the crack is one of the most influential damages of surface wave propagation. The analytical model of a vertical subsurface crack [3] is obtained by establishing an integral representation for the scattered field expressed in terms of the fundamental potentials first introduced by Lapwood [6]. Furthermore, numerical methods have been devised to solve the singular integral equations of the analytical model. For example, the Gauss-Chebyshev quadrature rule as a numerical method was adopted by Erdogan and Gupta [7]. Understanding the numerical solution of a singular integration is vital to providing an accurate and reliable solution to the analytical wave scattering model. Such studies Figure 1. Illustration of subsurface crack on the two-dimensional − half-plane. is the distance between the ground and the top of the crack tip, is the distance between the ground and the bottom of the crack tip and and are components of the polar coordinate. The image shows the directions of the incident waves and scattered waves influenced by different and distances.

Derivation of Displacement Potentials.
The displacement field as the final theoretical model can be expressed as two potential functions, a scalar potential ( , ) and a vector potential ( , ) by Helmholtz decomposition [14]: where is the displacement field of scattered waves and is the Nabla or Del operator. Computing the displacement potentials is the main part of defining the theoretical scattering model. The displacement field is defined in the -invariant condition. The waves in the half-plane are zinvariant waves where all functions are dependent on the z-direction. The displacement potentials are defined in the scattered field and satisfy a radiation condition of the elliptic boundary-value problem by Wickham [15] as the equation including the exponential and small ' ' functions: where and are constant; , and are the Rayleigh wavenumbers, longitudinal waves and transverse waves, respectively. and are the components of the polar coordinate ( , ). The ± signs denote forward direction, + and backward direction, −. The part including the polar coordinate in the Equation (2) is ignored since the incident wave angle is 90 degrees. The coefficient and satisfy the equation and are given as: The constant B can be obtained by Equation (3) when the constant can be calculated as an explicit formula given Achenbach's approach [16], which can be expressed as: Where Figure 1. Illustration of subsurface crack on the two-dimensional x − y half-plane. a is the distance between the ground and the top of the crack tip, b is the distance between the ground and the bottom of the crack tip and r and θ are components of the polar coordinate. The image shows the directions of the incident waves and scattered waves influenced by different a and b distances.

Derivation of Displacement Potentials.
The displacement field as the final theoretical model can be expressed as two potential functions, a scalar potential ϕ(x, y) and a vector potential ψ(x, y) by Helmholtz decomposition [14]: where u is the displacement field of scattered waves and ∇ is the Nabla or Del operator. Computing the displacement potentials is the main part of defining the theoretical scattering model. The displacement field is defined in the z-invariant condition. The waves in the half-plane are z-invariant waves where all functions are dependent on the z-direction. The displacement potentials are defined in the scattered field and satisfy a radiation condition of the elliptic boundary-value problem by Wickham [15] as the equation including the exponential and small 'o' functions: where A and B are constant; k R , k L and k T are the Rayleigh wavenumbers, longitudinal waves and transverse waves, respectively. r and θ are the components of the polar coordinate (r, θ). The ± signs denote forward direction, + and backward direction, −. The part including the polar coordinate in the Equation (2) is ignored since the incident wave angle is 90 degrees. The coefficient A and B satisfy the equation and are given as: The constant B can be obtained by Equation (3) when the constant A can be calculated as an explicit formula given Achenbach's approach [16], which can be expressed as: dy ±d y (y)P(y) + id x (y)Q(y) Where The components of integration in Equation (4) are given as: b a dy ±d y (y)P(y) + id x (y)Q(y) (8) where, d y (y) and d x (y) are the displacement functions in y-direction and x-direction, respectively. The displacement functions can be obtained by Hooke's law, d = F/K, where d is a displacement matrix, F is a force matrix and K is a stiffness matrix. The force matrix can be computed by the stress field given by Achenbach [3] as: where µ is Poisson's ratio.
The stiffness matrix is obtained by the equations representing the Hankel function, which is a Bessel function of the third kind.
where H 1 is the first order of the first kind of Hankel function and where ρ is a density and ψ(k + 1) is a digamma function of a complex variable obtained by differentiating the logarithm of a gamma function.

Numerical Integral Solutions of an Analytical Model
The ultimate objective of this study is to calculate and verify displacement variations, which is defined by the ratio of displacements as scattered wave and incident wave occurrences. The displacements are calculated from the potentials derived by the singular integral Equation (4) Notably, numerical integration methods are key to solve the analytical model to provide reliable and accurate displacement variation. Five Gaussian quadrature numerical integration methods are employed: GLEQ, GCQ, GJQ, GHQ and GLAQ.
The key to the Gaussian quadrature method is to define the node (x) and weight (w), which is the additional coefficient known as "weight" at the element and given by the weight function, to calculate integration in numerically as described in Equation (12). The rule of integration is stated as where f (x) is a continuous function, f (x i ) is a dependent variable at a discrete node x i and w i is the weight function at the discrete node x i . To define the node or the reference point in the x coordinate, an orthogonal polynomial is utilized. Depending on the applied orthogonal polynomial, the five numerical integration methods are classified. The Gauss-Legendre quadrature is the simplest integration problem by Legendre polynomials on [−1, 1], and the weight function is: where P n (x) is the Legendre polynomial, which is the solution of the Legendre differential equation The Gauss-Chebyshev quadrature uses Chebyshev polynomials with respect to the weight function, which is given as: The Gauss-Jacobi quadrature is for an approximate integral with the Jacobi polynomials and weight function expressed as: where f is a smooth function on [−1, 1] and α, β > −1. The interval [−1, 1] can be replaced by any other interval by a linear transformation. The Gauss-Hermite quadrature is a form of Gaussian quadrature for approximating the value of integrals of the following kind, +∞ −∞ e −x f (x)dx and the weight function is given as: where H n (x) is a Hermite polynomial, which is the solution to Hermite's differential equation, y − 2xy + 2ny = 0. The Gauss-Laguerre quadrature is an extension of the Gaussian quadrature method for approximating the value of integrals for the following kind +∞ 0 e −x f (x)dx and the weight function is: where L n (x) is the Laguerre polynomial, which gives the solution to Laguerre's equation of a second-order linear differential equation, xy + (1 − x)y + ny = 0. Weight functions of each numerical integration method are described in Table 1. Table 1. Weight function of integral methods.

Computational Simulations Model
The computational simulations, especially the finite element (FE) analysis, is the most widely used method to solve field problems using a numerical approach. In this study, the FE modeling is performed to verify analytical models solved with different numerical integral methods. The model will simulate the wave propagation in the solid media, including the vertical subsurface crack with different frequencies of incident waves. In detail, the cracks are designed by the different 'a' values and a/b ratios. Three values of the parameter 'a' to simulate subsurface crack in those depths are 12 mm, 18 mm and 24 mm cracks, as denoted in Table 2. The selected parameters are typical sizes in the field area. In addition, two a/b ratios with 0.1 and 0.2 are considered to design the internal crack. The increase in the a/b ratio indicates the decrease in crack sizes. These ratios are selected to compare with the existing analytical solution [3]. In addition, the proposed analytical model will be verified with the FE simulation. Each case has a different wavenumber range. A total of 177 FE models are simulated. A total of 150 FE models are designed for 6 different crack geometry cases simulated varying different frequency ranges. For example, there are 27 models for C12-0.1 case for simulating different frequencies described in Table 2. Twenty-seven FE models out of 177 are designed for no crack case.

Model Description
Typically, an FE simulation requires a large-sized simulation model to avoid unexpected waves reflected from the boundary; otherwise, it will show considerable computation time. Several researchers [10] have conducted the FE simulation with damper and energy-absorbing elements (i.e., infinite elements) to reduce the size of the model. The study of FE simulation with damper and infinite element recommends 10 damping zones in an 0.4-0.6/length (distance between wave source and boundary) with 1000 to 2000 of step damping factor increase. Since the simulation model of this study has a shorter distance between the wave source and boundary, 1000 of the initial damping factor, and 5000 of the incremental damping factor are chosen. The FE models in this study compose of three groups: the solid medium group, damper group and wave absorbing infinite element group. The solid medium size is 4000 mm (width) × 1000 mm (height). The 500 mm damper zone is at the boundary of the solid medium group and is represented in the rainbow color of Figure 2. The infinite element is placed at the dampers' end. The applied material properties are assumed for the typical solid material (ρ = 2400 kg/m 2 , E = 35 GPa and υ = 0.2), which is the largest component in the model. Both the analytical model and the FE model have applicability limitations for typical material cases having material nonlinearities that may affect wave propagation. We assume linear elastic wave propagation. The explicit FE simulation does not enforce the structural equilibrium [17], therefore the structural boundary condition is not necessary for the structural equilibrium of the internal structure forces under the external load. The model of the ABAQUS solver is designed using a two-dimensional (2-D) four-node plane strain element (CPE4R) with 2 mm mesh size for the solid medium group and dampers. 1 is the listening node for the backward scattered waves, 2 is the listening node for the forward scattered waves, − 1 and − 1 are the distances between the crack and nodes, 1 and 2, respectively, − 1 is the distance between the wave source and 1. The figure also includes ten dampers in rainbow colors on the three sides of the solid medium and an infinite element in the thin red layer at the end of dampers. The illustration shows the required component and geometry to simulate the wave propagation in the cracked model.
The artificially attenuated dampers are designed by gradually increased alpha massproportional damping factors, as described in Table 3. Typically, the damper zone allows attenuating low-frequency waves, while the infinite element absorbs high-frequency waves [10]. The type of the element is an infinite element, CINPE4.

Discussion of the Wave Separation in FE Simulation
The model based on the wave separation in the P-wave and surface waves is studied in this section. As illustrated in Figure 2, the − 1, 2 represents a distance between the crack location and listening nodes in the x-direction, while − 1 represents the distance between the source point and backward listening node in the x-direction. Typically, each propagated wave (e.g., P-wave, S-wave and surface wave) has a different arrival time due to each wave's different speed. For example, in the simulation with given material property, the speed of the P-wave ( ) and surface wave ( ) are 4084 m/s and 2265 m/s, respectively. However, the model requires a minimum distance between the wave source and the listening nodes to obtain the pure waveforms (e.g., surface wave) so as to avoid overlapped waves. Figure 3 shows four different wave groups excited by an incident sine wave. The two-wave groups (P-wave and surface wave) are incident waves in the time range between 0.2 and 0.8 x 10 −3 s, while the other two-wave groups are the reflected waves at the crack in the time range between 1.2 and 1.6 x 10 −3 s. The three following groups must be considered to obtain a clear P-wave and surface wave separation: i) Reflected P-waves by the crack should not pass when the incident surface waves pass through at 1, as − 1 + − 1 + − 2 ≠ − 1 . ii) the distance, − 1 should be enough to have a wave separation in P-and Surface-waves at 1, (see Region A in Figure 3), − 1 − − 1 > , where T is period of the incident wave. − 1, 2 should be enough to have a wave separation between the reflected P-and Surface-waves by the crack at the 1 and 2, (see region B in Figure 3), − 1, 2 − − 1, 2 > . In addition, the simulations are designed to avoid the S-wave effect to obtain the pure The artificially attenuated dampers are designed by gradually increased alpha mass-proportional damping factors, as described in Table 3. Typically, the damper zone allows attenuating low-frequency waves, while the infinite element absorbs high-frequency waves [10]. The type of the element is an infinite element, CINPE4.

Discussion of the Wave Separation in FE Simulation
The model based on the wave separation in the P-wave and surface waves is studied in this section. As illustrated in Figure 2, the D C−N1, N2 represents a distance between the crack location and listening nodes in the x-direction, while D C−N1 represents the distance between the source point and backward listening node in the x-direction. Typically, each propagated wave (e.g., P-wave, S-wave and surface wave) has a different arrival time due to each wave's different speed. For example, in the simulation with given material property, the speed of the P-wave (C P ) and surface wave (C R ) are 4084 m/s and 2265 m/s, respectively. However, the model requires a minimum distance between the wave source and the listening nodes to obtain the pure waveforms (e.g., surface wave) so as to avoid overlapped waves. Figure 3 shows four different wave groups excited by an incident sine wave. The two-wave groups (P-wave and surface wave) are incident waves in the time range between 0.2 and 0.8 × 10 −3 s, while the other two-wave groups are the reflected waves at the crack in the time range between 1.2 and 1.6 × 10 −3 s. The three following groups must be considered to obtain a clear P-wave and surface wave separation: i) Reflected P-waves by the crack should not pass when the incident surface waves pass through at N1, as ii) the distance, D C−N1 should be enough to have a wave separation in P-and Surface-waves at N1, (see Region A in Figure 3), where T is period of the incident wave. D C−N1, N2 should be enough to have a wave separation between the reflected P-and Surface-waves by the crack at the N1 and N2 , (see region B in Figure 3), In addition, the simulations are designed to avoid the S-wave effect to obtain the surface wave by receiving only the horizontal displacement, which is in a longitudinal direction, thus the S-wave was eliminated.

Results and Discussions
The FE simulation is performed to verify the analytical model introduced in Section 2. The numerically calculated analytical model and FE simulation results are described to the displacement variation of the ratio of the backward scattered wave and incident wave, ⁄ and the ratio of the forward scattered wave and incident wave ⁄ by two parameters of 'a' and ' ', which is related to the frequency ( = 2 / : k is a wavenumber, f is a frequency and C is the wave speed).

Results of the Analytical Model and FE Simulation
The analytical model with five different integral methods is computed using MATLAB (2019, MathWorks, Natick, MA) in accordance with the procedures of Section 2. The FE simulation is designed based on the details described in Section 3. The analytical model and FE simulation results are normalized with the highest magnitude of the variation ratio to compare the shape and tendency of each case. The wavenumbers are considered to have the displacement variation curves in the < 4 range by the given ' ' values (e.g., 12 mm, 18 mm and 24 mm). In the simulation result of the forward scattered case, the pure scattered waves are obtained at a higher frequency range than 10 kHz. Since the lower frequency (e.g., 10 kHz or below) has a longer wavelength, the scattered P-wave overlaps the scattered surface wave, as discussed in Section 3.2. Figure 4 includes two graphs of the backward and forward scattered cases of the C18-0.2 model, which has a 72-mm crack size. In this geometry, the higher amplitude of scattered energy variation shows in the forward scattered direction (see the amplitude of the ⁄ in Figure 4) than in the backward direction (see the amplitude of the ⁄ in Figure 4).

Results and Discussions
The FE simulation is performed to verify the analytical model introduced in Section 2. The numerically calculated analytical model and FE simulation results are described to the displacement variation of the ratio of the backward scattered wave and incident wave, u bs /u in and the ratio of the forward scattered wave and incident wave u f s /u in by two parameters of 'a' and 'k R ', which is related to the frequency (k = 2π f /C: k is a wavenumber, f is a frequency and C is the wave speed).

Results of the Analytical Model and FE Simulation
The analytical model with five different integral methods is computed using MATLAB (2019, MathWorks, Natick, MA, USA) in accordance with the procedures of Section 2. The FE simulation is designed based on the details described in Section 3. The analytical model and FE simulation results are normalized with the highest magnitude of the variation ratio to compare the shape and tendency of each case. The wavenumbers are considered to have the displacement variation curves in the k R a < 4 range by the given 'a' values (e.g., 12 mm, 18 mm and 24 mm). In the simulation result of the forward scattered case, the pure scattered waves are obtained at a higher frequency range than 10 kHz. Since the lower frequency (e.g., 10 kHz or below) has a longer wavelength, the scattered P-wave overlaps the scattered surface wave, as discussed in Section 3.2. Figure 4 includes two graphs of the backward and forward scattered cases of the C18-0.2 model, which has a 72-mm crack size. In this geometry, the higher amplitude of scattered energy variation shows in the forward scattered direction (see the amplitude of the u f s /u in in Figure 4) than in the backward direction (see the amplitude of the u bs /u in in Figure 4).
The displacement variation changes depending on distance between the ground and top crack, denoted as 'a' and the crack size are described in Figure 5. The three models considered are C12-0.2, C18-0.2 and C24-0.2, which have depths of 12 mm, 18 mm and 24 mm, respectively and a crack size of 48 mm, 72 mm and 96 mm. The models have the same a/b ratio, 0.2. Regardless of the crack depth and size, a similar amount of wave scattering occurs in the same a/b ratio.  The displacement variation changes depending on distance between the ground and top crack, denoted as ' ' and the crack size are described in Figure 5. The three models considered are C12-0.2, C18-0.2 and C24-0.2, which have depths of 12 mm, 18 mm and 24 mm, respectively and a crack size of 48 mm, 72 mm and 96 mm. The models have the same / ratio, 0.2. Regardless of the crack depth and size, a similar amount of wave scattering occurs in the same / ratio. The displacement variations based on crack size change is described in Figure 6, including C12-0.1 and C12-0.2 case, which have a 108-mm and 48-mm crack size. The graph result showing 0.1 of the a/b ratio (see Figure 6a) is obtained by a longer crack than the 0.2 case (see Figure 6b) and the model shows a sharper triangular shape of curves in the lower k R a value. The small crack size case has its peak value at the relatively higher k R a. Since the lower frequency waves can pass through the smaller length or size of the crack, it shows lesser wave scattering energy in the lower k R a. Thus, the peak displacement variation is formed in the higher k R a as part of a smaller crack model.  Figure 6a) is obtained by a longer crack than the 0.2 case (see Figure 6b) and the model shows a sharper triangular shape of curves in the lower value. The small crack size case has its peak value at the relatively higher . Since the lower frequency waves can pass through the smaller length or size of the crack, it shows lesser wave scattering energy in the lower . Thus, the peak displacement variation is formed in the higher as part of a smaller crack model.

Representative Correlation Values of Curves
The correlation coefficient is used for quantifying the similarity between two variables. In this study, Pearson's correlation coefficient [18] is adopted to evaluate the similarity of the displacement variation curves. The method of Pearson's correlation coefficient gives the coefficient (r) value between −1 and 1 to define the similarity of the two curves. When r > 0, it means that the two variables are positively correlated and when r < 0, the two variables are negatively correlated. When |r| = 1, it means that the two variables are completely and linearly correlated, that is, they have a functional relationship. When r = 0, it refers to the nonlinear correlation between the two variables. When 0 < |r| < 1, it means there is a certain degree of linear correlation between the two variables, as described in Table 4. And, the closer |r| is to 1, the closer the linear relationship is between two variables. If r is closer to 0, the linear correlation is weaker between two variables, as described in Table 5. Table 4. Condition of correlation coefficient and its relationship between variables.

Condition of Correlation Coefficient
Relationship between Variables r > 0 or r < 0 Positive correlation or negative correlation |r| = 1 Completely linear correlation (functional relationship) r = 0 Nonlinear correlation 0 < |r| <1 Certain degree of linear correlation Table 5. Criteria for a correlation coefficient.

Rank of Correlation Coefficient Meaning
Significance correlation 0.7 ≤ |r| < 1.0 Highly linear correlation The Pearson's correlation coefficient (r) is calculated by Equation (18): where µ X and µ Y is the mean of variable X and Y. In the discussion section, the results of the FE simulation will be compared with the analytical models using five different numerical integration approaches. The calculated correlation coefficients are described in Table 6. The averages of the obtained correlation coefficients from the five integral methods range from 0.79 to 0.94. Therefore, all five integral methods are in the highly linear correlation condition in Table 5. Among the applied integration methods, the Gauss-Legendre quadrature method shows a higher correlation, e.g., 0.94 and obtained an averaged correlation coefficient. The distance of the correlation coefficient between the analytical model and FE simulation can be explained with considerable limitations: 1) the assumption of the analytical model which is not considered in FE simulation and 2) methodological difference between the analytical solutions and FE simulation.
The analytical model is derived based on many assumptions to make the mathematical logic simple. For example, the analytical model uses Green's function expressed in terms of the Hankel functions to obtain the stiffness function (see G I and G II in Equation (10)). Since Green's function adopts a point source effective solution based on the Dirac delta function in mathematical terms, the analytical model limits the non-localized explicit solution. Thus, the results show a difference from the FE simulation results, which is the explicit time variable solution. In addition, the FE simulation results are affected by the distance (D C−N1, N2 and D S−N1 ) and incident waveform; however, the analytical model has no variables pertaining to the horizontal distance and incident waveform.

Conclusions
The presented study mainly focuses on singular integral solutions to investigate the existing analytical model and to provide a new analytical solution. The study also presents an efficient dynamic FE analysis with an optimized model configuration for wave separation. The simulation that uniquely offers the wave scattering phenomena is designed to obtain the displacement value under the same given conditions as with the analytical model. The data from the simulation results are used to demonstrate the viability of the wave scattering model. The obtained analytical solutions with displacement variation curves are investigated by defining the Pearson's correlation of the FE results. Finally, the result shows the significantly correlated displacement variation of the proposed model with the FE results. The obtained results are discussed in terms of the curved shape and the correlation coefficient between the analytical model and FE simulation. The following conclusions are drawn based on the analytical solutions and FE simulation results and discussion presented in the study.

•
The computation time of the analytical model is much shorter than the FE simulation. In the study, 177 FE models are created and simulated in order to provide an adequate resolution of the displacement variation curves. Each simulation computing time is about 20 minutes, and post-processing is required, while the numerical solution for analytical model takes a few seconds without designing and post-processing time.

•
Regardless of crack depth and size, a similar amount of wave scattering energy occurs in the same a/b ratio.

•
The smaller crack sized model has a peak value at the relatively higher k R a. Since the lower frequency waves can pass through the smaller length of the crack, it shows lesser wave scattering in the lower k R a. Therefore, the peak displacement variation is formed in the higher k R a in the case of the smaller crack model.

•
The Gauss-Legendre quadrature integration method shows the highest correlation coefficient values, with an averaged value of 0.94 compared to all considered models, between its analytical model and the FE simulation. • The analytical model can provide an efficient and reliable approach to solving the parametric relationship of wave scattering. Further studies are needed to investigate the limits of applicability of the developed analytical formulation, such as nonlinearity.