A Non-Ordinary State-Based Peridynamic Formulation for Failure of Concrete Subjected to Impacting Loads

Strain hardening and strain rate play an important role in dynamic deformation and failure problems such as high-velocity impact cases. In this paper, a non-ordinary state-based peridynamic model for failure and damage of concrete materials subjected to impacting condition is proposed, taking the advantages of both damage model and nonlocal peridynamic method. The Holmquist-Johnson-Cook (HJC) model describing the mechanical character and damage of concrete materials under large strain, high strain rate and high hydrostatic pressure was reformulated in the framework of non-ordinary statebased peridynamic theory, and the corresponding numerical approach was developed. The proposed model and numerical approach were validated through simulating typical impacting examples and comparing the results with available experimental observations and results in literature.


Introduction
Concrete is a kind of widely used material in protective engineering and infrastructures, which attracts interest of researchers from nuclear industry and fortification installations [Corbett, Reid and Johnson (1996) ;Børvik, Langseth, Hopperstad et al. (2002)]. Since the progressive failure of internal defects will lead to the weakness of concrete structures and probably will bring a sudden destruction without visible damage, it has become a long-standing and rough challenge for researchers to study the mechanism of damage evolution, crack nucleation and propagation in concrete materials and structures. In order to ensure the structural safety, especially under impact and fatigue loadings, one should accurately describe the destruction mechanism of concrete materials and structures. Many researches on the mechanical behavior of concrete under dynamic condition have been conducted and some simplified analytical models are determined by experiments [Kennedy (1976); Ben-Dor, Dubinsky and Elperin (2005) ;Li, Reid, Wen et al. (2005)]. On the other hand, numerical simulations with different constitutive models of concrete were used in penetration and impact problems [Teng, Chu, Chang et al. (2005); Leppä nen (2006); Tai and Tang (2006)]. In addition, different constitutive models for concrete under dynamic loadings were investigated in Bićanić et al. [Bićanić and Zienkiewicz to the solid mechanics problems under dynamic condition. The viscoplasticity model with non-ordinary state-based peridynamic theory was proposed by Foster et al. [Foster, Silling and Chen (2010) ;Foster, Silling and Chen (2011)]. The stress-based failure criterion was implemented into the non-ordinary state-based peridynamic theory by Zhou et al. [Zhou, Wang and Xu (2016)] to simulate the initiation, propagation and coalescence process of cracks under quasi-static and dynamic loads. Oterkus et al. [Oterkus, Madenci and Agwai (2014)] studied heat conduction problem using state-based peridynamic theory and employed the state-based peridynamics in the formulation for thermoplastic fracture [Amani, Oterkus, Areias et al. (2016)]. Lai et al. [Lai, Liu, Li et al. (2018)] formulated the fracture of quasi-brittle materials under dynamic condition. Fan et al. [Fan, Bergel and Li (2016); Fan and Li (2017)] performed soil fragmentation simulation under blast loads of buried explosive using hybrid PD-SPH method, and Yaghoob et al. [Yaghoobi and Mi (2017)] developed numerical methods in order to improve the stability of peridyamic calculation by suppressing zero-energy mode. For the damage of concrete materials and structures, Kilic et al. [Kilic and Madenci (2009)] analyzed the influence of interface size and boundary conditions on the critical destabilization load for concrete target. Huang et al. [Huang, Lu and Qiao (2015)] proposed a new method for force loading and a more accurate constitutive model by introducing a quasi-static solution algorithm for local damping and unbalanced force convergence criteria. Gerstle et al. [Gerstle, Sau and Silling (2005); Gerstle, Sau and ] simulated the progressive failure process of concrete and reinforced concrete structures under tension, compression, shear and combined loading conditions. As to dynamic deformation and simulation of concrete, the effects of strain hardening and strain rate play an important role. Nevertheless, these factors have seldom been taken into account in previous peridynamic simulations. In this paper, the specific purpose is to implement a traditional constitutive model into the non-ordinary state-based peridynamics framework, thus making use of its advantages in simulating impact problems of concrete, considering the effect of strain hardening and strain rate. To accurately characterize the damage progress of concrete, the Holmquist-Johnson-Cook (HJC) model [Holmquist, Johnson and Cook (1993) ;Johnson, Beissel, Holmquist et al. (1998)] was employed in the present work. The HJC model shows advantages on concrete damage description and large-scale computations, and has been successfully implemented into LS-DYNA, a general-purpose finite element program developed by LSTC (Livermore Software Technology Corporation), for penetration simulations before. The remainder of this paper is organized as following. In Section 2 the framework of HJC constitutive model is introduced. The brief concept of the peridynamics and the nonordinary state-based peridynamic theory are presented in Section 3. Section 4 presents the numerical approaches including the discretization of the equation of motion, short range force and the failure criterion. The proposed model and approach are validated in Section 5 with two benchmark numerical examples. Finally, some concluding remarks are drawn in Section 6. , vol.118, no.3, pp.561-581, 2019 2 The HJC model In the HJC constitutive model, the normalized equivalent stress is determined by the constitutive relation, as shown in Fig. 1. ( ) ( )

CMES
in which * σ is the normalized equivalent stress, * =/ c σ σ f .  and c f is the actual equivalent stress and the uniaxial compressive strength respectively. * P is the normalized pressure, and P is the hydrostatic pressure, * / c P P f = . * 0 / ε ε ε = is the dimensionless strain rate, in which ε is the actual strain rate and -5 1 0 =10 s ε − is the reference strain rate. In addition, A , B , C and N are material parameters, which are determined by test data. A is the cohesive parameter, B is the pressure hardening coefficient, C is the strain-rate sensitivity coefficient and N is the pressure hardening exponent. Material degradation is described by the damage variable D , leading to reduction of the cohesive strength.
where P ε  and P μ  are the effective plastic strain increment and plastic volumetric strain during a cycle of integration, respectively.
ff pp ε μ + is the plastic strain until concrete fractures under a constant pressure, and EFMIN is a value of the minimum plastic strain causing the fracture of the material, which the value of 0.01 is chosen in our work. 1 D and 2 D represent the damage constants. 1 D is the date from unconfined compression test, while 2 D is chosen to be 1.0, which assumes that plastic fracture strain increases linearly with increasing pressure.   CMES, vol.118, no.3, pp.561-581, 2019 where lock P is the compacting pressure, and lock μ is the locking volumetric strain.
In the third phase (BC), all air voids in concrete are crushed out of it. Concrete is locked and cannot be compressed any further in either plastic deformation or void collapse. Therefore, this region can be assumed completely non-linear elastic.
is the modified volumetric strain, 1 k , 2 k and 3 k are the material constants.

Basic theory
A brief description of peridynamics by Silling [Silling (2000)] is presented below. The peridynamics describes the dynamics process of a body from its reference configuration to the current configuration. A schematic of the body is shown in Fig. 4. In the peridynamics, the motion equation for any material point x in the reference configuration at time t ( 0 t  ) is given as where  and b are the mass density and external applied body force density respectively,

( )
,t ux is the displacement field, the initial conditions ( ) ( ) Lx is a function of displacement, which represents the internal force density (per unit volume) that is exerted on x by other body-points [Silling and Lehoucq (2008)]. In the peridynamic theory, the motion of the body is presented by considering the interaction of any material point, x , with the other material points,  x , within a horizon x H . The size of horizon of a given point x is finite, defined as , where 0   . The relative position vector between points x and  x is referred to as the bond, which is denoted by ξ , and is defined as where  is the relative displacement of two interacting material points x and  x , which is defined as: However, in the state-based peridynamic theory, the deformation of the bond is described by the deformation vector state Y .
The equation of motion in state-based peridynamics can be written as: Where   ,t Tx is the force vector state representing the relationship between material points at time t . To make the notation more concise, (6) will be abbreviated as:

Constitutive model of non-ordinary state-based peridynamics
To calculate the force-vector state T mentioned above, the non-local deformation gradient tensor F of point x is calculated firstly as the following expression Where () ω ξ is the influence function of the bond and K is the non-local shape tensor defined by It should be noted that the nonlocal shape tensor K is a symmetric and positive definite second order tensor [Silling, Epton, Weckner et al. (2007) And then, the velocity gradient tensor is given by which can be decomposed into symmetric and skew-symmetric parts. The further one is referred to as the rate of deformation tensor D and is shown as follows: If the polar decomposition theory is applied to F , it can be expressed by

F = VR = RU
(18) where R is an orthogonal tensor proposed by Flanagan et al. [Flanagan and Taylor (1987)], representing a rigid-body rotation. And then, V and U are the left and right stretch tensors respectively. The velocity gradient tensor can be defined in another way by substituting the right polar decomposition from Eq. (18) into Eq. (16) T The term T RR in Eq. (19) is skew symmetric and describes a rate of rotation. And for completeness, the unrotated rate of deformation tensor d will be employed.
where t R describes the rigid-body rotation at current time t , calculated by the incremental formulation as follow: . ikj e is the permutation tensor, and the axial vector ω is given as V is the left stretch tensor. The material is firstly assumed to be elastic, then elastic strain increment tensor and deviatoric strain increment tensor are calculated as The trial unrotated Cauchy stress at time t is defined as According to the von-Mises plasticity theory, we can get the trial deviatoric stress tensor and the von-Mises yield stress as follows: ( ) Then, the rotated Cauchy stress tensor σ is based on the unrotated Cauchy stress as follow where n is the number of time steps, 3 = p V x stands for the involved volume of particle where t  is the time step, and in order to obtain a stable numerical result, the time step t  should be satisfied the following inequation [Warren, Silling, Askari et al. (2009) where δ denotes the size of horizon, represents the dilatational wave speed,  and  are the Lame's elastic constants of the material.

Short range force
When simulating the movement of the material point in the target substance, a shortrange force model [Silling and Askari (2005)] is employed to prevent target material points from penetrating each other, as shown in Fig. 5.

Failure criterion
In this paper, damage is given as the ratio shown in Eq. (36), which is relevant to the amount of broken bonds and the total amount of bonds in the horizon as: where s denotes the stretch of a bond, , where 0 G is energy release rate [Silling and Askari (2005)].

Numerical results
In this section, two benchmark examples are considered to validate and demonstrate the proposed peridynamic model. In the first example, the fracture of a single-edge notched, three-point-bending concrete beam under impact loading is analyzed. The experimental results were reported by John et al. [John and Shah (1986)] and John [John (1988)]. In the second example, the penetration experiment conducted by Hanchak et al. [Hanchak, Forrestal, Young et al. (1992) CMES, vol.118, no.3, pp.561-581, 2019

Crack propagation in a three-point-bending specimen under dynamic condition
In the three-point-bending impact test, the concrete specimen with an initial notch is subjected to impact loading (see Fig. 6). The experiment reported the crack propagation paths for several test specimens with the offset notch at different locations, and Belytschko et al. predicted the crack propagation angles and time data by element-free Galerkin (EFG) method [Belytschko, Organ and Gerlach (2000)]. The location of the notch in the concrete specimen is described by a normalized parameter γ , which is the ratio of the distance from the notch to the midspan to the distance from one support to the midspan. And the height of the notch is decided by a parameter β , which is the ratio between notch height and specimen height.
The following parameters are same as the experiments: the distance from the midspan to the support is 101 mm, the distance from the support to the edge is 13 mm, the height of the specimen is 76 mm, and the thickness of the specimen is 25.4 mm. To effectively decrease the inertial oscillations, the impact velocity is increased linearly to the maximum value 1 v and then held unchanged  CMES, vol.118, no.3, pp.561-581, 2019 Figure 9: Crack propagation paths of the three-point-bending specimen (experimental and LEFM results [Belytschko, Organ and Gerlach (2000)]) To determine the relationship between crack propagation path and different location of notch, the results of cracking angle and crack initiation time with different γ are listed in Tab. 2. It shows that no matter where the pre-existed notch located, all the cracks propagate from the notch tip towards the loading point. As the γ increases, the initial angle and crack initiation time will increase.

Hanchak penetration experiment
In the Hanchak penetration experiment, a rectangular block of reinforced concrete suffers from projectile impact. Fig. 10 and Fig. 11 show the geometry of the target and location of the reinforcement. The impact point is not the center of the specimen so as to avoid the collision between the projectile and steel bars. For the sake of numerical calculation, the models are discretized into particles with the uniform grid spacing 5 mm x = . In the numerical model, the specimen is discretized to 570972 material points. The stable time step is adopted of 8.0×10 -7 s and the horizon size is three times of grid spacing, =15 mm  Fig. 13 and Fig. 14, with damage maps of the impact surface and damage maps of the exit surface respectively. It shows that the numerical results on the pattern of damage distribution agree well with the experimental observations. , vol.118, no.3, pp.561-581, 2019 (a) impact surface (b) exit surface  Fig. 16. It shows that the numerical results match well with experimental data, indicating that the proposed peridynamic approach is capable of analyzing this kind of impacting problems.

Figure 16:
Comparison of residual velocities of the projectile with various striking velocities between the numerical results and experimental results [Hanchak, Forrestal, Young et al. (1992)]

Conclusions
In this work, a non-ordinary state-based peridynamic model for concrete failure process simulation under impact loading is presented by reformulating the HJC constitutive model for concrete in continuum mechanics theory. In the proposed model, the strain hardening, strain-rate effects, and pressure-dependence is taken into account, to characterize the damage and dynamic fracture of concrete under impacting loads. Two benchmark problems have been studied to verify the proposed model and approach.
In the dynamic fracture of a three-point-bending beam, the crack propagation paths of the specimen subjected to impact loading were investigated, and crack propagation, 578 Copyright © 2019Tech Science Press CMES, vol.118, no.3, pp.561-581, 2019 orientation and cracking time with offset notches at different locations were analyzed. The results show that when the pre-existing notch is closer to edge of the specimen, the cracking angle and initiation time will increase. This agrees well with available experimental observations and results from other numerical methods. In the Hanchak penetration test, the evolution of damage and stress wave in the target at a nominal striking velocity of 749 m/s was investigated. The residual velocities of the projectile with different impacting velocities were predicted by using the proposed model, and which keeps good accordance with experimental results. This indicates the capability of the proposed non-ordinary state-based peridynamic model to simulate the deformation, damage and fracture of concrete structures subjected to dynamic loads such as impacting.