Rotation of Stresses in French Wheel Tracking Test

The main function of a pavement is to distribute the traffic-induced load over its different layers. While the flexible pavement design methods are based on a linear elastic calculation, the real behavior of the different layers is highly nonlinear and elastic. They can also, in some cases, be plastic and viscous. This research aims to develop a three-dimensional numerical model that is closely similar to the test FWTT conditions. The model will have a real geometry wheel footprint (rather than a rectangular shape). As a substitute for incremental loading, the wheel movement during its passage over the specimen will be simulated by a horizontal displacement. These important characteristics of the model represent the novelty and the major difference between the current research and previous studies. The current model, which is based on the finite elements method, uses Abaqus software and a viscoelastic constitutive model. The materials' viscoelastic properties have been described by the Prony series, also called the relaxation modulus, which is a function of time. This parameter can be defined in most computer-aided engineering (CAE) software. The procedure for calculating the Prony series from experimental data is explained. The results obtained agree with the stress signal amplitude, the stress rotation principal, and the total displacement rotation when the load approaches the node considered and located in the middle of the specimen.


Introduction
In the design of flexible pavements, the limitation of rutting is one of the most important phenomena that must be taken into consideration. It may be easy to measure the development of a rut, but it is more complex to predict its development. In most cases, knowledge of the material used alone is not enough; the environmental conditions and the distribution of stresses throughout the life of the pavement are also given consideration.
The result of the test carried out by Lekarp et al. (2000) [1] showed the presence of more permanent deformations in a hollow cylinder with shear stress than one without stress. Research by Bilodeau (2008) [2] also shows that the permanent deformations are doubled when a rotation of the stresses is applied during the triaxial test. This means rutting cannot be quantified with a conventional triaxial test. The cyclic loading tests carried out in the laboratory are only a simplified representation of the loading that solicits the unbound layers during road traffic, in which the applied loads are mobile. The movement of the wheel may cause a rotation of the main directions of the stresses. However, the triaxial test with repeated loading, which is mostly used for the study of permanent deformations, does not take into consideration this rotation of stresses and, until today, little is known about the parameters affecting permanent deformations. Gidel et al. (2001) [3] worked on improving the characterization of the behavior of granular materials based on TCR tests. The study, which was carried out both in the laboratory and on a field test board, takes into account the initial state of the material and its long-term evolution. Perret (2004) [4] highlighted the load factor by determining the related parameters that influence permanent deformations. Habiballah's (2003) [5] study, which was inspired by the works of   [6], relied on adaptation theory to predict the long-term behavior of a pavement structure. Its model was implemented in a code for finite element calculation. Cast3M, used in this case, showed the contribution of the different layers in the appearance of the rutting phenomenon. Dongmo-Engeland (2005) [7] highlighted this phenomenon in the cyclic tension/compression test by adding the temperature factor and its influence on the rutting with a proposal for a law of accumulation of permanent deformations. Bassem (2006) [8] studied the mechanisms responsible for the development of rutting by proposing a numerical model that is based on elastoplastic behavior for untreated seat layers and an elasto-viscoplastic model for surface layers. There was also a similar effect on Dawson (1996) [9] studied, where the temperature was taken into account. The study showed the influence of the type of tires and that of traffic.
Abd (2006) [10] experimented with rutting at the fatigue merry-go-round where he proposed a simplified rutting prediction method based on a finite element model (ORNI module of the CESAR-LCPC software). This was calculated in situ through measurements. On the one hand, what was studied was the mechanical behavior of granular materials by cyclic triaxial tests in the laboratory. Lv et al. (2019) [11] utilized three grades of modified asphalt mixture in the Hamburg using the Hamburg Wheel-Tracking Device test and Multiple Stress Creep Recovery test, with different content. The study found that optimal content exists for most additives. Lin et al. (2021) [12] study modified the Japanese pavement design method. The study considered the modification of the commonly used MEPDG model as well as the effect of the rotation of the principal stress axis. In Liu et al. (2021) [13] study, the microstructural changes in asphalt mixtures before and after wheel tracking test (WTT) were investigated using X-ray tomography (CT). A tomographybased aggregate tracking method is proposed to observe the spatial movement of coarse aggregates during the WTT. Fedakar et al. (2021) [14] investigated the deformation behavior of consolidated sand-clay mixtures by means of cyclic triaxial (CT) and hollow cylinder (CHC) tests. The test results clearly showed that the effect of a rotation of the principal stress must be taken into account to better estimate the deformation behavior of sand-clay mixtures under repetitive traffic loads. Alnedawi et al. (2019) [15] recent study, which uses the repeated load triaxial test (RLTT), examines the influence of loading frequencies on the rutting of unbound granular material. The test significantly showed that by increasing the loading frequency, the RLTT test duration could be significantly reduced.
The interaction between vehicles and pavements has been the subject of some studies in recent years, but these have not had a significant impact on the design of flexible pavements [16][17][18]. These studies have focused on the effects of traffic dynamics on accelerated pavement degradation due to axle aggressiveness or the effects of pavement roughness on truck tyres deterioration. Considering the huge budgets spent each year on road infrastructure maintenance and transport activities, there is no doubt about the importance of striking the right balance between the maximum axle load that should be allowed and the cost of restoring road infrastructure, all with the aim of maximizing the overall benefits to taxpayers.
One of the advantages of empirical mechanical methods is that they rely on one or more of the fundamental mechanical properties of pavements and soil layers to determine the stress state and hence predict pavement performance [19]. One of the most important properties studied in this context is the modulus of elasticity. The modulus of elasticity has many advantages over other index values such as the AASHTO layer coefficients, R-value, and CBR (California Bearing Ratio) in that it directly influences the analytical models used to predict the deformation state and consequently the stress state. Despite this important advantage, the use of the modulus of elasticity is also accompanied by some important problems. Firstly, the materials used in asphalt pavements are not elastic. Consequently, a surrogate for the modulus of elasticity -the modulus of resilience -is used to characterize the flexural strength that material in a given layer would possess under in-situ stress condition. Another problem with this method is the difficulty of measuring the resilient modulus accurately in the laboratory. While improvements in laboratory methods for testing the modulus are still expected, another method which uses non-destructive testing combined with back analysis is also more promising. In this method, surface deflection measurements are obtained non-destructively in the field and then evaluated using a mechanistic procedure (in this case a computerized back-calculation process) to determine the in-situ modulus of resilience of each layer. This process is particularly useful for the design of rehabilitation strategies and can also be applied to the design of a new pavement, provided that non-destructive test data are available along the intended route of the road.
Multiple mathematical relationships has been developed to correlate the stress state of pavement with its overall performance. The main transfer functions used in current flexible pavement design methods include: firstly, the maximum tensile elongation of hot mix surface layers under traffic loading as well as their resistance to fatigue cracking, and secondly, the compressive stresses caused by traffic loading at the top of pavement sub-bases and their resistance to rutting. These models are commonly derived from statistical correlations between pavement responses and observed performance in laboratory specimen tests and/or comprehensive road section experiments [20]. Transfer functions are the most important component of any mechanistic method and significant efforts are still required to develop more efficient and realistic models for predicting their performance.
In the category of simulation tests is the rutting test. The principle of this test is to subject a slab of the asphalt mixture to be tested to the forward and return motion of a surface loaded wheel [21]. During the test, the slab is kept as constant as possible in temperature. The change in rut depth is recorded as a function of the number of times the load is passed. The acceptance or not of the material tested is based on the rut depth obtained after a certain number of cycles. In particular, the LCPC rutting machine developed in France, standardized "NF EN 12697-22" is used in the study of the formulation of materials for pavement design.
The success of the approach lies in the good representativeness of the rutting test in relation to in-situ situation. In particular, it allows the material to be subjected to mechanical stresses that are relatively similar to those experienced in pavements (stress level, rotations, etc.) and makes it possible to adequately detect the ability of the material to resist or not to permanent deformations. However, the test alone, carried out according to these standards, does not make it possible to predict in detail the evolution of rutting in a pavement, in so far as it does not make it possible to take into account the good boundary conditions as well as the variability of the real conditions encountered on pavements, i.e. variations in temperature, wheel load, sweeping of the pavements, real stacking of the bituminous layers with their thicknesses, …etc. We will therefore no longer be interested in simulation tests, but in rheological tests in homogeneous stress and strain conditions.
There are several models that can be used for the stress state under rolling load. Most of these models are based on multilayer pavement elasticity theory and/or finite element analysis. Models based on the elasticity of multi-layered pavements are considered satisfactory for predicting the responses of flexible pavements to wheel load stresses and are also relatively easy to apply. However, these models are not able to predict the responses of pavements to environmental stresses (i.e. stresses associated with daily temperature and moisture variations, temperature gradients, etc.). In contrast, models based on finite element analysis are able to take into account both environmental and traffic constraints, but they are quite complex to implement and time consuming.
The cyclic loading tests carried out in the laboratory are only a simplified representation of the loading of materials during road traffic, during which the applied loads are mobile [22,23]. The displacement of the loads causes a rotation of the principal stress directions, however, the triaxial repeated loading tests (which is the most commonly used test for permanent deformations) does not take this stress rotation into account. Little is therefore known up today about this parameter affecting permanent deformations.

The Mesh Model
The French rutting test is modelled in 3D using the finite element code Abaqus [24]. The specimen is usually parallelepiped with dimensions of 500×180×100 mm. During the test, the wheel moves from one side to the other. The footprint of the wheel is simulated by a shell-type elliptical surface, with dimensions of 135 mm for the major axis and 80 mm for the minor axis (see Figure 2). As the problem is symmetrical along the wheel path, only half of the slab is simulated (Figures 3 and 4). The mesh of the slab is made up of 57600 C3D8R cubic elements and 64923 nodes, the footprint of the wheel that will carries the load is discretized into 890 shell elements of types S3 and S4R.

Loading
The loads can be simulated in the form of concentrated forces, pressures or displacements; the forces are applied to the nodes, the pressure loads are transformed into nodal forces and applied directly to the nodes in the loaded area. The research carried out in this context simulate the wheel in motion by moving the footprint of the tyre represented by a rectangular or circular surface. The movement of the wheel on the road is done incrementally by loading a set of elements, at each increment, adding a new line of loaded elements at the front end of the tyre and unloading a line of elements at the back end are added. In my opinion this process cannot model the real stresses induced by the tyre pressure, in the current, study the movement of the wheel is simulated by a horizontal displacement U3 of the elliptical surface between the two extremities of the asphalt concrete specimen, at a speed of 10 km/h. Its stroke is 365 mm for a duration of 0.123 s. The tyre footprint was modelled as a pressure load of 0.58 MPa, corresponding to a load of 5 kN on an elliptical area of 8523 mm 2 .
In order to be as accurate as possible in reproducing the stress states in the laboratory tests, we also sought to model the stress states due to the passage of a load in pavement structures, by taking into consideration the actual footprint of the wheel, the pressure on the specimen as well as the speed of movement of the wheel.

The Behavior of Material
Christensen (1982) [25] indicates that the time-dependent stress Sij(t) and strain εij(t) can be expressed by Boltzmann superposition integral, Equation 1. The concept of the linear viscoelastic model described in this paper is the generalized maxwell model represented by a set of elements connected in parallel, each element is composed of one spring for elastic behavior and one dashpot for the viscous behavior to which we add a spring in parallel Figure 5. The material of the specimen is assumed to be linear viscoelastic. (1)

Figure 5. Generalized Maxwell model
The simplest mechanical model, simulating linear viscoelastic behavior is a spring combined in series with a dashpot. The constitutive equation of this element is expressed as follows: After solving the above differential equations, the time-related expression for the stress was obtained.
The generalized model, shown in Figure 5 consisting of a spring and several Maxwell elements in parallel can be used. The relaxation modulus is then given by: where E0 is the equilibrium modulus, E ∞ is the instantaneous relaxation modulus, E i is the Maxwell element modulus, and i is the relaxation time. In the case of tension or compression stress: and in the case of shearing stress: The material is considered to be isotropic linear viscoelastic. For the specimen, the Poisson's ratio is considered as a constant, equal to 0.35. The values of the complex modulus were chosen for an asphalt mixture. They are given in Table  1 at a temperature of 20°C. For the elliptical surface it is elastic, its Young's modulus is E = 3500 MPa, and its poisson's ratio is =0.3.

Prony Series Fitting
The Prony series is utilized in this paper to describe the relaxation behavior of the bitumen mixture, which is subjected to the rolling load [26]. The Maxwell i th element has the shear relaxation modulus Gi, the viscosity ηi and the relaxation time τi=ηi/Gi. While the single spring that plays the role of the elastic behavior of the material after unloading it has stiffness G∞. Equation 7 shows the relationship between the relaxation modulus of the model and that of each model term, where n is the number of elements in the model. The more chains the more accurately the model describes the material. Fitting the experimental behavior of bitumen by the Prony series means minimizing the error between the predicted and the experimental values by adjusting the number n and coefficients. Instead, the Prony series can also be defined by G0, gi.
Since flexible pavements are subjected to time-depend temperatures and frequencies, the rheological properties of asphalt concrete are highly affected by its complex modulus. In the case of a viscoelastic material, the determination of the complex modulus G*(dynamic modulus) Equation 8 of asphalt mixtures is obtained by dynamic tests in which a sinusoidal load of frequency f is generally applied in shear, flexion or tension-compression. It is characterized by its norm |G*| and the phase angle, these two parameters depend on both the temperature and the frequency of loading. This modulus is composed of a real part G', also called the modulus of elasticity, which quantifies the stored elastic energy, and an imaginary part G'', called the loss modulus, which quantifies the energy dissipated by internal friction in the material. The expressions of G' and G'' In the frequency domain are given by the Equations 10 and 11 [27]. * = ′ + ′′ Using two terms (g1, g2, 1 , 2 and G∞) and three terms (g1, g2, g3, 1, 2 ,3 and G∞) of the Prony series, the fit of the model for G'(w) to the experimental data is given by the Equation 16 in the frequency domain and the equation of relaxation modulus Equation 15 in the time domain. Figures 6 and 7 show the relaxation and storage modulus variation curves, respectively. The Prony series parameters of the asphalt mixture used in this study are listed in Table 2.   Figure 8 shows the stress states of the specimen under rolling load at node N9539 in the middle of its top surface. It can be clearly seen that the evolution of the tensor components is similar to previous research, in particular Lekarp et al. (2000) and Di Benedetto & Corté (2005) [1,28]. The intensity of the vertical component S22 (Figure 9), responsible for the densification of the material, is much higher than the longitudinal and transverse components S33 and S11 ( Figures  10 and 11), respectively. Comparing these results with those found by the aforementioned authors [1,28], there is a difference between the pressure exerted by the wheel on the specimen and the amplitude of the vertical stress, for the node considered, this difference is greater than 15% compared to the pressure exerted by wheel, while for the present study they are almost identical. Figure 12 illustrates the state of stress of the sample under rolling load, in node N9551, located at mid-depth of the sample, the evolution of the components during the progression of the wheel, is similar to that of the previous node. The amplitude of the vertical stress is equal to 95% compared to the model of Di Benedetto and Corté (2005) [28], while in the present study, it is equal to 75% of the pressure applied on the surface. The perturbations due to surface friction are attenuated at depth. It can also be seen that the shear components S13 and S23 are closer in intensity, the maximum values of the three normal components S11, S22 and S33 are located at the plumb of the wheel passage, while at this same position, the shear stresses have a value of zero on the surface, S12 is slightly negative.  16 illustrate the principal stresses rotation evolution SI, SII and SIII at the node N9539 during the backward and forward movement of the wheel. The angles have been taken in absolute values, and the important variation of the rotation angle is clearly highlighted when the intensity is close to the maximum value, just before and after its passage on the considered node (Figures 17 and 18). The principal stress SII has been subjected to an important total rotation α22 which is close to 90°. This is thought to be the effect rolling that may have generated this great variation. It was also observed that the jump in the intensity has caused a rotation of 20° of the stress which is in the range 43 to 46 degrees (depending on the intensity of the stresses). The curves representing these rotations are angular peak curves. The three angles of rotation cancel when the wheel passes over the node and are therefore collinear to the reference frame (o, x1, x2, x3) ( Figure 19). Regarding the other two principal stresses SI and SII of respective rotation angles α11 and α33 their values increase rapidly as the wheel approaches the node taken into consideration and they have values of 49°. Furthermore, after its passage this rotation increases to 42°and the total rotation angle is then 91°. Figure 20 shows the principal stresses rotation direction. From the present study, it was found that the principal stresses rotation direction is directly related to the wheel movement direction on the test specimen, when the wheel moves from left to right, the rotation direction is clockwise, and when it moves from right to left, the rotation is counter-clockwise.  Figure 21 associates the displacement vector rotation α2 and the principal stress rotation SII. It establishes the relationship between these two magnitudes, which increase rapidly as the wheel approaches the node considered, in this case N9539, It should be noted that the rotation α22 of the displacement vector as the wheel approaches is three times greater than the rotation α2 of the displacement vector, and a variation of the angle of rotation of the constraint SII with a value of 45 degrees while in the same range. On the other hand, the displacement vector increases by 15 degrees, while the intensity of the constraint reaches maximum. They present the different components signals of the deformation tensor, when the load approaches the surface node. These graphs indicate that the shear component ε23 is substantially larger than the other components. The contact pressure induces, during its movement, shear forces at the surface, mainly at the edges of the wheel. This deformation is added to the vertical component ε22 will cause more rutting. This effect is essentially due to the shape of the load distribution on the specimen as well as its movement. However, concerning the node N9551, the vertical component is more important than the tangential one. This is essentially due to the shear attenuation in depth.

Conclusion
An accurate description of the behavior of bituminous mixtures is necessary to adequately predict and evaluate the time-dependent characteristics and the evolution of stresses and strains in the bituminous layers of flexible pavement. The viscoelastic characteristics of the materials implemented on the ACE can provide more accurate results, as they showcase a better approximation of the real behavior of the material. The values of the Prony series constants were obtained by a regression curve using the least squares method. A good fit was obtained for all measured values with 7 parameters, which is better than that of 5 parameters. The relaxation modulus of the analysed bituminous mixture showed a significant time dependence as the modulus values changed from the maximum value of G0 = 3898 MPa to the minimum value of G∞ = 297 MPa. The results of the model used show that the evolution of the signals representing the stress and strain states are accurate and are similar to reality when referring to previous studies. It was also noticed that the rotation of the vertical principal stress SII made a total angle of 90 degrees around the vertical axis by following the movement of the wheel during its passage. The same finding was made regarding the displacement vector U by making a sudden increase in its rotation angle around the vertical axis with a total value of 30 degrees. The other two main horizontal constraints, SI and SII, have also been subjected to a rotation around the vertical axis in the same direction, with a total angle of 90 degrees. The use of the proposed model allowed us to realistically evaluate the effects of the wheel movement when it is close to the vertical of the studied element. The elastoviscoplastic models of pavement structures and the effect of the interface between bituminous and granular layers will be dealt with in future work.

Data Availability Statement
The data presented in this study are available in the article.

Funding
This Project was supported by Civil Engineering Laboratory-Risks and Structures in Interactions (LGC-ROI), University of Batna 2, Algeria.

Conflicts of Interest
The authors declare no conflict of interest.