Stochastic Extension of NMMD Model for Fracture Behaviors of Concrete Materials

1. Introduction

Up to now, concrete is the most widely used composite material in civil engineering. The two remarkable features of concrete material, i.e., non-linearity and randomness (Li and Chen 2009 [1]; Li et al. 2014 [2]), make it still a challenging topic in safety analysis and reliability design of architecture structures.

The non-linearity of concrete material mainly originates from the microcrack nucleation and the macroscopic crack growth. Ever since the groundbreaking work of Griffith (1921 [3]), fracture mechanics and damage mechanics have become the two mainstream approaches in dealing with crack evolution. In the framework of fracture mechanics (Anderson 2005 [4]), the hypothetic pre-existing crack as a geometric entity is extracted from the cracked solid body, and then the eccentric equilibrium equation of crack tip becomes the primary task for crack propagation. Along this line, many important concepts of driving force attached to the crack tip, such as energy release rate (Irwin 1956 [5]), stress intensity factor (Irwin 1957 [6]), J-integral (Rice 1968 [7]) and crack-tip opening displacement (Wells 1963 [8]), have been successively proposed. Correspondingly, the accompanied resistances, e.g., fracture energy, fracture toughness and the critical crack-tip opening displacement, regarded as the constant material properties for a crack to resist fracture, are left to the experimental measurement. However, the ongoing experimental observations reveal that a realistic crack would arise at the location far from the pre-existing crack tip (Haeri et al. 2014 [9]; Cao et al. 2023 [10]; etc.) and the aforementioned fracture parameters (especially the fracture energy) depend on the crack size (instead of being the constant properties) for a broad category of engineering materials including concrete (Malvar and Warren 1988 [11]; Refai and Swartz 1988 [12]; Gao et al. 2022 [13]; etc.), which challenge the two fundamental hypotheses in fracture mechanics. On the other hand, damage mechanics (Lemaitre 2012 [14]), without introducing the hypothesis of pre-existing crack, has the potential advantages in overcoming the drawbacks of fracture mechanics. Although the solid mathematical foundation has been already established in the earlier days (Simó and Ju 1987 [15]; Ju 1989 [16]), traditional local damage models remain far from being popular in simulations of cracking problems due to the lack of appropriate crack representation. Such an embarrassed situation sustains a relatively long period of time perhaps until the appearances of peridynamics (Silling 2000 [17]; Silling et al. 2007 [18]; Sheikhbahaei et al. 2023 [19]) and phase field model (Bourdin et al. 2000 [20]; Borden et al. 2012 [21]; Wu 2017 [22]; Feng and Li 2023 [23]).

With taking advantages of both peridynamics and phase field model, the nonlocal macro-meso-scale consistent damage (NMMD) model (Lu and Chen 2020 [24]; Chen et al. 2021 [25]) opens up a new avenue in simulation of crack evolution from nucleation to propagation. In NMMD model, meso-structure of material point-pairs adopted from peridynamics (Silling 2000 [17]; Silling et al. 2007 [18]) is employed to represent crack topology in geometry. By advocating the geometry-energy translational mechanism, which is linked with an energetic degradation factor also existing in phase field model (Bourdin et al. 2000 [20]; Wu 2017 [22]), the governing equations of NMMD model are constructed on the basis of continuum damage mechanics (Li et al. 2014 [2]). Theoretically speaking, the NMMD model can be deemed as the communication bridge connecting peridynamics and phase field model. Compared to the two popular approaches, nevertheless, numerical evidences (Ren et al. 2023 [26]; Chen et al. 2024 [27]) have demonstrated that the NMMD model owns the highest computational efficiency, for instance, the CPU time consuming of NMMD model is merely one-tenth of phase field model and only one-thirtieth of peridynamics for the same problems with considerable degrees of freedom and similar hardware. Besides, issues of the fixed Poisson’s ratio in bond-based scheme (Silling 2000 [17]), the zero-energy modes of state-based formulation (Silling et al. 2007 [18]) and the mesh-size (or particle distance) sensitivity of strain-softening regime (Sheikhbahaei et al. 2023 [19]) suffered by peridynamics due to the reformulation of differential equations, can be overcome spontaneously by NMMD model (Lu and Chen 2020 [24]; Ren et al. 2024 [28]). Benefiting from the intrinsic macro-meso-scale consistent mechanism, the NMMD model can capture the experimentally-observed non-constant fracture parameters (Lu et al. 2024 [29]), which cannot be predicted by phase field model with logical consistency. In recent years, some representative applications of NMMD model include various quasi-static fractures with different mechanisms (Ren et al. 2023 [26]; Xia et al. 2024 [30]; Ren et al. 2024 [31]), dynamic crack branching with and without rate-dependent effect (Lu and Chen 2021 [32]; Zhao et al. 2024 [33]) and even thermo-mechanical coupling cracking problems (Lv et al. 2024 [34]). However, almost all of these works focus on the deterministic perspective with neglecting the ubiquitous randomness (or stochasticity) existing in engineering materials and structures.

For heterogeneous materials like concrete, randomness plays an important role in the mechanical behaviors (Li and Chen 2009 [1]; Li et al. 2014 [2]), including stochastic fluctuations of bearing capacity and fracture toughness, random walk of crack path with uncertainty transition of fracture modes. To this end, this paper focuses on material randomness, and extends the deterministic NMMD model to the stochastic version equipped with stochastic harmonic function method of second kind. The effects of material randomness in the stochastic NMMD model are investigated in behaviors of load-displacement responses of the entire cracking process, crack patterns and energy dissipations (especially fracture energy).

The rest sections are arranged as follows. In Section 2, the basic idea and numerical implementation of deterministic NMMD model are briefly introduced. On that basis, the stochastic extension of NMMD model for concrete materials are depicted in Section 3 with sufficient details. Numerical examples of three point bending beams without and with various sizes of initial crack, concerning the influence of random fields with different correlation lengths and the stochastic fracture failure analysis by using the stochastic NMMD model, are discussed in Section 4. Conclusions and problems to be studied in the future are given in Section 5.

2. Deterministic NMMD Model

Let $B\subset {ℝ}^d$ ($d=1,2,3$) be the reference configuration of a solid body, and $\partial B$ denote the external boundary with the outward normal vector $n$. The density of body force vector $\overline{b}$ is distributed in the interior of the solid body, and on the exterior, the surface traction vector $\overline{t}$ and the prescribed displacement vector $\overline{u}$ are respectively applied on Neumann part $\partial B_t$ and Dirichlet part $\partial B_u=\partial B\backslash \partial B_t$. The resulting displacement vector and effective stress tensor of any material point $x\in B$ are denoted by $u(x)$ and ${\sigma }_0(x)$, correspondingly.

The governing equation of NMMD model (Lu and Chen 2020 [24]; Chen et al. 2021 [25]) are established on the framework of continuum damage mechanics (Li et al. 2014 [2]) in the manner of articulation, accuracy and effectiveness, i.e.,

$$\mathrm{div}[(1-D_{\mathcal{l}}\circ {\Omega }_{\mathcal{l}}){\sigma }_0]+\overline{b}=0\mathrm{in}B$$

(1)

with the following Neumann and Dirichlet conditions

$$\left\{\begin{array}{ll}(1-D_{\mathcal{l}}\circ {\Omega }_{\mathcal{l}}){\sigma }_0\cdot n=\overline{t}& \mathrm{on}\partial B_t\\ u=\overline{u}& \mathrm{on}\partial B_u\end{array}\right.$$

(2)

where $\mathrm{div}(\cdot )$ is the divergence operator, and $D_{\mathcal{l}}\circ {\Omega }_{\mathcal{l}}$ denotes the physical translation from geometry damage ${\Omega }_{\mathcal{l}}$ to the energy one $D_{\mathcal{l}}$.

In NMMD model, the geometry damage ${\Omega }_{\mathcal{l}}$ at macroscale is expressed by the ratio of damaged mesoscale point-pairs to all the point-pairs within an influence domain (see Figure 1 for the intuitive illustration) with a radius of $\mathcal{l}>0$, similar to peridynamics (Silling 2000 [17]; Silling et al. 2007 [18]), that is

$${\Omega }_{\mathcal{l}}(x)={\displaystyle \frac{1}{{\displaystyle {\int }_{H_{\mathcal{l}}(x)}\varphi \mathrm{d}V(x^{\prime })}}}{\displaystyle {\int }_{H_{\mathcal{l}}(x)}\omega (x,x^{\prime })\varphi \mathrm{d}V(x^{\prime })},\forall x\in B$$

(3)

where the influence domain $H_{\mathcal{l}}(x)=\{x^{\prime }\in {ℝ}^d:0\text{<}‖x^{\prime }-x‖\le \mathcal{l}\}\cap B$ is actually the $d$-dimensional ball intersecting with the reference configuration, $\varphi$ is the influence function with the default convention $\varphi =1$ for isotropic materials, and $\omega \in [0,1]$ is the mesoscale damage of point-pair $(x,x^{\prime })$. It should be noted that the point set shown in Figure 1 is the simplest case of the $n$-point meso-structures (Chen et al. 2024 [27]).

In dealing with the evolution of mesoscale damage, the natural deformation (also called predominant elongation) is employed to depict the separation of a material point-pair, namely

$$\lambda (x,x^{\prime })=[u(x^{\prime })-u(x)]\cdot (x^{\prime }-x)/‖x^{\prime }-x‖$$

(4)

and hence the rate-independent evolution of mesoscale damage can be expressed as

$$\omega (x,x^{\prime })=\varpi \circ \kappa (x,x^{\prime }),\kappa =\underset{\tau \in [0,t]}{\max }[\lambda (x,x^{\prime },\tau ),{\lambda }_{\mathrm{c}}]$$

(5)

where: $\kappa$ denotes the maximum deformation of a material point-pair during the possible loading and unloading history; ${\lambda }_{\mathrm{c}}={\sigma }_{\mathrm{t}}\mathcal{l}/E$is the critical threshold relating to tensile strength ${\sigma }_{\mathrm{t}}$ and Young’s modulus $E$ of materials; $\varpi$ is the evolution law of mesoscale damage that can be determined from the Ab-initio molecular dynamics simulation, but herein the specified formulation is adopted (Lu and Chen 2020 [24])

$$\varpi \circ \kappa =1-\exp (-\gamma \kappa )$$

(6)

in which $\gamma$ is left to be calibrated by some kind of standard tests of material.

Up to now, there are two approaches to establish the relationship between damage variables in geometric and in mechanical (or energetic) senses: the macroscopic (Lu and Chen 2020 [24]) and the mesoscopic (Ren et al. 2024 [28]). In this paper, the macroscopic approach is preferred due to the simplicity, and it reads

$$D_{\mathcal{l}}\circ {\Omega }_{\mathcal{l}}+g({\Omega }_{\mathcal{l}})=1$$

(7)

where $g\in [0,1]$ is called energetic degradation factor, which is utilized from phase field model (Bourdin et al. 2000 [20]; Wu 2017 [22]). It is physically identified that the energetic degradation factor must be a nonincreasing and convex function of geometry damage (Chen et al. 2021 [25]), which reveals the nonlinear translation mechanism from geometry damage ${\Omega }_{\mathcal{l}}$ to energy damage $D_{\mathcal{l}}$ that the irreversible breakage of mesoscale point-pairs in geometry leads to the non-commensurate energetic degradation. For quasi-brittle materials like concrete, the energetic degradation factor can be empirically taken as (Lu et al. 2024 [29])

$$g({\Omega }_{\mathcal{l}})={\displaystyle \frac{1}{1+q{\Omega }_{\mathcal{l}}}}{(1-{\Omega }_{\mathcal{l}})}^p,p>1,q\ge 0$$

(8)

where $p$ and $q$ are the two degradation parameters, which can be physically identified by the standard test specimens.

Once the formulation of effective stress tensor ${\sigma }_0$ is determined for the specified material, the solutions of NMMD model (i.e., the displacement vector $u$ and the geometry damage ${\Omega }_{\mathcal{l}}$) can be directly obtained by the following weak form

$${\displaystyle {\int }_B{\nabla }^{\mathrm{S}}w\text{:}g({\Omega }_{\mathcal{l}})}{\sigma }_0\mathrm{d}V(x)-{\displaystyle {\int }_{B}w\cdot \overline{b}\mathrm{d}V(x)}-{\displaystyle {\int }_{\partial B_t}w\cdot \overline{t}dS(x)}=0$$

(9)

with the discretization of finite element method, where ${\nabla }^{\mathrm{S}}w$ denotes the symmetric gradient of test function $w$. It should be noted that the arc-length method (e.g., May and Duan 1997 [36]) is necessary for capturing the post-peak behavior of strain-softening. In the solution process, cracks represented by the geometry damage ${\Omega }_{\mathcal{l}}$ naturally form as a result of deformation without any extra crack tracking strategy.

3. Stochastic Extension of NMMD Model for Concrete Materials

In this section, the deterministic NMMD model described in Section 2 is extended by incorporating random spatial fluctuations in the material properties to describe stochastic fracture behavior in heterogeneous materials like concrete. Due to incomplete observation and imperfect control (Li and Chen 2009 [1]) during the manufacturing process, the constituents (e.g., aggregate configurations and locations, inhomogeneous distribution of water and holes with various shapes), in addition to non-uniform drying and hardening, make the physical parameters of concrete materials hard to be measured precisely in practice. The method of packing aggregates and generating interface layers (Yang et al. 2019 [37]; Xia et al. 2021 [38]; etc.), though proved to be feasible to depict the concrete heterogeneity, suffers expensive computational costs and even sometimes un-convergence problems, excepting the complicated pre-processing. In this paper, the characterization of uncertain fluctuations is therefore preferred to be included in the formulation of random fields.

For concrete materials considered herein, the infinitesimal strain tensor $\epsilon ={\nabla }^{\mathrm{S}}u$ is sufficient to measure the resulting deformation intensity, and hence the effective stress tensor in Eq. (9) is expressed as

$${\mathbf{\sigma }}_0(\mathbf{x})=\mathbf{E}(\mathbf{x})\text{:}\epsilon (\mathbf{x}),\forall \mathbf{x}\in B$$

(10)

where $\mathbf{E}$ is the fourth-order elastic tensor with the essential concentration of material properties. If isotropic property is further considered, we have

$$\mathbf{E}=E\left[{\displaystyle \frac{1}{1+\nu }}{(\mathbf{I}\otimes \mathbf{I})}^{\mathrm{S}}+{\displaystyle \frac{\nu }{(1+\nu )(1-2\nu )}}\mathbf{I}\odot \mathbf{I}\right]$$

(11)

which involves Young’s modulus $E$ and Poisson's ratio $\nu$, where ${(\mathbf{I}\otimes \mathbf{I})}^{\mathrm{S}}$ and $\mathbf{I}\odot \mathbf{I}$ are respectively the super-symmetric identity tensor and the trace projection tensor, or alternatively in the index form, i.e., ${[{(\mathbf{I}\otimes \mathbf{I})}^{\mathrm{S}}]}_{\mathbf{I}jkl}=({\delta }_{\mathbf{I}k}{\delta }_{jl}+{\delta }_{\mathbf{I}l}{\delta }_{jk})/2$ and$[\mathbf{I}\odot \mathbf{I}{]}_{\mathbf{I}jkl}={\delta }_{\mathbf{I}j}{\delta }_{kl}$, by the representation of Kronecker delta function ${\delta }_{\mathbf{I}j}$.

In the framework of stochastic representation and computer simulation, the uncertain fluctuation of material parameter of at each material point is assigned to a non-zero Gaussian random field, that is

$$Y(\mathit{x})=\mu +Y_0(\mathit{x})$$

(12)

where $\mu >0$ is the mean value, and $Y_0(\mathit{x})$ is the zero-mean Gaussian random field sharing the same standard deviation with $Y(\mathit{x})$.

To generate the zero-mean Gaussian random field $Y_0(\mathit{x})$, Karhunen-Loève expansion (KL) method (Kosambi 1943 [39]) and the classical spectral representation method (Shinozuka and Deodatis 1991 [40]; 1996 [41]) are available choices. However, these two common methods work adequately for low-dimensional problems. The challenge of large computational cost is encountered by both methods when simulating high-dimensional random fields (Chen et al. 2018 [42]). To solve this issue, the concept of stochastic harmonic function (SHF) was proposed by Chen and Li (2011 [43]). Similar to the classical spectral representation method, random fields are still represented as the combination of a series of harmonic components in SHF method. The difference is that, in the SHF method, the wave numbers of the harmonic components are also random, not just the phase angles. Precisely because of this, the target wave-number power spectral density (PSD) function can be reproduced by SHF method with only a limited number of harmonic components. According to the probability distribution of the wave number within each wave-number subdomain, this method can be further divided into the SHF of the first kind (SHF-I; Chen and Li 2011 [43]) and of the second kind (SHF-II; Sun et al. 2011 [44]).

Herein the SHF-II method endowed with tensor products expansion (Liang et al. 2012 [45]; Chen et al., 2018 [42]) is adopted, and for two-dimensional cases, i.e., $\mathit{x}=(x_1,x_2)$, it reads

$$\begin{array}{c}{Y}_0^{\mathrm{SHF}-\mathrm{II}}\left(\mathit{x}\right)={\displaystyle \sum _{\mathbf{I}=1}^{N_1}{\displaystyle \sum _{j=1}^{N_2}A\left(k_{1,\mathbf{I}},k_{2,j}\right)\left[\cos \left(k_{1,\mathbf{I}}{x}_1+k_{2,j}{x}_2+{\phi }_{\mathbf{I}j}^{(1)}\right)\right.}}\\ \left.+\cos \left(-k_{1,\mathbf{I}}{x}_1+k_{2,j}{x}_2+{\phi }_{\mathbf{I}j}^{(2)}\right)\right]\end{array}$$

(13)

where ${\phi }_{\mathbf{I}j}$ is the independently random phase distributed uniformly over $\left[0,2\pi \right]$, $k_{1,\mathbf{I}}$ and $k_{2,j}$ are the random wave numbers required to be positive, $N_1,N_2\in {\text{ℕ}}^{+}$ represent the two truncation numbers. According to the relationship between autocorrelation function and PSD function, namely the Wiener-Khinchin formula (Li and Chen 2009 [1]), the amplitude $A\left(k_{1,\mathbf{I}},k_{2,j}\right)$ in Eq. (13) can be given by

$$A\left(k_{1,\mathbf{I}},k_{2,j}\right)=\sqrt{4S_Y\left(k_{1,\mathbf{I}},k_{2,j}\right)\left(k_{1,\mathbf{I}}^U-k_{1,\mathbf{I}}^L\right)\left(k_{2,j}^U-k_{2,j}^L\right)}$$

(14)

where $S_Y\left(k_1,k_2\right)$ represents the target PSD function of random field $Y(x)$ in the first-quadrant.

For the stationary Gaussian random field usually employed to characterize the condensed state of material properties, the form of $S_Y\left(k_1,k_2\right)$ is taken as (Vanmarcke 2010 [46])

$$S_Y\left(k_{1,\mathbf{I}},k_{2,j}\right)={\displaystyle \frac{{\zeta }^2{a}^2}{4\pi }}\exp \left[-{\displaystyle \frac{a^2}{4}}\left(k_{1,\mathbf{I}}^2+k_{2,j}^2\right)\right]$$

(15)

where: $\zeta$ is the standard deviation of the zero-mean random field $Y\left(\mathit{x}\right)$, and can be determined by the mean value $\mu$ and the coefficient of variation $\beta$, i.e., $\zeta =\mu \beta$; the scale parameter $a$ can be expressed by $a=\chi /\sqrt{\pi }$, and $\chi$ is the correlation length of $Y\left(\mathit{x}\right)$.

In practical implementation of SHF-II method, once a material property is specified to the non-zero Gaussian random field with positive wave-number domain $K=\left[k_1^{\mathrm{L}},k_1^{\mathrm{U}}\right]\times \left[k_2^{\mathrm{L}},k_2^{\mathrm{U}}\right]$, where $k_1^{\mathrm{L}},k_2^{\mathrm{L}}\in {\text{ℝ}}^{+}$ and $k_1^{\mathrm{U}},k_2^{\mathrm{U}}\in {\text{ℝ}}^{+}$ are the lower and upper bounds of the wave numbers, the wave-number domain $K$ can be discretized a series of small regions, i.e., $K={\displaystyle {\cup }_{\mathbf{I}=1}^{N_1}{\displaystyle {\cup }_{j=1}^{N_2}{K}_{\mathbf{I}j}}}$, and $K_{\mathbf{I}j}=\left[k_{1,\mathbf{I}}^{\mathrm{L}},k_{1,\mathbf{I}}^{\mathrm{U}}\right]\times \left[k_{2,j}^{\mathrm{L}},k_{2,j}^{\mathrm{U}}\right]$. And then, in each small wave-number region $\left[k_{1,\mathbf{I}}^{\mathrm{L}},k_{1,\mathbf{I}}^{\mathrm{U}}\right]\times \left[k_{2,j}^{\mathrm{L}},k_{2,j}^{\mathrm{U}}\right]$, $k_{1,\mathbf{I}}$ and $k_{2,j}$ in Eq. (13) can be deemed as uniform distributions.

4. Numerical Examples

Let us consider the three-point bending concrete beam (TPBB) tests conducted by Grégoire et al. (2013 [47]), the geometry and boundary conditions are shown in Figure 2, in which the thickness $b$ and height $H$ are $50\mathrm{mm}$ and $400\mathrm{mm}$, respectively. Each beam is simply supported at the bottom edge, and the vertical force $P$ is applied at the midpoint of the top edge of the beam. During the loading process, the force $P$ at loading point, the displacement of loading point (DLP) and the crack mouth opening displacement (CMOD) are recorded. The preliminary analyses of these beams assumed to plane stress problems were performed previously by Lu et al. (2024 [29]) using the deterministic NMMD model, and the mechanical parameters (including Young’s modulus $E$, Poisson’s ratio $\nu$, radius of influence domain $\mathcal{l}$, critical elongation ${\lambda }_{\mathrm{c}}$ and the two degradation parameters $p$ and $q$) listed in Table 1 were determined according to the experiments (Grégoire et al. 2013 [47]).

In this section, consistent with Lu et al. (2024 [29]), different notch-to-depth ratios, i.e., $\xi =\left\{0,0.02,0.05,0.1,0.2,0.5,0.7,0.9\right\}$ with the notch width of $10\mathrm{mm}$, are considered for investigating stochastic fluctuations of bearing capacity, fracture energy and the random walk of crack path with uncertainty transition of fracture modes by using the proposed stochastic NMMD model. Based on the intuition sharpened by hundreds of numerical examples of deterministic NMMD model, only the spatial variability of Young's modulus is considered in the following analyses.

4.1. Influence of Random Fields with Different Correlation Lengths

Without loss of generality, the configuration of TPBB is might as well divided into three parts according to material properties, as shown in Figure 3. The Young’s modulus in the region of initial notch is set to zero, and the random filed of Young’s modulus is generated only in the region where the crack is expected to propagate. Apart from the two cases, the mean value of Young’s modulus, i.e., $E=3.7\times {10}^4\mathrm{MPa}$, is adopted in the remaining regions. Correspondingly, the mesh is refined only within the range of random field and a coarse mesh is used in the remaining region. Additionally, the material points are only placed at the Gaussian integration points of elements within the mesh refined region. The aforementioned mesh discretization strategy and arrangement of material points are solely aimed at reducing unnecessary computations.

Shown in Figure 4 is the comparison of random field samples with three different correlation lengths, i.e., $\chi =\left\{1.5\times {10}^{-3}\mathrm{m},1.5\times {10}^{-2}\mathrm{m},1.5\times {10}^{-1}\mathrm{m}\right\}$. In the three cases, the coefficient of variation $\beta$ is $0.15$ and the truncation number of the wave number in each dimension of wave-number domain is taken to be $25$, i.e., $N_1=N_2=25$. From Figure 4 (though only displayed ten samples in each case), it is clearly found that the larger the correlation length, the more pronounced homogeneity of random field. Visually speaking, these samples of Young’s modulus are very similar to the realistic concrete materials containing aggregates and holes with various shapes. In other words, the method of packing aggregates and generating interface layers (Yang et al. 2019 [37]; Xia et al. 2021 [38]) can be alternatively replaced by the formulation of random field in the numerical generation of concrete materials.

After computation of stochastic NMMD model, the load-CMOD curves and load-DLP curves for the case $\xi =0.5$, which are obtained from the experiments (Grégoire et al. 2013 [47]), the deterministic NMMD model (Lu et al. 2024 [29]) and the proposed stochastic NMMD model with numerical samples concerning different correlation lengths, are depicted in Figure 5 and Figure 6. It can be concluded from Figure 5 and Figure 6 that, as the correlation length decreases, the dispersion of the curves gradually decreases, which is consistent with the conclusions reached when considering the spatial variability of the critical elongation (Ren et al. 2021 [48]). Especially, when the correlation length $\chi$ is $1.5\times {10}^{-3}\mathrm{m}$, the results of the numerical samples closely match the results of the deterministic analysis. The reason is that the topologic damage defined in Eq. (3) itself represents the mean of the mesoscale damage within the influence domain.

When the correlation length becomes as small as the mesh size, the Young’s modulus in the random field tends to become a series of independent and identically distributed random variables. Since the Young’s modulus characterizes the ability of material to resist deformation, this also means, from a qualitative perspective, that the deformations of point pairs are independent and identically distributed random variables. According to the law of large numbers, as the number of independent and identically distributed random variables approaches infinity, the variance of their mean tends to zero. Therefore, the topologic damage calculated from the numerical samples is close to the result of the deterministic analysis, leading to the reduction in the variability of the curves. On the other hand, for the appropriate choice of correlation length, the fluctuation of load-displacement responses in experimental observation (Grégoire et al. 2013 [47]) can be completely covered by the stochastic NMMD model.

The crack patterns calculated by numerical samples with different correlation lengths are shown in Figure 7. From Figure 7 it can be observed that when the correlation length $\chi$ of the random field is $1.5\times {10}^{-2}\mathrm{m}$, the zig-zag crack patterns will emerge with the instantaneous transition from mode-I to mixed-mode, which is very closed to the experimental results (Yin et al. 2023 [51]). Nevertheless, in the other two cases, the cracks are straighter. The reason for the occurrence of the above phenomenon with small correlation length has already been explained when investigating the variation in the dispersion of the curves aforementioned.

Based on the above analysis of numerical results, it can be concluded that when the correlation length is comparable to the size of influence domain, the range of load-CMOD curves of numerical samples matches well with the experimental results, and the crack patterns is also consistent with physical phenomena.

4.2. Stochastic Fracture Failure Analysis

In this subsection, the cases with the correlation length $\chi =1.5\times {10}^{-2}\mathrm{m}$ are investigated. Shown in Figure 8 and Figure 9 are the load-CMOD curves and load-DLP curves of the experimental results (Refai and Swartz 1988 [12]), deterministic results (Lu et al. 2024 [29]) and numerical samples obtained by the proposed stochastic NMMD model. From Figure 8, it can be seen that the load-CMOD curves of numerical samples with $\xi =\left\{0,0.2,0.5\right\}$ fall within the range of experimental results, indicating that the nonlinear characteristics of the fracture failure process can be captured by NMMD model. It is also of interest to see that in Figure 9 as the notch-to-depth $\xi$ varies from zero to $0.9$, the characteristics of snap-back gradually decrease.

Crack patterns of specimens with different notch-to-depth are plotted in Figure 10. As expected and consistent with the experimental phenomena, the crack still initiates from the tip of the initial notch (or from the midspan in the case of $\xi =0$), propagates tortuously in the vertical direction, and eventually reaches the upper edge of the TPBB. Thus, the numerical results of NMMD model can be more closely align with physical reality by introducing the spatial variability of material properties. That is to say, the stochastic NMMD model with single random field (namely Young’s modulus) can pretty well capture not only the stochastic fluctuations of bearing capacity but also the random walk of crack path. This discovery shakes up the accustomed cognition (Le et al. 2016 [52]; Gorgogianni et al. 2022 [53]) that stochastic simulations of quasi-brittle fracture should contain at least two mechanical properties with spatial randomness.

In the published paper (Lu et al. 2024 [29]), it is found that the NMMD model does not require the input of fracture energy $G_f$ but can still yield the conclusion that fracture energy $G_f$ is not a constant but rather depends on the notch-to-depth $\xi$. This phenomenon can also be observed from experiments. However, a slight drawback of this work is that there is only one result for each notch-to-depth case. Therefore, to further strengthen the foundation of the conclusion mentioned before, the fracture energy $G_f$ obtained from stochastic analysis is also calculated here. For the TPBB tests, the fracture energy $G_f$ can be expressed as

$$G_f={\displaystyle \frac{1}{b\left(\text{ℒ}-\xi H\right)}}{\text{ℰ}}_f\left(u_{\mathrm{cut}}\right)$$

(16)

where $\text{ℒ}$ is the total length of the crack that can be measured by the topologic damage in Eq. (3), and ${\text{ℰ}}_f\left(u_{\mathrm{cut}}\right)$ represents the energy dissipation at the cutting displacement $u_{\mathrm{cut}}$ of loading point. As depicted in Figure 11, ${\text{ℰ}}_f\left(u_{\mathrm{cut}}\right)$ can be determined by

$${\text{ℰ}}_f\left(u_{\mathrm{cut}}\right)={\displaystyle {\int }_0^{u_{\mathrm{cut}}}P\left(s\right)\mathrm{d}s}-{\displaystyle \frac{1}{2}}{P}_f{u}_{\mathrm{cut}}$$

(17)

The comparison of fracture energy of the experiment (Refai and Swartz 1988 [12]), deterministic analysis (Lu et al. 2024 [29]) and stochastic simulation is plotted in Figure 12 with the cutting displacement tends to infinity, i.e., $u_{\mathrm{cut}}\to \infty$. It can be seen that the experimental results can be enveloped by the numerical sample results. Moreover, when the notch-to-depth $\xi$ is less than $0.5$, it is of interest to see that the mean value of numerical samples is larger than that of deterministic results (Lu et al. 2024 [29]). The reason is that fracture energy represents the material's ability to resist fracture from the physical perspective, i.e., the greater the fracture energy, the more energy is dissipated as the crack propagates within the material.

Recalling the stochastic results shown in Figure 10, the crack paths are tortuous with instantaneous transition of fracture modes, while without considering the spatial variability of the material properties, the crack paths are always straight along the direction the initial crack (i.e., mode-I) as shown in Lu et al. (2024 [29]). As a result, the crack length in the former case is significantly greater than in the latter, leading to higher energy dissipation. In contrast, when the notch-to-depth $\xi$ exceeds $0.5$, the influence of the random field on the propagation of crack decreases, resulting in the mean value of fracture energy becoming closer to the deterministic results. Anyway, the logical inconsistency with physical phenomena in models that use constant fracture energy as an input parameter, such as the crack band model (Jirásek and Bauer 2012 [49]), the gradient-enhanced damage model (Geers et al. 1998 [50]), and the phase-field model (Bourdin et al. 2000 [20]; Wu et al. 2017 [22]), is once again validated by non-constant fracture energy observed experimentally and numerically.

5. Conclusions

In this paper, the deterministic nonlocal macro-meso-scale consistent damage (NMMD) model is extended to analyze stochastic fracture behaviors of typical composite materials like concrete with the equipment of stochastic harmonic function method of second kind. In the proposed stochastic NMMD model, non-linearity caused by crack evolution and randomness resourced from uncertain distribution of concrete materials are unified in the manner of articulation. Numerical examples of three point bending beams without and with various sizes of initial crack, also compared with the experimental observations, demonstrate the effectiveness of stochastic NMMD model. The main conclusions are given as follows:

(1) Though directly inherited the deterministic NMMD model, the proposed stochastic version still manifests different capacities of capturing fracture behaviors of concrete materials as observed in experiments, including stochastic fluctuations of peak load and random walk of crack path with instantaneous transition of fracture modes.

(2) It is widely believed that stochastic simulations of quasi-brittle fracture should contain at least two mechanical properties with spatial randomness. However, the preliminary success of stochastic NMMD model only assigned with a single random field (i.e., Young’s modulus) implies that the previous understanding might be inappropriate.

(3) It is found that in quasi-brittle fracture problems the mean value of fracture energy obtained by the stochastic NMMD model as an output always greater than that obtained by the deterministic NMMD model with the mean value of random field as an input parameter. Such an unexpected result is attributed to the nonlinear energy-dissipating capacity of crack, that is, the extension of zig-zag crack appearing in the stochastic model requires more energy than that of straight crack occurring in the deterministic model.

Although proved to be feasible and effective, the presented stochastic NMMD model still needs further improvements in the future, including the revelation of relationship between the correlation length of the random field and the size of influence domain, the extension of anisotropic randomness and the analysis of dynamic cracking problems with rate-dependent effect.