#### Abstract

In this study, a contact surface constitutive model with zero-thickness unit and variable shear stiffness was improved based on the statistical damage constitutive model. The model parameters were derived by means of the shear stress-shear displacement curve of the pile-soil contact surface, and the model parameters under different normal stresses were obtained by the linear interpolation method. At the same time, the influence of the interpolation interval range on the model parameters was explored. The shear stiffness adjustment factor was introduced, and the improved pile-soil contact surface constitutive model was applied to the numerical simulation of pile-soil contact surface shear calculations using the fish language embedded in FLAC 3D, and the variation of monopile bearing characteristics and ultimate bearing capacity were investigated and analyzed. The results show that the improved contact surface constitutive model is able to reflect the nonlinear variation of shear stiffness, and that different normal stresses correspond to different fitting parameters, demonstrating the depth effect of the model. The accuracy of the model parameters reduces as the interpolation interval increases, and the interpolation results are more accurate when the interval range is smaller. The numerical model accurately simulates the pile-soil contact surface shear calculation and the monopile bearing calculation, and the simulation results of the ultimate pile bearing capacity are closer to the results computed by the equations in the Chinese code. At the same time, the variation law of pile axial force and pile lateral frictional resistance along the depth direction and the variation of pile ultimate bearing capacity under different working conditions are reasonable, which shows the validity of the contact surface principal structure model and the reasonableness of numerical calculation in this study.

#### 1. Introduction

Pile foundations have a long history and a wide range of applications in the field of civil engineering. Whether a bridge erection is being built, a foundation pit is being supported, or a housing is being constructed, a pile foundation will be used. Scholars currently employ a variety of research methods, including field tests [1, 2], model tests [3, 4], numerical simulations [5–7], and others, to investigate numerous macroscopic problems of pile foundations, including settlement [8–10], deformation [11, 12], and bearing characteristics [13–15], and various explanations are put forward for pile-soil interactions [16–18]. Due to the considerable differences in soil properties, the mechanical properties of the contact surface are also complicated and variable. At the same time, since the pile-soil contact surface has a significant effect on pile foundation bearing capacity and bearing characteristics, it is required to conduct an in-depth analysis of the pile-soil contact surface.

A contact surface implies the presence of two different materials in touch and interacting, so that their mechanical properties are determined by the combined influence of the two materials [19]. Coulomb proposed an Earth pressure theory as early as the 18th century, based on the friction between soil and walls in practical engineering, which essentially represents the strength attributes of the contact surface within the soil. Since piles are made of concrete, early pile-soil contact surface studies were mainly carried out through soil-concrete contact surface shear tests. Potyondy [20] conducted hundreds of tests to determine the magnitude of contact surface friction and discovered that soil type, moisture content, roughness, and normal stress significantly affect the contact surface. Shi et al. [21] utilized large cyclic load direct shear tests to evaluate the shear characteristics of the contact surfaces of red clay and concrete under a different number of cycles. Xiong et al. [22] conducted a shear test at the concrete-frozen soil interface and found that the initial moisture content and temperature had a substantial effect on the shear behavior of the frozen pile-soil interface. In addition, experimental studies on pile-soil contact surface shear have been undertaken to determine the soil particle size [23, 24], roughness [25], shear rate [26], and shear path [27].

According to the stress-strain relationship during shear, several conventional contact surface instantiation models have been developed, including the hyperbolic model [28], the elasto-plastic model [29], the rigid-plastic model [30], and the damage model [31]. In response to the results of contact surface experimental studies, the contact surface element theory was proposed to better explain the mechanical behavior of contact surface shear. The contact surface element theory is mainly divided into two types: contact surface element without thickness [32] and element with thickness [33]. Numerous scholars have proposed the contact surface constitutive models under a variety of conditions based on the early investigations [34–37], and numerous studies on contact surface units have been done [38–40]. With the help of continuous strength theory, Yang [41] developed the statistical damage constitutive model of the soil-structure interface from the randomness of the internal damage distribution. This study will expand on the research using this model.

With the advancement of measurement technology and test equipment, related scholars have conducted in-depth studies of the pile-soil contact surface [42–45]. Aldaeef and Rayhani [46] introduced a roughness factor to study the characteristics of the pile-soil interface in permafrost and found that the residual strength of the interface was primarily due to residual interfacial friction and that the roughness factor decreased with decreasing temperature. Zhang et al. [47] employed direct shear tests to determine the mechanical properties of the pile-soil interface in clay soils. They found that roughness, water content, and shear rate were the main influencing factors. For the precast concrete pile-hydraulic soil interface friction problem, Zhao et al. [48] found that the interface friction capacity is highly dependent on the hydraulic soil strength. Additionally, numerical simulations have become a critical tool for investigating pile-soil contact surfaces [49–52]. For example, Wang et al. [53] developed a finite element model based on shear test data to explore the effect of thermal loads and the pile-soil interface on the thermo-mechanical behavior of piles. González et al. [54] developed a BEM-FEM equivalent linear model to analyze the horizontal load response of piles in the case of pile-soil interface degradation.

The majority of the literature on the above-mentioned research focuses on the factors influencing the mechanical properties of the pile-soil interface [46–48, 53, 54], whereas the thickness of the contact surface itself receives less attention. In addition, it is difficult to determine the thickness of the contact surface, which complicates practical application. At the same time, to facilitate modeling and computation, the model and contact surfaces are typically simplified by using the software’s built-in model and setting the shear stiffness of the pile-soil interface to a constant value [50–52]. However, the contact surface shear stiffness is not constant in practice, which significantly impacts the accuracy of the pile-soil numerical simulation. Therefore, it is necessary to develop new contact surface intrinsic structure models.

In order to accurately simulate the pile-soil contact surface unit and its interaction, the effect of contact surface thickness is eliminated, which accounts for the nonlinear variation of shear stiffness at the pile-soil interface. This study develops a new contact surface constitutive relationship independent of the contact surface thickness using the statistical damage constitutive model proposed by Yang and Liu [41]. The first section of the study reviews the literature and analyzes pile-soil contact surface tests, as well as the constitutive model and the pile-soil interface. In the second section, mathematical treatment is used to get the improved contact surface constitutive equations. The following section describes the calculation of the model parameters. The fourth section programs the improved contact surface constitutive model using the FISH language in FLAC3D and applies it to pile-soil contact surface shear numerical simulations. Ultimately, the case study in Section 5 simulates the monopile bearing and analyzes the pile bearing characteristics as well as the final pile bearing capacity under a variety of working conditions.

#### 2. Improved Constitutive Model of Pile-Soil Interface

##### 2.1. Contact Surface Shear Band and Thickness

Under the extrusion of the adjacent soil, normal stress will exist between the pile surface and the soil particles, and the two will be in close contact. After the load is applied to the top of the pile, it tends to move downward, generating static friction at the pile-soil contact surface. When the pile top load gradually increases and exceeds the maximum static friction, the pile will generate downward displacement with respect to the soil; at that point, the static friction is transformed into sliding friction and shear occurs between the two. During the shearing process, shear stress is applied to the soil along the contact surface, casing, and the soil particles to be compressed, dislocated, and slipped amongst one another, before being reorganized and finally reaching a stable condition. Therefore, the soil region adjacent to the contact surface is referred to as the shear band, as illustrated schematically in Figure 1.

The contact surface thickness is difficult to measure since the shear mechanical behavior of the contact surface is affected by numerous factors, including normal stress roughness soil properties. On the relatively smooth surface, shear stress is minimal, and the range of soil influence is small, resulting in a thin contact surface. However, for rough surfaces, the contact surface thickness is typically greater. The schematic diagram is shown in Figure 2.

**(a)**

**(b)**

##### 2.2. Improved Contact Surface Constitutive Model

Yang and Liu [41] developed a statistical damage constitutive model for the soil contact surface based on the continuous strength and statistical damage theory, whose expression is as follows:where *τ* is the shear stress at the pile-soil contact, *G* is the shear modulus, *γ* is the shear strain, and *m* and *F* are the fitted parameters. Further improvement of this constitutive model is achieved by considering the relatively small thickness of the pile-soil contact surface and assuming that the shear stress and shear strain are uniformly distributed along the contact surface and the shear strain is linearly related to the shear displacement, that is, *γ* *=* *λ*Δ, where *λ* *=* 1/*t* and *t* is the thickness of the contact surface, allowing equation (1) to be expressed as follows:

Taking the logarithm of both sides of the above equation, we obtain

Then, considering the displacement Δ of the derivative of equation (2), we can get the shear stiffness *k*_{s} and shear displacement of the relationship as follows:

Let the displacement Δ = 0 and the initial shear stiffness be

Then, *Y* in equation (3) can be expressed as follows:

Therefore, all that is required is to find the initial stiffness *k*_{si} and to combine it with the points on the contact surface shear test-curve to determine *C* and *D* in equation (3). The parameters *m* and *F* can be expressed in terms of *C* and *D*, respectively:

Substituting equations (5) and (7) into equations (2) and (4) yields

From equations (8) and (9), it can be seen that *k*_{si} will have a large effect on the fitting effect of the shear process, so in order to ensure the validity of the fitting results, an appropriate solution must be chosen to calculate the initial shear stiffness *k*_{si}. Considering the presence of slip between the pile and soil and the fact that the shear stiffness at the pile-soil interface does not vary uniformly, a solution approach proposed by Alonso [55] was borrowed. As shown in Figure 3, in the *τ*−Δ curve, the slope of the initial tangent line is the initial shear stiffness: , and in order to calculate , the parameter is introduced to establish the link between and so that the relationship between the initial shear stiffness *k*_{si} and the ultimate shear stress at the pile-soil interface can be established . In the above equation, the ultimate shear stress and the ultimate shear displacement can be derived from the experimental shear curve.

The proposed procedure successfully obtains the initial pile-soil shear stiffness *k*_{si}, which in turn results in an improved constitutive model of the pile-soil contact surface and the shear stiffness *k*_{s} during shear. This method produces *k*_{si} directly by deriving the displacement Δ in equation (2) and calculating its value without regard for the contact surface thickness. This means that the contact surface constitutive model established in this section is not directly related to the contact surface thickness, and because the contact surface thickness has no direct effect on the shear stress and shear stiffness, this study referred to it as the “zero-thickness element.”

According to equation (9), shear stiffness varies continuously with shear displacement, which can be utilized to characterize the nonlinear fluctuation of contact surface shear stiffness, and is thus will more accurate. In addition, the shear stiffness-shear displacement relationship in equation (9) is programmatically calculable, which forms the basis for the numerical computation of the monopile bearing characteristics that will be performed in this study.

#### 3. Methodology for Calculating Model Parameters

##### 3.1. Pile-Soil Contact Surface Shear Test

To evaluate the accuracy of the improved pile-soil contact surface intrinsic model presented in the previous subsection, this section uses the data from Yang et al.’s contact surface shear test [56], which includes the following pertinent information:

The test apparatus was modified based on a powered single shear apparatus. The lower box of the shear apparatus was a precast concrete slab to simulate the concrete pile, and the shear stress was directly applied through the right weight. The contact circle diameter between the soil sample and the concrete slab was raised from 61.8 to 88 mm, and the shear area was doubled, reducing stress concentration at the edge of the soil sample and the effect of oblique shear. Figure 4 illustrates the schematic diagram of shear test.

Soil samples were taken from yellow clay excavated from the foundation pit of a new university building. Prior to the test, the soil samples were air dried and then sieved to a fineness of 0.5 mm. The redesigned soil samples were then taken according to the geotechnical test procedure, and the dry density was maintained at 1.63 g/cm^{3} for 1–2 days by stratified compaction. The loading control standard is as follows: (1) the shear displacement is stable under the action of this stage load; (2) the displacement is not stable, but the load of this stage has been applied for 2 minutes and the next load is applied at this time; (3) if the shear displacement continues to develop under the current load, the soil sample is regarded as destroyed. The contact surface shear stress-shear displacement curves for various normal stresses are shown in Figure 5.

##### 3.2. Calculation of Model Parameters under Different Normal Stresses

First, the initial shear stiffness *k*_{si} is determined from , where Δ_{u} is not greatly affected by the normal stress and is uniformly taken to be 2.5 mm. To ensure the accuracy of the simulation results, the error analysis coefficient *R*^{2} is introduced to determine the parameter *χ*, and the expression is as follows:where *y _{i}* is the experimental stress value,

*f*is the simulated stress value, is the average value of the experimental stress. By fitting the experimental data from the literature [56], selecting suitable parameters

_{i}*χ*to make the fitting results closer to the experimental data, and then deriving the initial shear stiffness

*k*

_{si}at different normal stresses, the results are obtained, as listed in Table 1.

Then, the data *k*_{si}, *τ*, and Δ are substituted into equation (3) to obtain the model parameters under different normal stress and the results are listed in Table 2.

Tables 1 and 2 demonstrate that as the normal stress increases, the initial shear stiffness *k*_{si} of the contact surface increases; additionally, different pile depths correspond to different normal stresses, and changes in normal stresses lead to changes in the values taken for the parameters, indicating that the constitutive model has a depth effect that accurately describes the nonlinear variation in contact surface shear stiffness.

After obtaining the values of the parameters at the four nodes, the next step is to request the values of the parameters at remaining normal stresses. The normal stress range of 22–107 kPa is divided into three intervals: 22–55 kPa, 55–88 kPa, and 88–107 kPa, with the respective intervals of 33 kPa, 33 kPa, and 19 kPa. Because the range of each interval has a minor difference, we may obtain the relevant model parameter values under various normal stress conditions by performing linear interpolation computation on the node parameters.

To prove the effectiveness of the linear interpolation method, this section examines the normal stress of 42 kPa that is included in the range 22–55 kPa, and the corresponding model parameter values are listed in Table 3. Also, the simulation curve for 42 kPa normal stress may be obtained, and the comparison with the shear test curve obtained in reference [56] is displayed in Figure 6.

Figure 6 reveals that the shear stress-shear displacement curves derived from the calculation of the model parameters using the linear interpolation approach fit very well and meet the accuracy requirements, indicating that this method is feasible for calculating the model parameter values.

##### 3.3. Influence of the Range of Normal Stress Interpolation Intervals on the Parameters

The interpolation interval range of the linear interpolation method may affect the accuracy of interpolation results. Interpolation calculations are performed on the model parameters corresponding to the normal stress 42 kPa, in the interval ranges 22–55 kPa, 22–88 kPa, and 22–107 kPa, and values are 85 kPa, 66 kPa, and 33 kPa, respectively. The model parameters obtained using different interpolation regions for the normal stress 42 kPa are summarized in Table 4.

From the data in Table 4, the corresponding simulation curve can be obtained. The comparison of the simulation and the test curves at 42 kPa normal stress calculated in various interpolation areas is illustrated in Figure 7.

As can be seen from Figure 7, the accuracy of the simulation curve varies with the interpolation interval ranges; from high to low, the accuracy of the simulation curve is 33 kPa, 66 kPa, and 85 kPa, and its accuracy falls as the interpolation area increases. To ensure the accuracy of the model parameters, it is recommended that the interpolation interval should be kept below 40 kPa in order to obtain a reasonable simulation curve depicting the mechanical behavior of the contact surface.

#### 4. Numerical Realization of Improved Constitutive Model

##### 4.1. Numerical Calculation Process Design

Due to the ideal elastoplastic model embedded in FLAC3D having a fixed pile-soil contact surface stiffness, it cannot accurately reflect the nonlinear characteristics. Considering the characteristics of variable shear stiffness and the depth effect of contact surface as expressed in equation (9), and in order to adapt to the FLAC3D solution method and apply it more accurately to FLAC3D, the shear stiffness adjustment factor is added to equation (9) to obtain the following equation:where *k*_{s,j} is the contact surface shear stiffness of the *j*th node of the contact surface, Δ_{j} is the shear displacement of the *j*th node of the contact surface, *k*_{si,j} is the initial shear stiffness of the *j*th node of the contact surface, *C*_{j} and *D*_{j} are model parameter value of the *j*th node of the contact surface, and *φ* is the adjustment coefficient for shear stiffness.

The contact surface fish function embedded in FLAC3D may be used to acquire the normal stress value corresponding to the contact surface nodes, and the model parameters under the normal stress can be derived using the approach mentioned above for calculating model parameters. The initial contact surface shear stiffness can be determined by substituting the model parameter values into equation (9), which is then corrected and the contact surface shear stiffness parameters assigned. The specific implementation and solution process of the shear stiffness adjustment factor are depicted in Figure 8.

##### 4.2. Model Building and Parameter Setting

On the basis of the foregoing, the numerical calculation process was successfully applied using an improved pile-soil contact surface principal structure model. Due to the symmetry of the pile structure and applied loads, a half-pile-soil FLAC3D model is developed here. The diameter of the pile is 0.6 m, and the length is 5.1 m in this model (exposed ground 0.1 m). The soil specifications are as follows: 5 m in thickness, 16 m in length, and 8 m in width; grid division: 2196 units and 3000 nodes (not excluding empty units and nodes). The model is shown in Figure 9.

The numerical calculation model boundary cases are as follows: the top surface is free, the bottom surface is constrained in the *Z* direction only (except for solid pile), the *X* direction is constrained in the displacements at the two boundary surfaces at *X* = −8 m and *X* = 8 m, and the *Y* direction is constrained in the displacements at the two boundary surfaces at *Y* = 0 m and *Y* = 8 m. The pile was modeled as an elastic model with a bulk modulus of 13.9 GPa and a shear modulus of 10.4 GPa, while the surrounding soil was modeled as a Mohr–Coulomb model. The values of the parameters are listed in Table 5.

According to the shear test results [56], the internal friction angle of the pile-soil contact surface was set at 29.4° and the cohesive force was set to 14.8 kPa. It is sufficient that the normal stiffness is not negligible, and twice the maximum normal stiffness of the surrounding unit body was considered to be 292.9 MPa/m based on the numerical calculation experience. Since there is a constant variation in shear stiffness, the contact surface shear stiffness in the interval 22–107 kPa can be computed using the fitted data from Section 3 and the interpolation method. The following solution is proposed for shear stiffnesses in the 0–22 kPa normal stress interval; it is assumed that when the normal stress value is small, the shear stiffness in the 0–22 kPa normal stress interval varies similarly to the 22–55 kPa normal stress range, such that the initial shear stiffness in the 0–22 kPa normal stress interval can be obtained via linear extrapolation.

##### 4.3. Implementation of Numerical Simulation and Comparison with Theoretical Calculation Results

Vertical displacement of the pile can be achieved through the command flow in the software by applying vertical displacement to the pile unit to simulate pile-soil shear, and the greater the vertical displacement, the greater the shear displacement between the pile and the soil. The track of the shear displacement and shear stress is kept during the shear process in order to plot the shear curve. In this study, the ratio of the maximum imbalance force to the typical internal force is less than 10^{−6} as the convergence criterion of the calculation, and the calculation is considered complete when the criterion is met.

However, the range of values for the normal stress is rather scattered, resulting in insufficient data to be used to determine constant normal stress. Therefore, to simulate the shear stress-shear displacement curves of the nodes at the contact surface of the pile-soil interface under constant normal stress, the normal stress increments are applied to each node using fish language programming to ensure that the nodal normal stress is a certain value, ensuring that the majority of the nodal normal stresses are concentrated around the simulated normal stress value, thereby satisfying the requirement for data acquisition. Four different normal stresses of 16 kPa, 22 kPa, 25 kPa, and 27 kPa are selected below, and the associated shear stress-shear displacement curves are obtained through simulation using the shear stiffness adjustment factors for the four cases, as listed in Table 6. The shear stress-shear displacement curves obtained numerically were then compared to those derived using the improved contact surface constitutive structure model, which is displayed in Figure 10.

In this section, the improved pile-soil contact surface was put to the numerical calculation software FLAC3D, and the contact surface shear calculation was successfully completed, yielding a numerical simulation graph of shear stress-shear displacement. As can be seen from the comparison graph of the results, the numerical and theoretical results are quite similar, and the fitting degree is quite high. Additionally, when the normal stress is in the range of 0–22 kPa, the linear extrapolate approach can be used to solve the shear stiffness more accurately. Using equation (10) to calculate the correlation coefficient of the two, when the normal stress is 16 kPa, 22 kPa, 25 kPa, and 27 kPa, the correlation coefficients are 0.9832, 0.9826, 0.9844, and 0.9851, respectively, all of which are close to 1. Therefore, the numerical calculation can be regarded as precise.

#### 5. Case Study

##### 5.1. Numerical Calculation of Ultimate Bearing Capacity of Single Pile

In this section, the simulation and analysis of monopile load bearing will be introduced with the pile-soil shear model. Considering that the shear stress-shear displacement curves for the selected pile-soil indoor shear tests exhibit lateral resistance hardening characteristics, this section will further optimize the implementation path of the model reflecting lateral resistance hardening in FLAC3D. As shown in equation (9), when Δ = (*D*exp*C*)^{−1/D}, the shear stiffness *k*_{s} = 0 and the shear stress is at its maximum value, it can be assumed that when the shear stiffness is less than or equal to 0, the shear stress is maintained constant to reflect the lateral resistance hardening characteristics of the pile-soil contact surface during the shear process, and the above idea can be implemented in FLAC3D through fish language programming. According to the data in Table 6, the shear stiffness adjustment factor can be assumed as 1 during the numerical calculation process. In the line of the above ideas, the calculation program of the previous section for the pile-soil shear process can be modified to perform the numerical calculation for single pile bearing.

Due to the symmetry of the pile foundation structure and the added load, only half of the model needs to be built. The model specifies the following pile parameters: diameter 0.6 m and length 5.1 m (exposed ground 0.1 m). The soil specifications are as follows: 8 m in soil thickness, 16 m in length, 8 m in width, and 3492 cells and 4205 nodes of the grid. The model is depicted in Figure 11. The numerically calculated boundary cases of the model are as follows: the top surface is the free surface; in the *X* direction, the displacements of the two boundary surfaces at *X* = −8 m and *X* = 8 m are constrained; in the *Y* direction, the displacements of the two boundary surfaces at *Y* = 0 m and *Y* = 8 m are constrained (except for the solid pile); in the *Z* direction, the surface at *Z* = −8 m is constrained. The pile and soil model parameters are the same as in the numerical pile-soil shear simulation approach described previously.

The pile-soil interface constitutive model is based on the lateral resistance hardening model given previously, with the remaining parameters identical to those in Figure 9. To execute the numerical simulation of the bearing capacity of a single pile, the pile top is loaded incrementally at a rate of 40 kPa. The ultimate bearing capacity of a single pile is determined by observing the dramatic shift of pile foundation settlement, and the load-settlement curve is seen in Figure 12.

As displayed in Figure 12, when the load applied to the top of the pile is minimal, the load-settlement curve changes essentially linearly. However, when the applied load reaches a critical value, the settlement displacement of the top of the pile reduces steeply as the load increases. The settlement curve presented in Figure 8 follows a pattern similar to that observed during the static load test of the monopile, which has a bearing capacity of approximately 760 kPa under such conditions.

##### 5.2. Applied Research on Pile Foundation Bearing Characteristics

###### 5.2.1. Numerical Analysis of Axial Stresses in Piles

Figure 13 describes the variation of axial stresses within the pile body as the function of pile top load. As can be seen from the figure, the pile is subjected to vertical loads of 120 kPa, 240 kPa, 360 kPa, 480 kPa, 600 kPa, and 720 kPa, respectively. As the load increases, the axial stress in the pile body increases accordingly, resulting in a compressed stress state throughout the pile foundation. Within the ultimate bearing capacity of the pile foundation, the axial stress at the same position on the pile body grows essentially linearly with the vertical load.

Axial pile stresses exhibit the same variation rule under different vertical loads, which can be summarized as follows: as the depth of the pile foundation increases, the axial pile stresses gradually decrease and the pile end resistance values drop, indicating that the pile side frictional resistance carries the majority of the pile top load. The axial stress of pile diminishes nonlinearly with depth at all levels of load. The slope of the curve is greater in the upper soil part, revealing a faster rate of change, whereas the slope of the curve declines in the lower soil part, revealing a slower rate of change. This suggests that the axial force of the upper part of the pile drops rapidly, while the lower part slows down and tends to remain stable.

###### 5.2.2. Numerical Analysis of Pile Side Friction Resistance

Figure 14 depicts the variation curves for pile lateral friction resistance under various pile top loads. As the pile top load increases, the pile side frictional resistance increases proportionately, and the pile side frictional resistance caused by the soil at the upper part of the pile perimeter is greater. During the loading process, the upper part of the pile is compressed, resulting in downward displacement relative to the soil, while the pile side is also subjected to the upward frictional resistance of the soil; thus, the pile top load is transferred to the soil around the pile through the frictional resistance, resulting in decreasing pile axial stress with depth.

When the load value on the pile top is increased, the load is transferred from top to bottom, increasing the compression and displacement of the pile body. The relative displacement between the pile and the lower soil occurs, and the friction of the lower soil layer of the pile body is gradually exerted. This demonstrates that the lateral friction of the upper and lower soil layers of the pile is not synchronized, but that the friction of the upper soil layer of the pile plays a role first. If the load continues to increase, the resistance at the pile end will begin to play out until it reaches the bearing limit of the bearing layer, and the pile top displacement will increase significantly, resulting in damage.

From the above analysis, it is clear that the simulation results for the variation curves of axial stress and lateral frictional resistance of the pile under vertical load are consistent with the conventional pile foundation load transfer law, which more accurately simulates the top-down pile top load transfer.

##### 5.3. Theoretical Calculation of Monopile Bearing Capacity

There are numerous requirements and methods in various codes for calculating the vertical bearing capacity of a single pile; but in general, the bearing capacity of the pile is divided into two parts: pile side friction and pile end resistance. The pile side friction is determined by the different rock layers traversed by the pile, and different coefficients are chosen, but the method of calculation is based on friction theory, whereas the pile end resistance is primarily determined by the lithology of the bearing layer and the load transfer mechanism.

According to the Chinese code—“Code for Design on Subsoil and Foundation of Railway Bridge and Culvert” [57], the axial compressive bearing capacity of a single pile can be computed using the following equation:where U is the perimeter of the pile section, *f*_{i} is the ultimate frictional resistance of the soil layer on the side of the pile, *l*_{i} is the thickness of the soil layer on the side of the pile, *m*_{0} is the discount factor for the support force at the base of the pile, *A* is the area of the base of the pile, and [*σ*] is the allowable bearing capacity of the foundation soil at the base of the pile. Based on the basic parameters of the pile foundation model in this section and the provisions of the code on the relevant coefficients, the above parameters are obtained, as listed in Table 7.

By substituting the coefficients from Table 7 into equation (12), the theoretical value of the ultimate bearing capacity of the single pile can be obtained. After calculation, the final bearing capacity is about 233.36 kN, which is then divided by the pile area to get around 825 kPa. In comparison to the numerical calculation result of 760 kPa, the theoretical calculation value is 825 kPa, and the discrepancy between the two values is stable at 10%. This reveals that the theoretical calculation results provide a more robust verification of the numerical simulation’s correctness.

##### 5.4. Analysis of Influencing Factors of Ultimate Bearing Capacity of Single Pile

To facilitate comprehension of the application of the improved contact surface constitutive model and to further verify the validity of the model, this subsection will simulate the ultimate bearing capacity of a single pile with varying pile lengths and soil moduli.

###### 5.4.1. Effect of Different Pile Lengths

In Section 5.1, the length of the pile in the model is 5 m in length and 0.6 m in diameter. Keeping the pile diameter constant and changing the length of the pile, models with pile lengths of 7.5 m, 10 m, and 12.5 m, were established to investigate the effect of varying pile lengths on the ultimate bearing capacity of the monopile. To minimize the influence of boundary effects, the model height is adjusted to double the length of the pile, and the remaining material characteristics and boundary conditions remain unchanged. The model loading scheme is still loaded incrementally with 40 kPa, and the settlement of pile tops is recorded for each load step. Figure 15 illustrates the load-settlement curves for various pile lengths. The inflection point of the curve represents the ultimate bearing capacity of the single pile.

As shown in Figure 15, the pile top settlement increases linearly with increasing load at the commencement of loading. After the inflection point of the curve, the settlement accumulated rapidly and the pile foundation is damaged. The longer the pile length is, the less settlement there will be under the same load. At the same time, the ultimate bearing capacity of a pile foundation grows with the pile length increases. This is because as the pile length increases, the pile side frictional resistance increases bearing capacity, and the pile axial force is transferred more to the pile perimeter soil, which reduces pile body compression and pile bottom soil compression. Therefore, increasing the pile length can reduce settlement and improve the pile bearing capacity.

###### 5.4.2. Effect of Different Soil Moduli

The pile length of 5 m and the pile diameter of 0.6 m are maintained, but the soil modulus around the pile is altered. To test the effect of different soil deformation moduli on the ultimate bearing capacity of the monopile, the soil deformation modulus was set to 10 MPa, 16.67 MPa, and 20 MPa. The remaining material characteristics and boundary conditions remain unchanged, as does the loading method. The load-settlement curves for the pile with varying soil moduli around the pile are plotted in Figure 16.

As depicted in Figure 16, the variation pattern of pile Q-S curves with varying soil moduli is similar to that of Q-S curves with varying pile lengths, where settlement initially increases linearly and then rapidly increases once the inflection point appears. Increases in the soil modulus surrounding the pile result in a decrease in the settlement at the top of the pile and an increase in ultimate bearing capacity. This corresponds to the impact of lengthening the pile. When the soil modulus around the pile increases, the pile lateral frictional resistance increases as well. Because the soil at the bottom of the pile gets stronger, its bearing capacity increases and the compression deformation decreases, and increasing the soil modulus can reduce settlement and enhance the pile bearing capacity.

As can be observed from the above analysis, the simulation results for the ultimate bearing capacity of a single pile under various pile lengths and soil moduli are consistent with the actual pile bearing capacity variation law, which more accurately simulates the pile bearing capacity variation law.

This section applies an improved constitutive model of the pile-soil interface to the simulation of single pile bearing capacity, verifies the ultimate bearing capacity of a single pile, analyzes the variation law of pile axial force and pile side friction, and investigates the influence of pile length and soil modulus around the pile on the ultimate bearing capacity of a single pile. The case study proves the effectiveness of the improved constitutive model, which has a certain engineering application value.

#### 6. Conclusion

The constitutive model of the pile-soil interface with variable shear stiffness and zero thickness element is improved using the original contact surface model. At the same time, the numerical simulation of the new constitutive model is realized, as well as a case study of a single pile is carried out. The following conclusions are drawn:(1)The improved constitutive model is capable of describing the nonlinearity of the contact surface shear stiffness, and it is beneficial to program computation and numerical simulation research, laying the foundation for further contact surface and single pile bearing numerical calculations.(2)Combined with the data from the shear test, the parameter calculation method is proposed. The parameter values for various normal stresses can be determined by linear interpolation, which reflects the depth effect of the model. Because the accuracy of the model parameters reduces as the interpolation area increases, it is recommended that the interpolation interval be kept small.(3)Implementing an improved constitutive model of the pile-soil interface in FLAC 3D is proposed, as is the usage of a linear extrapolation method for the shear stiffness in the normal stress range of 0–22 kPa. The model and calculation procedure successfully simulate pile-soil shear, and the numerical results are in good agreement with the theoretical results, which accurately depict the nonlinear stiffness of the contact surface.(4)The improved constitutive model is used in the bearing simulation of a single pile. The case study demonstrates that during pile top load transfer, the upper soil part of the pile provides more pile lateral friction resistance than the lower soil part because the pile and upper soil body are the first to be displaced relative to one another. By increasing the pile length and the soil modulus around the pile, settlement can be reduced and the pile bearing capacity increased. The feasibility of the model is verified by comparing the calculated results to the equations in the Chinese code, which serves as a guide for engineering practice.

#### Data Availability

The data used to support the findings of this study are included in the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors would like to thank “SCI Capability Sharing United Group” for the English editing during the preparation of this manuscript. This work was supported by the Natural Science Foundation of China (grant nos. 51308552 and 51778633).