1 Introduction

The study of wave propagation and dynamic response of poroelastic soils is greatly significant for geotechnical engineering, vibration science, and seismology. For example, structures exposed to dynamic loadings such as the large power equipment and machine foundations are common in practical engineering. Moreover, railways built on the viaduct also vibrate the surroundings at the time of train operation. Therefore, the deformation caused by the dynamic vibration should be limited to an acceptable level, and researches of this subject have attracted wide attention from scholars and engineers in many fields. Biot ^[1-3] built wave propagation rules in the three-dimensional consolidation, which laid the foundation for further studies in poroelastic soil media. Since then, other researchers have carried out numerous investigations on dynamic responses of poroelastic soils with various analytical and numerical methods based on this theory. Paul ^[4-5] firstly used Biot's theory to solve Lamb's problem ^[6] in a saturated half-space by using the Helmholtz decomposition, the Hankel transform, and the modified Cagniard method ^[7]. Halpern and Christiano ^[8-9] studied vibrations of saturated soils due to time-harmonic loadings on the surface. Philippacopoulos ^[10-11] used the Helmholtz decomposition and Hankel transform to tackle with Lamb's problem in poroelastic soil media and took the soil's viscosity into consideration. Senjuntichai and Rajapakse ^[12] presented an analytical and numerical treatment of the dynamic response of a dissipative poroelastic half-plane. Rajapakse and Senjuntichai ^[13] used an exact matrix method to evaluate the displacements and stresses of a multilayered poroelastic medium subject to time-harmonic loads. Zhang and Huang ^[14] studied the non-axisymmetrical dynamic response of transversely isotropic saturated poroelastic media by means of Fourier expansion and Hankel transform, and then the work was extended to general cases in the rectangular coordinate by Wang and Huang ^[15]. Cai et al. ^[16-17] considered the dynamic response of a single soil layer with subjacent rock-stratum subject to an axisymmetric harmonic excitation, and obtained fundamental solutions for general stresses and displacements under various types of excitations by means of Laplace-Hankel transform. Ding et al. ^[18] presented Lamb's integral formulae with drainage for a two-phase saturated medium. Zhou et al. ^[19] and He et al. ^[20] derived the dynamic 2.5-D Green's function for a poroelastic half-space and a layered poroelastic half-space, respectively. Zheng et al. ^[21] presented the three-dimensional dynamic Green's functions for a poroelastic half-space due to an internal point load or a fluid source. Zheng et al. ^[22] derived the three-dimensional dynamic Green's functions in a multilayered poroelastic half-space with the propagator matrix method. Lu and Hanyga ^[23] obtained the fundamental solution for a layered porous half-space due to a vertical point force or a point fluid source using the transmission and reflection matrices (TRM) method. Liu and Zhao ^[24] used the generalized transfer matrix method to analyze the dynamic response of multilayered poroelastic media. Liu et al. ^[25] presented a modified stiffness matrix method for the dynamic response of a poroelastic layered half-space subject to a three-dimensional concentrated load.

Recently, the first author and his coworkers developed a stable and efficient computational approach named as the analytical layer-element method for the quasi-static analysis of multilayered thermoelastic media ^[26-27], poroelastic media ^[28-29], and porous thermoelastic media ^[30-31] as well as the dynamic analysis of multilayered transversely isotropic elastic media ^[32-35]. The analytical layer-element method ^[26-35] only involves negative exponentials functions in the exact stiffness matrices of media, which not only simplifies the computational processes but also improves the numerical efficiency and stability. Therefore, this study is aimed at extending the analytical layer-element method to solve the multilayered poroelastic soils due to an axisymmetric harmonic excitation and investigate effects of soil properties on dynamic responses of soil media. By combining equilibrium equations, fluid movement equations and the seepage continuity equation, the analytical layer-elements for a single layer and the underlying half-space are derived with the aid of algebraic manipulation and Hankel transform. A global matrix for multilayered soil media can be established based on the adjacent continuity conditions of general stresses and displacements in the transform domain. Solutions in the frequency domain are further obtained with the inverse Hankel transform. Numerical examples are used to confirm the precision of the present method and demonstrate effects of soil parameters and load conditions on dynamic responses of poroelastic soil media.

2 Analytical layer-element for single poroelastic soil layer

The problem in this study is simulated by an isotropic poroelastic soil system subject to a harmonic excitation Pe^iωt with a radius of r₀ on the surface, in which P is the amplitude of external load, ω is the circular frequency, and . In this section, a single layer with a finite thickness is investigated firstly. Before we start with the analytical derivation, some parameters should be denoted. G, λ, ρ, h, k, and n are the shear modulus, Lame's constant, density, soil thickness, fluid permeability, and porosity of the bulk material, respectively. g is the gravitational acceleration, and ρ_f is the density of pore water.

The governing equations of a poroelastic soil medium for an axisymmetric problem in a cylindrical coordinate system are used in this paper. With the assumption of negligible body forces and incompressibility of soil skeleton and pore water, the equilibrium equations in rectangular coordinates are given by Biot ^[2] as follows:

(1a)

(1b)

(1c)

where σ_x, σ_y, and σ_v denote normal stresses components in the x-, y-, and z-directions, respectively. σ_xy, σ_xz, and σ_zy are the shear stresses in the xy-, xz-, and zy-planes, respectively. u_s, v_s, and w_s are the displacement components of soil skeleton in the x-, y-, and z-directions, respectively, and v_x, v_y, and v_z are the displacement components of fluid in the x-, y-, and z-directions, respectively. The over dot is the derivation with respect to t.

The equilibrium equations in cylindrical coordinates can be obtained from the equilibrium equations (1a)-(1c) in rectangular coordinates, and the relationships of expression in rectangular coordinates and cylindrical coordinates at the same position are written as

(2a)

(2b)

(2c)

The relationships of stress components in rectangular coordinates and cylindrical coordinates are written as ^[36]

(3a)

(3b)

(3c)

(3d)

(3e)

(3f)

where σ_r and σ_θ denote normal stress components in the r-and θ-directions, respectively, and σ_rθ, σ_rz, and σ_θz are the shear stresses in the rθ-, zr-, and zθ-planes, respectively.

Based on Eqs. (1) -(3), the equilibrium equations in cylindrical coordinates can be obtained as follows:

(4a)

(4b)

(4c)

where u, v, and w are the displacement components of soil skeleton in the r-, θ-, and z-directions, and v_r, v_θ, and v_z are the displacement components of fluid in the r-, θ-, and z-directions.

For axisymmetric problems, the displacement components of soil grain and pore fluid in the θ-direction v=v_θ=0, and the shear stresses in the zθ-and θr-planes σ_zθ =σ_rθ=0. Therefore, we have

(5a)

(5b)

In accordance with the principle of effective stress, it is convenient to separate the total stress σ into the effective stress σ' and the pore pressure p^[37]. Define the tension stress as positive, and we have

(6a)

(6b)

(6c)

Substitution of Eqs. (6a)-(6c) into Eqs. (5a)-(5b) yields

(7a)

(7b)

According to the constitutive equations and the principle of effective stress, we obtain

(8a)

(8b)

(8c)

(8d)

in which .

Combined with elasticity equations, the effective stresses can be replaced by the displacement of soil grain. Thus, we have

(9a)

(9b)

The balance equations of fluid motion are given as

(10a)

(10b)

in which n denotes the porosity. The porosity n is defined as the ratio of the volume of voids V_v to the total volume of soil V.

The seepage continuity equation is given as

(11)

where α and M are Biot's parameters accounting for compressibility of soil grain and pore fluid. According to Ref. [1], α =1 and M=∈∞ stand for soils with incompressible constituents. Therefore, Eq. (11) can be simplified as

(12)

In the analytical derivation of steady state, the motion under consideration is set as time-harmonic with a factor e^iωt. That is to say, σ_r, σ_θ, σ_v, σ_rz, p, u, w, v_r, and v_z are expressed as σ_re^iωt, σ_θe^iωt, σ_ze^iωt, σ_rz e^iωt, pe^iωt, ue^iωt, we^iωt, v_r e^iωt, and v_z e^iωt, respectively. Then, substitution of these components into Eqs. (9a) and (9b), the fluid motion equations (10a) and (10b), and the seepage continuity equation (12) leads to

(13a)

(13b)

(13c)

(13d)

(13e)

where .

Similarly, by substituting the expressions of v_r and v_z into Eqs. (13a), (13b), and (13e) and applying the Hankel transform ^[38] to the achieved equations, we have

(14a)

(14b)

(14c)

in which a=ω²ρ + ω⁴ρ_f²c, and b=1+ω²ρ_fc. u, w, and p are the Hankel transform expressions of u, w, and p, respectively, and ξ is the Hankel transform parameter.

For convenience, two matrix vectors are defined as follows:

(15)

(16)

Thus, Eqs. (14a)-(14c) can be expressed in a matrix form as

(17)

where , 0_3×3 and I_3×3 are the zero matrix and the identity matrix of corresponding order, and

Let W(ξ, z)=X(ξ)e^Λz, and then Eq. (17) can be written as

(18)

Its characteristic equation is expressed as

(19)

By means of mathematical computation, the eigenvalues for the above equation are calculated as

(20)

By solving the eigenvalue problem, the solution for the matrix equation can be further obtained as follows:

(21)

Substituting Eq. (21) into Eq. (17), we obtain

(22)

where , , and l₁, …, l₆ denote arbitrary constants. Here, , and .

Substituting Eq. (22) into Eq. (21) yields

(23)

in which , and

In the meantime, the flow quantity through the unit cross section area Q (which is simplified from the time-harmonic form Qe^iωt) is expressed as

(24)

Applying the Hankel transform to Eqs.(8a), (8c), and (24), the following equations can be obtained:

(25a)

(25b)

(25c)

Substitution of Eq. (23) into Eqs. (25a)-(25c) yields

(26)

in which , and

With the aid of Eqs. (23) and (26), the relationship between general stresses and displacements in the transform domain for a single layer, which is shown in Fig. 1, can be written in a matrix form,

(27)

where , and Φ is a matrix with an order of 6× 6 and indicates the relationship between the general stress vector and the general displacement vector.

For the underlying half-space, the relationship between V(ξ, z) and U(ξ, z) can be established by considering z→∞ in Eqs. (23) and (26),

(28)

where Φ₁ is a matrix of order 3× 3, which indicates the relationship between the general stress vector and the general displacement vector for the underlying half-space in the transform domain.

Elements of Φ and Φ₁ in Eqs. (27) and (28) are provided in Appendix A and Appendix B, respectively.

3 Solutions of multilayered poroelastic soils

Due to the natural sedimentation, the soil media often show the stratification characteristic. It is therefore relevant to the study of the dynamic response of multilayered soils subject to harmonic loads. Here, an n-layered soil media system with a depth H is plotted in Fig. 2 to model the stratified ground. A harmonic load Pe^iωt with a radius of r₀ is applied at a certain depth of the soil media. The thickness of the i th layer is written as Δ H_i =H_i -H_i-1, where H_i and H_i-1 denote the depths from the surface of the i th and the (i-1) th layers to the surface of the soil media, respectively. The interfaces at the depth H_i (i=1, …, n-1) are assumed to be in perfect bonded contact. Before the analytical derivation of the solutions, several important boundary and continuity conditions should be stressed.

If we suppose that the surface of the multilayered soil media is free and permeable, the basic stress components and pore pressure will be

(29a)

(29b)

(29c)

For the problem in Fig. 2(a), the bottom of the soil media system is assumed to be fixed and impermeable. Then,

(30)

The adjacent continuity conditions are expressed as

(31a)

(31b)

where [V(ξ, H_j)_j and [U(ξ, H_j)_j are the general stress vector and the general displacement vector of the j th layer when the depth z=H_j, respectively. [V(ξ, H_j)_j+1 and [U(ξ, H_j)_j+1 are the general stress vector and displacement vector of the (j+1) th layer when the depth z=H_j, respectively. F(ξ, H_j) is the vertical harmonic load vector acting at a certain depth H_j of soil media in the Hankel transform domain. For a circular uniform loading, F(ξ, H_j)=. Here, J₁ (ξr) denotes the first order Bessel function of the first kind. If there is no force acting at the depth H_j, F(ξ, H_j)=[0, 0, 0]^T.

Taking into account the boundary and continuity conditions mentioned above, the relationship between the general stress vector and the general displacement vector connected by the global matrix for the multilayered soil can be established as

(32)

Similarly, for a half-space configuration, the global matrix can be assembled as follows:

(33)

Once the computation of general displacements is accomplished in the transform domain, solutions for multilayered soil media in the physical domain can be obtained with the inverse Hankel transform. To accomplish the computation of inverse Hankel transform, we compile the corresponding FORTRAN programs based on the Gauss-Legendre quadrature suggested by Ai. et al. ^[39], where the Gauss-Legendre quadrature of 32 points and 8 points is adopted for the first integrated interval and the other integral intervals of the Bessel function.

4 Verification

Numerical programs based on the derivation process are implemented to study the dynamic response of multilayered soil media subject to a harmonic excitation. In this section, a numerical example is presented to verify the accuracy of the proposed method. Before we start with the analyses of the numerical examples, some parameters should be defined firstly, i.e., a₀ =, and , in which f=nρ_fg/k. Moreover, with the purpose of reducing the influence of singularities on the precision of numerical integration, Lame parameters are complex in program computation ^[32-35], that is to say, λ₀ =λ (1+iδ₀), and G₀ =G(1+iδ₀), in which δ₀=0.03.

In order to confirm the correctness of the present method, comparisons of the numerical solutions for vertical displacements along the radius direction obtained from the present study and those given by Cai et al. ^[16] are made. The parameters in the example are r₀ = 1 m, a₀ =0.346, b₀ =42.5, P=10 kPa, and H=5 m. As shown in Fig. 3, the current solutions are in good agreement with the existing ones for the single soil layer with a subjacent rock-stratum subject to a harmonic vertical load on the surface.

5 Parameter study

For soil subsidence, there are various aspects of influencing factors in practical engineering, such as soil properties and load conditions. In this section, representative parameters are chosen to discuss their effects on dynamic responses of a multilayered saturated half-space due to the vertical excitation. Furthermore, numerical examples are performed to demonstrate effects of circular frequency, soil permeability, force depth, and soil stratification.

5.1 Circular frequency

A case is carried out here to investigate the influence of circular frequency ω. Solutions are presented for five different frequencies (ω =10 rad/s, 20 rad/s, 40 rad/s, 80 rad/s, and 120 rad/s). The parameters of poroelastic soil are G=λ =30 000 kPa, ρ =2 000 kg/m³, ρ_f =1 000 kg/m³, k=1× 10^-6 m/s, n=0.4, r₀ =1 m, P=10 kPa, H is a half-space, and the external vertical load is applied on the ground surface.

Variations of amplitude and real parts and imaginary parts of vertical displacements along the r-direction for different circular frequencies are plotted in Figs. 4-6, respectively. It is evident from these solutions that the response of the half-space depends significantly on the circular frequency of excitation. For each curve of a circular frequency ω, the amplitude of displacement along the r-direction shows a decreasing tendency. However, in the cases of ω = 80 rad/s and ω =120 rad/s, the amplitudes, real parts and imaginary parts of displacements show an apparent wave propagation phenomenon. As for the real parts of displacements, each case of a circular frequency ω≤40 rad/s shows a steady variation trend of decrease with the distance from the load center, while the imaginary parts of displacements show a steady variation trend of increase. This variation pattern is mainly related to the physical property of soil, showing that, when the circular frequency is close to the natural vibration frequency of soil, the propagation phenomenon is more obvious. It is observed that, at the same distance from the load center, the increase in the circular frequency does not lead to a monotonic increase in the vertical displacement. From Fig. 4, we can find that at the same position within a radius of 5 m, the amplitude of surface subsidence increases when the circular frequency increases from 10 rad/s to 20 rad/s, whereas the amplitude of displacement decreases when the circular frequency increases from 20 rad/s to 40 rad/s. Therefore, the circular frequency exerts a prominent effect on the soil deformation, and the reasonable numerical analysis should be performed for reference in the practical engineering.

5.2 Permeable coefficient

The parameters of poroelastic soil are G=λ =30 000 kPa, ρ =2 000 kg/m³, ρ_f = 1 000 kg/m³, ω =40 rad/s, n=0.4, r₀ =1 m, P=10 kPa, and H is a half-space. The load is applied on the ground surface. In order to investigate the influence of permeability on the displacement variation along the r-direction, six cases with various permeable coefficients are chosen. It can be found in Fig. 7 that, within a radius of 4 m near the load center, the case with a strong permeability presents a larger surface settlement value under the same loading condition. With the distance from the load center increasing, the influence of applied loads on the soil is weakened. Therefore, their differences are reduced. Moreover, when the permeable coefficient becomes small enough, such as the cases of k=10^-6 m/s and k=10^-7 m/s, the displacement variations tend to be the same. Figures 8 and 9 are the real and imaginary parts of vertical displacements along the r-direction for different permeable coefficients, and the similar rule can be concluded like Fig. 7.

5.3 Force depth

In this section, several cases of various force depths are studied. The basic parameters of poroelastic soil are G=λ =30 000 kPa, ρ=2 000 kg/m³, ρ_f=1 000 kg/m³, k=1× 10^-6 m/s, ω =40 rad/s, n=0.4, r₀ =1 m, P=10 kPa, and H is a half-space. Figures 10-12 show the amplitude, real parts and imaginary parts of displacements varying with z due to vertical harmonic excitations at four depths (h=0 m, 4 m, 10 m, and 14 m). We can observe that the displacement amplitude reaches the peak value at the depth of applied load (see Fig. 10). In cases of force depths h=4 m, 10 m, and 14 m, the displacements exhibit an approximately symmetrical variation on both sides of the peak value. The same varying pattern can also be found in Fig. 11. When the load is applied on the ground surface, the amplitude of displacement decreases first and then shows an upward increase with the depth. For the imaginary parts of displacements, in cases of force depths h =4 m, 10 m, and 14 m, displacements decrease before the load depth and increase after that depth. When the load is applied on the soil surface, the imaginary parts of displacements tend to increase with the depth.

5.4 Stratification characteristic

With the purpose of studying the influence of stratification characteristic on the soil subsidence, six cases are used to analyze displacements along the r-direction and the z-direction. Here, we choose three types of soil to independently combine and form six types of soil media. The basic parameters for the three types of soils are as follows: S1: G=1 500 kPa, λ = 7 875 kPa, ρ=1 800 kg/m³, k=10^-4 m/s; S2: G=4 500 kPa, λ =23 625 kPa, ρ =1 800 kg/m³, k=10^-5 m/s; S3: G=30 000 kPa, λ =30 000 kPa, ρ =2 000 kg/m³, k=10^-6 m/s. The investigated soil media consist of three layers, whose thickness is set as h₁ /r₀ =10, h₂ /r₀ =10, h₃ is a half-space. The arrangements of free combinations of three types of soils are as follows: Case 1: S1S2S3; Case 2: S2S1S3; Case 3: S2S3S1; Case 4: S1S3S2; Case 5: S3S1S2; Case 6: S3S2S1. Additionally, the dynamic load P = 10 kPa is applied on the ground surface.

Figures 13-15 show the displacement variation along the r-direction under various layer conditions. Through the results plotted in Fig. 13, the six cases are divided into three variation trends, which are separately Case 1 and Case 4, Case 2 and Case 3, Case 5 and Case 6. We denote these three variation trends as three groups. The similarity between the two cases in the same group is that their first layers are the same. Therefore, it is easy to observe that the property of the first layer is of great significance for soil subsidence compared with the other underlying layers. Based on the comparison of deformation of those three groups, soil media with a softer first layer show a much larger displacement near the load center. Additionally, as shown in Fig. 14, the order of real parts of displacement fluctuation is Group 1 > Group 2 > Group 3. That is to say, soil media with a softer first layer have a more obvious fluctuation variation along the r-direction than soil media with a stiffer first layer. The similar varying pattern can be found in Fig. 15.

Figures 16-18 show the displacement variation along the z-direction under various layer conditions. According to the results plotted in Fig. 16, displacements under different soil conditions show the similar variation trend to Fig. 13. In general, the amplitude of displacement in each case decreases with the depth. Meanwhile, curves of Group 1 and Group 2 show an obvious fluctuation near the ground surface. When the depth reaches 10 m, the deformation tends to be stable. By contrast, curves of Group 3 change smoothly and gently. That is because Group 3 consists of a stiffer first layer, and the load is applied on the ground surface.

5.5 Layer thickness

In order to investigate the influence of layer thickness on the soil subsidence, three cases are selected for comparisons of displacements along the r-direction. The soil types and parameters are consistent with the numerical examples in Subsection 5.4. According to their thicknesses, three cases for the research are chosen to carry out. Case Ⅰ: h₁ /r₀ =2, h₂ /r₀ =8, and h₃ is a half-space. Case Ⅱ: h₁ /r₀ =5, h₂ /r₀ =5, and h₃ is a half-space. Case Ⅲ: h₁ /r₀ =5, h₂ /r₀ =10, and h₃ is a half-space. Additionally, the dynamic load P = 10 kPa is applied on the ground surface. From Fig. 19, it is obvious that the thickness of soil layer also exerts a noticeable impact on the surface subsidence, especially near the load center. Figure 19 also reveals that the thickness of the first layer affects the surface subsidence more evidently than the underlying layers. Therefore, in the practical engineering, the first layer of soil foundation should be improved to a reasonable depth to satisfy the deformation requirement.

6 Conclusions

In this study, fundamental solutions for multilayered poroelastic soils subject to an axisymmetric harmonic vertical excitation are presented. The analytical layer-element method and Hankel transform are used in the derivation process. With the aid of boundary and adjacent continuity conditions of general stresses and displacements, the global matrix for multilayered soil media can be developed. Finally, solutions in the physical domain are obtained with the inverse Hankel transform. Selected numerical examples are carried out to verify the accuracy of the present method and demonstrate effects of soil parameters and load conditions on the dynamic responses of the multilayered poroelastic saturated soils. In view of the above investigations, it is reasonable to draw that the circular frequency exerts a prominent effect on the soil deformation. The influence of the soil permeability on the surface subsidence is more apparent near the load center. According to studies on the stratification effect, it is easy to observe that the property of the first layer is of great significance for the soil surface subsidence. The soil media with a softer first layer show a much larger displacement near the load center. Moreover, soil media with a softer first layer have a more obvious fluctuation variation along the r-direction than soil media with a stiffer first layer, and the force depth and the layer thickness also have an important impact on the dynamic response of saturated soils.

In conclusion, with the validated accuracy and efficiency, investigations of the present method and strategy can be used to solve the dynamic problem of stratified poroelastic soils, and it will provide a valuable theoretical guide for practical engineering in the emerging computational geotechnical field.

Appendix A

Elements of Φ for a finite layer are

where T=8εξ (Fγδεa₁a₃f₅ +ξa₂d₅h₁)-m₅θ₁²+m₆θ₂²+m₇θ₃²+m₈θ₄², a₁ =e^-ξz, a₂ =e^-εz, a₃ =e^-γz, b₁ =G(ε²+ξ²), b₂ =2Gξ²-δ, b₃ =2Gγ²-F+λ (γ²-ξ²), f₁ =1+a₁², f₂ =1+a₂², f₃ =1+a₃², f₄ =1-a₁², f₅ =1-a₂², f₆ =1-a₃², f₇ =a₁²+a₂², f₈ =a₁²+a₃², f₉ =a₂²+a₃², f₁₀ =a₁²-a₂², f₁₁ =a₁²-a₃², f₁₂ =a₂²-a₃², d₁ =γδ +Fξ, d₂ =γδ -Fξ, d₃ =F(ξ -ε)-δξ, d₄ =F(ξ +ε)-δξ, d₅ =F-δ, d₆ =F+δ, d₇ =γ +ξ, d₈ =γ -ξ, d₉ =γε +ξ², d₁₀ =γε -ξ², d₁₁ =ε +ξ, d₁₂ =ε -ξ, θ₁ =γδε +ξd₃, θ₂ =γδε -ξd₃, θ₃ =γδε +ξd₄, θ₄ =γδε -ξd₄, m₁ =1+a₁²a₂²a₃², m₂ =a₁²a₃²+a₂², m₃ =a₂²a₃²+a₁², m₄ =a₁²a₂²+a₃², m₅ =a₁²a₂²a₃²-1, m₆ =a₁²a₂²-a₃², m₇ =a₂²a₃²-a₁², m₈ =a₁²a₃²-a₂², m₉ =1+a₁²a₂², m₁₀ =1+a₂²a₃², m₁₁ =1+a₁²a₃², m₁₂ =a₁²a₂²-1, m₁₃ =a₂²a₃²-1, m₁₄ =a₁²a₃²-1, g₁ =-δd₁₀ +Fγd₁₂, g₂ =δd₉ -Fγd₁₂, g₃ =-δd₁₀ +Fγd₁₁, g₄ =δd₉ -Fγd₁₁, g₅ =d₉f₁₂ +d₁₀m₁₃, g₆ =d₇f₁₁ +d₈m₁₄, g₇ =d₉²f₉ -d₁₀²m₁₀, g₈ =-d₁₁²f₇ +d₁₂²m₉, g₉ =-d₇²f₈ +d₈²m₁₁, j₁ =2Gεd₂ +b₁d₅, j₂ =2Gεd₁ -b₁d₅, j₃ =2Gεd₁ +b₁d₅, j₄ =2Gεd₂ -b₁d₅, j₅ =Fb₂ +δb₃ -2Gξ²d₆, j₆ =(F-2δ)b₂ +δb₃ -2Gξ²d₅, j₇ =Fb₂ +(δ -2F)b₃ +2Gξ²d₅, h₁ =-γδa₃f₄ +Fξa₁f₆, h₂ =a₃f₁ -a₁f₃, h₃ =8Fγδξa₁a₃ -d₁²f₈ +d₂²m₁₁, h₄ =d₉f₉ +d₁₀m₁₀, h₅ =-d₁d₇f₈ +d₂d₈m₁₁, h₆ =2Gεξ²d₈ -b₂d₁₀ +ξb₃d₁₂, h₇ =2Gεξ²d₇ -b₂d₉ -ξb₃d₁₂, h₈ =2Gεξ²d₇ -b₂d₁₀ -ξb₃d₁₁, and h₉ =2Gεξ²d₈ -b₂d₉ +ξb₃d₁₁.

Appendix B

Elements of Φ₁ for the underlying half-space are

where b₁ =G(ε²+ξ²), b₂ =2Gξ²-δ, b₃ =2Gγ²-F+λ (γ²-ξ²), d₂ =γδ -Fξ, d₃ =F(ξ -ε)-δξ, d₆ =F+δ, d₈ =γ -ξ, d₁₀ =γε -ξ², d₁₂ =ε -ξ, j₁ =2Gεd₂ +b₁d₅, j₅ =Fb₂ +δb₃ -2Gξ²d₆, h₆ =2Gεξ²d₈ -b₂d₁₀ +ξb₃d₁₂, g₁ =-δd₁₀ +Fγd₁₂, and θ₁ =γδε +ξd₃.