# A Study on Integrity Testing for Platform-Pile Systems considering Wave Propagation

Authors: Minghui Fu; Meihong Lin
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Platform-pile systems are typically used in bridge foundations and seaport structures. These systems are difficult to analyze using a conventional low-strain integrity test (LST) due to periodic oscillations. A staggered-grid finite-difference (SFD) method is proposed for a platform-pile system in soil and verified with ABAQUS and measured data. Different impact locations on the top surface of the platform and the lateral surface of the pile are analyzed, and parallel velocities are obtained at different locations. Choosing an appropriate impact location reduces the influence of the periodic oscillation signal; the impact locations at the center of symmetry are the best choices. At the same impact energy, the impact location should not be outside the cross section of the pile when the impact occurs on the top surface of the platform to achieve a meaningful damage signal. When the impact locations are on the lateral surface of the pile, the interferences are relatively small, and the damage vibration signal can be detected. The material and size of the hammer determine the impact duration and should be chosen to reduce the high-frequency influences and obtain reasonable data for analysis.


---

## Body

## 1. Introduction

The need to evaluate the integrity of superstructure pile foundations has significantly increased since an increasing number of wharves and bridges require repair, maintenance, or expansion. As a versatile nondestructive test (NDT) method for the condition assessment of concrete structures such as pile foundations, the impulse-response test has been widely employed in the NDT industry and is recommended for the quantitative condition assessment of concrete structures [1]. The impulse-response test can be used to estimate the flaw size of damaged pile structures [2]. The low-strain integrity test (LST) is an NDT method for assessing the structural integrity of piles in soil. It has the advantages of low cost and straightforward implementation and is a suitable choice for foundation engineering inspections and evaluations of concrete or timber piles. The method has been used in many countries for foundation quality assurance [3]. The optimal distance between the striking point and the receiver in an LST is 0.5 R–0.7 R when the impact location is on the top surface of the pile [4, 5]. The effectiveness of the impulse-response test for the condition assessment of piles with caps has been analyzed using 2D and 3D finite element models without considering the surrounding soil. Typically, the interpretation of the high-frequency wave is very complex due to the reflection from the boundaries of the pile cap and its asymmetry [6]. A triangular wedge was used to provide a longitudinal excitation signal for the integrity assessment of bridge piles and estimate the pile length [7].Experimental studies of the integrity assessment of groups of piles with and without caps were conducted using the impulse-response test, parallel seismic (PS) surveys, and sonic echo tests. The magnitude of the wave reflected by the pile toe could not be determined due to interference by the boundary of the pile cap and the wave energy absorbed by the pile cap, making it difficult to identify the toe signals [8]. The theoretical capabilities of the PS test were investigated to determine the length of the piles and characterize possible defects [9, 10]. The double-velocity symmetrical superposition method (DVSSM) was proposed to prevent high-frequency interference at the pile head without increasing the predominant period of the impact pulse [11]. Three-dimensional characteristics cause high-frequency interference at the head of large solid piles [12]. Thus, it is challenging to determine the optimum excitation position and sensor position. The interference wave can be suppressed when the ratio of the characteristic wavelength to the pile radius is sufficiently large [13].In superstructures, the secondary reflection can be confused with the defect reflection, which reduces the accuracy of defect identification. The asymmetry of the impact load produces bending vibration and bending waves, resulting in high-frequency interference near the pile head. The frequency range of high-frequency interference is affected by the pile radius and modulus [14]. Johnson and Rausche (1996) proposed an algorithm based on two acceleration signals to separate the upward traveling wave signals in the pile direction [15]. The algorithm has been implemented in equipment and software. However, due to the complexity of superstructures and a lack of guidance and experience, it is difficult to determine the optimum excitation point and sensor position to identify the upward traveling wave in site inspections.A coupled finite element and scaled boundary finite element method (FEM) was used to model the pile and the surrounding soil. The results were used as input into a genetically optimized neural network to identify structural defects in concrete piles [16]. The spectral FEM was proposed for the damage detection of structures [17, 18]; this approach can detect cracks and small delaminations. The boundary element method (BEM) can deal with irregular surfaces but is not suitable for inhomogeneous materials [19–21]. The FEM has been used to model wave propagation problems [22, 23]. The elastodynamic finite integration technique (EFIT) considers stress-free boundary conditions and has been used to analyze structures containing flaws [24]. Stress imaging was suggested for P-SV wave propagation using a finite-difference scheme and free surface boundary conditions [25]. The staggered-grid finite-difference (SFD) method was used to model a platform-pile system, considering zero stress on the free surface and absorbing boundary conditions [26]. Zero-stress formulation was used to deal with a planar free surface problem in SFD simulation [27]. H and W formats have been proposed to solve single free surface boundary for 3D fourth-order SFD modeling of seismic motion. Stress imaging is commonly used in seismic wave analysis to simulate the free surface of the earth but is less used for simulating the free surfaces of structures. Compared with other methods, the SFD algorithm has the advantages of high computational efficiency and few computational requirements and has been successfully applied for wave modeling. In addition, the SFD method can deal with composite computational problems.There is a need to determine the optimum excitation positions and sensor locations for integrity testing of platform-pile systems and clarify the formation mechanism of the high-frequency interference generated by the system to eliminate this interference. An SFD numerical solution is developed considering the free surface, heterogeneity, and absorption boundary layers. The rationality and accuracy of the SFD method are verified by comparing the results with ABAQUS simulation results and measurements. The optimized impact locations and sensor locations are provided for integrity testing of platform-pile systems.


## 2. Computational Model and Elastic Wave Equations

### 2.1. Computational Model and Hypothesis

A numerical model of two damaged piles with a cap and a numerical model of a damaged pile in soil were established to understand the influences of the platform and soil on stress wave propagation. A test model of the first model was created for verification. The first model was used to identify the 3D vibration effects of an asymmetric superstructure, and the second model was used to identify the 3D vibration effects of a symmetric superstructure and the surrounding soil. The platform-pile systems were assumed to be a linear-elastic isotropic medium, and body forces were not considered. The integrity coefficientβ in the pile length direction was used for damage establishment. According to the definition of the pile foundation integrity coefficient, β=Z2/Z1, where Zi=ρiCiAi is the cross section impedance, ρi is the density, Ci is the elastic wave velocity, and Ai is the section area. According to the 1D wave theory in the pile length direction,(1)∂2uz,t∂t2=c2∂2uz,t∂z2,where the elastic wave velocity c=E/ρ is used. The relationship between the reflected wave velocity VR (the upward wave) and the incoming wave velocity VI (the downward wave) by the impact hammer is defined as follows:(2)VR=Z1−Z2Z1+Z2VI=1−β1+βVI,where VI and VR are the W velocities in the pile length direction. The material parameters of the pile cross section are altered, and a reflected wave is generated and is monitored by the sensors to analyze the influence of the damage. The velocity can be used for the integrity analysis using the traditional LST method.The computational model is defined using a series of 3D domains (the pile cap, pile body, damage, and soil), and the material parameters are defined separately in each domain. The independent definitions of the material parameters and computational domain ensure that the damaged area is accurately determined. Figure1 shows the computational domains of the two-pile model and the pile-soil model.Figure 1
The computational domains of (a) the two-pile model and (b) the pile-soil model.
(a)(b)


### 2.2. 3D Elastic Wave Equations and Discretization

The computational model of the pile-soil system was based on linear-elastic materials subjected to an LST. According to Hooke’s law and Cauchy’s equation of motion, the 3D elastic wave equations are obtained as follows:(3)ρ∂u∂t=∂σxx∂x+∂τxy∂y+∂τxz∂z+fx,(4)ρ∂v∂t=∂τxy∂x+∂σyy∂y+∂τyz∂z+fy,(5)ρ∂w∂t=∂σxz∂x+∂τyz∂y+∂σzz∂z+fz,(6)∂σxx∂t=λ+2μ∂u∂x+λ∂v∂y+λ∂w∂z,(7)∂σyy∂t=λ∂u∂x+λ+2μ∂v∂y+λ∂w∂z,(8)∂σzz∂t=λ∂u∂x+λ∂v∂y+λ+2μ∂w∂z,(9)∂τxy∂t=μ∂u∂y+∂v∂x,(10)∂τxz∂t=μ∂u∂z+∂w∂x,(11)∂τyz∂t=μ∂v∂z+∂w∂y.


### 2.3. Initial Conditions

AZ-direction impact force p(t) is applied to the impact surface as follows:(12)pt=It01−cos2πt0t,t≤t0,0,t>t0,where t0 is the contact time and I is the unit impulse of the force. The τzz values in the impact area are determined by the pressure P=pt/πr2, and r is the radius of the impact surface.


### 2.4. Free Surface Boundary Condition

Zero stress must be satisfied on the free surface boundaries.(13)σxx=τxy=τxz=0,σyy=τxy=τyz=0,σzz=τxz=τyz=0.


### 2.5. Perfectly Matched Layer Technique

The velocity-stress of a perfectly matched layer (PML) is used to eliminate the reflection of waves at the boundaries [28], and the maximum damping value d0 is calculated as follows:(14)d0=log1RτVsnbΔh,where R is the reflection coefficient and Vs is the shear velocity.


### 2.6. Convergence Condition

The Courant–Friedrichs–Lewy criterion must be satisfied to obtain numerical stability.(15)Δt≤1c1/Δx2+1/Δy2+1/Δz2,where c is the fastest longitudinal wave speed. When an unequal step is used, Δz=aΔx=aΔy=aΔs, and the convergence condition is as follows:(16)Δt≤ac2a2+1Δs,where Δs must meet the requirement of approximately 10 grid points per shear wavelength [24].


## 2.1. Computational Model and Hypothesis

A numerical model of two damaged piles with a cap and a numerical model of a damaged pile in soil were established to understand the influences of the platform and soil on stress wave propagation. A test model of the first model was created for verification. The first model was used to identify the 3D vibration effects of an asymmetric superstructure, and the second model was used to identify the 3D vibration effects of a symmetric superstructure and the surrounding soil. The platform-pile systems were assumed to be a linear-elastic isotropic medium, and body forces were not considered. The integrity coefficientβ in the pile length direction was used for damage establishment. According to the definition of the pile foundation integrity coefficient, β=Z2/Z1, where Zi=ρiCiAi is the cross section impedance, ρi is the density, Ci is the elastic wave velocity, and Ai is the section area. According to the 1D wave theory in the pile length direction,(1)∂2uz,t∂t2=c2∂2uz,t∂z2,where the elastic wave velocity c=E/ρ is used. The relationship between the reflected wave velocity VR (the upward wave) and the incoming wave velocity VI (the downward wave) by the impact hammer is defined as follows:(2)VR=Z1−Z2Z1+Z2VI=1−β1+βVI,where VI and VR are the W velocities in the pile length direction. The material parameters of the pile cross section are altered, and a reflected wave is generated and is monitored by the sensors to analyze the influence of the damage. The velocity can be used for the integrity analysis using the traditional LST method.The computational model is defined using a series of 3D domains (the pile cap, pile body, damage, and soil), and the material parameters are defined separately in each domain. The independent definitions of the material parameters and computational domain ensure that the damaged area is accurately determined. Figure1 shows the computational domains of the two-pile model and the pile-soil model.Figure 1
The computational domains of (a) the two-pile model and (b) the pile-soil model.
(a)(b)


## 2.2. 3D Elastic Wave Equations and Discretization

The computational model of the pile-soil system was based on linear-elastic materials subjected to an LST. According to Hooke’s law and Cauchy’s equation of motion, the 3D elastic wave equations are obtained as follows:(3)ρ∂u∂t=∂σxx∂x+∂τxy∂y+∂τxz∂z+fx,(4)ρ∂v∂t=∂τxy∂x+∂σyy∂y+∂τyz∂z+fy,(5)ρ∂w∂t=∂σxz∂x+∂τyz∂y+∂σzz∂z+fz,(6)∂σxx∂t=λ+2μ∂u∂x+λ∂v∂y+λ∂w∂z,(7)∂σyy∂t=λ∂u∂x+λ+2μ∂v∂y+λ∂w∂z,(8)∂σzz∂t=λ∂u∂x+λ∂v∂y+λ+2μ∂w∂z,(9)∂τxy∂t=μ∂u∂y+∂v∂x,(10)∂τxz∂t=μ∂u∂z+∂w∂x,(11)∂τyz∂t=μ∂v∂z+∂w∂y.


## 2.3. Initial Conditions

AZ-direction impact force p(t) is applied to the impact surface as follows:(12)pt=It01−cos2πt0t,t≤t0,0,t>t0,where t0 is the contact time and I is the unit impulse of the force. The τzz values in the impact area are determined by the pressure P=pt/πr2, and r is the radius of the impact surface.


## 2.4. Free Surface Boundary Condition

Zero stress must be satisfied on the free surface boundaries.(13)σxx=τxy=τxz=0,σyy=τxy=τyz=0,σzz=τxz=τyz=0.


## 2.5. Perfectly Matched Layer Technique

The velocity-stress of a perfectly matched layer (PML) is used to eliminate the reflection of waves at the boundaries [28], and the maximum damping value d0 is calculated as follows:(14)d0=log1RτVsnbΔh,where R is the reflection coefficient and Vs is the shear velocity.


## 2.6. Convergence Condition

The Courant–Friedrichs–Lewy criterion must be satisfied to obtain numerical stability.(15)Δt≤1c1/Δx2+1/Δy2+1/Δz2,where c is the fastest longitudinal wave speed. When an unequal step is used, Δz=aΔx=aΔy=aΔs, and the convergence condition is as follows:(16)Δt≤ac2a2+1Δs,where Δs must meet the requirement of approximately 10 grid points per shear wavelength [24].


## 3. Discretization Schemes

### 3.1. Discretization Architecture

The stress image technique is typically used in seismic wave analysis where only a single free surface exists. Since the velocity components and the stress components of the staggered grid are not saved in the same location, the discretization of the free surface becomes more complex, and a half-grid deviation from the original free surface position will occur during the computation [29]. Due to the existence of faces, edges, and corner points, few studies focused on the implementation of the stress image technique in composite structures. The velocity components are located on the physical surface of the computational model, the shear stress terms are also on the surface so that they can be set to zero, and guided elastic waves in structures are analyzed at MHz frequencies [24]; thus, these free surface treatments limit the SFD modeling and generalization ability. The SFD method has been used for the analysis of platform-pile systems without considering the free surface of 3D structures [26, 30, 31]. The computational domain is a regular cuboid domain. Therefore, it is necessary to extend the SFD method for use in platform-pile systems with free surfaces.An SFD for the platform-pile system of composite structures is proposed to provide a systematic and widely applicable solution for the structures. We extend the stress image technique [32] to 3D structures and use multidirectional imaging on the free surfaces. The SFD system for the free surface treatment (720 grid types), heterogeneity, and absorption boundary layers (72 grid types) is implemented using group batch processing, and the calculation of the SFD model is straightforward. It must be pointed out that the overlap of the image grid points of two very close free surfaces (including faces, edges, and corner points) is handled by hash map implementation of the grid point indices to improve the applicability of the SFD scheme. The SFD system is highly adaptable to complex composite structures and is suitable for the simulation of a platform-pile system with a small impact block.


### 3.2. Discretization of Elastic Wave Equations

The (2M)-th-order finite-difference scheme of the first-order spatial derivative of the function ux can be expressed as follows [33]:(17)∂u∂x=1h∑r=1Marux+rh−0.5h−ux−rh+0.5h,where h is the grid size, which can be replaced by Δx, Δy, and Δz, and arr=1,2,…,M are the SFD coefficients of the spatial derivatives.The configuration of the grid points for the finite-difference discretization is shown in Figure2. All stress and velocity components are located in the small cell. S represents the σxx, σyy, and σzz grid points. Txy, Tyz, and Txz, respectively, represent τxy, τyz, and τxz grid points. U, V, and W are the wave velocities in the X-, Y-, and Z-directions.Figure 2
The configuration of the stress and velocity grid points.The details of the discretization equations are as follows:(18)um+1/2=um−1/2+ΔtρDxτxxm+Dyτxym+Dzτxzm,vm+1/2=vm−1/2+ΔtρDxτxym+Dyτyym+Dzτyzm,wm+1/2=wm−1/2+ΔtρDxτxzm+Dyτyzm+Dzτzzm,τxxm+1=τxxm+Δtλ+2μDxum+1/2+λDyvm+1/2+Dzwm+1/2,τyym+1=τyym+ΔtλDxum+1/2+λ+2μDyvm+1/2+Dzwm+1/2,τzzm+1=τzzm+ΔtλDxum+1/2+λDyvm+1/2+λ+2μDzwm+1/2,τxym+1=τxym+ΔtμDxvm+1/2+Dyum+1/2,τyzm+1=τyzm+ΔtμDywm+1/2+Dzvm+1/2,τxzm+1=τxzm+ΔtμDxwm+1/2+Dzum+1/2,where Dx, Dy, and Dz are the x, y, and z derivatives that can be calculated by (17).


### 3.3. Discretization of Free Surface Boundaries

Few studies analyzed free surface treatments of 3D structures using the stress imaging method. Thus, the SFD for the free surface treatment of composite structures is proposed. The types of grid points are analyzed in detail and are divided into 9 cases to clarify the discretization types of the free surface boundaries. Finally, the SFD system can be solved by discretizing all grid points according to their type. It must be pointed out that the complex problem of free surface boundaries is simplified, resulting in a straightforward implementation.The free surfaces are defined as 7 types to identify the types of boundaries: faces (X, Y, Z), edges (XY, YZ, XZ), and eight corner points (XYZ), as exhibited in Figure 3. Moreover, we use one-directional imaging on the boundary faces (Figure 4), two-directional imaging on the boundary edges (Figure 5), and three-directional imaging on the boundary corner points (Figure 6). The types of image grid points for the boundary faces, edges, and corner points are defined as 1 (one-directional imaging) × 3 (faces X, Y, Z) + 2 (two-directional imaging) × 3 (edges XY, YZ, XZ) + 3(three-directional imaging) × 1 (corner points XYZ) = 12 types.Figure 3
The definitions of the free surface boundary.Figure 4
Free surface of boundariesX, Y, and Z. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Figure 5
Free surface of boundariesXY, XZ, and YZ. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Figure 6
Free surface of boundaryXYZ. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Further, three boundary layers outside the free surface are required to generate the free surface treatment, namely, the first, second, and third layers. Leth be the mesh step corresponding to the X-, Y-, and Z-directions; h can be replaced with Δx, Δy, and Δz. Let ±h/2,±h, and ±3h/2 be the distances between the first, second, and third boundary layers and the free surface; there are 6 types of boundary layers.Finally, there are 3 (U, V, W) + 6 (Txx, Txy, Txz, Tyy, Tyz, Tzz) =  9 components. Thus, there are (1 + 7) × 9 = 72 grid types in the internal area and on the free surfaces. Further, the number of image grid types is 6 (the type count of boundary layers) × 12 (the types of image points) × 9 = 648. Thus, a total of 720 grid types must be considered for the free surface treatment in the SFD calculation.Note. (+) and (−) values indicate the normal face points in the positive and negative directions of the related axis.The stress componentsσij of the boundary grid points outside the free surface are computed by the interior image components [32]. The details of the SFD discretization of the boundary velocity components and the stress components on the free surface are listed as follows.


#### 3.3.1.X, Y, Z Boundaries

We use boundaryZ as an example; there are two types of boundary grid configurations on the free surface, as shown in Figure 4. The first type is the S grid points on the free surface Z, and the second type is the velocity grid points on the free surface. The discretization details are described below.Case 1.
Stress grid points on the free surface.
We useTxx on the boundary Z as an example, as shown in Figure 4(a). When the S grid points are on the free surface Z, the stress component Tzz=0, and Txx (τxx) is computed by (6) and (8).(19)Txxm+1=Txxm+AΔtDU,x+BΔtDV,y,where A=4λμ+4μ2/λ+2μ and B=2λμ/λ+2μ. DU,x is the x derivative of U according to (17). Equation (27) can be expressed as Ti,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,y; similarly, Tyy (τyy) and Tzz (τzz) of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 1.Table 1
Summary of stress components on the free surface.
BoundaryComponentFormulaXTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DW,zXTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DV,yYTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DW,zYTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DU,xZTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,yZTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DU,xNote.i, j, k are the x, y, z indices of grid points on the free surface; m is the time step.Case 2.
Velocity grid points at the first boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
We useW−Δz/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(a). W−Δz/2 is calculated according to the second-order finite-difference approximation of the condition Tzz,t0=0 (Tzz,t is the time derivative of Tzz) and Hooke’s law for Tzz:(20)Wi,j,k−1/2m+1/2=Wi,j,k+1/2m+1/2+Δzλλ+2μDU,x+Δzλλ+2μDV,y.
Equation (20) can be expressed as Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 2.Table 2
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f2Ui+1/2,j,km+1/2,−Δx,DV,y,DW,zX+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f2Ui−1/2,j,km+1/2,Δx,DV,y,DW,zY−Vi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f2Vi,j+1/2,km+1/2,−Δy,DU,x,DW,zY+Vi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f2Vi,j−1/2,km+1/2,Δy,DU,x,DW,zZ−Wi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,yZ+Wi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f2Wi,j,k−1/2m+1/2,Δz,DU,x,DV,yCase 3.
Velocity grid points at the second layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
We useU−h at the first layer outside the boundary Z as an example, as shown in Figure 4(a). U−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzx−h/2=−Tzxh/2 and Hooke’s law for Tzx:(21)U−h=U0+hwx−h/2−hTxz,t−h/2/μ,where wx−h/2=DxW−h/2. Equation (21) can be expressed as Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xz; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 3.Table 3
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Vi−1,j,km+1/2Vi−1,j,km+1/2=f3Vi,j,km+1/2,−Δx,Di−1/2,j,kU,y,Ti−1/2,j,kxyX−Wi−1,j,km+1/2Wi−1,j,km+1/2=f3Wi,j,km+1/2,−Δx,Di−1/2,j,kU,z,Ti−1/2,j,kxzX+Vi+1,j,km+1/2Vi+1,j,km+1/2=f3Vi,j,km+1/2,Δx,Di+1/2,j,kU,y,Ti+1/2,j,kxyX+Wi+1,j,km+1/2Wi+1,j,km+1/2=f3Wi,j,km+1/2,Δx,Di+1/2,j,kU,z,Ti+1/2,j,kxzY−Ui,j−1,km+1/2Ui,j−1,km+1/2=f3Ui,j,km+1/2,−Δy,Di,j−1/2,kV,x,Ti,j−1/2,kxyY−Wi,j−1,km+1/2Wi,j−1,km+1/2=f3Vi,j,km+1/2,−Δy,Di,j−1/2,kV,z,Ti,j−1/2,kyzY+Ui,j+1,km+1/2Ui,j+1,km+1/2=f3Ui,j,km+1/2,Δy,Di,j+1/2,kV,x,Ti,j+1/2,kxyY+Wi,j+1,km+1/2Wi,j+1,km+1/2=f3Vi,j,km+1/2,Δy,Di,j+1/2,kV,z,Ti,j+1/2,kyzZ−Ui,j,k−1m+1/2Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xzZ−Vi,j,k−1m+1/2Vi,j,k−1m+1/2=f3Vi,j,km+1/2,−Δz,Di,j,k−1/2W,y,Ti,j,k−1/2yzZ+Ui,j,k+1m+1/2Ui,j,k+1m+1/2=f3Ui,j,km+1/2,Δz,Di,j,k+1/2W,x,Ti,j,k+1/2xzZ+Vi,j,k+1m+1/2Vi,j,k+1m+1/2=f3Vi,j,km+1/2,Δz,Di,j,k+1/2W,y,Ti,j,k+1/2yzCase 4.
Velocity grid points at the second boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use W−h as the second layer outside the boundary Z as an example, as shown in Figure 4(d). W−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2 and Hooke’s law for Tzz:(22)W−h=W0+−Δzλ+2μ×Tzz,t−h/2−λux−h/2−λvy−h/2.
Equation (22) can be expressed as Wi,j,k−1m+1/2=f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 4.Table 4
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Ui−1,j,km+1/2f4Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,t,Di−1/2,j,kV,y,Di−1/2,j,kW,zX+Ui+1,j,km+1/2f4Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,t,Di+1/2,j,kV,y,Di+1/2,j,kW,zY−Vi,j−1,km+1/2f4Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,t,Di,j−1/2,kU,x,Di,j−1/2,kW,zY+Vi,j+1,km+1/2f4Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,t,Di,j+1/2,kU,x,Di,j+1/2,kW,zZ−Wi,j,k−1m+1/2f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,yZ+Wi,j,k+1m+1/2f4Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,t,Di,j,k+1/2U,x,Di,j,k+1/2V,yCase 5.
Velocity grid points at the first layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use U−h/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(d). U−h/2 is calculated according to the second-order finite-difference approximation of the condition Tzx0=0 and Hooke’s law for Tzx. Moreover, V−h/2 is calculated in the same way.(23)U−h/2=Uh/2+h×wx0,where wx0=DxW0. Equation (23) can be expressed as Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,x; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 5.Table 5
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Vi−1/2,j,km+1/2Vi−1/2,j,km+1/2=f5Vi+1/2,j,km+1/2,−Δx,Di,j,kU,yX−Wi−1/2,j,km+1/2Wi−1/2,j,km+1/2=f5Wi+1/2,j,km+1/2,−Δx,Di,j,kU,zX+Vi+1/2,j,km+1/2Vi+1/2,j,km+1/2=f5Vi,j,k−1/2m+1/2,Δx,Di,j,kU,yX+Wi+1/2,j,km+1/2Wi+1/2,j,km+1/2=f5Wi,j,k−1/2m+1/2,Δx,Di,j,kU,zY−Ui,j−1/2,km+1/2Ui,j−1/2,km+1/2=f5Ui,j+1/2,km+1/2,−Δy,Di,j,kV,xY−Wi,j−1/2,km+1/2Wi,j−1/2,km+1/2=f5Wi,j+1/2,km+1/2,−Δy,Di,j,kV,zY+Ui,j+1/2,km+1/2Ui,j+1/2,km+1/2=f5Ui,j+1/2,km+1/2,Δy,Di,j,kV,xY+Wi,j+1/2,km+1/2Wi,j+1/2,km+1/2=f5Wi,j+1/2,km+1/2,Δy,Di,j,kV,zZ−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,xZ−Vi,j,k−1/2m+1/2Vi,j,k−1/2m+1/2=f5Vi,j,k+1/2m+1/2,−Δz,Di,j,kW,yZ+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f5Ui,j,k−1/2m+1/2,Δz,Di,j,kW,xZ+Vi,j,k+1/2m+1/2Vi,j,k+1/2m+1/2=f5Vi,j,k−1/2m+1/2,Δz,Di,j,kW,y


#### 3.3.2.XY, YZ, XZ Boundaries

The grid point arrangements of boundariesXY, XZ, and YZ are shown in Figure 5. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 6.
Velocity grid points at the first boundary layer.
When theS grid points are on the free surface boundary Z, as shown in Figure 5(a), U−h/2 at the first layer of boundary XY is obtained from the second-order finite-difference approximation of the condition Txx0=0, Tyy0=0 and Hooke’s law for Txx, Tyy.(24)U−h/2=Ui−1/2,j,km+1/2=Ui+1/2,j,km+1/2−−Δxλ2λ+μwz.
Equation (24) can be expressed as Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,z; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 6.Table 6
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaXYUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,zXYUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kW,zXYVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kW,zXYVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kW,zYZVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kU,xYZVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kU,xYZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kU,xYZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kU,xXZUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kV,yXZUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kV,yXZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kV,yXZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kV,yCase 7.
Velocity grid points at the second layer.
When theW grid points are on the free surface boundary Z, as shown in Figure 5(d), W−h at the second layer of boundary XZ is calculated by (25) according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2=−Tzzh/2 and Hooke’s law for Tzz.(25)W−h=W0+−Δz4μ2+4λμ×λ+2μTzz,t−h/2−2λμvy−h/2.
Equation (25) can be expressed as Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 7.Table 7
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaXYUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kW,zXYUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,m,Di+1/2,j,kW,zXYVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kW,zXYVi,j+1,km+1/2=Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kW,zYZVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kU,xYZVi,j+1,km+1/2Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kU,xYZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2U,xYZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2U,xXZUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,yXZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2V,y


#### 3.3.3.XYZ Boundaries

The grid point arrangements of boundaryXYZ are shown in Figure 6. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 8.
Velocity grid points at the first layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0, and in the first layer, U−h/2=Uh/2, V−h/2=Vh/2, and W−h/2=Wh/2. The boundary grid points are shown in Figure 6(a).Case 9.
Velocity grid points at the second layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0; in the second layer, U−h=Uh, V−h=Vh, and W−h=Wh. The boundary grid points are shown in Figure 6(b)–6(d).


### 3.4. Discretization of Heterogeneity Conditions

The densityρ is calculated using the arithmetic average and stored with the velocity grid points. The Lamé parameters λ and μ are calculated using the volume harmonic average and stored with the stress grid points. In this work, the Lamé parameters of the boundary grid points outside the free surface are obtained by the image’s interior grid points. The detailed implementations of ρ and μ are as follows:(26)ρi,j,kx=ρi−1,j,k+ρi+1,j,k2,ρi,j,ky=ρi,j−1,k+ρi,j+1,k2,ρi,j,kz=ρi,j,k−1+ρi,j,k+12,μi,j,kxy=41/μi−1,j,k+1/μi+1,j,k+1/μi,j−1,k+1/μi,j+1,k,μi,j,kyz=41/μi,j−1,k+1/μi,j+1,k+1/μi,j,k−1+1/μi,j,k+1,μi,j,kxz=41/μi−1,j,k+1/μi+1,j,k+1/μi,j,k−1+1/μi,j,k+1,where the superscripts xy, yz, and xz, respectively, represent the stresses Txy, Tyz, and Txz, and the subscripts i, j, and k, respectively, represent the location of the current grid point.


### 3.5. Discretization of PML Computational Domain

The types of grid points in the PML computational domain are identified, and the number of PML grid types is 8 (inner grid type and 7 free surface types) × 9 = 72. In each calculation time step, the PML difference schemes of the internal grid points and free surface grid points are first calculated. Then, the difference schemes of the imaging grid points outside the free surface are calculated. The time derivative terms of the formulas in Section3.3 are required to deal with the PML technique to ensure the seamless combination of the difference scheme for the PML boundary conditions and free surface boundary conditions. Equations (3)–(11) and (19) are discretized with the PML technique.


### 3.6. Algorithm Design for Parallel Computing

Six indicesi,j,k, m,n,l are used to identify each grid point to describe the complex combination of the computational domain, where i,j,k represent the location x,y,z; e.g., x=x0+iΔx/2, y=y0+jΔy/2, and z=z0+kΔz/2; x0,y0,z0 is the starting point used for grid point generation. m,n,l represent the boundary grid point location. Similarly, xm=xb+mΔx/2, yn=yb+nΔy/2, and zl=zb+lΔz/2, xb,yb,zb are the grid points on the free surface boundary. The advantage of using the indices m,n,l is that they prevent the overlap of the indices i,j,k caused by corners and small structural cracks between two very close surfaces. Staggered-grid points are generated automatically by defining the computational domains and determining if a point is within a domain. Moreover, the material parameters of each grid point are defined by the domain properties, and the six pointers of the neighboring points are updated if they exist. The grid generation begins at the minimum XYZ coordinate and ends at the maximum XYZ coordinate of all boundaries.The grid points in the computational domain are divided into seven groups for the convenience and efficiency of parallel and multithread computing. According to the SFD scheme, the seven groups are the following:(1)
Interior and free surface stress grid points.(2)
Stress grid points on the first boundary layer.(3)
Stress grid points on the second boundary layer.(4)
Stress grid points on the third boundary layer.(5)
Interior and free surface velocity grid points.(6)
Velocity grid points on the first boundary layer.(7)
Velocity grid points on the second boundary layer.


## 3.1. Discretization Architecture

The stress image technique is typically used in seismic wave analysis where only a single free surface exists. Since the velocity components and the stress components of the staggered grid are not saved in the same location, the discretization of the free surface becomes more complex, and a half-grid deviation from the original free surface position will occur during the computation [29]. Due to the existence of faces, edges, and corner points, few studies focused on the implementation of the stress image technique in composite structures. The velocity components are located on the physical surface of the computational model, the shear stress terms are also on the surface so that they can be set to zero, and guided elastic waves in structures are analyzed at MHz frequencies [24]; thus, these free surface treatments limit the SFD modeling and generalization ability. The SFD method has been used for the analysis of platform-pile systems without considering the free surface of 3D structures [26, 30, 31]. The computational domain is a regular cuboid domain. Therefore, it is necessary to extend the SFD method for use in platform-pile systems with free surfaces.An SFD for the platform-pile system of composite structures is proposed to provide a systematic and widely applicable solution for the structures. We extend the stress image technique [32] to 3D structures and use multidirectional imaging on the free surfaces. The SFD system for the free surface treatment (720 grid types), heterogeneity, and absorption boundary layers (72 grid types) is implemented using group batch processing, and the calculation of the SFD model is straightforward. It must be pointed out that the overlap of the image grid points of two very close free surfaces (including faces, edges, and corner points) is handled by hash map implementation of the grid point indices to improve the applicability of the SFD scheme. The SFD system is highly adaptable to complex composite structures and is suitable for the simulation of a platform-pile system with a small impact block.


## 3.2. Discretization of Elastic Wave Equations

The (2M)-th-order finite-difference scheme of the first-order spatial derivative of the function ux can be expressed as follows [33]:(17)∂u∂x=1h∑r=1Marux+rh−0.5h−ux−rh+0.5h,where h is the grid size, which can be replaced by Δx, Δy, and Δz, and arr=1,2,…,M are the SFD coefficients of the spatial derivatives.The configuration of the grid points for the finite-difference discretization is shown in Figure2. All stress and velocity components are located in the small cell. S represents the σxx, σyy, and σzz grid points. Txy, Tyz, and Txz, respectively, represent τxy, τyz, and τxz grid points. U, V, and W are the wave velocities in the X-, Y-, and Z-directions.Figure 2
The configuration of the stress and velocity grid points.The details of the discretization equations are as follows:(18)um+1/2=um−1/2+ΔtρDxτxxm+Dyτxym+Dzτxzm,vm+1/2=vm−1/2+ΔtρDxτxym+Dyτyym+Dzτyzm,wm+1/2=wm−1/2+ΔtρDxτxzm+Dyτyzm+Dzτzzm,τxxm+1=τxxm+Δtλ+2μDxum+1/2+λDyvm+1/2+Dzwm+1/2,τyym+1=τyym+ΔtλDxum+1/2+λ+2μDyvm+1/2+Dzwm+1/2,τzzm+1=τzzm+ΔtλDxum+1/2+λDyvm+1/2+λ+2μDzwm+1/2,τxym+1=τxym+ΔtμDxvm+1/2+Dyum+1/2,τyzm+1=τyzm+ΔtμDywm+1/2+Dzvm+1/2,τxzm+1=τxzm+ΔtμDxwm+1/2+Dzum+1/2,where Dx, Dy, and Dz are the x, y, and z derivatives that can be calculated by (17).


## 3.3. Discretization of Free Surface Boundaries

Few studies analyzed free surface treatments of 3D structures using the stress imaging method. Thus, the SFD for the free surface treatment of composite structures is proposed. The types of grid points are analyzed in detail and are divided into 9 cases to clarify the discretization types of the free surface boundaries. Finally, the SFD system can be solved by discretizing all grid points according to their type. It must be pointed out that the complex problem of free surface boundaries is simplified, resulting in a straightforward implementation.The free surfaces are defined as 7 types to identify the types of boundaries: faces (X, Y, Z), edges (XY, YZ, XZ), and eight corner points (XYZ), as exhibited in Figure 3. Moreover, we use one-directional imaging on the boundary faces (Figure 4), two-directional imaging on the boundary edges (Figure 5), and three-directional imaging on the boundary corner points (Figure 6). The types of image grid points for the boundary faces, edges, and corner points are defined as 1 (one-directional imaging) × 3 (faces X, Y, Z) + 2 (two-directional imaging) × 3 (edges XY, YZ, XZ) + 3(three-directional imaging) × 1 (corner points XYZ) = 12 types.Figure 3
The definitions of the free surface boundary.Figure 4
Free surface of boundariesX, Y, and Z. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Figure 5
Free surface of boundariesXY, XZ, and YZ. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Figure 6
Free surface of boundaryXYZ. (a) S grid points on the free surface Z. (b) U grid points on the free surface X. (c) V grid points on the free surface Y. (d) W grid points on the free surface Z.
(a)(b)(c)(d)Further, three boundary layers outside the free surface are required to generate the free surface treatment, namely, the first, second, and third layers. Leth be the mesh step corresponding to the X-, Y-, and Z-directions; h can be replaced with Δx, Δy, and Δz. Let ±h/2,±h, and ±3h/2 be the distances between the first, second, and third boundary layers and the free surface; there are 6 types of boundary layers.Finally, there are 3 (U, V, W) + 6 (Txx, Txy, Txz, Tyy, Tyz, Tzz) =  9 components. Thus, there are (1 + 7) × 9 = 72 grid types in the internal area and on the free surfaces. Further, the number of image grid types is 6 (the type count of boundary layers) × 12 (the types of image points) × 9 = 648. Thus, a total of 720 grid types must be considered for the free surface treatment in the SFD calculation.Note. (+) and (−) values indicate the normal face points in the positive and negative directions of the related axis.The stress componentsσij of the boundary grid points outside the free surface are computed by the interior image components [32]. The details of the SFD discretization of the boundary velocity components and the stress components on the free surface are listed as follows.


### 3.3.1.X, Y, Z Boundaries

We use boundaryZ as an example; there are two types of boundary grid configurations on the free surface, as shown in Figure 4. The first type is the S grid points on the free surface Z, and the second type is the velocity grid points on the free surface. The discretization details are described below.Case 1.
Stress grid points on the free surface.
We useTxx on the boundary Z as an example, as shown in Figure 4(a). When the S grid points are on the free surface Z, the stress component Tzz=0, and Txx (τxx) is computed by (6) and (8).(19)Txxm+1=Txxm+AΔtDU,x+BΔtDV,y,where A=4λμ+4μ2/λ+2μ and B=2λμ/λ+2μ. DU,x is the x derivative of U according to (17). Equation (27) can be expressed as Ti,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,y; similarly, Tyy (τyy) and Tzz (τzz) of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 1.Table 1
Summary of stress components on the free surface.
BoundaryComponentFormulaXTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DW,zXTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DV,yYTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DW,zYTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DU,xZTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,yZTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DU,xNote.i, j, k are the x, y, z indices of grid points on the free surface; m is the time step.Case 2.
Velocity grid points at the first boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
We useW−Δz/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(a). W−Δz/2 is calculated according to the second-order finite-difference approximation of the condition Tzz,t0=0 (Tzz,t is the time derivative of Tzz) and Hooke’s law for Tzz:(20)Wi,j,k−1/2m+1/2=Wi,j,k+1/2m+1/2+Δzλλ+2μDU,x+Δzλλ+2μDV,y.
Equation (20) can be expressed as Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 2.Table 2
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f2Ui+1/2,j,km+1/2,−Δx,DV,y,DW,zX+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f2Ui−1/2,j,km+1/2,Δx,DV,y,DW,zY−Vi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f2Vi,j+1/2,km+1/2,−Δy,DU,x,DW,zY+Vi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f2Vi,j−1/2,km+1/2,Δy,DU,x,DW,zZ−Wi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,yZ+Wi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f2Wi,j,k−1/2m+1/2,Δz,DU,x,DV,yCase 3.
Velocity grid points at the second layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
We useU−h at the first layer outside the boundary Z as an example, as shown in Figure 4(a). U−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzx−h/2=−Tzxh/2 and Hooke’s law for Tzx:(21)U−h=U0+hwx−h/2−hTxz,t−h/2/μ,where wx−h/2=DxW−h/2. Equation (21) can be expressed as Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xz; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 3.Table 3
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Vi−1,j,km+1/2Vi−1,j,km+1/2=f3Vi,j,km+1/2,−Δx,Di−1/2,j,kU,y,Ti−1/2,j,kxyX−Wi−1,j,km+1/2Wi−1,j,km+1/2=f3Wi,j,km+1/2,−Δx,Di−1/2,j,kU,z,Ti−1/2,j,kxzX+Vi+1,j,km+1/2Vi+1,j,km+1/2=f3Vi,j,km+1/2,Δx,Di+1/2,j,kU,y,Ti+1/2,j,kxyX+Wi+1,j,km+1/2Wi+1,j,km+1/2=f3Wi,j,km+1/2,Δx,Di+1/2,j,kU,z,Ti+1/2,j,kxzY−Ui,j−1,km+1/2Ui,j−1,km+1/2=f3Ui,j,km+1/2,−Δy,Di,j−1/2,kV,x,Ti,j−1/2,kxyY−Wi,j−1,km+1/2Wi,j−1,km+1/2=f3Vi,j,km+1/2,−Δy,Di,j−1/2,kV,z,Ti,j−1/2,kyzY+Ui,j+1,km+1/2Ui,j+1,km+1/2=f3Ui,j,km+1/2,Δy,Di,j+1/2,kV,x,Ti,j+1/2,kxyY+Wi,j+1,km+1/2Wi,j+1,km+1/2=f3Vi,j,km+1/2,Δy,Di,j+1/2,kV,z,Ti,j+1/2,kyzZ−Ui,j,k−1m+1/2Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xzZ−Vi,j,k−1m+1/2Vi,j,k−1m+1/2=f3Vi,j,km+1/2,−Δz,Di,j,k−1/2W,y,Ti,j,k−1/2yzZ+Ui,j,k+1m+1/2Ui,j,k+1m+1/2=f3Ui,j,km+1/2,Δz,Di,j,k+1/2W,x,Ti,j,k+1/2xzZ+Vi,j,k+1m+1/2Vi,j,k+1m+1/2=f3Vi,j,km+1/2,Δz,Di,j,k+1/2W,y,Ti,j,k+1/2yzCase 4.
Velocity grid points at the second boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use W−h as the second layer outside the boundary Z as an example, as shown in Figure 4(d). W−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2 and Hooke’s law for Tzz:(22)W−h=W0+−Δzλ+2μ×Tzz,t−h/2−λux−h/2−λvy−h/2.
Equation (22) can be expressed as Wi,j,k−1m+1/2=f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 4.Table 4
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Ui−1,j,km+1/2f4Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,t,Di−1/2,j,kV,y,Di−1/2,j,kW,zX+Ui+1,j,km+1/2f4Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,t,Di+1/2,j,kV,y,Di+1/2,j,kW,zY−Vi,j−1,km+1/2f4Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,t,Di,j−1/2,kU,x,Di,j−1/2,kW,zY+Vi,j+1,km+1/2f4Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,t,Di,j+1/2,kU,x,Di,j+1/2,kW,zZ−Wi,j,k−1m+1/2f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,yZ+Wi,j,k+1m+1/2f4Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,t,Di,j,k+1/2U,x,Di,j,k+1/2V,yCase 5.
Velocity grid points at the first layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use U−h/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(d). U−h/2 is calculated according to the second-order finite-difference approximation of the condition Tzx0=0 and Hooke’s law for Tzx. Moreover, V−h/2 is calculated in the same way.(23)U−h/2=Uh/2+h×wx0,where wx0=DxW0. Equation (23) can be expressed as Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,x; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 5.Table 5
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Vi−1/2,j,km+1/2Vi−1/2,j,km+1/2=f5Vi+1/2,j,km+1/2,−Δx,Di,j,kU,yX−Wi−1/2,j,km+1/2Wi−1/2,j,km+1/2=f5Wi+1/2,j,km+1/2,−Δx,Di,j,kU,zX+Vi+1/2,j,km+1/2Vi+1/2,j,km+1/2=f5Vi,j,k−1/2m+1/2,Δx,Di,j,kU,yX+Wi+1/2,j,km+1/2Wi+1/2,j,km+1/2=f5Wi,j,k−1/2m+1/2,Δx,Di,j,kU,zY−Ui,j−1/2,km+1/2Ui,j−1/2,km+1/2=f5Ui,j+1/2,km+1/2,−Δy,Di,j,kV,xY−Wi,j−1/2,km+1/2Wi,j−1/2,km+1/2=f5Wi,j+1/2,km+1/2,−Δy,Di,j,kV,zY+Ui,j+1/2,km+1/2Ui,j+1/2,km+1/2=f5Ui,j+1/2,km+1/2,Δy,Di,j,kV,xY+Wi,j+1/2,km+1/2Wi,j+1/2,km+1/2=f5Wi,j+1/2,km+1/2,Δy,Di,j,kV,zZ−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,xZ−Vi,j,k−1/2m+1/2Vi,j,k−1/2m+1/2=f5Vi,j,k+1/2m+1/2,−Δz,Di,j,kW,yZ+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f5Ui,j,k−1/2m+1/2,Δz,Di,j,kW,xZ+Vi,j,k+1/2m+1/2Vi,j,k+1/2m+1/2=f5Vi,j,k−1/2m+1/2,Δz,Di,j,kW,y


### 3.3.2.XY, YZ, XZ Boundaries

The grid point arrangements of boundariesXY, XZ, and YZ are shown in Figure 5. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 6.
Velocity grid points at the first boundary layer.
When theS grid points are on the free surface boundary Z, as shown in Figure 5(a), U−h/2 at the first layer of boundary XY is obtained from the second-order finite-difference approximation of the condition Txx0=0, Tyy0=0 and Hooke’s law for Txx, Tyy.(24)U−h/2=Ui−1/2,j,km+1/2=Ui+1/2,j,km+1/2−−Δxλ2λ+μwz.
Equation (24) can be expressed as Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,z; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 6.Table 6
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaXYUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,zXYUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kW,zXYVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kW,zXYVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kW,zYZVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kU,xYZVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kU,xYZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kU,xYZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kU,xXZUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kV,yXZUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kV,yXZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kV,yXZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kV,yCase 7.
Velocity grid points at the second layer.
When theW grid points are on the free surface boundary Z, as shown in Figure 5(d), W−h at the second layer of boundary XZ is calculated by (25) according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2=−Tzzh/2 and Hooke’s law for Tzz.(25)W−h=W0+−Δz4μ2+4λμ×λ+2μTzz,t−h/2−2λμvy−h/2.
Equation (25) can be expressed as Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 7.Table 7
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaXYUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kW,zXYUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,m,Di+1/2,j,kW,zXYVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kW,zXYVi,j+1,km+1/2=Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kW,zYZVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kU,xYZVi,j+1,km+1/2Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kU,xYZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2U,xYZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2U,xXZUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,yXZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2V,y


### 3.3.3.XYZ Boundaries

The grid point arrangements of boundaryXYZ are shown in Figure 6. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 8.
Velocity grid points at the first layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0, and in the first layer, U−h/2=Uh/2, V−h/2=Vh/2, and W−h/2=Wh/2. The boundary grid points are shown in Figure 6(a).Case 9.
Velocity grid points at the second layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0; in the second layer, U−h=Uh, V−h=Vh, and W−h=Wh. The boundary grid points are shown in Figure 6(b)–6(d).


## 3.3.1.X, Y, Z Boundaries

We use boundaryZ as an example; there are two types of boundary grid configurations on the free surface, as shown in Figure 4. The first type is the S grid points on the free surface Z, and the second type is the velocity grid points on the free surface. The discretization details are described below.Case 1.
Stress grid points on the free surface.
We useTxx on the boundary Z as an example, as shown in Figure 4(a). When the S grid points are on the free surface Z, the stress component Tzz=0, and Txx (τxx) is computed by (6) and (8).(19)Txxm+1=Txxm+AΔtDU,x+BΔtDV,y,where A=4λμ+4μ2/λ+2μ and B=2λμ/λ+2μ. DU,x is the x derivative of U according to (17). Equation (27) can be expressed as Ti,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,y; similarly, Tyy (τyy) and Tzz (τzz) of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 1.Table 1
Summary of stress components on the free surface.
BoundaryComponentFormulaXTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DW,zXTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DV,yYTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DW,zYTzzTi,j,kzz,m+1=f1Ti,j,kzz,m,DW,z,DU,xZTxxTi,j,kxx,m+1=f1Ti,j,kxx,m,DU,x,DV,yZTyyTi,j,kyy,m+1=f1Ti,j,kyy,m,DV,y,DU,xNote.i, j, k are the x, y, z indices of grid points on the free surface; m is the time step.Case 2.
Velocity grid points at the first boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
We useW−Δz/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(a). W−Δz/2 is calculated according to the second-order finite-difference approximation of the condition Tzz,t0=0 (Tzz,t is the time derivative of Tzz) and Hooke’s law for Tzz:(20)Wi,j,k−1/2m+1/2=Wi,j,k+1/2m+1/2+Δzλλ+2μDU,x+Δzλλ+2μDV,y.
Equation (20) can be expressed as Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 2.Table 2
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f2Ui+1/2,j,km+1/2,−Δx,DV,y,DW,zX+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f2Ui−1/2,j,km+1/2,Δx,DV,y,DW,zY−Vi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f2Vi,j+1/2,km+1/2,−Δy,DU,x,DW,zY+Vi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f2Vi,j−1/2,km+1/2,Δy,DU,x,DW,zZ−Wi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f2Wi,j,k+1/2m+1/2,−Δz,DU,x,DV,yZ+Wi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f2Wi,j,k−1/2m+1/2,Δz,DU,x,DV,yCase 3.
Velocity grid points at the second layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
We useU−h at the first layer outside the boundary Z as an example, as shown in Figure 4(a). U−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzx−h/2=−Tzxh/2 and Hooke’s law for Tzx:(21)U−h=U0+hwx−h/2−hTxz,t−h/2/μ,where wx−h/2=DxW−h/2. Equation (21) can be expressed as Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xz; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 3.Table 3
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Vi−1,j,km+1/2Vi−1,j,km+1/2=f3Vi,j,km+1/2,−Δx,Di−1/2,j,kU,y,Ti−1/2,j,kxyX−Wi−1,j,km+1/2Wi−1,j,km+1/2=f3Wi,j,km+1/2,−Δx,Di−1/2,j,kU,z,Ti−1/2,j,kxzX+Vi+1,j,km+1/2Vi+1,j,km+1/2=f3Vi,j,km+1/2,Δx,Di+1/2,j,kU,y,Ti+1/2,j,kxyX+Wi+1,j,km+1/2Wi+1,j,km+1/2=f3Wi,j,km+1/2,Δx,Di+1/2,j,kU,z,Ti+1/2,j,kxzY−Ui,j−1,km+1/2Ui,j−1,km+1/2=f3Ui,j,km+1/2,−Δy,Di,j−1/2,kV,x,Ti,j−1/2,kxyY−Wi,j−1,km+1/2Wi,j−1,km+1/2=f3Vi,j,km+1/2,−Δy,Di,j−1/2,kV,z,Ti,j−1/2,kyzY+Ui,j+1,km+1/2Ui,j+1,km+1/2=f3Ui,j,km+1/2,Δy,Di,j+1/2,kV,x,Ti,j+1/2,kxyY+Wi,j+1,km+1/2Wi,j+1,km+1/2=f3Vi,j,km+1/2,Δy,Di,j+1/2,kV,z,Ti,j+1/2,kyzZ−Ui,j,k−1m+1/2Ui,j,k−1m+1/2=f3Ui,j,km+1/2,−Δz,Di,j,k−1/2W,x,Ti,j,k−1/2xzZ−Vi,j,k−1m+1/2Vi,j,k−1m+1/2=f3Vi,j,km+1/2,−Δz,Di,j,k−1/2W,y,Ti,j,k−1/2yzZ+Ui,j,k+1m+1/2Ui,j,k+1m+1/2=f3Ui,j,km+1/2,Δz,Di,j,k+1/2W,x,Ti,j,k+1/2xzZ+Vi,j,k+1m+1/2Vi,j,k+1m+1/2=f3Vi,j,km+1/2,Δz,Di,j,k+1/2W,y,Ti,j,k+1/2yzCase 4.
Velocity grid points at the second boundary layer:U grid points outside the boundary X, V grid points outside the boundary Y, and W grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use W−h as the second layer outside the boundary Z as an example, as shown in Figure 4(d). W−h is calculated according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2 and Hooke’s law for Tzz:(22)W−h=W0+−Δzλ+2μ×Tzz,t−h/2−λux−h/2−λvy−h/2.
Equation (22) can be expressed as Wi,j,k−1m+1/2=f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 4.Table 4
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaX−Ui−1,j,km+1/2f4Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,t,Di−1/2,j,kV,y,Di−1/2,j,kW,zX+Ui+1,j,km+1/2f4Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,t,Di+1/2,j,kV,y,Di+1/2,j,kW,zY−Vi,j−1,km+1/2f4Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,t,Di,j−1/2,kU,x,Di,j−1/2,kW,zY+Vi,j+1,km+1/2f4Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,t,Di,j+1/2,kU,x,Di,j+1/2,kW,zZ−Wi,j,k−1m+1/2f4Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,t,Di,j,k−1/2U,x,Di,j,k−1/2V,yZ+Wi,j,k+1m+1/2f4Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,t,Di,j,k+1/2U,x,Di,j,k+1/2V,yCase 5.
Velocity grid points at the first layer:V, W grid points outside the boundary X; U, W grid points outside the boundary Y; U, V grid points outside the boundary Z.
When theW velocity grid points are on the free surface Z, we use U−h/2 at the first layer outside the boundary Z as an example, as shown in Figure 4(d). U−h/2 is calculated according to the second-order finite-difference approximation of the condition Tzx0=0 and Hooke’s law for Tzx. Moreover, V−h/2 is calculated in the same way.(23)U−h/2=Uh/2+h×wx0,where wx0=DxW0. Equation (23) can be expressed as Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,x; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is given in Table 5.Table 5
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaX−Vi−1/2,j,km+1/2Vi−1/2,j,km+1/2=f5Vi+1/2,j,km+1/2,−Δx,Di,j,kU,yX−Wi−1/2,j,km+1/2Wi−1/2,j,km+1/2=f5Wi+1/2,j,km+1/2,−Δx,Di,j,kU,zX+Vi+1/2,j,km+1/2Vi+1/2,j,km+1/2=f5Vi,j,k−1/2m+1/2,Δx,Di,j,kU,yX+Wi+1/2,j,km+1/2Wi+1/2,j,km+1/2=f5Wi,j,k−1/2m+1/2,Δx,Di,j,kU,zY−Ui,j−1/2,km+1/2Ui,j−1/2,km+1/2=f5Ui,j+1/2,km+1/2,−Δy,Di,j,kV,xY−Wi,j−1/2,km+1/2Wi,j−1/2,km+1/2=f5Wi,j+1/2,km+1/2,−Δy,Di,j,kV,zY+Ui,j+1/2,km+1/2Ui,j+1/2,km+1/2=f5Ui,j+1/2,km+1/2,Δy,Di,j,kV,xY+Wi,j+1/2,km+1/2Wi,j+1/2,km+1/2=f5Wi,j+1/2,km+1/2,Δy,Di,j,kV,zZ−Ui,j,k−1/2m+1/2Ui,j,k−1/2m+1/2=f5Ui,j,k+1/2m+1/2,−Δz,Di,j,kW,xZ−Vi,j,k−1/2m+1/2Vi,j,k−1/2m+1/2=f5Vi,j,k+1/2m+1/2,−Δz,Di,j,kW,yZ+Ui,j,k+1/2m+1/2Ui,j,k+1/2m+1/2=f5Ui,j,k−1/2m+1/2,Δz,Di,j,kW,xZ+Vi,j,k+1/2m+1/2Vi,j,k+1/2m+1/2=f5Vi,j,k−1/2m+1/2,Δz,Di,j,kW,y


## 3.3.2.XY, YZ, XZ Boundaries

The grid point arrangements of boundariesXY, XZ, and YZ are shown in Figure 5. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 6.
Velocity grid points at the first boundary layer.
When theS grid points are on the free surface boundary Z, as shown in Figure 5(a), U−h/2 at the first layer of boundary XY is obtained from the second-order finite-difference approximation of the condition Txx0=0, Tyy0=0 and Hooke’s law for Txx, Tyy.(24)U−h/2=Ui−1/2,j,km+1/2=Ui+1/2,j,km+1/2−−Δxλ2λ+μwz.
Equation (24) can be expressed as Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,z; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is presented in Table 6.Table 6
Summary of velocity components at the first boundary layer.
BoundaryComponentFormulaXYUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kW,zXYUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kW,zXYVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kW,zXYVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kW,zYZVi,j−1/2,km+1/2Vi,j−1/2,km+1/2=f6Vi,j+1/2,km+1/2,−Δy,Di,j,kU,xYZVi,j+1/2,km+1/2Vi,j+1/2,km+1/2=f6Vi,j−1/2,km+1/2,Δy,Di,j,kU,xYZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kU,xYZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kU,xXZUi−1/2,j,km+1/2Ui−1/2,j,km+1/2=f6Ui+1/2,j,km+1/2,−Δx,Di,j,kV,yXZUi+1/2,j,km+1/2Ui+1/2,j,km+1/2=f6Ui−1/2,j,km+1/2,Δx,Di,j,kV,yXZWi,j,k−1/2m+1/2Wi,j,k−1/2m+1/2=f6Wi,j,k+1/2m+1/2,−Δz,Di,j,kV,yXZWi,j,k+1/2m+1/2Wi,j,k+1/2m+1/2=f6Wi,j,k−1/2m+1/2,Δz,Di,j,kV,yCase 7.
Velocity grid points at the second layer.
When theW grid points are on the free surface boundary Z, as shown in Figure 5(d), W−h at the second layer of boundary XZ is calculated by (25) according to the second-order finite-difference approximation of the anti-symmetry condition Tzz−h/2=−Tzzh/2 and Hooke’s law for Tzz.(25)W−h=W0+−Δz4μ2+4λμ×λ+2μTzz,t−h/2−2λμvy−h/2.
Equation (25) can be expressed as Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,y; similarly, the velocity components of the other boundaries can be calculated in the same way. A summary of all the formulas for this case is provided in Table 7.Table 7
Summary of velocity components at the second boundary layer.
BoundaryComponentFormulaXYUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kW,zXYUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti+1/2,j,kxx,m,Di+1/2,j,kW,zXYVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kW,zXYVi,j+1,km+1/2=Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kW,zYZVi,j−1,km+1/2Vi,j−1,km+1/2=f7Vi,j,km+1/2,−Δy,Ti,j−1/2,kyy,m,Di,j−1/2,kU,xYZVi,j+1,km+1/2Vi,j+1,km+1/2=f7Vi,j,km+1/2,Δy,Ti,j+1/2,kyy,m,Di,j+1/2,kU,xYZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2U,xYZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2U,xXZUi−1,j,km+1/2Ui−1,j,km+1/2=f7Ui,j,km+1/2,−Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZUi+1,j,km+1/2Ui+1,j,km+1/2=f7Ui,j,km+1/2,Δx,Ti−1/2,j,kxx,m,Di−1/2,j,kV,yXZWi,j,k−1m+1/2Wi,j,k−1m+1/2=f7Wi,j,km+1/2,−Δz,Ti,j,k−1/2zz,m,Di,j,k−1/2V,yXZWi,j,k+1m+1/2Wi,j,k+1m+1/2=f7Wi,j,km+1/2,Δz,Ti,j,k+1/2zz,m,Di,j,k+1/2V,y


## 3.3.3.XYZ Boundaries

The grid point arrangements of boundaryXYZ are shown in Figure 6. Similarly, the configuration of the free surface is based on the S grid points and velocity grid points on the free surface. The discretization details are described below.Case 8.
Velocity grid points at the first layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0, and in the first layer, U−h/2=Uh/2, V−h/2=Vh/2, and W−h/2=Wh/2. The boundary grid points are shown in Figure 6(a).Case 9.
Velocity grid points at the second layer.
According to the free surface boundary conditionsσxx=0, σyy=0, and σzz=0, there exists ux=vy=wz=0; in the second layer, U−h=Uh, V−h=Vh, and W−h=Wh. The boundary grid points are shown in Figure 6(b)–6(d).


## 3.4. Discretization of Heterogeneity Conditions

The densityρ is calculated using the arithmetic average and stored with the velocity grid points. The Lamé parameters λ and μ are calculated using the volume harmonic average and stored with the stress grid points. In this work, the Lamé parameters of the boundary grid points outside the free surface are obtained by the image’s interior grid points. The detailed implementations of ρ and μ are as follows:(26)ρi,j,kx=ρi−1,j,k+ρi+1,j,k2,ρi,j,ky=ρi,j−1,k+ρi,j+1,k2,ρi,j,kz=ρi,j,k−1+ρi,j,k+12,μi,j,kxy=41/μi−1,j,k+1/μi+1,j,k+1/μi,j−1,k+1/μi,j+1,k,μi,j,kyz=41/μi,j−1,k+1/μi,j+1,k+1/μi,j,k−1+1/μi,j,k+1,μi,j,kxz=41/μi−1,j,k+1/μi+1,j,k+1/μi,j,k−1+1/μi,j,k+1,where the superscripts xy, yz, and xz, respectively, represent the stresses Txy, Tyz, and Txz, and the subscripts i, j, and k, respectively, represent the location of the current grid point.


## 3.5. Discretization of PML Computational Domain

The types of grid points in the PML computational domain are identified, and the number of PML grid types is 8 (inner grid type and 7 free surface types) × 9 = 72. In each calculation time step, the PML difference schemes of the internal grid points and free surface grid points are first calculated. Then, the difference schemes of the imaging grid points outside the free surface are calculated. The time derivative terms of the formulas in Section3.3 are required to deal with the PML technique to ensure the seamless combination of the difference scheme for the PML boundary conditions and free surface boundary conditions. Equations (3)–(11) and (19) are discretized with the PML technique.


## 3.6. Algorithm Design for Parallel Computing

Six indicesi,j,k, m,n,l are used to identify each grid point to describe the complex combination of the computational domain, where i,j,k represent the location x,y,z; e.g., x=x0+iΔx/2, y=y0+jΔy/2, and z=z0+kΔz/2; x0,y0,z0 is the starting point used for grid point generation. m,n,l represent the boundary grid point location. Similarly, xm=xb+mΔx/2, yn=yb+nΔy/2, and zl=zb+lΔz/2, xb,yb,zb are the grid points on the free surface boundary. The advantage of using the indices m,n,l is that they prevent the overlap of the indices i,j,k caused by corners and small structural cracks between two very close surfaces. Staggered-grid points are generated automatically by defining the computational domains and determining if a point is within a domain. Moreover, the material parameters of each grid point are defined by the domain properties, and the six pointers of the neighboring points are updated if they exist. The grid generation begins at the minimum XYZ coordinate and ends at the maximum XYZ coordinate of all boundaries.The grid points in the computational domain are divided into seven groups for the convenience and efficiency of parallel and multithread computing. According to the SFD scheme, the seven groups are the following:(1)
Interior and free surface stress grid points.(2)
Stress grid points on the first boundary layer.(3)
Stress grid points on the second boundary layer.(4)
Stress grid points on the third boundary layer.(5)
Interior and free surface velocity grid points.(6)
Velocity grid points on the first boundary layer.(7)
Velocity grid points on the second boundary layer.


## 4. Methodologies for Parallel Velocity Analysis of the Platform-Pile System

The LST has been widely used to monitor the integrity of pile foundations. However, due to the three-dimensional effects, further research and improvements of the LST are required for platform-pile systems of bridge foundations and seaport wharves. The DVSSM has been proposed to eliminate the high-frequency interference during pipe pile integrity testing [11]. The PS method was proposed for detecting piles in layered soil [10]. A triangular wedge was used to provide a longitudinal excitation signal for the integrity assessment of bridge piles and estimate the pile length [7]. Due to the reflection waves generated by the superstructures such as caps and beams, the vibration wave signals acquired by the sensors tend to oscillate, exhibiting coupling of the upward and downward waves. Therefore, the upward waves generated by the damaged area are difficult to identify. The upward algorithm for parallel velocity analysis was proposed to extract the downward wave signals and obtain meaningful upward reflected signals from the vibration signal [15].(27)V1↓t=V1t−V2t−T+V1↓t−2T,V1↑t=V1t−V1↓t,where V1t and V2t are the parallel velocities, V1↓t is the downward velocity component measured at location A1, and V1↑t is the upward velocity of sensor 1. The locations of the sensors and the impact hammer must be considered to prevent 3D vibration and obtain meaningful signals for damage identification using this equation. The locations considered in this study are shown in Figures 7-8.Figure 7
Illustration of parallel wave testing. (a) The two-pile model. (b) The pile-soil model.Figure 8
Impact locations. (a) The two-pile model. (b) The pile-soil model.
(a)(b)


## 5. Results and Discussion

### 5.1. Model Information

The two-pile model consists of two square piles and a cap. The size of the pile section is400×400mm, the pile length is 8.0 m, and the size of the pile cap is 1.6×0.8×0.8m. According to (2), a wood partition is used for the damage simulation of a solid structure, and the reflected wave is generated by changing the material parameters in the pile length direction. The damage is located at z = 5.3–5.34 m, and the thickness of the damaged area is 0.04 m, as shown in Figures 9-10. The impact locations at the top surface of the pile cap are (0.0, 0, 0), (−0.2, 0, 0), and (−0.4, 0, 0), as shown in Figure 8; the impact location at the lateral surface of the pile is (0.0, 0.2, 1.3). The computational models are shown in Figure 1. The parameters of the two-pile model are listed in Table 8, and the theoretical integrity coefficient β is 0.36. The SFD numerical model of the two-pile model was constructed with the following parameters: Δx=Δy=0.01 m, Δz=0.02 m, and Δt=1μs. The impact parameters of (12) are I=1.0, t0=1 ms, and r=0.05 m. The W velocities of the points A1 (0.0, 0.2, 2.1), A2 (0.0, 0.2, 3.1), B1 (–0.2, 0.0, 2.1), and B2 (–0.2, 0.0, 3.1) are monitored and analyzed.Figure 9
The damage of the two-pile model (the size of the middle hole is13×13 cm).Figure 10
The arrangement of the concrete impact block and sensors of the two-pile model.Table 8
The parameters of the two-pile model.
MaterialE (GPa)Poisson’s ratioρ (kg/m3)Concrete350.242400Wood80.30850The pile-soil model has the following soil parameters:S wave velocity = 100 m/s, Poisson’s ratio ν = 0.35, and density =  1600kg/m3. The damage is located at z = 5.0–5.2 m. The soil location is x = –1 m, 1 m; y = –1 m, 1 m; z = 4–10 m. The other parameters of the pile-soil model are listed in Table 9, and the theoretical integrity coefficient β is 0.49. The computational domain is exhibited in Figure 1. The impact locations at the top surface of the pile cap are (0.0, 0, 0), (0.0, 0.2, 0.0), and (0.0, 0.4, 0.0) as shown in Figure 8, and the impact location at the lateral surface of pile is (0.0, 0.2, 1.3). The W velocities of points B1 and B2 are monitored and analyzed. The step size is Δx=Δy=0.02 m, Δz=0.04 m, and t=2μs. The thickness of the PML is 0.4 m in the X- and Y-directions and 0.8 m in the Z-direction to meet the grid requirements of the PML.Table 9
The parameters of the pile-soil model.
MaterialE (GPa)Poisson’s ratioρ (kg/m3)Pile and cap350.242400Damage100.302000The computational models of the undamaged domains (the pile cap, pile body, and soil), the damaged domains (the wood), and the PML domains are defined using the material parameters in Tables8 and 9. The SFD models are shown in Figure 11.Figure 11
The SFD models. (a) The two-pile model. (b) The pile-soil model.
(a)(b)


### 5.2. Verification of SFD Simulation

The SFD results of the two-pile model were compared with the results of the ABAQUS software to verify the accuracy and reliability of the numerical simulation [34]. A PIT-VV pile integrity tester with two accelerometers was used to measure the velocity data at B1 and B2. The impact location at the top surface of the pile cap is (0.0, 0, 0).Figure12 shows the ABAQUS model. Figures 13-14 show that the results of ABAQUS and SFD are consistent, and the velocities at A1, A2, B1, and B2 are the same. Figures 15–19 reveal that the peak values and phases are the same. Figures 20-21 indicate that the upward velocities of the test data and SFD model are consistent. These results verify the accuracy of the SFD algorithm and the SFD differential system for the free surface treatment.Figure 12
The pile model in ABAQUS software.Figure 13
TheW field of the ABAQUS model at t = 2.1 ms.Figure 14
TheW field of the SFD model at t = 2.1 ms.Figure 15
The velocity componentW at A1.Figure 16
The velocity componentW at A2.Figure 17
The velocity componentW at B1.Figure 18
The velocity componentW at B2.Figure 19
The upward velocityW at A1.Figure 20
The velocity componentW at B1 and B2 (impact location at y = 0.0 m).Figure 21
The upward velocity componentW at B1 and B2 (impact location at y = 0.0 m).


### 5.3. Results of Two-Pile Model

When the nylon hammer struck the top surface of the pile cap, downward stress waves were generated and reflected by the cap. The upward and downward stress waves scattered and propagated in the longitudinal direction of the pile structure, as shown in Figures22–25. The wave velocities of locations A1, A2, B1, and B2 were monitored to determine the most suitable excitation position of the impact (Figures 26–32). Due to the 3D effects and asymmetric excitation, interference was caused by the pile cap; therefore, the two signals have a large difference, as shown in Figures 26–28. The identification of the defect signal was difficult due to the coupling of the interference signal propagating along the pile. When the impact location was on the symmetry axis of the top surface, the impact signal and toe signal of the velocities received by A1 and A2 were consistent, as shown in Figures 29 and 31. The damage signal first arrived at A2 and B2 and then at A1 and B1. When the sensor locations were on the bending plane relative to the impact location, the peak values of the downward velocities at A1 and A2 were the same. When the damage signals were monitored at the same time, the downward and upward waves were superimposed at sensor location A2, resulting in an increase in the peak value relative to the peak value received at A1; thus, the damage signal could be identified (Figures 29–32). The damage signals at the two monitoring points had similar frequencies and amplitudes, as shown in the results of the lateral impact. The interference was less pronounced when the impact location was on the symmetry axis of the platform-pile system.Figure 22
The profile of the velocity componentW of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (impact location x = −0.2).Figure 23
The profile of the velocity componentW of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (impact location y = 0.2).Figure 24
TheW profile of model 1 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (lateral impact).Figure 25
The profile of the velocity componentW at time t = 2.25 ms (lateral impact).Figure 26
The velocityW at A1 and A2 (impact location y = 0.2).Figure 27
The velocityW at B1 and B2 (impact location y = 0.2).Figure 28
The velocityW at B1 and B2 (impact location x = −0.2).Figure 29
The velocityW at A1 and A2 (impact location x = 0.0).Figure 30
The velocityW at A1 and A2 (impact location x = −0.2).Figure 31
The velocityW at A1 and A2 (impact location x = −0.4).Figure 32
The velocityW of model 1 at B1 and B2 (lateral impact).


### 5.4. Results of Pile-Soil Model

The pile-soil model was biaxially symmetric. The results of B1 and B2 were significantly different for the same impact duration of 0.5 ms. Figures33 and 34 show that the frequency of the received signals from the top impact was higher than that from the lateral impact. The signal showed attenuation because of the surrounding soil. Figures 35–37 show that the reflected waves generated by the damage cannot be identified due to the periodic reflection generated by the platform.Figure 33
The velocityW of model 2 at B1 and B2 (impact from top surface).Figure 34
The velocityW of model 2 at B1 and B2 (impact from lateral surface).Figure 35
TheW profile of model 2 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (top impact).Figure 36
TheW profile of model 2 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (lateral impact).Figure 37
The profile of the velocity componentW at time t = 2.20 ms (lateral impact).


### 5.5. Parallel Wave Analysis of the Platform-Pile Systems

The upward velocity results of the different impact locations are shown in Figures38–44. Due to the interference signals generated by the pile cap, the upward wave velocity curves in Figures 38–40 show large oscillations. Although the oscillation amplitude was high at the damage locations and the pile toe, it was difficult to identify the location and size of the damage. The material and size of the hammer determine the impact duration and should be considered to reduce the high-frequency interferences and obtain more reasonable data. Figures 41 and 42 show that the peak value of the damage signal was smaller when the impact location deviated from the center point of the pile cap. The impact location should not be outside the boundary of the cross section of the pile to obtain meaningful damage signals. Figures 42–44 show that the upward wave curves were not affected by oscillation and the damage locations could be detected. Furthermore, for the same impact parameters, the peak value of the excitation position 0.2 was relatively higher, and the amplitudes of the damage and pile toe were larger. Therefore, the impact location should be on the symmetry axis, and the monitoring points should be placed on the neutral bending plane relative to the impact location. Figures 43 and 44 show that the upward velocity signal of the damage is easily identified because the damage signals are less affected by the pile cap. Therefore, the impact should be placed on the biaxial symmetry axis of the pile cap, and the sensors should be placed on the bending plane relative to the impact location and symmetry axis of the pile cap.Figure 38
The upward velocitiesW of the two-pile model at A1 (impact location: y = 0, 0.2, 0.4).Figure 39
The upward velocitiesW of the two-pile model at B1 (impact location: y = 0, 0.2, 0.4).Figure 40
The upward velocitiesW of the two-pile model at B1 (impact location: x = 0, −0.2, −0.4).Figure 41
The upward velocitiesW of the two-pile model at A1 (impact location: x = 0, −0.2, −0.4).Figure 42
The upward velocitiesW of the pile-soil model at A1 (top impact).Figure 43
The upward velocitiesW of the two-pile model at B1 (lateral impact).Figure 44
The upward velocitiesW of the pile-soil model at B1 (lateral impact).


## 5.1. Model Information

The two-pile model consists of two square piles and a cap. The size of the pile section is400×400mm, the pile length is 8.0 m, and the size of the pile cap is 1.6×0.8×0.8m. According to (2), a wood partition is used for the damage simulation of a solid structure, and the reflected wave is generated by changing the material parameters in the pile length direction. The damage is located at z = 5.3–5.34 m, and the thickness of the damaged area is 0.04 m, as shown in Figures 9-10. The impact locations at the top surface of the pile cap are (0.0, 0, 0), (−0.2, 0, 0), and (−0.4, 0, 0), as shown in Figure 8; the impact location at the lateral surface of the pile is (0.0, 0.2, 1.3). The computational models are shown in Figure 1. The parameters of the two-pile model are listed in Table 8, and the theoretical integrity coefficient β is 0.36. The SFD numerical model of the two-pile model was constructed with the following parameters: Δx=Δy=0.01 m, Δz=0.02 m, and Δt=1μs. The impact parameters of (12) are I=1.0, t0=1 ms, and r=0.05 m. The W velocities of the points A1 (0.0, 0.2, 2.1), A2 (0.0, 0.2, 3.1), B1 (–0.2, 0.0, 2.1), and B2 (–0.2, 0.0, 3.1) are monitored and analyzed.Figure 9
The damage of the two-pile model (the size of the middle hole is13×13 cm).Figure 10
The arrangement of the concrete impact block and sensors of the two-pile model.Table 8
The parameters of the two-pile model.
MaterialE (GPa)Poisson’s ratioρ (kg/m3)Concrete350.242400Wood80.30850The pile-soil model has the following soil parameters:S wave velocity = 100 m/s, Poisson’s ratio ν = 0.35, and density =  1600kg/m3. The damage is located at z = 5.0–5.2 m. The soil location is x = –1 m, 1 m; y = –1 m, 1 m; z = 4–10 m. The other parameters of the pile-soil model are listed in Table 9, and the theoretical integrity coefficient β is 0.49. The computational domain is exhibited in Figure 1. The impact locations at the top surface of the pile cap are (0.0, 0, 0), (0.0, 0.2, 0.0), and (0.0, 0.4, 0.0) as shown in Figure 8, and the impact location at the lateral surface of pile is (0.0, 0.2, 1.3). The W velocities of points B1 and B2 are monitored and analyzed. The step size is Δx=Δy=0.02 m, Δz=0.04 m, and t=2μs. The thickness of the PML is 0.4 m in the X- and Y-directions and 0.8 m in the Z-direction to meet the grid requirements of the PML.Table 9
The parameters of the pile-soil model.
MaterialE (GPa)Poisson’s ratioρ (kg/m3)Pile and cap350.242400Damage100.302000The computational models of the undamaged domains (the pile cap, pile body, and soil), the damaged domains (the wood), and the PML domains are defined using the material parameters in Tables8 and 9. The SFD models are shown in Figure 11.Figure 11
The SFD models. (a) The two-pile model. (b) The pile-soil model.
(a)(b)


## 5.2. Verification of SFD Simulation

The SFD results of the two-pile model were compared with the results of the ABAQUS software to verify the accuracy and reliability of the numerical simulation [34]. A PIT-VV pile integrity tester with two accelerometers was used to measure the velocity data at B1 and B2. The impact location at the top surface of the pile cap is (0.0, 0, 0).Figure12 shows the ABAQUS model. Figures 13-14 show that the results of ABAQUS and SFD are consistent, and the velocities at A1, A2, B1, and B2 are the same. Figures 15–19 reveal that the peak values and phases are the same. Figures 20-21 indicate that the upward velocities of the test data and SFD model are consistent. These results verify the accuracy of the SFD algorithm and the SFD differential system for the free surface treatment.Figure 12
The pile model in ABAQUS software.Figure 13
TheW field of the ABAQUS model at t = 2.1 ms.Figure 14
TheW field of the SFD model at t = 2.1 ms.Figure 15
The velocity componentW at A1.Figure 16
The velocity componentW at A2.Figure 17
The velocity componentW at B1.Figure 18
The velocity componentW at B2.Figure 19
The upward velocityW at A1.Figure 20
The velocity componentW at B1 and B2 (impact location at y = 0.0 m).Figure 21
The upward velocity componentW at B1 and B2 (impact location at y = 0.0 m).


## 5.3. Results of Two-Pile Model

When the nylon hammer struck the top surface of the pile cap, downward stress waves were generated and reflected by the cap. The upward and downward stress waves scattered and propagated in the longitudinal direction of the pile structure, as shown in Figures22–25. The wave velocities of locations A1, A2, B1, and B2 were monitored to determine the most suitable excitation position of the impact (Figures 26–32). Due to the 3D effects and asymmetric excitation, interference was caused by the pile cap; therefore, the two signals have a large difference, as shown in Figures 26–28. The identification of the defect signal was difficult due to the coupling of the interference signal propagating along the pile. When the impact location was on the symmetry axis of the top surface, the impact signal and toe signal of the velocities received by A1 and A2 were consistent, as shown in Figures 29 and 31. The damage signal first arrived at A2 and B2 and then at A1 and B1. When the sensor locations were on the bending plane relative to the impact location, the peak values of the downward velocities at A1 and A2 were the same. When the damage signals were monitored at the same time, the downward and upward waves were superimposed at sensor location A2, resulting in an increase in the peak value relative to the peak value received at A1; thus, the damage signal could be identified (Figures 29–32). The damage signals at the two monitoring points had similar frequencies and amplitudes, as shown in the results of the lateral impact. The interference was less pronounced when the impact location was on the symmetry axis of the platform-pile system.Figure 22
The profile of the velocity componentW of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (impact location x = −0.2).Figure 23
The profile of the velocity componentW of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (impact location y = 0.2).Figure 24
TheW profile of model 1 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (lateral impact).Figure 25
The profile of the velocity componentW at time t = 2.25 ms (lateral impact).Figure 26
The velocityW at A1 and A2 (impact location y = 0.2).Figure 27
The velocityW at B1 and B2 (impact location y = 0.2).Figure 28
The velocityW at B1 and B2 (impact location x = −0.2).Figure 29
The velocityW at A1 and A2 (impact location x = 0.0).Figure 30
The velocityW at A1 and A2 (impact location x = −0.2).Figure 31
The velocityW at A1 and A2 (impact location x = −0.4).Figure 32
The velocityW of model 1 at B1 and B2 (lateral impact).


## 5.4. Results of Pile-Soil Model

The pile-soil model was biaxially symmetric. The results of B1 and B2 were significantly different for the same impact duration of 0.5 ms. Figures33 and 34 show that the frequency of the received signals from the top impact was higher than that from the lateral impact. The signal showed attenuation because of the surrounding soil. Figures 35–37 show that the reflected waves generated by the damage cannot be identified due to the periodic reflection generated by the platform.Figure 33
The velocityW of model 2 at B1 and B2 (impact from top surface).Figure 34
The velocityW of model 2 at B1 and B2 (impact from lateral surface).Figure 35
TheW profile of model 2 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (top impact).Figure 36
TheW profile of model 2 of the YZ plane at x=0.0 m and the XZ plane at y=0.0 m at time t=2.25 ms (lateral impact).Figure 37
The profile of the velocity componentW at time t = 2.20 ms (lateral impact).


## 5.5. Parallel Wave Analysis of the Platform-Pile Systems

The upward velocity results of the different impact locations are shown in Figures38–44. Due to the interference signals generated by the pile cap, the upward wave velocity curves in Figures 38–40 show large oscillations. Although the oscillation amplitude was high at the damage locations and the pile toe, it was difficult to identify the location and size of the damage. The material and size of the hammer determine the impact duration and should be considered to reduce the high-frequency interferences and obtain more reasonable data. Figures 41 and 42 show that the peak value of the damage signal was smaller when the impact location deviated from the center point of the pile cap. The impact location should not be outside the boundary of the cross section of the pile to obtain meaningful damage signals. Figures 42–44 show that the upward wave curves were not affected by oscillation and the damage locations could be detected. Furthermore, for the same impact parameters, the peak value of the excitation position 0.2 was relatively higher, and the amplitudes of the damage and pile toe were larger. Therefore, the impact location should be on the symmetry axis, and the monitoring points should be placed on the neutral bending plane relative to the impact location. Figures 43 and 44 show that the upward velocity signal of the damage is easily identified because the damage signals are less affected by the pile cap. Therefore, the impact should be placed on the biaxial symmetry axis of the pile cap, and the sensors should be placed on the bending plane relative to the impact location and symmetry axis of the pile cap.Figure 38
The upward velocitiesW of the two-pile model at A1 (impact location: y = 0, 0.2, 0.4).Figure 39
The upward velocitiesW of the two-pile model at B1 (impact location: y = 0, 0.2, 0.4).Figure 40
The upward velocitiesW of the two-pile model at B1 (impact location: x = 0, −0.2, −0.4).Figure 41
The upward velocitiesW of the two-pile model at A1 (impact location: x = 0, −0.2, −0.4).Figure 42
The upward velocitiesW of the pile-soil model at A1 (top impact).Figure 43
The upward velocitiesW of the two-pile model at B1 (lateral impact).Figure 44
The upward velocitiesW of the pile-soil model at B1 (lateral impact).


## 6. Conclusions

A 3D SFD method was proposed for wave propagation analysis of a platform-pile system considering the free surface, material heterogeneity, and soil environment. The accuracy and feasibility of the algorithm were verified by ABAQUS results and testing. Different impact locations and sensor configurations were considered. The following conclusions were obtained:(1)
The comparison of the proposed SFD method with the ABAQUS and testing result demonstrated the accuracy of damage identification, indicating the suitability of the method for the analysis of similar problems.(2)
The reflected waves generated by the platform exhibited periodic oscillations. Choosing an appropriate impact location reduced the influence of the periodic oscillation signal. Excitation points at the center of symmetry were the best locations.(3)
At the same impact energy, the impact location should not be outside the boundary of the cross section of the pile when the impact occurs on the top surface of the platform to achieve meaningful damage signals.(4)
When the impact locations were on the lateral surface of the pile, the interferences were relatively small, and the damage signal and secondary damage signal could be distinguished.(5)
The material and size of the hammer determine the impact duration and should be chosen to reduce the high-frequency interferences and obtain reasonable data.



---

*Source: 1431774-2022-05-16.xml*