1 Introduction

Heavy-haul railway has been a worldwide transportation mode that improves the train's weight substantially. A heavy-haul railway has been constructed by Shenhua Group in Northwest China, where subgrade filling materials are typically subjected to freeze-thaw cycles. Therefore, railways deteriorate with time due to traffic loading and environmental exposure(Pylkkänen et al., 2012). Embankment deformation is one of the important research topics in the railway field, where deformation not only increases the cost of railway operation, but also influences transportation safety(Ma et al., 2000; Li GY et al., 2009; Thomas et al., 2009). Research into improvements of subgrade filling materials and deformation of railway evaluation subjected to freeze-thaw cycles, including laboratory testing and numerical simulation, can extend the railway life and lead to significant cost savings in railway maintenance and rehabilitation.

The techniques available to predict deformations in heavy-haul railway's embankment include multilayer elastic theory, and analytical, finite difference, finite element, boundary element, and hybrid methods. The development of numerical tools and methodologies that can accurately predict the deformation performance of high embankment in heavy-haul railway when subjected to freeze-thaw cycles has become one of the main issues of research in the past few decades(Belforte et al., 2008; Geng et al., 2008; Xiao et al., 2011; Mehdipour et al., 2013; Paixão et al., 2015). Although using three dimensional(3D)finite element models can solve all problems that can be solved with two-dimensional(2D)models, it is very expensive to use 3D models in terms of data preparation and computational time. The railway structure has a large longitudinal dimension, which makes it suitable for using 2D plane strain models when studying transverse profiles. Therefore, 2D plane strain models were appropriate for this study. A plane strain model assumes that thickness in the horizontal plane is infinite. Geometry in the horizontal plane does not change so that strain is nil. Finite element analysis methods have been widely used to predict railway embankment deformation, and improve guidelines by investigating the response of railway track to various factors affecting its performance(Shahin and Indraratna, 2006). The developed finite element track model is also used to demonstrate the capability of different soil stabilization techniques in enhancing track performance(Dahlberg, 2001; Puzavac and Popović, 2010).

An important step in numerical analysis of subgrade filling materials is the selection of a suitable constitutive equation to model the behavior of a railway under an external load. Despite the large attention that numerical models of the track has had for the past few decades, very seldom has the soil's nonlinear response been taken into account in the simulation of track under heavy-haul railway subjected to freeze-thaw cycles. In conventional railway analysis methods, subgrade filling materials on the surface of the embankment are assumed to be linear and elastic. However, a linear elastic model cannot always properly model the behavior of subgrade filling materials( LiSY et al., 2009; Wang and Al-Qadi, 2012; Gomes and Cunha, 2014).

Therefore, a nonlinear model has been developed in this paper that is suitable for railway analysis subjected to freeze-thaw cycles with different moisture content. The applicability and implementation of the nonlinear constitutive model are presented. In this study, taking the subgrade structure system of Ba-quasi heavy-haul railway as research background, the deformation characteristics of the subgrade under various moisture contents subjected to freeze-thaw cycles are studied by using a 2D finite element model considering nonlinear subgrade properties.

2 Numerical simulation 2.1 Constitutive model

The material model that is used for strain consideration of subgrade filling materials is a softening-type stress-strain nonlinear model based on the Duncan-Chang model(Duncan and Chang, 1970). It is developed to analyze the effects of freeze-thaw cycles with different moisture content according to results of the static triaxial tests.

Figure 1 shows the stress-strain curves of the unfrozen and thawed filling materials after the freeze-thaw cycles with different moisture content. It is clearly seen in Figure 1 that the freeze-thaw processes did not affect the patterns of the stress-strain curves. Therefore, a softening-type stress-strain model with a high degree of accuracy can be used to approximate the nonlinear stress-strain curves of filling materials after being subjected to freeze-thaw cycles as shown in Figure 2. These parameters(i.e., R_σ, R_ε, R_σt, and R_εt)are defined by Equations (1a)-(1d) to describe the relationships between stress and strain as follows:

${R_\sigma } = \frac{{{\sigma _1} - {\sigma _3}}}{{{{\left( {{\sigma _1} - {\sigma _3}} \right)}_p}}}$

(1a)

${R_\varepsilon } = \frac{\varepsilon }{{{\varepsilon _p}}}$

(1b)

${R_{\sigma t}} = \frac{{{{\left( {{\sigma _1} - {\sigma _3}} \right)}_t}}}{{{{\left( {{\sigma _1} - {\sigma _3}} \right)}_p}}}$

(1c)

${R_{\varepsilon t}} = \frac{{{\varepsilon _t}}}{{{\varepsilon _p}}}$

(1d)

where(σ₁−σ₃)is the deviator stress; σ₁ and σ₃ are the major principal stress and minor principal stress, respectively; ε is the axial strain;(σ₁−σ₃)_p, ε_pare the deviator stress and axial strain, respectively, at the peak point(i.e., point P in Figure 2); and (σ₁−σ₃)_t and ε_t are the deviator stress and axial strain, respectively, at the inflection point(i.e., point T in Figure 2).

After careful inspection, the piecewise continuous functions for the proposed constitutive model can be expressed as follows:

${R_\sigma } = D\frac{{1 + {b_2}{R_\varepsilon }}}{{1 + {b_3}{R_\varepsilon }}} \cdot \frac{{{R_\varepsilon }}}{{1 + {b_1}{R_\varepsilon }}}\begin{array}{*{20}{c}} {}&{} \end{array}\left( {{R_\varepsilon } \le {R_{\varepsilon t}}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} $

(2a)

${R_\sigma } = {R_{\sigma t}} - \frac{{{R_\varepsilon } - {R_{\varepsilon t}}}}{{{b_4} + {b_5}\left( {{R_\varepsilon } - {R_{\varepsilon t}}} \right)}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {{R_\varepsilon } > {R_{\varepsilon t}}} \right){\kern 1pt} {\kern 1pt} $

(2b)

where D, b₁, b₂, b₃, b₄ and b₅ are undetermined coefficients that can be obtained from experimental results, and D = E_max/E_p, where E_max and E_p(i.e., E_p =(σ₁−σ₃)_p/ε_p)are the initial tangent modulus at the original point(i.e., point O in Figure 2) and the secant modulus at the peak point, respectively. The determination of E_max is equivalent to the determination of E_max in the Duncan-Chang model(Duncan and Chang, 1970).

The stress-strain curves obtained from static triaxial tests with moisture content of 6.5%, 8.0% and 9.5% subjected to zero freeze-thaw cycles(e.g., unfrozen filling materials)are fitted by Equations (2a) and (2b), and then the parameters D, b₁, b₂, b₃, b₄, and b₅ can be attained. It is clear that all parameters are dependent on moisture content, as shown in Figure 3. For convenience, a unified equation f_x(w)was selected for the three moisture contents to fit the undetermined parameters in Equations (2a) and (2b). Meanwhile, to directly express the relationships between b₁, b₂, b₃, b₄, b₅, D and the freeze-thaw cycles, N_FT(0≤ N_FT ≤10). R_x(N_FT)is introduced as reduced factors of six parameters. These parameters can be evaluated by the computed values of b₁, b₂, b₃, b₄, b₅, and D, which are normalized with the corresponding values for the unfrozen filling materials(Figure 4).

Afterwards, the parameters of different moisture content and freeze-thaw cycles can be evaluated by Equations (3a) to (3f) based on the results of the static triaxial tests with moisture content of 6.5% subjected to zero freeze-thaw cycles. Research in this way leads to considerable reduction of the duration without influencing results.

$\left\{ \begin{array}{l} {f_b}_1\left( w \right){\rm{ }} = {\rm{ }}3.343{{\rm{e}}^{ - 0.02911w}}\\ {R_{b1}}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = 0.9158{{\rm{e}}^{ - 0.1151}}^{{N_{{\rm{FT}}}}} \end{array} \right.$

(3a)

$\left\{ \begin{array}{l} {f_b}_2\left( w \right) = {\rm{ }} - 0.2385{{\rm{e}}^{ - 0.003976w}}\\ {R_{b2}}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = {\rm{ }}1.019{{\rm{e}}^{0.04431{N_{{\rm{FT}}}}}} \end{array} \right.$

(3b)

$\left\{ \begin{array}{l} {f_b}_3\left( w \right) = {\rm{ }} - 0.1513{{\rm{e}}^{ - 0.2247w}}\\ {R_{b3}}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = {\rm{ }}0.7374{{\rm{e}}^{ - 0.05758{N_{{\rm{FT}}}}}} \end{array} \right.$

(3c)

$\left\{ \begin{array}{l} {f_{b4}}\left( w \right) = {\rm{ }}3.07{{\rm{e}}^{ - 0.03084w}}\\ {R_{b4}}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = {\rm{ }}0.8882{{\rm{e}}^{ - 0.07142{N_{{\rm{FT}}}}}} \end{array} \right.$

(3d)

$\left\{ \begin{array}{l} {f_b}_5\left( w \right) = {\rm{ }}9.447{{\rm{e}}^{ - 0.1528w}}\\ {R_{b5}}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = {\rm{ }}1.038{{\rm{e}}^{ - 0.08935{N_{{\rm{FT}}}}}} \end{array} \right.$

(3e)

$\left\{ \begin{array}{l} {f_D}\left( w \right) = {\rm{ }}5.206{{\rm{e}}^{ - 0.01481w}}\\ {R_D}\left( {{N_{{\rm{FT}}}}} \right){\rm{ }} = {\rm{ }}0.9252{{\rm{e}}^{ - 0.05567{N_{{\rm{FT}}}}}} \end{array} \right.$

(3f)

2.2 Numerical implementation of the constitutive model

The parameters of nonlinear constitutive model for filling materials include: b₁, b₂, b₃, b₄, b₅, R_εt, ε_p, E_max. The value of all the above mentioned parameters should be properly determined by using the method as described in Section 2.1. Triaxial test results are used to calculate parameters D, b₁, b₂, b₃, b₄, and b₅ of unfrozen filling materials with 6.5% moisture content for Equations (2a) and (2b) in the model, as shown in Table 1. Then, the parameters of filling materials with other moisture content subjected to various freeze-thaw cycles can be evaluated by multiplying by f_x(w) and R_x(N_FT), as shown in Equations (3a) to (3f), respectively.

ABAQUS software offers an interface program called UMAT for secondary development. Here, the user subroutine UMAT is adopted to complete the numerical implementation of the nonlinear constitutive model of filling materials. In this process, the tangent modulus is the core problem. Therefore, a UMAT is developed to calculate the tangent modulus of filling materials corresponding to the stress state which is initially obtained from the main program. The Jacobian Matrix defined as Equation (4) of integration points of all elements should be required, and the stress tensor and state variables for integration points should be updated. In the end, the updated stress state returns to the main program to calculate the deformation.

$\left[ D \right] = \frac{{\partial \Delta \bar \sigma }}{{\partial \Delta \bar \varepsilon }}$

(4)

where [D] = Jacobian matrix(incremental stiffness matrix); $\partial \Delta \bar \sigma $= increment in stress; and $\partial \Delta \bar \varepsilon $= increment in strain. The UMAT program requires eight parameters(shown in Table 1)for calculating the nonlinear tangent modulus. The tangent modulus before and after the inflection point in the curve can be obtained from Equations (5a) and (5b), respectively.

$\begin{array}{l} \overline E = \frac{{{\rm{d}}\left( {{\sigma _1} - {\sigma _3}} \right)}}{{{\rm{d}}\varepsilon }}\\ \quad \;\, = \frac{{{E_{\max }}}}{{\left( {1 + {b_1}{R_\varepsilon }} \right)\left( {1 + {b_3}{R_\varepsilon }} \right)}}\left[ {\frac{{\left( {{b_1} - {b_2}} \right){R_\varepsilon }}}{{1 + {b_3}{R_\varepsilon }}} + \frac{{1 + {b_2}{R_\varepsilon }}}{{1 + {b_1}{R_\varepsilon }}}} \right]\\ \quad \quad \quad \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad ({R_\varepsilon } &1t; {R_{\varepsilon t}}) \end{array}$

(5a)

$\begin{array}{l} \overline E = \frac{{{\rm{d}}\left( {{\sigma _1} - {\sigma _3}} \right)}}{{{\rm{d}}\varepsilon }}\; = \frac{{\overline {{E_t}} }}{{{{\left[ {1 + \frac{{{b_5}}}{{{b_4}}}\left( {{R_\varepsilon } - {R_{\varepsilon t}}} \right)} \right]}^2}}}\\ \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad ({R_\varepsilon } > {R_{\varepsilon t}}) \end{array}$

(5b)

Based on the above-mentioned theory, the nonlinear constitutive model of filling materials is embedded in ABAQUS. Then a numerical example is presented and analyzed in the following. A series of static triaxial numerical tests for a specimen of filling materials whose confining pressure is 200 kPa in different moisture content and freeze-thaw cycles are conducted. Moisture content starts from 6.5% to 9.5% with step by 1.5% and the numbers of freeze-thaw cycles are 0, 1, 3, 6, and 10, respectively. Figure 5 shows the sample in the static triaxial and the calculation results of the corresponding model in ABAQUS. The model parameters used in this analysis is calculated as the above mentioned.

Results of static triaxial numerical tests and comparison with the experimental data are shown in Figure 6. The points represent experimental results, meanwhile, the curve is obtained from the static triaxial numerical tests conducted in ABAQUS. The comparisons in this step show that the largest difference between the calculated stress-strain curve and measured in static triaxial tests is about 15.9% with the moisture content of 9.5% subjected to three freeze-thaw cycles. By this value, it is concluded that results of the experiment and prediction of model inosculate, and the correctness of the proposed program development using the proposed model can be concluded.

2.3 Finite element model

The section DK97+454 on Bazhun heavy-haul railway is taken as a typical calculation model with embankment height of 21 m, and embankment top surface width of 12 m. The stratigraphic structure, as shown in Figure 7, is determined according to the site survey.

The commercial finite element software ABAQUS is used to create the model and conduct the analysis. Plane strain analysis can be carried out due to the fact that the railway line has semi-infinite characteristic. Taking advantage of the symmetry condition, only half of the model is meshed. The quasi-static finite element embankment-field coupling model is shown in Figure 8. The full model is 57.5 m in the longitudinal direction and 63 m in depth. Therein, the 4-node bilinear plane strain quadrilateral is used to simulate the high embankment, which is composed of 1, 098 nodes and 1, 020 elements. The nonlinear constitutive model developed in this paper is employed in the top layer of 2.5 m to simulate the properties of the filling materials(Group A)which is in the range of freeze-thaw zone, while the remaining depth of embankment fill and strata layer is based on Mohr Coulomb criterion. The physical and mechanical parameters are adopted mainly based on the site survey in Bazhun heavy-haul railway, as well as other relevant references concerning soil properties(Chu et al., 1999). The parameters of embankment fill and soil strata are shown in Table 2. The nodes at the bottom of the model are fixed in the vertical and horizontal direction, while only horizontal displacement is constrained at the left of the model. At the symmetry plane only translations in the plane and rotations around the normal vector to the plane are allowed.

2.4 Quasi-static loads

The train live load pattern C80 is demonstrated by four 250 kN axles. For simplicity, the train dynamic load is simulated through the definition of stationary loads. This load pattern should be distributed uniformly across the sleeper and ballast on the railway embankment. Furthermore, the whole concentrated load of the C80 pattern is assumed to be applied to the adjacent sleepers and to be distributed on the track subgrade through a ballast layer with 35 cm thickness and 1:1(H:V)distribution pattern according to fundamental code for design on railway bridge and culvert(CRFSDIG, 2005). The load pattern distributed on the track subgrade level is illustrated in Figure 9. According to the illustrated load pattern, the train moving load could be substituted with a rectangular soil column whose unit weight is identical to the substratum. Therefore, an equivalent uniform live load of 61.8 kPa provided in the code for design heavy-haul railway(CRFSDIG, 2011)is used on the surface of the high embankment to simulate the load passing from superstructure(rail, sleeper and ballast)to substructure.

3 Numerical results for the deformation of railway embankment

By using ABAQUS software, a two-dimensional railway track with plane strain conditions was analyzed. To summarize the output of the numerical models, the results are presented using 2D plots regarding the deformation of railway embankment, with respect to the moisture content, freeze-thaw cycles and depth. It presents the deformation of railway subgrade predicted in ABAQUS at the moisture content of 6.5%, 8.0% and 9.5% subjected to different freeze-thaw cycles in Figure 10. It can be concluded that the vertical deformation of embankment central line(as shown in Figure 7)may sharply decrease with increased depth of the railway subgrade. The results highlight the fact that the freeze-thaw processes and moisture content did not affect the patterns of the deformation curves with the depth.

Considering all the calculated results along the vertical direction, it is interesting to compare the deformation of the surface of subgrade bed and embankment fill with different moisture content subjected to freeze-thaw cycles(Figure 11). It is shown that deformation may increase with increased freeze-thaw cycles in the moisture content of 6.5%, 8.0% and 9.5%. Figure 11a indicates that embankment deformation with moisture content of 6.5%(the optimum moisture content)reaches significant deformation variations subjected to the first freeze-thaw cycle, and becomes stable after six freeze-thaw cycles. However, embankment deformation with moisture content of 8.0% and 9.5%(Figures 11b and 11c)increases continuously even after ten freeze-thaw cycles. The freeze-thaw process is a process by which soil in an unstable state develops into a dynamic stable state, and repeated freeze-thaw cycles will change the soil's structure towards a new dynamic stable balance state(Yang et al., 2004). After six cycles of freeze-thaw, embankment deformation with moisture content of 6.5% arrives at a constant value, reflecting that a dynamic equilibrium is reached in the specimen internal fabric. Those cases with higher moisture content, however, may reach a stable state with more freeze-thaw cycles, owing to a mass of pore water frozen and migration. Therefore, additional experiments should be conducted to investigate micro mechanism, such as appearance of ice crystals, migration of pore water, and reorientation of the filling materials, to make sure how and when filling materials become stable. From the results in this paper, it is seen that the optimum moisture content used in construction can not only reduce deformation, but also stabilize the deformation in a changing climate condition.

4 Conclusions

This study provides a simple nonlinear finite element approach for predicting embankment deformation subjected to freeze-thaw cycles with different moisture content. To achieve this goal, the ABAQUS finite element program, with two-dimensional plane strain analysis, is used for simulation of the railway subgrade, and a nonlinear constitutive model based on static triaxial tests is developed and used in the finite element model. The following conclusions can be drawn.

1)A practical constitutive model is developed to represent the stress-dependent and nonlinear stress-strain behaviors of the granular materials with different moisture content subjected to freeze-thaw cycles. A simple technique is used to interpret the experimental results and to conveniently determine the parameters(b₁, b₂, b₃, b₄, b₅, R_εt, ε_p, E_max)in the model.

2)Static triaxial tests are simulated by ABAQUS to calibrate the nonlinear constitutive model. Results show that there is a good relationship between the stress-strain curves predicted in the finite element model according to the nonlinear constitutive model and the curves measured in the static triaxial tests. It is reasonable to use the nonlinear constitutive model to predict the deformation of railway subgrade in changing climate conditions.

3)The deformation of the subgrade gradually increases as the number of freeze-thaw cycles increases. The most significant variation in deformation is observed after the initial cycle. After six cycles, embankment deformation with content moisture of 6.5% remained nearly constant regardless of additional freeze-thaw cycles, however, embankment deformation with content moisture of 8.0% and 9.5% will not achieve a stable state with increased freeze-thaw cycles. These results indicated that the optimum moisture content can be used for engineering design in regions that seasonally freeze.

Acknowledgments:

This research is supported by the National Natural Science Foundation of China(Grant No. 41430634), the Foundation Project Program 973 of China(No. 2012CB026104), the Foundation Project Program of SHENHUA BAOSHEN Railway Corporation Limited(No. 201212240384), Technology Research and Development Plan Program of Heilongjiang Province, China(No. GZ13A009), and State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining & Technology(Grant No. SKLGDUEK1209).