Riccardo Artoni
email: riccardo.artoni@ifsttar.fr

Patrick Richard
email: patrick.richard@ifsttar.fr

Average balance equations, scale dependence, and energy cascade for granular materials

Keywords: numbers: 45.70.-n, 47.57.Gc, 45.70.-n, 83.10.Ff

A new averaging method linking discrete to continuum variables of granular materials is developed and used to derive average balance equations. Its novelty lies in the choice of the decomposition between mean values and fluctuations of properties which takes into account the effect of gradients. Thanks to a local homogeneity hypothesis, whose validity is discussed, simplified balance equations are obtained. This original approach solves the problem of dependence of some variables on the size of the averaging domain obtained in previous approaches which can lead to huge relative errors (several hundred percentages). It also clearly separates affine and nonaffine fields in the balance equations. The resulting energy cascade picture is discussed, with a particular focus on unidirectional steady and fully developed flows for which it appears that the contact terms are dissipated locally unlike the kinetic terms which contribute to a nonlocal balance. Application of the method is demonstrated in the determination of the macroscopic properties such as volume fraction, velocity, stress, and energy of a simple shear flow, where the discrete results are generated by means of discrete particle simulation.

I. INTRODUCTION

The behavior of granular materials under deformation and flow is often studied by means of numerical methods that solve the equations of motion for every single body and account for the interactions between particles through appropriate models (see, e.g., Refs. [1][2][3]). These methods, commonly grouped under the name of discrete element methods (DEM) [START_REF] Smit | Understanding Molecular Simulation[END_REF][START_REF] Schäfer | [END_REF][6], can be extremely useful for simulating discrete media with the purpose of establishing micro-macro relations, and in this perspective they need consistent averaging procedures [7][8][9][10][11][12]. An averaging method computes continuum-like variables (e.g., a velocity field) starting from their discrete, particle-scale counterparts (e.g., particle velocities). When using such an approach, two issues have to be considered with utmost attention: the representativity of the average and the physical meaning of the obtained estimates. We can say that an average performed at a given point is representative if a sufficiently large domain in space and time can be defined, where particles share the same properties. Such a volume, usually called "representative volume element," may be defined with respect to only some variables.

An important step in the development of averaging techniques was made by Babic [7]. In that work the author developed a general framework for weighted space-time averages, applicable to a wide range of conditions, and giving a large set of self-consistent continuum balance equations. In Babic's method, the definition of a representative volume is intrinsically that of a zone where affine velocity fields are locally uniform, i.e., no gradients occur. The assumption that such a volume can be defined is called by Babic the "continuum assumption." Now, the lack of scale separation typical of granular flows (the scale of spatial variation of variables has the same order of magnitude of the particle size) implies that in a three-dimensional (3D) flow we cannot in practice define such a volume. As a consequence, the physical meaning of averages obtained using Babic's approach in a 3D flow is generally questionable. At least some average estimates will then probably depend on the size of the averaging volume [8,13].

We can see a simple example of this by considering a uniform shear flow; if all particles follow the mean flow with no frustration, we will tend to say that fluctuations are negligible. But this depends on the definition of fluctuation. If we calculate the fluctuational kinetic energy ε T [which is related to granular temperature [14]; see also Eq. [START_REF] Woods | The Thermodynamics of fluid Systems[END_REF]] by considering fluctuations as Babic does, we find ε T ∼ ( γD) 2 , where D is the averaging scale. Thus the fluctuations depend on the averaging scale which casts doubt on such a method. Also the kinetic contribution to the stress tensor will suffer the same problem. Thus, as we will show below, such an approach may lead to errors that can reach several hundred percentages.

In this work we focus on this problem, with a particular attention on the mechanical energy equation. We discuss both the issue of dependence on the size of the averaging domain and the energy cascade picture associated with the averaging method.

It is clear that the dependence of some average variables on the averaging domain size is intrinsic to Babic's definition of fluctuations. In the following we show how a different fluctuation decomposition taking into account gradients yields a new derivation of continuum balance equations. These balance equations have the merit of clearly separating affine and nonaffine fields and identifying the terms responsible for averaging scale dependence. With this decomposition it is possible to attribute a more sound physical meaning to the terms appearing in the final balance equations. We show how these balance equations may be greatly simplified by a "local homogeneity assumption," which is strongly related to the fluctuation decomposition itself and which yields less strict requirements than Babic's "continuum assumption." Applying this method to a simple class of flows, we find that two separate paths for fluctuating energy dissipation are implied: while kinetic stress power enters a nonlocal balance, the contact stress power seems to be dissipated locally. This result, as well as the robustness of the method, is verified with discrete element simulations of a simple shear flow.

II. A NEW DERIVATION OF BALANCE EQUATIONS A. Weighted average, fluctuation decomposition and local homogeneity assumption

Given a particle property ψ p (x p , t ), Babic defined the averaged property ψ at point x and time t as a weighted spacetime average:

ρ ψ = ∞ -∞ p w(x p -x, t -t)m p ψ p dt , (1)
where w(x p -x, t -t) = w p is a normalized weighting function, ρ the average density, and m p the mass of a grain. The fluctuation was simply defined by Babic with respect to the center of the averaging volume:

ψp = ψ p -ψ(x). (2)
This may seem harmless, but in fact-as shown in Ref. [8]-the dependence of some terms of the continuum averaged equations on the size of the averaging domain is strictly related to this choice for the fluctuation decomposition.

We propose to define the fluctuation not with respect to the average value at point x (as was done in Ref. [7]), but with respect to the particle center, ψp = ψ p -ψ(x p ).

(

In the following we will redevelop continuum-averaged equations in the spirit of Babic's approach, but using this new fluctuation decomposition.

In order to simplify the result, we will make a constitutive assumption. We suppose that a scale exists where the gradients of ψ are smooth and that if this scale is the averaging scale intrinsic to the weighting function w p , we can approximate ∇ ψ as a constant near the averaging point. We call this assumption the local homogeneity assumption. Together with this assumption we also assume a homogeneous distribution of particle centers near the averaging point and the decorrelation of positions and velocities. These three assumptions are strictly related.

While the absence of scale separation usually prevents identification of a scale where gradients are zero, such a first-order approximation is reasonable and-as we will see-very fruitful.

The local homogeneity assumption, involving gradients instead of variables, is less strict than Babic's continuum hypothesis and therefore more likely to be valid under a proper choice of the weighting function. Based on this assumption, we can therefore approximate ψ(x p ) by

ψ(x p ) = ψ(x) + (x p -x) • ∇ ψ(x). ( 4
)
Following Babic's approach, it is possible to develop a generic balance equation for the particle property ψ(x p ). Details are given in Appendix A; here we show the main findings, with a particular regard to the continuity, momentum, and translational energy equations. When performing the derivation, many new terms appear (with respect to Babic's treatment), but most of them vanish due to the local homogeneity hypothesis.

B. Continuum balances

The local homogeneity hypothesis leaves the continuity equation unchanged:

∂ ∂t ρ + ∇ • ρv = 0, (5)
where average density and velocity are defined by:

ρ = ∞ -∞ p w p m p dt , (6) ρv
= ∞ -∞ p w p m p v p dt . (7)
As for the momentum balance, we end up with the classical equation:

∂ ∂t ρv + ∇ • ρvv = ∇ • T + ρg, (8)
where the total stress tensor is composed by three contributions:

T = T c + T k + T γ .
The first two terms are the same as in Babic. The first contribution is due to contact forces (both long-lasting contacts and collisions) and is given by

T c = ∞ -∞ p q>p w F pq l pq f pq dt , (9)
where w F pq is a weight function related to the fraction of the branch vector joining the two particles p and q lying within the averaging function. Due to the new decomposition, Babic's kinetic stress tensor is decomposed in a truly kinetic part,

T k = - ∞ -∞ p w p m p ṽp ṽp dt , (10)
and a new term,

T γ = -ρ(D • ∇v)(D • ∇v), (11)
which contains all the dependence on averaging domain size through the vector D. It is useful to recall that the emergence of a part of the stress tensor depending on velocity gradients has nothing to do with constitutive relations but is the joint product of coarse graining and scale dependence. Due to their definition, both T γ and T k are symmetric. Evidently, in absence of velocity gradients, T γ will vanish; in that case, the kinetic stress tensor as defined by Babic will no longer depend on the coarse-graining length. The components of D, related to the characteristic size of the domain along each direction, are defined by the following equation:

ρD 2 i = ∞ -∞ p w p m p (x pi -x i ) 2 dt . ( 12
)
where x i is the i -th coordinate of the averaging point.

If the local homogeneity hypothesis holds, D i can be shown to scale with D m,i , the size of the averaging domain in the i direction (see Appendix B). The scaling will depend also on the shape of the weighting function.

Moreover, if the assumption of local homogeneity holds, T γ contains all the dependence on the coarse-graining width of the stress tensor while the kinetic stress tensor T k defined by means of the present fluctuation decomposition appears to be independent on averaging domain size.

The translational kinetic energy balance is only slightly more complicated than the one obtained by Babic. Once more, most of the new terms are negligible under the local homogeneity assumption.

The total translational kinetic energy is found to be the sum of three contributions, ρE T + ρE γ T + ρε T , which are respectively the translational kinetic energy of the mean flow, the translational kinetic energy related to the velocity gradients, and the truly fluctuational translational kinetic energy. These three variables are defined as

ρE T = ρ 1 2 v • v, (13)
ρE γ T = 1 2 ρ(D • ∇v) • (D • ∇v), (14)
ρε T = 1 2 ∞ -∞ p w p m p ṽp • ṽp dt . ( 15
)
The translational kinetic energy balance equation is derived as:

∂ ∂t (ρE T + ρE γ T + ρε T ) + ∇ • (ρE T + ρE γ T + ρε T )v = ∇ • (q k + q c ) +∇ • ((T k + T c + T γ ) • v) +ρg • v -γ k -γ γ , ( 16
)
where q k and q c are Babic's energy fluxes (see Appendix A).

The last two terms correspond to the rate of conversion of kinetic energy to other forms of energy (dissipation due to friction, collisions, plasticity, but also to storage and restitution of elastic energy). In the rigid limit they represent dissipation through friction and inelastic collisions and are therefore positive. The first one contains the contribution of fluctuations:

γ k = - ∞ -∞ p q>p f pq • ( ṽp -ṽq )w S pq dt , (17)
and the second one the rate of conversion of mechanical energy due to affine deformations:

γ γ = - ∞ -∞ p q>p f pq • ((x p -x q ) • ∇v)w S pq dt , (18)
where w S pq is the mean of w p and w q . Both terms are independent of the size of the averaging domain. We can see here again the benefit of our approach: On one hand the domain size-dependent terms are clearly identified in the equation (ρE γ T and T γ ); on the other hand, the contribution of affine deformations is separated from the contribution of the fluctuations.

In fact, by superficially interpreting Babic's approach, one could conclude that all the energy is dissipated by means of fluctuations. This is incorrect, since affine deformations do also dissipate energy; the reason for this misunderstanding is that affine deformations are considered as fluctuations in Babic's method.

We can remove the mean flow kinetic energy terms by manipulating the momentum equation in a classical way (for example, see Ref. [START_REF] Woods | The Thermodynamics of fluid Systems[END_REF]). The result is an equation expressing the balance for the fluctuating energy related to fluctuations and to affine deformations:

∂ ∂t (ρE γ T + ρε T ) + ∇ • (ρE γ T + ρε T )v = (T c + T γ + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ , ( 19
)
where † stands for transposition. Let us do some more manipulation. By looking at the term accounting for the rate of conversion of mechanical energy due to affine deformations, manipulation of Eq. ( 18), making use of l pq = x q -x p , allows us to rewrite γ γ in the following form:

γ γ = T c * † : ∇v, (20)
where

T c * = ∞ -∞ p q>p w S pq l pq f pq dt , (21)
is a new contact stress tensor, differing from the T c by the weight used. Indeed, as mentioned above, w S pq is the simple mean of w p and w q , whereas w F pq is related to the fraction of the branch vector l pq lying within the averaging region.

III. ONE DIRECTIONAL, STEADY, FULLY DEVELOPED FLOW

In this section we will apply our method to a simplified flow situation in order to gather some information on the effect of the new terms on force balance and on energy cascade.

Let us consider a one directional, steady, fully developed flow [defined by v = (v x , 0, 0), ∂ x = 0, ∂ t = 0]. In such a case, the continuity equation is identically zero. The momentum equation becomes:

∇ • (T c + T γ + T k ) + ρg = 0, ( 22
)
and the fluctuating energy equation

(T c + T γ + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ = 0. ( 23
)
According to the definition of T γ , it is easy to see that its only nonzero component is:

T γ xx = -ρ(D y ∂ y v x + D z ∂ z v x ) 2 , ( 24
)
and it is therefore straightforward to see that

(T γ ) † : ∇v = 0, (25)
that is, the contribution of affine deformations to the stress tensor does not perform any work. Moreover, it is obvious that T γ directly affects only the x component of the force balance. Given that the flow is fully developed, we can safely remove it out from the balances, obtaining:

∇ • (T c + T k ) + ρg = 0, (26)
(T c + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ = 0. ( 27
)
According to Eq. ( 20):

(T c ) † : ∇v -γ γ = (T c -T c * ) : ∇v. ( 28
)
It is clear that T c and T c * are not in principle the same, particularly when the averaging domain is smaller than a particle diameter. However, if the averaging domain spans more than one particle diameter, we can assume that T c ≈ T c * . Equation (27) becomes therefore:

T k : ∇v + ∇ • (q k + q c ) -γ k = 0. ( 29
)
This means that the work provided by the contact part of the stress tensor is converted locally by affine deformations, while only the kinetic part of the stress tensor enters into a nonlocal energy balance. This result has important consequences, since it shows that separate paths may exist for the dissipation of the work of kinetic and contact stresses. However, it does not cast doubt on the validity of hydrodynamic theories [START_REF] Bocquet | [END_REF][17][18][19] since our results do not deal with constitutive relations but suggest that a revision of the balances used for such theories could lead to a simplified description in some cases. Concerning the size of the averaging domain, it may seem strange to talk about domains smaller than a particle diameter, because it may seem obvious that such domains are inappropriate: Few particles are indeed contained in a snapshot of such a domain. However, we are always dealing with space-time averages, so even if the averaging domain is small, in many cases (e.g., slowly evolving or steady-state processes) the domain indeed contains a lot of particles since many of them traveled through it during their history without affecting local homogeneity. In other cases, periodic boundary conditions may help too, allowing use of small domains along some directions. That is one of the reasons why in the next section we will employ the averaging method to some DEM simulation results pertaining to this class of simplified flows.

IV. DEM SIMULATIONS A. Method & parameters

We use our own 2D implementation [2,20] of the classical discrete element method where Newton's equations of motion for a system of N "soft" disks are integrated. Such a technique, which is able to reproduce successfully the experimental results in many configurations (e.g., gravity-driven flows [1,2,21,22], sheared systems [3,23], granular materials close to jamming [24], silos [25], and rotating drums [26,27]), requires giving an explicit expression for the inter-particle forces. The discrete element method is classical and well known and can found in the aforementioned references. Therefore, we just present here the forces used in this work. For the normal force between two overlapping disks we use a standard linear spring-dashpot interaction model [6],

F n = k n δ n -γ n v n ,
where δ n is the normal overlap, k n is the spring constant, γ n the damping coefficient, and v n the normal relative velocity. The damping models the dissipation characteristic of granular materials. Likewise, the tangential force is modeled as a linear elastic and linear dissipative force in the tangential direction F t = k t δ t -γ t v t , where k t is the tangential spring constant, δ t the tangential overlap, γ t the tangential damping, and v t the tangential velocity at the contact point. The magnitude of δ t is truncated as necessary to satisfy the Coulomb law, |F t | ≤ µ |F n |, where µ is the grain-grain friction coefficient. The simulated system is two dimensional (Fig. 1). The granular material is a dense assembly of N dissipative disks of average diameter d and average mass m. A small polydispersity of ±20% is considered to prevent crystallization.

The granular material is submitted to a plane shear, without gravity, leading to a uniform stress distribution. The material is sheared between two parallel rough walls, separated by a distance H. One of the walls is fixed, while the other moves at the prescribed velocity V . The flow and transverse directions are respectively called x and y. Periodic boundary conditions are applied along the flow direction, and the length of the simulation box, L x , is set to 60 grain sizes. The wall roughness is made of disks sharing the characteristics of the flowing grains (same polydispersity and mechanical properties, no rotation). Their centers are equally spaced by a distance equal to the largest disk diameter 1.2d. Along the y axis, the position y = 0 corresponds to the center of the glued grains on the fixed wall. The normal stress applied on the moving wall, P 0 , is controlled. The vertical position of the moving wall is thus not fixed and, using the method described in Ref. [28], the height of the system H obeys Ḣ = (P 0 -P w )L x /g p , where g p is a viscous damping parameter and P w is the normal stress exerted by the grains on the moving wall.

The following values of the parameters are used: k n /P 0 = 10 5 , k t = 2k n /7, γ t = 0, and µ = 0.5 and g p = 100 √ k n m. The value of γ n is adjusted to obtain a normal restitution coefficient e n = 0.5 [START_REF] Schäfer | [END_REF]. The equations of motion for the translational and rotational degrees of freedom are integrated with a velocity-Verlet scheme with a time step 10 -4 m/P 0 .

The number of grains, N , is adjusted so the size of the system in the y direction, H, is roughly equal to 80d. The initial state of the system is a randomly diluted hexagonal lattice with disorder both on grain positions and velocities. The attainment of a steady state is verified by observing timeinvariant total kinetic energy and distance between the two walls.

B. Averaging procedure

For the simulation set up described above, we implemented an averaging procedure following the method developed in the earlier sections of the paper.

Several snapshots of the state of the system (positions, velocities, contact forces) at given times were captured. For each time a space average was performed; given that the system was periodic, averages were computed at different y values. A Heaviside step function was chosen as the weighting function. The only parameter of the weighting function was therefore the averaging window size. Some authors suggested [9,10] smoother weighting functions. Here we preferred a step function for its simplicity and to avoid smoothing the profiles too much. Given that the system was simple and periodic, we felt that a nonsmooth function was the best choice, so as not to lose information about layering, localization, and so on.

Given that the simulations displayed a stationary state, the time average was performed as a simple mean over the results obtained for each snapshot. Standard errors of the time averages were estimated by the blocking method [29] in order to take into account time correlation of the data.

Note that, in order to calculate the full set of variables, two steps were necessary: a first step to calculate average values of solid fraction, velocity, and contact stresses and a second step to calculate terms containing velocity fluctuations. Two steps are needed since velocity fluctuations are defined with respect to the mean field.

C. Results

In the following we will apply the method to a simulation with a velocity of the top wall V top = 100d P/m; as stated above, the flow is at constant pressure, so the vertical position of the top wall evolved freely to a stationary value, y top ≈ 83.9d.

It is convenient to characterize such a shear using the socalled inertial number I that compares the typical time scale of microscopic rearrangements with the typical time scale of macroscopic deformations: I = γd/ P/ρ S . Note that I corresponds to the square root of the Savage number [14] and is also called the Coulomb number [30]. Our system corresponds to a rapid flow with a global value of the inertial number corresponding to I > 1; such a flow regime was chosen because kinetic and contact stresses have the same order of magnitude, allowing us to test the results of the previous sections. A snapshot of the system at steady state is given in Fig. 1. All profiles shown in the following are calculated using an averaging window of 2d unless otherwise stated.

Figure 2(a) shows the velocity profile at steady state. As highlighted in the figure, the profile is clearly linear except from two zones (∼ 10d wide) near the bumpy walls, where the shear rate increases approaching the boundary.

The solid fraction profile [Fig. 2(b)] displays a central zone with a constant value, and a decreasing behavior approaching the walls. The state of the system, characterized by a lower compaction near the walls, can be appreciated also from Fig. 1. In fact, all the variables suggest that the system is in a simple shear condition (uniform shear, stresses, no energy flux) in the center and deviates from this state only in two narrow bands (∼ 10d wide) near the bumpy walls. Another feature of the system that can be appreciated in Fig. 2(b) is a slight asymmetry of the profiles near the walls. As it could be seen also in Fig. 1, in the shear zone near the upper wall the material is denser and less agitated than in the shear zone near the bottom wall. This slight difference is related to the asymmetry of boundary conditions: while the bottom wall is fixed, the upper one is free to move vertically in order to impose a constant normal stress. This yields slightly asymmetric profiles for all the variables. Figure 3(a) displays the fluctuating energy profile calculated with the new procedure developed in this work for an averaging window of 2d. To compare the two approaches (the present one and Babic's one) we report in Fig. 3(b) the quantity k -1, where k is the ratio of the fluctuating energy profile calculated using one of the methods with different averaging sizes to a reference profile obtained with the present method and an averaging size equal to D = 1d. It clearly shows that the Babic's method leads to an estimate of the fluctuating energy which increases with the domain size. On the contrary, the present method leads to an estimate which does not depend on the size of the averaging window. However, if the latter size is too large it influences the averages close to the walls due to the inhomogeneity of the shear in such an averaging window. To illustrate this point, we have reported in Fig. 3(b) results obtained with D = 5d. Note also that, using the present method, the fluctuating energy profiles for D = 1d and D = 2d are almost equal for any position y. Figure 4(a) displays the profiles of the components of the contact part of the stress tensor (here mainly due to collisions, since the flow is quite rapid). Again, a plateau is found at the center; when approaching the walls all the components decrease in absolute value.

At the same time [Fig. 4(b)], the components of the kinetic part of the stress tensor increase when approaching the walls. This is a well-known direct consequence of the decreasing of solid fraction. In our simulation, near the walls, kinetic stresses are larger than collisional ones.

The sum of the two parts gives the total stress tensor, shown in Fig. 4(c). As predicted by the stress balance, the (xy) and (yy) components of the total stress tensor are uniform along y. Concerning the third part of the stress tensor, T γ , it is clear that according to Eq. ( 24) its only nonzero component is the (xx) one. For the reference simulation it was verified that, apart for a small difference in a thin zone (∼ 2d) close to the walls (where slight deviations from local homogeneity appear due to the shape of the wall), D y ≈ D m,y / √ 12, confirming the scaling given above.

Figure 5 collects contact and kinetic flux profiles. All the components of these fluxes differ from zero only when approaching the wall; as regards the y component of the fluxes, fluxes are negative when the fluctuating energy decreases with y and positive otherwise. This agrees with the sign convention used for q in the derivation, and supports a Fourier-like expression for the flux dependence on fluctuating energy gradient. On the other hand, the horizontal components of the fluxes are not negligible near the walls. This contrasts with a linear transport theory since there is no gradient of fluctuating energy in the x direction. Figure 6 displays the estimates of the stress power related to contact and kinetic contributions and of the dissipative terms related to affine and non affine motions. Firstly, it can be seen that far from the walls, the algebraic relations hold

T c † : ∇v = γ γ , T k † : ∇v = γ k , (30)
therefore justifying the separate dissipation paths for kinetic and contact stress powers. Then, near the walls, the energy fluxes seem to influence only the kinetic part, therefore supporting the result obtained in the previous section: only kinetic terms enter in a nonlocal energy balance, while contact stress power is dissipated locally.

V. DISCUSSION

In the previous sections we developed an averaging procedure for granular materials which was based on a weighted average plus an original fluctuation decomposition. Through numerical simulation results we proved the strength of the approach which gives consistent predictions for kinetic terms, and-more importantly-gives a new picture of the energy cascade in granular flows, at least for simple flows. In this section we discuss the validity of such an approach. It is clear that the assumption of local homogeneity holds only when the solid fraction and the shear rate profiles do not display large gradients with respect to the averaging domain size. For the simulations presented above, two zones were identified near the walls where such gradients exist. The independence of the fluctuating energy profile on the size of the averaging window and the respect of the stress balance are two proofs that the assumption is still acceptable.

Figure 7 allows us to discuss another important point: the physical meaning of the variables obtained by the method, with a particular regard to solid fraction. In the figure, the solid fraction profile was computed by two other methods: estimating the fraction of (1) a horizontal line and (2) of a horizontal slice (1d wide) intersecting the particles. These two methods have a clear geometrical meaning, the second being a moving average of the first.

For the rapid flow simulation which was taken as the reference [Fig. 7(a)], where no layering occurs, the three methods agree very well. Figure7(b) displays the results of the three methods for a denser flow (V top = 10, y top ≈ 50,⇒ I ≈ 0.2). Please note that the quite ordered bumpy wall was explicitly prepared to display layering. In this case the estimates on lines strongly fluctuate near the wall, and the average on slices is a smoother profile which still fluctuates, while the estimate produced by the present averaging method seems to be intermediate between the two. Far from the wall results from the methods coincide. Therefore, in the presence of layering, the solid fraction calculated by use of the present method may not have the direct geometrical meaning that its name implies. This problem may be reduced by choosing a smoother weighting function (in this paper we chose on purpose a very abrupt weighting function); we also suspect that in three dimensions this rather pathological issue may be less important.

Another issue caused by layering is that the strong spatial fluctuation of variables may imply that some of the assumptions of the local homogeneity hypothesis do not hold: The distribution of particles around the averaging point is not in principle homogeneous. Fluctuations are evidently associated with strong local gradients of solid fraction and shear rate. Therefore, when strong layering occurs, care has to be taken in applying the present method due to the approximations therein.

The inhomogeneous distribution of particles inside the averaging domain yields also a deviation of the virtual center of mass X from the center of the domain. This is not, in principle, a problem, because the averaging procedure is related to a volume and not to a point; average variables may be harmlessly referred to the virtual center of mass instead of the simple geometric center of the averaging domain. In other words the relation X = x which seemed to be a consequence of the local homogeneity assumption can be taken as a definition of x, the point where averages are supposed to be related.

In addition, results from simulations at V top = 1 and V top = 10 (which show quite strong layering), not shown here for the sake of brevity, show that stress balance is respected and confirm locality of the dissipation of the contact stress power. However, if layering occurs, local homogeneity may be difficult to assume, and therefore the simplified balances obtained in this work have always to be critically analyzed when dealing with regions very close to the walls.

VI. CONCLUSIONS

In this work, we propose a new derivation of continuum balance equations for granular materials which is inspired by the work of Babic [7]. We introduce a new decomposition of particle scale velocities into mean and fluctuations, taking into account velocity gradients. A very reasonable local homogeneity hypothesis allows us to simplify the balances, ensuring at the same time representativity and physical meaning of the continuum estimates obtained.

Two important results are obtained. First, the method solves the problem of dependence of some variables (e.g., fluctuating energy, kinetic stresses) on the averaging domain size which is intrinsic to previous methods [7] by clearly identifying scale-dependent terms. Second, the development suggests a new physical picture for the energy cascade in granular flows: kinetic terms enter in a nonlocal energy balance, while contact terms (containing both collisions and long-lasting contacts) seem to be dissipated locally. Both results are verified against discrete element simulations of granular flows, with respect to which the method is shown at work. This point is particularly important with regards to the hope of obtaining a statistical mechanics for granular flows whose predictions can be faithfully tested in physical experiments and numerical simulations. Note also that the question of the size of the averaging window close to a boundary [11] is problematic since it is bounded by the distance between the center of a grain and that boundary.

A last note has to be added on the domain of application of the present method. Discrete element methods are sometimes divided into smooth and nonsmooth ones [START_REF] Jean | Micromechanics of Granular Materials[END_REF]. The present averaging method (as Babic's one) applies rigorously only to smooth methods, since the idea of smooth (differentiable) particle fields was implicit in the expression of the equation of motion. The adaptation of such a method to non smooth dynamics is a very important point which will be the subject of future works.

FIG. 1 :

 1 FIG.1: Snapshot of the simulated system. A 2D granular material is sheared between two bumpy walls. The lower wall is fixed, the upper one moves horizontally at a constant velocity Vtop = 100 and is submitted to a constant vertical stress P = 1.

FIG. 2 :FIG. 3 :

 23 FIG.2:(Color online) (a) x velocity and (b) solid fraction profiles along the y direction (Vtop = 100). The system appears to be in a simple shear condition (linear velocity profile, constant solid fraction) in the center of the cell. The size of the averaging window is 2d.

FIG. 4 :

 4 FIG. 4: (Color online) Components of the stress tensor [contact part (a), kinetic part (b), and total tensor (c)] along the y direction. The upper wall velocity is Vtop = 100 and the size of the averaging window is 2d.

FIG. 5 :

 5 FIG. 5: (Color online) Kinetic (a) and contact (b) fluxes of fluctuating energy along the y direction (Vtop = 100). The size of the averaging window is 2d.

FIG. 7 :

 7 FIG. 7: (Color online) Layering close to a bumpy wall: solid fraction computed using the present method with an averaging window of 1d, on lines and on slices (1d wide) for (a) Vtop = 100, and (b) Vtop = 10.

Acknowledgments

We are indebted to M. Y. Louge for fruitful discussions and to S. McNamara for a critical reading of the manuscript. The numerical simulations were carried out at the CCIPL (Centre de Calcul Intensif des Pays de la Loire) under the project MTEEGD.

Appendix A: Derivation of balance equations

The derivation follows closely that of Babic, except for the choice of the fluctuation decomposition. The evolution equation for a generic particle property ψ p is:

We multiply by w p , sum on p, and integrate on time:

The left hand side term can be developed into [7]:

As regards the right hand side of Eq. (A2), Babic also showed that the interaction term could be developed as:

where w F pq and w S pq are two weighting functions respectively defined as:

and

Now let us see the effect of the fluctuation decomposition on Eq. (A3). Decomposing the velocity field with the new decomposition leads to:

The divergence term on Eq. (A3) can be written as:

Decomposing also ψ p yields for the second term in Eq. (A8):

where X is the average center of mass of the averaging domain, the integral in the third term of Eq. (A8) becomes

Following the homogeneity assumption, we neglect correlations of fluctuations and of particle positions, and we consider that the average center of mass coincides with the averaging point. This means that terms under the divergence operator in Eq. (A8) simplify to

Developing the k component of the third term in Eq. (A12) [using the definition ∆x pi = (x pi -x i )] yields:

Assumption of local homogeneity also implies decorrelation of the components of particle position with respect to the averaging point. This means that:

) where D i is related to the scale of the averaging volume in the i direction. If D is the vector with components D i , the third term of Eq. (A12) is

This term renders explicit the scale dependence which was otherwise implicit in Babic's method. Equation (A3) is finally reduced to:

As regards to the effect of the fluctuation decomposition on the right hand side of Eq. (A2), the result depends on the particular property ψ p considered. Thus in the following we will analyze each property separately. If mass conservation is considered (ψ p = 1, P pq = 0, g p = 0), the right-hand-side term of Eq. ( A2) is zero and it is straightforward to see that the classical continuity equation is obtained [Eq. [START_REF] Schäfer | [END_REF].

If force balances are considered (ψ p = v p , P pq = f pq , g p = g), Eq. (A3) becomes

The last two terms can be considered as related to two components of an effective stress tensor as defined by Eqs. (11) and (10). On the other hand, the right-hand side of Eq. (A2) is not affected by the nature of the fluctuation decomposition. Substitution of the P pq = f pq and g p = g yields to the standard momentum equation, Eq. ( 8), with the third component of the effective stress tensor given by Eq. ( 9).

Let us derive now the translational kinetic energy equation. In this case the reference particle property is

The average translational kinetic energy is

and the fluctuation decomposition allows us to split this quantity into six terms, three of which disappear due to the local homogeneity assumption, yielding:

where, with respect to Babic's development, ρε T is likely to be scale independent, and a new term appears, which is the translational kinetic energy related to affine deformation within the average volume ρE γ T . The divergence term in Eq. (A3) corresponds, for the kinetic energy, to:

When developing this term, the fluctuation decomposition gives 27 terms, most of which can be deleted as in the previous development through (1) the postulate X = x, (2) the assumption of decorrelation among components of x p -x, and (3) the assumption of decorrelation between components of x p -x and velocity fluctuations. It can be shown that the irreducible terms are

where

Let us treat the terms coming from the right hand side term of Eq. (A2). First, we must identify the particle scale interaction terms. Taking the product of v p and Newton's equation,

Therefore, the particle interaction term is

and the particle source term is

Concerning the source term, local homogeneity implies:

(A27) Developing the interaction terms yields:

Using the fluctuation decomposition, the relation f pq = -f qp , and the assumption of local homogeneity this becomes:

which correspond respectively to:

where the first two terms are the conversion rates of translational kinetic energy into other forms of energy respectively due to velocity fluctuations (-γ k ) and to affine deformations (-γ γ ); q c can be seen as a translational kinetic energy flux related to contact forces. The translational kinetic energy balance equation is therefore given by Eq. ( 16), which can be further simplified to Eq. ( 19) when subtracting the mean flow energy balance.

Appendix B: Derivation of a scaling for D

In the following, a simplified scaling for the vector D which enters in the definition of T γ is derived using the local homogeneity hypothesis. Let us assume equal-sized spheres, and a step function on space and time as a weighting function. Let us consider variables computed at a point (x, t) in space and time. The average density is therefore:

where w p = 1/(V T ) for particles residing in the averaging volume V during the period (t -T /2, t + T /2) and zero otherwise. If we call N p the time average of the number of particles lying in the averaging volume, it is clear that ρ = m p N p /V . Developing the definition of the vector D, we find:

where x pi is the i component of the particle position. The local homogeneity hypothesis says that X = x and that x pi is uniformly distributed around x i . This yields for large N p T to:

Therefore if the local homogeneity holds, D i is likely to scale on the size of the averaging domain in the i direction. The proportionality will then depend on the choice of the weighting function, on polydispersity, on shape, and so on.


Juliette Blanc
email: juliette.blanc@ifsttar.fr

-Thomas Gabet
email: thomas.gabet@ifsttar.fr

Pierre Hornych
email: pierre.hornych@ifsttar.fr

- Jean-Michel Piau*-Hervé
email: jean-michel.piau@ifsttar.fr

D I Benedetto
email: herve.di-benedetto@entpe.fr

Cyclic triaxial tests on bituminous mixtures

Keywords:

In the framework of a study on rutting of bituminous layers, thermo-controlled triaxial tests under constant confining pressure on homogenous asphalt concrete specimen have been performed. These tests are cyclic sinusoidal tests, with the aim to reproduce in some way the loading and unloading phases due to traffic on roads. The influences of confinement and load frequency have been studied. Tests have highlighted that confinement notably decreases the strain rate as well as the level of strain reached at the end of the test. The tests did not allow to clearly observe the influence of frequency on the response of the material, due to problems of repeatability. In order to avoid such problems, another test procedure was adopted. It consisted in changing either confinement or frequency during the test, with the purpose to see an evolution of strain rate, and to establish correlations. These tests have clearly shown the influence of cyclic loading and the relative influence of loading frequency on the response of asphalt concrete. The main objective of this triaxial program is to develop a numerical model , based on rheological laws, able to predict the behaviour of bituminous layers.

Introduction

Mechanical damage of bituminous pavement layers is due to two main mechanisms: rutting and cracking. This damage is mainly induced by traffic, and especially traffic of heavy and slow vehicles. Our study focuses on rutting, which is characterized by the development of longitudinal ruts due to permanent deformations of pavement layers, which affect safety and comfort of road users. Two types of rutting can be encountered on road structures: large radius rutting generally related to permanent deformation of granular base layers or subgrades, and small radius rutting due to permanent deformations appearing in the bituminous layers, mainly in the wheel path and close to it. Small radius rutting is the subject of this work.

When studying permanent deformations (rutting) appearing in bituminous layers, three main difficulties are encountered. The first problem is the complexity of loading, which is dynamic, local and characterized by a complex geometry and a non-homogeneous distribution pressure. The second difficulty is the complexity of the rheological behaviour of bituminous mixtures, which have a heterogeneous structure, made of a viscous material, the bitumen, and a granular skeleton. From a rheological point of view, bituminous mixtures present characteristics of viscoelastic, viscoplastic, thermosuceptible materials, and their behaviour is also dependent on the triaxial loading path and cyclic loading, like for granular materials. The third difficulty is related to the scales encountered, in terms of cyclic strains (very small), time and number of cycles (very large), which complicate the study of bituminous concrete in the laboratory.

For studying the rutting behaviour of bituminous mixtures in the laboratory, determining representative triaxial stress paths is a key point. Loading paths generated in pavements by moving wheel loads are difficult to reproduce in the laboratory, as shown by [START_REF] Perraton | Rutting of bituminous mixtures : wheel tracking tests campaign analysis[END_REF], and [START_REF] Gabet | French wheel tracking round robin test on a polymer modified bitumen mixture[END_REF]. They induce cyclic, three-dimensional loading paths, with stress rotations and evolve with time. Due to the difficulty to simulate accurately such loadings with a typical press, like a thermo-controlled triaxial apparatus, more empirical, simulative tests like the French Wheel tracking test are used for designing bituminous mixes. Although wheel tracking tests apply more realistic loading conditions, they do not allow to determine the rheological properties of the material tested, due to the complex, non homogeneous, stress conditions in the specimen. Thus, triaxial tests, on homogenous specimen are more appropriate for a rheological approach.

The most common method for studying the rheological behaviour of bituminous mixtures consists in using thermo-controlled uniaxial tests (monotonic or cyclic). They allow determining the viscoelastic properties of bituminous mixtures, but in terms of rutting, several authors have emphasized the importance of considering triaxial behaviour with confining pressure [START_REF] Visscher | The permanent deformation law of asphalt mixtures: investigation of the effect of mix composition and material properties[END_REF][START_REF] Ebels | Determination of materials properties of bitumen stabilised materials using tri-axial testing[END_REF][START_REF] Taherkhani | Characterization of repeated creep recovery behaviour of asphaltic mixtures[END_REF][START_REF] Taherkhani | Compressive uniaxial and triaxial testing of asphaltic mixtures[END_REF].

Cyclic triaxial tests on bituminous mixtures 3

At IFSTTAR, a triaxial apparatus has been developed for studying bituminous mixes [START_REF] Sohm | Cyclic triaxial apparatus for the study of permanent deformations of bituminous mixtures[END_REF]. It is able to perform a large range of tests: monotonic creep tests, complex modulus tests and cyclic permanent deformation tests at small and large strain levels (from 50 µstrain up to 5 %) and at different temperatures (from 5 to 60°C), different frequencies (up to 10 Hz) and different levels of confining pressure. First tests performed with this apparatus have been presented by [START_REF] Sohm | Creep tests on bituminous mixtures and modelling[END_REF]. This first study focused on the monotonic creep behaviour of bituminous mixtures, in order to examine the effect of different levels of vertical stress, confining pressure and temperature. These tests have shown, in particular, that the level of confinement has a very significant effect on the response of asphalt concretes. When confinement increases, the strain rate decreases and the volumetric behaviour of bituminous mixtures is modified. For a given level of axial strain, volumetric strain is dilatant (volume increase) at low confining pressures and may become contractant for high confining pressures. It has also been shown that increasing temperature during the creep test leads to higher strain rates, as it is the case in the viscoelastic domain.

Another interesting result of [START_REF] Sohm | Creep tests on bituminous mixtures and modelling[END_REF] was that the time temperature superposition principle (TTSP), usually considered as valid in the linear domain for most asphalt concretes, could be extended to the non linear domain even under confinement. This result confirmed previous results presented in the literature, but without confinement [START_REF] Di Benedetto | General overview of the timetemperature superposition principle for materials containing bituminous binder[END_REF][START_REF] Di Benedetto | Time-temperature superposition principle for bituminous mixtures: three dimensional approach and extension in the non linear domain[END_REF] and [START_REF] Nguyen | Time-Temperature Superposition Principle for bituminous mixtures[END_REF]), or in the very high frequency domain [START_REF] Sohm | Cyclic triaxial apparatus for the study of permanent deformations of bituminous mixtures[END_REF]. Finally a viscoelastic-viscoplastic model has been developed in order to simulate the experimental response. It shows both good kinetic and volumetric responses, and in particular reproduces well the dilatant/contractant aspects, despite some limitations for high deviatoric stresses without confinement. This paper presents results of cyclic sinusoidal triaxial tests, which simulate the stress paths generated in pavements by moving wheel loads. To reproduce such moving loads, cyclic sinusoidal triaxial tests with cyclic confining pressure ( [START_REF] Brown | Essais triaxiaux sur enrobés bitumineux et chargement répété ou en fluage[END_REF] and [START_REF] Brown | The measurement of recoverable and irrecoverable deformations in the repeated load triaxial test[END_REF]) are considered most representative. However, in this work, for technical limitations of the test apparatus, we focus on cyclic triaxial tests with constant confining pressure. Influence of the confining pressure and loading frequency on the evolution of permanent strains and on viscoelastic characteristics of a bituminous mixture are studied.

A first tests program (T1) was performed to evaluate the repeatability of the results and to evaluate the influence of the confining pressure and frequency. In a second test program (T2), the test procedure was modified, to overcome problems of repeatability. With the second procedure, different load sequences, with variable confinement and/or frequency, were applied successively on a single specimen, to avoid the scatter due to the use of several different specimens. However, with this second procedure, tests interpretation is more complex, due to stress history effects (each stress sequence can possibly modify the subsequent results).

BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLOR AND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

Conventions and notations

In general, the mean stress p and the deviatoric stress q are defined as:

( ) ( ) 2 3 3 2 p tr q tr S  = S   =   [1] and [2]
Where S is the deviatoric part of S, the stress tensor, defined by ( )

1 3 S tr I = ∑ - ∑ .
In the case of the triaxial test, the mean stress p and the deviatoric stress q are: 2 3

H V V H p q σ σ σ σ +  =    = -  [3] and [4]
Where H σ is the horizontal stress and V σ is the vertical stress.

In this paper, the soil mechanics sign convention is used (positive stress/strain in compression/contraction).

In the particular case of the tests presented in this paper, 0

V H σ σ ≥ ≥ so : V H q σ σ = - [5]
For the following, the axial or vertical stress is named σ x and the radial or horizontal stress is named σ r .

During a triaxial cyclic test, a cyclic sinusoidal deviatoric stress q and a constant confining pressure (or radial stress) σ r are applied to the specimen at the start of the test, and the evolutions of the axial and radial strains, ε x and ε r , are measured.

Volumetric strains and deviatoric strains are deduced from axial and radial strains, using equations 6 and 7.

( )

2 2 3 V x r d x r ε ε ε ε ε ε = +    = -   [6] and [7]
Cyclic triaxial tests on bituminous mixtures 5

Calculation of stress paths in pavement layers and determination of triaxial test conditions

To define representative stress paths for the triaxial tests, pavement structure calculations have been performed with the pavement design software Alize LCPC (http://www.lcpc.fr/en/produits/alize/index.dml). Computations are based on a linear elastic model. The pavement layers are supposed isotropic, and infinite in the horizontal plane. Additional calculations have also been performed using a viscoelastic model, (Viscoroute software, [START_REF] Chabot | Viscoroute 2.0: a tool for the simulation of moving load effects on asphalt pavement[END_REF][START_REF] Chupin | Influence of sliding interfaces on the response of a visco-elastic multilayered medium under a moving load[END_REF]. However, differences between the results were not significant, in terms of stress paths induced in the pavement. For brevity, only results of elastic simulations will be presented. The simulations have been performed for a pavement structure consisting of 24 cm of dense-graded asphalt mix (DGAM) and 16 cm of road base asphalt (RBA) at 20°C. 24 cm is not a realistic thickness for a bituminous layer, but it was used to highlight the influence of the thickness of the structure. The tire load is simulated by a circular surface, with a diameter of 25 cm and a contact pressure of 0.662 MPa (see Figure 1).

Results are presented on Figure 1 and Figure 2. Figure 1 shows the values and signs of the vertical stress σ V and horizontal (or radial) stress σ H under the centre of the wheel:

-On the top of the pavement (A1), σ V and σ H are both in compression.

-At mid-depth in the pavement (A2), there is compression in the vertical direction and a small tension in the horizontal direction. -At the bottom of the pavement (A3), there is compression in the vertical direction and tension in the horizontal direction (as shown by [START_REF] Heck | Modélisation des déformations réversibles et permanentes des enrobés bitumineux : application à l'orniérage des chaussées[END_REF].

Figure 1. Stress states in a pavement layer under wheel loading

Figure 2 shows the stress paths produced, at different depths, by the moving load, in terms of mean stress p and deviatoric stress q (equations 1 and 2). The stress paths induced in the structure have a complex shape, and range from compression at the top of the layer to tension, at the bottom, with a large range of q/p slopes (for this example from -2 to 3).

Figure 2. Stress paths at different depths Z in bituminous layers and experimental stress path in the triaxial tests (solid line)

To simulate this range of stress paths during a cyclic triaxial test, it would be necessary to cycle both the axial load and the cell pressure. However, this requires a complex loading system, and for the IFSTTAR apparatus, it was decided to work, first, with a constant confining pressure, and to add the possibility of a variable confining pressure in a second phase. The experimental stress paths were thus limited to a constant slope q/p = 3, and to stress states in compression (p>0), as presented on figure 2.

Test equipment and material tested

Triaxial apparatus with temperature control

The triaxial apparatus has been developed using an existing servo-hydraulic loading frame. The pressure cell, the temperature control system, and the instrumentation have been developed specifically for triaxial tests on bituminous mixtures (see figure 3).

The apparatus is designed for specimens with a diameter of 80 mm and a height of 160 mm. The maximum axial load is 25 kN, the maximum confining pressure is 700 kPa, and the temperature range is from 5°C to 60°C. Monotonic and cyclic axial loading can be applied, keeping a constant confining pressure. The maximum loading frequency is 10 Hz.

Cyclic triaxial tests on bituminous mixtures 7

Triaxial cell and loading frame

The loading frame used to apply the axial load is a SCHENCK frame, with a servohydraulic actuator, with a load capacity of ±100 kN and a ± 100 mm axial stroke. The loading frame is equipped with a control system, which allows to perform a large range of loading tests (monotonic or cyclic) with load or displacement control. The triaxial cell for bituminous materials has been developed by GDS Instruments (see figure 3). It is made of aluminium, has an external diameter of 300 mm, and can accept 80 mm or 100 mm diameter specimens. The confining fluid that is also used for temperature regulation can be water or air. Water ensures a better homogeneity of the temperature. The cell is equipped with an internal load transducer of 25 kN capacity, and a 1000 kPa pressure transducer, placed outside the cell, on the confining fluid supply circuit. The temperature in the triaxial cell is regulated using a double heating and cooling system (figure 3). The first system consists of a coil, situated inside the triaxial cell, in which circulates the heating and cooling fluid (water with anti-freeze). This coil is connected with an external, removable heating and cooling system. To improve the control of the temperature in the cell, a thermal housing is added around the triaxial cell. A thermocouple placed inside the cell, at mid-height of the specimen is used to control the temperature. To ensure stabilization of temperature before testing, specimens are stored in the triaxial cell at the testing temperature for at least 4 hours prior to testing.

BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLOR AND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

Instrumentation

Two kinds of sensors, placed on the central part of the specimen, are used in the tests:

-Strain gages are used to measure small strains, in a range of about 10 to 500 microstrains -LVDT sensors are used to measure larger strains, corresponding to the permanent deformations of the material, which can reach up to 5 per cent of deformation.

Strain gages

The axial and radial small strains of the specimen are measured using four strain gages: two axial and two radial (see figure 4). The axial strains are measured on the central part of the specimen. The two radial gages are placed at mid-height of the specimen, along the circumference. The two gages are placed vertically, on opposite sides of the specimen. The average of the signals of the two gages is used to calculate the axial strain. The homogeneity of the strain field is also checked by comparing the values given by each sensor. The gages have a measurement range of 1000 μdef and a resolution of 5 μdef. The procedure to glue the gages on the specimen is described in [START_REF] Sohm | Prédiction des déformations permanentes des matériaux bitumineux[END_REF] and [START_REF] Gabet | Comportement triaxial du béton sous fortes contraintes : influence du trajet de chargement[END_REF].

Axial and radial displacement transducers

The axial and radial strains of the specimen are also measured by means of LVDT sensors. The axial strains are measured on the central part of the specimen using two LVDTs, placed vertically, on opposite sides of the specimen (see Figure 5). The

Figure 4. Axial and radial strain gages glued on a specimen (at the end of the test)

Cyclic triaxial tests on bituminous mixtures 9 average of the signals of the two sensors is used to calculate the axial strain. The homogeneity of the strain field is also checked by comparing the values given by each sensor. The radial strains are measured using an articulated ring, equipped with an LVDT, placed at mid-height of the specimen. The radial LVDT measures the opening of the ring, and thus the variations of the specimen diameter. The axial LVDTs and the articulated ring are glued on the specimen, using metallic clamps. The LVDTs have a measurement range of 5 mm and a resolution of 5 μm.

Figure 5. View of specimen with the axial and radial LVDT sensors

Experimental procedures

Material and sample preparation

The asphalt concrete tested in this study (using the triaxial apparatus) was a French O/10 mm dense-graded bituminous mix, named BBSG or "Béton Bitumineux Semi-Grenu". This material is typically used for wearing courses in France. The characteristics of this mix are given in Table 1. It has been designed with a void percentage of 4.5%. The void percentage was measured by a nuclear (gamma ray) method (EN 12 697-7). The mean of the void percentage, on 60 specimens, was 4.54% of voids and the standard deviation was 0.62% of voids (13.6% in relative value). Once mixed, the material is compacted using a slab compactor with rubber-tired wheels (EN 12697-33 + A1). The dimensions of the plates are 600 mm by 400 mm by 120 mm thick. After a minimum of two weeks of storage, 6 cylindrical specimens (80 mm in diameter and 160 mm high) were cored horizontally, from the central part of each slab. It would have been more relevant to core them vertically, in the compaction direction, but the thickness of the plates did not allow this. The specimens were stored in a room at a controlled temperature of 18°C before testing.

Preparation of triaxial test specimen

To perform tests in compression and tension, a metallic platen is glued to each end of the specimen. These platens are attached to the cell base and to the piston. To ensure homogeneity of the stress states in the specimen, a good alignment of the specimen and end platens is required. A specific stand, for gluing the end platens, has therefore been developed.

Then the LVDT clamps and the four strain gages are glued on the specimen. A latex membrane is placed around the specimen, above the clamps and finally, the LVDT sensors are fixed on the clamps (see figure 5).

Standard mechanical characterisation of the asphalt concrete

Two standard tests were made to characterize the bituminous mixture: the French wheel tracking test (WTT, large size device, T=60°C, pressure of tyre is 6 bars and wheel frequency is 1Hz) (EN 12697-22) and the complex modulus test. The WTT test gives a rut depth of 7.2% at 30000 cycles. Rutting resistance of the chosen bituminous mixture is considered as medium. The complex modulus test is a twopoint bending test on a parallelepiped. The value of the norm of the complex modulus at 15°C and 10 Hz is 14400 MPa.

Testing program and test data analysis

Description of the tests

The aim of the experimental investigation is to study the influence of confining pressure and frequency on permanent deformations and viscoelastic properties of our material under cyclic loading. During a cyclic test, a cyclic sinusoidal deviatoric Cyclic triaxial tests on bituminous mixtures 11 stress q and a constant confining pressure (or radial stress) σ r are applied to the specimen at the start of the test, and the evolutions of the axial and radial strains, ε x and ε r , are measured (see Figure 6).

The stress state is described using the mean stress p and the deviatoric stress q. For axisymmetric triaxial test conditions, stresses are defined by equations 3 and 5. -Then cyclic sinusoidal test without rest periods, but with several successive loading phases at different frequencies and different confining pressures. These tests are named T2.

Description of test procedure T1

In this test procedure, cyclic compression and rest periods are successively applied leading to a high level of accumulated strain, outside the linear viscoelastic domain.

The rest periods are used to measure the relaxation behaviour of the material after unloading. Applied stresses are presented in Figure 7 and table 2. With the test procedure T1 presented above, studying the influence of the load frequency or the confinement presents some difficulties in terms of repeatability, because each experimental condition is applied on a different specimen. As two specimens are never exactly identical, mainly due to manufacturing (in particular void content variations), it is sometimes difficult to know if differences between results are due to changing parameters between tests or to specimen variability. The aim of the test procedure T2 was to overcome repeatability problems due to manufacturing, by applying variable experimental conditions (different confining pressures or loading frequencies) in a single test. Three tests were performed according to procedure T2. In the first test (T2-1), the variable parameter was the Cyclic triaxial tests on bituminous mixtures 13 confining pressure. The temperature, the cyclic deviatoric load and the frequency were kept constant all along the test. In the tests T2-2 and T2-3, the variable parameter was the frequency. In each test, several phases were set up. Tests are presented in detail in paragraph 5.2.2.

Experimental program

Tests using procedure T1

A synthesis of all the tests T1 performed in this study is presented in Table 3. For all T1 tests, the temperature is 20°C. The values of q min and q max are : q min =60 kPa and q max =460 kPa. So the mean deviatoric stress q m is 260 kPa and the amplitude of the sinusoidal deviatoric stress is 200 kPa.

Tests using procedure T2

The test T2-1 included 3 phases. Each phase consisted of 5 loading sequences, with increasing confining pressures (σ r = 0, 50, 100, 200 or 400 kPa). The stress states applied during this test are shown in figure 8. The tests T2-2 and T2-3 included 4 phases. Each phase consisted of 5 loading sequences, with increasing load frequencies (f = 0 (creep), 1, 2, 5 and 10 Hz). The confining pressure was different in the two tests. The stress states applied during these tests are presented in figure 9. The characteristics of each test are summarized in table 4.

For the three tests using procedure T2, the values of q min and q max are : q min =60 kPa and q max =460 kPa. So the mean deviatoric stress q m is 260 kPa and the amplitude of the sinusoidal deviatoric stress is 200 kPa. For the two tests with creep phases, q m is also 260 kPa. In test T2-2, the main parameter studied was the loading frequency, but in phase 4, the confinement was increased to 100 kPa, to have also an idea of the influence of variable frequencies under confinement. As this change was made at the end of the test, it was then decided to perform another test T2-3, with confinement applied at the beginning of the test.

The three samples have been cored in the same plate. The void ratio are respectively 4.4% , 4.7% and 5.0% for the 3 samples.

The analysis of the results of these tests is made in terms of cyclic strain amplitude (determination of complex modulus and Poisson ratio) and in terms of mean strain rate (evolution of viscoplastic strains). In paragraph 4.3, the sinusoidal regression used for analyzing data in terms of cyclic strain amplitude is presented. In paragraph 4.4, the method for studying strain rates discontinuities is presented.

Sinusoidal regression for analyzing test data

Each cyclic test represents a large quantity of data. About 10 signals are saved at the same time (load, vertical displacement, confining pressure, local axial and radial strains). A test at 10 Hz, lasting 10000 seconds (around 3 hours) is considered. If 100 points per cycle are recorded, this will lead to a file with 10 million lines. It was decided to limit the amount of data by modelling solicitations and responses. As both solicitation and response are supposed to be sinusoidal, a sinusoidal regression was used to fit the experimental data. This allowed to decrease the amount of stored data and to calculate for each cycle the mean value and amplitude of each signal. From these values, complex modulus and permanent deformations could be determined.

The evolution with time of the deviatoric stress was fitted with the following function:

( ) ( ) ( ) ( ) ( ) ( ) ( ) { } 0 0 , sin Im i m m
q t q N q N q N q N q N e ωt t ωt = = + = + [START_REF] Airey | Viscoelastic linear limits for bituminous materials[END_REF] Where :

-T is the period of the cyclic loading and 2 T ω π = the angular frequency -N is the cycle number -t is the local time within the cycle number N, ( )

1 t N T t =- - , 0 T t ≤ ≤ -( ) 0
q N is the amplitude of the sinusoidal deviatoric stress -( ) m q N is the mean deviatoric stress

The measured axial and radial strains oscillate around a mean non-linear deformation curve (see figure 10). Evolutions with time (for two consecutive cycles) of these strains were fitted using the following functions ( [START_REF] Dongmo | Caractérisation des déformations d'orniérage des chaussées[END_REF], [START_REF] Neifar | Comportement thermomécanique des enrobés bitumineux : expérimentation et modélisation[END_REF] and [START_REF] Kappl | Enhanced algorithms for the derivation of material parameters from triaxial cyclic compression tests on asphalt specimen[END_REF]):

Cyclic triaxial tests on bituminous mixtures 17 The equations for the complex modulus and the complex Poisson ratio are:

( ) { } ( ) { } 0 * 0 0 * 0 ( ) ( , ) ( ) ( ) ( )sin ( ) ( ) ( ) Im ( ) ( ) ( , ) ( ) ( ) ( )sin ( ) ( ) ( ) Im ( ) x x xm x x x i xm x x r r rm r r r i rm r r t N N N N N N N N e t N N N N N N N N e ωt ωt ε ε t ε α t ε ωt ϕ ε α t ε ε ε t ε α t ε ωt ϕ ε α t ε  = = + + +  = + +   = = + + +   = + +  
( ) ( ) ( ) 0 * * * 0 i x q t E E e t ϕ ω ε = = [11]
( ) ( ) ( ) ( ) * 0 * * 0 * 0 0 r i r x i r x x N e e N ν ϕ ϕ π ϕ ε ε ν ω ν ε ε -- = = = [12]
Where :

-* E and * ν are the norms of the complex modulus and Poisson ratio.

-ϕ and ν ϕ are the phase angles of the complex modulus and Poisson ratio.

-0 q is the amplitude of the sinusoidal deviatoric stress -0 x ε and 0 r ε are the amplitudes of the sinusoidal axial and radial strains x ϕ and r ϕ are the phase lags between the axial and radial strains and the axial stress.

It is important to emphasize that for the stress levels applied, the amplitude of the sinusoidal axial strain is more than 10 -4 m/m. According to [START_REF] Airey | Viscoelastic linear limits for bituminous materials[END_REF] and [START_REF] Airey | Linear rheological behavior of bituminous paving materials[END_REF], beyond this level, the behaviour of the bituminous mixture cannot be considered linear viscoelastic anymore. So, the modulus calculated in the following is not a true complex modulus.

Regression for strain rate discontinuity analyses

The principle of the strain rate analysis is presented in Figure 11. In this figure, a part of the axial response ε m of the test T2-2 is presented. The mean axial strain is a function of time. At time t=1500s, the frequency is changed from f=0 Hz (creep) to f=1 Hz. The experimental points of the axial strain clearly show a discontinuity in the slope at this time. The analysis consists in assessing the mean strain rate discontinuity at this time. The relative strain rate discontinuity ti xm ε ∆ at time ti is defined as follows:

( )

100 2 1 (%) ⋅ + - = ∆ - + - + ti xm ti xm ti xm ti xm ti xm ε ε ε ε ε      [13] + ti xm

Results

Repeatability of triaxial cyclic sinusoidal compression tests (T1)

A triaxial cyclic sinusoidal compression test, made at T=20°C, q m =260 kPa, σ r =100 kPa and f=1 Hz, was repeated 3 times.

Test 2 and test 3 were made on two samples coming from the same manufactured plate. Test 5 was made on a sample cored from another plate. Their void contents were respectively 4.5%, 4.6% and 5.7%. Figure 12b presents the volumetric strain as a function of axial strain. For the three tests, volumetric strain is slightly contractant at low strain level, and then strongly dilatant at higher strain levels. These results tend to show that a denser bituminous mixture (test 2 and 3) is less deformable (referring to axial strains) and more dilatant than a less dense one (test 5).

Presently, it is difficult to assess if the differences observed are due only to different densities or more generally to the different origins of the specimens. More repeatability tests should clearly be performed, but the results show that using specimens cored in the same plate reduces the scatter of the results. During these three tests, complex modulus and Poisson ratio were also measured. For the three tests, mean values were very close, with a standard deviation of 5 %, except for the phase angle of the Poisson ratio, which value is very low (about -3°) (see table 5). Repeatability problems seem to be less important for complex modulus than for permanent deformations. It is also interesting to note that the density variations between the specimens did not affect the viscoelastic properties.

Table 5. Standard deviation of complex modulus and Poisson ratio and their phase angle

Loading and unloading behaviour (T1)

In order to separate the recoverable and irrecoverable strains, rest periods were applied at the end of each loading phase. During the rest periods, the deviatoric stress q was not set to 0 (due to constraints of the servo-hydraulic loading system) but to 60 kPa. The strains were measured during these "recovery" phases.

An example of strain response during loading and unloading, for the test at q m =260 kPa, σ r =100 kPa, T=20°C and f=1 Hz (test 2) is shown in figure 13.

Figure 83. a) Loading and unloading phases, for test 2; b) Part of irrecoverable and recoverable strains

It can be seen that the resilient strain after unloading is very small compared to the irrecoverable strain. The results were similar in all tests, the recoverable strain never exceeds 6% of the total strain. [START_REF] Dongmo | Caractérisation des déformations d'orniérage des chaussées[END_REF] has already presented similar conclusions, without confining pressure.

Influence of the confining pressure on cyclic response

Influence on permanent deformations : test procedure T1

Three tests with procedure T1 have been performed at 20°C, 1Hz, at a mean deviatoric stress of q m =260 kPa and at three different confining pressures: 0, 100 and 200 kPa. The three specimens were cored in the same plate, in order to limit problems of repeatability. The results of these tests are presented in Figure 14a, b andc. They show that increasing the confining pressure reduces the axial and radial strains. Figure 14d represents the volumetric strain function of axial strain. The test BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLOR AND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293 without confining pressure leads to the highest negative volumetric strains. The two tests with confining pressure lead to smaller negative strains. This effect of the confining pressure clearly confirms the benefits of triaxial tests for understanding and modelling the permanent deformation behaviour. Behaviour observed in triaxial tests is clearly different from uniaxial behaviour: when a confinement is applied, it leads to a significant decrease of permanent strains, and a change of the volumetric behaviour; less dilatancy is observed when the confinement increases. In test T2-1, one sample was submitted to sinusoidal creep tests at different and successive confining pressures, to avoid problems of repeatability. The temperature (20°C), the frequency (1 Hz) and the mean level of the deviatoric load (q m =260 kPa) were maintained constant during the test.

Cyclic triaxial tests on bituminous mixtures 23

Results of test T2-1 are summarized in Table 6. Because the confining pressure was changed several times during this test, the analysis of the results is focussed on the strain rate discontinuity at each change of confining pressure. According to the results of tests T1, it is expected that an increase of confining pressure leads to a decrease of the strain rate. This trend is observed at each change of confinement in the first phase of the test and for most confinement levels in phases 2 and 3, except 3 contradictory results. These contradictory results can be observed for low confinement changes, in phase 2, from 0 to 50 and from 50 to 100 kPa, and in phase 3, from 0 to 50 kPa. These contradictory results are difficult to explain, they may be due to an imprecise data fitting. But they could also be due to the fact that increasing the confinement also leads to increase the axial stress (σ x = σ r + q). And it is thought that increasing both confinement and axial stress leads to opposite trends in term of strain rates. In conclusion, we consider that results of tests T1 are confirmed because increasing significantly the confinement leads to decrease significantly the axial strain rate. Complex modulus is measured in laboratory, using a two point bending test on a parallelepiped. The value at 20°C and 1Hz is 6630 MPa. That is close to the values measured during the cyclic sinusoidal test (T2-1). Values of complex modulus, Poisson ratio and phase angle remain practically the same for the 5 levels of confining pressures (from 0 to 400 kPa). There is no significant influence of the confining pressure on the viscoelastic properties of the tested bituminous mixture. These results confirm those obtained by [START_REF] Clec | Multidirectional behaviour of bituminous mixtures[END_REF].

Influence of loading frequency

6.4.1. Influence on permanent deformations: test procedure T1.

Two tests have been performed at 20°C, at a mean deviatoric stress of q m =260 kPa, at a confining pressure of σ r =100 kPa, and at two different frequencies: 1 Hz and 10 Hz. The two specimens were taken from the same plate.

Figures 15a and15b show the axial and radial strains as functions of time. For the two frequencies, the axial strains are the same. Radial strains are slightly more important at 1 Hz than at 10 Hz. These results were not really expected.

Figures 15c and15d show the axial and radial strains as function of the number of load cycles. In this scale, axial and radial strains are much higher at 1Hz than at 10 Hz. This seems to indicate that the development of permanent strains depends much more on the loading rate, or on the total loading time than on the number of cycles (on figure 15a

Influence on permanent deformation: test procedure T2

Tests T2-2 and T2-3 were performed at 20°C, at a mean deviatoric stress of 260 kPa and a confining pressure of 0 or 100 kPa, depending on the testing phase (see Table 8). The objective was to analyse the effect of the frequency on the strain rate discontinuity. Results are presented in Table 8. Whatever the test and loading phase, changing from a creep (static) solicitation to a cyclic solicitation leads to an BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLOR AND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293 important increase of the strain rate. Then, under cyclic loading, increasing the frequency leads to a decrease of the strain rate. It is true when the frequency is increased from 1 Hz to 5 Hz and from 5 Hz to 10 Hz. Decreasing the frequency, here from 10 Hz to 1 Hz, always leads to a strong increase of the strain rate. All these strain rate variations are summarized in Table 9. A test at a frequency of 1 Hz leads to an axial strain rate 45% higher than a creep test; a test performed at a frequency of 1 Hz leads to an axial strain rate 30% higher than a test performed at 10 Hz. According to these results, it is clear that changing frequency has an effect on permanent strain rate. This may be contradictory with the result presented in figure 15 a). Several explanations maybe proposed. As a first possible explanation, it is thought that the effect of changing frequency is temporary only, so that permanent deformation strains would be finally the same after a long loading time. However, this assumption is not in accordance with the results of the literature, without confining pressure [START_REF] Dongmo | Caractérisation des déformations d'orniérage des chaussées[END_REF]. The second explanation considers that the two results presented in figure 15 a) are accidentally the same. The problems of repeatability of both tests may have lead to higher permanent strains value than usual at 10 Hz and/or lower permanent strains than usual at 1 Hz, so that the expected influence of changing frequency could not be highlighted. We are convinced that some other tests would validate this hypothesis, which would more suit with the literature. Explaining why cyclic loading leads to higher permanent deformations than creep loading is not obvious. The authors try to explain this by two reasons. Firstly, it is thought that cycling loads facilitates the relative movements of the grains (rotations, displacements) of the granular skeleton. It is secondly thought that cycling around a mean value necessarily lead to reach higher stress levels (in our tests, a mean value of q at 200 kPa leads to reach 400 kPa). More plastic phenomenon may appear due to reached higher values and passed thresholds. It would also be interesting to explain why lower frequencies lead to higher permanent deformations than higher frequencies. The authors think that for a given test duration, lower frequencies allows to remain at higher stress levels during longer periods (the time of a half cycle), which may more affect the material than reaching more often the higher level. It is thought that at lower frequency, higher viscoelastic strains are reached (see viscoelastic results), then relative movements of the grains of the granular skeleton are more important; this may facilitate the setting of irreversible (permanent strains). -41% -52% -81% -58% -43% -41% -70% -51% 1Hz -> 5Hz -15% -28% -32% -64% -35% -13% -25% -31% -64% -33% 5Hz -> 10Hz -9% -13% -11% -10% -11% -6% -12% -15% -10% Table 10 shows that:

Cyclic triaxial tests on bituminous mixtures 27

• For each frequency (1, 5 and 10 Hz), complex modulus, Poisson ratio, and phase angle are the same with or without confining pressure. • With or without confining pressure, the complex modulus increases and the phase angle of the complex modulus decreases when frequency increases.

Cyclic triaxial tests on bituminous mixtures 29

• For the Poisson ratio, the conclusions are less clear, due to some experimental scatter, but it seems that increasing the frequency leads to a decrease of the Poisson ratio and an increase of the phase angle of this coefficient. However, this phase angle remains small (between 2° and 5°). Di [START_REF] Di Benedetto | Time-temperature superposition principle for bituminous mixtures: three dimensional approach and extension in the non linear domain[END_REF] carried out similar studies, without confining pressure, and have come to the same conclusions. It seems that variations of viscoelastic properties with frequency are not affected by the confining pressure, for the range of pressures applied in this study (0 kPa to 200 kPa), which is representative of pavement loading conditions.

Conclusions and future developments

Conclusions

In this paper cyclic triaxial thermo-regulated tests on bituminous mixtures are presented. These tests were performed at different loading frequencies and confining pressures, to evaluate the influence of these parameters on the viscoelastic properties and permanent deformations of the material.

Two types of test procedures were developed and applied to our specimens. The test procedure T1 was a classical cyclic triaxial test, which consisted in applying a single loading level and frequency all along each test. However, some problems of repeatability where emphasized. They were attributed to the manufacturing process and to scatter of specimen densities and structures. The second procedure named T2 was developed to overcome problems of repeatability. This procedure consisted in applying successively different loading conditions (confining pressure or loading frequency) on a single specimen, to avoid sample preparation effects. On these tests, a focus was made on the strain rate variations.

Our experimental program highlighted several interesting results. Firstly, it showed that slight variations of specimen density and structure due to the manufacturing process can have a significant effect on permanent strain behaviour. Once overtaken these problems of repeatability, tests have clearly shown that :

• Increasing the confinement reduces significantly permanent deformations.

• Increasing the loading frequency reduces permanent deformation rate

• Creep tests lead to less deformation than sinusoidal cyclic tests at low frequency, at the same mean value. This emphasises the influence of cyclic solicitations on the response of bituminous mixtures.

• In terms of viscoelastic properties, results have confirmed the influence of the frequency of loading on the viscoelastic response of the material, classically observed in complex modulus tests. The tests have also shown no significant effect of the confining pressure on the axial response of the material.

It is clear that both the triaxial aspect of loading, and cyclic effects cannot be neglected when studying bituminous mixtures. From a rheological point of view, BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLOR AND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293 this highlights that their behaviour results from the contributions of the viscous binder (viscoelastic properties changing with temperature and frequency) and of the granular skeleton (influenced by the triaxial loading path). From these observations, it can be concluded that a good rheological model for a bituminous mixture needs to take into account all these characteristics : viscoelasticity, viscoplasticity, thermosuceptibility, cycling effects..

Outlook

The final objective of the presented study is to develop efficient rheological models for FEM simulations for the prediction of the permanent deformations of pavement structures, under monotonic and cyclic loading. Expected applications are for instance :

• Rutting of road pavements under severe loading conditions (channelized traffic, slopes, etc..)

• Evaluation of the effects of high temperatures on pavement rutting (impacts of climate change)

• Rutting of industrial platforms under static or moving loads

• Effect of very high loads or tyre pressures on pavements (aircraft loads…).

Cyclic triaxial tests on bituminous mixtures 31

Figure 3 .

 3 Figure 3. Scheme of the triaxial cell

Figure 6 .

 6 Figure 6. Principle of a cyclic triaxial test performed in this study. Definition of the variables q and σ r . The soil mechanics sign convention is used. Two types of test procedures were used during the experimental program: -First cyclic sinusoidal tests with successive load periods and rest periods, with fixed experimental conditions (frequency and confining pressure) all along the tests. These tests are named T1.

Figure 7 .

 7 Figure 7. Applied stress sequences for tests T1

Figure 8 .Figure 9 .

 89 Figure 8. Applied stress paths for test T2-1

  lags between the axial and radial strains and the axial stress at cycle N, which are not considered as equal.

Figure 10 .

 10 Figure 10. Analysis of experimental data. Sinusoidal regression.

  from polynomial regressions fitted on the values of the axial mean strain function of time, before and after the discontinuity respectively. The axial mean strain equation[9]. It has to be noticed that it is necessary to assess this discontinuity without regression because of some scatter in the experimental data.The authors would point out that in tests T2, considering the viscoplastic strain, it has no sense to consider the total accumulated strains (from the beginning to the end of the test) because loading modalities change several times during the test. At Cyclic triaxial tests on bituminous mixtures 19 each modality change (confinement or frequency), a focus is made on the strain rate discontinuity, defined in equation[10].

Figure 6 .

 6 Figure 6. Strain rate analysis. Definition of the strain rate discontinuity at the time t=1500 sec, where frequency is changed from f=0 Hz (creep test) to f=1 Hz.

  Figure 7. a) Axial strain function of time and b) volumetric strain function of axial strain, for the three repeatability tests

  a) Axial strains as a function of time; b) Radial strains as a function of time; c) Axial strains as a function of number of loads and d) volumetric strain as a function of axial strain, for tests at 20°C, 1Hz, q m =260 kPa and different confinements 6.3.2. Influence on strain rate : test procedure T2

  Figures15c and 15dshow the axial and radial strains as function of the number of load cycles. In this scale, axial and radial strains are much higher at 1Hz than at 10 Hz. This seems to indicate that the development of permanent strains depends much more on the loading rate, or on the total loading time than on the number of cycles (on figure15a and 15b, multiplying the number of cycles by 10 does not change the results, provided that the loading time is the same).

  BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLORAND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

Table 1 .

 1 Characteristics of the tested bituminous mixtureBLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLORAND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

	Size of aggregates (mm)	Percentage (by mass)
	0/2	27%

Table 2 .

 2 Load history for tests T1

				frequency
				1 Hz	10 Hz
		sequence	Cycle	Time (s)	Cycle	Time (s)
	1	load rest	30 0	30 120	300 0	30 120
	2	load rest	100 0	100 400	1000 0	100 400
	3	load rest	300 0	300 1200	3000 0	300 1200
	4	load rest	1000 0	1000 4000	10000 0	1000 4000
	5	load rest	3000 0	3000 12000	30000 0	3000 12000
	6	load rest	10000 0	10000 40000	100000 0	10000 40000
	5.1.2. Description of test procedure T2

Table 3 .

 3 Synthesis of tests T1

	Test	σr	Frequency Material plate	Void content
		kPa	Hz		%
	1	0	1	11	5.3
	2	100	1	12	4.5
	3	100	1	12	4.6
	4	200	1	11	4.9
	5	100	1	11	5.7
	6	100	10	12	4.2

  BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLORAND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

  BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLORAND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

Table 6 .

 6 Test conditions and results for test T2-1Mean values of complex modulus, Poisson ratio and phase angle were calculated for each confining pressure and each loading phase. Results are presented in table 7.BLANC, Juliette, GABET, Thomas, HORNYCH, Pierre, PIAU, Jean Michel, DI BENEDETTO, Hervé, 2015, Cyclic triaxial tests on bituminous mixtures , Road materials and pavement design, 1, 16, TAYLORAND FRANCIS, pp.46-69, DOI: 10.1080/14680629.2014.964293

			test T2-1
		confinement changing	test time (sec)	strain rate discontinuity at test time
		start	0	//
		0 -> 50 MPa	225	-6%
	phase 1	50->100 MPa 100 -> 200 MPa	450 675	-10% -11%
		200 -> 400 MPa	900	-15%
		400 -> 0 MPa 1125	15%
		0 -> 50 MPa 1625	4%
		50->100 MPa 2125	9%
	phase 2	100 -> 200 MPa 2625	-5%
		200 -> 400 MPa 3125	-18%
		400 -> 0 MPa 3625	41%
		0 -> 50 MPa 4625	5%
		50->100 MPa 5625	-20%
	phase 3	100 -> 200 MPa 6625	-18%
		200 -> 400 MPa 7625	-23%
		stop 8625	//
	6.3.3. Influence of confining pressure on viscoelastic properties

Table 7 .

 7 Values of complex modulus, Poisson ratio and phase angles, in function of confining pressure at 20°C, 1 Hz.

	Phase	σ rad	|E*|	φ	|ν*|	φ ν
		0 kPa	5790 MPa	30.2°0.31	-2.6°5
		0 kPa 5840 MPa	29.4°0.32	-2.6°1
		00 kPa 5860 MPa	29.3°0.32	-3.2°2
		00 kPa 5860 MPa	28.8°0.32	-3.1°4
		00 kPa 5870 MPa	28.8°0.32	-3.1°0
		kPa	5850 MPa	29.2°0.32	-2.5°5
		0 kPa 5830 MPa	29.1°0.32	-2.4°1
	2	00 kPa 5830 MPa	29.1°0.31	-2.5°2
		00 kPa 5830 MPa	28.8°0.32	-2.4°4
		00 kPa 5880 MPa	28.6°0.31	-2.2°0
		kPa	5860 MPa	29.2°0.31	-2.2°5
		0 kPa 5830 MPa	29.2°0.31	-2.3°1
	3	00 kPa 5820 MPa	29.2°0.31	-2.3°2
		00 kPa 5860 MPa	29.0°0.31	-2.3°4
		00 kPa 5950 MPa	28.6°0.31	-2.2°1

Table 8 .

 8 Test conditions and results of testsT2-2 and T2-3

			test T2-2	test T2-3
	frequency changing	test time(sec)	strain rate gap at test time
	start	0	//	//
	creep -> 1Hz	200	19%	36%
	1Hz -> 5Hz	400	-15%	-13%
	5Hz -> 10Hz	600	-9%	-6%
	10 -> 1Hz	800	13%	28%
	1Hz -> creep	1000	-41%	-43%
	creep -> 1Hz	1500	36%	53%
	1Hz -> 5Hz	2000	-28%	-25%
	5Hz -> 10Hz	2500	-13%	-12%
	10 -> 1Hz	3000	26%	28%
	1Hz -> creep	3500	-52%	-41%
	creep -> 1Hz	4500	51%	45%
	1Hz -> 5Hz	5500	-32%	-31%
	5Hz -> 10Hz	6500	-11%	-15%
	10 -> 1Hz	7500	41%	41%
	1Hz -> creep	8500	-81%	-70%
	creep -> 1Hz	9500	59%	59%
	1Hz -> 5Hz	10500	-64%	-64%
	5Hz -> 10Hz	11500	-10%	-10%
	stop	12500	//	//

Table 9 .

 9 Summary of the test conditions and results for testT2-2 and T2-3

	frequency changing	mean	mean
	creep -> 1Hz 19% 36% 51% 59% 41% 36% 53% 45% 59% 48%
	1Hz -> creep

Table 10 .

 10 Values of complex modulus, Poisson ratio and phase angle, from triaxial tests, as function of frequency, for test T2-2

	Phase	σ rad	|E*|	φ	|ν*|	φ ν
		creep	-	-	-	-
	1 σ r = 0 kPa	1 Hz Hz 0 Hz	6140 MPa 9840 MPa 11600 MPa	30°0.37 23°0.35 22°0.34	-3°5 -4°1 -5°1
		Hz	7140	29°0.38	-4°c
		reep	-	-	-	-
		1 Hz	6100 MPa	29°0.38	-3.8°5
		Hz	-	-	-	-
		10 Hz	11640 MPa	21°0.34	-5.1°1
		Hz	6070 MPa	29°0.39	-4.3°c
		reep	-	-	-	-
		1 Hz	6100 MPa	29°0.31	-2.2°5
		Hz	9900 MPa	23°0.31	-2.3°1
		0 Hz	11720 MPa	21°0.31	-2.3°1
		Hz	6030 MPa	29°0.31	-2.2°c
		reep	-	-	-	-
		1 Hz	6050 MPa	29°0.4	-4°5
		Hz	9910 MPa	24°0.37	-4.9°1
		0 Hz	-
		1 Hz	-	-	-	-