Concrete fracture caused by corros ion 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 fail s to captur e 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 product s is enforced here as preset, incremental radial displacement s. The computed fracture patterns and the order in which various cracks develop match what is seen in experiment s. The model ’s robust ness is tested under different stochastic realizations and discretization grid types . Keywords Peridynamics ; reinforced c oncrete ; fracture ; brittle fracture; rebar c orrosion 1. 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 rebar s [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 manufac turing or early loading stage s [3, 4]. Small amounts of a ggressive elements then reach the rebar s 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 product s (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 , accelerat ing the corrosion process [5]. The corrosion of the steel rebar decrease s its effective cross -sectional area and break s the bond between the rebar and concrete , causing performance deg radation of the structure [6]. Moreover, changes in the layer * Corresponding author. E-mail address: fbobaru2@unl.edu (F. Bobaru). 2 affected by corrosion trigger significant reduction s 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]). Non- uniformities 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-12]. Experiment al investigations of corrosion -induced fracture in concr ete can offer some insights into the process but they are expensive and time-consuming . Usually , these are performed using external current s to accelerate th e 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 th is 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 specia lized 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 use d in t hese 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 stud ies were not presented for comparing variability of crack patterns with those from experiments . A 2D mecha no-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, t he active zone on the rebar evolve s automatically. However, for the fracture model, using the crack band theory , no details (material properties, boundary condition s, 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 experimental ly, but fracture patterns obtained for the 3 - rebar case do not capture the experimental observation very well. One of the limitations of the above -mentioned work s 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 c reate 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 mode ls have shown their potential s for such problems [18]. However, issues such as the selection of the geometric shape of the aggregates can significantly affect the fracture behavior [19-21]. It is possible to use a meso -scale structure acquired from X-ray tomography images [22], but it is computational ly costly to extend me so-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 involve s some of concrete’s microscale features , including their random ness. 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 . 3 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 sensitivit y in fracture problems , while nonlocal ones can prevent it. Another r eason 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 t heory which has received considerable a ttention since its introduction almost two decades ago [25, 26]. PD reformulate s the classical continuum mechanics by eliminating spatial derivatives to model mechanical [27-37], diffusion/corrosion [38-43], or mech ano-chemical [44-46] etc., behavior s in materials involving damage. Using spatial integration rather than differentiation lead s to a mathematically consistent formulation that works naturally for problems in which discontinuit ies in the domain (such as cracks) appear . In a PD mode l, cracks /damage can initiate and propagate auto nomously [47]. One of the first application s of the PD theory to concrete structure s 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 formulat ion of pressure -dependent PD plasticity model was shown to work well for compression, impact , and spallation of concrete structure s [51]. The trilinear peridynamic model introduced in [52] has shown good results in terms of the load -CMOD curves for three - point -bending tests. This model, h owever, being a homogeniz ed 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 -Hom ogenized 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 resul ts 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 . 4 This paper is organized as follows : in Section 2 we give a brief review of t he 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 discuss es 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. 2. 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 model s 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 analys es 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 condition s (Poisson’s ratio equals 1/4 ) everywhere in the paper . The equation s of motion for the BBPD can be written as [25]; 𝜌(𝒙)𝒖̈(𝒙,𝑡)=∫𝒇 𝐻𝒙(𝒖(𝒙̂,𝑡)−𝒖(𝒙,𝑡),𝒙̂−𝒙)𝑑𝑉𝒙̂+𝒃(𝒙,𝑡) (1) 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 force s exist between 𝒙 and 𝒙̂, an arbitrary point locate d inside 𝐻𝒙 (see Fig. 1 ). 𝑉𝒙̂ is the volume (area in 2D, length in 1D) occupied by the materia l 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 𝛿. Fig. 1. Nonlocal interaction between point 𝒙 and an arbitrary point 𝒙̂ located in its horizon 𝐻𝒙. The pairwise force for a proto type microelastic brittle material [47] is defined as: 𝒇(𝝃,𝑡)=𝑔(𝑠(𝝃,𝑡))ℎ(𝝃,𝑡)𝝃+𝜼 ‖𝝃+𝜼‖ (2) 5 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 : 𝑔(𝑠)={ 𝑐(𝜉)𝑠 𝜉≤𝛿 0 𝜉>𝛿 (3) ℎ(𝜉,𝑡)={ 1 if 𝑠(𝜉,𝑡′)<𝑠0 for all 0≤𝑡′≤𝑡 0 otherwise (4) where 𝑐(𝜉) is the micro -modulu s function or the el astic 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 condition s and the “conical ” micromodulus function [60]: 𝑐(𝜉)=𝑐1(1−𝜉 𝛿)=192 𝐸 5𝜋𝛿3(1−𝜉 𝛿) (5) 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 irreversibl y 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]: 𝑠0=√25𝜋𝐺0 48𝐸𝛿 (6) 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 : 𝑑(𝒙,𝑡)=1−∫ℎ(𝜉,𝜂)d𝑉𝒙̂ 𝐻𝒙 ∫d𝑉𝒙̂ 𝐻𝒙 (7) 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 . 3. 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 combina tion of both in which the FEM is used far from cracks , and the meshfree discretization is used where damage happens [63-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-70] grids are possible. Fig. 2 (a) shows a 2D uniform discretization with grid spacing 𝛥𝑥 around a node 𝐱𝒊. Non- 6 uniform grid s conform better for shapes with rounded boundaries . Such grids can be easily obtained based on finite element mesh es generated , for example, with ANSYS and a simple APDL code (see the appendix in [55]). As shown in Fig. 2 (b), the PD nodes are centroids of element s, 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. Comparison s 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 f or a peridynamic model . The circular region is the horizon region of node 𝐱𝑖. At time 𝑡, the spatially discretized form of Eq. (1) is written as: 𝜌𝒖̈𝑖𝑡=∑ 𝑐𝑖𝑗𝑠𝑖𝑗𝑡𝑉𝑖𝑗𝒏 𝑗,𝑗∈𝐻𝑖+𝒃𝑖𝑡 (8) where 𝒖𝑖𝑡 is the displacement of node 𝑖, 𝐻𝑖 is the horizon r egion 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 integra tion 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 (se e 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 si mulations 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 s0 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 t he 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 7 softening part, due to loss of positive definiteness. Here we are using the linear elastic (with brittle failure) force -bond strain mod el, and the NCG is a good option. The overall simulation process is shown in Fig. 3 . At a given load step, t he NCG solver is called to find the equilibrium displacement field . On these displacements, t he bond -breaking subroutine is called to check if any bond s exceeded their critical strain. If there are any such bond break s, the NCG solver is called again (at the same load step) until no more bonds br eak and the physical system reaches equilibrium . If too many bonds break a fter one NCG solution , numerical instabilities may result [32]. To p revent 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 step s 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. 4. 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 homogen eous 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 (FGM s) and porous mater ials, like rock, in [54, 58]. Here we adapt it for a two-phase composite , in which we define three types of inter mingled sets of bonds: one bond -type for each of the phases (aggregates and matrix), and one for interfacial 8 (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 position s fall in one of the two phases, A and B. We assum e 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 interfac ial properties, d epend s 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 interfac ial properties will be called an AB -bond. Fig. 5 shows an exam ple 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 fraction s at the two end nodes , but other choices could also be made . If the volume fraction s of phas e A are 𝑅 and 𝑅′ at 𝒙 and 𝒙̂, respectively , the chance for this bond t o be an A-bond , B-bond, or an AB-bond is 𝑅𝑅′, (1−𝑅)(1−𝑅′), or 1−(1−𝑅)(1−𝑅′)−𝑅𝑅′, 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 t he 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 interfac ial (aggregate -cement) bonds. Fig. 4. Determining the properties of a PD bond in the IH-PD model , based on the local volume fraction s of phase A at point s 𝒙 (𝑅), and 𝒙′(𝑅′). 9 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. The pre-processing step for instantiating bond properties in the IH-PD model . The mechanical properties for the interface can be chosen as the arithmetic or the harmonic average s [58] of the two phases A and B . In a recent paper on failure in solder joint s [77], an area -weighted harmonic average method was introduced for computing elastic properties of PD interfac ial 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]: 𝑀𝐴𝐵=2𝑀𝐴𝑀𝐵 𝑀𝐴+𝑀𝐵 (9) 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 (s0 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 propert y, we take into account that in con crete, 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 10 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 c orresponding 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 simulation s 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, difference s between such solutions become smaller the finer the grid (the larger the ratio of horizon to grid spacing , or the “𝑚-value ” [80]) is. 5. 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 p assive film locally , which results in non-uniform pitting corrosion on the rebar surface . Corrosion can be more non- uniform when one considers the heterogeneities (defects and pores ) at the concrete -rebar interface. Non- uniform corrosion leads to the non-uniform ly distributed expansion of corrosion products . Various methods have been developed to mimic this expansion, rather t han 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 model ing the rebar . The von Mises model for corrosion pattern is simple to implement and has shown good accuracy compared to available experiment al data. The parameters in the formula have physical meaning s 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 ): 𝑇𝑑(𝜃)=(𝑎−1)𝜋𝑅𝜂𝑒𝑘cos(𝜃−𝜇) 2𝜋𝐼0(𝑘) (10) where 𝐼0(𝑘) is the modified Bessel function of order 0 written as 𝐼0(𝑘)=∑(1 4𝑘2)𝑖 (𝑖!)2∞ 𝑖=0 . The meaning of all parameters can be found in Table 1 . The v alues 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 expan sion pus hes 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 11 region instead of only to the domain nodes located nearest to the surface of the hole . This is done to red uce 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 p arameters used in Eq. (10). Parameter Meaning 𝑇𝑑 Thickness of the expanded layer 𝑎 Ratio of rust expansion to corroded rebar 𝑅 Radius of the rebar 𝜂 Corrosion degree 𝜇 Location of maximum thickness 𝑘 Non-uniform coefficient (a) Schematic d istribution of radial displacement s caused by corrosion product expansion (redrawn from [19]). (b) Thickness of the expanded layer aroun d the rebar , for different parameter k values. Fig. 7. The von Mises radial displacement model . There are different ways to implement the displacement boundary condition s 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 inc rease 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, i n 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 condition s 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 . 12 Fig. 8. Imposing displacement boundary condition s at the rebar hole surface using a fictitious node layer . Notice that the rebar is not included in the model. Other factors such as t he 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 pattern s and its evolution in a quasi -static setting. 6. Fracture in concrete due to corrosion of single/multiple rebar s 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 th is case, only the region inside the red d ashed 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