The assumptions made in design codes can result in unconservative predictions of shear strength for reinforced concrete members. The limitations of empirical methods have prompted the development and use of numerical techniques. A three-dimensional bond-based peridynamic framework is developed for predicting shear failure in reinforced concrete members. The predictive accuracy and generality of the framework is assessed against existing experimental results. Nine reinforced concrete beams that exhibit a wide range of failure modes are modelled. The shear-span-to-depth ratio is systematically varied from 1 to 8 to facilitate a study of di erent load-transfer mechanisms and failure modes. A comprehensive validation study such as this has until now been missing in the peridynamic literature. A bilinear constitutive law is employed, and the sensitivity of the model is tested using two levels of mesh reﬁnement. The predictive error between the experimental and numerical failure loads ranges from +3% to -57%, highlighting the importance of validation against a series of problems. The results demonstrate that the model captures many of the factors that contribute to shear and bending resistance. New insights into the capabilities and deﬁciencies of the peridynamic model are gained by comparing the expected load-transfer mechanisms with the predictive error. Keywords: Bond-based peridynamics, Reinforced concrete, Shear failure, Model validation, Three-dimensional fracture propagation 1. Introduction There is a pressing need to address the overdesign of reinforced concrete structures. The utilisation of structural concrete members is often low and structural material wastage in the order of 50% is common [ 1]. There is signiﬁcant scope for reducing material usage and associated carbon emissions by combing smarter structural design with innovative construction methods. Concrete structural elements are predominantly prismatic (cross-sectional area is constant along the member length). The reasons are twofold: prismatic reinforced concrete members are simple and economical to construct, and decades of experimental work has provided a comprehensive understanding of their behaviour. Using innovative construction methods, such as fabric formwork, it is possible to design and build e cient non-prismatic concrete members that use up to 40% less concrete than prismatic members of equivalent strength [ 2]. These material savings have been achieved by creating structural members whose geometry reﬂects the requirements of their loading envelope. Ensuring the safety and reliability of highly e cient structural members requires accurate analysis methods. Predicting the shear behaviour of reinforced concrete members is notoriously di cult and codiﬁed design methods are generally based on empirical formulas derived from the testing of prismatic beams. Non-prismatic members fall outside of conventional design codes and current methods for determining their structural response and ultimate limit state (ULS) behaviour are inadequate. These di culties have been highlighted in a number of papers. Rombach et al. [3]investigated the shear design of concrete members without shear reinforcement and compared results from existing experimental shear databases to design codes and numerical models. The authors concluded that no consistent Corresponding author Email addresses: mch61@cam.ac.uk (Mark Hobbs), gh465@cam.ac.uk (Gabriel Hattori), jjo33@cam.ac.uk (John Orr) Preprint submitted to Computers &Structures May 22, 2020 mechanical model exists for the design of reinforced concrete (RC) members without shear reinforcement. There are large discrepancies between experimental results and analytical models, indicating that existing approaches for shear design have a large amount of uncertainty. In an earlier paper, Rombach and Nghiep [4], investigated the shear design of variable depth beams (haunched beams) and noted that there was only a very limited number of investigations conducted regarding the shear capacity of variable depth beams. Orr et al. [5]further investigated the shear behaviour of non-prismatic RC beams, testing nineteen beams designed using three di erent approaches. The results showed that the assumptions made in some design codes can lead to unconservative predictions of shear strength for non-prismatic members. Shear failure is typically a sudden brittle failure, underlying the importance of correctly predicting the load carrying capacity and failure mechanism [ 6]. As our understanding of shear behaviour in reinforced concrete has improved, it has become clear that older structures do not always conform with modern standards. Extensions to the functional lifespan of a structure and increased loading requirements necessitate accurate structural assessments. It is crucial for the safety and serviceability of reinforced concrete structures that reliable design models are available. Empirical formulas are only applicable to the range of problems for which they were developed. Numerical models o er higher ﬂexibility and can overcome the severe limitations of empirical methods. A wide range of numerical methods are available for modelling the failure behaviour of reinforced concrete members. Computational approaches for modelling fracture can be broadly divided into three categories: discrete crack models, smeared crack models (continuum damage mechanics), and lattice and particle models [ 7]. It is di cult to determine the leading method due to a lack of validation studies that methodically examine a range of failure modes. For instance, Rabczuk and Belytschko [8], Rabczuk et al. [9,10]have developed meshfree approaches based on the cracking particle method. The accuracy of the developed model has been validated against a range of problems, including: plain concrete members, prestressed reinforced concrete beams, and concrete under explosive loading. Very good agreement with experimental data is shown. This validation approach is not su cient to make robust conclusions on the generality of the model, as according to Thacker et al. [11], validation is the process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended use of the model. The goal of validation is to measure the agreement between model predictions and experimental data from appropriately chosen tests. The intention of this paper is to develop a model that can predict the shear strength of reinforced concrete members; therefore, a systematic evaluation of the multiple parameters that inﬂuence shear failure is needed. Collins et al. [12] presented a di cult benchmark problem to evaluate the blind predictive accuracy of reinforced concrete models. Engineers were invited to predict the shear strength of a very deep beam. 66 entries were received from academia and industry. The overall winner of the contest used a smeared crack approach and a fracture mechanics-based cohesive crack model, implemented in the commercial ﬁnite element code ATENA. The shear strength predicted by the winner was 745 kN and the experimental value was 685 kN, a prediction error of +9%. ˇCervenka et al. [13] provide further details of the winning model. Shear strength predictions from the 66 participants were evenly distributed from 250 kN to 3773 kN, further demonstrating the challenge of predicting shear strength. Whilst numerical models have been applied successfully, accurately predicting the load capacity and failure mode for complex reinforced concrete members remains an open problem. Robust numerical models, capable of accurately capturing the structural behaviour of reinforced concrete elements of any geometry, under any loading condition are needed. The peridynamic theory of solid mechanics, introduced by Silling [14], provides a promising theoretical framework for developing robust numerical models, capable of capturing the complex ultimate limit state behaviour of reinforced concrete. Fracture is an emergent behaviour of the peridynamic governing equation, and no additional assumptions or techniques are required. The original bond-based theory limits the Poisson’s ratio to a ﬁxed value [ 14]. Silling et al. [15] later introduced a generalised state-based theory that overcomes the limitations of the original theory. Peridynamics has been successfully applied to a wide range of material failure problems and Javili et al. [16] provide a thorough review of the peridynamic theory and its many applications. There are a number of published papers that address concrete and reinforced concrete problems using a peridynamic model. Gerstle and Sau [17] were the ﬁrst to model concrete and reinforced concrete using a two-dimensional bond-based peridynamic model. Gerstle et al. [18] later introduced a micropolar peridynamic model that accounts for the bending sti ness of bonds. The model is applied to a reinforced concrete cantilever failing in shear [ 19] and a simply supported reinforced concrete beam, with and without shear reinforcement [ 20]. Numerical results are compared against predictions from design codes and a good agreement is observed. Signiﬁcant work remains to quantify the conﬁdence and predictive accuracy of peridynamic models. Diehl et al. [21] reviewed benchmark experiments for the validation of peridynamic models. They identiﬁed 39 publications 2 that compare numerical predictions from peridynamic simulations against experimental data. Peridynamic models have been assessed against four benchmark experiments that address the behaviour of concrete: reinforced concrete lap-splice problem, concrete beams in 3-point bending, Brazilian disk under compression, and anchor bolt pullout tests. Diehl notes that the majority of papers calibrate model parameters against a single experiment and only a few addressed further validation. For a thorough review of the development of peridynamics and applications to concrete and reinforced concrete problems, readers are referred to the work of Hattori et al. [22]. A three-dimensional bond-based peridynamic framework is developed for predicting shear failure in reinforced concrete members. It is a deliberate choice not to consider the more complex state-based theory and there are a number of reasons that justify this decision: (1) The simplest model should be examined ﬁrst. Only if the model is found to be deﬁcient should the complexity be increased. (2) The developed framework serves as a good baseline from which to measure any improvements obtained from more sophisticated models. (3) Bond-based models are computationally cheaper. This paper di erentiates itself by quantitatively examining the accuracy and generality of the model against a series of nine reinforced concrete beams tested by Leonhardt and Walther [23]. The chosen series of tests has been carefully selected to provide a robust and methodical validation of the bond-based peridynamic framework. The shear-span-to-depth ratio is systematically varied from 1 to 8 to facilitate a study of di erent load-transfer mechanisms and failure modes. The change in load-transfer behaviour and the important inﬂuence of fracture behaviour on the redistribution of stresses makes this a challenging benchmark. To the best of the authors knowledge, a robust quantitative examination such as this has until now been missing in the peridynamic literature. The paper is organised as follows: Section 2 reviews the bond-based peridynamic theory and the di erent constitutive models that have been proposed. Section 3 describes the numerical implementation and the displacement-controlled loading scheme. Section 4 provides a comprehensive validation study against a series of nine beams. Four validation metrics are used: load-deﬂection response, ultimate load capacity, failure mode, and fracture behaviour. Section 5 provides a discussion of the results and deﬁciencies within the model are identiﬁed. Section 6 concludes the paper and summarises the ﬁndings. Future research and development needs are listed. 2. Bond-based peridynamic theory Peridynamics is a non-local theory of solid mechanics, introduced by Silling [14]. The peridynamic model deﬁnes the 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 deﬁned by partial di erential equations. The classical theory is only valid if the body under analysis has a spatially continuous and di erentiable displacement ﬁeld. Spatial derivatives are not deﬁned across discontinuities and additional techniques are required for modelling fracture. The peridynamic theory does not include spatial derivatives and remains valid across discontinuities, allowing for the natural inclusion of fracture behaviour. This section brieﬂy presents the bond-based peridynamic theory. Detailed explanations can be found in the literature [14, 24]. 3 2.1. Peridynamic continuum model Figure 1: Peridynamic body and kinematics of particle pair and bond-based pairwise force function Assuming that a body occupies a spatial region R, for any material point x2R, a pairwise force function fcan be deﬁned to describe the interaction between particles within a ﬁnite distance ofx, at any time t, where urepresents the displacement of a material point (see Fig. 1). f=f(x;x0;u(x;t);u(x0;t);t);x02R:jjx0xjj6 (1) The peridynamic equation of motion for a single material point xat a point in time tis given by Newton’s second law of motion, and is deﬁned by Eq. (2). ¨u(x;t)=Z Hxf(u(x0;t)u(x;t);x0x)dVx0+b(x;t) (2) is mass density, ¨uis particle acceleration, bis body force per unit volume, and Hxis the neighbourhood of material point x. The size of the neighbourhood is deﬁned 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. Hx=Hx(x;)=fx02R:jjx0xjj6g (3) The pairwise force function frepresents the force that particle x0exerts on particle xand 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 fis deﬁned by Eq. (4). f(;)=f(j+j;)+ j+j(4) In a bond-based model, the force vector fis parallel to the deformed bond and the vector magnitude fis proportional to the bond stretch s. The initial relative position vector of a pair of particles is denoted by =x0x, and the relative displacement vector is denoted by =u0u. The current relative position vector is given by +. 2.2. Constitutive model The constitutive behaviour of a material within the bond-based peridynamic framework is expressed as a relationship between bond force and stretch. Bonds are essentially springs and they can possess linear or non-linear characteristics. Silling and Askari [24] introduced the ﬁrst peridynamic constitutive law, known as the prototype microelastic brittle (PMB) model. This model represents the behaviour of a linear elastic material and is depicted by Fig. 2. A brief 4 description of the linear elastic model is provided before discussion on its applicability to concrete problems and more suitable models are identiﬁed. The force magnitude fis deﬁned by Eq. (5). fis dependent only on the scalar distance between nodes in the deformed conﬁguration and does not depend on bond direction. The described material model is therefore isotropic. f(j+j;)=cs (5) The micromodulus constant cdeﬁnes the sti ness of a pair-wise bond and is often referred to as the bond sti ness or spring constant. ccontains all material speciﬁc information and can be correlated with classical material properties by equating the strain energy density of a homogeneous body under isotropic extension using the peridynamic theory and classical elasticity theory. For a pure shear state there is also an equivalence of strain energy, and a second expression forccan be determined. By equating the two expressions for cit is found that the Poisson’s ratio is limited to 1 /4 for the 3D case and the bond sti ness cfor 3D problems is deﬁned by Eq. (6). The value of the micromodulus constant is identical for all bonds within the horizon of a node and reduces to zero when the length of a bond exceeds the horizon length. c=12E 4(6) The scalar bond stretch sis deﬁned by Eq. (7). Bond stretch sis positive when a peridynamic bond is under tension and negative when under compression. s=j+jjj jj(7) A method for capturing failure must be included in the material model to enable the simulation of fracture. Damage is introduced by allowing bonds to break when a critical stretch value scis exceeded. Once a bond has broken, it can no longer sustain any tensile force. Damage is permanent and bonds cannot be repaired. The material model is consequently history dependent. The failure of bonds is captured by a history dependent scalar-valued function . When bond stretch is less than sc,is equal to 1, and when bond stretch is greater than sc,is equal to 0. (t;)=8>><>>:1if s (t;)<sc 0if s (t;)>sc(8) Silling and Askari [24] derived the critical stretch value in terms of the material fracture energy GF. A peridynamic body is completely separated into two halves and all the bonds that initially crossed the fracture surface are considered to have broken. By determining the energy required to break a single bond, the energy per unit fracture area for complete separation of a peridynamic body can be determined. The energy release rate GFis a material parameter and can be determined experimentally. The critical stretch value scfor a three-dimensional peridynamic material is deﬁned by Eq. (9). sc=r 5GF 6E(9) The damage 'of a peridynamic material at node xis deﬁned as the ratio of broken bonds to total number of bonds attached to node x. '(x;t)=1R Hx(x;t;)dVR HxdV(10) The restriction on the Poisson’s ratio has been addressed by the state-based peridynamic theory [ 15]. The bond- based model has also been revised to overcome the Poisson’s ratio restriction [ 18]. The Poisson’s ratio of concrete is generally between 0.1 to 0.2 and the limitation in the bond-based theory is unlikely to play a crucial role in the modelling of the ultimate limit state behaviour of concrete. This work uses the standard bond-based theory and the Poisson’s ratio is limited to 0.25. It is reasonable to trade o small improvements in predictive accuracy for the need to retain simplicity. 5 2.2.1. Concrete constitutive model This section presents and discusses the di erent constitutive models that have been proposed within the bond-based peridynamic framework and applied to the simulation of quasi-brittle materials. It should be noted that the compressive failure of bonds is not considered in any of the following models. The fracture behaviour of concrete at a micro and meso level is essentially a tensile phenomenon, even whilst under compression [ 25]. The linear elastic damage model, introduced by Silling and Askari [24], was proposed for capturing the behaviour of ideal-brittle materials, such as glass and ceramics, and has been used to accurately capture complex fracture patterns [ 26,27]. The linear model has also been applied to the modelling of concrete structures and has achieved reasonable results [ 28,29], but the use of linear damage models for concrete is primarily suited to the analysis of large structures (e.g. large dams with aggregate sizes up to 80 mm). Due to the size e ect law, the fracture behaviour of large concrete structures can be considered as brittle. This assumption can be made when the size of the fracture process zone (FPZ) is insigniﬁcant compared to the size of the structure [30]. Figure 2: Linear Damage Model In small and medium sized structures, concrete acts as a quasi-brittle material and exhibits strain-softening behaviour when in tension. The linear damage model does not consider softening behaviour, and this limits its ability to capture the true behaviour of concrete structures. Gerstle et al. [28] recognised this limitation and proposed a bilinear damage model, illustrated in Fig. 3. No details were provided for determining the value of the bond stretch at the linear elastic limit s0, or the value of the bond stretch at failure sc. Zaccariotto et al. [31] determined the parameters of the bilinear damage model with material fracture energy GFand experimental load-displacement curves. The energy required to start the fracture process is deﬁned as G0andkris the ratio sc=s0. The stretch at the linear elastic limit s0is deﬁned by Eq. (11) and kris determined from experimental load-displacement curves or through sensitivity studies. Examples of a concrete like material in three-point bending and a concrete four-point shear beam test are provided. For the concrete member, kris set at 20. All examples are two-dimensional and further validation of the bilinear damage model is required. s0=r 5G0 6Ekr=sc s0=GF G0(11) A bilinear damage model capable of reproducing fracture energy GFand crack nucleation strain independently is introduced in Chapter 2 of [ 32]. According to the original linear damage model, bonds break when their stretch s exceeds a critical value scthat has been calibrated to the material fracture energy GF. This implies that scis the crack nucleation strain. The problem, as pointed out in [ 32], is that the critical stretch scdepends strongly on the horizon and cannot be speciﬁed independently of GF. The fracture energy GFand crack nucleation strain can be deﬁned independently by adding a softening tail to the constitutive law. The crack nucleation strain is the linear elastic limit of a bond s0, and can be approximated by ft=E, where ftis the tensile strength, and Eis the elastic modulus. The critical stretch sccan then be deﬁned so that the energy release rate calculation agrees with GF(see Eq. 12). This model assumes that sc>s0, which imposes a constraint on the selection of the horizon and thus restricts the mesh spacing. sc=5GF 6Es0s0=ft E(12) 6 Ideally, the parameters in the constitutive model should be deﬁned independently of the horizon and the numerical results should remain consistent when the mesh discretisation is altered. Rossi Cabral et al. [33] recently proposed a bilinear model that introduces a material horizon related to the characteristic length of the material, and a numerical horizon that is chosen for convenience. The proposed model improves the ﬂexibility of the model calibration and reduces the computational cost. Validation is performed on three-dimensional sandstone specimens subjected to uniaxial tension. Figure 3: Bilinear Model Yang et al. [34] proposed the trilinear damage model, illustrated in Fig 4. They noted that the tension softening part of the bilinear model was not based on the true softening behaviour of concrete. A generalised method for determining the parameters s0,s1, and scwas proposed and only requires elastic modulus E, tensile strength ft, and fracture energy GF. The position of the kink point in the softening curve is determined with the initial fracture energy Gf, as proposed by Bažant [35]. Validation of the constitutive model is performed on two-dimensional concrete members in three-point bending and numerical results are compared to experimental results in the literature. Concrete members with di erent strengths, dimensions, and crack-depth ratios are examined. The bilinear and trilinear model show close agreement with experimental results. The linear damage mode signiﬁcantly overestimates the failure load of concrete members. This behaviour is observed in quasi-brittle materials (such as concrete), because fracture energy is predominantly the integral of the work consumed by the bond during the softening stage. This is e ectively the opposite of an ideal-brittle material (such as glass), where all the fracture energy is consumed during the elastic stage. Figure 4: Trilinear Model The critical stretch of a peridynamic bond is highly dependent on fracture energy GF. The selection of GFis a major source of uncertainty in the peridynamic model and the value chosen will have a signiﬁcant impact on the numerical results. Concrete material properties are typically determined through the testing of concrete compressive strength fc. Additional material parameters such as elastic modulus E, tensile strength ft, and fracture energy GFare rarely tested and must be determined from empirical formulas. At present, there is no standard theory relating concrete compressive strength to fracture energy and Neville [36] states that fracture mechanics has not succeeded in developing parameters that adequately quantify the resistance of concrete to cracking. Furthermore, to paraphrase the ﬁb[37]: deﬁning and measuring the fracture energy of concrete is still a highly debated issue, despite the availability of many test results. ﬁbModel Code 2010 does attempt to relate fracture energy to mean compressive strength fcmthrough the 7 empirical Eq. (13). This equation does not account for additional factors that are likely to play an important role such as maximum aggregate size. GF=73f0:18 cm (13) This work uses the bilinear damage model proposed by Zaccariotto et al. [31]. This model was chosen as it accounts for softening behaviour and is simple to implement. Validation on three-dimensional problems and mesh sensitivity studies are needed. 2.2.2. Steel constitutive model The microplastic material model introduced by Macek and Silling [38] is used to capture the yielding behaviour of steel reinforcing bars, see Fig. 5. Plasticity is introduced by assuming that a bond has a yielding stress. At an individual bond level the material is elastic perfectly plastic. At a macroscopic level, the material shows a strain hardening e ect due to bonds yielding at di erent values of bulk strain. In the microplastic model, the force magnitude fis deﬁned by Eq. (14), where sp(t) is the plastic stretch history deﬁned by Eq. (15), and syis the yielding stretch. f=8>><>>:c(ssp(t))ifjx0xj6 0 otherwise(14) sp(0)=0;d dtsp=8>><>>:0 ifjsspj<sy d dts otherwise(15) The bond yield stretch sycan be related to the ultimate yield strength ultby considering that all bonds will have yielded when the ultimate strength is reached. By following the same process used to derive the critical bond stretch in the original linear elastic material model, the bond yield stretch is deﬁned by Eq. (16), where Kis bulk modulus. syult 3K(16) Figure 5: Microplastic Model 2.2.3. Interface constitutive model The bond between the steel reinforcing bars and concrete plays an important role in the behaviour of a reinforced concrete member. Within the peridynamic framework there are a limited number of papers that explicitly discuss the steel-concrete interface model. Gerstle et al. [39] modelled a reinforced concrete lap splice problem and treated the steel-to-concrete interface bonds as modiﬁed concrete bonds. The calculated pairwise force and pre-deﬁned critical stretch value in the modiﬁed bonds are multiplied by a factor of three. This factor accounts for the toughening e ect of the ribs on the reinforcing bars. Lu et al. [40] simulate anchor bolt pullout in concrete by introducing a short-range interaction force to reproduce the interaction between the concrete and steel anchor bolt. Li and Guo [41] recently used a dual-horizon peridynamic model to simulate the debonding process in ﬁbre-reinforced polymer (FRP)-to-concrete bonded joints. A multi-scale discretisation scheme is used and the mesh is reﬁned at the bond interface. Using a linear elastic damage model, the captured interfacial slip is in accordance with experimental results. 8 The interface behaviour has been investigated by many researchers using conventional numerical methods. Jendele and ˇCervenka [42] presented a bond slip model for simulating the interface between reinforcing bars and concrete within a non-linear ﬁnite element model. They found that when bond slip is neglected, a ﬁner mesh is required to capture the correct behaviour. More recently, Xenos and Grassl [43] modelled the failure behaviour of a reinforced concrete member with non-local and crack band approaches using a damage-plasticity model and they assumed perfect bond between the steel reinforcement and concrete. This idealisation of the interface behaviour leads to reasonable results and simpliﬁes the problem. This work also assumes perfect bond between the reinforcement and concrete and bond slip is not accounted for. It is not the purpose of this paper to study the interface in detail. The sti ness of steel-to-concrete bonds are set equivalent to concrete bonds and bond failure can not occur. 3. Numerical framework The peridynamic equation of motion contains both integrals and derivatives. Equations of this type are known as integro-di erential equations and are generally hard to solve analytically. Most problems of interest to engineers concern complex geometries and consist of multiple materials. Numerical techniques must be used to obtain solutions. A bond-based peridynamic analysis code has been developed in MATLAB and C. The code makes use of shared-memory parallelism using OpenMP. It is the intention of the authors to make the code open source in the near future. At a high level, a peridynamic analysis code consists of several distinct components: a routine for spatial integration /calculating bond forces, a time integration routine, a method for applying initial and boundary conditions, a damping method for obtaining steady-state solutions, a routine for applying surface corrections, a routine for bond damage evaluation and tracking, and a method for reading input data and saving output data [ 44]. Details of the chosen methods are described brieﬂy. 3.1. Spatial integration The meshfree approach of Silling and Askari [24] has been utilised in the majority of peridynamic simulations and has been implemented in this work. The original meshfree method is popular due to its simplicity and relatively low computational cost. In the meshfree method, the peridynamic body is discretised into a regular grid of nodes and every node is assigned a corresponding cell. The discretised form of the peridynamic equation of motion (Eq. 2) is given by Eq. (17). Where un i=u(xi;t=nt),nsigniﬁes the nth time step, tis the time step size, and prepresents any material particle in the horizon of node i. ¨un i=X pf(un pun i;xpxi)dVp+bn i (17) Spatial integration over the horizon of a node is implemented numerically with a summation of integrals over all cells contained within the horizon boundary. A one-point quadrature method is applied over every cell. The centre point of every cell is used as a quadrature point, and the area (2D problems) or volume (3D problems) of a cell is used as a weight. For some nodes, a fraction of the cell volume will sit across the horizon boundary, and using the whole cell volume will lead to errors in the spatial integration. A volume correction procedure can be implemented to improve the accuracy of the spatial integration. Seleson [45] provides a detailed analysis of the various methods available for volume correction. This work has implemented the volume correction procedure introduced by Parks et al. [46], referred to as the PA-PDLAMMPS algorithm. This method was chosen for its simplicity. 3.2. Time integration The peridynamic equation of motion can be solved using explicit or implicit time integration methods. Explicit schemes are robust and have been successfully applied in a wide range of peridynamic simulations involving large deformations and material failure. The main limitation of explicit methods is that they become numerically unstable for large time steps t. To ensure numerical stability, small time steps are required and this can signiﬁcantly increase simulation run time. This work has implemented the Euler-Cromer method, deﬁned by Eq. (18), and the stable time step condition derived by Silling and Askari [24] is used. ˙un+1 i=˙un i+ t¨un i un+1 i=un i+ t˙un+1 i(18) 9 3.3. Damping The peridynamic equation of motion is in a dynamic form and direct application to static /quasi-static problems is not possible. All materials exhibit damping to some extent and damping should be included within the peridynamic model. Without damping the solution will exhibit large oscillations. An additional local damping term is introduced into the governing motion equation, Eq. (19), where the local damping coe cient Cis a positive real constant with units of kg=m3s. ¨u(x;t)=Z Hxf(u(x0;t)u(x;t);x0x)dVx0+b(x;t)C˙u(x;t) (19) This method is sometimes referred to as dynamic relaxation. Local damping can cause unphysical e ects in fully dynamic simulations; however, dynamic e ects play a small role in the quasi-static behaviour of materials and structures. The choice of local damping coe cient Ccan have a large e ect on the convergence speed for quasi-static problems. Determining the most e ective damping coe cient Ccan be di cult in non-linear problems and this led to the introduction of an adaptive dynamic relaxation (ADR) method for quasi-static simulations using the peridynamic theory [ 47]. In ADR, the damping coe cient is changed adaptively in every time step and this can lead to a speed up in convergence. Preliminary studies found that standard dynamic relaxation works well for composite material problems and the use of ADR was not necessary. 3.4. Boundary conditions This work uses a displacement-controlled loading scheme. Displacement boundary conditions cannot be imposed on the boundary surface of a peridynamic body and must be applied directly to a material point. As noted by Mehrmashhadi et al. [48], applying imposed displacements to a single node will result in a large surface e ect. Imposed displacements should be distributed over a number of a nodes but this can also be problematic. Nodes subjected to the same applied displacement move together and forces in the connecting bonds are zero. To limit these issues, Mehrmashhadi applies imposed displacements over a two by two node region. In this work, loading is applied using a displacement-controlled cylindrical rigid impactor. The contact algorithm described by Madenci and Oterkus [49] is used to model the interaction between the rigid impactor and deformable body. This method appears to solve the aforementioned issues. In quasi-static problems, the rate of applied loading should be smooth to reduce dynamic e ects and improve accuracy. Discontinuities in the loading rate produce stress waves that introduce noise and errors into the solution. The applied displacement is increased incrementally using a smoothstep function. This work uses a ﬁfth-order smoothstep function deﬁned by Eq. (20) and illustrated in Fig. 6, where Atrepresents the applied displacement at time t. The function deﬁnes a smooth loading amplitude that ramps up from zero at the start of loading ( t0=0s,A0=0mm) to the ﬁnal value ( t1;A1). The ﬁrst and second derivatives of the curve are continuous and are zero at t0andt1. There is minimal change in acceleration from one increment to the next. This method is based on that used by Abaqus, a commercial non-linear ﬁnite element analysis tool. To ensure that quasi-static conditions are met, a general rule is to check that the kinetic energy in the system is less than 5% of the total internal energy. This approach is used in Abaqus for explicit analysis [50]. At=A0+(A1A0) 3(1015 +6 2) =tt0 t1t0(20) 10 Figure 6: Applied displacement against time using a smooth amplitude curve 3.5. Surface correction factors The derivation of the peridynamic bond sti ness cis based on the assumption that the horizon of a material point Hxis contained completely within the bulk of a body. This assumption does not hold true for nodes within a distance  of the body edge. Nodes within of the body edge do not possess a full non-local neighbourhood. Consequently the material properties of nodes contained within the main bulk and nodes within a distance of the body edge di ers. This is known as the peridynamic surface e ect. Nodes that do not possess a full non-local neighbourhood have less interactions (bonds) and this leads to a reduction in material sti ness. Material behaviour is not accurately captured at the free surface of a peridynamic body, and displacements are larger than that predicted by classical methods. Surface correction factors can be used to correct sti ness softening behaviour. This work has implemented the V olume Method, ﬁrst proposed in Chapter 2 of [ 32]. Le and Bobaru [51] provide a detailed comparison of the most common surface correction techniques. They found that the volume method is one of the most e ective techniques for correcting surface e ects and is relatively simple to implement. Sti ness softening behaviour is corrected by increasing the micromodulus of bonds connected to nodes within a distance of the body edge. This is done by ensuring that nodes under homogeneous deformation have the same value of strain energy density when located near the surface or within the main bulk. The micromodulus of a bond connecting node xandx0is corrected using Eq. (21). The horizon area /volume of node xisV(x)and the horizon area /volume of node x0isV(x0). The horizon area /volume of a node within the main bulk is V0. For 3D problems, V0=43=3; and for 2D problems, V0=2. The dimensionless sti ening factor will be greater than or equal to 1 because V0will always be larger than or equal to V(x) and V(x0). ccorrected =c=2V0 V(x)+V(x0)(21) 4. Model validation The goal of this paper is to develop and benchmark the capabilities of a three-dimensional bond-based peridynamic model for the simulation of failure behaviour in reinforced concrete members. Whilst peridynamic models have been studied fairly extensively, a rigorous validation study for reinforced concrete problems is missing. Limitations of previous studies include: (1) a dependency on qualitative validation through visual comparisons of fracture behaviour, (2) simulations are restricted to two-dimensions but the complex interaction between reinforcing bars and concrete is three-dimensional, (3) validation is limited to a narrow range of problems and fails to provide su cient insight into the generality of the model. Proving that a model can correctly capture all possible scenarios is not feasible but a suciently wide range of problems should be addressed so that the limitations of a model can be identiﬁed. The inﬂuence of the discretisation on the structural response and fracture behaviour is also investigated. The predictive capabilities of the peridynamic model are quantiﬁed by comparison with experimental results. The model is evaluated against the following criteria: load-deﬂection curve, ultimate load capacity, failure mode, and fracture behaviour. The generality of the model is assessed by simulating a comprehensive range of reinforced concrete members that exhibit distinct failure modes. Members failing in ﬂexure and shear are considered. Shear 11 failure is distinguished from ﬂexural failure by the development of inclined cracks. The ﬂexural behaviour of reinforced concrete members is well understood and the ultimate ﬂexural strength can be accurately predicted using analytical methods based on simpliﬁed rectangular stress blocks. Failure loads predicted by analytical methods match closely with experimental results and form the basis of modern design codes. Correctly capturing ﬂexural failure behaviour still remains an important test for any numerical model. The shear behaviour of reinforced concrete members is much more complex and predictions of shear strength rely on empirical methods that lack a robust theoretical basis [ 52]. The ultimate shear capacity is governed by the combined resistance to shear force o ered by (1) the uncracked compressive zone, (2) arch or direct strut action, (3) aggregate interlock, (4) dowel action, and (5) the residual tensile strength in the fracture process zone. The contribution of each action to the overall shear resistance is related to parameters such as shear span, beam depth, and reinforcement ratio [ 53,54]. There is no consensus on the relational theory between shear failure and the many inﬂuencing parameters. This uncertainty is highlighted by the variability that exists within di erent design codes. Bentz et al. [55] state that predictions of ultimate shear strength can vary by a factor of more that two when using di erent design codes, whilst the ﬂexural strength predicted by the same codes does not vary by more than 10%. Correctly capturing di erent modes of shear failure is an essential test for any numerical model due to the sudden and catastrophic nature of members failing in shear. 4.1. Stuttgart shear tests The peridynamic model is validated against a series of reinforced concrete beams tested by Leonhardt and Walther [23], commonly referred to as the Stuttgart Shear Tests. Ten beams subjected to four-point loading were tested to investigate the inﬂuence of the shear-span-to-depth ratio ( av=d) upon the shear strength of simply supported prismatic beams. Shear reinforcement was not considered so that the development of diagonal shear cracks could be studied. Extensive testing has shown that the shear-span-to-depth ratio is one of the most important factors a ecting the failure behaviour of reinforced concrete beams. Kani [56] proved experimentally that the load transfer behaviour changes at a value of av=dapproximately equal to 2.5. If the value of av=dis less than 2.5, beams develop an internal arch and shear capacity increases. For values of av=dgreater than 2.5, beam action becomes dominant and sudden and brittle failure in diagonal-tension is expected. Full ﬂexural capacity, Mf l, is not obtained for values of av=dbetween 1.5 to 5, The largest reduction in ﬂexural capacity is observed for av=d=2.5, and the moment at failure, Mu, is approximately 50% lower than the theoretical capacity. The observed reduction in ﬂexural strength due to shear is often referred to as ‘Kani’s valley of shear failure’ and is illustrated in Fig. 7. The wide range of failure modes, the multiple inﬂuencing parameters, the complex load transfer mechanisms, and the important role of fracture behaviour on stress redistribution make this series of tests a suitable challenge for model validation. Figure 7: Kani’s shear valley Four broad categories of failure mode can be deﬁned: shear-web failure ( av=d<1), shear-compression failure (1<av=d<2:5), shear-ﬂexural failure ( 2:5<av=d<6), and ﬂexural failure ( av=d>6) [57,58]. For RC members within the same category, the expected failure behaviour (sequence of events and fracture characteristics) will be approximately the same. The results section is divided into these categories and a detailed description of each failure 12 mode is provided within the discussion. It should be noted that providing an exact deﬁnition of the di erent shear failure modes that may occur in a reinforced concrete beam is not trivial and the terminology used within the literature is often conﬂicting. Fig. 10 illustrates the experimental fracture behaviour from the series of tests by Leonhardt and Walther [23]. This is one of the best examples available that clearly demonstrates the inﬂuence of av=don fracture behaviour. Beam 4, 5, and 6 were repaired on one side after the occurrence of shear failure and further loaded until failure occurred on the other side. A schematic diagram of the experimental setup is illustrated in Fig. 8 and member dimensions and material properties can be found in Table 1. The shear span, av, is deﬁned as the distance between the applied point load and the point of nearest restraint. The shear force remains constant throughout the shear span. The e ective depth, d, is deﬁned as the distance from the top compression ﬁbre to the centre of the tensile reinforcement and is constant across all members ( d=270mm). All members have the same cross sectional area (320 mm 190 mm) and reinforcement ratio (2.05%). Longitudinal reinforcement is provided in the form of two ribbed steel bars of 26 mm diameter (As=1062 mm2). The reinforcing bars were manufactured from high yield steel and possess the following properties: modulus of elasticity Es=208,000 MPa, yield strength fy=465 MPa, and yield strain sy=0:2%. Input parameters for the concrete constitutive model have been determined using empirical formulas from ﬁbModel Code 2010 [ 37], that relate the experimentally measured cubic compressive strength to common concrete properties. For every test, the mean cube strength fcm;cubewas determined from twenty test cubes 200 mm in size. Model Code 2010 relates the compressive strength of cubes 150 mm in size to cylindrical compressive strength, thus introducing a degree of error. The value of kr, the ratio between the critical stretch and the linear elastic limit, has been set at 25, and was determined using experimental load-deﬂection curves from Barr et al. [59]. The value of 25 is approximate and there is a large degree of uncertainty. Yang et al. [34] note that it is not easy to determine the exact point at which linearity is lost in a load-deﬂection curve and this is a weakness of the calibration method proposed by Zaccariotto et al. [31]. Figure 8: Experimental Setup. L1is the shear span av. Adapted from Leonhardt and Walther [23]. Table 1: Stuttgart Shear Tests - Experimental Parameters Designation L(m) L1(m) L2(m) av=d f cm;cube(MPa) fcm;cyl(MPa) GF(N/m) Ec(MPa) 1 0.90 0.27 0.36 1.0 35 28.6 133 30,500 2 1.15 0.40 0.35 1.5 35 28.6 133 30,500 3 1.45 0.54 0.37 2.0 35 28.6 133 30,500 4 1.70 0.67 0.36 2.5 35 28.6 133 30,500 5 1.95 0.81 0.33 3.0 35 28.6 133 30,500 6 2.55 1.10 0.35 4.0 35 28.6 133 30,500 7 3.10 1.35 0.40 5.0 36.5 29.6 134 30,900 8 3.60 1.62 0.36 6.0 36.6 29.7 134 30,900 10 4.70 2.16 0.38 8.0 35.4 28.9 134 30,600 13 All members have been modelled with a constant peridynamic horizon =3:14xand regular grid spacing x =11.5 mm. To study mesh sensitivity, Beam 3, 5, and 7 have also been modelled using a ﬁner mesh with regular grid spacing x=7.67 mm. The grid spacing has been chosen so that the correct area of steel can be modelled. Experimental load-deﬂection curves are available for Beam 3, 5 and 7. The concrete constitutive model parameters for the two levels of mesh reﬁnement are deﬁned in Table 2. The loading conﬁguration of Beam 9 is not consistent with the remainder of the test series and has been omitted from the simulations. Peridynamic model parameters remain the same across all tests. This is essential in a robust validation procedure. Adjusting model parameters to improve the agreement between numerical and experimental measurements does not constitute validation [11]. Table 2: Simulation Parameters. Stretch is a dimensionless scalar. MeshGrid Spacing Bilinear Model Parameters x(mm) s0 sc Normal 11.5 6.34E-5 1.60E-3 Fine 7.67 7.78E-5 1.95E-3 Figure 9: Discretised cross-section. The nodal spacing for the normal mesh is 11.5 mm, and for the ﬁne mesh is 7.67 mm. The area of steel Asis approximately 1060 mm2in both meshes. 14 Figure 10: Failure behaviour in beams tested by Leonhardt and Walther [23]. Refer to Table 1 for dimensions. 4.2. Results This section presents the numerical results. The experimental and predicted failure loads are summarised in Table 3. Fig. 11 plots the experimental and numerical shear force and moment at failure against av=d. The numerical fracture pattern at failure for all beams can be found in Fig. 12. Numerical fracture patterns using the ﬁne mesh are provided in Fig. 13. Load-deﬂection results for Beam 3, 5, and 7 are plotted in Fig. 14, Fig. 15, and Fig. 16. Interpretation of the results is saved for the discussion. 15 Table 3: Comparison of experimental results, design code predictions, and numerical predictions Specimen 1 2 3 4 5 6 7 8 10 Experimental P(kN) 777 520 294 161 118 118120 1311126 126194 1031 Failure Mode2s s s s s s s s f Eurocode 2 P(kN) 287 194 144 144 144 144 144 125 93 Failure Mode s s s s s s s f f % Error -63% -63% -51% -10% +22% +22% +20% -1% +1% Normal Mesh P(kN) 335 253 170 131 112 85 87 76 56 Failure Mode s s s s s s s f f % Error -57% -51% -42% -19% -5% -28% -28% -40% -41% Fine Mesh P(kN) - - 181 - 122 - 92 - - Failure Mode - - s - s - s - - % Error - - -38% - +3% - -23% - - 1Beam 7, 8, and 10 were tested twice using identical material properties 2Failure mode - shear (s), ﬂexural (f) Figure 11: Shear force and moment at failure against av=d 16 Figure 12: Numerical fracture pattern at failure using the normal mesh. The cross-section is taken through a reinforcing bar 17 Figure 13: Numerical fracture pattern at failure for Beam 3, 5, and 7 using the ﬁne mesh 4.2.1. Shear-web failure: Beam 1 Shear-web failure, sometimes referred to as as true shear failure, deep beam failure, or splitting failure, occurs when the shear span is less than the e ective depth, 0 <av=d<1. Diagonal shear cracks form in the web and propagate in an approximately direct line between the applied load and support. This is due to the splitting action of the compressive force between the load and support. The shear force cannot transfer across the crack, and members behave like a truss. The compression zones act as struts and the longitudinal reinforcement ties the struts together. Final failure can be caused by splitting, compression failure at the supports, or reinforcement anchorage failure due to the large tensile forces generated. Flexural cracking is limited. Beam 1 has a shear-span-to-depth ratio of 1 and lies at the boundary of shear-web and shear-compression failure. The experimental fracture behaviour is typical of shear-web failure. Beam 1 failed experimentally at a load of 777 kN. Failure is predicted at a load of 335 kN and a mid-span deﬂection of approximately 1.1 mm. The experimental mid-span deﬂection is not known. The percentage error between the actual and predicted failure load is -57%; the largest error observed across all tests. The potential reasons for this error are numerous and analysis is provided within the discussion. The numerical and experimental fracture behaviour are in reasonable agreement but there are clear di erences that might provide insight into the reasons for premature failure. Uniformly spaced ﬂexural cracks form in the span of constant bending. As the applied load is increased, symmetric cracks propagate vertically towards the point of load application. The numerical and experimental fracture behaviour then start to diverge. In the numerical model, splitting between the steel-concrete interface occurs. This behaviour was not observed in the experimental work. Final failure occurs with the formation of a diagonal crack in the web that propagates suddenly between the load and support. In the numerical model, propagation of the diagonal crack is stopped by the reinforcement. The numerical model correctly captures that no concrete crushing occurs in the compression zone. 4.2.2. Shear-compression failure: Beam 2, 3 Shear-compression failure typically occurs in the range 1 <av=d<2.5 and is characterised by crushing failure of the concrete at the point of load application. This mode of failure can occur suddenly and at a substantially lower load than that seen in shear-web failure. Flexural cracks will form initially before a diagonal crack starts from the last ﬂexural crack and turn gradually towards the applied load. With few exceptions, diagonal cracks are extensions of ﬂexural cracks. High vertical compressive stresses directly beneath the applied load halt the propagation of the diagonal crack. Internal forces are then redistributed and the member is able to carry additional load before compression failure ﬁnally occurs adjacent to the applied loads. Similarly, high vertical compressive stresses near the supports prevent the further propagation of splitting failure between the reinforcing bars and concrete. Fracture behaviour is generally symmetric with two diagonal cracks forming. The experimental failure behaviour observed in Beam 2 and 3 is typical of shear-compression failure. 18 The experimental failure load for Beam 2 is 520 kN. Failure is predicted to occur at a load of 253 kN and a mid-span deﬂection of 1.9 mm. The prediction error of -51% is smaller than the error observed in Beam 1. The experimental and numerical fracture behaviour are in reasonable agreement. The numerical model captures the crushing failure between the applied loads and symmetrical diagonal cracks form as expected. The location and angle of the diagonal cracks is incorrect, initiating closer to the applied load and at a steeper angle than that observed experimentally. Again, splitting between the steel-concrete interface occurs and propagates to the supports. This behaviour was not seen experimentally. High compressive forces beneath the applied load halt the direct propagation of the crack toward the point of load application. The experimental failure load for Beam 3 is 294 kN. The numerical model predicts that failure will occur at a load of 170 kN and a mid-span deﬂection of 1.8 mm. The prediction error is -42%. The predicted deﬂection at failure is similar for Beam 2 and 3. Experimental load-deﬂection results are available for Beam 3 and the numerical load-deﬂection behaviour in Fig. 14 is in good agreement with the experimental results until premature failure occurs. The peridynamic model accurately reproduces the structural sti ness in the initial stages of loading. The model captures the sti ness softening observed at the start of crack initiation and there is very good agreement between the numerical and experimental results until premature failure occurs. The mesh size does not inﬂuence the structural response during the initial loading period. The predicted structural response diverges moderately in the later stages of loading. Using the ﬁne mesh, the predicted failure load is 181 kN and mid-span deﬂection is approximately 1.8 mm. In Beam 3, crushing occurs between the applied loads but only a single diagonal crack forms. The degree of crushing in the model is less than that observed experimentally. Using the normal mesh, the location and angle of the diagonal crack is incorrect. The ﬁne mesh is more accurate in capturing the correct location and angle but errors remain. Splitting at the steel-concrete interface occurs when using the normal and ﬁne mesh. Beam 3 and 5 have been modelled by Malm [60] using two commercial ﬁnite element codes, ATENA and Abaqus. The best correlation was achieved with Abaqus and the results have been included in Fig. 14. The results have been included for completeness but a fair comparison is dicult. The material properties for Beam 3 and 5 are identical but Malm adjusts the simulation parameters to improve agreement with the experimental data. This is calibration not validation. Figure 14: Beam 3 load-displacement curve ( av=d=2) 4.2.3. Shear-ﬂexural failure: Beam 4, 5, 6, 7, 8 Shear-ﬂexural failure occurs in the range 2.5 <av=d<6.0. Flexural cracks form during the early stages of loading. As the loading process progresses, a critical diagonal shear crack initiates from a ﬂexural crack and travels towards the applied load. The diagonal crack rotates as the fracture propagates further into the compression zone. For members subjected to shear and bending, cracks will vary in slope from 90at the extreme ﬁbre (pure bending) to 45at the neutral axis (pure shear). For beams in the range 2.5 <av=d<4, the diagonal crack meets resistance as it propagates 19 upwards into the compression zone and progression of the crack may be halted. With a further increase in loading, the fracture gradually extends along a shallow slope until sudden failure occurs when the arch mechanism is no longer capable of sustaining the applied load, leading to failure. This mode of failure is sometimes referred to as diagonal tension failure. For beams that lie within the range 4 <av=d<6, additional ﬂexural cracks might form between the critical diagonal crack and the support and failure is typically very sudden. Leonhardt and Walther [23] remarked that the failure of Beam 7 and 8 was sudden. Fracture behaviour is typically asymmetrical and failure occurs at a load approximately half that observed in shear-web failure. Shear-ﬂexural failure is also characterised by splitting between the concrete and longitudinal reinforcement. The experimental failure behaviour exhibited in Beam 4, 5, 6, 7, and 8 is characteristic of shear-ﬂexural failure. Beam 4 has a shear-span-to-depth ratio of 2.5 and lies at the boundary of shear-compression and shear-ﬂexural failure. The experimental failure load is 161 kN, and the numerical model predicts that failure will occur at a load of 131 kN and mid-span deﬂection of 2.24 mm. The prediction error is -19%. The experimental fracture pattern demonstrates characteristics of shear-compression and shear-ﬂexural failure. The slope of the critical diagonal crack is approximately constant with limited crack rotation. Concrete crushing does not occur. The numerical model correctly captures that crushing failure does not occur in Beam 4 but the location and angle of the diagonal crack is incorrect. The diagonal crack initiates closer to the applied load and at a steeper angle than that exhibited experimentally. The crack travels under the applied load and into the span of constant bending, this behaviour was not observed in the experiment. Beam 5 failed experimentally at a load of 118 kN and the numerical model predicted the peak load to be 112 kN at a mid-span deﬂection of 2.5 mm. Final failure occurs at a mid-span deﬂection of 2.8 mm. The numerical prediction error is -5%. This is the best prediction achieved across all tests using the normal mesh. The di erence between the actual and predicted failure load reduces as the internal behaviour switches from arch action to beam action. The numerical structural response, plotted in Fig. 15, is in good agreement with the experimental load-displacement curve. The initial structural sti ness is captured correctly. As loading progresses, structural sti ness is moderately overestimated. Concrete softening and the initiation of damage reduce the sti ness of the structural member and good agreement is observed near the peak load. The mesh size does not impact on the structural response until the later stages of loading. The failure load predicted using the ﬁne mesh is 122 kN and mid-span deﬂection is approximately 2.9 mm. Beam 5 has previously been modelled by Most [61], Malm [60], Jendele and ˇCervenka [42],ˇCervenka and Papanikolaou [62], and more recently by Xenos and Grassl [43]. Xenos and Grassl [43] modelled the failure of Beam 5 with non-local and crack band approaches using the damage-plasticity model CDPM2, originally introduced by Grassl et al. [63]. The results for the non-local model with a ﬁne mesh from Xenos and Grassl [43] are included in Fig. 15. The member was modelled in two-dimensions and planes of symmetry were exploited, with only half of the beam considered. This is not strictly valid due to the asymmetrical failure behaviour of reinforced concrete members in the range 2.5 <av=d< 6.0. The two models show good agreement until approximately 100 kN where the results diverge and the model of Xenos and Grassl [43] underestimates the peak load. The peridynamic model reproduces the expected asymmetrical failure mode and the fracture pattern is in reasonable agreement with the experimentally observed behaviour. Using the normal mesh, the location and angle of the diagonal crack is incorrect. The distance between the fracture origin and applied load is only 0.2 m. Experimentally this distance is measured to be approximately 0.5 m. Using the ﬁne mesh improves the fracture results. As the diagonal crack propagates into the compression zone, the characteristic rotation behaviour is seen and the distance between the fracture origin and applied load is over 0.3 m. 20 Figure 15: Beam 5 load-displacement curve ( av=d=3) Beam 6 fails experimentally at a load of 118 kN and the numerical model predicts that failure will occur at a load of 85 kN and mid-span deﬂection of 4.5 mm. The prediction error is -28%. This is a notable increase in the error of the numerical prediction. The numerical model captures the asymmetrical fracture behaviour but the location of the diagonal crack is wrong. The diagonal crack initiates closer to the applied load than that observed in the physical test. This behaviour has been observed in Beam 2, 3, 4, and 5. Beam 7 was tested twice using identical material parameters. The ﬁrst specimen failed experimentally at a load of 120 kN and the second specimen failed at 131 kN. The results demonstrate the inherent variability of ultimate strength in identical reinforced concrete elements. Herbrand [54] notes that the exact replication of ﬂexural shear tests is unreliable and the coe cient of variation of repeated tests is typically between 6 to 12%. Failure is predicted at a load of 87 kN and a mid-span deﬂection of approximately 8.1 mm. Using the lower failure load of 120 kN, the prediction error is -28%. The numerical and experimental load-deﬂection behaviour, plotted in Fig. 16, are in good agreement during the early stages of loading. As the applied load is increased, the model predicts a sti er response than the experimental results. After damage, the numerical model overestimates the reduction in member sti ness. The member carries additional load until premature failure occurs. The predicted structural response using the normal and ﬁne mesh is in good agreement until the peak. The peak load is mesh dependent and the failure load predicted using the ﬁne mesh is 92 kN and mid-span deﬂection is 9.1 mm. Herbrand et al. [64] modelled the failure of Beam 7 using a plastic damage material model with Abaqus. Full bond is assumed between the concrete and reinforcement. The two models show good agreement during the initial stages of loading. The model of Herbrand et al. [64] better captures the structural sti ness after the onset of damage. The peridynamic model captures the expected asymmetrical failure behaviour and the fracture behaviour is in good agreement with the experimental results. Fig. 17 illustrates the three-dimensional fracture path. Using the normal and ﬁne mesh, the location and angle of the critical diagonal crack is approximately correct. The predicted distance between the fracture origin and applied load is 0.45 m. The characteristic rotation behaviour of the shear crack is captured and the fracture propagates at a shallow angle beneath the applied load. The ﬂexural crack width spacing also shows good agreement with the experimental results. In the experimental test, additional ﬂexural cracks form between the critical diagonal crack and the support. This behaviour is only observed when using the ﬁne mesh. Fig. 18 displays the deformed shape of Beam 7 at failure and illustrates the opening of the diagonal crack. 21 Figure 16: Beam 7 load-displacement curve ( av=d=5) Figure 17: Fracture paths in Beam 7 using the ﬁne mesh Figure 18: Deformed shape of Beam 7 at failure. Displacement is scaled up by a factor of 10 Beam 8 has a shear span-to-depth ratio av=dof 6 and lies at the boundary of shear-ﬂexural and ﬂexural failure. The experimental failure behaviour is typical of shear-ﬂexural failure. Beam 8 was tested twice using identical material 22 parameters and both specimens failed at 126 kN. Failure is predicted at a load of 76 kN and a mid-span deﬂection of 11 mm. The prediction error is -40%. The numerical model fails to capture the correct failure mode and the fracture behaviour is typical of ﬂexural failure. Uniformly spaced ﬂexural cracks form and grow proportionally with increases in applied load. Compressive failure of the concrete occurs between the applied loads resulting in a complete loss of load carrying capacity. Yielding of steel reinforcement does not occur. 4.2.4. Flexural failure: Beam 10 Members with a shear-span-to-depth ratio greater than 6 are expected to reach their full theoretical ﬂexural capacity. The failure behaviour of a slender beam under constantly increasing load is well known. Narrow ﬂexural cracks will begin to appear well before the ultimate load. If proper bond is provided, the steel and concrete will attain the same strain. As the loading is increased, the width and length of the ﬂexural cracks will increase and the area of the compression zone will decrease. This behaviour will become more pronounced after yielding of the reinforcing steel. The area of the compression zone will continue to decrease and the compressive stresses will eventually reach a level where crushing failure of the concrete occurs [65]. Beam 10 was also tested twice. The ﬁrst specimen failed at 94 kN and the second specimen failed at 103 kN. Both specimens failed in ﬂexure. The numerical model predicts failure at 56 kN and a mid-span deﬂection of 18.4 mm. Using the lower failure load of 94 kN, the prediction error is -41%. The numerical fracture behaviour is in good agreement with the experimental results. The average crack spacing in the experimental test is approximately 110 mm and the maximum ﬂexural crack length is 230 mm. In the numerical model, the average crack spacing is 120 mm and the maximum ﬂexural crack length is 218 mm. Leonhardt and Walther [23] note that Beam 10 fails due to the crushing of concrete before the yield point was reached in the reinforcement. This behaviour is observed in the numerical model and yielding of the reinforcement does not occur. 5. Discussion 5.1. Shear-transfer actions The numerical results are in fair agreement with the experimental results but the ultimate load capacity is under- estimated in all cases. Insight into the behaviour of the numerical model can be gained through examination of the variation in error for di erent values of av=d. The prediction error between the numerical and experimental results decreases as av=dchanges from 1 to 3 (-57% !-5%). Unexpectedly the prediction error then increases as av=dchanges from 3 to 8 (-5%!-41%). The potential reasons for these di erences are numerous and it is probable that no single issue is solely responsible for the observed errors. The results suggest that the peridynamic model captures some but not all of the factors that contribute to shear resistance. Isolating sources of discrepancy between experimental and numerical results is challenging. By considering how the contribution of di erent shear-transfer actions changes for a range of av=dvalues, and relating to variations in prediction error, it is possible to identify speciﬁc areas of deﬁciency within the model. A number of recent papers have examined the contribution of the various shear-transfer actions in members without web reinforcement [ 53,54,66,67]. These works are used to guide the discussion but it should be noted that there is no clear consensus on the exact relationship between the individual shear-transfer mechanisms and inﬂuencing parameters. The largest prediction errors are observed in members where load transfer occurs predominantly through arch action ( av=d<2.5). This is not completely unexpected as this behaviour is observed within existing design codes and numerical models. The degree of uncertainty in shear strength predictions increases as av=dapproaches 0 and the ultimate shear strength of a deep beam ( av=d<1) can be four times greater than that predicted by Eurocode 2 (EC2) [3].ˇCervenka et al. [68] used non-linear ﬁnite element analysis to simulate 33 experimental tests of reinforced concrete elements and found that model uncertainty is higher for members failing in shear. Fig. 11 shows that the numerical model does capture an increase in shear capacity for members in the range av=d<2.5 but the degree of increase is smaller than expected. It is hypothesised that the failure to capture the correct ultimate load and fracture behaviour could be due to (1) limitations of the concrete damage model, and /or (2) the assumption of perfect bond between the steel reinforcement and concrete. Potential shortcomings of the concrete damage model include the the absence of compressive softening, and diculty capturing mode II fracture. It might be necessary to include compressive softening behaviour in the concrete 23 constitutive model. Fig. 19 illustrates the principal compressive strain distribution in Beam 1 before the onset of damage. Bottle-shaped compressive struts appear to be forming correctly but damage does not form within the strut as expected. Claus [69] found that the inclusion of compressive softening is important for the accurate representation of the formation and degradation of the compressive struts within a non-linear ﬁnite element model. In shear-compression and shear-ﬂexural failure, the critical diagonal crack initiates from a ﬂexural crack (mode I type fracture). In shear-web failure, the critical diagonal crack is expected to initiate within the web and the exact fracture mode or combination of modes is not clear but mode II type fracture potentially plays a signiﬁcant role. The prediction error between the numerical and experimental results is highest for Beam 1 ( av=d=1), and this might imply that the concrete constitutive model does not properly account for mode II fracture behaviour. Ren et al. [70] proposed a new damage rule based on the maximal deviatoric bond strain. They report that this damage rule is compatible with the three modes of brittle failure, and compared with the classic maximum bond stretch rule, the proposed damage model can predict in-plane and out-of-plane shear failure with ease. This requires further investigation. Figure 19: Principal compressive strain distribution in Beam 1. The purpose of this ﬁgure is to illustrate the formation of compressive struts and the inclusion of strain values is considered to be redundant. The load-deﬂection results for Beam 3 in Fig. 14 provide useful insight into the numerical behaviour of a beam under arch action. The peridynamic model accurately reproduces the structural response in the initial stages of loading and correctly captures the reduction in sti ness observed at the start of crack initiation. The results are in good agreement until premature failure occurs. This suggests that the model fails to capture an important aspect of structural behaviour. Fenwick and Pauley [71] found that arch action can only occur at the expense of slip (complete loss of bond transfer). Swamy et al. [72] also showed that arch action action can be developed in a beam only when the steel is fully unbonded between strong end anchorages. It was found that arch action enhances the shear resistance of an unbonded beam beyond that of a bonded beam failing in shear but also increases deﬂection. This evidence suggests that the assumption of perfect bond between the reinforcement and concrete is a likely source of error for beams that transfer load primarily through arch action. Work is needed to develop a constitutive model for bonds that interface between reinforcement and concrete. Beam 5 ( av=d=3) is the most brittle member tested (lowest moment capacity) and lies at the bottom of Kani’s shear valley. Fig. 11 plots the ultimate moment at failure from the experimental and numerical tests. It is interesting to observe that the numerical shear valley is much shallower than the experimental shear valley. The predicted ultimate moment capacity ranges from 44 kNm to 62 kNm, whereas the experimental ultimate moment capacity ranges from 48 kNm to 105 kNm. The predicted failure load for Beam 5 is also the most accurate across all tests and the key question is why? The shear force resistance of the inclined compressive strut and aggregate interlock is expected to be limited. The contribution of the inclined compressive strut is vanishingly small in members with a shear-span-to-depth ratio greater than 3. This is well known. The exact contribution of aggregate interlock to the shear strength is much less certain but the following points suggest that it is comparatively low: (1) Large aggregate interlock stresses are expected when the angle of the diagonal crack is greater than 45[66]. For Beam 5, the angle of the diagonal crack is approximately 28at the steepest point. (2) Aggregate interlock is much more pronounced in deep beams [ 53]. It is assumed that the shear force is resisted primarily by a combination of the uncracked compressive zone, residual tensile stresses at ﬂexural crack tips, and dowel action. For Beam 7, Herbrand [54] determined that approximately 30% of the shear force at failure is resisted by the uncracked compressive zone. As the contribution of the uncracked compressive zone is mainly governed by the reinforcement ratio , a similar value is expected for Beam 5. The contribution of the residual tensile stresses at ﬂexural crack tips to the shear strength is estimated to be approximately 20 - 40% of the 24 shear capacity [ 53,54], and the contribution of the dowel action to the shear strength is estimated to be approximately 10 - 20% of the shear capacity [ 53,67]. The high predictive accuracy suggests that the peridynamic model does a good job of capturing these shear-transfer actions. According to ˇCervenka et al. [68], the uncertainty in a non-linear ﬁnite element model of a reinforced concrete member decreases as the failure mode changes from shear to ﬂexural. Unexpectedly, the prediction error in the peridynamic model increases as the failure mode transitions from shear to ﬂexural. It is theorised that for members in the range 3 <av=d<6, discrepancies between the experimental and numerical results can be partly accounted for by the neglect of crack friction within the numerical model. Aggregate interlock is known to have an important role in the transfer of shear forces. As a crack opens, sliding occurs and contact between the faces produces normal and tangential stresses. Huber et al. [53] found that the crack opening and sliding behaviour is directly correlated with the shape of the critical shear crack, and aggregate interlock is more pronounced in members with steeply inclined shear cracks. The steep inclination of the shear crack at mid-span in Beam 7 and 8 (see Fig. 10) produces a large aggregate lock e ect. Cavagnis et al. [67] used digital image correlation to investigate the role of various shear-transfer actions in members without web reinforcement. It was found that for members governed by beam shear-transfer actions, the contribution of aggregate interlock is dominant. Herbrand [54] developed a mechanical model to predict the shear strength of reinforced concrete members without web reinforcement. Herbrand used the model to determine the contribution of each shear transfer action to the overall shear resistance of Beam 7 at every load step. At failure, aggregate interlock carries 43% of the shear force. Unbroken bonds that cross between fracture faces introduce a degree of resistance to crack sliding and the numerical prediction error is -28%. Contact between crack faces sti ens the structural response and the reduced numerical sti ness observed in Beam 7 (see Fig. 16) further supports the theory that crack friction must be accounted for. The numerical fracture pattern predicted by the peridynamic model (see Fig. 17 and 18) is superior to that of Herbrand et al. [64] and it is interesting to note that the asymmetrical failure is captured with homogeneous material properties. Correctly predicting the ultimate ﬂexural capacity of a reinforced concrete member is a key test when evaluating model robustness. The principles governing ﬂexural strength are clearly understood and accurately predicting the bending capacity of a reinforced concrete member is trivial. The size of the prediction error between the actual and predicted failure load for Beam 10 is somewhat unexpected. Unfortunately there is no experimental load-displacement graph from which insight into the numerical behaviour might be gained. The numerical fracture behaviour is in good agreement with the experimental results and the failure is governed by concrete crushing in the compression zone. This suggests that early compressive failure might be responsible for the disagreement between the numerical and experimental results. According to Eurocode 2 (EC2) [ 73], concrete reaches its maximum compressive stress when compressive strain is approximately 0.002 and crushing failure occurs at an ultimate strain of 0.0035. The maximum compressive strain in Beam 10 reaches 0.0013 shortly before failure occurs. This value is consistent with the prediction error. Beam 8 failed experimentally in shear but the numerical fracture displays the characteristics of ﬂexural failure. Crushing failure of concrete in the compression zone occurs before a diagonal shear crack can form. The numerical failure mechanism and prediction error in Beam 8 and 10 is almost identical. 5.2. Loading system Di erences in obtained structural behaviour and failure load can also be partly explained by the loading system. Leonhardt and Walther [23] tested all members with a discrete load-controlled system. Beams were loaded in incremental stages of approximately 1 /10th of the predicted failure load. In discrete load-controlled testing, the precise failure load and corresponding displacement are not measured. In the numerical simulations, all members were loaded with a forced displacement. Displacement-controlled testing has some important advantages. The exact failure load is identiﬁed and the structural response of the member can be studied continuously, including the post-peak behaviour. To explain the problem; the ﬁnal loading step prior to failure for Beam 3 was 235 kN, and shear failure occurred between this load step and the ﬁnal load step of 294 kN. The exact failure load is not known and the potential prediction error between the experimental and numerical results ranges from a minimum of -28% to a maximum of -42%. 5.3. Discretisation The mesh discretisation is another source of error and uncertainty. Mesh dependency has been investigated using two levels of mesh reﬁnement. The load-deﬂection curves are independent of the mesh reﬁnement until the peak load 25 is approached (see Fig. 14, 15, and 16). The structural response at the peak load is mesh dependent, and the ﬁne mesh provides better predictions of load capacity. The mesh reﬁnement inﬂuences the fracture patterns. The average crack spacing between ﬂexural cracks remains approximately the same for both meshes, but the location and path of the diagonal shear crack improves as the mesh is reﬁned. The improved shear crack location and path could explain the increased prediction accuracy. Bobaru and Hu [74] and Zhao et al. [75] note that the size of the horizon should be deﬁned by the smallest relevant geometrical feature. In this work the maximum horizon size should ideally be smaller than the diameter of the steel reinforcing bars. This was not reasonably practicable due to the high computational expense. This is a clear limitation and the interaction between the reinforcing bars and concrete might not be captured correctly. Computational expense and time limits prevented further examination of the mesh dependency. Ren et al. [76,70,77]recently introduced a dual-horizon peridynamic formulation which allows for an irregular discretisation. Using a multi-scale discretisation scheme, it would be possible to reﬁne the mesh at critical locations such as the bond interface, similar to the work of Lu et al. [40]. This is a promising area of future research. 5.4. Material parameters There is a signiﬁcant degree of uncertainty in the selection of material parameters and this will inﬂuence the reliability of the predicted failure load and behaviour. This uncertainty can be investigated by considering a probabilistic model of material properties. ˇCervenka et al. [13] modelled fracture energy GFas a random ﬁeld over the domain of a concrete member. The introduction of randomly distributed material properties reduced discrepancies between numerical and experimental results. Future work should consider the micro-structure and spatial variability of concrete. 5.5. Comparison with design codes Load capacity and failure mode predictions using EC2 [ 73] are provided in Table 3. All safety factors have been set to 1. The poor prediction accuracy demonstrates the uncertainty that surrounds shear failure. The peridynamic model o ers an appreciable improvement in accuracy for members failing in shear, and in all cases, the model provides safe predictions of load capacity. The design code predictions are unconservative for Beam 5, 6, and 7, and the predicted load capacity is up to 22% higher than that observed experimentally. EC2 predicts that Beam 8 will fail in ﬂexure at 125 kN. The predicted load capacity is correct but the failure mode is wrong. It is interesting that the failure mode in the numerical analysis agrees with EC2. Predicting the correct failure type for members that lie at the boundary of distinct failure modes is challenging. As has already been stated, EC2 provides accurate predictions of ﬂexural capacity. 5.6. Importance of further validation As pointed out by Diehl et al. [21], the majority of papers calibrate model parameters against a single experiment and only a few addressed further validation. This work has demonstrated the importance of further validation and areas of future research have emerged. The experimental test series by Leonhardt and Walther [23], modelled in this paper, is proposed as a benchmark problem for modelling the failure of reinforced concrete elements. To the best of the authors knowledge, this is the ﬁrst paper to address the full series of tests. Good benchmark problems help facilitate the fair comparison of di erent numerical models and provide the required credibility to make deﬁnitive conclusions. Validation against the complete test series will provide a robust investigation of the ability of a numerical model to capture a range of load-transfer and failure mechanisms. When viewing the load-deﬂection graphs in isolation it can be di cult to discern how the structural response di ers between members. Fig. 20 combines the experimental and numerical load-displacement curves for Beam 3, 5, and 7 using the normal mesh. The horizontal line indicates the load capacity predicted by EC2. This ﬁgure illustrates the capability of the peridynamic model to predict a range of structural responses. 26 Figure 20: Comparison of experimental and numerical load-displacement curves. The horizontal line indicates the load capacity predicted by EC2. 5.7. Summary The preceding results discussion has identiﬁed a number of potential areas of deﬁciency within the model and provided evidence to support the interpretation of results. There are multiple other factors that may introduce errors and detailed uncertainty quantiﬁcation studies are needed. The presented approach has demonstrated signiﬁcant promise for predicting the shear capacity of reinforced concrete members, and the peridynamic model possesses a number of features that justify further development: 1.Crack initiation, propagation and coalescence naturally emerge from the governing equations. Asymmetrical failure modes (as illustrated in Fig. 18) are predicted without explicitly modelling the microstructure. 2.Non-planar three-dimensional crack paths are captured with ease (see Fig. 17). This is important when considering more complex problems such as reinforced slabs. 3.The constitutive model that has been employed is very simple and has few parameters, and in some cases, as illustrated in Fig. 15, the results are better than that of well-established numerical methods. Xenos and Grassl [43] used the damage-plasticity model CDPM2, and they note that this model requires a large number of input parameters. V orel et al. [78] 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 justiﬁcation. The simplicity of the presented approach should be regarded as an advantage. 6. Conclusions Determining the shear capacity of reinforced concrete members is notoriously di cult. In this work, a three- dimensional bond-based peridynamic framework has been developed and validated for predicting shear failure in reinforced concrete beams. This is the ﬁrst study to comprehensively validate the predictive capability of a peridynamic model against a series of problems. The experimental test series by Leonhardt and Walther [23] was chosen as a benchmark due to the wide range of failure modes exhibited. It was hypothesised that this series of beams would be dicult to simulate due to the change in load-transfer behaviour as shear slenderness is increased, and the important inﬂuence of fracture behaviour on the redistribution of stresses. The predictive error between the experimental and numerical results ranges from +3% to -57%, highlighting the importance of validation against a series of problems. The results demonstrate that the model captures many, but not all, of the factors that contribute to shear and bending resistance. The main ﬁndings are listed: 1.The model successfully captures many of the key features of shear and ﬂexural failure. This is achieved using a simple bilinear constitutive model, a regular mesh discretisation, and without the need for prior knowledge of 27 fracture behaviour. The predicted fracture patterns are generally in good agreement with the experimental results. Using the ﬁne mesh, the predicted fracture behaviour is equivalent to, if not better than that found elsewhere in the literature. 2.The model accurately predicts the load capacity of Beam 5 and the numerical structural response is in good agreement with the experimental load-deﬂection curve. This suggests that the model captures the following load-transfer actions: shear transfer in the uncracked compressive zone, residual tensile stresses at ﬂexural crack tips, and dowelling action of the tensile reinforcement. 3.The assumption of perfect bond between the steel reinforcement and concrete is a likely source of error for beams that transfer load primarily through arch action. Arch action can only occur at the expense of slip and development of an interface constitutive model is needed. Factors such as bond-slip need to be considered. 4.The results for short members suggest that there are limitations with the concrete constitutive model. Potential shortcomings of the concrete damage model include the absence of compressive softening, and di culty capturing mode II fracture. It remains unclear to which degree the prediction error can be attributed to these factors but they are worthy of further investigation. The early compressive failure of concrete in slender beams could be linked to the identiﬁed shortcomings. It is also possible that there are unknown factors at play. 5.Discrepancies between the experimental and numerical results can be partly accounted for by the neglect of crack friction within the numerical model. For members governed by beam shear-transfer actions the contribution of aggregate interlock is dominant, and the numerical model should be expanded to explicitly account for crack friction. In addition, there a number of areas that have not been addressed by this work but likely play an important role: (1) the inﬂuence of size e ect on structural strength, and (2) the inﬂuence of concrete micro-structure on fracture initialisation and growth. Further challenges include addressing computational e ciency with a multi-scale mesh discretisation, and bounding predictions with upper and lower limits using uncertainty quantiﬁcation methods. The results presented in this study help to elucidate where future avenues of research lie for the application of peridynamics to reinforced concrete problems. With further research and development, it is conceivable that the peridynamic model might enable a better understanding of shear failure and the multitude of inﬂuencing factors. Acknowledgements This work was funded by the UK Engineering and Physical Sciences Research Council (EPSRC), grant no. EP/L016095 /1 - University of Cambridge Centre for Doctoral Training in Future Infrastructure and Built Environment, and grant no. EP /M020908 /1 - Concrete Modelled using Random Elements. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP /P020259 /1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk). Data access statement Data underlying this article can be accessed at https: //doi.org /10.17863 /CAM.52041. The dataset includes the raw numerical output used to create Figures 12 to 18.