# Experimental investigation and multi-scale Mori-Tanaka modeling of viscoelastic asphalt mastic with imperfect interfaces

K. Lakshmi Roja, George Chatzigeorgiou, Fodil Meraghni, Eyad Masad

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

## To cite this version:

K. Lakshmi Roja, George Chatzigeorgiou, Fodil Meraghni, Eyad Masad. Experimental investigation and multi-scale Mori-Tanaka modeling of viscoelastic asphalt mastic with imperfect interfaces. Construction and Building Materials, 2026, 506, pp.144715. 10.1016/j.conbuildmat.2025.144715. hal-05423461

HAL Id: hal-05423461

https://artsetmetiers.hal.science/hal-05423461v1

Submitted on 18 Dec 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Experimental investigation and multi-scale Mori–Tanaka modeling of viscoelastic asphalt mastic with imperfect interfaces

K. Lakshmi Roja<sup>a ID,*</sup>, G. Chatzigeorgiou<sup>b</sup>, F. Meraghni<sup>b ID</sup>, Eyad Masad<sup>a</sup>

## ARTICLE INFO

**Keywords:**

- Asphalt mastic
- Imperfect interface
- Size effect
- Multi-scale model
- Viscoelasticity

## ABSTRACT

Asphalt mastics, comprising asphalt binder and mineral filler, play a critical role in determining the performance of asphalt mixtures. This study examines the viscoelastic behavior of mastics composed of two asphalt binders and three fillers (gabbro, quartz, and hydrated lime) with varying particle sizes. Repeated Creep and Recovery Tests (RCRT) were conducted to evaluate the influence of filler type and concentration on mastic strain response. Then, a modified multi-scale Mori–Tanaka model was successfully employed to predict the viscoelastic properties of the mastics. The model parameters were found to be influenced by the filler type and size, as well as by imperfections at the binder–filler interface. Results show that mastics containing hydrated lime exhibit higher moduli, particularly at high filler content (40% by volume), with significantly lower final strain in RCRT compared to those with gabbro or quartz. This enhanced performance is attributed to the high surface area of hydrated lime and its chemical interactions with the asphalt binder, which promote the formation of a network structure within the mastic. While the proposed model effectively captures the physical interactions in asphalt mastics, it requires further refinement to account for the additional influence of chemical reactions between fillers, such as hydrated lime, and the asphalt binder.

## 1. Introduction

Asphalt mastics consist of asphalt binder and mineral filler. The interactions between these two components influence asphalt mixture properties, and they are critical in determining the performance of asphalt mixtures [1–3]. Numerous studies have shown that the performance of asphalt mastic is primarily governed by three key mechanisms: mechanical reinforcement, physico-chemical interactions, and chemical reactions [4–6]. Mechanical reinforcement involves two key processes: volume filling and particle structurization. Volume filling refers to the ability of filler particles to occupy the voids within the asphalt mixture, thereby improving the overall mechanical properties. At lower filler contents, the asphalt mastic behaves as a dilute suspension, where the filler particles replace small portions of the asphalt binder, thus resulting in a modest increase in stiffness. As the filler content increases, the interactions between the filler particles intensify as they form an interconnected network structure [7]. This transition marks the onset of particle structurization, a process in which the fillers form a rigid network, significantly increasing the mastic stiffness [4,8,9]. The particle size distribution of mineral fillers plays a crucial role in determining the behavior of asphalt mastics. Finer fillers increase the

total surface area available for interaction with the asphalt binder, which can strengthen adhesion but also increase the likelihood of imperfect or weak interfaces when high filler content is used. These size-dependent effects directly influence the viscoelastic properties and long-term durability of asphalt mastics. Therefore, accounting for the filler size effect in micromechanical modeling is essential to capture the behavior of asphalt mastics [10].

Physico-chemical interactions in asphalt mastic refer to weak interfacial forces, such as Van der Waals, electrostatic, and hydrogen bonding, between the binder and filler, which enhance adhesion without the formation of new chemical compounds. On the other hand, chemical reactions involve the formation of new chemical bonds and compounds at the interface. The fillers investigated in this study (gabbro, quartz, and hydrated lime) are commonly used in asphalt mastics. The mechanical interactions of these fillers with the asphalt binder are influenced by their physical characteristics, such as particle size and surface texture [4,11,12]. Among these fillers, hydrated lime exhibits more pronounced physico-chemical interactions and chemical reactions with the binder, in contrast to the predominantly physical interactions observed with gabbro and quartz. When hydrated lime is blended

<footnote><sup>a</sup> College of Science and Engineering, Hamad Bin Khalifa University, Doha, Qatar</footnote>

<footnote><sup>b</sup> Arts et Métiers Institute of Technology, CNRS, LEM3-UMR 7239 F-57070 Metz, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: lkakumanu@hbku.edu.qa (K.L. Roja), georges.chatzigeorgiou@ensam.eu (G. Chatzigeorgiou), fodil.meraghni@ensam.eu (F. Meraghni), emasad@hbku.edu.qa (E. Masad).</footnote>
with asphalt binder, the calcium hydroxide (Ca(OH)₂) reacts with the acidic/polar components of asphalt, such as carboxylic acids (2R-COOH), forming calcium carboxylate salts (Ca(RCOO)₂). This reaction is irreversible and removes the polar components from the asphalt matrix through their adsorption onto the surface of the hydrated lime particles [13]. The functional group transformation resulting from the lime-bitumen interaction is evidenced in the Fourier Transform Infrared Spectroscopy (FTIR) spectra by a reduction in the intensity of the carbonyl (C=O) peak at approximately 1700 cm⁻¹, coupled with the enhancement of characteristic peaks corresponding to ionic carboxylate salts in the 1560–1580 cm⁻¹ region [14].

Since the late 1990s, researchers have developed micromechanical models to predict the viscoelastic behavior of asphalt mastics and mixtures based on the properties and concentrations of asphalt binder and fillers [15–17]. For example, researchers [17,18] have shown that classical micromechanical models, such as the Nielsen model, exhibit a reasonable agreement with experimental measurements of asphalt mastics at low filler volume concentrations (less than 10%), but underestimate stiffness at higher concentrations (greater than 25%). The difference between the model and experimental results was more pronounced when hydrated lime was used, suggesting that the classical micromechanical models cannot capture the physico-chemical and chemical interactions. Consequently, Underwood et al. [19] expanded on these findings by developing a four-phase micromechanical model that accounts for the physico-chemical interactions to predict the stiffening of asphalt mastics across a range of concentrations. This model successfully predicted experimental data up to a filler volume concentration of about 40%. It is postulated that beyond this concentration, an aggregate structural skeleton begins to form, which is not accounted for in the model. Abbas et al. [20] employed the discrete element method (DEM) and various micromechanical models to describe the exponential increase in mastic stiffness with increasing filler concentration. The authors noted that classical micromechanical models underestimated the stiffening effect of mineral fillers. Their work demonstrated that filler content and type significantly affect the mastic properties, with hydrated lime having a more substantial impact on binder stiffness.

This study builds on early theories by Eshelby [21], who showed how elastic solutions for inclusions can be extended to viscoelastic materials using the elastic–viscoelastic correspondence principle. Hashin [22] developed a model that relates the elastic moduli of mastic to its relaxation moduli and creep compliance of the viscoelastic mastic. Then, Christensen [23] established two theoretical bounds for predicting the effective viscoelastic behavior of two-phase composite materials. Classic approaches such as the self-consistent method [24], the Hashin–Shtrikman bounds [22], and the Mori–Tanaka model [25] were later adapted to describe the time- and frequency-dependent behavior of asphalt binders and mastics. While these classic models laid the foundation for micromechanical predictions, each has limitations when compared to the Mori–Tanaka model. The Hashin–Shtrikman bounds [22] provide only upper and lower bounds for the effective properties rather than a single, explicit prediction, which can be less practical for design and analysis. The self-consistent method [24] can handle complex microstructures but is often more mathematically demanding and can overestimate stiffness when the filler phase dominates. Christensen's bounds [23] similarly give only approximate limits, not exact values. Hashin's viscoelastic extension [22] is conceptually useful but does not directly incorporate the effects of filler shape and interaction as effectively as Mori–Tanaka. Among these models, the Mori–Tanaka approach is more straightforward and computationally efficient while still capturing key effects of filler shape, stiffness, and concentration. Unlike purely empirical models, it directly connects the microstructure to the mastic's viscoelastic response and can be extended to account for frequency and temperature effects. These features make the Mori–Tanaka model well suited for predicting the realistic behavior of asphalt mastics [26–29].

A well-established methodology in the literature to model finite-thickness interphases in composite materials involves replacing them with a hypothetical zero-thickness interface model. The latter is characterized by certain traction and/or displacement jumps. The origins of this idea can be traced back to the pivotal works of researchers [30]. This modeling strategy has been extended in the literature [31–34]. Two main types of non-perfect interfaces appear in the literature: the coherent interface and the elastic interface. In the former case, the zero-thickness region between two bulk materials experiences a displacement jump, while the traction remains continuous. For the latter case, the interface displays a traction jump; however, the displacement remains continuous. Coherent interface models frequently appear in various studies, dating back to the seminal works by Barenblatt and Dugdale [35,36] and their subsequent extensions in the literature [37–44]. Elastic-type interfaces originate from the theory of Gurtin and Murdoch [45], according to which the interface is characterized by its thermodynamic potential. Such an assumption results in a traction jump across the interface due to the Young–Laplace equation [46]. Composites with elastic interfaces have been studied extensively in the literature [47–52]. A combination of elastic and cohesive interfaces also appears in theoretical studies, leading to the concept of generalized interfaces [53–57].

## 2. Study objectives and scope

The objective of this study is to characterize and model the viscoelastic responses of asphalt mastics comprising different concentrations of fillers (gabbro, quartz, and hydrated lime) with varying particle sizes, as well as two types of asphalt binders. The paper also examines the interfacial interactions between the fillers and the asphalt binder. The study combines experimental characterization and mean-field modeling to explicitly account for the influence of filler size and imperfect interfaces. The following tasks were conducted:

1. Examine the influence of filler type and concentration on the viscoelastic response of asphalt mastics through repeated creep and recovery tests (RCRT).

2. Utilize the multi-scale Mori–Tanaka model to describe the viscoelastic response of mastics while accounting for binder–filler interface effects. This is achieved by analyzing the elastic and viscous components of the mastic's response in the RCRT.

3. Investigate the relationship between the Mori–Tanaka model parameters and asphalt–filler interactions, focusing on how filler type and particle size distribution affect these parameters and contribute to the overall performance of the mastic.

The manuscript is organized as follows: Section 3 presents the adopted multiscale model. Section 4 is devoted to the description of asphalt mastic. Section 5 presents experimental characterization, modeling simulations, and predictions, along with a relevant discussion comparing both experimental and modeling results. The main conclusions drawn are discussed in the last section. The article also includes an Appendix that describes the computational details of the Composite Cylinders Assemblage.

## 3. Multiscale strategy

### 3.1. Microscale: Equivalent particle

Mean-field based micromechanics methods rely on the Eshelby solution of the elasticity problem considering an inhomogeneity embedded in an infinite medium [58]. For matrix-inhomogeneity type of mastics with imperfect interfaces, the aim is to identify global dilute elastic and inelastic concentration tensors for the system inhomogeneity/interface. The idea of providing dilute concentration tensors is motivated by similar techniques in the literature for coated particles or fibers [59–61].
<figure><img src="image_4.png" /><figcaption>Fig. 1. Spherical inhomogeneity with imperfect interface inside an infinite matrix undergoing inelastic deformation. The zoomed image is schematic illustration of the occurring displacement jump upon loading.</figcaption></figure>

In Fig. 1, the spherical inhomogeneity of radius $r_i$ with elasticity modulus $L_i$ occupies the space $\Omega_i$, is surrounded by the imperfect interface $I$ and is embedded in the infinite matrix with elastic modulus $L_0$. The matrix occupies the space $\Omega_0$ and is subjected to i) uniform inelastic stress $\sigma_0^p$ and ii) linear displacement $u_0 = \varepsilon_0 \cdot x$ at the far field $\partial\Omega^\infty$.

In the above description, it has already been assumed that only the matrix phase presents inelastic response. The uniformity of the inelastic stress tensor is a strong assumption that can lead to stiff responses [62, 63]. Under this hypothesis, the constitutive relation between stresses and strains is written as

$$ \sigma(x) = \begin{cases} L_i : \varepsilon(x), & x \in \Omega_i, \\ L_0 : \varepsilon(x) + \sigma_0^p, & x \in \Omega_0, \end{cases} \qquad (1) $$

for the bulk phases. The kinematics and the equilibrium (ignoring body and inertia forces) read

$$ \varepsilon(x) = \frac{1}{2} [I \cdot \text{grad}u + [\text{grad}u]^T \cdot I], \quad \text{div}\sigma = 0, \qquad (2) $$

where $I$ is the second order identity tensor, grad is the gradient operator and div is the divergence operator.

For the cohesive imperfect interface the constitutive relation takes into account the displacement jump. The latter requires a second order tensor $\hat{\Lambda}$ to be introduced for this law. The constitutive relation is written

$$ t = \sigma \cdot n = \hat{\Lambda} \cdot [u], \quad x \in I. \qquad (3) $$

In the above expressions $[\cdot]$ is the jump operator [54], based on the methodology of Duan [64], proposed to substitute the inhomogeneity/interface system with an equivalent particle, using energy principles on the problem of Fig. 1. Instead of that approach, the goal of this contribution is to identify the dilute strain concentration tensors and dilute stress-strain concentration tensors for the same system, according to the following formulas:

$$ \langle \varepsilon \rangle_{\Omega_i}^{+} = T : \varepsilon_0 + T^p : \sigma_0^p, \quad \langle \sigma \rangle_{\Omega_i}^{+} = H : \varepsilon_0 + H^p : \sigma_0^p, \qquad (4) $$

with

$$ \langle \varepsilon \rangle_{\Omega_i}^{+} = \frac{1}{2|\Omega_i|} \int_I [u^+ \otimes n + n \otimes u^+] \, dS, \\ \langle \sigma \rangle_{\Omega_i}^{+} = \frac{1}{2|\Omega_i|} \int_I [\sigma^+ \cdot n] \otimes x + x \otimes [\sigma^+ \cdot n] \, dS. \qquad (5) $$

In the above equations, $n$ is the unit normal vector of the interface $I$. Combining the two expressions of Eq. (4) yields

$$ \langle \sigma \rangle_{\Omega_i}^{+} = L_{\text{eq}} : \langle \varepsilon \rangle_{\Omega_i}^{+} + \sigma_{\text{eq}}^{p}. $$

$$ \begin{gathered} L_{\text{eq}} = H \cdot T^{-1}, \\ \sigma_{\text{eq}}^{p} = [H^{p} - L_{\text{eq}} : T^{p}] : \sigma_{0}^{p}. \end{gathered} \qquad (6) $$

For isotropic phases and spherical particles the dilute concentration tensors are isotropic. Hill's notation permits to represent these fourth order tensors in the compact form [65]

$$ \begin{gathered} T = (3T_b, 2T_s), \quad H = (3H_b, 2H_s), \\ T^p = (3T_b^p, 2T_s^p), \quad H^p = (3H_b^p, 2H_s^p). \end{gathered} \qquad (7) $$

Consequently,

$$ L_{\text{eq}} = \left( \frac{H_b}{T_b}, \frac{H_s}{T_s} \right), \quad H^p - L_{\text{eq}} : T^p = \left( 3H_b^p - 3T_b^p \frac{H_b}{T_b}, 2H_s^p - 2T_s^p \frac{H_s}{T_s} \right). \qquad (8) $$

Moreover,

$$ K_{\text{eq}} = \frac{H_b}{3T_b}, \quad \mu_{\text{eq}} = \frac{H_s}{2T_s}. \qquad (9) $$

The dilute concentration tensors are computed using an approach similar to the one introduced by Hashin et al. [66] for the Composite Sphere Assemblage methodology. Two types of boundary value problems are solved analytically for the material system of Fig. 1: one for the hydrostatic response and one for the shear response. The elastic and inelastic dilute concentration tensors are computed through explicit formulas. The computational details of this approach are provided in the Appendix. The final formulas are

$$ 3T_b = \frac{3K_0 + 4\mu_0}{3\tilde{K}_i + 4\mu_0}, \quad 3H_b = 3\tilde{K}_i, \quad 3T_b, \qquad (10) $$

$$ 3T_b^p = \frac{1}{3\tilde{K}_i + 4\mu_0}, \quad 3H_b^p = 3\tilde{K}_i, \quad 3T_b^p, \qquad (11) $$

$$ 2T_s = \frac{5\mu_0[3K_0 + 4\mu_0]}{\mu_0[9K_0 + 8\mu_0] + 6\tilde{\mu}_i[K_0 + 2\mu_0]}, \quad 2H_s = 2\tilde{\mu}_i, \quad 2T_s, \qquad (12) $$

$$ 2T_s^p = \frac{3[K_0 + 2\mu_0]}{\mu_0[9K_0 + 8\mu_0] + 6\tilde{\mu}_i[K_0 + 2\mu_0]}, \quad 2H_s^p = 2\tilde{\mu}_i, \quad 2T_s^p, \qquad (13) $$

where

$$ \tilde{K}_i = K_i \frac{1}{1 + 3K_i / [\hat{\lambda}r_i]}, \quad \tilde{\mu}_i = \mu_i \frac{1}{1 + 2\mu_i / [\hat{\lambda}r_i]}. \qquad (14) $$

In the above expressions, $K_0$, $\mu_0$ are the bulk and shear modulus of the matrix, $K_i$, $\mu_i$ are the bulk and shear modulus of the particle, $\hat{\lambda}$ is the interface elastic parameter (assuming that $\hat{\Lambda} = \hat{\lambda}I$) and $r_i$ is the particle's radius.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Mesoscale RVE consisting of a matrix phase and N spherical particles with generalized interfaces.</figcaption></figure>

**Table 1**

User Material (UMAT) Algorithm for the numerical implementation of the multiscale model in ABAQUS.

<table><tbody><tr><td>Step 1</td><td>At time n everything is known</td></tr><tr><td></td><td>At time n + 1 the increment Δε<sup>meso</sup> is known</td></tr><tr><td></td><td>Set stress σ<sub>0</sub><sup>(old)</sup> = σ<sub>0</sub><sup>(n)</sup> and internal variables IV<sub>0</sub><sup>(old)</sup> = IV<sub>0</sub><sup>(n)</sup></td></tr><tr><td>Step 2</td><td>Use matrix UMAT to correct all matrix fields. Compute T, T<sup>p</sup>, H, H<sup>p</sup> from the dilute problem solutions for all particle sizes. Compute A<sub>0</sub>, A<sub>0</sub><sup>p</sup> from Eq (19). Compute Δε<sub>0</sub>. If prediction is ok, go to Step 6. If not, then:</td></tr><tr><td>Step 3</td><td>Set σ<sub>0</sub> = σ<sub>0</sub><sup>(old)</sup>, IV<sub>0</sub> = IV<sub>0</sub><sup>(old)</sup>, Δε<sub>0</sub><sup>(pr)</sup> = Δε<sub>0</sub></td></tr><tr><td>Step 4</td><td>Use matrix UMAT to correct all matrix fields. Compute T, T<sup>p</sup>, H, H<sup>p</sup> from the dilute problem solutions for all particle sizes. Compute A<sub>0</sub>, A<sub>0</sub><sup>p</sup> from Eq (19). Compute L<sup>meso</sup> and B<sup>p</sup> meso from Eq (24).</td></tr><tr><td>Step 5</td><td>Compute Δε<sub>0</sub>. Check residual: err = ||Δε<sub>0</sub> - Δε<sub>0</sub><sup>(pr)</sup>||. If err > tolerance, return to Step 3 Else move to next step.</td></tr><tr><td>Step 6</td><td>Update ε<sub>0</sub>, σ<sub>0</sub>, σ<sub>0</sub><sup>p</sup>. Compute σ<sup>meso</sup> at time n + 1 from Eq (24). Update all internal variables.</td></tr><tr><td>Step 7</td><td>Compute mesoscopic tangent modulus. Proceed to next time increment.</td></tr></tbody></table>

## 3.2. Mesoscale: Concentration tensors

The dilute concentration tensors is the first step towards computing the behavior of the mesoscale RVE.

Let us assume an RVE consisting of a matrix and N spherical particles with different radii and generalized interfaces (Fig. 2). The matrix has elasticity tensor L₀ and is subjected to inelastic stress σ₀^p. The Jₜh particle has radius rᵢ^j, is characterized by equivalent elastic properties Lₐeq and is subjected to equivalent inelastic stress σₐeq^p. The equivalent properties and inelastic stresses are defined according to the procedure discussed in Section 3.1. Each equivalent particle has volume fraction cⱼ and the matrix has volume fraction c₀. Consequently,

$$c_0 + \sum_{j=1}^{N} c_j = 1. \qquad (15)$$

Under these conditions, the constitutive relation between stresses and strains is written as

$$\sigma(x) = \begin{cases} L_{\text{eq}}^{j} : \epsilon(x) + \sigma_{\text{eq}}^{p j}, & x \in \Omega_{i}^{j}, \\ L_{0} : \epsilon(x) + \sigma_{0}^{p}, & x \in \Omega_{0}. \end{cases} \qquad (16)$$

It is noted that the equivalent particle properties depend on the matrix elasticity tensor.

To apply the Mori–Tanaka/TFA method, the following steps are required: The dilute problem at each equivalent particle gives its average strain and stress as

$$\begin{aligned} \epsilon_{\text{eq}}^{j} &= T^{j} : \epsilon_{0} + T^{p j} : \sigma_{0}^{p}, \\ \sigma_{\text{eq}}^{j} &= H^{j} : \epsilon_{0} + H^{p j} : \sigma_{0}^{p}. \end{aligned} \qquad (17)$$

Consequently,

$$\begin{aligned} \epsilon^{\text{meso}} &= c_0 \epsilon_0 + \sum_{j=1}^{N} c_j [T^j : \epsilon_0 + T^{p j} : \sigma_0^p] \Rightarrow \\ \epsilon^{\text{meso}} &= \left[ c_0 \mathcal{I} + \sum_{j=1}^{N} c_j T^j \right] : \epsilon_0 + \sum_{j=1}^{N} c_j T^{p j} : \sigma_0^p \Rightarrow \end{aligned} \qquad (18)$$

$$\begin{aligned} \epsilon_0 &= A_0 : \epsilon^{\text{meso}} + A_0^p : \sigma_0^p, \\ A_0 &= \left[ c_0 \mathcal{I} + \sum_{j=1}^{N} c_j T^j \right]^{-1}, \end{aligned} \qquad (19)$$

$$A_0^p = -A_0 : \sum_{j=1}^{N} c_j T^{p j}. \qquad (19)$$

Also,

$$\begin{aligned} \sigma^{\text{meso}} &= c_0 [L_0 : \epsilon_0 + \sigma_0^p] + \sum_{j=1}^{N} c_j [H^j : \epsilon_0 + H^{p j} : \sigma_0^p] \Rightarrow \\ \sigma^{\text{meso}} &= \left[ c_0 L_0 + \sum_{j=1}^{N} c_j H^j \right] : \epsilon_0 + \left[ c_0 \mathcal{I} + \sum_{j=1}^{N} c_j H^{p j} \right] : \sigma_0^p. \end{aligned} \qquad (20)$$

The latter relation leads to

$$\sigma^{\text{meso}} = L^{\text{meso}} : \epsilon^{\text{meso}} + \sigma^{p} \epsilon^{\text{meso}}, \qquad (21)$$

$$L^{\text{meso}} = \left[ c_0 L_0 + \sum_{j=1}^{N} c_j H^j \right] : A_0, \qquad (22)$$

$$\sigma^p \text{ meso} = B^p \text{ meso} : \sigma_0^p, \qquad (23)$$

$$B^p \text{ meso} = c_0 \mathcal{I} + \sum_{j=1}^{N} c_j H^{p j} - L^{\text{meso}} : \sum_{j=1}^{N} c_j T^{p j}. \qquad (24)$$

This concludes the micromechanics scheme. To complete the model presentation, the constitutive law of the matrix is discussed in the next subsection.

## 3.3. Viscoelastic matrix

The binder material, which in the multiscale model plays the role of the matrix, is considered a viscoelastic medium with N Kelvin–Voigt viscoelastic branches, whose Helmholtz energy potential depends on the total strain ε, the viscoelastic strain of all branches ε<sub>vᵢ</sub> (i = 1, ..., N). The potential takes the form [67]

$$\Psi \equiv \Psi(\epsilon, \epsilon_{v_i}) = \underbrace{\frac{1}{2}\epsilon_e : C : \epsilon_e}_{\psi^{\text{el}}} + \underbrace{\sum_l \frac{1}{2}\epsilon_{v_l} : C_{v_l} : \epsilon_{v_l}}_{\psi^{\text{v}}}, \qquad (25)$$

where

$$\epsilon_e = \epsilon - \sum_l \epsilon_{v_l}. \qquad (26)$$
<table><caption>**Table 2**

Mineral composition of fillers.</caption><thead><tr><th>S.No</th><th>Filler type</th><th>Mineral composition</th><th>weight%</th></tr></thead><tbody><tr><td rowspan="5">1</td><td rowspan="5">Gabbro</td><td>Plagioclase</td><td>26.90</td></tr><tr><td>Pyroxene</td><td>59.11</td></tr><tr><td>Quartz</td><td>3.58</td></tr><tr><td>Zeolite</td><td>9.75</td></tr><tr><td>Clinochlore</td><td>0.70</td></tr><tr><td>2</td><td>Quartz</td><td>Quartz</td><td>100.00</td></tr><tr><td>3</td><td>Hydrated lime</td><td>Portlandite</td><td>92.22</td></tr><tr><td></td><td></td><td>Calcite</td><td>7.78</td></tr></tbody></table>

In addition, C is the elastic fourth order tensor and $C_{v_i}$ is the viscoelastic fourth order tensor of the branch i. These tensors present minor and major symmetries.

The energy potential leads to the following constitutive relations

$$ \frac{\partial \Psi}{\partial \epsilon} = C : \epsilon_e = \sigma, \qquad (27) $$

$$ -\frac{\partial \Psi}{\partial \epsilon_{v_i}} = \sigma - C_{v_i} : \epsilon_{v_i} = \sigma_{v_i}. \qquad (28) $$

For the inelastic mechanisms one needs to identify evolution laws for the internal variables and activation criteria. These are given in the form of the functions

$$ \Phi_{v_i} = V_{v_i}^{-1} : \sigma_{v_i} - \dot{\epsilon}_{v_i} = 0. \qquad (29) $$

In the above expressions $V_{v_i}$ is the viscous fourth order tensor that presents minor and major symmetries. The elastic part of the material is isotropic, with Young's modulus E and the Poisson's ratio v. Thus, the elasticity is expressed in Hill's notation as

$$ C = \left( \frac{E}{1 - 2v}, \frac{E}{1 + v} \right). \qquad (30) $$

This elastic fourth order tensor corresponds to the fourth order tensor $L_0$ in the micromechanics framework. Each viscoelastic branch is isotropic, with the same Poisson's ratio v, elastic modulus $E_{v_i}$ and characteristic time $\tau_{v_i}$. Thus, the tensors $C_{v_i}$ and $V_{v_i}$ is expressed are Hill's notation as

$$ C_{v_i} = \left( \frac{E_{v_i}}{1 - 2v}, \frac{E_{v_i}}{1 + v} \right), \qquad (31) $$

$$ V_{v_i} = \left( \frac{\tau_{v_i} E_{v_i}}{1 - 2v}, \frac{\tau_{v_i} E_{v_i}}{1 + v} \right). \qquad (32) $$

The numerical implementation of this constitutive law follows the procedure of a return mapping algorithm. Details about this technique is provided in the literature [68,69].

It should be noted that in the studied composite, the binder is several orders of magnitude softer than the fillers. When there is large contrast in the material parameters, it is well known that the analytical micromechanics techniques are not very accurate. In addition, experimentally it is observed that as the filler concentration increases, the binder tends to become stiffer. To address these two effects, a new parameter a is introduced in the micromechanics scheme. This parameter is a multiplier (greater than 1) in the binder's elastic, E, and inelastic $E_{v_i}$ properties, whose value depends on the filler concentration. Thus, for the binder's constitutive law, the micromechanics scheme considers the following substitutions:

$$ C \rightarrow aC, \quad C_{v_i} \rightarrow aC_{v_i}. \qquad (33) $$

### 3.4. Numerical implementation

The integration of the above scheme into a Finite Element (FE) code follows the standard procedure of an ABAQUS UMAT subroutine: All variables at previous time step and the mesoscopic strain increment are provided as input. The mesoscopic stress, all the internal variables and the tangent modulus are computed as output. The algorithm for

the numerical implementation of the reinforced binder is presented in Table 1.

For the mesoscopic tangent modulus (necessary only for FE calculations), the computation follows the same steps with those of the $L^{meso}$. The only difference is that for the matrix, the instantaneous tangent modulus, $L_0^{tang}$, is given as input. In case that $L_0^{tang}$ is not isotropic, it can be “isotropized” to follow the CSA requirements. A simple way is the following:

$$ K_0^{\text{tang}} = \frac{1}{9} \sum_{i=1}^{3} \sum_{j=1}^{3} L_{0_{ij}}^{\text{tang}}, \quad \mu_0^{\text{tang}} = \frac{1}{3} \sum_{i=4}^{6} L_{0_{ii}}^{\text{tang}}. \qquad (34) $$

### 4. Materials description

In this study, the mastics were produced using two types of binder: the first binder is 100% Pen 60–70 (denoted as 0RB), while the second binder is a blend of 70% Pen 60–70 binder and 30% reclaimed asphalt pavement (RAP) binder (denoted as 30RB). The Performance Grade (PG) of the virgin binder (0RB) is PG 64S-22, whereas the 30RB has a grade of PG 70H-16. The use of 30% RAP increased material stiffness, resulting in an increase in the high-temperature grade from 64 °C to 76 °C and an increase in the low-temperature grade from –22 °C to –16 °C.

The mineral composition of the fillers (gabbro, quartz, and hydrated lime) is presented in Table 2. The gabbro filler is primarily composed of plagioclase and pyroxene, while the hydrated lime filler contains calcium hydroxide, and the quartz filler consists of silicon dioxide. Photos of these fillers are presented in Fig. 3. Hydrated lime is a chemically reactive filler that undergoes surface reactions with the asphalt binder, affecting its rheological properties and performance [14]. In contrast, gabbro and quartz are non-reactive fillers that interact with the binder primarily through physical mechanisms.

Fig. 4 shows the particle size distribution of gabbro, quartz, and hydrated lime. For the gabbro (G) sample, 90% of the particles were smaller than 107 µm, with an average particle diameter of 48 µm. For quartz (Q), 90% of the particles were smaller than 70 µm, with an average particle diameter of 34 µm. The hydrated lime (HL) filler had 90% of the particles smaller than 15 µm and an average particle diameter of 7 µm.

Two volume concentrations of fillers (10% and 40%) were used to produce asphalt mastics with ORB and 30RB binders. The asphalt mastic composition and designations are shown in Table 3.

### 5. Characterization of material properties and response

#### 5.1. Linear viscoelastic test

To study the linear viscoelastic behavior of asphalt mastics, frequency sweep tests were conducted over a temperature range of 25 °C to 64 °C, with a frequency range of 1 to 20 Hz. The storage (G') and loss (G") moduli of the asphalt mastics, along with those of the base binders (0RB, 30RB), were measured. The results were used to determine the dynamic shear moduli (|G*|) and phase angle (δ). A minimum of three
<figure><img src="image_6.png" /><figcaption>Fig. 3. Photos showing different fillers- 5 grams of each material.</figcaption></figure>

<table><caption>**Table 3**
Material description and composition of asphalt mastic.</caption><thead><tr><th rowspan="2">Sample ID</th><th colspan="2">Binder percentage</th><th rowspan="2">% Filler relative to binder by volume</th></tr><tr><th>Virgin binder (Pen 60–70)</th><th>RAP</th></tr></thead><tbody><tr><td>0RB</td><td>100</td><td>0</td><td>0</td></tr><tr><td>30RB</td><td>70</td><td>30</td><td>0</td></tr><tr><td>0R-0.1G</td><td>100</td><td>0</td><td>10% gabbro</td></tr><tr><td>0R-0.4G</td><td>100</td><td>0</td><td>40% gabbro</td></tr><tr><td>30R-0.1G</td><td>70</td><td>30</td><td>10% gabbro</td></tr><tr><td>30R-0.4G</td><td>70</td><td>30</td><td>40% gabbro</td></tr><tr><td>0R-0.1HL</td><td>100</td><td>0</td><td>10% hydarted Lime</td></tr><tr><td>0R-0.4HL</td><td>100</td><td>0</td><td>40% hydarted Lime</td></tr><tr><td>30R-0.1HL</td><td>70</td><td>30</td><td>10% hydarted Lime</td></tr><tr><td>30R-0.4HL</td><td>70</td><td>30</td><td>40% hydarted Lime</td></tr><tr><td>0R-0.1Q</td><td>100</td><td>0</td><td>10% quartz</td></tr><tr><td>0R-0.4Q</td><td>100</td><td>0</td><td>40% quartz</td></tr><tr><td>30R-0.1Q</td><td>70</td><td>30</td><td>10% quartz</td></tr><tr><td>30R-0.4Q</td><td>70</td><td>30</td><td>40% quartz</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Volume percentage vs particle diameter for three types of fillers.</figcaption></figure>

tests per material were conducted to ensure the repeatability of tests, and their average was considered in the analysis.

As expected, the dynamic moduli |G*| increased with an increase in filler content. An example is Fig. 5(a), which shows asphalt mastic containing different proportions of gabbro filler, along with its control binders (0RB and 30RB). Examining the results in Fig. 5(b) reveals two findings. First, at low filler volume concentrations, there are slight differences among the mastic moduli. This could be because, at low

filler concentrations, the asphalt mastic behaves like a dilute suspension with little interaction with the binder and among the particles [7]. The second observation in Fig. 5(b) is that the increase in hydrated lime content had the most pronounced effect on increasing the dynamic modulus |G*| (0R-0.4HL vs 0R-0.1HL in Fig. 5(b)). For instance, at 1 Hz frequency, for a 10% filler volume concentration, the modulus of mastic with hydrated lime was 10% higher than that of Gabbro mastic, whereas at a 40% filler volume concentration, it was 211% higher than that of Gabbro mastic. This is attributed to the chemical reaction of hydrated lime and its high surface area, resulting in the formation of a network or rigid aggregate structure.

## 5.2. Repeated creep and recovery test (RCRT)

To characterize the viscoelastic response of the binders and the mastic, a RCRT test protocol was utilized. This protocol consists of 50 loading cycles in total. Every ten cycles, pre-defined shear stress levels are applied: 100 Pa, 320 Pa, 960 Pa, 1600 Pa and 3200 Pa. Each step lasts 10 seconds: during the first second, stress is maintained at the predetermined level, and the specimen is unloaded at practically zero stress for the remaining 9 s. The stress-time curve for the RCRT is illustrated in Fig. 6. The strain accumulation of asphalt mastic prepared with different types of fillers was recorded. Furthermore, the interaction between the binder and filler, as well as the influence of filler size and volume proportions on the performance of asphalt mastic, was modeled using the Mori–Tanaka model. A minimum of three test trials were conducted to ensure repeatability, and their average was considered for analysis.
<figure><img src="image_8.png" /><figcaption>(a) An example of dynamic moduli master curves constructed at a reference temperature of 25 °C for mastics containing gabbro filler.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Experimental master curves of dynamic moduli values at a reference temperature of 25 °C for different filler types and concentrations.</figcaption></figure>

<caption>Fig. 5. Dynamic moduli of asphalt binders and mastics.</caption>

### 5.2.1. Asphalt binder viscoelastic response

As shown earlier in Table 3, two different binders were used in this study (0RB and 30RB). The two binders were subjected to the loading conditions as shown in Fig. 6. The experimentally obtained stress-strain curves are shown in Fig. 7. These experimental curves identified the shear viscoelastic properties of the two different binders through a Python script that combines the least squares method with the trust

region reflective algorithm (see Eq. (31), (32), (35)). The obtained parameters are summarized in Table 4.

Due to the isotropic nature of the mastic, the shear and Young moduli for the elastic part and the viscoelastic branches are connected through the usual relation

$$E = \frac{G}{2[1 + \nu]}, \quad E_{\nu_i} = \frac{G_{\nu_i}}{2[1 + \nu]} \quad (35)$$
<figure><img src="image_10.png" /><figcaption>Fig. 6. Illustration of stresses applied in the RCRT protocol.</figcaption></figure>

<table><caption>**Table 4**

Viscoelastic properties of the two binders, ORB and 30RB, at 25 °C (see Equations (31), (32), (35)).</caption><thead><tr><th>Property</th><th>ORB</th><th>30RB</th></tr></thead><tbody><tr><td>G (MPa)</td><td>3.002</td><td>4.575</td></tr><tr><td>G<sub>v1</sub> (MPa)</td><td>1.986</td><td>3.542</td></tr><tr><td>τ<sub>v1</sub> (s)</td><td>0.493</td><td>0.495</td></tr><tr><td>G<sub>v2</sub> (MPa)</td><td>0.4939</td><td>0.3404</td></tr><tr><td>τ<sub>v2</sub> (s)</td><td>1.744</td><td>3.65</td></tr><tr><td>G<sub>v3</sub> (MPa)</td><td>6.11 · 10<sup>-4</sup></td><td>1.43 · 10<sup>-3</sup></td></tr><tr><td>τ<sub>v3</sub> (s)</td><td>1586.6</td><td>2307.4</td></tr><tr><td>G<sub>v4</sub> (MPa)</td><td>1.43 · 10<sup>-5</sup></td><td>3 · 10<sup>-5</sup></td></tr><tr><td>τ<sub>v4</sub> (s)</td><td>38780.4</td><td>136 834</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 7. The calibration results to determine the model viscoelastic properties using the RCRT experimental data for two types of binders, ORB and 30RB, at 25 °C.</figcaption></figure>

The Poisson's ratio is assumed to be unique for all branches of the viscoelastic model. Its value is equal to 0.35 for both binders. The obtained viscoelastic mastic parameters yield a model that accurately captures the experimental curves, as shown in Fig. 7. It is worth noting that the maximum strain (20%) at the ORB slightly exceeds the typical limit for the small deformation theory assumption. However, this assumption will still be used here, bearing in mind that the masticics operate at smaller strains.

Numerous studies have shown that the performance of asphalt mastic is primarily governed by three key mechanisms: mechanical reinforcement, physico-chemical interactions, and chemical reactions. At this point, it is essential to note that using the viscoelastic material properties derived from the frequency sweep tests at small strains, led to an underestimation of the binder behavior under the more complex repeated loading path depicted in Fig. 6. It has been demonstrated in other types of materials, like semi-crystalline polymers, that proper calibration from creep or relaxation tests can better capture the overall viscoelastic response under complex cyclic loading [67,70,71]. For this reason, the viscoelastic properties of both ORB and 30RB binders are calibrated using the RCRT (Fig. 7 and Table 4).

### 5.2.2. Strain response of asphalt mastic

The RCRT protocol outlined in Section 5.2 was applied to several types of masticics. As discussed earlier, these masticics consist of the two binders (ORB and 30RB) and three types of fillers (gabbro, quartz, and hydrated lime). Two filler contents by volume have been examined: 10% and 40% (see Table 3). The obtained strain-time curves are shown in Fig. 8. The asphalt masticics with gabbro filler exhibited the highest strain values, while those with hydrated lime had the least strain values. As discussed earlier, the small concentration of fillers (10%) did not have much effect on the mastic viscoelastic moduli (see Fig. 5(b)). The same applies to the strain results of the RCRT. However, when 40% filler was used, the strain values of the asphalt significantly. This is especially true for the hydrated lime due to its chemical reactions with the asphalt binders.

Fig. 9 shows the strain accumulation (%) at the end of the RCRT tests for all asphalt masticics. In both Figs. 9a and b, the plain binders (ORB and 30RB) exhibited higher strain values, and as the filler concentration increased in the 0% RAP binder, the strain accumulation
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Experimental strain curves of gabbro, quartz and hydrated lime asphalt mastics considering two binder types (0RB and 30RB) and two filler contents (10% and 40%).</figcaption></figure>

decreased. For the mastics prepared with 30RB, the strain values further decreased because of the use of an aged asphalt binder. The DSR was unable to achieve the required strain values to conduct the RCRT tests on the asphalt mastic prepared with a combination of 30% RAP binder and 40% quartz (30R-0.4Q). It was noticed that during the tests, slippage occurred between the parallel plates, leading to a loss of contact between the rigid surface of the asphalt mastic and the plates.

### 5.3. Modeling results

For modeling, the filler elastic properties are assumed to be the same: bulk modulus 60 GPa and shear modulus 30 GPa. While the actual values may vary and differ between gabbro, quartz, and hydrated

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** The strain accumulation of asphalt mastics contained gabbro, quartz and hydrated lime, (a) 10% filler, (b) 40% filler.</figcaption></figure>

lime, from a modeling perspective, such variation is not expected to cause a noticeable impact on the overall mastic behavior. This is due to the large gap between the binder and the filler's elastic properties.

For the simulation of the mastics, the material parameters presented in the previous subsections should be accompanied by two additional coefficients: the interface parameter $\lambda$ which accounts for the binder-filler interaction, and the binder viscoelastic properties multiplier "a". These two coefficients are sensitive to the type of mastic under study and cannot be determined from isolated experiments on the binder or the fillers. Numerical investigation has demonstrated that $\lambda$ depends on the types of the binder and filler, but it is practically insensitive to the filler's content.

The value of $\lambda$ mainly controls how soft or stiff the overall strain-time response is. On the other hand, "a" strongly depends on the filler type and volume content but seems to be independent of the binder type. The main impact of the value of "a" resides in the oscillation height per cycle in the strain-time curve. Fig. 10 demonstrates the micromechanics model simulations and the experimental findings for gabbro mastics with 0RB at 25 °C. The $\lambda$ for both 10% and 40% filler content is equal to $1.5 \cdot 10^{-3}$ MPa/µm. On the other hand, the "a" is equal to 1 for 10% gabbro (Fig. 10a) and 4 for 40% gabbro (Fig. 10b). Similar analysis, but for mastics prepared with 30RB is illustrated in Fig. 11. For 10% filler content (Fig. 11a), the a is equal to 1 (as with Fig. 10a case), while $\lambda$ is calibrated and found equal to $6 \cdot 10^{-3}$
<figure><img src="image_14.png" /><figcaption>(a) 10% gabbro (0R-0.1G). Model parameters: $\hat{\lambda} = 1.5 \cdot 10^{-3}$ MPa/µm, $a=1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) 40% gabbro (0R-0.4G). Model parameters: $\hat{\lambda} = 1.5 \cdot 10^{-3}$ MPa/µm, $a=4$.</figcaption></figure>

<caption>**Fig. 10.** Modeling simulations vs experimental results for strain-time curves for gabbro mastic with 0RB binder at 25 °C.</caption>

MPa/µm. Using the latter value for the interface coefficient and “a=4” from the case in Fig. 10b, the model predicts the strain-time response quite accurately for 30R-0.4G (Fig. 11b).

Fig. 12 shows the model simulations vs the experimental strain-time curves for quartz mastics with 0RB at 25 °C. The $\hat{\lambda}$ for both 10% and 40% filler content is equal to $2 \cdot 10^{-3}$ MPa/µm. On the other hand, the “a” is equal to 2 for 10% quartz (Fig. 12a) and 4 for 40% quartz (Fig. 12b). It is worth noticing that the model simulation for 40% quartz shows a slightly stiffer response than the experimental curve. A similar analysis, but for a 30R-0.1Q, is illustrated in Fig. 13. The calibration of $\hat{\lambda}$ in this study gives a value equal to $1 \cdot 10^{-2}$ MPa/µm, while the value of “a” is considered equal to 2, similar to the one from Fig. 12a case.

Fig. 14 shows the multi-scale simulations and the experimental results for hydrated lime mastics with 0RB at 25 °C. For both 10% and 40% filler content, $\hat{\lambda}$ is equal to $2.5 \cdot 10^{-2}$ MPa/µm. On the other hand, the “a” is equal to 1.5 for 10% filler (Fig. 14a) and 4 for 40% filler (Fig. 14b). Similar analysis, but for 30RB, is illustrated in Fig. 15. For 10% filler content (Fig. 15a), the “a” is equal to 1.5 (as with Fig. 14a case), while $\hat{\lambda}$ is found to be equal to $2.5 \cdot 10^{-1}$ MPa/µm. Using the latter value for the interface coefficient and “a = 12” from the case of Fig. 14b, the model predicts the strain-time response for 30RB and 40% filler content (Fig. 15b). It is observed that the response of hydrated lime mastics is captured less accurately compared to the other two types of mastics.

The model parameters of all mastics are listed in Table 5. As can be seen, the model parameters $K_p$ and $\mu_p$ were kept constant for all the samples. In most cases, the parameter $\hat{\lambda}$ decreased after introducing the RAP binder. The model parameter “a” both the shift in strain response and its magnitude. The value of “a” increased with an increase in filler volume concentration. However, for hydrated lime samples, even with very high “a” values, the model did not satisfactorily predict the experimental data. This is attributed to the chemical reactions between the hydrated lime and binder, especially given the high surface area of this filler due to its size distribution. Therefore, one can conclude that the Mori-Tanaka model successfully captures the physical interaction between gabbro-asphalt and quartz-asphalt, whereas the chemical interaction between hydrated lime and asphalt is not perfectly captured.

The parameters, “a” (Eq. (33)) and $\hat{\lambda}$ (Eqs. (3) and (14)), capture two distinct interactions between the fillers and binder matrix. As shown in Table 5, the $\hat{\lambda}$ value depends on the materials (binder and filler types) pairing at the interface. A higher $\hat{\lambda}$ value indicates stronger adhesive bonding, which is reflected in the modeling framework as a smaller displacement jump at the interface (Eqs. (3) and (14)). As shown in Table 5, hydrated lime has a higher $\hat{\lambda}$ value due to its physico-chemical interactions and strong chemical bonds with the binder (see Fig. 16a). In contrast, gabbro and quartz primarily exhibit mechanical
<figure><img src="image_16.png" /><figcaption>(a) 10% gabbro (30R-0.1G). Model parameters: $\hat{\lambda} = 6 \cdot 10^{-3}$ MPa/μm, $a=1$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) 40% gabbro (30R-0.4G). Model parameters: $\hat{\lambda} = 6 \cdot 10^{-3}$ MPa/μm, $a=4$.</figcaption></figure>

<caption>**Fig. 11.** Modeling simulations vs experimental results for strain-time curves for gabbro mastic with 30RB at 25 °C.</caption>

interactions with the binder (see Fig. 16b). These weak bonds are associated with displacement jumps at the interface, which are reflected in low $\hat{\lambda}$ values in the model. For all fillers, the use of RAP leads to an increase in the $\hat{\lambda}$.

The “a” parameter captures the effect of filler type and concentration on stiffening the binder matrix. As shown in Table 5, the “a” value increases with an increase in filler content. This might be explained by the development of an interphase, which exhibits unique properties depending on the interactions of the two materials. This interphase has higher moduli than the bulk binder (“a” value higher than one) for the three fillers. Interestingly, hydrated lime has a greater effect on stiffening the binder phase, as indicated by a higher “a” parameter value, than gabbro and quartz. This is attributed to the strong and possibly wide interphase formed by the physico-chemical interactions and strong chemical bonds in the case of hydrated lime (see Fig. 16a) and the small interphase formed primarily by mechanical interactions between the binder matrix and quartz or gabbro (see Fig. 16b).

## 6. Conclusions

The study examined the interactions of various mineral fillers with asphalt binders and their effects on the viscoelastic properties of asphalt mastics. This was achieved through the use of repeated creep and recovery tests (RCRT) to test combinations of materials, and then adopting the appropriate Multi-Scale Mori–Tanaka model. The mineral

fillers used for mastic preparation were gabbro, quartz, and hydrated lime. The hydrated lime had significantly smaller particles than the other two fillers. The RCRT strain response data over a wide range of loading cycles were used to determine the model parameters. Then, the model parameters were related to the binder type, filler type, and the interactions between the two materials. The following conclusions can be drawn from the study.

1. At low filler concentrations (10%), the fillers had a minimal impact on increasing the mastic modulus compared to the plain asphalt binder, and the mastic behaved more like a diluted suspension. In contrast, at higher concentrations (40%), a substantial increase in modulus was observed. Specifically, mastics containing 40% hydrated lime exhibited the highest moduli, followed by those with quartz and gabbro. This behavior was attributed to the physico-chemical and chemical interactions between hydrated lime and the asphalt binder. Furthermore, the high surface area of hydrated lime contributed to enhancing these interactions.

2. The multi-scale Mori–Tanaka model successfully captured the RCRT viscoelastic strain response of the mastics. However, the model predictions were closer to the experimental strain values for samples with low filler content than for those with high filler content. The model showed better agreement for gabbro and quartz mastics, as it implicitly accounted for their physical
<figure><img src="image_18.png" /><figcaption>(a) 10% quartz (0R-0.1Q). Model parameters: $\hat{\lambda} = 2 \cdot 10^{-3}$ MPa/µm, $a=2$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) 40% quartz (0R-0.4Q). Model parameters: $\hat{\lambda} = 2 \cdot 10^{-3}$ MPa/µm, $a=4$.</figcaption></figure>

<caption>**Fig. 12.** Modeling simulations vs experimental results for strain-time curves for quartz mastic with ORB binder at 25 °C.</caption>

<table><caption>**Table 5**
Mori–Tanaka model parameters.</caption><thead><tr><th>Material</th><th>K<sub>p</sub> (MPa)</th><th>μ<sub>p</sub> (MPa)</th><th>λ̂ (MPa/μm)</th><th>a (-)</th></tr></thead><tbody><tr><td>0R-0.1G</td><td>60 000</td><td>30 000</td><td>1.50E-03</td><td>1</td></tr><tr><td>0R-0.4G</td><td>60 000</td><td>30 000</td><td>1.50E-03</td><td>4</td></tr><tr><td>30R-0.1G</td><td>60 000</td><td>30 000</td><td>6.00E-03</td><td>1</td></tr><tr><td>30R-0.4G</td><td>60 000</td><td>30 000</td><td>6.00E-03</td><td>4</td></tr><tr><td>0R-0.1HL</td><td>60 000</td><td>30 000</td><td>2.50E-02</td><td>1.5</td></tr><tr><td>0R-0.4HL</td><td>60 000</td><td>30 000</td><td>2.50E-02</td><td>12</td></tr><tr><td>30R-0.1HL</td><td>60 000</td><td>30 000</td><td>2.50E-01</td><td>1.5</td></tr><tr><td>30R-0.4HL</td><td>60 000</td><td>30 000</td><td>2.50E-01</td><td>12</td></tr><tr><td>0R-0.1Q</td><td>60 000</td><td>30 000</td><td>2.00E-03</td><td>2</td></tr><tr><td>0R-0.4Q</td><td>60 000</td><td>30 000</td><td>2.00E-03</td><td>4</td></tr><tr><td>30R-0.1Q</td><td>60 000</td><td>30 000</td><td>1.00E-02</td><td>2</td></tr></tbody></table>

interactions with the binder. In contrast, the model was less effective in capturing the chemical interactions associated with hydrated lime, particularly at higher filler concentrations.

3. The model parameters, “a” and $\lambda$ capture two distinct phenomena in the mastic. The $\lambda$ characterizes the quality of the binder-filler interface. It reflects the degree of adhesion between the two materials (binder and filler) by accounting for displacement jumps at the interface: weaker bonding allows larger relative displacements, while stronger bonding restricts them. Thus, a higher $\lambda$ value corresponds to stronger adhesive interactions and manifests in the model as a smaller displacement discontinuity across the interface. The hydrated lime has

a higher $\lambda$ value due to its physico-chemical interactions and strong chemical bonds with the binder. On the other hand, gabbro and quartz primarily exhibit mechanical interactions with the binder, which are associated with larger displacement jumps at the interface. The “a” parameter captures the effect of filler type and concentration on stiffening the binder matrix. The “a” value increases with an increase in filler content. Hydrated lime has a greater effect on stiffening the binder phase, as indicated by a higher “a” parameter value, than gabbro and quartz.

Future research will expand the testing matrix to include various filler types, particle size distributions, asphalt binder sources, and
<figure><img src="image_20.png" /><figcaption>**Fig. 13.** Modeling simulations vs experimental results for strain-time curves for 10% quartz mastic with 30RB binder (30R-0.1Q) at 25 °C. Model parameters: $\hat{\lambda} = 1 \cdot 10^{-2}$ MPa/µm, $a=2$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(a) 10% hydrated lime (0R-0.1HL). Model parameters: $\hat{\lambda} = 2.5 \cdot 10^{-2}$ MPa/µm, $a=1.5$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) 40% hydrated lime (0R-0.4HL). Model parameters: $\hat{\lambda} = 2.5 \cdot 10^{-2}$ MPa/µm, $a=12$.</figcaption></figure>

<caption>**Fig. 14.** Modeling simulations vs experimental results for strain-time curves for hydrated lime mastic with ORB binder at 25 °C.</caption>

asphalt mixture-level modulus evaluations. Such efforts will advance the understanding of filler effects on mastic behavior and overall asphalt mixture performance. Additionally, broadening experimental measurements will help link model parameters to intrinsic material properties and interfacial interactions.

<figure><img src="image_23.png" /><figcaption>(a) 10% hydrated lime (30R-0.1HL). Model parameters: $\hat{\lambda} = 2.5 \cdot 10^{-1}$ MPa/µm, $a=1.5$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) 40% hydrated lime (30R-0.4HL). Model parameters: $\hat{\lambda} = 2.5 \cdot 10^{-1}$ MPa/µm, $a=12$.</figcaption></figure>

<caption>**Fig. 15.** Modeling simulations vs experimental results for strain-time curves for hydrated lime mastic with 30RB binder at 25 °C.</caption>

CRediT authorship contribution statement

**K. Lakshmi Roja:** Writing – original draft, Resources, Methodology, Investigation, Conceptualization. **G. Chatzigeorgiou:** Writing – original draft, Validation, Formal analysis, Conceptualization. **F. Meraghni:** Writing – review & editing, Validation, Methodology, Formal analysis, Conceptualization. **Eyad Masad:** Writing – review & editing, Supervision, Resources, Methodology, Investigation, Conceptualization.

Declaration of competing interest

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.

Acknowledgments

The research reported in this paper was supported by the Qatar Research Development and Innovation Council [ARG01-0429-230026]. The content is solely the responsibility of the authors and does not necessarily represent the official views of Qatar Research Development and Innovation Council.
<figure><img src="image_25.png" /><figcaption>Fig. 16. Schematic illustrating the interphase formed by (a) the physico-chemical interaction and chemical reaction between hydrated lime and binder, and (b) the mechanical interaction between gabbro or quartz with binder.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. A.17. Spherical coordinates.</figcaption></figure>

Appendix. Composite sphere assemblage strategy for equivalent particle

A.1. Problem definition in spherical coordinates

We consider the axes $r, \theta, \phi$ (Fig. A.17) to be represented with the numbers 1, 2, 3 respectively.

For inelastic isotropic constituents we have in Voigt notation [69]

$$ \begin{bmatrix} \sigma_{rr} \\ \sigma_{\theta\theta} \\ \sigma_{\phi\phi} \\ \sigma_{r\theta} \\ \sigma_{r\phi} \\ \sigma_{\theta\phi} \end{bmatrix} = \begin{bmatrix} K + \frac{4\mu}{3} & K - \frac{2\mu}{3} & K - \frac{2\mu}{3} & 0 & 0 & 0 \\ K - \frac{2\mu}{3} & K + \frac{4\mu}{3} & K - \frac{2\mu}{3} & 0 & 0 & 0 \\ K - \frac{2\mu}{3} & K - \frac{2\mu}{3} & K + \frac{4\mu}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix} \begin{bmatrix} \varepsilon_{rr} \\ \varepsilon_{\theta\theta} \\ \varepsilon_{\phi\phi} \\ 2\varepsilon_{r\theta} \\ 2\varepsilon_{r\phi} \\ 2\varepsilon_{\theta\phi} \end{bmatrix} + \begin{bmatrix} \sigma_{rr}^p \\ \sigma_{\theta\theta}^p \\ \sigma_{\phi\phi}^p \\ \sigma_{r\theta}^p \\ \sigma_{r\phi}^p \\ \sigma_{\theta\phi}^p \end{bmatrix} \quad (A.1) $$

It is noted that

$$ \begin{align} \varepsilon_{rr} &= \frac{\partial u_r}{\partial r}, \quad \varepsilon_{\theta\theta} = \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} + \frac{u_r}{r}, \quad \varepsilon_{\phi\phi} = \frac{1}{r} \frac{\partial u_\phi}{\partial \phi} + \frac{u_r}{r} + \frac{u_\theta \cos \theta}{r \sin \theta}, \nonumber \\ 2\varepsilon_{r\phi} &= \frac{\partial u_\phi}{\partial r} + \frac{1}{r} \frac{\partial u_r}{\partial \phi} - \frac{u_\phi}{r}, \quad 2\varepsilon_{\theta\phi} = \frac{1}{r} \frac{\partial u_\phi}{\partial \theta} + \frac{1}{r} \frac{\partial u_\theta}{\partial \phi} - \frac{u_\phi \cos \theta}{r \sin \theta}, \nonumber \\ 2\varepsilon_{r\theta} &= \frac{\partial u_\theta}{\partial r} + \frac{1}{r} \frac{\partial u_r}{\partial \theta} - \frac{u_\theta}{r}. \tag{A.2} \end{align} $$

The equilibrium equations are written as

$$ \begin{align*} \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\sigma_{r\theta} \cos \theta}{r \sin \theta} + \frac{2\sigma_{rr} - \sigma_{\theta\theta} - \sigma_{\phi\phi}}{r} + \frac{1}{r \sin \theta} \frac{\partial \sigma_{r\phi}}{\partial \phi} &= 0, \\ \frac{\partial \sigma_{r\theta}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta\theta}}{\partial \theta} + \frac{3\sigma_{r\theta}}{r} + \frac{[\sigma_{\theta\theta} - \sigma_{\phi\phi}] \cos \theta}{r \sin \theta} + \frac{1}{r \sin \theta} \frac{\partial \sigma_{\theta\phi}}{\partial \phi} &= 0, \end{align*} $$

$$ \frac{\partial \sigma_{r\phi}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta\phi}}{\partial \theta} + \frac{3\sigma_{r\phi}}{r} + \frac{2\sigma_{\theta\phi} \cos \theta}{r \sin \theta} + \frac{1}{r \sin \theta} \frac{\partial \sigma_{\phi\phi}}{\partial \phi} = 0. \quad (A.3) $$

The particle is considered to have radius $r = r_i$. The interface is assumed to be isotropic, i.e.

$$ \hat{\Lambda} = \hat{\lambda} I, \quad t = \hat{\lambda} \|u\|. \quad (A.4) $$

The necessary vectors and tensors for the computations are represented in tensorial form as

$$ \begin{align} \boldsymbol{u} &= u_r \boldsymbol{n}_r + u_\theta \boldsymbol{n}_\theta + u_\phi \boldsymbol{n}_\phi, \nonumber \\ \boldsymbol{\sigma} \cdot \boldsymbol{n}_r &= \sigma_{rr} \boldsymbol{n}_r + \sigma_{r\theta} \boldsymbol{n}_\theta + \sigma_{r\phi} \boldsymbol{n}_\phi, \nonumber \\ \tilde{\boldsymbol{\epsilon}} &= \frac{1}{2} [\boldsymbol{u} \otimes \boldsymbol{n}_r + \boldsymbol{n}_r \otimes \boldsymbol{u}], \nonumber \\ \tilde{\boldsymbol{\sigma}} &= \frac{1}{2} [[\boldsymbol{\sigma} \cdot \boldsymbol{n}_r] \otimes \boldsymbol{x} + \boldsymbol{x} \otimes [\boldsymbol{\sigma} \cdot \boldsymbol{n}_r]], \tag{A.5} \end{align} $$

where

$$ \boldsymbol{n}_r = \begin{bmatrix} \sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \end{bmatrix}, \quad \boldsymbol{n}_\theta = \begin{bmatrix} \cos \theta \cos \phi \\ \cos \theta \sin \phi \\ -\sin \theta \end{bmatrix}, \quad \boldsymbol{n}_\phi = \begin{bmatrix} -\sin \phi \\ \cos \phi \\ 0 \end{bmatrix}, \quad (A.6) $$

$$ \boldsymbol{x} = \begin{bmatrix} r \sin \theta \cos \phi \\ r \sin \theta \sin \phi \\ r \cos \theta \end{bmatrix} = r \boldsymbol{n}_r. \quad (A.7) $$

In Voigt notation and spherical coordinates one has

$$ \tilde{\boldsymbol{\epsilon}}^{\text{sph}} = \begin{bmatrix} u_r \\ 0 \\ 0 \\ u_\theta \\ u_\phi \\ 0 \end{bmatrix}, \quad \tilde{\boldsymbol{\sigma}}^{\text{sph}} = r \begin{bmatrix} \sigma_{rr} \\ 0 \\ 0 \\ \sigma_{r\theta}/2 \\ \sigma_{r\phi}/2 \\ 0 \end{bmatrix}. \quad (A.8) $$

Passing from spherical (local) to Cartesian (global) coordinates is a procedure that requires the introduction of the rotator tensor

$$ \mathbf{R} = \begin{bmatrix} R_{11} & R_{12} & R_{13} \\ R_{21} & R_{22} & R_{23} \\ R_{31} & R_{32} & R_{33} \end{bmatrix} = \begin{bmatrix} \sin \theta \cos \phi & \cos \theta \cos \phi & -\sin \phi \\ \sin \theta \sin \phi & \cos \theta \sin \phi & \cos \phi \\ \cos \theta & -\sin \theta & 0 \end{bmatrix}. \quad (A.9) $$

When dealing with Voigt notation, the above rotator should be transformed to the following 6 × 6 matrices (see Box I).

Then the transformation of stress and strain fields between global ($\sigma, \epsilon$) and local ($\sigma', \epsilon'$) is obtained through the relations

$$ \varepsilon' = Q_{\varepsilon} \cdot \varepsilon, \quad \sigma' = Q_{\sigma} \cdot \sigma, \quad \varepsilon = Q_{\varepsilon}^{T} \cdot \varepsilon', \quad \sigma = Q_{\sigma}^{T} \cdot \sigma'. \quad (A.11) $$

The particle has outer radius $r_i$. The average total fields of the equivalent particle are given by the expressions

$$ \begin{align} \langle \varepsilon \rangle_{\Omega_i}^{+} &= \frac{3}{4\pi r_i} \int_0^{2\pi} \int_0^{\pi} \mathbf{Q}_{\sigma}^{T} \cdot \tilde{\boldsymbol{\epsilon}}^{\text{sph+}}(r_i, \theta, \phi) \sin \theta d\theta d\phi, \nonumber \\ \langle \sigma \rangle_{\Omega_i}^{+} &= \frac{3}{4\pi r_i} \int_0^{2\pi} \int_0^{\pi} \mathbf{Q}_{\varepsilon}^{T} \cdot \tilde{\boldsymbol{\sigma}}^{\text{sph+}}(r_i, \theta, \phi) \sin \theta d\theta d\phi. \tag{A.12} \end{align} $$
$$
\begin{align}
Q_{\epsilon} &= \begin{bmatrix}
R_{11}^2 & R_{21}^2 & R_{31}^2 & R_{11} R_{21} & R_{11} R_{31} & R_{21} R_{31} \\
R_{12}^2 & R_{22}^2 & R_{32}^2 & R_{12} R_{22} & R_{12} R_{32} & R_{22} R_{32} \\
R_{13}^2 & R_{23}^2 & R_{33}^2 & R_{13} R_{23} & R_{13} R_{33} & R_{23} R_{33}
\end{bmatrix}, \nonumber \\
Q_{\sigma} &= \begin{bmatrix}
2R_{11} R_{12} & 2R_{21} R_{22} & 2R_{31} R_{32} & R_{12} R_{21} + R_{11} R_{22} & R_{12} R_{31} + R_{11} R_{32} & R_{22} R_{31} + R_{21} R_{32} \\
2R_{11} R_{13} & 2R_{21} R_{23} & 2R_{31} R_{33} & R_{13} R_{21} + R_{11} R_{23} & R_{13} R_{31} + R_{11} R_{33} & R_{23} R_{31} + R_{21} R_{33} \\
2R_{12} R_{13} & 2R_{22} R_{23} & 2R_{32} R_{33} & R_{13} R_{22} + R_{12} R_{23} & R_{13} R_{32} + R_{12} R_{33} & R_{23} R_{32} + R_{22} R_{33}
\end{bmatrix} \nonumber \\
&\qquad \begin{bmatrix}
R_{11}^2 & R_{21}^2 & R_{31}^2 & 2R_{11} R_{21} & 2R_{11} R_{31} & 2R_{21} R_{31} \\
R_{12}^2 & R_{22}^2 & R_{32}^2 & 2R_{12} R_{22} & 2R_{12} R_{32} & 2R_{22} R_{32} \\
R_{13}^2 & R_{23}^2 & R_{33}^2 & 2R_{13} R_{23} & 2R_{13} R_{33} & 2R_{23} R_{33}
\end{bmatrix} \nonumber \\
&\qquad \begin{bmatrix}
R_{11} R_{12} & R_{21} R_{22} & R_{31} R_{32} & R_{12} R_{21} + R_{11} R_{22} & R_{12} R_{31} + R_{11} R_{32} & R_{22} R_{31} + R_{21} R_{32} \\
R_{11} R_{13} & R_{21} R_{23} & R_{31} R_{33} & R_{13} R_{21} + R_{11} R_{23} & R_{13} R_{31} + R_{11} R_{33} & R_{23} R_{31} + R_{21} R_{33} \\
R_{12} R_{13} & R_{22} R_{23} & R_{32} R_{33} & R_{13} R_{22} + R_{12} R_{23} & R_{13} R_{32} + R_{12} R_{33} & R_{23} R_{32} + R_{22} R_{33}
\end{bmatrix}. \tag{A.10}
\end{align}
$$

Box I.

A.2. Dilute concentration terms for bulk modulus

A.2.1. Elastic

For these properties we apply the following displacement-strain field at far distance:

$$
\mathbf{u}_{0}^{(x,y,z)} = \begin{bmatrix} x \\ y \\ z \end{bmatrix}, \quad \varepsilon_{0} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \qquad (\text{A.13})
$$

In spherical coordinates, the displacement is transformed to the vector

$$
\mathbf{u}_{0}^{(r,\theta,\phi)} = \begin{bmatrix} r \\ 0 \\ 0 \end{bmatrix}. \qquad (\text{A.14})
$$

For these boundary conditions, the important displacements and the stresses for the calculations at each phase are given by the general expressions

$$
u_r^{(q)} = \Xi_1^{(q)} r + \Xi_2^{(q)} \frac{r_i^3}{r^3}, \quad q = 0, i, \qquad (\text{A.15})
$$

$$
\sigma_{rr}^{(q)} = 3K_q \Xi_1^{(q)} - 4\mu_q \Xi_2^{(q)} \frac{r_i^3}{r^3}, \quad q = 0, i, \qquad (\text{A.16})
$$

The unknowns that need to be defined are the $\Xi_1^{(i)}$, $\Xi_2^{(i)}$, $\Xi_1^{(0)}$, $\Xi_2^{(0)}$. The boundary and interface conditions that we need to satisfy are

$$
u_r^{(i)} \text{finite at } r = 0 \rightarrow \Xi_2^{(i)} = 0, \qquad (\text{A.17})
$$

$$
u_r^{(0)}(r \to \infty) = r \to \Xi_1^{(0)} = 1, \qquad (\text{A.18})
$$

$$
t_r = \hat{\lambda} \llbracket u_r \rrbracket \rightarrow \sigma_{rr}^{(i)}(r_i) = \hat{\lambda} [u_r^{(0)}(r_i) - u_r^{(i)}(r_i)], \qquad (\text{A.19})
$$

$$
\llbracket \sigma_{rr}(r_i) \rrbracket = 0 \rightarrow \sigma_{rr}^{(0)}(r_i) = \sigma_{rr}^{(i)}(r_i). \qquad (\text{A.20})
$$

This linear system takes the form

$$
\mathbf{K}_b \cdot \mathbf{\Xi}_b = \mathbf{F}_b, \qquad (\text{A.21})
$$

where

$$
\mathbf{K}_b = \begin{bmatrix} 1 + \frac{3K_i}{\hat{\lambda}r_i} & -1 \\ 3K_i & 4\mu_0 \end{bmatrix}, \qquad (\text{A.22})
$$

$$
\Xi_b = \begin{bmatrix} \Xi_1^{(i)} \\ \Xi_2^{(0)} \end{bmatrix}, \qquad (A.23)
$$

$$
F_b = \begin{bmatrix} 1 \\ 3K_0 \end{bmatrix}. \qquad (A.24)
$$

After solving this linear system, the average stain and stress in the particle/interface system are equal to

$$
\langle \epsilon \rangle_{\Omega_i}^{+} = \frac{1}{r_i} u_r^{(0)}(r_i) \epsilon_0, \quad \langle \sigma \rangle_{\Omega_i}^{+} = \sigma_{rr}^{(0)}(r_i) \epsilon_0. \qquad (\text{A.25})
$$

Consequently, the bulk dilute concentration terms are equal to

$$
3T_b = 1 + \Xi_2^{(0)}, \quad 3H_b = 3K_0 - 4\mu_0\Xi_2^{(0)}, \qquad (A.26)
$$

or

$$
3T_b = \frac{3K_0 + 4\mu_0}{3\tilde{K}_i + 4\mu_0}, \quad 3H_b = 3\tilde{K}_i 3T_b, \quad \tilde{K}_i = K_i \frac{1}{1 + 3K_i / [\hat{\lambda}r_i]}. \qquad (\text{A.27})
$$

A.2.2. Inelastic

For these properties we apply zero displacement field at far distance.
We also consider the following inelastic stress tensor

$$
\sigma_0^p = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \qquad (\text{A.28})
$$

The displacement field is described by the same general forms with the elastic case. The matrix phase stress is expressed as

$$
\sigma_{rr}^{(0)} = 3K_0 \Xi_1^{(0)} - 4\mu_q \Xi_2^{(0)} \frac{r_i^3}{r^3} + 1. \qquad (\text{A.29})
$$
The unknowns that need to be defined are the $\Xi_1^{(i)}$, $\Xi_2^{(i)}$, $\Xi_1^{(0)}$, $\Xi_2^{(0)}$. The boundary and interface conditions that we need to satisfy are

$$u_r^{(i)} \text{ finite at } r = 0 \to \Xi_2^{(i)} = 0, \qquad (A.30)$$

$$u_r^{(0)}(r \to \infty) = 0 \to \Xi_1^{(0)} = 0, \qquad (A.31)$$

$$t_r = \lambda \llbracket u_r \rrbracket \to \sigma_{rr}^{(i)}(r_i) = \lambda [u_r^{(0)}(r_i) - u_r^{(i)}(r_i)], \qquad (A.32)$$

$$\llbracket u_{rr}(r_i) \rrbracket = 0 \to \sigma_{rr}^{(0)}(r_i) = \sigma_{rr}^{(i)}(r_i). \qquad (A.33)$$

This linear system takes the form

$$K_b \cdot \Xi_b = F_b^p, \qquad (A.34)$$

where $K_b$ and $\Xi_b$ are the same with those of the elastic case and

$$F_b^p = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \qquad (A.35)$$

After solving this linear system, the average stain and stress in the particle/interface system are equal to

$$\langle \epsilon \rangle_{\Omega_i}^{+} = \frac{1}{r_i} u_r^{(0)}(r_i) \sigma_0^p, \quad \langle \sigma \rangle_{\Omega_i}^{+} = \sigma_{rr}^{(0)}(r_i) \sigma_0^p. \qquad (A.36)$$

Consequently, the bulk inelastic dilute concentration terms are equal to

$$3T_b^p = \Xi_2^{(0)}, \quad 3H_b^p = -4\mu_0\Xi_2^{(0)} + 1. \qquad (A.37)$$

or

$$3T_b^p = \frac{1}{3\tilde{K}_i + 4\mu_0}, \quad 3H_b^p = 3\tilde{K}_i \cdot 3T_b^p, \quad \tilde{K}_i = K_i \frac{1}{1 + 3\tilde{K}_i / [\hat{\lambda}r_i]}. \qquad (A.38)$$

### A.3. Dilute concentration terms for shear modulus

#### A.3.1. Elastic

For these properties we apply the following displacement-strain field at far distance [65]:

$$u_0^{(x,y,z)} = \begin{bmatrix} x \\ -y \\ 0 \end{bmatrix}, \quad \varepsilon_0 = \begin{bmatrix} 1 \\ -1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \qquad (A.39)$$

In spherical coordinates, the displacement is transformed to the vector

$$u_0^{(r,\theta,\phi)} = \begin{bmatrix} r \sin^2 \theta \cos 2\phi \\ r \sin \theta \cos \theta \cos 2\phi \\ -r \sin \theta \sin 2\phi \end{bmatrix}. \qquad (A.40)$$

For these boundary conditions, the important displacements and stresses at each phase are given by the general expressions

$$u_r^{(q)} = U_r^{(q)}(r) \sin^2 \theta \cos 2\phi, \qquad (A.41)$$

$$u_{\theta}^{(q)} = U_{\theta}^{(q)}(r) \sin \theta \cos \theta \cos 2\phi, \qquad (A.42)$$

$$u_{\phi}^{(q)} = -U_{\phi}^{(q)}(r) \sin \theta \sin 2\phi, \qquad (A.43)$$

$$\sigma_{rr}^{(q)} = \Sigma_{rr}^{(q)}(r) \sin^2 \theta \cos 2\phi, \qquad (A.44)$$

$$\sigma_{r\theta}^{(q)} = \Sigma_{r\theta}^{(q)}(r) \sin \theta \cos \theta \cos 2\phi, \qquad (A.45)$$

$$\sigma_{r\phi}^{(q)} = -\Sigma_{r\phi}^{(q)}(r) \sin \theta \sin 2\phi, \qquad (A.46)$$

for $q=0, i$, where

$$U_r^{(q)}(r) = \Xi_1^{(q)} r + \frac{9K_q - 6\mu_q}{15K_q + 11\mu_q} \Xi_2^{(q)} \frac{r^3}{r_i^2} - \frac{3}{2} \Xi_3^{(q)} \frac{r_i^5}{r^4} + \frac{3K_q + 3\mu_q}{2\mu_q} \Xi_4^{(q)} \frac{r_i^3}{r^2}, \qquad (A.47)$$

$$U_{\theta}^{(q)}(r) = \Xi_1^{(q)} r + \Xi_2^{(q)} \frac{r^3}{r_i^2} + \Xi_3^{(q)} \frac{r_i^5}{r^4} + \Xi_4^{(q)} \frac{r_i^3}{r^2}, \qquad (A.48)$$

$$\Sigma_{rr}^{(q)}(r) = 2\mu_q\Xi_1^{(q)} - \mu_q\frac{9K_q - 6\mu_q}{15K_q + 11\mu_q}\Xi_2^{(q)}\frac{r^2}{r_i^2} + 12\mu_q\Xi_3^{(q)}\frac{r_i^5}{r^5} - [9K_q + 4\mu_q]\Xi_4^{(q)}\frac{r_i^3}{r^3}, \qquad (A.49)$$

$$\Sigma_{r\theta}^{(q)}(r) = 2\mu_q\Xi_1^{(q)} + \mu_q\frac{48K_q + 10\mu_q}{15K_q + 11\mu_q}\Xi_2^{(q)}\frac{r^2}{r_i^2} - 8\mu_q\Xi_3^{(q)}\frac{r_i^5}{r^5} + 3K_q\Xi_4^{(q)}\frac{r_i^3}{r^3}. \qquad (A.50)$$

The unknowns that need to be defined are the $\Xi_1^{(i)}$, $\Xi_2^{(i)}$, $\Xi_3^{(i)}$, $\Xi_4^{(i)}$, $\Xi_1^{(0)}$, $\Xi_2^{(0)}$, $\Xi_3^{(0)}$, $\Xi_4^{(0)}$. The boundary and interface conditions that we need to satisfy are

$$u_r^{(i)}, u_{\theta}^{(i)} \text{ finite at } r = 0 \to \Xi_3^{(i)} = \Xi_4^{(i)} = 0, \qquad (A.51)$$

$$u_r^{(0)}(r \to \infty, \theta) = r \sin^2 \theta \cos(2\phi) \to \Xi_1^{(0)} = 1, \quad \Xi_2^{(0)} = 0, \qquad (A.52)$$

$$t_r = \lambda \llbracket u_r \rrbracket \to \Sigma_{rr}^{(i)}(r_i) = \lambda [U_r^{(0)}(r_i) - U_r^{(i)}(r_i)], \qquad (A.53)$$

$$t_{\theta} = \lambda \llbracket u_{\theta} \rrbracket \to \Sigma_{r\theta}^{(i)}(r_i) = \lambda [U_{\theta}^{(0)}(r_i) - U_{\theta}^{(i)}(r_i)], \qquad (A.54)$$

$$\llbracket \sigma_{rr}(r_i, \theta, \phi) \rrbracket = 0 \to \sigma_{rr}^{(0)}(r_i, \theta, \phi) = \sigma_{rr}^{(i)}(r_i, \theta, \phi), \qquad (A.55)$$

$$\llbracket \sigma_{r\theta}(r_i, \theta, \phi) \rrbracket = 0 \to \sigma_{r\theta}^{(0)}(r_i, \theta, \phi) = \sigma_{r\theta}^{(i)}(r_i, \theta, \phi). \qquad (A.56)$$

This linear system takes the form

$$K_s \cdot \Xi_s = F_s, \qquad (A.57)$$

where

$$K_s = \begin{bmatrix} K_{11} & K_{12} & K_{13} & K_{14} \\ K_{21} & K_{22} & K_{23} & K_{24} \\ K_{31} & K_{32} & K_{33} & K_{34} \\ K_{41} & K_{42} & K_{43} & K_{44} \end{bmatrix}, \qquad (A.58)$$

$$K_{11} = 1 + \frac{2\mu_i}{\lambda r_i}, \quad K_{12} = \frac{9K_i - 6\mu_i}{15K_i + 11\mu_i} - \frac{9K_i - 6\mu_i}{15K_i + 11\mu_i} \frac{\mu_i}{\lambda r_i}, \qquad (A.59)$$

$$K_{13} = \frac{3}{2}, \quad K_{14} = -\frac{3K_0 + 3\mu_0}{2\mu_0}, \qquad (A.60)$$

$$K_{21} = 1 + \frac{2\mu_i}{\lambda r_i}, \quad K_{22} = 1 + \frac{48K_i + 10\mu_i}{15K_i + 11\mu_i} \frac{\mu_i}{\lambda r_i}, \qquad (A.61)$$

$$K_{23} = -1, \quad K_{24} = -1, \qquad (A.62)$$

$$K_{31} = -2\mu_i, \quad K_{32} = \frac{9K_i - 6\mu_i}{15K_i + 11\mu_i} \mu_i, \qquad (A.63)$$

$$K_{33} = 12\mu_0, \quad K_{34} = -9K_0 - 4\mu_0, \qquad (A.64)$$

$$K_{41} = -2\mu_i, \quad K_{42} = -\frac{48K_i + 10\mu_i}{15K_i + 11\mu_i} \mu_i, \qquad (A.65)$$

$$K_{43} = -8\mu_0, \quad K_{44} = 3K_0, \qquad (A.66)$$

$$\Xi_s = \begin{bmatrix} \Xi_1^{(i)} \\ \Xi_2^{(i)} \\ \Xi_3^{(i)} \\ \Xi_4^{(i)} \end{bmatrix}, \qquad (A.67)$$

$$F_s = \begin{bmatrix} 1 \\ 1 \\ -2\mu_0 \\ -2\mu_0 \end{bmatrix}. \qquad (A.68)$$
After solving this linear system, the average strain and stress inside the particle are equal to

$$ 
\begin{aligned}
\langle \epsilon \rangle_{\Omega_i}^{+} &= \frac{1}{5r_i} \left[ 2U_r^{(0)}(r_i) + 3U_{\theta}^{(0)}(r_i) \right] \epsilon_0, \\
\langle \sigma \rangle_{\Omega_i}^{+} &= \frac{1}{5} \left[ 2\Sigma_{rr}^{(0)}(r_i) + 3\Sigma_{r\theta}^{(0)}(r_i) \right] \epsilon_0.
\end{aligned}
\quad (A.69) $$

Consequently, the shear dilute concentration terms are equal to

$$ 2T_s = 1 + \frac{3K_0 + 6\mu_0}{5\mu_0} \Xi_4^{(0)}, \quad 2H_s = 2\mu_0 - \frac{9K_0 + 8\mu_0}{5} \Xi_4^{(0)}, \quad (A.70) $$

or

$$ 
\begin{aligned}
2T_s = \frac{5\mu_0[3K_0 + 4\mu_0]}{\mu_0[9K_0 + 8\mu_0] + 6\tilde{\mu}_i[K_0 + 2\mu_0]}, \quad 2H_s = 2\tilde{\mu}_i \cdot 2T_s, \\
\tilde{\mu}_i = \mu_i \frac{1}{1 + 2\mu_i/[\hat{\lambda}r_i]}.
\end{aligned}
\quad (A.71) $$

### A.3.2. Inelastic

For these properties we apply zero displacement field at far distance. We also consider the following inelastic stress tensor

$$ \sigma_0^p = \begin{bmatrix} 1 \\ -1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \quad (A.72) $$

The displacement field is described by the same general forms with the elastic case. The matrix phase stress is expressed as

$$ 
\begin{aligned}
\sigma_{rr}^{(0)} &= [\Sigma_{rr}^{(0)}(r) + 1] \sin^2 \theta \cos 2\phi, \\
\sigma_{r\theta}^{(0)} &= [\Sigma_{r\theta}^{(0)}(r) + 1] \sin \theta \cos \theta \cos 2\phi, \\
\sigma_{r\phi}^{(0)} &= -[\Sigma_{r\phi}^{(0)}(r) + 1] \sin \theta \sin 2\phi.
\end{aligned}
\quad (A.73) $$

The unknowns that need to be defined are the $\Xi_1^{(i)}$, $\Xi_2^{(i)}$, $\Xi_3^{(i)}$, $\Xi_4^{(i)}$. The boundary and interface conditions that we need to satisfy are

$$ 
\begin{aligned}
u_r^{(i)}, u_\theta^{(i)} \text{ finite at } r=0 \rightarrow \Xi_3^{(i)} = \Xi_4^{(i)} = 0, \\
u_r^{(0)}(r \to \infty, \theta) = r \sin^2 \theta \cos(2\phi) \rightarrow \Xi_1^{(0)} = \Xi_2^{(0)} = 0,
\end{aligned}
\quad (A.74) $$

$$ 
\begin{aligned}
t_r &= \hat{\lambda} \llbracket u_r \rrbracket \rightarrow \Sigma_{rr}^{(i)}(r_i) = \hat{\lambda} [U_r^{(0)}(r_i) - U_r^{(i)}(r_i)], \\
t_\theta &= \hat{\lambda} \llbracket u_\theta \rrbracket \rightarrow \Sigma_{r\theta}^{(i)}(r_i) = \hat{\lambda} [U_\theta^{(0)}(r_i) - U_\theta^{(i)}(r_i)],
\end{aligned}
\quad (A.75) $$

$$ 
\begin{aligned}
\llbracket \sigma_{rr}(r_i, \theta, \phi) \rrbracket &= 0 \rightarrow \sigma_{rr}^{(0)}(r_i, \theta, \phi) = \sigma_{rr}^{(i)}(r_i, \theta, \phi), \\
\llbracket \sigma_{r\theta}(r_i, \theta, \phi) \rrbracket &= 0 \rightarrow \sigma_{r\theta}^{(0)}(r_i, \theta, \phi) = \sigma_{r\theta}^{(i)}(r_i, \theta, \phi).
\end{aligned}
\quad (A.76) $$

This linear system takes the form

$$ K_s \cdot \Xi_s = F_s^p, \quad (A.77) $$

where $K_s$ and $\Xi_s$ are the same with those of the elastic case and

$$ F_s^p = \begin{bmatrix} 0 \\ 0 \\ -1 \\ -1 \end{bmatrix}. \quad (A.78) $$

After solving this linear system, the average strain and stress inside the particle are equal to

$$ 
\begin{aligned}
\langle \epsilon \rangle_{\Omega_i}^{+} &= \frac{1}{5r_i} \left[ 2U_r^{(0)}(r_i) + 3U_{\theta}^{(0)}(r_i) \right] \sigma_0^p, \\
\langle \sigma \rangle_{\Omega_i}^{+} &= \frac{1}{5} \left[ 2\Sigma_{rr}^{(0)}(r_i) + 3\Sigma_{r\theta}^{(0)}(r_i) + 5 \right] \sigma_0^p.
\end{aligned}
\quad (A.79) $$

Consequently, the shear inelastic dilute concentration terms are equal to

$$ 2T_s^p = \frac{3K_0 + 6\mu_0}{5\mu_0} \Xi_4^{(0)}, \quad 2H_s^p = -\frac{9K_0 + 8\mu_0}{5} \Xi_4^{(0)} + 1, \quad (A.80) $$

or

$$ 2T_s^p = \frac{3[K_0 + 2\mu_0]}{\mu_0[9K_0 + 8\mu_0] + 6\tilde{\mu}_i[K_0 + 2\mu_0]}, \quad 2H_s^p = 2\tilde{\mu}_i \cdot 2T_s^p, \\
\tilde{\mu}_i = \mu_i \frac{1}{1 + 2\mu_i/[\hat{\lambda}r_i]}. \quad (A.81) $$

## Data availability

Data will be made available on request.

## References

[1] W.G. Buttlar, D. Bozkurt, G.G. Al-Khateeb, A.S. Waldhoff, Understanding asphalt mastic behavior through micromechanics, Transp. Res. Rec. 1681 (1) (1999) 157-169, http://dx.doi.org/10.3141/1681-19.

[2] P.S. Kandhal, C.Y. Lynn, F. Parker, Characterization tests for mineral fillers related to performance of asphalt paving mixtures, Transp. Res. Rec. 1638 (1) (1998) 101-110, http://dx.doi.org/10.3141/1638-12.

[3] Y. Bi, F. Guo, J. Zhang, J. Pei, R. Li, Correlation analysis between asphalt binder/asphalt mastic properties and dynamic modulus of asphalt mixture, Constr. Build. Mater. 276 (2021) 122256, http://dx.doi.org/10.1016/j.conbuildmat.2021.122256.

[4] C. Clopotel, R. Velasquez, H. Bahia, Measuring physico-chemical interaction in mastic using glass transition, Road Mater. Pavement Des. 13 (sup1) (2012) 304-320, http://dx.doi.org/10.1080/14680629.2012.657095.

[5] M. Kim, W.G. Buttlar, Stiffening mechanisms of asphalt-aggregate mixtures: From binder to mixture, Transp. Res. Rec. 2181 (1) (2010) 98-108, http://dx.doi.org/10.3141/2181-11.

[6] B.S. Underwood, Y.R. Kim, Experimental investigation into the multiscale behaviour of asphalt concrete, Int. J. Pavement Eng. 12 (4) (2011) 357-370, http://dx.doi.org/10.1080/10298436.2011.574136.

[7] W. Xu, X. Qiu, S. Xiao, H. Hong, F. Wang, J. Yuan, Characteristics and mechanisms of asphalt-filler interactions from a multi-scale perspective, Materials 13 (12) (2020) 2744, http://dx.doi.org/10.3390/ma13122744.

[8] D.A. Anderson, W.H. Goetz, Mechanical behavior and reinforcement of mineral filler-asphalt mixtures, Tech. rep, Purdue University, 1973, http://dx.doi.org/10.5703/1288824313845.

[9] A.F. Faheem, H.U. Bahia, Modelling of asphalt mastic in terms of filler-bitumen interaction, Road Mater. Pavement Des. 11 (sup1) (2010) 281-303, http://dx.doi.org/10.1080/14680629.2010.9690335.

[10] L.J. Walpole, A coated inclusion in an elastic medium, Math. Proc. Cambridge Philos. Soc. 83 (1978) 495-506, http://dx.doi.org/10.1017/S0305004100054773.

[11] A. Faheem, H. Wen, L. Stephenson, H. Bahia, Effect of mineral filler on damage resistance characteristics of asphalt binders, Asph. Paving Technol.-Proc. 77 (2008) 885.

[12] J. Zhang, Q. Shen, P. Ji, J. Wang, Y. Xiao, et al., Influence of filler type and rheological properties of asphalt mastic on the asphalt mastic-aggregate interaction, Materials 16 (2) (2023) 574, http://dx.doi.org/10.3390/ma16020574.

[13] D.N. Little, J.C. Petersen, Unique effects of hydrated lime filler on the performance-related properties of asphalt cements: Physical and chemical interactions revisited, J. Mater. Civ. Eng. 17 (2) (2005) 207-218, http://dx.doi.org/10.1061/(ASCE)0899-1561(2005)17:2(207).

[14] R.M. Alfaqawi, A. Fareed, S.B.A. Zaidi, G.D. Airey, A. Rahim, Effect of hydrated lime and other mineral fillers on stiffening and oxidative ageing in bitumen mastic, Constr. Build. Mater. 315 (2022) 125789, http://dx.doi.org/10.1016/j.conbuildmat.2021.125789.

[15] W.G. Buttlar, R. Roque, Evaluation of empirical and theoretical models to determine asphalt mixture stiffnesses at low temperatures, 1996, pp. 99-141, Asphalt Paving Technology: Association of Asphalt Paving Technologists- Proceedings of the Technical Sessions, 65.

[16] N. Shashidhar, A. Shenoy, On using micromechanical models to describe dynamic mechanical behavior of asphalt mastics, Mech. Mater. 34 (10) (2002) 657-669, http://dx.doi.org/10.1016/S0167-6636(02)00166-7.

[17] Y.-R. Kim, D. Little, Linear viscoelastic analysis of asphalt mastic, J. Mater. Civ. Eng. 16 (2) (2004) 122-132, http://dx.doi.org/10.1061/(ASCE)0899-1561(2004)16:2(122).

[18] Y.-R. Kim, D.N. Little, I. Song, Effect of mineral fillers on fatigue resistance and fundamental material characteristics: mechanistic evaluation, Transp. Res. Rec. 1832 (1) (2003) 1-8, http://dx.doi.org/10.3141/1832-01.

[19] B.S. Underwood, Y.R. Kim, A four phase micro-mechanical model for asphalt mastic modulus, Mech. Mater. 75 (2014) 13-33, http://dx.doi.org/10.1016/j.mechmatr.2014.04.001.
[20] A. Abbas, E. Masad, T. Papagiannakis, A. Shenoy, Modelling asphalt mastic stiffness using discrete element analysis and micromechanics-based models, Int. J. Pavement Eng. 6 (2) (2005) 137-146, http://dx.doi.org/10.1080/10298430500159040.

[21] J.D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 252 (1271) (1959) 561-569, http://www.jstor.org/stable/100791.

[22] Z. Hashin, Complex moduli of viscoelastic composites - I. General theory and application to particulate composites, Int. J. Solids Struct. 6 (5) (1970) 539-552, http://dx.doi.org/10.1016/0020-7683(70)90029-6.

[23] R. Christensen, Viscoelastic properties of heterogeneous media, J. Mech. Phys. Solids 17 (1) (1969) 23-41, http://dx.doi.org/10.1016/0022-5096(69)90011-8.

[24] R. Hill, A self-consistent mechanics of composite materials, J. Mech. Phys. Solids 13 (4) (1965) 213-222.

[25] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta Metall. 21 (5) (1973) 571-574.

[26] X. Gao, Z. Fan, J. Zhang, S. Liu, Micromechanical model for asphalt mixture coupling inter-particle effect and imperfect interface, Constr. Build. Mater. 148 (2017) 696-703, http://dx.doi.org/10.1016/j.conbuildmat.2017.05.015.

[27] Y. Jia, G. Liu, Y. Gao, S. Wang, Z. Li, X. Wang, L. Yan, Micromechanics prediction of dynamic modulus for asphalt mastic considering filler distribution characteristic: A new perspective, Constr. Build. Mater. 364 (2023) 129916, http://dx.doi.org/10.1016/j.conbuildmat.2022.129916.

[28] J. Xu, Z. Fan, G. Lu, D. Wang, P. Liu, Effects of bitumen origins and filler types on the micromechanics prediction of complex modulus of asphalt mastic considering interparticle interactions, J. Mater. Civ. Eng. 36 (1) (2024) 04023521, http://dx.doi.org/10.1061/JMCEE7.MTENG-16609.

[29] Y. Jia, Y. Li, Z. Wu, C. Si, J. Zhang, S. Wang, Y. Gao, Z. Li, Z. Wei, Evaluation of filler particle size effect of dynamic modulus for asphalt mastic with a new combined micromechanics approach, Constr. Build. Mater. 477 (2025) 141335.

[30] H. Pham Huy, E. Sanchez-Palencia, Phénomènes de transmission à travers des couches minces de conductivité élevée, J. Math. Anal. Appl. 47 (2) (1974) 284-309, http://dx.doi.org/10.1016/0022-247X(74)90023-7.

[31] Z. Hashin, Thermoelastic properties of particulate composites with imperfect interface, J. Mech. Phys. Solids 39 (6) (1991) 745-762, http://dx.doi.org/10.1016/0022-5096(91)90023-H.

[32] A. Klarbring, Derivation of a model of adhesively bonded joints by the asymptotic expansion method, Internat. J. Engrg. Sci. 29 (4) (1991) 493-512, http://dx.doi.org/10.1016/0020-7225(91)90090-P.

[33] P.B. Bövik, On the modelling of thin interface layers in elastic and acoustic scattering problems, Q. J. Mech. Appl. Math. 47 (1994) 17-42, http://dx.doi.org/10.1093/qjmam/47.1.17.

[34] Y. Benveniste, T. Miloh, Imperfect soft and stiff interfaces in two-dimensional elasticity, Mech. Mater. 33 (6) (2001) 309-323, http://dx.doi.org/10.1016/S0167-6636(01)00055-2.

[35] G.I. Barenblatt, The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks, J. Appl. Math. Mech. 23 (3) (1959) 622-636, http://dx.doi.org/10.1016/0021-8928(59)90157-1.

[36] D. Dugdale, Yielding of steel sheets containing slits, J. Mech. Phys. Solids 8 (2) (1960) 100-104, http://dx.doi.org/10.1016/0022-5096(60)90013-2.

[37] X.P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (1994) 1397-1434, http://dx.doi.org/10.1016/0022-5096(94)90003-5.

[38] M.G.A. Tijssens, B.L.J. Sluys, E. Van der Giessen, Numerical simulation of quasi-brittle fracture using damaging cohesive surfaces, Eur. J. Mech. A Solids 19 (2000) 761-779, http://dx.doi.org/10.1016/S0997-7538(00)00190-X.

[39] M.J. van den Bosch, P.J.G. Schreurs, M.G.D. Geers, An improved description of the exponential xu and needleman cohesive zone law for mixed-mode decohesion, Eng. Fract. Mech. 73 (9) (2006) 1220-1234, http://dx.doi.org/10.1016/j.engfracmech.2005.12.006.

[40] M. Fagerström, R. Larsson, Theory and numerics for finite deformation fracture modelling using strong discontinuities, Internat. J. Numer. Methods Engrg. 66 (6) (2006) 911-948, http://dx.doi.org/10.1002/nme.1573.

[41] K. Park, G.H. Paulino, J.R. Roesler, A unified potential-based cohesive model of mixed-mode fracture, J. Mech. Phys. Solids 57 (6) (2009) 891-908, http://dx.doi.org/10.1016/j.jmps.2008.10.003.

[42] J. Mosler, I. Scheider, A thermodynamically and variationally consistent class of damage-type cohesive models, J. Mech. Phys. Solids 59 (8) (2011) 1647-1668, http://dx.doi.org/10.1016/j.jmps.2011.04.012.

[43] R. Dimitri, M. Trullo, L. De Lorenzis, G. Zavarise, Coupled cohesive zone models for mixed-mode fracture: A comparative study, Eng. Fract. Mech. 148 (2015) 145-179, http://dx.doi.org/10.1016/j.engfracmech.2015.09.029.

[44] N. Despringre, Y. Chemisky, K. Bonnay, F. Meraghni, Micromechanical modeling of damage and load transfer in particulate composites with partially debonded interface, Compos. Struct. 155 (2016) 77-88, http://dx.doi.org/10.1016/j.compstruct.2016.06.075.

[45] M.E. Gurtin, A.I. Murdoch, A continuum theory of elastic material surfaces, Arch. Ration. Mech. Anal. 57 (4) (1975) 291-323, http://dx.doi.org/10.1007/BF00261375.

[46] A. Javili, A. McBride, P. Steinmann, Thermomechanics of solids with lower-dimensional energetics: On the importance of surface, interface and curve structures at the nanoscale, a unifying review, Appl. Mech. Rev. 65 (1) (2013) 010802, http://dx.doi.org/10.1115/1.4023012.

[47] H.L. Duan, J. Wang, Z.P. Huang, B.L. Karihaloo, Size-dependent effective elastic constants of solids containing nano-inhomogeneities with interface stress, J. Mech. Phys. Solids 53 (7) (2005) 1574-1596, http://dx.doi.org/10.1016/j.jmps.2005.02.009.

[48] R. Dingreville, J. Qu, Interfacial excess energy, excess stress and excess strain in elastic solids: Planar interfaces, J. Mech. Phys. Solids 56 (5) (2008) 1944-1954, http://dx.doi.org/10.1016/j.jmps.2007.11.003.

[49] P. Chhapadia, P. Mohammadi, P. Sharma, Curvature-dependent surface energy and implications for nanostructures, J. Mech. Phys. Solids 59 (10) (2011) 2103-2115, http://dx.doi.org/10.1016/j.jmps.2011.06.007.

[50] A. Javili, A. McBride, J. Mergheim, P. Steinmann, U. Schmidt, Micro-to-macro transitions for continua with surface structure at the microscale, Int. J. Solids Struct. 50 (16-17) (2013) 2561-2572, http://dx.doi.org/10.1016/j.ijsolstr.2013.03.022.

[51] G. Chatzigeorgiou, A. Javili, P. Steinmann, Multiscale modeling for composites with energetic interfaces at the micro- or nanoscale, Math. Mech. Solids 20 (9) (2015) 1130-1145, http://dx.doi.org/10.1177/1081286513516122.

[52] Q. Chen, G. Chatzigeorgiou, F. Meraghni, A. Javili, Homogenization of size-dependent multiphysics behavior of nanostructured piezoelectric composites with energetic surfaces, Eur. J. Mech. A Solids 96 (2022) 104731, http://dx.doi.org/10.1016/j.euromechsol.2022.104731.

[53] Y. Benveniste, A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media, J. Mech. Phys. Solids 54 (4) (2006) 708-734, http://dx.doi.org/10.1016/j.jmps.2005.10.009.

[54] S.T. Gu, J.T. Liu, Q.C. He, Size-dependent effective elastic moduli of particulate composites with interfacial displacement and traction discontinuities, Int. J. Solids Struct. 51 (13) (2014) 2283-2296, http://dx.doi.org/10.1016/j.ijsolstr.2014.02.033.

[55] G. Chatzigeorgiou, F. Meraghni, A. Javili, Generalized interfacial energy and size effects in composites, J. Mech. Phys. Solids 106 (2017) 257-282, http://dx.doi.org/10.1016/j.jmps.2017.06.002.

[56] S. Firooz, G. Chatzigeorgiou, F. Meraghni, A. Javili, Homogenization accounting for size effects in particulate composites due to general interfaces, Mech. Mater. 139 (2019) 103204, http://dx.doi.org/10.1098/rspa.1957.0133.

[57] S. Firooz, G. Chatzigeorgiou, F. Meraghni, A. Javili, Bounds on size effects in composites via homogenization accounting for general interfaces, Contin. Mech. Thermodyn. 32 (1) (2020) 173-206, http://dx.doi.org/10.1007/s00161-019-00796-w.

[58] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 241 (1226) (1957) 376-396, http://dx.doi.org/10.1098/rspa.1957.0133.

[59] Y. Benveniste, G.J. Dvorak, T. Chen, Stress fields in composites with coated inclusions, Mech. Mater. 7 (1989) 305-317, http://dx.doi.org/10.1016/0167-6636(89)90021-5.

[60] Z. Wang, R.J. Oelkers, K.C. Lee, F.T. Fisher, Annular coated inclusion model and applications for polymer nanocomposites - part I: Spherical inclusions, Mech. Mater. 101 (2016) 107-184, http://dx.doi.org/10.1016/j.mechmat.2016.07.004.

[61] Z. Wang, R.J. Oelkers, K.C. Lee, F.T. Fisher, Annular coated inclusion model and applications for polymer nanocomposites - part II: Cylindrical inclusions, Mech. Mater. 101 (2016) 50-60, http://dx.doi.org/10.1016/j.mechmat.2016.07.005.

[62] D.C. Lagoudas, A.C. Gavazzi, H. Nigam, Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme, Comput. Mech. 8 (1991) 193-203, http://dx.doi.org/10.1007/BF00372689.

[63] J. Chaboche, P. Kanoute, A. Ross, On the capabilities of mean-field approaches for the description of plasticity in metal matrix composites, Int. J. Plast. 21 (2005) 1409-1434, http://dx.doi.org/10.1016/j.ijplas.2004.07.001.

[64] H.L. Duan, X. Yi, Z.P. Huang, J. Wang, A unified scheme for prediction of effective moduli of multiphase composites with interface effects. Part I: Theoretical framework, Mech. Mater. 39 (1) (2007) 81-93, http://dx.doi.org/10.1016/j.mechmat.2006.02.009.
[65] G. Chatzigeorgiou, F. Meraghni, N. Charalambakis, Multiscale Modeling Approaches for Composites, Elsevier, Amsterdam, 2022, http://dx.doi.org/10.1016/C2019-0-05214-4.

[66] Z. Hashin, The elastic moduli of heterogeneous materials, J. Appl. Mech. 29 (1) (1962) 143–150, http://dx.doi.org/10.1115/1.3636446.

[67] F. Praud, G. Chatzigeorgiou, J. Bikard, F. Meraghni, Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation, Mech. Mater. 114 (2017) 9–29, http://dx.doi.org/10.1016/j.mechmat.2017.07.001.

[68] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer-Verlag, New York, 1998, http://dx.doi.org/10.1007/b98904.

[69] G. Chatzigeorgiou, N. Charalambakis, Y. Chemisky, F. Meraghni, Thermomechanical Behavior of Dissipative Composite Materials, ISTE Press - Elsevier, London, 2018, http://dx.doi.org/10.1016/C2017-0-01481-7.

[70] B. Miled, I. Doghri, L. Delannay, Coupled viscoelastic-viscoplastic modeling of homogeneous and isotropic polymers: Numerical algorithm and analytical solutions, Comput. Methods Appl. Mech. Engrg. 200 (47–48) (2011) 3381–3394, http://dx.doi.org/10.1016/j.cma.2011.08.015.

[71] S. Satouri, R. Chekkour, G. Chatzigeorgiou, F. Meraghni, G. Robert, Numerical-experimental approach to identify the effect of relative humidity on the material parameters of a rate-dependent damage model for polyamide 66, Mech. Mater. 184 (2023) 104735, http://dx.doi.org/10.1016/j.mechmat.2023.104735.