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

Authors: Jiezhi Lu; Yaoting Zhang; Habib Muhammad; Zhijun Chen
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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 peridynamic method was 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.


---

## Body

## 1. 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–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–11]. One of the early methods [7] includes modeling of concrete as the four-node 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–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–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–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.


## 2. 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 domainR 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 t is given by Silling [22] as follows:(1)ρu¨x,t=∫Ωfη,ξdVx′+bx,t,(2)Ω=x′∈R:ξ=x-x′<δ,where ρ denotes the mass density; b is the prescribed external body force density; ux,t and u¨x,t 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:(3)η=u′-u=ux′,t-ux,t,ξ=x′-x.Figure 1
Pairwise interaction betweenx and x′ at time t.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(4)fξ,η=fξ,ηξ+ηζifξ∈Ω,where f(ξ,η) 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(5)fξ,η=∂ωξ,η∂η.The general form of the linear microelastic potential is obtained as follows:(6)ωξ,η=cξs2ξ,ηξ2,where c(ξ)=c0κ(ξ) represents the bond elastic stiffness known as “micromodulus” function, κ(ξ) is the scalar-valued Boolean influence function, c0 is a constant that depends on κ(ξ), and s denotes the relative elongation of a bond and it is given as(7)sξ,η=ζ-ξξ.If the relative elongation s=0, then the pairwise force f between the particles does not exist. According to the above expressions, the corresponding pairwise force becomes(8)fξ,η=c0sξ,ηκξμt,ξξ+ηζifξ∈Ω0otherwise,where μ(t,ξ) is a history-dependent function which depends on the value of s and it can be written as(9)μt,ξ=1ifst′,ξ,η<s0,0≤t′≤t0otherwise.s0 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(10)ϕx,t=1-∫Ωμt,ξdVx′∫ΩdVx′,where ϕ(x,t) 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 c(ξ) 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(11a)Un,e=Eε021-ν,Us,e=Eε021+ν,where E is the elastic modulus, ν is the Poisson ratio, and ε0 denotes the strain. On the other hand, from the two-dimensional peridynamic theory, s=ε0, the strain energy density can be calculated as [34](11b)Un,pd=c0ε02πtδ36,Us,pd=c0ε02πtδ312,where t denotes the thickness. By solving the equations Un,e=Un,pd and Us,e=Us,pd, the constant c0 can be obtained as(12)c0=6E1-νtπδ3ifξ∈Ω.The critical stretch s0 is associated with the fracture energy GF and the bond would break when the elongation goes beyond the critical value s0. In 2D plane stress conditions, the work GF required to break all the bonds per unit fracture area can be derived as in [23]:(13)GF=2∫0δ∫zδ∫0cos-1z/ξcξs02ξ2ξdφdξdz,where z is the distance between the point x and the crack surface as shown in Figure 3 and s0 can be obtained by substituting c(ξ) into (13) as follows: (14)s0=4πGF9Eδ.Figure 2
A two-dimensional plate subjected to uniform deformation.Figure 3
Evaluation of fracture energy.


## 3. Peridynamic Model for Anchor Bolts

### 3.1. 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.Figure 4
Constitutive model of concrete and steel.
(a)
Concrete (b)
SteelAs shown in Figure4(a), the yield compressive and tensile stretch limits of concrete are denoted as sc and st, respectively. For small bond stretches, sc<s<st, the bond remains in the linear elastic range. However, if s is less than the compressive stretch limit sc, the pairwise force of concrete remains constant. On the other hand, an abrupt drop will happen in pairwise force if s reaches the tensile stretch limit st and the pairwise force will become αft until s reaches the critical stretch βst. 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 s reaches the tensile stretch limit of steel sy. st, sc, and sy can be given as(15)sy=fy′E,sc=-fc′E,st=ft′E.To provide required fracture energy in tension, according to (13),(16)GF=2∫0δ∫zδ∫0cos-1z/ξcξst2ξ22+αftβ-1stξ2dθdξdz.Substituting (11a) and (11b) into (13), the relationship between α and β can be computed as(17)1+2αβ-1=4πGF9Est2δ.The value of α here is 0.5 as advised by Cheng et al. [35], so β value can be obtained as(18)β=4πGF9Est2δ.


### 3.2. 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(19)fξ,η=ξ+ηζmin⁡0,c-ζΔx-1,where Δx denotes the distance between two types of particles and the constant c- denotes the average micromodulus which is assumed to be(20)c-=cconcrete+csteel2.


### 3.3. 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:(21)ρu¨ni=∑j=1Nifηn,ξVij+bni,where u¨ni denotes the acceleration of the point xi at time step n, Ni is the total number of particles within the horizon of the point xi, f(ηn,ξ) and bni are the pairwise forces and the body force density at time step n, and Vij 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 n, as follows: (22)u¨ni+cnu˙ni=Fniλi,where Fni is the resultant force density vector; cn is the damping coefficient at the nth iteration and λi is the modified density at the point xi and they can be given as(23)cn=2uniTKnii1uniuniTuni,λi≥14Δt2∑j=1Niξ·e·cξξ2.Here,Δt is the time step size, e is a unit vector along x, y, or z, and Knii1 is the diagonal stiffness matrix of the system, given as(24)Knii1=-Fni/λi-Fn-1i/λiΔtu˙n-1/2i.Finally, with the assumptions thatu0i≠0 and u˙0i=0, velocities and displacement at point xi for the next time step can be obtained by central difference explicit integration as follows:(25)u˙1/2i=ΔtF0i2λi,u˙n+1/2i=2-cnΔtu˙n-1/2i+2ΔtFni/λi2+cnΔt,un+1i=uni+Δtu˙n+1/2i.Both the numerical algorithm and constitutive modeling are implemented in Fortran-90 language based on Visual Studio using an in-house peridynamic code.


## 3.1. 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.Figure 4
Constitutive model of concrete and steel.
(a)
Concrete (b)
SteelAs shown in Figure4(a), the yield compressive and tensile stretch limits of concrete are denoted as sc and st, respectively. For small bond stretches, sc<s<st, the bond remains in the linear elastic range. However, if s is less than the compressive stretch limit sc, the pairwise force of concrete remains constant. On the other hand, an abrupt drop will happen in pairwise force if s reaches the tensile stretch limit st and the pairwise force will become αft until s reaches the critical stretch βst. 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 s reaches the tensile stretch limit of steel sy. st, sc, and sy can be given as(15)sy=fy′E,sc=-fc′E,st=ft′E.To provide required fracture energy in tension, according to (13),(16)GF=2∫0δ∫zδ∫0cos-1z/ξcξst2ξ22+αftβ-1stξ2dθdξdz.Substituting (11a) and (11b) into (13), the relationship between α and β can be computed as(17)1+2αβ-1=4πGF9Est2δ.The value of α here is 0.5 as advised by Cheng et al. [35], so β value can be obtained as(18)β=4πGF9Est2δ.


## 3.2. 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(19)fξ,η=ξ+ηζmin⁡0,c-ζΔx-1,where Δx denotes the distance between two types of particles and the constant c- denotes the average micromodulus which is assumed to be(20)c-=cconcrete+csteel2.


## 3.3. 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:(21)ρu¨ni=∑j=1Nifηn,ξVij+bni,where u¨ni denotes the acceleration of the point xi at time step n, Ni is the total number of particles within the horizon of the point xi, f(ηn,ξ) and bni are the pairwise forces and the body force density at time step n, and Vij 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 n, as follows: (22)u¨ni+cnu˙ni=Fniλi,where Fni is the resultant force density vector; cn is the damping coefficient at the nth iteration and λi is the modified density at the point xi and they can be given as(23)cn=2uniTKnii1uniuniTuni,λi≥14Δt2∑j=1Niξ·e·cξξ2.Here,Δt is the time step size, e is a unit vector along x, y, or z, and Knii1 is the diagonal stiffness matrix of the system, given as(24)Knii1=-Fni/λi-Fn-1i/λiΔtu˙n-1/2i.Finally, with the assumptions thatu0i≠0 and u˙0i=0, velocities and displacement at point xi for the next time step can be obtained by central difference explicit integration as follows:(25)u˙1/2i=ΔtF0i2λi,u˙n+1/2i=2-cnΔtu˙n-1/2i+2ΔtFni/λi2+cnΔt,un+1i=uni+Δtu˙n+1/2i.Both the numerical algorithm and constitutive modeling are implemented in Fortran-90 language based on Visual Studio using an in-house peridynamic code.


## 4. Peridynamic Results for Anchor Bolt Pullout in Concrete

### 4.1. Problem Setup

Because of the distinctive advantages in solving crack propagation problems, the peridynamic method is adopted in the present study to model the anchor bolt pullout in concrete. The anchor bolt pullout experiment conducted by Vervuurt et al. [6] is considered here. Table 1 lists the geometrical parameters of the experiment. Figure 5 shows the corresponding geometrical features of the prediction model.Table 1
Geometrical parameters used in the simulation.
Group ID W (mm) L (mm) a (mm) d (mm) b (mm) t (mm) 1 300 300 100 50 15 5 2 600 600 200 100 30 10 3 900 900 300 150 45 15Figure 5
Specimen geometry and the pullout test setup.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 short-range 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.Table 2
Material parameters used in the simulation.
Material type Concrete Steel Young’s modulusE (GPa) 30 200 Yield or tensile strengthft′ (MPa) 3 400 Compressive strengthfc′ (MPa) 30 - Fracture energyGF (N/m) 100 - Densityρ (kg/m3) 2500 7890Figure 6
Numerical model.


### 4.2. 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∗dx, 4∗dx, and 5∗dx) 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 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∗dx is adopted; therefore, the default material horizon is set to be three times the grid spacing in the following simulation.Figure 7
Load-displacement response with different horizon sizes.Figure 8
Comparison of the numerical load-displacement curve with past research work.Figure8 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 load-displacement 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].Table 3
Comparison of the peak loads.
Group ID This paper (KN) Experimental results (KN) Relative error (%) 1 17.04 13.4 27.2% 2 28.86 24.5 17.8% 3 37.48 33.6 11.6%Figure 9
Load-displacement response for different groups.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.Figure 10
Comparison of crack patterns (PD versus EFGM).
(a)
Group 1 (d=50mm, displacement × 100) (b)
Group 2 (d=100mm, displacement × 100) (c)
Group 3 (d=150mm, displacement × 100)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].Figure 11
Comparison of failure modes (PD versus EXP).
(a)
Group 1 (d=50mm, displacement × 50) (b)
Group 2 (d=100mm, displacement × 50) (c)
Group 3 (d=150mm, displacement × 50)


## 4.1. Problem Setup

Because of the distinctive advantages in solving crack propagation problems, the peridynamic method is adopted in the present study to model the anchor bolt pullout in concrete. The anchor bolt pullout experiment conducted by Vervuurt et al. [6] is considered here. Table 1 lists the geometrical parameters of the experiment. Figure 5 shows the corresponding geometrical features of the prediction model.Table 1
Geometrical parameters used in the simulation.
Group ID W (mm) L (mm) a (mm) d (mm) b (mm) t (mm) 1 300 300 100 50 15 5 2 600 600 200 100 30 10 3 900 900 300 150 45 15Figure 5
Specimen geometry and the pullout test setup.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 short-range 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.Table 2
Material parameters used in the simulation.
Material type Concrete Steel Young’s modulusE (GPa) 30 200 Yield or tensile strengthft′ (MPa) 3 400 Compressive strengthfc′ (MPa) 30 - Fracture energyGF (N/m) 100 - Densityρ (kg/m3) 2500 7890Figure 6
Numerical model.


## 4.2. 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∗dx, 4∗dx, and 5∗dx) 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 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∗dx is adopted; therefore, the default material horizon is set to be three times the grid spacing in the following simulation.Figure 7
Load-displacement response with different horizon sizes.Figure 8
Comparison of the numerical load-displacement curve with past research work.Figure8 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 load-displacement 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].Table 3
Comparison of the peak loads.
Group ID This paper (KN) Experimental results (KN) Relative error (%) 1 17.04 13.4 27.2% 2 28.86 24.5 17.8% 3 37.48 33.6 11.6%Figure 9
Load-displacement response for different groups.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.Figure 10
Comparison of crack patterns (PD versus EFGM).
(a)
Group 1 (d=50mm, displacement × 100) (b)
Group 2 (d=100mm, displacement × 100) (c)
Group 3 (d=150mm, displacement × 100)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].Figure 11
Comparison of failure modes (PD versus EXP).
(a)
Group 1 (d=50mm, displacement × 50) (b)
Group 2 (d=100mm, displacement × 50) (c)
Group 3 (d=150mm, displacement × 50)


## 5. 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.



---

*Source: 8412620-2018-03-21.xml*