Three-Dimensional Probabilistic Semi-Explicit Cracking Model for Concrete Structures

Abstract

This paper introduces a three-dimensional (3D) semi-explicit probabilistic numerical model for simulating crack propagation within the framework of the finite element method. The model specifically addresses macrocrack propagation using linear volume elements. The criteria governing the macrocrack propagation is based on the softening behavior observed in concrete under uniaxial tension. This softening behavior corresponds to a dissipated cracking energy that is equal to the mode I critical fracture energy $\left(G_{IC}\right)$ used in the Linear Elastic Fracture Mechanics theory (LEFM). The probabilistic nature of this model revolves around the random distribution of two mechanical properties: tensile strength $\left(f_t\right)$ and fracture energy, which varies based on the volume of finite elements. The scattering of the fracture energy increases as the volume of finite elements decreases in order to consider the strong heterogeneity of the material. This work primarily aims to estimate the relationship between the standard deviation of $G_{IC}$ and the volume of finite elements through the development of the numerical model. For this purpose, an inverse analysis is conducted based on a fracture mechanical test simulation. This test involves macrocrack propagation in a large Double Cantilever Beam (DCB) specimen with a crack length exceeding two meters. The proposed inverse analysis procedure yields highly significant results, indicating that the numerical model effectively evaluates both crack length and crack opening during propagation.

1. Introduction

The accurate modeling of crack initiation and propagation holds significant importance in the numerical analysis of concrete structures, leading to the development of various constitutive models within this domain. Among them, the following approaches can be mentioned: damaged models [1,2,3,4]; the smeared crack model [5,6,7]; plasticity [8]; the cohesive crack model or fictitious crack model [9,10]; and lattice models [11]. These modeling approaches mainly follow deterministic principles. However, there exist also probabilistic models that utilize statistical techniques to directly address material heterogeneity, offering promising alternatives for numerically simulating the nonlinear behavior of concrete structures [12,13,14,15,16]. Among the probabilistic models, two distinct formulations have emerged, differing in their approach to cracking analysis (whether it focuses on crack initiation, localization, or propagation) and their treatment of kinematic discontinuities associated with cracks, namely, explicit and semi-explicit models.

The first formulation [12] adopts a highly localized approach. This method finely models cracking by employing interface elements to simulate, with a lot of precision, all cracks opening, including very small cracks. In contrast, the second model employs linear volume elements to simulate only macrocrack propagation in large concrete structures. The utilization of interface elements in the first approach can lead to significant computation times, limiting the applicability of the explicit method. In contrast, the second model, employing volume elements for macrocrack modeling, ensures reasonable computation times, rendering it more suitable for analyzing large structures.

The concept behind the semi-explicit probabilistic cracking numerical model has been developed in recent years [14,15]. In this model, the macrocrack propagation is governed by two criteria: (a) the macrocrack initiation is associated with the uniaxial tensile strength ($f_t$) and (b) the macrocrack propagation is linked to the mode I critical fracture energy ($G_{IC}$) used in the well-known linear elastic fracture mechanics (LEFM) theory.

The propagation of a macrocrack occurs when all $G_{IC}$ is fully dissipated, leading to the representation of a macrocrack as a succession of fully damaged elements converging to form the macrocrack structure. The originality of the proposed approach lies in its characterization as a type of LEFM model, wherein volume elements are employed instead of singular elements usually used in this type of numerical model. This modification of finite element type results in smoothed stresses concentration at the macrocracks tip.

Therefore, the model presented in this paper is not a classical damage model or smeared crack model [1,2,6,17,18,19,20,21]. In these conventional models, damaged zones are treated as physically realistic entities, and crack openings are deduced from situations involving these zones of damage. On the contrary, in the proposed model, physically existing damaged zones are not considered. The introduction of a damage parameter in the theoretical approach is arbitrary and serves solely as a method to dissipate $G_{IC}$.

The mechanical properties $f_t$ and $G_{IC}$ are considered as random parameters depending on the size of the mesh elements employed. In the case of $f_t$, previous experimental and numerical investigations [15,22] have determined and validated its probabilistic properties in terms of mean and standard deviation. Estimating these values relies on factors such as the size of finite elements (FE) and the maximum aggregate size (further described in Section 2.3). Concerning $G_{IC}$, its mean value represents an intrinsic characteristic of concrete that remains unaffected by the size of the mesh elements. However, its standard deviation exhibits size-dependent variation, as observed for $f_t$. This assumption is a key aspect of the proposed model. Consequently, the challenge lies in obtaining the standard deviation value for $G_{IC}$, with the quantification of this size-dependent variation being the primary objective of the present paper.

Another contribution of the present work is the improvement in the implementation of the Monte Carlo (MC) method. This method is used due to the probabilistic nature of the numerical model, and involves conducting numerous numerical simulations of a structural problem using distinct spatial distributions of mechanical material properties. However, executing multiple finite element analyses leads to significant computational cost. To address this issue, a parallelization strategy is adopted in the implementation of the method, providing a significant decrease in the execution time of three-dimensional simulations. This approach overcomes the computational time constraints associated with the intensive computational cost related to the inverse analysis process, ensuring more efficient and timely results in large-scale simulations.

Given these considerations, this study showcases significant progress in both the scale of application, particularly for large-scale concrete structures, and the methodology for determining the probabilistic properties of the 3D semi-explicit model. To achieve this, a comprehensive explanation of the inverse analysis procedure used to ascertain the standard deviation of $G_{IC}$ based on element size is provided. Additionally, the results of three-dimensional simulations conducted on a large Double Cantilever Beam (DCB) specimen, depicting macrocrack propagation extending up to 2 m, are presented.

2. Development and Implementation of the 3D Model

2.1. Probabilistic Model Formulation

The three-dimensional semi-explicit probabilistic model focuses on crack propagation, specifically addressing failure in mode I due to tensile stresses within concrete structures. This model is developed within the framework of the finite element method (FEM), incorporating considerations of heterogeneity and volume effects through a probabilistic approach. Implemented within a platform continually refined by researchers from the Civil Engineering Department at COPPE/UFRJ, the code is written in FORTRAN language. This platform has served as a foundational tool for numerous studies across various applications over the years [16,23,24].

A summary of the model’s formulation along with its primary features is depicted in Figure 1. This modeling approach utilizes volume elements, specifically tetrahedral solid elements employing linear interpolation functions. To account for material heterogeneity, each finite element represents a volume of heterogeneous material (Figure 1a), quantified by a heterogeneity degree ($r_e$). This degree is evaluated as the ratio of the finite element’s volume ($V_e$) to that of the maximum aggregate ($V_a$) (Figure 1b).

The model’s probabilistic nature is established through the random distribution of tensile strength and fracture energy among the mesh elements. These properties follow Weibull and lognormal distributions, respectively. Notably, the randomness of these mechanical properties directly correlates with the material’s heterogeneity degree. The selection of the Weibull distribution [25,26] for the tensile strength was supported by experimental evidence as demonstrated by Rossi et al. [22]. The rationale behind this choice was rooted in the observation that concrete exhibits a highly heterogeneous nature and displays a brittle behavior under uniaxial tension. Hence, the adoption of the Weibull theory aligns well with these characteristics. On the other hand, the preference for a lognormal distribution to represent fracture energy is somewhat more arbitrary but remains a classical choice in such contexts.

At the elementary level, specifically in each finite element volume, the cracking process induces energy dissipation, governed by an isotropic damage law. The progression of damage is characterized by the $(\sigma \times \epsilon )$ diagram (Figure 1c). This constitutive law accounts for both the tensile strength and volumetric density of dissipated energy ($g_{IC}$). The value of $g_{IC}$ is determined using an energetic regularization technique proposed by Bazant and Oh [27]. It relates to the material fracture energy as follows: $g_{IC}=G_{IC}/l_e$, where $l_e$ denotes the elementary characteristic length evaluated as: $l_e={\left(V_e\right)}^{1/3}$. Therefore, given the random distribution of $f_t$ and $g_{IC}$ over the finite element mesh, each finite element follows distinct laws characterizing its own dissipation process.

The model’s global structural responses are obtained using a Monte Carlo approach (Figure 1d). This methodology involves conducting numerous simulations of the structural problem, each employing different spatial distributions of mechanical material properties, to achieve statistically significant results.

The initiation of the elementary dissipative process occurs when the maximum principal stress (${\sigma }_1$) at a specific Gauss point reaches the randomly assigned material tensile strength ($f_t$). Subsequently, when the total available energy for the finite element is dissipated, the element is considered cracked, and its elementary stiffness matrix is then nullified  [14]. This dissipation mechanism is mathematically expressed through an isotropic damage law, as depicted in Equation (1), where ${\tilde{\epsilon }}_0$, ${\tilde{\epsilon }}_{fi}$, and ${\tilde{\epsilon }}^k$, respectively, denote the damage initialization strain, the maximum critical strain, and the equivalent strain.

$$D=1-\frac{{\tilde{\epsilon }}_0}{{\tilde{\epsilon }}^k}\left[1-\frac{({\tilde{\epsilon }}^k-{\tilde{\epsilon }}_0)}{({\tilde{\epsilon }}_{fi}-{\tilde{\epsilon }}_0)}\right]$$

(1)

The maximum critical strain is defined as ${\tilde{\epsilon }}_{fi}=\frac{2g_{IC}}{f_t}$ and the equivalent strain is evaluated according to Equation (2), where ${〈{\epsilon }_i〉}_{+}$ is the positive part of the main strain in the i direction, which means that ${〈{\epsilon }_i〉}_{+}={\epsilon }_i$, if ${\epsilon }_i\ge 0$; otherwise, ${〈{\epsilon }_i〉}_{+}=0$. However, for a state of uniaxial tension in the direction 1, this expression is reduced to ${\tilde{\epsilon }}^k={\epsilon }_1$.

$${\tilde{\epsilon }}^k=\sqrt{\sum _{i=1}^3{〈{\epsilon }_i〉}_{+}^2}$$

(2)

The stress–strain relationship of the material during a stage of damage can be expressed in relation to the undamaged stress–strain relationship, as depicted in Equation (3). Here, $\tilde{E}$ and $E_0$ represent the elastic modulus of the damaged and undamaged material, respectively, while D denotes the damage variable.

$$\sigma =\tilde{E}\epsilon ,\tilde{E}=E_0(1-D),(0\le D\le 1)$$

(3)

2.2. Generation of Random Mechanical Properties

2.2.1. Tensile Strength Random Generation

The Weibull distribution, established by Weibull [25,26], characterizes the heterogeneity of the tensile strength of the material. The probability density function $f_w(x,b,c)$, cumulative function $F(x,b,c)$, and inverse cumulative function $F^{-1}$ for a random variable $x\ge 0$, are outlined in Equations (4)–(6). Here, the terms $b>0$ and $c>0$ denote the shape and scale parameters governing this distribution.

$$f_w(x,b,c)=\frac{b}{c}{\left(\frac{x}{c}\right)}^{b-1}{e}^{-{\left(\frac{x}{c}\right)}^b}$$

(4)

$$F(x,b,c)=1-e^{-{\left(\frac{x}{c}\right)}^b}$$

(5)

$$F^{-1}{(u,b,c)=c(-ln(1-u))}^{1/b}$$

(6)

The mean ${\mu }_w$ and standard deviation ${\sigma }_w$ of the distribution can be evaluated, respectively, by Equations (7) and (8).

$${\mu }_w=c\Gamma \left(1+\frac{1}{b}\right)$$

(7)

$${\sigma }_w=\sqrt{c^2\Gamma \left(1+\frac{2}{b}\right)-{\mu }^2}$$

(8)

where $\Gamma$ is the gamma function, given by Equation (9),

$$\Gamma \left(\eta \right)={\int }_0^{\infty }{t}^{\eta -1}{e}^{-t}dt;\eta >0$$

(9)

For positive integer values of $\eta$, this function can be evaluated as $\Gamma \left(n\right)=(n-1)!$.

The numerical generation of random values for tensile strength relies on the inverse transformation method. The essential steps of this procedure are summarized in Algorithm 1. In this process, random numbers $\left(u\right)$, derived from the uniform distribution $\left(U\right(0,1\left)\right)$ within the range of 0 to 1, are generated. Subsequently, these generated values are applied to Equation (6), providing random numbers following the Weibull distribution.

Algorithm 1 Generation of random tensile strength by inversion method

1: for each finite element do 2:       Generate a random number u from $U(0,1)$; 3:       Given the parameters $(b,c)$, evaluate $f_t=F^{-1}{\left(u\right)=c(-ln(1-u))}^{1/b}$; 4:       Return $f_t$; 5: end for

2.2.2. Fracture Energy Random Generation

The lognormal distribution characterizes the heterogeneity within the material’s fracture energy. Its probability density function, denoted as $f_L(x,{\mu }_L,{\sigma }_L)$, is defined in Equation (10), where ${\mu }_L$ represents the mean and ${\sigma }_L$ denotes the standard deviation of the variable’s natural logarithm. Furthermore, the expected value $E_L\left(X\right)$ and variance $Va{r}_L\left(X\right)$ of this distribution can be calculated using Equations (8) and (11).

$$f_L(x,{\mu }_L,{\sigma }_L)=\frac{1}{{\mu }_L{\sigma }_L\sqrt{2\pi }}{e}^{-\frac{{(ln\left(x\right)-{\mu }_L)}^2}{2{\sigma }_L^2}}$$

(10)

$$E_L\left(X\right)=e^{{\mu }_L+\frac{{\sigma }_L^2}{2}}$$

(11)

$$Va{r}_L\left(x\right)=\left(e^{{\sigma }_L^2}-1\right)e^{2{\mu }_L+{\sigma }_L^2}$$

(12)

The numerical generation of random fracture energy relies on the convolution method, with the fundamental steps summarized in Algorithm 2. The approach involves initially generating random numbers from $U(0,1)$ in order to after obtain random variables conforming to the standard normal distribution $N(0,1)$ (a normal distribution with ${\mu }_L=0$ and ${\sigma }_L=1$). Subsequently, through a straightforward algebraic operation, lognormally distributed random variables are obtained. Additional details regarding these random number generation procedures can be found in the work by Rita [16].

Algorithm 2 Generation of random fracture energy by the convolution method

1: for each finite element do 2:       Generate a sequence of 12 random values $u_1,u_2,...,u_{12}$ from $U(0,1)$; 3:       Evaluate $X=\left({\displaystyle \sum _{i=1}^n}{u}_i\right)-6$; // X follows N(0, 1) 4:       Evaluate $G_{IC}=e^{{\mu }_L+{\sigma }_{L}X}$; 5:       Return $G_{IC}$; // $G_{IC}$ follows the lognormal distribution 6: end for

2.3. Estimation of the Model’s Parameters

The estimation of Weibull distribution parameters relies on the correlation of equations for ${\mu }_w$ and ${\sigma }_w$ with the analytical scale law proposed by Rossi et al. [22]. This scale law, derived from extensive experimental investigations that correlate concrete heterogeneity and scale effects, plays a significant role in validating probabilistic models. The expressions for these correlations are detailed in Equations (13)–(16). Through these equations, the mean $\left({\mu }_{f_t}\right)$ and coefficient of variation of tensile strength $({\sigma }_{f_t}/{\mu }_{f_t})$ for a given volume of concrete can be estimated, taking into account parameters such as specimen volume $\left(V_s\right)$, maximum aggregate volume $\left(V_a\right)$, and material compressive strength $\left(f_c\right)$.

$${\mu }_{f_t}=a{\left(\frac{V_s}{V_a}\right)}^{-\gamma }$$

(13)

$$\frac{{\sigma }_{f_t}}{{\mu }_{f_t}}=A{\left(\frac{V_s}{V_a}\right)}^{-B}$$

(14)

where $a=6.5$ MPa, $\gamma$ is given by Equation (15), $A=0.35$, and B is given by Equation (16). Moreover, the constant $c_1=1$ MPa.

$$\gamma =0.25-3.6\times {10}^{-3}\left(\frac{f_c}{c_1}\right)+1.3\times {10}^{-5}{\left(\frac{f_c}{c_1}\right)}^2$$

(15)

$$B=4.5\times {10}^{-2}+4.5\times {10}^{-3}\left(\frac{f_c}{c_1}\right)-1.8\times {10}^{-5}{\left(\frac{f_c}{c_1}\right)}^2$$

(16)

These analytical expressions are applicable to various concrete types with compressive strengths ranging from 35 to 130 MPa, excluding fiber-reinforced and lightweight concretes. They are specifically utilized in the model at the elemental level, i.e., considering the finite element volume $\left(V_e\right)$. Consequently, the combination of Equations (7), (8), (13) and (14) results in a nonlinear system of equations.

To solve this system, an iterative numerical procedure is employed. Through this iterative process, a pair of $(b,c)$ values is determined for each finite element based on its $V_e$, as well as the concrete’s $V_a$ and $f_c$. Detailed implementation of this iterative procedure can be found in Rita [16].

The lognormal distribution characterizes the heterogeneity of fracture energy, typically governed by mean and standard deviation parameters for estimation. However, experimental findings by Rossi [28] suggest that fracture energy can be regarded as an intrinsic material property, whereby its mean value $\mu \left(G_{IC}\right)$ remains constant and aligns with experimentally obtained values. Consequently, in this context, the sole parameter necessitating determination is the standard deviation of the distribution, denoted as ${\sigma }_L\left(G_{IC}\right)$.

To estimate ${\sigma }_L\left(G_{IC}\right)$, an inverse analysis procedure is adopted. This method involves conducting multiple simulations of a macrocrack propagation test on a large double cantilever beam specimen (DCB). Utilizing the outcomes from this procedure, a function relating ${\sigma }_L\left(G_{IC}\right)$ to the heterogeneity degree of mesh elements is proposed. A comprehensive description of this inverse analysis procedure is detailed in Section 3.

2.4. Parallel Monte Carlo Implementation

Working with probabilistic models often presents a significant challenge due to the necessity of conducting numerous analyses for the same problem, resulting in high computational costs. This challenge becomes even more pronounced when conducting inverse analyses for three-dimensional models, significantly escalating computational expenses. Consequently, optimizing computational time becomes imperative to ensure the feasibility of the proposed methodology. In addressing this issue, the Monte Carlo (MC) procedure has been implemented with a parallelization strategy. This strategy serves as a crucial resource, enabling the execution of three-dimensional simulations within a viable computational execution time.

The implementation of parallelization in this work employs OpenMP (Open Multi-Processing), an application programming interface (API) supporting multi-platform shared-memory parallel programming in C/C++ and Fortran [29,30,31]. The parallel algorithm for the Monte Carlo method operates based on the principles of the sequential algorithm but adopts a master–slave parallel model. This model facilitates the management of the program, enabling the simultaneous execution of several independent finite element method (FEM) analyses.

A pseudo-code outlining the Parallel Monte Carlo (PMC) procedure is presented in Algorithm 3. The selection of this parallelization strategy stems from its ease of implementation and compatibility with the proposed methodology, which necessitates conducting numerous independent analyses. However, an in-depth discussion regarding the intricacies of parallelization extends beyond the scope of this work and is available in Rita [16] for further insight.

Algorithm 3 Parallel Monte Carlo Procedure

1: Variables initialization 2: Read maxmc and stop criteria; //maxmc is the maximum number of samples 3: Read nthreads; //nthreads is the number of threads 4: j = 0;//counter of performed samples 5: while stopping criterion is not reached 6:    j = j + nthreads 7:    Initialize the parallel environment 8:          For each block with nthreads samples do 9:                Perform the FEM analysis 10:                Save $(P,\delta )$ curve and relevant results of the sample 11:          end for 12:    Finalize the parallel environment 13: end while 14: Perform the statistical analysis

3. Estimation of the Coefficients of Variation of Fracture Energy

The term “coefficients of variation” is employed in the plural form in this context because, within the framework of the current probabilistic model, the coefficient of variation of $G_{IC}$ varies based on the volume of each finite element. Consequently, determining these coefficients of variation stands as a primary challenge concerning the model’s practical application and operational aspects.

To achieve this objective, an inverse analysis necessitates an experimental study providing detailed information on macrocrack propagation, including the evolution of macrocrack length and opening. This experimental study must be conducted at a sufficient distance to yield representative data concerning the dispersion associated with this propagation. Additionally, it should facilitate the determination of the intrinsic $G_{IC}$. However, determining these toughness characteristics experimentally poses a significant challenge.

Previous studies have demonstrated that experimental tests need to be conducted on large concrete specimens to accurately determine these characteristics [9,28,32,33,34,35,36,37,38]. This necessity arises due to the high level of concrete heterogeneity associated with larger aggregate sizes. As a result, the process zone at the tip of the propagating macrocrack can reach approximately 30 cm [28].

The inverse analysis in this study relies on an experimental investigation conducted on a large-scale Double Cantilever Beam (DCB) specimen detailed in [28,38]. The dimensions of the DCB specimen used in this experiment were 3.5 m in length, 1.1 m in width, and 0.3 m in depth.

This particular DCB test conducted on such a substantial specimen stands as a unique case in current concrete studies. It is noteworthy because it provided genuine and intrinsic values of $K_{IC}$ (mode I critical stress intensity factor) or $G_{IC}$ (mode I specific fracture energy). It is important to emphasize that most other experimental studies have yielded values that are sensitive to scale effects rather than being intrinsic [28,38]. Due to the computational demands associated with simulating this type of DCB test specimen, only one concrete composition was studied in our analysis.

3.1. A Brief Presentation of the Experimental Study Performed on the DCB Specimen

The geometry and loading conditions of the DCB test are depicted in Figure 2. The test is governed by the application of a constant notch displacement rate set at 25 $\mathsf{\mu }$m/min. Throughout the test, crack propagation occurs from the bottom to the top of the structure. The load is applied 17.5 cm from the lower side of the beam. Crack opening is measured by averaging the results from two extensometers—one on each side. To guide and maintain the crack within the median plane, the section was thinned and a longitudinal prestressing (via post-tensioning) was introduced using multiple cables.

During the test, the load/notch opening curve is recorded. Upon reaching the peak load, signifying the initiation of the macrocrack, the specimen undergoes several cycles of notch closing and reopening, all at the same frequency. These cycles enable the tracking of the specimen’s stiffness evolution, a critical factor used to calculate $K_{IC}$ and $G_{IC}$ employing the classical compliance method, a well-established technique in linear fracture mechanics theory [28,38,39,40]. A more detailed description of the entire test procedure can be found in [28,38].

The analysis of the Double Cantilever Beam (DCB) test using the compliance method results in the derivation of a curve depicting $K_{IC}$ versus the equivalent crack length (ECL). This ECL represents an idealized crack length—essentially a straight and smooth crack within a linear elastic material—that results in the same compliance as experimentally obtained from the specimen. Conceptually, this ECL can be visualized as an average crack length, encompassing the real macrocrack along with its associated process zone—a microcracked region at the macrocrack tip during propagation. Typically, the evolution of the specimen compliance with the idealized crack length (analytical relation) is obtained by performing linear elastic finite element analyses [28,38].

The experimental test campaign involved various types of concrete, including plain concrete with and without additional prestressing loads, concrete with fibers, and reinforced concrete. However, this numerical study focuses solely on the cases involving plain concrete. The mix design details for plain concrete are provided in

Table 1. Mix design of the plain concrete.

Constituent	Quantity
Aggregate (4–12 mm)	1105 (kg/m3)
Sand (0–5 mm)	700 (kg/m3)
Cement	400 (kg/m3)
Water	190 (L/m3)

as reported in [28,38]. The values for compressive strength ($f_c$), splitting tensile strength ($f_t$), the Young’s modulus (obtained from average values during compressive tests), and the mean values of $G_{IC}$ for plain concrete are provided in

Table 2. Mechanical properties of the plain concrete.

Mechanical Property	Value
Compressive strength	54 MPa
Splitting tensile strength	4.1 MPa
Young modulus	35.5 GPa
Fracture energy ($G_{IC}$)	1.25 $\times {10}^{-4}$ MN/m

.

During the tests conducted with the DCB specimen composed by plain concrete, three prestressing force (F) conditions were applied: case 1 with $F_1=1230$ KN, case 2 with $F_2=682$ KN, and case 3 with no prestressing force ($F_3=0$ KN). To focus on using the most rectilinear crack propagation, only the test with the highest prestressing load ($F_1$) is utilized for the numerical simulations. The purpose of prestressing load was to guide the crack propagation to the central part of the specimen (in the thinnest zone) and prevent any premature bifurcation, as observed in the case with no prestressing load. Figure 3 displays the load–notch opening curve obtained from the considered experimental test. It is noteworthy that the curve exhibits a softening phase, characterized by a continuous decrease, which correlates with the very rectilinear propagation of the crack.

3.2. Inverse Analysis Methodology

In accordance with the model’s hypotheses, the parameters need to be determined relative to the finite element volume. Hence, the estimation is conducted by considering three different finite element mesh refinements to determine the relationship between the standard deviation of $G_{IC}$ and the ratio $V_e/V_a$, where $V_e$ represents the volume of the finite elements and $V_a$ corresponds to the volume of the larger aggregate in the concrete composition. Specifically, three meshes, referred to as Mesh 1, Mesh 2, and Mesh 3, comprising linear tetrahedral elements, are taken into account to perform this inverse analysis. Their respective three-dimensional representations are illustrated in Figure 4. The finite element meshes were produced using the software Coreform Cubit [41].

It is noticeable that the elements comprising the meshes exhibit considerable variations in size. As the focus lies on assessing the behavior of macrocrack propagation, predominantly occurring within the central part of the specimen, an average volume of the elements situated in this area will be considered. Hence, the mean heterogeneity degree for each mesh refinement, denoted as $r_e^{mean}$ = $V_e^{mean}$/$V_a$, is established. This assumption seems rational as it concentrates solely on the elements in this specific region, undergoing the damage process and being subject to random variations in $G_{IC}$.

Additional details regarding the three meshes are presented in

Table 3. Description of the meshes.

	Neq	Nnode	Numel		Heterogeneity Degree
	Total	Total	Total	Center	${\mathit{r}}_{\mathit{e}}^{\mathit{mean}}$	${\mathit{r}}_{\mathit{e}}^{\mathit{sd}}$	${\mathit{r}}_{\mathit{e}}^{\mathit{min}}$	${\mathit{r}}_{\mathit{e}}^{\mathit{max}}$
Mesh 1	3786	1370	5222	3953	114.7	50.7	34.1	332.3
Mesh 2	11,848	4135	19,564	15,143	29.3	8.1	8.2	78.4
Mesh 3	23,273	8038	39,600	34,009	10.8	2.8	3.1	31.5

, indicating: the total number of equations (Neq) and nodes (Nnode) encompassing the entire specimen; the count of elements (Numel); and the mean heterogeneity degree at the center of the specimen. This includes the standard deviation ($r_e^{sd}$), as well as its minimum ($r_e^{min}$) and maximum ($r_e^{max}$) values. The refinement of Mesh 3 was specifically designed to accounts for elements with volumes exceeding that of the maximum aggregate in the concrete composition, ensuring that the elements represent a heterogeneity degree $r_e=V_e/V_a\ge 1$.

In Figure 5, the frontal view of Mesh 3 illustrates the locations where the prestressing forces and imposed displacements are applied. The boundary conditions employed to accurately represent the DCB test behavior are as follows: restriction of displacements along the X-axis within the YZ plane, along the Y-axis within the XZ plane, and along the Z-axis at the central nodes in the XZ plane. Additionally, Z-axis displacements at nodes with prescribed displacements are restricted. Prestressing force is simulated by applying prescribed compression forces in the Y-direction on the elements situated on the specimen’s bottom surface at the initial load step.

Each Monte Carlo simulation involves the execution of one hundred independent Finite Element Method (FEM) analyses, which proves to be a sufficient number for ensuring consistent results. An in-depth analysis about the influence of the Monte Carlo number of samples on the global response of this structural problem is detailed in [16], assuring the convergence of the numerical results in terms of mean numerical curve. These simulations encompass varied values of the standard deviation within the lognormal distribution of $G_{IC}$. To achieve this, ten distinct values of ${\sigma }_L\left(G_{IC}\right)$ were selected within the range of $[0,100\times G_{IC}]$. The configuration of ${\sigma }_L\left(G_{IC}\right)$ values is presented as follows:

$$\begin{array}{cccc}S1\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=0\times G_{IC}\hfill & S6\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=40\times G_{IC}\hfill \\ S2\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=5\times G_{IC}\hfill & S7\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=50\times G_{IC}\hfill \\ S3\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=10\times G_{IC}\hfill & S8\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=70\times G_{IC}\hfill \\ S4\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=20\times G_{IC}\hfill & S9\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=90\times G_{IC}\hfill \\ S5\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=30\times G_{IC}\hfill & S10\hfill & ⟶{\sigma }_L\left(G_{IC}\right)=100\times G_{IC}\hfill \end{array}$$

It is noteworthy to observe that simulation S1 represents a scenario where no dispersion exists concerning fracture energy, implying that all elements within the mesh receive an identical value of $G_{IC}$ as a parameter. Hence, simulation S1 exhibits a deterministic aspect in relation to fracture energy. However, it is crucial to note that tensile strength remains probabilistic across all cases. This decision stems from the intention to demonstrate that within the framework of the proposed numerical model, the probabilistic approach (regarding both parameters) stands as the sole pertinent solution. Throughout this numerical study, a total of three thousand finite element analyses were conducted, all featuring an imposed displacement increment value of 25 $\mathsf{\mu }$m.

3.3. Strategy to Choose the Best Value of ${\sigma }_L\left(G_{IC}\right)$ for Each Mesh Refinement

The approach for selecting the optimal ${\sigma }_L\left(G_{IC}\right)$ value for each mesh refinement relies on two mechanical outcomes from the experimental study: (1) the load–notch opening displacement (load–NOD) curve, and (2) the equivalent crack propagation length (ECL). To determine the appropriate ${\sigma }_L\left(G_{IC}\right)$ values, the strategy involves the following two steps:

(1)

First step: Comparison between numerical and experimental load–NOD curve. During this step, the following criteria are considered:

(a)

The experimental curve must lie within the range of numerical curves.

(b)

Minimization of the variability within the bundle of the numerical curves as possible. This emphasis on minimizing variability is crucial since the load–NOD curve characterizes a structural behavior (global) rather than material behavior (local), which typically exhibits more variability.

At the conclusion of this initial, two ${\sigma }_L\left(G_{IC}\right)$ values are selected for each mesh refinement. These values align most closely with the two aforementioned criteria and will proceed to the subsequent analysis in the second step.

(2)

Second step: Comparison between numerical and experimental macrocrack propagation length (in terms of ECL). Determining the macrocrack propagation length is the primary objective of the proposed model, thus serving as the final criterion for the optimal solution. In this comparison, it is essential to clarify some points:

(a)

The comparison is made concerning the last notch opening displacement (NOD) for which an equivalent crack length has been calculated. In the experimental work of Rossi [28], the last NOD is 1560 $\mathsf{\mu }$m, and the corresponding last ECL is 2.09 m.

(b)

Since the model defines the macrocrack as a series of fully damaged elements, the numerical macrocrack should generally be smaller than the experimental one (represented by the ECL) for identical notch openings. This size difference is expected to decrease with mesh refinement.

(c)

Numerically, the macrocrack tip is identified by the last fully damaged finite element. However, in experimental studies, surface cracks might appear larger than those within the specimen due to different stress situations. This discrepancy might lead to the equivalent crack being smaller than the visible surface crack. Additionally, numerical cracks could be larger than the equivalent crack for meshes comprising multiple elements in the specimen’s thickness. This scenario is particularly observed in Mesh 3 within the scope of this study.

4. Numerical Results

The outcomes are displayed in terms of load–notch opening displacement (NOD) curves, representing global behavior. Each figure includes the experimental curve alongside the calculated ones. It is crucial to note that NOD is directly measured in experiments and it is numerically computed at specific loading points.

This paper does not present all numerical outcomes due to space limitations. However, a comprehensive collection of these results is available in Rita [16]. The presentation herein initiates after the completion of the first step in the strategy, where two ${\sigma }_L\left(G_{IC}\right)$ values have been retained for each mesh refinement. These chosen values will be the focus of the subsequent analysis.

4.1. Importance of the Probabilistic Approach and the Adequate Choice of Random Properties

In Figure 6, a comparison is presented between the three numerical mean curves related to simulation S1 for the three meshes, representing the deterministic aspect of fracture energy (S1) as explained in Section 3.2. The comparison is made against the experimental curve. These simulations appear to inadequately represent the behavior observed in the DCB test. This outcome underscores the significance of the probabilistic approach and justifies the need to randomly distribute the proposed mechanical properties, tensile strength, and fracture energy within the model’s framework.

4.2. First Step of the Strategy

This section presents the results for the two best cases associated with each mesh refinement. In Figure 7, the Monte Carlo responses for the best cases related to M1 are displayed. These cases correspond to S4 with ${\sigma }_L\left(G_{IC}\right)=20\times G_{IC}$ and S5 with ${\sigma }_L\left(G_{IC}\right)=30\times G_{IC}$. Additionally, Figure 8 illustrates the Monte Carlo responses for the best cases of M2, which are associated with S5 having ${\sigma }_L\left(G_{IC}\right)=30\times G_{IC}$ and S6 with ${\sigma }_L\left(G_{IC}\right)=40\times G_{IC}$. Finally, Figure 9 showcases the global responses for the best cases regarding M3, linked to S6 with ${\sigma }_L\left(G_{IC}\right)=40\times G_{IC}$ and S7 with ${\sigma }_L\left(G_{IC}\right)=50\times G_{IC}$.

The observation of Figure 7, Figure 8 and Figure 9 leads to the following comments:

1.

The post-peak section of the curves exhibits the experimental global response consistently within the range of numerical responses.

2.

Except for Mesh 1, the peak load in the numerical simulations are lower compared to those observed in the experimental test. Moreover, this difference increases from Mesh 2 to Mesh 3. Several reasons can explain this behavior:

(a)

The peak load is correlated with the process of macrocrack localization (initiation) rather than the macrocrack propagation, which the proposed model is specifically designed to address. It focuses on the softening behavior of the specimen’s global response during macrocrack propagation, not the microcracking to macrocracking transition.

(b)

The numerical modeling of this macrocrack propagation employs linear volume elements, i.e., three-dimensional finite elements with linear interpolation functions. However, typically, strong stress concentrations at the narrow notch tip require special crack elements with different formulations or elements utilizing more complex interpolation functions, often nonlinear, to accurately represent such stress concentrations.

(c)

Using linear elements (as in this case) results in higher stress concentrations at the notch tip for Meshes 2 and 3 compared to Mesh 1, considering the same notch opening. Consequently, for Meshes 2 and 3, the tensile stresses at the notch tip are higher, leading to smaller notch openings and loads initiating cracking (D = 1) in the first elements at the notch tip compared to Mesh 1. Thus, the mean number of cracked elements at the moment of maximum peak load increases with mesh refinement. Specifically, the mean number of fully cracked elements at peak load for the presented simulations in Figure 7, Figure 8 and Figure 9 are: 6.6 (Mesh 1); 75.1 (Mesh 2); and 140.4 (Mesh 3) [16]. Consequently, the peak load for Meshes 2 and 3 is lower compared to Mesh 1 due to a higher number of cracked elements at the peak-load stage.

4.3. Second Step of the Strategy

Figure 10 showcases surface views of three macrocrack patterns associated with the three meshes during the final stage of the simulation, characterized by an imposed displacement of 1560 $\mathsf{\mu }$m. It is apparent that Mesh 3 exhibits a straighter crack pattern with a longer macrocrack length.

Table 4. Results of the mean numerical crack length (m).

	Mesh 1		Mesh 2		Mesh 3
	S4	S5	S5	S6	S6	S7
$C{L}_{NUM}$	1.975	2.011	1.97	2.026	2.075	2.113

outlines the mean values of the numerical crack lengths corresponding to the three mesh refinements, denoted as $C{L}_{NUM}$. These numerical mean values are to be juxtaposed with the experimental ECL, which stands at 2.09 m.

Analyzing the findings presented in Table 4, the following conclusions emerge: for Mesh 1 and Mesh 2, the optimal configurations are, respectively, S5 (${\sigma }_L\left(G_{IC}\right)=30\times G_{IC}$) and S6 (${\sigma }_L\left(G_{IC}\right)=40\times G_{IC}$). In these instances, the numerical crack length consistently appears slightly smaller but remains in close proximity to the equivalent crack length. Conversely, for Mesh 3, drawing a definitive conclusion is less straightforward. As mentioned in Section 3.3, in this scenario, the mean numerical crack length might be larger if the crack length at the specimen’s surface surpasses that in its interior. Evidence supporting faster crack propagation at the specimen’s surface is evident in Figure 11. This figure provides a comparison between the macrocrack pattern obtained for the two surfaces and a cut view, illustrating the interior of the specimen.

Considering both Figure 10 and Figure 11, the most suitable configuration for Mesh 3 appears to be S7 (${\sigma }_L\left(G_{IC}\right)=50\times G_{IC}$). The observations from Table 4 and Figure 11 imply that the macrocrack pattern is better defined with the most refined mesh. This inference aligns with the logical progression that as mesh refinement increases, precision in the information also improves. However, this enhancement in precision comes at the expense of increased computational cost in terms of simulation runtime.

4.4. Analysis of the Results

At the conclusion of the second step of the strategy, it was observed that ${\sigma }_L\left(G_{IC}\right)$ increases with the mesh refinement. This outcome aligns with the evolution of the concrete tensile strength standard deviation correlated to the heterogeneity degree [15,22]. Utilizing the optimal parameters (Section 4.3) determined for each $r_e^{mean}$ value (Table 3), a function has been formulated to estimate ${\sigma }_L\left(G_{IC}\right)$ corresponding to the heterogeneity degree of each finite element. This function, described in Equation (17), has been extrapolated to encompass all elements included in the three mesh refinements.

$${\sigma }_L\left(G_{IC}\right)=\left(Aln\left(r_e\right)+B\right)\times {\mu }_L\left(G_{IC}\right)$$

(17)

where $A=-8.538$, $B=70.88$, and ${\mu }_L\left(G_{IC}\right)$ is the mean value of $G_{IC}$.

Figure 12 illustrates the behavior of the estimated function concerning the coefficient of variation of $G_{IC}$, represented as ${\sigma }_L\left(G_{IC}\right)/{\mu }_L\left(G_{IC}\right)$. A focused view of the initial section of the function domain (encompassing the heterogeneity degree of all elements at the central part of the DCB specimen, $r_e\in [3,350]$) is presented in Figure 12a, while Figure 12b showcases the complete domain of the proposed function.

Analyzing the extrapolated function from Equation (17) raises caution due to the inherent risks associated with extrapolating from only three data points. Therefore, it is imperative to approach this function with caution, emphasizing the need for future validation. To initiate the validation process using available experimental information, the following pertinent points merit examination:

1.

If Equation (17) is taken into account, it evident that ${\sigma }_L\left(G_{IC}\right)$ equals zero at a value of $r_e\approx 4030$.

2.

In the work developed by Rossi [28], acoustic emission techniques were employed to assess the dimensions of the process zone at the macrocrack tip. The findings revealed a process zone volume approximately 3600 cm³ (10 cm high, 10 cm wide, 30 cm long). Considering the maximum aggregate volume of 1.13 cm³ (diameter of 12 mm, Table 1), the ratio of the process zone volume ($V_{FPZ}$) to the larger aggregate volume ($V_a$) is $V_{FPZ}/V_a\approx 3185$.

3.

Additionally, in the study of Rossi [28], during the macrocrack propagation (exceeding 2 m in length), it was evaluated that $\sigma \left(G_{IC}\right)=0.073\times {10}^{-4}$ MN/m for $\mu \left(G_{IC}\right)=1.25\times {10}^{-4}$ MN/m.

Hence, the deduction that the critical fracture energy’s standard deviation reaches zero when the finite element heterogeneity degree equals 4030 appears to be logical. Therefore, it could be reasonable to assume that ${\sigma }_L\left(G_{IC}\right)$ continues to remain at zero for a mesh refinement related to a value of $r_e>4300$. Additionally, it is important to highlight that Equation (17) is specifically linked to the proposed numerical model and is not an inherent property of the concrete investigated in this study. Consequently, if different formulations of finite elements or alternative softening curves in tension are employed, this relationship might undergo alterations.

A comprehensive review of Figure 7, Figure 8 and Figure 9 yields the following insights and comments: (a) The experimental global response consistently resides within the bundle of numerical responses concerning a notch opening displacement (NOD) around 0.4 mm. This implies that within the context of the DCB study, the numerical model becomes pertinent from this particular NOD value. (b) It is crucial to note that this NOD value is not equivalent to the crack opening. The observed value is higher due to a certain level of flexibility inherent in the DCB specimen. (c) To assess the domain of crack opening where the model remains acceptable—indicating a lack of significant mesh dependency—the relationship between the numerical mean crack opening at the crack tip and the notch opening has been computed for the three mesh refinements. This analysis is presented in detail in Figure 13.

5. Conclusions

In this paper, a three-dimensional semi-explicit probabilistic cracking numerical model was introduced which was developed within the finite element method framework. This model focuses solely on macrocrack propagation employing linear volume elements. The probabilistic nature of the model arises from the distribution of mechanical properties that depends on the volume of finite elements. Consequently, the primary objective of this study was to determine this dependency, crucial for the model’s development.

The parameter estimation process in our proposed strategy primarily involves determining the standard deviation of fracture energy. This estimation is achieved through an inverse analysis procedure, conducted via a fracture mechanical test simulation. The test, features macrocrack propagation spanning over two meters in a vast Double Cantilever Beam (DCB) specimen. The inverse approach yielded promising results. It demonstrated the model’s capability to accurately evaluate crack length and opening during propagation.

However, it is imperative to stress the necessity of validating the results obtained from our proposed inverse analysis and establishing the relationship between ${\sigma }_L\left(G_{IC}\right)$ and finite element volume. Future work should involve numerical simulations of large structures using this relation for distinct mesh refinements.

Furthermore, a key limitation of our study is that the ${\sigma }_L\left(G_{IC}\right)$ relationship is specific to a single concrete mix design. Therefore, to enhance the practical applicability of our proposed model in engineering settings, efforts must be made to generalize this relationship across different concrete mix designs—an essential challenge that needs addressing.

Lastly, the utilization of parallel Monte Carlo played a pivotal role in executing the inverse analysis procedure, encompassing three thousand finite element analyses. This optimization significantly reduced computational execution time, enabling the exploration of three-dimensional problems.