#### Abstract

Weaknesses of the subgrade structure induce the asphalt surface diseases and shorten the service life of flexible pavement. However, the resilient modulus (Mr) of subgrade soils is difficult to be evaluated directly since the subgrade is hidden and covered by the granular or asphalt layer. This study aimed to establish a numerical approach to predict the dynamic behavior of flexible pavements considering the stress sensitivity and moisture variation of subgrade soils. Firstly, 2D FEM simulations of flexible pavements were performed with half-sine loadings. A constitutive model of subgrade soils was proposed to incorporate soil suction and octahedral shear stress. It was validated using the laboratory triaxial test data of 3 selected soils. Then, the developed model was programmed by the user-defined material subroutine (UMAT) in the software ABAQUS. Subsequently, the validity of FEM model was verified by the laboratory tank model. Finally, the effect of moisture contents on the dynamic response of pavement structures was studied by tensile stress and vertical compressive strain. Results show that the surface deflection of the FEM model is similar to that of the actual pavement structure with the *R*^{2} of 98.44%. The developed UMAT program is reliable since the distribution of Mr in the FEM model is influenced by the stress and moisture condition of subgrade soils. When the moisture content is increased by 63%, the average Mr of subgrade soils is decreased by 18.7%. Meanwhile, the stiffness softening of subgrade soils increases vertical compressive strain at the top of the subgrade and the tensile stress at the bottom of the surface layer. It is interesting that the developed model can be applied to analyze the fatigue cracking of both subgrade and surface layers in the future.

#### 1. Introduction

Flexible pavements are widely used in the world and in general consist of three layers: subgrade soils, base course, and asphalt surface [1]. The asphalt surface layer provides the comfortable driving environment, and the subgrade layer ensures the structural deformation stability. As time goes on, as a result of long-term exposure to the external environment, the asphalt surface weakens gradually due to adverse temperature and traffic loads. Also, a large number of studies have been conducted to analyze the mechanism and causes of the asphalt surface damage [2, 3]. However, dynamic properties and permanent deformation of the subgrade structure are difficult to be directly measured during the service period, since subgrade soils are hidden and covered by granular or asphalt base [4]. Few studies have been proposed to describe the dynamic behavior of the subgrade structure. Meanwhile, the entire section of the subgrade is seriously wetted when the humidity of subgrade soils tends to balance. The balanced moisture content of the subgrade is up to 71% higher than the maximum moisture content of the initial design in the south China area [5]. The stiffness of the subgrade decreases significantly and the deformation increases continuously due to the effect of humidity and cumulative load [6]. These weaknesses of the subgrade structure become major reasons for the failure of flexible pavement, except for asphalt materials ageing.

Subgrade soils distribute nonuniform dynamic stresses caused by the traffic loading, and both moisture content and matric suction change continuously from top to bottom [7]. Most flexible pavements are usually static designed by treating the subgrade layer as a linear-elastic material with constant mechanical parameters. However, the dynamic response of flexible pavements is different from that of the static design, since the tensile strain at the bottom of the asphalt surface is 50% less than that of the static response while the shear stress of base course and asphalt surface is significantly increased [8]. Furthermore, several research studies have concluded that the soil subgrade exhibits nonlinear stress-dependent behavior related to soil physical properties [9]. Recognizing these deficiencies, the Mechanistic Empirical Pavement Design Guide (MEPDG) 2004 recommended the resilient modulus (Mr) of subgrade soils as a basic input to calculate the dynamic response of flexible pavements [10].

The Mr describes the stress-dependent elastic behavior of subgrade soils and is defined as the ratio of stress to strain in cyclic loading under the elastic behavior phase of granular materials [11]. Variations of the Mr often manifest as multiple nonlinearities with the change of test conditions in repeated triaxial tests. The laboratory test is an accurate approach to evaluate the Mr of soils in different conditions of stresses and moisture [12]. In repeated load triaxial tests, the confining stresses are applied to simulate the overburden pressure and wheel load but it requires well-trained professional and relatively expensive laboratory device. Considering the complexity and inconvenience of the laboratory testing, many researchers attempted to build Mr-predicting models [13]. They found the Mr of subgrade soils is obviously dependent on the stress and moisture status [14]. At the same time, these physical characteristics of soils are easy to be obtained and embedded in the predicting models. Therefore, it is highly appropriate to employ a Mr-predicting model to characterize mechanical behaviors of subgrade soils relating to physical properties.

The numerical methods are regarded as a beneficial supplement of the mechanistic empirical method [15]. Two constitutive models were always used to describe the elastic behavior of soils, such as the linear elastic model and nonlinear elastic model [16]. Based on the linear elastic model, the pavement structure is considered to be a multilayer elastic system, which provides a good analytic approach to explain its static mechanical behavior [17]. However, the model is not suitable for the pavement design under high-speed heavy-duty traffic conditions. Under dynamic loading, the theoretical response is smaller than the realistic measure, since the nonlinear material properties of subgrade soils are overlooked [18]. Recently, the finite element method (FEM) has been used widely to provide better solutions for the pavement dynamic analysis [19]. More specifically, the nonlinear behaviors of subgrade soils can be simulated through a user-defined material subroutine (UMAT), which is programmed according to the user-defined modulus and stress iterative process [20]. In addition, the modulus of a nonlinear UMAT can be calculated using a Mr-predicting model developed by repeated load triaxial tests. Therefore, to further understand the actual response of subgrade soils under traffic loads, a FEM model with the UMAT of subgrade soils should be established through a Mr-predicting model that is based on easily obtained parameters and that could adequately quantify the effect of the stress and moisture on the Mr of soils.

The objective of this study is to establish a numerical approach to predict the dynamic behavior of flexible pavements considering the stress sensitivity and moisture variation of subgrade soils. Specifically, a constitutive model of subgrade soils considering moisture variation was developed by the UMAT program. The predicting result for constitutive model was fitted and verified with laboratory results of repeated load triaxial tests. Then, the difference of surface deflection between FEM model and laboratory tank model was also discussed. Subsequently, the stress sensitivity of subgrade soils was analyzed by calculating the distribution of Mr and stress state in the subgrade layer. Finally, the effect of moisture contents on the dynamic response of pavement structures was studied by tensile stress and vertical compressive strain.

#### 2. FEM Simulation for Dynamic Response of Flexible Pavement

##### 2.1. FEM Model of Flexible Pavement

2D FEM simulations were performed on a flexible pavement under dynamic loadings. In general, the size of the finite domain should be greater than 20 times the width of the loading range (0.15 m) [21]. The axisymmetric model is used to simplify the flexible pavement structure. Preliminary trials indicated that the road vertical deformation error was less than 0.45, if the calculation depth was more than 3 m. Considering the validity and accuracy of the computational resources, the size of the finite element model is set to be 5 m × 3 m (in width and height), as shown in Figure 1(a).

**(a)**

**(b)**

The pavement structure consists of asphalt mixture layer, unsaturated granular base course, and soil subgrade from top to bottom [22]. Meanwhile, the interfaces between each layer are assumed to be fully bonded. The FEM model was meshed with different section and material properties. The size and mesh of pavement model are shown in Figure 1(b). The model mesh and nodes are defined as a series of 8-node biquadratic axisymmetric elements. The element mesh is first divided by global coefficients, and then the mesh density is increased in the load action area to ensure that the dynamic response calculation is accurate and precise. The subsidence of the simulated surface of the road was caused by the vehicle’s axial load, without considering the influence of self-gravity. The bottom and the side of the model were applied to the constraints, respectively.

##### 2.2. Material Parameters and Models

All the input parameters of pavement materials are displayed in Table 1. The mechanical behavior of the asphalt layer and the base layer (unsaturated granular material) is characterized by the linear elastic model [23]. Linear elastic moduli of all mentioned materials were derived from design parameters by the Specifications for Design of Highway Asphalt Pavement in China [24] in order to ensure that the FEM model is typical and accurate. In addition, the linear elastic modulus of subgrade soils is usually applied to response analysis of flexible pavements as a constant parameter. However, this research noticed that the influence of humidity and stress on subgrade soils is important, but its behavior is not described clearly by the linear elastic model. Therefore, the constitutive model of subgrade soils needs to be developed by the UMAT program.

##### 2.3. Load History of Traffic Loading

The pavement structure is subjected to dynamic loading cycles, which is always assumed to be uniformly distributed on the contact area. The position of the load is shown in Figure 1(a). To simulate the double circular loading, the loading area is simplified to be a circular rigid plate with a radius of 15 cm [18]. For the BZZ-100 standard axle load in China, axle load is 100 kN and contact pressure is 0.7 MPa. The dynamic loading sequence uses a half-sinusoidal loading waveform to simulate the traffic loads acting on the road surface. Both the load time and the unloading time are 0.1 s. The rest period for each loading cycle is 0.8 s.

**(a)**

**(b)**

**(c)**

#### 3. Development of Constitutive Model of Subgrade Soils

##### 3.1. Strain-Stress Relationship of Materials

The relationship between the stress matrix and the corresponding strain matrix is still a function, although the elastic behavior of subgrade soils is nonlinear. The stress-strain relationship of elastic materials on an axisymmetric plane can be described by generalized Hooke’s law, as shown in the following equation:where *E* is the Mr, *G* is the shear modulus, and *μ* is Poisson’s ratio.

In ABAQUS, the Jacobian matrix of the UMAT program is defined by elastic increment matrix, which is used to solve the unknown stress increments with the stress and displacement increments at the beginning of an analysis step. The constitutive relation in equation (1) needs to transform the position of the elastic matrix, as shown in equation (2), to obtain a new corresponding matrix. This Jacobian matrix simplifies the complex mechanical behavior of materials to the iterative approximation of the Mr with time steps. Therefore, the stress-strain relationship of subgrade soils will be equivalent to the nonlinear characterization of the Mr with stress and moisture conditions, if Poisson’s ratio remains constant.

##### 3.2. Material Model considering Stress State

The strain-stress relationship, especially the Mr, is a key parameter to compute the incremental elastic stress for a given node displacement. Over the years, many material models have been established to obtain the Mr of granular materials at different stress levels [25]. These generalized models were all validated by fitting well the experimental data from the repeated load triaxial tests under various levels of confining pressures. It is favorable evidence that the mechanical behavior of subgrade soils is nonlinear and stress-dependent. The *K*-*θ* model has been well documented as a commonly used nonlinear model of soils (equation (3)) [26]. A three-parameter model evolved from the *K*-*θ* model, since the octahedral shear stress is another factor that affects the triaxial test results (equation (4)) [27]. Then, the MEPDG model was adopted widely since it is consistent with the actual stress situation by bringing in the shear stress (equation (5)) [7, 10].where *θ* is the bulk stress; is the octahedral shear stress; is the atmospheric pressure, assumed to be 100 kPa; and *k*_{1}, *k*_{2}, and *k*_{3} are the regression coefficients.

##### 3.3. Developed Model considering Moisture Variation

In the MEPDG model, the confining effect is characterized by the bulk stress, which is calculated by adding together confining pressure and deviator stress. However, results of triaxial tests in different stress levels show that the confining pressure is positively correlated with the Mr, while the deviator stress is negatively correlated with the Mr, reported by Yao et al. [5]. As a result, with the variation of the bulk stress, the evolvement rule of the Mr is not clear.

To avoid no convergence when the deviator stress is dynamic, the minimum bulk stress was proposed as the factor in characterizing the confining effect and defined as the bulk stress under no deviator stress [5, 25]. Meanwhile, to integrate the moisture-dependent factor, a developed predicted model of the Mr is displayed in equation (6), which incorporates the soil suction and minimum bulk stress into the MEPDG model.where *ψ* is the soil suction; is the minimum bulk stress; and *k*_{0}, *k*_{1}, *k*_{2}, and *k*_{3} are the regression coefficients.

This equation contains three stress terms, which reflect the moisture state, the confinement effect, and the shear effect, to consider the influencing factors of the Mr [25]. In order to validate the developed model, the predicted values of equation (3) and the measured values of the three selected soils are shown in Figure 2. Among them, the soil suction of the same soil is the same, and the soil suction of different soils is different. With the change of stress level, the prediction results obtained in a certain soil suction condition agree well with the result of repeated load triaxial tests. The result shows that the material model can predict the dynamic characteristics of subgrade soils correlated with stress and moisture, since the statistical *R*^{2} between predicted and measured values is up to 90%.

The low plasticity clay from Guangdong, China, was selected to determine the Mr of subgrade soils with different compaction degrees. The optimum moisture content of low plasticity clay from Guangdong, China, is 12.5%, and the saturated moisture content is 22.02%. According to the developed model considering moisture variation (equation (5)), fitting parameters were calculated and are shown in Table 2. It is noticed that the degree of compaction in this research is slightly different from the status of subgrade soils in field construction for the convenience and distinction of the repeated load triaxial test.

#### 4. User-Defined Material Subroutine for Subgrade Soils

The moisture-dependent behavior of subgrade soils was defined by a UMAT subroutine. The flowchart of the user-defined material subroutine for subgrade soils is shown in Figure 3. The moisture of subgrade was added as a custom variable, and the different layers have different regression coefficients (*k*_{0}, *k*_{1}, *k*_{2}, and *k*_{3}). During the loading process, the distribution of the Mr is a function of the strain state, so the finite element calculation is performed using the tangential stiffness method [29].

The tangent stiffness method solves the nonlinear stress-strain behavior by updating the tangent stiffness at each iteration until the stress increment converges. When the subgrade structure is subjected to a dynamic load, the calculated stress of each unit node changes with the loading time history, and the total stress state of the model is updated by the Euler (display integral) method through the stiffness Jacobian matrix. In order to improve the convergence of the iteration, the Mr in each iteration is calculated by the following equation [18]:where is the Mr output from the *i*^{th} iteration; is the Mr output from the (*i*−1)^{th} iteration; *λ* is the damping factor (e.g., initial *λ* is 0.95); and is the computed modulus from equation (3) in the *i*^{th} iteration.

In the subgrade model, the computed stresses of element nodes are different along with different dynamic loading. The Mohr–Coulomb failure theory is applied to adjust the initially stresses of model, and the cumulative error for the entire model is controlled to be less than 0.5%. Moreover, the individual convergence criterion for each node is shown in the following equation:where is the iteration error.

#### 5. Results and Discussion

##### 5.1. Model Validation by Laboratory-Measured Deflection

The critical pavement responses were commonly described by the surface deflection. Figure 4 compares the surface deflection of the developed pavement model and the laboratory tank test when they are subjected to a 565 kPa load on a circular area with a radius of 15 cm [20]. The size of circular tank model in laboratory is 2.4 m in diameter and 1.8 m in height, including 0.15 m asphalt mixture layer, 0.25 m base course layer, and 1.40 m subgrade layer with 95% degree of compaction. Figure 4(a) illustrates the location of the sensors within the tank model, where “S1–S5” represent the linear displacement sensors.

**(a)**

**(b)**

In Figure 4(b), it is seen that as the distance to the center of the load increases, the measured and predicted road surface deflection gradually decreases. Meanwhile, the result of this model is similar to the actual measured response of pavement structure with the *R*^{2} of 98.44%. The main cause of measurement error is that the strengths of pavement materials are enhanced by the cyclohoop effect of the steel tank. In addition, the good agreements indicate the feasibility and correctness of the developed UMAT subroutine.

##### 5.2. Effect of Stress Distribution on Resilient Modulus

The distribution of Mr within the subgrade is plotted in Figure 5(a). The subgrade soils were considered as a stress-dependent material. Along with the depth, the bulk stress of soils increases and the shear stress is weakened. Under the center of dynamic loading, the Mr of the surface is the smallest, while the Mr at the bottom of the subgrade is the largest away from the surface. The distribution of the modulus is nonuniform in the vertical and horizontal directions because the internal stress state of the subgrade changes.

**(a)**

**(b)**

The vertical stress is a significant factor affecting the response of pavement structures under traffic load [30]. Figure 5(b) shows the distribution of vertical stresses within the subgrade layer. This is consistent with the state of the predicted modulus, i.e., the vertical stress near the load center changes drastically. The horizontal distribution of stress is limited, and the bulk stress of soils becomes the major factor influencing the variation of the modulus.

##### 5.3. Effect of Moisture Content on Dynamic Responses

Soil suction is an input parameter in the user-defined material model of subgrade soils (equation (5)). In order to study the effect of soil moisture content on the rebound modulus and road surface response, the influence law was analyzed by changing the soil suction of the model. In this study, three different moisture states of subgrade soils are 9.5%, 12.5%, and 15.5%, corresponding to low moisture content, optimum moisture content, and high moisture content. The soil water characteristic curve (SWCC) depicts the relationship between matric suction and water content, as shown in Figure 6. The modified Fredlund and Xing model was applied to describe the SWCC of the low plasticity clay, in equation (9) [31]. The soil suction can be back-calculated by the corresponding moisture content through the Fredlund and Xing model:where is the volume moisture content; is the saturated volume moisture content; *ψ* is the soil suction; and *a*, *n*, and *m* are the fitting parameters; the result is shown in Table 3.

Figure 7 compares the average modulus of three sublayers at different moisture conditions. The result shows that the Mr of soils increases layer by layer of the subgrade. Moisture content is an important factor affecting the Mr of subgrade, since the increase of the moisture content significantly reduces the Mr of each subgrade layer. From low moisture to high moisture, the moisture content is increased by 63%, the average Mr of subgrade soils is decreased by 18.7%.

The variations of Mr in turn change the dynamic response of the pavement construction. As shown in Figure 8, the moisture condition of subgrade soils affects the tensile stress at the bottom of the surface layer. Peaks of time history curves of tensile stress appear near 0.1 s, since the response of the asphalt surface layer is directly affected by the traffic loading in terms of waveform and period. Meanwhile, as the moisture content increases, the tensile stress of the surface layer increases from 451.46 kPa to 499.72 kPa. It indicates that the increase in the moisture content of soils will accelerate the damage of the road surface, since the change makes a difference in the stiffness softening of subgrade soils and the increase of the tensile stress at the bottom of the surface layer.

The vertical compressive strain at the top of subgrade is an important performance indicator of subgrade. As shown in Figure 9, time-history curves of vertical strain on the top surface of the subgrade at the load centerline were displayed with different moisture content. The plus sign (+) represents the tensile strain and the minus sign (−) represents the compressive strain in the FEM calculation. The compressive stress firstly increases and then decreases with the increase of the loading time. The dynamic deformation within subgrade soils is nonlinear since the curve peak is delayed by about 0.05 s and the residual strain cannot be neglected after stopping the loading. Meanwhile, as the moisture content increases, the compressive strain increases remarkably from 19.50 × 10^{−5} to 22.51 × 10^{−5}.

#### 6. Conclusions

This study proposed a numerical approach to predict the dynamic response of flexible pavements by considering the stress sensitivity and moisture variation of subgrade soils. The constitutive equation of subgrade soils was verified by repeated triaxial test results of three soils. Meanwhile, the proposed stress-strain relationship was embedded in the finite element model to characterize the nonlinear behavior of subgrade soils. The validity of the FEM model was verified by the laboratory tank model. Then, the stress sensitivity of subgrade soils was analyzed by calculating the distribution of Mr and stress state in the subgrade layer. Finally, the influence of moisture condition of subgrade soils on dynamic responses of flexible pavements was analyzed. The main contributions and conclusions are as follows:(1)The constitutive model of subgrade soils was established by transition of the Jacobian matrix and iteration of the Mr. A material model of subgrade soils was developed by incorporating the soil suction and minimum bulk stress into the MEPDG model, to take into account the moisture variation. The result shows that predicted values have a good agreement with results of repeated load triaxial tests, since *R*^{2} of the fitting result is greater than 90%.(2)A convergence principle for the UMAT program was presented by the tangent stiffness method and Mohr–Coulomb failure theory. The surface deflection of the FEM model is similar to that of the actual pavement structure with *R*^{2} of 98.44%. It was also found that the distribution of Mr is nonuniform in the vertical and horizontal directions because the internal stress state of the subgrade changes. It indicates that the FEM model embedded by the developed UMAT program is accurate and reliable.(3)The moisture content of subgrade soils is characterized by the matric suction based on the modified Fredlund and Xing model. From low moisture to high moisture, the moisture content is increased by 63% and the average Mr of subgrade soils is decreased by 18.7%. The stiffness softening of subgrade soils increases vertical compressive strain at the top of the subgrade and then increases the tensile stress at the bottom of the surface layer.(4)According the pavement design theory, these changes of design indicators with the rise of moisture content lead to the shorter serve life of flexible pavements. It is interesting that the pavement life can be accurately estimated by taking into account the stress and moisture condition of subgrade soils. The research will be carried out in the future by the developed model to analyze the fatigue cracking of both subgrade and surface layers.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors gratefully acknowledge the financial support of the National Key Research and Development Program of China (2017YFC0805307), National Natural Science Foundation of China (51478054, 51878078, 51608057, and 5181102194), Excellent Youth Foundation of Natural Science Foundation of Hunan Province (2018JJ1026), Key Project of Education Department of Hunan Province (17A008), Hunan Provincial Innovation Foundation for Postgraduate (CX2018B521), and Open Research Fund of Science and Technology Innovation Platform of State Engineering Laboratory of Highway Maintenance Technology, Changsha University of Science & Technology (kfj150103). This research was also supported by the Youth Scientific Research Foundation, Central South University of Forestry and Technology (QJ2018008B), and Open Fund of Engineering Research Center of Catastrophic Prophylaxis and Treatment of Road & Traffic Safety of Ministry of Education (Changsha University of Science & Technology) (kfj170405).