An examination of the size eﬀect in quasi-brittle materials using a bond-based peridynamic model

This paper examines the size eﬀect in quasi-brittle materials using a three-dimensional bond-based peridynamic model. This is the ﬁrst time that the capability of a peridynamic model to capture the size eﬀect in quasi-brittle materials has been examined. Correctly reproducing the size eﬀect is an essential check on the validity of any computational model and it is demonstrated that a bond-based peridynamic model can accurately capture the failure stress of geometrically identical structures over a range of sizes. A systematic examination of geometrically similar notched and unnotched beams provides new insights into the predictive capabilities of the model, and mode I and mixed-mode problems are considered to examine the generality of the model. The model is validated using published experimental data, and the predictive accuracy is equivalent, if not superior, to well-established numerical methods whilst oﬀering several beneﬁts that justify further research and development. This study provides evidence that the length scale in a peridynamic model is a numerical parameter and not a material property.


Introduction
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 on 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 [1]. Quasi-brittle materials exhibit a size effect where large elements fail at lower stresses than geometrically identical but smaller-scale elements. Accurately capturing the structural size effect is essential for safe predictions of load capacity.
According to the strength-of-materials theory, the maximum stress that a material can resist is independent of size when geometrically similar specimens are considered. Any deviation from predictions made using stress failure criteria is known as the structural size effect. The maximum stress that a quasi-brittle material can resist deviates significantly from predictions made using stress failure criteria, and large specimens fail at significantly lower stresses than small specimens.
There are two primary sources of size effect in quasi-brittle materials [2]: (1) release of stored energy (deterministic size effect), and (2) statistical variability in material properties (statistical size effect). The release of stored energy is by far the most important factor influencing the size effect on structural strength, and the statistical size effect is of secondary importance. The deterministic size effect is governed by the size of the fracture process zone (zone of energy dissipation) relative to the size of the structure. The statistical size effect is a result of the randomness of material properties and defects. The probability that a specimen contains a defect from which failure will initiate increases as the size of the specimen increases.
Understanding the structural size effect is essential for safe predictions of load capacity. The design of large scale structures relies on the properties of quasi-brittle materials measured in small scale laboratory tests. These properties must be extrapolated to structures that are one or two orders of magnitude larger than laboratory tests. Full scale tests are prohibitively expensive and reliable models that can capture structural size effect are needed. Accurately capturing the size effect is of utmost importance for the design of efficient structural forms with small margins of safety. Existing analytical formulas for predicting size effect in geometrically similar structures are extremely limited, and they generally employ empirical parameters that must be determined through best-fit procedures [3]. Numerical methods offer the only practicable approach for addressing complex industry motivated problems. This paper examines the size effect in quasi-brittle materials using a deterministic threedimensional bond-based peridynamic model. The capability of a peridynamic model to capture the size effect in quasi-brittle materials has not been investigated yet and remains a major question. Understanding the predictive capabilities of the peridynamic model is essential for further application to more complex industry-motivated problems. As stated by Bažant [4], the capability to correctly reproduce the size effect is an important check on the validity of any computational model.
There are relatively few papers in the peridynamic literature that compare numerical results against experimental data, and when a comparison is made, the majority of papers calibrate model parameters against a single experiment, and only a few provide further validation. Diehl et al. [5] provided a review of benchmark experiments for the validation of peridynamic models and noted the same issue. Hobbs [6] also discusses this issue at length and expands on the key elements of model validation. The lack of rigorous validation studies is a concern that is limiting the application of peridynamic models. This study quantifies the predictive capabilities of the peridynamic model against carefully selected benchmark tests. Geometrically similar beams of different sizes are modelled to examine the size effect on structural strength. Section 2 provides a review of structural size effect and the mechanisms behind size effect are examined. Section 3 introduces the peridynamic theory and discusses the advantages of using a bond-based peridynamic model. Section 4 briefly outlines the numerical framework. Mode I problems are examined in section 5 and the effect of boundary types is explored. Mixed-mode problems are examined in section 6. Section 7 provides a discussion of the results and areas of future work are identified. Section 8 concludes the paper.

What is size effect?
According to the classical failure theories, such as elasticity with a material strength limit and plasticity theory, the nominal strength σ n of a structure is independent of size when geometrically similar structures are considered. The classical theories express the material failure criterion in terms of stress and strain and the nominal strength σ n is defined as the maximum stress that a structure can resist [4]. For a concrete beam in bending, the limiting nominal strength σ n is given by the ultimate tensile strength f t . Any deviation from predictions made using stress failure criteria is known as the structural size effect. A size effect is observed in many materials.
To clarify further, the basic theory of strength of materials predicts that a small and large beam made of the same material will fail at the same value of stress. In reality, the larger beam will fail at a lower stress. The post-peak behaviour of quasi-brittle materials also exhibits a size effect. At a small scale, concrete behaves like a ductile material and the post-peak load deflection curve descends slowly. At a sufficiently large scale, concrete becomes perfectly brittle, and the load-deflection curve descends rapidly and may exhibit snap-back behaviour [7]. Clearly, a proper understanding of structural size effect is essential for safe predictions of load capacity.
For a detailed examination of size effect in quasi-brittle materials, see the work of Bažant and Planas [2] and Bažant [7,4].

Size effect mechanisms
In brittle and quasi-brittle materials, the size effect can primarily be explained by two mechanisms [2,4]: (1) release of stored energy, and (2) statistical variability in material properties. These effects are often referred to as the deterministic and statistical size effect.
Secondary factors that influence the size effect include the boundary layer effect, diffusion phenomena, and hydration heat [8]. The influence of these secondary factors is insignificant and will not be considered any further. The release of stored energy is by far the most important factor influencing the size effect on structural strength. The influence of statistical variability in material properties is less significant. But until the late 1980's, the primary cause of size effect was generally believed to be of statistical origin [9].
When a fracture initiates and begins to propagate, the elastic strain energy stored in the structure is consumed by the fracture propagation. In quasi-brittle materials, energy is dissipated in a non-linear zone of micro-cracking formed at the crack tip. This damage zone is known as the fracture process zone (FPZ). The ratio between the elastic energy stored in the structure and the energy consumed by the fracture process varies with the size of the structure. This is the main mechanism responsible for structural size effect. The mechanism by which the statistical size effect occurs is simpler to comprehend. As the size of a structure increases, so does the probability that a defect will be present from which a fracture will initiate.
Studying the mechanisms that govern size effect is limited by the physical size of the experiments that can be performed within a reasonable timescale and budget. Many structures of interest, such as dams and reactor containment buildings, are far beyond the scale of laboratory experiments. Numerical methods provide new ways to explore the structural size effect without the aforementioned limitations.

Size effect laws
Numerous analytical expressions have been developed to describe the size effect of quasibrittle materials. Size effect laws provide predictions of nominal strength as a function of the characteristic size of the structure for geometrically similar structures, with and without notches. These expressions have been developed from experimental observations and theoretical considerations.
The central characteristic of quasi-brittle materials is that there exists a sizeable fracture process zone and area of softening damage. The size effect is governed by the size of the fracture process zone relative to the size of the structure and depending on the structure size, different theories are applicable for predicting strength. At a sufficiently small scale, the theory of plasticity, which assumes no size effect, is applicable. The theory of plasticity is valid when the size of the fracture process zone is comparable to the structural size. At a sufficiently large scale, the theory of linear elastic fracture mechanics, which predicts the strongest possible deterministic size effect, is applicable. The theory of linear elastic fracture mechanics is valid when the size of the fracture process zone is negligible with respect to the structural size. Many attempts have been made to bridge the gap between the theory of plasticity and the theory of linear elastic fracture mechanics.
Numerous expressions have been proposed to combine the different theories and provide a general size effect law. The analytical expressions defined by Eq. 1 were proposed by Bažant [10] and Bažant and Planas [2] to describe the deterministic size effect. Two types of size effect law are defined: Type 1 -structures with no notches or pre-existing cracks (fracture initiates from a smooth surface); Type 2 -structures with a notch. Notched problems are representative of reinforced concrete structures, where due to the reinforcement, the stable growth of cracks occurs prior to the ultimate load.
Type 2 (notched case) (1) σ N is the nominal strength and D is the characteristic size of the structure. For Type 1 structures, f ∞ r , D b , l p , and r are positive constants representing unknown empirical parameters.  An extended analytical expression was proposed by Bažant et al. [12] to describe the combined deterministic-stochastic size effect for Type 1 structures.
l s is a stochastic scaling length, m is the dimensionless Weibull modulus (shape parameter of the Weibull distribution), and n is the number of spatial dimensions in which the structure is scaled. The proposed expression accounts for the stochastic and deterministic component of size effect.
The practical applicability of size effect laws is limited. This has been noted by Moallemi et al. [3]. Many of the parameters are empirical and there is no systematic method to determine the parameter values. A precise definition for the characteristic length is lacking and this ambiguity further undermines the general applicability of the size effect laws. This is discussed further in the following section. The nominal strength is influenced by multiple factors including the geometry of the structure, boundary conditions, and material properties. It is unlikely that it will be possible to develop a general size effect law that can be reliably applied to a wide range of non-trivial problems. Numerical methods offer the only practicable approach for addressing complex industry motivated problems.

Numerical methods
A significant body of work has been published on the numerical simulation of size effect in quasi-brittle materials, and models can be broadly categorised as deterministic or stochastic. In a deterministic model, the material properties are constant over the spatial domain. Deterministic models have been employed by many researchers including Moallemi et al. [3], Feng and Wu [13] and Barbat et al. [8]. In a stochastic model, the spatial variability of the material properties is considered. Stochastic models have been used by Gutiérrez and De Borst [14], Vorechovský [15], Bažant et al. [12], Yang and Frank Xu [16], Bobiński et al. [17], Grassl et al. [18], Syroka-Korol et al. [11], Syroka-Korol et al. [19], Zhou and Chen [20], Eliáš et al. [21] and Eliáš and Vorechovský [22]. These references are not exhaustive. A stochastic model is required for a complete examination of the mechanisms that govern the size effect, and allows the magnitude of the deterministic and stochastic components to be determined.
Capturing the correct energy dissipation as a crack forms and propagates is essential to accurately model structural size effect [23,8]. Many existing numerical methods use a localisation limiter that governs the width of the softening zone. Examples include crack band, gradient enhanced, non-local, and phase-field models. For further details, the reader is referred to [2,24,25,8]. The use of a localisation limiter is necessary to prevent damage from localising into a zone of zero volume [26]. The localisation limiter is chosen to ensure correct energy dissipation when a crack forms and is determined by an internal characteristic length l c that is a property of the material being modelled. Note that the correlation length in a probabilistic framework plays a role equivalent to the characteristic length [27]. Bažant and Pijaudier-Cabot [26] define the characteristic length as a material property that governs the minimum possible width of a zone of strain-softening damage. The characteristic length l c is the ratio of the fracture energy (energy dissipated per unit area) to the energy dissipated per unit volume. For concrete, the characteristic length is determined to be approximately 2.7 times the maximum aggregate size [26]. Nguyen et al. [28] states that a clear physical interpretation and direct link between the length parameter in the numerical model and the characteristic length of the material is debatable. Moallemi et al. [3] also note that the definition of the characteristic length is ambiguous. In different works it has loosely been interpreted as the size of the fracture process zone or the size of material inhomogeneities. The peridynamic model is non-local and material particles interact over a length-scale determined by the horizon δ. The peridynamic model differs from existing approaches as the size of the peridynamic horizon δ is not explicitly related to material properties.

Peridynamic theory
The limitations of existing methods are a result of the mathematical deficiencies of the governing partial differential equations. Silling [29] proposed a non-local theory of solid mechanics, known as the peridynamic theory, 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. 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 [30].
There are two primary formulations of the peridynamic theory: bond-based [29] and state based theory [31]. It is not the purpose of this work to explain the peridynamic theory in detail and the reader is referred to [29,30,32] for a detailed treatment of the theory. A mechanically intuitive but less rigorous way of obtaining the governing equations can be found in [33]. For a review of the peridynamic theory and its major applications, the reader is referred to work of Javili et al. [34].
For a review of the development of peridynamics for quasi-brittle materials, the reader is referred to Hobbs [6] and Hattori et al. [35]. Peridynamic models posses a number of advantageous features for modelling the failure of quasi-brittle materials: 1. The classical theory assumes that all forces are contact forces that act across zero distance (local theory). The peridynamic theory is a non-local theory in which material points interact with each other directly over finite distances. The fracture of quasi-brittle materials is a process in which non-locality is known to be important [36]. The failure of quasi-brittle materials is characterised by the accumulation of damage over a significant volume of material before localising into a discontinuity that generally propagates along a complex three-dimensional path.
2. The fundamental feature of any model capable of correctly capturing the structural size effect is the presence of some form of characteristic length (length scale) [7].

3.
A key advantage of using a bond-based peridynamic model is the simplicity of the constitutive model. In a standard continuum-based framework, constitutive models generally require numerous parameters that have no physical basis and must be determined using a calibration procedure. Vorel et al. [37] recently compared four state of the art concrete constitutive models in a standard continuum-based framework. The average number of model parameters is 25, and many of these parameters have no credible justification.

Numerical framework
A deterministic three-dimensional bond-based peridynamic model is employed in this paper.
It is not the purpose of this work to describe the numerical framework in detail, and the reader is referred to Hobbs [6] for a complete description of the model. The code used to generate the results in this study is available at https://github.com/mhobbs18/BB_PD.
All the results presented in this paper were obtained using an explicit scheme (as outlined in Fig. 4.14 of Hobbs [6]). Explicit schemes are the most popular and widely used method in the peridynamic literature due to their easy implementation and computational efficiency. All problems are discretised using the meshfree method proposed by Silling and Askari [30]. Note that the physical domain is discretised using a regular grid of nodes. Partial volume correction factors (to improve spatial integration accuracy) are calculated using the PA-PDLAMMPS algorithm introduced by Parks et al. [38]. Surface correction factors (to correct the perdynamic surface effect) are calculated using the volume method first proposed in Chapter 2 of Bobaru et al. [33].
The critical stretch correction scheme, first proposed by Hobbs [6], is employed for all computations. Hobbs [6] demonstrated that the application of surface correction factors produces a toughening effect. When surface correction factors are applied, the energy required to initiate a fracture at the free surface of a peridynamic body is greater than the material fracture energy G F and consequently leads to an overestimation of the peak load. For an illustration of the problem, the reader is referred to Fig. 5.11 in Hobbs [6]. The critical stretch s c should not be a constant across all bonds, and must be uniquely assigned depending on the applied surface correction factor. This can be easily achieved by using the corrected micro-modulus c (bond stiffness) to calculate s c .
Boundary conditions (supports and loading) are applied using a contact model. The contact algorithm described by Madenci and Oterkus [39] is used to model the interaction between the deformable body and a rigid support/impactor. The loading is applied using a displacementcontrolled cylindrical rigid impactor, and to reduce dynamic effects, the applied displacement is increased incrementally using a fifth-order smooth step function. To ensure that quasi-static conditions are met (i.e., dynamic effects are negligible), a general rule is to check that the kinetic energy in the system is less than 5% of the total internal energy. This approach has been employed for all simulations. Note that there is no damping in the system.

Constitutive model
The non-linear softening law, first proposed by Hobbs [6] and illustrated in Fig. 2, is used for all simulations. Note that the force-stretch (f -s) relationship of a peridynamic bond should be consistent with the macroscopic material response. The reader is referred to Hobbs [6] for a detailed examination of the different constitutive laws that have been proposed within the bond-based peridynamic framework and applied to the simulation of quasi-brittle materials.
The inclusion of softening behaviour is essential when modelling strain-softening materials. to linear decay, and k controls the rate of exponential decay. Note that the value of d will range from 0 to 1, where 0 indicates that the bond is still in the elastic range, and 1 represents a bond that has failed.
The critical stretch s c is defined by Eq. (5). For the full derivation, the reader is referred to Hobbs [6].
The proposed model provides an unambiguous approach for determining constitutive parameters and can be easily calibrated by adjusting k and α.

Mode I problem
Validation is performed against the full set of experimental results published by Grégoire et al. [40]. Geometrically similar notched and unnotched concrete members were tested in three-point bending to investigate size and boundary effects. The aim of the experiments was (1) to provide suitable experimental data for the validation of constitutive models, and (2) to provide a benchmark for assessing the accuracy of non-linear finite element models in capturing size effect. It should be noted that Hoover et al. [41] independently performed a similar experimental study that has also widely been used for validation. A schematic diagram of the experimental setup is illustrated in Fig. 3. Four different sizes of geometrically similar specimens were considered: depth d = 50 mm, 100 mm, 200 mm, and 400 mm; constant thickness of 50mm; and fixed span-to-depth ratio of 2.5. Member dimensions can be found in Table 1 for all test cases. This is essential in a robust validation procedure [43].  [8]. The authors employ a crack band model in a smeared crack framework. This is the most widely used approach for modelling the failure of concrete [46,8]. Comparing the predictive accuracy of the peridynamic model against established numerical methods is an essential step in validation. The above discussion is not exhaustive and numerical results have been reported by many other authors.

Results: notched cases
The predicted load-CMOD response for the full set of notched tests is illustrated in Fig. 5 (a) and (b). Results from Barbat et al. [8] are included for comparison. The numerical and experimental results are in good agreement. The peak-load and softening response is accurately captured for almost all notched specimens, and the correct fracture behaviour is predicted for all cases. The computed fracture paths for the half-notched specimens are illustrated in Fig. 4.
Results for the fifth-notched specimens have not been included to avoid repetition of similar results. For all problems, the fracture initiates from the mid-span of the beam and propagates vertically towards the applied load. The results demonstrate that a bond-based peridynamic model can capture a deterministic size effect. The statistical size effect is negligible in notched specimens [17]. The influence of randomness in material properties is lessened due to the high concentration of stress at the notch tip.

Results: unnotched case
The predicted load-CMOD response for the full set of unnotched tests is illustrated in Fig.   5 (c). As the structure size increases, so does the overprediction of structural strength for unnotched members. For Beam 3, the predictive error between the experimental and predicted failure load is approximately +10%; for Beam 2, the predictive error is approximately +20%; and for Beam 1, the predictive error is approximately +40%. To some extent, this behaviour was expected. For notched problems, statistical variability in material properties can largely be ignored, but for unnotched problems, the randomness of material properties has a significant effect on structural strength [11,22]. For unnotched specimens, the probability that a defect is present in the region of highly stressed material is much higher. Vorechovský [15], Bobiński et al. [17], Syroka-Korol et al. [11], and Eliáš et al. [21] found Syroka-Korol et al. [11] compared the stochastic and deterministic numerical results with size effect laws proposed in the literature [2]. The deterministic-stochastic size effect law of Bažant et al. [12] was fitted to the stochastic numerical results using a Weibull modulus of 48 (see Fig. 21 in [11]). The Weibull modulus m is a dimensionless parameter of the Weibull distribution, used to describe the variability in material properties. For concrete, the Weibull modulus is generally in the range of 5 -24, where a higher value signifies less variability [47]. Eliáš et al. [21] modelled a similar set of unnotched beams using a discrete meso-scale model.
Deterministic calculations were performed using a random spatial arrangement of particles and a small size effect was predicted, but the reduction in nominal strength was significantly less than that observed experimentally.
The above analysis does not provide any definitive conclusions and clearly the next step is to model the spatial variability in material properties. This is beyond the scope of the current work but implementation details are provided in the discussion. The influence of an unstructured discretisation also requires further study.
Computed crack paths for the unnotched specimens have not been included to avoid repetition of similar results. It is illustrated in the next section (Section 5.3 Energy dissipation) that a wide zone of softening damage develops in the tensile region prior to the peak load (see Fig.   7 (c) and Fig. 8). At the peak load, the damage localises, and a fracture initiates from the mid-span of the beam and propagates vertically towards the applied load.

Energy dissipation
This section provides an examination and discussion of the damage energy dissipation.
Stress redistribution and the associated energy release caused by micro-cracking is the primary mechanism influencing the size effect on structural strength [9]. Investigating the energy dissipation might help to explain why the numerical model fails to predict a deterministic size effect for the unnotched problems.
For every specimen, the width of the region of energy dissipation at peak load has been determined and the results are presented in Table 2. Note that the region of energy dissipation represents the distribution of micro-cracking. Similar to the work of Grassl et al. [18], the region of energy dissipation is determined using the mid-point of bonds that have exceeded the elastic limit. Fig. 7 illustrates the region of damage energy dissipation in Beam 4 (d = 50 mm) at peak load using a fine mesh (∆x = 1.25 mm). Each point represents the mid-point of a bond that has exceeded the elastic limit. The size of the region of energy dissipation at peak load depends on the depth of the beam and the boundary type.    [18] were later validated against acoustic emission data and a good agreement was observed [48].
For the unnotched beams, the width of the region of energy dissipation is strongly influenced by the depth of the beam, increasing from 62.5 mm for the smallest beam (d = 50 mm) to 553 mm for the largest beam (d = 400 mm). For the largest beam (d = 400 mm), Grassl et al. [18] determined that the region of energy dissipation is approximately 200 mm wide. This prediction was later validated against acoustic emissions data [48]. As illustrated in Fig. 8, the region of energy dissipation scales linearly with the depth. The width of the region of energy dissipation to beam depth ratio is approximately constant for all unnotched specimens and this is believed to be the source of the overprediction of peak load. It is hypothesised that a source of randomness must be introduced to trigger the localisation of damage in unnotched specimens.
This hypothesis requires further study. The reader is referred to Hobbs [6] for further discussion of energy dissipation in quasi-brittle materials.

Mixed-mode problem
To test the generality of the model, the mixed-mode fracture tests of García-Álvarez et al. [49] are simulated. Feng and Wu [13] noted that the overwhelming majority of numerical studies only considerer mode I problems. Accurately predicting mixed-mode fracture behaviour is generally considered to be a harder problem than predicting mode I fracture behaviour, and most problems of practical relevance will be subject to mixed-mode loading conditions. Eccentrically notched concrete beams in three-point bending were experimentally tested.
A schematic diagram of the experimental setup is illustrated in Fig. 9. Three different sizes of geometrically similar specimens were considered: depth d = 80 mm, 160 mm, and 320 mm.
The notch-to-depth ratio λ was fixed at 0.25, and three different notch eccentricities µd were considered: 0.0d, 0.3125d, and 0.625d. The material properties used for computations are as follows: cylindrical compressive strength f cm,cyl = 20 MPa; modulus of elasticity E = 33.8 GPa; tensile strength f t = 3.5 MPa; material fracture energy G F =125.2 N/m. The cylindrical compressive strength f cm,cyl and modulus of elasticity E were determined experimentally. The value of tensile strength f t was determined numerically by García-Álvarez et al. [49]. The specimens were cured for 720 days in a fog room. The long curing period and high moisture content explains the high tensile strength f t [42]. The material fracture energy G F was estimated using the following empirical equation [42]: G F = 73f 0.18 cm . Numerical results have been reported by García-Álvarez et al. [49] who used a non-linear finite element model with interface elements and a cohesive crack model to simulate the experimental tests. Interface elements are placed along the crack path, requiring pre-existing knowledge of the fracture behaviour and limiting this method to simplistic problems. Numerical results have also been reported by Feng and Wu [13], Mendonça et al. [50] and Barbat et al. [8]. The methods used in these works are more generally applicable and do not require pre-existing knowledge of the fracture path. Mendonça et al. [50] used the boundary element method to model the beam of depth 160 mm. Feng and Wu [13] and Barbat et al. [8] modelled the full set of experiments using a phase-field regularised cohesive zone model and a smeared crack approach respectively.
All previously reported results are two-dimensional. Results reported in this work were obtained using a three-dimensional model. Fig. 15 illustrates the predicted load-CMOD response for the mixed-mode fracture tests of García-Álvarez et al. [49]. The numerical results are in very good agreement with the experimental data. The peak load is correctly predicted in almost all cases and the softening response shows good agreement with the experimental response. Results from Barbat et al. [8] are included for comparison. It would be fair to say that the predictive accuracy of the peridynamic model is equivalent, if not superior, to the smeared crack approach employed by Barbat. The smeared crack approach, first introduced over 50 years ago by Rashid [51], in combination with the crack band theory [52], is the most widely used method for modelling the failure of concrete. The relative immaturity of the peridynamic model and the high predictive accuracy support the argument that this is a method worthy of further investigation.

Results
The correct fracture behaviour is predicted for all cases and the computed crack paths are illustrated in Fig. 10 and Fig. 11. Results for the mode I case (µd = 0.0d) have not been included to avoid repetition. For all mode I problems, the fracture initiates from the notch (mid-span) and propagates vertically towards the applied load. For the mixed-mode problems, the fracture initiates from the notch and propagates at an angle towards the applied load.    a single realisation of the random field. In many problems of practical interest, the mean failure load is of limited concern and a probability distribution of the load capacity is required.
Calculating structural failure probabilities requires sampling the entire probability distribution using Monte Carlo methods. Monte Carlo methods require thousands of simulations and are currently prohibitively expensive. The number of samples can be reduced by using variance reduction methods such as Latin hypercube sampling. Syroka-Korol et al. [11] examined the statistical size effect in concrete beams and used Latin hypercube sampling to reduce the number of samples required from thousands to less than a hundred. Multilevel Monte Carlo schemes, similar to those proposed by Blondeel et al. [54] and Dodwell et al. [55], offer another approach.
These schemes employ a hierarchy of mesh discretisations, with many samples taken on cheaper coarse meshes, and relatively few samples on expensive fine meshes.
Chen et al. [56] recently proposed a new model, referred to as the Intermediately Homogenized Peridynamic (IH-PD) model. The explicit modelling of the material microstructure is avoided and only essential details of the microstructure, such as void volume, are required by applying material properties directly to bonds. The authors modelled the failure behaviour of porous materials by deleting bonds to achieve a defined pore volume. Bonds are selected for deletion using a random uniform distribution. Mehrmashhadi et al. [57] adopted the same method and modelled fibre-reinforced composites by stochastically assigning bond properties to achieve a specific fibre volume fraction. This method has also been used by Zhao et al. [58] and Wu et al. [59] to model concrete. Bonds are randomly assigned properties to produce a defined volume fraction of aggregates and cement. Oterkus and Madenci [60] used a comparable method to model nuclear fuel pellets. Randomness was introduced into the model by assigning a Young's Modulus value to bonds using a random Gaussian distribution. These methods are computationally cheap but they lack a robust theoretical basis. By introducing randomness, realistic crack paths and crack tortuosity are predicted, but the failure probability distribution cannot be determined. Jones et al. [61] 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.
After addressing computational expense, there are further issues. Jones et al. [61] found that naively randomising the properties of bonds affects the results in unexpected ways, for example producing very narrow failure distributions. Accurately reproducing strength distributions, such as the Weibull distribution, requires a correction scheme. The authors proposed a scheme that successfully generates a Weibull strength distribution in a one-dimensional peridynamic body. The proposed method is only applicable to one-dimensional problems, and expansion to two and three-dimensional peridynamic problems is still required. The preceding discussion demonstrates that the inclusion of spatially varying material properties is not a trivial problem, and further work is required to reduce computational expense and to accurately reproduce strength distributions.

Load path following
The direct displacement-controlled loading scheme fails to capture the post-peak structural response for large unnotched beams, see Fig. 5 (c). This behaviour was expected, as direct displacement control is not capable of passing limit points of vertical tangent ('snap-back') [62]. The post-peak response of large structures is steeper than small structures, and for very large structures snap-back behaviour will occur [7,11,8]. Path-following techniques, such as the arc-length method [63], can be used to capture the full equilibrium path of materials that exhibit strain-softening behaviour. Sun et al. [64] have implemented an arc-length method within a peridynamic framework. Note that this requires the use of an implicit solver. Due to the simplicity, computational efficiency, and robustness of explicit schemes, implicit methods are rarely found in the peridynamic literature. The reader is referred to Hobbs [6] for a discussion of the advantages and disadvantages of implicit schemes.

Localisation limiter
The localisation of damage into a band that is one element wide is a major issue when using a finite element model. Results from finite element models lack objectivity, because upon mesh refinement, damage localises into a infinitely small volume and the energy dissipation tends to zero. Furthermore, capturing the correct energy dissipation as a crack forms and propagates is essential to accurately model structural size effect. A localisation limiter, defined by a characteristic length of the material, is used to enforce a minimum width of the strain localisation band [25]. This prevents the damage from localising into a zone of zero volume upon mesh refinement.
As stated by Bažant [7], the fundamental feature of any model capable of bridging the gap between the theory of plasticity and the theory of linear elastic fracture mechanics is the presence of some form of characteristic length. Due to the non-locality of the peridynamic model, a size effect is naturally captured. In the peridynamic model, the horizon δ defines the range over which particles interact. The peridynamic horizon δ is not an intrinsic material property.
This differentiates the peridynamic model from common numerical approaches, and the use of an ambiguous characteristic length parameter is avoided. The non-locality of the peridynamic model prevents both damage localisation into a zone of zero volume and zero energy dissipation upon mesh refinement. Note that a crack band model prevents zero energy dissipation upon refinement, but damage localisation still occurs in a single row of elements, even if the size of the elements becomes infinitely small [25]. The peridynamic model is physically consistent, and therefore is more fundamental than a crack band model.

Additional comments
The high computational expense prevented a mesh sensitivity study for all members. For completeness, the influence of the mesh resolution on the predictive accuracy should be examined.
Detailed mesh sensitivity studies will become increasingly practicable as solvers advance and computational power increases. The influence of an unstructured mesh discretisation also requires examination. It is hypothesised that an unstructured mesh might influence the localisation behaviour and improve the predictive accuracy for unnotched problems. Molinari et al. [65] used a cohesive zone model to show that by introducing randomness into the mesh structure, the total energy dissipated converges significantly faster. Finally, the mode II energy dissipation component in the mixed-mode fracture tests is minor [49], and further problems should be tested to improve the robustness of the validation.

Conclusions
Bažant [36] stated that correct modelling of the size effect on material strength should be adopted as the basic criterion of acceptability of any numerical model. This paper has demonstrated the capability of a bond-based peridynamic model to capture the size effect in quasi-brittle materials. This is the first time that a peridynamic model has been used to examine the size effect, and provides an important check on the validity of the numerical model. The key findings are as follows: 1. The peridynamic model accurately captures the deterministic size effect in notched problems and the predicted load-CMOD response and crack paths are in close agreement with the observed experimental data. Comparisons with leading numerical models are provided, and the predictive accuracy of the peridynamic model is equivalent, if not superior.
2. In conventional numerical models, capturing size effect requires a localisation limiter that governs the width of the softening zone. The localisation limiter is determined by an internal characteristic length that is a property of the material. The definition of the characteristic length is highly ambiguous and its value can not be determined experimentally. The non-local length scale in the peridynamic model is not related to material properties and avoids the aforementioned problems.
3. The peridynamic model fails to accurately capture the structural strength for all unnotched specimens. As the structure size increases, so does the overprediction of structural strength.
To some extent, this behaviour was expected. For notched problems, statistical variability in material properties can largely be ignored, but for unnotched problems, the randomness of material properties has a significant effect on structural strength.
4. The peridynamic model does not predict any deterministic size effect for the unnotched problems, and whilst it is expected that the inclusion of statistical variability in material properties would significantly improve the predictive accuracy, a small deterministic size effect should have been observed.

5.
To trigger the localisation of damage in unnotched specimens, it is hypothesised that a source of randomness must be introduced.
In 1999, Bažant [7] stated that discrete element models will be needed to help explain the mechanics of size effect and to separate the important processes from the unimportant ones.
A complete understanding of size effect is still lacking and peridynamic models might help to clarify the underlying mechanics. The bond-based peridynamic model could be used to inform a greater understanding of the physical basis of size effect, and extending the model to account for spatial variability in material properties is the obvious next step.