Yue’s solution of classical elasticity in n-layered solids: Part 2, mathematical verification

Abstract

This paper presents a detailed and rigorous mathematical verification of Yue’s approach, Yue’s treatment, Yue’s method and Yue’s solution in the companion paper for the classical theory of elasticity in n-layered solid. It involves three levels of the mathematical verifications. The first level is to show that Yue’s solution can be automatically and uniformly degenerated into these classical solutions in closed-form such as Kelvin’s, Boussinesq’s, Mindlin’s and bi-material’s solutions when the material properties and boundary conditions are the same. This mathematical verification also gives and serves a clear and concrete understanding on the mathematical properties and singularities of Yue’s solution in n-layered solids. The second level is to analytically and rigorously show the convergence and singularity of the solution and the satisfaction of the solution to the governing partial differential equations, the interface conditions, the external boundary conditions and the body force loading conditions. This verification also provides the easy and executable means and results for the solutions in n-layered or graded solids to be calculated with any controlled accuracy in association with classical numerical integration techniques. The third level is to demonstrate the applicability and suitability of Yue’s approach, Yue’s treatment, Yue’s method and Yue’s solution to uniformly and systematically derive and formulate exact and complete solutions for other boundary-value problems, mixed-boundary value problems, and initial-boundary value problems in layered solids in the frameworks of classical elasticity, boundary element methods, elastodynamics, Biot’s theory of poroelasticity and thermoelasticity. All of such applications are substantiated by peer-reviewed journal publications made by the author and his collaborators over the past 30 years.

Keywords

elasticity / boundary element method / elastodynamics / poroelasticity / thermoelasticity

Introduction

In the companion paper [ 1], details of the Yue’s approach, Yue’s treatment, Yue’s method and Yue’s solution have been presented for the mathematical formulation of the solutions in n-layered solids in both transform and physical domains and in both Cartesian and cylindrical coordinate systems, where n is an arbitrary nonnegative integer. This paper presents a detailed and rigorous mathematical verification of the approach, the treatment, the method and the solution. To achieve this objective, the author has used the following three levels of the mathematical verifications.

Since Yue’s solution for the n-layered solid is a logical extension of these classical solutions in homogeneous solid, it shall be automatically degenerated into these classical solutions (or basic solutions) such as Kelvin’s, Boussinesq’s and Mindlin’s solutions. If the material properties and boundary conditions are the same, they shall be exactly the same. The fundamental singular solutions in exact closed-form are systematically and uniformly presented for the basic and classical boundary-value problems in either homogeneous or bi-homogeneous solids. This is the first level of the mathematical verification of Yue’s solutions in the n-layered solid. It shows that Yue’s approach gives all the classical solutions in closed-form by different authors with different mathematical tools and also many complete sets of new fundamental singular solutions in closed-form. In addition, the mathematical properties and singularities of the classical fundamental singular solutions in both the transform and physical domains are examined analytically and in closed-form. They give and serve a clear and concrete understanding on the mathematical properties and singularities of Yue’s solution in n-layered solids.

Secondly, Yue’s solution of the displacement vector \mathbf{u}(x,y,z) , the vertical stress vector {\mathbf{T}}_z(x,y,z) , and the plane strains {\mathbf{\Gamma }}_p(x,y,z) in the n-layered solid of either transverse isotropy or isotropy due to the internal loading concentrated on a horizontal plane ( \mathbf{f}(x,y,z)=\mathbf{f}(x,y)\delta (z-d) ), can be exactly and uniformly expressed as follows in the Cartesian coordinate systems.

(1a) \mathbf{u}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}\mathbf{\Pi }\mathbf{\Phi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(1b) {\mathbf{T}}_z(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }\mathbf{\Psi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(1c) {\mathbf{\Gamma }}_p(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}_p\mathbf{\Phi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

where , n is a positive integer and the body force vector \tilde{\mathbf{f}}(\xi ,\eta ) in the transform domain is expressed as follows.

(1d) \tilde{\mathbf{f}}\left(\xi ,\eta \right)={\displaystyle \frac{1}{2\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{f}\left(x,y\right)K^{*}\mathrm{d}x\mathrm{d}y}

The solution in Eqs. (1) is in the forms of improper integrals of infinite intervals either over the entire horizontal plane or along the radial axis . The improper integrals have many depending parameters including \mathrm{5}\times (n+\mathrm{2}) elastic constants (c_{ij},i=\mathrm{1},\mathrm{.2},\mathrm{3},\mathrm{4},\mathrm{5};j=\mathrm{0},\mathrm{1},\cdots ,n,n+\mathrm{1}) , n layer thicknesses (h_j,j=\mathrm{1},\mathrm{2},\cdots ,n) , the three independent variables (x,y,z) and the applied loading vectors \tilde{\mathbf{f}}(\xi ,\eta ) . Therefore, the second level of the mathematical verification of Yue’s solution is presented by analytically and rigorously examining the following three questions: a) the convergence of the solution; b) the singularity of the solution; c) satisfaction of the solution to the governing partial differential equations, the interface conditions, the external boundary conditions and the body force loading conditions. This verification also provides the easy and executable means and results for the solutions in n-layered or graded solids to be calculated with any controlled accuracy in association with classical numerical integration techniques. It is noted that the other three stresses {\mathbf{T}}_p and three strains {\mathbf{\Gamma }}_z can be obtained and examined using the Hooke’s law Eqs. (6) and the solution of {\mathbf{T}}_z and {\mathbf{\Gamma }}_p in Eqs. (1).

The third level of the mathematical verification of Yue’s solution is to demonstrate the applicability and suitability of Yue’s approach, Yue’s treatment and Yue’s method to uniformly and systematically derive and formulate exact and complete solutions for other boundary-value problems, mixed-boundary value problems, and initial-boundary value problems in the n-layered solids. Therefore, the applications to other problems in layered solids are further briefly presented in the frameworks of classical elasticity, boundary element methods, elastodynamics, Biot’s theory of poroelasticity and thermoelasticity. All of such applications are substantiated by peer-reviewed journal publications made by the author and his collaborators over the past 30 years. They show that the solutions for other problems can also be derived and formulated similarly and systematically and in the form of Yue’s matrix operations in both Cartesian and cylindrical coordinates.

Basic solutions in closed-form

General

To better understand and verify the solution derived above for the general layered solid, this section examines the basic solutions and their properties and singularities for some simplified material cases. They include

(a) The Kelvin’s case: one isotropic homogeneous solid of infinite extent subject to \mathbf{f}=f(x,y)\delta (z-d) ;

(b) The Boussinesq’s case: one isotropic homogeneous solid of semi-infinite extent subject to boundary loading \mathbf{f}=-f(x,y)\delta (z) ;

(c) The Mindlin’s case: one isotropic homogenous solid of semi-infinite extent subject to \mathbf{f}=f(x,y)\delta (z-d) ;

(d) The bi-material case: two perfectly bonded isotropic homogenous solid of infinite extent subject to \mathbf{f}=f(x,y)\delta (z-d) ;

(e) The extended Kelvin’s case: one transversely isotropic homogenous solid of infinite extent subject to \mathbf{f}=f(x,y)\delta (z-d) ;

(f) The extended bi-material case: two perfectly bonded transversely isotropic homogenous solids of infinite extent subject to \mathbf{f}=f(x,y)\delta (z-d) .

Since \mathbf{w}(\xi ,\eta ,z)=\mathbf{\Phi }(\rho ,z)\mathbf{g}(\xi ,\eta ) and {\mathbf{Y}}_z(\xi ,\eta ,z)=\mathbf{\Psi }(\rho ,z)\mathbf{g}(\xi ,\eta ) , \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) are only dependent on the solid materials and their occupied regions. They are derived and examined at first. The loading vector \mathbf{g}(\xi ,\eta ) is independent on the solid materials and is examined afterward for the closed-form solutions of the elastic field variables ( \mathbf{u} , {\mathbf{T}}_z and {\mathbf{\Gamma }}_p ) in physical domain.

Basic solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z)

The Kelvin’s case

In this case, the isotropic homogenous solid of infinite extent is subjected to the concentrated body force \mathbf{f}=f(x,y)\delta (z-d) . Its solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) is given below.

(2a) {\mathbf{\Phi }}_k(\rho ,z)=e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Phi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{k\mathrm{1}}\right)

(2b) {\mathbf{\Psi }}_k(\rho ,z)=e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Psi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{k\mathrm{1}}\right)

where , , and

(2c)| {\mathbf{\Phi }}_{k\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mu }}\left(\begin{array}{lll}\mathrm{1}+\alpha & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}+\alpha \end{array}\right);{\mathbf{\Phi }}_{k\mathrm{1}}={\displaystyle \frac{(\mathrm{1}-\alpha )}{\mathrm{4}\mu }}\left(\begin{array}{lll}-\mathrm{1}& \mathrm{0}& -{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}& \mathrm{0}& \mathrm{1}\end{array}\right)

(2d) {\mathbf{\Psi }}_{k\mathrm{0}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{lll}{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}& \mathrm{0}& -\alpha \\ \mathrm{0}& {\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}& \mathrm{0}\\ -\alpha & \mathrm{0}& {\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}\end{array}\right);{\mathbf{\Psi }}_{k\mathrm{1}}={\displaystyle \frac{(\mathrm{1}-\alpha )}{\mathrm{2}}}\left(\begin{array}{lll}{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{0}& -{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}\end{array}\right)

The Boussinesq’s case

In this case, the isotropic homogenous solid of semi-infinite extent is subjected to the concentrated body force \mathbf{f}=-f(x,y)\delta (z) on the external boundary z=\mathrm{0} . Its solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) is given below.

(3a) {\mathbf{\Phi }}_b(\rho ,z)=e^{-\rho z}({\mathbf{\Phi }}_{b\mathrm{0}}+\rho z{\mathbf{\Phi }}_{b\mathrm{1}})

(3b) {\mathbf{\Psi }}_b(\rho ,z)=e^{-\rho z}({\mathbf{\Psi }}_{b\mathrm{0}}+\rho z{\mathbf{\Psi }}_{b\mathrm{1}})

where , , and

(3c) {\mathbf{\Phi }}_{b\mathrm{0}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mu (\mathrm{1}-\alpha )}}\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \alpha \\ \mathrm{0}& \mathrm{2}(\mathrm{1}-\alpha )& \mathrm{0}\\ \alpha & \mathrm{0}& \mathrm{1}\end{array}\right);{\mathbf{\Phi }}_{b\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mu }}{\mathbf{I}}_{a\mathrm{1}};{\mathbf{\Psi }}_{b\mathrm{0}}={\mathbf{I}}_{\mathrm{3}};{\mathbf{\Psi }}_{b\mathrm{1}}=-{\mathbf{I}}_{a\mathrm{1}}

where

(3d) {\mathbf{I}}_{\mathrm{3}}=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right);{\mathbf{I}}_{a\mathrm{1}}=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{0}& -\mathrm{1}\end{array}\right)

The Mindlin’s case

In this case, the isotropic homogenous solid of semi-infinite extent is subjected to the concentrated body force \mathbf{f}=f(x,y)\delta (z-d) . Its solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) is given below.

(4a) {\mathbf{\Phi }}_m(\rho ,z)=e^{-\rho (z+d)}({\mathbf{\Phi }}_{m\mathrm{0}}+\rho z{\mathbf{\Phi }}_{m\mathrm{1}z}+\rho d{\mathbf{\Phi }}_{m\mathrm{1}d}+{\rho }^{\mathrm{2}}dz{\mathbf{\Phi }}_{m\mathrm{2}})+e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Phi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{k\mathrm{1}}\right)

(4b) {\mathbf{\Psi }}_m(\rho ,z)=e^{-\rho (z+d)}({\mathbf{\Psi }}_{m\mathrm{0}}+\rho z{\mathbf{\Psi }}_{m\mathrm{1}z}+\rho d{\mathbf{\Psi }}_{m\mathrm{1}d}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{m\mathrm{2}})+e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Psi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{k\mathrm{1}}\right)

where , , and

(4c) {\mathbf{\Phi }}_{m\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mu (\mathrm{1}-\alpha )}}\left(\begin{array}{lll}\mathrm{1}+{\alpha }^{\mathrm{2}}& \mathrm{0}& \mathrm{2}\alpha \\ \mathrm{0}& \mathrm{2}(\mathrm{1}-\alpha )& \mathrm{0}\\ \mathrm{2}\alpha & \mathrm{0}& \mathrm{1}+{\alpha }^{\mathrm{2}}\end{array}\right);{\mathbf{\Psi }}_{m\mathrm{0}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \alpha \\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \alpha & \mathrm{0}& \mathrm{1}\end{array}\right)

(4d) {\mathbf{\Phi }}_{m\mathrm{1}z}=-{\displaystyle \frac{(\mathrm{1}+\alpha )}{\mathrm{4}\mu }}{\mathbf{I}}_{a\mathrm{1}};{\mathbf{\Phi }}_{m\mathrm{1}d}=-{\displaystyle \frac{(\mathrm{1}+\alpha )}{\mathrm{4}\mu }}{\mathbf{I}}_{a\mathrm{2}};{\mathbf{\Phi }}_{m\mathrm{2}}={\displaystyle \frac{(\mathrm{1}-\alpha )}{\mathrm{2}\mu }}{\mathbf{I}}_{a\mathrm{3}}

(4e) {\mathbf{\Psi }}_{m\mathrm{1}z}={\displaystyle \frac{(\mathrm{1}+\alpha )}{\mathrm{2}}}{\mathbf{I}}_{a\mathrm{1}};{\mathbf{\Psi }}_{m\mathrm{1}d}={\displaystyle \frac{(\mathrm{1}-\alpha )}{\mathrm{2}}}{\mathbf{I}}_{a\mathrm{2}};{\mathbf{\Psi }}_{m\mathrm{2}}=(\alpha -\mathrm{1}){\mathbf{I}}_{a\mathrm{3}}

where

(4f) {\mathbf{I}}_{a\mathrm{2}}=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& -\mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& -\mathrm{1}\end{array}\right);{\mathbf{I}}_{a\mathrm{3}}=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& -\mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{0}& \mathrm{1}\end{array}\right)

It can be easily shown that the Eqs. (4) are reduced to those in Eqs. (3) if assuming d=\mathrm{0} . In other words, the solution for the Mindlin’s case includes that for the Boussinesq’s case.

The bi-material case

In this case, the two perfectly bonded isotropic homogenous solid of infinite extent is subjected to the concentrated body force \mathbf{f}=f(x,y)\delta (z-d) . The solid is perfectly bonded at z={\mathrm{0}}^{\pm } and. It is assumed that . Its solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) is given below.

In the solid of upper halfspace region with the two elastic constants {\mu }_{\mathrm{1}} and {\alpha }_{\mathrm{1}} , we have

(5a) {\mathbf{\Phi }}_t(\rho ,z)=e^{-\rho (z+d)}({\mathbf{\Phi }}_{t\mathrm{01}}+\rho z{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{1}}+\rho d{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{1}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Phi }}_{t\mathrm{21}})+e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Phi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{k\mathrm{1}}\right)

(5b) {\mathbf{\Psi }}_t(\rho ,z)=e^{-\rho (z+d)}({\mathbf{\Psi }}_{t\mathrm{01}}+\rho z{\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{1}}+\rho d{\mathbf{\Psi }}_{t\mathrm{0}z\mathrm{1}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{t\mathrm{21}})+e^{-\rho \left|\overline{z}\right|}\left({\mathbf{\Psi }}_{k\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{k\mathrm{1}}\right)

where , , {\mathbf{\Phi }}_{t\mathrm{01}} , {\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{1}} , {\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{1}} , {\mathbf{\Phi }}_{t\mathrm{21}} , {\mathbf{\Psi }}_{t\mathrm{01}} , {\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{1}} , {\mathbf{\Psi }}_{t\mathrm{1}z\mathrm{1}} and {\mathbf{\Psi }}_{t\mathrm{21}} are eight constant square matrices and defined in Appendix A.

In the solid of lower halfspace region with the two elastic constants {\mu }_{\mathrm{2}} and {\alpha }_{\mathrm{2}} , we have

(6a) {\mathbf{\Phi }}_t(\rho ,z)=e^{-\rho (d-z)}({\mathbf{\Phi }}_{t\mathrm{02}}+\rho z{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{2}}+\rho d{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{2}})

(6b) {\mathbf{\Psi }}_t(\rho ,z)=e^{-\rho (d-z)}({\mathbf{\Psi }}_{t\mathrm{02}}+\rho z{\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{2}}+\rho d{\mathbf{\Psi }}_{t\mathrm{1}z\mathrm{2}})

where , , {\mathbf{\Phi }}_{t\mathrm{02}} , {\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{2}} , {\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{2}} , {\mathbf{\Psi }}_{t\mathrm{02}} , {\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{2}} and {\mathbf{\Psi }}_{t\mathrm{1}z\mathrm{2}} are six constant square matrices and defined in Appendix A.

The extended Kelvin’s case

In this case, the transversely isotropic homogenous solid of infinite extent is subjected to \mathbf{f}=f(x,y)\delta (z-d) . Its five material constants are c_{j\mathrm{1}} , where j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5} . Its solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) is given below in Cases 1 and 2 for {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} , respectively, where {\mathrm{\Delta }}_{\mathrm{1}}=\sqrt{c_{\mathrm{11}}{c}_{\mathrm{31}}}-c_{\mathrm{21}}-\mathrm{2}{c}_{\mathrm{41}} .

Case 1: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0}

(7a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})

(7b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})

where , , and

(7c) {\mathrm{\Phi }}_v={\displaystyle \frac{1}{2c_{41}{\gamma }_{01}}}\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right);{\mathrm{\Phi }}_u\left(\chi \right)={\displaystyle \frac{1}{2\left({\gamma }_{11}^2-{\gamma }_{21}^2\right)\chi }}\left(\begin{array}{ccc}\frac{{\chi }^2}{c_{41}}-\frac{1}{c_{31}}& 0& \frac{\overline{z}}{\left|\overline{z}\right|}\frac{c_{21}+c_{41}}{c_{31}{c}_{41}}\chi \\ 0& 0& 0\\ -\frac{\overline{z}}{\left|\overline{z}\right|}\frac{c_{21}+c_{41}}{c_{31}{c}_{41}}\chi & 0& \frac{{\chi }^2}{c_{41}}-\frac{1}{c_{31}}\end{array}\right)

(7d) {\mathrm{\Psi }}_v={\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}0& 0& 0\\ 0& -\frac{\overline{z}}{\left|\overline{z}\right|}& 0\\ 0& 0& 0\end{array}\right);{\mathrm{\Psi }}_u\left(\chi \right)={\displaystyle \frac{-1}{2\left({\gamma }_{11}^2-{\gamma }_{21}^2\right)\chi }}\left(\begin{array}{ccc}\frac{\overline{z}}{\left|\overline{z}\right|}\left({\chi }^3+\frac{c_{21}}{c_{31}}\chi \right)& 0& \frac{c_{21}}{c_{31}}{\chi }^2+\frac{c_{11}}{c_{31}}\\ 0& 0& 0\\ -{\chi }^2-\frac{c_{21}}{c_{31}}& 0& \frac{\overline{z}}{\left|\overline{z}\right|}\left({\chi }^3+c_{q1}\chi \right)\end{array}\right)

Case 2: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0}

(8a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Phi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{u\mathrm{1}}\right)

(8b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Psi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{u\mathrm{1}}\right)

where , , and

(8c) {\mathbf{\Phi }}_{u\mathrm{0}}={\displaystyle \frac{c_{\mathrm{21}}+\mathrm{3}{c}_{\mathrm{41}}}{\mathrm{4}{c}_{\mathrm{31}}{c}_{\mathrm{41}}{\gamma }_{\mathrm{31}}^{\mathrm{3}}}}\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\gamma }_{\mathrm{31}}^{\mathrm{2}}\end{array}\right);{\mathbf{\Phi }}_{u\mathrm{1}}={\displaystyle \frac{c_{\mathrm{21}}+c_{\mathrm{41}}}{\mathrm{4}{c}_{\mathrm{31}}{c}_{\mathrm{41}}{\gamma }_{\mathrm{31}}^{\mathrm{3}}}}\left(\begin{array}{ccc}-\mathrm{1}& \mathrm{0}& -{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}{\gamma }_{\mathrm{31}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}{\gamma }_{\mathrm{31}}& \mathrm{0}& {\gamma }_{\mathrm{31}}^{\mathrm{2}}\end{array}\right)

(8d) {\mathbf{\Psi }}_{u\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{ccc}-{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}& \mathrm{0}& {\displaystyle \frac{c_{\mathrm{41}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{31}}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\displaystyle \frac{c_{\mathrm{41}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{31}}^{\mathrm{3}}}}& \mathrm{0}& -{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}\end{array}\right);{\mathbf{\Psi }}_{u\mathrm{1}}={\displaystyle \frac{c_{\mathrm{21}}+c_{\mathrm{41}}}{\mathrm{2}{c}_{\mathrm{31}}{\gamma }_{\mathrm{31}}^{\mathrm{2}}}}\left(\begin{array}{ccc}{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}{\gamma }_{\mathrm{31}}& \mathrm{0}& {\gamma }_{\mathrm{31}}^{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{0}& -{\displaystyle \frac{\overline{z}}{\left|\overline{z}\right|}}{\gamma }_{\mathrm{31}}\end{array}\right)

It is noted that the following limits are valid.

(9a) \underset{{\gamma }_{\mathrm{11}}\to {\gamma }_{\mathrm{21}}}{\lim }\left(e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)=e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Phi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{u\mathrm{1}}\right)

(9b) \underset{{\gamma }_{\mathrm{11}}\to {\gamma }_{\mathrm{21}}}{\lim }\left(e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})\right)=e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Psi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{u\mathrm{1}}\right)

It can also be shown that the solution for Case 2 ( {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} ) in Eqs. (8) can be degenerated to the solution for the Kelvin’s case in Eqs. (2) if the transversely isotropic material becomes to the isotropic material.

The extended bi-material case

In this case, the two perfectly bonded transversely isotropic homogenous solids of infinite extent is subjected to \mathbf{f}=f(x,y)\delta (z-d) . The two solids are perfectly bonded at z={\mathrm{0}}^{\pm } . The upper halfspace region is occupied by the solid 1 with the five material constants are c_{j\mathrm{1}} , where j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5} . The lower halfspace region is occupied by the solid 2 with the five material constants c_{j\mathrm{2}} , where j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5} . The loading plane is located within the solid 1, where . There are four combinations of {\mathrm{\Delta }}_{\mathrm{1}} and {\mathrm{\Delta }}_{\mathrm{2}} as follows: a) {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0} ; b) {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0} ; c) {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0} ; d) {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0} , where {\mathrm{\Delta }}_{\mathrm{2}}=\sqrt{c_{\mathrm{21}}{c}_{\mathrm{32}}}-c_{\mathrm{22}}-\mathrm{2}{c}_{\mathrm{42}} . The solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) for each of the four cases is specifically given below.

Case a: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

In the solid of upper halfspace region ,

(10a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Phi }}_{aJ\mathrm{1}}}+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})

(10b) {\mathrm{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathrm{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathrm{\Psi }}_{aJ\mathrm{1}}}+e^{-\rho {\gamma }_{\mathrm{01}}\left|\bar{z}\right|}{\mathrm{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\bar{z}\right|}{\mathrm{\Psi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\bar{z}\right|}{\mathrm{\Psi }}_u({\gamma }_{\mathrm{21}})

where z_{\mathrm{01}}={\gamma }_{\mathrm{01}}(d+z) , z_{a\mathrm{11}}={\gamma }_{\mathrm{11}}(d+z) , z_{a\mathrm{21}}={\gamma }_{\mathrm{11}}d+{\gamma }_{\mathrm{21}}z , z_{a\mathrm{31}}={\gamma }_{\mathrm{21}}d+{\gamma }_{\mathrm{11}}z , z_{a\mathrm{41}}={\gamma }_{\mathrm{21}}(d+z) ; The constant square matrices {\mathbf{\Phi }}_{\mathrm{01}} , {\mathbf{\Psi }}_{\mathrm{01}} , {\mathbf{\Phi }}_{aJ\mathrm{1}} and {\mathbf{\Psi }}_{aJ\mathrm{1}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in Refs. [ 2, 3].

In the solid of lower halfspace region ,

(11a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Phi }}_{\mathrm{02}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{2}}}{\mathbf{\Phi }}_{aJ\mathrm{2}}}

(11b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Psi }}_{\mathrm{02}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{2}}}{\mathbf{\Psi }}_{aJ\mathrm{2}}}

where z_{\mathrm{02}}={\gamma }_{\mathrm{01}}d-{\gamma }_{\mathrm{02}}z , z_{a\mathrm{12}}={\gamma }_{\mathrm{11}}d-{\gamma }_{\mathrm{12}}z , z_{a\mathrm{22}}={\gamma }_{\mathrm{11}}d-{\gamma }_{\mathrm{22}}z , z_{a\mathrm{32}}={\gamma }_{\mathrm{21}}d-{\gamma }_{\mathrm{12}}z , z_{a\mathrm{42}}={\gamma }_{\mathrm{21}}d-{\gamma }_{\mathrm{22}}z ; The constant square matrices {\mathbf{\Phi }}_{\mathrm{02}} , {\mathbf{\Psi }}_{\mathrm{02}} , {\mathbf{\Phi }}_{aJ\mathrm{2}} and {\mathbf{\Psi }}_{aJ\mathrm{2}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in Refs. [ 2, 3].

Case b: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

In the solid of upper halfspace region ,

(12a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Phi }}_{b\mathrm{11}}+\rho z{\mathbf{\Phi }}_{b\mathrm{21}}+\rho d{\mathbf{\Phi }}_{b\mathrm{31}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Phi }}_{b\mathrm{41}})+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Phi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{u\mathrm{1}}\right)

(12b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Psi }}_{b\mathrm{11}}+\rho z{\mathbf{\Psi }}_{b\mathrm{21}}+\rho d{\mathbf{\Psi }}_{b\mathrm{31}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{b\mathrm{41}})+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Psi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{u\mathrm{1}}\right)

where z_{b\mathrm{1}}={\gamma }_{\mathrm{31}}(d+z) , {\mathbf{\Phi }}_{bJ\mathrm{1}} and {\mathbf{\Psi }}_{bJ\mathrm{1}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in Refs. [ 2, 3].

In the solid of lower halfspace region ,

(13a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Phi }}_{\mathrm{02}}+e^{-\rho z_{b\mathrm{12}}}({\mathbf{\Phi }}_{b\mathrm{12}}+\rho d{\mathbf{\Phi }}_{b\mathrm{22}})+e^{-\rho z_{b\mathrm{22}}}({\mathbf{\Phi }}_{b\mathrm{32}}+\rho d{\mathbf{\Phi }}_{b\mathrm{42}})

(13b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Psi }}_{\mathrm{02}}+e^{-\rho z_{b\mathrm{12}}}({\mathbf{\Psi }}_{b\mathrm{12}}+\rho d{\mathbf{\Psi }}_{b\mathrm{22}})+e^{-\rho z_{b\mathrm{22}}}({\mathbf{\Psi }}_{b\mathrm{32}}+\rho d{\mathbf{\Psi }}_{b\mathrm{42}})

where z_{b\mathrm{12}}={\gamma }_{\mathrm{31}}d-{\gamma }_{\mathrm{12}}z , z_{b\mathrm{22}}={\gamma }_{\mathrm{31}}d-{\gamma }_{\mathrm{22}}z , {\mathbf{\Phi }}_{bJ\mathrm{2}} and {\mathbf{\Psi }}_{bJ\mathrm{2}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in Refs. [ 2, 3].

Case c: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

In the solid of upper halfspace region

(14a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Phi }}_{cJ\mathrm{1}}}+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})

(14b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Psi }}_{cJ\mathrm{1}}}+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{11}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})-e^{-\rho {\gamma }_{\mathrm{21}}\left|\overline{z}\right|}{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})

where z_{a\mathrm{11}}={\gamma }_{\mathrm{11}}(d+z) , z_{a\mathrm{21}}={\gamma }_{\mathrm{11}}d+{\gamma }_{\mathrm{21}}z , z_{a\mathrm{13}}={\gamma }_{\mathrm{21}}d+{\gamma }_{\mathrm{11}}z , z_{a\mathrm{14}}={\gamma }_{\mathrm{21}}(d+z) ; {\mathbf{\Phi }}_{cJ\mathrm{1}} and {\mathbf{\Psi }}_{cJ\mathrm{1}} are given in Refs. [ 2, 3].

In the solid of lower halfspace region ,

(15a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Phi }}_{\mathrm{02}}+e^{-\rho z_{c\mathrm{12}}}({\mathbf{\Phi }}_{c\mathrm{12}}+\rho z{\mathbf{\Phi }}_{c\mathrm{22}})+e^{-\rho z_{c\mathrm{22}}}({\mathbf{\Phi }}_{c\mathrm{32}}+\rho z{\mathbf{\Phi }}_{c\mathrm{42}})

(15b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Psi }}_{\mathrm{02}}+e^{-\rho z_{c\mathrm{12}}}({\mathbf{\Psi }}_{c\mathrm{12}}+\rho z{\mathbf{\Psi }}_{c\mathrm{22}})+e^{-\rho z_{c\mathrm{22}}}({\mathbf{\Psi }}_{c\mathrm{32}}+\rho z{\mathbf{\Psi }}_{c\mathrm{42}})

where z_{c\mathrm{12}}={\gamma }_{\mathrm{11}}d-{\gamma }_{\mathrm{32}}z , z_{c\mathrm{22}}={\gamma }_{\mathrm{21}}d-{\gamma }_{\mathrm{32}}z , {\mathbf{\Phi }}_{cJ\mathrm{2}} and {\mathbf{\Psi }}_{cJ\mathrm{2}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in Refs. [ 2, 3].

Case d: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

In the solid of upper halfspace region

(16a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Phi }}_{d\mathrm{11}}+\rho z{\mathbf{\Phi }}_{d\mathrm{21}}+\rho d{\mathbf{\Phi }}_{d\mathrm{31}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Phi }}_{d\mathrm{41}})+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Phi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Phi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Phi }}_{u\mathrm{1}}\right)

(16b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Psi }}_{d\mathrm{11}}+\rho z{\mathbf{\Psi }}_{d\mathrm{21}}+\rho d{\mathbf{\Psi }}_{d\mathrm{31}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{d\mathrm{41}})+e^{-\rho {\gamma }_{\mathrm{01}}\left|\overline{z}\right|}{\mathbf{\Psi }}_v+e^{-\rho {\gamma }_{\mathrm{31}}\left|\overline{z}\right|}\left({\mathbf{\Psi }}_{u\mathrm{0}}+\rho \left|\overline{z}\right|{\mathbf{\Psi }}_{u\mathrm{1}}\right)

where z_{b\mathrm{1}}={\gamma }_{\mathrm{31}}(d+z) , {\mathbf{\Phi }}_{dJ\mathrm{1}} and {\mathbf{\Psi }}_{dJ\mathrm{1}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in [ 2, 3].

In the solid of lower halfspace region ,

(17a) {\mathbf{\Phi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Phi }}_{\mathrm{02}}+e^{-\rho z_{d\mathrm{2}}}({\mathbf{\Phi }}_{d\mathrm{12}}+\rho z{\mathbf{\Phi }}_{d\mathrm{22}}+\rho d{\mathbf{\Phi }}_{d\mathrm{32}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Phi }}_{d\mathrm{42}})

(17b) {\mathbf{\Psi }}_{kt}(\rho ,z)=e^{-\rho z_{\mathrm{02}}}{\mathbf{\Psi }}_{\mathrm{02}}+e^{-\rho z_{d\mathrm{2}}}({\mathbf{\Psi }}_{d\mathrm{12}}+\rho z{\mathbf{\Psi }}_{d\mathrm{22}}+\rho d{\mathbf{\Psi }}_{d\mathrm{32}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{d\mathrm{42}})

where z_{d\mathrm{2}}={\gamma }_{\mathrm{31}}d-{\gamma }_{\mathrm{32}}z , {\mathbf{\Phi }}_{dJ\mathrm{2}} and {\mathbf{\Psi }}_{dJ\mathrm{2}} ( J=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} ) are given in [ 2, 3].

Constant limits of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) at loading plane

The Kelvin’s case

For this case, the solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) has the following constant limits at loading plane z=d .

(18a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_k(\rho ,z)={\mathbf{\Phi }}_{k\mathrm{0}}

(18b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_k(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{\mathrm{0}}

where , and

(18c) {\mathbf{\Phi }}_{k\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mu }}\left(\begin{array}{ccc}\mathrm{1}+\alpha & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}+\alpha \end{array}\right);{\mathbf{\Psi }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{lll}\mathrm{0}& \mathrm{0}& \alpha \\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \alpha & \mathrm{0}& \mathrm{0}\end{array}\right);{\mathbf{I}}_{\mathrm{3}}=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)

The Boussinesq’s case

For this case, the solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) has the following constant limits at loading plane z=\mathrm{0} .

(19a) \underset{z\to {\mathrm{0}}^{+}}{\lim }{\mathbf{\Phi }}_b(\rho ,z)={\mathbf{\Phi }}_{b\mathrm{0}}

(19b) \underset{z\to {\mathrm{0}}^{+}}{\lim }{\mathbf{\Psi }}_b(\rho ,z)={\mathbf{I}}_{\mathrm{3}}

where .

The Mindlin’s case

For this case, the solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) has the following constant limits at loading plane z=d , where d\ge \mathrm{0} .

(20a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_m(\rho ,z)={\mathbf{\Phi }}_{k\mathrm{0}}+e^{-\mathrm{2}\rho d}({\mathbf{\Phi }}_{m\mathrm{0}}+\rho d{\mathbf{\Phi }}_{m\mathrm{1}z}+\rho d{\mathbf{\Phi }}_{m\mathrm{1}d}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Phi }}_{m\mathrm{2}})

(20b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_m(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{\mathrm{0}}+e^{-\mathrm{2}\rho d}({\mathbf{\Psi }}_{m\mathrm{0}}+\rho d{\mathbf{\Psi }}_{m\mathrm{1}z}+\rho d{\mathbf{\Psi }}_{m\mathrm{1}d}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Psi }}_{m\mathrm{2}})

where .

In particular, if d=\mathrm{0} , we have, \underset{z\to d^{+}}{\lim }{\mathbf{\Phi }}_m(\rho ,z)={\mathbf{\Phi }}_{b\mathrm{0}} and \underset{z\to d^{+}}{\lim }{\mathbf{\Psi }}_m(\rho ,z)=-{\mathbf{I}}_{\mathrm{3}} , which is the Boussinesq’s case.

The bi-material case

For this case, the solution of \mathrm{\Phi }(\rho ,z) and \mathrm{\Psi }(\rho ,z) has the following constant limits at loading plane z=d , where the two solids are bonded at z={\mathrm{0}}^{\pm } and .

(21a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_t(\rho ,z)={\mathbf{\Phi }}_{k\mathrm{0}}+e^{-\mathrm{2}\rho d}({\mathbf{\Phi }}_{t\mathrm{01}}+\rho d{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{1}}+\rho d{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{1}}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Phi }}_{t\mathrm{21}})

(21b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_t(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{\mathrm{0}}+e^{-\mathrm{2}\rho d}({\mathbf{\Psi }}_{t\mathrm{01}}+\rho d{\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{1}}+\rho d{\mathbf{\Psi }}_{t\mathrm{0}z\mathrm{1}}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Psi }}_{t\mathrm{21}})

where .

In particular, if d=\mathrm{0} , we have \underset{z\to {\mathrm{0}}^{+}}{\lim }{\mathbf{\Phi }}_t(\rho ,z)={\mathbf{\Phi }}_{t\mathrm{01}}+{\mathbf{\Phi }}_{k\mathrm{0}}={\mathbf{\Phi }}_{t\mathrm{02}} and \underset{z\to {\mathrm{0}}^{+}}{\lim }{\mathbf{\Psi }}_t(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{\mathrm{0}}+{\mathbf{\Psi }}_{t\mathrm{01}} .

The extended Kelvin’s case

For this case, the solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) has the following cases of constant limits at loading plane z=d .

Case 1: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0}

(22a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{1}}

(22b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{1}}

where , and

(22c) {\mathbf{\Phi }}_{L\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{c}_{\mathrm{41}}(\gamma +_{\mathrm{11}}{\gamma }_{\mathrm{21}})}}\left(\begin{array}{ccc}\mathrm{1}+{\displaystyle \frac{c_{\mathrm{41}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{11}}{\gamma }_{\mathrm{21}}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\gamma +_{\mathrm{11}}{\gamma }_{\mathrm{21}}}{{\gamma }_{\mathrm{01}}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}+{\displaystyle \frac{c_{\mathrm{41}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{11}}{\gamma }_{\mathrm{21}}}}\end{array}\right)

(22d) {\mathbf{\Psi }}_{L\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}({\gamma }_{\mathrm{11}}+{\gamma }_{\mathrm{21}})}}\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& {\displaystyle \frac{c_{\mathrm{11}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{11}}{\gamma }_{\mathrm{21}}}}-{\displaystyle \frac{c_{\mathrm{21}}}{c_{\mathrm{31}}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}-{\displaystyle \frac{c_{\mathrm{21}}}{c_{\mathrm{31}}{\gamma }_{\mathrm{11}}{\gamma }_{\mathrm{21}}}}& \mathrm{0}& \mathrm{0}\end{array}\right)

Case 2: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0}

(23a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{2}}

(23b)| \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{2}}

where and

(23c) {\mathbf{\Phi }}_{L\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{c}_{\mathrm{41}}}}\left(\begin{array}{ccc}{\displaystyle \frac{c_{\mathrm{21}}+\mathrm{3}{c}_{\mathrm{41}}}{\mathrm{2}{c}_{\mathrm{31}}{\gamma }_{\mathrm{31}}^{\mathrm{3}}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{{\gamma }_{\mathrm{01}}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{c_{\mathrm{21}}+\mathrm{3}{c}_{\mathrm{41}}}{\mathrm{2}{c}_{\mathrm{31}}{\gamma }_{\mathrm{31}}}}\end{array}\right)

(23d) {\mathbf{\Psi }}_{L\mathrm{2}}={\displaystyle \frac{c_{\mathrm{41}}}{\mathrm{2}{c}_{\mathrm{31}}{\gamma }_{\mathrm{31}}^{\mathrm{3}}}}\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& {\gamma }_{\mathrm{31}}^{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right)

It is noted that \underset{{\gamma }_{\mathrm{11}}\to {\gamma }_{\mathrm{21}}}{\lim }{\mathbf{\Phi }}_{L\mathrm{1}}={\mathbf{\Phi }}_{L\mathrm{2}} and \underset{{\gamma }_{\mathrm{11}}\to {\gamma }_{\mathrm{21}}}{\lim }{\mathbf{\Psi }}_{L\mathrm{1}}={\mathbf{\Psi }}_{L\mathrm{2}} for {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} .

The extended bi-material case

For this case, the solution of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) has the following constant limits at loading plane z=d , where the two solids are bonded at z={\mathrm{0}}^{\pm } and .

Case a: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

(24a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{1}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Phi }}_{aJ\mathrm{1}}}

(24b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{1}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Psi }}_{aJ\mathrm{1}}}

where z_{\mathrm{01}}=\mathrm{2}{\gamma }_{\mathrm{01}}d , z_{a\mathrm{11}}=\mathrm{2}{\gamma }_{\mathrm{11}}d , z_{a\mathrm{21}}={\gamma }_{\mathrm{11}}d+{\gamma }_{\mathrm{21}}d , z_{a\mathrm{31}}={\gamma }_{\mathrm{21}}d+{\gamma }_{\mathrm{11}}d , and z_{a\mathrm{41}}=\mathrm{2}{\gamma }_{\mathrm{21}}d.

Case b: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

(25a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{2}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Phi }}_{b\mathrm{11}}+\rho d{\mathbf{\Phi }}_{b\mathrm{21}}+\rho d{\mathbf{\Phi }}_{b\mathrm{31}}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Phi }}_{b\mathrm{41}})

(25b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{2}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Psi }}_{b\mathrm{11}}+\rho d{\mathbf{\Psi }}_{b\mathrm{21}}+\rho d{\mathbf{\Psi }}_{b\mathrm{31}}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Psi }}_{b\mathrm{41}})

where z_{b\mathrm{1}}=\mathrm{2}{\gamma }_{\mathrm{31}}d

Case c: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

(26a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{1}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Phi }}_{cJ\mathrm{1}}}

(26b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{1}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{e}^{-\rho z_{aJ\mathrm{1}}}{\mathbf{\Psi }}_{cJ\mathrm{1}}}

where z_{a\mathrm{11}}=\mathrm{2}{\gamma }_{\mathrm{11}}d , z_{a\mathrm{21}}={\gamma }_{\mathrm{11}}d+{\gamma }_{\mathrm{21}}d , z_{a\mathrm{13}}={\gamma }_{\mathrm{21}}d+{\gamma }_{\mathrm{11}}d , and z_{a\mathrm{14}}=\mathrm{2}{\gamma }_{\mathrm{21}}d .

Case d: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

(27a) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{2}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Phi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Phi }}_{d\mathrm{11}}+\rho d{\mathbf{\Phi }}_{d\mathrm{21}}+\rho d{\mathbf{\Phi }}_{d\mathrm{31}}+{\rho }^{\mathrm{2}}{d}^{\mathrm{2}}{\mathbf{\Phi }}_{d\mathrm{41}})

(27b) \underset{z\to d^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{2}}+e^{-\rho z_{\mathrm{01}}}{\mathbf{\Psi }}_{\mathrm{01}}+e^{-\rho z_{b\mathrm{1}}}({\mathbf{\Psi }}_{d\mathrm{11}}+\rho z{\mathbf{\Psi }}_{d\mathrm{21}}+\rho d{\mathbf{\Psi }}_{d\mathrm{31}}+{\rho }^{\mathrm{2}}zd{\mathbf{\Psi }}_{d\mathrm{41}})

where z_{b\mathrm{1}}=\mathrm{2}{\gamma }_{\mathrm{31}}d .

In particular, if d=\mathrm{0} , we have

Case a: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

(28a) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{1}}={\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{\Phi }}_{aJ\mathrm{1}}}

(28b) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{1}}+{\mathbf{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{\Psi }}_{aJ\mathrm{1}}}

Case b: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

(29a) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{2}}+{\mathbf{\Phi }}_{\mathrm{01}}+{\mathbf{\Phi }}_{b\mathrm{11}}

(29b) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{2}}+{\mathbf{\Psi }}_{\mathrm{01}}+{\mathbf{\Psi }}_{b\mathrm{11}}

Case c: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

(30a) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{1}}+{\mathbf{\Phi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{\Phi }}_{cJ\mathrm{1}}}

(30b) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{1}}+{\mathbf{\Psi }}_{\mathrm{01}}+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{\Psi }}_{cJ\mathrm{1}}}

Case d: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

(31a) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Phi }}_{kt}(\rho ,z)={\mathbf{\Phi }}_{L\mathrm{2}}+{\mathbf{\Phi }}_{\mathrm{01}}+{\mathbf{\Phi }}_{d\mathrm{11}}

(31b) \underset{z\to {\mathrm{0}}^{\pm }}{\lim }{\mathbf{\Psi }}_{kt}(\rho ,z)=\mp {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{I}}_{\mathrm{3}}+{\mathbf{\Psi }}_{L\mathrm{2}}+{\mathbf{\Psi }}_{\mathrm{01}}+{\mathbf{\Psi }}_{d\mathrm{11}}

where .

The properties of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z)

It is evident that the above specific solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) have only the exponential functions in the form of (a_{\mathrm{0}}+a_{\mathrm{1}}d+a_{\mathrm{2}}\rho z+a_{\mathrm{3}}{\rho }^{\mathrm{2}}dz)e^{-a\rho z} (where a,a_{\mathrm{0}},a_{\mathrm{1}},a_{\mathrm{2}},a_{\mathrm{3}} are material constants (or zero) and \mathrm{Re}(a)>\mathrm{0} ). Their mathematical properties can be listed as follows.

All of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) are continuous functions of \rho , where .

\mathbf{\Phi }(\rho ,z) in Eq. (2a) for the Kelvin’s case and in Eq. (7a) for the extended Kelvin’s case are continuous functions of z , where . \mathbf{\Phi }(\rho ,z) in Eq. (3a) for the Boussinesq’s case and in Eq. (4a) for the Mindlin’s case are continuous functions of z , where . \mathbf{\Phi }(\rho ,z) in Eqs. (5a) and (6a) for the bi-material case and in Eqs. (10a)−(17a) for the extended bi-material case are continuous functions of z for .

\mathbf{\Psi }(\rho ,z) in Eq. (2b) for the Kelvin’s case and in Eq. (7b) for the extended Kelvin’s case are continuous functions of z for either or . \mathbf{\Psi }(\rho ,z) in Eq. (3b) for the Boussinesq’s case are continuous functions of z for . \mathbf{\Psi }(\rho ,z) in Eq. (4b) for the Mindlin’s case are continuous functions of z for either {\mathrm{0}}^{+}\le z\le d^{-} or . \mathbf{\Psi }(\rho ,z) in Eqs. (5b) and (6b) for the bi-material case and in Eqs. (10a)–(17b) for the Eq. extended bi-material case are continuous functions of z for either or {\mathrm{0}}^{+}\le z\le d^{-} or .

At the loading plane z=d , \mathbf{\Psi }(\rho ,z) has the unit step discontinuity as follows.

(32) \mathbf{\Psi }(\rho ,d^{+})-\mathbf{\Psi }(\rho ,d^{-})=-{\mathbf{I}}_{\mathrm{3}}

The partial differentiations of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) with respect to z are continuous functions for any z and \rho within each region of the homogeneous solids, except at the loading plane z=d or at the interface z=\mathrm{0} of the two dissimilar solids.

As \rho \to +\infty and z\ne d , \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) converge to zero with the rate of e^{-\rho \omega \left|z-d\right|} , where \omega >\mathrm{0} . Specifically, \omega =\mathrm{1} for isotropic solid; \omega ={\gamma }_{\mathrm{01}} , \omega ={\gamma }_{\mathrm{21}} for {\mathrm{\Delta }}_{\mathrm{1}}>\mathrm{0} ; \omega =c_{a\mathrm{1}} for ; and \omega ={\gamma }_{\mathrm{31}} for {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} . At the loading plane z=d , \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) have the constant limits for . These constant limits are given in Eqs. (18a)–(31a) for \mathbf{\Phi }(\rho ,d) and Eqs. (28b)–(31b) for \mathbf{\Psi }(\rho ,d^{\pm }) . They are the singularities of the solutions in the physical domain, which are discussed in the ensuing.

The Green’s function matrices

Definitions in Cartesian coordinate system

At first, ten basic functions q_{km} are defined as follows.

(33a) q_{\mathrm{00}}=\mathrm{1}

(33b) q_{\mathrm{10}}={\displaystyle \frac{i\xi }{\rho }};\hspace{1em}{q}_{\mathrm{01}}={\displaystyle \frac{i\eta }{\rho }}

(33c) q_{\mathrm{20}}={\displaystyle \frac{{\xi }^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}};\hspace{1em}{q}_{\mathrm{11}}={\displaystyle \frac{\xi \eta }{{\rho }^{\mathrm{2}}}};\hspace{1em}{q}_{\mathrm{02}}={\displaystyle \frac{{\eta }^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}}

(33d) q_{\mathrm{30}}={\displaystyle \frac{i{\xi }^{\mathrm{3}}}{{\rho }^{\mathrm{3}}}};\hspace{1em}{q}_{\mathrm{21}}={\displaystyle \frac{i{\xi }^{\mathrm{2}}\eta }{{\rho }^{\mathrm{3}}}};\hspace{1em}{q}_{\mathrm{12}}={\displaystyle \frac{i\xi {\eta }^{\mathrm{2}}}{{\rho }^{\mathrm{3}}}};\hspace{1em}{q}_{\mathrm{03}}={\displaystyle \frac{i{\eta }^{\mathrm{3}}}{{\rho }^{\mathrm{3}}}}

Secondly, the following harmonic functions g_{LIJ} are defined,

(34) g_{LIJ}=g_{LIJ}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{L-\mathrm{1}}{e}^{-\rho z}{q}_{IJ}\tilde{f}K\mathrm{d}\xi \mathrm{d}\eta }

where L, I and J are nonnegative integers, \mathrm{0}\le I+J\le \mathrm{3} ; for L=\mathrm{0} \mathrm{or} \mathrm{1} , z\ge \mathrm{0} ; for L\ge \mathrm{2} , z>\mathrm{0} ; and K=e^{i(x\xi +y\eta )} .

It can be shown that the harmonic functions g_{LIJ}\left(x,y,z,\tilde{f}\right) satisfy the Laplace equations as follows:

(35) {\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}{g}_{LIJ}\left(x,y,z,\tilde{f}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}{g}_{LIJ}\left(x,y,z,\tilde{f}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}{g}_{LIJ}\left(x,y,z,\tilde{f}\right)=\mathrm{0.}

They also have the following recursive relations.

(36a) g_{(L+\mathrm{1})IJ}\left(x,y,z,\tilde{f}\right)=-{\displaystyle \frac{\partial }{\partial z}}{g}_{LIJ}\left(x,y,z,\tilde{f}\right)

(36b) g_{\mathrm{121}}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{\partial }{\partial y}}{g}_{\mathrm{020}}\left(x,y,z,\tilde{f}\right)

(36c) g_{\mathrm{103}}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{\partial }{\partial y}}{g}_{\mathrm{002}}\left(x,y,z,\tilde{f}\right)

(36d) g_{\mathrm{112}}\left(x,y,z,\tilde{f}\right)=g_{\mathrm{121}}\left(y,x,z,\tilde{f}\right)={\displaystyle \frac{\partial }{\partial x}}{g}_{\mathrm{002}}\left(x,y,z,\tilde{f}\right)

(36e) g_{\mathrm{130}}\left(x,y,z,\tilde{f}\right)=g_{\mathrm{103}}\left(y,x,z,\tilde{f}\right)={\displaystyle \frac{\partial }{\partial x}}{g}_{\mathrm{020}}\left(x,y,z,\tilde{f}\right)

Accordingly, the Green’s function matrices {\mathbf{G}}_v[L,z,\mathbf{\phi }] and {\mathbf{G}}_p[L,z,\mathbf{\phi }] ( L=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}; z\ge \mathrm{0} ) are defined as follows.

(37a) \mathrm{4}\mathrm{\pi }{\mathbf{G}}_v(L,z,\mathbf{\phi })={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{L-\mathrm{1}}{e}^{-\rho z}\mathbf{\Pi }\mathbf{\phi }{\mathbf{\Pi }}^{*}K\mathrm{d}\xi \mathrm{d}\eta }={\phi }_{\mathrm{22}}\left(\begin{array}{ccc}{g}_{L\mathrm{02}}& -g_{L\mathrm{11}}& \mathrm{0}\\ -g_{L\mathrm{11}}& g_{L\mathrm{20}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)+\left(\begin{array}{ccc}{\phi }_{\mathrm{11}}{g}_{L\mathrm{20}}& {\phi }_{\mathrm{11}}{g}_{L\mathrm{11}}& {\phi }_{\mathrm{13}}{g}_{L\mathrm{10}}\\ {\phi }_{\mathrm{11}}{g}_{L\mathrm{11}}& {\phi }_{\mathrm{11}}{g}_{L\mathrm{02}}& {\phi }_{\mathrm{13}}{g}_{L\mathrm{01}}\\ -{\phi }_{\mathrm{31}}{g}_{L\mathrm{10}}& -{\phi }_{\mathrm{31}}{g}_{L\mathrm{01}}& {\phi }_{\mathrm{33}}{g}_{L\mathrm{00}}\end{array}\right)

(37b) \mathrm{4}\mathrm{\pi }{\mathbf{G}}_p(L,z,\mathbf{\phi })={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{L-\mathrm{1}}{e}^{-\rho z}{\mathbf{\Pi }}_p\mathbf{\phi }{\mathbf{\Pi }}^{*}K\mathrm{d}\xi \mathrm{d}\eta }={\phi }_{\mathrm{22}}\left(\begin{array}{ccc}{g}_{L\mathrm{12}}& -g_{L\mathrm{21}}& \mathrm{0}\\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(g_{L\mathrm{03}}-g_{L\mathrm{21}})& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(g_{L\mathrm{30}}-g_{L\mathrm{12}})& \mathrm{0}\\ -g_{L\mathrm{12}}& g_{L\mathrm{21}}& \mathrm{0}\end{array}\right)+\left(\begin{array}{ccc}{\phi }_{\mathrm{11}}{g}_{L\mathrm{30}}& {\phi }_{\mathrm{11}}{g}_{L\mathrm{21}}& -{\phi }_{\mathrm{13}}{g}_{L\mathrm{20}}\\ {\phi }_{\mathrm{11}}{g}_{L\mathrm{21}}& {\phi }_{\mathrm{11}}{g}_{L\mathrm{12}}& -{\phi }_{\mathrm{13}}{g}_{L\mathrm{11}}\\ {\phi }_{\mathrm{11}}{g}_{L\mathrm{12}}& {\phi }_{\mathrm{11}}{g}_{L\mathrm{03}}& -{\phi }_{\mathrm{13}}{g}_{L\mathrm{02}}\end{array}\right)

where {\phi }_{\mathrm{11}} , {\phi }_{\mathrm{22}} , {\phi }_{\mathrm{33}} , {\phi }_{\mathrm{13}} , and {\phi }_{\mathrm{31}} are the five elements of the material constant matrix \mathbf{\phi } defined as follows.

(37c) \mathbf{\phi }=\left(\begin{array}{ccc}{\phi }_{\mathrm{11}}& \mathrm{0}& {\phi }_{\mathrm{13}}\\ \mathrm{0}& {\phi }_{\mathrm{22}}& \mathrm{0}\\ {\phi }_{\mathrm{31}}& \mathrm{0}& {\phi }_{\mathrm{33}}\end{array}\right)

Closed-form results for point loading

The point body force f_{po\mathrm{int}}(x,y) and its Fourier integral transform {\tilde{f}}_{po\mathrm{int}} can be expressed as follows.

(38a) f_{po\mathrm{int}}(x,y)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}\delta (x)\delta (y)

(38b) {\tilde{f}}_{po\mathrm{int}}={\tilde{f}}_{po\mathrm{int}}(\xi ,\eta )=1

The six basic harmonic functions g_{\mathrm{0}IJ} for L=\mathrm{0} can be integrated in closed-forms in terms of elementary functions [ 2, 4, 5]. The other harmonic functions g_{LIJ} for L\ge \mathrm{1} and \mathrm{0}\le I+J\le \mathrm{3} can be found using the recursive relations Eq. (36). Expressions of those harmonic functions are specifically given in Appendix B.

Closed-form results for rectangular loading

The rectangular loading is the uniformly distributed for over a horizontal rectangular area. It and its Fourier integral transform can be expressed as follows.

(39a) f_{rect}(x,y)=\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathrm{4}ab}}& \mathrm{i}\mathrm{f} \left|x\right|\le a \mathrm{and} \left|y\right|\le b\\ \mathrm{0}& \mathrm{i}\mathrm{f} \left|x\right|>a \mathrm{or} \left|y\right|>b\end{array}.

(39b) {\tilde{f}}_{rect}={\tilde{f}}_{rect}(\xi ,\eta )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }} {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(x,y) K^{*} \mathrm{d}x \mathrm{d}y}}={\displaystyle \frac{\sin (a\xi )\sin (b\eta )}{\mathrm{\pi }ab\xi \eta }}

where a>0; b>0.

The six basic harmonic functions g_{\mathrm{0}IJ} for L=\mathrm{0} can be integrated in closed-forms in terms of elementary functions [ 3, 6]. The other harmonic functions g_{LIJ} for L\ge \mathrm{1} and \mathrm{0}\le I+J\le \mathrm{3} can be found using the recursive relations Eq. (36). Expressions of those harmonic functions are specifically given in Appendix C.

Closed-form result for circular ring loading

The circular ring loading is the body force vector uniformly concentrated on the circular ring and can be expressed as follows in cylindrical coordinate system.

(40) f_{ring}(r,\theta )={\displaystyle \frac{\delta (r-r_{\mathrm{0}})}{\mathrm{2}\mathrm{\pi }r}}

Four basic functions q_{km} are defined as follows in terms of products of the Bessel functions of orders of zero and one.

(41a) q_{\mathrm{00}}=J_{\mathrm{0}}(\rho r)J_{\mathrm{0}}(\rho r_{\mathrm{0}})

(41b) q_{\mathrm{01}}=J_{\mathrm{0}}(\rho r)J_{\mathrm{1}}(\rho r_{\mathrm{0}})

(41c) q_{\mathrm{10}}=J_{\mathrm{1}}(\rho r)J_{\mathrm{0}}(\rho r_{\mathrm{0}})

(41d) q_{\mathrm{11}}=J_{\mathrm{1}}(\rho r)J_{\mathrm{1}}(\rho r_{\mathrm{0}})

Secondly, the following harmonic functions g_{LIJ} are defined,

(42) q_{LIJ}=q_{LIJ}(r,r_{\mathrm{0}},z)={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\rho }^{L-\mathrm{1}}{e}^{-\rho z}{q}_{IJ}\mathrm{d}\rho }

where L, I and J are nonnegative integers, I=\mathrm{0} \mathrm{or} \mathrm{1},J=\mathrm{0} \mathrm{or} \mathrm{1} ; for L=\mathrm{0} \mathrm{or} \mathrm{1} , z\ge \mathrm{0} ; for L\ge \mathrm{2} , z>\mathrm{0} . It can be shown that the harmonic functions g_{LIJ}\left(z,\tilde{P}\right) satisfy the Laplace equations as follows:

(43) {\displaystyle \frac{\partial }{r\partial r}}\left(r{\displaystyle \frac{\partial }{\partial r}}{q}_{LIJ}(r,r_{\mathrm{0}},z)\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}{q}_{LIJ}(r,r_{\mathrm{0}},z)=\mathrm{0}

They also have the following recursive relations.

(44) g_{(L+\mathrm{1})IJ}\left(x,y,z,\tilde{f}\right)=-{\displaystyle \frac{\partial }{\partial z}}{g}_{LIJ}\left(x,y,z,\tilde{f}\right)

Accordingly, the Green’s function matrices {\mathbf{G}}_v[L,z,\mathbf{\phi }] and {\mathbf{G}}_p[L,z,\mathbf{\phi }] ( L=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3};z\ge \mathrm{0} ) are defined in the following.

(45a) 4\mathrm{\pi }{\mathbf{G}}_v\left(L,z,\mathbf{\phi }\right)={\displaystyle \underset{0}{\overset{+\infty }{\int }}{\rho }^L{e}^{-\rho z}{\mathbf{\Pi }}_{c0}\left(\rho r\right)\mathbf{\phi }{\mathbf{\Pi }}_{c0}\left(\rho r_0\right)d\rho }=\left(\begin{array}{ccc}{\phi }_{11}{q}_{L11}& 0& {\phi }_{13}{g}_{L10}\\ 0& {\phi }_2{g}_{L11}& 0\\ -{\phi }_{31}{q}_{L01}& 0& {\phi }_{33}{g}_{L00}\end{array}\right)

(45b) 4\mathrm{\pi }{\mathbf{G}}_p\left(L+1,z,\mathbf{\phi }\right)={\displaystyle \frac{1}{2\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{L+1}{e}^{-\rho z}{\mathbf{\Pi }}_{cp0}\left(\rho r\right)\mathbf{\phi }{\mathbf{\Pi }}_{c0}\left(\rho r_0\right)\mathrm{d}\rho }=\left(\begin{array}{ccc}{\phi }_{11}{q}_{\left(L+1\right)01}& 0& -{\phi }_{13}{q}_{\left(L+1\right)00}\\ 0& \frac{1}{2}{\phi }_{22}{q}_{\left(L+1\right)01}& 0\\ 0& 0& 0\end{array}\right)-{\displaystyle \frac{1}{r}}\left(\begin{array}{ccc}{\phi }_{11}{q}_{L11}& 0& -{\phi }_{13}{g}_{L10}\\ 0& {\phi }_{22}{g}_{L11}& 0\\ -{\phi }_{11}{g}_{L11}& 0& {\phi }_{13}{g}_{L10}\end{array}\right)

The four basic harmonic functions q_{0IJ} for L=0 can be integrated in closed-forms in terms of the complete elliptic integrals of the first, second and third kind [ 2, 4, 5]. The other harmonic functions q_{LIJ} for L\ge 1 can be found using the recursive relations Eq. (44). Expressions of those harmonic functions are specifically given in Appendix D.

The basic solutions of u , T_z , and {\mathrm{\Gamma }}_p

Accordingly, the solution of the displacement vector \mathbf{u} , vertical stresses {\mathbf{T}}_z and the plane strains {\mathbf{\Gamma }}_p induced by the concentrated loading f\left(x,y,z\right)=f\left(x,y\right)\delta \left(z-d\right)L or f\left(r,r_0,z\right)=\frac{\delta \left(r-r_0\right)}{2\pi r}\delta \left(z-d\right)L can be presented in the following unified matrix forms. For the three particular loading in Eqs. (38), (39) and (40), the solutions are presented in closed-form as follows.

The solutions for the three Kelvin’s cases

(46a) \mathbf{u}=\left\{{\mathbf{G}}_v\left(\mathrm{0},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{1}}\right)\right\}\mathbf{L}

(46b) {\mathbf{T}}_z=\left\{{\mathbf{G}}_v\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Psi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{2},\left|\overline{z}\right|,{\mathbf{\Psi }}_{k\mathrm{1}}\right)\right\}\mathbf{L}

(46c) {\mathbf{\Gamma }}_p=\left\{{\mathbf{G}}_p\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{1}}\right)\right\}\mathbf{L}

where and \overline{z}=z-d.

In particular, the Kelvin’s solution for the point loading can be specifically expressed as follows,

(47a) \mathbf{u}={\displaystyle \frac{1}{8\mathrm{\pi }\mu R}}\left\{\left(1+\alpha \right)\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)+{\displaystyle \frac{1-\alpha }{R^2}}\left(\begin{array}{ccc}{x}^2& xy& x\overline{z}\\ xy& y^2& y\overline{z}\\ x\overline{z}& y\overline{z}& {\overline{z}}^2\end{array}\right)\right\}\mathbf{L}

(47b) {\mathbf{T}}_z= {\displaystyle \frac{-1}{4\mathrm{\pi }{R}^3}}\left\{\alpha \left(\begin{array}{ccc}\overline{z}& 0& x\\ 0& \overline{z}& y\\ -x& -y& \overline{z}\end{array}\right)+{\displaystyle \frac{3\left(1-\alpha \right)\overline{z}}{R^2}}\left(\begin{array}{ccc}{x}^2& xy& x\overline{z}\\ xy& y^2& y\overline{z}\\ x\overline{z}& y\overline{z}& {\overline{z}}^2\end{array}\right)\right\}\mathbf{L}

(47c) T_p={\displaystyle \frac{-1}{4\mathrm{\pi }{R}^3}}\left\{\alpha \left(\begin{array}{ccc}x& -y& -\overline{z}\\ y& x& 0\\ -x& y& -\overline{z}\end{array}\right)+{\displaystyle \frac{3\left(1-\alpha \right)}{R^2}}\left(\begin{array}{ccc}{x}^3& x^{2}y& x^2\overline{z}\\ x^{2}y& x{y}^2& xy\overline{z}\\ x{y}^2& y^3& y^2\overline{z}\end{array}\right)\right\}L

where R=\sqrt{x^2+y^2+{\overline{z}}^2} .

The solutions for the three Boussinesq’s cases

(48a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z,{\mathbf{\Phi }}_{b\mathrm{0}})+z{\mathbf{G}}_v(\mathrm{1},z,{\mathbf{\Phi }}_{b\mathrm{1}})\}\mathbf{L}

(48b) {\mathbf{T}}_z=\{G_v(\mathrm{1},z,{\mathbf{\Psi }}_{b\mathrm{0}})+z{\mathbf{G}}_v(\mathrm{2},z,{\mathbf{\Psi }}_{b\mathrm{1}})\}\mathbf{L}

(48c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z,{\mathbf{\Phi }}_{b\mathrm{0}})+z{\mathbf{G}}_p(\mathrm{2},z,{\mathbf{\Phi }}_{b\mathrm{1}})\}\mathbf{L}

where and . In particular, the solution for the point loading is the Boussinesq’s solution.

The solutions for the three Mindlin’s cases

(49a) \mathbf{u}=\left\{{\mathbf{G}}_v\left(0,z_1,{\mathbf{\Phi }}_{m0}\right)+d{\mathbf{G}}_v\left(1,z_1,{\mathbf{\Phi }}_{m1d}\right)+z{\mathbf{G}}_v\left(1,z_1,{\mathbf{\Phi }}_{m1z}\right)+dz{\mathbf{G}}_v\left(2,z_1,{\mathbf{\Phi }}_{m2}\right)+{\mathbf{G}}_v\left(0,\left|\overline{z}\right|,{\mathbf{\Phi }}_{k0}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(1,\left|\overline{z}\right|,{\mathbf{\Phi }}_{k1}\right)\right\}\mathbf{L}

(49b) {\mathbf{T}}_z=\left\{{\mathbf{G}}_v\left(1,z_1,{\mathbf{\Psi }}_{m0}\right)+d{\mathbf{G}}_v\left(2,z_1,{\mathbf{\Psi }}_{m1d}\right)+z{\mathbf{G}}_v\left(2,z_1,{\mathbf{\Psi }}_{m1z}\right)+dz{\mathbf{G}}_v\left(3,z_1,{\mathbf{\Psi }}_{m2}\right)+{\mathbf{G}}_v\left(0,\left|\overline{z}\right|,{\mathbf{\Psi }}_{k0}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(1,\left|\overline{z}\right|,{\mathbf{\Psi }}_{k1}\right)\right\}\mathbf{L}

(49c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{1}},{\mathbf{\Phi }}_{m\mathrm{0}})+d{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Phi }}_{m\mathrm{1}d})+z{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Phi }}_{m\mathrm{1}z})+dz{\mathbf{G}}_p(\mathrm{3},z_{\mathrm{1}},{\mathbf{\Phi }}_{m\mathrm{2}})+{\mathbf{G}}_p\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{1}}\right)\}\mathbf{L}

where , , z_{\mathrm{1}}=d+z , and \bar{z}=z-d . The solution has a part of the boundary effect and a part of Kelvin solution. In particular, the point loading case is the Mindlin’s solution.

The solutions for the three bi-material cases

In the rock solid of the semi-infinite region k=\mathrm{1} ( ), we have

(50a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{01}})+d{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{1}})+z{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{1}})+dz{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{21}})+{\mathbf{G}}_v\left(\mathrm{0},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{1}}\right)\}\mathbf{L}

(50b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{1}},{\mathbf{\Psi }}_{t\mathrm{01}})+d{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{1}})+z{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Psi }}_{t\mathrm{1}z\mathrm{1}})+dz{\mathbf{G}}_v(\mathrm{3},z_{\mathrm{1}},{\mathbf{\Psi }}_{t\mathrm{21}})+{\mathbf{G}}_v\left(\mathrm{0},\left|\overline{z}\right|,{\mathbf{\Psi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Psi }}_{k\mathrm{1}}\right)\}\mathbf{L}

(50c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{01}})+d{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{1}})+z{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{1}})+dz{\mathbf{G}}_p(\mathrm{3},z_{\mathrm{1}},{\mathbf{\Phi }}_{t\mathrm{21}})+{\mathbf{G}}_p\left(\mathrm{1},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},\left|\overline{z}\right|,{\mathbf{\Phi }}_{k\mathrm{1}}\right)\}\mathbf{L}

where and z_{\mathrm{1}}=d+z . The solution has a part of the interface effect and the part of Kelvin solution.

In the rock solid of the semi-infinite region , we have

(51a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{02}})+d{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{2}})+z{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{2}})+dz{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{22}})\}\mathbf{L}

(51b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{2}},{\mathbf{\Psi }}_{t\mathrm{012}})+d{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{2}},{\mathbf{\Psi }}_{t\mathrm{1}d\mathrm{2}})+z{\mathbf{G}}_v(\mathrm{2},z_{\mathrm{2}},{\mathbf{\Psi }}_{t\mathrm{1}z\mathrm{2}})+dz{\mathbf{G}}_v(\mathrm{3},z_{\mathrm{2}},{\mathbf{\Psi }}_{t\mathrm{22}})\}\mathbf{L}

(51c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{02}})+d{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{1}d\mathrm{2}})+z{\mathbf{G}}_p(\mathrm{2},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{1}z\mathrm{2}})+dz{\mathbf{G}}_p(\mathrm{3},z_{\mathrm{2}},{\mathbf{\Phi }}_{t\mathrm{22}})\}\mathbf{L}

where and z_{\mathrm{2}}=d-z .

The solutions for the three extended Kelvin’s cases

Case 1: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0}

(52a) \mathbf{u}=\{{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

(52b) {\mathbf{T}}_Z=\{{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

(52c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_{{\mathrm{\Gamma }}_p}\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

Case 2: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0}

(53a) \mathbf{u}=\{{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(53b) {\mathbf{T}}_Z=\{{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{1}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(53c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

The solutions for the three extended bi-material cases

Case a: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

In the rock solid of the semi-infinite region k=\mathrm{1} ( ), we have

(54a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_v(\mathrm{0},z_{aJ\mathrm{1}},{\mathbf{\Phi }}_{aJ\mathrm{1}})}+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

(54b) {\mathbf{T}}_Z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Psi }}_{\mathrm{01}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_v(\mathrm{1},z_{aJ\mathrm{1}},{\mathbf{\Psi }}_{aJ\mathrm{1}})}+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

(54c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_p(\mathrm{1},z_{aJ\mathrm{1}},{\mathbf{\Phi }}_{aJ\mathrm{1}})}+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

In the rock solid of the semi-infinite region , we have

(55a) \mathbf{u}=\left\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_v(\mathrm{0},z_{nJ\mathrm{2}},{\mathbf{\Phi }}_{aJ\mathrm{2}})}\right\}\mathbf{L}

(55b) {\mathbf{T}}_Z=\left\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{02}},{\mathbf{\Psi }}_{\mathrm{02}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_v(\mathrm{1},z_{nJ\mathrm{2}},{\mathbf{\Psi }}_{aJ\mathrm{2}})}\right\}\mathbf{L}

(55c) {\mathbf{\Gamma }}_p=\left\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_p(\mathrm{1},z_{nJ\mathrm{2}},{\mathbf{\Phi }}_{aJ\mathrm{2}})}\right\}\mathbf{L}

Case b: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}\ne \mathrm{0}

In the rock solid of the semi-infinite region k=\mathrm{1} ( ), we have

(56a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\mathbf{G}}_v(\mathrm{0},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{11}})+z{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{21}})+d{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{31}})+zd{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{41}})+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(56b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Psi }}_{\mathrm{01}})+{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Psi }}_{b\mathrm{11}})+z{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Psi }}_{b\mathrm{21}})+d{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Psi }}_{b\mathrm{31}})+zd{\mathbf{G}}_v(\mathrm{3},z_{b\mathrm{1}},{\mathbf{\Psi }}_{b\mathrm{41}})+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(56c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\mathbf{G}}_p(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{11}})+z{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{21}})+d{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{31}})+zd{\mathbf{G}}_p(\mathrm{3},z_{b\mathrm{1}},{\mathbf{\Phi }}_{b\mathrm{41}})+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

In the rock solid of the semi-infinite region , we have

(57a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\mathbf{G}}_v(\mathrm{0},z_{b\mathrm{12}},{\mathbf{\Phi }}_{b\mathrm{12}})+d{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{12}},{\mathbf{\Phi }}_{b\mathrm{22}})+{\mathbf{G}}_v(\mathrm{0},z_{b\mathrm{22}},{\mathbf{\Phi }}_{b\mathrm{32}})+d{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{22}},{\mathbf{\Phi }}_{b\mathrm{42}})\}\mathbf{L}

(57b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{02}},{\mathbf{\Psi }}_{\mathrm{02}})+{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{12}},{\mathbf{\Psi }}_{b\mathrm{12}})+d{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{12}},{\mathbf{\Psi }}_{b\mathrm{22}})+{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{22}},{\mathbf{\Psi }}_{b\mathrm{32}})+d{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{22}},{\mathbf{\Psi }}_{b\mathrm{42}})\}\mathbf{L}

(57c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\mathbf{G}}_p(\mathrm{1},z_{b\mathrm{12}},{\mathbf{\Phi }}_{b\mathrm{12}})+d{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{12}},{\mathbf{\Phi }}_{b\mathrm{22}})+{\mathbf{G}}_p(\mathrm{1},z_{b\mathrm{22}},{\mathbf{\Phi }}_{b\mathrm{32}})+d{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{22}},{\mathbf{\Phi }}_{b\mathrm{42}})\}\mathbf{L}

Case c: {\mathrm{\Delta }}_{\mathrm{1}}\ne \mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

In the rock solid of the semi-infinite region k=\mathrm{1} ( ),

(58a) \mathbf{u}=\left\{{\mathbf{G}}_v\left(0,z_{01},{\mathbf{\Phi }}_{01}\right)+{\displaystyle \sum _{J=1}^4{\mathbf{G}}_v\left(0,z_{aJ1},{\mathbf{\Phi }}_{cJ1}\right)}+{\mathbf{G}}_v\left(0,{\gamma }_{01}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(0,{\gamma }_{11}\left|\overline{z}\right|,{\mathbf{\Phi }}_u\left({\gamma }_{11}\right)\right)-{\mathbf{G}}_v\left(0,{\gamma }_{21}\left|\overline{z}\right|,{\mathbf{\Phi }}_u\left({\gamma }_{21}\right)\right)\right\}\mathbf{L}

(58b) {\mathbf{T}}_Z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Psi }}_{\mathrm{01}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_v(\mathrm{1},z_{aJ\mathrm{1}},{\mathbf{\Psi }}_{cJ\mathrm{1}})}+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Psi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

(58c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\displaystyle \sum _{J=\mathrm{1}}^{\mathrm{4}}{\mathbf{G}}_p(\mathrm{1},z_{aJ\mathrm{1}},{\mathbf{\Phi }}_{cJ\mathrm{1}})}+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{11}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{11}})\right)-{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{21}}\left|\overline{z}\right|,{\mathbf{\Phi }}_u({\gamma }_{\mathrm{21}})\right)\}\mathbf{L}

In the rock solid of the semi-infinite region , we have

(59a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\mathbf{G}}_v(\mathrm{0},z_{c\mathrm{12}},{\mathbf{\Phi }}_{c\mathrm{12}})+z{\mathbf{G}}_v(\mathrm{1},z_{c\mathrm{12}},{\mathbf{\Phi }}_{c\mathrm{22}})+{\mathbf{G}}_v(\mathrm{0},z_{c\mathrm{22}},{\mathbf{\Phi }}_{c\mathrm{32}})+z{\mathbf{G}}_v(\mathrm{1},z_{c\mathrm{22}},{\mathbf{\Phi }}_{c\mathrm{42}})\}\mathbf{L}

(59b) T_z=\{G_v(\mathrm{1},z_{\mathrm{02}},{\mathrm{\Psi }}_{\mathrm{02}})+G_v(\mathrm{1},z_{c\mathrm{12}},{\mathrm{\Psi }}_{c\mathrm{12}})+z{G}_v(\mathrm{2},z_{c\mathrm{12}},{\mathrm{\Psi }}_{c\mathrm{22}})+G_v(\mathrm{1},z_{c\mathrm{22}},{\mathrm{\Psi }}_{c\mathrm{32}})+z{G}_v(\mathrm{2},z_{c\mathrm{22}},{\mathrm{\Psi }}_{c\mathrm{42}})\}L

(59c) {\mathrm{\Gamma }}_p=\{G_p(\mathrm{1},z_{\mathrm{02}},{\mathrm{\Phi }}_{\mathrm{02}})+G_p(\mathrm{1},z_{c\mathrm{12}},{\mathrm{\Phi }}_{c\mathrm{12}})+z{G}_p(\mathrm{2},z_{c\mathrm{12}},{\mathrm{\Phi }}_{c\mathrm{22}})+G_p(\mathrm{1},z_{c\mathrm{22}},{\mathrm{\Phi }}_{c\mathrm{32}})+z{G}_p(\mathrm{2},z_{c\mathrm{22}},{\mathrm{\Phi }}_{c\mathrm{42}})\}L

Case d: {\mathrm{\Delta }}_{\mathrm{1}}=\mathrm{0} and {\mathrm{\Delta }}_{\mathrm{2}}=\mathrm{0}

In the rock solid of the semi-infinite region k=\mathrm{1} ( )

(60a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\mathbf{G}}_v(\mathrm{0},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{11}})+z{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{21}})+d{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{31}})+zd{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{41}})+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{0},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(60b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Psi }}_{\mathrm{01}})+{\mathbf{G}}_v(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Psi }}_{d\mathrm{11}})+z{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Psi }}_{d\mathrm{21}})+d{\mathbf{G}}_v(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Psi }}_{d\mathrm{31}})+zd{\mathbf{G}}_v(\mathrm{3},z_{b\mathrm{1}},{\mathbf{\Psi }}_{d\mathrm{41}})+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Psi }}_v\right)+{\mathbf{G}}_v\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_v\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Psi }}_{u\mathrm{1}}\right)\}\mathbf{L}

(60c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\mathbf{G}}_p(\mathrm{1},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{11}})+z{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{21}})+d{\mathbf{G}}_p(\mathrm{2},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{31}})+zd{\mathbf{G}}_p(\mathrm{3},z_{b\mathrm{1}},{\mathbf{\Phi }}_{d\mathrm{41}})+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{01}}\left|\overline{z}\right|,{\mathbf{\Phi }}_v\right)+{\mathbf{G}}_p\left(\mathrm{1},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{0}}\right)+\left|\overline{z}\right|{\mathbf{G}}_p\left(\mathrm{2},{\gamma }_{\mathrm{31}}\left|\overline{z}\right|,{\mathbf{\Phi }}_{u\mathrm{1}}\right)\}\mathbf{L}

In the rock solid of the semi-infinite region , we have

(61a) \mathbf{u}=\{{\mathbf{G}}_v(\mathrm{0},z_{\mathrm{02}},{\mathbf{\Phi }}_{\mathrm{02}})+{\mathbf{G}}_v(\mathrm{0},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{12}})+z{\mathbf{G}}_v(\mathrm{1},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{22}})+d{\mathbf{G}}_v(\mathrm{1},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{32}})+zd{\mathbf{G}}_v(\mathrm{2},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{21}})\}\mathbf{L}

(61b) {\mathbf{T}}_z=\{{\mathbf{G}}_v(\mathrm{1},z_{\mathrm{02}},{\mathbf{\Psi }}_{\mathrm{02}})+{\mathbf{G}}_v(\mathrm{1},z_{d\mathrm{2}},{\mathbf{\Psi }}_{d\mathrm{12}})+z{\mathbf{G}}_v(\mathrm{2},z_{d\mathrm{2}},{\mathbf{\Psi }}_{d\mathrm{22}})+d{\mathbf{G}}_v(\mathrm{2},z_{d\mathrm{2}},{\mathbf{\Psi }}_{d\mathrm{32}})+zd{\mathbf{G}}_v(\mathrm{3},z_{d\mathrm{2}},{\mathbf{\Psi }}_{d\mathrm{42}})\}\mathbf{L}

(61c) {\mathbf{\Gamma }}_p=\{{\mathbf{G}}_p(\mathrm{1},z_{\mathrm{01}},{\mathbf{\Phi }}_{\mathrm{01}})+{\mathbf{G}}_p(\mathrm{1},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{12}})+z{\mathbf{G}}_p(\mathrm{2},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{22}})+d{\mathbf{G}}_p(\mathrm{2},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{32}})+zd{\mathbf{G}}_p(\mathrm{3},z_{d\mathrm{2}},{\mathbf{\Phi }}_{d\mathrm{42}})\}\mathbf{L}

Singularities of the basic solutions

The harmonic functions g_{LIJ} in terms of the improper integrations in Eqs. (34) and (42) are absolutely and uniformly convergent if \left|z-d\right|>\epsilon >\mathrm{0} , where \epsilon is a positive small value. However, at the loading plane z\to d^{\pm } ( d\ge \mathrm{0} ), \mathrm{\Phi }(\rho ,z) and \mathrm{\Psi }(\rho ,z) have the constant matrix terms for . The g_{LIJ} for z\to d^{\pm } at the loading plane are expressed in terms of the following singular integrals.

(62a) \underset{z\to d^{\pm }}{\lim }{g}_{\mathrm{0}IJ}=g_{\mathrm{0}IJ}\left(x,y,{\mathrm{0}}^{\pm },\tilde{f}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{-\mathrm{1}}{q}_{IJ}\tilde{f}K\mathrm{d}\xi \mathrm{d}\eta }

(62b) \underset{z\to d^{\pm }}{\lim }{g}_{\mathrm{1}IJ}=g_{\mathrm{1}IJ}\left(x,y,{\mathrm{0}}^{\pm },\tilde{f}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{\mathrm{0}}{q}_{IJ}\tilde{f}K\mathrm{d}\xi \mathrm{d}\eta }

They are convergent in the sense of Cauchy principal values. They cannot be directly integrated accurately and rapidly with normal numerical integration methods. They have to be isolated and integrated analytically in closed-forms.

For the point loading, using the closed-form results, we have the following limiting results of the singular integrations g_{\mathrm{0}IJ} and g_{\mathrm{1}IJ} as z\to {\mathrm{0}}^{+} .

(63a) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{000}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{\mathrm{1}}{r}}

(63b) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{010}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i\xi }{{\rho }^{\mathrm{2}}}}K\mathrm{d}\xi \mathrm{d}\eta }=-{\displaystyle \frac{x}{r^{\mathrm{2}}}}

(63c) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{020}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{2}}}{{\rho }^{\mathrm{3}}}}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{\mathrm{1}}{r}}\left(\mathrm{1}-{\displaystyle \frac{x^{\mathrm{2}}}{r^{\mathrm{2}}}}\right)

(63d) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{011}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\xi \eta }{{\rho }^{\mathrm{3}}}}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{xy}{r^{\mathrm{3}}}}

(63e) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{030}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{3}}}{{\rho }^{\mathrm{4}}}}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{x}{\mathrm{2}{r}^{\mathrm{2}}}}\left(\mathrm{3}-{\displaystyle \frac{\mathrm{2}{x}^{\mathrm{2}}}{r^{\mathrm{2}}}}\right)

(63f) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{021}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{2}}\eta }{{\rho }^{\mathrm{4}}}}K\mathrm{d}\xi \mathrm{d}\eta }=-{\displaystyle \frac{y}{\mathrm{2}{r}^{\mathrm{2}}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}{x}^{\mathrm{2}}}{r^{\mathrm{2}}}}\right)

and

(64a) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{100}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}K\mathrm{d}\xi \mathrm{d}\eta }=\mathrm{2}\mathrm{\pi }\delta (x)\delta (y)

(64b) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{110}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i\xi }{\rho }}K\mathrm{d}\xi \mathrm{d}\eta }=-{\displaystyle \frac{x}{r^{\mathrm{3}}}}

(64c) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{120}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{y^{\mathrm{2}}-x^{\mathrm{2}}}{r^{\mathrm{4}}}}+\mathrm{\pi }\delta (x)\delta (y)

(64d) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{111}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\xi \eta }{{\rho }^{\mathrm{2}}}}K\mathrm{d}\xi \mathrm{d}\eta }={\displaystyle \frac{\mathrm{2}xy}{r^{\mathrm{4}}}}

(64e) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{130}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{3}}}{{\rho }^{\mathrm{3}}}}K\mathrm{d}\xi \mathrm{d}\eta }=-{\displaystyle \frac{\mathrm{3}x}{\mathrm{2}{r}^{\mathrm{3}}}}\left(\mathrm{1}-{\displaystyle \frac{x^{\mathrm{2}}}{r^{\mathrm{2}}}}\right)

(64f) \underset{z\to {\mathrm{0}}^{+}}{\lim }{g}_{\mathrm{121}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{i{\xi }^{\mathrm{2}}\eta }{{\rho }^{\mathrm{3}}}}K\mathrm{d}\xi \mathrm{d}\eta }=-{\displaystyle \frac{y}{r^{\mathrm{3}}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}{x}^{\mathrm{2}}}{r^{\mathrm{2}}}}\right)

where r=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}>\mathrm{0.}

Similarly, the singular integrals for q_{LIJ} at the horizontal plane of the concentrated circular ring loading have the following results.

(65a) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{000}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\displaystyle \frac{J_{\mathrm{0}}(\rho r)J_{\mathrm{0}}(\rho r_{\mathrm{0}})}{\rho }}\mathrm{d}\rho }={\displaystyle \frac{\mathrm{2}}{\mathrm{\pi }(r+r_{\mathrm{0}})}}K({\kappa }_{\mathrm{1}})

(65b) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{001}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\displaystyle \frac{J_{\mathrm{0}}(\rho r)J_{\mathrm{1}}(\rho r_{\mathrm{0}})}{\rho }}\mathrm{d}\rho }={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{0}}}}H(r_{\mathrm{0}}-r)

(65c) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{011}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\displaystyle \frac{J_{\mathrm{1}}(\rho r)J_{\mathrm{1}}(\rho r_{\mathrm{0}})}{\rho }}\mathrm{d}\rho }={\displaystyle \frac{\mathrm{1}}{\mathrm{\pi }r{r}_{\mathrm{0}}}}\left({\displaystyle \frac{r^{\mathrm{2}}+r_{\mathrm{0}}^{\mathrm{2}}}{r+r_{\mathrm{0}}}}K({\kappa }_{\mathrm{1}})-(r+r_{\mathrm{0}})E({\kappa }_{\mathrm{1}})\right)

and

(66a) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{100}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{J}_{\mathrm{0}}(\rho r)J_{\mathrm{0}}(\rho r_{\mathrm{0}})\mathrm{d}\rho }={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{0}}}}\delta (r-r_{\mathrm{0}})

(66b) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{101}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{J}_{\mathrm{0}}(\rho r)J_{\mathrm{1}}(\rho r_{\mathrm{0}})\mathrm{d}\rho }={\displaystyle \frac{\mathrm{1}}{\mathrm{\pi }{r}_{\mathrm{0}}}}\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{0}}-r}}E({\kappa }_{\mathrm{1}})+{\displaystyle \frac{\mathrm{1}}{r+r_{\mathrm{0}}}}K({\kappa }_{\mathrm{1}})\right)

(66c) \underset{z\to \mathrm{0}}{\lim }{q}_{\mathrm{111}}={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{J}_{\mathrm{0}}(\rho r)J_{\mathrm{0}}(\rho r_{\mathrm{0}})\mathrm{d}\rho }={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{0}}}}\delta (r-r_{\mathrm{0}})

For general body force \tilde{f} , in particular, the rectangular body force, the above singular integrals can be integrated in their closed-forms using the expressions in Appendix C.

Summary notes

In the above, exact solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) have been given for specific boundary-value problems in either a homogenous solid of infinite or semi-infinite extent or a bi-material solid of infinite extent. The materials are either isotropic or transversely isotropic. The solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) are independent to the applied loading \mathbf{f}=f(x,y)\delta (z-d) on a horizontal plane and can be directly used to systematically formulate the solutions of \mathbf{u} , {\mathbf{T}}_z and {\mathbf{\Gamma }}_p in matrix form for any specific load distribution f(x,y) .

The three specific forces f(x,y) (i.e., point, rectangular and ring force distributions) are considered in details. Their solutions in physical domain are given exactly and systematically in closed-forms in terms of some harmonic functions. The closed-form results include those given by Kelvin in 1848, Boussinesq in 1885, Mindlin in 1936 and some others with various methods.

Furthermore, it can be shown that the closed-form solutions can be systematically and easily obtained for any other specific force f(x,y) if the following basic harmonic integral can be integrated in closed-form.

(67a) Q\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{e^{-\rho z}}{{\rho }^{\mathrm{3}}}}\tilde{f}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

or

(67b) Q(x,y,z,f)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}[z\ln (z+R)-R]f(s,t)\mathrm{d}s\mathrm{d}t}

where z>\mathrm{0} and R=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+z^{\mathrm{2}}} .

It is because of the following recursive relations of differentiations for the four basic harmonic functions and other harmonic functions.

(68a) g_{\mathrm{000}}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}Q\left(x,y,z,\tilde{f}\right).

(68b) g_{\mathrm{010}}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial z}}Q\left(x,y,z,\tilde{f}\right)

(68c) g_{\mathrm{011}}\left(x,y,z,,\tilde{f}\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial y}}Q\left(x,y,z,\tilde{f}\right)

(68d) g_{\mathrm{020}}\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}Q\left(x,y,z,\tilde{f}\right)

The mathematical properties of the basic solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) have been examined exactly and clearly. For all the cases, the improper integrals of \mathbf{u} , {\mathbf{T}}_z and {\mathbf{\Gamma }}_p are uniformly and absolutely convergent if \left|z-d\right|>\epsilon >\mathrm{0} . At the loading plane, the improper integrals become the singular integrals and have been analytically integrated for the point, ring and rectangular forces.

The solutions for the bi-solid cases are given under the perfectly bonded interface condition. Other interface bonding conditions for the bi-solids and other boundary conditions for a homogeneous solid of semi-infinite can also be examined similarly and systematically. In particular, the smooth bonded interface condition and the horizontally inextensible bonded interface condition have been examined in Refs. [ 4, 5]. Their solutions in closed-forms have been derived. The solutions in a homogeneous transversely isotropic solid of semi-infinite extent with traction free boundary condition have been examined in Ref. [ 6]. Their closed-formed solutions of the complete elastic field variables for all the boundary-value problems in the homogeneous solids of infinite or semi-infinite extent due to the three forces in Eqs. (38a), (39a) and (40) have been given exactly and similarly with the above method. The three concentrated loadings can be re-expressed as follows.

(69a) \mathbf{f}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}\delta (x)\delta (y)\delta (z-d)\mathbf{L}

(69b) \mathbf{f}(r,\theta ,z)={\displaystyle \frac{\delta (r-r_{\mathrm{0}})}{\mathrm{2}\mathrm{\pi }r}}\delta (z-d)\mathbf{L}

(69c) \mathbf{f}(x,y,z)=\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathrm{4}ab}}\delta (z-d)\mathbf{L},& if \left|x\right|\le a \mathrm{and} \left|y\right|\le b;\\ \mathrm{0},& if \left|x\right|>a \mathrm{o}\mathrm{r} \left|y\right|>b,\end{array}

where \mathbf{L} is a constant force vector.

Mathematical properties of the solution

General

The basic solutions given in the above section show the basic properties of the basic solutions and their common singularities. Similarly, they can be used to show the solutions given in the Sections 4 and 5 for the layered solids in Ref. [ 1] because they are also systematically expressed in matrix forms in terms of either inverse 2-D Fourier integral transforms or inverse Hankel integral transforms. The integrals are improper integrals with depending parameters and 2D-infinite or semi-infinite integration intervals. Therefore, the following issues have to be examined and answered in this section:

1) the mathematical properties of the solution matrices \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) (or similarly {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) ) in the transform domain;

2) the convergence of the improper integrals of the solution;

3) the singularities of the solution;

4) the interchangeability of the integration limits and the integrations and the partial differentiations for the solution of displacements, stresses and strains expressed in the forms of the improper integrals; and

5) Satisfaction of the solution to the governing equations, interfacial conditions and the imposed boundary and internal loading conditions.

The properties of {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z)

Using the principal of mathematical induction and/or numerical techniques and under the positive strain energy constraints of the five elastic parameters of each layer (7) in Ref. [ 1], it can be shown that the determinants of the four coefficient matrices {\mathbf{M}}_{Ap} , {\mathbf{M}}_{Aq} , {\mathbf{M}}_{Qp} and {\mathbf{M}}_{Qq} respectively in Eqs. (55c), (60c), (73c) and (78c) as well as the Appendices A to D in Ref. [ 1] do not have any zero value for any \rho in [ 7, 8]. i.e.,

(70a) \det {\mathbf{M}}_{Ap}(\rho )\ne \mathrm{0}

(70b) \det {\mathbf{M}}_{Aq}(\rho )\ne \mathrm{0}

(70c) \det {\mathbf{M}}_{Qp}(\rho )\ne \mathrm{0}

(70d) \det {\mathbf{M}}_{Qq}(\rho )\ne \mathrm{0}

where .

Accordingly, it is evidently that {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) have the following properties, which are similar to those of the basic solutions of \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) in the above section.

(a) {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) are continuous functions of \rho , where .

(b) {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) are continuous functions of z , where or .

(c) Within every single layer, {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) are partial differentiable of any orders with respect to z except at the loading plane z=d .

(d) At z=d , the leading diagonal elements of the six matrix solutions have the discontinuity of the first kind while the non-diagonal elements are continuous. The difference between the diagonal elements as z=d^{+} and z=d^{-} is unity, i.e.,

(71a) {\mathbf{\Psi }}_V(\rho ,d^{+})-{\mathbf{\Psi }}_V(\rho ,d^{-})=-{\mathbf{I}}_{\mathrm{2}}

(71b) {\mathbf{\Psi }}_U(\rho ,d^{+})-{\mathbf{\Psi }}_U(\rho ,d^{-})=-{\mathbf{I}}_{\mathrm{4}}

where I_{\mathrm{2}}=\mathrm{unit}\mathrm{matrix}\mathrm{of}\mathrm{2}\times \mathrm{2} and I_{\mathrm{4}}=\mathrm{unit}\mathrm{matrix}\mathrm{of}\mathrm{4}\times \mathrm{4} .

(e) In particular, as \rho \to M (M is a very large value approaching +\infty ), {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) have asymptotic representations and can be expressed as follows.

(i) For and \mathrm{0}\le j\le k ,

(72a) \underset{\rho \to M}{\lim }{\mathbf{\Psi }}_V(\rho ,z)\approx e^{-{\gamma }_{\mathrm{0}j}\rho (H_j-z)-{\gamma }_{\mathrm{0}(j+\mathrm{1})}\rho h_{j+\mathrm{1}}-\cdots -{\gamma }_{\mathrm{0}(k-\mathrm{1})}\rho h_{k-\mathrm{1}}-{\gamma }_{\mathrm{0}k}\rho (d-H_{k-\mathrm{1}})}{\mathbf{A}}_j^p(z-H_{j-\mathrm{1}}){\mathbf{\Gamma }}_{Ap}

(72b) \underset{\rho \to M}{\lim }{\mathbf{\Psi }}_U(\rho ,z)\approx e^{-{\gamma }_{aj}\rho (H_j-z)-{\gamma }_{a(j+\mathrm{1})}\rho h_{j+\mathrm{1}}-\cdots -{\gamma }_{a(k-\mathrm{1})}\rho h_{k-\mathrm{1}}-{\gamma }_{ak}\rho (d-H_{k-\mathrm{1}})}{\mathbf{Q}}_j^p(z-H_{j-\mathrm{1}}){\mathbf{\Gamma }}_{Qp}

(ii) For and k\le j\le n+\mathrm{1},

(72c) \underset{\rho \to M}{\lim }{\mathbf{\Psi }}_V(\rho ,z)\approx e^{-{\gamma }_{\mathrm{0}j}\rho (z-H_{j-\mathrm{1}})-{\gamma }_{\mathrm{0}(j-\mathrm{1})}\rho h_{(j-\mathrm{1})}-\cdots -{\gamma }_{\mathrm{0}(k+\mathrm{1})}\rho h_{(k+\mathrm{1})}-{\gamma }_{\mathrm{0}k}\rho (H_k-d)}{\mathbf{A}}_j^q(z-H_j){\mathbf{\Gamma }}_{Aq}

(72d) \underset{\rho \to M}{\lim }{\mathbf{\Psi }}_U(\rho ,z)\approx e^{-{\gamma }_{aj}\rho (z-H_{j-\mathrm{1}})-{\gamma }_{a(j-\mathrm{1})}\rho h_{(j-\mathrm{1})}-\cdots -{\gamma }_{a(k+\mathrm{1})}\rho h_{(k+\mathrm{1})}-{\gamma }_{ak}\rho (H_k-d)}{\mathbf{Q}}_j^q(z-H_j){\mathbf{\Gamma }}_{Qq}

where {\mathbf{\Gamma }}_{Ap} , {\mathbf{\Gamma }}_{Qp} , {\mathbf{\Gamma }}_{Aq} and {\mathbf{\Gamma }}_{Qq} are dependent on the material parameters.

(f) Furthermore, all the elements of {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) and their mth differential derivatives with respect to z have the following bounds.

(i) For and \mathrm{0}\le j\le k ,

(73a) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}{\mathbf{\Psi }}_V(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{\mathrm{0}j}\rho (H_j-z)-{\gamma }_{\mathrm{0}(j+\mathrm{1})}\rho h_{j+\mathrm{1}}-\cdots -{\gamma }_{\mathrm{0}(k-\mathrm{1})}\rho h_{k-\mathrm{1}}-{\gamma }_{\mathrm{0}k}\rho (d-H_{k-\mathrm{1}})}{\mathbf{\Psi }}_{V\max }^m\le {\rho }^m{e}^{-{\gamma }_{\mathrm{0}\min }\rho (d-z)}{\mathbf{\Psi }}_{V\max }^m

(73b) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}{\mathbf{\Psi }}_U(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{aj}\rho (H_j-z)-{\gamma }_{a(j+\mathrm{1})}\rho h_{j+\mathrm{1}}-\cdots -{\gamma }_{a(k-\mathrm{1})}\rho h_{k-\mathrm{1}}-{\gamma }_{ak}\rho (d-H_{k-\mathrm{1}})}{\mathbf{\Psi }}_{U\max }^m\le {\rho }^m{e}^{-{\gamma }_{a\min }\rho (d-z)}{\mathbf{\Psi }}_{U\max }^m

(ii) For and k\le j\le n+\mathrm{1} ,

(73c) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}{\mathbf{\Psi }}_V(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{\mathrm{0}j}\rho (z-H_{j-\mathrm{1}})-{\gamma }_{\mathrm{0}(j-\mathrm{1})}\rho h_{(j-\mathrm{1})}-\cdots -{\gamma }_{\mathrm{0}(k+\mathrm{1})}\rho h_{(k+\mathrm{1})}-{\gamma }_{\mathrm{0}k}\rho (H_k-d)}{\mathbf{\Psi }}_{V\max }^m\le {\rho }^m{e}^{-{\gamma }_{\mathrm{0}\min }\rho (z-d)}{\mathbf{\Psi }}_{V\max }^m

(73d) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}{\mathbf{\Psi }}_U(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{aj}\rho (z-H_{j-\mathrm{1}})-{\gamma }_{a(j-\mathrm{1})}\rho h_{(j-\mathrm{1})}-\cdots -{\gamma }_{a(k+\mathrm{1})}\rho h_{(k+\mathrm{1})}-{\gamma }_{ak}\rho (H_k-d)}{\mathbf{\Psi }}_{U\max }^m\le {\rho }^m{e}^{-{\gamma }_{a\min }\rho (z-d)}{\mathbf{\Psi }}_{U\max }^m

where , {\gamma }_{\mathrm{0}\min }=\min ({\gamma }_{\mathrm{00}},{\gamma }_{\mathrm{01}},\cdots ,{\gamma }_{\mathrm{0}n},{\gamma }_{\mathrm{0}(n+\mathrm{1})}) , {\gamma }_{a\min }=\min ({\gamma }_{a\mathrm{0}},{\gamma }_{a\mathrm{1}},\cdots ,{\gamma }_{an},{\gamma }_{a(n+\mathrm{1})}) , m=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\cdots ; {\mathbf{\Psi }}_{V\max }^m and {\mathbf{\Psi }}_{U\max }^m are two positive square matrices and dependent on the material parameters only.

(g) Correspondingly, \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) and their mth differential derivatives with respect to z have the following bounds.

(74a) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}\mathbf{\Phi }(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{\min }\rho \left|z-d\right|}{\mathbf{\Phi }}_{\max }^m

(74b) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}\mathbf{\Psi }(\rho ,z)\right|\le {\rho }^m{e}^{-{\gamma }_{\min }\rho \left|z-d\right|}{\mathbf{\Psi }}_{\max }^m

where , {\gamma }_{\min }=\min ({\gamma }_{\mathrm{0}\min },{\gamma }_{a\min }) , m=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\cdots ; {\mathbf{\Phi }}_{\max }^m and {\mathbf{\Psi }}_{\max }^m are two positive square matrices and dependent on the material parameters only.

(h) Finally, as \left|z-d\right|\to \mathrm{0} and \rho \to M (M is a very large value approaching +\infty ), {\mathbf{\Psi }}_V(\rho ,z) and {\mathbf{\Psi }}_U(\rho ,z) have the following asymptotic solutions.

(75a) \underset{\left|z-d\right|\to \mathrm{0}}{\lim }{\mathbf{\Psi }}_V(\rho ,z)\approx {\mathbf{\Psi }}_V^{two}(\rho (z-d))

(75b) \underset{\left|z-d\right|\to \mathrm{0}}{\lim }{\mathbf{\Psi }}_U(\rho ,z)\approx {\mathbf{\Psi }}_U^{two}(\rho (z-d))

where {\mathbf{\Psi }}_V^{two}(\rho (z-d)) and {\mathbf{\Psi }}_U^{two}(\rho (z-d)) are the matrix solutions for the body force vector acting in the interior of two perfectly bonded elastic halfspace with the material properties of the (k-\mathrm{1})\mathrm{th} and k\mathrm{th} layers for d-H_{k-\mathrm{1}}\ge H_k-d or the k\mathrm{th} and (k+\mathrm{1})\mathrm{th} layers for . They can also be obtained from the solution in Eqs. (10)–(17) using the following method. (a) for d-H_{k-\mathrm{1}}\ge H_k-d , letting c_{lj}=c_{l(k-\mathrm{1})} where l=\mathrm{1}\mathrm{,}\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5} and j=\mathrm{0},\mathrm{1},\cdots ,k-\mathrm{2} and c_{lj}=c_{lk} where l=\mathrm{1},2,\mathrm{3},\mathrm{4},\mathrm{5} and j=n+\mathrm{1},n,\cdots ,k+\mathrm{1} ; (b) for , letting c_{ij}=c_{ik} where i=\mathrm{1},2,\mathrm{3},\mathrm{4},\mathrm{5} and j=\mathrm{0},\mathrm{1},\cdots ,k and c_{ij}=c_{i(k+\mathrm{1})} where i=\mathrm{1},2,\mathrm{3},\mathrm{4},\mathrm{5} and j=n+\mathrm{1},n,\cdots ,k+\mathrm{2} .

Convergence of the solution

Without loss of generality, the internal loading of the body force vector in the transform domain \tilde{\mathbf{f}}(\xi ,\eta ) can be assumed to be a normal function of the independent variables (\xi ,\eta ) . It has no singularity for any values of (\xi ,\eta ) in and approaches to constant values as \rho \to +\infty . In other words,

(76) \left|\tilde{f}(\xi ,\eta )\right|\le F_{\mathrm{0}}{I}_a

where ; F_{\mathrm{0}}=a\mathrm{constant} ; I_a={(\mathrm{1},\mathrm{1},\mathrm{1})}^{\mathrm{T}} .

Accordingly, for a given small positive value \epsilon (\epsilon >\mathrm{0}) , under the condition \left|z-d\right|\ge \epsilon >\mathrm{0} , it can be shown that the solutions of \mathbf{u}(x,y,z) , {\mathbf{T}}_z(x,y,z) and {\mathbf{\Gamma }}_p(x,y,z) as well as their partial derivatives with respect to x, y and z can be expressed in the following double improper integrals depending on parameters. The integrals are uniformly and absolutely convergent.

(77a)

(77b)

(77c)

where the integers k,l,m\ge \mathrm{0} , , and if z=H_j , the partial derivatives have the values in the sense of either z=H_j^{-} or z=H_j^{+} , J=\mathrm{0},\mathrm{1},\mathrm{2},\cdots ,n,n+\mathrm{1.}

Secondly, under the condition , the convergence of the solutions and their partial derivatives with respect to x, y and z can be examined using the asymptotic solutions in Eqs. (72). In general, the solutions can be re-expressed as follows.

(78a) \mathbf{u}(x,y,z)={\mathbf{u}}^a(x,y,z)+{\mathbf{u}}^{two}(x,y,z)

(78b) {\mathbf{T}}_z(x,y,z)={\mathbf{T}}_z^a(x,y,z)+{\mathbf{T}}_z^{two}(x,y,z)

{\mathbf{\Gamma }}_p(x,y,z)={\mathbf{\Gamma }}_p^a(x,y,z)+{\mathbf{\Gamma }}_p^a(x,y,z)

(where {\mathbf{u}}^a(x,y,z) , {\mathbf{T}}_z^a(x,y,z) , and {\mathbf{\Gamma }}_p^a(x,y,z) are expressed as follows.

(79a) {\mathbf{u}}^a(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}\mathbf{\Pi }[\mathbf{\Phi }(\rho ,z)-{\mathbf{\Phi }}^{two}(\rho (z-d))]{\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(79b) {\mathbf{T}}_z^a(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }[\mathrm{\Psi }(\rho ,z)-{\mathrm{\Psi }}^{two}(\rho (z-d))]{\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(79c) {\mathbf{\Gamma }}_p^a(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}_p[\mathrm{\Phi }(\rho ,z)-{\mathrm{\Phi }}^{two}(\rho (z-d))]{\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

It can be shown that [\mathbf{\Phi }(\rho ,z)-{\mathbf{\Phi }}^{two}(\rho ,z)] and [\mathbf{\Psi }(\rho ,z)-{\mathbf{\Psi }}^{two}(\rho ,z)] have the following bounds

(80a) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}[\mathbf{\Phi }(\rho ,z)-{\mathbf{\Phi }}^{two}(\rho (z-d))]\right|\le {\rho }^m{e}^{-{\gamma }_{\min }\rho \epsilon }{\mathbf{\Phi }}_{\mathrm{M}\mathrm{a}\mathrm{x}}^m

(80b) \left|{\displaystyle \frac{{\partial }^m}{\partial z^m}}[\mathbf{\Psi }(\rho ,z)-{\mathbf{\Psi }}^{two}(\rho (z-d))]\right|\le {\rho }^m{e}^{-{\gamma }_{\min }\rho \epsilon }{\mathbf{\Psi }}_{\mathrm{M}\mathrm{a}\mathrm{x}}^m

where , m=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\cdots , {\mathbf{\Phi }}_{\mathrm{M}\mathrm{a}\mathrm{x}}^m and {\mathbf{\Psi }}_{\mathrm{M}\mathrm{a}\mathrm{x}}^m are two positive square matrices and dependent on the material parameters only.

Accordingly, {\mathbf{u}}^a(x,y,z) , {\mathbf{T}}_z^a(x,y,z) and {\mathbf{\Gamma }}_p^a(x,y,z) as well as their partial derivatives with respect to x, y and z in the double improper integrals depending on parameters are uniformly and absolutely convergent. The uniform and absolute convergence of the double improper integrals in Eqs. (77) and (79) can be shown using the following result of the bounds:

(81) {\displaystyle \frac{F_{\mathrm{0}}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rho }^{k+l+m}{e}^{-{\gamma }_{\min }\epsilon \rho }\mathrm{d}\xi \mathrm{d}\eta =F_{\mathrm{0}}}{\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\rho }^{k+l+m}{e}^{-{\gamma }_{\min }\epsilon \rho }\mathrm{d}\rho }={\displaystyle \frac{(k+l+m)!F_{\mathrm{0}}}{{({\gamma }_{\min }\epsilon )}^{k+l+m}}}

Thirdly, {\mathbf{u}}^{two}(x,y,z) , {\mathbf{T}}_z^{two}(x,y,z) , and {\mathbf{\Gamma }}_p^{two}(x,y,z) are expressed as follows.

(82a) {\mathbf{u}}^{two}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}\mathbf{\Pi }{\mathbf{\Phi }}^{two}(\rho (z-d)){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(82b) {\mathbf{T}}_z^{two}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }{\mathbf{\Psi }}^{two}(\rho (z-d)){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

(82c) {\mathbf{\Gamma }}_p^{two}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}_p{\mathbf{\Phi }}^{two}(\rho (z-d)){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

where {\mathbf{\Phi }}^{two}(\rho (z-d)) and {\mathbf{\Psi }}^{two}(\rho (z-d)) are the matrix solutions for the body force vector acting in the interior of two bonded elastic halfspace with the material properties of the (k-\mathrm{1})\mathrm{th} and k\mathrm{th} layers for d-H_{k-\mathrm{1}}\ge H_k-d or the k\mathrm{th} and (k+\mathrm{1})\mathrm{th} layers for .

The improper integrals in Eq. (82) are convergent in the sense of the Cauchy principal value. It can be shown that as \rho \to M (M is a very large value approaching +\infty and \left|z-d\right|\to \mathrm{0} , \mathbf{\Phi }(\rho ,z) and \mathbf{\Psi }(\rho ,z) have asymptotic solutions and can be expressed as follows.

(83a) \underset{\rho \to M}{\lim }\underset{\left|z-d\right|\to \mathrm{0}}{\lim }\mathbf{\Phi }(\rho ,z)={\mathbf{\Phi }}^{two}(\rho (z-d))

(83b) \underset{\rho \to M}{\lim }\underset{\left|z-d\right|\to \mathrm{0}}{\lim }\mathbf{\Psi }(\rho ,z)={\mathbf{\Psi }}^{two}(\rho (z-d))

Closed-form singularity of the solution

According to Eqs. (79) – (82), it is evident that the singularity of the solutions in the layered heterogeneous solid is exactly the same as that of the solutions for the body force vector \mathbf{f}(x,y)\delta (z-d) concentrated at a horizontal plane in the interior of two perfectly bonded elastic halfspaces. The upper halfspace has the material properties of the (k-\mathrm{1})\mathrm{th} layer and the lower halfspace has the material properties of the k\mathrm{th} layer, under the condition of d-H_{k-\mathrm{1}}\ge H_k-d (i.e., the loading plane z=d is closer to the interface of the (k-\mathrm{1})\mathrm{th} and k\mathrm{th} layers). On the other hand, the upper halfspace has the material properties of the k\mathrm{th} layer and the lower halfspace has the material properties of the (k+\mathrm{1})\mathrm{th} layer, under the condition of (i.e., the loading plane z=d is closer to the interface of the k\mathrm{th} and (k+\mathrm{1})\mathrm{th} layers).

{\mathbf{\Phi }}^{two}(\rho (z-d)) and {\mathbf{\Psi }}^{two}(\rho (z-d)) have been given exactly given in (10) to (17) and Yue [ 2]. The solutions of {\mathbf{u}}^a(x,y,z) , {\mathbf{T}}_z^a(x,y,z) and {\mathbf{\Gamma }}_p^a(x,y,z) in Eq. (82) have been derived systematically and exactly in the closed-form in terms of elementary harmonic functions or the complete elliptic integrals of the first, second and third kinds for the three concentrated body force vectors (see the above Section 2 and Refs. [ 2, 3, 6]).

Satisfaction of the required conditions

Under the condition of \left|z-d\right|\ge \epsilon >\mathrm{0} , the improper integrals of the solution and their partial differentiations with respect to x, y and z are uniformly and absolutely convergent. Consequently, the integration limits and the integrations and the partial differentiations for the solution in the forms of the improper integrals are interchangeable. Accordingly, it can be shown that the solution given in Eqs. (1) satisfies the governing partial differential equations using the following method. Based on the geometric Eqs. (2), the static equilibrium Eqs. (4) and the Hooke’s law (6) in Ref. [1], the governing partial differential equations for the j-th layer can be re-expressed in terms of the displacement vector u as follows.

(84a) {\mathbf{L}}_j\mathbf{u}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}{\mathbf{L}}_j[\mathbf{\Pi }\mathbf{\Phi }(\rho ,z)K]\mathbf{g}(\xi ,\eta )\mathrm{d}\xi \mathrm{d}\eta }=\mathbf{0}

where and for j=\mathrm{0},\mathrm{1},\cdots ,n,n+1\mathrm{and}j\ne k or and for j=k ;

(84b) {\mathbf{L}}_j=\left(\begin{array}{ccc}{c}_{\mathrm{1}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+c_{\mathrm{5}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}+c_{\mathrm{4}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}& (c_{\mathrm{1}j}-c_{\mathrm{5}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial y}}& (c_{\mathrm{2}j}+c_{\mathrm{4}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial y}}\\ (c_{\mathrm{1}j}-c_{\mathrm{5}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial y}}& c_{\mathrm{5}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+c_{\mathrm{1}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}+c_{\mathrm{4}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}& (c_{\mathrm{2}j}+c_{\mathrm{4}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y\partial z}}\\ (c_{\mathrm{2}j}+c_{\mathrm{4}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial z}}& (c_{\mathrm{2}j}+c_{\mathrm{4}j}){\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y\partial z}}& c_{\mathrm{4}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+c_{\mathrm{4}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}+c_{\mathrm{3}j}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\end{array}\right)

The Eq. (84a) is valid because of the following identify.

(85) {\mathbf{L}}_j[\mathbf{\Pi }\mathbf{\Phi }(\rho ,z)K]=\mathbf{0}

Secondly, the solution in Eqs. (1) satisfies Hooke’s law and the geometric equations within each single layer because of the following two identifies.

(86a) {\mathbf{L}}_{pj}\mathbf{u}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}{\mathbf{L}}_{pj}[\mathbf{\Pi }\mathbf{\Phi }(\rho ,z)K]\mathbf{g}(\xi ,\eta )\mathrm{d}\xi \mathrm{d}\eta }={\mathbf{\Gamma }}_p(x,y,z)

(86b) {\mathbf{L}}_{zj}\mathbf{u}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}{\mathbf{L}}_{zj}[\mathbf{\Pi }\mathbf{\Phi }(\rho ,z)K]\mathbf{g}(\xi ,\eta )\mathrm{d}\xi \mathrm{d}\eta }={\mathbf{T}}_z(x,y,z)

where and for j=\mathrm{0},\mathrm{1},\mathrm{..},n,n+1\mathrm{and}j\ne k or and for j=k ;

(86c) {\mathbf{L}}_{pj}=\left(\begin{array}{ccc}{\displaystyle \frac{\partial }{\partial x}}& \mathrm{0}& \mathrm{0}\\ {\displaystyle \frac{\partial }{\mathrm{2}\partial y}}& {\displaystyle \frac{\partial }{\mathrm{2}\partial x}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\partial }{\partial y}}& \mathrm{0}\end{array}\right);{\mathbf{L}}_{zj}=\left(\begin{array}{ccc}{c}_{\mathrm{4}j}{\displaystyle \frac{\partial }{\partial z}}& \mathrm{0}& c_{\mathrm{4}j}{\displaystyle \frac{\partial }{\partial x}}\\ \mathrm{0}& c_{\mathrm{4}j}{\displaystyle \frac{\partial }{\partial z}}& c_{\mathrm{4}j}{\displaystyle \frac{\partial }{\partial y}}\\ c_{\mathrm{2}j}{\displaystyle \frac{\partial }{\partial x}}& c_{\mathrm{2}j}{\displaystyle \frac{\partial }{\partial y}}& c_{\mathrm{3}j}{\displaystyle \frac{\partial }{\partial z}}\end{array}\right)

The two identifies in Eqs. (86a) and (86b) are valid because of the following two identifies.

(87a) {\mathbf{L}}_{pj}[\mathbf{\Pi }K]=\rho {\mathbf{\Pi }}_{p}K

(87b) {\mathbf{L}}_{zj}[\mathbf{\Pi }\mathbf{\Phi }(\rho ,z)K]=\rho \mathbf{\Pi }\mathbf{\Psi }(\rho ,z)K

Thirdly, the solution of \mathbf{u}(x,y,z) , {\mathbf{T}}_z(x,y,z) and {\mathbf{\Gamma }}_p(x,y,z) in Eqs. (1) satisfies the perfectly bonded interface condition. The solution of \mathbf{u}(x,y,z) and {\mathbf{\Gamma }}_p(x,y,z) in Eqs. (1) is continuous functions of x,y \mathrm{and} z , where . The solution of {\mathbf{T}}_z(x,y,z) in Eq. (1) is continuous functions of x,y\mathrm{and}z , where and or .

Fourthly, the solution of {\mathbf{T}}_z(x,y,z) satisfies the internal loading condition of the concentrated body force \mathbf{f}(x,y)\delta (z-d) because the following identity. In other word, the discontinuity Eqs. (71) is due to the body force loading concentrated at z=d .

(88) \underset{\mathrm{\Delta }z\to +\mathrm{0}}{\lim }[{\mathbf{T}}_z(x,y,z+\mathrm{\Delta }z)-{\mathbf{T}}_z(x,y,z-\mathrm{\Delta }z)]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }\underset{\mathrm{\Delta }z\to +\mathrm{0}}{\lim }[\mathbf{\Psi }(\rho ,z+\mathrm{\Delta }z)-\mathbf{\Psi }(\rho ,z-\mathrm{\Delta }z)]{\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }=\{\begin{array}{ll}\mathrm{0}& for z\ne d\\ -\mathbf{f}(x,y)& for z=d\end{array}

Finally, the solution of \mathbf{u}(x,y,z) , {\mathbf{T}}_z(x,y,z) and {\mathbf{\Gamma }}_p(x,y,z) in Eqs. (1) satisfies the natural regularity conditions as R\to +\infty , where R=\sqrt{r^{\mathrm{2}}+{(z-d)}^{\mathrm{2}}} and r=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}} In other words, it can be shown the following limiting results are valid.

(89a) \underset{R\to +\infty }{\lim }\mathbf{u}(x,y,z)=\{\begin{array}{l}\underset{\left|z-d\right|\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}\mathbf{\Pi }\mathbf{\Phi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }\right]\\ \underset{r\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }r}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{\mathrm{1}}{\rho }}\mathbf{\Pi }\mathbf{\Phi }\left({\displaystyle \frac{\rho }{r}},z\right){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}\left({\displaystyle \frac{\xi }{r}},{\displaystyle \frac{\eta }{r}}\right)K_{\mathrm{0}}\mathrm{d}\xi \mathrm{d}\eta }\right]\end{array}=\{\begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}

(89b) \underset{R\to +\infty }{\lim }{\mathbf{T}}_z(x,y,z)=\{\begin{array}{l}\underset{\left|z-d\right|\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }\mathbf{\Psi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }\right]\\ \underset{r\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }{r}^{\mathrm{2}}}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathbf{\Pi }\mathbf{\Psi }\left({\displaystyle \frac{\rho }{r}},z\right){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}\left({\displaystyle \frac{\xi }{r}},{\displaystyle \frac{\eta }{r}}\right)K_{\mathrm{0}}\mathrm{d}\xi \mathrm{d}\eta }\right]\end{array}=\{\begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}

(89c) \underset{R\to +\infty }{\lim }{\mathbf{\Gamma }}_p(x,y,z)=\{\begin{array}{l}\underset{\left|z-d\right|\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}_p\mathbf{\Phi }(\rho ,z){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }\right]\\ \underset{r\to +\infty }{\lim }\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }{r}^{\mathrm{2}}}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}_p\mathbf{\Phi }\left({\displaystyle \frac{\rho }{r}},z\right){\mathbf{\Pi }}^{*}\tilde{\mathbf{f}}\left({\displaystyle \frac{\xi }{r}},{\displaystyle \frac{\eta }{r}}\right)K_{\mathrm{0}}\mathrm{d}\xi \mathrm{d}\eta }\right]\end{array}=\{\begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}

where K_{\mathrm{0}}=e^{i(\xi \cos \theta +\eta \sin \theta )} and the following inequalities can be shown.

(90a)

(90b)

(90c)

where F_{\mathrm{0}}>\mathrm{0} .

It can be shown that the following bounds are valid with the constant parameters,

(91a)

(91b)

(91c)

where {\mathbf{\Phi }}_{\max }^{\mathrm{0}} and {\mathbf{\Psi }}_{\max }^{\mathrm{0}} are constant square matrices with positive elements depending on the material parameters.

On the other hand, it can be shown that the following integrations are convergent in the sense of the Cauchy principal value.

(92a)

(92b)

(92c)

where \mathbf{\Phi }(\mathrm{0},z) and \mathbf{\Psi }(\mathrm{0},z) are two constant matrices.

Summary notes

The solution given in Eqs. (82) is uniformly and absolutely convergent provided that or . If z=d^{-} or z=d^{+} , the solution is convergent in the sense of Cauchy principal value, where z=d is the loading plane of the concentrated body force vector. The solution satisfies the governing partial differential equations, the interfacial connection conditions, the internal loading condition and the external regularity conditions as R\to +\infty . The singularity of the solution occurs only at the loading plane location and can be exactly presented in closed-form in terms of elementary harmonic functions and the special functions. It is also noted that the above mathematical properties of the solution can be examined and are also valid in the cylindrical coordinate system.

Applications to other solutions

General

The classical elasticity theory is the common foundation of many other continuum mechanics. They include contact mechanics, fracture mechanics, inclusion mechanics, elastodynamics, thermo-elasticity, thermo-elastodynamics, poroelasticity and viscoelasticity. The matrix Fourier integral transform approach presented above can be applied to derive and formulate analytical solutions for many other boundary-value and initial-boundary-value problems of linear continuum mechanics [ 9− 29]. For the initial-boundary-value problems, the Laplace integral transform method has to be applied to the independent variable of time t . Some cases and examples are further discussed in this section.

Other boundary-value problems in n-layered solids

Solutions in n-layered solids of semi-infinite extent

Similar to the solution in the n-layered solid of infinite extent, the solutions can be derived and formulated systematically for the boundary-value problems in a general layered elastic solid of semi-infinite extent . The general layered solid has (n + 1) dissimilar layers. For H_{j-\mathrm{1}}^{+}\le z\le H_j^{-} , it is the jth homogeneous elastic layer with the layer thickness h_j=H_j^{+}-H_{j-\mathrm{1}}^{-} (j=\mathrm{1},\mathrm{2},\mathrm{3},\cdots ,n) . For , it is the (n + 1)th homogeneous elastic halfspace. The five elastic constants for the jth layer are c_{\mathrm{1}j},c_{\mathrm{2}j},c_{\mathrm{3}j},c_{\mathrm{4}j},c_{\mathrm{5}j} ( j=\mathrm{1},\mathrm{2},\mathrm{3},\cdots ,n,n+\mathrm{1} ).

In addition to the internal loading of the body force vector \mathbf{f}(x,y,z)=\mathbf{f}(x,y)\delta (z-d) ( ), the external boundary conditions on the external boundary z=H_{\mathrm{0}} can be described as one of the following four conditions.

(93a) {\mathbf{u}}^a=\mathbf{u}(x,y,H_{\mathrm{0}})=\left(\begin{array}{l}{u}_x(x,y,H_{\mathrm{0}})\\ u_y(x,y,H_{\mathrm{0}})\\ u_z(x,y,H_{\mathrm{0}})\end{array}\right)

(93b) {\mathbf{T}}_z^a={\mathbf{T}}_z(x,y,H_{\mathrm{0}})= \left(\begin{array}{l}{\sigma }_{xz}(x,y,H_{\mathrm{0}})\\ {\sigma }_{yz}(x,y,H_{\mathrm{0}})\\ {\sigma }_{zz}(x,y,H_{\mathrm{0}})\end{array}\right)

(93c) {\mathbf{T}}_{zu}^a={\mathbf{T}}_{zu}(x,y,H_{\mathrm{0}})= \left(\begin{array}{l}{u}_x(x,y,H_{\mathrm{0}})\\ u_y(x,y,H_{\mathrm{0}})\\ {\sigma }_{zz}(x,y,H_{\mathrm{0}})\end{array}\right)

(93d) {\mathbf{T}}_{uz}^a={\mathbf{T}}_{uz}(x,y,H_{\mathrm{0}})= \left(\begin{array}{l}{\sigma }_{xz}(x,y,H_{\mathrm{0}})\\ {\sigma }_{yz}(x,y,H_{\mathrm{0}})\\ u_z(x,y,H_{\mathrm{0}})\end{array}\right)

The corresponding loading conditions in the transform domain can be expressed by one of the following four loading vectors at z=H_{\mathrm{0}} and plus the internal body force vector \mathbf{g} .

(94a) {\mathbf{w}}^a=\left(\begin{array}{l}{w}_{\mathrm{1}}^a\\ w_{\mathrm{2}}^a\\ w_{\mathrm{3}a}^a\end{array}\right), {\mathbf{Y}}_z^a=\left(\begin{array}{l}{\tau }_{\mathrm{1}}^a\\ {\tau }_{\mathrm{2}}^a\\ {\tau }_{\mathrm{3}}^a\end{array}\right), {\mathbf{Y}}_{zu}^a=\left(\begin{array}{l}{w}_{\mathrm{1}}^a\\ w_{\mathrm{2}}^a\\ {\tau }_{\mathrm{3}}^a\end{array}\right), {\mathbf{Y}}_{uz}^a=\left(\begin{array}{l}{\tau }_{\mathrm{1}}^a\\ {\tau }_{\mathrm{2}}^a\\ w_{\mathrm{3}}^a\end{array}\right)

where

(94b) {\mathbf{w}}^a={\mathbf{w}}^a(\xi ,\eta ,H_{\mathrm{0}})={\displaystyle \frac{\rho }{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}^{*}{\mathbf{u}}^a(x,y)K^{*}\mathrm{d}x\mathrm{d}y}

(94c) {\mathbf{Y}}_z^a={\mathbf{Y}}_z^a(\xi ,\eta ,H_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{\Pi }}^{*}{\mathbf{T}}_z^a(x,y)K^{*}\mathrm{d}x\mathrm{d}y}

(94d) {\mathbf{Y}}_{zu}^a={\mathbf{Y}}_{zu}^a(\xi ,\eta ,H_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{I}}_a(\rho ){\mathbf{\Pi }}^{*}{\mathbf{T}}_{zu}^a(x,y)K^{*}\mathrm{d}x\mathrm{d}y}

(94e) {\mathbf{Y}}_{uz}^a={\mathbf{Y}}_{uz}^a(\xi ,\eta ,H_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\mathbf{I}}_b(\rho ){\mathbf{\Pi }}^{*}{\mathbf{T}}_{uz}^a(x,y)K^{*}\mathrm{d}x\mathrm{d}y}

where

(94f) \mathbf{I}(\rho )_a=\left(\begin{array}{lll}{\displaystyle \frac{\mathrm{1}}{\rho }}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{\rho }}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right),{\mathbf{I}}_b(\rho )=\left(\begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{\rho }}\end{array}\right)

For each of the four boundary loading cases, the solution can be obtained using the same equations and the formulation procedure given in Sections 4 and 5 in [ 1] and the above Sections 2 and 3 except the three regularity boundary conditions at z=H_{\mathrm{0}} in Eqs. (50) and (68) in Ref. [ 1]. They have to be replaced by the given three boundary conditions at z=H_{\mathrm{0}} as follows, respectively.

(95a) {\mathbf{q}}_a\mathbf{V}(H_{\mathrm{0}})={\mathbf{q}}_a{\mathbf{V}}^a

(95b) {\mathbf{P}}_{qa}\mathbf{U}(H_{\mathrm{0}})={\mathbf{P}}_{qa}{\mathbf{U}}^a

where

(95c) {\mathbf{V}}^a=\left(\begin{array}{l}{w}_{\mathrm{2}}^a\\ {\tau }_{\mathrm{2}}^a\end{array}\right);{\mathbf{q}}_a=\{\begin{array}{ll}\left[\begin{array}{ll}\mathrm{1}& \mathrm{0}\end{array}\right],& \mathrm{for} {\mathbf{u}}^a \mathrm{or} {\mathbf{T}}_{zu}^a \mathrm{is} known \\ \left[\begin{array}{ll}\mathrm{0}& \mathrm{1}\end{array}\right],& \mathrm{for} {\mathbf{T}}_z^a \mathrm{or} {\mathbf{T}}_{uz}^a \mathrm{is} \mathrm{known}\end{array}

(95d) {\mathbf{U}}^a=\left(\begin{array}{l}{w}_{\mathrm{1}}^a\\ w_{\mathrm{3}}^a\\ {\tau }_{\mathrm{3}}^a\\ {\tau }_{\mathrm{1}}^a\end{array}\right);{\mathbf{P}}_{qa}=\{\begin{array}{ll}\left[\begin{array}{llll}\begin{array}{l}\mathrm{1}\\ \mathrm{0}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{1}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right]& \mathrm{for} {\mathbf{u}}^a \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]& \mathrm{for} {\mathbf{T}}_{zu}^a \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]& \mathrm{for} {\mathbf{T}}_{uz}^a \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]& \mathrm{for} {\mathbf{T}}_z^a \mathrm{is} \mathrm{known}\end{array}

For each boundary condition case, two sets of the solution can be obtained. One set is for the internal loading of the body force vector \mathbf{g} in the transform domain or \mathbf{f} in the physical domain. The other set is for each of {\mathbf{w}}^a , {\mathbf{Y}}_z^a , {\mathbf{Y}}_{zu}^a and {\mathbf{Y}}_{uz}^a in the transform domain or each of {\mathbf{u}}^a , {\mathbf{T}}_z^a , {\mathbf{T}}_{zu}^a and {\mathbf{T}}_{uz}^a in the physical domain. More details can be found in Yue [ 8], Yue and Yin [ 12] and Yue et al. [ 13]

Solutions in n-layered solid of finite thickness

Similarly, the solutions in general layered solid of finite thickness (H_{\mathrm{0}}\le z\le H_n) can be systematically derived and formulated. The solid has n dissimilar layers. For H_{j-\mathrm{1}}^{+}\le z\le H_j^{-} it is the jth homogeneous elastic layer with the layer thickness h_j=H_j^{+}-H_{j-\mathrm{1}}^{-} and the five elastic constants c_{\mathrm{1}j},c_{\mathrm{2}j},c_{\mathrm{3}j},c_{\mathrm{4}j},c_{\mathrm{5}j} ( j=\mathrm{1},\mathrm{2},\mathrm{3},\cdots ,n ). In addition to the internal loading of the body force vector and the external boundary conditions on the external boundary z=H_{\mathrm{0}} in Eq. (93), the other boundary conditions on the external boundary z=H_n can be described as one of the following four conditions.

(96a) {\mathbf{u}}^b=\mathbf{u}(x,y,H_n)=\left(\begin{array}{l}{u}_x(x,y,H_n)\\ u_y(x,y,H_n)\\ u_z(x,y,H_n)\end{array}\right)

(96b) {\mathbf{T}}_z^b={\mathbf{T}}_z(x,y,H_n)= \left(\begin{array}{l}{\sigma }_{xz}(x,y,H_n)\\ {\sigma }_{yz}(x,y,H_n)\\ {\sigma }_{zz}(x,y,H_n)\end{array}\right)

(96c) {\mathbf{T}}_{zu}^b={\mathbf{T}}_{zu}(x,y,H_n)= \left(\begin{array}{l}{u}_x(x,y,H_n)\\ u_y(x,y,H_n)\\ {\sigma }_{zz}(x,y,H_n)\end{array}\right)

(96d) {\mathbf{T}}_{uz}^b={\mathbf{T}}_{uz}(x,y,H_n)= \left(\begin{array}{l}{\sigma }_{xz}(x,y,H_n)\\ {\sigma }_{yz}(x,y,H_n)\\ u_z(x,y,H_n)\end{array}\right)

The corresponding loading conditions in the transform domain can be expressed by one of the following four loading vectors at z=H_n and plus the internal body force vector \mathbf{g} .

(97) {\mathbf{w}}^b=\left(\begin{array}{l}{w}_{\mathrm{1}}^b\\ w_{\mathrm{2}}^b\\ w_{\mathrm{3}}^b\end{array}\right), {\mathbf{Y}}_z^b=\left(\begin{array}{l}{\tau }_{\mathrm{1}}^b\\ {\tau }_{\mathrm{2}}^b\\ {\tau }_{\mathrm{3}}^b\end{array}\right), {\mathbf{Y}}_{zu}^b=\left(\begin{array}{l}{w}_{\mathrm{1}}^b\\ w_{\mathrm{2}}^b\\ {\tau }_{\mathrm{3}}^b\end{array}\right), {\mathbf{Y}}_{uz}^b=\left(\begin{array}{l}{\tau }_{\mathrm{1}}^b\\ {\tau }_{\mathrm{2}}^b\\ w_{\mathrm{3}}^b\end{array}\right)

where the given boundary conditions {\mathbf{w}}^b, {\mathbf{Y}}_z^b, {\mathbf{Y}}_{zu}^b, {\mathbf{Y}}_{uz}^b can be obtained by changing the superscript a with b in Eq. (97).

For each of the four boundary loading cases, the solution can be obtained using the same equations and the formulation procedure given in Sections 4 and 5 in [ 1] and the above Sections 2 and 3 except the following changes. The three regularity boundary conditions at z=H_{\mathrm{0}} in Eqs. (50) and (68) in Ref. [ 1] have to be replaced by Eq. (95). The three regularity boundary conditions at z=H_n in Eqs. (51) and (69) in Ref. [ 1] have to be replaced by the following three given boundary conditions at z=H_n .

(98a) {\mathbf{p}}_b\mathbf{V}(H_n)={\mathbf{p}}_b{\mathbf{V}}^b

(98b) {\mathbf{P}}_{pb}\mathbf{U}(H_n)={\mathbf{P}}_{pb}{\mathbf{U}}^b

where

(98c) {\mathbf{V}}^b=\left(\begin{array}{l}{w}_{\mathrm{2}}^b\\ {\tau }_{\mathrm{2}}^b\end{array}\right);{\mathbf{p}}_b=\{\begin{array}{ll}\left[\begin{array}{ll}\mathrm{1}& \mathrm{0}\end{array}\right]& \mathrm{for} {\mathbf{u}}^b \mathrm{or} {\mathbf{T}}_{zu}^b \mathrm{is} known \\ \left[\begin{array}{ll}\mathrm{0}& \mathrm{1}\end{array}\right]& \mathrm{for} {\mathbf{T}}_z^b \mathrm{or} {\mathbf{T}}_{uz}^b \mathrm{is} \mathrm{known}\end{array}

(98d) U^b=\left(\begin{array}{l}{w}_{\mathrm{1}}^b\\ w_{\mathrm{3}}^b\\ {\tau }_{\mathrm{3}}^b\\ {\tau }_{\mathrm{1}}^b\end{array}\right);P_{pb}=\{\begin{array}{ll}\left[\begin{array}{llll}\begin{array}{l}\mathrm{1}\\ \mathrm{0}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{1}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}& \begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right],& \mathrm{for} u^b \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right],& \mathrm{for} T_{zu}^b \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],& \mathrm{for} T_{uz}^b \mathrm{is} \mathrm{known}\\ \left[\begin{array}{llll}\mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],& \mathrm{for} T_z^b \mathrm{is} \mathrm{known}\end{array}

Each of the two boundaries at z=H_{\mathrm{0}} and z=H_n has four types of three given conditions. Hence, there are 16 cases of the boundary-value problems. Each case has three sets of solutions in terms of the body force vector, the three given conditions at z=H_{\mathrm{0}} and the given three conditions at z=H_n . At z=H_{\mathrm{0}} , the three given conditions are either {\mathbf{w}}^a , {\mathbf{Y}}_z^a , {\mathbf{Y}}_{zu}^a or {\mathbf{Y}}_{uz}^a in the transform domain and either {\mathbf{u}}^a , {\mathbf{T}}_z^a , {\mathbf{T}}_{zu}^a or {\mathbf{T}}_{uz}^a in the physical domain. At z=H_n , the three given conditions are either {\mathbf{w}}^b , {\mathbf{Y}}_z^b , {\mathbf{Y}}_{zu}^b or {\mathbf{Y}}_{uz}^b in the transform domain and either {\mathbf{u}}^b , {\mathbf{T}}_z^b , {\mathbf{T}}_{zu}^b or {\mathbf{T}}_{uz}^b in the physical domain.

Mix-boundary value problems

Mix-boundary value problems in elasticity are usually called contact mechanics problems and fracture mechanics. Two methods can be used to derive and formulate solutions for the mixed-boundary value problem. The first method is the integral equation method [ 9, 30, 31]. The second method is the boundary element method [ 14− 19]. They are discussed below.

Integral equation method

An example is given to illustrate how to use the matrix Fourier integral approach to solve the mix-boundary value problems in the general layered solids [ 9]. This mixed boundary value problem is a layered solid of semi-infinite extent subjected to the eccentric indentation of a rigid circular and smooth plate. The mixed boundary conditions at the surface z=H_{\mathrm{0}} of the layered solid are written in the cylindrical coordinate system as follows.

(99a)

(99b)

(99c)

(99d)

where a is the radius of the rigid circular plate; D_z is the axial translation of the rigid plate along the z-axis; {\mathrm{\Omega }}_y is the central rotation of the rigid plate about the y-axis.

The axial load P_z and its associated moment M_y acting on the rigid plate have the following integral relations with the contact normal stress {\sigma }_{zz}(r,\theta ,H_{\mathrm{0}}) .

(100a) P_z=-{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\pi }}{\int }}{\sigma }_{zz}(r,\theta ,H_{\mathrm{0}})r\mathrm{d}\theta \mathrm{d}r}}

(100b) M_y=P_{z}b=-{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\pi }}{\int }}{\sigma }_{zz}(r,\theta ,H_{\mathrm{0}})r^{\mathrm{2}}\cos \theta \mathrm{d}\theta \mathrm{d}r}}

where b is the eccentricity of the external load P_z acting on the rigid circular plate.

Using the approach presented above, the \mathbf{w}(z) and {\mathbf{Y}}_z(z) can be expressed as follows in terms of the boundary loading vector {\mathbf{Y}}_z^a .

(101a) \mathbf{w}(\rho ,\phi ,z)=\mathbf{\Phi }(\rho ,z){\mathbf{Y}}_z^a

(101b) {\mathbf{Y}}_z(\rho ,\phi ,,z)=\mathbf{\Psi }(\rho ,z){\mathbf{Y}}_z^a

Because of Eqs. (99c) and (99d), so {\tau }_{\mathrm{1}}^a=\mathrm{0} , and {\tau }_{\mathrm{2}}^a=\mathrm{0} . Then, the following results are valid.

(102a) w_{\mathrm{1}}(\rho ,\phi ,z)={\mathrm{\Phi }}_{\mathrm{13}}(\rho ,z){\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}})

(102b) w_{\mathrm{2}}(\rho ,\phi ,z)=\mathrm{0}

(102c) w_{\mathrm{3}}(\rho ,\phi ,z)={\mathrm{\Phi }}_{\mathrm{33}}(\rho ,z){\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}})

(102d) {\tau }_{\mathrm{1}}(\rho ,\phi ,z)={\mathrm{\Psi }}_{\mathrm{13}}(\rho ,z){\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}})

(102e) {\tau }_{\mathrm{2}}(\rho ,\phi ,z)=\mathrm{0}

(102f) {\tau }_{\mathrm{3}}(\rho ,\phi ,z)={\mathrm{\Psi }}_{\mathrm{33}}(\rho ,z){\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}})

where {\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}}) is the unknown and to be found using the mixed boundary condition Eqs. (99a) and (99b).

Secondly, because of Eqs. (102c) and (102f), the mixed boundary condition Eqs. (99a) and (99b) can be expressed as follows.

(103)

Using the Fourier series expansions of {\tau }_{\mathrm{3}}(\rho ,\phi ,H_{\mathrm{0}})={\tau }_{\mathrm{30}}(\rho )-\mathrm{2}i\sin \phi {\tau }_{\mathrm{31}}(\rho ) , the Eq. (103) can be decoupled into the following two sets of dual integral equations.

(104a)

and

(104b)

Furthermore the following solution representations are used.

(105a) {\tau }_{\mathrm{30}}(\rho )={\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{\phi }_{\mathrm{0}}(s)\cos (\rho s)\mathrm{d}s}

(105b) {\tau }_{\mathrm{31}}(\rho )={\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{\phi }_{\mathrm{1}}(s)\sin (\rho s)\mathrm{d}s}

Then, the Eqs. (104a) and (104b) can be respectively reduced to the following two sets of Fredholm integral equations of the second kind.

(106a) {\phi }_{\mathrm{0}}(s)+{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{K}_a(s,t){\phi }_{\mathrm{0}}(t)\mathrm{d}t={\displaystyle \frac{\mathrm{2}}{\mathrm{\pi }A}}{D}_z}

where \mathrm{0}\le s\le a , A=\underset{\rho \to +\infty }{\lim }{\mathrm{\Phi }}_{\mathrm{33}}(\rho ,\mathrm{0}) and

(106b) K_a(s,t)={\displaystyle \frac{\mathrm{2}}{\mathrm{\pi }}}{\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}\left({\displaystyle \frac{\mathrm{1}}{A}}{\mathrm{\Phi }}_{\mathrm{33}}(\rho ,\mathrm{0})-\mathrm{1}\right)\cos (\rho s)\cos (\rho t)\mathrm{d}\rho }

(106c) P_z=-\mathrm{2}\mathrm{\pi }{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{\phi }_{\mathrm{0}}(s)\mathrm{d}s}

and

(107a) {\phi }_{\mathrm{1}}(s)+{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}{K}_b(s,t){\phi }_{\mathrm{1}}(t)\mathrm{d}t={\displaystyle \frac{\mathrm{2}s}{\mathrm{\pi }A}}{\mathrm{\Omega }}_y}

where \mathrm{0}\le s\le a , and

(107b) K_b(s,t)={\displaystyle \frac{\mathrm{2}}{\mathrm{\pi }}}{\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}\left({\displaystyle \frac{\mathrm{1}}{A}}{\mathrm{\Phi }}_{\mathrm{33}}(\rho ,\mathrm{0})-\mathrm{1}\right)\sin (\rho s)\sin (\rho t)\mathrm{d}\rho }

(107c) M_y=-\mathrm{4}\mathrm{\pi }{\displaystyle \underset{\mathrm{0}}{\overset{a}{\int }}s{\phi }_{\mathrm{1}}(s)\mathrm{d}s}

The above integral equations can be calculated numerically and accurately. The numerical results can be used to find the solution of \mathbf{u} , {\mathbf{T}}_z and {\mathbf{\Gamma }}_p in Eqs. (93) or (100) in Ref. [ 1] with the following solution expressions.

(108a) u(s,\theta ,z_a)={\displaystyle \frac{-P_z}{\mathrm{2}\mathrm{\pi }a}}\left[{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\phi (t){\mathbf{K}}_{uz}(s,z_{\mathrm{1}},t)\mathrm{d}t}+{\displaystyle \frac{b}{\mathrm{2}a}}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\psi (t){\mathbf{K}}_{uy}(s,z_{\mathrm{1}},t)\mathrm{d}t}{\mathbf{\Pi }}_a(\theta )\right]

(108b) {\mathbf{T}}_z(s,\theta ,z_a)={\displaystyle \frac{-P_z}{\mathrm{2}\mathrm{\pi }{a}^{\mathrm{2}}}}\left[{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\phi (t){\mathbf{K}}_{zz}(s,z_{\mathrm{1}},t)\mathrm{d}t}+{\displaystyle \frac{b}{\mathrm{2}a}}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\psi (t){\mathbf{K}}_{zy}(s,z_{\mathrm{1}},t)\mathrm{d}t}{\mathbf{\Pi }}_a(\theta )\right]

(108c) {\mathbf{\Gamma }}_p(s,\theta ,z_a)={\displaystyle \frac{-P_z}{\mathrm{2}\mathrm{\pi }{a}^{\mathrm{2}}}}\left[{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\phi (t){\mathbf{K}}_{pz}(s,z_{\mathrm{1}},t)\mathrm{d}t}+{\displaystyle \frac{b}{\mathrm{2}a}}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\psi (t){\mathbf{K}}_{py}(s,z_{\mathrm{1}},t)\mathrm{d}t}{\mathbf{\Pi }}_a(\theta )\right]

where , , , \phi (s)=-\mathrm{2}\mathrm{\pi }a{\phi }_{\mathrm{0}}(r)/P_z , \psi (s)=-\mathrm{4}\mathrm{\pi }{a}^{\mathrm{2}}{\phi }_{\mathrm{1}}(r)/M_y . {\mathbf{K}}_{uz} , {\mathbf{K}}_{uy} , {\mathbf{K}}_{pz} and {\mathbf{K}}_{py} are functions of {\mathbf{\Phi }}_{\mathrm{13}}(\rho ,z) and {\mathbf{\Phi }}_{\mathrm{33}}(\rho ,z) . {\mathbf{K}}_{zz} and {\mathbf{K}}_{zy} are functions of {\mathbf{\Psi }}_{\mathrm{13}}(\rho ,z) and {\mathbf{\Psi }}_{\mathrm{33}}(\rho ,z) ., which are exactly given in Ref. [ 9].

(108d) {\mathbf{\Pi }}_a(\theta )=\left(\begin{array}{lll}\cos \theta & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sin \theta & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \cos \theta \end{array}\right)

Details of the analytical formulation and numerical solution can be found in Ref. [ 9].

Boundary element methods

The classical closed-form fundamental singular solutions given by Kelvin in 1848 [ 32] and Mindlin in 1936 [ 33] have been used in the formulation and development of the powerful boundary element methods (BEM) [ 19]. These Kelvin’s and Mindlin’s solutions based BEMs have some drawbacks in solving both boundary-value and mixed boundary-value problems in layered and functionally graded materials. Since 2002, Yue’s fundamental singular solution given in Eqs. (84) to (89) in Ref. [ 1] has been used to replace the Kelvin’s and Mindlin’s solutions and formulate the Yue’s solution based BEMs [ 14− 19]. They have been used to find the solutions for many specific problems of interests in science and technology.

The singular terms of the Yue’s solutions are presented in exact closed-form in terms of the elementary functions given in Appendix B and can be analytically isolated and exactly calculated in the new BEMs. The remaining parts of the Yue’s solutions (say, {\mathrm{\Phi }}_r for example) are presented in the forms of the inverse Hankel transform integrals with the Bessel functions J_K(\rho r) ( K=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3} ) and can be accurately calculated with any numerical integration techniques. For example, the following proceeding limit technique has been used [ 7, 8, 9, 14, 19].

(109a) {\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho }\approx {\displaystyle \underset{\mathrm{0}}{\overset{A_{\mathrm{1}}}{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho }+{\displaystyle \underset{A_{\mathrm{1}}}{\overset{A_{\mathrm{2}}}{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho }+\cdots +{\displaystyle \underset{A_m}{\overset{A_{m+\mathrm{1}}}{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho }

where is a sequence of numbers that approaches infinity and can be selected using the rule A_l={\lambda }_c{A}_{l-\mathrm{1}} , where l=\mathrm{2},\mathrm{3},\cdots ,m,m+\mathrm{1} and in particular, A_{\mathrm{1}}=\mathrm{2} and {\lambda }_c=\mathrm{1.5} . Each finite integral on the right-hand is a proper integral and can be calculated using the Simpson’s quadrature based adaptively iterative integration with a pre-specified allowable error {\delta }_c as follows.

(109b) {\displaystyle \frac{{\displaystyle \left|\underset{A_m}{\overset{A_{m+\mathrm{1}}}{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho \right|}}{\mathrm{1}+{\displaystyle \left|\underset{A_m}{\overset{A_{m+\mathrm{1}}}{\int }}{\mathrm{\Phi }}_r(\rho ,z,d)J_K(\rho r)\mathrm{d}\rho \right|}}}\le {\delta }_c

Elastodynamics

The governing equations

The governing equations of the mathematical theory of elastodynamics are exactly the same as those given in Eqs. (1) − (3) and Eqs. (5) − (9) in Refs. [ 1] for the elastostatics [ 20, 34− 36]. The equations of static equilibrium (4) in Ref. [ 1] are replaced by the equations of equilibrium of motion at any point within the solid. They take the following form by adding the inertial force vector on the left side of Eq. (4) in Ref. [ 1] along the x, y and z directions, respectively.

(110a) {\displaystyle \frac{\partial {\sigma }_{xx}}{\partial x}}+{\displaystyle \frac{\partial {\sigma }_{xy}}{\partial y}}+{\displaystyle \frac{\partial {\sigma }_{xz}}{\partial z}}+f_x=\gamma {\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_x}{\partial t^{\mathrm{2}}}}

(110b) {\displaystyle \frac{\partial {\sigma }_{yx}}{\partial x}}+{\displaystyle \frac{\partial {\sigma }_{yy}}{\partial y}}+{\displaystyle \frac{\partial {\sigma }_{yz}}{\partial z}}+f_y=\gamma {\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_y}{\partial t^{\mathrm{2}}}}

(110c) {\displaystyle \frac{\partial {\sigma }_{zx}}{\partial x}}+{\displaystyle \frac{\partial {\sigma }_{zy}}{\partial y}}+{\displaystyle \frac{\partial {\sigma }_{zz}}{\partial z}}+f_z=\gamma {\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_z}{\partial t^{\mathrm{2}}}}

where \gamma is the density of the solid material and t is the independent variable of time.

Two sets of governing ordinary differential equations

The matrix Fourier integral transform in Eqs. (14) – (17) in Ref. [ 1] are used on the two horizontal axes x and y. Then the classical Laplace transform is used on the time t. The two unknown vectors \mathbf{w}(\xi ,\eta ,z,t) and {\mathbf{Y}}_z(\xi ,\eta ,z,t) and the internal body force loading vector \mathbf{g}(\xi ,\eta ,z,t) have to be changed in the forms of Laplace transform as follows.

(111a) \widehat{\mathbf{w}}(\xi ,\eta ,z,s)={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}\mathbf{w}(\xi ,\eta ,z,t)}{e}^{-st}\mathrm{d}t

(111b) {\widehat{\mathbf{Y}}}_z(\xi ,\eta ,z,s)={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}{\mathbf{Y}}_z(\xi ,\eta ,z,t)}{e}^{-st}\mathrm{d}t

(111c) \widehat{\mathbf{g}}(\xi ,\eta ,z,s)={\displaystyle \underset{\mathrm{0}}{\overset{+\infty }{\int }}\mathbf{g}(\xi ,\eta ,z,t)}{e}^{-st}\mathrm{d}t-\gamma s{\displaystyle \frac{\partial }{\partial t}}\mathbf{w}(\xi ,\eta ,z,\mathrm{0})-\gamma \mathbf{w}(\xi ,\eta ,z,\mathrm{0})

where s is the independent variable corresponding to the time t in Laplace transform. \mathrm{Re}(s)>\mathrm{0} . \mathbf{w}(\xi ,\eta ,z,\mathrm{0}) and {\displaystyle \frac{\partial }{\partial t}}\mathbf{w}(\xi ,\eta ,z,\mathrm{0}) are the initial displacement and initial velocity vectors in the Fourier transform domain, respectively.

Accordingly, the governing equations can be reduced two sets of first order ordinary differential equations [ 20]. The first set is due to the anti-symmetry about the z-axis and has two linear ordinary differential equations with two field variables and variable coefficients with z. It can be expressed as follows.

(112a) {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\widehat{\mathbf{V}}(z)=\rho {\widehat{\mathbf{C}}}_v(z)\widehat{\mathbf{V}}(z)+{\widehat{\mathbf{G}}}_v(z)

where a\le z\le b , , and

(112b) \widehat{\mathbf{V}}(z)=\left(\begin{array}{l}{\widehat{w}}_{\mathrm{2}}\\ {\widehat{\tau }}_{\mathrm{2}}\end{array}\right),{\widehat{\mathbf{G}}}_v=\left(\begin{array}{l}\mathrm{0}\\ {\widehat{g}}_{\mathrm{2}}\end{array}\right),{\widehat{\mathbf{C}}}_v(z)=\left[\begin{array}{ll}\mathrm{0}& \mathrm{1}/c_{\mathrm{4}}\\ c_{\mathrm{5}}+\gamma s^{\mathrm{2}}/{\rho }^{\mathrm{2}}& \mathrm{0}\end{array}\right]

The second set is due to the axial symmetry about the z-axis and has four linear ordinary differential equations with four field variables and variable coefficients with z. It can be expressed as follows.

(113a) {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\widehat{\mathbf{U}}(z)=\rho {\widehat{\mathbf{C}}}_u(z)\widehat{\mathbf{U}}(z)+{\widehat{\mathbf{G}}}_u(z)

where a\le z\le b , , and

(113b) \widehat{\mathbf{U}}(z)=\left(\begin{array}{l}{\widehat{w}}_{\mathrm{1}}\\ {\widehat{w}}_{\mathrm{3}}\\ {\widehat{\tau }}_{\mathrm{3}}\\ {\widehat{\tau }}_{\mathrm{1}}\end{array}\right),{\widehat{\mathbf{G}}}_u=\left(\begin{array}{l}\mathrm{0}\\ \mathrm{0}\\ {\widehat{g}}_{\mathrm{3}}\\ {\widehat{g}}_{\mathrm{1}}\end{array}\right),{\widehat{\mathbf{C}}}_u(z)=\left[\begin{array}{cccc}\mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{1}/c_{\mathrm{4}}\\ c_{\mathrm{2}}/c_{\mathrm{3}}& \mathrm{0}& \mathrm{1}/c_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \gamma s^{\mathrm{2}}/{\rho }^{\mathrm{2}}& \mathrm{0}& \mathrm{1}\\ c_p+\gamma s^{\mathrm{2}}/{\rho }^{\mathrm{2}}& \mathrm{0}& -c_{\mathrm{2}}/c_{\mathrm{3}}& \mathrm{0}\end{array}\right]

The five elastic parameters in Eqs. (112b) and (113b) can be arbitrary functions of the depth z, i.e., c=_i{c}_i(z),i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5} . Most importantly, the matrix approach eliminates the two independent variables \xi and \eta in the six governing ordinary differential equations and preserves only the radial distance \rho of the material axial symmetry about the z-axis. The two coefficient matrices {\mathbf{C}}_v(z) and {\mathbf{C}}_u(z) contain only the five material parameters c(=c_i(z),i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5})_i and the combined variable \gamma s^{\mathrm{2}}/{\rho }^{\mathrm{2}} .

The general solution of \widehat{\mathbf{V}}(\mathit{z})

The basic solution for the first set with constant elastic parameters can be obtained as follows [ 20],

(114a) \widehat{\mathbf{V}}(z)=\widehat{\mathbf{A}}(z-z_{\mathrm{1}})\widehat{\mathbf{V}}(z_{\mathrm{1}})-{\displaystyle \underset{z_{\mathrm{1}}}{\overset{z}{\int }}\widehat{\mathbf{A}}(z-\varsigma ){\widehat{\mathbf{G}}}_v(\varsigma )\mathrm{d}\varsigma }

where z\ge z_{\mathrm{1}} or z\le z_{\mathrm{1}} . The first basic square matrix \widehat{\mathbf{A}}(z) is defined as follow,

(114b) \widehat{\mathbf{A}}(z)=\mathbf{B}\left({\widehat{\gamma }}_{\mathrm{0}}\right)e^{{\widehat{\gamma }}_{\mathrm{0}}\rho z}+\mathbf{B}\left(-{\widehat{\gamma }}_{\mathrm{0}}\right)e^{-{\widehat{\gamma }}_{\mathrm{0}}\rho z}

where the material characteristic root {\widehat{\gamma }}_{\mathrm{0}} and the material square matrix \mathbf{B}(\chi ) are defined as follows,

(114c) {\widehat{\gamma }}_{\mathrm{0}}=\sqrt{{\displaystyle \frac{c_{\mathrm{5}}}{c_{\mathrm{4}}}}+{\displaystyle \frac{S}{c_{\mathrm{4}}}}}>\mathrm{0}

(114d) \mathbf{B}(\chi )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{ll}\mathrm{1}& {\displaystyle \frac{\mathrm{1}}{c_{\mathrm{4}}\chi }}\\ c_{\mathrm{4}}\chi & \mathrm{1}\end{array}\right]

where S={\displaystyle \frac{\gamma s^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}} .

It is evident that the above basic solution matrix \widehat{\mathbf{A}}(z) automatically reduces to the basic solution matrix \mathbf{A}(z) in (20) in [ 1] once S=\mathrm{0} or s=\mathrm{0} . It also has the following properties.

(115a) \det \widehat{\mathbf{A}}(z)=\mathrm{1}

(115b) \widehat{\mathbf{A}}(\mathrm{0})={\mathbf{I}}_{\mathrm{2}}

(115c) \widehat{\mathbf{A}}(z)\widehat{\mathbf{A}}(z_{\mathrm{1}})=\widehat{\mathbf{A}}(z+z_{\mathrm{1}})

(115d) \widehat{\mathbf{A}}{(z)}^{-\mathrm{1}}=\widehat{\mathbf{A}}(-z)

The general solution of \widehat{\mathbf{U}}(z)

Similarly, the general matrix solution can be obtained as follows for the second set of four linear ordinary differential equations with constant coefficients [ 20],

(116a) \widehat{\mathbf{U}}(z)=\widehat{\mathbf{Q}}(z-z_{\mathrm{1}})\widehat{\mathbf{U}}(z_{\mathrm{1}})-{\displaystyle \underset{z_{\mathrm{1}}}{\overset{z}{\int }}\widehat{\mathbf{Q}}(z-\varsigma ){\widehat{\mathbf{G}}}_u(\varsigma )\mathrm{d}\varsigma }

where z\ge z_{\mathrm{1}} or z\le z_{\mathrm{1}} . The second basic square matrix \widehat{\mathbf{Q}}(z) is defined as follow.

(116b) \widehat{\mathbf{Q}}(z)=\mathbf{C}\left({\widehat{\gamma }}_{\mathrm{1}}\right)e^{{\widehat{\gamma }}_{\mathrm{1}}\rho z}+\mathbf{C}\left(-{\widehat{\gamma }}_{\mathrm{1}}\right)e^{-{\widehat{\gamma }}_{\mathrm{1}}\rho z}-\mathbf{C}\left({\widehat{\gamma }}_{\mathrm{2}}\right)e^{{\widehat{\gamma }}_{\mathrm{2}}\rho z}-\mathbf{C}\left(-{\widehat{\gamma }}_{\mathrm{2}}\right)e^{-{\widehat{\gamma }}_{\mathrm{2}}\rho z}

where {\widehat{\gamma }}_{\mathrm{1}} and {\widehat{\gamma }}_{\mathrm{2}} are the roots of the following material characteristic equation,

(116c) c_{\mathrm{3}}{c}_{\mathrm{4}}{\chi }^{\mathrm{4}}-[c_{\mathrm{1}}{c}_{\mathrm{3}}-c_{\mathrm{2}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{4}}+(c_{\mathrm{3}}+c_{\mathrm{4}})S]{\chi }^{\mathrm{2}}+(c_{\mathrm{1}}+S)(c_{\mathrm{4}}+S)=\mathrm{0}

They can be expressed as follows.

(116d) {\widehat{\gamma }}_{\mathrm{1}}=\sqrt{{\displaystyle \frac{c_{\mathrm{1}}{c}_{\mathrm{3}}-c_{\mathrm{2}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{4}}+(c_{\mathrm{3}}+c_{\mathrm{4}})S+\sqrt{{[c_{\mathrm{1}}{c}_{\mathrm{3}}-c_{\mathrm{2}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{4}}+(c_{\mathrm{3}}+c_{\mathrm{4}})S]}^{\mathrm{2}}+\mathrm{4}{c}_{\mathrm{3}}{c}_{\mathrm{4}}(c_{\mathrm{1}}+S)(c_{\mathrm{4}}+S)}}{\mathrm{2}{c}_{\mathrm{3}}{c}_{\mathrm{4}}}}}

(116e) {\widehat{\gamma }}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{c_{\mathrm{1}}{c}_{\mathrm{3}}-c_{\mathrm{2}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{4}}+(c_{\mathrm{3}}+c_{\mathrm{4}})S-\sqrt{{[c_{\mathrm{1}}{c}_{\mathrm{3}}-c_{\mathrm{2}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{4}}+(c_{\mathrm{3}}+c_{\mathrm{4}})S]}^{\mathrm{2}}+\mathrm{4}{c}_{\mathrm{3}}{c}_{\mathrm{4}}(c_{\mathrm{1}}+S)(c_{\mathrm{4}}+S)}}{\mathrm{2}{c}_{\mathrm{3}}{c}_{\mathrm{4}}}}}

The material square matrix \mathbf{C}(\chi ) is defined as follows.

(117a) \mathbf{C}(\chi )=\left(\begin{array}{llll}{C}_{\mathrm{11}}(\chi )& C_{\mathrm{12}}(\chi )& C_{\mathrm{13}}(\chi )& C_{\mathrm{14}}(\chi )\\ C_{\mathrm{21}}(\chi )& C_{\mathrm{22}}(\chi )& C_{\mathrm{23}}(\chi )& -C_{\mathrm{13}}(\chi )\\ C_{\mathrm{31}}(\chi )& C_{\mathrm{32}}(\chi )& C_{\mathrm{22}}(\chi )& -C_{\mathrm{12}}(\chi )\\ C_{\mathrm{41}}(\chi )& -C_{\mathrm{21}}(\chi )& -C_{\mathrm{31}}(\chi )& C_{\mathrm{11}}(\chi )\end{array}\right)

where

(117b) C_{\mathrm{11}}(\chi )=B{c}_{\mathrm{4}}[c_{\mathrm{3}}{\chi }^{\mathrm{2}}+(c_{\mathrm{2}}-S)]

(117c) C_{\mathrm{12}}(\chi )=-B[c_{\mathrm{4}}{c}_{\mathrm{3}}\chi +c_{\mathrm{2}}(c_{\mathrm{4}}+S)/\chi ]

(117d) C_{\mathrm{13}}(\chi )=-B(c_{\mathrm{2}}+c_{\mathrm{4}})

(117e) C_{\mathrm{14}}(\chi )=B[c_{\mathrm{3}}\chi -(c_{\mathrm{4}}+S)/\chi ]

(117f) C_{\mathrm{21}}(\chi )=B{c}_{\mathrm{4}}[c_{\mathrm{2}}\chi +(c_{\mathrm{1}}+S)/\chi ]

(117g) C_{\mathrm{22}}(\chi )=B[c_{\mathrm{4}}{c}_{\mathrm{3}}{\chi }^{\mathrm{2}}+c_{\mathrm{2}}(c_{\mathrm{2}}+c_{\mathrm{4}})-c_{\mathrm{3}}(c_{\mathrm{1}}+S)]

(117h) C_{\mathrm{23}}(\chi )=B[c_{\mathrm{4}}\chi -(c_{\mathrm{1}}+S)/\chi ]

(117i) C_{\mathrm{31}}(\chi )=B{c}_{\mathrm{44}}[(c_{\mathrm{3}}{c}_{\mathrm{1}}-c_{\mathrm{2}}^{\mathrm{2}})+(c_{\mathrm{3}}+c_{\mathrm{2}})S]

(117j) C_{\mathrm{32}}(\chi )=B\{c_{\mathrm{4}}{c}_{\mathrm{3}}^{\mathrm{2}}{\chi }^{\mathrm{3}}+c_{\mathrm{3}}[c_{\mathrm{2}}(c_{\mathrm{2}}+\mathrm{2}{c}_{\mathrm{4}})-c_{\mathrm{3}}(c_{\mathrm{1}}+S)]\chi +c_{\mathrm{2}}^{\mathrm{2}}(c_{\mathrm{4}}+S)/\chi \}

(117k) C_{\mathrm{41}}(\chi )=B{c}_{\mathrm{4}}^{\mathrm{2}}[c_{\mathrm{3}}{\chi }^{\mathrm{3}}+(\mathrm{2}{c}_{\mathrm{2}}-S)\chi +(c_{\mathrm{1}}+S)/\chi ]

where B={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{3}}{c}_{\mathrm{4}}\left({\widehat{\gamma }}_{\mathrm{1}}^{\mathrm{2}}-{\widehat{\gamma }}_{\mathrm{2}}^{\mathrm{2}}\right)}} .

It is also evident that the above basic solution matrix \widehat{\mathbf{Q}}(z) automatically reduces to the basic solution matrix \mathbf{Q}(z) in Eq. (22) in Ref. [1] once S=\mathrm{0} or s=\mathrm{0} . It also has the following properties.

(118a) \det \widehat{\mathbf{Q}}(z)=\mathrm{1}

(118b) \widehat{\mathbf{Q}}(\mathrm{0})=I_{\mathrm{4}}

(118c) \widehat{\mathbf{Q}}(z)\widehat{\mathbf{Q}}(z_{\mathrm{1}})=\mathbf{Q}(z+z_{\mathrm{1}})

(118d) \widehat{\mathbf{Q}}{(z)}^{-\mathrm{1}}=\widehat{\mathbf{Q}}(-z).

Formulation of solutions

The above general solutions can be used to derive and formulate solutions for specific initial-boundary value problems in elastodynamics in layered solids. The mathematical formulation procedures are similar to those given above for the boundary-value problems in elastostatics in the transform domain. The remaining issue is the inverse Hankel and Laplace transforms for the solution in spatial and time domain. Such inverse integral transforms can have singularities because of the dynamic behavior. Their accurate results in the spatial and time domain are difficult and need further studies with both analytical and numerical methods. Furthermore, it can be shown that the results given in Eqs. (112) – (117) can be reduced to those given in Ref. [ 35] for isotropic layered solid in frequency domain.

Poroelasticity

The theory of three dimensional linear poroelasticity for a saturated porous material (say, soil) was formulated by Biot in 1941 and 1956 [ 37, 38] to model the behavior of saturated soils within the working load range. The saturated soils are modeled as deformable, linear, porous, elastic materials saturated with compressible fluids. A set of partial differential equations was established to describe the coupled behavior of saturated soils and make the poroelasticity become a completely self-consistent initial-boundary value problem. Many people have derived solutions for poroelasticity.

In particular, the above matrix Fourier and Laplace transform approach has been used to derive analytical solutions for the initial-boundary value problems of the poroelasticity [ 21− 28]. The set of governing partial differential equations in the physical domain has reduced into two sets of first order ordinary differential equations in the transform domain. Their general solutions have been derived for a homogeneous poroelastic layer. They are further used for solving specific initial-boundary value problems in layered poroelastic solids. These solutions have been used to examine and predict the ground settlements and porewater pressure in saturated poroelastic layers due to flexible and rigid foundation loadings and the time-dependent behavior of rigid disc in a saturated poroelastic material. Details can be found in the publications [ 21− 28].

Thermoelasticity

Similarly, the above matrix Fourier and Laplace transform approach can be used to derive analytical solutions for the initial-boundary value problems of the mathematical theory of linear thermoelasticity. The Hooke’s law shall include the temperature expansion effect for the linear relationship of the three normal stresses with the three normal strains plus the thermal expansion terms. In addition, there are four time-dependent heat conduction equations. The following material parameters shall be added: the coefficient of linear heat expansion, the specific heat, the heat conduction coefficient, the thermoelastic coefficient. The set of governing partial differential equations in the physical domain can be reduced into two sets of first order ordinary differential equations in the transform domain. Their general solutions can be derived for a homogeneous thermoelastic layer. They can be further used for solving specific initial-boundary value problems in layered or functionally graded thermoelastic solids. Details of the mathematical formulations and equations can be found in Yue [ 29] and Ai et al. [ 39].

Concluding remarks

In the above, the author has verified Yue’s approach, Yue’s treatment, Yue’s method and Yue’s solution for complete set of exact and analytical solutions for boundary-value problems in n-layered elastic solids of either transverse isotropy or isotropy. Three levels of the mathematical verification have made and presented. They are the degeneration of Yue’s solution to the basic solutions in closed-form, the convergence, singularity and satisfaction of Yue’s solution, and the systematic and uniform applications of Yue’s approach, treatment and method to derive and formulate new solutions of other problems with wide interests in science and engineering.

In particular, for the three cases of loading distribution f(x,y) including the point loading, the concentrated ring loading and the concentrated rectangular loading, the singularities of the solutions have been isolated and integrated analytically in exact closed-forms in terms of elementary harmonic functions and special functions. Consequently, the solutions at any point in n-layered or graded solids can be calculated with any controlled accuracy in association with any classical numerical integration techniques. Furthermore, the author has shown that the closed-form singularity can be systematically and easily obtained if the following basic harmonic integral in Eq. (67a) or (67b) can be integrated analytically. They are re-expressed as follows.

(119a) Q\left(x,y,z,\tilde{f}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \frac{e^{-\rho z}}{{\rho }^{\mathrm{3}}}}\tilde{f}(\xi ,\eta )K\mathrm{d}\xi \mathrm{d}\eta }

or

(119b) Q(x,y,z,f)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }} }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}[z\ln (z+R)-R]f(s,t)\mathrm{d}s\mathrm{d}t}

where z>\mathrm{0} and R=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+z^{\mathrm{2}}} .

It is further noted that, the mathematical tools used by the author are classical tools of more than 150 years history. They include the Fourier integral transform, Hankel transforms, Laplace transforms, Fourier series, matrix operations, linear algebra equations, integrations, differentiations, partial differential equations, ordinary differential equations, improper integrals with depending parameters, harmonic functions and special functions. Using these classical mathematical tools, the author has systematically and uniformly presented the mathematical formulation and verification of the solutions in n-layered solid in matrix forms. Accordingly, from 1984 to present, the author has found and presented many new solutions with known mathematical properties and singularities in homogeneous or layered solid. Because the classical elasticity theory is the common foundation of many other continuum or field theories, the approach, the treatment, the method and the solutions have been applied to systematically and uniformly derive and formulate many new solutions in elastodynamics, poroelasticity and thermoelasticity and with boundary element methods.

Hence, the researchers at Research Centre Jülich and Massachusetts Institute of Technology have shortly named these mathematical formulation, verification and solutions as Yue’s approach, Yue’s treatment, Yue’s method and Yue’s solution.

References

[1] Yue Z Q. Yue’s solution of classical elasticity in n-layered solids: Part 1, mathematical formulation. Frontiers of Structural and Civil Engineering, 2015, 9(3): 215–249

[2] Yue Z Q. Elastic fields in two joined transversely isotropic solids due to concentrated forces. International Journal of Engineering Science, 1995, 33(3): 351–369

[3] Xiao H T, Yue Z Q. Elastic fields in two joined transversely isotropic media of infinite extent as a result of rectangular loading. International Journal for Numerical & Analytical Methods in Geomechanics, 2013, 37(3): 247–277

[4] Yue Z Q. Closed-form Green’s functions for transversely isotropic bi-solids with a slipping interface. Structural Engineering and Mechanics. An International Journal, 1996, 4(5): 469–484

[5] Yue Z Q, Yin J H. Closed-form fundamental solutions for transversely isotropic bi-materials with inextensible interface. Journal of Engineering Mechanics, ASCE, 1996, 122(11): 1052–1059

[6] Yue Z Q, Xiao H T, Tham L G, Lee C F, Yin J H. Stresses and displacements of a transversely isotropic elastic halfspace due to rectangular loadings. Engineering Analysis with Boundary Elements, 2005, 29(6): 647–671

[7] Yue Z Q. On generalized Kelvin solutions in multilayered elastic media. Journal of Elasticity, 1995, 40(1): 1–44

[8] Yue Z Q. On elastostatics of multilayered solids subjected to general surface traction. Quarterly Journal of Mechanics and Applied Mathematics, 1996, 49(3): 471–499

[9] Yue Z Q. 1996. Elastic field for an eccentrically loaded rigid plate on multilayered solids. International Journal of Solids and Structures, 33(27): 4019–4049

[10] Yue Z Q, Svec O C. Effect of tire-pavement contact pressure distribution on the response of asphalt concrete pavements. Canadian Journal of Civil Engineering, 1995, 22(5): 849–860

[11] Yue Z Q, Yin J H. Backward transfer-matrix method for elastic analysis of layered solid with imperfect bonding. Journal of Elasticity, 1998, 50(2):109–128

[12] Yue Z Q, Yin J H. Layered elastic model for analysis of cone penetration testing. International Journal for Numerical & Analytical Methods in Geomechanics, 1999, 23(8): 829–843

[13] Yue Z Q, Yin J H, Zhang S Y. Computation of point load solutions for geo-materials exhibiting elastic non-homogeneity with depth. Computers and Geotechnics, 1999, 25(2): 75–105

[14] Yue Z Q, Xiao H T, 2002. Generalized Kelvin solution based boundary element method for crack problems in multilayered solids. Engineering Analysis with Boundary Elements, 26(8): 691–705

[15] Yue Z Q, Xiao H T, Tham L G. Boundary element analysis of crack problems in functionally graded materials. International Journal of Solids and Structures, 2003, 40(13–14): 3273–3291

[16] Yue Z Q, Xiao H T, Tham L G. Elliptical crack normal to functionally graded interface of bonded solids. Theoretical and Applied Fracture Mechanics. 2004, 42: 227–248

[17] Yue Z Q, Xiao H T, Tham L G, Lee C F, Pan E. Boundary element analysis of three-dimensional crack problems in two joined transversely isotropic solids. Computational Mechanics, 2005, 36(6): 459–474

[18] Yue Z Q, Xiao H T, Pan E. Stress intensity factors of square crack inclined to interface of transversely isotropic bi-material. Engineering Analysis with Boundary Elements, 007, 31(1): 50–65

[19] Xiao H T, Yue Z Q. Fracture Mechanics in Layered and Graded Solids: Analysis Using Boundary Element Methods. Berlin: De Gruyter & Higher Education Press, 2014, 305

[20] Yue Z Q. Solutions of transversely isotropic elastodynamics in a vertically heterogeneous halfspace. In: Proceedings of the 2nd Chinese Congress of Earthquake Engineering. Wuhan, China, 1987 (in Chinese)

[21] Yue Z Q, Selvadurai A P S. The role of Poisson’s ratios on the consolidation response of soils. In: Proceedings of the 45th Conference of Canadian Geotechnical Engineering. Toronto, Canada, <month>October</month> 1992, 11–1 to 11–11

[22] Yue Z Q, Selvadurai A P S. Eccentric settlement of a rigid foundation on a consolidating thin layer. Vertical and Horizontal Deformation of Foundation and Embankments, ASCE Geotechnical Special Publication, 1994, 40 (2): 612–627

[23] Yue Z Q, Selvadurai A P S. On the asymmetric indentation of a consolidating poroelastic halfspace. International Journal of Applied Mathematical Modeling, 1994, 18: 170–185

[24] Yue Z Q, Selvadurai A P S, Law K T. Excess pore pressure in a poroelastic seabed saturated with a compressible fluid. Canadian Geotechnical Journal, 1994, 31: 989–1003

[25] Selvadurai A P S, Yue Z Q. On the indentation of a poroelastic layer. International Journal of Numerical & Analytical Methods for Geomechanics, 1994, 18: 161–175

[26] Yue Z Q, Selvadurai A P S. On the mechanics of a rigid disc inclusion in a fluid saturated poroelastic medium. International Journal of Engineering Science, 1995, 33(11): 1633–1662

[27] Yue Z Q, Selvadurai A P S. Contact problem for saturated poroelastic solid. Journal of Engineering Mechanics, ASCE, 1995, 121(4): 502–512

[28] Yue Z Q, Selvadurai A P S. Axisymmetric indentation of a multilayered poroelastic solid. Mechanics of Poroelastic Media, ASME Special Publication, Kluwer Academic Publishers. The Netherlands: Dordrecht, 1996, 235–241

[29] Yue Z Q. Solutions for the thermoelastic problems in vertically inhomogeneous media. Acta Mechnica Sinica, 1988, 4(2): 182–189

[30] Sneddon I N. Mixed Boundary Value Problems in Potential Theory. Amsterdam: North Holland Publ. Co., 1966

[31] Sneddon I N, Lowengrub M. Crack Problems in the Classical Theory of Elasticity. New York: Wiley, 1969, 221

[32] Thompson, W (Lord Kelvin). Note on the integration of the equations of equilibrium of an elastic solid. Cambridge and Dublin Mathematical Journal, 1848, 1: 97–99

[33] Mindlin R D. Force at a point in the interior of a semi-infinite solid. Physics, 1936, 7: 195–202

[34] Kupradze V D. Three-dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity. Amsterdam: North Holland Publ. Co., 1979, 929

[35] Aki K, Richards P. G. Quantitative Seismology. 2nd ed. Sansalito, CA, University Science Books, 2002, 700

[36] Ai Z Y, Li Z X. Time-harmonic response of transversely isotropic multilayered half-space in a cylindrical coordinate system. Soil Dynamics and Earthquake Engineering, 2014, 66: 69–77

[37] Biot M A. General theory of three-dimensional consolidation. Journal of Applied Physics, 1941, 12(2): 155–164

[38] Biot M A. General solution of the equations of elasticity and consolidation for a porous material, Journal of Applied Mechanics, 1956, 23: 91–95

[39] Ai Z Y, Wang L J, Zeng K. Analytical layer-element method for 3D thermoelastic problem of layered medium around a heat source. Meccanica 2015, 50(1): 49–59

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg