Nonlocal Peridynamic Modeling and Simulation on Crack Propagation in Concrete Structures

An extended peridynamic approach for crack propagation analysis in concrete structures was proposed. In the peridynamic constitutive model, concrete material was described as a series of interacting particles, and the short-range repulsive force and anisotropic behavior of concrete were taken into account in the expression of the interactive bonding force, which was given in terms of classical elastic constants and peridynamic horizon. The damage of material was defined locally at the level of pairwise bond, and the critical stretch of material bond was described as a function of fracture strength in the classical concrete failure theory. The efficiency and accuracy of the proposed model and algorithms were validated by simulating the propagation of mode I and I-II mixed mode cracks in concrete slabs. Furthermore, crack propagation in a double-edge notched concrete beam subjected to four-point load was simulated, in which the experimental observations are captured naturally as a consequence of the solution.


Introduction
Even though the widely used concrete materials have been studied for so many years, their failure mechanism has not been fully understood yet.In particular, the damage accumulation and progressive cracking in concrete materials cannot be predicted accurately with satisfactory effectiveness due to the incapability of current available simulation tools, such as finite element method (FEM) and finite differential method (FDM).The essential reason lies in their mathematical framework, which is based on the classical continuum mechanics theory.In the classical continuum mechanics, it is assumed that the material remains continuous as it deforms, and the spatial partial derivatives are employed to define stress and strain in the body.However, the derivatives will become undefined or ill-suited when discontinuities, such as interfaces or cracks, are present.Various remedial methods (e.g., extended finite element method [1,2], virtual crack closure technique [3], cohesive elements [4][5][6], etc.) have been applied to deal with the discontinuous situation.For example, Cendón et al. [7], Unger et al. [8], and Ooi et al. [9] have modeled the crack growth in concrete using cohesive crack model, extended finite element method, and polygon scaled boundary finite elements, respectively.Nevertheless, numerical methods incorporated with these special techniques still more or less suffer from the unsatisfactory simulation accuracy and low efficiency when handling complicated fracture problems.This motivates a reformulation of the basic equations of classical continuum mechanics and results in an alternative of solid mechanics theory referred to as peridynamics based on a nonlocal mathematical framework, in which the identical spatial integral equations rather than derivatives are applied both on and off of a discontinuity and no special remedial techniques are required whenever a crack or interface appears.
The peridynamic theory was first introduced by Silling [10] as an alternative formulation of solid mechanics for handling discontinuities and long-range forces, which presents a unique capability of analyzing damage and fracture process of materials and structures by predicting the displacements, crack nucleation, and propagation with arbitrary paths directly without any extra numerical techniques or criteria.Since then, peridynamics has been successfully utilized to analyze various problems including progressive failure [11,12], phase transformation [13], and stationary and dynamic fracture [14][15][16] at both macro-and microscales, and it shows great potential analogous but advantageous to both classical meshfree and molecular dynamic methods.
Although peridynamic method has been employed successfully to simulate various failure and fracture problems at different scales, the application of peridynamics to concrete materials and structures is comparatively fewer.Gerstle et al. [17][18][19] proposed a simple bond-based peridynamic model for concrete with a microelastic stiffness  and a cut-off stretch , and this model was updated by Gerstle et al. [20][21][22] for several times to describe the mechanical behavior of concrete in both tension and compression regime more accurately.The latest updated model has been used to predict the failure and cracking of plain and reinforced concrete structures, and the simulation results agree well with experimental observation.The structural stability of a concrete column was investigated by Kilic and Madenci [23] and Kilic [24] by describing concrete material as a microelastic brittle peridynamic material.Huang et al. [25] and Shen et al. [26] simulated the damage accumulation, crack propagation, and failure process of a concrete plate subjected to in-plane uniaxial tension, compression, and out-of-plane impact load.
The present paper aims to investigate the capability of peridynamic model to predict the paths accurately along which the crack will grow and to simulate the fracture process of concrete structures at the bond level.We present here a very simple bond-based peridynamic model for concrete, considering the short-range repulsive force and anisotropic behavior of concrete materials.Propagation of mode I and modes I-II mixed cracks in concrete slabs is simulated by using the proposed model, and the paths along which the cracks grow match the theoretical prediction reasonably closely.To further demonstrate the capabilities of the proposed model, crack propagation in a double-edge notched concrete beam subjected to four-point load was simulated, in which the experimental observations are captured naturally as a consequence of the solution.

Peridynamic Model and Numerical Method
Like a coarse-grain molecular dynamics, in peridynamics, the material concerned in spatial domain  is described as numbers of interactive particles, and a pairwise force function f is defined to describe the interactions of particles within some surrounding finite distance  at any time  (as shown in Figure 1): where u denotes the displacement of a material point.The peridynamic equation of motion at point x and time  according to Newton's law is given as where  is the mass density, b is the prescribed loading force density field which represents the external force per unit reference volume, and  is a neighborhood of material point x within the horizon size of : From the peridynamic equation of motion, it can be found that the material points interact with each other through the pairwise force function f, which contains all the constitutive information of the material but eliminates any spatial derivatives or assumption on the continuity of displacement, analogous to the interatomic potential function of a molecular dynamics model.For a general homogeneous material without memory, the function f can be reduced to where  = x  − x and  = u  − u denote the relative position and relative displacement in the reference configuration, respectively.Hence, the vector ( + ) = (x  + u  ) − (x + u) gives the current relative position between points x and 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000 0.00e + 000  x  in the deformed configuration, and it keeps in the same direction as the interactive force function f; that is, Thus, a general form for the peridynamic constitutive force function f can be written as where (, ) is a scalar-valued function.Considering concrete as a prototype microelastic brittle material in peridynamics, the scalar-valued function  is given by Silling [10]  (, ) = , where  With  referring to the peridynamic microelastic constant,  is the elastic modulus, ] is the Poisson ratio, and  is the bond stretch, which is defined as The break of the bond is described by the factor  as follows: where  0 is defined as the critical stretch or failure of the bond between pairs of material points.Take the different behavior and strength of concrete when subjected to tension and compression into account; here we define   and   as tensile and compressive critical stretch of concrete, respectively, where   and   represent the uniaxial tensile and uniaxial compressive strength coming from experimental results and  is the elastic modulus.Accounting for the number of broken bonds, the damage of material can be defined locally as To avoid the overlap of concrete particles with each other or invasion of one particle into the boundary of other particles, we here introduced a repulsive force or short-range to the force function of concrete, as Parks et al. defined earlier in another work [27], Mathematical Problems in Engineering  where   is the short-range interaction distance between particles and   is a multiple of constant .They are given as follows: where   was defined as the particle radius, which was chosen to be half the lattice constant Δ in the present work.
To solve the peridynamic equations of motion, the region where a peridynamic material occupies is often discretized into particles (material points), and a numerical approximation is usually implemented by reducing the spatial integration into a finite sum.In this study, the concrete peridynamic model is discretized into particles forming a cubic lattice with lattice constant Δ, and the integral in ( 2) is replaced by where  denotes the number of time step, u   = u(x  ,   ), and   = |Δx| 3 stands for the involved volume of particle x  in a cubic lattice.An explicit central difference formula is used for acceleration in the present work The boundary conditions are applied as follows:

Mode I and Mixed Mode Cracks Propagation in Concrete
Slabs.A peridynamic model for a concrete slab with the size of 500 mm * 400 mm is established as shown in Figure 2.
To model a typical concrete, here we select Young's modulus  = 25 GPa, density  = 2400 kg/m 3 , uniaxial tensile strength   = 3.5 Mpa, and uniaxial compressive strength   = 28 MPa.The concrete slab is discretized into particles forming a simple cubic lattice with lattice constant Δ = 1 mm and with material horizon  = 4mm (four times the lattice constant).The time step is set to be 1 × 10 −7 s, based on the combined optimization test of both integration stability and computation efficiency.
A crack with the size of 40 mm * 4 mm is preset by deleting the particles at the center of the slab.To model the propagation process of typical mode I crack and modes I-II mixed crack, the angle between the orientation of the crack and the direction of the uniform distributed load on the boundary, that is,  in Figure 2, is set to be 90 ∘ and 45 ∘ separately in two examples.The loading rate is 2200 kN/s, and the force is applied to each particle on left and right side gradually within each time step.
Figures 3 and 4 show the propagation process of a central mode I crack and the contours of horizontal displacements in concrete slab.It appears that the crack will propagate naturally as an outgrowth of the peridynamic equations of motion.When  = 0.92 ms, there is no damage in the slab.It means no bond is broken and the displacement of material points in the slab keeps continuously distributed, as Figure 4(a) shows.At the time of  = 0.94 ms, several bonds at the crack tips are broken, which leads to damage initiation locally.The damage accumulates as more and more bonds are broken near the crack tip, which results in the crack propagation, as shown in Figures 3(c) and 3(d).The crack will rapidly propagate to the boundary of the slab, and the slab will be broken thoroughly, as Figure 3(e) shows.From Figures 4(d) and 4(e) we can see that the displacement field in the slab obviously becomes discontinuously near the crack, which will be complicated for numerical methods based on classical continuum mechanics theory to handle.
Figures 5 and 6 show the crack propagation process and the contours of horizontal displacement field in the concrete slab with a preexisted I-II mixed mode crack.Similarly, it can be seen that, at the time of  = 1.12 ms, no damage still occurs; that is, no bond is broken in the slab.Pairwise bonds at the crack tip begin to break when  = 1.14 ms.After that time, the crack will propagate rapidly along some definite orientation and grow to the boundary of the slab, as Figures 5(c)-5(e) show.Also, we can see from Figure 6 that the displacement field in the slab keeps continuously distributed before the crack grows, and it gets visibly discontinuously distributed near the crack in Figures 6(d)-6(e).
For the linear elastic fracture problems concerned in this work as shown in Figure 2, the theory on fracture mechanics of maximum circumferential stress can be used to deduce the cracking angle.The circumferential stress at point near the crack tip in Figure 2 can be written as where (, ) denotes the position of points in local polar coordinate. is shown in Figure 2.
The maximum circumferential stress should meet the condition of   / = 0 and  2   / 2 < 0. The first equation can be extended as The solution for cos(/2) = 0 is  = ±, which does not agree with experimental results; thus the condition which the cracking angle will meet can be written as follows: For the problems described in Figure 2, the stress intensity factors are Substituting (21) to (20), the condition becomes where  0 is the cracking angle and  is the angle between the orientation of the preexisted crack and the remote loading direction.
Mathematical Problems in Engineering 0.000e + 000 0.000e + 000 0.000e + 000 0.000e + 000 0.000e + 000 0.000e + 000 0.000e + 000 0.000e + 000 Thus, by solving (22), the cracking angle  0 , that is, the paths along which the preexisted crack will propagate initially, can be obtained.For mode I crack,  = 90 ∘ , the solution is  0 = 0 ∘ .For I-II mixed mode crack shown in Figure 2 with  = 45 ∘ , the analytical cracking angle is  0 = 53.13∘ .From Figures 3 and 5, we can see that, by solving the peridynamic equations of motion with the present model, the crack will propagate spontaneously with the cracking angle of 0 ∘ and 53.1 ∘ separately, which acts in significant accordance with theoretical prediction.

Crack Propagation in a Double-Edge Notched Concrete
Beam.Crack propagation in a double-edge notched concrete beam subjected to four-point load is a classical failure problem and has been investigated by various researchers using different experimental and numerical approaches.For example, Bocca et al. [29] and Carpinteri et al. [30] have made four-point shear tests with different-sized concrete specimens and modeled the fracture with a cohesive crack model as well.What is more is that Saleh and Aliabadi [28] have modeled the same crack growth problem using boundary element method.To further demonstrate the capabilities of the proposed model and compare with other approaches, in this section, we simulate the crack propagation in a doubleedge notched concrete beam subjected to four-point load.
The beam is discretized into 40381 material points with a constant grid length size of 2 mm.The external load is applied on the beam step by step with a load increment of 10 kN, and in each loading step the effect of the external load can be transmitted from the loading point to the whole model and a new equilibrium state can be obtained.When damage occurs around the crack tip, the load increment will be reduced to 1 kN.Finally when the external load increases to some critical value, the system cannot reach a new equilibrium again, and the crack will become instable and propagate rapidly.
Figure 8 shows the evolution of damage maps (crack paths) in the beam as external load increases.For more clear display of the damage contours around the crack tips, the deformation of the beam has been enlarged for 1000 times.It is shown that the beam may hold an external load of 70 kN without any damage appearance (see Figure 8(a)).While the external load increases to 80 kN, damage will occur at the crack tips with broken material bonds.However, the structure can still hold the external load without crack propagation, even if the load increases to 90 kN.However, when the external load increases to 99 kN, the beam system cannot reach a new equilibrium state again, the crack will propagate rapidly along the crack tips (as shown in Figures 8(d)-8(f)).
The deformation and crack paths in the beam which resulted from the presented PD model are compared with the numerical results by using boundary element method (BEM) [28] and experimental observations [29], as shown in Figures 9 and 10.We may find that the deformation diagrams by using different numerical methods are similar to each other reasonably, while the crack paths by PD simulation are closer to the experimental results, especially at the end of the crack paths.What is more is that the failure load of the beam (99 kN) obtained in the present PD simulation is close to the BEM simulation result (92 kN) as well.

Conclusions
As a reformation of the classical continuum mechanics theories, the recently developed peridynamic theory has shown outstanding capabilities of modeling various discontinuous problems without any ill-suited situation.In this paper, the updated bond-based peridynamic model for concrete materials has been presented, in which the different mechanical behavior of concrete under tensile and compressive is taken into account, as well as the short-range repulsive forces between material points.
The present model is applied successfully to simulate the crack propagation process in a concrete slab with mode I and I-II mixed mode central cracks and is further employed to simulate the classical crack growth problems in a doubleedge notched concrete beam subjected to four-point shear load.Numerical results show that the present peridynamic model can simulate the entire fracture process in concrete structures effectively with great potential analogous

Figure 3 :
Figure 3: Mode I crack propagation process as an outgrowth of peridynamic equations of motion, contours of damage shown with different colors, and damage is defined at the level of pairwise bond locally.(a) No damage occurs with elastic deformation, (b) damage initiation at the crack tip, (c)-(d) damage accumulates, and crack grows symmetrically along the crack tip as a typical mode I crack will behave, and (e) crack propagates to the boundary of the slab.

Figure 4 :
Figure 4: Contours of horizontal displacement in concrete slab with mode I crack.(a)-(c) Material deforms continuously; (d)-(e) displacement field near the crack gets discontinuously distributed as crack grows.

Figure 5 :
Figure 5: I-II mixed mode crack propagation in concrete slab, with damage defined at the level of pairwise bond locally.(a) No damage occurs with elastic deformation, (b) damage initiation at the crack tip, (c)-(d) damage accumulates, and crack grows along definite orientation, and (e) crack propagates to the boundary.

Figure 6 :
Figure 6: Contours of horizontal displacement in concrete slab with I-II mixed mode crack.(a)-(b) Material deforms continuously; (c)-(e) displacement field around the crack gets discontinuously distributed as crack grows.

Figure 8 :
Figure 8: Damage maps (crack paths) in the beam with increasing external load (deformation has been enlarged for 1000 times for more clear display of the crack tip).(a) No damage appearance; (b)-(c) damage accumulation without crack propagation; (d)-(f) crack propagation.
(a) Results by using the proposed PD model (b) Results by using boundary element method (BEM)[28]

Figure 9 :
Figure 9: Comparison of the deformation of the beam.

Figure 10 :
Figure 10: Comparison of the crack paths.