Resilient moduli of subgrade materials are important parameters in the mechanistic-empirical pavement design guide (MEPDG). The resilient moduli values are generally obtained through repeated load triaxial (RLT) tests. Since the test procedure is complex and time-consuming, extensive efforts have been made to develop predictive models based on observations from laboratory tests. Several empirical resilient modulus predictive equations have been proposed by incorporating state variables, such as confining stress, bulk stress, deviator stress and octahedral shear stress[1-2]. Some classical models, such as the k-θ model and the Uzan-Witczak model, have been widely used in numerical modeling of flexible pavements due to their simplicity[3-5]. Some modifications of these models have taken cross-anisotropy into account[6], especially in the geosynthetics reinforced layer[7].
However, the aforementioned models are stress-dependent types.Recently, the importance of environmental factors, in particular, the effect of moisture content, is emphasized in pavement performance[8]. Experimental studies demonstrated that the resilient modulus of unsaturated soil is sensitive to changes of matric suction, which directly relates to moisture change[9], and the resilient modulus-suction relationships were proposed[10-12]. Nevertheless, numerical implements in FE models of these models remains vacant.
Considering the water effect in geotechnical analyses, the coupled hydro-mechanical method is generally employed for simulating the interaction of solid skeleton with flow through unsaturated soils[13]. A simplified approach was adopted by Saad[14], which combined the effective stress principle of unsaturated soil with the modified cam-clay model. However, no sufficient experimental evident supported this simplified approach. Another approach attempted to capture the coupling effect between water flow and soil skeleton from constitutive laws. A coupled two-phase fluid flow and elastoplastic model proposed by Hu et al.[15] incorporated matric suction into the modified cam-clay model to develop a 3D yield surface for unsaturated soils. Nevertheless, the complexities of these models made their numerical implementation cumbersome and impractical for routine pavement analyses. Gu et al.[16] integrated the moisture-dependent resilient modulus constitutive model into the FE pavement model. However, it can only simulate specific situations with uniform moisture distribution, which limited its application in coupled hydro-mechanical analyses.
To address the aforementioned problems, this study aimed to integrate one suction-dependent resilient modulus constitutive model for unsaturated soils into FE models based on the code Abaqus. Specifically, a UMAT subroutine was developed to characterize suction-dependent nonlinear stress-strain behavior of unsaturated granular material. A comparison of measured and simulated triaxial test results was carried out to validate the developed model. Finally, the developed model was applied to a 3D pavement model subjected to three moisture states, and the coupled hydro-mechanical responses under FWD load were analyzed to demonstrate the detrimental effects of moisture intrusion on the resilient moduli of unsaturated soils.
1 Methodology
1.1 Suction in unsaturated soils
Soil suction is comprised of matric suction and osmotic suction. Herein,the matric suction represents attraction due to capillary and surface adsorptive forces in the solid-liquid interphase of unsaturated soils; the osmotic suction is stimulated by the presence of dissolved salts in the pore fluid. Generally, the effect of osmotic suction can be negligible. Therefore, generalized suction in unsaturated soils refers to the matric suction, which is expressed as
ψ=ua-uw
(1)
where ua and uw are the pore air and pore water pressure, respectively.
The matric suction has a vital impact on the effective stress of unsaturated soils according to the effective stress principle proposed by Bishop[17].
σ′=(σ-ua)+χ(ua-uw)
(2)
where σ′ and σ are the effective and total normal stress, respectively; ua-uw presents the matric suction (ψ); χ is the Bishop’s effective stress parameter, which is a suction/saturation-dependent factor, representing the contribution of matric suction to effective stress.
For fully saturated soils, χ is equal to 1.0; and for unsaturated soils, the values of χ ranges from 0 to 1.0. A commonly used relationship between the matric suction and χ was noted by Khalili and Khabbaz[18], which is expressed as
(3)
where ψb presents the air entry value, namely, matric suction where air begins to enter pores in the soils. In other literature, χ is assumed to be equal to the saturation Sr of the medium for simplicity[19]. Besides, when ua is a constant and equals external atmospheric pressure, it can be removed from the equation of effective stress. Consequently, the effective stress principle in Abaqus is simplified to
σ′=σ+Sr(-uw)
(4)
In this case, the matric suction can be regarded as negative pore water pressure (ψ=-uw) and integrated into the concept of pore pressure with pore water pressure uw.
The soil water characteristic curve (SWCC) depicts the relationship between matric suction and water content (gravimetric, volumetric, or degree of saturation). In other words, the value of matric suction implies the moisture content in unsaturated soils. Many efforts were made to obtain a formulation that captures the configuration of SWCC. The Fredlund and Xing[20] model is one of widely accepted models to describe SWCC with the mathematical expression,
(5)
where the subscripts s and r indicate the saturated and residual values of the water content, respectively; a, n and m are the fitting parameters; ψ denotes the matric suction; and θr is assumed to be zero.
1.2 Resilient modulusconstitutive models
Resilient modulus Mr is the ratio of the deviator stress σd to the elastic axial strain εr in the repeated load triaxial test (RLT), and it is expressed as
(6)
Various models were proposed to characterize the stress-dependency of resilient modulus. Among them, the earliest one was the so-called k-θ model[21], which indicates that the resilient modulus is an exponential function of bulk stress.
(7)
where θ presents the bulk stress, and herein, a positive stress means compression; pa presents the atmospheric pressure (normalizing factor), assumed to be 100 kPa generally.
Uzan[2] proposed a classical model including the impact of shear stress,
(8)
where τoct is the octahedral shear stress.
To avoid numerical singularity when τoct is equal to zero, MEPDG adopts the Witczak and Uzan model as[22]
(9)
Recently, two resilient modulus predictive models incorporating matric suction were developed by Liang et al.[11-12] and they are expressed as
(10)
(11)
where χL is given as χ in Eq.(3); ψ represents the matric suction; θnet denotes the net bulk stress, θnet=θ-3ua; Δuw-sat represents the excessive pore water pressure under saturated conditions (ψ=0); ψ0 denotes the initial matric suction and Δψ represents the variation of matric suction due to the pore water pressure build-up under unsaturated conditions (Δuw-sat=0). The coefficients k1,k2,k3,k4 are the regression parameters which represent the effect extent of factors on the resilient modulus.
Aforementioned models were validated by accurately predicting resilient modulus from RLT tests under specified conditions. However, few attempts were made to integrate suction-dependent resilient modulus models into coupled hydro-mechanical numerical analysis. Since the accuracy of Eq.(10) strongly depends on the correct χ estimation and no commonly adopted estimation method of χ is available, the constitutive model Eq.(11) is utilized to model the hydro-mechanical coupling behaviors of unsaturated soil in this paper.
It should be noted that the stress used in code Abaqus is the effective stress expressed in Eq.(4), which incorporates the pore-water pressure and matric suction into pore pressure (matric suction is negative). In addition, the increment of the pore pressure is calculated before stress computation. Therefore, the constitutive equation for numerical analysis in code Abaqus is adjusted as
(12)
where θ′ represents the effective bulk stress, θ′=θnet-3Srψ, incorporating the impact of the pore-water pressure.
2 Model Development
2.1 Development of suction-dependent resilient modulus model
A user-defined material (UMAT) subroutine is programmed to characterize suction-dependent nonlinear elastic behaviors that are not included in the code Abaqus material library. Since the Newton-Raphson method is used to solve the nonlinear problems in code Abaqus, a tangent stiffness matrix is required for iteration. A fully consistent tangent stiffness matrix can be obtained through the partial derivative[4], but it results in a non-symmetric matrix, lengthy mathematical expressions, and consequently, the increase of computational cost. Young’s elasticity modulus is replaced by state-dependent resilient modulus for simplicity, and gives the stress-strain relationship in tensor notation,
(13)
where σij,εij represent the stress and strain tensor, respectively; δij denotes the Kronecker delta; and v is the Poisson’s ratio. The convergence problems of this algorithm were reported by Hjelmstad et al[3]. The modified Euler method with error control was adopted as the integration algorithm of UMAT subroutine to determine the nonlinear stress solution in each iteration. The convergence criterion is
(14)
where Er is the iteration error; △σ1=D(σ0)△ε;△σ2=D(σ0+△σ1)△ε;△σaver=(△σ1+△σ2)/2;D(σ) represents the stiffness matrix, which is dependent on state variables; △ε is the strain increment for the current iteration; and the subscripts 0, 1, 2 stand for the initial, the first trial and the second trial conditions, respectively; the notation “
” is the 2-norm operator.
In order to obtain the computational outputs of the pore pressure and degree of saturation, the user-defined field variables (USDFLD) subroutine is utilized. The outputs can be stored in state variable array and transmitted to UMAT subroutine.
The flowchart of the developed UMAT subroutine is shown in Fig.1.
Fig.1 Flowchart of the developed UMAT subroutine
2.2 Model validation
To validate the developed model (Eq.(12)), the results of triaxial compression tests in consolidation-drainage (CD) condition under constant suction for expansive clay conducted by Chen[23] were used. The confining pressure of tests was 350 kPa, and the controlled suctions were 0, 200, 350 kPa, respectively.
The relationship between the deviator stress σd and axial strain εa was used to calibrate parameters k1,k2,k3 and k4 of the proposed model. The variation of deviator stress with axial strain can be divided into multiple linear elastic processes, as shown in Fig.2, in which the elastic modulus is replaced by resilient modulus Mr determined by the state of current incremental step in code Abaqus. Under the triaxial condition, each state variable involved in the proposed model can be determined for each step,
(15)
(16)
Mr=Es ,i(1+v)
(17)
where σc represents the confining stress, which is a constant (350 kPa) for the entire process; σd ,i denotes the deviator stress of initial state in the current step; v is the Poisson ratio, specified as 0.3; Es ,i=△σd ,i/△εa ,i, which is the secant slope in σd-εa curve.
Fig.2 Schematic of multiple linear elastic processes
The effectiveness of data was checked in advance, for instance, points with discontinuous slopes were excluded. Nonlinear fitting was carried out to determine the developed model parameters k1,k2,k3 and k4. Parameters for the simulation of the triaxial test are shown in Tab.1. Besides, the SWCC was fitted by the Fredlund and Xing model (Eq.(5)) with parameters a=141.43 kPa, m=1.238, and n=0.193.
Tab.1 Parameters for simulation of triaxial test
k1/kPak2k3k4Poisson s ratioSaturated permeability/(m·s-1)2.4862.538-4.680.7880.31.8×10-8
To verify the feasibility of the proposed resilient modulus constitutive model, triaxial compression tests were simulated under the CD condition. A whole specimen with a diameter of 39.1 mm and a height of 80 mm was meshed with 8-node trilinear pore pressure-displacement brick element (C3D8P), as shown in Fig.3. A coupled hydro-mechanical analysis, namely, soil step in code Abaqus, was conducted by applying displacement load at the top end of the specimen, and the loading velocity was small enough (0.01%/min) for the dissipation of excessive pore pressure. The other boundary conditions were consistent with laboratory tests. The simulations started with a geostatic step, in which the specimen achieved an initial equilibrium state under the confining pressure and suction.
Fig.3 FE mesh of a triaxial test specimen
Stress-strain curves obtained from laboratory tests and numerical simulations are plotted in Fig.4. Apparently, numerical results are very close to the measured ones, indicating that the developed constitutive model is capable of capturing nonlinear behavior of granular materials at different suctions. In addition, the good agreements indicate the feasibility and correctness of the developed UMAT subroutine.
Fig.4 Measured and computed stress-strain curves at different matric suctions
2.3 Model application in coupled hydro-mechanical pavement analysis
To demonstrate the effectiveness of the proposed constitutive model, 3D FE simulations were performed on a flexible pavement under different moisture states. The pavement was composed of a 127 mm hot mixture asphalt (HMA) surface layer and a 305 mm aggregate base layer. A FWD test load with 150 mm radius was applied on the pavement surface. This load was a half-sine impact load of 0.03 s duration with 53 kN maximum amplitude.
A linear viscoelastic constitutive model in terms of Prony series is used to characterize the HMA mechanical behavior,
(18)
(19)
where G(t) and K(t) are the relaxation shear modulus and bulk modulus, respectively; G0 and K0 are the instantaneous shear modulus and bulk modulus; t is the relaxation time; Gi, Ki and τi are the Prony series parameters. The fitting Prony series for the HMA relaxation modulus at a reference temperature of 25 ℃ according to Ref.[6] are shown in Tab.2, thus herein, the instantaneous modulus is 4 170 MPa and Poisson’s ratio is 0.35.
Tab.2 Viscoelastic prony series for HMA
igi or kiτi10.3380.00220.2410.130.150140.1087.150.06039.460.044153
The proposed resilient modulus model was used as the constitutive model of granular base material (GB) and subgrade material (SG). The parameters used and SWCC for both materials are adopted from Cary et al.[12], as shown in Tab.3 and Fig.5.
Tab.3 Material properties of the GB and SG materials
MaterialGBSGDry density/(g·cm-3)2.261.91Poisson s ratio0.350.35Initial void ratio0.220.5Saturated permeability coefficient/(m·s-1)8×10-35.5×10-10k1/kPa1 992.891 480.32k20.6280.42k3-0.294-2.982k40.181.65
Fig.5 SWCCs of the GB and SG materials[12]
A quarter of the cube model was used in the analysis due to symmetry. The horizontal remote boundaries were truncated at a distance of 6 m away from the load center, and the total depth of the model was 3 m. The bottom of the model was assumed to be fixed in all directions, whereas displacements were restrained only in the normal direction at the side.
With respect to the hydraulic boundary conditions, no flow was allowed through the bottom and sides of the model in the FWD loading step, and the HMA surface was treated as impervious. Three moisture states were inputted as initial conditions of the model, which are as follows.
1) Reference moisture state: This state was a rather dry state in newly constructed pavement. Moisture contents of the base layer and subgrade were around optimum (The corresponding matric suction was 81 kPa), and the groundwater table was below the bottom of the model.
2) Moisture state 1: The groundwater table was 2 m below the subgrade top surface.
3) Moisture state 2: The groundwater table was at the same level as the subgrade top surface.
The FE model is shown in Fig.6. An 8-node brick element (C3D8) is used for the HMA layer, and an 8-node trilinear pore pressure-displacement brick element (C3D8P) is used for both the base layer and subgrade. Fine mesh is generated in the region adjacent to the loading area, while coarser mesh is generated as the distance to the load center increases.
Fig.6 The FE model of a 3D pavement
The simulations started with a geostatic step to establish an initial equilibrium state under the gravity load and boundary conditions, so that the initial stiffness matrix of nonlinear materials in the base and subgrade can be formed correctly. Subsequently, a transient soil step was created to simulate the coupled deformation-pore pressure response of the pavement under the FWD load.
3 Results and Discussion
Fig.7 illustrates the distributions of the calculated modulus (SDV1) near the load center in moisture state 2 before and after the load application. It is found that the calculated modulus near the load center in the base layer increases since net bulk stress grows, representing the stress-hardening behavior of base granular materials. However, the moduli beneath the load center in top subgrade decrease due to the shear stress and excessive pore water pressure build-up, revealing the stress-softening behavior of subgrade soils. Actually, these results can be forecasted with the relative values of regression parameters (k2, k3, k4), which determine the extent that different factors affect the modulus.
Fig.7 Isosurface of calculated modulus for pavement subjected to moisture state 2.(a) Before the load application (t=0 s). (b) At the peak load time (t=0.15 s).
(a)
(b)
Fig.8 In-depth distributions.(a) Modulus; (b) Pore pressure beneath the load centerline at the peak load time
Fig.8 presents the in-depth distributions of the calculated modulus and pore pressure beneath the load centerline at peak load time under three specific moisture states. It can be seen that matric suction (negative pore pressure) above groundwater is nearly linear with distance from the groundwater table for moisture states 1 and 2 in Fig.8(b), which indicates that matric suction distribution will affect the value of the resilient modulus. Excessive pore pressure generation is observed at the groundwater table level under the FWD load, which substantially decreases the modulus of soil. Also, the increase in the groundwater table results in pore pressure build-up and reduction of the modulus consequently. The modulus at the top of the base and subgrade decreases by 48% and 68% from the reference moisture state to state 2.
The surface deflections under FWD loading are widely used for evaluating the structural bearing capacity of the pavement structure.Here, the surface deflection profiles at different moisture states are plotted in Fig.9. An increase of the deflection basin area with the increase of the groundwater table can be clearly observed. The maximum deflections from the reference state to state 1 and from state 1 to state 2, increase by 61% and 25%, respectively. Similar results were observed by Salour et al.[24], who conducted field measurements by FWD and found that the increase in the groundwater table level from 2.5 to 1.0 m below the pavement surface can increase the maximum deflection by 36%.
Fig.9 Surface deflection at the peak load time
Some representative mechanical responses related to major distress encountered in flexible pavements are summarized in Tab.4 at positions below the load centerline. The results show that in moisture state 2, the tensile strain at the bottom of the HMA and the compressive strain on the top of subgrade are 53% and 64% greater than those in the reference state, respectively. It implies that high groundwater table level accelerates failure of the pavement structure, such as premature fatigue cracking in the HMA layer and excessive permanent deformation (rutting) of the subgrade. Salour et al.[24] drew a similar conclusion by back-calculating the critical strain above based on the FWD data measured at different groundwater levels. In summary, the detrimental effect on pavement performance due to the groundwater table rising is remarkable.
4 Conclusions
Tab.4 Summary of predicted pavement responses at different moisture states
Moisture stateTensile strain at bottom of HMA/10-6Compressive strain on top of subgrade/10-6Maximum deflection at pavement surface/mmCalculated modulus on top of subgrade/MPaReference state197.51 043.10.343161.6State 1253.21 333.80.55280.9State 2302.31 715.60.68951.3
1) A suction-dependent resilient modulus constitutive model (Eq.(12)) is successfully implemented in numerical simulation by means of UMAT subroutine in code Abaqus. The developed model can characterize the nonlinear stress-strain relationship under different suction conditions, and good agreements are achieved by comparing predicted and measured results of triaxial tests.
2) The developed model is integrated into a 3D pavement FE model and hydro-mechanical analyses are conducted to investigate the effect of moisture on the pavement response. It is found that the groundwater table rising decreases the modulus of granular material in pavement due to suction reduction, and then increases the maximum surface deflection.
3) Based on the results of numerical simulations, the critical strains used in pavement design are greater under a higher groundwater table condition, leading to shorter estimated pavement life when subgrade rutting and fatigue cracking are considered.
Only one of the suction-dependent resilient modulus models is implemented in this study. Comparisons between different models are recommended for future studies.
[1]Hicks R G, Monismith C L. Factors influencing the resilient response of granular materials [J]. Highway Research Record, 1971(345): 15-31.
[2]Uzan J. Characterization of granular material [J]. Transportation Research Record, 1985, 1022: 52-59.
[3]Hjelmstad K D, Taciroglu E. Analysis and implementation of resilient modulus models for granular solids [J]. Journal of Engineering Mechanics, 2000, 126(8): 821-830. DOI:10.1061/(asce)0733-9399(2000)126:8(821).
[4]Kuo C M, Huang C W. Three-dimensional pavement analysis with nonlinear subgrade materials [J]. Journal of Materials in Civil Engineering, 2006, 18(4): 537-544. DOI:10.1061/(asce)0899-1561(2006)18:4(537).
[5]Wang H, Imad L A Q. Importance of nonlinear anisotropic modeling of granular base for predicting maximum viscoelastic pavement responses under moving vehicular loading [J]. Journal of Engineering Mechanics, 2013, 139(1): 29-38. DOI:10.1061/(asce)em.1943-7889.0000465.
[6]Li M Y, Wang H, Xu G J, et al. Finite element modeling and parametric analysis of viscoelastic and nonlinear pavement responses under dynamic FWD loading [J]. Construction and Building Materials, 2017, 141: 23-35. DOI:10.1016/j.conbuildmat.2017.02.096.
[7]Gu F, Luo X, Luo R, et al. Numerical modeling of geogrid-reinforced flexible pavement and corresponding validation using large-scale tank test [J]. Construction and Building Materials, 2016, 122: 214-230. DOI:10.1016/j.conbuildmat.2016.06.081.
[8]Claudia Z, Yugantha P, William H. Matric suction prediction model in new aashto mechanistic-empirical pavement design guide [J]. Transportation Research Record, 2009, 2101: 53-62. DOI:10.3141/2101-07.
[9]Khoury N, Zaman M. Correlation between resilient modulus, moisture variation, and soil suction for subgrade soils [J]. Transportation Research Record, 2004, 1874: 99-107. DOI:10.3141/1874-11.
[10]Han Z, Vanapalli S K. Model for predicting resilient modulus of unsaturated subgrade soil using soil-water characteristic curve [J]. Canadian Geotechnical Journal, 2015, 52(10): 1605-1619. DOI:10.1139/cgj-2014-0339.
[11]Liang R Y, Rabab’Ah S, Khasawneh M. Predicting moisture-dependent resilient modulus of cohesive soils using soil suction concept [J]. Journal of Transportation Engineering, 2008, 134(1): 34-40. DOI:10.1061/(asce)0733-947x(2008)134:1(34).
[12]Cary C E, Zapata C E. Resilient modulus for unsaturated unbound materials [J]. Road Materials and Pavement Design, 2011, 12(3): 615-638. DOI:10.1080/14680629.2011.9695263.
[13]Wang L J, Liu S H, Fu Z Z, et al. Coupled hydro-mechanical analysis of slope under rainfall using modified elasto-plastic model for unsaturated soils [J]. Journal of Central South University, 2015, 22(5): 1892-1900. DOI:10.1007/s11771-015-2708-2.
[14]Saad B. Analysis of excess water impact on the structural performance of flexible pavements [J]. International Journal of Pavement Engineering, 2014, 15(5): 409-426. DOI:10.1080/10298436.2013.790546.
[15]Hu R, Chen Y F, Liu H H, et al. A coupled two-phase fluid flow and elastoplastic deformation model for unsaturated soils: Theory, implementation, and application [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(7): 1023-1058. DOI:10.1002/nag.2473.
[16]Gu F, Luo X, Zhang Y, et al. Modeling of unsaturated granular materials in flexible pavements [C]// E3S Web of Conferences. Paris. France, 2016:20002. DOI:10.1051/e3sconf/20160920002.
[17]Bishop A W. The principle of effective stress [J]. Teknisk Ukeblad, 1959, 39: 859-863.
[18]Khalili N, Khabbaz M. Application of effective stress concept to unsaturated soils [C]// Proceedings of the 8th Australia New Zealand Conference on Geomechanics: Consolidating Knowledge. Hobart, Australia, 1999: 849-854.
[19]Komolvilas V, Kikumoto M. Fully undrained cyclic loading simulation on unsaturated soils using an elastoplastic model for unsaturated soils [C]// E3S Web of Conferences. Paris, France, 2016, 9: 17008. DOI:10.1051/e3sconf/20160917008.
[20]Fredlund D G, Xing A. Equations for the soil-water characteristic curve [J]. Canadian Geotechnical Journal, 1994, 31(4): 521-532. DOI:10.1139/t94-061.
[21]Seed H, Mitry F, Monosmith C, et al. Prediction of pavement deflection from laboratory repeated load tests. NCHRP Report No.35 [R]. Washington, DC: National Cooperative Highway Research Program, 1967.
[22]Inc. ARA. Guide for the mechanistic empirical design of new and rehabilitated pavement structures. Final report. NCHRP 1-37A [R]. Washington, DC: Transportation Board of the National Academies, 2004.
[23]Chen R. Experimental study and constitutive modelling of stress-dependent coupled hydraulic hysteresis and mechanical behaviour of an unsaturated soil [D]. Hong Kong: Civil Engineering Department,Hong Kong University of Science and Technology, 2007.
[24]Salour F, Erlingsson S. The influence of groundwater level on the structural behaviour of a pavement structure using fwd [C]//The Ninth International Conference on the Bearing Capacity of Roads, Railways, and Airfields. Trondheim, Norway, 2013: 25-27.