A Rate-Dependent Peridynamic Model for the Dynamic Behavior of Ceramic
Materials

In this study, a new bond-based peridynamic model is proposed to describe the dynamic properties of ceramics under impact loading. Ceramic materials show pseudo-plastic behavior under certain compressive loadings with high strain-rate, while the characteristic brittleness of the material dominates when it is subjected to tensile loading. In this model, brittle response under tension, softening plasticity under compression and strain-rate effect of ceramics are considered, which makes it possible to accurately capture the overall dynamic process of ceramics. This enables the investigation of the fracture mechanism for ceramic materials, during ballistic impact, in more detail. Furthermore, a bond-force updating algorithm is introduced to perform the numerical simulation and solve the derived equations. The proposed model is then used to analyze the dynamic response of ceramics tiles under impact loading to assess its validity. The results of damage development in ceramic materials are calculated and compared with the experimental results. The simulation results are consistent with the experiments, which indicates that the proposed rate-dependent peridynamic model has the capability to describe damage propagation in ceramics with good accuracy. Finally, based on a comparison between simulation and experimental results, it can be concluded that the damage results are in better agreement with experimental results than non-ordinary state-based peridynamic method.


Introduction
Ceramic materials are known to work well in protective armor systems such as armored vehicles, bodyarmor for soldiers, and spacecraft modules. Specifically, they provide effective protection against the penetration of high-speed projectiles. Their unique low-density and high compressive-strength make them also suitable for lightweight protective structures [1]. However, ceramics are inherently brittle and exhibit low fracture-toughness compared to metals. In addition, extensive fragmentation of the ceramic material occurs when subjected to impact from a projectile at high speed. It is necessary to investigate this fragmentation process.
The failure events are not necessarily localized to the immediate impact zone [2]. When a projectile penetrates a ceramic plate, the compressive shock waves are generated at the projectile/ceramic interface, and propagate from the front surface to the rear surface of the tile. When the compressive waves reach the rear surface of the ceramic tile, they will turn into tensile waves due to the reflection at the free surface, and initial cracks occur because of the interaction of high tensile stress. Then, the ceramic tile fractures due to nucleation, propagation, and interaction of the cracks during the penetration. Consequently, there is a comminuted fracture zone in the ceramic tile, which has the shape of a conoid [2]. To design and optimize armor applications, it is essential to have a deep understanding of the dynamic responses of ceramics under impact loading, especially their fracturing behavior.
Until now, many studies that involved experimental observation, theoretical considerations and numerical simulations were conducted to improve the understanding of the dynamic behavior of brittle materials. However, these experimental observations were limited by the accuracy of measurements over very short periods as well as the ability to link experimental observations with meaningful interpretations [1]. Furthermore, theoretical considerations can currently only solve certain simple scenarios and models. Therefore, several researchers developed different numerical techniques to model and predict the dynamic response of brittle materials. Several computational methods can now capture the characteristics of dynamic fracture in brittle materials-but each method has their benefits and drawbacks [3].
Hughes et al. [4] used the finite element method to solve a class of contact-impact problems, which requires finding the solution to a nonlinear algebraic system of equations. In classical continuum mechanics, the material properties are assumed to remain continuous. Crack propagation plays an important role in the failure analysis of brittle materials [5]. However, the spatial partial derivatives become unsuitable along the crack, where the displacement field is discontinuous [6]. Currently-available numerical continuum methods, which use a mesh or grid structure (such as Finite Difference Methods (FDM), Finite Element Method (FEM)) cannot handle discontinuous problems effectively. Thus, modified methods like the Extended Finite Element Method (XFEM) were developed to deal with discontinuous problems. The XFEM, developed by Belytschko et al. [7] and Moёs et al. [8], allows cracks to pass through the finite element. This is a relatively accurate method to model cracks of arbitrary geometries using the finite method without remeshing. However, it fails to predict experimentally-observed crackpropagation speeds without adjustment of the energy release rate. The mesh-free methods to solve such problems, e.g., the Smoothed Particle Hydrodynamics (SPH) and Element-free Galerkin (EFG), are more efficient but may suffer from instabilities and low efficiency when handling more complicated problems. Although significant progress has been made in computational mechanics, the simulation of the dynamic response of ceramics has been a great challenge for a long time. This is due to the complexity of the involved failure mechanism.
In recent years, a non-local continuum theory, called peridynamics (PD), was proposed by Silling [9] to handle discontinuous problems such as the crack propagation and damage. In contrast to the classical continuum theories, peridynamics reformulates the basic equation of continuum in the form of integration rather than differential equation, which can be better suited to simulate discontinuities. In PD theory, the solid domain is discretized into material points, and the material point interacts with other material points within a finite distance called horizon-a domain associated to the material point. Peridynamic formulations can be divided, generally, into two distinct branches, bond-based peridynamics and statebased peridynamics.
Most past research effort concentrated on the state-based formulation because it is possible to use very general material models and accommodate any constitutive relations, using a classic continuum framework. A non-ordinary state-based peridynamic model was developed for the fracture analysis of brittle materials, in which the modified Johnson Holmquist (JH-2) model [10] was used to characterize the plastic softening behavior of ceramic materials [11]. Wu et al. [12] proposed a non-ordinary state-based peridynamic model to investigate the failure and damage of concrete materials under impact loadings. An improved stated-based peridynamic lattice model (SPLM) which coupled elasticity, plasticity and damage was proposed [13]. The improved SPLM could simulate concrete structures accurately and efficiently with less computational effort. Gu et al. [14] presented an effective approach to integrate peridynamics with an open source software. This integrated method can capture the complex process of crack propagation of elastoplastic material by combining the state-based peridynamic program with 3D plastic constitutive models. Statebased peridynamics was also adopted to study the thermal shock problems in ceramics, including crack propagation [15]. However, in the case of both quasi-static and dynamic loading, the non-ordinary statebased peridynamic method reveals numerical oscillations, especially when severe deformation gradients were present due to cracks [16]. The zero-energy mode is also associated with the non-ordinary statebased peridynamic method, which would lead to the inaccurate calculation of the deformation gradient.
The bond-based peridynamic method is proved to be stable when simulating the crack-growth problems. While the bond-based peridynamic formulation is generally more suitable for brittle fracture analysis, it has also been successfully used to study dynamic fracture in brittle materials [17]. Gerstle et al. [18] developed a new model using bond-based peridynamics by adding pairwise peridynamic moments to simulate linear elastic materials with a varying Poisson's ratios. The new model is called the "micropolar peridynamic model" and can simulate damage and cracking in concrete structures using an implicit rather than an explicit algorithm. This model can efficiently explain a number of different microcracking (damage) and fracture mechanisms observed in concrete. The peridynamic theory has also been applied successfully to simulate impact damage and failure analysis in composite laminates [19,20]. Modeling of fracture problems related to crack growth in brittle materials represents an active ongoing challenge in computational mechanics [21]. It is required that the computation algorithm has the capacity to deal with the fracture pattern and crack-propagation paths in the ballistic application to produce accurate results [11]. Ha et al. [3] reproduced various characteristics of dynamic brittle fractures like crack-branching, crack-path instability, and asymmetries of crack paths, using the bond-based peridynamic model. In many applications, it is necessary to discrete the computational model with changing horizon sizes. Spurious wave reflections occur in original peridynamics when use different horizon sizes [22]. A new peridynamic approach called dual-horizon peridynamics (DH-PD) was developed to solve the issue of varying horizons and ghost force [22,23]. The capabilities of DH-PD to deal with the issue of spurious wave reflection, multiple material problems and crack stability problems on random particle arrangement were proved [22]. Later, a dual-support smoothed particle hydrodynamics (DS-SPH) inspired by the DH-PD is developed [24]. The DS-SPH can compute the unbalanced interactions between the particles with different smoothing strengths correctly. Researchers also investigated the ballistic problems by combining the peridynamics with other numerical methods, and good progress has been made so far. Houfu [25] accomplished coupling the peridynamic model with a modified SPH model to simulate the fragmentation for buried explosive loads, and the simulation results generally agree with the experimental data. Bo et al. [26] proposed an explicit dynamics implementation of the bond-based peridynamic formulation to simulate the dynamics of the fracture-process via the Discontinuous Galerkin method, using the classic peridynamic governing equation. This model was capable to capture the 3D dynamic crack process in brittle materials effectively. Sun et al. [27] presented an approach which coupled the peridynamic theory with the numerical substructure method (NSM) to model structures with local discontinuities. This approach takes advantage of the strong capacity of the peridynamics for dealing with discontinuities and high computational efficiency of the NEM. Giannakeas et al. [28] proposed a numerical model by combining the finite element method and peridynamics to simulate cracks in ceramics after severe thermal shock. A mesh-free method known as peridynamic differential operator (PDDO), of which the concept comes from peridynamic theory, was adopted to investigate incompressible inviscid fluid flow with moving boundaries [29]. There are many attempts to overcome the limitation of Poisson's ratio for original bond-based peridynamics [30,31]. For example, Wang et al. [30] proposed a numerical method that takes into account the interacting forces between two material points along the horizon, which are not only related to bond stretching but also to the rotation of the conjugated bond angles. The author also developed a coupled thermo-mechanical bond-based peridynamic model to simulate thermal fracturing in ceramic nuclear systems [32]. Furthermore, fatigue is a major cause of failure in many material structures and prediction of fatigue crack propagation is a challenging problem. Bazazzadeh et al. [33] adopted bond-based peridynamics to simulate fatigue crack propagation. Yang et al. [34] proposed a new damagemodel based on the bond-based peridynamic formulation to investigate crack propagation in concrete. Cheng et al. [4] proposed a bond-based peridynamic model to analyze the dynamic fracture of shale material. The bond-based peridynamic model was also introduced to simulate dynamic brittle fracture in functionally graded materials [35]. It can be seen that bond-based peridynamics has the capacity to capture the dynamic fracture of brittle Materials accurately.
The capacity of the numerical method to simulate impact problems of brittle materials also depends heavily on the constitutive models. Suitable constitutive relations are needed to describe the fracture and damage properties of brittle materials accurately. In the last two decades, a few constitutive relations have been proposed to simulate and predict the fracture characteristic of brittle materials. Fahrenthold [36] provided a constitutive model that includes a continuum damage mechanics description to simulate the impact of brittle materials. The model was used to estimate the effects of both damage accumulation and brittle fracture on the penetration of ceramic by simulating the penetration of spherical projectiles of a steel plate. Rajendran et al. [37] proposed a constitutive model based on an internal state variable to describe the shock and high strain rate properties of some ceramics under impact loading conditions. Simha [38] developed a computation model based on the bar impact and plate impact data for a penetration-response test of AD-99.5 alumina. Deshpande et al. [39] presented a mechanism-based constitutive model for the inelastic deformation and fracture of ceramics, considering the effects of microstructural parameters. One of the most widely used constitutive models to simulate the fracture response of brittle materials under impact loads was proposed by Johnson et al. [10]. This model, called JH-2 [10], considers strength softening, damage accumulation and rate-dependence of ceramic materials, and was fairly accurate in predicting the dynamic behavior of ceramics. In the classical JH-2 model, ceramics are treated as elastic materials before damage occurs, and the material begins to soften when the damage begins to accumulate and they are regarded as intact material whose strength varies with the cumulative damage. Because the bond-based peridynamic method can naturally solve the discontinuity associated with cracks, the aim of this work is to develop a new bond-based peridynamic model that takes into account strength softening, damage accumulation and the rate-dependence, which were described in classical JH-2 model, of brittle materials for the fracture analysis.
In this paper, we first construct a rate-dependent bond-based peridynamic formulation for ceramic materials. The Prototype Microelastic Brittle (PMB) model, which was introduced in bond-based peridynamic theory by Silling [40], assumes that the bond between material points breaks, when stretched beyond a predefined elastic limit. In this work, the PMB model is modified to take into account the unique pseudo-plastic properties of each bond under compression. In this model, the different responses under tension and compression of ceramics are considered, and tensile damage and compressive damage are shown respectively. Finally, numerical examples are presented to verify the implementation with the analytical solutions. The proposed rate-dependent bond-based peridynamic formulation is used to simulate the crack evolution and damage process of the ceramics. This paper is organized as follows: in Section 2, we briefly review the mathematical formulation of bond-based peridynamics, especially the PMB model. In Section 3, the dynamic responses during each stage of the brittle materials are described in detail. Furthermore, the rate-dependent bond-based peridynamic model is introduced to describe the dynamic response of ceramics under impact. Then, the numerical techniques that include the constitutive update algorithm and contact algorithm are given in Section 4. In Section 5, the numerical implementation and examples are provided to verify the proposed peridynamic model for ceramics, and the results are compared with the experiment. Finally, the conclusions are summarized in Section 6.

Formulation of Bond-Based Peridynamics
Bond-based peridynamics is a nonlocal numerical method, which reformulates the basic equations of motion in integral form, instead of a differential formulation. Thus, it can be used to analyze the dynamic mechanical behavior of brittle materials, which is inherently discontinuous. In bond-based peridynamic theory, a continuum solid Ω is separated by a set of individual material points x i with associated mass density ρ i and volume V i , where i = 1, 2, …, n is the index of each material particle. A material particle x i only interacts with particles x j within a finite distance δ, which makes up a region called Horizon H x i of x i , see Fig. 1. Interaction between material points is described directly by a pairwise force function f of the bond, which contains all of the constitutive information associated with the material. The equation of motion for a material point x i at time t in the reference configuration is given as denote the displacement vector field and accelerating vector field of the material point x i at time t respectively, and b[x i , t] is the body force density. V j denotes the volume of material point x j within the horizon of x i , and the pairwise bond force vector is a function of the relative positive vector ξ in the reference configuration and the relative displacement vector η in the current configuration, which are expressed, respectively, as Here, ξ + η is the current relative position vector of the material point x i and x j , and, according to bondbased peridynamic theory, for the PMB model, the pairwise force function f(ξ, η) can be expressed as Here, ξ represents the original bond-length in the reference configuration, and η the current length of the bond after deformation. Eq. (3) implies that the pairwise force vector is parallel to the current relative position vector. c(ξ, η) is the micro-elastic modulus function and represents the stiffness of the pair-wise bond, which can be expressed as c(ξ, δ) = c(0, δ)g(ξ, δ). The accuracy of the simulation is highly dependent on the micromodulus function c(ξ, δ), which contains all information about the material. In the original PMB model, the function c(ξ, δ) is reduced to a constant c(0, δ) = 18 K/πδ 4 for 3D case, in which K is the bulk modulus of the material, by simplifying the kernel function g(ξ, δ) to [17] gðn; dÞ Moreover, s(η, ξ) is the scalar bond stretch, which is defined as For the PMB model, it is assumed that the bond breaks when the corresponding stretch s exceeds the critical stretch s 0 and no longer sustains a force. The linear damage model is shown in Fig. 2. The critical stretch s 0 of the PMB material is determined by the fracture energy of material G 0 . It is given by The damage at a material point is defined as [40] 'ðx; tÞ Here, μ(x, t, ξ) is a history-dependent scalar function that takes on values of either 1 or 0, and expressed as It is worth noting that 0 ≤ μ ≤ 1, with 0 representing the intact material, while 1 represents complete disconnection of a point with all of the points within its horizon [40].

Basic Formulation of the Rate-Dependent Peridynamic Model
In this Section, we first summarize the damage and fracture process for ceramic materials under dynamic loading. Then, the rate-dependent bond-based peridynamic model for ceramics is proposed to describe the dynamic response.

Dynamic Response of Ceramics During Ballistic Penetration
Composite armor consists mainly of a ceramic front plate and a ductile back plate such as steel, aluminum nitride, or fiber reinforced composite. The ceramic tile plays an important role in making armor bulletproof. The reason why a ceramic tile can be bulletproof is because it can absorb most of the kinetic energy of an armor-piercing projectile.
The energy absorption mechanism can be divided into three stages: Firstly, the initial stage, when an armor-piercing projectile hits a ceramic tile, it produces a strong compression wave on the projectile/ ceramic interface area, which generates compressive stress both inside the projectile and the ceramic plate. The compressive stress increases quickly due to the high hardness and compressive strength of ceramics. Plastic deformation occurs in ceramic tile when the compressive stress exceeds the yield strength. At this time, initial conical and radial cracking forms in the ceramic tile, which absorbs much energy. Next, during the erosion stage, the blunt or broken bullet-projectile continues to penetrate into the ceramic tile and is eroded by the ceramic target simultaneously. Cracks in the ceramic tile begin to propagate and a successive fragmentation layer forms, which will further absorb the kinetic energy of the projectile. Finally, during the comminuted conoid stage, with the penetration of the projectile, the cracks in the ceramic tile continue to propagate and form a ceramic comminuted conoid, the apex of which is the tip of the projectile. The comminuted conoid pushes on the ductile back plate and distributes the impact load over a large area of the ductile back plate, provided the ceramic tile is supported by a ductile panel. The associated fracture pattern is illustrated in Fig. 3.
Considering the process described above, it clear that both the initiation and evolution of the cracks in the ceramics play an important role during the penetration. Ceramics are inherently brittle and have very high compressive strength, especially under confinement and a high loading-rate, but they have low tensile strength [41]. The failure mechanism of brittle materials is significantly different under quasi-static conditions and dynamic loads. In fact, the fracture strength of the brittle materials is highly rate-sensitive, which increases with increasing strain-rate under compression. Under quasi-static loads, the brittle material is fractured, or crushed due to micro-cracks, which propagate to form macro-cracks, and almost no plastic deformation occurs. However, this process represents a unique inelastic response characteristic that is different from the plastic flow of ductile materials. The unique inelastic behavior of ceramics under compression can be explained by the micro-crack formation and crystal plasticity [39,42].
To capture this process accurately, we summarize the damage features according to the damage patterns for different loadings and stages, based on existing research as follows: The penetration mechanism under impact loading of brittle materials is complex and can be generally divided into three successive loading stages [43].The three successive stages described above are shown in Fig. 4. The first stage, the elastic Figure 3: Schematic of the damage for a ceramic target, during penetration by a projectile stage, is marked as OB in Fig. 4. The shock waves, which are compressive waves initially, propagate through the target, and dynamic compressive loading is generated simultaneously. The pressure increases in the ceramic and may exceed the Hugoniot Elastic Limit (HEL), point B in Fig. 4, which defines the second stage, the damage stage (see BC in Fig. 4). Important pseudo-plastic behavior, which is different from ductile materials, can be observed during this period [44]. Then, the material begins to soften with increasing plastic strain, which allows for gradual decreasing in strength. The shock waves are reflected and spread, along the thickness of the ceramic plate, to the projectile. The waves become tensile waves upon reaching the rear surface of the ceramic target and radial cracks occur. Dynamic tensile loads, caused by the radial cracks, produce extensive fragmentation. A fractured (cone shaped) zone is generated due to the interaction of the radial cracks with the conical cracks. Different from the compressive damage, the material fails immediately because the tensile elastic limit is exceeded and cannot withstand any tensile loads (only compressive loads). The last stage, the fractured stage, is shown in the CD Section in Fig. 4. The projectile penetrates the damaged ceramic, which fails completely during this period. The fractured ceramic can withstand compressive loads with constant strength.
To capture the damage process at each stage accurately, it is critical to develop an appropriate and effective numerical model for brittle materials. Several constitutive models were proposed to describe the damage characteristics and fracture pattern of brittle materials-see Section 1. In this work, the properties described above, which were considered in the Johnson-Holmquist (JH-2) model, are taken into account to develop a rate-dependent bond-based peridynamic model for ceramics. This model considers that the material begins to soften when the damage starts to accumulate. The proposed model is developed to capture the pseudo-plastic properties of a ceramic under compression.

Rate-Dependent Peridynamic Model for Ceramics
The original linear damage model called PMB model in bond-based peridynamic theory is suitable for brittle materials describing the materials characterized brittle. Ceramic materials are inherently brittle, and it is known to have a brittle response under tensile loading, but they have an intrinsic strain-softening behavior under compression. Ceramics can have significant strength after failure under pressure, but the compressive strength of damaged ceramic is lower than that of the intact ceramic. This characteristic of ceramic materials under compression can be described as plastic softening behavior. During projectile impact a ceramic material experiences both tensile loading and compression [42]. Therefore, the plasticity of ceramics such as strength softening and damage accumulation should be taken into account. These mechanical properties, which were introduced above, are considered in the bond-based peridynamic constitutive model (see Fig. 5). Ceramics are generally brittle materials and exhibit brittle behavior under tensile loading. In other words, the material fails, when the tensile elastic limit is reached-see Section OA in  Fig. 5. Because the effect of Poisson's ratio is not significant for problems that involve dynamic cracks, the response of the materials is mainly governed by the fracture mechanism here [34], and the correction of Poisson's ratio is not considered in this model. The constitutive relations of the bond, under tensile and compressive loading, are given, individually, as follows:

Tensile Loading
When the bond is stretched, the constitutive relation describing the bonds follows the original PMB model: fðg; nÞ ¼ f ðg; nÞ n þ g n þ g k k (9) f ðg; nÞ ¼ csðg; nÞ 0 sðg; nÞ s 0 0 sðg; nÞ > s 0 (10) Here, f(η, ξ) represents the scalar-valued pairwise force, c and s 0 are parameters that relate to the material introduced in Section 2. The form of damage, when the bond stretches, is considered in the PMB model-see Eqs. (7) and (8).

Compressive Loading
When under compressive load, the form of the constitutive relation refers to the JH-2 model, and both material softening and rate-dependence are considered. The strength of the bond between two material points is considered via where D is the accumulate damage of the bond under compression, which is given by in which Δs p is the relative plastic deformation increment during a time step, and s p ¼ P Ds p represents the accumulate plastic deformation. s 1 is the critical fracture deformation of the bond and s e is the elastic where p 0 is the static strength without considering rate-dependence, and _ s is the deformation rate of the bond without considering rate-dependence. Furthermore, β and γ are material constants. A method to identify these parameters is shown in Section 3.3. Using Eq. (5), the _ s is given by [20] Because the pseudo-plastic behavior of the bond is considered under compression, the pairwise force is given by

Yield Criterion
Note that the force of the bond can never be higher than the current yield strength. However, it can be lower if unloading occurs. Furthermore, the function is defined to evaluate whether the bond force reaches the maximum value. This is called the yield condition-referring to the continuum mechanics. It is given by Using Eq. (17), it can be formulated as follows 'ðf ; s p ; _ sÞ ¼ f À 1 À ð1 À bÞs p s 1 À s e p i (18) Referring to the consistency condition, we can write Here, the second derivative term of s is neglected, implying that the bond force must remain at the yield strength value to consider any decrease due to softening. Using the time derivative of Eq. (16) as _ f ¼ cð_ s À _ s p Þ, the plastic deformation rate can be found as where H ¼ ð1 À bÞp i s 1 À s e .

Determination of the Parameters
In this Section, we determine the parameters p 0 , β and γ, which relate to the pseudo-plastic properties proposed in Section 3.2.2, are introduced. It is assumed that a homogeneous sphere, whose radius equals to the horizon δ, is subjected to isotropic pressure p (see Fig. 6). During the elastic stage, the bond-force is f = cs, which satisfies Then, we can obtain the relation for bond-relative deformation s and the pressure P In the original JH-2 model, intact strength and fracture strength of the materials are given, respectively, as where A, B, N, M, σ HEL and T are material constants (see [10]). For this condition, we assumed previously that relative deformation of the bond equals to the compressive strain during the elastic stage. Hence, the elastic compressive limit s e can be derived as Here, E is Young's modulus. Combining Eqs. (22)-(25), the elastic compressive limit s e can be obtained, when reaching the elastic limit. Consequently the compressive strength is For the proposed model, the reduction factor β, due to the fracture of the material can be calculated using Because the relative bond deformation is similar to the strain, the deformation-rate factor γ is equal to the original JH-2 model.

Solving Method
In order to better implement numerical formulations in the program introduced in Section 3, the computational body is normally discretized into material particles. After the spatial discretization, the motion equation, Eq. (1), can be replaced by a discretized form: where f is the pair-wise bond force, n is the time-step number, and the subscripts denote the node number which can be expressed as u n i ¼ u x i ; t n ð Þ. V p is the volume of node p, which represents the nodes in the horizon of node i. An explicit central difference-technique is used for the acceleration: Thus, the displacement of material point i at the next time step (n + 1) can be obtained by When the bond stretches, the pair-wise force f can be calculated using Eq. (10). When the bond shortens, the plasticity of the bond is included, in which case the bond-force updating algorithm is needed to calculate the compressive bond force. The bond-force updating algorithm will be introduced in Section 4.3. The flow chart for the calculation is shown in Fig. 7. This flow chart helps clarify the numerical implementation of the proposed model.

Kernel Function
Since the kernel function g(ξ, δ), mentioned in Section 2, describes the effect of a spatial distribution of the materials on the bond force, the form for the kernel function g(ξ, δ), which is used in this work was proposed by Huang [17]:

Bond-Force Updating Algorithm
In this part, the bond-force updating algorithm is given to solve the numerical formulations, which takes into account the pseudo-plastic behavior, when the bond under compression.
To ensure consistency with the yield condition for each incremental time step, the incremental nature of the computational process is taken into account in this algorithm. Such a procedure is similar to the returnmapping algorithm. To calculate the real bond force at each stage, we first assume the deformation for a time step is purely elastic, which can be formulated as: where the subscript n and n + 1 also represent the time step number at times t and t + 1, respectively. l n+1 is the relative distance of the material points connected by a bond during time increment Δt, and l p n is the permanent length with existing plastic deformation. The force of the bond at time t + 1 is given by Eq. (33). The superscript trial is given because of the assumption of pure elasticity. If s n+1 > 0, the bond tends to become tensile and no plastic deformation occurs. As a result, the actual force is equal to the trial force, i.e., f t ¼ f trial tþ1 . Otherwise the force should be updated as plastic deformation occurs: Whether or not plastic deformation occurs depends on substituting the trial bond force into the yield criterion-as given in Eq. (18). If 'ðf trial tþ1 ; s p ; _ sÞ < 0, the deformation of the bond is elastic. Otherwise, we need to calculate the real elastic-deformation increment, the plastic-deformation increment, and the true force of the bond. Referring to the return-mapping algorithm, the current force of the bond is computed as below.
Using Eq. (20), the plastic stretch increment Δs p during time increment Δt, can be obtained by Then, the accumulated plastic damage is given by Now, we can calculate the current strength associated with the damage The force of the bond is updated according to Eqs. (37) and (38) until 'ðf ðkÞ tþ1 ; s p ; _ sÞ < 0 is satisfied, and the real bond force is given Eq. (39) A flow diagram for bond-force updating algorithm is shown in Fig. 8.

Contact Algorithm
For ballistic applications, impact and penetration are the most common dynamic processes [11]. The contact algorithm used in this work was introduced in [45]. There is initially an interpenetration of material points after contacting occurs, between the projectile and the target material. The projectile is assumed to be rigid and is not deformable at any distance, when moving at its own velocity, and the target material is deformable and governed by the peridynamic equation of motion. The material points inside the projectile are relocated to their new positions outside the projectile. The new locations are selected to satisfy the condition for the closest distance to the surface of the projectile. This process is illustrated in Fig. 9, and will create a contact surface between the projectile and material points at the target at time t.
The velocity of such a material point x i , in its new location at the next time step t+Δt, can be calculated as where u tþDt i is the modified displacement vector at time t+Δt, and u t i is the displacement vector at time t, and Δt is the time increment. At time t+Δt, the contribution of the material point x i to the reaction force, from the target material to the projectile, F tþDt i , can be computed as where v tþDt i is the velocity vector at the time before relocating the material point x i . Furthermore, ρ i and V i represent density and volume, respectively. Summation of the contributions of all material points inside the projectile produces the total reaction-force on the projectile at time t + Δt, which can be expressed as

Numerical Simulations
In this Section, the damage process of the ceramic under different impact loadings is investigated using the proposed rate-dependent bond-based peridynamic model for brittle material. Firstly, the sensitivity of the results to variations of the main discretization parameters is analyzed. Then, a numerical example is presented to validate the feasibility and accuracy of the proposed model for brittle material. Simulations have been conducted to demonstrate the capability of the proposed rate-dependent bond-based peridynamic model. The damage characteristic of the ceramics under different impact velocities is provided to analyze the fracture mechanism of ceramics in more detail. Finally, the numerical result, using the proposed bond-based peridynamic model, is compared with the experiments. The results of the proposed bond-based peridynamic model are in better agreement with the experiment than for the nonordinary state-based peridynamic method.

Analysis on Convergence
In this section, a drop-ball test using the proposed rate-dependent peridynamic model is simulated to analyze two types of convergence: the m-convergence and the δ-convergence. The configuration of the Al 2 O 3 ceramic plate is 10 × 10 × 1 mm 3 , see Fig. 10, whose density ρ and Young's modulus E are 3740 kg/m 3 and 310 GPa respectively. The drop ball is steel and assumed to be rigid with a diameter of 2.5 mm and a density of 7800 kg/m 3 . Initial velocity of 300 m/s is assigned to the drop ball. In the m-convergence the horizon δ is a constant, and take different values of m = δ/Δx. In this example, we consider the horizon δ = 2.0 × 10 -4 m fixed, and three different values of m and Δx are used as m = 2, 3, 4 and Δx = 1.0 × 10 -4 , 6.67 × 10 -5 , 5.0 × 10 -5 m respectively, see Fig. 11. The crack pattern on rear surface of ceramic plate is shown in Fig. 12. For a fixed horizon, when m increases, the crack path is more obvious and more consistent with the phenomenon described in literature [46]. Since the cone Figure 9: Relocation of material points inside a target material to represent contact with the impact [45] cracks and radial cracks can be captured when m = 3 or 4, m = 3 is employed for the rest examples in terms of the computational efficiency.

Verification
The validation example studied here is presented by Johnson et al. [10] in the JH-2 model. A cube with sides of 1.0 m is considered, which is characterized as brittle, and displacement in three directions of the bottom face is constrained, while the four faces on the side are constrained with respect to their normal directions-see Figs. 15 and 16. A normal displacement is slowly applied to the top face, until it reaches 0.05 m. Then it is slowly released until a zero stress-state is reached. In the peridynamic implementation in this work, the cube is discretized into 5 × 5 × 5 particles, as given in Fig. 17. The displacement boundary condition is imposed by prescribing constraints through the material points on the displacement field of a "fictitious material layer" R c according to [45]: The time step Δt = 10 -6 s was chosen. For case C, the material was defined as having fractured strength and allowed to accumulate plastic strain. Furthermore, the comparison of the damage result between this work and the original JH-2 model is shown for a force vs. displacement, as displayed in Fig. 18. The response is complex. From points 1 to 2 the material loads elastically until the intact strength is encountered at point 2. From points 2 to 3 the material flows plastically from the intact strength at point 2 to the fractured strength at point 3. At the same time, the damage goes from D ¼ 0 at point 2 to D ¼ 1:0 at point 3. From points 3 to 4 the material continues to flow plastically. The loading direction is reversed at point 4 and material unloads elastically from points 4 to 5. From points 5 to 6 the elastic unloading continues and unloads plastically from points 6 to 7. The simulation is generally consistent with the solution in the JH-2 model. Since case C includes the entire process of damage evolution for brittle materials, considering the rate sensitivity and strain softening, the consistency of the results   confirms the validity of the proposed rate-dependent bond-based peridynamic model and its ability to capture the complete process of damage during loading and unloading. In addition, there is still a slight discrepancy between simulation and theory. This may be for two reasons: One reason could be that it is due to the fixed Poisson ratio in the bond-based peridynamic model. Because the structure is under quasi-static loading in this example, the effect of Poisson's ratio is clearer. The second reason may be that the identification of the parameters in Section 3.3 is relatively simple and requires further optimization in the future.

Drop-Ball Test
In this section, a drop-ball test is modeled using the proposed rate-dependent peridynamic model to investigate the dynamic response of the ceramic plate under impact loading. Firstly, the crack evolution of the Al 2 O 3 ceramic plate under impact loading with time is derived to analyze the damage and fracture pattern of ceramics. Then, different damage patterns under different impact velocities are investigated.

Initiation and Propagation of Cracks
As shown in Fig. 10, the configuration of the Al 2 O 3 ceramic plate is 10 × 10 × 1 mm 3 , whose density ρ and bulk modulus K are 3700 kg/m 3 and 231 GPa respectively. The drop ball is steel and assumed to be rigid with a diameter of 2.5 mm and a density of 7800 kg/m 3 . Impact velocity of the drop ball is assigned to be 300 m/s. In the following example, the target plate is discretized into 200 × 200 × 20 material points, for a particle spacing of Δx = 0.05 mm. The peridynamic horizon is chosen proportionate to the grid spacing (i.e., δ = mΔx, for some m ≥ 1). In these example, the horizon size of the material point is chosen to be three times the particle spacing and therefore m = δ/Δx = 3. The time-step size used in this example is Δt = 1.0 × 10 -9 s. The time-step size is chosen to ensure that details are being captured during the simulation without loss of efficiency, considering the minimum discretization size of the system and the maximum velocity of the particles. We also take the effect of the material intensity into account, using the kernel function in Eq. (31). The time evolution of conical and radial cracking at the front surface, rear surface, and cross-section of the ceramic tiles is shown in Tabs. 1-3. The particle distribution under impact loading is shown in Tab. 4.
Using the crack changes as a function of time, it can be seen that the conical cracks form at the front face firstly, due to the propagation of the shock waves. At the front face of the tile, no radial cracks formed until t = 4 μs, which is caused by the tensile waves that were reflected from the rear surface. There is almost no compressive damage on the back of the tile. The radial cracks, which are caused by tensile stress, are     initiated at the rear face of the ceramic tile and propagated from the center toward the edges [47]. Based on the damage feature at the front and rear surfaces shown in Tabs. 1 and 2, the cracks on the rear surface appear later than those on the front surface. This is because it takes time for the shock waves to propagate from the impact point to the rear surface of the tile. When the shock waves propagate to the rear surface, radial cracks appear due to tensile stress. From the damage distribution of the cross-section shown in Tab. 3, fractured cone forming in the ceramic tile under impact loading can be observed clearly. In this work, we can obtain tensile damage and compressive damage separately, which facilitates a more detailed analysis of the damage and fracture mechanism for ceramic materials. The initiation and evolution of cracks in the tiles can be captured using the proposed rate-dependent bond-based peridynamic model.

Damage Pattern under Different Impact Velocities
In this part, the different damage features under different impact velocities are investigated to capture the characteristics of rate sensitivity, which is considered in the proposed model based on the fracture mechanism of brittle materials. The development of conical and radial cracks, as a function of impact velocity, at the rear face of the ceramic target is shown in Tabs. 5 and 6. It is compared to the experimental results given in [46]. By comparing, it can be seen that the proposed bond-based peridynamic model in this work can capture the cone and radial cracks of ceramics, for low to high impact velocities, well. In [46], three numerical models were used to simulate the evaluation of damage and compared with experimental results. A comparison of radial and conical cracking between the simulation and experiment is shown in Tab. 5. It can be seen that the crack evolution, which is observed in the experiment, is captured accurately using the proposed ratedependent bond-based peridynamic model from low to high impact-velocities. The radial cracks form from the periphery of the cone crack and extend outwards. The damage results, using the proposed ratedependent peridynamic model, suggest that the number of radial cracks in the plate increases with increasing impact velocity. This tendency was observed in the experiments. The radial cracks are more pronounced for higher impact velocities. The damage results for the cross-section are given in Tab. 6. It can be seen that a fragment of the ceramic tile can be captured by the proposed rate-dependent PD model in this work.

Edge-On Impact
The Edge-On impact (EOI) technique can observe the dynamic fracture of brittle materials directly. It is convenient to compare the numerical simulation with the results from EOI results [11]. In the EOI test, a  projectile impacts the sample edge, as illustrated in Fig. 19, and damage propagation is observed with a highspeed camera. In this section, the experimental set-up is numerically simulated using the proposed bondbased peridynamic model. The simulation result is compared with the experimental result and the numerical result, which was calculated using the non-ordinary state-based peridynamic model [11]. The properties of the sample and the projectile are as described in [11]. The configuration of the target is 100 × 100 × 10 mm 3 . The projectile is a solid cylinder and assumed to be rigid with a mass density of 8060 kg/m 3 . The numerical results are shown in Fig. 20. By comparing the damage results between the simulation and experimental results, we can see that the result using the proposed bond-based peridynamic model agree better with the experiments than the non-ordinary state-based peridynamic model.  Table 6: Comparison of damage results for the cross-section of the tile given in [47] under different impact velocities Figure 19: Configuration of edge-on impact

Conclusion
A bond-based PD model that considers strength-softening and rate-dependence of ceramic materials is proposed in this work. A bond-force updating algorithm is introduced to solve the pseudo-plastic feature of ceramic materials in the numerical model. This rate-dependent bond-based PD model was used to simulate the complete damage-and fracture-process under impact loading for ceramic materials. Using this model, it is possible to capture, accurately, both the brittleness (tensile loading) and pseudo-plastic behavior (compressive loading). The damage features under different impact velocities were simulated and compared with the experiment and other numerical models. The results, using the proposed model in this work, are in good agreement with the experiments. The Edge-On impact test of glass was simulated in this work and compared with the experimental result. Based on our simulations, the following conclusions can be drawn: 1. The proposed bond-based PD model can capture the entire damage process, from loading to unloading, including the fracture feature of ceramic materials. This was validated by implementing the example given in the original JH-2 model. The results of the proposed model are consistent with the results giving in the JH-2 model. 2. This model can accurately capture the dynamic response during ballistic impact. The entire process, including the initiation and evolution of the cracks, can be observed clearly with the proposed model. Both tensile-and compressive-damage are described in this model, which enables a deep analysis of the damage and fracture mechanism in ceramic tiles. 3. This model can capture the damage and fracture process of ceramic materials under impact from low-to high-velocity. The results are in good agreement with the experiment. The simulation results also show that the splash of the ceramics can be observed clearly. 4. The Edge-On impact test was simulated using the proposed bond-based peridynamic model in this work.
A comparison between the simulation and experiments indicates that the bond-based peridynamic model is in better agreement with the experimental results than the non-ordinary state-based peridynamic model.
In general, the rate-dependent bond-based peridynamic model, which is proposed in this work, can be used to simulate the damage and fracture behavior of ceramic materials.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.