Numerical modelling of mode-I cracking in asphalt concrete

This paper presents a numerical model to simulate mode-I fracture in asphalt concrete. The phase field method is used to simulate diffuse damage in asphalt concrete. In addition to the displacement degrees of freedom, a scalar damage degree of freedom is added at each node of a finite element. Evolution of the damage variable is governed by a crack potential which regularizes a sharp crack over a finite volume. The phase field method does not suffer from mesh dependency problem like the continuum damage models and therefore preferred for damage analysis. A single edge notch beam is analysed to show the performance of the model. The damage pattern as well as load versus crack mouth opening displacement graph is compared with the experimental data. The numerical results show good agreement with the experimental observations.


Introduction
Asphalt pavements are subjected to damage over their service life. The damage in asphalt pavements can be due to environmental effects and/or traffic loads. A damaged pavement increases the maintenance cost and reduces the service life of the pavement. For an efficient design of pavement structures, development of an accurate and reliable numerical tool is pivotal.
Different computational models have been developed in the past to study the damage mechanics of asphalt mixtures. Cohesive zone model in combination with interface elements is a highly popular method to simulate damage in asphalt mixtures. [1,2] used a cohesive zone model with interface elements to simulate mode-I fracture in asphalt concrete at low temperatures. A potential based exponential cohesive constitutive law was used to model the inelastic material response near a crack tip region. The model successfully simulated mode-I fracture response of asphalt concrete. [3] performed a numerical analysis to study reflective cracking in asphalt pavements using a cohesive zone model with a linear softening constitutive law. [4] used a cohesive zone model with interface elements to simulate reflective and thermal induced cracking in asphalt pavements. The advantage of interface elements is that it can simulate fracture behaviour of a material however, the limitation is that the cracks are only allowed to propagate along finite element boundaries. Therefore, realistic damage behaviour cannot be obtained. The limitation of a crack to grow along element boundaries is eliminated using advanced finite element methods like the extended finite element method. [5] used extended finite element method to simulate multiple cracks in an asphalt pavement. Other researchers e.g. [6,7] successfully used extended finite element method in combination with visco-elastic bulk constitutive law to simulate damage in asphalt concrete and hot-mix asphalt mixtures. However, the damage behaviour of asphalt mixture was assumed to be brittle and inelastic zone ahead of crack tip was assumed to be very small. Moreover, even though extended finite element method can simulate multiple mesh independent cracking in a solid, simulation of multiple interacting cracks and its computer implementation is complicated and burdensome.
Apart from a fracture mechanics-based approach like the cohesive zone model, [8] presented a continuum mechanics-based approach to simulate damage in asphalt mixtures. A damage phase field variable was defined to separate a fully damaged and undamaged part of the material through a diffuse interface. The method does not require an algorithm to track the crack path unlike interface element and extended finite element method and is easy to implement in a computer program. A double-well potential function was used to define the crack driving force. However, such a function may give unrealistic damage. Moreover, the damage may grow perpendicular to the actual crack path [9]. [10] presented a phase field model which uses a generalized potential function to simulate brittle and quasi-brittle material failure.
The aim of this work is to present a diffuse damage model for mode-I fracture in asphalt concrete. The phase field model of [10] is used to simulate damage. A damage degradation function with linear and exponential softening behaviour is used to simulate quasi-brittle fracture in asphalt concrete. It is also the aim of this work to study the effect of different softening curves on material response. Within the framework of finite element method, a damage phase field variable in addition to the displacement field is defined to simulate damage in the bulk material.
Remainder of the manuscript is organized as follows: Section 2 briefly discusses the phase field damage model. Details of the analysed numerical problem, i.e. a notched beam, are discussed in section 3 whereas section 4 presents analysis results. Section 5 briefly summarizes main conclusions of the manuscript.

Phase field model
Consider a solid body with domain Ω, subjected to prescribed displacements ū on the boundary ∂Ω u and prescribed tractions ̅ on the boundary ∂Ω t, figure 1 . Volumetric forces is acting over the domain Ω. Additionally, the body also contains an internal boundary of a sharp crack S, which is regularized over a localization band B of width b. The spaces of displacement field U u and its small variation V u are defined as The spaces of damage phase field variable U d and its small variation V d are defined as (2) in which u and d are the unknown displacement and damage variables, respectively. The damage variable d varies within a localization band B such that d=0 and d=1 represents undamaged and fully damaged material.
The total potential energy Π of the system is given as (3) in which c 0 is the model parameter which is required to regularize a sharp crack into a diffuse crack, G f is the material fracture toughness, α is the crack geometric function. In this work α = 2d -d 2 is used. It is worth mentioning that [9] used α = 2d for simulating damage in brittle materials whereas [11] used α = d 2 . The advantage of using function α = 2d -d 2 is that damage develops within a finite domain which is numerically useful. Another, advantage of finiteness of the phase field domain is that we can apply the phase field model to a small portion of the whole structure/body where the damage is more likely to grow and rest of the body can be modelled as an elastic material. This reduces the computational cost. Ψ 0 is the elastic strain energy density function. The Cauchy stress tensor σ is defined as σ = ωσ o with . D is the elastic material stiffness tensor. П ext is the potential energy of external forces. ω is the damage degradation function which is given as [10] (4) in which the parameters , , are determined to simulate a particular softening behaviour similar to the cohesive constitutive laws used in cohesive zone models.
Minimization of the total potential, equation (1), yields the following two coupled variational equations which are the solved using finite element method (5) (6) in which is the effective crack driving force given as (7) Equation (7) basically defines the damage initiation criteria. Accordingly, damage initiates when local tensile principal stress reaches the tensile strength of the material. Equations (5) and (6) are discretised using standard finite element shape functions and are solved for the unknown displacement and damage phase field .

Three-point bend test on a notched beam
A three-point bending test is performed to see the performance of the phase field model. The beam also contains a notch at its mid-length to simulate mode-I fracture. Geometry and boundary conditions of the notched beam model is shown in figure 2. The beam was experimentally tested by [12]. The beam is made of asphalt concrete and is tested at a temperature of -10 o C to study fracture behaviour of an asphalt concrete at low temperatures. In the numerical model a notch of 19mm height is introduced at the mid-length of the beam. Notch width is taken to be 2mm. The geometry of the beam is discretised with four node quadrilateral elements with full integration. A fine mesh with a maximum element size of 1mm is used near the notch region where the damage is more likely to grow whereas a coarse mesh with a maximum element size of 5mm is used elsewhere. A load control analysis is performed by linearly increasing the load with an increment of 0.001mm. An in-house C++ finite element code with the phase field model is used for analysis whereas MATLAB is used for post-processing. Material properties used for the analysis are taken from [13]. Young's modulus E o =14.2GPa, Poisson ratio ν=0.35, material tensile strength f t =3.4MPa, fracture toughness G f =0.344N/mm. In order to investigate the effect of softening law, the analysis is performed with two softening laws, i.e. the bilinear and exponential softening laws, figure 3. Phase field model parameters to simulate the two softening behaviours are given in table 1.

Analysis results and discussion
The phase field analysis results are compared with the experimental load versus crack mouth opening displacement (CMOD) of [13] in Figure 4. It is observed from the figure that the initial stiffness of the single edge notch beam is well predicted by the phase field model both with the linear and exponential softening laws. Moreover, the post peak response is also in good agreement with the experimental result.    [13]. [13] used interface elements with a potential based exponentially decaying cohesive law to simulate model-I cracking in a single edge notch beam. The initial pre-peak branch of the curve shows a compliant behaviour compared to the experimental curve whereas the post peak branch of the curve is in good agreement with the experimental results. This mismatch in the pre-peak branch of the curve is due to the complaint initial stiffness of the underlaying cohesive law. On the other hand, the phase field model is able to capture the initial high stiffness branch of the load versus CMOD curve both with the linear and exponential softening law. The reason is that in the phase field model the material remains intact until the failure stress is approached. Whereas in the interface element method an initial stiffness of the cohesive law simulates rigid material behaviour before the initiation of damage. It is for this reason analysis results with interface elements are dependent upon the initial stiffness of the cohesive constitutive law. Comparing the phase field analysis results of linear and exponential softening laws, it is evident from figure 4 that the peak load is underestimated when exponential softening law is used. This is due to rapid softening of the material after the peak stress when using exponential softening law. Consequently, stress free crack develops quickly.
The effect of load drop is also evident from figure 5. Figure 5 shows a zoom-in view of the damaged zone and presents damage progression in a notch beam at different levels of prescribed displacements u. At u=0.14mm, the overall damaged zone for the case of exponential softening is (a) narrow near the tip (b) more elongated and (c) shows more damage value (shown by the orange colour in the centre) compared to the case of linear softening. This initial increased damage in case of exponential softening results in a lower peak load compared to the analysis with linear softening law using the same material parameters. This is obvious because in case of exponential softening the stress decays exponentially and faster after damage initiation. Moreover, the width of the damaged zone is more for the case of analysis with linear softening compared to the analysis with exponential softening, figure 5. This is again due to the nature of softening curves. Exponential softening causes rapid material degradation and therefore material does not have a time for stress redistribution. On the other hand, linear softening causes relatively slow material degradation therefore stress redistribution occurs and more area surrounding the crack is stressed resulting in broader width of localization/damage zone. This is also evident from the damage profiles at u=0.49mm. The fully damaged zone is broader in case of linear softening.

Conclusions
A phase model to simulate mode-I fracture in asphalt concrete is presented. Fracture is modelled as a diffuse damage like in continuum damage models. A damage evolution function considering linear and exponential softening behaviour is used to simulate quasibrittle behaviour of asphalt concrete at low temperature. The numerical results show that the phase field model can simulate model-I fracture in asphalt concrete. Moreover, the choice of softening curve in a modelling process is crucial as the shape of the softening curve affect the post-peak behaviour of the structure. Based on the numerical results it is concluded that the type of softening law influences the peak load and width of the localization zone.