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.

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 Figure 1, 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 , and Poisson’s ratio of the medium .

As shown in Figure 1, 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 aswhere P and Q are, respectively, the source and field points of the fundamental solution; and are the kernel functions of the displacements and stresses of a layered medium presented by Yue [10]; 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 area S is discretized into ne quadrilateral elements . 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.

(a)

(b)

(a)
(b)

Figure 2 

The isoparametric element. (a) Linear element with 4 nodes. (b) The quadrilateral element with 8 nodes.

In global and local coordinate systems, there are the following transform relationships of coordinate and traction values within an element:where is the shape function of the element at node α and and 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:where J is the Jacobian determinant.

The interpolation function of the quadrilateral element iswhere , and 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 and 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

The triangle area is transformed into a square area and by the coordinate change of the linear element shown in Figure 2(a). Since the transformation is linear, the relationship between and is established by linear interpolation function:

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 Figure 3.

(a)

(b)

(a)
(b)

Figure 3 

Singularity at node 1or node 5.

When the singular point is at node 1, the element is divided into two triangles: to triangle 1-3-4,and to triangle 1-2-3,

When the singular point is at node 5, the element is divided into three triangles: to triangle 5-4-1,to triangle 5-2-3,and to triangle 5-3-4,

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 by n 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 and 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 (radius a) 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 f_z 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 and h/a = 1.

Figure 5 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 , the largest absolute error of is 0.001 and appears at z = a. Figure 6 illustrates the vertical stress at the points close to the loading plane for different meshes. It is known that the loading causes a jump discontinuity of the stress across the loading plane and the jump at the center of the loading area is equal to. 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.9934 and 0.99743, respectively, and have the absolute errors of 0.0066 and 0.00257, 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 displacement u_z(0, 0, z) for different meshes.

Figure 6 

Variations of 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:where ,, , , 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, and . 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 (radius a) is located at the depth h = a and is subjected to an inclined load p₀, 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.

Figure 7 illustrates the displacements and at the point induced by the inclined load p₀ for different loading angles. Because of symmetry, the displacement at the point 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 u_z. Furthermore, the horizontal displacement u_x becomes small with the loading angle θ increasing and is equal to zero for θ = 90°. The vertical displacements u_z increase with the loading angle θ increasing and arrive at the maximum values at θ = 90°.

(a)

(b)

(a)
(b)

Figure 7 

Variations of the displacements u_x and u_z for d = a and h = a for different θ. (a) u_xE₀/p₀a; (b) u_zE₀/p₀a.

Figure 8 illustrates the variations of the stresses , , , and at the point induced by the inclined load from different loading angles. Because of symmetry, the stresses and at the point 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 , , and .(3)With the loading angle increasing, the absolute values of decrease and the ones of and 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 decrease with the loading angle increasing and for θ = 90°.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 8 

Variations of the four stress components for d = a and h = a for different θ. (a) σ_xx/p₀; (b) σ_yy/p₀; (c) σ_zz/p₀; (d) σ_xz/p_0.

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

It is assumed that the circular loading area (radius a) is located at different depths (h = 0, 2a, 3a, and 5a) and is subjected to an inclined load, 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.

Figure 9 illustrates the displacements and at the point induced by the inclined load . It can be found that, for each loading depth, the maximum values of u_x and u_z appear at the loading plane. And the maximum values of u_x and u_z decrease with the loading depth increasing. Furthermore, the displacements u_x and u_z for h = 0 are much larger than the ones for other loading depths. At z = 5a, the displacements of u_x and u_z for h = 5a are larger than the ones for other loading depths.

(a)

(b)

(a)
(b)

Figure 9 

Variations of the displacements u_x and u_z for d = a and θ = 45° for different h. (a) u_xE₀/pa; (b) u_zE₀/pa.

Figure 10 illustrates the variations of the stresses , , , and at the point 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.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 10 

Variations of the four stress components for d = a and θ = 45° for different h. (a) σ_xx/p; (b) σ_yy/p; (c) σ_zz/p; (d) σ_xz/p.

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.

Abbreviations

Data Availability

All data analyzed during this study are included in the figures in this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

The work presented in this paper was funded by the National Natural Science Foundation of China (Grant no. 41672291), SDUST Research Fund (2015TD-JH104), and the Opening Foundation of Shandong Key Laboratory of Civil Engineering Disaster Prevention and Mitigation (CDPM2019ZR07). The authors deeply thank Professor A.P.S. Selvadurai from McGill University, Canada, for supplying the data of their research in Figure 4.