# Elastic Fields of a Nonhomogeneous Half-Space Subject to an Inclined Circular Load

Authors: Yunyue Xie; Hongtian Xiao; Zhongqi Quentin Yue
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

The paper examines the elastic fields of displacements and stresses for a nonhomogeneous elastic half-space where the elastic parameters have a linear variation over a finite depth beyond which it is constant. The circular loading area is subjected to a uniform inclined load. The numerical method is developed by applying the fundamental solution of a layered elastic solid and integrating numerically it over the loading area. As a result, only the loading area needs to be discretized in using the proposed numerical method. Numerical examples of calculation of displacements are conducted, and excellent agreement with the existing closed-form solutions is obtained. The results obtained are used to understand the elastic fields induced by different types of loads in a nonhomogeneous medium.


---

## Body

## 1. Introduction

The nonhomogeneity of materials is a common phenomenon that can be observed in many natural and engineering materials. Naturally, soil and in situ rock media exhibit strong spatial variations in their properties because of their natural formation process. The nonhomogeneity problem in geomechanics has applications to many problems of technological importance. There have been numerous efforts attempting to understand the influence of nonhomogeneity on the response of an elastic half-space induced by various loads [1–3]. Due to the mathematical difficulties associated with the geomechanics of nonhomogeneity, the majority of the relevant investigations available in the literature make some assumptions on the elastic properties of materials and the distribution forms of loads.Considering the influence of gravity on the natural formation process of soils and rocks, the elastic properties of geomaterials are assumed to vary in the depth direction and keep constant along the other two directions perpendicular to the depth direction. It is generally assumed that the shear modulus varies continuously with depth in power law, linear, hyperbolic, and exponential functions [4–8]. For easy analysis, it is further assumed that Poisson’s ratio keeps constant with depth. Actually, depth-dependent nonhomogeneity is distributed in a much complex form and an effective analyzing approach can model arbitrary variations with depth.In this paper, the elastic fields induced by an inclined circular load in a nonhomogeneous elastic half-space, shown in Figure1, are analyzed. The problem of interior loading in a nonhomogeneous half-space has applications in geomechanics. Selvadurai and Katebi [9] analyzed the problem of the application of an axisymmetric circular load at the interior of a nonhomogeneous isotropic half-space. Here, we extend the study to include the influence of a nonaxisymmetric circular load, that is, an inclined load at the interior of a nonhomogeneous isotropic half-space. In particular, the variations of elastic modulus and Poisson’s ratio with depth are considered. The numerical method used for analysis is developed by applying the fundamental solution of layered elastic solids [10–17] and integrating it over the loading area. The adaptive numerical quadrature and parallel computation techniques are used to evaluate the integrals over the elements on the discretized area. Numerical results are presented in order to show the influence of nonhomogeneity in an elastic half-space on the elastic fields induced by different types of loads.Figure 1
Inclined load at the interior of a nonhomogeneous elastic half-space.


## 2. Numerical Method of Analyzing a Nonhomogeneous Half-Space

### 2.1. Basic Equations

The fundamental solution of a multilayered solid developed by Yue [10] is the analytical solution for the elastostatic field in a layered solid of infinite extent due to the action of concentrated point loads. The dissimilar homogeneous layers adhere to an elastic solid of upper semi-infinite extent and another elastic solid of lower semi-infinite extent. The interface between any two connected dissimilar layers is fully bonded and the layer number is an arbitrary nonnegative integer. Since 2000, Xiao and Yue [18] have incorporated this fundamental solution into the BEMs for the analysis of the fracture mechanics in layered solids and found the solutions for many specific problems of interests in science and technology. The fundamental solution of infinite multilayered media is also suitable for the layered media of semi-infinite extent. In this case, the elastic modulus of the upper semi-infinite medium has an extremely small value, such as E0=1.0×10−15MPa, and Poisson’s ratio of the medium ν0=0.3.As shown in Figure1, a nonhomogeneous half-space is subjected to internal loads in the x, y, and/or z directions. Using the fundamental solution [10], the displacements and stresses at any point of the nonhomogeneous medium are described as(1)uiQ=∫Suik∗Q,PtkPdSP,(2)σijQ=∫Sσijk∗Q,PtkPdSP,where P and Q are, respectively, the source and field points of the fundamental solution; uik∗Q,P and σijk∗Q,P are the kernel functions of the displacements and stresses of a layered medium presented by Yue [10]; tkP is the traction of the source point P; and the integral domain S is the loading area.


### 2.2. Numerical Scheme

#### 2.2.1. The Discretization of the Loading Area

The loading areaS is discretized into ne quadrilateral elements S=∑e=1neSe. For the convenience of calculation, isoparametric element is used. All discrete elements in the loading domain correspond to a standard element. On the boundary, the variable node element is adopted for the convenience of discretization. Here, variable node isoparametric elements with 4–8 nodes are selected, and the form is shown in Figure 2. Only four corner elements are linear elements, and eight nodes are quadratic elements.Figure 2
The isoparametric element. (a) Linear element with 4 nodes. (b) The quadrilateral element with 8 nodes.
(a)(b)In global and local coordinate systems, there are the following transform relationships of coordinate and traction values within an element:(3)xi=∑α=14to8Nαξ,ηxiα,ti=∑α=14to8Nαξ,ηtiα,where Nα is the shape function of the element at node α and xiα and tiα are, respectively, the coordinate and traction values of the element at node α.Considering equations (1) and (2) with (3), we have the following discretized forms:(4)uiQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11uik∗Q,PNαJdξdη,(5)σijQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11σijk∗Q,PNαJdξdη,where J is the Jacobian determinant.The interpolation function of the quadrilateral element is(6)Nαξ,η=121−ξ21+ηαη,α=5,7,Nαξ,η=121−η21+ξαξ,α=6,8,N1ξ,η=N1∗−12N8+N5,N2ξ,η=N2∗−12N5+N6,N3ξ,η=N3∗−12N6+N7,N4ξ,η=N4∗−12N7+N8,where Nα∗ξ,η=1/41+ξαξ1+ηαη,α=1,2,3,4, and ξa,ηa is the local coordinate component at node α. For the variable node element, the nodes in the four edges can be available or not. In this case, the interpolation function corresponding to the missing nodes in the edge can be taken as 0.


#### 2.2.2. Numerical Integration

The integrals of equations (4) and (5) may be calculated by using the regular Gaussian quadrature. If the field point Q is located at the integral element and Q = Pα, the singular integrals of equations (4) and (5) appear and can be computed by, respectively, applying a linear coordinate transformation and an indirect method. When the field point Q is not located on the loading area and the distance r between the points P and Q has a small value, the kernel functions uik∗ and σijk∗ in equations (4) and (5) have sharp variations. An adaptive integration approach for this so-called nearly singular integral is developed through dividing the integrating element into subelements according to the ratio of the length of an element to the distance from a field point to the element.In the following, the linear coordinate transformation for the weakly singular integral in equation (4) is presented.The coordinates of any point in the triangle can be expressed as(7)ξ=∑α=14Lαξ1,η1ξα,η=∑α=14Lαξ1,η1ηα.The triangle area is transformed into a square area−1≤ξ1≤1 and −1≤η1≤1 by the coordinate change of the linear element shown in Figure 2(a). Since the transformation is linear, the relationship between ξ,η and ξ1,η1 is established by linear interpolation function:(8)L1ξ1,η1=141−ξ11−η1,L2ξ1,η1=141+ξ11−η1,L3ξ1,η1=141+ξ11+η1,L4ξ1,η1=141−ξ11+η1,Taking the singular point appearing in corner node 1 and middle node 5 in the edge as an example, the transformation relationship between coordinate systems is established, as shown in Figure3.Figure 3
Singularity at node 1or node 5.
(a)(b)When the singular point is at node 1, the element is divided into two triangles: to triangle 1-3-4,(9a)ξ=L1+L2ξ1+L3ξ3+L4ξ4=121+ξ11+η1−1,(9b)η=L1+L2η1+L3η3+L4η4=η1,(10)J1=∂ξ,η∂ξ1,η1=∂ξ∂ξ1∂η∂ξ1∂ξ∂η1∂η∂η1=121+η10121+ξ11=121+η1=οr,and to triangle 1-2-3,(11a)ξ=L1+L2ξ1+L3ξ2+L4ξ3=η1,(11b)η=L1+L2η1+L3η2L4η3=121−ξ11+η1−1,(12)J1=121+η1=οr.When the singular point is at node 5, the element is divided into three triangles: to triangle 5-4-1,(13a)ξ=L1+L2ξ5+L3ξ4+L4ξ1=−121+η1,(13b)η=L1+L2η5+L3η4+L4η1=121+ξ11+η1−1,(14)J1=141+η1=οr,to triangle 5-2-3,(15a)ξ=L1+L2ξ5+L3ξ2+L4ξ3=121+η1,(15b)η=L1+L2η5+L3η2+L4η3=121−ξ11+η1−1,(16)J1=141+η1=οr,and to triangle 5-3-4,(17a)ξ=L1+L2ξ5+L3ξ3+L4ξ4=12ξ11+η1,(17b)η=L1+L2η5+L3η3+L4η4=η1,(18)J1=121+η2=οr.In the same way, it is not difficult to write out the coordinate transformation relationship when the singular point appears in other nodes.


### 2.3. Parallel Computing and Discretization of the Nonhomogeneous Layer

Since the integrals in equations (4) and (5) for each element are independent of the other, parallel computing can be implemented by the proposed numerical method. Here, a straight approach is to use OpenMP directives to parallelize the internal loop, which controls element iterations. For details of OpenMP method and implementation, please refer to [19].By using the proposed numerical method for the analysis of a nonhomogeneous half-space exhibiting arbitrary variations in depth, the half-space is discretized into a large number of homogeneous layers. The nonhomogeneous layer in a half-space is closely approximated byn bonded layers of elastic homogeneous media. Each layer has the elastic modulus and Poisson’s ratio equal to the values at the lower depth of the layer. A homogeneous material bonded with the nonhomogeneous layer is considered as a semi-infinite domain.


## 2.1. Basic Equations

The fundamental solution of a multilayered solid developed by Yue [10] is the analytical solution for the elastostatic field in a layered solid of infinite extent due to the action of concentrated point loads. The dissimilar homogeneous layers adhere to an elastic solid of upper semi-infinite extent and another elastic solid of lower semi-infinite extent. The interface between any two connected dissimilar layers is fully bonded and the layer number is an arbitrary nonnegative integer. Since 2000, Xiao and Yue [18] have incorporated this fundamental solution into the BEMs for the analysis of the fracture mechanics in layered solids and found the solutions for many specific problems of interests in science and technology. The fundamental solution of infinite multilayered media is also suitable for the layered media of semi-infinite extent. In this case, the elastic modulus of the upper semi-infinite medium has an extremely small value, such as E0=1.0×10−15MPa, and Poisson’s ratio of the medium ν0=0.3.As shown in Figure1, a nonhomogeneous half-space is subjected to internal loads in the x, y, and/or z directions. Using the fundamental solution [10], the displacements and stresses at any point of the nonhomogeneous medium are described as(1)uiQ=∫Suik∗Q,PtkPdSP,(2)σijQ=∫Sσijk∗Q,PtkPdSP,where P and Q are, respectively, the source and field points of the fundamental solution; uik∗Q,P and σijk∗Q,P are the kernel functions of the displacements and stresses of a layered medium presented by Yue [10]; tkP is the traction of the source point P; and the integral domain S is the loading area.


## 2.2. Numerical Scheme

### 2.2.1. The Discretization of the Loading Area

The loading areaS is discretized into ne quadrilateral elements S=∑e=1neSe. For the convenience of calculation, isoparametric element is used. All discrete elements in the loading domain correspond to a standard element. On the boundary, the variable node element is adopted for the convenience of discretization. Here, variable node isoparametric elements with 4–8 nodes are selected, and the form is shown in Figure 2. Only four corner elements are linear elements, and eight nodes are quadratic elements.Figure 2
The isoparametric element. (a) Linear element with 4 nodes. (b) The quadrilateral element with 8 nodes.
(a)(b)In global and local coordinate systems, there are the following transform relationships of coordinate and traction values within an element:(3)xi=∑α=14to8Nαξ,ηxiα,ti=∑α=14to8Nαξ,ηtiα,where Nα is the shape function of the element at node α and xiα and tiα are, respectively, the coordinate and traction values of the element at node α.Considering equations (1) and (2) with (3), we have the following discretized forms:(4)uiQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11uik∗Q,PNαJdξdη,(5)σijQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11σijk∗Q,PNαJdξdη,where J is the Jacobian determinant.The interpolation function of the quadrilateral element is(6)Nαξ,η=121−ξ21+ηαη,α=5,7,Nαξ,η=121−η21+ξαξ,α=6,8,N1ξ,η=N1∗−12N8+N5,N2ξ,η=N2∗−12N5+N6,N3ξ,η=N3∗−12N6+N7,N4ξ,η=N4∗−12N7+N8,where Nα∗ξ,η=1/41+ξαξ1+ηαη,α=1,2,3,4, and ξa,ηa is the local coordinate component at node α. For the variable node element, the nodes in the four edges can be available or not. In this case, the interpolation function corresponding to the missing nodes in the edge can be taken as 0.


### 2.2.2. Numerical Integration

The integrals of equations (4) and (5) may be calculated by using the regular Gaussian quadrature. If the field point Q is located at the integral element and Q = Pα, the singular integrals of equations (4) and (5) appear and can be computed by, respectively, applying a linear coordinate transformation and an indirect method. When the field point Q is not located on the loading area and the distance r between the points P and Q has a small value, the kernel functions uik∗ and σijk∗ in equations (4) and (5) have sharp variations. An adaptive integration approach for this so-called nearly singular integral is developed through dividing the integrating element into subelements according to the ratio of the length of an element to the distance from a field point to the element.In the following, the linear coordinate transformation for the weakly singular integral in equation (4) is presented.The coordinates of any point in the triangle can be expressed as(7)ξ=∑α=14Lαξ1,η1ξα,η=∑α=14Lαξ1,η1ηα.The triangle area is transformed into a square area−1≤ξ1≤1 and −1≤η1≤1 by the coordinate change of the linear element shown in Figure 2(a). Since the transformation is linear, the relationship between ξ,η and ξ1,η1 is established by linear interpolation function:(8)L1ξ1,η1=141−ξ11−η1,L2ξ1,η1=141+ξ11−η1,L3ξ1,η1=141+ξ11+η1,L4ξ1,η1=141−ξ11+η1,Taking the singular point appearing in corner node 1 and middle node 5 in the edge as an example, the transformation relationship between coordinate systems is established, as shown in Figure3.Figure 3
Singularity at node 1or node 5.
(a)(b)When the singular point is at node 1, the element is divided into two triangles: to triangle 1-3-4,(9a)ξ=L1+L2ξ1+L3ξ3+L4ξ4=121+ξ11+η1−1,(9b)η=L1+L2η1+L3η3+L4η4=η1,(10)J1=∂ξ,η∂ξ1,η1=∂ξ∂ξ1∂η∂ξ1∂ξ∂η1∂η∂η1=121+η10121+ξ11=121+η1=οr,and to triangle 1-2-3,(11a)ξ=L1+L2ξ1+L3ξ2+L4ξ3=η1,(11b)η=L1+L2η1+L3η2L4η3=121−ξ11+η1−1,(12)J1=121+η1=οr.When the singular point is at node 5, the element is divided into three triangles: to triangle 5-4-1,(13a)ξ=L1+L2ξ5+L3ξ4+L4ξ1=−121+η1,(13b)η=L1+L2η5+L3η4+L4η1=121+ξ11+η1−1,(14)J1=141+η1=οr,to triangle 5-2-3,(15a)ξ=L1+L2ξ5+L3ξ2+L4ξ3=121+η1,(15b)η=L1+L2η5+L3η2+L4η3=121−ξ11+η1−1,(16)J1=141+η1=οr,and to triangle 5-3-4,(17a)ξ=L1+L2ξ5+L3ξ3+L4ξ4=12ξ11+η1,(17b)η=L1+L2η5+L3η3+L4η4=η1,(18)J1=121+η2=οr.In the same way, it is not difficult to write out the coordinate transformation relationship when the singular point appears in other nodes.


## 2.2.1. The Discretization of the Loading Area

The loading areaS is discretized into ne quadrilateral elements S=∑e=1neSe. For the convenience of calculation, isoparametric element is used. All discrete elements in the loading domain correspond to a standard element. On the boundary, the variable node element is adopted for the convenience of discretization. Here, variable node isoparametric elements with 4–8 nodes are selected, and the form is shown in Figure 2. Only four corner elements are linear elements, and eight nodes are quadratic elements.Figure 2
The isoparametric element. (a) Linear element with 4 nodes. (b) The quadrilateral element with 8 nodes.
(a)(b)In global and local coordinate systems, there are the following transform relationships of coordinate and traction values within an element:(3)xi=∑α=14to8Nαξ,ηxiα,ti=∑α=14to8Nαξ,ηtiα,where Nα is the shape function of the element at node α and xiα and tiα are, respectively, the coordinate and traction values of the element at node α.Considering equations (1) and (2) with (3), we have the following discretized forms:(4)uiQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11uik∗Q,PNαJdξdη,(5)σijQ=∑e=1ne∑α=14to8tkαPα∫−11∫−11σijk∗Q,PNαJdξdη,where J is the Jacobian determinant.The interpolation function of the quadrilateral element is(6)Nαξ,η=121−ξ21+ηαη,α=5,7,Nαξ,η=121−η21+ξαξ,α=6,8,N1ξ,η=N1∗−12N8+N5,N2ξ,η=N2∗−12N5+N6,N3ξ,η=N3∗−12N6+N7,N4ξ,η=N4∗−12N7+N8,where Nα∗ξ,η=1/41+ξαξ1+ηαη,α=1,2,3,4, and ξa,ηa is the local coordinate component at node α. For the variable node element, the nodes in the four edges can be available or not. In this case, the interpolation function corresponding to the missing nodes in the edge can be taken as 0.


## 2.2.2. Numerical Integration

The integrals of equations (4) and (5) may be calculated by using the regular Gaussian quadrature. If the field point Q is located at the integral element and Q = Pα, the singular integrals of equations (4) and (5) appear and can be computed by, respectively, applying a linear coordinate transformation and an indirect method. When the field point Q is not located on the loading area and the distance r between the points P and Q has a small value, the kernel functions uik∗ and σijk∗ in equations (4) and (5) have sharp variations. An adaptive integration approach for this so-called nearly singular integral is developed through dividing the integrating element into subelements according to the ratio of the length of an element to the distance from a field point to the element.In the following, the linear coordinate transformation for the weakly singular integral in equation (4) is presented.The coordinates of any point in the triangle can be expressed as(7)ξ=∑α=14Lαξ1,η1ξα,η=∑α=14Lαξ1,η1ηα.The triangle area is transformed into a square area−1≤ξ1≤1 and −1≤η1≤1 by the coordinate change of the linear element shown in Figure 2(a). Since the transformation is linear, the relationship between ξ,η and ξ1,η1 is established by linear interpolation function:(8)L1ξ1,η1=141−ξ11−η1,L2ξ1,η1=141+ξ11−η1,L3ξ1,η1=141+ξ11+η1,L4ξ1,η1=141−ξ11+η1,Taking the singular point appearing in corner node 1 and middle node 5 in the edge as an example, the transformation relationship between coordinate systems is established, as shown in Figure3.Figure 3
Singularity at node 1or node 5.
(a)(b)When the singular point is at node 1, the element is divided into two triangles: to triangle 1-3-4,(9a)ξ=L1+L2ξ1+L3ξ3+L4ξ4=121+ξ11+η1−1,(9b)η=L1+L2η1+L3η3+L4η4=η1,(10)J1=∂ξ,η∂ξ1,η1=∂ξ∂ξ1∂η∂ξ1∂ξ∂η1∂η∂η1=121+η10121+ξ11=121+η1=οr,and to triangle 1-2-3,(11a)ξ=L1+L2ξ1+L3ξ2+L4ξ3=η1,(11b)η=L1+L2η1+L3η2L4η3=121−ξ11+η1−1,(12)J1=121+η1=οr.When the singular point is at node 5, the element is divided into three triangles: to triangle 5-4-1,(13a)ξ=L1+L2ξ5+L3ξ4+L4ξ1=−121+η1,(13b)η=L1+L2η5+L3η4+L4η1=121+ξ11+η1−1,(14)J1=141+η1=οr,to triangle 5-2-3,(15a)ξ=L1+L2ξ5+L3ξ2+L4ξ3=121+η1,(15b)η=L1+L2η5+L3η2+L4η3=121−ξ11+η1−1,(16)J1=141+η1=οr,and to triangle 5-3-4,(17a)ξ=L1+L2ξ5+L3ξ3+L4ξ4=12ξ11+η1,(17b)η=L1+L2η5+L3η3+L4η4=η1,(18)J1=121+η2=οr.In the same way, it is not difficult to write out the coordinate transformation relationship when the singular point appears in other nodes.


## 2.3. Parallel Computing and Discretization of the Nonhomogeneous Layer

Since the integrals in equations (4) and (5) for each element are independent of the other, parallel computing can be implemented by the proposed numerical method. Here, a straight approach is to use OpenMP directives to parallelize the internal loop, which controls element iterations. For details of OpenMP method and implementation, please refer to [19].By using the proposed numerical method for the analysis of a nonhomogeneous half-space exhibiting arbitrary variations in depth, the half-space is discretized into a large number of homogeneous layers. The nonhomogeneous layer in a half-space is closely approximated byn bonded layers of elastic homogeneous media. Each layer has the elastic modulus and Poisson’s ratio equal to the values at the lower depth of the layer. A homogeneous material bonded with the nonhomogeneous layer is considered as a semi-infinite domain.


## 3. Numerical Verification

Selvadurai and Katebi [9] analyzed the undrained behavior of a nonhomogeneous elastic medium induced by a uniform vertical load. It was assumed that the variation of shear modulus with the depth was described by the expressions Gz=G0eλzz≤d and Gz=G0eλdz≥d and Poisson’s ratio ν = 0.5. We reexamine the problems by verifying the accuracy of the proposed method and choosing the discretized mesh on the loading area.The circular loading area (radiusa) is discretized into the four meshes with eight-noded elements, which have, respectively, 390 elements and 1235 nodes, 1589 elements and 4896 nodes, 2439 elements and 7478 nodes, and 6172 elements and 18773 nodes. As shown in Figure 4, the nonhomogeneous half-space with λ = 0.1 is discretized into 100 homogeneous elastic layers and the displacements and stresses induced by the vertical uniform load fz at h = a are presented for the approximation by depth-dependent variations of the elastic parameters.Figure 4
Approximation of the continuous depth variation of the shear modulus by a layered system of 100 piecewise homogeneous layers forλ=0.1 and h/a = 1.Figure5 shows the vertical displacement of the points along the vertical axis passing through the center of the loading area for different meshes. It can be found that the vertical displacements obtained by using the four meshes are much close to each other. With comparison to the analytical results by Selvadurai and Katebi [9], the displacements above the loading plane z > a have small errors and the ones below the loading plane have relatively large errors. For 0≤z<5a, the largest absolute error of uz0,0,z/fz is 0.001 and appears at z = a. Figure 6 illustrates the vertical stress σzz at the points close to the loading plane for different meshes. It is known that the loading fz causes a jump discontinuity of the stress σzz across the loading plane and the jump at the center of the loading area is equal tofz. For two points with a distance of 0.003a from the loading surface for Meshes 3 and 4, the jumps at the center of the loading surface are 0.9934fz and 0.99743fz, respectively, and have the absolute errors of 0.0066fz and 0.00257fz, respectively. These calculation errors are induced by nearly singular integrals and it is difficult to completely eliminate the errors by the subelement method above. Thus, we choose Mesh 4 to perform the following analysis and present the results of points with a distance greater than or equal to 0.003a from the loading surface.Figure 5
Variation of the displacementuz(0, 0, z) for different meshes.Figure 6
Variations ofσzz at the points close to the loading plane for different meshes.


## 4. Numerical Evaluations and Examples

### 4.1. Elastic Parameters of a Nonhomogeneous Medium

In order to consider the variation of Poisson’s ratio with depth, a simple linear fit has been completed for the data by Pan [3], who investigated the depth variation of the geotechnical properties of sand soil, clay loam, clay, and soft soil. The variations of the elastic modulus and Poisson’s ratio for sand soil are presented as follows:(19)Ez=E01+m1z,νz=ν01+m2z,where E0=53.09MPa,m1=0.5065, ν0=0.3469, m2=0.0123, and the SI unit of z is meters (m). The thickness of the sand soil is d = 12 m. For z > d, the elastic properties keep constant; that is, Ez≥d=E01+m1d and νz≥d=ν01+m2d. It is considered that the depth of the nonhomogeneous medium is equal to the radius of the circular loading area; that is, d = a.


### 4.2. The Elastic Fields Induced by Inclined Loads Located at a Given Depth

It is assumed that the circular loading area (radiusa) is located at the depth h = a and is subjected to an inclined load p0, which has an angle θ with the loading plane. Here, we present only the results along the vertical axis passing the center of the loading area. For comparison, two cases of the depth-dependent and depth-independent variations of Poisson’s ratio are analyzed.Figure7 illustrates the displacements ux and uz at the point x,y,z=0,0,z induced by the inclined load p0 for different loading angles. Because of symmetry, the displacement uy at the point x,y,z=0,0,z is equal to zero. It can be found that the variation of Poisson’s ratio with depth exerts an influence on the distribution of the displacements. The influence of Poisson’s ratio is more obvious on the vertical displacement uz. Furthermore, the horizontal displacement ux becomes small with the loading angle θ increasing and is equal to zero for θ = 90°. The vertical displacements uz increase with the loading angle θ increasing and arrive at the maximum values at θ = 90°.Figure 7
Variations of the displacementsux and uz for d = a and h = a for different θ. (a) uxE0/p0a; (b) uzE0/p0a.
(a)(b)Figure8 illustrates the variations of the stresses σxx, σxz, σyy, and σzz at the point x,y,z=0,0,z induced by the inclined load p0 from different loading angles. Because of symmetry, the stresses σxy and σyz at the point x,y,z=0,0,z are equal to zero. From Figure 6, we have the following observations:(1)
All the stresses tend to zero with the distance from the loading plane increasing and have jumps across the loading plane.(2)
The variation of Poisson’s ratio with depth exerts an influence on the distribution of the stresses. The influence of Poisson’s ratio is more obvious on the normal stressesσxx, σyy, and σzz.(3)
With the loading angle increasing, the absolute values ofσxx decrease and the ones of σyy and σzz increase below the loading plane. Above the loading plane, the variations of the three normal stresses with the loading angle are different in different depths. The absolute values of σxz decrease with the loading angle increasing and σxz=0 for θ = 90°.Figure 8
Variations of the four stress components ford = a and h = a for different θ. (a) σxx/p0; (b) σyy/p0; (c) σzz/p0; (d) σxz/p0.
(a)(b)(c)(d)


### 4.3. The Elastic Fields Induced by Inclined Loads Located at Different Depths

It is assumed that the circular loading area (radiusa) is located at different depths (h = 0, 2a, 3a, and 5a) and is subjected to an inclined loadp0, which has an angle θ = 45° with the loading plane. Here, we use the depth-independent model of elastic modulus and Poisson’s ratio and present only the results along the vertical axis passing through the center of the loading area.Figure9 illustrates the displacements ux and uzat the point x,y,z=0,0,z induced by the inclined load p0. It can be found that, for each loading depth, the maximum values of ux and uz appear at the loading plane. And the maximum values of ux and uz decrease with the loading depth increasing. Furthermore, the displacements ux and uz for h = 0 are much larger than the ones for other loading depths. At z = 5a, the displacements of ux and uz for h = 5a are larger than the ones for other loading depths.Figure 9
Variations of the displacementsux and uz for d = a and θ = 45° for different h. (a) uxE0/pa; (b) uzE0/pa.
(a)(b)Figure10 illustrates the variations of the stresses σxx, σxz, σyy, and σzz at the point x,y,z=0,0,z with the loading depth. It can be found that all the jumps of the stresses appear at the loading planes. Apart from the loading plane, all the stresses approach zero. Furthermore, except for h = 0, the jumps of the four stresses move to the right of this figure with the loading depth increasing. This may be related to the thickness of the medium above the loading plane and the loading direction.Figure 10
Variations of the four stress components ford = a and θ = 45° for different h. (a) σxx/p; (b) σyy/p; (c) σzz/p; (d) σxz/p.
(a)(b)(c)(d)


## 4.1. Elastic Parameters of a Nonhomogeneous Medium

In order to consider the variation of Poisson’s ratio with depth, a simple linear fit has been completed for the data by Pan [3], who investigated the depth variation of the geotechnical properties of sand soil, clay loam, clay, and soft soil. The variations of the elastic modulus and Poisson’s ratio for sand soil are presented as follows:(19)Ez=E01+m1z,νz=ν01+m2z,where E0=53.09MPa,m1=0.5065, ν0=0.3469, m2=0.0123, and the SI unit of z is meters (m). The thickness of the sand soil is d = 12 m. For z > d, the elastic properties keep constant; that is, Ez≥d=E01+m1d and νz≥d=ν01+m2d. It is considered that the depth of the nonhomogeneous medium is equal to the radius of the circular loading area; that is, d = a.


## 4.2. The Elastic Fields Induced by Inclined Loads Located at a Given Depth

It is assumed that the circular loading area (radiusa) is located at the depth h = a and is subjected to an inclined load p0, which has an angle θ with the loading plane. Here, we present only the results along the vertical axis passing the center of the loading area. For comparison, two cases of the depth-dependent and depth-independent variations of Poisson’s ratio are analyzed.Figure7 illustrates the displacements ux and uz at the point x,y,z=0,0,z induced by the inclined load p0 for different loading angles. Because of symmetry, the displacement uy at the point x,y,z=0,0,z is equal to zero. It can be found that the variation of Poisson’s ratio with depth exerts an influence on the distribution of the displacements. The influence of Poisson’s ratio is more obvious on the vertical displacement uz. Furthermore, the horizontal displacement ux becomes small with the loading angle θ increasing and is equal to zero for θ = 90°. The vertical displacements uz increase with the loading angle θ increasing and arrive at the maximum values at θ = 90°.Figure 7
Variations of the displacementsux and uz for d = a and h = a for different θ. (a) uxE0/p0a; (b) uzE0/p0a.
(a)(b)Figure8 illustrates the variations of the stresses σxx, σxz, σyy, and σzz at the point x,y,z=0,0,z induced by the inclined load p0 from different loading angles. Because of symmetry, the stresses σxy and σyz at the point x,y,z=0,0,z are equal to zero. From Figure 6, we have the following observations:(1)
All the stresses tend to zero with the distance from the loading plane increasing and have jumps across the loading plane.(2)
The variation of Poisson’s ratio with depth exerts an influence on the distribution of the stresses. The influence of Poisson’s ratio is more obvious on the normal stressesσxx, σyy, and σzz.(3)
With the loading angle increasing, the absolute values ofσxx decrease and the ones of σyy and σzz increase below the loading plane. Above the loading plane, the variations of the three normal stresses with the loading angle are different in different depths. The absolute values of σxz decrease with the loading angle increasing and σxz=0 for θ = 90°.Figure 8
Variations of the four stress components ford = a and h = a for different θ. (a) σxx/p0; (b) σyy/p0; (c) σzz/p0; (d) σxz/p0.
(a)(b)(c)(d)


## 4.3. The Elastic Fields Induced by Inclined Loads Located at Different Depths

It is assumed that the circular loading area (radiusa) is located at different depths (h = 0, 2a, 3a, and 5a) and is subjected to an inclined loadp0, which has an angle θ = 45° with the loading plane. Here, we use the depth-independent model of elastic modulus and Poisson’s ratio and present only the results along the vertical axis passing through the center of the loading area.Figure9 illustrates the displacements ux and uzat the point x,y,z=0,0,z induced by the inclined load p0. It can be found that, for each loading depth, the maximum values of ux and uz appear at the loading plane. And the maximum values of ux and uz decrease with the loading depth increasing. Furthermore, the displacements ux and uz for h = 0 are much larger than the ones for other loading depths. At z = 5a, the displacements of ux and uz for h = 5a are larger than the ones for other loading depths.Figure 9
Variations of the displacementsux and uz for d = a and θ = 45° for different h. (a) uxE0/pa; (b) uzE0/pa.
(a)(b)Figure10 illustrates the variations of the stresses σxx, σxz, σyy, and σzz at the point x,y,z=0,0,z with the loading depth. It can be found that all the jumps of the stresses appear at the loading planes. Apart from the loading plane, all the stresses approach zero. Furthermore, except for h = 0, the jumps of the four stresses move to the right of this figure with the loading depth increasing. This may be related to the thickness of the medium above the loading plane and the loading direction.Figure 10
Variations of the four stress components ford = a and θ = 45° for different h. (a) σxx/p; (b) σyy/p; (c) σzz/p; (d) σxz/p.
(a)(b)(c)(d)


## 5. Conclusions

The problem of the inclined loadings within a nonhomogeneous elastic half-space can serve as a much useful model for examining the response of geological media. In the past analysis, it is generally assumed that the elastic shear modulus of a nonhomogeneous half-space varies with depth and Poisson’s ratio keeps constant with depth. In this paper, it is emphasized that Poisson’s ratio of a nonhomogeneous medium also varies with depth and the loadings have any angle with the loading plane. The results illustrate that the nonhomogeneity of an elastic half-space exerts an obvious influence on the stress and displacement fields. The proposed numerical method can also be used to examine more complicated variations of the elastic parameters with depth.



---

*Source: 1402839-2020-08-24.xml*