A probabilistic peridynamic framework with an application to the study of the statistical size effect

Mathematical models are essential for understanding and making predictions about systems arising in nature and engineering. Yet, mathematical models are a simplification of true phenomena, thus making predictions subject to uncertainty. Hence, the ability to quantify uncertainties is essential to any modelling framework, enabling the user to assess the importance of certain parameters on quantities of interest and have control over the quality of the model output by providing a rigorous understanding of uncertainty. Peridynamic models are a particular class of mathematical models that have proven to be remarkably accurate and robust for a large class of material failure problems. However, the high computational expense of peridynamic models remains a major limitation, hindering outer-loop applications that require a large number of simulations, for example, uncertainty quantification. This contribution provides a framework to make such computations feasible. By employing a Multilevel Monte Carlo (MLMC) framework, where the majority of simulations are performed using a coarse mesh, and performing relatively few simulations using a fine mesh, a significant reduction in computational cost can be realised, and statistics of structural failure can be estimated. The results show a speed-up factor of 16x over a standard Monte Carlo estimator, enabling the forward propagation of uncertain parameters in a computationally expensive peridynamic model. Furthermore, the multilevel method provides an estimate of both the discretisation error and sampling error, thus improving the confidence in numerical predictions. The performance of the approach is demonstrated through an examination of the statistical size effect in quasi-brittle materials.


Introduction
Design approaches remain broadly the same across different engineering disciplines (e.g. aerospace, structural and mechanical). All disciplines depend heavily upon empirical formulas and large safety factors. This approach leads to highly conservative designs with low material utilisation. The benefits of improving material utilisation are clear (e.g. lighter vehicles achieve greater fuel efficiency) but certifying the safety and reliability of novel structural forms requires expensive programmes of testing. As the demand for more efficient structures increases, the need for new design approaches becomes more pressing. When experimental data is incomplete, a better approach to examine structural reliability might be provided through numerical simulations and stochastic methods.
Uncertainties in structural analysis arise from multiple sources, for example, material properties, loading conditions and geometry [1]. Current design approaches generally rely on deterministic models, and large safety factors must be applied to account for the inherent uncertainties. Uncertainties can be examined by using stochastic simulation methods, where uncertain input parameters, such as material properties, are treated as random variables.
Methods for the forward propagation of uncertainty (where sources of uncertainty are propagated through a model to evaluate the uncertainty in the output) can be broadly classified as intrusive or non-intrusive [2]. Intrusive uncertainty quantification (UQ) methods reformulate the original deterministic governing equations that describe the physical process. Non-intrusive UQ methods sample uncertain input parameters from a probability distribution and the deterministic governing equations are solved for each sample. The output is a distribution of the quantity of interest (QoI) from which various statistics, such as the mean and variance, can be computed. In this work, we employ the multilevel Monte Carlo (MLMC) method. The aim of MLMC is attain the same solution error as MC but at a significantly reduced computational cost. The standard MC estimator is computationally expensive as all samples must be computed using a fine mesh that guarantees a small discretisation error. A significant reduction in computational cost can be realised by taking the majority of samples on a coarse mesh (low accuracy but computationally cheap), and taking relatively few samples on a fine mesh (high accuracy but computationally expensive). This is made possible by isolating the error sources in the estimator: (1) sampling error (variance) and (2) discretisation error (deterministic error). The sampling error is controlled by using a low accuracy but computationally cheap model to take a large number of samples, and the discretisation error is reduced to a defined tolerance by employing a sufficiently fine mesh.
Multilevel techniques were first introduced by Heinrich and Sindambiwe [4] and Heinrich [5] and later popularised by Giles [6] for option pricing in computational finance. Cliffe et al. [7] were the first to apply multilevel methods in the field of engineering, motivated by the study of uncertainty in groundwater flows. Since Cliffe et al. [7] recognised the potential of multilevel methods, there has been a wide range of applications in engineering and scientific fields, for example, Dodwell et al. [8] employed MLMC to estimate the probability of failure of composite materials, and Clare et al. [9] assessed the risk of coastal flooding. For a detailed review of multilevel Monte Carlo methods, the reader is referred to the work of Giles [10].
The peridynamic theory of solid mechanics, introduced by Silling [11], is an integral-type non-local theory of solid mechanics that provides a robust theoretical framework for developing numerical models capable of simulating the failure behaviour of a wide range of materials. The peridynamic model defines material behaviour at a point in a continuum body as an integral equation of the surrounding displacement. This is in contrast to the classical theory of solid mechanics, where the material behaviour at a point is defined by partial differential equations. We focus on quasi-brittle materials because the range of experimental data is greater than that of any other material and quasi-brittle materials exhibit a significant size effect. A stochastic model is required for a complete examination of the mechanisms that govern the structural size effect. This work also provides new insights into the convergence behaviour of bond-based peridynamic models. To the best of the authors knowledge, convergence studies of the predicted structural response are missing from the peridynamics literature. Existing convergence studies only consider static elastic problems [12].
The paper is organised as follows: Section 2 introduces the peridynamic theory and Section 3 briefly describes the numerical model (a bond-based peridynamic model). Section 4 details the standard and multilevel Monte Carlo methodology. Section 5 describes the modelling of material properties as spatially correlated random fields. Section 6 presents two case studies.
The presented problems have been selected as examples where a deeper understanding of the physical behaviour can be gained by considering uncertainty. Section 7 discusses the results and Section 8 concludes the paper.

Peridynamic theory
The peridynamic theory, introduced by Silling [11] in 2000, is a non-local theory of solid mechanics that is formulated in terms of integral equations rather than partial differential equations. The governing equations do not require a spatially continuous and differentiable displacement field and damage localisation and fracture naturally emerge. No additional assumptions or techniques are required for modelling damage and fracture.
There are two primary formulations of the peridynamic theory: bond-based [11] and state based theory [13]. In the original bond-based theory, the Poisson's ratio is limited to a fixed value. Silling et al. [13] later introduced a generalised state-based theory that overcomes the limitations of the original theory. This paper employs the bond-based theory due to its lower computational expense and proven predictive capabilities.

Peridynamic continuum model
The bond-based peridynamic theory is briefly presented. It is not the purpose of this work to explain the peridynamic theory in detail and the reader is referred to [11,14,15] for a comprehensive treatment of the theory. A mechanically intuitive but less rigorous way of obtaining the governing equations can be found in [16]. Assuming that a body occupies a spatial region R, for any material point x ∈ R, a pairwise force function f can be defined to describe the interaction between particles within a finite distance δ of x, at any time t, where u represents the displacement of a material point (see Fig.   1) .
The peridynamic equation of motion for a single material point x at a point in time t is given by Newton's second law of motion, and is defined by Eq. (2).
ρ is mass density,ü is particle acceleration, b is body force per unit volume, and H x is the neighbourhood of material point x. The size of the neighbourhood is defined by the horizon length δ. For a 3D problem, the material point neighbourhood will be a sphere, and for a 2D problem, the neighbourhood will be circular.
The pairwise force function f represents the force that particle x exerts on particle x and contains all the constitutive information of the material under analysis. This interaction is commonly referred to as the peridynamic bond force. Particles separated by a distance greater than the horizon length δ do not interact. The pairwise force function f is defined by Eq. (4).
In a bond-based model, the force vector f is parallel to the deformed bond and the scalar bond force f (vector magnitude) is proportional to the bond stretch s. The initial relative position vector of a pair of particles is denoted by ξ = x − x, and the relative displacement vector is denoted by η = u − u. The current relative position vector is given by ξ + η.
To make a distinction between the peridynamic theory and other non-local theories, note that most non-local theories average some measure of strain within a neighbourhood of a material particle. The peridynamic theory dispenses with the concept of strain, which by its definition, requires the evaluation of partial derivatives of displacement [14].

Non-locality
The peridynamic theory is a non-local theory in which material points interact with each other directly over finite distances. This is in contrast to the classical theory of solid mechanics, where it is assumed that all forces are contact forces that act across zero distance (local theory).
Physical justification of non-locality was provided by Bažant [17], and further discussion on the origins of non-locality (with a focus on the peridynamic theory), can be found in Chapter 1 of Bobaru et al. [16] and Hobbs [18].
At the macroscale, the peridynamic horizon δ is a numerical constant with no physical meaning. This differentiates the peridynamic model from many numerical approaches, and the use of an ambiguous characteristic length parameter is avoided. For a given value of δ, the parameters in a peridynamic model can be chosen to match a given set of physically measurable material properties. Therefore, an optimum value of δ must be chosen that provides high accuracy whilst balancing computational expense. Section 3.2 discusses the selection of an optimum value of δ.
The reader should note the distinction between the non-local length scale in the peridynamic model (horizon δ), and the non-local length scale in a spatially correlated random field (correlation length l c ). The correlation length l c is generally considered to be a material parameter reflecting the internal length scale of the microstructure. This will be discussed throughout the paper.

Numerical model
To illustrate the framework, this work employs a two-dimensional bond-based peridynamic model. The aim of this work is to demonstrate the multilevel framework and a detailed treatment of the numerical model is not provided. All the results presented in this paper were obtained using an explicit scheme (as outlined in Fig. 4.14 of Hobbs [18]). The reader is referred to Hobbs [18] for implementation details. The main distinction of the model used in this work is the existence of two length scales: (1) the peridynamic horizon δ and (2) the correlation length l c in the random field. The generation of spatially correlated random fields is discussed in Section 5.

Constitutive model
It is generally assumed that the force-stretch (f -s) relationship of a peridynamic bond should be consistent with the macroscopic material response, and a failure mechanism is introduced into the model by eliminating the interaction between particle pairs when the stretch of the connecting bond exceeds a critical value. The stress-strain response of quasi-brittle materials is characterised by strain softening behaviour in the post-peak stage, and hence we employ the non-linear softening law, illustrated in Fig. 10, first proposed by Hobbs [18]. The derivation of the model parameters for the two-dimensional plane stress and plane strain case are provided in Appendix A.

Numerical convergence
The accuracy and convergence behaviour of a peridynamic model is complicated by the presence of a length scale. To determine an optimum value of δ, an additional parameter m must be introduced. m is the ratio between the horizon radius and grid resolution (m = δ/∆x). Bobaru et al. [19] and Ha and Bobaru [20] define and discuss two fundamental types of convergence: (1) m-convergence: δ is fixed and m → ∞. This can also be stated as δ is fixed and ∆x → 0. (2) δ-convergence: m is fixed and δ → 0. This can also be stated as m is fixed and ∆x → 0. See Fig. 2 for a graphical representation of the types of convergence. A third type of convergence can be defined: δm-convergence. This is a combination of δ-and m-convergence. See Bobaru et al. [19] for details.
In this work, we consider δ-convergence, as it is generally agreed that m should be close to 3. Madenci and Oterkus [21] investigated the choice of m for macroscale problems and it was found that values of m = 1 and m = 3 achieved the highest accuracy when compared to the classical analytical solution for the displacement of a one-dimensional bar subjected to a defined initial strain. Values of m much larger than 3 lead to excessive wave dispersion and become extremely computationally expensive. When fracture behaviour is also considered, values of m less than 3 lead to grid dependence on crack propagation [22,21]. Hu et al. [23] and Seleson [24] examined the m-convergence behaviour for two-dimensional models and Hobbs [18] examined the m-convergence behaviour for three-dimensional models. Higher values of m improve the spatial integration accuracy but m ≈ 3 provides an acceptable approximation. A value of m = π (δ = π∆x) is generally recommended for macroscale problems and is found extensively throughout the literature. The m-ratio is set to π for all problems in this paper. and ω represents a vector of random variables that takes values in R M . ω represents sources of uncertainty in the problem, in this case, the material properties. Note that the quantity of interest (Q) could be a function, for instance, the load-deflection response of a structure. For the presented case studies, the quantity of interest is the failure load, and the objective is to compute the expected value of Q, denoted E[Q], with a quantified level of uncertainty. However, for many real world applications, the probability distribution of Q is of more interest. Methods for obtaining the probability distribution of Q will be discussed.

Standard Monte Carlo simulation
In a standard Monte Carlo (MC) simulation, a large number (N ) of independent random realisations (or samples) of the parameters are generated. For every sample, the solution is computed using a numerical solver (finite element model, particle model etc). The accuracy of the solution is directly proportional to the resolution of the discretisation, and it is assumed An advantage of quantifying the accuracy of the estimator in this way is that the mean square error can be expanded and two distinct sources of error can be isolated: (1) the bias error and (2) the sampling error.
The first term in Eq. (7) is the bias error (sometimes referred to as the discretisation or numerical error). This arises as we are actually interested in the expected value E[Q] of Q, the unobtainable random variable corresponding to the exact solution without any numerical error.
If we assume that the numerical model converges to the exact solution as the mesh is refined, as M → ∞, then we can state the following where α is the order of convergence and α > 0 1 . The value of α is problem dependent and depends on numerous factors, such as, the chosen numerical model, the material model and the smoothness of the random field. By making M sufficiently large, the discretisation error can be reduced to any tolerance value b .
The second term in Eq. (7) is the sampling error and represents the variance of the estimator and decays inversely with the number of samples N . To ensure that the sampling error is less than a defined tolerance s , it is reasonable to determine the number of samples N using Eq.
To reduce the total error to a defined tolerance, the number of degrees of freedom M and the number of samples N must both be increased. This can be prohibitively computationally expensive when the cost to compute each sample to the required level of accuracy is high. The cost C M to compute a single sample of Q M is dependent on the computational complexity of the solver. The computational cost will grow as follows for some γ ≥ 1. The rate at which the computational cost grows (γ) is dependent on a number of factors, such as, the dimension of the problem and the chosen solver (explicit/implicit).
Standard MC estimators are proven to be robust and accurate for many problems, however the slow convergence rate limits applications to problems where the QoI can be computed cheaply. For problems that require the solution of computationally expensive models it is not possible to achieve reasonable estimations in an acceptable time. Different strategies have been examined to accelerate MC estimators, and all are based on the idea of reducing the sampling error.

Multilevel Monte Carlo simulation (MLMC)
The Multilevel Monte Carlo method (MLMC) was introduced by Giles [6] in 2008, but the first work on multilevel methods was by Heinrich and Sindambiwe [4] for parametric integration.
Further details on the origins of MLMC are provided in Giles [10]. Multilevel methods have been widely applied to engineering and scientific problems (i.e. solving partial differential equations).
Examples include the computation of the failure probability of composite structures by Dodwell et al. [8], the study of the travel time of particles through random heterogeneous porous media by Crevillén-García and Power [25], and the study of flood risk by Clare et al. [9].
The standard MC estimator is too costly as the quantity of interest for every sample must be computed to the level of accuracy required to ensure that the discretisation error is less than a defined tolerance. The key idea of MLMC is to compute a sequence of estimates of the quantity of interest using a hierarchy of nested meshes, as illustrated in Fig. 3. A significant reduction in computational cost can be realised by taking the majority of samples on a coarse mesh (low accuracy but computationally cheap), and taking relatively few samples on a fine mesh (high accuracy but computationally expensive). corresponds to a level 0 ≤ ≤ L in the multilevel method with M 0 < · · · < M < · · · < M L degrees of freedom.
We restrict ourselves to the case of uniform mesh refinement where the node spacing (∆x) is halved every time.
Because of the linearity of the expectation operator, the expected value of Q on the finest can be expressed as a telescopic sum of the expectation of Q on the coarsest mesh plus a sum of correction terms that account for the difference between evaluations on consecutive mesh levels.
where Y is the discrepancy between the QoI at successive mesh resolutions and is defined as The multilevel estimator for E[Q] is given by Eq. (13).
The number of samples N on each level is determined such that the total computational cost of the estimator is minimised for a defined sampling error (see Eq. (23)). It is important to highlight that the same random sample ω (i) is used to compute the quantity Q Fig. 4 for clarification). Note that the resolution of the coarsest level ∆x 0 must be smaller than the correlation length l c in the random field. Cliffe et al. [7] states that the optimal choice for the resolution of the coarsest mesh is such that ∆x 0 is slightly smaller than l c .
As all the expectations E[Y ] are estimated independently, the variance of the multilevel The accuracy of the estimator can be quantified by considering the mean square error.
Much like the standard MC estimator, the mean square error is composed of two terms, the bias error and the sampling error. The bias error is exactly the same as in the MC estimator (see Eq. (7)), and the number of degrees of freedom on the finest level (M L ) must be sufficiently large to satisfy Eq. (8), and thus ensuring that the bias error is less than b .
The multilevel estimator is cheaper than the standard MC estimator as the number of samples N on every level can be chosen to ensure that the sampling error is less than s , whilst minimising the total computational cost of the estimator. The computational cost of the multilevel Monte Carlo estimator is given by the following where C is the cost to compute a single sample of Y on level ≥ 1 or Q M 0 on level 0. Note that taking a sample of Y requires the numerical approximation of Q on two consecutive mesh levels M −1 must be computed). The determination of the optimal sample allocation is detailed in Section 4.2.2.
To achieve a RMSE of , it can be asserted that the multilevel estimator is computationally cheaper than the standard MC estimator due to the significant reduction in variance [7]. As the MLMC estimator is unbiased, the variance of the estimator is equal to The variance of the multilevel estimator is reduced as both numerical approximations Q M and Q M −1 converge to Q and consequently It is assumed that there exists a β > 0, where β is the order of convergence of the sampling error, such that By the central limit theorem, it is clear that fewer samples will be required to accurately approximate the expectation of the difference Q − Q −1 as → ∞. Consequently, the majority of samples will be taken on level 0 (computationally cheap), and relatively few samples will be required at the finest level L (computationally expensive).

Error estimation
The aim is to estimate E[Q] such that the RMSE is below a defined tolerance , whilst minimising the total computational cost of the estimator C(Q M L M ). The RMSE, defined by Eq. (14), is comprised of two parts: (1) the bias error and (2) the sampling error. To ensure that the RMSE is less than , it is sufficient to bound each term by 2 /2. To estimate the bias error, it is assumed that M is sufficiently large so that the decay in E[Q M − Q] is in the asymptotic region and satisfies the following Following the derivation of Dodwell et al. [8], for uniform mesh refinement, where the number of degrees of freedom on level is given by M ≈ m M 0 , the bias error on level can be over-estimated as follows where r is set to 1. This is equivalent to the assumption that M is sufficiently large so that the decay in E[Q M − Q] is in the asymptotic region. The user may wish to select a more conservative values for r, for example 0.7 or 0.9. If the bias error is greater than the tolerance, then M L must be increased.
To ensure that the sampling error is less than or equal to the sample tolerance s , the following constraint is enforced As the number of samples increases, the variance of the sample mean decreases and hence precision increases. The sample variance is estimated in the standard way [8]

Sample allocation
The optimal sample allocation (number of samples per level N ) is determined by solving a constrained optimisation problem that minimises C(Q M L M ) with respect to N , subject to the constraint that the sampling error of the multilevel estimator is less than or equal to the defined tolerance s .
The computational cost of the MLMC estimator grows as follows: The rate at which the computational cost grows with respect to the number of degrees of freedom M is given by Eq. (25), for some γ ≥ 1.
The reader is referred to Cliffe et al. [7] and Giles [10] for a full proof of the MLMC computational complexity theorem with bounds on the RMSE. converged == True  [26]. For a detailed treatment of the theory of random fields and further applications, the reader is referred to Hristopulos [27].

MLMC implementation
In this contribution, matrix decomposition and KL expansion are employed because of their practical simplicity. The matrix decomposition method generates accurate spatial random fields, but the computational expense is prohibitive for large-scale problems. KL expansion produces less accurate random fields but due to the lower computational cost and ease of implementation, this method was employed for all considered problems. Examining large problems, such as three-dimensional models, is prohibitively expensive. To overcome this issue, the spatial domain can be split into several smaller sub-domains, and a sample of the random field is generated for each sub-domain. A sample for the entire domain is obtained by using an overlapping technique [28,29].

Covariance function and length scale
The value of a random variable at two adjacent points in space is correlated. Conversely, there is negligible correlation between the value of a random variable at two distant points.
Many choices exist for the covariance function. Popular choices include Matérn covariance, exponential covariance, Gaussian covariance, spherical covariance and many others [26]. In this work we have employed an exponential covariance function, as defined by Eq. (26), where ρ is the correlation coefficient between a random variable at point x i and x j , σ 2 is the variance (set to 1), l c is the correlation length and x i − x j 2 is the Euclidean distance between two material points. This form has been selected due to its popularity in the literature.
There is limited guidance in the literature for selecting a suitable covariance function for different material types. The Joint Commission of Structural Safety (JCSS) Probabilistic Model Code [30] provides guidance on the probabilistic assessment of concrete structures. The correlation coefficient ρ ij between a random variable at point x i and x j is defined by Eq. (27).
The covariance function defined in the JCSS probabilistic model code is unusual as the function contains a threshold value for ρ. According to the JCSS probabilistic model code the default threshold value is 0.5. To the best of our knowledge, this approach is not seen elsewhere in the literature. By setting the threshold value to 0, the exponential covariance function is obtained.
The distance over which a correlation exists is determined by the length scale. The correlation length l c is a highly uncertain parameter that has a significant influence on the final results. For quasi-brittle materials, Grassl and Bažant [31] suggested that the correlation length must, at a minimum, be as large as the fracture process zone (FPZ). For concrete, the size of the FPZ is approximately two to three times the maximum aggregate size [32]. The JCSS probabilistic model code recommends a correlation length of 5 m. This value is significantly higher than that found elsewhere in the literature and there is no clear rationale. Fig. 5 illustrates the influence of the correlation length l c on the generated random field. It is not the purpose of this paper to examine the correlation length in detail but our studies suggest that a shorter correlation length improves convergence. Note that l c must be greater than the resolution of the coarsest level ∆x 0 or the value of a random variable at two adjacent points in space will be uncorrelated (white noise). For a detailed examination of the influence of the correlation length l c , the reader is referred to the work of Syroka-Korol et al. [33]. that the correlation length l c must be greater than the resolution of the coarsest level ∆x 0 or the value of a random variable at two adjacent points in space will be uncorrelated (white noise).

Material strength distribution (probability distribution function)
The chosen material strength distribution plays a key role in the predicted results and convergence of the model. In the literature, normal (Gaussian), log-normal, Gauss-Weibull and Weibull distributions have all been employed for modelling quasi-brittle materials.
The choice of a normal distribution is generally made for convenience as opposed to physical reasons [34]. In particular, material parameters are usually bounded (values must be positive), but negative values are possible when using a normal distribution. Our preliminary studies determined that a normal distribution is not suitable for modelling quasi-brittle materials. This is discussed further in Section 6.3. The JCSS Probabilistic Model Code recommends that the properties of quasi-brittle materials are modelled using a log-normal distribution [30]. Van der Have [35] provides a detailed study of random field generation and the differences between the use of normal and log-normal distributions are explored. The log-normal distribution guarantees that the material parameters are positive, is easy to implement and is widely employed throughout the literature. However, it has been demonstrated that on the scale of a representative volume element (RVE), the probability distribution of strength of quasi-brittle materials is best approximated by a Gaussian distribution onto which a far-left Weibull tail is grafted [36,37,38]. Eliáš et al. [39] and Eliáš and Vorechovský [40] modelled the size effect in quasi-brittle materials using a lattice discrete particle model (LDPM) and the cumulative distribution function of the random field was assumed to be Gaussian with a left Weibullian tail. The far-left tail of the strength distribution has a huge influence on the the failure load when considering small failure probabilities. For example, for a failure probability of 10 −6 (structures are generally designed for a failure probability lower than 10 −6 per lifetime [41]), the difference between the failure load and the mean strength will almost double when the strength distribution changes from Gaussian to Weibull (with the same mean and coefficient of variation) [36]. It should be noted that the modelling of quasi-brittle materials is complicated by the transition of the strength distribution from Gaussian to Weibullian as the structure size increases [37].
It is not the purpose of this paper to examine different strength distributions in detail but we explored the use of a normal, log-normal and Weibull distribution. The Weibull distribution provided the best agreement with experimental data and improved the rate of convergence of the discretisation error. This is discussed further in Section 7.4.
To easily generate a random field, where the probability distribution function of a material parameter at a given location is a univariate Weibull distribution, we follow the approach of Rappel et al. [42] and Rappel et al. [43]. In a Gaussian random field, the probability density function of a material parameter at a given location is a univariate Gaussian distribution. Using The following subsection briefly discuss the structural size effect. For a detailed review of the structural size effect, the reader is referred to Bažant and Planas [45] and Bažant [46].

Structural size effect
Based on the strength-of-materials theory, structural failure is assumed to occur when the maximum stress in a structure exceeds some limiting value of stress that can be determined from small scale tests of representative material samples. Simple fundamental tests such as uniaxial tension, uniaxial compression and flexural tests are used to establish the limiting stress for different loading conditions. This simplistic view does not suffice for quasi-brittle materials [47]. Quasi-brittle materials exhibit a size effect where large elements fail at lower stresses than geometrically identical but smaller-scale elements.
In brittle and quasi-brittle materials, the size effect can primarily be explained by two mechanisms [45,46]: (1)  Hobbs et al. [49] previously examined the size effect in quasi-brittle materials using a deterministic bond-based peridynamic model. The model did not consider the spatial variability in material properties and the magnitude of the statistical size effect remains to be established.
Due to the high computational expense of peridynamic simulations, examining the statistical size effect was impracticable but the presented framework allows us to overcome the aforementioned issues. Hobbs et al. [49] validated the deterministic model against the full set of experimental results published by Grégoire et al. [50]. This work only considers two members from the test series as the aim of this study is to demonstrate the possible computational savings that can be realised using the MLMC framework, and to demonstrate the importance of examining uncertainty. Future work will use the MLMC framework to examine the full series of tests and provide a comprehensive study of the statistical size effect.
6.2. Case study 1: Statistical size effect in quasi-brittle materials (Type 2) The first problem that we consider is a notched concrete beam in three-point bending, tested experimentally by Grégoire et al. [50]. A schematic diagram of the experimental setup is illustrated in Fig. 6. The mean compressive strength f cm,cyl = 42.3 MPa is used to generate a realisation of the random field. The Young's modulus E, tensile strength f t and fracture energy G F are then estimated using empirical equations [44]. The density of the concrete mixture was 2346 kg/m 3 and the maximum aggregate diameter was 10 mm. The correlation length l c is set to 20 mm.
Please refer back to Section 5.1 for a discussion of the length scale. The Weibull modulus m is set to 3. This is an uncertain parameter with high sensitivity and a wide range of values can be found in the literature. According to the Weibull theory, the modulus m is a material property that is independent of the geometry and scale of the structure, however Syroka-Korol et al. [33] found that the Weibull modulus m does depend on the size of the structure and length scale l c . All the presented results have been obtained using a constant peridynamic horizon δ = 3.14∆x and regular grid spacing. Table 1

Results
We start by taking 100 samples on all levels and estimate α, β and γ. The first step is to estimate how the computational cost scales as M increases. The time to compute each sample is recorded and it is determined that the computational costs grows linearly. The computational cost is given by Eq. (28), where γ = 1. Note that the performance of PeriPy scales linearly with the number of nodes ∴ γ = 1 [51].
The next step is to estimate the parameters α and β for the QoI, which is taken to be the peak load. Fig. 7 illustrates the log-log plots of the estimated means and variances of Q and Y = Q − Q −1 , for = 0, ..., 4, with respect to the number of degrees of freedom M on each level. The rate of convergence of the discretisation error is given by Eq. (29), where α is approximately 0.528.
The rate of convergence of the sampling error is given by Eq.      [39] who found that considering spatial variability in material properties does not significantly influence the mean failure load, but does lead to an increase in the variance of the structural response. Using the finest mesh ( = 4), the deterministic model of Hobbs et al. [49] predicts that the specimen will fail at approximately 1800 N. Setting the sampling tolerance s to 10 N, the mean stochastic strength is predicted to be approximately 1790 N. Note that the bias error is approximately 0.75 N. The observed results are in agreement with theory, which predicts that the difference between the deterministic strength and mean stochastic strength will be small [52,48]. Note that the experimental failure load ranged between 1580 N and 1710 N.

Case study 2: Statistical size effect in quasi-brittle materials (Type 1)
The second problem that we consider is an unnotched concrete beam in three-point bending, tested experimentally by Grégoire et al. [50]. We consider Specimen 3 (illustrated in Fig. 6) again but with no notch (λ = 0). Beyond demonstrating the computational savings that can be achieved using the MLMC framework, the presented example provides insight into the following areas: Statistical size effect -Hobbs et al. [49] showed that a deterministic bond-based model accurately captures the structural size effect for Type 2 (notched) problems, but fails to capture the correct response for Type 1 (unnotched) problems. This was expected as it is well known that the randomness of material properties has a significant effect on the structural strength of Type 1 problems [53,40]. In Type 1 problems, the volume of highly stressed material is much larger than that observed in Type 2 problems, and the probability that a defect is present in the stressed region is consequently higher. In Type 2 problems, the presence of a notch results in a localised region of highly stressed material, and the influence of randomness in material properties is consequently lessened. It is expected that the inclusion of statistical variability in the material properties will improve the predictive accuracy of the peridynamic model. [18] demonstrated that a deterministic bond-based model fails to converge for Type 1 problems (the predicted strength is coupled with the mesh resolution).

Convergence -Hobbs
It was hypothesised that accounting for randomness in the material properties is required to initiate the localisation of damage and improve convergence.

Results
Again we start by taking 100 samples on all levels and estimate α, β and γ. As per the previous example, the computational cost grows linearly (γ = 1). Taking 100 samples on every level, α is estimated to be 0.337 and β is estimated to be 0.682 (refer to Fig. 9). The rate of convergence of the discretisation error and sampling error is slower than that observed in problem 1 (Type 2). Using the estimated values of α, β and γ, Eq. (24) predicts that the cost of the MLMC simulations will grow proportionally to −2.94 , whilst the cost of the standard MC simulations will grow proportionally to −4.97 . Table 4 presents the optimal number of samples N across the mesh levels for different values of sampling tolerance ( s = 10, 50 and 100 N), plus the number of samples required when using the standard MC estimator (N ). Due to the higher variance of the estimator, the number of samples required is considerably higher than that required for the Type 2 problem. Type 1 problems are subject to a high degree of natural variability and consequently the computational cost is higher as significantly more samples are required.  By including uncertainty in the material properties, the bond-based model converges for Type 1 problems (α ≈ 0.337). This is the first time that this behaviour has been demonstrated, but the convergence behaviour is significantly worse than that observed for Type 2 problems Weibull distribution with a low Weibull modulus) but this might not be physically realistic for the considered problem.
As the size of a structure increases, so does the probability that a defect will be present from which a fracture will initiate. Syroka-Korol et al. [53] determined numerically that the deterministic and mean stochastic strength start to diverge when the beam depth is greater than 50-60 mm. Specimen 3 is 100 mm deep and the magnitude of the statistical size effect is expected to be non-negligible. Setting the sampling tolerance s to 50 N, the mean stochastic strength is estimated to be approximately 6250 N. Note that the bias error is approximately 200 N. Using the finest mesh ( = 4), the deterministic model predicts that the specimen will fail at approximately 9200 N. The experimental failure load ranged between 7620 N and 8770 N. The numerical results are consistent with the theory, i.e., the difference between the deterministic strength and mean stochastic strength is much larger than that observed for Type 2 problems.
However, the deterministic model does not converge for Type 1 problems and the prediction of strength is therefore unreliable, and a rigorous comparison is not possible.
The objective of the multilevel framework is to estimate the expectation of an output variable, in this case, the peak load. However, for many industrial applications, engineers are more concerned with the probability of an output variable exceeding a specific value and the cumulative distribution function (CDF) is needed. It is possible to obtain the CDF by following the method outlined in Gregory and Cotter [54]. The reader is also referred to Clare et al. [9] for further information.

Statistical size effect
A key aim of this study was to select case studies where uncertainty must be considered to gain a comprehensive understanding of the physical behaviour. We focussed our studies on the structural size effect in quasi-brittle materials. Bažant [17] stated that the correct modelling of the size effect on material strength should be adopted as the basic criterion of acceptability of any model. The results demonstrate that a bond-based peridynamic model can be used to examine both the statistical and deterministic component of the structural size effect. The intention of this study was never to provide a detailed examination of the statistical size effect, and further studies on a wider range of problems are required to improve confidence in the models predictive capabilities.
By employing the presented MLMC framework, studying the statistical component of the structural size effect using a peridynamic model becomes computationally feasible. Future work aims to employ the presented MLMC framework to study the full series of tests published by Grégoire et al. [50] and provide a detailed examination of influential factors, such as the shape of the material strength distribution and the correlation length l c . Grassl and Bažant [31] state that the ratio of the correlation length l c to the size of the fracture process zone (FPZ) is the main parameter that influences the statistical size effect.

Convergence
Numerical results should be independent of the mesh resolution. This is a basic test of the adequacy of any numerical model. To the best of the authors knowledge, Hobbs [18] was the first to consider the effect of mesh refinement (δ-convergence) on the predicted peak load and load-deflection response for Type 1 and 2 problems. Hobbs found that a deterministic bond-based peridynamic model fails to converge for Type 1 problems. Note that Niazi et al. [55] also published a convergence study that considered the complete structural response. The study of Niazi et al. [55] is limited as Type 1 problems were not considered.
The results in this study confirm that, as previously hypothesised, a source of randomness must be introduced to trigger the localisation of damage in Type 1 specimens and eliminate problems of mesh dependence that occur in peridynamic models. Niazi et al. [55] reported that the convergence behaviour is improved by randomly deleting 1% of all bonds, as first suggested by Chen et al. [56]. Whilst the method of Chen et al. [56] is computationally cheap and does improve convergence behaviour, it is an oversimplified approach that lacks a robust theoretical basis (heuristic) and does not consider the spatial correlation of material properties. Jones et al. [57] note that these methods are generally used to avoid problems related to symmetry, and they do not attempt to capture the true material behaviour by implementing an experimentally measured probability distribution of material properties.

Length scales
The correlation length l c was set to be 20 mm for all considered problems. This value was selected after running a number of preliminary simulations. However, the aim of this contribution was not to identify the parameters that describe the spatial fields. It is important to note that a theoretically grounded probabilistic framework based on Bayesian inference (see [42,43,58]) is essential to identify the parameters of the spatial fields (e.g. length scale l c ) rigorously.
Furthermore, the interaction between the two length scales (peridynamic horizon δ and the correlation length l c in the random field) requires further examination. It remains uncertain how the ratio of the two length scales influences the predictive accuracy of the model.

Probability distribution
The material strength distribution plays an important role in the predicted results and convergence of the model. Three distribution were considered (normal, log-normal and Weibull) and it was determined that the Weibull distribution provides the best predictions of mean strength for quasi-brittle materials. This was expected and has been extensively discussed in the literature. A more novel observation is that the selected probability distribution influences the convergence rate of the bias error. Extreme values in the left-tail are required to initiate the localisation of damage and eliminate problems of mesh dependence. Note that the model failed to converge for Type 1 problems when using a normal distribution.

Model calibration
Many of the model parameters are impossible to determine exactly and are subject to significant uncertainties, for example: the length scale l c and the Weibull modulus (shape parameter). Future work will examine the integration of the multilevel method with experimental data in a Bayesian setting to quantify modelling uncertainties as proposed by Dodwell et al. [59,60]. This will be an important step in the validation of peridynamic models, enabling the identification of model discrepancy and measurement bias, and providing better estimates of model parameters.

Conclusions
Peridynamic models are computationally expensive, thus preventing the use of standard Monte Carlo methods for the assessment of uncertainties in model outputs propagated from uncertain inputs. The aim of this study was to demonstrate the possible computational savings that can be realised using the MLMC framework. The results show a speed-up factor of 16× over a standard Monte Carlo estimator, enabling the forward propagation of uncertain parameters in a computationally expensive peridynamic model. Beyond demonstrating the computational savings that can be achieved using the multilevel framework, the results presented in this paper are of interest for two further reasons: 1. Deterministic bond-based models suffer from a strong mesh dependency when simulating Type 1 problems. It has been demonstrated that by including uncertainty in the material properties, the bond-based peridynamic model converges for both Type 1 and Type 2 problems. The need to consider uncertainty is essential for robust and accurate predictions.
Furthermore, the multilevel method provides an estimate of the discretisation error, thus improving the interpretability of numerical predictions.

2.
A secondary aim was to select case studies where uncertainty must be considered to gain a comprehensive understanding of the physical behaviour. We examined the structural size effect in quasi-brittle materials as the random variability of material properties is known to play an important role. Bažant [17] stated that the correct modelling of the size effect on material strength should be adopted as the basic criterion of acceptability of any model. The results demonstrate that a bond-based peridynamic model can be used to study the statistical size effect but further studies on a wider range of problems are required to improve confidence in the models predictive capabilities. Future work will consider the full series of tests published by Grégoire et al. [50] and provide a detailed study of the statistical size effect.
We have motivated the use of the MLMC framework by studying the statistical size effect in quasi-brittle materials. But forward uncertainty quantification is equally important for cases where a high degree of reliability is required, as is common in many aerospace and power generation applications.
derived for quasi-brittle materials from the experimental work of Cornelissen et al. [61]. Note that the area under the σ-w curve is a measure of the material fracture energy G F .
The σ-w relationship is described by an exponentially decaying model with a term that forces the curve to intersect with the horizontal axis at w c . If the softening relationship is asymptotic with the horizontal axis, and thus never intersects, a unique value for the critical stretch of a bond cannot be determined.
The bond stiffness constant for different problem types is defined by Eq. (32), where t is the thickness of the domain under analysis.
The evolution of the non-linear bond softening parameter d is defined by Eq. (33). This function describes an exponentially decaying curve with a linear term. As the bond stretch s approaches the critical stretch s c , the linear term forces the softening curve to decay linearly and intersect with s c . α controls the position of the transition from exponential to linear decay, and k controls the rate of exponential decay. The linear elastic limit s 0 is defined empirically as f t /E. This definition of s 0 is not objective but it has been shown to provide good results.
The energy required to break a bond is defined by Eq. (34). Only the energy consumed during the softening stage is considered (between the limits s 0 and s c ). It is important that the softening curve intersects with s c so that the integral in Eq. (34) can be evaluated.
The proposed model provides an explicit definition of the critical stretch s c and an unambiguous relationship between s c , k and α.