#### 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 . 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 . 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.

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

##### 2.1. Basic Equations

The fundamental solution of a multilayered solid developed by Yue  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  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 , 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 ; is the traction of the source point P; and the integral domain S is the loading area.

##### 2.2. Numerical Scheme

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.

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.

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 .

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  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 fz at h = a are presented for the approximation by depth-dependent variations of the elastic parameters.

#### 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 , 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 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.

Figure 7 illustrates the displacements and at the point induced by the inclined load p0 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 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°.

##### 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.

#### 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.