COMPARISON OF NONLINEAR ANALYSIS ALGORITHMS FOR TWO TYPICAL ASPHALT PAVEMENT ANALYSIS PROGRAMS

Two representative programs, MICH-PAVE and KENLAYER, are selected and compared to many key aspects of their analysis algorithms to achieve an in-depth understanding of the features of the Finite Element Method and elastic layered system theory in nonlinear material analysis of the structure of asphalt pavement. Furthermore, by conducting a case study, the impact of using different analysis methods on the calculation results is presented. Moreover, the feasibility of the equivalent resilient modulus obtained by the Finite Element Method is discussed. The results show that the difference among the nonlinear analysis algorithms used by the two software packages is mainly reflected in the determination of the initial resilient modulus, the stress


Introduction
The mechanistic-empirical analysis is the mainstream approach used to design of asphalt pavements. The first task in this approach is to obtain the real mechanical response of the asphalt pavement structure under the wheel load as accurately as possible. Boussinesq (1885), Burmister, Palmer, Barber, & Middlebrooks (1944), and Burmister (1945) provided analytical solutions for the mechanical response of a semiinfinite homogeneous body, a two-layered and a three-layered system under a circular uniform load, respectively. However, with an increase in the number of layers in the pavement structure, it is impossible and unrealistic to obtain a relatively closed theoretical solution for the mechanical response through cumbersome analytical derivation. Instead, it becomes necessary to introduce computational techniques and to exploit relevant programs to solve for the mechanical responses of pavement structure with increasing numbers of layers. Therefore, several computational programs have been developed based on the theory of elastic layered system, including BISAR (Bitumen Business Group, 1998), KENLAYER (Huang, 2004), CIRCLY (MINCAD Systems Pty. Ltd., 2009) and NonPAS (Ghanizadeh & Ziaie, 2015).
On the other hand, the Finite Element Method (FEM) has become an essential analytical tool in many different fields including pavement engineering, due to its wide applicability and excellent performance in solving many complex problems. Duncan, Monismith, & Wilson (1968) successfully applied FEM to analyse the mechanical response of flexible pavements for the first time. Since then, many FEM-based programs, such as axisymmetric Finite Element (FE) programs GT-PAVE (Tutumluer, 1995), MICH-PAVE (Harichandran & Baladi, 2000), ILLI-PAVE (Thompson, 1982) and 3D FE program EverStressFE (Davids & Clapp, 2009) have been developed to analyse the structure of asphalt pavement exclusively and applied to the calculation of mechanical responses of pavements.
Previous literature studies have shown that the theory of the elastic layered system and the FEM are the major approaches for analysing the mechanical response of asphalt pavement. In particular, Kruntcheva, Collop, & Thom (2005) studied the critical response and life of the pavement with different interlayer bonding condition by the elastic layered system and FEM. Park, Fernando, & Leidy (2005a) and Park, Martin, & Masad (2005b) measured the real tire contact stresses under three loading condition, which was used to predict the rutting and fatigue life of the pavement based on the FEM. Ziari & Khabiri (2007) analysed the impact of different interlayer bonding condition and the service life of the pavement by the programs based on the elastic layered system. Jiang, Zeng, Gao, Liu, & Qiu (2019) applied the program developed by the FEM to investigated the mechanical responses of flexible base asphalt pavements under six nonuniform contact stress distribution loads. However, each pavement structure layer was treated as linear elastic to simplify the analysis or improve the computational efficiency in most studies mentioned above. It has been demonstrated in former studies that both the granular materials and fine-grained soil materials exhibit significant nonlinear mechanical properties. As such, the modulus of these materials is strongly dependent on their stress state (Hicks & Monismith, 1971). They were simplifying the granular and fine-grained soil materials as linear elastic bodies fail to reflect their real behaviours and therefore yields erroneous results.
Hence, the nonlinear mechanical properties of granular and finegrained soil materials, instead of ideal linear elastic simplification, must be fully considered to get more accurate responses. You & Ling (2015) performed a FEM analysis of the mechanical response of asphalt pavements with either semi-rigid base layer or flexible base layer under three types of loads. Their research considered three types of material properties, including granular nonlinearity, subgrade nonlinearity, and a combination of them. The results showed that the structural performance is affected by the nonlinear properties of the materials for different levels under different structures and load conditions. Sun, Zhuo, & Liao (2014) used the FEM software ABAQUS to analyse the fatigue damage of asphalt pavement under repeated loads in consideration of the nonlinear properties of the materials. Kim & Tutumluer (2010) also utilised the same software to analyse the mechanical response of asphalt pavement with a flexible base layer under multi-wheel loads, considering the nonlinearities of both the subgrade soil and the granular material. In all these studies, the analysis was directly performed using a well-developed program without a deep understanding of the underlying nonlinear analysis algorithm. A detailed OF ROAD AND BRIDGE ENGINEERING 2 02 0/1 5 (4) discussion on the influence of different nonlinear analysis algorithms used in the theory of the elastic layered system and the FEM on the mechanical response of pavement structure is still lacking.
Among the special-purpose computer programs described above, KENLAYER and MICH-PAVE are two representative procedures, which are employed to perform nonlinear analysis. They are embedded with a nonlinear material model that fully considers the nonlinearities (i.e., the dependence of the resilient modulus on the stress state) of the granular material and the fine-grained soil. Due to their continuous update and improvement, both typical programs developed earlier still play significant roles in asphalt pavement structure nonlinear analysis field and show strong vitality even to this day. Many software products currently available have learned from their research and development experience on nonlinear analysis algorithms. Huang (2004) compared the solutions obtained by KENLAYER with MICH-PAVE for nonlinear layers to check KENLAYER correctness. However, comparative work focusing on nonlinear analysis algorithms still needs further comprehensive and in-depth. Hence, these two typical programs are selected to analyse and compare their nonlinear analysis algorithms in this paper. Besides, the impact of different analytical methods on the mechanical responses such as pavement surface deflection, radial strain, and vertical compressive strain at the top of the subgrade is also analysed in this study by using the asphalt pavement structure with the granular base as an example. These analyses facilitate a better understanding of nonlinear analysis algorithms and promote the application of the programs as mentioned above to the further research and design of asphalt pavement.

1.
Comparison of the algorithms used for material nonlinearity analysis in the programs

Overview
In both programs, nonlinear analysis is performed using an iterative method in which the stiffness of the structure is gradually changed to approach the final solution. The iteration processes used in the two programs are similar in general, but the determination of initial resilient modulus, stress correction, and convergence conditions are different. The similarities, differences, and the characteristics of the two nonlinear analysis algorithms used in the two programs are discussed in detail in the following sections.

Nonlinear constitutive relations
As shown in Figure 1, the nonlinear behaviour of granular materials is described by the K -θ model in both the KENLAYER and MICH-PAVE programs. This model was employed early in nonlinear analysis of pavement and laid a substantial foundation for followed research. The K -θ model is described as Eq. (1): where θ is the first invariant of the stress that is equal to the sum of three principal stresses or normal stresses, i.e., θ = σ 1 + σ 2 + σ 3 = σ x + σ y + σ z , MPa; M R is the resilient modulus, MPa; K 1 (in MPa) and K 2 are the material parameters obtained from experimental measurements, which are affected by the gradation, humidity, density, type and other factors of the material (Hicks & Monismith, 1971). The reference values of K 1 and K 2 for some typical materials provided by the two programs are listed in Table 1. The nonlinear behaviour of fine-grained soil materials is described by the classical bilinear model in both the KENLAYER and MICH-PAVE programs, as shown in Figure 2. In the KENLAYER program, the Eq. (2) of the bilinear model is given by In the MICH-PAVE program, the Eq. (3) of the bilinear model is given by  where K 1 (in kPa), K 2 (in kPa), K 3 , and K 4 − all material parameters obtained from experiments, which are affected by multiple factors including the size of the sample, compaction, density, and lateral confining pressure of the material (Thompson & Robnett, 1976); σ d is the deviatoric stress that is equal to the difference between the maximum and minimum principal stresses, i.e. σ d = σ 1 -σ 3 , kPa; and M R − the resilient modulus, kPa. As shown in Figure 2, K 1 and K 2 are defined in the opposite manner in the two programs. Also, in the KENLAYER program, the slope of the line is negative when σ d > K 2 . As such, the resilient modulus decreases with increasing deviatoric stress. However, the slope of the line for the same segment is positive in the MICH-PAVE program, which indicates that the resilient modulus increases with increasing deviatoric stress. This difference results in a significant deviation of the resilient modulus calculated at a point with high deviatoric stress among the two programs. Furthermore, to avoid unreasonable values of the resilient modulus, the user is required to provide the material parameters, i.e. K 2 , K 3 , K 4 , together with K 1 , the maximum and minimum resilient modulus based on the type of the subgrade, including stiff, medium, soft, and very soft when using the KENLAYER program. Based on these parameters, the relationship between the resilient modulus and the deviatoric stress was determined, as shown in Figure 3. When the modulus used for calculation is higher than the maximum value or lower than the minimum value input to the program, it resets to the maximum or minimum value, respectively. In the MICH-PAVE program, however, there is non-restriction on the maximum or minimum value of the modulus. The modulus is calculated directly from Eq. (3) based on the four material parameters.   (Thompson & Robnett, 1976) psi Note: 1 psi = 6.895 kPa 2 02 0/1 5 (4)

The gravity stress and the residual stress
When analysing a nonlinear behaviour that is highly dependent on stress, all the associated stresses, including the additional stress, gravity stress, and the residual stress, is considered in the process of calculating the modulus. Among these stresses, the gravity stress of the calculation point is computed as the accumulation of the layer thicknesses multiplied by the appropriate unit weights. σ g is denoted as the vertical gravity stress calculated at an arbitrary point in the layer of the pavement, then the lateral stress σ h induced by the material at the same point, is given by Eq. (4): where K 0 − the coefficient of earth pressure at rest, which generally ranges from 0.4 (untreated) to 3.0 (fully compacted). Apart from the gravity stress, the repeated loading and the compaction during the construction produce horizontal residual stress in the pavement structure. This stress is also known as locked-in stress. Although this stress is relatively low, it can partially balance the additional horizontal stress generated by the load (Tutumluer, 1995). Therefore, the residual stress should still be considered during nonlinear analysis. In both programs, the residual stress is taken into consideration by setting the value of K 0 to be slightly higher than the coefficient of earth pressure at rest.

The initial resilient modulus
The initial resilient modulus of each layer must be determined before the iterative calculation. For a linear layer, the resilient modulus is kept as a constant during the iteration process in both programs. However, the initial resilient modulus for a nonlinear layer is determined differently in the two cases. In the MICH-PAVE program, the user is not required to provide any initial resilient modulus for the nonlinear layer. Instead, the programs assume that the load applied on the pavement surface spread with a slope of 2:1 (vertical direction: radial direction). Based on this assumption, the initial resilient modulus of each nonlinear element is then calculated by substituting the gravity stress and additional stress induced by the assumption into the constitutive model. This resilient modulus is used as the initial value in the iterative calculation. In the KENLAYER program, however, the user must provide an initial resilient modulus for the nonlinear layer, which is then used as the initial value in the iterative calculation process. Hanyan Gu, Zhenkun Li, Yanjun Qiu

Comparison of Nonlinear Analysis Algorithms for Two Typical Asphalt Pavement Analysis Programs
The determination of the initial resilient modulus as the initial value in the iteration process has a significant impact on the number of iteration steps required during the nonlinear analysis. In the MICH-PAVE program, the initial resilient modulus is determined based on an embedded algorithm and is, therefore, more stable. Thus, the number of iteration steps required in the MICH-PAVE program gets less. On the contrary, the KENLAYER program allows for a higher degree of freedom and requires a reasonable initial input value to ensure high accuracy of the calculation result. Therefore, the program recommends that the users use the parameter K 1 in the constitutive model as the initial resilient modulus of the nonlinear layer.

Stress points
The MICH-PAVE program calculates the resilient modulus of each element based on the stress condition (including the gravity stress) at the centre of the element. This method considers the resilient modulus variation along the radial and vertical direction because of the stress conditions change.
The KENLAYER program, which was developed based on the theory of the elastic layered system assumes uniform mechanical properties in each layer. However, stress states at diverse locations are different in the same layer. Therefore, it is necessary to select a stress point to reflect the nonlinear properties of the layer fully. The resilient modulus calculated at this point (considering the gravity stress) is used to represent the modulus for the entire layer. The user requires to define a point on the pavement surface, the slope of load distribution, and a z coordinate, as shown in Figure 4, to determine the location of the stress point. The z coordinate of stress point is determined differently for the granular base and the subgrade. The z coordinate of the stress point in the granular base is determined based on the selected stress correction method, which is discussed in detail in Section 1.6. The stress point in subgrade is selected at 2.54 cm beneath the top will yield a lower resilient modulus than at 60.96 cm beneath the top. Moreover, the deflection at locations far away from the load is highly dependent on the resilient modulus of the subgrade. Therefore, yields significantly larger deflection at these locations using a stress point at 2.54 cm beneath the top of the subgrade. Therefore, the stress point in subgrade is selected at 60.96 cm beneath the top when calculating the deflection basin, but 2.54 cm beneath the top for all other conditions (Huang, 2004).
Furthermore, a point on the pavement surface and the slope of load distribution were determined based on the output type. If only the maximum mechanical response is required, then the stress point is selected beneath the centre of the wheel for a single wheel load, or beneath the centre of the gap between two wheels for a two-wheel load condition. However, if the average mechanical response under a single wheel load is required such as the deflection basin, then the point on the pavement surface is set at the edge of load, and the slope of load distribution is set as 0.5:1 (the vertical direction: the radial direction).

Stress correction
In the MICH-PAVE program, the Mohr-Coulomb failure theory is used to correct the stress of each element during each iteration. Considering the vertical stress as principal stress, the upper limit of the maximum principal stress σ 1 and the less limit of the minimum principal stress σ 3 of an element are determined by the following the Eq. (5): where σ v is the vertical stress, c is the cohesion of the material, and φ is the angle of internal friction. At the end of each iteration, the stress, calculated in each element, should not exceed the strength of the element. As such, the maximum principal stress σ 1 should be lower than σ 1,max , and the minimum principal stress σ 3 should be higher than σ 3,min . Also, the maximum principal stress σ 1 should also be lower than σ′ 1 , which is given by Eq. (6):   Figure 5 describes the stress correction process based on the Mohr-Coulomb failure theory. The circle C shown in the figure represents the Mohr circle associated with the stress state of the element before (after) the correction. Circle A and B represent the Mohr circle formed by the vertical stress σ v with the calculated limits of the minimum and maximum principal stresses, respectively. Figure 5a represents the stress state 1 before the correction, where the minimum principal stress σ 3 is lower than σ 3,min . In this case, the stress circle intersects with the straight line of the failure criterion. Therefore, the stress must be corrected by setting σ 3 = σ 3 ,min and σ 1 = σ′ 1 = σ y . Figure 5b shows the corrected stress states for those shown in Figure 5a. Figure 5c represents stress state 2 before the correction where both the minimum principal stress σ 3 and the maximum principal stress σ 1 satisfy the requirement but σ 1 ˃ σ′ 1 . In this case, the stress circle intersects with the Figure 5. Plots of the stress correction used in the MICH-PAVE program (Raad & Figueroa, 1980) a) the stress state 1 before the correction b) the stress state 1 after the correction c) the stress state 2 before the correction d) the stress state 2 after the correction straight line of the failure criterion. Therefore, stress is also corrected by setting σ 1 = σ′ 1 . The corrected stress state 2 is shown in Figure 5d. Given that the horizontal stress generated by the load at the stress point in the subgrade is usually positive and satisfies the safety criterion by the Mohr-Coulomb failure theory, it is unrequired that stress correction in the subgrade is required in the KENLAYER program  (Huang, 2004). For granular materials, there are three types of stress correction methods, as shown in Table 2. The desired correction method and the corresponding z coordinate of the stress point are determined based on the value of an input parameter PHI.
It has been shown in previous studies that the first method allows for the most accurate correction among all three stress correction methods for granular materials. The nonlinear behaviour of the mechanical properties in the base layer and the subgrade is much more significant for pavement structures with a thin asphalt layer. Moreover, using the first method better simulate the change of the resilient modulus in the vertical direction in the granular base layer and therefore yield more accurate results. As a result, this method should be chosen for a thin asphalt layer where the thickness of each sublayer should be set as 5.08 cm. For a pavement structure with a thick asphalt layer, the impact of the nonlinear material property is very subtle. Consequently, the less accurate second and third methods also are used for correcting the stress in this case.

Convergence condition of the iteration process
A single indicator controls the iteration process in the MICH-PAVE program. Notably, the iteration is stopped when the error in the resilient modulus is lower than an allowable value that is set as 0.001 by default and cannot be changed. However, the iteration process used in the KENLAYER program is controlled by two indicators, including an allowable error and the maximum iteration steps. The iteration process is stopped whenever one of these two criteria is satisfied. Also, both values can be customised by the user to control the accuracy and the required time for the calculation. By default, the allowable error and the maximum iteration step are set as 0.01 and 15, respectively, in the KENLAYER program. Given that the KENLAYER program offers a higher degree of freedom with many adjustable parameters, the following issues are encountered: • divergence of the resilient modulus during the iteration process; • the accuracy criterion cannot be satisfied at the maximum iteration step. When the resilient modulus is too low for the subgrade, this value diverges during the iteration process. In this case, the relaxation coefficient is modified in the program from a default value of 0.5 to 0.25 or 0.125. When the accuracy criterion cannot be satisfied at the maximum iteration step, the user modifies either the allowable maximum iteration steps or the allowable error value.

Equivalent resilient modulus applicable for linear analysis
The MICH-PAVE program also outputs the equivalent resilient modulus for the nonlinear layer. This value is the average of the resilient modulus of each element in the region with an assumed load distribution of 1:2 (the radial direction: the vertical direction). As shown in Figure 6, the equivalent resilient modulus of the layer 1, 2, and 3 is equivalent to the average resilient modulus of all the elements in regions ABGH, BCFG, and CDEF, respectively. The equivalent resilient modulus reflects the nonlinear mechanical properties of the layer to a certain extent in the linear elastic analysis.  (Harichandran, Yeh, & Baladi, 1990)

Description of the case
This paper picks the asphalt pavement structure with the granular base as an example to analyse using both programs to compare the effect of the nonlinear analysis algorithm on the mechanical response of the pavement structure in parallel. Specifically, the parameters in nonlinear calculations, including pavement structure combination, material models and relevant parameters, and single circular load, are identical to a previous study (Harichandran & Baladi, 2000), shown in Figure 7.
The mechanical properties of the granular base and the subgrade are described by the Kθ and the bilinear models, respectively. The parameters, used in MICH-PAVE program, during the modelling stage are the same as those listed in Figure 7. In the KENLAYER program, however, the values of K 1 and K 2 for subgrade are assigned oppositely compared to those listed in Figure 7. Besides, the thickness of the subgrade is set as infinite in the KENLAYER program.

Model construction
As shown in Figure 8, an axisymmetric model is constructed using the MICH-PAVE program for the analysis. The radius of the loading circle is a, and the length from the symmetric axis to the side boundary is 10a. The distance from the top to the bottom boundary is set as 127 cm. The mesh used in the model is the same as that used in the literature study (Harichandran & Baladi, 2000). Specifically, the entire model is divided into 156 rectangular elements. In the radial direction, the model is divided into four regions with each region covering a radial distance from 0-a, a-3a, 3a-6a, and 6a-10a, respectively. Each of the four radial regions is further divided into 4, 4, 3, and 2 elements, respectively. In the vertical direction, each of the three layers is divided into 4, 5, and 3 elements, respectively. The left boundary of the model is the ote: ρ -density, v -Poisson ratio, E -resilient modulus a -radius of the load, p -level of load concentration symmetric axis. A zero radial displacement constraint is imposed on the right boundary of the model. The top is set as a free surface, and the bottom boundary is set as a flexible boundary. Furthermore, the form of output in the program includes the vertical section and the horizontal section. Because the rectangular shape of the element implies that the stress and strain calculated via FE analysis are most accurate at the centre of each element, the vertical output section is taken at the centre of each element and origin in the radial direction. This method allows one to obtain the mechanical response at the centre of each element and the interface among neighbouring layers. Figure 9 shows the model constructed using the KENLAYER program considered the symmetric feature. The first stress correction method, with the highest accuracy, is used in this model. The thickness of each sublayer in the granular layer is set as 5.08 cm. The stress point is selected at the intermediate height of each sublayer. The radial coordinate of the point on the pavement surface is set as 13.59 cm. The slope of the load distribution is set as 0.5:1 (the radial direction: the vertical direction). When calculating the deflection basin, the stress The initial resilient modulus for the granular layer and the fine-grained soil layer is set to K 1 (62.05 MPa and 20.82 MPa, respectively) in their constitutive models. In addition, to compare the results with those obtained using the MICH-PAVE program in parallel, the same locations are used for the calculation.

Analysis and discussion of the calculation results
The results show that for this example, convergence is reached after 8 and 3 iterations in the KENLAYER and MICH-PAVE programs, respectively. The Finite Element Method is more advantageous than the theory of the elastic layered system in terms of the iteration speed. Besides, the algorithm used to obtain the initial value of the iteration in the MICH-PAVE program is found to be stable. This paper has examined the calculation data generated by both programs to compare the similarities and differences among the output results from the two programs in detail. It is important to note that the sign of the mechanical response is defined differently in the MICH-PAVE and KENLAYER programs. In the former, the output tensile stress and strain are positive, but compressive are negative. Furthermore, the deflection is negative along the direction of the increasing depth of the pavement. On the contrary, compressive stress and strain are positive, but tensile are negative in the latter. Moreover, the deflection is positive along the direction towards the depth of the pavement. When the results from the two programs are compared, the sign needs to be converted for one set results.
2.3.1. The resilient modulus at the end of the iteration Figure 10 shows the distribution of the resilient modulus in the granular base at the end of the iteration calculated by both programs. It is determined that when using the KENLAYER program, a maximum resilient modulus of 114.1 MPa for the granular base is obtained in the top sublayer. With increasing depth, the resilient modulus decreases until a minimum value of 101.8 MPa is reached at the bottom of the granular base,  Figure 10a. When using the MICH-PAVE program, a maximum resilient modulus of 132.0 MPa for the granular base is obtained for the top central element directly under load. The resilient modulus continues to decrease with increasing distance from the centre of the load. However, at the bottom of the granular base, the resilient modulus increases slightly at locations far away from the load centre. A minimum resilient modulus of 99.2 MPa is obtained at the rightmost element on the top of the granular base. The elements subjected to Mohr-Coulomb failure or tensile failure at the end of the iteration are concentrated near the wheel load (the regions enclosed by the red dashed line in Figure 10b). Figure 11 shows the distribution of the resilient modulus in the subgrade at the end of the iteration calculated using both programs. The entire subgrade in the KENLAYER model shares the same resilient modulus of 50.2 MPa because there exists only one stress point in this case, as shown in Figure 11a. In the results calculated using the MICH-PAVE program, however, a maximum resilient modulus of 54.8 MPa is obtained at the second rightmost element on the top of the subgrade. In addition, the resilient modulus tends to decrease when moving away from this element. A minimum resilient modulus of 49.0 MPa is obtained at the rightmost element on the bottom of the subgrade. Furthermore, no element in the subgrade is subjected to Mohr-Coulomb failure or tensile failure, as shown in Figure 11b.
It is pointed out in the above results that the nonlinear analysis based on the theory of the elastic layered system only consider the change in the resilient modulus with a varying thickness within the layer. In contrast, the FEM consider the variation of the resilient modulus induced by the change in the stress condition for both the radial and the depth direction, simultaneously. Therefore, the FEM provides a more accurate calculation. However, the calculation process for this approach is affected by multiple factors such as the element type, boundary condition, and the mesh. Besides, this calculation method is much more complex than that based on the theory of the elastic layered system. Consequently, one should be cautious when analysing the nonlinear problem using the FEM.
The equivalent resilient modulus for the granular base and the subgrade provided by MICH-PAVE is 116.98 MPa and 50.19 MPa, respectively. This paper input these values to the MICH-PAVE and KENLAYER programs to analyse the feasibility of applying the equivalent resilient modulus obtained using the FEM for linear elastic analysis. The same pavement structure and load parameters as those shown in Figure 7 are used for the analysis. The mesh and boundary conditions used in the MICH-APVE model are identical to those shown in Figure 8. Finally, the results obtained based on the linear and nonlinear analysis are compared and discussed in detail in the following sections. Figure 12 shows the pavement surface deflection curve. It is determined that the distributions of the two nonlinear curves are similar where the maximum deflection is located at the centre of the load, and the deflection decreases with the increase in the distance from the load centre. Also, the deflection curve is depressed as a basin shape in the region directly affected by the load. At locations far away from the load, the deflection is attenuated linearly with a minimal change in slope in the results obtained from KENLAYER program. However, the change in the deflection magnitude is nonlinear for the results obtained using MICH-PAVE where the slope varies substantially. Furthermore, a cross-over among the two curves is observed at 76 cm away from the load centre approximately. The above phenomenon is possibly caused by the fewer number of elements in the MICH-PAVE model, and a very coarse mesh at region far away from the load, which leads to insufficient calculation accuracy in there. The results obtained using the KENLAYER program exhibited a more significant deflection where the maximum value is 0.4288 mm. This value is 3.9% larger than the maximum deflection obtained via MICH-PAVE.

Pavement surface deflection
Also, the results obtained from the linear analysis using the equivalent resilient modulus well describe the behaviour reflected by the nonlinear analysis. The two deflection curves obtained from the linear analysis intersect at 23 cm away from the load centre approximately. The essence of nonlinear analysis is to determine the actual resilient modulus of a nonlinear material under an applied load by repeated linear calculations in an iterative process. The resulting resilient modulus then is used to perform linear analysis for the results. It is also concluded in Figure 12 that the results obtained from the linear analysis are all smaller than those obtained from the nonlinear analysis. In particular, the maximum deflection obtained for the KENLAYER and MICH-PAVE programs based on linear analysis is 7.8% and 1.9% smaller compared to those obtained from the nonlinear analysis, respectively. Figure 13 shows the contour map of the radial strain distribution, among which Figure 13a and radial strain is located at the bottom of the asphalt layer under the load centre. Besides, the pavement structure exhibits a compressive state near the regions that are directly affected by the load and are far away from the load. The distribution of the radial strain obtained using linear analysis, and the equivalent resilient modulus has similar features to that obtained via nonlinear analysis.  Figure 14 shows the distribution of radial strain on the bottom of the asphalt layer. As shown in the figure, the distribution of the two curves obtained from the nonlinear analysis is very similar to each other where the maximum radial strain is located at the load centre and the radial strain decrease drastically near the region directly affected by the load. With the increase in the distance from the load centre, the slope of the curve becomes gentler. Furthermore, the radial strain calculated by the KENLAYER results in a transformation from a tensile state to a compressive state at 41 cm away from the load centre approximately. Increasing the distance from the load centre yields a slight change in the compressive region.

Radial strain
Similarly, the radial strain calculated based on MICH-PAVE also results in a transformation from a tensile state to a compressive state at 37 cm away from the load centre approximately. With a further increase in the distance from the load centre, the compressive radial strain increases firstly and then decreases. Finally, the KENLAYER program provides a more considerable maximum tensile strain of 129.9⋅10 −6 , which is 16.0% larger than the maximum tensile strain obtained using the MICH-PAVE program.
The results of the radial strain obtained from linear analysis using the equivalent resilient modulus shares remarkably similar features with those obtained via the nonlinear analysis. There is only a subtle difference between the results obtained from the linear and nonlinear analysis. The differences between the maximum radial strains obtained from the linear and nonlinear analysis are 2.4⋅10 −6 and 0.3⋅10 −6 for the KENLAYER and MICH-PAVE program, respectively.  Figure 15 shows the distribution of the vertical compressive strain on the top of the subgrade. As shown in the figure, the distribution of the two curves obtained from the nonlinear analysis is very similar to each other. Specifically, both curves exhibit a compressive state on the top surface of the subgrade and a maximum vertical compressive strain at the centre of the load. The magnitude of the vertical compressive strain decreases with increasing distance from the load centre. The change of the vertical compressive strain is relatively small in the region directly affected by the load. In the regions away from the loading zone, the vertical compressive strain calculated based on the KENLAYER program tends to decrease linearly with a minor change in the slope. However, the vertical compressive strain calculated by the MICH-PAVE program decreases nonlinearly away from the loading zone during which the slope initially increases, then decreases, before finally becoming constant. The KENLAYER program results in the most considerable vertical compressive strain of 189.0⋅10 −6 , which is 70.2% larger than the maximum value calculated using MICH-PAVE. The above result is possible because the mechanical responses at the interface among the intermediate layers are obtained by linear extrapolation in the FEM. Therefore, the vertical compressive strain calculated at these locations is less accurate than that calculated at the centre of the element. Furthermore, the fewer number of elements in the MICH-PAVE model would also reduce the accuracy of the calculation.  The results obtained from the linear analysis based on the equivalent resilient modulus exhibit very similar features compared to the results obtained using nonlinear analysis. There is only a slight difference between the results obtained from the linear and nonlinear analysis. The differences between the maximum vertical compressive strain obtained from the linear and nonlinear analysis are 2.2⋅10 −6 and 0.3⋅10 −6 for KENLAYER and MICH-PAVE program, respectively.

Conclusions
1. The differences in the algorithms for nonlinear analysis in the case of KENLAYER and MICH-PAVE programs are mainly reflected in three aspects: • the initial resilient modulus of the nonlinear layer is calculated using the MICH-PAVE program automatically without any user input, but the value must be provided by the user in the case of the KENLAYER program; • the MICH-PAVE program only uses the Mohr-Coulomb failure criterion to correct the stress, but KENLAYER program offers three stress correction methods; • the convergence condition for the iteration process in MICH-PAVE program is controlled by only one indicator, which is the allowable error. In contrast, the convergence condition in the KENLAYER program is controlled by two indicators, including the allowable error and the maximum iteration steps. 2. The essence of nonlinear analysis is first to determine the actual resilient modulus of nonlinear material using an iterative method and then to use this modulus to perform linear analysis for the final results. The theory of elastic layered system has a slower iteration speed. It only considers the variation of the resilient modulus in the depth direction. In comparison, the FEM exhibits a faster iteration speed and considers the variation of the resilient modulus in both the radial and the depth directions simultaneously. Each of the two analysis methods exhibits its characteristics. The desired analysis methods and corresponding procedures need to be appropriately selected depending on the calculation conditions. 3. Using different nonlinear analysis methods have different levels of impact on multiple mechanical responses. The vertical compressive strain calculated at the top of the subgrade within the loading zone differs substantially for the two different programs. However, both programs yield quite similar results for the pavement surface deflection and the radial strain.
4. The equivalent resilient modulus obtained from the nonlinear analysis can be used in the linear analysis of the radial strain and the vertical compressive strain at the top of the subgrade. However, the nonlinear analysis and the linear analysis based on the equivalent resilient modulus yield considerably different results for the pavement surface deflection.