# Calibration of non-local damage models from full-field measurements: Application to discrete element fields

Louis Védrine, Flavien Loiseau, Cécile Oliver-Leblond, Rodrigue Desmorat

► To cite this version:

Louis Védrine, Flavien Loiseau, Cécile Oliver-Leblond, Rodrigue Desmorat. Calibration of non-local damage models from full-field measurements: Application to discrete element fields. European Journal of Mechanics - A/Solids, 2025, 112, pp.105611. 10.1016/j.euromechsol.2025.105611. hal-05005942

HAL Id: hal-05005942

https://hal.science/hal-05005942v1

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Full length article

# Calibration of non-local damage models from full-field measurements: Application to discrete element fields

Louis Védrine<sup>a,b, D,*</sup>, Flavien Loiseau<sup>a,c, D</sup>, Cécile Oliver-Leblond<sup>a, D</sup>, Rodrigue Desmorat<sup>a, D</sup>

## ARTICLE INFO

Keywords:

Quasi-brittle materials

Eikonal equation

Regularization

Evolving non-local interactions

Beam-Particle model

Full fields measurements

## ABSTRACT

Continuous damage models are increasingly used in numerical simulations to design structures, but their local formulations are sensitive to mesh size and present localization of strains in an infinitely thin region. To overcome these problems non-local damage models and related regularization methods, have been developed introducing a characteristic/internal length, thus avoiding pathological mesh dependence. Those methods make the damage evolution depends on the mechanical quantities at the current material point (local) and its neighborhood (non-local).

Existing approaches to calibrate the internal length use global quantities in the calibration process, although local data is now becoming accessible (e.g., using digital image correlation). In this study, we investigate the use of full-field displacement measurements and propose a new methodology for calibrating a non-local damage model based on local field measurements and we apply it to calibrate the Eikonal Non-local Gradient (ENL-G) approach from the measured strain and damage field. After detailing the calibration procedure, we then apply it on a simple ideal case. We illustrate, and analyze the robustness of the calibration procedure with respect to the choice of evolution law and measurement noise of the proposed calibration method. To confront the procedure to a more realistic case, we employed a 2D beam-particle model. This discrete model is first identified with respect to the size and shape effect based on one of the comprehensive experimental data sets available in the literature, including four shapes with three sizes each. Then, it is used to generate a “reference” evolution of the damage and strain fields in beams of different sizes subjected to uniaxial tension.

The parameters of the discrete model used have been calibrated to represent the scale and size effects, giving a very good representation of the experiments. We also illustrate the evolution of non-local interactions in the Eikonal approach using Green functions. Finally, the application of the calibration procedure shows that it is possible to determine the internal length of the non-local problem studied as well as the damage evolution law and its parameters. The outcomes of this study contribute to shed light on a new methodology to identify non-local damage models based on full-field measurements, and call for experimental size effect campaign with displacement field.

## Introduction

Quasi-brittle materials are massively used in civil engineering; modeling and predicting their mechanical degradation is necessary to guarantee the integrity of structures (e.g., nuclear power plants, and dams). Hence, it is crucial to be able to predict the mechanical behavior of this structure. To model the material degradation and its impact on the behavior, local damage models have been developed (Kachanov, 1958; Rabotnov et al., 1970; Lemaitre, 1971). These models are often

used in structural design. Nevertheless, local damage models lead to strain localization into an infinitely thin region due to their softening nature and their dissipated energy tend to zero. Thus, the numerical simulations using local damage model are mesh-dependent (Hillerborg et al., 1976; Bažant, 1986; Benallal et al., 1989). To overcome the spurious mesh sensitivity of continuous damage models and to describe the size effect at the structural level, non-local damage models were introduced (Pijaudier-Cabot and Bažant, 1987; Peerlings et al., 1996).

<footnote><sup>a</sup> Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS — Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Université Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Etudes de la Neige, Grenoble, France</footnote>

<footnote><sup>c</sup> IMSIA, ENSTA Paris, CNRS, EDF, Institut Polytechnique de Paris, 91120 Palaiseau, France</footnote>

<footnote>* Corresponding author at: Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS — Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France.</footnote>

<footnote>E-mail addresses: louis.vedrine@ens-paris-saclay.fr (L. Védrine), flavien.loiseau@ensta-paris.fr (F. Loiseau), cecile.oliver@ens-paris-saclay.fr (C. Oliver-Leblond), rodrigue.desmorat@ens-paris-saclay.fr (R. Desmorat).</footnote>
Among the great diversity of non-local damage models found in the literature (Bažant and Jirásek, 2002b; Jirásek, 2004; Pijaudier-Cabot and Grégoire, 2014), a main feature is the introduction of a characteristic length into the constitutive relations. To be more precise, in the damage evolution law, the damage growth is driven by a non-local quantity being a weighted spatial average of the local damage driving quantity. The characteristic length of a non-local continuum was initially estimated experimentally using an energetic approach, defined as the ratio of the fracture energy to the energy dissipated per unit volume of the material (Bažant and Pijaudier-Cabot, 1989). Their validations are mainly based on criteria of “good” damage behavior (Jirásek, 1998), its ability to reproduce experimental results and to be insensitive to mesh size.

To properly reproduce experimental results, the constitutive equations for quasi-brittle material must reproduce the size effect (Bažant and Jirásek, 2002a; Alam et al., 2014). This size effect is induced by the existence of a Fracture Process Zone (FPZ) of finite size. The non-local feature induces a diffusion of the damage over a region of its characteristic length, which can be interpreted as the modeling of the process zone. Different non-local damage models will generate different shapes of the process zone as well as different spatial distributions of damage inside it, leading to a different size effect (Havlásek et al., 2016; Zhang et al., 2021). While a relation between the internal length(s) of a model and the size of the process zone can be found (Bažant and Jirásek, 2002a; Pijaudier-Cabot and Grégoire, 2014), the main issue for modeling is finding the internal length (a material parameter) of any non-local damage model best-fitting the reference data. Consequently, the calibration of these models has emerged as a significant challenge in the application of this damage modeling approach (Bažant et al., 1990).

Indeed, direct calibration is not feasible, the primary focus of calibration work on non-local models has predominantly been on global approaches (Grégoire et al., 2013). Based on the sensitivity of the size effect to the internal length, Le Bellégo et al. (2003) proposed a calibration procedure based on a three-point bending test, taking into account the size effect on notched specimens to calibrate the internal length of the non-local model. However, this method addresses only the macroscopic behavior of the beam while neglecting its local behavior. Moreover, this methodology does not enable the decoupling of the validation of the evolution law from the calculation of the non-local variable. Furthermore, it has been demonstrated that this approach results in non-unique solutions due to the ill-posed nature of the inverse problem (Iacono et al., 2006). This limitation necessitates the incorporation of additional local experimental data into the inverse procedure, such as the spatial distribution of dissipated energy, fracture patterns and/or strain field.

A common challenge in this context is obtaining local information, such as the micro-cracking state (Mindess, 1991). Acoustic emission is a widely utilized method for assessing the cracking state of quasi-brittle materials. However, linking these emissions to material damage remains challenging (Rossi et al., 1990; Landis, 1999; Alam et al., 2014; Thirumalaiselvi and Sasmal, 2021). Nevertheless, acoustic measurements from fracture events can be employed to evaluate energy dissipation. These emissions are spatially localized, and their intensity can be used to estimate local energy dissipation (Mihashi et al., 1991; Otsuka and Date, 2000; Haidar et al., 2005; Muralidhara et al., 2010; Grégoire et al., 2015).

Since the 2000s, displacement field measurements, such as those obtained using Digital Image Correlation (DIC) techniques (Choi and Shah, 1997; Sutton et al., 2009), have become accessible and opened new avenues for identifying local damage properties (Geers et al., 1999; Claire et al., 2002, 2003, 2004; Geymonat et al., 2002). Subsequently, full-field measurements facilitated the identification of elastic properties and damage fields on material surfaces by coupling elasticity and damage states (Lemaitre and Marquis, 1992). These methods have advanced the detection of micro-cracking (Ben Azzouna et al., 2011;

Chateau et al., 2018) and improved the calibration and validation of local damage models (Hild et al., 2015; Roux and Bouchard, 2015; Nguyen et al., 2015), incorporating local quantities rather than relying on single strain measurements.

Based on these local field measurements, various methods have been proposed to integrate them into the calibration of non-local damage models. For example, the size of the FPZ can be directly determined from displacement field measurements (Skarżyński et al., 2011; Alam et al., 2013; Doll et al., 2017; Tang et al., 2020; Li et al., 2020). Similarly, dissipated energy can also be used to determine the size of the FPZ of a quasi-brittle material (Grassl et al., 2012; Grégoire et al., 2015). Xenos et al. (2015) calibrated a non-local damage model by matching experimental and numerical standard deviations of the spatial distributions of dissipated energy densities. This approach allows for the incorporation of additional averaged local quantities into the calibration process by assuming a predefined profile of dissipated energy.

Geers et al. (1999) were among the first to use the full displacement field to calibrate the internal length of a non-local damage model used for composite materials. Digital Image Correlation (DIC) was utilized to obtain the displacement field, and the model was calibrated by minimizing the pixel-by-pixel error between the simulated and experimental displacement fields. Recently, Kosin et al. (2024) proposed calibrating the internal length of a phase-field fracture model using displacement field measurements. In this method, the phase-field model generates a series of displacement fields, and the internal length is identified through integrated digital image correlation.

The aim of this study is to propose a general methodology for calibrating any non-local damage model based on local displacement fields. For a homogeneous material, as the damage evolution should be the same for every point of the domain, the relation between damage growth ($\dot{D}$) and the damage driving quantity should always be the same. Moreover, the damage driving variable is non-local.

Yet, choosing the wrong non-local damage driving quantity leads to a non-unique relation between damage growth and the (non-local) damage driving quantity $\dot{\epsilon}^{NL}$ (i.e., the point in the space $(\dot{\epsilon}^{NL}, D)$ will not follow a unique evolution curve). Therefore, the evolution law can be considered to be unique over the domain. Based on this property, a minimization problem will be developed to propose a calibration procedure. This methodology is a priori size-invariant; thus, it should not require to carry experiments with multiple structures sizes. This method can be made more robust by incorporating data from multiple sources, allowing information from various experiments — potentially involving different structure sizes — to be combined to account for size effects.

Section 1 introduces the constitutive equations in the framework of the continuum damage mechanics, as well as non-local regularization methods. The general procedure for calibrating the internal length and identifying non-local damage models is presented in Section 2. To illustrate and analyze the method, the procedure is applied to a fictive ideal case (Section 3). Moreover, the ideal case is also employed to explore the sensitivity of the procedure to an (evolution) model error and to measurement noise. Then, in order to validate the procedure, we chose to generate reference fields by performing virtual experiments. The latters are carried out using a discrete model that explicitly represent heterogeneities and micro-cracking. Indeed, the discrete lattice models are well recognized in the literature for modeling concrete fracture (Nikolić et al., 2018) due to their simplicity, a faithful representation of the cracking mechanism and a computational efficiency to generate large data sets (Loiseau et al., 2023). The beam-particle model introduced and its parameters are identified with respect to the size and shape effects in Section 4. Finally, based on this virtual testing machine, we apply the proposed procedure to calibrate a non-local damage model with evolving internal length (Desmorat et al., 2015) (Section 5).
**Remark 1.** We use the term 'calibration procedure' to refer to the calibration of the non-local damage model and the term 'identification' to refer to the identification of the parameters of the considered discrete model.

**Remark 2.** Most data and codes used is this article are openly available in the Zenodo repository (Védrine et al., 2025) at Zenodo repository.

# 1. Local and non-local damage models

## 1.1. Local damage model

### 1.1.1. State law

A classic (local) damage model can be separated into different equations. In this study, we are only interested in the behavior of an isotropic material, described by a single isotropic damage variable D. The stress-strain ($\sigma - \epsilon$) relation for this type of models reads,

$$ \sigma = \bar{E}(D) : \epsilon, \quad \bar{E}(D) = (1 - D)E, \qquad (1) $$

where E and $\bar{E}$ are respectively the initial elasticity and the effective (damaged) elasticity tensors. The scalar value D = 0 represents the undamaged material with the initial stiffness and D = 1 represents a state of complete loss of stiffness, in which no stresses can be transferred.

### 1.1.2. Damage evolution laws

Whether damage growth is possible is decided on the basis of a yield function f,

$$ f(\epsilon) = \dot{\epsilon} - \kappa \le 0 \qquad (2) $$

where $\dot{\epsilon}$ is the Mazars equivalent strain, and $\kappa$ is the consolidation function. For $f < 0$ there can be no growth of damage and the response thus remains linear elastic, for $f = 0$ & $f' = 0$ damage growth is possible and follows an evolution law. Let us introduce $\kappa = \max_i \dot{\epsilon}$ as the maximum value up to the time t of the equivalent strain, that will be the damage-driving variable. Three evolution laws in tension are studied in their integrated form:

* The (Mazars, 1984) damage model in tension,

$$ D = g_{\text{Mazars}}(\kappa) = 1 - \frac{\kappa_0(1 - A)}{\kappa} - A \exp(-B(\kappa - \kappa_0)) \qquad (3) $$

* The modified Mazars damage model in tension (Eq. (4)), which correspond to the Mazars model with $A = 1$,

$$ D = g_{\text{Mazars\_mod.}}(\kappa) = 1 - \exp(-B(\kappa - \kappa_0)) \qquad (4) $$

* The (Feenstra, 1993) evolution law,

$$ D = g_{\text{Feenstra}}(\kappa) = 1 - \frac{\kappa_0}{\kappa} \exp(-B(\kappa - \kappa_0)) \qquad (5) $$

All these models are valid for $\kappa > \kappa_0$, otherwise D = 0. The parameter $\kappa_0$ is the damage threshold and A, B are two material damage evolution parameters.

## 1.2. Non-local damage model

The main idea of non-local damage models is to take into account the state of the neighborhood of a material point, by replacing a local variable with its non-local counterpart. Following the work of Bažant and Pijaudier-Cabot (1989) and Jirásek (1998), the damage driving variable (the equivalent strain $\dot{\epsilon}$ in our case) is replaced by a non-local variable $\dot{\epsilon}^{\text{NL}}$ (that will act as a localization limiter) in the damage criterion and the damage evolution law.

### 1.2.1. Non-local integral (NLI)

Initially proposed by Pijaudier-Cabot and Bažant (1987), the main idea of the NLI damage theory was to fix the dissipated energy by spreading damage over a non-zero and mesh-independent area. Therefore, it replaces the local field $\dot{\epsilon}(x)$ by a weighted average over the whole volume V as follows:

$$ \begin{gathered} \dot{\epsilon}^{\text{NL}}(x) = \frac{1}{V_r(x)} \int_V \Psi(x, \xi, l_c) \dot{\epsilon}(\xi) \, dV_{\xi} \\ \text{and} \quad V_r(x) = \int_V \Psi(x, \xi, l_c) \, dV_{\xi} \end{gathered} \qquad (6) $$

where $l_c$ is the characteristic length of the material and $\Psi(x, \xi, l_c)$ is a weight function, which generally follows a Gaussian type distribution: $\Psi(x, \xi, l_c) = \exp(-\frac{\|x-\xi\|^2}{2l_c^2})$. However, this type of model needs to detect each point's neighbors and evaluate the distance between them to compute the non-local field, which is computationally expensive, as an integral over the domain must be calculated at each Gauss point.

### 1.2.2. Non-local gradient (NLG)

Another approach was proposed by Peerlings et al. (1996), based on the work done by Aifantis (1987). The idea here is to get an equivalent expression of the non-local equivalent strain by replacing the equivalent strain with its expansion to a Taylor series, after additional adjustments to obtain a stable numerical scheme (for more details see Peerlings et al., 2001), we obtain the implicit gradient-type formulation (called NLG in this study):

$$ \begin{gathered} \dot{\epsilon}^{\text{NL}}(x) - \frac{1}{2} l_c^2 \nabla^2 \dot{\epsilon}^{\text{NL}}(x) = \dot{\epsilon}(x) \quad \text{on } \Omega \\ \text{and} \quad \nabla \dot{\epsilon}^{\text{NL}}.n = 0 \quad \text{on } \partial \Omega \end{gathered} \qquad (7) $$

where $\Omega$ is the domain under consideration, $\partial \Omega$ denotes its boundary, n is the outward unity vector to $\partial \Omega$, and $l_c$ is a (constant) characteristic length.

### 1.2.3. Eikonal non-local gradient (ENL-G)

The previous formulations, with constant internal length, successfully address the issue of spurious mesh dependency, and were shown to induce convergence toward meaningful solutions with a finite non-zero dissipated energy. Nevertheless, those models introduce some non-physical interactions of material points close to free boundaries or across damaged bands and holes (Geers et al., 1998). To solve this issue, non-local interactions evolving with mechanical fields (e.g., stress, strain, damage) were introduced by different authors (Giry et al., 2011; Pijaudier-Cabot et al., 2004; Desmorat et al., 2015; Poh and Sun, 2017; Vandoren and Simone, 2018) to better model strain localization, thus indirectly describing a progressive damage-to-fracture transition. In this study, we will mainly focus on the Eikonal approach, proposed by Desmorat et al. (2015), which accounts for damage-dependent non-local interactions. This model represents non-local interactions which are described by a function of the damage state. The resulting differential problem to be solved can be written as:

$$ \begin{gathered} \dot{\epsilon}^{\text{NL}} - \frac{1}{2} l_c^2 \frac{1}{\sqrt{\det g}} \nabla \cdot (\sqrt{\det g} g^{-1} \nabla \dot{\epsilon}^{\text{NL}}) = \dot{\epsilon} \quad \text{on } \Omega \\ \text{and} \quad \nabla \dot{\epsilon}^{\text{NL}}.n = 0 \quad \text{on } \partial \Omega \end{gathered} \qquad (8) $$

where g is the damage-dependent Riemannian metric: $g = (1 - D1)^{-1}$. Note that this metric is reduced to a Euclidean metric ($g = 1$) for an undamaged material, which gives back the NLG model.

**Remark 3.** Although all these models introduce a characteristic length $l_c$, its value may differ between models for the same material.

# 2. Procedure to calibrate the internal length and to identify the non-local damage models

In this section, we propose a procedure to calibrate the internal length of a non-local model and to identify the parameters of the evolution law.
## 2.1. General procedure

For an isotropic homogeneous material, the damage evolution law is a relationship between damage growth and the damage driving variable (Section 1.1.2), here $κ = max_t \hat{\epsilon}^{NL}$. This relationship should be unique for all the material points and independent with the loading, sample geometry and sample size. The proposed calibration procedure is based on this property.

It means that there is an unique curve representing the evolution of damage $D$ with the damage-driving variable $κ$. Note that this curve is independent with specimen size and geometry. The proposed calibration procedure consists in different steps. First, both the displacement field $u(x, t)$ and the damage field $D(x, t)$ are measured. Then, non-local damage model (e.g., Eqs. (6)–(8)) is applied to the measured displacement field $u(x, t)$ to compute the damage-driving quantity $κ$ with a specific choice of internal length $l_c$.

**Remark 4.** To illustrate the calibration procedure, we have chosen to calibrate non-local damage models with a single parameter: the internal length. However, this method can also be used to distinguish two different non-local damage models (and its parameter(s)) thanks to the residual model error.

The next step consists in plotting the curve of damage $D$ vs. the damage-driving variable $κ$. As the damage evolution law is supposed to be unique, the point should align along a single curve if the internal length $l_c$ is well-chosen. However, if the internal length $l_c$ is not the adequate one, the points will not follow an unique curve.

In order to make this procedure systematic and more robust (in particular to measurement noise and discrepancies), an arbitrary evolution model has to be introduced. The idea is to fit the evolution law on the damage $D$ vs. the damage-driving variable $κ$ data points. Therefore, the calibration problem becomes finding the best set of parameters (characteristic length of the non-local damage model $l_c$ and the parameters of the evolution law $p(κ_0, A, B$ in our case)) to minimize the error,

$$ \mathcal{L}(l_c, \mathbf{p}) = \left( \sum_{i \text{ in dataset}} \frac{(D(\hat{\epsilon}_i^{\text{NL}}, l_c, \mathbf{p}) - D_i)^2}{D_i^2} \right)^{1/2} \quad (9) $$

where $(\hat{\epsilon}_i^{\text{NL}}, D_i)$ are the measured data sets. Algorithm 1 details the procedure to calibrate a non-local damage model.

**Algorithm 1: Calibrate a non-local damage model**

Perform experiments (potentially on different geometries/sizes);
**foreach experiments do**

Measure the damage field (see Section 2.2);

Measure the strain field (from the displacement field);

Compute the non-local strain $\hat{\epsilon}^{\text{NL}}$ (Equations. (6)–(8)) from strain with a specific $l_c$;

Compute the damage-driving variable $κ$ from $\hat{\epsilon}^{\text{NL}}$;

**end**

Gather results and plot the evolution of $D(κ)$;

Solve the calibration problem by finding the characteristic length $l_c$ and the parameters of the evolution law that minimizes the error (Equation (9));

**Remark 5.** For proportional loading, the analysis can be restricted to the relationship between the damage variable, $D$, and the evolution driving variable, $κ$. However, in the case of non-proportional loading, the damage evolution must also take into account the rate of damage, $\hat{D}$, meaning the study should be extended to the space $(κ, D, \hat{D})$.

**Remark 6.** Note that this procedure not only provides a mean to calibrate the non-local damage model parameters, it also provides a mean to verify and to calibrate the proposed damage evolution law.

## 2.2. Generate a reference set of tests on different beam sizes

Algorithm 1 requires monitoring the evolution of damage and strain during mechanical testing. The most challenging aspect is accurately measuring the damage. Various techniques have been employed in mechanical testing; for instance, acoustic emission can provide insights into the cracking state, but linking this information to material damage remains difficult (Rossi et al., 1990; Landis, 1999). However, with the development of Digital Volume Correlation (DVC) and advanced post-processing techniques (Zhang and Zhou, 2022), determining damage from displacement fields has become more straightforward (Hild et al., 2015; Nguyen et al., 2015).

In a first step, we propose to apply the procedure to an ideal case to illustrate it and show its capability. Then, we propose to use a beam-particle model, considered as a virtual testing machine, to generate a more complex set of reference mechanical tests to apply the proposed calibration procedure (Algorithm 1). Since the measured fields are discrete, an interpolation phase is necessary.

## 3. Application of the calibration method on a 1D ideal case

This section is dedicated to the illustration and analysis of the proposed calibration method. The robustness of the calibration procedure with respect to the choice of evolution law and measurement noise is also assessed and discussed. To this aim, an multi-dimensional tensile test on a bar serves as our reference. Based on this simple test case, we show how a miscalibration of the internal length affects the damage $D$ vs damage driving variable $κ$ curve. Then, we show how the errors in this space can be used to calibrate the non-local damage model, while also analyzing the sensitivity to both model errors and the measurement set. To conclude this section, we also show how that the $(κ, D)$ is invariant under size changes, making this method ideal for incorporating various geometries and sizes into the calibration process.

### 3.1. Ideal reference problem

Let us consider a bar of length $L_b = 1$ m. The bar has a Young modulus $E = 46$ GPa. The damage evolution is governed by the Feenstra model in Eq. (5) with $κ_0 = 10^{-4}$, $A = 1$, and $B = 1/9 \times 1e4$. The gradient non-local damage model from Eq. (7) is applied with the characteristic length $l_c = l_{c,ref} = 0.04$ m. Those properties corresponds to the properties of a concrete. To ensure the localization of the damage field, a thickness reduction of 1% is applied in the out-of-plane direction on the central element. The strain and damage fields are computed from this non-local damage model and are considered to be the only known quantities to apply the calibration procedure. To use as a reference result, the damage driving variable $κ$ used in the reference to compute the damage evolution is also computed (but kept unknown in the calibration procedure). In this section, this damage-driving variable is called the reference in contrast with the test damage-driving variable (the guess).

### 3.2. Error due to a miscalibration of the internal length

We consider that the strain and damage field in the bar are both known. To illustrate the principle of the calibration procedure, for each time step and each point of the domain we

1. compute the test equivalent strain from the known strain field, since we are in uniaxial tension $\hat{\epsilon} = \epsilon$,

2. compute the test non-local equivalent strain via Eq. (7) for different values of characteristic length $l_c$ ranging from $l_{c,ref}/8 = 0.005$ m to $8 \times l_{c,ref} = 0.32$ m, and

3. compute the test damage-driving variable $κ = max_t \hat{\epsilon}^{\text{NL}}$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Damage evolution as a function of the damage-driving variable $\kappa$ where the test non-local equivalent strain is calculated with different characteristic length $l_c < l_{c,ref}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Damage evolution as a function of the damage-driving variable $\kappa$ where the non-local equivalent strain is calculated with different characteristic length $l_c > l_{c,ref}$.</figcaption></figure>

To properly describe the proposed method, let us first illustrate how a miscalibration of the internal length $l_c$ affects the damage $D$ vs. damage-driving variable $\kappa$.

Fig. 1 shows the data points $(\kappa, D)$ when the characteristic length $l_c$ is smaller than the reference characteristic length $l_{c,ref}$. The reference evolution model is also displayed along the data points. In this case, we observe vertical branches in the data points $(\kappa, D)$; thus, the data points do not align along a single curve. Those horizontal lines are typical of an under-estimation of the characteristic length $l_c$. They are due to the test strain diffusion region being too small compared to the reference. Thus, the (test) strain diffusion region is smaller than the reference damage diffusion region. As a result, the 'high' reference damage values are associated with 'small' test damage-driving variable values (as the strain diffusion region is smaller). This error appears as vertical branches in this plots.

Fig. 2 shows the same type of plot when the characteristic length is larger than the reference characteristic length. In this case, we observe horizontal branches in the data points $(D, \kappa)$; thus, the data points do not align along a single curve. Indeed, damage peaks can be observed when damage reaches values between 0.8 and 1, due to the under-estimation of the characteristic length $l_c$. Those peaks are due to points where the damage-driving variable $\kappa$ is under-estimated. They are due to the strain diffusion region being too large compared to the reference. Thus, the (test) strain diffusion region is larger than the reference damage diffusion region. As a result, some 'small' reference damage values are associated with 'large' test damage-driving variable values (as the strain diffusion region is smaller). This error appears as horizontal branches in this plots.

In this ideal case, the evolution (relation between damage $D$ and this choice of damage driving variable $\kappa$) is already known to be unique over both space and time in the bar. When the internal length $l_c$ is well-chosen, the data points should be perfectly aligned along the reference

<figure><img src="image_4.png" /><figcaption>Fig. 3. Damage evolution as a function of the damage-driving variable $\kappa$ where the non-local equivalent strain is calculated with the reference characteristic length $l_{c,ref}$.</figcaption></figure>

model. Due to the wrong choice of characteristic length, the reference evolution model does not fit perfectly to the data points in both Figs. 1 and 2. Fig. 3 depicts the case where the internal length used to compute $\kappa$ is correctly chosen. In this case, the data points are aligned along a curve with correspond to the reference model.

### 3.3. Application of the procedure

This subsection is dedicated to applying the procedure to the ideal reference data. We fit the damage evolution models in Eqs. (3)–(5) to the data points $(\kappa, D)$ shown in Figs. 1–3.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Error between the fitted model and the reference data ($D, \kappa$) for different values of the characteristic length.</figcaption></figure>

Fig. 4 shows the error $L$ (see Eq. (9)) depending on the characteristic length $l_c$. It shows that the calibration error for the Feenstra model goes to zero ($L \approx 10^{-5}$) only when the characteristic length $l_c$, used to calculate the damage-driving variable $\kappa$, is equal to the reference characteristic length $l_{c,ref}$.

### 3.3.1. Influence of model selection on residual errors in internal length estimation

It is important to note that the other models also reach their minimum for the reference internal length. Yet, those other models still display a residual error for the reference internal length. This residual error is the so-called model error, due to the specific choice of damage evolution law. This residual error could be exploited in a more generic approach. Indeed, if one had no idea of the model structure, interpretable machine learning methods could be used to reduce the model error (e.g., symbolic regression Wang et al., 2019 or sparse regression Tibshirani, 1996; Loiseau et al., 2023). The error reaches a plateau when the characteristic length is too short and should converge to the error in the local evolution model. When the characteristic length is too high, the error highly grows. Note that this error seems to be more sensitive to an overestimation of the characteristic length than to an underestimation.

**Remark 7.** Interestingly, note that the error in Fig. 4 admit only one minimum and that the error function seems to be convex in $l_c$, facilitating in minimization. A strategy to perform the minimization would be to employ a Newton-Raphson algorithm with an initial guess that over-estimate the internal length.

### 3.3.2. Influence of measurement noise on damage field and internal length estimation

The aim of this section is to assess the consequences of uncorrelated noise corrupting the measured displacement field, which results in perturbations of the damage field and the damage driving variable field. We employ the same method as Avril et al. (2004) and Amiot et al. (2007), by adding Gaussian white noise characterized by a zero mean and a specified standard deviation $\sigma_n$, represented as $N'(0, \sigma_n^2)$, to all fields $Z$ (where $Z = \kappa$ or $D$). We quantify the noise on the field $Z$ by the Noise-to-Signal Ratio (NSR$_Z = \frac{\mu_Z}{\sigma_n Z}$), where $\mu_Z$ is the field's mean value.

We study different pairs of Noise-to-Signal ratio for the pair $(\kappa, D)$ in order to quantify their impact on the identification of the internal

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of the error between the fitted model and reference data under varying levels of white noise on $(D, \kappa)$, evaluated for different values of the characteristic length.</figcaption></figure>

<caption>Table 1</caption>

Calibrated values of parameters $\kappa_0$ and $B$ as well as their errors compared to the reference.

<table><thead><tr><th>(NSR<sub>D</sub>, NSR<sub>κ</sub>)</th><th>κ<sub>0</sub></th><th>B</th><th>Error</th></tr></thead><tbody><tr><td>(0,0)</td><td>9.960 · 10<sup>-5</sup></td><td>2738</td><td>0</td></tr><tr><td>(0.01,0.01)</td><td>1.010 · 10<sup>-4</sup></td><td>2907</td><td>0.0027</td></tr><tr><td>(0.05,0.05)</td><td>1.116 · 10<sup>-4</sup></td><td>3673</td><td>0.0196</td></tr><tr><td>(0.1,0.1)</td><td>1.255 · 10<sup>-4</sup></td><td>3750</td><td>0.0433</td></tr></tbody></table>

length. The results shown in Fig. 5 illustrate the error $L$ for different levels of noise in the data set. As the Noise-to-Signal ratio approaches zero, the error for the reference internal length also approaches zero, and the identification error vanishes, thereby demonstrating the consistency of the calculations. The plot also shows that noising the damage field has a small impact on the error compared to noising the strain field: the results indicate that the error $L$ is more sensitive to noise in $\kappa$ than in $D$. This can be explained by the saturating shape of the curve $D(\kappa)$. Similarly, an excessively large internal length, as well as an error in the $\kappa$ field, leads to the diffusion of strains for intermediate damage and significant errors, whereas the diffusion of damage has a lesser effect. This effect results in an increased error for small internal lengths. Our results show that the identified internal length does not deviate if $NSR_\kappa < 10\%$ regardless of the noise level on the damage field. This threshold is relatively high compared to the tolerance level found in Amiot et al. (2007).

We also studied the sensitivity of the calibrated evolution parameters to the Gaussian noise (applied to the input fields). The calibrated values of damage parameters $\kappa_0$ and $B$ as well as their errors compared to the reference are reported in Table 1. Note that for an NSR of 10% on the fields $\kappa$ and $D$, we obtain a relative error of 26% on $\kappa_0$ and of 37% on $B$. To check the influence of this error on the evolution law, Fig. 6 compares the calibrated evolution law for the different NSR to the reference one. The mean-square error between the calibrated evolution law and the reference are also reported in Table 1. It shows that the error on the evolution is minor despite seemingly large errors on the parameters.

### 3.4. Discussion on combining different structures in the calibration procedure

As we integrate local data in the calibration procedure, the method can easily incorporate data from multiple sources. For instance, it
<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution law with a parametrization identified based on different Noise to Signal ratio on the (D, κ) field.</figcaption></figure>

would enable to combine the results from a large experimental campaign with the distinctive experiments proposed by Bažant and Nguyen (2023). Indeed, the damage D vs damage-driving variable κ curve is expected to be invariant with the geometry and size of the structure making it an ideal candidate for this type of calibration. To illustrate this claim, we perform simulations in the ideal case with different beam length $L_b \in \{0.75, 1, 2\}$ m.

**Remark 8.** Due to the well-known boundary effect of the gradient non-local model (Geers et al., 1998), smaller beam sizes could not be reached.

The spatio-temporal evolution the damage field is represented on the left part of Fig. 7. On the right part of Fig. 7, we show the damage-driving variable κ vs. damage D plot for each size. As expected, this curve is invariant with the size of the structure making it an ideal candidate for this type of calibration.

**Remark 9.** We show that different structural sizes can be considered but as we considered only 1D cases, they does not allow to conclude on the size effect.

# 4. Introduction of a virtual testing machine

To generate more realistic reference damage fields, we have chosen to use the beam-particle model developed by Delaplace (2008) and Vassaux et al. (2015, 2016). This model has already been used to: investigate the formation and impact of drying shrinkage crack (Chan et al., 2018), correctly reproduce cracking patterns in concrete (Oliver-Leblond, 2019), model an anchor system (Rima et al., 2021), study the effect of micro-cracking on elasticity (Oliver-Leblond et al., 2021) and model it with an anisotropic damage model (Loiseau et al., 2023). However, only few tests have been performed to reproduce size effect (Giry et al., 2014; Vassaux et al., 2016).

Recent studies have shown the ability of lattice model to represent size effects (Grassl and Bažant, 2009; Grassl et al., 2012; Murali and Deb, 2018; Suchorzewski et al., 2019), a major criterion in the validation of numerical damage modeling. These capacities are linked to the fulfillment of two conditions: the introduction of heterogeneities through the meso-structure and its mechanical properties, and the introduction of a characteristic size (Grassl and Bažant, 2009). Once the discrete model is calibrated for scale effects, it becomes a valuable tool for numerically investigating the size effect, such as its impact on

fracture process zones in concrete structures (Grassl et al., 2012). In this part, we will verify that the beam-particle program (Discrete Element Analysis Program) can represent the size effect on the fracture process in quasi-brittle materials and identify the material parameters based on the (Grégoire et al., 2013)'s experiments campaigns.

## 4.1. Beam-particles model

In the beam-particle model used, the medium is considered as a random assembly of polygonal particles based on the procedure proposed by Moukarzel and Herrmann (1992). The model can be described as follows.

**Geometry.** Firstly, a regular grid is defined with a constant spacing $l_b$, which corresponds to the average particle size. Then, a node is allocated randomly into each grid cell, these nodes correspond to the centroids of the particles. Finally, Delaunay's triangulation of the set of nodes is performed. The edges of the Delaunay triangles represent the cohesive beams linking neighboring particles, and the polygonal particles are computed from Voronoi's tessellation.

**Elasticity.** The elastic Euler-Bernoulli beams connecting the particle nodes ensure cohesive interactions between the particles. Their elastic behavior is defined by a Young modulus $E_b$ and a coefficient of inertia $\alpha = (64I_b\pi)/A_b^2$, with $I_b$ and $A_b$ the second moment of area and the section of the beam. Note that, these elastic quantities are assumed to be the same for all the beams.

**Fracture.** The quasi-brittle fracture behavior is obtained by considering a perfectly brittle law for the beams with different breaking thresholds to reproduce the spatially variable tensile strength and account for the material heterogeneity. The failure criteria of beam ij linking two particles i and j is

$$P_{ij} = \frac{\epsilon_{ij}}{\epsilon_{ij}^{cr}} + \frac{|\theta_i - \theta_j|}{\theta_{ij}^{cr}} > 1, \quad (10)$$

where $\epsilon_{ij} = \frac{\|u_i - u_j\|}{l_{ij}}$ is the beam strain, $|\theta_i - \theta_j|$ is the beam curvature, and $\epsilon_{ij}^{cr}$ and $\theta_{ij}^{cr}$ are respectively the breaking thresholds in tension and rotation. As proposed by Oliver-Leblond (2019), we assumed that the breaking threshold follows a three-parameter Weibull distribution,

$$f(x) = \frac{k}{\lambda} \left(\frac{x}{\lambda} - m\right)^{k-1} e^{\left(-\frac{x}{\lambda} - m\right)^k} \quad (11)$$

with $\lambda$ the scale factor, k the shape factor and m the position factor. The position factor corresponds to the ratio between the minimum resistance in tension and the scale factor.

**Remark 10.** Note that the friction between the crack faces induces irreversible strains (Andrieux et al., 1986) which is modeled in Vassaux et al. (2015, 2016). However, to ease the modeling process, the crack closure effects, such as the permanent strain at unloading, will be omitted for now and left to future work.

## 4.2. Identification of beam-particle model parameters with respect to size and shape effects

### 4.2.1. Identification procedure

As Murali and Deb (2018) and Suchorzewski et al. (2019), we use four 3-point bending tests to verify that the beam-particle model can represent the size effect on the fracture process in quasi-brittle materials. We compare the numerical results with Grégoire et al. (2013) experiments for four specimen sizes and three geometries (illustrated in Fig. 8): unnotched (UN), fifth notched (FN) and half notched (HN). To the best of our knowledge, this is one of the most important experimental campaigns to study size effect with Hoover et al. (2013). Four sizes of geometrically similar prismatic specimens with four depths $D = 50$ mm, 100 mm, 200 mm, and 400 mm are tested. In addition, different shapes are explored, a notch-to-depth ratios $\alpha = a/D$ of 0, 0.2 and 0.5 are tested, as shown on Fig. 8a. A crack width of 2 mm is set for all notched beams. The beams, of length 3.5D, has a span-to-depth ratio $S/D$ of 2.5,
<figure><img src="image_8.png" /><figcaption>Fig. 7. Illustration of the invariance with the structure size. The left part shows the spatio-temporal evolution of damage for each beam length $L_b$. The right part shows the damage-driving variable $\kappa$ vs. damage $D$.</figcaption></figure>

and the out-of-plane thickness is 50 mm. The samples are tested under Crack Mouth Opening Displacement (CMOD) control to obtain a stable post-peak response. The resolution method used in the beam-particle model consists of progressively removing the beams from the network when they are considered to be broken, and then updating the general stiffness matrix of the system to solve a new elasticity problem. So even if the numerical experiments are not controlled with the CMOD, the snap-back can be correctly described. The CMOD is measured with an extensometer attached to the bottom fiber of the beam at a distance from the mid-span of half the beam depth, to ensure crack initiation occurred between the two legs of the extensometer.

The boundary conditions and geometric parameters for the simulation, shown on Fig. 8b, are chosen to correspond to the experimental conditions (Grégoire et al., 2013). The simulations are performed under displacement control conditions and the CMOD is computed by calculating the displacement of two points separated by a distance equal to one beam depth $D$, as in Grégoire et al. (2013). To take into account the heterogeneity of the behavior, the displacement is averaged over a square of 4 particles's width.

### 4.2.2. Identification of beam-particle model parameters

Mesh size. First, we set the mesh size by imposing an average particle size of $l_b = 2$ mm to ensure both a controlled computational cost and convergence of elastic properties, as well as a correct crack pattern. This way, there are at least 25 particles along the smallest side of the

<figure><img src="image_9.png" /><figcaption>Fig. 8. (a) Geometry tested. (b) Numerical setup and two-dimensional spatial discretization of the specimens. Half-notched (HN) beam of size depth $D = 100$ mm.</figcaption></figure>

smallest specimen. In addition, the largest simulation, consisting of 35000 particles, is performed in almost 5 hours.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Comparison of simulated and experimental load-CMOD curves for different geometry. The gray zone corresponds to experimental variability, and the cyan/blue dashed line corresponds to the specimen used for identification.</figcaption></figure>

<table><caption>**Table 2**
Parametrization of the beam-particle model.</caption><thead><tr><td>Parameter</td><td>Values</td></tr></thead><tbody><tr><td>Young modulus <i>E</i> (GPa)</td><td>46</td></tr><tr><td>Inertia coefficient <i>α</i></td><td>0.83</td></tr><tr><td>Average particle size <i>l</i><sub>b</sub> (mm)</td><td>2</td></tr><tr><td>Shape factor <i>k</i></td><td>0.5</td></tr><tr><td>Position factor <i>m</i></td><td>0.3</td></tr><tr><td>Extension scale factor <i>λ</i><sub>c</sub></td><td>4.494 × 10<sup>-4</sup></td></tr></tbody></table>

**Material parameters.** Firstly, for the elastic phase, we use the elastic parameters proposed by Vassaux (2015) for the same experiments, which correctly represent the elastic phase. A beam Young modulus of *E* = 46 GPa and an inertia coefficient of *α* = 0.83 are used for the rest of the study, reproducing the elastic properties of the concrete (*E*<sub>m</sub> = 37 GPa and *ν* = 0.21). As proposed by Oliver-Leblond (2019), the failure thresholds in extension and in rotation are supposed to be correlated, therefore the shape factor is shared for both thresholds and a ratio of 0.073 is fixed between the extension and rotation scale parameters.

$$k_c = k_{\theta} = k \quad m_c = m_{\theta} = m. \qquad (12)$$

In practice, a threshold are obtained in one draw of the Weibull distribution and scaled for extension and rotation. It induces a spatial correlation between both failure thresholds in extension and rotation. The three parameters related to the tensile failure threshold (i.e., $\lambda_c, k, m$) are calibrated based on two statistical indicators: the relative error on peak loads and the Mean Absolute Percentage Error (MAPE), which corresponds to the relative error made in the local energy calculation (for more details, see Appendix A). The identification of the fracture parameters is performed on a single specimen, an unnotched beam with depth *D* = 100 mm. The experiments are repeated 32 times and at each time step the mean force and CMOD are calculated. The impact of parameters *k* and *m* to reproduce the scale and shape effects are thoroughly investigated. Building on the work of Oliver-Leblond (2019), we selected $k \in \{0.5, 1, 1.5, 2, 2.5, 3\}$ and $m \in \{0, 0.1, 0.2, 0.3, 0.4, 0.5\}$ as our parameter ranges. For each combination of these parameters, the peak error is calculated and the shape parameter in extension ($\lambda_c$) is determined to minimize the peak error. Then, the optimal pair ($k, m$) is then chosen to minimize the MAPE error.

### 4.2.3. Representation of size effect with the beam-particle model

The procedure proposed in Section 4.2.2 is followed, and we retain one set of material parameters: $k = 0.5, m = 0.3$ and $\lambda_c = 4.494 \times 10^{-4}$ (Table 2).

The model is then used to predict the load-CMOD curves for all sizes and notch lengths of the Grégoire et al. (2013)'s experiments, based on the chosen material parameters (i.e., without adjusting the model parameters). Fig. 9 shows the predicted force-CMOD curves for all the sizes, all the notch-to-depth ratios. The gray area represents the experimental results range, and the solid line is the mean response of the Monte Carlo simulations. To achieve the statistical convergence, respectively 32 and 16 simulations are required for the specimen of size 50 mm, 100 mm, and 200 mm, 400 mm (for more details see Appendix A.2). Firstly, it can be observed that the numerical simulations predict well the mechanical behavior in the elastic, near-peak, and post-peak regimes for the different geometries and sizes and all parameters. We can note that the behavior of the two unnotched beams with the largest sizes can be reproduced, contrary to the recent study of Pathirage et al. (2023), since our model can take into account the snapback.

The model identified on a single specimen (100 mm depth unnotched beam) is then assessed based on blind predictions for different beam sizes and notch lengths. Fig. 10(a) shows the scatter between predicted and experimental peak forces for all geometries and sizes. The cloud of points is close to the 1:1 line. An average absolute relative error on peak load of 16.38% is obtained in this study using a coarse identification procedure with only 7 material parameters. Incorporating additional phenomena into the model may help to minimize its error. As demonstrated by Pathirage et al. (2023), a more sophisticated model was able to reduce the error to 9.69%. This is close to the error of 10.50%, achieved by Grassl et al. (2012). The load points associated with the smallest half-notched specimen, and therefore the smallest peak forces, exhibit the greatest deviation from the experimental results. The reason for this could be that the model introduces a minimum fracture strength, which increases the fracture strength for small notched specimens that have a low fracture strength. In addition, the number of particles of the smallest specimen may not be sufficient, making the model too coarse (Oliver-Leblond, 2019). The quality of the prediction can also be assessed by looking at the evolution of the peak force with beam depth (Fig. 10(b)). Based on the experimental campaigns of Grégoire et al. (2013) and compared with the numerical study of Pathirage et al. (2023), we can observe a good ability of our model to predict the peak forces of large unnotched beams. This ability could be related to the fact that the model has been calibrated for an unnotched beam and to the randomness in the material properties introduced in our model, which allows for representing the statistical size effect.

## 5. Application of the proposed calibration procedure to virtual tests

In this section, we propose to use Algorithm 1 to calibrate the ENL-G model detailed in Section 1.2.3. We use the beam-particle model, with
<figure><img src="image_11.png" /><figcaption>(a) Scatter of the predicted versus measured values of the peak loads. The dashed lines correspond to the 1:1 line.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Peak loads versus size for experimental data and simulation results from Pathirage et al. (2023), Grassl et al. (2012), as well as our simulation results. The red dots in Pathirage et al. (2023) correspond to the maximum load reached before the simulation stopped due to snapback.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Methodology to generate reference test: (a) cracking test, (b) measurement test with an elastic loading, (c) effective discrete young modulus field, (d) discrete damage field.</figcaption></figure>

the material parameter identified with respect to the size and shape effects (Section 4), as a virtual testing machine to generate reference damage and strain fields.

## 5.1. Creation of a database of reference damage field

### 5.1.1. Virtual testing setup

The beam-particle model is used as a virtual testing machine to generate damage field based on the explicit micro-cracking of the material. Three beam sizes have been studied, with a constant height of 32 mm (16 particles) and three different lengths $L$: 256 mm, 512 mm, 1024 mm, which correspond respectively to 128, 256 and 512 particles. The beams are loaded with a 1D uniaxial tensile test until failure (Fig. 11(a)). The loading are repeated for 16 different particles distributions for each size to represent the statistical variability of the material.

### 5.1.2. Damage measurement and post process

Let us now describe how to determine the strain and damage fields during a virtual test. The proposed method is described in Fig. 11. Since we are using a discrete tool, we must model our 1D problem in 2D. Therefore, the challenge in this section is to transition from discrete 2D

measurements to continuous 1D measurements. After performing the virtual test (see Fig. 11(a)), we define different observation windows as shown in Fig. 11(b). These windows are used to compute the average strains and damages along the $x$ direction. Given the orientation of micro-cracks, we assume that the damage mainly affects the elastic properties in the loading direction, allowing us to consider the problem as uni-dimensional. Knowing both strain and damage in each window gives a discrete field for both mechanical quantities. For every cracking state (i.e., at each step during which a new beam breaks) during the virtual test, the average strain is measured in each window by measuring the average displacement along the $x$ direction (in a wide band of two particles) on each side of the window. Similarly, for every cracking state, the effective Young's modulus $\bar{E}$ is measured using an uni-axial elastic tensile test as shown in Fig. 11(bc). Knowing the initial stiffness $E_0$ of this observation window, the associated damage can be calculated as $D(x) = 1 - \bar{E}(x)/E_0$. An example of damage field is provided in Fig. 11(d).

**Remark 11.** Following the sensitivity analysis in Appendix B, the size of the observation windows is set to 8 particles along the $x$ direction and 16 particles in height (the all height). There is no overlap of observation windows. Appendix B.

Thus, for each observation window and each state of damage, we have access to the average value per window of the longitudinal strain and the damage. However, a transition from these discrete data to a continuous field is necessary to calculate the non-local strains. For the sake of simplicity, and as we have no information on the typical damage field in a beam during a tensile test, we have chosen to use a piecewise linear interpolation, as shown in Fig. 12.

## 5.2. Application of the calibration methodology

We now propose to apply the proposed calibration procedure (Algorithm 1) to the ENL-G nonlocal damage model (described in Eq. (8)). The different steps of this application are:

(1) Calculate and interpolate damage fields and local strain for different geometries, which is explained in Section 5.1.2 using the beam-particle (discrete) model.

(2) Compute the non-local damage driving quantity, $\kappa = \max_t \bar{e}^{NL}$.

(3) Calibrate the ENL-G non-local model and estimate the parameters of the evolution law within the framework of continuous non-local damage mechanics. Different characteristic lengths
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Mechanical fields in a beam 128 particles long after the full failure of the macroscopic beam (Geometry Fig. 11a).</figcaption></figure>

$l_c/l_b$ are explored, where $l_c$ is the characteristic lengths of the ENL-G model and $l_b$ is the mean beam length of the discrete model, used for the non-dimensionalization of the problem.

### 5.2.1. Non-local interactions and non-local strains induce by the ENL-G approach

The non-local strain is calculated by solving the PDE in Eq. (8). The only unknown parameter of the Eikonal equation to compute the non-local strain is its internal length and potentially the expression of the metric.

**Remark 12.** In a more generic approach, the metric in Eq. (8) could also be considered as an unknown in the calibration problem.

To understand and illustrate the calibration, we propose to calculate and analyze the weighting function. We have shown that the Green's function (G) solution of the Eikonal boundary value problem corresponds to the weighting function for calculating the non-local variable (Eq. (13)). Details of the equivalence and the finite-element approach to solving the problem are given in Appendix C.

$$ \hat{\epsilon}^{\text{NL}}(x) = \int_{\Omega} G(x, y)\epsilon(y)d\Omega \qquad (13) $$

Based on the damage field, it is possible to calculate the Green's functions associated with the boundary value problem. The function obtained for a correlation length $l_c = 5$ and a fully cracked beam is shown in Fig. 11(a). This function can be interpreted as a set of weighting functions (Fig. 13(ab)). For a fixed Dirac perturbation

at coordinate y, we obtain the weighting function used to calculate the non-local variable at coordinate x. Fig. 13(c) shows the weight functions at three different locations, away from the damage peak ●, and on each side near the peak ●●. We can clearly see that damage cuts off the interactions as the shape functions becoming asymmetric. Next, the effect of damage and internal length will be studied by 2D representations of Green's function. Also note that the Green function G depends on the damage field $D(x)$ through the metric g.

The Green functions are plotted in Fig. 14 for given geometries at different crack states: uncracked, intermediate cracking and fully cracked; and for different internal lengths. A dashed cross indicates the damage peak. Without damage, we can see the effect of internal length on the interactions. The greater the internal length, the further the non-local interactions. As the damage increases, two phenomena can be observed. First, it acts as a break in the interactions, and the points on each side of the damage peak no longer interact. Indeed, we observe a break in the interactions in quadrants 2 and 4, in other words, the parts to the left and right of the damage peak do not interact with each other. On the other hand, quadrants 1 and 3, which represent respectively the interactions between the points to the right and left of the damage peak, are less affected. Secondly, we can see that as the damage increases, the Green's functions tend towards a Dirac at the level of the interaction peak. This relocalization phenomenon is desired in the Eikonal approach, by introducing a singularity in the metric for damage of one. It makes the damage tends to 1 for infinitely thin ban, i.e. a crack.

Using Green's functions, it is possible to calculate non-local strains as a function of local strains (Eq. (13)). The integral is calculated solely on the basis of observation points, in order to preserve the discrete approach of the methodology used in this study and not to depend too heavily on interpolation choices. The non-local strains are plotted for different characteristic lengths for a cracking state of 57% (i.e. the ratio of the micro-beam of the beam-particle model broken) and a geometry given in Figure 15. It can be seen that a long internal length spreads the strains. In addition, it can be seen that for a short interaction length, the interactions are cut off close to the damage peak and non-local strain values increase at the damage peak.

### 5.2.2. Calibration of the internal length and identification of the damage evolution law

For each characteristic length tested, we calculate the non-local strain and plot the evolution of damage as a function of the driving variable $\kappa = \max_t \hat{\epsilon}^{\text{NL}}$ (Fig. 16). As expected, damage is only concentrated in one part of the bar, while the other parts are only slightly damaged by diffuse micro-cracking.

**Remark 13.** Note that the point set generated by this virtual experiment is unbalanced, as observed in Fig. 16, with low damage being more prevalent than high damage. This imbalance will likely exist in real full-field measurements, as only a part of the specimen will be fully damaged. To address this imbalance, under-sampling or over-sampling techniques could be applied (Spelmen and Porkodi, 2018) but is out of the scope of this study. Moreover, the choice of loss function could also be discussed to mitigate this issue.

For a short internal length, $l_c/l_b = 2$, represented in Fig. 16(a), we observe that the damage starts to grow after a threshold on the driving variable. Then, it increases with a high degree of variability, especially for high damage values. For higher internal lengths, such as $l_c/l_b = 100$ represented in Fig. 16(b), the curve ($\kappa, D$) is almost vertical with low variability, the strains have been averaged around the mean value of strains in the beams.

Let us now compare how the different damage evolution laws introduced in Section 1.1.2 fit the data points. A least square regression has been performed to calibrate the parameters of the evolution models. The error $L(l_c, p)$, proposed in Eq. (9), between the calibrated evolution models and the data sets are plotted in Fig. 17. We observe that the
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Illustration of the link between Green function and weight function for a characteristic length $l_c^{1/2} = 5$. (a) Green function associated to a fracture beam show on the top, (b) Green function show as surface plot, the gray lines corresponds to the damage peak, (c) Weight functions for different position (y).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Green's function for the left-hand side of the beam, for different characteristic lengths and different damage states. The corresponding beam is shown in the upper part of the Figure.</figcaption></figure>

Mazars model provides the best representation of the damage evolution law obtained in this numerical experiment for all characteristic lengths. We also note that the error is lowest for a characteristic length of 2 $l_b$,

suggesting that this is the best fitting internal length for our non-local damage model. This result is consistent with the fact that the mesh size of the discrete model introduces a characteristic length.
<figure><img src="image_18.png" /><figcaption>Fig. 15. Local and non-local strain for different characteristic lengths. Specimen shown in Fig. 11, at 57% of final crack stage ($l_b$ = 2 mm).</figcaption></figure>

**Remark 14.** Note that several measurement points are required in the process zone. However, as in experimental studies, only a few measurement points are accessible with the virtual testing machine. In the case of discrete methods, introducing a spatially correlated field could be foreseen. This approach removes the mesh-dependency of the behavior. With spatially correlated parameters fields, the internal length of the material is proportional to the correlation length (Grassl and Bažant, 2009; Eliáš et al., 2015; Zhang et al., 2018). It would enable to refine the mesh without changing the internal length, leading to a better description of the process zone. To apply the procedure in real experiments, we suppose that the resolution of the field measurement must be fine enough to capture the process zone. Thus, the spatial resolution should be several time smaller than the internal length.

Interestingly, solving the minimization problem provides access to the parameters of the damage evolution law. For a characteristic length of $2 l_b$, the identified parameters of the Mazars model are $\kappa_0 = 2.5 \times 10^{-4}$, $A = 4.2 \times 10^{-1}$ and $B = 2.2 \times 10^2$.

Looking more closely at the evolution of the damage as a function of $\kappa$, we observe little variability in the non-local strains for the calibrated length (Fig. 16). Nevertheless, for an overestimated internal length the non-local strains are significantly more diffuse for low damage states and underestimated for high damage states, associated with an attenuation and spreading of Green's functions. It is interesting to note that the same results were observed for all sizes tested. It is important to emphasize that discrete approaches represent the behavior on average and introduce the statistical variability of the material; as observed in experiments. Therefore, part of the dispersion of the damage-$\kappa$ data is intrinsically related to the variability of the discrete model, and a deterministic model will most likely not be able to model this dispersion.

**Remark 15.** It is important to note that the errors observed in the real case are greater than the model error obtained in the ideal case without any errors (see Section 3.3). Indeed, the combination of both model error and measurement noise and lead to larger residual error. In the context of the virtual testing machine, measurement error manifests as variability in the material properties, which introduces a systemic error, as well as the choice to interpolate the discrete fields.

Fig. 18 plots the error $L(l_c, p)$ (Eq. (9)) calculated separately for each beam sizes $L$ (256 mm, 512 mm, and 1024 mm). Interestingly, we observe

<figure><img src="image_19.png" /><figcaption>Fig. 16. Evolution of the damage with respect to $\kappa$ in a beam submitted to uniaxial tension, obtained for 3 different beam sizes and 16 random draws of the geometry for each size, for two different non-local damage models.</figcaption></figure>

that the same minimum is reached at $l_c = 2 l_b$, whatever the beam size. However, it should be noted that a single size is insufficient to confirm the accurate representation of the size effect. Indeed, multiple non-local models with different size effects may fit a single size equally well (Havlášek et al., 2016; Zhang et al., 2021). Since these models exhibit varying size effects, they may fail to adequately represent other sizes. Therefore, incorporating multiple structure sizes into the data set is essential for ensuring robust modeling of the size effect. Once the appropriate non-local model has been identified, however, only one structure size might be sufficient, as the size effect is inherently embedded in the non-local model.

## 6. Discussion and conclusion

In this study, we propose an advanced method to calibrate a non-local damage model coupled with an evolution law. This method has been applied to an ideal case for illustration purposes, and to virtual
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Evolution of the error (Eq. (9)) with the characteristic length for different evolution laws (Eqs. (3)–(5)).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Evolution of the error (Eq. (9)) with the characteristic for the Mazars evolution law, effect of the beam size considered.</figcaption></figure>

tests using a discrete model to assess its feasibility in a more realistic setting.

We initially examined an ideal case. Using local fields extracted from a non-local 1D FEM gradient model. We sought the characteristic length introduced into the gradient model with the proposed calibration procedure. This analysis emphasized the sensitivity of the procedure to the characteristic length. Interestingly, the error is greater when the characteristic length is overestimated than when it is underestimated. Therefore, it seems preferable to start with a large characteristic length and then decrease this value until a plateau is reached in terms of error. Additionally, the sensitivity of both the choice of damage evolution model and measurement noise has been evaluated in the calibration process. Interestingly, the calibration procedure demonstrates robustness against incorrect choices of evolution model and white noise with a variation coefficient of up to 5%. Furthermore, a study has been proposed to identify errors in the evolution model. Moreover, it shows that a wrong choice of evolution law leads to a systematic error that could be exploited to figure out the evolution law expression.

In order to test the calibration procedure in a more realistic settings, we used the beam-particle model as a virtual testing machine. The first

step was to identify the material parameters of the beam-particle model with respect to the size and shape on the basis of the experimental campaign of Grégoire et al. (2013). We observed that a randomly drawn geometry and a random distribution of the failure parameters with no spatial correlation are sufficient to represent size effects. In addition, the performance of the tested model is compared with the recent (Pathirage et al., 2023) numerical identification, suggesting a similar predictive power, with the ability to reproduce all geometries, in particular unnotched geometries.

Before applying the calibration procedure to the ENL-G model using data from the virtual testing machine, an analysis of the ENL-G model have been proposed. An analysis of the effect of the characteristic length on the Green's function shows (1) how the characteristic length affects non-local interactions and (2) how damage cuts non-local interactions. Its effect on non-local strains has also been illustrated. After this analysis, the calibration procedure has been applied to the ENL-G model using the fields obtained from the discrete model. Three different evolution law have been considered. In this case, the Mazars evolution gave the best fit. Concerning the size effect, which is discussed by Guinea et al. (1992) and Kosin et al. (2024), we considered only 1D cases which does not allow to conclude on the size effect. As this method can easily gather results from various experiments, it should be possible include tests for different sizes. The damage evolution law should still be the same, then the damage *D* vs damage-driving variable *κ* curve should still superimpose for different sizes. To summarize, the examination of both the ideal case (Section 3) and the real case (Section 5) demonstrated the capabilities of the method to identify the characteristic length and the evolution law of a non-local model, the error *L*(*l*<sub>c</sub>, *p*) (Eq. (9)) being highly sensitive to these two parameters.

This study proposes a methodology for calibrating a non-local model based on full-field measurement which enables the integration of local quantities into the calibration process. Incorporating local quantities, that captured nonlocality, is essential to ensure that the problem is well-posed (Iacono et al., 2006), and dissipated energy has often been proposed as a key parameter (Xenos et al., 2015; Grégoire et al., 2015). Our approach builds on advances in displacement field measurements, using not only the measurement of the Fracture Process Zone (FPZ) (Skarżyński et al., 2011; Doll et al., 2017) but the full displacement fields in the calibration process. This comprehensive integration allows for the identification of the internal length and the evolution law. Therefore, this work offer new perspectives for future research. One potential avenue for further investigation is the application of this procedure to real experiments and different non-local damage models. In this work, we used a virtual testing machine, based on the beam-particle model, to generate realistic damage fields. In practice, Digital Image Correlation (DIC) or Digital Volume Correlation (DVC) could be employed to access full-field data during an experiment (Hild et al., 2015; Nguyen et al., 2015). To go even further, an integrated DIC approach could even be developed based on this procedure. Therefore, it would be valuable to apply this calibration methodology to various non-local damage models within an experimental campaign, similar to those conducted by Grégoire et al. (2013) and Hoover et al. (2013), adding displacements tracked using image correlation. This approach would allow for verifying the ability of non-local damage models to accurately reproduce size effects. Another perspective is to determine the evolution law expression in the calibration problem (through symbolic regression Wang et al., 2019 or sparse regression Tibshirani, 1996; Loiseau et al., 2023), allowing to mitigate the model error induced by the arbitrary choice of evolution law.

Despite the assumptions and limitations of this study, it highlights the potential of the proposed procedure to calibrate a non-local damage model. Indeed, this method is highly sensitive to the characteristic length, enabling a straightforward calibration. The result of this work can be applied to the calibration of the parameters of other non-local damage models and also to determine the damage evolution law.
### CRediT authorship contribution statement

**Louis Védrine:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **Flavien Loiseau:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Cécile Oliver-Leblond:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software. **Rodrigue Desmorat:** Writing – review & editing, Writing – original draft, Validation, Supervision.

### Declaration of competing interest

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

### Acknowledgment

We acknowledge Breno Ribeiro Nogueira from LMPS for insightful discussions on the topic of non-local damage approach. The authors gratefully acknowledge the anonymous reviewers for their comments and advice, which helped improve the quality of the final paper.

### Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

## Appendix A. Identification of beam-particle model parameters

### A.1. Statistical indicators

Two statistical indicators were used to calibrate the parameters of the discrete-element model. For a given test j (i.e. a given size and shape) the experimental results are averaged to obtain the mean data $f^j$ and the Monte Carlo simulations are averaged to obtain $F^j$. The CMOD (the predictor variable) $CMOD^j$ is discretized into equally spaced intervals $CMOD_i^j$, and $f_i^j$ and $F_i^j$ are the linearly interpolated value of $f_i^j$ and $F_i^j$ on $CMOD_i^j$, where i is the point index and n is the number of points.

Our first indicator is the relative error on peak loads ($\delta F^j$) (Eq. (A.1)).

$$ \delta F^j = \frac{\max_i(f_i^j) - \max_i(F_i^j)}{\max_i(f_i^j)} \qquad (A.1) $$

The second indicator is the Mean Absolute Percentage Error (MAPE), it is interesting to note that this error corresponds to the relative error made in the instantaneous energy calculation (Eq. (A.2)). The advantage of this error is that it allows to take into account an energetic criterion to identify the parameters of the model and its local nature allows avoiding the phenomenon of energy balance/compensation as absorbed on the total dissipated energy criterion (Thierry, 2021).

$$ \text{MAPE}^j = \frac{1}{n} \sum_{i=1}^{n} \frac{|f_i^j - F_i^j|}{f_i^j} \qquad (A.2) $$

Fig. A.19(a) shows the experimental force-crack mouth opening displacement curve and the simulated mean response for the initial parameters of the identification phase. The peak error (Eq. (A.1)) and MAPE error (Eq. (A.2)) are defined in this Figure. Fig. A.19(b) shows the experimental curve and the simulated mean response, together with the associated errors. Note that we obtain a small peak error (1.6%), and a MAPE error of 5.5%.

<figure><img src="image_22.png" /><figcaption>(a) Initial material paramters.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Calibrated material paramters.</figcaption></figure>

<caption>Fig. A.19. Simulated and experimental load-CMOD curves for the calibrated geometry. The gray zone corresponds to experimental variability.</caption>

### A.2. Monte Carlo simulations

Fig. A.20(a) underlines the evolution of the mean peak stress over the peak stress for 128 simulations as the function of the simulation number, we can note that all the specimen size, the mean peak force converge quickly. In addition, Fig. A.20(b) shows the evolution of the coefficient of variation for the different specimens size, one can observe that the smaller is the specimen the larger is the statistical variability. To estimate the number of simulations required to achieve the statistical convergence, we use a criteria based on the mean peak stress (relative error lower than 3% compared to the overkill solution), we obtain that 32 simulations are required for the specimen of size 50 mm and 100 mm and 16 simulations are required for the specimen of size 200 mm and 400 mm.

## Appendix B. Choice of observation window parameters to measure discrete damage field

A sensitivity study of the damage measurement technique has been conducted on the beam with a size of L = 256 mm. A sensitivity analysis of the damage measurement technique has been carried out. The first idea is to refine the size of the observation window, to increase the number of observation points. Fig. B.21(a) shows the damage field in the beam for different observation window width. A window width $l_W$ of 8 particles allows a much more accurate description of the damage field than when coarser windows are chosen. It is important to note
<figure><img src="image_24.png" /><figcaption>(a) Mean peak force over mean peak force obtained for 128 simulations versus sample number.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b) Coefficient of variation versus sample number.</figcaption></figure>

<caption>Fig. A.20. Effect of the sample number on the numerical response for the different sizes of the unnotched specimens.</caption>

that it is not possible to carry out tests on smaller volumes, as it is difficult to apply the boundary conditions due in particular to the rotation of the particles and the short length of the sample. For the selected window width, the assumption that the damage mainly affects the elastic properties in the loading direction is reasonable, allowing us to consider the problem as uni-dimensional.

A moving average measurement type has been also studied. Fig. B.21(b) shows, for a given observation window size, the damage field for different total numbers of windows $N_w$. For example, the first corresponds to non-overlapping windows, while the last corresponds to windows successively shifted by a $\frac{l_w}{8}$ step. We can note that the fact of carrying out moving averages contributes to artificially diffuse the damage, because a localization zone will be counted several times. Consequently, this solution will not be adopted, and we have chosen non-overlapping observation windows of eight particles wide.

<figure><img src="image_26.png" /><figcaption>(a) Damage field as a function of observation window size $l_w$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) Damage field for different total numbers of observation windows $N_w$.</figcaption></figure>

<caption>Fig. B.21. Sensitivity analysis on the choice of observation windows.</caption>

## Appendix C. Equivalence between the Green's function of the Eikonal equation and the weight function, numerical computation and approximation

This section details the theoretical approach considered in this study. All calculations and results are carried out in dimension one to facilitate their formulation and visualization, but the whole approach can easily be applied on a higher dimension using the thermodynamic derivation of Ribeiro Nogueira et al. (2023) based on differential geometry concepts.

### C.1. Boundary value problem governing the Green's function

Let us define the Boundary Value Problem associated to the ENL-G approach (Desmorat et al., 2015):

$$ \hat{\epsilon}^{\text{NL}} - \frac{1}{2} l_c^2 \frac{1}{\sqrt{\det g(x)}} \frac{d}{dx} \left( \sqrt{\det g(x)} g^{-1}(x) \frac{d}{dx} \hat{\epsilon}^{\text{NL}} \right) = \hat{\epsilon} \quad \text{on } \Omega \quad (C.1) $$

with the boundary conditions:

$$ g^{-1}(x) \nabla \hat{\epsilon}^{\text{NL}}.n = 0 \quad \text{on } \partial \Omega \quad (C.2) $$
where '∇' denotes the divergence operator and n is the outward unity vector to δΩ. To replace our problem in the 3D case, the metric g is defined as:

$$g = \begin{bmatrix} g(x) & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \qquad (C.3)$$

Any arbitrary expression of the metric could be used, for instance:

$$g(x) = (1 - d(x))^{-1} \qquad (C.4)$$

$$g(x) = (1 - d(x))^{-2} \qquad (C.5)$$

Let us simplify the equation

$$\hat{\epsilon}^{\text{NL}} - \frac{1}{2} l_c^2 \frac{1}{\sqrt{g(x)}} \frac{d}{dx} \left( \frac{1}{\sqrt{g(x)}} \frac{d}{dx} \hat{\epsilon}^{\text{NL}}(x) \right) = \hat{\epsilon}(x) \qquad (C.6)$$

with the boundary conditions

$$\left\{ \begin{array}{l} g^{-1}(0) \frac{dG(x,y)}{dx} \Big|_{x=0} = 0 \\ g^{-1}(L) \frac{dG(x,y)}{dx} \Big|_{x=L} = 0 \end{array} \right. \qquad (C.7)$$

To obtain the Green's function associated with this problem, the right-hand side must be set to a Dirac function representing a perturbation in y. In addition, the equation is multiplied by $\sqrt{g(x)}$ to simplify the calculation of the weak formulation, and we obtain a problem similar to the Sturm-Liouville problem (Bonnet, 1999), which leads to:

$$\sqrt{g(x)} G(x, y) - \frac{1}{2} l_c^2 \frac{d}{dx} \left( \frac{1}{\sqrt{g(x)}} \frac{d}{dx} G(x, y) \right) = \sqrt{g(x)} \delta(x - y) \qquad (C.8)$$

with the boundary conditions (). Note that the Green's function is also symmetric: $G(x, y) = G(y, x)$. There are also two conditions:

* Continuity of the Green's function:

$$G(y^+, y) - G(y^-, y) = 0 \qquad (C.9)$$

* Jump in the derivative of the Green's function. To obtain it, we calculate the integral of the problem:

$$\int_{x=y-\epsilon}^{x=y+\epsilon} \left( \sqrt{g(x)} G(x, y) - \frac{1}{2} l_c^2 \frac{d}{dx} \left( \frac{1}{\sqrt{g(x)}} \frac{d}{dx} G(x, y) \right) \right) dx \\ = \int_{x=y-\epsilon}^{x=y+\epsilon} \sqrt{g(x)} \delta(x, y) dx \qquad (C.10)$$

We then take the limit when $\epsilon \to 0$, and we obtain the following derivative jump:

$$\frac{dG(x, y)}{dx} \Big|_{x=y^+} - \frac{dG(x, y)}{dx} \Big|_{x=y^-} = \frac{2 g(x)}{l_c^2} \qquad (C.11)$$

## C.2. Equivalence between the Green's function of the eikonal equation and the weight function

The right-hand side of the original Eq. (C.6) is considered as a superposition of Dirac functions. Therefore, the solution to this problem is then obtained by the same superposition of the Green's functions $G(x, y)$, which can be rewritten as

$$\hat{\epsilon}^{\text{NL}}(x) = \int_{\Omega} G(x, y) \hat{\epsilon}(y) d\Omega \qquad (C.12)$$

This expression is of exactly the same form as Eq. (6) of the one-dimensional non-local integral model. Thus, the ENL-G damage model based on the differential Eq. (C.6) is a special case of the class of non-local model, in which the weight function $\Psi(x, y)$ is defined as $\Psi(x, y) = G(x, y)$.

<figure><img src="image_28.png" /><figcaption>(a) Numerical and symbolic comparison of the Green's functions (top). Damage field in the beam (bottom).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b) Numerically calculated Green's function.</figcaption></figure>

<caption>Fig. C.22. Green's function (ENL-G) for an internal length $c^{1/2} = 1$ m and an 1D undamaged beam.</caption>

## C.3. Numerical formulation and validation

The solution can be computed analytically using a symbolic solver for simple cases (no damage or damage with a polynomial field). Nevertheless, to solve the problem in a general case, a numerical solver needs to be used. In this study, we suggest using the Finite Element (FE) solver Freefem (Hecht, 2012).

### C.3.1. Numerical formulation

To compute the Green's function for any damage problem, we will solve the problem numerically using FE methods. Therefore, the problem associated with the BVP has been rewritten in weak form, so that it can be implemented in a FE solver. Then, the 1D problem (Eqs. (C.8) & ) was implemented in the FreeFem solver using the weak formulation to compute the Green's function.

### C.3.2. Numerical comparison

To validate the numerical implementation of the resolution of the Eikonal equation, the Green's function obtained numerically and symbolically (with Mathematica) are compared.
Fig. C.22 shows the Green's function obtained for a 1D undamaged beam, we observe that the Green's function is well approached compared to the symbolic solution, even on the border. It is interesting to note that we obtain the same results as Peerlings et al. (2001), since the ENL-G and NLG models are equivalent for undamaged material. In addition, on the basis of a few simple examples, different metrics and damage fields, where it is possible to calculate symbolically the Green's functions associated with the boundary condition problem, we were able to validate its numerical approximation by a FE calculation. This numerical method will be able to solve a wider range of problems, including different damage fields, different metrics and different boundary conditions.

## Data availability

The data set and software generated for this work is available on Védrine et al. (2025).

## References

Aifantis, E.C., 1987. The physics of plastic deformation. Int. J. Plast. 3 (3), 211-247. http://dx.doi.org/10.1016/0749-6419(87)90021-0.

Alam, S.Y., Kotronis, P., Loukili, A., 2013. Crack propagation and size effect in concrete using a non-local damage model. Eng. Fract. Mech. 109, 246-261. http://dx.doi.org/10.1016/j.engfracmech.2013.07.009.

Alam, S.Y., Saliba, J., Loukili, A., 2014. Fracture examination in concrete through combined digital image correlation and acoustic emission techniques. Constr. Build. Mater. 69, 232-242. http://dx.doi.org/10.1016/j.conbuildmat.2014.07.044.

Amiot, F., Hild, F., Roger, J.P., 2007. Identification of elastic property and loading fields from full-field displacement measurements. Int. J. Solids Struct. 44 (9), 2863-2887. http://dx.doi.org/10.1016/j.ijsolstr.2006.08.031.

Andrieux, S., Bamberger, Y., Marigo, J., 1986. A model of micro-cracked material for concretes and rocks. J. Mec. Theor. Appl. 5 (3), 471-513.

Avril, S., Grédiac, M., Pierron, F., 2004. Sensitivity of the virtual fields method to noisy data. Comput. Mech. 34 (6), 439-452. http://dx.doi.org/10.1007/s00466-004-0589-6.

Bažant, Z.P., 1986. Mechanics of distributed cracking. Appl. Mech. Rev. 39 (5), 675-705. http://dx.doi.org/10.1115/1.3143724.

Bažant, Z.P., Jirásek, M., 2002a. Nonlocal integral formulations of plasticity and damage: survey of progress. J. Eng. Mech. 128 (11), 1119-1149.

Bažant, Z.P., Jirásek, M., 2002b. Nonlocal integral formulations of plasticity and damage: Survey of progress. J. Eng. Mech. 128 (11), 1119-1149. http://dx.doi.org/10.1061/(ASCE)0733-9399(2002)128:11(1119), Publisher: American Society of Civil Engineers.

Bažant, Z.P., Nguyen, H.T., 2023. Proposal of m-index for rating fracture and damage models by their ability to represent a set of distinctive experiments. J. Eng. Mech. 149 (8), 04023047.

Bažant, Z.P., Pijaudier-Cabot, G., 1989. Measurement of characteristic length of nonlocal continuum. J. Eng. Mech. 115 (4), 755-767. http://dx.doi.org/10.1061/(ASCE)0733-9399(1989)115:4(755), Publisher: American Society of Civil Engineers.

Bažant, Z.P., Tabbara, M.R., Kazemi, M.T., Pijaudier-Cabot, G., 1990. Random particle model for fracture of aggregate or fiber composites. J. Eng. Mech. 116 (8), 1686-1705. http://dx.doi.org/10.1061/(ASCE)0733-9399(1990)116:8(1686).

Ben Azzouna, M., Périé, J.-N., Guimard, J.-M., Hild, F., Roux, S., 2011. On the identification and validation of an anisotropic damage model using full-field measurements. Int. J. Damage Mech. 20 (8), 1130-1150. http://dx.doi.org/10.1177/1056789510395555.

Benallal, A., Billardon, R., Geymonat, G., 1989. Conditions de bifurcation à l'intérieur et aux frontière pour une classe de matériaux non-standards. (bifurcation conditions inside and at the boundary for a class of non-standard materials). C. R. Acad. Sci. Sér. 2 Publisher: Gauthier-Villars.

Bonnet, M., 1999. Boundary Integral Equation Methods for Solids and Fluids, first ed. Wiley, Chichester ; New York.

Chan, N., Oliver-Leblond, C., Benboudjema, F., Ragueneau, F., Soleilhet, F., 2018. Prediction of the impact of shrinkage on concrete strength with the beam-particle approach. In: SynerCrete'18. Funchal, Portugal.

Chateau, C., Nguyen, T.T., Bornert, M., Yvonnet, J., 2018. DVC-based image subtraction to detect microcracking in lightweight concrete. Strain 54 (5), e12276. http://dx.doi.org/10.1111/str.12276.

Choi, S., Shah, S.P., 1997. Measurement of deformations on concrete subjected to compression using image correlation. Exp. Mech. 37 (3), 307-313. http://dx.doi.org/10.1007/BF02317423.

Claire, D., Hild, F., Roux, S., 2002. Identification of damage fields using kinematic measurements. C. R. Méc. 330 (11), 729-734. http://dx.doi.org/10.1016/S1631-0721(02)01524-3.

Claire, D., Hild, F., Roux, S., 2004. A finite element formulation to identify damage fields: the equilibrium gap method. Internat. J. Numer. Methods Engrg. 61 (2), 189-208. http://dx.doi.org/10.1002/nme.1057.

Claire, D., Roux, S., Hild, F., 2003. Identification of conductivités thermiques et de propriétés élastiques locales par analyse de champs. Méc. Ind. 4 (6), 655-665. http://dx.doi.org/10.1016/j.mecind.2003.09.007.

Delaplace, A., 2008. Modélisation discrete appliquée au comportement des matériaux et des structures. Mém. Habilit. Dir. Rech. Éc. Norm. Supér. Cachan.

Desmorat, R., Gatuingt, F., Jirásek, M., 2015. Nonlocal models with damage-dependent interactions motivated by internal time. Eng. Fract. Mech. 142, 255-275. http://dx.doi.org/10.1016/j.engfracmech.2015.06.015.

Doll, B., Ozer, H., Rivera-Perez, J., Al-Qadi, I.L., Lambros, J., 2017. Damage zone development in heterogeneous asphalt concrete. Eng. Fract. Mech. 182, 356-371. http://dx.doi.org/10.1016/j.engfracmech.2017.06.002.

Eliáš, J., Vořechovský, M., Skoček, J., Bažant, Z.P., 2015. Stochastic discrete meso-scale simulations of concrete fracture: Comparison to experimental data. Eng. Fract. Mech. 135, 1-16. http://dx.doi.org/10.1016/j.engfracmech.2015.01.004.

Feenstra, P.H., 1993. Computational aspects of biaxial stress in plain and reinforced concrete. Publisher: Delft University Press.

Geers, M.G.D., de Borst, R., Brekelmans, W.A.M., Peerlings, R.H.J., 1998. Strain-based transient-gradient damage model for failure analyses. Comput. Methods Appl. Mech. Engrg. 160 (1), 133-153. http://dx.doi.org/10.1016/S0045-7825(98)00011-2.

Geers, M.G.D., de Borst, R., Peijs, T., 1999. Mixed numerical-experimental identification of non-local characteristics of random-fibre-reinforced composites. Compos. Sci. Technol. 59 (10), 1569-1578. http://dx.doi.org/10.1016/S0266-3538(99)00017-2.

Geymonat, G., Hild, F., Pagano, S., 2002. Identification of elastic parameters by displacement field measurement. C. R. Méc. 330 (6), 403-408. http://dx.doi.org/10.1016/S1631-0721(02)01476-6.

Giry, C., Dufour, F., Mazars, J., 2011. Stress-based nonlocal damage model. Int. J. Solids Struct. 48 (25), 3431-3443. http://dx.doi.org/10.1016/j.ijsolstr.2011.08.012.

Giry, C., Oliver-Leblond, C., Ragueneau, F., Kishta, E., 2014. Non-local model and global-local cracking analysis for the study of size effect. In: Bí caní, N., Mang, H., Meschke, G., de Borst, R. (Eds.), Comput. Model. Concr. Struct. 255-262. http://dx.doi.org/10.1201/b16645-28.

Grassl, P., Bažant, Z.P., 2009. Random lattice-particle simulation of statistical size effect in quasi-brittle structures failing at crack initiation. J. Eng. Mech. 135 (2), 85-92. http://dx.doi.org/10.1061/(ASCE)0733-9399(2009)135:2(85).

Grassl, P., Grégoire, D., Rojas Solano, L., Pijaudier-Cabot, G., 2012. Meso-scale modelling of the size effect on the fracture process zone of concrete. Int. J. Solids Struct. 49 (13), 1818-1827. http://dx.doi.org/10.1016/j.ijsolstr.2012.03.023.

Grégoire, D., Rojas-Solano, L., Pijaudier-Cabot, G., 2013. Failure and size effect for notched and unnotched concrete beams. Int. J. Numer. Anal. Methods Geomech. 37 (10), 1434-1452. http://dx.doi.org/10.1002/nag.2180.

Grégoire, D., Verdon, L., Lefort, V., Grassl, P., Saliba, J., Regoin, J.-P., Loukili, A., Pijaudier-Cabot, G., 2015. Mesoscale analysis of failure in quasi-brittle materials: comparison between lattice model and acoustic emission data: Mesoscale analysis of failure in quasi-brittle materials. Int. J. Numer. Anal. Methods Geomech. 39 (15), 1639-1664. http://dx.doi.org/10.1002/nag.2363.

Guinea, G.V., Planas, J., Elices, M., 1992. Measurement of the fracture energy using three-point bend tests: Part 1-Influence of experimental procedures. Mater. Struct. 25 (4), 212-218. http://dx.doi.org/10.1007/BF02473065.

Haidar, K., Pijaudier-Cabot, G., Dubé, J.F., Loukili, A., 2005. Correlation between the internal length, the fracture process zone and size effect in model materials. Mater. Struct. 38 (2), 201-210. http://dx.doi.org/10.1007/BF02479345.

Havlásek, P., Grassl, P., Jirásek, M., 2016. Analysis of size effect on strength of quasi-brittle materials using integral-type nonlocal models. Eng. Fract. Mech. 157, 72-85. http://dx.doi.org/10.1016/j.engfracmech.2016.02.029.

Hecht, F., 2012. New development in FreeFem++. J. Numer. Math. 20 (3-4), 251-265.

Hild, F., Bouterf, A., Roux, S., 2015. Damage measurements via DIC. Int. J. Fract. 191 (1), 77-105. http://dx.doi.org/10.1007/s10704-015-0007-7.

Hillerborg, A., Modéer, M., Petersson, P.E., 1976. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem. Concr. Res. 6 (6), 773-781. http://dx.doi.org/10.1016/0008-8846(76)90007-7.

Hoover, C.G., P. Bažant, Z., Vorel, J., Wendner, R., Hubler, M.H., 2013. Comprehensive concrete fracture tests: Description and results. Eng. Fract. Mech. 114, 92-103. http://dx.doi.org/10.1016/j.engfracmech.2013.08.007.

Iacono, C., Sluys, L., van Mier, J., 2006. Estimation of model parameters in nonlocal damage theories by inverse analysis techniques. Comput. Methods Appl. Mech. Engrg. 195 (52), 7211-7222. http://dx.doi.org/10.1016/j.cma.2004.12.033.

Jirásek, M., 1998. Nonlocal models for damage and fracture: Comparison of approaches. Int. J. Solids Struct. 35 (31), 4133-4145. http://dx.doi.org/10.1016/S0020-7683(97)00306-5.

Jirásek, M., 2004. Non-local damage mechanics with application to concrete. Rev. Fr. Génie Civ. 8 (5-6), 683-707. http://dx.doi.org/10.1080/12795119.2004.9692625, Publisher: Taylor & Francis eprint: http://dx.doi.org/10.1080/12795119.2004.9692625.

Kachanov, L.M., 1958. On creep rupture time 8. pp. 26-31.

Kosin, V., Fau, A., Jailin, C., Hild, F., Wick, T., 2024. Parameter identification of a phase-field fracture model using integrated digital image correlation. Comput. Methods Appl. Mech. Engrg. 420, 116689. http://dx.doi.org/10.1016/j.cma.2023.116689.
Landis, E.N., 1999. Micro-macro fracture relationships and acoustic emissions in concrete. Constr. Build. Mater. 13 (1), 65-72. http://dx.doi.org/10.1016/S0950-0618(99)00009-4.

Le Bellégo, C., François Dubé, J., Pijaudier-Cabot, G., Gérard, B., 2003. Calibration of nonlocal damage model from size effect tests. Eur. J. Mech. A Solids 22 (1), 33-46. http://dx.doi.org/10.1016/S0997-7538(02)01255-X.

Lemaitre, J., 1971. Evaluation of Dissipation and Damage in Metals Submitted to Dynamic Loading. Evaluation of the Dissipation and de Lendommagement Dans les Metaux Soumis à Des Sollicitations Dynamiques. Technical Report N7211801, Office National Detudes Et de Recherches Aérospatiales, Paris (France), Num Pages: 20.

Lemaitre, J., Marquis, D., 1992. Modeling complex behavior of metals by the "state-kinetic coupling theory". J. Eng. Mater. Technol. 114 (3), 250-254. http://dx.doi.org/10.1115/1.2904169.

Li, S., Fan, X., Chen, X., Liu, S., Guo, Y., 2020. Development of fracture process zone in full-graded dam concrete under three-point bending by DIC and acoustic emission. Eng. Fract. Mech. 230, 106972.

Loiseau, F., Oliver-Leblond, C., Verbeke, T., Desmorat, R., 2023. Anisotropic damage state modeling based on harmonic decomposition and discrete simulation of fracture. Eng. Fract. Mech. 293, 109669. http://dx.doi.org/10.1016/j.engfracmech.2023.109669.

Mazars, J., 1984. Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure. These de docteur sœces présentée à l'Université Pierre et Marie Curie -Paris 6.

Mihashi, H., Nomura, N., Niiseki, S., 1991. Influence of aggregate size on fracture process zone of concrete detected with three dimensional acoustic emission technique. Cem. Concr. Res. 21 (5), 737-744. http://dx.doi.org/10.1016/0008-8846(91)90168-1-H.

Mindess, S., 1991. The fracture process zone in concrete. In: Shah, S.P. (Ed.), Toughening Mechanisms in Quasi-Brittle Materials. In: NATO ASI Series, Springer Netherlands, Dordrecht, pp. 271-286. http://dx.doi.org/10.1007/978-94-011-3388-3_17.

Moukarzel, C., Herrmann, H.J., 1992. A vectorizable random lattice. J. Stat. Phys. 68 (5), 911-923. http://dx.doi.org/10.1007/BF01048880.

Murali, K., Deb, A., 2018. Effect of meso-structure on strength and size effect in concrete under tension. Int. J. Numer. Anal. Methods Geomech. 42 (1), 181-207. http://dx.doi.org/10.1002/nag.2721.

Muralidhara, S., Prasad, B.K.R., Eskandari, H., Karihaloo, B.L., 2010. Fracture process zone size and true fracture energy of concrete using acoustic emission. Constr. Build. Mater. 24 (4), 479-486. http://dx.doi.org/10.1016/j.conbuildmat.2009.10.014.

Nguyen, M.T., Desceliers, C., Soize, C., Allain, J.M., Gharbi, H., 2015. Multiscale identification of the random elasticity field at mesoscale of a heterogeneous microstructure using multiscale experimental observations. Int. J. Multiscale Comput. Eng. 13 (4), http://dx.doi.org/10.1615/IntJMultCompEng.2015011435, Publisher: Begel House Inc.

Nikolić, M., Karavelić, E., Ibrahimbegovic, A., Miščević, P., 2018. Lattice element models and their peculiarities. Arch. Comput. Methods Eng. 25 (3), 753-784. http://dx.doi.org/10.1007/s11831-017-9210-y.

Oliver-Leblond, C., 2019. Discontinuous crack growth and toughening mechanisms in concrete: A numerical study based on the beam-particle approach. Eng. Fract. Mech. 207, 1-22. http://dx.doi.org/10.1016/j.engfracmech.2018.11.050.

Oliver-Leblond, C., Desmorat, R., Kolev, B., 2021. Continuous anisotropic damage as a twin modelling of discrete bi-dimensional fracture. Eur. J. Mech. A Solids 89, 104285. http://dx.doi.org/10.1016/j.euromechsol.2021.104285.

Otsuka, K., Date, H., 2000. Fracture process zone in concrete tension specimen. Eng. Fract. Mech. 65 (2), 111-131. http://dx.doi.org/10.1016/S0013-7944(99)00111-3.

Pathrage, M., Tong, D., Thierry, F., Cusatis, G., Grégoire, D., Pijaudier-Cabot, G., 2023. Discrete modeling of concrete failure and size-effect. Theor. Appl. Fract. Mech. 124, 103738. http://dx.doi.org/10.1016/j.tafmec.2022.103738.

Peerlings, R.H.J., De Borst, R., Brekelmans, W.A.M., De Vree, J.H.P., 1996. Gradient enhanced damage for quasi-brittle materials. Internat. J. Numer. Methods Engrg. 39 (19), 3391-3403. http://dx.doi.org/10.1002/(SICI)1097-0207(19961015)39:19<3391::AID-NME7>3.0.CO;2-D.

Peerlings, R.H.J., Geers, M.G.D., de Borst, R., Brekelmans, W.A.M., 2001. A critical comparison of nonlocal and gradient-enhanced softening continua. Int. J. Solids Struct. 38 (44), 7723-7746. http://dx.doi.org/10.1016/S0020-7683(01)00087-7.

Pijaudier-Cabot, G., Bažant, Z.P., 1987. Nonlocal damage theory. J. Eng. Mech. 113 (10), 1512-1533. http://dx.doi.org/10.1061/(ASCE)0733-9399(1987)113:10(1512), Publisher: American Society of Civil Engineers.

Pijaudier-Cabot, G., Grégoire, D., 2014. A review of non local continuum damage: Modelling of failure? Netw. Heterog. Media 9 (4), 575-597. http://dx.doi.org/10.3934/nhm.2014.9.575, Publisher: Networks and Heterogeneous Media.

Pijaudier-Cabot, G., Haidar, K., Dubé, J.-F., 2004. Non-local damage model with evolving internal length. Int. J. Numer. Anal. Methods Geomech. 28 (7-8), 633-652, https://doi.org- sid2nomade-2-grenet.fr/10.1002/nag.367.

Poh, L.H., Sun, G., 2017. Localizing gradient damage model with decreasing interactions. Internat. J. Numer. Methods Engrg. 110 (6), 503-522. http://dx.doi.org/10.1002/nme.5364, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5364.

Rabotnov, Y.N., Leckie, F.A., Prager, W., 1970. Creep problems in structural members. J. Appl. Mech. 37 (1), 249. http://dx.doi.org/10.1115/1.3408479.

Ribeiro Nogueira, B., Rastiello, G., Giry, C., Gatingt, F., Callari, C., 2023. Differential geometry-based thermodynamics derivation of isotropic and anisotropic eikonal non-local gradient (ENLG) damage models using a micromorphic media framework. Eng. Fract. Mech. 109670. http://dx.doi.org/10.1016/j.engfracmech.2023.109670.

Rima, A., Ragueneau, F., Oliver-Leblond, C., Guillet, T., Bonhomme, M.H., Nahas, G., 2021. Study of the behavior of post installed fasteners in civil engineering structures of french nuclear installations.

Rossi, P., Robert, J.L., Gervais, J.P., Bruhat, D., 1990. The use of acoustic emission in fracture mechanics applied to concrete. Eng. Fract. Mech. 35 (4), 751-763. http://dx.doi.org/10.1016/0013-7944(90)90158-D.

Roux, E., Bouchard, P.O., 2015. On the interest of using full field measurements in ductile damage model calibration. Int. J. Solids Struct. 72, 50-62. http://dx.doi.org/10.1016/j.ijsolstr.2015.07.011.

Skaržynski, L., Syroka, E., Tejchman, J., 2011. Measurements and calculations of the width of the fracture process zones on the surface of notched concrete beams. Strain 47 (s1), e319-e332. http://dx.doi.org/10.1111/j.1475-1305.2008.00605.x, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1475-1305.2008.00605.x.

Spelmen, V.S., Porkodi, R., 2018. A review on handling imbalanced data. In: 2018 International Conference on Current Trends Towards Converging Technologies. ICCTCT, IEEE, pp. 1-11.

Suchorzewski, J., Tejchman, J., Nitka, M., Bobiński, J., 2019. Meso-scale analyses of size effect in brittle materials using DEM. Granul. Matter 21 (1), 9. http://dx.doi.org/10.1007/s10035-018-0862-6.

Sutton, M.A., Orteu, J.J., Schreier, H., 2009. Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts,Theory and Applications. Springer Science & Business Media, Google-Books-ID: AlkqMxpQMLsC.

Tang, J.H., Chen, X.D., Dai, F., 2020. Experimental study on the crack propagation and acoustic emission characteristics of notched rock beams under post-peak cyclic loading. Eng. Fract. Mech. 226, 106890. http://dx.doi.org/10.1016/j.engfracmech.2020.106890.

Thierry, F., 2021. Modélisation de la localisation des déformations dans les milieux adoucissants par une approche eikonal (Ph.D. thesis). université Paris-Saclay.

Thirumalaiselvi, A., Sasmal, S., 2021. Pattern recognition enabled acoustic emission signatures for crack characterization during damage progression in large concrete structures. Appl. Acoust. 175, 107797. http://dx.doi.org/10.1016/j.apacoust.2020.107797.

Tibshirani, R., 1996. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B Stat. Methodol. 58 (1), 267-288.

Vandoren, B., Simone, A., 2018. Modeling and simulation of quasi-brittle failure with continuous anisotropic stress-based gradient-enhanced damage models. Comput. Methods Appl. Mech. Engrg. 332, 644-685. http://dx.doi.org/10.1016/j.cma.2017.12.027.

Vassaux, M., 2015. Comportement mécanique des matériaux quasi-fragiles sous sollicitations cycliques : de l'expérimentation numérique au calcul de structures. (Ph.D. thesis). École normale supérieure de Cachan - ENS Cachan.

Vassaux, M., Oliver-Leblond, C., Richard, B., Ragueneau, F., 2016. Beam-particle approach to model cracking and energy dissipation in concrete: Identification strategy and validation. Cem. Concr. Compos. 70, 1-14. http://dx.doi.org/10.1016/j.cemconcomp.2016.03.011.

Vassaux, M., Richard, B., Ragueneau, F., Millard, A., Delaplace, A., 2015. Lattice models applied to cyclic behavior description of quasi-brittle materials: advantages of implicit integration. Int. J. Numer. Anal. Methods Geomech. 39 (7), 775-798. http://dx.doi.org/10.1002/nag.2343, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.2343.

Védrine, L., Loiseau, F., Oliver-Leblond, C., Desmorat, R., 2025. Calibration of non-local damage models from full-field measurements: application to discrete element fields. http://dx.doi.org/10.5281/zenodo.13835397.

Wang, Y., Wagner, N., Rondinelli, J.M., 2019. Symbolic regression in materials science. MRS Commun. 9 (3), 793-805. http://dx.doi.org/10.1557/mrc.2019.85, Publisher: Cambridge University Press.

Xenos, D., Grégoire, D., Morel, S., Grassl, P., 2015. Calibration of nonlocal models for tensile fracture in quasi-brittle heterogeneous materials. J. Mech. Phys. Solids 82, 48-60. http://dx.doi.org/10.1016/j.jmps.2015.05.019.

Zhang, H., Šavija, B., Xu, Y., Schlangen, E., 2018. Size effect on splitting strength of hardened cement paste: Experimental and numerical study. Cem. Concr. Compos. 94, 264-276. http://dx.doi.org/10.1016/j.cemconcomp.2018.09.018.

Zhang, Y., Shedbale, A.S., Gan, Y., Moon, J., Poh, L.H., 2021. Size effect analysis of quasi-brittle fracture with localizing gradient damage model. Int. J. Damage Mech. 30 (7), 1012-1035. http://dx.doi.org/10.1177/1056789520983872, Publisher: SAGE Publications Ltd STM.

Zhang, J.Z., Zhou, X.P., 2022. Fracture process zone (FPZ) in quasi-brittle materials: Review and new insights from flawed granite subjected to uniaxial stress. Eng. Fract. Mech. 274, 108795.