Complex Variable Solution for Stress and Displacement of Layered Soil with Finite Thickness

: The static problem of a layered isotropic elastic body is a very useful research subject in relation to the analysis and design of foundation works. Due to the complexity of the problem, there is no analytical solution to the problem so far. This study provides an efﬁcient analytical approach to accurately calculate the displacement and stress ﬁelds of the soil. The constraints of bedrock on soil, different soil layer thickness and the shear stress of the foundation on soil were all taken into account in the analysis. In this study, each layer is regarded as an isotropic elastomer with inﬁnite width, and the layers are in complete contact. By using conformal mapping, each layer is mapped to a unit circle, and the two complex potential functions are expanded into Taylor series with unknown coefﬁcients. These unknown coefﬁcients are obtained by satisfying boundary conditions and continuity conditions. The boundary and continuity conditions were veriﬁed in this paper. As a validation step, we compared the analytical results for the settlement with the results of the ANSYS numerical simulations and found good agreement. Parametric analyses were also carried out to investigate the inﬂuence of different distribution forms of base pressure on surface settlement, and the effects of layered properties on the surface settlement and stress ﬁeld.


Introduction
In many places, it is common to see compressible soil overlying rigid bedrock. In this case, the stress and displacement of soil are obviously different from that of a semi-infinite body medium [1]. Natural soils tend to be horizontally layered during deposition, the density and hardness of each layer are generally different, and the deformation characteristics may vary greatly [2]. When the stress in the soil is too large, the deformation of the soil will cause unacceptable settlement of the building and even the instability of the entire soil. Therefore, an accurate calculation of soil stress and displacement is an important basis for the stability analysis of building foundations and geotechnical structures [3].
Finite element method can provide useful solutions to many complex geological conditions. However, extensive computational time is required, especially when complete parametric analyses need to be performed [4]. Analytical solutions provide an efficient and quick approach to gain insight into the nature of the problem [5]; in addition, they enable analyses of a wide range of parameter values so that the physics of the problem can be better understood [4]. In the past, researchers mainly used the integral transformation method to analyze the elastic static problems of layered soil. Burmister [6] provided the integral expression of stress and displacement of two layers of soil under a flexible circular foundation based on the integral transformation method, and also established an empirical formula for surface settlement when both layers of soil consist of incompressible materials. In many studies in the literature, the soil surface is subjected to a circular uniformly distributed load, which is considered a spatial axisymmetric problem. Both layers of soil

Problem Statement and Method Presentation
This paper deals with the case of complete contact between layers of multi-layer soil and complete constraint with bedrock below it. When the soil surface is subjected to distributed load, the deformation of bedrock is considered to be negligible compared with the deformation of soil, that is, the displacement component on the interface between the soil and bedrock is equal to zero. Figure 1a shows a problem in which the upper boundary is the stress boundary condition, the lower boundary is the displacement boundary condition, and the contact surface of each layer is the stress and displacement continuity condition.
In Figure 1a, the four end points (A, C, D, F) of each layer of soil are points at infinity. t is any point on the contact surface. There is a vector relation about t: Appl. Sci. 2022, 12, 766

of 22
As shown in Figure 1, the j-th layer (j = 1,2,· · · ,n) soil is mapped to the inner domain of the unit circle in the ζ j plane. The mapping function is: where z j is any point in the j-th layer soil, z j = x j + iy j , i = √ −1; H j is the thickness of the j-th layer soil. ζ j is the point where z j maps to the image plane. ζ j = ρ j e iθ j , ρ j = 1 corresponds to the upper and lower boundaries in the j-th layer soil, and σ j is used to denote ζ j at this time. In Figure 1a, the four end points (A, C, D, F) of each layer of soil are points at infinity. is any point on the contact surface. There is a vector relation about : As shown in Figure 1, the j-th layer (j = 1,2,⋯,n) soil is mapped to the inner domain of the unit circle in the ζj plane. The mapping function is: where is any point in the j-th layer soil, = + , = √−1; is the thickness of the j-th layer soil. ζj is the point where maps to the image plane. ζj = e , = 1 corresponds to the upper and lower boundaries in the j-th layer soil, and is used to denote ζj at this time.
As shown in Figure 1, in the ζ1 plane, the circumference of the unit circle (│ │ = 1) corresponds to the upper and lower boundaries of the first layer of soil. A , B , C , D , As shown in Figure 1, in the ζ 1 plane, the circumference of the unit circle (σ 1 = 1) corresponds to the upper and lower boundaries of the first layer of soil. A 1 , B 1 , C 1 , D 1 , E 1 , F 1 are mapped to A 1 , B 1 , C 1 , D 1 , E 1 , F 1 , respectively. The upper boundary (A 1 B 1 C 1 ) is mapped to the upper semicircle ( A 1 B 1 C 1 ). The lower boundary (D 1 E 1 F 1 ) is mapped to the lower semicircle ( D 1 E 1 F 1 ). The coordinate origin (O 1 ) is mapped to the center of the unit circle (O 1 ) of the ζ 1 plane. On the upper and lower boundaries, all points of x 1 → +∞ are mapped to the same point (σ 1 = 1) of the ζ 1 plane, and all points of x 1 → −∞ are mapped to the same point (σ 1 = −1) of the ζ 1 plane. The coordinate axes x 1 and y 1 are mapped to ξ 1 and η 1 axes. The points a and b are mapped to a and b . The rest of the layers can be deduced by analogy.
By substituting Equation (2) into Equation (1), the relationship between θ 1 and θ 2 at the same point t on the contact surface of the 1st layer and the 2nd layer can be obtained: In the same way, on the contact surface of the j-th layer and the j + 1-th layer: σ j,x , σ j,y and τ j,xy are the horizontal normal stress, vertical normal stress and shear stress of the j-th layer soil, respectively, and u j and v j are the horizontal displacement and vertical displacement of the j-th layer soil, respectively. The three stress components σ j,x , σ j,y , τ j,xy and two displacement components u j , v j of point z j can be expressed as [21]: where κ j = 3 − 4µ j , and Re[ . . . ] means to take the real part of [ . . . ], the imaginary part represented by Im[ . . . ] will also be involved later. G j and µ j are the shear modulus and Poisson's ratio of the j-th layer soil, respectively. The single valued analytic functions of the j-th layer are ϕ 1,j z j and ψ 1,j z j . The superscripts (.)' and (.)" denote the first and second derivatives of (.) with respect to z j . For the problem analyzed in this paper, the stress boundary conditions can be obtained from Equations (5) and (6): where z 1 is the point on the boundary of A 1 B 1 C 1 , p x (x 1 ) is the horizontal load and p y (x 1 ) is the vertical load.
The continuity condition of displacement and stress on the contact surface between the lower boundary of the j-th layer and the upper boundary of the j + 1-th layer can be obtained from Equations (5)-(7).
In Equations (9) and (10), j = 1, 2, . . . , n − 1. z j and z j+1 are the same point on the contact surface, where z j is on D j E j F j , and z j+1 is on A j B j C j .
Assuming that the bottom boundary is completely constrained, i.e., the horizontal and vertical displacements are equal to zero, the displacement boundary condition on D n E n F n can be obtained from Equation (7): Substituting Equation (2) into ϕ 1,j z j and ψ 1,j z j , we can obtain ϕ j ζ j and ψ j ζ j . In the unit circle, they can be expressed as follows: where a j,k and b j,k are the complex coefficients to be solved: In Equation (13), a j,k,R and a j,k,I are the real part and imaginary part of a j,k respectively; b j,k,R and b j,k,I are the real part and imaginary part of b j,k respectively. This can be solved by Equations (8)- (11).
On the boundary, ζ j is denoted by σ j , Equations (8)- (11) in the image plane can be rewritten as 2Re In Equation (14) In Equations (15) and (16), j = 1, . . . , n − 1: In Equation (17) In this paper, the power series method is used to solve ϕ j ζ j and ψ j ζ j .
This paper adopts the power series method, ϕ j ζ j and ψ j ζ j are taken as finite terms, and the highest powers are n j,1 and n j,2 respectively.
Equations (31) and (32) can be obtained from the real and imaginary parts of Equation (27).
Equations (33) and (34) can be obtained from the real and imaginary parts of Equation (28).
Re l j,ks b j,k,I + κ j G j a j,0,I + 1 Re g j+1,ks + h j+1,ks a j+1,k,I + Equations (35) and (36) can be obtained from the real and imaginary parts of Equation (29).
Re 2c j,ks + d j,ks − e j,ks a j,k,R + n j,1 ∑ k=1 Im −2c j,ks − d j,ks + e j,ks a j,k,I + Re d j,ks − e j,ks a j,k,I + Equations (37) and (38) can be obtained from the real and imaginary parts of Equation (30).
Im[g n,ks + h n,ks ]a n,k,I − Im[g n,ks − h n,ks ]a n,k,R + n n,1 ∑ k=1 Re[g n,ks + h n,ks ]a n,k,I − Re[l n,ks ]b n,k,I + κ n G n a n,0,I + 1 G n b n,0,I = 0 Since a j,0 and b j,0 do not affect the stress, they represent the translation of the rigid body. If a j,0 or b j,0 is equal to zero, it will not affect the calculation result of displacement, so in this paper, b j,0 = 0.
The system of linear equations of a j,k,R , a j,k,I , b j,k,R , b j,k,I can be obtained from Equations (31)-(38).
where: X = · · · · · · a j,1,R a j,1,I · · · a j,n j,1 ,R a j,n j,1 ,I b j,1,R b j,1,I · · · b j,n j,2 ,R b j,n j,2 ,I a j,0,R a j,0,I 2n j,1 +2n j,2 +2 coefficients in the j−th layer There are ∑ n j=1 2 m j + m j+1 + 2 terms in matrix F. The elements in coefficient matrix D are coefficients of Equations (29)-(36) for a j,k,R , a j,k,I , b j,k,R , b j,k,I , and the size is row ∑ n j=1 2 m j + m j+1 + 2 and column ∑ n j=1 2 n j,1 + n j,2 + 1 . In order to ensure the accuracy of calculation, n j needs to take a large value. After calculation, it is found that the result obtained by Equation (39) may be very unsatisfactory. In order to obtain satisfactory results, this paper uses ordinary least squares. Taking ∑ n j=1 2 m j + m j+1 + 2 > ∑ n j=1 2 n j,1 + n j,2 + 1 : The analytic functions ϕ j ζ j and ψ j ζ j can be obtained from Equation (40). The displacement and stress of any point in the j-th layer can be calculated by Equations (41) and (42), respectively.
From Equations (41) and (42), we have where where t j,1 t j,8 are functions of ζ j ; the details are provided in Appendix B.

Analysis and Discussion
In order to ensure accuracy, the values of m j , n j,1 , n j,2 in this paper are large, which means matrix D is too large, thus we took two layers of soil as an example. The example only discusses the effect of a smooth foundation on soil, which means p x (x 1 ) = 0.
According to the given known conditions, the solution process in Section 3 can be used to calculate the ϕ j ζ j and ψ j ζ j of each layer (j = 1, 2), and then the stress components σ 1,y and τ 1,xy of the surface can be obtained by Equation (42). If σ 1,y = p y (x 1 ), τ 1,xy = p x (x 1 ) = 0, the stress boundary condition is satisfied. It can be seen from the calculation that the value of τ 1,xy on the surface is very small, basically equal to zero, and σ 1,y is also basically equal to p y (x 1 ). Figure 2 shows the distribution curve of σ 1,y and p y (x 1 ) at the surface. It can be seen that the vertical stress curve of the surface basically coincides with the load curve, which indicates that the stress boundary conditions on the surface are well satisfied. Appl. Sci. 2022, 12, x FOR PEER REVIEW 11 of 24 According to the calculated ( ) and ( ), the displacement and stress of each layer on the contact surface can also be obtained through Equations (41) and (42), therefore the continuity condition of the displacement and stress can be verified. Figure 3 shows the verification of the displacement continuity condition. It can be seen that ≈ , ≈ , where and are the x-direction displacement and y-direction displacement of the lower boundary of the first layer of soil, respectively, while and are the x-direction displacement and y-direction displacement of the same position of the upper boundary of the second layer of soil, respectively. Figure 4 shows the verification of the stress continuity condition, which shows that , ≈ , , , ≈ , . It can be seen from Figure 3 and Figure 4 that the displacement and stress continuity condition can also be satisfied very well. The maximum displacement in Figure 3 is −0.255 mm, and compared with this, the maximum displacement of the lower boundary of the second layer is 1.3791294 × 10 mm, which means the displacement boundary conditions of the lower boundary of the second layer are also very well satisfied. According to the calculated ϕ j ζ j and ψ j ζ j , the displacement and stress of each layer on the contact surface can also be obtained through Equations (41) and (42), therefore the continuity condition of the displacement and stress can be verified. Figure 3 shows the verification of the displacement continuity condition. It can be seen that u 1 ≈ u 2 , v 1 ≈ v 2 , where u 1 and v 1 are the x-direction displacement and y-direction displacement of the lower boundary of the first layer of soil, respectively, while u 2 and v 2 are the x-direction displacement and y-direction displacement of the same position of the upper boundary of the second layer of soil, respectively. Figure 4 shows the verification of the stress continuity condition, which shows that σ 1,y ≈ σ 2,y , τ 1,xy ≈ τ 2,xy . It can be seen from Figures 3 and 4 that the displacement and stress continuity condition can also be satisfied very well. The maximum displacement in Figure 3 is −0.255 mm, and compared with this, the maximum displacement of the lower boundary of the second layer is 1.3791294 × 10 −4 mm, which means the displacement boundary conditions of the lower boundary of the second layer are also very well satisfied.

Comparison with Numerical Method
In Section 4.1, the surface stress boundary conditions, the displacement boundary condition of the lower boundary of the second layer and the stress and displacement continuity conditions between layers were verified. It was found that the boundary condition and the continuity conditions were very well satisfied; however, this is not enough to establish that the whole derivation process presented in this paper is correct. In order to verify the correctness of the whole derivation process, it is necessary to compare the calculation results of this paper with those of the numerical method. When the results of the two methods are in good agreement, the whole derivation process and the calculation program used in this paper are correct.
The schematic diagram of the ANSYS [30] model is shown in Figure 5. The model is divided into upper and lower soil layers. The width of the upper soil model is 100 m, the thickness is 30 m, the load action interval is −2m ≤ x 1 ≤ 2m, y = 0, the element type is plane42, the elastic modulus is 2.4 MPa, the Poisson's ratio is 0.2, The element size in most areas is 0.5 (shown in Figure 6), and the element size is 0.05 in the area of −2m ≤ x 1 ≤ 2m, −10m ≤ y ≤ 0 (shown in Figure 7),. The width of the lower soil model is 100 m, the thickness is 20 m, the element type is plane42, the elastic modulus is 5.  Figure 8 shows the ANSYS deformation diagram. Figure 9 shows a comparison between the analytical solution and numerical solution of the surface settlement. ANSYS version 15.0 was used in this paper.

Comparison with Numerical Method
In Section 4.1, the surface stress boundary conditions, the displacement boundary condition of the lower boundary of the second layer and the stress and displacement continuity conditions between layers were verified. It was found that the boundary condition and the continuity conditions were very well satisfied; however, this is not enough to establish that the whole derivation process presented in this paper is correct. In order to verify the correctness of the whole derivation process, it is necessary to compare the calculation results of this paper with those of the numerical method. When the results of the two methods are in good agreement, the whole derivation process and the calculation program used in this paper are correct.
The schematic diagram of the ANSYS [30] model is shown in Figure 5. The model is divided into upper and lower soil layers. The width of the upper soil model is 100 m, the thickness is 30 m, the load action interval is −2 ≤ ≤ 2 , = 0, the element type is plane42, the elastic modulus is 2.4 , the Poisson's ratio is 0.2, The element size in most areas is 0.5 (shown in Figure 6), and the element size is 0.05 in the area of − 2 ≤ ≤ 2 , − 10 ≤ ≤ 0 (shown in Figure 7),. The width of the lower soil model is 100 m, the thickness is 20 m, the element type is plane42, the elastic modulus is 5.  Figure 8 shows the ANSYS deformation diagram. Figure 9 shows a comparison between the analytical solution and numerical solution of the surface settlement. ANSYS version 15.0 was used in this paper.

Influence of Base Pressure Distribution on Surface Settlement
As shown in Figure 10a, three types of base pressure distribution under vertical action are discussed: (1) The concentration of base pressure is large at the edge: p y1 (x 1 ) = −3x 2 1 − 12 /16; (2) The concentration of base pressure is large at the center: p y2 (x 1 ) = −3 4 − x 2 1 /8; (3) Uniform distribution of base pressure: p y3 (x 1 ) = −1.

Influence of Base Pressure Distribution on Surface Settlement
As shown in Figure 10a, three types of base pressure distribution under vertical tion are discussed: In order to compare the influence of load distribution on surface settlement, the sultant forces of the three types of base pressure distribution are the same, and they all 4 . The other parameters are the same as those described in Section 4.1, and surface settlement that they caused are shown in Figure 10b,c. It can be seen from Figure 10 that the surface settlement caused by the three type base pressure distribution only shows certain differences near the foundation, and maximum surface settlement occurs in the center of the foundation action area. When load concentration at the center of the foundation is large, the surface settlement is largest, followed by uniform distribution, and the surface settlement is smallest when load concentration at the edge of the foundation is large. At the same time, it can be s from Figure 10c that the surface settlement caused by the first distribution is appr mately a horizontal line in the foundation action area. Only when the strip foundatio absolutely rigid and it is subjected to axial compression, can the bottom of the founda sink evenly. Therefore, we can infer that the distribution form of the base pressure un the absolutely rigid strip foundation is close to this distribution form. In order to compare the influence of load distribution on surface settlement, the resultant forces of the three types of base pressure distribution are the same, and they are all 4 kPa. The other parameters are the same as those described in Section 4.1, and the surface settlement that they caused are shown in Figure 10b,c.
It can be seen from Figure 10 that the surface settlement caused by the three types of base pressure distribution only shows certain differences near the foundation, and the maximum surface settlement occurs in the center of the foundation action area. When the load concentration at the center of the foundation is large, the surface settlement is the largest, followed by uniform distribution, and the surface settlement is smallest when the load concentration at the edge of the foundation is large. At the same time, it can be seen from Figure 10c that the surface settlement caused by the first distribution is approximately a horizontal line in the foundation action area. Only when the strip foundation is absolutely rigid and it is subjected to axial compression, can the bottom of the foundation sink evenly. Therefore, we can infer that the distribution form of the base pressure under the absolutely rigid strip foundation is close to this distribution form.

Influence of Layered Soil on Surface Settlement and Additional Stress Distribution in Soil
When H 1 = H 2 = 30 m, G 1 = 1.0 MPa, µ 1 = µ 2 = 0.2, the distribution of base pressure is p y (x 1 ) = −3 4 − x 2 1 /8, and other parameters are the same as those described in Section 4.1. Figure 11 shows the surface settlement curve when G 2 /G 1 is 0.1, 0.5, 1.0, 2.0 and 10.0, respectively.

Influence of Layered Soil on Surface Settlement and Additional Stress Distribution in Soil
When = = 30 , = 1.0 , = = 0.2, the distribution of base pressure is ( ) = − 3(4 − ) 8 ⁄ , and other parameters are the same as those described in Section 4.1. Figure 11 shows the surface settlement curve when / is 0.1, 0.5, 1.0, 2.0 and 10.0, respectively. It can be seen from Figure 11 that the smaller the shear modulus ratio / , the greater the surface settlement. The influence of the change in / on the surface settlement can be divided into two scenarios: when / > 1, the change in / has little influence on the surface settlement; when / < 1, the change in / has a great influence on the surface settlement. Figure 12 shows the , curves (solid line) of the lower boundary of the first layer when / is 0.1, 0.5, 1.0, 2.0 and 10.0 and the , curves (dotted line) of the upper boundary of the second layer when / is 0.1, and 10.0. Figure 13 shows the , curve on the contact surface, and Figure 14 shows the , curve on the contact surface. It can be seen from Figure 11 that the smaller the shear modulus ratio G 2 /G 1 , the greater the surface settlement. The influence of the change in G 2 /G 1 on the surface settlement can be divided into two scenarios: when G 2 /G 1 > 1, the change in G 2 /G 1 has little influence on the surface settlement; when G 2 /G 1 < 1, the change in G 2 /G 1 has a great influence on the surface settlement. Figure 12 shows the σ 1,x curves (solid line) of the lower boundary of the first layer when G 2 /G 1 is 0.1, 0.5, 1.0, 2.0 and 10.0 and the σ 2,x curves (dotted line) of the upper boundary of the second layer when G 2 /G 1 is 0.1, and 10.0. Figure 13 shows the σ 1,y curve on the contact surface, and Figure 14 shows the τ 1,xy curve on the contact surface. Appl. Sci. 2022, 12, x FOR PEER REVIEW 18 of 24      It can be seen from Figure 12 that with the decrease in the shear modulus ratio / of the two layers of soil, the maximum value of , (absolute value) of the lower boundary of the first layer of soil decreases at first and then increases, and the tensile stress occurs. Moreover, with the decrease in / , the tensile stress increases. Compared with the first layer, , on the upper boundary of the second layer is very small and the change in / has little effect on , . Therefore, this paper only shows the , curves when / is 0.1, 10.0 respectively. As can be seen from Figure 13, when the soft soil layer is covered by the hard soil layer ( / > 1), near = 0, the , (absolute value) on the contact surface is larger than that of homogeneous soil ( / = 1), but in the distance, it is smaller than that of homogeneous soil, that is, the stress concentration phenomenon occurs near = 0. When the hard soil layer covers the soft soil layer ( / < 1), the , on the contact surface is smaller near = 0 than that of the homogeneous soil layer, but larger in the distance than that of the homogeneous soil layer, which is known as the stress diffusion phenomenon. At this time, the vertical normal stress distribution on the contact surface is relatively uniform, so the vertical deformation of the lower soil will be relatively uniform, but it can be seen from Figure 12 that the horizontal tensile stress of the lower boundary of the upper soil will be relatively large.
It can be seen from Figure 14 that with the decrease in / , the shear stress , (absolute value) on the contact surface gradually decreases. When / > 1, the change in / has little effect on the shear stress distribution on the contact surface, while when / < 1, the change in / has a greater effect on the shear stress. Taking the same parameters as described in Section 4.1, since the base pressure is symmetrical about the y-axis, there are only two normal stress components on the y-axis, and the shear stress component is 0. Figure 15 shows the variation in the two normal stresses on the y-axis with depth ℎ. At the surface, the horizontal normal stress is basically equal to the vertical normal stress. With the increase in ℎ, the horizontal and vertical It can be seen from Figure 12 that with the decrease in the shear modulus ratio G 2 /G 1 of the two layers of soil, the maximum value of σ 1,x (absolute value) of the lower boundary of the first layer of soil decreases at first and then increases, and the tensile stress occurs. Moreover, with the decrease in G 2 /G 1 , the tensile stress increases. Compared with the first layer, σ 2,x on the upper boundary of the second layer is very small and the change in G 2 /G 1 has little effect on σ 2,x . Therefore, this paper only shows the σ 2,x curves when G 2 /G 1 is 0.1, 10.0 respectively.
As can be seen from Figure 13, when the soft soil layer is covered by the hard soil layer (G 2 /G 1 > 1), near x = 0, the σ 1,y (absolute value) on the contact surface is larger than that of homogeneous soil (G 2 /G 1 = 1), but in the distance, it is smaller than that of homogeneous soil, that is, the stress concentration phenomenon occurs near x = 0. When the hard soil layer covers the soft soil layer (G 2 /G 1 < 1), the σ 1,y on the contact surface is smaller near x = 0 than that of the homogeneous soil layer, but larger in the distance than that of the homogeneous soil layer, which is known as the stress diffusion phenomenon. At this time, the vertical normal stress distribution on the contact surface is relatively uniform, so the vertical deformation of the lower soil will be relatively uniform, but it can be seen from Figure 12 that the horizontal tensile stress of the lower boundary of the upper soil will be relatively large.
It can be seen from Figure 14 that with the decrease in G 2 /G 1 , the shear stress τ 1,xy (absolute value) on the contact surface gradually decreases. When G 2 /G 1 > 1, the change in G 2 /G 1 has little effect on the shear stress distribution on the contact surface, while when G 2 /G 1 < 1, the change in G 2 /G 1 has a greater effect on the shear stress.
Taking the same parameters as described in Section 4.1, since the base pressure is symmetrical about the y-axis, there are only two normal stress components on the y-axis, and the shear stress component is 0. Figure 15 shows the variation in the two normal stresses on the y-axis with depth h. At the surface, the horizontal normal stress is basically equal to the vertical normal stress. With the increase in h, the horizontal and vertical normal stress decay rapidly, and the decay speed of the horizontal normal stress is faster than that of the vertical normal stress. The horizontal normal stress at the surface is about −1.5, and it basically decays to zero at the depth of h = 4 m. At the same time, it can be seen that due to the different material properties of the two layers of soil, the horizontal normal stress curve jumps at h = 30 m (contact surface).
normal stress decay rapidly, and the decay speed of the horizontal normal stress is faster than that of the vertical normal stress. The horizontal normal stress at the surface is about −1.5, and it basically decays to zero at the depth of ℎ = 4 . At the same time, it can be seen that due to the different material properties of the two layers of soil, the horizontal normal stress curve jumps at ℎ = 30 (contact surface).

Conclusions
The method proposed in this paper can be used to calculate the stress and displacement of multi-layer soil above bedrock under the action of a strip foundation. In fact, it is an analytical method, and the accuracy of its calculation results is related to the number of terms taken by the analytic function. Through the programming calculation, it was found that when the number of terms of the analytic function is 550, the accuracy of the calculation results is quite high, which can well meet the stress boundary condition, displacement boundary condition, and the stress and displacement continuity conditions. The results are in good agreement with those obtained by the finite element method.
In this paper, the effects of different distribution forms of base pressure and layered soil on the surface settlement and stress distribution inside the soil were analyzed by numerical examples. The results show that the surface settlement curves corresponding to the three types of basement pressure distribution are basically the same, but there are some differences near the foundation action area. When the distribution of the base pressure reaches the maximum in the center of the base, the displacement is at its maximum. The displacement caused by uniform distribution takes second place. When the base pressure is largest on both sides of the base, the displacement is the smallest. In the foundation action area, the surface settlement curve is approximately a horizontal straight line, which

Conclusions
The method proposed in this paper can be used to calculate the stress and displacement of multi-layer soil above bedrock under the action of a strip foundation. In fact, it is an analytical method, and the accuracy of its calculation results is related to the number of terms taken by the analytic function. Through the programming calculation, it was found that when the number of terms of the analytic function is 550, the accuracy of the calculation results is quite high, which can well meet the stress boundary condition, displacement boundary condition, and the stress and displacement continuity conditions. The results are in good agreement with those obtained by the finite element method.
In this paper, the effects of different distribution forms of base pressure and layered soil on the surface settlement and stress distribution inside the soil were analyzed by numerical examples. The results show that the surface settlement curves corresponding to the three types of basement pressure distribution are basically the same, but there are some differences near the foundation action area. When the distribution of the base pressure reaches the maximum in the center of the base, the displacement is at its maximum. The displacement caused by uniform distribution takes second place. When the base pressure is largest on both sides of the base, the displacement is the smallest. In the foundation action area, the surface settlement curve is approximately a horizontal straight line, which proves that the base pressure of a vertical sinking rigid foundation is similar to this distribution form.
For two-layer soil, the surface settlement decreases with the increase in the ratio of shear modulus (G 2 /G 1 ), G 1 and G 2 are the shear modulus of the upper and lower soil, respectively. When the G 2 /G 1 increases: (1) The maximum horizontal normal stress of the lower boundary of the upper soil first decreases and then increases, and when G 2 /G 1 reaches a certain value, the tensile stress begins to appear; (2) The maximum vertical normal stress on the contact surface increases, while the vertical normal stress decreases at a distance from the foundation; (3) The shear stress on the contact surface only changes significantly when G 2 < G 1 .
Author Contributions: A.L. conceived and designed the study. X.S. derived the formulas and wrote the programs. X.S. wrote the paper. H.C. and C.Y. helped with drawing and translation. A.L., H.C. and C.Y. reviewed the edited manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
The study is supported by the National Natural Science Foundation of China (Grant no.51974124).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data in the present study are available from the corresponding author upon reasonable request.