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

State Key Laboratory of Mining Disaster Prevention and Control Cofounded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266590, China Shandong Key Laboratory of Civil Engineering Disaster Prevention and Mitigation, College of Civil Engineering & Architecture, Shandong University of Science and Technology, Qingdao 266590, China Department of Civil Engineering, +e University of Hong Kong, Hong Kong, China


Introduction
e 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. e nonhomogeneity problem in geomechanics has applications to many problems of technological importance.
ere have been numerous efforts attempting to understand the influence of nonhomogeneity on the response of an elastic half-space induced by various loads [1][2][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][5][6][7][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. e 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. e numerical method used for analysis is developed by applying the fundamental solution of layered elastic solids [10][11][12][13][14][15][16][17] and integrating it over the loading area. e 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.

Basic Equations.
e 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. e dissimilar homogeneous layers adhere to an elastic solid of upper semi-infinite extent and another elastic solid of lower semi-infinite extent. e 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. e 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 E 0 � 1.0 × 10 −15 MPa, and Poisson's ratio of the medium ] 0 � 0.3.
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 as where P and Q are, respectively, the source and field points of the fundamental solution; u * ik (Q, P) and σ * ijk (Q, P) are the kernel functions of the displacements and stresses of a layered medium presented by Yue [10]; t k (P) is the traction of the source point P; and the integral domain S is the loading area. e loading area S is discretized into ne quadrilateral elements S � ne e�1 S e . 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.

Numerical
In global and local coordinate systems, there are the following transform relationships of coordinate and traction values within an element: where N α is the shape function of the element at node α and x α i and t α i 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. e interpolation function of the quadrilateral element is 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.

Numerical Integration.
e 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 u * ik 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.
e coordinates of any point in the triangle can be expressed as e 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: 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, Advances in Materials Science and Engineering 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. (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].

Parallel Computing and Discretization of the Nonhomogeneous Layer. Since the integrals in equations
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. e nonhomogeneous layer in a halfspace 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.

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 G(z) � G 0 e λz (z ≤ d) and G(z) � G 0 e λd (z ≥ 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. e 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 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 0 ≤ z < 5a, the largest absolute error of u z (0, 0, z)/f z 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 f z causes a jump discontinuity of the stress σ zz across the loading plane and the jump at the center of the loading area is equal tof z . 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.9934f z and 0.99743f z , respectively, and have the absolute errors of 0.0066f z and 0.00257f z , respectively. ese calculation errors are induced by nearly singular integrals and it is difficult to completely eliminate the errors by the subelement method above.
us, 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.

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. e variations of the elastic modulus and Poisson's ratio for sand soil are presented as follows: where E 0 � 53.09 MPa, m 1 � 0.5065, ] 0 � 0.3469, m 2 � 0.0123, and the SI unit of z is meters (m). e thickness of the sand soil is d � 12 m. For z > d, the elastic properties keep constant; that is, It is considered that the depth of the nonhomogeneous medium is equal to the radius of the circular loading area; that is, d � a.

e 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 = E 0 (1 + m 1 z) and Poisson's ratio keeps constant = E 0 (1 + m 1 z) and Poisson's ratio keeps constant    Advances in Materials Science and Engineering inclined load p 0 , 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 u x and u z at the point (x, y, z) � (0, 0, z) induced by the inclined load p 0 for different loading angles. Because of symmetry, the displacement u y 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. e 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°. e vertical displacements u z increase with the loading angle θ increasing and arrive at the maximum values at θ � 90°. Figure 8 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 p 0 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.

e 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 loadp 0 , 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 u x and u z at the point (x, y, z) � (0, 0, z) induced by the inclined load p 0 . 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. Figure 10 illustrates the variations of the stresses σ xx , σ xz , σ yy , and σ zz at the point (x, y, z) � (0, 0, z) with the loading = E 0 (1 + m 1 z) and Poisson's ratio keeps constant   Figure 10: Continued.   Figure 9: Variations of the displacements u x and u z for d � a and θ � 45°for different h. (a) u x E 0 /pa; (b) u z E 0 /pa. 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. is may be related to the thickness of the medium above the loading plane and the loading direction.

Conclusions
e 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. e results illustrate that the nonhomogeneity of an elastic half-space exerts an obvious influence on the stress and displacement fields. e proposed numerical method can also be used to examine more complicated variations of the elastic parameters with depth.

P, Q:
Source and field points of the fundamental solution u * ik (Q, P): Kernel functions of the displacements of a layered medium σ * ijk (Q, P): Kernel functions of the stresses of a layered medium t k (P): Traction of the source point P S: Integral domain, loading area N α : Shape function of the element at node α x α i , t α i : Coordinate and traction values of the element at node α ξ a , η a : Local coordinate component at node α J: Jacobian determinant a: Radius of the circular loading area ]: Poisson's ratio E: Elastic modulus θ: Load angle with the loading plane p 0 : Inclined load (σ xz , σ yz , σ xz ): Stresses induced by the inclined load p 0 (u x , u y , u z ): Displacements induced by the inclined load p 0 G 0 : Shear modulus λ: Constant.

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