A Simplified Calculation Model for Determining the Permanent Deformation of Soils Under Flexible Pavements Due to Repeated Traffic Load

1. Introduction

In flexible pavement structures, ruts can be also formed besides flexural and fatigue cracks, which can be caused by inadequate asphalt mixture, improperly compacted asphalt, or subsequent compaction of the soil caused by traffic. Conventional design procedures [1,2] do not consider permanent deformations from the compaction of the soil under traffic loads, and rutting is prevented by an appropriately chosen asphalt mixture and its careful compaction on site. Some analytical design procedures [3,4] assume that rutting is of soil origin only, and can be reduced by limiting the vertical deformation at the top of the subgrade, but does not quantify the rutting depth. In his Ph.D. thesis, [5] concludes based on in-situ measurements that permanent settlements may not only be of soil origin, but the unbound granular base course can also significantly contribute to the rutting depending on the thickness of the asphalt course. [6,7,8] developed calculation procedures for the permanent deformation of unbound granular materials under repeated traffic loads.

There is no generally accepted, widely recognized and documented procedure for calculating permanent deformations in unbound granular layers (base and subbase courses and subgrade) beneath the asphalt and thus for determining the depth of ruts. Some methods [8] only focus on the deformation of the granular base course, which contradicts the findings of [3]. In other cases [7], the material constants or the calibration of the parameters are not documented, therefore these models cannot be used in different boundary conditions. When studying rutting, priority should be given to methods that can determine permanent deformation of each layer under repeated traffic loads.

One of the best documented, most advanced and most reliable such material model to determine permanent deformations under repeated loads that is also supported by numerous laboratory tests is the High-Cycle Accumulation (HCA) model or the high-cycle accumulation model presented in [9], which was verified by several model tests in [10]. The HCA model was developed for sandy soils under low intensity (ε^ampl < $3·{10}^{-3}$), high-cycle (even 2 million cycles) loading. The model was calibrated with hundreds of cyclic triaxial laboratory tests presented in [10,11,12,13,14,15,16].

2. Explicit Calculation Method and HCA Model

2.1. Explicit and Implicit Methods

Determination of permanent deformations has two fundamental approaches. Implicit methods directly calculate each cycle using a material model based on stress-strain increments $(\dot{\sigma }-\dot{\epsilon })$ and the residual deformation appears as a "by-product" of the calculation. The biggest disadvantage of these methods is that the calculation (integration) and material model errors are gradually accumulated with each load cycle, therefore they can only be used reliably for N< 50 cycles [17].

Due to the above reasons, explicit calculation methods come to the fore in case of high-cycle loading. Instead of calculating the deformation trajectory at each load step and each cycle, the accumulated permanent deformation (ε^acc) is calculated in pseudo-time. The procedure is similar to creep and viscoplastic models, but instead of time (t) the number of cycles (N) is used. Explicit methods are empirical in nature, as the observed behavior of the soil in laboratory is described with mathematical formulae.

2.2. HCA Model

The basic principle of the HCA model is to combine the increment-based $(\dot{\sigma }-\dot{\epsilon })$ advanced implicit material model and the cycle number–permanent deformation $(N-{\epsilon }^{acc})$ explicit model. As a final result, the HCA model gives the accumulated strain ε^acc as a vector quantity under a cyclic loading with N-number of constant stress amplitudes. The HCA model requires the void ratio, the mean effective normal stress, the stress state and the cyclic strain amplitude (ε^ampl) as input parameters. The latter is implicitly calculated using a conventional material model by gradually increasing the cyclic stress and calculating the strain amplitude as the sum of the strain increments. Then, the strain amplitude is used as input to determine explicitly the permanent deformation under a cyclic load with N-number of constant stress amplitudes. The model is presented in detail in the habilitation of [18]; hereinafter only the most important calculation steps are discussed.

In the first step, the strain amplitude ε^ampl is determined due to q^ampl in a traditional implicit way, which is to be interpreted as the vector sum of the vertical and horizontal strains (1). The strain amplitude is defined as the mean value of the strains calculated at q^max − q^min points.

$${\epsilon }^{ampl}=\sqrt{{\left({ϵ}_1^{ampl}\right)}^2+{\left({ϵ}_2^{ampl}\right)}^2+{\left({ϵ}_3^{ampl}\right)}^2}$$

(1)

Then, the calculation switches to the explicit part of the HCA model, in which the permanent strains due to ΔN cycle package are determined with the above-defined ε^ampl. Since the soil structure and state parameters change in the long term, the explicit calculation is interrupted from time to time and the implicit analysis is repeated in control cycles in order to continue the explicit calculation for another ΔN cycle package with the actual value of ε^ampl.

By many constitutive models, when a stress path touches the yield surface, a plastic strain occurs. If the model takes hardening into account, then the yield surface shifts. For the HCA calculation any suitable constitutive relationship as plastic yield surface or strain hardening might be considered. The direction of the strain (e.g. volumetric or deviatoric) will be described with a flow rule. The difficulty by high cycle accumulations models is that the intensity of loading is somewhat moderate, a cyclic plastic strain develops over the years of service due to the thousands and millions of load repetition although the static yield requirements are not satisfied. For the general case, the strains are determined using the basic constitutive equation of the HCA model Equation (2):

$$\dot{\mathit{\sigma }}=\mathit{E}:(\dot{\mathit{\epsilon }}-{\dot{\mathit{\epsilon }}}^{\mathit{a}\mathit{c}\mathit{c}}-{\dot{\mathit{\epsilon }}}^{\mathit{p}\mathit{l}})$$

(2)

where:

$\dot{\mathit{\sigma }}$ is the Cauchy stress tensor,

$\dot{\mathit{\epsilon }}$ is the total strain rate,

${\dot{\mathit{\epsilon }}}^{acc}$ is the accumulated strain rate,

${\dot{\mathit{\epsilon }}}^{pl}$ is the plastic strain rate, which starts to accumulate when the loading reaches the yield surface,

E is the stress-dependent elastic stiffness matrix.

Equation (2) gives the change in the stress (e.g. increase in pore pressure in undrained conditions, $\dot{\mathit{\sigma }}\ne 0,\dot{\mathit{\epsilon }}=0$) or strain (e.g. accumulated plastic strain in drained conditions, $\dot{\mathit{\sigma }}=0,\dot{\mathit{\epsilon }}\ne 0$) depending on the boundary conditions. In the analysis performed for this research, the latter boundary conditions are valid, since drained conditions and sufficiently deep groundwater level were assumed. In this research it is assumed, that the stress path might only touch the yield surface during construction or by the first regular load cycle. Therefore no static plastic strain develops during the repeated cyclic loading, so that ${\dot{\mathit{\epsilon }}}^{pl}=0$. Thus, Equation (2) is simplified to ${\dot{\mathit{\epsilon }}}^{acc}=\dot{\mathit{\epsilon }}$. [10,18] provide further information on the plastic constitutive relationships of the HCA model.

The HCA model gives the rate of accumulation as a vector quantity, so in addition to the scalar magnitude of the strain, a flow rule is also needed to specify the direction of accumulation (m), i.e. to specify the deviatoric ε_q and volumetric ε_v parts of the strain. These can be used to determine the vertical strain components needed for the settlement analysis.

$${\dot{\mathit{\epsilon }}}^{\mathit{a}\mathit{c}\mathit{c}}={\dot{\epsilon }}^{acc}\mathit{m}$$

(3)

The scalar part of the rate of accumulation can be determined using the other basic constitutive equation of the HCA model Equation (4), which is given as the product of 5 empirical factors. The $f_{ampl}$, ${\dot{f}}_N$, $f_p$, $f_Y$, and $f_e$ account for the influence of strain amplitude, the number of cycles, the mean normal stress, the mean stress ratio and the void ratio, respectively. In some cases, a sixth term is also included in the equation, which describes the polarization of the cyclic loading $f_{\pi }$. However, this factor can be neglected in most practical applications and can be taken as 1.0 based on the recommendation of [18].

$${\dot{\epsilon }}^{acc}=f_{ampl}{\dot{f}}_N{f}_e{f}_p{f}_Y{f}_{\pi }$$

(4)

The above-mentioned 5 f_i functions can be described with 7 C_i material constants. Table 1 provides details of the constitutive relationships of the HCA model.

where:

${\epsilon }^{ampl}$ is the cyclic strain amplitude

$g^A$ is the preloading variable

${\overline{Y}}^{av}$ is the average normalized stress state according to Matsuoka and Nakai

$e$ is the void ratio

${\mathrm{e}}_{max}$ is the maximal void ratio at loosest state

In addition, the critical friction angle ${\phi }_{cc}$ is also needed to define the flow rule and $f_Y$, which is not necessarily the same as the critical friction angle ${\phi }_c$ that can be calculated from the monotonic CU triaxial test. In this research, the recommendations of [11] were used to determine the critical friction angle ${\phi }_{cc}$ of the cyclic flow rule.

2.3. Direction of Accumulation

Based on laboratory test results, [11,19] found that the flow rule of MCC is approximately valid for sandy soils in the case of cyclic loading, so the vector of cyclic flow is mainly governed by the ratio of deviatoric and hydrostatic stresses ${\eta }^{av}=q^{av}/p^{av}$ describing the average stress state. The relationship between the number of cycles and the strain rate ratio for large number of cycles (10⁵–2$·$10⁶) was investigated by [20] on 22 different particle size distribution curves of a clean quartz sand. Laboratory tests showed that the findings in [19] are somewhat inaccurate, as the plastic strain accumulation vector inclines towards the p axis, i.e. towards the volumetric strain, as the number of cycles increases. However, the MCC model is acceptable as an approximation for isotropic soils, since basically there are no shear strains in the case of ${\eta }^{av}$ = 0, while only shear strains develop in the case of ${\eta }^{av}$ = M_c. In this research, the MCC flow rule was applied to the sandy subgrade described in Section 4. Tests carried out on crushed stone base curse by [21] showed that, unlike sand, the material behaves anisotropically under compaction, so the flow rule of the MCC does not apply. Therefore an anisotropic flow rule presented in [21] were used in this research for the base and subbase layers.

2.4. Determination and Calibration of the HCA Parameters

The calibration of C_i parameters is well-documented for the HCA model [13,22,23]. Based on these, there are three possible methods for determining these parameters. In the first case, at least 11 drained cyclic triaxial tests shall be performed, which is quite time and work consuming. In the simplified procedure, C_e, C_p, C_Y and C_ampl are determined with correlations and C_Ni parameters are calibrated using a single cyclic drained triaxial test. In the simplest approximation procedure, all parameters are determined using correlations based on particle size distribution tests (d₅₀ and C_U).

Because no additional laboratory tests were performed for this research, the latter simplified procedure is presented in detail below. The simplified calibration procedure discussed in the papers of [22,24] was developed on clean subangular quartz sand and gravel soils with no fines content. The samples had a coefficient of uniformity C_U between 1.5 and 8.0 and a d₅₀ between 0.1 mm and 3.5 mm, and the maximum cycle number was N_max = 2·10⁶. The maximum number of cycles of 2 million is an outstanding upper value compared to similar tests found in the literature. The running time of a single cyclic triaxial test in this case is approx. 23 days, during which technical problems may occur even despite the greatest precautions, e.g. water leaks from the triaxial cell into the sample via diffusion [24]. The 2 million cycles used in the laboratory can therefore be considered the upper limit of technical reliability. However, road pavements can be subjected to a higher number of load cycles. Since behaviour of the soil is not known in the range greater than $2·{10}^6$ cycles, the assumptions corresponding to the range of N > $2·{10}^6$ should be treated cautiously and with appropriate criticism as the validity of the C_i parameters is uncertain.

The boundary conditions of the particle size distribution discussed above are not valid for the crushed stone base and the subbase courses, therefore the HCA parameters should be determined directly by laboratory testing of these layers [13]. The tests were performed on unbound crushed stone and the particle size distribution curve of the material Figure 2 is in accordance with the Colombian Road Specifications. The obtained HCA parameters are presented in Table 3. In the absence of literature data, it was assumed that the HCA parameters are also valid for the subbase course.

3. The Developed Layered Calculation Method

3.1. Basic Principles

The HCA model was originally implemented in the finite element program ABAQUS. Implicit calculations are generally time and work consuming, and implementation of the HCA model also requires high level of programming skills and mechanical knowledge. A number of technical difficulties are documented in [25] when implementing the HCA-model in the geotechnical finite element software Plaxis. In order to make the application of the HCA model easier for practicing engineers, simplified layered models were developed in MS Excel for shallow foundations under vertical cyclic loading and pile foundations under horizontal cyclic loading [26,27,28]. These layered models do not require finite element calculations, so the results can be quickly obtained without the need of any special software. The layered model is very useful for parametric studies or scheme design since the effect of different parameters on the settlement can be quickly investigated.

In this research, a simplified calculation model was developed in MS Excel on the analogy of the layered models, which can determine implicitly the strain amplitude (ε^ampl) in the layers (crushed stone base course, subbase course, subgrade) beneath flexible pavements based on the theory of elasticity, the accumulated permanent strain (ε^acc) explicitly using the HCA model, and the rutting depth under cyclic traffic loads.

The simplified calculation model works similarly to the layered models that are generally used for settlement calculations as the subsoil is divided into several sublayers. Then, stress distribution under the axle/wheel load is determined and the additional vertical load Δσ_1,i in each sublayer is calculated. The strain amplitude in the center line of each sublayer ε_i^ampl (mean of the values at the top and at the base of the sublayer) is determined from the mean normal stress, relative density and strain-dependent nonlinear stiffness. Then, the accumulated permanent strain ε_i^acc is calculated in each sublayer using the ε_i^ampl as input parameter, and compression of each sublayer is obtained by multiplying it with the thickness and the associated cyclic flow rule: s_1,i^acc = m₁ ∗ ε_i^acc d_i. Direction of the accumulation was determined following the MCC flow rule in the subgrade, and using the anisotropic flow rule of [21] for the base and subbase courses. The total settlement or rutting depth is then obtained as the sum of the compression of each sublayer Σs_1,i^acc. The principles of the layered model are summarized in Figure 1.

3.2. Calculation Steps

The steps of the implicit calculation part can be summarized as follows. Subscript “i” refers to the number of the sublayer, and subscript “j” refers to the load step:

• the vertical and horizontal stresses due to the own weight of the soil are determined:
• ${\sigma }_{1,i}=\mathsf{\Sigma }{\gamma }_i{h}_i$
• ${\sigma }_{3,i}={\sigma }_{1,i}{K}_{xy,0,i}$, where $K_{xy,0,i}$ is the coefficient of earth pressure at rest by also considering overconsolidation;
• the additional vertical load in each sublayer ${\mathsf{\Delta }\sigma }_{1,i}$ due to the wheel load is determined as the function of the geometry and stiffness of the sublayer;
• the additional horizontal load is determined from the additional vertical load ${\mathsf{\Delta }\sigma }_{3,i,j}={\mathsf{\Delta }\sigma }_{1,i,j}{K}_{xy,i,j}$, where $K_{xy,i,j}$ is the coefficient of earth pressure considering overconsolidation and the incremental loading;
• the additional vertical strain is calculated from the additional stresses, mean normal stress, relative density and strain-dependent stiffness: ${\mathsf{\Delta }ϵ}_{1,i,j}=\frac{{\mathsf{\Delta }\sigma }_{i,j}}{E_{i,j}(p_{i,j},e,{\gamma }_{i,j}^{ampl})}$. Because the incremental strain is the function of stiffness, and stiffness is the function of strain, iteration shall be used to determine the strain increment;
• total vertical strain and strain amplitude are calculated as the sum of increments in each sublayer: ${{2\epsilon }_{1,i}^{ampl}=\mathsf{\Delta }ϵ}_{1,i}=\mathsf{\Sigma }{ϵ}_{1,i,j}$;
• total horizontal strain and strain amplitude are calculated from the vertical strain and the ratio of horizontal strains $\vartheta$ in each sublayer: ${2\epsilon }_{3,i}^{ampl}=\mathsf{\Delta }{\epsilon }_{1,i}\vartheta$ = $\mathsf{\Sigma }{ϵ}_{1,i,j}\vartheta$;
• Then, the elastic strain amplitude ε_i^ampl is obtained using Equation (1).
• For the base and subbase courses, the above implicit calculation method is simplified to the extent that the stiffness is taken into account with the M_R resilient modulus, which is independent of the void ratio and strain level.

The explicit calculation part can be summarized as follows, where subscript “i” refers to the number of the sublayer, and subscript “k” refers to the packages of cycles:

11.

the mean normal stress and stress state ($p_i^{av},{\eta }_i^{av}$) is determined in the center line of sublayer „i” then factors independent of the number of cycles f_p, f_Y and ${\dot{f}}_N^B$ are calculated;

12.

in the next step, factors of Equation (4) that depend on the number of cycle are obtained: f_ampl from the actual ε_i^ampl, ${\dot{f}}_{N,k}^A$ from the actual g^A_k and f_e from the actual e void ratio;

13.

the accumulated strain rate ${\dot{\mathit{\epsilon }}}_{\mathit{i}}^{\mathit{a}\mathit{c}\mathit{c}}$ is calculated using Equation (4). In the case of drained conditions, Equation (2) becomes ${\dot{\epsilon }}_i^{acc}={\dot{\epsilon }}_i$;

14.

the strain increment due to $\mathsf{\Delta }N$ cycles is determined by $\mathsf{\Delta }{\epsilon }_i={\dot{\epsilon }}_i\mathsf{\Delta }N$, then ${\epsilon }_i(N+\mathsf{\Delta }N)={\epsilon }_i\left(N\right)+\mathsf{\Delta }{\epsilon }_i$;

15.

calculating the rate of the preloading variable ${\dot{g}}^A$ with the actual value of g^A_k;

16.

the actual value of g^A_k is updated due to $\mathsf{\Delta }N$ load cycles by $\mathsf{\Delta }{g_i}^A={\dot{g}}_i^A\mathsf{\Delta }N, \mathrm{then} g_i^A(N+\mathsf{\Delta }\mathrm{N})=g_i^A\left(N\right)+\mathsf{\Delta }{g}_i^A$;

17.

vertical strains are calculated using the flow rule ${\epsilon }_{1,i,k}^{acc}={\epsilon }_{i,k}^{acc}{m}_1$, then compression of each sublayer by $s_{1,k}^{acc}={\mathsf{\Sigma }\epsilon }_{1,i,k}^{acc}{h}_i$;

18.

volumetric strain increment is determined due to $\mathsf{\Delta }N$ cycles, then the new void ratio is obtained by ${\epsilon }_{v,i,k}^{acc}={\epsilon }_{i,k}^{acc}\left(m_1+2m_3\right)$, $e_{i,N+\mathsf{\Delta }N}=e_{i,N}-{\epsilon }_{v,i,k}^{acc}\left(1+e_{i,N}\right)$;

19.

the $e_{i,N+\mathsf{\Delta }N}$ characterizing a denser state results in a higher stiffness $E_{i,N+\mathsf{\Delta }N}(p,e,{\gamma }_{i,j}^{ampl},{\vartheta }_i)$, thus points d)–g) of the implicit calculation are repeated;

20.

steps b)–j) of the explicit calculation are repeated until the goal value of number of cycles.

4. Analyzed Pavement Structures and Material Models

4.1. Analyzed Pavement Structures

The goal of this research was to determine the accumulated permanent strains and rutting depths in the standardized pavement cross-section according to [1] shown in Table 1 due to repeated passes of a standard axle. The cross-sections are composed of 4 layers (sand and gravelly sand subgrade course, granular subbase course, crushed stone base course and asphalt course). The subgrade comprises quartz sand and gravelly sand of subangular particles; particle size distribution is shown in Figure 2, and it was assumed to have a relative compaction of R.C. = 93%. The subgrade is overlain by a subbase course of 20 cm thickness. The material of the subbase is discussed in detail in the thesis of [29], and its particle size distribution is presented in Figure 2. The subbase is overlain by an unbound well graded crushed stone base course of 20 cm thickness (Table 1); this material is described in detail in the works of [21,30]. Particle size distribution of the pavement layers is shown in Figure 2, whereas their properties are summarized in Table 2 The base course is overlain by an asphalt course of traffic load class-dependent thickness.

Figure 2. Particle size distribution of the pavement layers.

Figure 2. Particle size distribution of the pavement layers.

4.2. Materials and Their Properties

The asphalt course was considered as a homogeneous layer, which can be characterized by its traffic load class-dependent thickness, Young’s modulus and Poisson's ratio. Stiffness of the asphalt was taken to be constant (E_a = 4.000 MPa, υ_a = 0.35) according to the equivalent temperature [31]. Poisson’s ratio of the subbase and base courses was taken from [31]. The Poisson’s ratio in the subgrade is not constant either in the sublayers or during loading, and its value can be calculated based on the relationships of the theory of elasticity from the shear and constrained (oedometer) modulus.

It is assumed that the subbase course has a relative compaction of R.C. = 95%. Bearing capacity (E₂) modulus of the subbase material based on finite element calculation [32] is approx. 150 MPa, thus an E_2,sb of 60 MPa can be expected on top of the subbase. It is assumed that the crushed stone base course has a relative compaction of R.C. = 96%. Bearing capacity (E₂) modulus of the crushed stone material based on finite element calculation [32] is approx. 230 MPa, thus an E_2,b of 125 MPa can be expected on top of the base course, which complies with the international practice of flexible pavement construction. Because there was no Proctor test available for the subgrade, it was assumed that ρ_d^max corresponds to a void ratio of e = 0.95 e_min. Particle size distribution of the analyzed layers is shown in Figure 2, whereas their properties are summarized in Table 2.

The HCA parameters of the subgrade were calculated using the relationships of [18,33,34] from d₅₀ and C_U. The input parameters for the subbase and base courses were determined considering the tests of [21]; the only exception is the C_e parameter as it was calculated from the actual e_min of the subbase using the difference of e_min − C_e in the crushed stone base. The HCA parameters used in the analysis are summarized in Table 3. The notations in brackets for the subgrades correspond the same notation as the examined soils in [18].

Table 3. Geotechnical parameters of the modelled layers.

Layer	d50 (mm)	CU (-)	emin (-)	emax (-)	ρdmax (g/cm3)	e0 (-)	υ (-)
Subgrade 1	0.2	3.0	0.540	0.920	1.75	0.627	0.33–0.36
Subgrade 2	2.0	3.0	0.491	0.783	1.81	0.577	0.33–0.36
Subgrade 3	0.6	3.0	0.474	0.829	1.83	0.559	0.33–0.36
Subgrade 4	0.6	5.0	0.394	0.749	1.93	0.478	0.35–0.38
Subgrade 5	0.6	8.0	0.356	0.673	1.98	0.439	0.36–0.40
Subbase	2.0	11.9	0.364	0.513	2.06	0.340	0.40
Base	6.3	100.0	0.230	0.440	2.30	0.188	0.40

Table 3. Geotechnical parameters of the modelled layers.

Layer	d50 (mm)	CU (-)	emin (-)	emax (-)	ρdmax (g/cm3)	e0 (-)	υ (-)
Subgrade 1	0.2	3.0	0.540	0.920	1.75	0.627	0.33–0.36
Subgrade 2	2.0	3.0	0.491	0.783	1.81	0.577	0.33–0.36
Subgrade 3	0.6	3.0	0.474	0.829	1.83	0.559	0.33–0.36
Subgrade 4	0.6	5.0	0.394	0.749	1.93	0.478	0.35–0.38
Subgrade 5	0.6	8.0	0.356	0.673	1.98	0.439	0.36–0.40
Subbase	2.0	11.9	0.364	0.513	2.06	0.340	0.40
Base	6.3	100.0	0.230	0.440	2.30	0.188	0.40

4.3. Material Model of the Implicit Calculation

The HCA model is sensitive to the input of ε^ampl, so it is important to calculate it as accurately and carefully as possible. Generally, the hypoplastic material model with small strain stiffness is used for the implicit calculation part [10], however in [28] an isotropic elastic material model according to [35] is used in addition. In [28] it is noted that ε^ampl can be determined by basically any material model however the (elastic) deformations should be modelled as precisely as possible, therefore the material model used should take into account the nonlinearity of stiffness depending on relative density, mean normal stress and small strains.

In the implicit calculations of this research, the subgrade was modelled as a nonlinearly elastic material with a stiffness depending on the void ratio, the mean normal stress and the strain level. Small strain shear modulus G_max and constrained (oedometer) modulus of the subgrade were calculated using the relationships of [18,33,34] from d₅₀ and C_U. The degradation curve was determined according to the relationships presented by [33]. As stiffness (G, E, M, υ) depends on the mean normal stress and strain, it changes constantly and should be determined for each sublayer in each package of cycles in every load increment. Moreover, the actual strain increment and the corresponding actual stiffness should be determined by iteration for each load step.

For unbound granular base and subbase materials, the resilient modulus M_R can be calculated using the following nonlinear elastic relationship [36]:

$$M_R=k_1{\theta }^{k_2}$$

(5)

where:

M_R is the resilient modulus

Θ is the bulk stress

k₁ and k₂ are material constants.

Table 4 shows the typical range of k₁ and k₂ parameters. In the absence of laboratory tests, the k values should be determined by engineering judgement based on the drainage conditions and material characteristics.

Stiffness of the crushed stone base and subbase were determined using Equation (5). The "k" parameters for the crushed stone base were taken as the mean of the ranges given in Table 4 (k₁ = 5500 psi, k₂ = 0.6), whereas they were taken slightly above the mean for the subbase course (k₁ = 5000 psi, k₂ = 0.55). For the subbase, the slightly higher values are justified by the fact that the particle size distribution and type of the material (50% crushed stone) makes the subbase of higher quality than conventional subbase materials. It is also noted that the resilient moduli calculated by the above assumptions (M_R,sub~150 MPa, M_R,base~210 MPa show good agreement with the back-analysis of plate load tests presented in [32]

5. Implicit Calculation Model and Comparative Study

5.1. Vertical Stress Distribution

Traffic load is characterized by the standard axle load, which is taken as 10 tons. The axle load is distributed evenly on the two circular wheels, each with a diameter of 30 cm, so the additional vertical stress will be $\mathsf{\Delta }{\sigma }_1=707 kPa$ on top of the asphalt course. This stress is distributed proportionally to the thickness of the layers and their relative stiffnesses (E, v). The stress distribution and the additional stress $\mathsf{\Delta }{\sigma }_{1,i}$ in each sublayer are determined based on the theory of multi-layered elastic isotropic infinite half-space derived by Boussineq. Since the model includes the small strain behaviour of the subgrade, the layer’s stiffness significantly depends on the vertical stress increment. However, it is also true that the stress distribution in the layered elastic half-space is determined by the stiffness of the individual layers as an input parameter. The stiffer layers attract proportionally higher additional stress, as a result of which the stiffness of them starts to decrease due to the increasing strain level. Consequently, the final stress distribution should be determined by iteration.

Figure 3 shows the distribution of the additional vertical stress with depth for traffic load classes A-C-E-R. For traffic load classes B-D-K, the distribution curves fall between the adjacent load classes, so they are not displayed in the figure for better visibility. As the figure also shows, thickness of the asphalt course has a significant influence on the stress distribution. In traffic load class A, the crushed stone base, the subbase and the subgrade are subjected to approx. 5 times, 4 times and 2.5 times higher additional vertical load, respectively, than in traffic load class R. The curves converge well, and, although the Figure 3 is somewhat deceptive, even at the depth of 1.5 m the difference is still around 50% between traffic load class A and R.

5.2. Horizontal Strain Ratio

Vertical strain of the i^th sublayer due to additional vertical stress Δσ_1,i can be calculated form the three-dimensional Hooke’s law Equation (6).

$$∆{\sigma }_{1,i}=\frac{E_i}{(1+\upsilon )(1-2\upsilon )}\left[(1-\upsilon )·{\epsilon }_{1,i}+\upsilon ·{\epsilon }_{2,i}+\upsilon ·{\epsilon }_{3,i}\right]$$

(6)

where:

E is the Young’s modulus,

$\upsilon$ is the Poisson’s ration.

Let’s determine the lateral strains, ${\epsilon }_2$ and ${\epsilon }_3$, in the function of vertical strain by assuming that the strains in the two lateral directions are the same due to the axisymmetric conditions Equation (7). Combining Equation (6) and (7) and rearranging it to $\epsilon$, gives Equation (8).

$${\epsilon }_2={\epsilon }_3=\vartheta {\epsilon }_1$$

(7)

where: $\vartheta$ is the ratio of lateral and vertical strain

$${\mathsf{\Delta }\epsilon }_{1,i}=\frac{∆{\sigma }_{1,i}}{E_i}\frac{(1+\upsilon )(1-2\upsilon )}{1-\upsilon +2\upsilon \vartheta }$$

(8)

The dependence of ϑ factor on depth is discussed by [37] for shallow foundations under static load. This relationship did not provide satisfactory agreement with the comparative analysis, so a new relationship was developed for ϑ(z) in [32] for each pavement layer. [32] presents a finite element analysis for traffic load classes A-R where the load was gradually increased in steps of 2 tons, and the ϑ(z) relationship was obtained by curve fitting on the results using the least squares method. In this research, the ratio of horizontal strains is described by the ϑ(z) relationship presented in there.

5.3. Coefficient of Earth Pressure

An increase in the horizontal earth pressure reduces the permanent strains. The higher horizontal stress increases the mean normal stress and thus the stiffness of the soil E = f(p), so the value of f_ampl also decreases due to the decreasing cyclic strain amplitude ε^ampl. Although its effect is smaller, the increasing mean normal stress also reduces the degree of accumulation through the f_p function. Accordingly, the horizontal stresses should be taken into account as precisely as possible. The soils beneath road pavements are considered to be overconsolidated, thus conventional coefficients of earth pressure are not valid for them. The overconsolidation is caused by the traffic and compaction loads during construction. The stresses generated by these actions are many times higher than the stresses occurring during operation of the road. Consequently, OCR >> 1 and thus K_xy ≠ K₀^NC. [32] performed finite element analysis for the determination of K_xy(z) relationship for traffic load classes A-R by increasing the load in 2-ton steps; the relationship for the crushed stone base, subbase and subgrade was obtained by curve fitting on the results using the least squares method. The horizontal coefficient of earth pressure in this research is calculated using the functions and parameters presented in that paper.

5.4. Number of Load Increments

Since the applied material model is nonlinearly elastic, the result depends on the number of load steps. In the model, the load was divided into j equal parts and then it was applied in uniform steps of $∆{\sigma }_1/j$. A parametric study was performed on soils with different thickness and loading and the number of load steps was changed between j = 2-3-4-5-10-20-40-80-200-1000, and then the strain $∆{\epsilon }_1$ was determined due to $∆{\sigma }_1$ load. The relative error, i.e. the ratio of the strain calculated with the actual "j" and belonging to j = 1000 converges to 1.0 over time, so the calculated strains are almost the same, and from an engineering point of view, there is basically no difference between the results. The relative error was calculated in 3 different depths for 5 load increments and 10 load steps. The relative error obtained for depths of h = 0.5 m and h = 2.5 m is presented in Figure 4 as the function of the number of load steps. The relative error quickly converges to 1.0; the difference is already negligible at j = 40. The value of relative error ranged from 0.97 to 1.00 in the entire analysed spectrum for j = 40. Due to the above reasons, the loading step of j = 40 was used in all the calculations, which already yields sufficiently accurate results from an engineering point of view.

5.5. Number and Arrangement of Layers

After clarifying the implicit part of the calculation, an answer had to be given to the question of how to discretize the soil, i.e., how the number and distribution of sublayers affect the settlement. The goal was to determine the minimum number of sublayers, with which the settlement no longer changes significantly. A parametric study was performed for traffic load class C and compression was evaluated in each layer with an increasing number of sublayers.

In the first step discretization of the subgrade was analysed. The sublayers can be distributed evenly with depth, i.e., all sublayers have the same thickness. Since the input parameters and thus the strain amplitude changes rapidly at the top of the layer, this leads to an overestimation of the settlement. Discretization of the layers was therefore optimized in such a way that, although the number of sublayers remained the same, the sublayers were thinner at the top of the layer, whereas they were thicker at the bottom; from top to bottom: 0.1 m - 0.1 m - 0.1 m - 0.15 m - 0.25 m - 0.25 m - 0.25 m - 0.25 m - 0.25 m - 0.25 m - 0.5 m. Figure 5 shows the cumulative settlements with depth and the compression of each sublayer for both types of discretization. The cumulative settlement is almost 8% higher in the case of uniform discretization, and this develops in the top approx. 0.5 m zone of the layer. The figures also show that the cumulative settlement is smaller in the deeper layers in the case of uniform discretization. The reason behind is that the value of ε^ampl is very close to 10⁻⁵, below which cumulative strains don’t develop. Since the uniform discretization describes the deeper layers with greater accuracy, settlement starts to accumulate from the depth of approx. 1.95 m, compared to the 2.1 m depth of the optimized model. Nevertheless, the optimized model describes the settlement more accurately in the most important depths.

In the next step, the discretization was optimized to determine to the necessary and sufficient number of sublayers. The developing settlements are shown in Figure 6 as the function of the number of sublayers (1-2-4-7-11-22). Since the sublayers were made thinner in the relevant places, an accurate result could already be obtained with a relatively small number of sublayers (from i = 4) and the correlation converges quickly. Based on the results, the subgrade was divided into i = 11 sublayers in the subsequent calculations, as above this further increase in accuracy was not experienced. The additional vertical stress, ϑ and K_xy coefficient of earth pressure changes rapidly with depth in the crushed stone base and subbase courses, so a finer discretization is required than for the subgrade. The developing settlements are shown in Figure 6 as the function of the number of sublayers (1-2-4-8). It can be seen that there is no significant change starting from i = 4, so the subbase and base courses were divided into 4 sublayers in the model. Since the sublayer thickness of 5 cm is quite small and the input parameters change continuously with depth, a uniform discretization was applied in these courses.

6. Discussion of Layered Model Results

6.1. Implicit Calculation Results

Using the simplified layered method presented in Section 3, Section 4 and Section 5, ε₁ and ε₃ strains were determined for the traffic load classes A-R and the above-presented pavement layers under a standard axle load, and the results were compared with the finite element analysis results presented in [32]. Figure 7a,b show the vertical and horizontal strains in the base and subbase courses for traffic load class E. The ε₁ values adequately approximate the vertical strain obtained by the finite element analysis, except for the boundary of the subbase and base courses. Shape of the functions is also satisfactory however the simplified model predicts that the vertical strain decreases with depth in the subbase layer, whereas the finite element analysis shows that it increases slightly with depth Figure 7a. The minor horizontal strain almost perfectly matches with the results of the finite element analysis (Figure 7b). Figure 7c,d show the results for the subgrade, where good match between the results can be observed. The vertical strain is slightly overestimated in the top of the layer, whereas it is slightly underestimated in the lower portions. The horizontal strain is always underestimated except for the topmost part of the subgrade.

Although each strain component is of interest, the focus is given to ε^ampl as the most important input parameter of the HCA model. Figure 8 shows the calculated ε^ampl for traffic load classes A-C-E-R and compares it with the results of the comparative finite element analysis. As expected, the greater the thickness of the asphalt course, which increases with traffic load class, the better the wheel load is distributed, thus the additional vertical stress acting on the underlying layers becomes significantly less. Accordingly, the strain amplitudes also decrease gradually in the higher traffic load classes (stronger pavements), and an almost three-fold difference is observed in the strain amplitude between the classes A and R. The relationship between ε^ampl and depth and the shape of the functions are the same regardless of the traffic load class. As the figure shows, the results obtained by the simplified model agrees well with the results of the finite element analysis.

6.2. Comparison of Permanent Deformation Results between Layered and FE Model

Then, the rutting depth predicted by the HCA model in traffic load class C was analysed in the cases when the input parameters, such as p^av, η^av, ε^ampl are determined using the simplified layered procedure and when they are obtained from the finite element model. Since the selected finite element material model (HS-Small) does not take into account the increase in stiffness with increasing relative density, and thus would overestimate ε^ampl and ε^acc for high number of cycles, the ε^ampl was reduced in the first cycle according to Equation (9). It was assumed that the strains are still proportional to the component of the G_max relationship that contains the void ratio, and the strains at the beginning of the i^th packages of cycles were calculated from the ratio of the void ratio at the beginning of the i^th packages of cycles and the void ratio at the first cycle.

$${\epsilon }_i^{ampl}={\epsilon }_{Plaxis}^{ampl}/\frac{\frac{({a_G-e_{N=i.})}^2}{1+e_{N=i.}}}{\frac{({a_G-e_{N=1.})}^2}{1+e_{N=1.}}}$$

(9)

Figure 9 illustrates the obtained rutting depths with the number of cycles for traffic load class C. The results given by the simplified procedure and the finite element model are basically the same regardless of the soil type of the subgrade as the difference between the rutting depths of the two approaches at the end of the design working life is approx. 3-7%. The difference was of a similar magnitude for the other traffic load classes as well. Accordingly, it can be concluded that the simplified layered procedure presented above describes well the short-term behaviour of the subgrade under wheel load, therefore it is suitable to provide an appropriate input for the HCA model.

6.3. Effect of Traffic Load Class on Permanent Deformations

This section presents the residual rutting depth for subgrade 1 for each traffic load class as a function of the number of cycles, and comparison is made between the results obtained by using the simplified layered model to determine the input parameters and by using the finite element model (Figure 10). As the figure shows, the rutting depth calculated from the input data of the simplified layered model approximates well the result of the finite element model regardless of the traffic load class. At the end of the design working life, the difference between the calculated settlements is approx. 2-5%.

The pavement structure of the traffic load classes was devised in such a way that the stresses resulting from the repeated passing of axles with different numbers of cycles cause the same fatigue in the asphalt course over the entire design working life. Based on Figure 10, this is not true for the rutting depths as at the end of the design working life, they increase with the traffic load class. However, the rutting depths are smaller in the higher load classes for the same number of cycles. This is logical since the thicker asphalt course of the higher traffic load class distributes the wheel load better, so smaller stresses reach the underlying layers, and therefore ε^ampl will also be smaller. Consequently, the accumulated strains are smaller in the higher load classes for the same number of cycles. This means that in load class A at the end of the design working life, i.e. after 65 000 axle passes, the rutting depth is approx. four times the rutting depth of the class R at the same number of cycles. However, traffic load class R suffers approx. 4.6 times greater rutting depth by the end of the design working life as it is subjected to approx. 460 times more axle passes.

Figure 10. The relationship between rutting depth and number of cycles for subgrade “Soil 1” for traffic load classes “A”-“C”-“E”-“R” with the simplified procedure and with finite element modelling.

Figure 10. The relationship between rutting depth and number of cycles for subgrade “Soil 1” for traffic load classes “A”-“C”-“E”-“R” with the simplified procedure and with finite element modelling.

7. Conclusions

A calculation algorithm was developed in MS Excel, which implicitly calculates the throughout the lifetime of the pavement continuously changing elastic strains (ε^ampl) developing in the different layers of a flexible pavement system, and by using the HCA model explicitly determines the accumulated strains (ε^acc) and rutting depths s^acc under repeated axle load in the axis of the wheel. The advantage of the presented method is that the rutting depth can be predicted based solely on literature data and basic particle size distribution parameters (d₅₀, C_u, k₁, k_2, e_min, e_max), therefore it is suitable for parametric studies and quick preliminary designs, where the time and cost intensive laboratory tests and FE calculations are not available. For the sake of reliability it is however necessary, that for the further pavement design the HCA material parameters and the resilient behaviour of the soil layers are determined with time consuming cyclic triaxial laboratory tests to get more accurate material input parameters.

The model considers the change in the horizontal earth pressure resulting from the preloading of the pavement during construction and the typical horizontal strain conditions of road pavements. The actual behaviour of the subgrade is taken into account by the stiffness depending on the relative density, the mean normal stress and strain level, while the behaviour of the crushed stone base course and subbase course is taken into account by a nonlinearly elastic k-Θ model depending on the normal stress.

Due to the nonlinear behaviour, the strains shall be determined in incremental form $(\dot{\sigma }-\dot{\epsilon })$. It was shown that the necessary and sufficient number of load steps is i = 40. Discretization of the layers, i.e. required number and distribution of the sublayers, was also analysed. A parametric study was performed in the subgrade, which showed that more accurate results will be obtained if the sublayers are thinner at the top and thicker at the bottom instead of a uniform thickness (Figure 5). By changing the number of sublayers, a recommendation was made on the necessary and sufficient number of sublayers in the base, subbase and subgrade courses, which are 4, 4 and 11, respectively. Finer discretization does not lead to a significant increase in accuracy.

The elastic strains calculated by the developed simplified layered model for 5 different subgrades and 7 pavement structures were compared with the results of a reference finite element model [32]. It was shown that the developed simplified layered model can provide sufficiently accurate results for the strain amplitude ε^ampl based solely on literature data and basic particle size distribution parameters (d₅₀, C_u, k₁, k₂) without the need of finite element calculations (Figure 8). Depending on, if the input parameters (ε^ampl, p^av, η^av) of the explicit model from the developed simplified procedure or from the finite element calculations originate, the differences in the accumulated rutting depths for the traffic load class C with 3 different subgrades are approx. 3-7% (Figure 9). Similarly, the calculated rutting depths with the explicit model also show a good agreement (2-5%) regardless of the input parameters originate from the developed layered model or from the reference finite element model (Figure 10).

Rutting depth in each traffic load class and subgrade type was calculated and plotted on a diagram as a function of the number of cycles (Figure 10). It was shown that the rutting depth at the end of the design working life increases with increasing traffic load class. However, for the same number of cycles the rutting depths are smaller in the higher load classes.

It was shown that the permanent rutting depths due to traffic loads can be determined with sufficient accuracy using the developed simplified layered model for flexible pavements with different sections and subgrades. The future goal of this research is to perform additional parametric studies with the model and to explore which factors have the greatest influence on permanent settlement.