A multiscale micromechanical progressive elastic-damage model for cementitious composites featuring superabsorbent polymer (SAP)

Abstract

A multiscale micromechanics-based progressive damage model is developed to investigate the overall mechanical behavior and the interfacial microcrack evolutions of the cementitious composites featuring superabsorbent polymer (SAP) under uniaxial tension. Elastic properties, progressive damage process, and homogenization procedure of cementitious composites are systematically integrated in this model. The effective elastic moduli of the composites are determined based on a multiscale micromechanical framework. According to the small strain assumption, the total strain tensor and the elastic-damage compliance tensor are additively decomposed into elastic and damage-induced components. The damage-induced strains and compliances are then deduced from micromechanics. To characterize the progressive elastic-damage induced by microcracks, stages of microcrack propagation are identified from the interface contact stress and the matrix cleavage stress. The complex potentials and stress intensity factors for kinked interface cracks are derived from the distributed dislocations method. By implementing the homogenization process, the macroscopic mechanical behavior is obtained from the micro/mesoscale. The results indicate that the material parameters have clear mechanical significance. Different parameters, such as the SAP addition ratio, aggregate content, initial interfacial crack size, and initial interfacial crack location, are revealed to be influential in the overall mechanical behavior of the composites. The proposed model can be generalized to other particle-reinforced composites with different constituent properties, which can potentially contribute to the design and optimization of durable composites.

Introduction

The applications of polymeric materials in cementitious composites have grown steadily over the past decades. Due to the excellent water absorption and retention capacity, superabsorbent polymers (SAPs) have been widely employed in cementitious composites to improve properties and broaden application areas (Jung et al., 2020; Mignon et al., 2017; Schröfl et al., 2022). When in service, the mechanical properties of composites are critical for structural deformation, development of internal defects, and life cycle determination. To gain insight into the mechanical properties of cementitious composites featuring SAP, considerable experimental studies have been conducted. It has been shown that SAP can effectively mitigate autogenous shrinkage and cracking (Justs et al., 2015; Wyrzykowski et al., 2011), enhance freeze-thaw resistance (Kim et al., 2021; Laustsen et al., 2015; Xu et al., 2022), and promote self-healing ability (Kua et al., 2019; Park and Choi, 2018; Snoeck et al., 2014). However, most of these studies focused on performance enhancement at the macroscale, ignoring the interconnections of complex microstructures. Cementitious composites, as a multiscale material, have various compositions from the nanoscale to the macroscale, and their macroscopic mechanical behaviors are determined by the constituent phases and their interactions at each scale. Therefore, it is necessary to utilize some multiscale methods to capture the material and mechanical characteristics at each scale and explore the relevant correlations among different scales to effectively obtain the overall macroscopic response of the composites.

Recently, based on continuum micromechanics, many studies have explored the different properties of cementitious composites using multiscale approaches, effectively realizing the correlations between different scales. Yan et al. (2019) proposed an equivalent elastoplastic damage model that can predict the macroscale response of hybrid fiber-reinforced composites from the behavior of single fibers at the mesoscale. Zhang et al. (2020) presented a four-scale micromechanical model that effectively estimated the damage of hybrid-fiber concrete at or after elevated temperatures. Chen et al. (2020) developed a multiscale framework for investigating the macroscopic equivalent elastic properties of ultra-high-performance concrete. Liu et al. (2022) simulated and characterized the freeze-thaw damage and overall macroscopic mechanical behavior of shotcrete at multiple scales. Not only cementitious composites (Luo et al., 2020; Park et al., 2022; Rivas Murillo et al., 2016; Voyiadjis et al., 2022; Zhu et al., 2015, 2016), multiscale models based on continuum micromechanics can also be of great help in studying and developing the properties of other composites (Chen et al., 2021; Görthofer et al., 2022; Ivica et al., 2022; Jia et al., 2015; Ju et al., 2012; Lou et al., 2020; Marcin et al., 2011; Onodera et al., 2022; Panamoottil et al., 2018; Satouri et al., 2022; Yanase and Ju, 2014; Wei et al., 2023; Wu and Ju, 2017; Zhang et al., 2021, 2022; Zhao et al., 2023).

Investigations into the material properties and damage mechanisms of cementitious composites have revealed that microcracks or microdefects develop internally during their fabrication, especially at the aggregate-matrix interface. When subjected to external loading, interfacial microcracks will generate damages at lower scales. From damage to failure, interfacial cracks will undergo different stages, including aggregates debonding from the matrix, interfacial cracks deviating into the matrix from the aggregate-matrix interface, and kinked cracks propagating in the matrix. For various damage processes in cementitious composites, experimental studies or full-scale simulations are costly. A multiscale progressive damage model that can describe the different damage stages is warranted to examine the damage mechanism of cementitious materials in depth. Currently, research on progressive damage modeling is mainly focused on other types of composites (Fu and Reiner, 2022; Haj-Ali, 2009; Ju and Ko, 2008; Ju and Wu, 2016; Khayyam Rayeni et al., 2020; Labeas et al., 2012; Lu et al., 2019; Tay et al., 2008; Wei et al., 2022; Zhang and Zhao, 2012), and there are few studies on the overall progressive damage process (debonding, cracking, propagation) of cementitious composites. To the best of the authors’ knowledge, most of the existing numerical studies on cementitious composites simplified the models to a certain extent. For example, only aggregate-matrix debonding or crack propagation in the matrix is modeled; characterization of the elastic properties of the composites is neglected; parameters of the initial microcracks are simplified (Dutta and Kishen, 2018; Huang et al., 2011; Li et al., 2022; Wu and Tang, 2019). These simplifications need to be improved. Thus, it is beneficial and necessary to develop a multiscale progressive damage model capable of capturing the damage evolutions at lower scales and predicting their effects on the overall macroscopic mechanical behaviors of composites.

In this work, a multiscale micromechanics-based progressive elastic-damage model is presented to study the overall mechanical behavior and the interfacial microcrack evolutions of the cementitious composites featuring SAP under uniaxial tension. To characterize the progressive damage, the total strain tensor and the elastic-damage compliance tensor are decomposed into an elastic component and a damage-induced component, assuming small strains. By systematically integrating the elastic properties, progressive damage processes, and homogenization processes of composites, this model can effectively predict the macroscopic mechanical responses of cementitious composites from the known properties of micro/mesoscale components and analyze the effects of different material parameters on uniaxial tensile strength. Based on this model, we can also propose adaptations and modifications depending on underlying material and damage characteristics to improve the model and extend it to other particle-reinforced composites with various constituent properties.

This paper is structured as follows: The micromechanical basis for formulating the effective macroscopic stress-strain law and the homogenization process is described in the next section. The progressive damage model used to describe the propagation evolutions of interfacial microcracks at various stress stages is detailed in Section ‘Progressive damage model’. Section ‘Computational algorithms’ outlines the computational algorithm for the proposed multiscale progressive damage model. The model is validated by comparison with experimental data in Section ‘Model validation’. A parametric study is performed in Section ‘Parametric study’ to examine the effect of key parameters on the uniaxial tensile strength. Finally, Section ‘Conclusions’ summarizes the conclusions of this work.

Micromechanics basis

In this section, the micromechanics basis is described to formulate the effective macroscopic stress-strain law and the homogenization process of cementitious composites, which includes the governing equations, Ju and Chen’s (1994a, 1994b) micromechanical framework, and the homogenization scheme for cementitious composites featuring SAP.

Governing equations

In this study, the ensemble-volume averaging process was carried out in a rational representative volume element (RVE) to represent the effective constitutive relationship and damage process of cementitious composites. For a composite RVE comprised of elastic inhomogeneities distributed in an elastic matrix, when it is subjected to specified far-field stress ${\sigma }_{ij}^0$ and strain ${\epsilon }_{ij}^0$ , the volume-averaged stress tensor $\overline{\mathbf{\sigma }}$ is (Ju and Chen, 1994a, 1994b):

$$\overline{\mathbf{\sigma }}\equiv \frac{1}{V}{\displaystyle \underset{V}{\int }\mathbf{\sigma }(\mathbf{x})}\text{d}\mathbf{x}=\frac{1}{V}\left[{\displaystyle \underset{V_m}{\int }\mathbf{\sigma }(\mathbf{x})}\text{d}\mathbf{x}+{\displaystyle \sum _{r=1}^n{\displaystyle \underset{V_r}{\int }\mathbf{\sigma }(\mathbf{x})}\text{d}\mathbf{x}}\right]={\varphi }_m{\overline{\mathbf{\sigma }}}_m\hspace{0.17em}+\hspace{0.17em}{\displaystyle \sum _{r=1}^n{\varphi }_r{\overline{\mathbf{\sigma }}}_r}$$

(1)

Further, the volume-averaged strain tensor $\overline{\mathbf{\epsilon }}$ is given by:

$$\overline{\mathbf{\epsilon }}\equiv \frac{1}{V}{\displaystyle \underset{V}{\int }\mathbf{\epsilon }(\mathbf{x})}\text{d}\mathbf{x}=\frac{1}{V}\left[{\displaystyle \underset{V_m}{\int }\mathbf{\epsilon }(\mathbf{x})}\text{d}\mathbf{x}+{\displaystyle \sum _{r=1}^n{\displaystyle \underset{V_r}{\int }\mathbf{\epsilon }(\mathbf{x})}\text{d}\mathbf{x}}\right]={\varphi }_m{\overline{\mathbf{\epsilon }}}_m\hspace{0.17em}+\hspace{0.17em}{\displaystyle \sum _{r=1}^n{\varphi }_r{\overline{\mathbf{\epsilon }}}_r}$$

(2)

where $V$ is the volume of the composite RVE, $V_m$ and $V_r$ are the volumes of the matrix and the r-th phase inhomogeneities, and $n$ is the number of inhomogeneities. ${\varphi }_i$ denotes the volume fractions, ${\overline{\mathbf{\sigma }}}_i$ and ${\overline{\mathbf{\epsilon }}}_i$ signify the volume-averaged stress and strain, respectively, and subscripts $\left(i=m,r\right)$ denote the matrix and r-th phase inhomogeneities.

Within the framework of the “homogenization” concept for the effective properties and constitutive relations of heterogeneous composites, the macroscopic stress-strain law and the overall elastic-damage compliance tensor S are given by (Ju, 1991):

$$\overline{\mathbf{\epsilon }}=\hspace{0.17em}\mathbf{\text{S}}:\overline{\mathbf{\sigma }};\mathbf{\text{ S}}\equiv \left[{\mathbf{\text{S}}}^{\mathbf{\text{e}}}\cdot \left(\mathbf{\text{I}}+\mathbf{\text{D}}\right)\right]={\mathbf{\text{S}}}^{\mathbf{\text{e}}}+{\mathbf{\text{S}}}^{\mathbf{\text{e}}}\cdot \mathbf{\text{D}}$$

(3)

where ${\mathbf{\text{S}}}^{\mathbf{\text{e}}}$ is the effective elastic (undamaged) compliance tensor of the composite, $\mathbf{\text{I}}$ is the fourth-rank unit tensor, and $\mathbf{\text{D}}$ denotes the damage tensor.

To characterize the damage evolution during the loading process, we introduce a fourth-rank tensor ${\mathbf{\text{S}}}^{*}$ to signify the damage-induced compliance, and the relationship between ${\mathbf{\text{S}}}^{*}$ and $\mathbf{\text{D}}$ can be formulated as follows:

$${\mathbf{\text{S}}}^{*}\equiv {\mathbf{\text{S}}}^{\mathbf{\text{e}}}\cdot \mathbf{\text{D}}$$

(4)

Given that the effective elastic (undamaged) compliance tensor ${\mathbf{\text{S}}}^{\mathbf{\text{e}}}$ is constant, it follows that ${\mathbf{\text{S}}}^{*}$ and $\mathbf{\text{D}}$ are in a one-to-one correspondence. Therefore, if the damage-induced compliance tensor ${\mathbf{\text{S}}}^{*}$ can be derived micromechanically, then the fourth-rank damage tensor $\mathbf{\text{D}}$ can be explicitly expressed by micromechanics, and we can characterize the evolution of the damage tensor $\mathbf{\text{D}}$ by tracking the damage-induced compliance tensor ${\mathbf{\text{S}}}^{*}$ .

During the damage process, the total strain tensor $\overline{\mathbf{\epsilon }}$ and the elastic-damage compliance tensor $\mathbf{\text{S}}$ are applicable to the additive decomposition based on the small strain assumption:

$$\overline{\mathbf{\epsilon }}={\overline{\mathbf{\epsilon }}}^e\hspace{0.17em}+\hspace{0.17em}{\overline{\mathbf{\epsilon }}}^{*};\mathbf{S}={\mathbf{S}}^{\text{e}}+{\mathbf{\text{S}}}^{*}$$

(5)

where ${\overline{\mathbf{\epsilon }}}^e$ and ${\overline{\mathbf{\epsilon }}}^{*}$ are the elastic strain and the damage-induced strain, respectively. ${\mathbf{\text{S}}}^{*}$ is the damage-induced compliance tensor.

Therefore, the relationship between macroscopic stress-strain can be defined as:

$$\overline{\mathbf{\epsilon }}=\hspace{0.17em}\mathbf{\text{S}}:\overline{\mathbf{\sigma }}=\left[{\mathbf{\text{S}}}^{\mathbf{\text{e}}}{\text{+\hspace{0.17em}S}}^{\text{*}}\right]:\overline{\mathbf{\sigma }}=\hspace{0.17em}{\mathbf{\text{S}}}^{\mathbf{\text{e}}}:\overline{\mathbf{\sigma }}\hspace{0.17em}+\hspace{0.17em}{\overline{\mathbf{\epsilon }}}^{*}$$

(6)

It is noted that this work assumes the matrix and the inclusions are isotropic and homogeneous, and they remain linear elastic in all phases of damage loading. The inelastic behavior results from damage, and the inelastic damage-induced strain ${\overline{\mathbf{\epsilon }}}^{*}$ is caused by displacement jumps between the crack faces of the microcracks.

With the mix proportions of the cementitious composites, the volume fractions of the aggregate and matrix can be determined. Once the stress and strain conditions of the aggregate and matrix are micromechanically derived, one can obtain the volume-averaged stress and strain. Then, if the effective elastic (undamaged) compliance and the damage-induced compliance are derived by micromechanics, the macroscopic stress-strain relationships will be achieved.

Ju and Chen’s micromechanical framework

To obtain the macroscopic stress-strain relationship, the overall elastic properties of the multiphase composites need to be calculated. Ju and Chen (1994a, 1994b) proposed a micromechanical framework to evaluate the effective elastic modulus of composites in an ensemble-volume averaged RVE. The framework considers direct interactions between particles and assumes all randomly distributed particles are of identical shape and impenetrable. Within this framework, the governing equations for determining the effective elastic moduli of composites with n-phase inhomogeneities can be formulated as:

$$\overline{\mathbf{\sigma }}=\hspace{0.17em}{\mathbf{C}}_0:\left(\overline{\mathbf{\epsilon }}-{\displaystyle \sum _{r=1}^n{\varphi }_r{\overline{\mathbf{\epsilon }}}_r^{\ast }}\right)$$

(7)

$$\overline{\mathbf{\epsilon }}=\hspace{0.17em}{\mathbf{\epsilon }}_0+{\displaystyle \sum _{r=1}^n{\varphi }_r\mathbf{B}:{\overline{\mathbf{\epsilon }}}_r^{\ast }}$$

(8)

$$\left(-{\mathbf{A}}_r-\mathbf{s}\right):{\overline{\mathbf{\epsilon }}}_r^{\ast }={\mathbf{\epsilon }}_0+{\overline{\mathbf{\epsilon }}}_r^{\prime p}$$

(9)

$$\left(1-{\displaystyle \sum _{r=1}^n{\varphi }_r}\right){\overline{\mathbf{\epsilon }}}_r^{\prime m}+{\displaystyle \sum _{r=1}^n{\varphi }_r}{\overline{\mathbf{\epsilon }}}_r^{\prime p}={\displaystyle \sum _{r=1}^n{\varphi }_r}\left(\mathbf{B}-\mathbf{s}\right):{\overline{\mathbf{\epsilon }}}_r^{\ast }=\mathbf{0}$$

(10)

where ${\mathbf{C}}_0$ is the stiffness tensor of the matrix, ${\overline{\mathbf{\epsilon }}}^{*}$ is the distributed eigenstrain, ${\mathbf{\epsilon }}_0$ is the uniform strain field induced by far-field loads for a homogeneous matrix material only, $\mathbf{B}$ is the fourth-order constant tensor (when all particles are unidirectionally aligned and identically shaped, $\mathbf{B}=\mathbf{s}$ ), $\mathbf{s}$ is the Eshelby tensor, ${\overline{\mathbf{\epsilon }}}_r^{\prime p}$ indicates the interaction effects of inter-particles, ${\overline{\mathbf{\epsilon }}}_r^{\prime m}$ corresponds to the ensemble-volume averaged perturbed strain in the matrix due to the existence of all particles. ${\mathbf{A}}_r$ is the elastic mismatch tensor, which can be written as:

$${\mathbf{A}}_r\equiv {\left({\mathbf{C}}_r-{\mathbf{C}}_0\right)}^{-1}\cdot {\mathbf{C}}_0$$

(11)

${\mathbf{C}}_r$ is the fourth order stiffness tensor of inclusion phase $r$ . When ${\overline{\mathbf{\epsilon }}}_r^{\prime p}$ is determined, the effective modulus of elasticity follows from the governing equations.

For two-phase composites with randomly dispersed spherical inhomogeneities, the effective stiffness tensor ${\mathbf{C}}_e$ reads (Ju and Chen, 1994a):

$${\mathbf{C}}_e={\mathbf{C}}_0\cdot \left\{\mathbf{I}-{\varphi }_1\mathbf{\Gamma }\cdot {\left(-{\mathbf{A}}_1-\mathbf{s}+{\varphi }_1\mathbf{s}\cdot \mathbf{\Gamma }\right)}^{-1}\right\}$$

(12)

The fourth-rank tensor ${\mathbf{A}}_1$ and $\mathbf{\Gamma }$ are defined as follows:

$${\mathbf{A}}_1={\left({\mathbf{C}}_1-{\mathbf{C}}_0\right)}^{-1}\cdot {\mathbf{C}}_0$$

(13)

$${\Gamma }_{\mathit{ijkl}}={\gamma }_1{\delta }_{ij}{\delta }_{kl}+{\gamma }_2\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)$$

(14)

where ${\mathbf{C}}_1$ is the stiffness tensor of the inclusion, other parameters in equation (14) are described in Appendix 1.

Regarding three-phase composites with two different types of randomly dispersed spherical inhomogeneities, the effective stiffness tensor ${\mathbf{C}}_e$ is derived as (Ko and Ju, 2012, 2013, 2021; Lin and Ju, 2009):

$${\mathbf{C}}_e={\mathbf{C}}_0\cdot \left\{\mathbf{I}-{\varphi }_1{\mathbf{\Gamma }}^1\cdot {\left({\mathbf{T}}^1\right)}^{-1}-{\varphi }_2{\mathbf{\Gamma }}^2\cdot {\left({\mathbf{T}}^2\right)}^{-1}\right\}$$

(15)

where

$${\mathbf{T}}^1=\left(-{\mathbf{A}}_1-\mathbf{s}+{\varphi }_1\mathbf{s}\cdot {\mathbf{\Gamma }}^1+{\varphi }_2\mathbf{s}\cdot {\mathbf{\Gamma }}^2\cdot {\left({\mathbf{A}}_2+\mathbf{s}\right)}^{-1}\cdot \left({\mathbf{A}}_1+\mathbf{s}\right)\right)$$

(16)

$${\mathbf{T}}^2=\left(-{\mathbf{A}}_2-\mathbf{s}+{\varphi }_2\mathbf{s}\cdot {\mathbf{\Gamma }}^2+{\varphi }_1\mathbf{s}\cdot {\mathbf{\Gamma }}^1\cdot {\left({\mathbf{A}}_1+\mathbf{s}\right)}^{-1}\cdot \left({\mathbf{A}}_2+\mathbf{s}\right)\right)$$

(17)

$${\text{Γ}}^1{\hspace{0.17em}}_{\mathit{ijkl}}=r_1{\delta }_{ij}{\delta }_{kl}+r_2\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)$$

(18)

$${\text{Γ}}^2{\hspace{0.17em}}_{\mathit{ijkl}}=r_3{\delta }_{ij}{\delta }_{kl}+r_4\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)$$

(19)

Other parameters in equations (16) to (19) are detailed in Appendix 2.

Therefore, it is possible to calculate the effective elastic properties of the composites based on the mechanical parameters, geometrical parameters, and volume fractions of the constituents.

Homogenization scheme for cementitious composites featuring SAP

For homogenization and elastic properties prediction, the microstructure of cementitious composites featuring SAP is divided into four scales, which are macroscale (concrete), mesoscale (mortar), microscale (cement paste), and nanoscale (calcium silicate hydrate crystals (C-S-H)), as depicted in Figure 1. The RVE for each scale is detailed as follows:

Nanoscale (C-S-H): C-S-H can be regarded as a two-phase composite on the nanoscale, with the high-density (HD) C-S-H embedded in the low-density (LD) C-S-H. From the perspective of composites, LD C-S-H is considered the matrix phase, and HD C-S-H is the inclusion.

Microscale (cement paste): At the microscale, cement paste mainly consists of cement hydration products (C-S-H, calcium hydroxide crystals (CH)), unhydrated cement clinkers, and microporous phases (capillary pores/water). As C-S-H has the highest content, cement paste is regarded as a four-phase composite consisting of C-S-H as the matrix, CH, unhydrated cement clinkers, and microporous phases as inclusions.

Mesoscale (mortar): It was observed that the scale lengths of saturated SAP particles and SAP voids are mesoscale (Schröfl et al., 2017). Therefore, mortar is viewed as a three-phase composite consisting of a homogenized cement paste matrix, sand particles, and SAP particles/voids. For simplification, the interfacial transition zone (ITZ) is not considered a separate mesoscale phase.

Macroscale (concrete): Concrete at the macroscale is considered as a two-phase composite with coarse aggregates embedded in the homogenized mortar matrix.

Figure 1.

Homogenization scheme for cementitious composites featuring SAP.

According to the microstructure division of cementitious composites featuring SAP, the effective elastic properties can be evaluated by a four-step homogenization procedure. It is assumed that the inclusions in the homogenization procedure are spherical.

At the nanoscale, equation (12) is employed to homogenize the two-phase composite of C-S-H. The homogenized C-S-H then becomes the matrix of the microscale RVE.

At the microscale, a sub-step is first performed to homogenize the hydration products by equation (12), and the homogenized hydration products is regarded as the matrix for cement paste. The three-phase cement paste can then be homogenized by equation (15) and becomes the matrix for mesoscale RVE.

At the mesoscale, the three-phase mortar is homogenized through equation (15) and becomes the matrix of the RVE at the macroscale.

At the macroscale, homogenization of concrete can be achieved using equation (12) to obtain the effective elastic properties of concrete featuring SAP.

The homogenization procedure is also presented in Figure 1 for a structured representation.

Progressive damage model

It was observed that the interface between matrix-inhomogeneities of cementitious materials is prone to the formation of microcracks. Under external loading, damage is initiated from interfacial microcracks at lower scales and transmitted to higher scales. The damage stages include debonding of aggregate particles from the matrix, deviation of microcracks from the interface to the matrix, and propagation of the deviated cracks in the matrix. In this section, based on micromechanics, the stress-strain relationships for an infinite matrix containing a single inhomogeneity with an interfacial arc crack are derived for different stages of far-field uniaxial tension. The results can subsequently be upscaled to higher scales by the homogenization procedure in Section ‘Governing equations’ to obtain the macroscopic stress-strain responses.

Damage-induced inelastic strains and compliances

For a single circular inhomogeneity with elastic properties $E^a$ and ${\nu }^a$ partially embedded in an infinite matrix with elastic properties $E^m$ and ${\nu }^m$ , Toya (1974) derived the stress and displacement solutions for an interfacial arc crack under far-field uniaxial tension. The problem for consideration is schematically exhibited in Figure 2. In the diagram, x−y and x′−y′ denote the global and local Cartesian coordinate systems, and the local polar coordinate system $\theta$ is measured from the x′-axis. In particular, $\alpha$ is the half-angle of the interfacial arc crack, $\varphi$ is the angle between applied tension and crack normal direction, and $a$ is the radius of the inhomogeneity. With remote uniaxial tension $p$ applied in the y-axis direction, the displacement components $u_r^{\prime }$ and $u_{\theta }^{\prime }$ of the crack faces are given in complex variable form by the solution of Toya (1974):

$$\begin{array}{l}{u}_r^{\prime }+i{u}_{\theta }^{\prime }=-A_{1}pa\left[\sin \frac{1}{2}(\alpha -\theta )\sin \frac{1}{2}(\alpha +\theta )\right]{}^{\frac{1}{2}}\\ \hspace{1em}\times \left\{G_0+i{H}_0-\frac{1}{k}-\frac{2\left(1-k\right)}{k}\exp \left[i\left(2\varphi -\theta \right)+2\lambda \left(\alpha -\pi \right)\right]\right\}\\ \hspace{1em}\times \exp \left[\lambda \left(\pi -\alpha \right)-i\left\{\frac{1}{2}\theta -\lambda \ln \left[\sin \frac{1}{2}(\alpha -\theta )/\sin \frac{1}{2}(\alpha +\theta )\right]\right\}\right]\end{array}$$

(20)

where

$$G_0=\frac{1-\left(\cos \alpha +2\lambda \sin \alpha \right)\exp \left[2\lambda \left(\pi -\alpha \right)\right]+\left(1-k\right)\left(1+4{\lambda }^2\right){\sin }^2\alpha \cos 2\varphi }{2-k-k\left(\cos \alpha +2\lambda \sin \alpha \right)\exp \left[2\lambda \left(\pi -\alpha \right)\right]}$$

(21)

$$H_0=-\frac{\left(1-k\right)\left(1+4{\lambda }^2\right){\sin }^2\alpha \sin 2\varphi }{k\left\{1+\left(\cos \alpha +2\lambda \sin \alpha \right)\exp \left[2\lambda \left(\pi -\alpha \right)\right]\right\}}$$

(22)

Figure 2.

Model of a single circular inclusion partially embedded in an infinite matrix with different elastic properties under far-field uniaxial tension.

Transforming the displacement components in the local polar coordinate system to the local Cartesian coordinate system:

$$u_x^{\prime }=u_r^{\prime }\cos \theta -u_{\theta }^{\prime }\sin \theta ;u_y^{\prime }=u_r^{\prime }\sin \theta +u_{\theta }^{\prime }\cos \theta$$

(23)

The inelastic strain components induced by a single arc microcrack are of the form:

$${\epsilon }_1^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_r^{\prime }{\cos }^2\theta -u_{\theta }^{\prime }\sin \theta \cos \theta \right)}d\theta$$

(24)

$${\epsilon }_2^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_r^{\prime }{\sin }^2\theta +u_{\theta }^{\prime }\sin \theta \cos \theta \right)}d\theta$$

(25)

$${\epsilon }_6^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_r^{\prime }\sin 2\theta -u_{\theta }^{\prime }\cos 2\theta \right)}d\theta$$

(26)

where $A$ is the two-dimensional surface area of the inhomogeneity.

To obtain the damage-induced compliance ${\mathbf{\text{S}}}^{*}$ , Ju (1991) derived the expressions for $u_{ri}^{\prime }$ and $u_{\theta i}^{\prime }$ :

$$u_{r1}^{\prime }=-B_1\left\{\begin{array}{l}\left[\begin{array}{l}\frac{1-\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}+(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{2-k-k\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}}\\ -\frac{1}{k}-\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \end{array}\right]\cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ \hspace{1em}+\left[\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\right]\end{array}\right\}$$

(27)

$$u_{\theta 1}^{\prime }=-B_1\left\{\begin{array}{l}-\left[\begin{array}{l}\frac{1-\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}+(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{2-k-k\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}}\\ -\frac{1}{k}-\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \end{array}\right]\sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ \hspace{1em}+\left[\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\right]\end{array}\right\}$$

(28)

$$u_{r2}^{\prime }=-B_1\left\{\begin{array}{l}\left[\begin{array}{l}\frac{1-\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}-(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{2-k-k\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}}\\ -\frac{1}{k}+\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \end{array}\right]\cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ \hspace{1em}-\left[\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\right]\end{array}\right\}$$

(29)

$$u_{\theta 2}^{\prime }=-B_1\left\{\begin{array}{l}-\left[\begin{array}{l}\frac{1-\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}-(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{2-k-k\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}}\\ -\frac{1}{k}+\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \end{array}\right]\sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ \hspace{1em}-\left[\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\right]\end{array}\right\}$$

(30)

$$u_{r6}^{\prime }=-2B_1\left\{\begin{array}{l}-\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ \hspace{1em}-\left[\frac{(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{k\left[1+\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}\right]}-\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \right]\\ \sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right))\end{array}\right\}$$

(31)

$$u_{\theta 6}^{\prime }=-2B_1\left\{\begin{array}{l}\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\sin \theta \sin \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\\ -\left[\frac{(1-k)(1+4{\lambda }^2){\sin }^2\alpha }{k\left[1+\left(\cos \alpha +2\lambda \sin \alpha \right)e^{2\lambda (\pi -\alpha )}\right]}-\frac{2(1-k)}{k}{e}^{2\lambda (\pi -\alpha )}\cos \theta \right]\cos \left(\frac{1}{2}\theta -\lambda \ln \left[\frac{\sin \frac{1}{2}(\alpha -\theta )}{\sin \frac{1}{2}(\alpha +\theta )}\right]\right)\end{array}\right\}$$

(32)

$$B_1=A_{1}a{\left[\sin \frac{1}{2}(\alpha -\theta )\sin \frac{1}{2}(\alpha +\theta )\right]}^{1/2}{e}^{\lambda \left(\pi -a\right)}$$

(33)

$A_1$ , $k$ , and $\lambda$ are material parameters defined as follows (Toya, 1974):

$$A_1=\frac{k}{4}\left(\frac{1+{\kappa }_m}{{\mu }_m}+\frac{1+{\kappa }_a}{{\mu }_a}\right);{\kappa }_m=3-4{\nu }_m;{\kappa }_a=3-4{\nu }_a$$

(34)

$$k=\frac{\beta }{1+\upsilon };\lambda =-\frac{\ln \upsilon }{2\pi }$$

(35)

$$\beta =\frac{{\mu }_m\left(1+{\kappa }_a\right)}{{\mu }_m+{\kappa }_m{\mu }_a};\upsilon =\frac{{\mu }_a+{\kappa }_a{\mu }_m}{{\mu }_m+{\kappa }_m{\mu }_a}$$

(36)

where ${\mu }_i\left(i=m,a\right)$ and ${\nu }_i$ are the shear modulus and Poisson's ratio, respectively, and subscripts $m$ and $a$ denote the matrix and inhomogeneity.

Therefore, the local damage-induced compliances by a single microcrack propagation are expressed as:

$$S_{1i}^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_{ri}^{\prime }{\cos }^2\theta -u_{\theta i}^{\prime }\sin \theta \cos \theta \right)}d\theta$$

(37)

$$S_{2i}^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_{ri}^{\prime }{\sin }^2\theta +u_{\theta i}^{\prime }\sin \theta \cos \theta \right)}d\theta$$

(38)

$$S_{6i}^{*\prime }=\frac{a}{A}{\displaystyle {\int }_{-\alpha }^{\alpha }\left(u_{ri}^{\prime }\sin 2\theta -u_{\theta i}^{\prime }\cos 2\theta \right)}d\theta$$

(39)

Moreover, the inelastic strain components ${\epsilon }_i^{*\prime }$ can be written in the form ( $i,j=1,2,6$ ):

$${\epsilon }_i^{*\prime }=S_{ij}^{*\prime }{\sigma }_j^{\prime }$$

(40)

${\sigma }_j^{\prime }$ are the local stress components, which can be calculated as:

$${\sigma }_1^{\prime }={\sigma }_{xx}^{\prime }=p{\cos }^2\varphi ;{\sigma }_2^{\prime }={\sigma }_{yy}^{\prime }=p{\sin }^2\varphi ;{\sigma }_6^{\prime }={\sigma }_{xy}^{\prime }=\frac{p}{2}\sin 2\varphi$$

(41)

To determine the average aggregate stress, we can integrate the traction along the interface, given the local stresses:

$${\sigma }_{xx}^{\prime a}=\frac{a}{A}{\displaystyle {\int }_{\alpha }^{2\pi -\alpha }{T}_1}\cos \theta d\theta$$

(42)

$${\sigma }_{yy}^{\prime a}=\frac{a}{A}{\displaystyle {\int }_{\alpha }^{2\pi -\alpha }{T}_2}\sin \theta d\theta$$

(43)

$${\sigma }_{xy}^{\prime a}=\frac{a}{2A}{\displaystyle {\int }_{\alpha }^{2\pi -\alpha }(T_1}\sin \theta +T_2\cos \theta )d\theta$$

(44)

where $T_i={\sigma }_{ij}^{\prime }{n}_j^{\prime }$ is the traction component.

Crack path determination

Under remote tensile loading, an interfacial microcrack can either expand along the interface, leading to more inclusion-matrix interfacial debonding, or it can veer from the interface and enter the matrix. The crack propagation path selection is mainly determined by the local stress conditions near the current crack tip. In this section, we utilize the interface contact stress ${\sigma }_{rr}$ and the matrix cleavage stress ${\sigma }_{\phi \phi }$ to determine the crack path. It is notable that in this study, the stiffness of the inhomogeneity is assumed to be higher than the matrix, and therefore crack penetration through the aggregate is neglected.

For crack propagation along the interface by $\Delta \alpha$ , the contact stress ${\sigma }_{rr}$ at the inclusion-matrix interface can be calculated as follows (Dutta and Kishen, 2018; Toya, 1974):

$$2\left({\sigma }_{rr}+i{\sigma }_{r\theta }\right)=-\frac{p}{4}ka\left(1-g\right)i\left(1+2i\lambda \right)\sin \alpha N_0{\left[a\exp i\theta -a\exp i\alpha \right]}^{\frac{1}{2}i\lambda }$$

(45)

where

$$\begin{array}{l}{N}_0=\left\{G_0-\frac{1}{k}-\frac{2}{k}\left(1-k\right)\upsilon \exp \left(2\alpha \lambda \right)\cos \left(2\varphi -\alpha \right)\right\}\\ +i\left\{H_0-\frac{2}{k}\left(1-k\right)\upsilon \exp \left(2\alpha \lambda \right)\sin \left(2\varphi -\alpha \right)\right\}\end{array}$$

(46)

$G_0$ and $H_0$ are the same as equations (21) and (22), $g$ is a material parameter defined in equation (72).

The matrix cleavage stress ${\sigma }_{\phi \phi }$ is calculated at a short distance $c$ from the interfacial crack tip. Figure 3 presents the deviated crack of length $c$ and angular location $\phi$ .

Figure 3.

Model of a single circular inclusion with a deviated crack of length $c$ and angle $\phi$ near the interfacial crack tip.

The matrix cleavage stress ${\sigma }_{\phi \phi }$ is given by (Toya, 1974):

$${\sigma }_{\phi \phi }=\frac{1}{2}\left({\sigma }_{xx}^{\prime }+{\sigma }_{yy}^{\prime }\right)+\frac{1}{2}\left({\sigma }_{xx}^{\prime }-{\sigma }_{yy}^{\prime }\right)\cos \left[2\left(\alpha -\phi \right)\right]+{\sigma }_{xy}^{\prime }\sin \left[2\left(\alpha -\phi \right)\right]$$

(47)

where ${\sigma }_{xx}^{\prime }$ , ${\sigma }_{yy}^{\prime }$ , and ${\sigma }_{xy}^{\prime }$ are local stress components. The quantities ${\sigma }_{xx}^{\prime }+{\sigma }_{yy}^{\prime }$ , ${\sigma }_{xx}^{\prime }-{\sigma }_{yy}^{\prime }$ , and ${\sigma }_{xy}^{\prime }$ can be calculated as follows:

$$\frac{1}{p}{\left(c/a\right)}^{\frac{1}{2}}({\sigma }_{xx}^{\prime }+{\sigma }_{yy}^{\prime })=\mathrm{Re}\left[\frac{ik{M}_0\exp \left[-\lambda \left(\alpha -\phi \right)\right]}{{\left(2\sin \alpha \right)}^{\frac{1}{2}}}\exp \left\{i\left(\frac{\phi -\alpha }{2}+\lambda \ln x_1\right)\right\}\right]$$

(48)

$$\frac{1}{p}{\left(c/a\right)}^{\frac{1}{2}}({\sigma }_{yy}^{\prime }-{\sigma }_{xx}^{\prime }+2i{\sigma }_{xy}^{\prime })=\frac{k}{2{\left(2\sin \alpha \right)}^{\frac{1}{2}}}\times \left[\begin{array}{l}-\left(1-2i\lambda \right)M_0\sin \phi \exp \left\{-\lambda \left(\alpha -\phi \right)+i{M}_1\right\}\\ +i{M}_0\exp \left\{-\lambda \left(\alpha -\phi \right)+i{M}_2\right\}\\ +\frac{i{\overline{M}}_0}{\upsilon }\exp \left\{-\lambda \left(\alpha +\phi \right)+i{M}_3\right\}\end{array}\right]$$

(49)

where

$$M_0=\left(i-2\lambda \right)\sin \alpha \left\{G_0+i{H}_0-\frac{1}{k}-\frac{2}{k}\left(1-k\right)\upsilon \exp \left[i\left(2\varphi -\alpha \right)+2\alpha \lambda \right]\right\}$$

(50)

$$M_1=-\frac{5}{2}\alpha +\frac{3}{2}\phi +\lambda \ln x_1$$

(51)

$$M_2=-\frac{5}{2}\alpha +\frac{1}{2}\phi +\lambda \ln x_1$$

(52)

$$M_3=-\frac{3}{2}\alpha +\frac{1}{2}\phi -\lambda \ln x_1$$

(53)

$$x_1=\frac{c}{2a\sin \alpha }$$

(54)

Once the interface contact stress ${\sigma }_{rr}$ and the matrix cleavage stress ${\sigma }_{\phi \phi }$ are obtained, we can compare them to the interface tensile strength $f_{\mathrm{int}}$ and the matrix tensile strength $f_m$ , respectively. If ${\sigma }_{rr}\le f_{\mathrm{int}}$ , the crack will be stationary. If ${\sigma }_{rr}>f_{\mathrm{int}}$ and ${\sigma }_{\phi \phi }\le f_m$ , the crack will propagate along the inclusion-matrix interface. Lastly, if ${\sigma }_{\phi \phi }>f_m$ , the crack will deviate into the matrix at a perpendicular angle relative to the maximum cleavage stress.

Kinked crack propagation

By comparing the stress conditions at the crack tip with the interface and matrix tensile strength, the criteria under which microcracks extend along the interface or deviate into the matrix are given. As a kinked interface crack enters the matrix and propagates, it is considered a normal sharp crack in the matrix, and the stress intensity factor will govern its propagation.

To analyze the kinked interface crack, it can be modeled as a distributed dislocation (Fang et al., 2003). Assume a circular inhomogeneity with an interfacial arc crack with an edge dislocation located in the matrix, as depicted in Figure 4. The elastic field relations within the medium can be described using Muskhelishvili’s complex potentials of $\Phi \left(z\right)$ and $\Psi \left(z\right)$ .

Figure 4.

Model of a single circular inclusion with an interfacial arc crack with an edge dislocation.

The stress and displacement fields expressed in complex potentials are given by:

$${\sigma }_{rr}+{\sigma }_{\theta \theta }=4\mathrm{Re}\left[\Phi \left(z\right)\right]$$

(55)

$${\sigma }_{rr}+i{\sigma }_{r\theta }=\Phi \left(z\right)+\overline{\Phi \left(z\right)}-\overline{z}\overline{\left\{\frac{d\left[\Phi \left(z\right)\right]}{dz}\right\}}-\frac{\overline{z}}{z}\overline{\Psi \left(z\right)}$$

(56)

$$2\mu \left(\frac{\partial u_x}{\partial \theta }+\frac{\partial u_y}{\partial \theta }\right)=iz\left[\kappa \Phi \left(z\right)-\overline{\Phi \left(z\right)}+\overline{z}\overline{\left\{\frac{d\left[\Phi \left(z\right)\right]}{dz}\right\}}+\frac{\overline{z}}{z}\overline{\Psi \left(z\right)}\right]$$

(57)

where $z=x+iy$ , $\mu$ is the shear modulus, and the over-bar indicates the complex conjugate. ${\sigma }_{rr}$ , ${\sigma }_{\theta \theta }$ , and ${\sigma }_{r\theta }$ are the components of the stress in polar coordinates, $u_x$ and $u_y$ are the components of the displacement in the Cartesian coordinates.

For the problem of an interfacial arc crack interacting with an edge dislocation positioned at an arbitrary point $z_0$ in a homogeneous matrix, the solutions for complex potentials of the matrix and inhomogeneity were derived by Fang et al. (2003) as follows:

$${\Phi }_m\left(z\right)=\frac{K}{1-g}\left[G\left(z\right)+D_0\right]-\frac{X_0\left(z\right)}{1-g}\left[G_0\left(z\right)+G_{\infty }\left(z\right)+G_{z_0}\left(z\right)+G_{z^{*}}\left(z\right)\right]$$

(58)

$${\Phi }_a\left(z\right)=-\frac{h}{1-g}\left[G\left(z\right)+D_0\right]+\frac{X_0\left(z\right)}{1-g}\left[G_0\left(z\right)+G_{\infty }\left(z\right)+G_{z_0}\left(z\right)+G_{z^{*}}\left(z\right)\right]$$

(59)

$${\Psi }_m\left(z\right)=\frac{a^2}{z^2}\left[{\Phi }_m\left(z\right)+\overline{{\Phi }_m}\left(\frac{a^2}{z}\right)-z\frac{d\left[{\Phi }_m\left(z\right)\right]}{dz}\right]$$

(60)

$${\Psi }_a\left(z\right)=\frac{a^2}{z^2}\left[{\Phi }_a\left(z\right)+\overline{{\Phi }_a}\left(\frac{a^2}{z}\right)-z\frac{d\left[{\Phi }_a\left(z\right)\right]}{dz}\right]$$

(61)

where

$$G\left(z\right)=\frac{{\gamma }_2}{z-z_0}+\frac{{\gamma }_2}{z}-\frac{{\gamma }_2}{z-z^{\ast }}+\frac{\overline{{\gamma }_2}{z}^{\ast }\left(z_0-z^{\ast }\right)}{{\overline{z}}_0{\left(z-z^{\ast }\right)}^2}\text{,}{z}^{\ast }=\frac{a^2}{{\overline{z}}_0}$$

(62)

$$G_0\left(z\right)=\frac{h{\gamma }_2}{z}\frac{1}{X_0\left(0\right)}$$

(63)

$$G_{\infty }\left(z\right)=\left[z-\frac{a}{2}\left(\exp \left[i{\alpha }_1\right]+\exp \left[i{\alpha }_2\right]\right)-ia\beta \left(\exp \left[i{\alpha }_1\right]-\exp \left[i{\alpha }_2\right]\right)\right]K{D}_0+h{\gamma }_2$$

(64)

$$G_{z_0}\left(z\right)=\frac{h{\gamma }_2}{X_0\left(z_0\right)}\frac{1}{z-z_0}$$

(65)

$$G_{z^{*}}\left(z\right)=\frac{h}{X_0\left(z^{\ast }\right)}\left\{\frac{\overline{{\gamma }_2}{z}^{\ast }\left(z_0-z^{\ast }\right)}{{\overline{z}}_0{\left(z-z^{\ast }\right)}^2}\left[1-\left(z-z^{\ast }\right)\frac{X_0^{\prime }\left(z^{\ast }\right)}{X_0\left(z^{\ast }\right)}\right]-\frac{{\gamma }_2}{z-z^{\ast }}\right\}$$

(66)

and

$$D_0=\frac{Q{Q}_1-\overline{Q_1}}{1-Q^2}$$

(67)

$$Q=\frac{-K}{1-g}{X}_0\left(0\right)a\left[\frac{1}{2}\left(\exp \left[i{\alpha }_1\right]+\exp \left[i{\alpha }_2\right]\right)+i\beta \left(\exp \left[i{\alpha }_1\right]-\exp \left[i{\alpha }_2\right]\right)\right]-\frac{K}{1-g}+1$$

(68)

$$\begin{array}{l}{Q}_1=\frac{h{\gamma }_2}{1-g}\frac{X_0^{\prime }\left(0\right)}{X_0\left(0\right)}+\frac{h{X}_0\left(0\right)}{1-g}\left\{\frac{\overline{{\gamma }_2}\left(z_0-z^{\ast }\right)}{{\overline{z}}_0{X}_0\left(z^{\ast }\right)}\left[\frac{{\overline{z}}_0}{a^2}+\frac{X_0^{\prime }\left(z^{\ast }\right)}{X_0\left(z^{\ast }\right)}\right]+{\gamma }_2-\frac{{\gamma }_2}{z_0{X}_0\left(z_0\right)}+\frac{{\gamma }_2}{z^{\ast }{X}_0\left(z^{\ast }\right)}\right\}\\ \hspace{1em}\hspace{0.17em}\hspace{1em}-\frac{h}{1-g}\left[\frac{{\gamma }_2}{z^{\ast }}+\frac{\overline{{\gamma }_2}\left(z_0-z^{\ast }\right)}{{\overline{z}}_0{z}^{\ast }}-\frac{{\gamma }_2}{z_0}\right]\end{array}$$

(69)

$$X_0\left(z\right)={\left(z-a\exp \left[i{\alpha }_1\right]\right)}^{-\frac{1}{2}-i\beta }\cdot {\left(z-a\exp \left[i{\alpha }_2\right]\right)}^{-\frac{1}{2}+i\beta },\beta =\frac{\ln \left|g\right|}{2\pi }$$

(70)

$${\gamma }_2=\frac{{\mu }_m}{\pi \left(1+{\kappa }_m\right)}\left(b_y-i{b}_x\right)$$

(71)

$\left(b_x,b_y\right)$ is Burger’s vector for the edge dislocation. The material parameters are defined as follows:

$$g=-\frac{{\mu }_m+{\kappa }_m{\mu }_a}{{\mu }_a+{\kappa }_a{\mu }_m},h=-\frac{{\mu }_a\left(1+{\kappa }_m\right)}{{\mu }_a+{\kappa }_a{\mu }_m},K=\frac{{\mu }_m\left(1+{\kappa }_a\right)}{{\mu }_a+{\kappa }_a{\mu }_m}$$

(72)

With the complex potentials being known, the components of stress and displacement may be determined.

To model the kinked crack propagation in the matrix as an ordinary sharp crack, the unknown dislocation distributions $b_x$ and $b_y$ are needed to obtain the complex potentials. As the crack propagation is determined by the stress intensity factor, the unknown dislocation distributions can be achieved after calculating the stress intensity factor at the kinked crack tip. The stress intensity factors at the kinked crack tip take the form:

$$\frac{K_I+i{K}_{II}}{p\sqrt{c}}=\frac{1}{\sqrt{2}}\left(b_{\theta }+i{b}_r\right)$$

(73)

where $c$ is the current crack length, $b_r$ and $b_{\theta }$ are the dislocation distributions in polar coordinates, $K_I$ and $K_{II}$ are the stress intensity factors corresponding to Mode I and Mode II loading, respectively. $K_I$ , $K_{II}$ are calculated as follows:

$$K_I=\left[{\sigma }_{xx}^{\prime }\sqrt{\left(\pi c\right)}\right]{\phi }_n{\sin }^2{\beta }_0+\left[{\sigma }_{yy}^{\prime }\sqrt{\left(\pi c\right)}\right]{\phi }_n\cos \left(\frac{\pi }{2}-{\beta }_0\right)\sin \left(\frac{\pi }{2}-{\beta }_0\right)$$

(74)

$$K_{II}=\left[{\sigma }_{yy}^{\prime }\sqrt{\left(\pi c\right)}\right]{\phi }_n{\sin }^2\left(\frac{\pi }{2}-{\beta }_0\right)+\left[{\sigma }_{xx}^{\prime }\sqrt{\left(\pi c\right)}\right]{\phi }_n\cos {\beta }_0\sin {\beta }_0$$

(75)

where

$${\phi }_n=0.707-0.18\eta +6.55{\eta }^2-10.54{\eta }^3+6.85{\eta }^4$$

(76)

with $\eta ={\left(1+\frac{c}{a}\right)}^{-1}$ and ${\beta }_0=\alpha +\frac{\pi }{2}-\phi$ .

For the kinked interface crack, an iteration process can be performed to capture the crack evolution for the current stress conditions and crack length. The corresponding stress conditions and crack length can be determined from the complex potentials and dislocation distributions.

Based on micromechanics, the formulations of the progressive damage model for different damage stages are derived. With these formulations, one can determine the stresses, strains, and damage-induced compliances at different damage stages from the material parameters, geometric parameters, and stress conditions of the matrix and aggregate, which will be further employed in the homogenization procedure.

Consequently, integrating the homogenization procedure of the macroscopic stress-strain relationship and effective elastic properties in Section ‘Micromechanics basis’ with the progressive damage model derived in Section ‘Progressive damage model’, it is possible to predict the macroscopic mechanical response of cementitious composites based on the known properties of the micro/mesoscale constituents without making assumptions for the essential properties of the composites.

Computational algorithms

In this section, the computational algorithm is detailed for the proposed multiscale micromechanical progressive elastic-damage model. First, the effective elastic compliance of the composite is obtained through the micromechanical framework in Section ‘Ju and Chen’s micromechanical framework’ for subsequent homogenization procedure calculations. Considering that cementitious composites experience various damage stages when subjected to external loading, the numerical results for different stages are derived from incremental analysis or an iteration process. At the stage of crack propagation along the inclusion-matrix interface, the stress-strain evolutions are obtained by incremental analysis, where the corresponding stress, strain, and compliance are calculated for different semi-debond crack angle increments. When the crack deviates into the matrix and propagates, numerical results are acquired by performing an iteration on the stress and crack length. For each step, the numerical results are scaled up by employing the homogenization procedure in Section ‘Governing equations’ to obtain macroscopic stress-strain relationships. The elaborated calculation schemes are as follows:

Determine the effective elastic compliance ${\mathbf{\text{S}}}^{\mathbf{\text{e}}}$ ( ${\mathbf{\text{S}}}^{\mathbf{\text{e}}}={\mathbf{C}}_e{\hspace{0.17em}}^{-1}$ ) through the micromechanical framework in Section ‘Ju and Chen’s micromechanical framework’.

Given the initial (current) value for the half-angle $\alpha$ of the interfacial arc crack.

Compute the interfacial crack-opening displacement components $u_{ri}^{\prime }$ and $u_{\theta i}^{\prime }$ for the current crack angle from equations (27) to (32).

Calculate the interface contact stress ${\sigma }_{rr}$ for the current crack angle from equation (45) and the matrix cleavage stress ${\sigma }_{\phi \phi }$ at a close point to the interfacial crack tip from equation (47).

Identify the crack path by comparing the interface contact stress ${\sigma }_{rr}$ and the matrix cleavage stress ${\sigma }_{\phi \phi }$ with the interface tensile strength $f_{\mathrm{int}}$ and the matrix tensile strength $f_m$ , respectively.

• If ( ${\sigma }_{rr}\le f_{\mathrm{int}}$ ) or ( ${\sigma }_{rr}>f_{\mathrm{int}}$ and ${\sigma }_{\phi \phi }\le f_m$ ), go to step (xi). Meanwhile, increase $\alpha$ by $\Delta \alpha$ (i.e., $\alpha =\alpha +\Delta \alpha$ ) and return to step (ii).

• If ( ${\sigma }_{\phi \phi }>f_m$ ), go to step (vi).

•    (vi) Obtain the current stress conditions and given the initial (current) crack length $c$ .

•   (vii) Compute the stress intensity factors $K_I$ and $K_{II}$ for the current stress conditions and crack length $c$ according to equations (74) and (75).

• (viii) Obtain the dislocation distributions ( $b_x$ , $b_y$ ) and the complex potentials ( $\Phi \left(z\right)$ , $\Psi \left(z\right)$ ) for the current crack propagation using equations (73) and (58) to (61).

•     (ix) Determine the current stress ( ${\sigma }_{rr}$ , ${\sigma }_{r\theta }$ ) and displacement ( $u_x$ , $u_y$ ) components using the complex potentials from equations (55) to (57).

•   (x) For current stress conditions and crack length, perform the following sub-steps simultaneously:

• Return current stress conditions and crack length to step (vi).

• Go to step (xi).

• (xi) Utilize the homogenization procedure in Section ‘Governing equations’ to acquire the macroscopic stress-strain relationship at each step.

For convenience, the computational algorithm is summarized in the flowchart below (Figure 5).

Figure 5.

Computational algorithm for the multiscale micromechanical progressive damage model.

Model validation

To validate the proposed multiscale micromechanical progressive elastic-damage model, experimental data from three independent studies were utilized for comparison with model predictions. Two sets are direct uniaxial tensile test results of plain concrete specimens (Gopalaratnam and Shah, 1985; López et al., 2008), and one is on SAP internally cured concrete (Xie et al., 2020). In addition, various incremental sizes, initial iteration parameters, and simulation input parameters were adopted for different experimental data in model predictions to assess and verify the convergence and stability of the proposed model.

In the simulations, concrete was regarded as a matrix-aggregate two-phase composite. First, the effective elastic properties of the equivalent concrete specimens were calculated from the multiscale micromechanical model detailed in Section ‘Ju and Chen’s micromechanical framework’. Then, the aggregates were considered as inhomogeneities distributed in the equivalent matrix, and the macroscopic stress-strain relationship of the composite was obtained using the progressive damage model in Section ‘Progressive damage model’. The material and mechanical parameters of the specimen components were given in the literature. For plain concrete specimens, the initial half-angle of the interfacial arc crack was taken as 10°, the half-angle increment $\Delta \alpha$ was set to $0.1\alpha$ , and the initial kinked crack length $c$ was assumed to be $0.1a$ . While the initial half-angle of the interfacial arc crack was also adopted to be 10° for the SAP internally cured concrete, the half-angle increment $\Delta \alpha$ was adjusted to $0.05\alpha$ , and the initial kinked crack length $c$ was considered to be $0.05a$ . The effects of SAP on the macroscopic mechanical response of composites were considered only for its role in the homogenization process of the equivalent matrix, and its influence on the damage or healing process was neglected at this stage, which will be the subject of our forthcoming papers.

Figure 6 displays the macroscopic stress-strain relationships obtained by the proposed multiscale micromechanical progressive damage model and the comparison with the experimental results (Gopalaratnam and Shah, 1985; López et al., 2008). It is evident that the proposed model can well simulate the trend of the stress-strain state at different stages and reasonably predict the constitutive relationships of plain concrete under uniaxial tension.

Figure 6.

Comparison with the experimental results for plain concrete under uniaxial tension.

The applicability of the presented model in predicting the macroscopic mechanical behavior of SAP-incorporated cementitious composites under uniaxial tension was verified by the experimental data obtained by Xie et al. (2020) for SAP internally cured concrete. This model is capable of estimating the variation of concrete tensile strength with SAP ratio satisfactorily, with the relative deviations between model predictions and the test data less than 10%, as displayed in Figure 7. As observed in Figure 7, the predicted tensile strengths are marginally greater than the experimental outcomes. This is mainly due to the model at this stage ignoring the effects of SAP particles on the damage or healing process, which will be addressed in our forthcoming studies.

Figure 7.

Comparison with the experimental results of SAP internally cured concrete with different SAP ratios (mass ratio of SAP to cement).

To verify the convergence of the prediction results and the numerical stability of the computational algorithm, the model-obtained stress-strain curves of different concrete specimens were analyzed. Figure 8 exhibits the macroscopic stress-strain responses of SAP internally cured concrete. Combined with the stress-strain curves of the plain concrete in Figure 6, it is observed that the numerical results obtained from the model predictions all converge to stable values at different increments and initial iteration parameters, which is in good agreement with the experimental data, indicating that the numerical results of the model are convergent. The numerical stability of the algorithm is demonstrated by analyzing the simulation results of SAP internally cured concrete. The results reveal that for concrete specimens with different SAP addition ratios, reliable predictions were obtained despite the variations in the simulation input parameters.

Figure 8.

Model-obtained macroscopic stress-strain responses of SAP internally cured concrete.

Conclusively, the comparisons between model predictions and published experimental results indicate that the proposed multiscale micromechanical progressive damage model can reasonably estimate the macroscopic mechanical responses of cementitious composites under uniaxial tensile loading.

Parametric study

This section investigates the effects of SAP addition ratio, aggregate content, initial interfacial crack size, and initial interfacial crack location on the uniaxial tensile strength of SAP-incorporated cementitious composites. In the parametric study, when examining one of these parameters, the others will remain fixed. The mechanical parameters employed in this study are listed in Table 1.

Table 1.

Mechanical parameters employed in this study.

Mechanical parameters	Value
Elastic modulus of coarse aggregate	70 GPa
Poisson’s ratio of coarse aggregate	0.23
Elastic modulus of fine aggregate	80 GPa
Poisson’s ratio of fine aggregate	0.22
Tensile strength of interface	3.0 MPa
Tensile strength of mortar	6.0 MPa
Tensile strength of cement paste	5.0 MPa
Coarse aggregate diameter	10 mm
Fine aggregate diameter	1.0 mm

Effect of SAP addition ratio

Experimental studies have demonstrated that the tensile strength of SAP-incorporated cementitious composites is affected by the ratio of SAP addition (Lam, 2005; Xie et al., 2020; Yao et al., 2012). The effect of various SAP addition ratios on the tensile strength of SAP-incorporated mortar and concrete is revealed in Figure 9. It can be seen from the figure that with a higher SAP addition ratio, the uniaxial tensile strength of the composites is lower. This is because for SAP-incorporated cementitious composites, the elastic properties of the matrix will decrease with the increase of SAP addition. Further, the elastic properties of matrix and aggregate will affect the mechanical behavior of the interfacial microcracks and further affect the overall mechanical properties of composites. With the lower matrix elastic modulus, the interfacial critical strength for the microcracks to start propagating will be lower. Once the crack deviates into the matrix, the lower elastic modulus will result in the crack expanding more easily in the matrix. Finally, this leads to a lower uniaxial tensile strength of composites.

Figure 9.

Effect of SAP addition ratio on the uniaxial tensile strength of SAP-incorporated cementitious composites.

Effect of aggregate content

Aggregate content is an essential parameter that affects the macroscopic mechanical properties of cementitious materials. In this study, four varying aggregate contents were investigated and their influences on the tensile strength of SAP-incorporated cementitious composites were evaluated. The results are rendered in Figure 10. It can be noticed that although the elastic moduli of aggregates are higher than those of the equivalent matrix, and increasing their volume fractions may improve the overall elastic properties of the composites, it does not directly increase the macroscopic tensile strength of the composites. As the aggregate content increases, the initial microcracks at the matrix-aggregate interface also increase during the manufacturing process of cementitious composites, leading to more distributed initial damages. Moreover, with increased aggregate volume fraction, the spaces between the aggregates become more compact, resulting in fewer crack propagation paths in the matrix. These reasons lead to an increase in the brittleness of composites, a decrease in the critical stress for crack propagation, and ultimately a reduction in uniaxial tensile strength.

Figure 10.

Effect of aggregate content on the uniaxial tensile strength of SAP-incorporated cementitious composites.

Effect of initial interfacial crack size

Cementitious composites are prone to initial damage or defects during production due to their intrinsic characteristics. When subjected to external loading, the initial damage or defects will affect the macroscopic mechanical response of the composites, which should be explored in the study. Figure 11 demonstrates the effects of various initial interfacial crack sizes on the tensile strength of SAP-incorporated cementitious composites, where the size of the initial interfacial crack is represented by the angle $\alpha$ . It is observed that the smaller the size of the initial crack, the higher the tensile strength of the composite. As the size of the initial crack grows, the initial defect or damage within the composites becomes more significant, and the damage-induced strain and compliance also increase. Moreover, as the degree of damage increases, the critical stress required for crack propagation is also reduced. Therefore, the uniaxial tensile strength of the composite decreases as the initial interfacial crack size increases.

Figure 11.

Effect of initial interfacial crack size on the uniaxial tensile strength of SAP-incorporated cementitious composites.

Effect of initial interfacial crack location

For cementitious composites, aggregates are randomly and uniformly dispersed inside the matrix, and the initial cracks at the aggregate-matrix interface are also randomly distributed with the aggregates. The location of the initial cracks that appear at the interface affects the stress state of the cracks, which in turn influences the macroscopic mechanical behavior of the cementitious composites. Hence, it is necessary to examine the effects of the initial interfacial crack location on the macroscopic mechanical properties of cementitious composites. The influence of the initial interfacial crack location on the tensile strength of SAP-incorporated cementitious composites was investigated by changing the angle $\varphi$ between the applied tension and the crack normal direction, and the results are rendered in Figure 12. It can be found that the initial interfacial crack location has a significant effect on the tensile strength of the composite. When the interfacial crack is perpendicular to the applied tensile loading, the crack is the easiest to propagate, and the tensile strength of the composite is the lowest. The reason is that when the angle between the applied tensile loading and the crack normal direction is larger, the local stress to which the crack propagation is subjected is weaker, which means that a larger uniaxial tension is needed to make the crack extend further.

Figure 12.

Effect of initial interfacial crack location on the uniaxial tensile strength of SAP-incorporated cementitious composites.

Conclusions

In this paper, a multiscale micromechanics-based progressive damage model is proposed for investigating the overall macroscopic mechanical behavior of cementitious composites featuring SAP under uniaxial tensile loading. A micromechanical framework is first presented for the homogenization process between different scales and for predicting the equivalent elastic properties of the composites. Then, a progressive elastic-damage model is developed to describe the propagation evolutions of interfacial microcracks at various stress stages. In this model, the damage-induced strain and compliance are determined based on micromechanics, the crack paths are identified according to the interface contact stress and matrix cleavage stress, and the crack propagation in the matrix is governed by the stress intensity factors. By comparing with the published experimental results, the model predictions show good agreement with the experimental data, which verifies the reliability and applicability of the proposed damage model in predicting the macroscopic tensile mechanical response of cementitious composites. Finally, the effects of SAP addition ratio, aggregate content, initial interfacial crack size, and initial interfacial crack location on the uniaxial tensile strength of cementitious composites featuring SAP are quantitatively investigated based on the proposed model. According to the study, the following main conclusions are drawn:

The proposed multiscale micromechanical progressive elastic-damage model systematically integrates the homogenization procedure of the macroscopic stress-strain relationship and effective elastic properties with the progressive damage processes (debonding, cracking, propagation) of composites. It is feasible to utilize the proposed model to predict the macroscopic mechanical response of cementitious composites featuring SAP from micro/mesoscopic scales. This provides an effective method for simulating and optimizing high-performance and multifunctional cementitious materials.

The effective elastic properties, overall progressive damage process, and initial microcrack parameters of composites are all accounted for in the model. In simulations, the required parameters/properties can be determined from the relevant material parameters of the micro/mesoscale constituents of the composites without making assumptions about the essential properties.

Quantitative investigations of cementitious composites featuring SAP reveal that the model can analyze the effects of different material parameters on uniaxial tensile strength. This model can also be extended to other particle-reinforced composites by appropriate adaptations and modifications depending on underlying material and damage characteristics.

Issues related to the initial damage of SAP-induced voids, the effect of SAP particles on the healing process, and the brittle behavior of cementitious materials warrant further studies to extend the proposed model and will be investigated in our forthcoming papers.