Peridynamic Model for the Numerical Simulation of Anchor Bolt Pullout in Concrete

Predictive simulation of anchor pullout from concrete structures is not only a serious problem in structural mechanics but also very important in structural design safety. In the finite element method (FEM), the crack paths or the points of crack initiation usually need to be assumed in advance. Otherwise, some special crack growth treatment or adaptive remeshing algorithm is normally used. In this paper, an extended peridynamicmethodwas introduced to avoid the difficulties found in FEM, and its application on anchor bolt pullout in plain concrete is studied. In the analysis, the interaction between the anchor bolt and concrete is represented by a modified short-range force and an extended bond-level model for concrete is developed. Numerical analysis results indicate that the peak pullout load obtained and the crack branching of the anchoring system agreed well with the experimental investigations.


Introduction
Anchor bolts are very important components of load transfer in a wide range of civil engineering structures such as dams, nuclear power plants, highways, and bridges.A better understanding of the pullout behavior of anchor bolts can contribute not only to the optimization of the design of the anchor system, but also to the improvement of the durability and stability of a structure.Therefore, the pullout behavior of anchor bolts in concrete structures has become a major concern in the past three decades and a lot of experimental studies have been performed [1][2][3][4][5][6].Likewise, various numerical methods have been adopted to analyze the failure mechanism and the progressive damage of anchor bolts in concrete structures.
Among previous works, most of the researchers focused on the finite element method (FEM) [7][8][9][10][11].One of the early methods [7] includes modeling of concrete as the fournode cracked element.In this approach, the tensile fracture behavior of the concrete for a solid body containing an internal discontinuous surface is formulated by deriving an incremental formulation.As a result, the crack is distributed in each grid and the localization of microcracking cannot be obtained.In order to better simulate the mode-I fracture, Alfaiate et al. [8] utilized the interface elements which are inserted along the interelement boundaries among the concrete elements to model the cracking path.Although multiple cracks can be obtained and no special remeshing technique is required with such approach, the crack direction is still restricted and ladder-shaped crack paths are formed.
Etse [9] predicted the distribution of equivalent fracture strain at peak load of the anchor system by adopting a fracture energy-based plasticity formulation.The propagation of crack is described in terms of an equivalent plastic softening process, but still it is difficult to obtain the final cracking pattern by this approach.Xu et al. [10] simulated the crack patterns and mechanical behavior of the anchor bolt pullout in concrete by using the Mohr-Coulomb criterion with tension cut-off.In this approach, the heterogeneity of the concrete material is modeled by randomly assigning strength and elastic modulus to the elements according to Weibull's distribution and the ongoing cracking process is represented by groups and alignments of failed elements.Feenstra [11] used the smeared crack approach to study two-dimensional pullout problems.In this approach, it is assumed that a crack can be distributed over a special band in the model and its influence on the mechanical behavior of a material can be represented by adjusting the constitutive matrix irrespective of the real displacement discontinuities within the band.It is obvious that these aforementioned approaches had overcome some deficiencies of crack propagation problems, but still these approaches require remeshing or redefining of the geometry to model the progressive crack growth.Furthermore, the accuracy of results heavily relies on the complex adaptive meshing algorithms [12].Therefore, the low efficiency and accuracy of simulation are still a major problem [12][13][14].
Recent developments of mesh-free (or meshless) methods such as diffuse element method (DEM), material point method (MPM), and element free Galerkin method (EFGM) are invented to circumvent the mesh-dependence problem and relieve the volumetric locking for suitable choice of support size of shape function [15,16].Soparat and Nanakorn [14] and Coetzee et al. [17] used mesh-free methods for pullout problems.Though many discussions on the advantages of the mesh-free methods have been reported [15][16][17][18], it should still be noted that the treatment of essential boundary conditions is not straightforward as the conventional FEM and the use of shape functions of any desired order of continuity may lead to computational difficulties and complexity in deriving the coefficients of the stiffness matrix [18,19].Other contributions can be found in the literature regarding anchors' modeling using the microplane model [20,21].
Silling and Askari [22,23] proposed peridynamics as an alternative technique for the solution of crack related problems where direct displacement is used in the formulations instead of displacement derivatives, and deformation continuity is not based on assumptions.Mainly, bonds containing constitutive information of the materials are used to reproduce the nonlocal interacting forces between particles over a certain distance.Moreover, in contrast to the partial differential equations used in the classical formulation, this theory uses spatial integral equations which permit spontaneous crack occurring at multiple sites and freely extend along an arbitrary path without extra remedial techniques [24].Therefore, the peridynamic theory has obvious advantages in handling cracks problems and predicting the progressive failure process in solid mechanics.Numerous achievements have been obtained using this model during the past few years including deformation of one-dimensional bar [25,26], progressive damage of composite laminates [27,28], damage and fracture of membranes and fibers [29], static and dynamic fracture of plain and reinforced concrete structure [30][31][32][33][34], dynamic fracture in functionally graded materials [35], and coupling with classical continuum mechanics [36].
In this paper, an extended bond-level model for concrete is proposed and the interaction between the anchor bolt and concrete is represented by a modified short-range force.The analysis of anchor bolt pullout problem is carried out using the bond-based peridynamic (PD) model.Moreover, the capabilities of the improved numerical method to capture the progressive damage process and the extreme load of the anchor bolt are validated.Finally, a parametric study was performed to investigate the influence of the size of the horizon and the embedded length.Comparison of the experiments and simulations with those in the literature is also carried out.

Introduction of the Peridynamic Theory (PD)
The PD theory may be viewed as a special version of particle method or mesh-free method.It is based on assumptions that an object possesses a spatial domain R modeled as a discrete set of particles and each particle x owns a subregion Ω within a certain radius  called the material horizon as shown in Figure 1.The peridynamic equation of motion at any time  is given by Silling [22] as follows: where  denotes the mass density; b is the prescribed external body force density; u(x, ) and ü (x, ) are the displacement vector and acceleration vector, respectively.Also, f is a pairwise force vector in the peridynamic bond that represents the nonlocal interactions between the particle x and the rest of the particles.Unlike MD, in peridynamics, the magnitude of f depends upon the initial reference configuration and the state of the bond, that is, the relative position vector  and the relative displacement vector , as follows: The direction of the pairwise force vector can be expressed as  + ; thus, a general expression for the peridynamic force function f can be written as where (, ) is a scalar-valued force function  = ‖‖ and  = ‖ + ‖.For a microelastic material [22], another expression form of the pairwise force vector can be derived from a pairwise potential function  such that The general form of the linear microelastic potential is obtained as follows: where () =  0 () represents the bond elastic stiffness known as "micromodulus" function, () is the scalar-valued Boolean influence function,  0 is a constant that depends on (), and  denotes the relative elongation of a bond and it is given as If the relative elongation  = 0, then the pairwise force f between the particles does not exist.According to the above expressions, the corresponding pairwise force becomes where (, ) is a history-dependent function which depends on the value of  and it can be written as 0 denotes the critical stretch of the bond and can be obtained by mathematical derivation and by processing the classical fracture parameters.The value of local damage at point x within a peridynamic material can be denoted by the percent of broken bonds as where (, ) can be viewed as a matrix representing the damage of material point.Silling [22] introduced an original constitutive model for quasi-brittle materials; the influence function () is assumed to be 1.0, and the micromodulus function () can be correlated with the classical elastic constants through the equivalent of elastic strain energy in elasticity and PD theory.For 2D plane stress problems, as illustrated in Figure 2, from the conventional theory of linear elasticity, the strain energy density due to a uniform principal strain state and uniform shear strain state can be calculated as where  is the elastic modulus, ] is the Poisson ratio, and  0 denotes the strain.On the other hand, from the twodimensional peridynamic theory,  =  0 , the strain energy density can be calculated as [34] where  denotes the thickness.By solving the equations  , =  ,pd and  , =  ,pd , the constant  0 can be obtained as The critical stretch  0 is associated with the fracture energy   and the bond would break when the elongation goes beyond the critical value  0 .In 2D plane stress conditions, the work   required to break all the bonds per unit fracture area can be derived as in [23]: where  is the distance between the point  and the crack surface as shown in Figure 3 and  0 can be obtained by substituting () into (13) as follows:

Peridynamic Model for Concrete and Steel.
As peridynamics do not need the continuous displacement field, there are no concepts of stress or strain required in the model.Thus, the constitutive model is defined through the relationship between the bond stretch and the pairwise force among material particles or, in other words, the material damage is introduced at the bond level.Analogous to the softening function proposed by Gerstle et al. [37], a simplified softening function is introduced in the present work and the bond force and bond stretch relationship are illustrated in Figure 4.As shown in Figure 4(a), the yield compressive and tensile stretch limits of concrete are denoted as   and   , respectively.For small bond stretches,   <  <   , the bond remains in the linear elastic range.However, if  is less than the compressive stretch limit   , the pairwise force of concrete remains constant.On the other hand, an abrupt drop will happen in pairwise force if  reaches the tensile stretch limit   and the pairwise force will become   until  reaches the critical stretch   .A brittle fracture will occur at this point and the pairwise force will drop to zero.As steel may yield in tension as shown in Figure 4(b), the pairwise force of steel remains constant if  reaches the tensile stretch limit of steel   .  ,   , and   can be given as To provide required fracture energy in tension, according to (13), Substituting (11a) and (11b) into (13), the relationship between  and  can be computed as The value of  here is 0.5 as advised by Cheng et al. [35], so  value can be obtained as

Interaction between Steel and Concrete.
Stress is transferred mainly by adhesion, mechanical interaction, and friction between steel and concrete.Adhesion comes from chemical bonding and stresses are generated during curing of concrete and experimentally it is very tough to measure them.Due to the limited experimental information, it is difficult to determine the friction coefficient.As stresses due to adhesion and friction are relatively small, therefore, only mechanical interaction is considered in the present simulation.In this paper, we introduce a short-range force [23,24,38] to reproduce the interaction between the concrete and anchor bolt as where Δ denotes the distance between two types of particles and the constant  denotes the average micromodulus which is assumed to be

Discretization and Numerical Implementation.
Peridynamic equations of motion (see ( 1)) could be solved by utilizing a numerical approximation method which involves the discretization of the reference configuration into particles with a certain volume.So, the integrals in (1) can be replaced by the finite sums: where ü   denotes the acceleration of the point   at time step ,   is the total number of particles within the horizon of the point   , f(  , ) and b   are the pairwise forces and the body force density at time step , and   denotes the equivalent calculation volume (see [38]).To solve quasi-static problems by applying peridynamic equations of motion, Kilic and Madenci [39] introduced artificial damping to attain a steady-state solution.The Adaptive Dynamic Relaxation (ADR) scheme proposed by Underwood [40] is used to determine the most effective damping coefficient at each time step , as follows: where F   is the resultant force density vector;   is the damping coefficient at the th iteration and   is the modified density at the point   and they can be given as  Here, Δ is the time step size, e is a unit vector along , , or , and 1 K   is the diagonal stiffness matrix of the system, given as Finally, with the assumptions that u  0 ̸ = 0 and u  0 = 0, velocities and displacement at point   for the next time step can be obtained by central difference explicit integration as follows: Both the numerical algorithm and constitutive modeling are implemented in Fortran-90 language based on Visual Studio using an in-house peridynamic code.

Peridynamic Results for Anchor Bolt
Pullout in Concrete  lists the geometrical parameters of the experiment.Figure 5 shows the corresponding geometrical features of the prediction model.W and L are the width and length of the concrete block as shown in Figure 5. a is the support span and it is the distance from the edge of the anchor to the support and d is the embedded depth.The breadth and thickness of the anchor bolt are b and t, respectively.All the concrete blocks and anchor bolts are 0.1 m thick, which is very small as compared to other dimensions, so the model is treated as a plane stress problem in the analysis.In most of the past research works [9, 10, 14], steel bolt was also not directly modeled and there were no experimental values of friction coefficients.Pullout load is modeled by controlling the vertical displacement of the top surface of the bolt head and for the convenience of simulations only the anchor head is assumed to be in contact with the concrete.As a result, both the friction resistance caused by the pullout movement and the vertical deformation of the anchor bar were ignored.In the present study, both the plain concrete and the steel anchor are modeled as discrete particles of corresponding volumes.In addition, a shortrange force is used to model the frictionless contact between concrete and the anchor bolt.To reduce the computation cost and modeling difficulty, a constant in-plane grid number of 240 is used for each prediction model so that the grid spacing of three groups is set to be 1.25 mm, 2.5 mm, and 3.75 mm, respectively.Representative of quasi-static loading, a very low velocity of 2.4 × 10 −9 m/s was applied on the top particles of the anchor bolt.The material properties are depicted in Table 2 and the corresponding numerical models are shown in Figure 6.

Results and Discussions.
To investigate the influence of the size of material horizon on the load-displacement response of the anchor bolt, three different material horizon sizes (i.e.,  = 3 * , 4 * , and 5 * ) were used.Figure 7 shows the relationship between the vertical displacement of the top particles and the pullout force in Group 2. It is observed that the results obtained with different material horizons approximately give the same peak load.For the  Experiment by Vervuurt et al. [6] Simulated by Vervuurt et al. [6] Simulated by Feenstra [11] Simulated by Soparat and Nanakorn [14] Displacement (휇m) This paper (Group 2: 훿 = ) 3 * dx larger material horizon, a flatter initial slope and greater peak displacement are witnessed.As the material horizon is related to critical stretch for the proposed concrete model (see (18)), a lower critical stretch shows that the bond is more brittle than the material.As shown in Figure 8, the curve is more close to the experiment results by Vervuurt et al. [6] when  = 3 *  is adopted; therefore, the default material horizon is set to be three times the grid spacing in the following simulation.
Figure 8 also shows the lattice model, FE, and EFG results by Vervuurt et al. [6], Feenstra [11], and Soparat and Nanakorn [14], respectively.The peridynamic model shows very good results in the upper section of the experimental   curve for load-displacement behavior.Figure 9 shows loaddisplacement curves for all groups.The predicted peak loads and experimental values comparison is given in Table 3.
Embedded depth has an influence on the pullout mechanism; however, the peak load does not increase proportionally with the embedded depth which shows that size has effect on pullout loads [6].
The values of the predicted peak loads in the present study are a little higher than the experimental values which may be because of the higher Poisson's ratio; 0.33 instead of normal 0.3, based on the peridynamic theory [22,23] which can be seen in Figure 8 and Table 3.
Although the precision obtained by the current model is not very high as compared to EFGM in peak load prediction, it was based on many simplifications [14].(a) Half of the model was considered due to symmetry.(b) Steel bolt was not directly modeled.(c) Supports were treated as points.(d) Crack propagation was assumed from the upper edge corner of the bolt head, among others.These simplifications may affect the reliability of simulation results.In peridynamics, both support and specimen are modeled as particles, so simulations are more reliable than in [14].Cracks can appear spontaneously as the material response includes damage without the need for any special criterion.Crack trajectories of the present model are compared with the EFGM results by Soparat and Nanakorn [14] for different groups as shown in Figure 10.The crack branching can be easily captured in PD theory while there is no localized branching in EFGM.
The crack pattern from experimental results is different from most of the previous numerical methods results [6,11,14].The reasons can be the assumption of a 2D model instead of the real 3D complex fracture model for simulations, heterogeneous nature of concrete [6,10], and experimental errors in the tests [10].Figure 11 compares the failure mode of the present study with the experimental observations [6].It is worth noting that both the crack direction and the crack branching obtained in the present study show good resemblance with the experimental result by Vervuurt et al. [6].

Conclusions
In this paper, the problem of anchor bolt pullout in plain concrete was investigated with an improved peridynamic model.An extended constitutive model was used to refine the behaviors of concrete material simulation.A short-range force was introduced to simulate anchor bolt and concrete interaction in 2D.The numerical discretization and iteration algorithms were implemented with an in-house peridynamic FORTRAN code.It was observed that the crack propagation of concrete was exposed in more detail by the proposed approach as compared to conventional FEM or EFGM.Compared with the results from the literatures and experiments, it can be concluded that the extreme failure load and the final failure mode of the anchor bolt by analysis of the peridynamic approach match well those of the experimental observations.From all the results and comparisons, this approach was proved to be a promising method for solving the problem of anchor bolt pullout in plain concrete.

Figure 1 :
Figure 1: Pairwise interaction between  and   at time .

Figure 2 :
Figure 2: A two-dimensional plate subjected to uniform deformation.

Figure 4 :
Figure 4: Constitutive model of concrete and steel.

Figure 5 :
Figure 5: Specimen geometry and the pullout test setup.

Figure 7 :
Figure 7: Load-displacement response with different horizon sizes.

Figure 8 :
Figure 8: Comparison of the numerical load-displacement curve with past research work.

Figure 9 :
Figure 9: Load-displacement response for different groups.

Table 1 :
Geometrical parameters used in the simulation.

Table 3 :
Comparison of the peak loads.