2D Mesoscale cracking simulation of partially saturated asphalt based on moisture diffusion and a cohesive zone model

ABSTRACT

The primary objective of this paper was to develop a combined model that incorporates moisture diffusion and a cohesive zone model, addressing anisotropic and loading-rate dependent cracking within partially saturated asphalt. Utilising X-ray CT scan, cross-sectional slices of asphalt were acquired and converted into vector images through Matlab and AutoCAD, forming a digital asphalt sample. Moisture concentration in the asphalt, after different immersion durations, was quantified by Fick's law. A sequentially coupled model of moisture diffusion and fracture investigated the effect of immersion duration, anisotropy, and loading rate on cracking performance of the asphalt during a digital indirect tensile strength test (DITST) at 5°C. Findings revealed that moisture evolution in partially saturated asphalt proceeds through two or three stages: near-zero growth (only applicable for locations far from the initial moisture-asphalt interface), rapid growth, and a plateau. Peak load, stiffness, and fracture work in DITST exponentially reduced with immersion duration, predominantly within the first four weeks. Anisotropy led to differential DITST results when varying loading direction. Moisture damage decreased crack resistance across all directions, while increasing loading rate enhanced it. Fracture stiffness and strength exhibited comparable impacts on cracking performance at a specific loading rate.

KEYWORDS:

• Asphalt
• cracking
• cohesive zone model
• moisture diffusion

1. Introduction

Moisture damage to asphalt pavements is one of the main distress types in the UK due to the regular occurrence of wet weather throughout the year. The durability of an asphalt pavement typically depends on the interface properties (i.e. adhesion and cohesion) between the aggregate and the bitumen or within the bituminous mastic (Cui et al. 2014, Azarhoosh et al. 2017). One of the documented critical factors contributing to the deterioration of these interface properties is the presence of water in the asphalt layer through diffusion or other transport methods (Huang and Luo 2022), although traffic volume, environmental temperature, construction quality etc. also affect the adhesive and cohesive failures in the pavement. A cohesive zone model (CZM) is commonly used to characterise the fracture properties of the asphalt due to its effective simplification of stress distribution at the crack tips, which practically avoids convergence issues in its application to finite element analysis. Caro et al. (2010a) investigated moisture-induced damage in a tiny element of asphalt (a rectangular sample consisting of asphalt mastic and two coarse aggregates) based on CZM. Ban et al. (2013) also employed CZM to evaluate the effects of air voids, diffusion coefficient, and degradation characteristics on the moisture damage of an asphalt beam. The adhesive and cohesive properties of asphalt manifest strong rate dependence (attributed to bulk and interface viscosity of the bitumen) and anisotropy (depending on the size, orientation, and sphericity of the aggregate), but how these properties influence the cracking resistance of asphalt is not yet well known. This uncertainty will be substantially added to when moisture damage is included.

At the mesoscale, asphalt is a type of heterogeneous material and can be regarded as a composition of fine asphalt mixture (simplified as a homogeneous viscoelastic material consisting of bitumen, fine aggregate, filler, and additives), air voids, and coarse aggregates. Chen et al. (2022a) used a high-resolution X-ray CT to characterise the fatigue damage of asphalt under the indirect tensile fatigue test. They found that fatigue damage changed the concentration of air voids, most of which were observed in the centre part of the asphalt sample. Fatigue cracks emerged in and around elongated aggregates, which implied that the aggregate angularity and shape had critical impacts on crack formation. Shi et al. (2021) developed a pore cellular structure model to assess the mesostructure of asphalt. They defined parameters (termed ‘interface coefficient’ and ‘skeleton rate’) that allowed effective mesostructure evaluation and found that asphalt gradation was crucial for the characterisation of contact points, and area of pore cell distribution. Zhao et al. (2020b) proposed a method to enhance the accuracy of 2D mesoscale modelling of an indirect tensile strength test of asphalt. Based on their correlation between the simulation results and aggregate content, they found that aggregate content is a good indicator allowing the selection of the most appropriate 2D CT scan image to best represent a 3D asphalt sample concerning its cracking performance. The focus of these referenced works primarily lay on the fracture features and traits of the asphalt under dry conditions. Given the application of asphalt in real-world environments, the pertinence of these fracture characteristics, particularly under the influence of environmental factors such as moisture, warrants further deliberation.

That there are two types of moisture damage mechanisms (adhesive damage and cohesive damage) is widely accepted. These are fundamentally due to (1) the emergence of diffusible moisture molecules at the interface between the bitumen and aggregate and within the bitumen; (2) chemical and physical reactions between the active alkali mineral components and the water molecules, and the inverted emulsion of water in the bitumen. Apeagyei et al. (2015) investigated certain aspects of asphalt mastic moisture diffusion characteristics using Fick’s law and non-Fick’ law, they found that both of them delivered excellent results with R² over 95%. Arambula et al. (2010b) numerically analysed the moisture diffusion in asphalt based on its digital image and Fick’s law. They concluded that Fick’s law worked very well for the prediction of moisture distribution in the asphalt by comparing it against the experimental measures. Huang and Luo (2022) proposed an accumulative water vapour diffusion test using a gravimetric sorption analyser to develop a two-phase diffusion model composed of free and bound water molecules in the asphalt. Based on this model, they provided a new method to calculate the diffusivity of water vapour and the maximum mass of free and bound water. Zaidi et al. (2022) quantified the effect of hydrated lime on the moisture susceptibility of asphalt by surface energy measurement, and modified saturation ageing and tensile stiffness tests. It was found that the moisture susceptibility results presented a close agreement between the surface energy measurement and the saturation ageing and tensile stiffness test technique. The asphalt with 1 wt.% hydrated lime was observed to have better moisture resistance than that with 2 wt.% hydrated lime in the aggregates (felsic intrusive, and argillaceous sandstone). Castillo et al. (2017) developed a coupled moisture-mechanical model to evaluate moisture damage in an asphalt mixture. A methodology for the random generation of mesostructure and a coupled moisture-mechanical continuum damage constitutive equation were integrated to evaluate the effect of material heterogeneity, air voids content, and moisture diffusion coefficient on the moisture susceptibility of the asphalt. They concluded that (1) the non-uniform distribution of moisture content and damage can be attributed to asphalt heterogeneity (characterised by the random nature of mesostructure); (2) both the air voids content and moisture diffusion coefficient have an impact on moisture damage, and the latter was founded to be more profound. On a broad scale, our understanding of moisture-induced damage to asphalt mixtures is fairly exhaustive. Yet, they showcase pronounced rate-dependence and anisotropy. The influence of these attributes on the susceptibility of asphalt mixtures to moisture damage remains unknown.

The adhesive and cohesive characteristics of asphalt show a pronounced dependence on rate, which is tied to the viscosity of the bitumen itself, both in bulk and at the interface. These properties also display a level of anisotropy, influenced by aspects such as the aggregate's size and positioning. However, the way these factors affect the asphalt's resilience against cracking is not entirely understood as of yet. This conundrum is poised to deepen significantly once the effects of moisture-induced damage are taken into account. This paper aims to characterise the cracking performance of partially saturated asphalt based on the coupled modelling of moisture diffusion and CZM crack modelling at the mesoscale. The ultimate goal is to simulate (and therefore design against) surface-originating asphalt damage. It is essential to clarify that the term ‘partially saturated’, as used in this context, dedicatedly describes a phase of the asphalt where its moisture content straddles the states of dryness and saturation. The precise degree of moisture saturation hinges on several variables, notably the duration of immersion and the mesostructural attributes of the asphalt. In this paper, a 2D geometrical model using the central slice (in the height direction) of a reconstructed three-dimensional asphalt based on an X-ray CT scan and image processing techniques has been developed. A digital immersion test has been used to generate the moisture profile in the asphalt followed by a digital indirect tensile strength test to quantify the moisture-induced reduction in the crack resistance, as well as the effects of anisotropy and loading rate.

2. Modelling methodology

2.1 Viscoelastic modelling of fine aggregate mixture with finite strain

Fine aggregate mixture (FAM) is a composite material consisting of four different phases: bitumen, filler particles, fine aggregates, and embedded tiny air voids that cannot be captured by the X-ray CT. Hence, the rheological properties of FAM depend on the behaviours of the bitumen, filler, and fine aggregates, and their mass percentage. Unlike the full asphalt mixture, homogeneity and isotropy are commonly acceptable assumptions for the mechanical characterisation of FAM. Herein, the following constitutive model (i.e. isotropic linear viscoelastic model) with finite strain is employed to quantify the rate-dependence of the FAM: (1) $\{\begin{array}{c}S\left(t\right)=S_0\left(t\right)+dev\left[{\int }_0^t\dot{G}\left(t\right){\overline{F}}_t^{-1}\left(t-\tau \right)\cdot e\left(t-\tau \right)\cdot {\overline{F}}_t^T\left(t-\tau \right)d\tau \right]\\ {\sigma }_{vol}\left(t\right)={\sigma }_{vol}^0\left(t\right)+{\int }_0^t\dot{K}\left(\tau \right){ϵ}_{vol}\left(t-\tau \right)d\tau \end{array}$(1) where S and e are the deviatoric stress and strain tensors; σ_vol and ϵ_vol are the volumetric stress and strain; S₀ and σ0 vol are the instantaneous shear stress tensor and volume stress at time t; F_t(t-τ) is the distortional deformation gradient of the state at t-τ relative to the state at t; G and K are the shear and bulk relaxation moduli of the FAM, respectively; t and τ are reduced present time and reduced time history, respectively; and dev(·) = (·) – 1/3(·):I ⊗ I, with I standing for the identity tensor.

The shear and bulk relaxation moduli for the FAM shown in Equation (1(1) $\{\begin{array}{c}S\left(t\right)=S_0\left(t\right)+dev\left[{\int }_0^t\dot{G}\left(t\right){\overline{F}}_t^{-1}\left(t-\tau \right)\cdot e\left(t-\tau \right)\cdot {\overline{F}}_t^T\left(t-\tau \right)d\tau \right]\\ {\sigma }_{vol}\left(t\right)={\sigma }_{vol}^0\left(t\right)+{\int }_0^t\dot{K}\left(\tau \right){ϵ}_{vol}\left(t-\tau \right)d\tau \end{array}$(1) ) can be further represented by the following Prony series model (Zhang et al. 2016): (2) $\{\begin{array}{c}G\left(t\right)=G_{\mathrm{\infty }}+{\sum }_{m=1}^M{G}_m\exp \left(-{\displaystyle \frac{t}{{\rho }_m}}\right)\\ K\left(t\right)=K_{\mathrm{\infty }}+{\sum }_{m=1}^M{K}_m\exp \left(-{\displaystyle \frac{t}{{\rho }_m}}\right)\end{array}$(2) where G_∞ and K_∞ are the long-term shear modulus and bulk modulus, respectively; G_m and K_m are the shear modulus and bulk modulus in the m^th Maxwell component; ρ_m is the relaxation time in the m^th Maxwell component; M is the total number of Maxwell components.

2.2 Rate-dependent traction-separation behaviour based on CZM

The time-domain viscoelasticity shown in Equations (1(1) $\{\begin{array}{c}S\left(t\right)=S_0\left(t\right)+dev\left[{\int }_0^t\dot{G}\left(t\right){\overline{F}}_t^{-1}\left(t-\tau \right)\cdot e\left(t-\tau \right)\cdot {\overline{F}}_t^T\left(t-\tau \right)d\tau \right]\\ {\sigma }_{vol}\left(t\right)={\sigma }_{vol}^0\left(t\right)+{\int }_0^t\dot{K}\left(\tau \right){ϵ}_{vol}\left(t-\tau \right)d\tau \end{array}$(1) ) and (2(2) $\{\begin{array}{c}G\left(t\right)=G_{\mathrm{\infty }}+{\sum }_{m=1}^M{G}_m\exp \left(-{\displaystyle \frac{t}{{\rho }_m}}\right)\\ K\left(t\right)=K_{\mathrm{\infty }}+{\sum }_{m=1}^M{K}_m\exp \left(-{\displaystyle \frac{t}{{\rho }_m}}\right)\end{array}$(2) ) can also be utilised to characterise the rate-dependent property of cohesive elements with traction-separation viscoelasticity (Zhao et al. 2021). Unlike the finite strain theory adopted for bulk material (e.g. the FAM in this paper), the small strain theory is commonly used in cohesive elements within the FAM and on the interfaces between the FAM and coarse aggregates. The evolution equations for the shear and normal nominal tractions in the two-dimensional case take the form: (3) $\{\begin{array}{c}{t}_s\left(t\right)=t_{s0}\left(t\right)+{\int }_0^t{\dot{E}}_{shear}\left(\tau \right){\delta }_{s0}\left(t-\tau \right)d\tau \\ t_n\left(t\right)=t_{n0}\left(t\right)+{\int }_0^t{\dot{E}}_{normal}\left(\tau \right){\delta }_{n0}\left(t-\tau \right)d\tau \end{array}$(3) where t_s(t) and t_n(t) represent the shear and normal nominal tractions, respectively. t_s0(t) and t_n0(t) are the instantaneous nominal tractions at time t in the local shear and normal directions, respectively. δ_s0(t) and δ_n0(t) stand for the nominal separations in the local shear and normal directions, respectively. E_shear and E_normal are the shear and normal relaxation moduli within the FAM or on the interfaces between the FAM and coarse aggregates.

2.3 Moisture diffusion and coupled moisture-mechanical damage model

The presence of moisture in the asphalt can introduce several types of distress by detrimentally reducing the bonds within the FAM and between the FAM and coarse aggregates. Moisture diffusion in the asphalt is commonly driven by gradients of temperature, pressure, and chemical potential of the moisture. In this paper, only the chemical potential is included in the governing equation to reduce the complexity of moisture diffusion in the asphalt. Equation (4(4) ${\int }_V\text{β}{\displaystyle \frac{d\varphi }{dt}}dV+{\int }_{S}n\cdot \left[-F\cdot \left(\text{β}{\displaystyle \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }x}+\varphi {\displaystyle \frac{\mathrm{\partial }\text{β}}{\mathrm{\partial }x}}}\right)\right]dS=0$(4) ) defines moisture diffusion using the extended Fick’s law (Crank 1979): (4) ${\int }_V\text{β}{\displaystyle \frac{d\varphi }{dt}}dV+{\int }_{S}n\cdot \left[-F\cdot \left(\text{β}{\displaystyle \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }x}+\varphi {\displaystyle \frac{\mathrm{\partial }\text{β}}{\mathrm{\partial }x}}}\right)\right]dS=0$(4) where ϕ is the normalised moisture concentration; β is the solubility of water; V is any volume, whose surface is S; n is the outward unit vector normal to S; F is the moisture diffusivity, and x is the location of volume V.

Moisture damage can typically be sorted into two mechanisms: (1) reduction of adhesion (characterised by damage variable D^a) and (2) loss of cohesion (characterised by damage variable D^c). The former is due to the moisture getting into the interfaces between the FAM and coarse aggregates resulting in adhesive damage. The latter is attributable to spontaneous emulsification of the FAM. The combined effect representing the integration of adhesion damage and cohesion damage was modelled by defining the moisture damage density (i.e. D_mois) as follows (Shakiba et al. 2013): (5) $\left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)=\left(1-D^{\mathrm{a}}\right)\left(1-D^{\mathrm{c}}\right)$(5) Najmeddine and Shakiba (2021) further proposed that the kinetics of D^a (or D^c) can be presented as a function of normalised moisture concentration $\varphi$ and moisture damage density D_mois: (6) ${\displaystyle \frac{d{D}^i}{dt}=k^i\varphi \left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}{)}^q\right(i=\mathrm{a}\mathrm{o}\mathrm{r}\mathrm{c})}$(6) where Dⁱ is adhesive damage density (i = a) or cohesive damage density (i = c); kⁱ is a material parameter characterising the rate of cracking performance (e.g. stiffness, strength, and fracture toughness in the CZM) degradation due to moisture; q is the exponent parameter quantifying the crack history.

The CZM framework mentioned above (e.g. Equation (3(3) $\{\begin{array}{c}{t}_s\left(t\right)=t_{s0}\left(t\right)+{\int }_0^t{\dot{E}}_{shear}\left(\tau \right){\delta }_{s0}\left(t-\tau \right)d\tau \\ t_n\left(t\right)=t_{n0}\left(t\right)+{\int }_0^t{\dot{E}}_{normal}\left(\tau \right){\delta }_{n0}\left(t-\tau \right)d\tau \end{array}$(3) )) will be modified by introducing moisture damage into the dry configuration, which will extend the applicability of the CZM to characterise the coupled moisture-mechanical performance of the asphalt. This paper mainly focuses on the scenario of digital moisture immersion of the asphalt followed by a DITST (digital indirect tensile strength test). Hence, only sequential coupling (one-way coupling) analysis (Castillo and Caro 2014) is adopted to quantify the crack growth within the predefined moisture field of the asphalt. These modifications update the effective tractions in the partially saturated asphalt as follows: (7) $\{\begin{array}{c}{t}_s^e=\left(1-D_{\mathrm{m}\mathrm{e}\mathrm{c}\mathrm{h}}\right)\left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)t_s\\ t_n^e=\{\begin{array}{c}\left(1-D_{\mathrm{m}\mathrm{e}\mathrm{c}\mathrm{h}}\right)\left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)t_n,t_n\ge 0\\ \left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)t_n,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\end{array}$(7) where t e s and t e n are the effective shear and normal nominal tractions in the cohesive elements. D_mech is the mechanical damage introduced by external loading.

3. Generation of two-dimensional mesostructure of the asphalt

Micromechanical models of the asphalt coupled with moisture diffusion are employed to predict the fracture performance of a partially saturated asphalt mixture. To facilitate the construction of these models, asphalt with a nominal maximum aggregate size of 19mm and asphalt binder content of 4.4% was fabricated. Details of the asphalt can be found in the authors’ previous work (Li et al. 2018b). An X-ray computer tomography (CT) scan (Discovery CT 750 HD (high definition) manufactured by General Electronic) is used herein as a non-destructive technique for the generation of the mesostructure of the asphalt. Vertical rectangular slices (i.e. sequences of vertical planar CT images of the selected asphalt) are first obtained following the reconstruction of a three-dimensional asphalt structure (100 mm in diameter and 50 mm in height) by combining these vertical slices spaced with preset intervals (i.e. 0.625mm). The resolution of the scanned rectangular image is 512 × 512 pixels for the asphalt specimen, which makes it clear enough to identify the coarse aggregates (particle size bigger than 2.36mm) embedded in the FAM (i.e. bitumen plus filler plus fine aggregates). To reduce the simulation time and effort, a two-dimensional microstructure (e.g. that digitally obtained from the central slice (in the height direction) of the reconstructed three-dimensional asphalt) is used to characterise the effects of moisture diffusion and content on the fracture performance of the asphalt.

Following CT image capture, MATLAB software is used to transform the original RGB image to grayscale, with subsequent denoising and smoothing processes. Morphological operators (i.e. opening and closing) enhance image quality by adjusting foreground region boundaries. Subsequently, pixel values at each location are calculated and organised into a histogram of pixel intensity values (shown in Figure 1), informing the optimal thresholds for segmenting the FAM, air voids, and coarse aggregates. The Otsu method utilises these thresholds for segmenting components. A threshold of 128 effectively segments air voids. Two potential thresholds (179 and 197) are identified for coarse aggregates and FAM segmentation. Additionally, pixel value 194, corresponding to a minimum count between 197 and 179, is used for comparison. It is found that 194 optimally segments the coarse aggregates and FAM, which aligns with the authors’ experience suggesting that, where two close potential thresholds exist, the pixel value corresponding to the minimum count between them should be the superior choice for image segmentation.

Figure 1. Histogram of a greyscale image of CT scan processed asphalt slice after the morphological operation (e.g. denoising, smoothing, opening, and closing).

In terms of very closely positioned coarse aggregates, another image processing approach (i.e. watershed segmentation) is additionally needed to strengthen the binary image. Specific steps of the watershed segmentation can be summarised as (1) background extraction, (2) foreground extraction, (3) finding the unknown area (neither sure foreground nor for background), and (4) applying the watershed algorithm. After that, vectorisation (using online Image Vectoriser) of these binary images is employed to assemble (using the software AutoCAD) coarse aggregates and air voids to yield a digital sample. Finally, a digital sample of the FAM is produced by Boolean operation (using software ABAQUS) between the above digital components (air voids and coarse aggregates) and the circumferential edge (e.g. a circle with a diameter of 100mm). Figure 2 presents the image processing for the asphalt from the CT scan image to FEM geometrical model.

Figure 2. Conversion from the CT image to FEM geometric model.

The diverse mechanisms by which moisture moves through materials have garnered considerable attention. These mechanisms can be classified into three main categories: water permeability, water capillary rise, and water diffusion. For this study, our focus is on exploring the 2D cracking characteristics of partially saturated asphalt mixtures. Visual representations, as depicted in Figure 2, demonstrate that moisture diffusion into the specimen's interior occurs exclusively through the circumferential boundary. As a result, our investigation excludes considerations of water permeability and water capillary rise as modes of moisture transport.

4. Material parameters and FEM geometry

The DITST at 5°C, chosen in accordance with Nottingham, UK's average annual low temperature on rainy days, is used to evaluate the impact of moisture on the asphalt's fracture performance. This test offers simplicity in terms of geometric modelling and loading conditions, and reproducibility through digitised specimens. Simulations have been carried out to create moisture profiles within the asphalt by soaking it for two, four, and eight weeks.

Mechanical damage to coarse aggregate is assumed to be negligible, and it is regarded as a linear elastic material. The selected values of elastic modulus and Poisson’s ratio are 6 × 10⁴ MPa and 0.15 (Kim and Buttlar 2009, Yin et al. 2012), respectively. In addition, no moisture damage is allowed in the coarse aggregate; hence, the elastic modulus and Poisson’s ratio are moisture independent.

Prony series parameters are used for the characterisation of FAM viscoelasticity. These parameters can be obtained from (1) dynamic mechanical testing of the FAM conducted at different temperatures; (2) construction of the master curve of the complex modulus (i.e. dynamic modulus and phase angle) of FAM at a reference temperature; (3) interconversion between complex modulus and Prony series representation of relaxation modulus. Calibration of Prony series parameters can be found in the authors’ previous work (Li et al. 2018a, Li et al. 2018b). Table 1 presents the Prony series parameters of the dry FAM.

Table 1. Fracture, moisture diffusion and damage, and viscoelastic parameters of the FAM and coarse aggregates at 5°C.

Material Parameter		Aggregate	FAM-aggregate interface	FAM
Fracture strength (MPa)	Normal	NA*	2.1	2.6
	Shear	NA	10.0	10.0
Fracture energy (kJ)	Normal	NA	1.3	1.3
	Shear	NA	13.0	13.0
Diffusion coefficient (m^2/s)		2.4 × 10^−10	2.8 × 10^−11	2.8 × 10^−11
Water solubility β		0.002	0.0013	0.0013
Moisture damage parameter k (s^−1)		NA	4.0 × 10^−6	3.1 × 10^−8
Moisture damage parameter q	NA	5		8
Prony series parameters	NA	τi (s)	Ki (MPa)	Gi (MPa)
		0	10210.2	3657.4
		1 × 10^−5	2603.0	932.4
		1 × 10^−4	2320.6	831.3
		1 × 10^−3	1166.9	418.0
		1 × 10^−2	645.3	231.2
		1 × 10^−1	197.0	70.6
		1 × 10^0	78.3	28.1
		1 × 10^1	30.4	10.9
		1 × 10^2	1.35	0.5
		1 × 10^3	16.9	6.0
		1 × 10^4	29.8	10.7

*NA stands for not applicable.

In this paper, E_normal is assumed to equal the relaxation modulus of fine asphalt mixture, which can be easily deduced from the G_i and K_i in Table 1. This assumption on E_normal is proved to be practical according to the work done by (Zhao et al. 2020a, Chen et al. 2022b). The selected water density in this paper is 1000kg/m³. Water is simplified as soft elastic material with Young’s modulus of 1.26MPa and Poisson’s ratio of 0.4999. The water solubility in the aggregate and FAM shown in Table 1 are referenced from (Kringos 2007) and (Caro et al. 2010a), respectively.

Some well-developed fracture experiments have been carried out to investigate the fracture behaviour of both FAM (i.e. cohesion) and the FAM-aggregate interface (i.e. adhesion). Espinosa et al. (2020) conducted Semi-Circular-Bending tests to evaluate the effect of loading rate on the fracture properties of FAM. Baldi-Sevilla et al. (2017) quantified the influence of bitumen and aggregate polarity on adhesive performance by adding different contents of diatoms and hydrophobic diatoms into the bitumen. This paper selected reliable fracture parameters (shown in Table 1) from the existing works of literature (Apeagyei et al. 2014, Gu 2019) to realistically model the fracture behaviour of the asphalt.

Many well-designed test protocols have been proposed to evaluate moisture diffusion and the resulting damage within FAM and on the coarse aggregate-FAM interface. Using gravimetric techniques, Apeagyei et al. (2014) quantified the moisture concentrations in asphalt mastics and aggregate substrates. Based on this work, they correlated the moisture content and moisture-induced strength degradation of adhesion (using a butt-jointed tensile test) and cohesion (with a dog bone-shaped tensile test). Regarding the moisture damage parameters (kⁱ and q) shown in Equation (6(6) ${\displaystyle \frac{d{D}^i}{dt}=k^i\varphi \left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}{)}^q\right(i=\mathrm{a}\mathrm{o}\mathrm{r}\mathrm{c})}$(6) ), Najmeddine and Shakiba (2021), Shakiba et al. (2015), and Castillo et al. (2017) explicitly presented values and numerically investigated the moisture damage of asphalt in the framework of continuum damage mechanics. This study selected the moisture diffusion and damage parameters (kⁱ and q) from these papers to quantify moisture diffusion and damage to the asphalt. These material parameters are listed in Table 1.

A micromechanical simulation was conducted using the commercial FEM software ABAQUS. It is noted that the distribution of mesh elements potentially affects the distribution of cohesive or adhesive damage. Hence, a mesh sensitivity (e.g. element size and shape) analysis was performed to assure model accuracy. It was found that the mesh size ranges from 10⁻³m to 10⁻⁴m is accurate enough for cracking performance prediction of the selected asphalt. A free meshing technique with mesh control was employed to allow flexibility in element generation (i.e. bulk and interface elements) due to the topological complexity of the coarse aggregates and air voids. The model mesh is presented in Figure 3. The selected boundary condition is identical to the real indirect tensile strength test (ASTM D6931-17). Vertical displacement-controlled load (with constant loading velocity) is applied to the model. The loading strips (12.7 ± 0.3mm) are simplified as rigid bodies to simplify their contact properties with the asphalt. Details of loading and boundary conditions can be found in Figure 3.

Figure 3. Mesh, material components, loading and boundary conditions of finite element model of the asphalt. (a) dry asphalt; (b) partially saturated asphalt; and RP stands for ‘reference point’ used for the definition of the rigid body of the loading strip.

5. Model verification using existing works of literature

The models employed in this paper include moisture diffusion (shown in Equation (4(4) ${\int }_V\text{β}{\displaystyle \frac{d\varphi }{dt}}dV+{\int }_{S}n\cdot \left[-F\cdot \left(\text{β}{\displaystyle \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }x}+\varphi {\displaystyle \frac{\mathrm{\partial }\text{β}}{\mathrm{\partial }x}}}\right)\right]dS=0$(4) )), cohesive zone element-based cracking in both dry (shown in Equations (1(1) $\{\begin{array}{c}S\left(t\right)=S_0\left(t\right)+dev\left[{\int }_0^t\dot{G}\left(t\right){\overline{F}}_t^{-1}\left(t-\tau \right)\cdot e\left(t-\tau \right)\cdot {\overline{F}}_t^T\left(t-\tau \right)d\tau \right]\\ {\sigma }_{vol}\left(t\right)={\sigma }_{vol}^0\left(t\right)+{\int }_0^t\dot{K}\left(\tau \right){ϵ}_{vol}\left(t-\tau \right)d\tau \end{array}$(1) )–(3(3) $\{\begin{array}{c}{t}_s\left(t\right)=t_{s0}\left(t\right)+{\int }_0^t{\dot{E}}_{shear}\left(\tau \right){\delta }_{s0}\left(t-\tau \right)d\tau \\ t_n\left(t\right)=t_{n0}\left(t\right)+{\int }_0^t{\dot{E}}_{normal}\left(\tau \right){\delta }_{n0}\left(t-\tau \right)d\tau \end{array}$(3) )) and partially saturated conditions (shown in Equations (1)–(3) and (5)–(7)). The authors admit that verifications of these models using laboratory-based experiments are very challenging because the inherent difference between the real asphalt specimen and the 2D model exists. Hence, the practical approach for these verifications adopted in this paper is comparing the results from the existing works of literature with that delivered by the models employed in this paper. It is stressed that the authors referenced the material parameters and mesoscale structures from the literature but predicted the performance of asphalt specimens using the equations mentioned above. Result comparisons can be found in Figure 4.

Figure 4. Comparisons between the literature results and those using the models in this paper. (a). Moisture diffusion model verification using literature 1 (Bozorgzad et al. 2018), literature 2 (Arambula et al. 2010a), and literature 3 (Caro et al. 2010b). (b). Cohesive zone model (dry condition) verification using literature 4 (Yin et al. 2012), literature 5 (Yin et al. 2015), and literature 6 (Wang et al. 2014). (c). Cohesive zone model-based moisture damage verification using literature 7 (Caro et al. 2010a), literature 3, and literature 1.

It is found that the moisture diffusion model, cohesive zone model, and cohesive zone model-based moisture damage model employed in this paper are reliable. Regarding the moisture diffusion, Figure 4(a) shows that the result calculated using the model in this paper matches very well with that in works of literature 1, 2 and 3. In Figure 4(b), the result predicted with this paper’s model is a little different from that in literature 4, but these two curves have the same trend with R² of 0.989, and RMSE of 18.652N; hence, its reliability is still acceptable. By comparing the measurement in literature 5 with the prediction in literature 5 and that delivered by this paper’s model, the latter seems to be in line with the experiment result better. Damage density in literature 6 agrees well with that predicted using the model in this paper. Regarding the sequential coupling of the cohesive zone model and moisture damage model, the result in literature 7 matches well with the prediction using this paper’s model, which proves its reliability. Comparisons between the modelling results in this paper and that, respectively, in works of literature 3 and 1 convey that the predictions using this paper’s model are still convincible due to the close curve trend, high R² and low RMSE.

Hence, the models developed in this paper are reliable and can be used to predict the moisture profile and fracture performance of dry or partially saturated asphalt specimens.

6. Results and discussion

6.1 Moisture diffusion

Figure 5 presents the geometric and mesh model, and boundary condition for the digital immersion test (Figure 5(a)), distribution of the moisture concentration at two weeks (Figure 5(b)), five selected points for the characterisation of moisture diffusion (Figure 5(c)), and moisture evolution at these points during the digital immersion test (Figure 5(d)).

Figure 5. Moisture diffusion into the asphalt in digital immersion test. (a), boundary conditions of digital immersion test with a unit normalised concentration on the circumferential edge. (b), moisture distribution after two weeks of diffusion (NT 11 stands for the normalised moisture concentration). (c), selected points in the FAM and coarse aggregate for the characterisation of moisture diffusion. (d), immersion time versus normalised concentration curves at the selected locations.

In the digital immersion test, the normalised moisture concentration (NMC) on the circumference (interface between the controlled humidity 100% and the asphalt) is 1 and the initial NMC in the asphalt is 0, which implies that the moisture will diffuse toward the centre of the asphalt. Due to the difference in the diffusion coefficient between the FAM and the coarse aggregate (shown in Table 1), the NMC is non-uniformly distributed within the partially saturated asphalt. Figure 5(b) demonstrates the contours of NMC after a two-week immersion duration. It clearly shows that the asphalt near the boundary edge is almost saturated, and the NMC near the centre part of the asphalt is substantially lower than that near the boundary edge. The lowest NMC (i.e. 0.3303) is at the centre part of the asphalt.

Figure 5(d) shows that the evolution of NMC can be divided into three stages (i.e. a near-zero growth stage, a rapid growth stage, and a plateau stage) for Points A and D; Points B, C, and E lack a clear near-zero growth stage. The near-zero growth stage represents the time before the moisture front arrives; when it arrives, this instigates the rapid growth stage; the final plateau stage represents near-saturation. It is noted that the NMC at Points B and C are not identical although the distances from Points B and C to the circumference are the same. This is attributable to moisture diffusion coefficient differences between FAM and coarse aggregate and the heterogeneity of the asphalt.

6.2 Moisture-induced reduction of crack resistance

Figure 6 shows the moisture damage quantified by the reduction of peak load in the DITST, stiffness (the slope of the initial part of the loading versus line displacement (LLD) curve), or work of fracture (i.e. the area formed by the loading versus LLD curve and abscissa) at four immersion times (i.e. zero, two weeks, four weeks, and eight weeks). All three measures decrease with immersion duration and can be fitted to exponential equations with R² over 0.98.

Figure 6. Moisture damage growth characterised by peak load, stiffness, and work of fracture in the digital indirect tensile strength test.

After an eight-week immersion duration, the peak load, stiffness, and work of fracture reduce to 63.3%, 74.8%, and 55.1% of the corresponding undamaged values, which represent a considerable amount of moisture damage. Reductions of peak load, stiffness, and work of fracture decline with immersion duration, and most of the moisture damage (81.5% for peak load, 82.4% for stiffness and 90.5% for work of fracture) happens during the soaking of the first four weeks. The main reasons for the above observations are as follows: (1) the increased moisture concentration during the second four weeks is substantially lower than that in the first four weeks (Figure 5(d)); (2) the relationship between normalised moisture diffusion and moisture damage density (characterised by Equations (5(5) $\left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)=\left(1-D^{\mathrm{a}}\right)\left(1-D^{\mathrm{c}}\right)$(5) ) and (6(6) ${\displaystyle \frac{d{D}^i}{dt}=k^i\varphi \left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}{)}^q\right(i=\mathrm{a}\mathrm{o}\mathrm{r}\mathrm{c})}$(6) )) gives very similar curves to that of the normalised moisture concentration versus immersion duration, which indicates that the moisture damage increases much more slowly in the second four weeks than in the first four weeks.

6.3 Anisotropic cracking of dry and partially saturated asphalt

Asphalt is an anisotropic material, which depends on the size, orientation, and sphericity of the coarse aggregate and air voids. The complexity of anisotropic cracking is further aggravated by the inclusion of adhesion (on the interface between the FAM and coarse aggregate) and cohesion (within the FAM). Figure 7 presents simulated cracking test data for dry and partially saturated asphalt (soaked for two weeks) under different directions of load. The anisotropic cracking was simulated by adopting displacement-controlled loading in different directions (shown in Figure 7(a)). Each loading scenario yields a specific simulation result (e.g. loading versus displacement curve) and some critical parameters (peak loading, stiffness, and fracture work) were recorded for the following anisotropic analysis, details of which can be found in Figure 7(b) and Figure 8. Equations (5(5) $\left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}\right)=\left(1-D^{\mathrm{a}}\right)\left(1-D^{\mathrm{c}}\right)$(5) ) and (6(6) ${\displaystyle \frac{d{D}^i}{dt}=k^i\varphi \left(1-D_{\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{s}}{)}^q\right(i=\mathrm{a}\mathrm{o}\mathrm{r}\mathrm{c})}$(6) ) detail how the moisture content introduces moisture damage, which goes ahead to be incorporated into the traditional cohesive zone mode (in dry condition) to calculate the performance reduction of the asphalt mixture exposed for a certain immersion duration.

Figure 7. Anisotropic cracking performance of the dry and partially saturated (two-week immersion) asphalt. (a), loading directions. (b), loading versus LLD curves.

Figure 8. Variations of peak load, stiffness, and work of fracture due to the asphalt anisotropy (in Figure 7).

In Figure 7(b) loading versus LLD curves are presented in directions 0°, 120°, and 180°, and it is noted that, as expected, the difference between the results at 0° and 180° (i.e. in the reverse direction) is negligible in both dry and partially saturated conditions. In addition, the LLD corresponding to the peak loads in the partially saturated and dry scenarios lies in a similar range.

One of the advantages of a digital cracking test over a real one is that the digital cracking test of asphalt mixture can be reused without any limitation, which avoids sample variation in multiple trials of the same experiment. Hence, the measured anisotropy in cracking performance is solely a function of the inherent anisotropy (distribution of coarse aggregate and voids, and consequent distribution of adhesion and cohesion). Figure 8 demonstrates the impact of the loading direction and moisture damage on the cracking performance of the asphalt. Results clearly show that the inclusion of two weeks of moisture damage effectively reduces the peak load, stiffness, and work of fracture except for the work of fracture in the loading direction of 60°. Compared with the failure modes in other loading directions, the one at 60° shows compressive rupture near the bottom plate due to the existence of a couple of air voids and small size aggregate.

Figure 8 reveals significant variations in DITST results due to the asphalt's anisotropy when loaded from different directions. Additionally, variability may be even more pronounced across actual asphalt specimens, despite using identical materials, design methods, and specifications. This emphasises the importance of widely used multiple repeat tests. Figure 9, showing averages, error bars, and the COV (coefficient of variation), highlights the effects of moisture damage with limited impact on variability. The subsequent section explores the specific impact of moisture on the crack growth pattern.

Figure 9. Statistical anisotropic information on peak load, stiffness, and fracture work of the asphalt mixture specimen (in Figure 7).

6.4 Loading-rate-dependent cracking performance of the dry and partially saturated asphalt

Equations (1(1) $\{\begin{array}{c}S\left(t\right)=S_0\left(t\right)+dev\left[{\int }_0^t\dot{G}\left(t\right){\overline{F}}_t^{-1}\left(t-\tau \right)\cdot e\left(t-\tau \right)\cdot {\overline{F}}_t^T\left(t-\tau \right)d\tau \right]\\ {\sigma }_{vol}\left(t\right)={\sigma }_{vol}^0\left(t\right)+{\int }_0^t\dot{K}\left(\tau \right){ϵ}_{vol}\left(t-\tau \right)d\tau \end{array}$(1) ) and (3(3) $\{\begin{array}{c}{t}_s\left(t\right)=t_{s0}\left(t\right)+{\int }_0^t{\dot{E}}_{shear}\left(\tau \right){\delta }_{s0}\left(t-\tau \right)d\tau \\ t_n\left(t\right)=t_{n0}\left(t\right)+{\int }_0^t{\dot{E}}_{normal}\left(\tau \right){\delta }_{n0}\left(t-\tau \right)d\tau \end{array}$(3) ) include the fact that both the bulk and interface properties of the asphalt are rate-dependent, and Figure 10 presents the effect on loading versus LLD curve, cumulative crack length at failure, LLD at crack initiation, and examples of the final cracking pattern (e.g. size and distribution) with 2.5mm LLD. It can be seen in Figure 10(a) that the stiffness, peak load, and work of fracture increase with loading rate and decrease with moisture damage, which matches well with the work done by Espinosa et al. (2020) and Caro et al. (2010a).

Figure 10. Results of digital indirect tensile strength test of the asphalt (dry and two-week immersion) under different loading rates at 5°C. (a), recorded loading versus LLD (LR stands for loading rate). (b), the cumulative crack length of all identified cracks at failure and the LLD at crack initiation. (c), crack patterns in the partially saturated asphalt under different loading rates when the LLD is 2.5mm.

Figure 10(b) demonstrates that the cumulative crack length of both the dry and partially saturated asphalt (e.g. two-week immersion) reduces with the increase in loading rate. The fundamental reason is that an increased loading rate effectively enhances fracture strength and toughness (Caro et al. 2010a, Espinosa et al. 2020). The cumulative crack length (calculated by summing the total length of the cracked cohesive zone element) in dry asphalt (ranging from 278mm to 320mm) is much longer than that in partially saturated asphalt (ranging from 141mm to 227mm). Although moisture damage lowers the fracture strength and toughness of the asphalt leading to reduced anti-cracking capability, it also decreases the stiffness, which makes the stress in the partially saturated asphalt smaller than that in the dry asphalt. In addition, heterogeneity, and moisture in the voids of the asphalt make this effect on the cracking performance more complex. For example, it is found that there are more tiny cracks in the dry asphalt sample than in the partially saturated scenario, which contributes to the cumulative crack length during this digital indirect tensile strength test. The fact that there is more cracking in the dry asphalt than in partially saturated asphalt indicates that the impact of the fracture stiffness of partially saturated asphalt on cracking performance is approximately comparable in significance to its fracture strength at a specified loading rate (i.e. 0.4mm/s, 0.83mm/s, 1mm/s, or 2mm/s).

Figure 10(b) also shows that the LLD at crack initiation increases with the loading rate. The fundamental reason is also that a high loading rate improves the fracture strength and toughness, which effectively delays the initiation of cracking in both dry and partially saturated cases. Only a very small amount of loading displacement (e.g. about 0.03mm) makes the crack initiation both in the dry and partially saturated asphalt sample take place. Compared with dry asphalt, partially saturated asphalt needs a smaller displacement (lower LLD) to get the crack initiated under loading rates of 0.4mm/s, 0.83mm/s, 1mm/s, or 2mm/s.

Figure 10(c) illustrates the pattern of cracking in partially saturated asphalt at a given LLD (2.5mm) under different loading rates. It demonstrates that cracking is resisted better at high loading rates, and the same effect is also seen in dry asphalt. This is well consistent with the observations in Figures 10(a,b).

The fact that the loading rate significantly affects the cracking performance of both dry and partially saturated asphalt can be attributed to the adhesive and cohesive viscosities of the asphalt. Equation (3(3) $\{\begin{array}{c}{t}_s\left(t\right)=t_{s0}\left(t\right)+{\int }_0^t{\dot{E}}_{shear}\left(\tau \right){\delta }_{s0}\left(t-\tau \right)d\tau \\ t_n\left(t\right)=t_{n0}\left(t\right)+{\int }_0^t{\dot{E}}_{normal}\left(\tau \right){\delta }_{n0}\left(t-\tau \right)d\tau \end{array}$(3) ) implies that the traditional time-temperature superposition principle is also applicable to the interface properties (e.g. bond failure) of the asphalt. Hence, the scenario (shown in Figure 10) with a higher loading rate at a given temperature is equivalent to that with a controlled loading rate at a reduced temperature. This is the fundamental reason that the cracking resistance of the asphalt decreases with a reduction in the loading rate.

7. Summary and conclusion

This paper has presented a 2D coupled modelling of moisture diffusion (i.e. Fick’s law) and fracture (i.e. cohesive zone model) at mesoscale to characterise the cracking performance of partially saturated asphalt. The geometrical components (e.g. fine aggregate asphalt mixture, coarse aggregate, and air voids) of the asphalt were accurately obtained by CT scanning and image processing of a cross-sectional slice of asphalt with a nominal maximum aggregate size of 19mm and asphalt binder content of 4.4%. The sequentially coupled modelling of moisture diffusion together with a cohesive zone model was developed to quantify the effects of immersion duration, anisotropy, and loading rate on cracking performance (e.g. peak load, stiffness, and work of fracture) using a digital immersion test and a digital indirect tensile strength test. The main findings and conclusions of this paper are summarised as follows:

1. The normalised moisture concentration is non-uniformly distributed in the partially saturated asphalt, and its evolution curve can be divided into a near-zero stage (only applicable for locations away from the initial moisture-asphalt interface) followed by a rapid growth stage and a final plateau stage.

2. The peak load, stiffness, and work of fracture obtained in a digital indirect tensile strength test all reduce with immersion duration, matched closely by exponential equations. Most of these reductions were completed during the immersion of the first four weeks for the selected asphalt.

3. Varying the loading direction on identical dry or partially saturated asphalt specimens yields different cracking results attributable to its anisotropy. Moisture damage detrimentally impacts the average crack performance in all directions but has a limited impact on its variation.

4. Cracking performance is enhanced by a high loading rate and degraded by moisture damage. The cumulative crack length in dry asphalt is much longer than that in partially saturated asphalt, which implies that the impact of the fracture stiffness of partially saturated asphalt on cracking performance is approximately comparable in significance to its fracture strength at a specified loading rate.

Data Accessibility Statement

The raw data used in this study can be made available based on reasonable requests subject to approval of the data owner.

Acknowledgements

The authors would like to acknowledge the financial support from the Engineering and Physical Sciences Research Council (grant number EP/T019506/1), and the National Natural Science Foundation of China (grant number 51978229 and No. 52278450).

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional information

Funding

This work was supported by Engineering and Physical Sciences Research Council: [Grant Number EP/T019506/1]; National Natural Science Foundation of China: [Grant Number 51978229].