# Continuum Damage Approach to Asphalt Concrete Fatigue Modelling

Didier Bodin, Gilles Pijaudier-Cabot, Chantal de La Roche, Jean-Michel Piau, Armelle Chabot

► To cite this version:

Didier Bodin, Gilles Pijaudier-Cabot, Chantal de La Roche, Jean-Michel Piau, Armelle Chabot. Continuum Damage Approach to Asphalt Concrete Fatigue Modelling. Journal of Engineering Mechanics - ASCE, 2004, 130 (6), pp.700-708. 10.1061/(ASCE)0733-9399(2004)130:6(700) . hal-01004693

HAL Id: hal-01004693

https://hal.science/hal-01004693

Submitted on 23 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuum Damage Approach to Asphalt Concrete Fatigue Modeling

Didier Bodin<fnref n="1" />, Gilles Pijaudier-Cabot<fnref n="2" />, Chantal de La Roche<fnref n="1" />, Jean-Michel Piau<fnref n="1" />, and Armelle Chabot<fnref n="1" />

**Abstract:** A nonlocal damage model is proposed to predict the behavior of pavement fatigue cracking. This constitutive relation has been implemented in a finite-element code, along with a self-adaptive jump-in-cycle procedure for high cycle fatigue computations. Strain localization analysis shows that during uniaxial fatigue tests, bifurcation due to strain softening occurs much later than in monotonic tests. The incorporation of an internal length into the constitutive model is advocated since the model should encompass loading histories with very different amplitudes of cycles, in which localization may still occur. The influence of the internal length on the fatigue life of bending beams is also investigated. Calibration of the damage model is performed after thermal effects have been evaluated and accounted for in a simplified way, uncoupled to damage. Parameter identification is performed in bending and uniaxial tests. The resulting calibrated constitutive relation is found to yield a good description of several different uniaxial tests.

**keywords:** Asphaltic concrete; Fatigue; Damage; Cracking; Pavements; Constitutive relations.

## Introduction

Repeated stress and strain induced by traffic loads can be very damaging to asphalt pavements. Cyclic mechanical loading leads to fatigue damage in the bottom layers. Accumulation of these small degradations can initiate macrocrack propagation through the pavement body from the bottom to the top. Design methods have to be developed to avoid the initiation and propagation of cracks, at least over the life of service of the structures. For asphalt pavements, design methods often rely on standardized material testing. In the French design method, the principle is to compare the strain level calculated at the bottom of the asphalt layer to admissible strain, which is strain leading to failure for 1 million cycles during a standardized fatigue test, modified by a shift factor that takes into account in situ conditions that are more complex than in laboratory test conditions (healing, temperature, frequency, etc.).

Most laboratory fatigue tests consist of applying constant amplitude sinusoidal displacement or force at the boundary of a bi-tuminous mixture sample (see, e.g., SHRP 1994; Di Benedetto et al. 1997; Partl et al. 2000). During the tests, the variation in global stiffness is monitored. It is defined as the ratio of the force amplitude applied to the displacement amplitude. Because it is difficult to relate laboratory test results to field data, such experiments provide, in most cases, groundwork for comparisons between several different types of materials only (mix, granulometry, etc.). Rest periods can also be integrated into continuous fatigue tests to investigate more realistic fatigue performance (Kim and Little 1990; Breysse et al. 2002) that include healing phenomena which are not dealt with here.

Computational prediction of fatigue life is quite a challenge because permanent strain and variations in stiffness are due to several related physical processes. There is dispersed microcracking in the binder (Lefeuvre 2001), and viscoelastic behavior induces an increase of temperature due to viscous dissipation during cyclic tests associated with important thermosensitivity of the binder (de La Roche and Marsac 1996). During a typical test, the global stiffness, calculated as the ratio of the force to the displacement amplitude, decreases following three regimes of evolution (Fig. 1): During phase I, a fast decrease in stiffness is observed. It is followed by phase II that corresponds to a quasilinear decay in stiffness. Fracture occurs during phase III, due to damage acceleration and to ultimately to the propagation of macrocracks.

There are several approaches by which to model asphalt concrete damage. It can be treated directly with viscoelasticity coupled to damage (Lefeuvre et al. 2000). Kim and Little (1990), Lee and Kim (1998), and Lee et al. (2000) have proposed models in which damage and viscoelasticity are uncoupled, with the help of a corresponding principle given by Shapery (1984) that transforms viscoelastic analysis into elastic analysis. For sinusoidal loading histories, Alimami (1987), Piau (1989), and Ullidtz et al. (1997) provided an interpretation of the variation in material stiffness due to damage as variation of the complex modulus. It is this variation of the complex modulus due to damage that is considered here. Therefore the viscoelasticity of the material will not be considered in this paper.

This work focuses on the description of dispersed microcracks in material during fatigue load histories. The constitutive model is

<footnote n="1">¹Laboratoire Central des Ponts et Chaussées, Route de Bouaye, B.P. 4129, F-44341 Bouguenais Cedex, France</footnote>

<footnote n="2">²R&DO-GeM, Institut de Recherche en Génie Civil et Mécanique, Ecole Centrale de Nantes, B.P. 92101, F-44321 Nantes Cedex 3, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Evolution of stiffness during typical fatigue test. Notice the three regimes: in the first damage evolution decelerates and the second is almost linear. In the third, damage growth accelerates.</figcaption></figure>

based on an elastic isotropic continuum damage model, which expresses the decrease of elasticity parameters with loading. A simple model has been chosen and implemented in a finite-element code in order to highlight performance trends.

The complete model is described next in this paper. Then the problem of high cycle fatigue computations in the context of finite-element analyses is given. The method of Peerlings et al. (2000) is implemented in the finite-element code Castem 2000. In order to capture ultimate failure and strain localization consistently, a nonlocal approach to the damage model is advocated. The influence of the internal length, which enters in the nonlocal approach, is investigated on a cantilever beam subjected to controlled displacement cycles. Full thermomechanical coupling is not considered, but the mechanical response of specimens is corrected in order to account for the variation in structural stiffness due to a rise of temperature in a simplified way. The principle behind the correction is outlined. Model calibration and comparisons with test data are also provided.

## Mechanical Damage Modeling

The mathematical model used to describe mechanical damage is an elasticity based damage model for fatigue. It is inspired by the constitutive relations of Mazars (1984) and Pijaudier-Cabot and Bazant (1987). Compared to existing proposals, e.g., by Lee et al. (2000), the constitutive relations are simpler since they rely on a scalar, isotropic, damage model which is widely used for concrete and other quasibrITTLE materials. It should be mentioned that for concrete at least, isotropy of damage may not be a major restriction in terms of modeling capabilities. Although microcracks are always pretty much oriented with respect to the loading system, the isotropic damage approach provides reasonable predictions when damage is due to tensile strain or stress especially (Fichant et al. 1999). The influence of microcracking is introduced via the damage variable *d* that ranges from 0 to 1. The stress strain relation reads

$$ 
\begin{aligned}
\varepsilon_{ij} &= \frac{1 + \nu}{E(1 - d)} \sigma_{ij} - \frac{\nu}{E(1 - d)} [\sigma_{kk} \delta_{ij}] \\
\text{or } \sigma_{ij} &= C_{ijkl} \varepsilon_{kl} = (1 - d) C_{ijkl}^{0} \varepsilon_{kl}
\end{aligned}
\quad (1) $$

*E* and *ν* = Young's modulus and Poisson ratio of the undamaged material, respectively; *ε*<sub>ij</sub> and *σ*<sub>ij</sub> = strain and stress components; *C*<sub>ijkl</sub> and *C*<sub>ijkl</sub><sup>0</sup> = damaged and initial (elastic) secant stiffness of

the material; and δ<sub>ij</sub> = Kronecker symbol. The evolution of damage depends on the amount of tensile stress that the material experiences during mechanical loading. It is expected that repeated tensile stress is the main source of local decohesion in asphalt concrete. The damage growth criterion is based on a modified Rankine criterion (Patzák and Jirásek 2001) with zero threshold damage growth, in which is the average ε̄ of equivalent strain ε̄ defined as

$$ \bar{\varepsilon} = \sqrt{\sum_{i=1}^{3} \left( \frac{\langle \sigma_i \rangle_{+}}{E(1-d)} \right)^2} \quad (2) $$

where ⟨·⟩<sub>+</sub> = Macauley bracket; and σ<sub>i</sub> = principal stress. The average equivalent strain ε̄ is

$$ \bar{\varepsilon}(x) = \frac{1}{V_r(x)} \int_{\Omega} \psi(x-s) \bar{\varepsilon}(s) ds, \quad \text{with } V_r(x) = \int_{\Omega} \psi(x-s) ds \quad (3) $$

where Ω = volume of the structure; V<sub>r</sub>(x) = representative volume at point *x*; and ψ(x-s) = weight function:

$$ \psi(x-s) = \exp\left(-\frac{4\|x-s\|^2}{l_c^2}\right) \quad (4) $$

and l<sub>c</sub> = internal length of the nonlocal continuum. The rate of damage growth is defined as a function of the equivalent strain rate (Paas et al. 1993):

$$ \dot{d} = f(d) \bar{\varepsilon}^{\beta} \langle \dot{\bar{\varepsilon}} \rangle_{+} \quad (5) $$

where *f*(*d*) = function of damage, which will be defined later on, and β = model parameter.

When the state of strain is homogeneous over the specimen, e.g., in uniaxial tension or compression, the nonlocal equivalent strain, Eq. (3), is, by definition, equal to the local equivalent strain, Eq. (2). It is possible to integrate Eq. (5) over the history of loading, which is restricted in this case to displacement controlled cycles of constant amplitude:

$$ \frac{dd}{f(d)} = \bar{\varepsilon}^{\beta} \langle d\bar{\varepsilon} \rangle_{+} \Rightarrow F(d) = \int_{0}^{d} \frac{dd}{f(d)} = \sum_{cycle \ 1}^{N} \frac{\bar{\varepsilon}_{a}^{\beta+1}}{\beta+1} \quad (6) $$

where ε̄<sub>a</sub> = amplitude of the equivalent strain over one cycle; and *F*(*d*) = scalar function of damage. Let us define fatigue life as the number of load cycles needed prior to reaching a critical value of damage *d*<sub>crit</sub>. Eq. (6) provides the corresponding allowable number of cycles *N*<sub>crit</sub> of constant amplitude ε̄<sub>a</sub>:

$$ N_{\text{crit}} = \frac{F(d_{\text{crit}})(\beta + 1)}{\bar{\varepsilon}_{a}^{\beta + 1}} \quad \text{with } F(d_{\text{crit}}) = \text{constant} \quad (7) $$

The value of *d*<sub>crit</sub> is understood as a parameter that is chosen by the user. It is related to the allowable relative decrease in stiffness of the material when subjected to uniaxial fatigue. In the French design method, this maximum decrease in stiffness could serve as an end of life criterion, with the help of corrective terms that take into account in situ conditions. From this expression, it can be also noticed that −(β + 1) corresponds to the slope of the Wöhler curve in a log–log diagram.

Fatigue tests at constant displacement amplitude have been carried out on a single bituminous mix (0/6 mm dense asphalt concrete containing 6.4% of 50/70 pure bitumen), part of a current interlaboratory test program within the framework of the RILEM TC 180 PEB (Partl et al. 2000). Tension compression tests have been performed on cylindrical specimens (height = 120 mm, diameter = 80 mm) with sinusoidal displacement
<figure><img src="image_3.png" /><figcaption>Fig. 2. Wöhler curve for typical bituminous mix loaded in uniaxial tension compression</figcaption></figure>

cycles corresponding to the following constant strain amplitudes: (80, 90, 100, 140, 150, 180 µm/m). Bending tests on cantilever trapezoidal samples (height=250 mm, large clamped base = 56 mm, small loaded base=25 mm and thickness=25 mm) have been also performed. These bending tests were displacement controlled. The displacement cycles were such that the maximum (absolute value) strain in the specimen, evaluated from beam theory, takes the following values: 140, 180, and 220 µm/m. Both types of tests have been performed at temperature of 10°C and frequency of 10 Hz. For each test, the decrease in structural stiffness of the specimen was recorded as a function of the number of cycles. Fig. 2 shows a typical Wöhler curve obtained for asphalt subjected to tension compression cycles of constant strain amplitude. In Fig. 2, the fatigue life is defined for a critical value, $d_{crit}=0.3$. The slope of the Wöhler curve is 5 approximately. The fact that a straight line is recovered in a log-log diagram is a justification of the mathematical form of the damage evolution law in Eq. (5). The function $F(d)$ needs, however, to be determined so that the evolution of stiffness with the number of cycles is recovered. This evolution of damage should be able to describe the three regimes of damage growth shown in Fig. 1. It is also quite important to devise a form of this function in which the duration of the three regimes of damage growth can be adjusted. The function used by Paas (1990) and Bodin et al. (2002),

$$f(d) = Cd^{\alpha} \qquad (8)$$

does not permit such control. Fig. 3 shows a fit of the same test data as in Fig. 2 with this equation. The first two regimes are captured very well. In bending tests, the occurrence of the third regime, with significant acceleration of damage growth, cannot be described properly (see Bodin et al. 2002). In fact, we observed that it was not possible to get a good description of the three phases of damage growth at the same time with the evolution of damage given by Eq. (8). In order to correct this weakness, a new function $F(d)$ is defined as

$$F(d) = \alpha_1 \left\{ 1 - \exp \left[ - \frac{d}{\alpha_2} \right]^{\alpha_3} \right\} \qquad (9)$$

where $\alpha_1$, $\alpha_2$, and $\alpha_3$= three model parameters. Fig. 4 shows the variation in material stiffness in a uniaxial fatigue test (sinusoidal displacement cycles with constant amplitude) for different values of these model parameters. These curves will be quite helpful later on for parameter calibration. The major influence of $\alpha_1$ is on

<figure><img src="image_4.png" /><figcaption>Fig. 3. Top: Wöhler curve for uniaxial tension-compression test and model predictions. Bottom: Calibration of the damage model for 80 × 10⁻⁶ constant strain amplitude. The model response should be within the gray area, between the upper and lower bounds (0–6 mm asphalt concrete containing 6.4% of 50/70 bitumen, 10 Hz and 10°C).</figcaption></figure>

fatigue life (here the number of cycles up to complete failure with zero stiffness). Parameter $\alpha_2$ essentially controls the duration of the three regimes, and keeps the fatigue life constant. The larger it is, the larger the growth of damage over a given number of cycles. The third parameter, $\alpha_3$, acts on the shape and steepness of the evolution of damage, in the second regime, where the rate of damage growth per cycle is almost linear. The larger $\alpha_3$, the lower the slope and the greater the damage growth per cycle.

The present constitutive relation has three parameters involved in evolution law for damage, in addition to the elastic constants and internal length that are in the nonlocal equivalent strain. Before describing calibration of the model, let us first discuss finite-element implementation of this model in the context of high cycle fatigue computations.

## Finite-Element Implementation

The above constitutive equations have been implemented into the finite-element code Castem 2000. The implementation does not differ from that of a standard, nonlocal isotropic damage model (see, e.g., Pijaudier-Cabot 1995). A secant stiffness algorithm is implemented, which avoids computation of the tangent material stiffness which has a large bandwith and is nonsymmetric (Jirásek and Patzák 2002). For high cycle fatigue computations, the fa-
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Sensitivity of parameters in damage model to uniaxial fatigue (α₁ top, α₂ middle, α₃ bottom)</figcaption></figure>

tigue life can be greater than 1 million cycles. Therefore, it does not seem to be realistic to perform a finite-element computation over the entire number of cycles. A minimum of 10 increments per cycle is needed to integrate the evolution of damage correctly over a single cycle (in the case of uniaxial tension compression) and the total number of increments for a typical complete fatigue computation would amount at least to 10 million, which is prohibitive for structural analysis. The so-called “jump-in-cycle” procedure devised by Lemaitre (1996) and improved by Peerlings et al. (2000) has been implemented.

## Jump-in-Cycle Procedure

The jump-in-cycle method is based on projection of the increment of damage computed within one cycle, onto a finite number of cycles. The number of cycles, extended to a continuous variable, is taken as a second time scale, in addition to the one that is used to perform computation within a single cycle incrementally. From an initial state defined at a number of cycles N-1 and for a state of damage d defined at each material point at this time, computation of the increment of damage over one cycle is performed. This increment of damage is then transformed into a rate of damage per cycle:

$$ \frac{\partial d}{\partial N} \approx \int_{\text{cycle N}} \dot{d} dt = G(d(N)) \quad (10) $$

Integration of damage growth between cycle N and cycle N + ΔN is approximated with the help of a trapezoidal rule:

$$ d(N + \Delta N) = d(N) + \frac{1}{2}[G(d(N)) + G(d(N + \Delta N))]\Delta N \quad (11) $$

In the present form, integration is implicit since it requires the computation of G[d(N + ΔN)]. This term is evaluated following Euler forward prediction,

$$ G(N + \Delta N) \approx G(N) + \frac{\partial G(N)}{\partial N} \Delta N \quad (12) $$

When this expression is substituted in Eq. (11), which is linearized to first order in ΔN, an explicit integration scheme is recovered along with evaluation of error e which is propagated from one increment of cycles ΔN to the next one (see Peerlings et al. 2000). The increment of cycles ΔN is adapted so as to limit this error. The increment of cycles is chosen so that the error is kept constant:

$$ \Delta N = \eta \frac{\partial d}{\partial G} \quad (13) $$

where η = fixed parameter in the computation (η = 5 × 10⁻²). The increment of cycles is chosen so that Eq. (13) is satisfied at the Gauss point of the finite-element mesh where growth of damage over a cycle is largest. Near failure, ΔN decreases greatly as the damage growth tends to become localized in a small region of the structure and eventually form a crack.

## Strain and Damage Localization

The constitutive relations described earlier incorporate an internal length. This additional parameter is classically introduced in strain softening models in order to control the occurrence of strain and damage localization and to avoid strain that localizes in a region of zero volume with a vanishing energy dissipation at failure (Pijaudier-Cabot and Bazant 1987; de Borst et al. 1993). In addition, nonlocal constitutive relations prevent spurious mesh dependence and exhibit a structural size effect, which in the case of monotonic loads provides a good indicator for obtaining the internal length from experiments (see Le Bellégo et al. 2003). Here we perform a simple, one-dimensional analysis of localization in order to show the influence of the internal length in fatigue problems.

Consider a one-dimensional specimen, in a state of homogeneous damage and strain, subjected to load cycles of amplitude εₐ. This initial state is defined by the value of damage d₀ and strain ε₀ at an arbitrary point of a specific load cycle. The stress strain relation, in the rate format, about this initial state follows from Eq. (1):
$$
\begin{align}
\dot{\sigma} &= (1-d_0)E\dot{\varepsilon} - \dot{d}E\varepsilon_0 = (1-d_0)E\dot{\varepsilon} \nonumber \\
&\quad - E\varepsilon_0 f(d_0) \bar{\varepsilon}_0^\beta \frac{1}{V_r(x)} \int_\Omega \psi(x-s) \dot{\bar{\varepsilon}}(s) ds \tag{14}
\end{align}
$$

We consider an infinite specimen so that the normalized weight function is independent of the coordinate $x$. Furthermore, damage is assumed to grow during the small perturbation (the strain and strain rates are positive). Equilibrium requires that the incremental stress is homogeneous so one arrives at the following integro-differential equation:

$$
(1-d_0)E \frac{d\dot{\varepsilon}(x)}{dx} - E\varepsilon_0 f(d_0) \bar{\varepsilon}_0^\beta \frac{1}{V_r(x)} \int_\Omega \psi(x-s) \frac{d\dot{\varepsilon}(s)}{dx} ds = 0 \quad (15)
$$

where the local equivalent strain rate has been replaced by the
local strain rate inside the integral because they are equal in the
one-dimensional case. Strain localization occurs when bifurcation
is observed in the rate problem. This means that Eq. (15) has
several, nontrivial solutions.

One may look for harmonic solutions to this equation in terms
of harmonic perturbations: $\dot{u}(x)=A \exp(-i\omega x)$ where
$A=\text{unknown constant; and } \omega=\text{angular frequency. Substitution of}$
>this expression of the displacement rate into Eq. (15) yields the
following characteristic equation:

$$
A \omega^2 \exp(-i\omega x) [(1-d_0)E - E\epsilon_0 f(d_0) \bar{\epsilon}_0^\beta \bar{\psi}(l_c, \omega)] = 0 \quad (16)
$$

in which $\bar{\psi}(l_c, \omega) = \text{Fourier transform of the normalized weight function. A similar result was obtained by Pijaudier-Cabot and Benallal (1993) in a more general [three-dimensional (3-D)] context. Eq. (16) shows that bifurcation occurs if}$

$$
\bar{\psi}(l_c, \omega) = \frac{(1-d_0)E}{E\varepsilon_0 f(d_0) \bar{\varepsilon}_0^\beta} \quad (17)
$$

For a fixed value of internal length, this equation provides the
angular frequency, or wavelength, of the localized mode. Note
that if the internal length vanishes, which is the case in a standard
local damage model, the wavelength of the localized mode is
arbitrary when Eq. (16) is satisfied. It follows that the boundary
value problem has an infinite number of solutions and is ill posed.
Same as in monotonic loading, this is a sufficient reason for in-
troducing an internal length in the model, since one cannot be
sure that the localization condition will ever be satisfied for any
arbitrary fatigue load history.

In the course of monotonic loading, bifurcation occurs at the onset of softening, at peak stress. The major difference between fatigue and monotonic loading is that the amplitude of the strain remains small during fatigue cycles. The denominator of the right hand-side term in Eq. (17) is always very small. The result is plotted in Fig. 5, which shows when localization may occur in the cases of fatigue and monotonic load histories. The values of the parameters in the evolution law for damage are provided in Table 1. We have considered in Fig. 5 the local damage case ($l_c=0$), which corresponds to the situation where localization occurs soonest in the course of the loading history (see Pijaudier-Cabot and Benallal 1993). As far as the bifurcation criterion is concerned, taking a nonzero value of the internal length would not change the comparison between monotonic and fatigue load histories substantially. The localized deformation modes would, however, be very different. Localization occurs when the load history (strain versus damage curve) intersects the bifurcation cri-

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Evolution of localization criterion during monotonic and fatigue uniaxial tests</figcaption></figure>

terion. For a monotonic increase in strain, it occurs for
ε₀ ≈ 8 × 10⁻³ and for a value of initial damage of d₀ ≈ 0.2. For
fatigue cycles with amplitude of εₐ = 140 × 10⁻⁶, bifurcation oc-
curs when damage is close to 1. If the amplitude increases, bifur-
cation occurs for smaller and smaller values of initial damage in
the fatigue problem. Strain localization, at least in uniaxial tests,
is almost avoided. Then, one might ask about the requirement of
introducing an internal length in the model whose aim is to re-
store well posedness of the boundary value problem upon local-
ization. Again, an internal length is still needed since it cannot be
proved that for any general loading localization will never occur.
The model should encompass fatigue loading histories with very
different amplitudes, including large ones for which bifurcation
due to softening will occur. In this last case, without an internal
length, failure without energy dissipation might be encountered
upon localization, the same as in usual monotonic loading situa-
tions. This feature, along with the ill posedness of boundary value
problems at the onset of localization would render finite-element
computations mathematically and physically meaningless.

Influence of Internal Length

Besides control of the localized mode at the onset of localization,
which results from Eq. (17), the internal length still has an influ-
ence on the growth of damage in boundary value problems where
the state of strain is not homogeneous. In order to exhibit this
influence, we consider the case of a cantilever beam subjected to
sinusoidal displacement of fixed amplitude at its free end. Fig. 6
shows the problem and the finite-element (FE) mesh that has been
used for the plane stress computations. The model parameters are
reported in Table 1. Fig. 7 shows the evolution of the overall
stiffness of the beam as a function of the number of cycles for
several values of internal length: l_c = 0, 9, 18, and 36 mm. When-
ever the state of strain is not homogeneous, there is a well known
sensitivity of the model predictions to the internal length (see,
e.g., Pijaudier-Cabot and Bazant 1987; de Borst et al. 1993). In a

<table><caption>**Table 1.** Model Parameters for Computation of Cantilever Beam for Several Internal Lengths</caption>
  <thead>
    <tr>
      <th>&alpha;<sub>1</sub></th>
      <th>&alpha;<sub>2</sub></th>
      <th>&alpha;<sub>3</sub></th>
      <th>&beta;</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>5 &times; 10<sup>-19</sup></td>
      <td>0.5</td>
      <td>3.0</td>
      <td>5</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Cantilever beam subjected to high cycling fatigue. Geometry of the specimen (a) and finite-element mesh (b). The specimen is clamped at the bottom and the upper end is subjected to sinusoidal horizontal displacement.</figcaption></figure>

monotonic loading problem, the peak load increases and the softening slope decreases when the internal length increases. In a fatigue problem, Fig. 7 shows that the larger the internal length, the larger the fatigue life time because the strain, which is not homogeneous, is averaged over a larger area as this length is increased (consequently the maximum nonlocal equivalent strain decreases).

It should be stressed that such sensitivity cannot be used for calibrating the internal length from test data on a single specimen size. In fact, it has been shown for monotonic loads that there is an infinite set of model parameters, with arbitrary values of inter-

<figure><img src="image_10.png" /><figcaption>Fig. 7. Influence of internal length on evolution of fatigue damage during cantilever test</figcaption></figure>

nal length that can fit a single structural response (three-point bend tests) with the same very good accuracy (Le Bellégo et al. 2003). To arrive at proper calibration of the model parameters, i.e., a unique set of model parameters, including the internal length, tests on several specimens subjected to different boundary conditions, with different geometries, or size effect tests must be available (Carmeliet 1999). Another possibility is to use an approximate value of the internal length. It is this second possibility that will be used next.

## Calibration and Comparisons with Test Data

Bituminous materials exhibit viscoelastic behavior, which leads to energy dissipation during each loading cycle of the fatigue tests. This dissipated energy is transformed into heat, leading to an increase of temperature inside the specimen. Dissipated energy, thermal exchange at specimen boundaries, and thermal properties of the material introduce a coupled thermomechanical effect. In uniaxial tension, compression fatigue tests, for instance, a 1°C rise may cause 5% loss of stiffness for usual conditions of temperature and frequency. We will evaluate the influence of thermal effects in a simplified uncoupled way. Although the viscoelastic response of the material is not explicitly modeled in this work (since we look at variation of the complex modulus), its consequences in terms of energy dissipation and thermomechanical effects will be taken into account for a proper comparison with experiments.

## Evaluation of the Temperature Increase during Fatigue Tests

Let us consider that the structural stiffness of the specimen $K(N)$ measured at cycle $N$ can be written as

$$ \frac{K(N)}{K(0)} = (1 - D_{\text{th}})(1 - D_{\text{mech}}) \quad (18) $$

where the variations in specimen stiffness due to thermal effects and mechanical damage are denoted $D_{\text{th}}$ and $D_{\text{mech}}$, respectively. Note that these variations result from a decrease of material stiffness due to the increase of temperature and from the distribution of damage $d$ in the specimen. Evaluation of the influence of thermal effects $D_{\text{th}}$ will now be performed independently from the calculation of mechanical damage.

This coupled thermomechanical problem without any damage has been computed with the finite-element program César-LCPC, and the results have been validated by infrared thermography measurements for bituminous mixtures in the past (de La Roche et al. 1998). The viscoelastic model used to describe bituminous mixture behavior is from work by Huet (1963) and Sayegh (1967). The model uses the “equivalence principle” between the frequency $f$ and the temperature $\theta$:

$$ E(f, \theta) = E[f, \tau(\theta)] = E_{\text{real}} + iE_{\text{imag}} = E_0 + \frac{E_{\infty} - E_0}{1 + \delta(ift)^{-k} + (ift)^{-h}} \quad (19) $$

$E_0$ and $E_{\infty}$ = limit of the complex modulus when $f$ tends towards 0 and infinity, respectively. $k$ and $h$ exponents (such that $1 > h > k > 0$) are related to the behavior of the ratio $E_{\text{imag}}/E_{\text{real}}$ when $f$ tends towards infinity (resp. $f$ tends towards 0). $\delta$ = a dimensionless constant. Function $\tau = \tau(\theta) = \text{the characteristic time, a function of the temperature, which Huet and Sayegh suggested approximating with Arrhenius or Eyring type law. } \log(\tau)$ is usually defined as the following parabolic function:
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Increase in sample surface temperature as a function of the number of cycles: experimental data and FE simulations for tension-compression tests (after Bodin et al. 2003)</figcaption></figure>

$$ \log(\tau) = A_0 + A_1 \theta + A_2 \theta^2 \quad \text{for } \theta \text{ in the interval } (-10, 50^{\circ}\text{C}) \qquad (20) $$

For the mix studied, the model parameters are the following: $E_{\infty}$=29,700 MPa, $E_0$=25 MPa, $h$=0.675, $k$=0.21, $\delta$=2.16, $A_0$ = 1.71, $A_1$ = -0.393 and $A_2$=0.00226. For the fatigue test conditions (10°C, 10 Hz), the complex modulus of the mix is 11,600 MPa and its phase angle $\phi$ is 15.1°. The Poisson ratio is considered to remain constant and =0.35. The thermal constants are as follows: heat conductivity 0.8 $10^6$ W/m/°C and volumetric heat capacity 1.6 J/m$^3$/°C. The exchange coefficients that give the heat flow at the boundaries of the specimen are determined in order to reach good agreement between measured and calculated temperatures at the surface of the specimen. The results of the simulation on tension compression tests are compared to thermal gauge measurements at the surface of the samples and presented in Fig. 8 (for more details also see Bodin et al. 2003).

Table 2 summarizes the relative decrease of structural stiffness computed for uniaxial tension compression tests and for bending tests on cantilever beams at several amplitudes. It can be observed that the variation in the structural stiffness of the specimen due to the thermal effect is almost negligible in bending tests, and slightly higher in tension compression tests. It has to be pointed out that the thermal calculations did not integrate the decrease in stiffness due to damage that occurred during fatigue tests. This assumption of noncoupled effects between damage and temperature is less and less true as damage grows and energy dissipation due to damage is converted into heat.

## Calibration of Fatigue Test Data

Calibration of the evolution law for damage, Eqs. (5), and (9), can be performed directly in tension compression tests and bending tests with small amplitudes ($\epsilon_a$=80×$10^{-6}$ and

<table><caption>**Table 2.** Maximum Correction Factors due to Thermal Effects</caption><thead><tr><th>Variable</th><th colspan="4">Tension compression</th><th colspan="3">Two-point bending</th></tr></thead><tbody><tr><td>ε<sub>max</sub> (10<sup>-6</sup>)</td><td>80</td><td>100</td><td>150</td><td>180</td><td>140</td><td>180</td><td>220</td></tr><tr><td>D<sub>th</sub> (10<sup>-2</sup>)</td><td>2.05</td><td>3.2</td><td>6.4</td><td>10.7</td><td>1.7</td><td>2.4</td><td>3.9</td></tr></tbody></table>

$140 \times 10^{-6}$ m/m, respectively). Strain localization analysis shows that it is possible to use the tension-compression tests under the assumption that the strain distribution is homogeneous throughout the specimen, since bifurcation occurs for values of damage that are close to 1. Unfortunately, the set of available experimental data in uniaxial tension does not exhibit the third regime of damage growth. The test control and the displacement transducers did not allow data acquisition during this very fast regime so experimental measurements are lacking. Hence, this regime of damage growth, which is important for calibration of parameter $\alpha_1$, is not accessible for the calibration from uniaxial tests and other experiments; here bending beams subjected to fatigue are required.

With the help of the Matlab standard least square fitting function, the fit of the uniaxial experiments, for the smallest amplitude after thermal correction is performed, yields an initial set of material parameters. Uniaxial tests are used because integration of the evolution of damage can be performed analytically [Eqs. (7) and (9)] and standard fitting techniques, which do not need any finite-element calculation, can be implemented without any difficulty.

This set of model parameters serves as initial input for finite-element simulation of bending tests and is corrected in order to capture the third regime of damage growth properly. To do this, the sensitivity analyses performed earlier (Fig. 4) are helpful. Parameter $\alpha_3$ is kept equal to its value obtained from tension compression tests. Since parameter $\alpha_1$ has a direct influence on fatigue life (Fig. 4), it is modified so that the computations provide the correct number of cycles up to complete failure. In order to balance the influence of the variation of $\alpha_1$, $\alpha_2$ is adjusted to keep the correct evolution of damage during phases I and II of damage growth (Bodin 2002). Fig. 9 shows the result of this calibration in the bending and uniaxial tests. The set of model parameters obtained is given in Table 3. For the asphalt mix studied, the internal length $l_c$ was chosen as equal to three times the maximum aggregate size, the same as for concrete (Bazant and Pijaudier-Cabot 1989). Since thermal effects are very small in the bending test (see Table 2), they are neglected.

In addition to this fit, the gray zone in Fig. 9 delineates lower and upper bound predictions that are a result of a fit of the extreme values in uniaxial tension compression [Fig. 9(b)]. Compared to the previous fit, which provides average specimen responses (Table 3), parameters $\alpha_1$ and $\alpha_2$ were changed significantly by keeping the third model parameter almost constant. These lower and upper bounds show how sensitive model responses in the bending problem are, depending on inherent dispersion of uniaxial tests. It appears that bending tests are more appropriate for model calibration. Dispersion of the test data is smaller and data up to complete failure are more easily accessible. Such calibration relies on inverse analysis and finite-element computations. It is therefore more difficult to perform, and it still requires an initial set of model parameters to start an inverse analysis which reaches an optimal fit that can be obtained from uniaxial tests.

The simulation of uniaxial tests at amplitude $100 \times 10^{-6}$ m/m is shown in Fig. 10, including the responses obtained for the upper and lower bound fits discussed above. The correction due to thermal effects was reported in order to exhibit the influence of the temperature which is rather limited in view of scattering of the experimental data. Overall, this simulation provides good agreement with test data.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Calibration of model parameters from cantilever bending tests (top) and upper and lower bounds from uniaxial test for amplitude of 80×10⁻⁶ m (bottom); gray zone shows lower and upper bound fits</figcaption></figure>

## Conclusion

An elasticity based nonlocal damage model for fatigue failure of asphalt concrete was presented. Compared to the existing model of Paas et al. (1993), a different form of evolution law for damage was implemented in order to capture more accurately the third regime of damage growth, just prior to complete failure. This constitutive relation was implemented in a finite-element code, along with a jump-in-cycle method for finite-element computations of high cycle fatigue tests.

Strain localization analysis in a one-dimensional setting shows the difference between the failure processes in monotonic loading and fatigue loading. In monotonic loading, bifurcation occurs at peak stress, while it occurs very close to complete failure (damage close to 1) in fatigue tests. Therefore, according to the proposed constitutive relations, the third regime of damage, with an increasingly growing rate of damage near failure, is not solely due to structural effects inherent to localization of damage.

<table><caption>**Table 3.** Set of Model Parameters Obtained from Calibration</caption><thead><tr><th>Limit</th><th>β</th><th>α<sub>1</sub> (10<sup>-19</sup>)</th><th>α<sub>2</sub> (10<sup>-1</sup>)</th><th>α<sub>3</sub></th></tr></thead><tbody><tr><td>Upper</td><td>5</td><td>4.73</td><td>5.04</td><td>3.07</td></tr><tr><td>Average fit</td><td></td><td>6.84</td><td>5.0</td><td>3.00</td></tr><tr><td>Bottom</td><td></td><td>15.8</td><td>5.43</td><td>3.10</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Simulation of tension-compression tests, including corrections due to the thermal effect; the gray zone indicates lower and upper bound predictions</figcaption></figure>

An internal length is still introduced to ensure the consistency of the constitutive relations upon localization of damage, which may still occur when the amplitude of loading is very high as, for instance, in the case of notched specimens. The internal length controls the size of the damage localization zone, and the larger this parameter, the longer the fatigue life of bending beams.

The model parameters were calibrated on bending beams and uniaxial test data. Prior to calibration, the thermal part of the decrease in stiffness due to thermomechanical effects was evaluated by a simplified approach. For the material studied and laboratory fatigue conditions tested (10°C, 10 Hz), the variation of stiffness due to thermal effects remains under 4% of the structural stiffness of the specimens when the maximum strain value is below 100×10⁻⁶ for tension compression tests and 220×10⁻⁶ for cantilever tests. The model provides reasonable simulations of tension-compression tests for several amplitudes of loading.

It should be pointed out that, in the calibration process, the internal length was chosen as a function of the maximum aggregate size, following experimental data obtained for concrete. This estimate of the internal length needs to be further justified, since the microcracking process in bituminous mixes may be quite different from that in concrete. Size effect tests on geometrically similar cantilever beams should be able to provide a better estimate. Since the internal length controls the size of the damage zone, the response of geometrically similar beams should not be geometrically similar. Inverse analysis of such size effect tests should provide a more accurate determination of this model parameter.

## Acknowledgments

The experimental results presented were performed within the framework of a scientific collaboration between the Laboratoire Central des Ponts et Chausseés (LCPC) and ENTPE, France. They are part of an international interlaboratory program for fatigue tests of bituminous mixtures conducted under the auspices of RILEM, TC 180 PEB. This work was part of LCPC research project Fatigue Damage of Bituminous Mixtures.
References

Alimami, M. (1987). "Contribution à l'étude de l'endommagement par fatigue des enrobés bitumineux." PhD thesis (Eng.), Univ. de Paris 6, Paris.

Bazant, Z. P., and Pijaudier-Cabot, G. (1989). "Measurement of characteristic length of nonlocal continuum." *J. Eng. Mech.*, 115(4), 755-767.

Bodin, D., de La Roche, C., Piau, J. M., and Pijaudier-Cabot, G. (2003). "Prediction of the Intrinsic damage during bituminous mixes fatigue tests." *Int. RILEM Symp. on Performance Testing and Evaluation of Bituminous Materials*.

Bodin, D., Pijaudier-Cabot, G., de La Roche, C. and Chabot, A. (2002). "A continuum damage approach of asphalt concrete fatigue tests," *Proc., 15th ASCE Engineering Mechanics Conf., Reston, Va.*

Breysse, D., de La Roche, C., Chauvin, J.-J., and Domec, V. (2002). "Fatigue tests on bituminous composites: Balance between damage and recovering." *Proc., 15th ASCE Engineering Mechanics Conf., Reston, Va.*

Carmeliet, J. (1999). "Optimal estimation of gradient damage parameters from localisation phenomena in quasibrITTLE materials." *Int. J. Mech. Cohes. Frict. Mater.*, 4, 1-16.

de Borst, R., Sluys, L. J., Muhlhaus, H. B., and Pamin, J. (1993). "Fundamental issues in finite element analysis of localisation of deformation." *Eng. Comput.*, 10, 99-121.

de La Roche, C., and Marsac, P. (1996). "Caractérisation expérimentale de la dissipation thermique dans un enrobé bitumineux sollicité en fatigue." *Proc., Euroasphalt and Eurobitume Congress*.

de La Roche C., Piau, J. M., and Dangla, P. (1998). "Thermal effects induced by viscoelastic dissipation during fatigue tests on bituminous mixtures." *Proc., 6th Int. Symp. On Creep and Coupled Processes*.

Di Benedetto, H., de La Roche, C., and Francken, L. (1997). "Fatigue of bituminous mixtures: Different approaches and RILEM interlaboratory tests." *Mechanical Tests for Bituminous Materials*—Recent Improvement and Future Prospects, *Proc., 5th RILEM Symp.* Balkema, Rotterdam, The Netherlands, 15-26.

Fichant, S., La Borderie, C., and Pijaudier-Cabot, G. (1999). "Isotropic and anisotropic descriptions of damage in concrete structures." *Int. J. Mech. Cohes. Frict. Mater.*, 4, 339-359.

Huet, C. (1963). "Etude par méthode d'impédance du comportement viscoélastique des matériaux hydrocarbonés." PhD thesis (Eng.), Faculté des Sciences de Paris, Paris.

Jirásek, M., and Patzák, B. (2002). "Consistent tangent stiffness for non-local damage models." *Comput. Struct.*, in press.

Kim, Y. R., and Little, D. N. (1990). "One-dimensional constitutive modeling of asphalt concrete." *J. Eng. Mech.*, 116(4), 751-772.

Le Bellégo, C., Dubé, J. F., Pijaudier-Cabot, G., and Gérard, B. (2003). "Calibration of nonlocal model from size effect tests." *Eur. J. Mech. A/Solids*, 22, 33-46.

Lee, H.-J., and Kim, Y. R. (1998). "Viscoelastic constitutive model for asphalt concrete under cyclic loading." *J. Eng. Mech.*, 124(1), 32-40.

Lee, H. J., Daniel, J. S., and Kim, Y. R. (2000). "Continuum damage mechanics-based fatigue model of asphalt concrete." *J. Mater. Civ.*

Eng., 12(2), 105-112.

Lefeuvre, Y. (2001). "Etude expérimentale et modélisation de la fissuration de matériaux bitumineux." PhD thesis, Ecole Nationale des Ponts et Chaussées, Paris, France.

Lefeuvre, Y., de La Roche, C., and Piau, J. M. (2000). "Asphalt material fatigue test under cyclic loading: An interpretation based on a viscoelastic model including unilateral damage." *Proc., 2nd Eurasphalt and Eurobitume Congress*.

Lemaitre, J. (1996). *A course on damage mechanics*, 2nd Ed., Springer, Berlin.

Mazars, J. (1984). "Application de la mécanique de l'endommagement au comportement nonlinéaire et à la rupture des bétons de structure." PhD thesis d'Etat, Univ. Pierre et Marie Curie, Paris.

Paas, R. H. J. W. (1990). "Continuum damage mechanics with an application to fatigue." PhD thesis, Technishe Universiteit Eindhoven, Eindhoven, The Netherlands.

Paas, R. H. J. W., Scheurs, P. J. G., and Brekelmans, W. A. M. (1993). "A continuum approach to brittle and fatigue damage: Theory and numerical procedures." *Int. J. Solids Struct.*, 30(4), 579-599.

Partl, M. N., Sybilsky, D., Sousa, J. B., Di Benedetto, H., and Piber, H. (2000). "Current research projects of the next RILEM TC PEB performance testing and evaluation of bituminous materials." *Proc., 2nd Eurasphalt and Eurobitume Congress*.

Patzák, B., and Jirásek, M. (2001). "Consistent tangent stiffness for non-local material models." *European Conf. on Computational Mechanics*, ECCM-2001.

Peerlings, R. H. J., Brekelmans, W. A. M., de Borst, R., and Geers, M. G. D. (2000). "Gradient-enhanced damage modelling of high-cycle fatigue." *Int. J. Numer. Methods Eng.*, 49(12), 1547-1569.

Piau, J. M. (1989). "Modélisation thermomécanique du comportement des enrobés bitumineux." *Bulletin de Liaison des Laboratoire des Ponts et Chaussées*, 163, 41-54.

Pijaudier-Cabot, G. (1995). Nonlocal damage, continuum models for materials with microstructure, Chap. 4, H. B. Muhlhaus, ed., Wiley, New York, 105-144.

Pijaudier-Cabot, G., and Bazant, Z. P. (1987). "Nonlocal damage theory." *J. Eng. Mech.*, 113(10), 1512-1533.

Pijaudier-Cabot, G., and Benallal, A. (1993). "Strain localization and bifurcation in a nonlocal continuum." *Int. J. Solids Struct.*, 30, 1761-1775.

Sayegh, G. (1967). "Viscoelastic properties of bituminous mixes." *Proc., 2nd Int. Conf. on the Structural Design of Asphalt Pavement*.

Shapery, R. A. (1984). "Correspondence principles and a generalised *J* integral for large deformation fracture analysis of viscoelastic media." *Int. J. Fract.*, 25, 195-223.

Strategic Highway Research Program (SHRP). (1994). "Fatigue response of asphalt aggregate mixes." National Research Council Rep. No. *SHRP-A-104*, Washington, D.C.

Ullidtz, P., Kieler, T. L., and Kargo, A. (1997). "Finite element simulation of asphalt fatigue tests, mechanical tests for bituminous materials—Recent improvement and future prospects," *Proc., 5th RILEM Sympo.*, Balkema, Rotterdam, The Netherlands, 233-239.