A stochastic multiscale peridynamic model for corrosion-induced fracture in reinforced concrete

Concrete fracture caused by corrosion of reinforcing bars may cause subsequent structure failure. To better predict this process, we introduce a partially-homogenized stochastic peridynamic model with the simplest constitutive relation (linear elastic with brittle failure). The model links microscale information (phase volume fractions of mortar, aggregates, interfaces) to macroscale fracture behavior, while costing the same as a fully homogenized model. We show, and explain why a fully-homogenized peridynamic model fails to capture the correct concrete fracture modes/patterns, while the new model succeeds. The multiscale model predicts the evolution of fracture in reinforced concrete caused by corrosion products expansion in samples with a single or multiple rebars. Non-uniform expansion of corrosion products is enforced here as preset, incremental radial displacements. The computed fracture patterns and the order in which various cracks develop match what is seen in experiments. The model’s robustness is tested under different stochastic realizations and discretization grid types.


Introduction
Reinforced concrete is one of the most commonly used construction materials.Reinforced concrete structures exposed to a corrosive environment can be greatly affected as penetration of aggressive substances (water, oxygen, chloride, carbon dioxide, etc.) can lead to degradation of the steel rebars [1,2].If the concrete cover is thick and free of defects, this process usually is slow, extending over decades.However, concrete may contain or develop small cracks, caused during the manufacturing or early loading stages [3,4].Small amounts of aggressive elements then reach the rebars and can initiate their corrosion.When the corrosion reaches a certain stage (less than 1% of cross-section area of the initial rebar), the expansion of corrosion products (volume of corrosion products can be 3-4 times the volume of consumed iron [2]) can lead to major cracks propagating through the concrete cover.These major cracks can then act as additional channels for penetration of aggressive agents, accelerating the corrosion process [5].The corrosion of the steel rebar decreases its effective cross-sectional area and breaks the bond between the rebar and concrete, causing performance degradation of the structure [6].Moreover, changes in the layer affected by corrosion trigger significant reductions in ductility [7], in addition to hydrogen embrittlement (see pp. 334-336 in [8]).Corrosion of the rebar is usually non-uniform (e.g.pitting corrosion [8]).Nonuniformities in corrosion damage of the rebars can reduce the concrete cracking pressure by more than 50% compared with uniform corrosion [9].Degraded rebars can fail, resulting in the collapse of the structure.During this entire process, concrete's fracture plays a key role.Corrosion-induced concrete fracture is the focus of our paper.
Analytical methods based on the thick-walled cylinder theory are available to estimate conditions that would lead to cracking of the concrete cover due to rebar corrosion, but these cannot model the actual failure process, and, in general, are limited to a single rebar [10][11][12].Experimental investigations of corrosion-induced fracture in concrete can offer some insights into the process but they are expensive and time-consuming.Usually, these are performed using external currents to accelerate the corrosion process, making corrosion patterns more uniform than those resulting from natural environmental conditions.Departures from uniformity in the corrosion process, and the complex evolution of concrete fracture induced by rebar corrosion require the use of computational modeling to obtain a more complete understanding of this phenomenon [13,14].
An important contribution to computational modeling in this field has appeared in [13,14], where a 3D chemo-hygro-thermo-mechanical model for concrete (with a specialized constitutive model for concrete) is used to simulate corrosion-induced damage and transport of corrosion products into cracks.These works assumed that the reinforcement bar was already depassivated.Several choices for corrosion sites along the longitudinal direction and around the cross-section of rebars were tried.For some of these choices, the obtained fracture patterns agree with experimental observations very well.However, although the crack band theory used in these publications can alleviate mesh size dependence for smeared crack approach, it cannot solve the mesh orientation dependence [15].For each example, only one mesh was employed, and convergence studies were not presented for comparing variability of crack patterns with those from experiments.A 2D mechano-chemical model ( [16]) coupled the ingress of chloride ions, carbonation, electrochemical reaction and mechanical damage for the prediction of rebar corrosion and concrete damage.In this model, the active zone on the rebar evolves automatically.However, for the fracture model, using the crack band theory, no details (material properties, boundary conditions, etc.) were provided.Recently, a 2D diffusion-mechanical model [17] studied depassivation of the steel surface due to chloride ingress in concrete and the subsequent corrosion of steel and crack propagation in concrete.The distribution of rust thickness obtained is close to what is measured experimentally, but fracture patterns obtained for the 3rebar case do not capture the experimental observation very well.
One of the limitations of the above-mentioned works is that they all use a homogenized model for the concrete structure.This may not work well in cases where the concrete microstructure does play a role in how cracks initiate and grow.Around corroding rebars that create pressures against the concrete, microcracks develop leading to major concrete cracks reaching the concrete cover.The evolution of such cracks can depend on certain microstructure characteristics.Meso-scale models have shown their potentials for such problems [18].However, issues such as the selection of the geometric shape of the aggregates can significantly affect the fracture behavior [19][20][21].It is possible to use a meso-scale structure acquired from X-ray tomography images [22], but it is computationally costly to extend meso-scale models to meter-scale samples.Thus, the general application of these models is rather limited.A more desirable approach would be to develop a partially-homogenized model which implicitly involves some of concrete's microscale features, including their randomness.Such models, spanning multiple scales, would allow for efficient simulation of fracture and failure at the macro-scale while accounting for the correct crack initiation at the micro-scale.
While most of the above-mentioned models are based on the classical (local) continuum mechanics, nonlocal models offer some important advantages over local ones in modeling fracture behavior of heterogeneous materials like concrete ( [23]).Classical (local) continuum-based models lead to spurious mesh sensitivity in fracture problems, while nonlocal ones can prevent it.Another reason for using nonlocality is the complex interactions between microcracks: these appear at scales too small to efficiently model with a local formulation, and their formation and growth are interlinked, resulting in an effectively nonlocal damage behavior.Nonlocality is also necessary in a macroscale framework to describe microstructural phenomena in concrete such as cohesion, friction and aggregate interlock [24].Peridynamics (PD) is a nonlocal theory which has received considerable attention since its introduction almost two decades ago [25,26].PD reformulates the classical continuum mechanics by eliminating spatial derivatives to model mechanical [27][28][29][30][31][32][33][34][35][36][37], diffusion/corrosion [38][39][40][41][42][43], or mechano-chemical [44][45][46] etc., behaviors in materials involving damage.Using spatial integration rather than differentiation leads to a mathematically consistent formulation that works naturally for problems in which discontinuities in the domain (such as cracks) appear.In a PD model, cracks/damage can initiate and propagate autonomously [47].
One of the first applications of the PD theory to concrete structures was [48].A micropolar PD model to better simulate damage in concrete was introduced in [49].Later, this model was employed in [50] for simulating fracture in short fiber-reinforced concrete.The authors of [50] introduced a semi-discrete method to represent the fiber-concrete interaction and considered the random nature of concrete by reducing the particle strength of some PD nodes.A formulation of pressure-dependent PD plasticity model was shown to work well for compression, impact, and spallation of concrete structures [51].The trilinear peridynamic model introduced in [52] has shown good results in terms of the load-CMOD curves for threepoint-bending tests.This model, however, being a homogenized one, cannot capture the rough and tortuous crack trajectories.Tortuous crack paths are an indication of considerable local mode-mixity, and this is lost in some models [53].Homogenized models may also fail to capture the observed fracture modes in porous/composite materials [54,55].A mesoscopic PD model for concrete (using explicit geometrical representation of aggregates) was shown in [56], but when comparing to experiments, different normalization schemes are used for the simulation results and the experimental data, which raises questions about the validity of the comparison.
In this paper, we introduce a multiscale stochastic peridynamic model that implicitly uses some information about the concrete meso-scale structure, to simulate fracture induced in reinforced concrete induced by the expansion of rebar corrosion-products.The model does not require the explicit geometrical representation of aggregates, for example, and in that sense is a partially-or Intermediately-Homogenized peridynamic (IH-PD) model.Notably, the model uses the simplest possible linear elastic with brittle failure constitutive relation.
A mathematical distribution function is used to mimic the expansion process of the corrosion product, which is the loading that induces fracture in this setting.The results from this multiscale peridynamic model are compared with experiments and with results from a "fully-homogenized" PD (FH-PD) model, to highlight deficiencies of complete homogenization in modeling failure in concrete and the need for preserving some information about material heterogeneity.We test the model for concrete structures with a single and multiple rebars, for which experimental data is available in the literature.We also perform parametric studies to show how the aggregates' fracture energy and various possible rebar corrosion patterns can affect the evolution of fracture in reinforced concrete.This paper is organized as follows: in Section 2 we give a brief review of the bond-based peridynamic theory; in Section 3 we present the numerical discretization of the peridynamic formulation; in Section 4 we show the IH-PD model for concrete and its implementation; Section 5 discusses the radial displacement model for imposing the effective expansion of corrosion product as a boundary condition on the concrete in the hole where the rebar is; numerical results are gathered in Section 6, where we test the IH-PD model in corrosion-induced fracture in concrete structures with a single (several configurations) and multiple rebars; conclusions are drawn in section 7.

Brief review of bond-based peridynamic theory
Bond-based PD (BBPD) is the original version of peridynamics, later generalized as state-based PD (SBPD) [25,47].BBPD leads to material models with a fixed Poisson's ratio (1/3 in 2D plane stress problems, and 1/4 in 2D plane strain and 3D problems).The focus of our paper is on concrete's fracture behavior and Poisson's ratio has little effect on such problems (see theoretical analyses in [57,58]).It is then reasonable to use the BBPD model, and the results and conclusions should not be affected by this choice.Given the geometry of the experimental samples used to compare our simulation results with, we employ the BBPD for plane strain conditions (Poisson's ratio equals 1/4) everywhere in the paper.
The equations of motion for the BBPD can be written as [25]; where  is the density field,  is the displacement vector field,  is the pairwise force in the peridynamic bond  −  ̂, and  is the body force field.  is called the "horizon" of , and is the region in which pairwise forces exist between  and  ̂, an arbitrary point located inside   (see Fig. 1).  ̂ is the volume (area in 2D, length in 1D) occupied by the material point  ̂, and  is the time (or a parameter tracking the loading step in quasi-static problems, in which case the acceleration in the equation above is zero).The horizon is usually taken to be a ball (sphere in 3D, disk in 2D, line segment in 1D) of radius .The pairwise force for a prototype microelastic brittle material [47] is defined as: where  =  ̂−  is the relative position of  ̂ and  in the reference configuration,  = ( ̂, ) − (, ) is the relative displacement with respect to the reference configuration,  = ‖+‖−‖‖ ‖‖ = −  is the relative deformation or bond strain ( = ‖‖ and  = ‖ + ‖). and ℎ are respectively given by: ℎ(, ) = { 1 if (,  ′ ) <  0 for all 0 ≤  ′ ≤  0 otherwise where () is the micro-modulus function or the elastic stiffness of the bond.The micro-modulus function can take different forms, depending on the required horizon-scale behavior [59].Here we only consider plane strain conditions and the "conical" micromodulus function [60]: For heterogeneous materials, micro-modulus depends on the location.This is discussed in the next section.
To simulate fracture and failure, peridynamics uses the notion of bond damage [25,47].Peridynamic bonds break irreversibly when they reach the critical relative deformation  0 , which can be related to material's fracture energy  0 .Note that other types of failure can be considered in PD models, including reforming bond connections [61].For the case of the conical micro-modulus function,  0 is given by [60]: Note that in heterogeneous materials  0 also depends on location.With the breakage of bonds, failure starts to accumulate, and cracks begin to initiate and propagate.The damage index  is used to measure the damage level: which, in the discrete version (see below) is the ratio of the number of broken bonds to that of total bonds connected to point  at time (or load step) .When all bonds connected to point  are broken, (, ) = 1 and point  becomes a free point.

Numerical discretization
Eq. ( 1) can be solved by any method that can solve integro-differential equations, including mesh-free direct discretization [47], the finite element method (FEM) [62,63], or a combination of both in which the FEM is used far from cracks, and the meshfree discretization is used where damage happens [63][64][65].Spectral methods can be alternative approaches to achieve efficient peridynamic computations [66].Here we use the meshfree discretization, which makes it easiest to handle damage and fracture [31,67].
We discretize the domain into cells with nodes in the center of those cells, effectively using the mid-point integration scheme to approximate the integral in Eq. ( 1).Both uniform [47] and non-uniform [68][69][70] grids are possible.Fig. 2 (a) shows a 2D uniform discretization with grid spacing  around a node   .Non-uniform grids conform better for shapes with rounded boundaries.Such grids can be easily obtained based on finite element meshes generated, for example, with ANSYS and a simple APDL code (see the appendix in [55]).As shown in Fig. 2  At time , the spatially discretized form of Eq. ( 1) is written as: where    is the displacement of node ,   is the horizon region of node ,  ∈   includes all the nodes covered by   (fully or partially),   is the micro-modulus of the bond  -,    is the relative stretch of bond  -,   is the area of node  covered by   ,  is the unit vector pointing from node  to node  in the current configuration and    is the body force at node .
Since node  may be only partially covered by the horizon of node , the "partial volume" integration scheme is used to improve the quadrature accuracy [31,71].Note that the partial volume integration scheme was developed for regular grids.We used the same scheme on irregular grids because the irregular grids used here do not depart much from uniform grids (most element sizes are about the same).Moreover, the use of the "conical" micromodulus (see Eq. ( 5)) helps with reducing the quadrature error since the influence of nodes near the edges of the horizon is smaller than that of nodes near the center of the horizon.
For the type of fracture we simulate, inertial effects are likely minor, thus all of the simulations performed in this work are quasi-static.The displacement-controlled loads are applied step-by-step, and at each step we solve the nonlinear (in displacements) system in Eq. ( 8) and use the criteria in Eq. ( 6) to determine which bonds need to break at this step.The values for the micro-moduli c  and the critical relative deformation s 0 for each bond are determined by the bond type (see section 4 below).We solve the equilibrium nonlinear system in Eq. ( 8) via the energy minimization method, using the nonlinear conjugate gradient (NCG) method with secant line search.The detail of the algorithm can be found in [32,35].Instead of the Polak-Ribiere formula [72], we use the hybrid Hu-Storey (HuS) formula [73], for a faster convergence.Note that other methods, such as the adaptive dynamic relaxation (ADR) method [74] and the direct sparse matrix solver [75], can also be used for quasi-static fracture problems in PD.Compared with a direct matrix solver, the NCG solver is faster and uses less memory.The ADR uses a variable artificial viscosity coefficient that sometimes can lead to some unphysical effects in the solution [76].The NCG does face convergence difficulties in problems in which the relation between bond force and bond stretch has a softening part, due to loss of positive definiteness.Here we are using the linear elastic (with brittle failure) force-bond strain model, and the NCG is a good option.
The overall simulation process is shown in Fig. 3.At a given load step, the NCG solver is called to find the equilibrium displacement field.On these displacements, the bond-breaking subroutine is called to check if any bonds exceeded their critical strain.If there are any such bond breaks, the NCG solver is called again (at the same load step) until no more bonds break and the physical system reaches equilibrium.If too many bonds break after one NCG solution, numerical instabilities may result [32].To prevent this, once the number of broken bonds at a step reaches a user-provided parameter  max , we need to go back half a step (re-compute the equilibrium for half of the load increment used in the previous iteration) and repeat the static solution.
In the examples shown in this paper, the crack growth is stable, and as long as a reasonable number of load steps are used, numerical instabilities are avoided.We have conducted tests with 100 and 1,000 load steps and found no significant difference between these splits.All the results shown here are therefore using 100 load steps.
Fig. 3. Flowchart for the simulation process of PD quasi-static fracture modeling.Here  is the load step.

The IH-PD model for concrete
In this paper, the regular homogenization approach is called the "fully-homogenized" peridynamic (FH-PD) model to distinguish it from the "intermediately-homogenized" peridynamic (IH-PD) model discussed below [58].In the FH-PD model, concrete is seen as a locally homogeneous material in terms of its mechanical properties (elasticity, density, and fracture energy).The properties used in the FH-PD model are the macro-scale properties obtained from direct experimental measurements.The multiscale IH-PD model has been originally introduced for functionally graded materials (FGMs) and porous materials, like rock, in [54,58].Here we adapt it for a two-phase composite, in which we define three types of intermingled sets of bonds: one bond-type for each of the phases (aggregates and matrix), and one for interfacial (aggregate-matrix) bonds.The model uses some meso-scale information (volume fraction of the phases) but does not preserve the topology of the microstructure phases.A discussion on conditions under which this approach is still sufficient to capture the fracture behavior accurately and efficiently can be found in [55].
Consider the two-phase (phases A and B) composite material shown on the left side of Fig. 4. At the microscale, an arbitrary PD bond connects points  and  ̂, whose geometrical positions fall in one of the two phases, A and B. We assume that the horizon, and therefore the size of the nodal volumes, are at a scale larger than that of the inclusions so that the composite volume fractions in a nodal volume are statistically representative of the macro-scale value.The probability for a bond to have the properties of phase A, B, or interfacial properties, depends on the volume fraction of the phases over the nodal volumes/areas of the two nodes.
A bond with properties of the A or B phases will be called an A-bond or a B-bond, respectively, while a bond with interfacial properties will be called an AB-bond.Fig. 5 shows an example of bond-type distribution in the IH-PD model around a particular node.
In the IH-PD model, we assume a linear relationship between the chance of the bond type and the local phase volume fractions at the two end nodes, but other choices could also be made.If the volume fractions of phase A are  and  ′ at  and  ̂, respectively, the chance for this bond to be an A-bond, B-bond, or an , respectively.While we will not use the specific distribution of the continuous or discontinuous phases (in order to end up with a computationally efficient model), the volume fraction information is included, and the model is, at the small scale, heterogeneous.For a discussion on when the topology of the phases is important, see [55].Since the explicit microstructure is not used here, and only the volume fraction information is input data, we generate the bond properties as a preprocessing step shown in Fig. 6.Notice that we only select bond types, not node types.The algorithm visits each node in the discretization, then considers each bond in that node's family (if it had not already been assigned its properties) and assigns its properties based on a random number generated from a uniform distribution (see Fig. 6).
For concrete, phase A is aggregate and phase B is cement.We assume that the concrete is homogeneous at the larger scale so that the phase volume-fractions are constants throughout the domain, i.e.  = ′.For example, if  = ′ = 40%, then 16% of all PD bonds end up as aggregate-bonds, 36% as cement-bonds and the remaining 48% as interfacial (aggregate-cement) bonds.The mechanical properties for the interface can be chosen as the arithmetic or the harmonic averages [58] of the two phases A and B. In a recent paper on failure in solder joints [77], an area-weighted harmonic average method was introduced for computing elastic properties of PD interfacial bonds in order to reduce/eliminate oscillations in strains observed at an interface when other options are used.For concrete, the harmonic average is a good option for the elasticity of the interface, according to the test results given in [78]: where   ,   , and   are the mechanical properties corresponding to aggregate, matrix, and interface, respectively.The micromoduli of the three bond-types will be computed to match (see Eq. ( 5)) the   ,   , and   moduli, respectively.
The fracture properties (s 0 values in Eq. ( 6)) for A-and B-bonds are computed based on the fracture energies of the two material phases: aggregate and mortar.In terms of the interfacial fracture property, we take into account that in concrete, the interface between aggregates and matrix is generally weaker than both of them.The fracture energy of the interface in concrete materials is found in experiments to be 4% -34% that of the mortar, with a value estimated to be between 2.5 -25.3 N/m [79].Surface roughness also affects the fracture energy, increasing fracture energy with increasing roughness [79].We choose the fracture energy of the interface to be 25% of mortar [19], and the corresponding property ( 0 ) of AB-bonds will be computed to match this value.(see Eq. ( 6)).
When the random number generator used in instantiating the bond properties in the model employs a different "seed", a different realization of the microstructure is obtained.Different simulations with the same IH-PD realization and same input data give, obviously, the same result.Different results are found using different microstructure realizations of the IH-PD model, even when the rest of the input data is the same.However, for a fixed horizon size, differences between such solutions become smaller the finer the grid (the larger the ratio of horizon to grid spacing, or the "-value" [80]) is.

Radial displacement model for effective expansion of corrosion product
Carbon dioxide and chloride ions from the environment can both depassivate the rebar surface [1,2].Carbon dioxide can neutralize the alkalinity of concrete and make the passive film unstable.Corrosion induced by carbonation usually happens uniformly.Chlorides, however, usually destroy the passive film locally, which results in non-uniform pitting corrosion on the rebar surface.Corrosion can be more nonuniform when one considers the heterogeneities (defects and pores) at the concrete-rebar interface.Nonuniform corrosion leads to the non-uniformly distributed expansion of corrosion products.Various methods have been developed to mimic this expansion, rather than solve for it: radial displacement [81], internal pressure [82], or thermal expansion [83].The purpose for these models is to use the distribution of displacements/pressures created by the expansion of corrosion products onto the concrete as a boundary condition, thus eliminating the need for explicitly modeling the rebar itself.
Here, we choose the radial displacement method and select the von Mises distribution model to approximate the corrosion pattern.We implement this distribution as displacement boundary condition on the inner surface of the rebar hole [19], without actually modeling the rebar.The von Mises model for corrosion pattern is simple to implement and has shown good accuracy compared to available experimental data.The parameters in the formula have physical meanings and are easy to manipulate.
Using the von Mises model, the expansion of corrosion product can be written in the form of a radial displacement as (see Fig. 7): where  0 () is the modified Bessel function of order 0 written as  0 () = ∑ . The meaning of all parameters can be found in Table 1.The values of  and  are fixed to be 4 and 0.0003, respectively, for all the following numerical problems, while the values for other parameters depend on each problem.
A sketch of the von Mises radial displacement model and displacement plots for different k values (level of non-uniformity, k = 0 means uniform corrosion) are given in Fig. 7.The thickness of the expanded corrosion product around the rebar surface (see Eq. ( 10)) will be used in our model as displacement boundary condition on the concrete inner hole surface, where the rebar expansion pushes against.To apply these conditions, considering the surface effect in peridynamic models [84], we use a fictitious layer of nodes outside of the domain (in the hole region) and enforce these displacements to all of the nodes in the fictitious region instead of only to the domain nodes located nearest to the surface of the hole.This is done to reduce the peridynamic surface effect.This layer is shown in Fig. 8, and its thickness equals the horizon size .
With the horizon size approaching zero, the PD boundary condition converges to the classical boundary condition [85].
Table 1 The meaning for parameters used in Eq. (10).
Parameter Meaning There are different ways to implement the displacement boundary conditions defined by Eq. (10).For the elastic problem (no damage is allowed), the radial displacement boundary condition at the rebar hole surface is assigned in one step.For cases involving fracture, we incrementally increase the imposed displacements.One way to reach the values provided by the formula in Eq. ( 10) is to split the total/final radial displacement distribution, at each point around the (initially) circular rebar hole surface, in multiple equal steps (e.g. 100 steps).This is a simple option.However, in reality, certain part of the rebar corrodes earlier than other parts of the rebar, as shown in Fig. A2 in Appendix A. Notice that the modeling of chloride diffusion given in Appendix A is important because it helps us in applying the von Mises boundary conditions on the rebar hole surface in Fig. 8. Therefore, another option is to control the sequence of radial displacement at different locations around the rebar hole surface, as described in Appendix B.
Fig. 8. Imposing displacement boundary conditions at the rebar hole surface using a fictitious node layer.Notice that the rebar is not included in the model.
Other factors such as the porous zone at the concrete-rebar interface [19] and the movement of corrosion products into cracks [13,86] also play a role in the cracking process, affecting the crack initiation time and the speed of crack growth.In the present work, we focus on the fracture patterns and its evolution in a quasi-static setting.

Fracture in concrete due to corrosion of single/multiple rebars 6.1. Concrete structure with a top-sided middle rebar
In this section, we study a particular reinforced concrete structure with one rebar located as shown in Fig. 9.For this case, only the region inside the red dashed contour will be shown in the numerical results.The bottom of the concrete structure is under roller support.The displacement boundary condition at the rebar hole surface was presented in section 5.The material properties, used in all of the following simulations are gathered from several references and given in Table 2. Concrete properties listed in Table 2 are measured directly, not derived from those of concrete components.The volume fraction of aggregates is 40%, which is a common value used in literature [19,56].We first verify the elastic solution by comparing with FEM results.Then we conduct a PD convergence study of the fracture patterns using different horizon sizes.We perform a parametric study to show the influence of aggregate fracture energy on the formation of cracks.We also show the evolution of fracture.After that, we compare the numerical simulation results for deformation of concrete surface with experiments.9)) 15.0 * See section 6.1.2.3 for a parametric study.

Verification for the elastic response
We first verify the PD model for the elastic response (damage is not active) using the concrete structure with a top-sided middle rebar shown in Fig. 9.We compare the FH-PD and IH-PD results with those from a FEM solution, when the displacements induced from the rebar corrosion are given by taking  = 5 in Eq. (10).It should be noticed that the FEM solves the local model while the PD-solutions are for corresponding nonlocal models.The local and nonlocal solutions, in general, are different, but the nonlocal model, in the limit of the horizon going to zero, should converge to the local elasticity solution [25,80].
For the FEM simulation, we generate a conforming map-mesh around the circular rebar with the total node number of 14,048 (see Fig. 10 (a) and (b)).The rebar is not included in the model, but the displacement boundary condition at the concrete inner hole surface is the same as that used in the PD simulation (but only applied at the interface nodes, not over a layer of nodes).Fig. 10 (c) and (d) show the displacements obtained with ANSYS using four-node plane strain elements (plane 182).For the PD simulations with a uniform grid (for horizon size 2 mm and node spacing 0.5 mm; total number of nodes 90,000), the displacements with both FH-PD and IH-PD models are shown in Fig. 11.The contour plots for both PD models (plotted with Tecplot) and the FEM ANSYS solution are close to identical.The same color legend was used in all PD results, and we tried to match with the ones produced by ANSYS (some colors may have slightly different nuances between Tecplot and ANSYS).Note that the horizon size must be smaller than the smallest relevant geometrical feature of the model [29], the rebar size in our case.Otherwise, stress concentrations and cracks initiating from the rebar may not be captured accurately.The horizon size may need also to be correlated to the damage process zone if damage is involved [29].The agreement between the PD and FEM solutions is very good, except for some small differences near the boundary, caused by the peridynamic surface effect [84].Note that a map-mesh cannot be used directly in the PD model due to the large size differences between elements near the rebar and those far from it.With an adaptive approach [68] or with the dual-horizon PD model [89], one could use such a mesh to generate the discretization nodes for the PD model.In the following PD simulations, we either use a uniform grid ( = / with  = 4) or free-mesh conforming grid generated in ANSYS (linear quad elements with edge-length equal to /4).See section 3 for how non-uniform grids generated with ANSYS are transformed into PD grids.

Fracture of concrete with a top-sided middle rebar
Fig. 12 shows the experimental results from [81], for fracture patterns in the case with a top-sided middle rebar.In our simulations we use the same geometry as used in the experiments.However, material properties for the aggregate, mortar and interface used in our simulation are from other sources (see Table 2) because they are not provided in [81].External current was used to speed up the corrosion process in this experiment, which may result in different corrosion profile compared to natural conditions.Also, experimental results for only two samples were given in [81], with significant variability between their fracture patterns.As a result, we can only use these experimental observations for qualitative comparisons with our simulation results.First, for the concrete structure with a top-sided middle rebar (see section 6.1), we compare the fracture pattern obtained by FH-PD and IH-PD models, using uniform and non-uniform grids, respectively.The fracture energy of aggregates is taken as 500 N/m and  = 5 (level of non-uniformity, see Eq. ( 10)).It should be noticed that fracture of aggregates is not considered in previous works [19,21], this being equivalent to selecting an infinite fracture energy for aggregates.In our model, however, we select a large but finite value to also allow failure of aggregate PD bonds because aggregates do, sometimes, fracture when concrete fractures.As seen in Fig. 13, with the FH-PD model we do not get convergence in terms of fracture patterns as we take the horizon size  to zero while keeping the ratio of  and  fixed (the convergence, see [80]).With either uniform or non-uniform grids, the fracture patterns change rather significantly when different horizon sizes are used.Although one of the fracture patterns (obtained using the larger horizon and the corresponding coarser grid) appears similar to one of the experimental observations, (see Fig. 12), the vertical crack found by the FH-PD model for this case initiates at the rebar hole surface and then propagates to the concrete surface.This is opposite to what is reported in the experiments [81,86] as well as previous numerical studies [17,19].In this fully homogenized model, with linear-elastic and brittle fracture, once the cracks start from the rebar hole surface, stress is redistributed uniformly through the material, and the tensile bending deformation at the concrete surface is relieved.Lower tensile stress at the surface means a crack cannot initiate there.The initial cracks continue to grow towards the concrete surface.Meanwhile, one or more cracks may initiate at the interface and propagate.
The reasons for the failure of the FH-PD model to converge come from changes in stress distribution caused by different grids around the rebar hole surface.Small changes in the location of highest tensile stress (strain) induced by use of different discretization grid sizes (especially with uniform grids), leads to variations in the location of crack initiation.These relative differences increase as the cracks grow.The differences in crack patterns are smaller when non-uniform conforming grids are used (see Fig. 13(b)) because the stress state is less dependent on the grid.Nevertheless, any small difference is still amplified in terms of crack paths because the material is homogeneous and brittle, and therefore, sensitive to any small numerical perturbations.10).
In the IH-PD models, the bond structures are different between the tests because different grid densities are used.When we use the same horizon size but create a new realization of the PD bonds structure, the fracture patterns obtained are remarkably similar to one another (see Appendix D).As can be seen from Fig. 14, the fracture patterns show similar features for different horizon sizes too, matching well those from experiments (see Fig. 12).Moreover, with the IH-PD models, the vertical crack happens before the horizontal ones fully split the sample, which is consistent with the experimental observations.To further explore the reasons for the apparent success of the IH-PD model in contrast with the FH-PD model, in the next section, we discuss in detail the evolution of the fracture process in the IH-PD simulations.
Since the fracture patterns obtained with  = 2 mm are not much different from those with  = 1 mm, for the remaining simulations we use the IH-PD model with  = 2 mm, unless otherwise stated.Also, for simplicity, we only use uniform grids for the rest of the simulations.

Fracture evolution with the multiscale PD model
To explain why  -convergence for crack patterns happens in IH-PD but not in FH-PD, it helps to recall that bond properties are randomly distributed (to match the volume fraction of the phases in terms of the elastic response) in the IH-PD case.This small-scale variability leads to a relative insensitivity to variability in the computational grid.The results in Fig. 14 (a) and Fig. 14 (b) show that crack patterns are about the same for the uniform grid and the non-uniform grid (conforming to the round rebar).
A typical evolution of fracture obtained with the IH-PD model is shown in Fig. 15.Micro-damage first accumulates around the rebar hole surface due to failure of weaker bonds (most of them should be ABbonds) around the interface where the displacement loading conditions are applied.Damage starts to localize into horizontal cracks.However, since the material ahead of the crack tip is composed of PD bonds with different failure resistance, the horizontal cracks propagate but may arrest in regions with higher crack growth resistance.When that happens, due to the continued loading, a vertical crack can initiate on the top surface and propagate across the concrete cover towards the rebar.The horizontal cracks may continue to grow and approach the sides of the concrete specimen.As the expansion of corrosion products continues, additional cracks may start from the rebar hole surface and propagate.
The stage when micro-damage around the rebar forms is difficult to detect experimentally.Nevertheless, the initiation and propagation of the vertical crack are both consistent with the post-mortem experimental observations [81,86,90] and other numerical simulations which utilized the meso-scale structure of concrete [19][20][21]81].It can be shown by both theoretical and FE analysis that, before fracture occurs, the maximum circumferential stress is located symmetrically about the vertical axis of the rebar [17].From these locations, horizontal cracks initiate and propagate first, but the propagation is constrained by "chains" of aggregate bonds.Once horizontal cracks have grown sufficiently, the loading is similar to bending of a beam: vertical cracks start from the concrete surface because the top concrete surface is under highest tensile horizontal loading, while the rebar top region is under bi-axial compression [91].Because of this, tensile horizontal stresses at the concrete surface build up and eventually lead to breakage of bonds.Only when the vertical crack reaches the rebar, can the horizontal cracks continue their propagation.In real problems, the microstructure randomness of concrete can lead to significantly different corrosion profiles around the rebar.We use our model to test how different "shapes" of corrosion distribution can affect the fracture patterns in the reinforced concrete.The results presented in Fig. C1 show that the shape of the corrosion product pressure function (see the parameter  in Eq. ( 10)) and the sequence of the imposed displacements (gradual or simultaneous corrosion around the rebar) affect the final fracture pattern but not in a significant way.Therefore, for the remaining simulations for the sample with a top-sided middle rebar (see Fig. 9), we only use the k = 5 value.The corrosion process is assumed to happen simultaneously around the rebar in all remaining simulations.It should also be noticed that, due to the embedded randomness, slightly different fracture patterns are obtained when using different bond-scale realizations of the IH-PD model.A brief study on this is included in Appendix D.

Parametric study in terms of the aggregates' fracture energy
As mentioned before, aggregate fracture is ignored in the available meso-scale models which use explicit representations of the aggregates.However, in physical tests, cracks do sometimes cut through aggregates.In order to understand the effect of aggregate-type bonds in the IH-PD model has on the overall failure behavior of the concrete-rebar structure, we perform a parametric study.Fracture patterns obtained by using different fracture for computing the critical strain aggregate bonds in the IH-PD model are shown in Fig. 16.We vary the fracture energy for such bonds from the small value equal to that of the mortar (not entirely realistic) to some arbitrary value (eight times larger than that of the mortar).When the fracture energy for aggregate bonds is too small, the vertical crack does not match the behavior seen in reality (Fig. 16 (a)).When it is too large, significant damage spreads sideways from the major crack paths (Fig. 16(d)), also not observed experimentally.The latter result is caused by the presence of the network of aggregatetype bonds in the IH-PD model, which is in contrast with the actual concrete microstructure (aggregates are inclusions, not long chains spanning the sample).The different topology at the micro-scale in the IH-PD model compared with the actual microstructure does not affect the material behavior in the elastic regime, but it can affect it once damage is considered [55].With an intermediate value of G = 500 N/m, the crack patterns match well the trends seen in experiments.This value will be used for the remaining simulations.10)) with  = 2 mm and  = 4.

Surface deformation
As the rebar corrodes and the expansion of the corrosion product leads to internal pressure build-up, the concrete surface is experiencing horizontal tensile stresses, similar to the case of beam bending.The surface crack initiates once these tensile stresses/strains reach a critical value.A comparison with experimental results can be made for the vertical deformation of the concrete surface [81].We also compare the width of the surface crack.For these simulations, we choose a horizon size equal to 1 mm because the smallest width of the surface crack measured in experiment is around 2 mm.
To estimate the surface crack width, we assume it to equal the smallest relative displacement between any two top-surface nodes on opposite sides of the crack with damage index bigger than 0.4.In Fig. 17 we plot the deformation of the concrete surface (vertical displacement) and give the values of the crack opening width in the legend.The experimental curves were measured from the sample given in Fig. 12 (a) at two different times.While the crack opening width is only qualitatively matched, the surface deformation curves are very similar to those measured in experiments.Some difference between experimental data and our numerical results is expected because our input data (the pressure profile) is only a rough approximation of the actual conditions.It should be noticed that in Fig. 17 (b) the maximum deformation in experimental curve is not located at the surface crack of the experimental fracture pattern (see Fig. 12 (a)), which is unexpected.No reason was given for this inconsistency, but it is possible that the deformation was not measured exactly at the location of the cross-section shown in Fig. 12 (a), but at a different location.

Concrete structure with a corner rebar
In this section, we study two cases in which the rebar is located at the corner of the concrete structure as shown in Fig. 18.The material properties are the same as those in the case with a top-sided middle rebar.The geometry data is given in Table 3.The horizon size is 1 mm, and node spacing is 0.25 mm.A smaller horizon size is used because the rebar size here is smaller than the concrete structure in the previous example (see section 6.1.1 for a discussion on selecting the horizon size).12 (a)) as well as the RBSM solutions in [81]).Data in the legends refers to the crack opening width at the top of the sample.Table 3 Geometry data for two samples of the reinforced concrete structure shown in Fig. 18 [92].
Geometry data Sample 1 (mm) Sample 2 (mm) The experimental results (and zoom-in images) for two different samples (see Table 3) are shown in Fig. 19.We note that the experimental setup used external current to accelerate corrosion [92].The chloride appears to penetrate through the left side of the concrete cover according to Fig. 19 (a) and (b).The zoom -in images show that the corrosion profiles vary significantly between the two samples, which could be, besides the slightly different geometrical setups, the main reason why the observed fracture patterns are also significantly different.
For our PD model, we assume the corrosion profile for the top-sided middle rebar is rotated by an angle , as shown in Fig. 20, to account for the chloride penetration from the left side.For sample 1 with  = π/2,  = 3π/8 and  = π/4, we obtain fracture patterns in Fig. 21  Because only one sample was provided in the experimental work [92] for each geometry, it is difficult to draw a definitive conclusion in these cases.Some differences between our simulation results shown in Fig. 21 and the experimental observations shown in Fig. 19 may be attributed to several factors, such as: the partially-homogenized PD bond structure is not an exact representation of the actual concrete microstructure; the boundary condition used to mimic the expansion of the corrosion products is an assumed approximate distribution instead of real values; some features could be reflections of 3D effects (possibly more significant for the corner rebar case than the previous symmetric rebar case), that our 2D model cannot be expected to replicate.
We noticed that references [13,14] (using a specialized constitutive model for concrete) obtained fracture patterns very similar to experimental observations using a 3D model.They obtained the radial expansion of corrosion product from the corrosion of rebar and applied it to the contact element between the rebar and the concrete.This coupling may produce more realistic boundary conditions for the corrosion expansion.However, the anodic regions around the rebar were determined beforehand in all simulations, by selecting just one pair of anode and cathode in each activated cross-section of the rebar.This may not be realistic since microcell corrosion in reality usually consists of many pairs of mixed anodic and cathodic areas [2], which would result in a relatively uniformly accumulated corrosion product along the depassivated rebar surface.The actual corrosion process is a combination of micro-and macro-cell corrosion [2].Macrocell corrosion involves pitting corrosion, leading to non-uniform corrosion patterns.Our way of applying radial displacement boundary condition (see Section 5) is a grossly simplified, but perhaps more reasonable approximation of real conditions because it can be seen as a weighted combination of pitting corrosion and uniform corrosion.Also, according to the pressure profile in [13,14], the corrosion product accumulates in zones that do not seem to agree with the experimental observations from [92] (shown in Fig. 19).Only one mesh was used for each example shown in [13,14].Since the crack band theory used to simulate the fracture process in [13,14] has mesh orientation dependence ( [15]), it would be interesting to see if any changes take place in the reported results if one uses a differently oriented mesh and whether possible changes in fracture patterns would match the fracture patterns variability seen experimentally.

Corrosion induced fracture from multiple rebars
A reinforced concrete structure with three rebars, shown in Fig. 22, is analyzed next.Only the region inside the red dash contour is shown in the following numerical results.

Conclusions
We introduced a 3-phase stochastic multiscale intermediately-homogenized peridynamic (IH-PD) model to study fracture in reinforced concrete due to non-uniform rebar corrosion.Different from traditional mesoscale heterogeneous models, our model only uses the volume fraction of different phases in the heterogeneous material.A simple constitutive model, linear-elastic with brittle failure, is used.We show and explain the reasons for which a fully homogenized peridynamic model leads to fracture patterns and failure evolution different from what is observed experimentally, while the new IH-PD model results match experimental observations very well.The model does not consider details such as aggregate sizes and shapes, making the analysis much simpler/cheaper than that of models with explicit representation of the microstructure geometry, matching the computational cost of a fully homogenous model.
The corrosion product expansion around the rebar was approximated here by the "von Mises model" from the literature.This simplified the analysis.For a more accurate representation of pressures induced by corrosion product formation around the rebar, the current model could be coupled with an explicit corrosion model.
We tested the IH-PD model on concrete structures with a single or multiple rebars.The numerical results match well experimentally observed crack patterns as well as the sequence/evolution of their growth.We performed computations using several different stochastic realizations of the peridynamic bond structures and found that fracture patterns remain within the variability of features observed experimentally.
The IH-PD model succeeds in balancing the accuracy of fracture prediction in concrete from (expensive) models that use an explicit representation of aggregates, with the efficiency of homogeneous models.It can be used for larger scale modeling without losing the influence the microscale has on the failure behavior of the material.

Appendix B. Applying radial displacement at the rebar hole surface
As mentioned in Section 5, the simplest way to apply the radial displacement boundary condition is to incrementally increase the imposed displacements at each point around the rebar hole surface in multiple equal steps, assuming that the corrosion around the interface happens simultaneously.However, considering the fact that certain part of the rebar corrodes earlier than other parts, as shown in Fig. A2, another option is to control the sequence of radial displacement at different locations around the rebar hole surface to mimic the evolution of the corrosion front (or corrosion products) around the rebar.To explain the idea, we only consider the right half of the rebar hole surface because of symmetry.As shown in Fig. B1, the right half of the interface is evenly divided into  sectors.The incremental steps (total number of 100) at which these sectors are applied the radial displacements are shown in Fig. B2.For boundary nodes in sector , the radial displacement is applied evenly from t −1 to t  .We consider both these options in Fig. C1.Because the total radial displacement is very small (the concrete is very stiff), we can assume that the shape of the rebar hole surface is close to circular at all times, and therefore the direction of the radial displacement is always normal to the rebar hole surface.

Fig. 1 .
Fig. 1.Nonlocal interaction between point  and an arbitrary point  ̂ located in its horizon   .
(b), the PD nodes are centroids of elements, and the PD nodal area of each node is the element area.To keep the quadrature error low, it is important to have meshes with relatively uniform element sizes.Comparisons between uniform and non-uniform grids for corrosion-induced fracture are shown in section 6.1.2.1.(a) Uniform grid.(b) Non-uniform grid.Fig. 2. Possible discretization types for a peridynamic model.The circular region is the horizon region of node   .

Fig. 4 .
Fig. 4. Determining the properties of a PD bond in the IH-PD model, based on the local volume fractions of phase A at points  (), and  ′ ( ′ ).

Fig. 5 .
Fig. 5.A possible distribution of bond properties at a node in the IH-PD model.Only bonds connecting the nearest eight neighbors of the central node are shown.

Fig. 6 .
Fig. 6.The pre-processing step for instantiating bond properties in the IH-PD model.

Fig. 9 .
Fig. 9. Geometry of the concrete structure with a top-sided middle rebar (mm).
(a) Map-mesh for the concrete with one rebar.(b)Zoom-in of the map-mesh for the concrete around the rebar (arrows schematically represent the imposed displacement boundary condition on the rebar hole surface in Fig.8).(c)Horizontal displacement.(d) Vertical displacement.

Fig. 10 .
Fig. 10.FEM mesh and displacements computed using ANSYS (only the highlighted region from Fig. 9 is shown).

Fig. 11 .
Fig.11.Displacements obtained with the PD models (only the highlighted region from Fig.9 isshown).For the PD simulations with a uniform grid (for horizon size 2 mm and node spacing 0.5 mm; total number of nodes 90,000), the displacements with both FH-PD and IH-PD models are shown in Fig.11.The contour plots for both PD models (plotted with Tecplot) and the FEM ANSYS solution are close to identical.The same color legend was used in all PD results, and we tried to match with the ones produced by ANSYS (some colors may have slightly different nuances between Tecplot and ANSYS).Note that the horizon size must be smaller than the smallest relevant geometrical feature of the model[29], the rebar size in our case.Otherwise, stress concentrations and cracks initiating from the rebar may not be captured accurately.The horizon size may need also to be correlated to the damage process zone if damage is involved[29].The agreement between the PD and FEM solutions is very good, except for some small differences near the boundary, caused by the peridynamic surface effect[84].Note that a map-mesh cannot be used directly in the PD model due to the large size differences between elements near the rebar and those far from it.With an adaptive approach[68] or with the dual-horizon PD model[89], one could use such a mesh to generate the discretization nodes for the PD model.In the following PD simulations, we either use a uniform grid ( = / with  = 4) or free-mesh conforming grid generated in ANSYS (linear quad elements with edge-length equal to /4).See section 3 for how non-uniform grids generated with ANSYS are transformed into PD grids.

Fig. 15 .
Typical fracture evolution using the IH-PD model.
(a) G = 100 N/m (b) G = 300 N/m (c) G = 500 N/m (d) G = 800 N/m Fig. 16.Fracture patterns produced by different choices for the aggregates' fracture energy (case with k = 5 in Eq. ( (a) step = 50 (b) step = 100 Fig. 17.Comparison of surface deformation at two different corrosion stages between IH-PD results, experimental observations (see the sample in Fig.

Fig. 18 .
Fig. 18.Concrete structure with a corner rebar: geometry and mechanical boundary conditions.The bottom and right sides are symmetry lines.
(a), (b) and (c), respectively, with two different k values.The crack patterns for sample 1 in the experiment, shown in Fig. 19 (a), have a cornertype symmetry, while the computed results less so.For sample 2, Fig. 21 (d) shows fracture patterns obtained with  = π/4.The fracture pattern for  = 10 shows diagonal symmetry, while the experimental fracture pattern, shown in Fig. 19 (b) does not.
(a) Sample 1 and the zoom-in at the rebar (b) Sample 2 and the zoom-in at the rebar Fig.19.Experimental results for cracking in reinforced concrete due to corrosion of a corner rebar (from[92]).

Fig. A1 .
Fig. A1.A zoom-in sketch of the concrete structure with a top-sided middle rebar under diffusion of chloride.A rebar node is "corroded" once half of its family nodes reach the critical chloride concentration.

Fig. A3 .
Fig. A3.Points in the rebars with a concentration larger than a threshold value (shown in red) obtained from chloride diffusion (from the top surface only) for the concrete structure with three rebars.

Fig. B2 .
Fig. B2.Applying radial displacements increments for different sectors around the rebar considering the evolution of the corrosion process.Radial displacement for each sector is applied incrementally during different time periods.

Fig. D1 .Fig. D2 .
Fig. D1.Fracture patterns obtained from different realizations of the bond-structure in the IH-PD model for the top-sided middle rebar geometry.