Next Article in Journal
Manufacture and Vibration-Damping Effect of Composites for Archery Carbon Fiber-Reinforced Polymer Limb with Glass Fiber-Reinforced Polymer Stabilizer
Next Article in Special Issue
Modeling and Simulation of the Hysteretic Behavior of Concrete under Cyclic Tension–Compression Using the Smeared Crack Approach
Previous Article in Journal
The Application of a Fluoride-and-Vitamin D Solution to Deciduous Teeth Promotes Formation of Persistent Mineral Crystals: A Morphological Ex-Vivo Study
Previous Article in Special Issue
Numerical Simulation of Failure Behavior of Reinforced Concrete Shear Walls by a Micropolar Peridynamic Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model for Impact Problems of Quasi-Brittle Materials

1
Hubei Key Laboratory of Theory and Application of Advanced Materials Mechanics, Wuhan University of Technology, Wuhan 430070, China
2
Department of Engineering Structure and Mechanics, Wuhan University of Technology, Wuhan 430070, China
*
Authors to whom correspondence should be addressed.
Materials 2023, 16(11), 4050; https://doi.org/10.3390/ma16114050
Submission received: 2 April 2023 / Revised: 13 May 2023 / Accepted: 26 May 2023 / Published: 29 May 2023

Abstract

:
In this work, we have developed a novel bond-associated non-ordinary state-based peridynamic (BA-NOSB PD) model for the numerical modeling and prediction of the impact response and fracture damage of quasi-brittle materials. First, the improved Johnson-Holmquist (JH2) constitutive relationship is implemented in the framework of BA-NOSB PD theory to describe the nonlinear material response, which also helps to eliminate the zero-energy mode. Afterwards, the volumetric strain in the equation of state is redefined by the introduction of the bond-associated deformation gradient, which can effectively improve the stability and accuracy of the material model. Then, a new general bond-breaking criterion is proposed in the BA-NOSB PD model, which is capable of covering various failure modes of quasi-brittle materials, including the tensile-shear failure that is not commonly considered in the literature. Subsequently, a practical bond-breaking strategy and its computational implementation are presented and discussed by means of energy convergence. Finally, the proposed model is verified by two benchmark numerical examples and demonstrated by the numerical simulation of edge-on impact and normal impact experiments on ceramics. The comparison between our results and references shows good capability and stability for impact problems of quasi-brittle materials. Numerical oscillations and unphysical deformation modes are effectively eliminated, showing strong robustness and bright prospects for relevant applications.

1. Introduction

In recent decades, ceramic materials have attracted great attention from researchers and engineers and have been extensively adopted for applications in the fields of aerospace, military equipment, biomedicine, and industrial manufacturing, due to its excellent performance, such as high hardness, high specific strength, high specific modulus, high temperature resistance and high abrasion resistance. In the initial stage, ceramics are applied in the field of armor protection with great success due to their excellent anti-penetration properties [1], which has great prospects. In military and protection applications, ceramics are often subjected to extreme loading conditions [2], such as shock waves from explosions or projectile penetration. Therefore, extensive research and testing have been devoted to the understanding of the fundamental nature and mechanism of the impact response and anti-penetration behavior of ceramic materials. Various approaches and numerical simulations have been carried out to study the damage and fracture of ceramics [3,4,5], including but not limited to the finite element method (FEM), cohesive zone model (CZM), and extended finite element method (XFEM). However, due to the mesh dependence, the damage or fracture of materials leads to grid distortion, and the calculation cannot continue in mesh-based approaches. On the other hand, the cohesive zone model only allows cracks to propagate along element boundaries, which may not be able to accurately describe the crack path. In XFEM, although cracks can propagate along any direction, the introduction of additional degrees of freedom brings considerable computational complexity. Damage or fracture will result in discontinuity of the displacement field, which creates inherent limitations for the above numerical simulation methods based on continuum mechanics. To date, it is still challenging for numerical approaches to give a precise prediction of the complex fracture and damage patterns on ceramics.
To naturally simulate the discontinuous behavior of materials such as damage and cracking, Silling [6] proposed the peridynamic (PD) theory with nonlocal interactions as the core. This theory is also called bond-based PD (BB PD). In the PD theory, a body is described by a set of discretized material points with mass and other physical properties. The description of macroscopic material properties is formed by nonlocal bonding and interaction between material points. The damage of materials is part of PD theory, which is described by breaking bonds and truncating interactions between material points. PD is a theory of nonlocal continuum mechanics, which uses integral equations instead of integral–differential equations and does not require the displacement field to satisfy the requirement of continuous derivability. Therefore, the PD theory is well suited to modeling discontinuity problems without the difficulties encountered by classical local theory. Further, Silling [7] successfully proved the ability of PD theory to model complex fracture problems using the Kalthoff–Winkler experiment. In addition, Silling and Askari [8] developed the PD prototype microelastic brittle (PMB) material model and its numerical implementation technology, and derived the critical stretch formulation. The BB PD theory proposed earlier does not include stress and strain measurements, encountered the restriction of fixed Poisson’s ratio, and cannot describe the shear deformation of materials, which had certain limitations. Hence, Silling et al. [9] proposed the state-based PD (SB PD) theory, which extends the PD theory and the material responses it can describe. The theory can be divided into the ordinary state-based PD (OSB PD) and the non-ordinary state-based PD (NOSB PD), depending on the form of interactions.
The PD theory has significant advantages in modeling discontinuity problems, and is widely used in simulating material damage, cracks, and fractures. Some promising applications of the PD theory in damage prediction and fracture modeling are outlined below. Ha and Bobaru [10,11] reproduced the crack growth and fracture behavior of glass materials in dynamic tensile experiments. Zhu and Ni [12] proposed the bond rotation effect in BB PD, which removed the limitation of Poisson’s ratio. Chu et al. [13] presented a rate-dependent BB PD model to simulate the anti-penetration behavior of ceramics. Liu et al. [14] proposed a comprehensive tensile-shear and compressive-shear failure criterion to evaluate the impact damage and fracture characteristics of ceramics. Erdogan and Oterkus [15] proposed an OSB PD model based on the Mises yield criterion. Zhang and Qiao [16] established a new bond-breaking criterion based on the critical relative rotation angle. Foster et al. [17] implemented an elastic viscoplastic constitutive relationship in NOSB PD and reproduced the Taylor impact test of 6061-T6 aluminum. O’Grady and Foster [18] developed a Euler–Bernoulli beam model. Lai et al. [19,20] implemented the Drucker–Prager and JH2 constitutive models in NOSB PD to model dynamic brittle fracture of quasi-brittle materials. Wu et al. [21] introduced the Holmquist–Johnson–Cook (HJC) constitutive theory to analyze the dynamic mechanical behavior of concrete during impact. Wang et al. [22] studied the thermo-viscoplastic responses of metals under impact loads. Zhu and Zhao [23] simulated single-hole rock blasting with the JH2 damage and tensile failure models. Yang et al. [24] modeled the impact spalling of concrete based on the BA-NOSB PD model for quasi-brittle materials. Li et al. [25] presented a PD model using the Mindlin–Reissner shell theory and simulated the brittle fracture of thin shells.
Although the capability of the NOSB PD model incorporating constitutive models of classical continuum mechanics to simulate complex material responses has been demonstrated, there are still urgent problems to be addressed. The zero-energy mode is one of the major challenges in the NOSB PD theory, which leads to numerical oscillations, affects calculation accuracy, and even produces incorrect calculation results. To date, a number of control schemes and stabilization strategies have been proposed. For example, a stable supplemental force state [26,27,28,29] has been added to the original force vector state to control the zero-energy mode. Yaghoobi and Chorzepa [30] suppressed numerical oscillations by introducing a higher-order version of the conventional deformation gradient. The stress-point method [31] and the higher-order stress-point method [32] were also developed to remove numerical oscillations. In particular, zero-energy mode oscillations are believed to be caused by non-unique mapping relationships from the deformation state to the force vector state, due to the conventional deformation gradient. Therefore, Chen [33] proposed the bond-associated deformation gradient to suppress or eliminate numerical oscillations. The BA-NOSB PD model has also been proved to meet the requirements for stability and convergence [34]. Gu et al. [35] further studied and discussed the possible causes for numerical oscillations in the model. On the other hand, the treatment of broken bonds is also a key problem that deserves attention and needs to be addressed. For the SB PD theory, the force vector state at a material point depends on the deformation behavior of all bonds in its family over time. Therefore, the treatment of broken bonds has a crucial effect on the calculation of deformation states for unbroken bonds. This further affects the stability of the material model, and needs to be taken into account. Silling [27] and Li et al. [36] introduced damage to the influence function so that a broken bond was no longer included in subsequent calculations. In order to make better use of NOSB PD to simulate complex material responses, the material model stability and the treatment of broken bonds will be part of our work.
Following the above introduction, in Section 2, we first briefly review the original PD models for quasi-brittle materials. After that, a modified BA-NOSB PD model is established, and its constitutive update scheme is presented. Then, the other four zero-energy mode control schemes are introduced. Subsequently, the numerical discretization and implementation used in this paper are discussed. In Section 3, the modified model is verified using two benchmark examples, and numerical simulations of edge impact and normal impact experiments on ceramics are carried out. Finally, the conclusions of this work are summarized in Section 4.

2. Methodology

2.1. Brief Review of Peridynamic Models for Quasi-Brittle Materials

To date, a variety of theoretical models [13,14,19,20,21,23,24,37,38,39] have been developed for damage prediction and fracture modeling of quasi-brittle materials. In this section, two PD models for quasi-brittle materials [20,24] are briefly reviewed. First, the original and bond-associated NOSB PD theory and their basic formulas are presented in Section 2.1.1 and Section 2.1.2, respectively. Then, the JH2 material model used for the constitutive update of the above two PD theories is reviewed in Section 2.1.3.

2.1.1. Original Non-Ordinary State-Based Peridynamics

In the PD theory, a body is modeled by a finite number of discretized material points. Further, the material point X interacts only with the material point X within a limited range, and the limited range of influence is called horizon δ . The material point X is called the neighbor of its center point X , and all the neighbors of the material point X form the family H X , i.e., H X = H X , δ = X B :   X X δ . In the reference configuration, the reference position state is defined as
X ¯ ξ = X X = ξ .
In the current configuration, the relative displacement state and the deformation state are defined as
U ¯ ξ = u u = η ,
Y ¯ ξ = y y = ξ + η ,
where y and y are the spatial positions of the two material points, η is a relative displacement between them, and ξ is called a bond.
The NOSB PD theory [9] is a generalization of BB PD, which provides generality in modeling material responses, such as volumetric strain or shear deformation. In NOSB PD, the equation of motion for the material point X at the instant of t can be expressed as
ρ 0 ( X ) u ¨ X , t = H X T ¯ X , t ξ T ¯ X , t ξ d V X + b X , t ,
in which ρ 0 ( X ) represents the mass density of materials, u ¨ X , t represents the acceleration, b X , t represents the external body force density, and T ¯ X , t represents the force vector state. In particular, the force state T ¯ X , t is a function of the deformation state Y ¯ .
In NOSB PD, the material response of a material point will depend on the deformation behavior of all bonds within its horizon. The mapping relationship between the deformation state and the force vector state is the constitutive model of materials in the PD framework. The force density vector at the material point X can be written as
t X , t = T ¯ X , t ξ = ω ¯ ξ P K 1 X ¯ ξ ,
in which ω ¯ ξ is the influence function on a bond.
The shape tensor K is defined as
K = H X ω ¯ ξ X ¯ ξ X ¯ ξ d V X
The first Piola–Kirchhoff stress tensor can be calculated by classical constitutive models, as
P = J σ F T ,
in which J = det F , T is the matrix transposition operator, and σ is the Cauchy stress. The approximate nonlocal deformation gradient tensor F is defined as
F = H X ω ¯ ξ Y ¯ ξ X ¯ ξ d V X K 1 .

2.1.2. Bond-Associated Non-Ordinary State-Based Peridynamics

The zero-energy mode can lead to instability of the material model, which is manifested in numerical oscillations of solutions or completely erroneous calculation results. To eliminate the numerical oscillations caused by the zero-energy mode, Chen [33] proposed the BA-NOSB PD theory and proved the convergence and stability of the PD model [34].
In the original NOSB PD theory, the point-associated deformation gradient results in a non-unique mapping relationship from the deformation vector state to the force vector state. The definition of “bond-associated horizon” is presented to meet the requirement of the unique mapping relationship. The bond-associated horizon is defined as H ξ = H X H X , as shown in Figure 1. In BA-NOSB PD, an individual bond corresponds to a unique deformation gradient tensor. Thus, the bond-associated deformation gradient tensor with respect to the bond ξ at the material point X is defined as
F b = H ξ ω ¯ ξ Y ¯ ξ X ¯ ξ d V X K b 1 ,
where the subscript b indicates that the physical quantity is bond-associated. The bond-associated deformation gradients ensure the desired unique mapping relationship, and the bond-associated shape tensor can be defined as
K b = H ξ ω ¯ ξ X ¯ ξ X ¯ ξ d V X .
The nonlocal deformation gradient tensor F is redefined based on Equation (9) by weighted averaging as
F = n = 1 N b ω ¯ ξ n F b n n = 1 N b ω ¯ ξ n .
where N b represents the number of bond-associated neighbors at the material point X .
The bond-associated force density vector is recalculated as
t b X , t = T b ¯ X , t ξ = H ξ 1 d V X H X 1 d V X ω ¯ ξ P b K b 1 X ¯ ξ ,
The bond-associated first Piola–Kirchhoff stress tensor is defined as
P b = J b σ b F b T ,
in which J b = det F b , and σ b is the bond-associated Cauchy stress tensor.
According to Gu’s work [35], the bond-associated deformation gradient not only satisfies a unique mapping relation from the deformation state to the force vector state, but also satisfies the kinematic constraints of the bonds, as follows:
ξ + η = F b X , ξ ξ ξ + η = F b X , ξ ξ .
Moreover, both the bond-associated shape tensors and the bond-associated deformation gradient tensors satisfy the symmetry, i.e., K b X , ξ = K b X , ξ and F b X , ξ = F b X , ξ . The bond-associated first Piola–Kirchhoff stress tensors satisfies symmetry as well:
P b X , ξ = P b X , ξ .
As shown in Figure 2, similar to the BB PD model with rotation effect [12], the two bond-associated force density vectors in a bond are also equal in magnitude and opposite in direction, but not parallel to the bond, which can be written as
T b ¯ X , t ξ = T b ¯ X , t ξ .
Finally, in BA-NOSB PD, the equation of motion for the material point X at time t can be expressed as
ρ 0 ( X ) u ¨ X , t = H X T b ¯ X , t ξ T b ¯ X , t ξ d V X + b X , t .

2.1.3. JH2 Constitutive Model

The quasi-brittle materials also have important engineering applications, such as armored ceramics, bulletproof glass, concrete dams, and rock materials. The mechanical behavior and material response of these materials in different applications is an important basis for evaluating and measuring their safety in use and working life. To date, some constitutive models [40,41,42,43,44] describing quasi-brittle materials have been developed and used. The original Johnson–Holmquist (JH1) model [45] does not allow gradual softening, may damage self-healing, and is too sensitive to material model parameters. The JH2 model [46] is an improved version of the JH1 model, which is a widely used constitutive model for brittle or quasi-brittle materials, including ceramics, concrete, and rocks.
In the JH2 constitutive model, the equivalent stress is normalized as
σ * = σ i * D JH 2 σ i * σ f * ,
in which σ i * is the normalized intact equivalent stress of materials, σ f * is the normalized fracture equivalent stress of materials, and D JH 2 is dimensionless damage ( 0 D JH 2 1.0 ).
The true equivalent stress can be normalized as
σ * = σ / σ HEL ,
in which σ represents the true equivalent stress, σ HEL represents the equivalent stress at HEL , and HEL represents the Hugoniot elastic limit.
The normalized intact equivalent strength of materials is expressed as
σ i * = A P * + T * N 1 + C ln ε ˙ * .
The normalized fracture equivalent strength of materials is expressed as
σ f * = B P * M 1 + C ln ε ˙ * .
Also, the fracture equivalent strength of the material should satisfy σ f * SFMAX . Moreover, A, B, C, M, N, and SFMAX are material parameters. ε ˙ * = ε ˙ / ε ˙ 0 is the dimensionless strain rate, in which ε ˙ is the true strain rate and ε ˙ 0 = 1.0   s 1 .
The normalized hydrostatic pressure is expressed as
P * = P / P HEL ,
where P is the true pressure, and P HEL is the pressure at HEL .
The maximum tensile hydrostatic pressure is normalized to
T * = T / P HEL ,
where T is the maximum hydrostatic tension that can be sustained.
Similar to the JH1 model [45], the damage in the JH2 model is accumulated as
D JH 2 = Δ ε P / ε f P ,
in which Δ ε P is the plastic strain accumulated under the pressure P . The fracture plastic strain can be expressed as
ε f P = D 1 P * + T * D 2 ,
in which D 1 and D 2 are material parameters. In addition, the material can no longer withstand any plastic strain when P * = T * , but ε f P can increase with the increase in P * .
Before materials fracture ( D JH 2 = 0 ), the energy loss can be negligible. At this moment, the hydrostatic pressure is
P = K 1 μ + K 2 μ 2 + K 3 μ 3 , If   μ 0 K 1 μ , otherwise .
After fractured ( D JH 2 > 0 ), damage begins to accumulate. As the volumetric strain increases, volume bulking can occur. Due to energy loss, incremental pressure Δ P is added to the pressure P as
P = K 1 μ + K 2 μ 2 + K 3 μ 3 + Δ P ,
in which K 1 , K 2 , and K 3 are material parameters, and μ is the volumetric strain. The material is in tension when μ 0 , otherwise it is in compression.
The elastic internal energy can be written as
U = σ 2 / 6 G ,
where σ is the equivalent stress, and G is the shear modulus.
The incremental energy loss can be calculated as
Δ U = U D JH 2 t Δ t U D JH 2 t ,
where U D JH 2 t Δ t and U D JH 2 t are the elastic internal energies calculated from Equation (28) for the corresponding time step.
The updated incremental pressure Δ P is calculated as
Δ P t = K 1 μ t + K 1 μ t + Δ P t Δ t 2 + 2 β K 1 Δ U ,
in which β is the material parameter, and Δ P is the incremental pressure.

2.2. Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model

In this section, the original BA-NOSB PD model based on the JH2 constitutive relationship is modified to more accurately model the impact problems of quasi-brittle materials. First, the concept of “bond-associated volumetric strain” and its new calculation scheme are proposed via the introduction of the bond-associated deformation gradient. Then, a general tensile-shear coupling bond-breaking criterion is proposed in the framework of the BA-NOSB PD theory. Finally, a new treatment strategy for broken bonds is also given.

2.2.1. Bond-Associated Volumetric Strain

It is worth noting that all physical quantities are bond-associated within the BA-NOSB PD theory. Some improvements are proposed in this section to implement a fully bond-associated JH2 constitutive model.
The deformation of an object under the action of an external force can generally be classified as volume change and shape change. In the theory of plasticity, it is generally considered that the volume change is caused by spherical stresses, while the shape change is caused by deviatoric stresses. Therefore, the stress state at a point can be decomposed into a spherical stress state and a deviatoric stress state. In addition, the stress state at the point is determined by the stress tensor σ , which can be expressed as
σ = P I + S ,
in which P I indicates the spherical stress tensor, P is the pressure, I is a second-order identity tensor, and S indicates the deviatoric stress tensor.
In the framework of BA-NOSB PD theory, the above equation can be rewritten as
σ b = P b I + S b ,
where S b represents the bond-associated deviatoric stress tensor, which is given in Section 2.3. According to Equation (26), the bond-associated hydrostatic pressure P b can be calculated as
P b = K 1 μ b + K 2 μ b 2 + K 3 μ b 3 , If   μ b 0 K 1 μ b , otherwise .
As seen in Equation (33), the bond-associated hydrostatic pressure P b depends on the bond-associated volumetric strain μ b . In the original JH2 constitutive model [45], the volumetric strain μ in the equation of state is defined as
μ = ρ / ρ 0 1 ,
where ρ 0 is the initial mass density, and ρ is the current mass density.
The volumetric strain μ in Equation (34) depends only on the mass density of a point, and can also be expressed as
μ = μ p = d V / d v 1 ,
in which the subscript b indicates that the physical quantity is point-associated, d V is the reference volume element, and d v is the current volume element.
Referring to classical continuum mechanics, and according to the third invariant of the bond-associated deformation gradient tensor F b in Section 2.1.2, J b = det F b , the transformation relationship between volume elements under different configurations can be expressed as
d V = J b d v .
The equation above holds for any shaped volume element since it can be approximated by an infinite number of parallel hexahedra. Substituting Equation (36) into Equation (35), the bond-associated volumetric strain μ b can be defined as
μ = μ b = 1 / J b 1 .
Obviously, μ p is a point-associated nonlocal physical quantity, while μ b is a bond-associated nonlocal physical quantity. Without losing generality, in addition to the JH2 model [45], the defined bond-associated volumetric strain μ b can also be used in other constitutive models [46,47,48] implemented in the framework of the BA-NOSB PD theory. Further details of the constitutive update using the JH2 model will be given in Section 2.3.

2.2.2. Bond-Breaking Criterion

If the external force exceeds the fracture strength of materials, local damage, cracks or fractures will occur. The PD can naturally model discontinuous behavior such as damage or cracks, and damage is also incorporated into its constitutive model. In the PD theory, the bonds are carriers that transmit non-local interactions, while damage is defined by breaking bonds. In addition, broken bonds will no longer transmit PD interaction. Accurately assessing and describing the discontinuity behavior of materials usually requires a robust bond-breaking criterion.
In most engineering problems, materials usually suffer shear failure, not just tensile and compressive failure. According to Equation (16), the bond-associated force density vector can describe not only tensile-compressive deformation, but also shear deformation. However, the critical stretch criterion [8] proposed earlier cannot evaluate material damage caused by shear deformation. Therefore, the assessment of shear damage needs to be refined in BA-NOSB PD. Based on the comprehensive failure criterion proposed by Liu et al. [14], a new general tensile-shear coupling bond-breaking criterion is proposed and implemented in the framework of BA-NOSB PD theory. As shown in Figure 3, under the condition of small deformation, the relative rotation angle vector [12] is approximated as
γ = η s ξ n ξ
where n = ξ + η / ξ + η is an identity direction vector of the bond ξ in the current configuration, and ξ = ξ .
The bond stretch in PD is defined as
s ξ , η = ξ + η ξ ξ .
Similarly, the bond-breaking threshold used to assess the tensile-shear failure of materials in this work is defined by the dimensionless quantity λ , which can be expressed as
λ = α 0 s s 0 2 + β 0 γ γ 0 2 ,
in which s 0 is the critical stretch, γ 0 is the critical relative rotation angle, and γ = γ is the relative rotation angle. The difference between the anti-tensile capability and anti-shear capability of materials is considered in this work by defining the anti-tensile weight α 0 and anti-shear weight β 0 respectively, as follows:
α 0 = s 0 s 0 + γ 0 β 0 = γ 0 s 0 + γ 0 .
In particular, the tensile-shear coupling bond-breaking criterion degrades to the critical stretch criterion when the anti-tensile weight and anti-shear weight are taken as α 0 = 1 and β 0 = 0 , respectively.
Further, the critical stretch s 0 and the critical relative rotation angle γ 0 are defined by the fracture energy of materials. In the PD theory, the energy release rate G 0 [8] and the shear fracture energy G s [14] are respectively expressed as
G 0 = π c s 0 2 δ 5 10 G S = π c γ 0 2 δ 5 10 .
For the three-dimensional case, the critical stretch is calculated as
s 0 = 5 G 0 / 9 K δ ,
where K represents the bulk modulus.
Similar to the derivation of the critical stretch s 0 by Silling and Askari [8], substituting the bond constant c = 18 K / π δ 4 of the PMB model into Equation (42), the critical relative rotation angle γ 0 is expressed as
γ 0 = 5 G s / 9 K δ .
In the PD theory, the local damage of the material point X at the instant of t can be defined as
φ X , t = 1 H X μ ξ , t d V X H X d V X ,
where the scalar function μ ξ , t is
μ ξ , t = 1 , λ < 1.0 0 , λ 1.0

2.2.3. Treatment Strategy of Broken Bonds

Unlike BB PD theory, the force density vectors at a material point depend on the deformation behavior of all bonds in SB PD. Therefore, the description of broken bonds in the state-based theory is very important, since it will affect the calculation update of the force vector state for unbroken bonds. The following two treatment strategies are given for a broken bond, which will be compared and discussed in the subsequent Section 3.2. If the bond ξ breaks, in each subsequent time-integral cycle, the relative position vector between two material points in the current configuration is considered as
Y ¯ B ξ = ξ ,
Y ¯ B ξ = ξ + η B ,
in which the subscript B indicates that the bond has been broken, and η B is the relative displacement vector between two material points when the bond ξ breaks.

2.3. Constitutive Update Scheme

In this section, the constitutive update scheme using the JH2 model is given in the framework of the BA-NOSB PD theory. First, based on the bond-associated deformation gradient tensor F b , the bond-associated spatial velocity gradient tensor L b can be calculated as
L b = F ˙ b F b 1 ,
where F ˙ b is the material derivative of the deformation gradient F b , which can be calculated as
F ˙ b = H ξ ω ¯ ξ Y ¯ ˙ ξ X ¯ ξ d V X K b 1 .
The spatial velocity gradient L b is decomposed into a bond-associated deformation rate tensor D b and a bond-associated rotation rate tensor W b :
D b = 1 2 L b + L b T and W b = 1 2 L b L b T .
To ensure that the resulting constitutive relationship satisfies the principle of coordinate invariance, the algorithm proposed by Rubinstein and Atluri [49] is used to update the bond-associated Cauchy stress σ b . Before this, the bond-associated unrotated deformation rate of the tensor d b can be expressed as
d b = R b t T D b R b t ,
where R b t is the bond-associated orthogonal rotation tensor that describes the rigid-body rotation, which can be calculated as
R b t = I + sin Δ t Ω Ω Ω 1 cos Δ t Ω Ω 2 Ω 2 R b t Δ t ,
in which the rigid-body rotation rate tensor Ω can be calculated as
z i = e i k j D b j m V b m k w i = 1 2 e i j k W b j k ϖ = w + tr V b I V b 1 z Ω i j = e i k j ϖ k ,
where e i k j is the permutation symbol, and the left stretch tensor is
V b t = V b t Δ t + Δ t V ˙ b Δ t ,
in which V ˙ b Δ t = L b V b t V b t Ω is the left stretch tensor rate.
As of now, we have calculated the bond-associated unrotated deformation rate of the tensor d b . In the JH2 constitutive relationship, the bond-associated unrotated Cauchy stress τ b is calculated using the bond-associated unrotated deformation rate of the tensor d b . At each time step, we will compute and store the rotated Cauchy stress σ b t Δ t at the previous step through the stored unrotated Cauchy stress τ b t Δ t , as
τ b t Δ t = R b t Δ t T σ b t Δ t R b t Δ t ,
where V b = I and R b = I when t = 0 .
We first assume that the material deformation is in the elastic stage, then the total elastic strain increment Δ e and the partial elastic strain increment Δ e dev are respectively expressed as
Δ e = d b Δ t , Δ e dev = Δ e 1 3 I .
At the current time t , the bond-associated unrotated trial Cauchy stress is
τ b t , trial = τ b t Δ t + λ 1 3 tr Δ e I + 2 G e dev
where λ and G are the Lamé constants of solids. Based on Von–Mises plasticity theory, the bond-associated trial deviatoric stress tensor can be calculated as
S b = τ b t , trial 1 3 tr τ b t , trial I
The bond-associated equivalent Von–Mises yield stress is
S b VM = 3 2 ( S b ) i j ( S b ) i j
The normalized strain rate is ε ˙ * = d b VM / ε ˙ 0 , in which d b VM = 2 3 d b i j d b i j is the total equivalent strain rate. Furthermore, the bond-associated hydrostatic pressure is calculated as
P b t = K 1 μ b t + K 2 μ b t 2 + K 3 μ b t 3 + Δ P b t Δ t , If   μ b t 0 K 1 μ b t + Δ P b t Δ t , otherwise
Based on Equation (60), the bond-associated yield function is expressed as
f σ b , ε ˙ b = S b VM σ b * σ HEL .
If f σ b , ε ˙ b < 0 , the material is in the elastic stage, and the bond-associated unrotated Cauchy stress is equal to the bond-associated unrotated trail Cauchy stress. Otherwise, the material yields, and we need to recalculate the bond-associated true elastic strain increment, plastic strain increment, and unrotated Cauchy stress.
For the three-dimensional stress state, the bond-associated plastic strain rate vector is calculated according to Prandtl–Reuss flow law as
ε ˙ b p = λ ˙ f σ b , ε ˙ b S b = λ ˙ a ,
in which λ ˙ is the bond-associated plastic strain rate multiplier, and a is an orthogonal vector. When the plastic deformation occurs, the stress rate or stress increment is perpendicular to the yield surface, as follows:
f σ b , ε ˙ b = a T σ ˙ b .
When the material yields, the total strain rate is decomposed into the elastic component and the plastic component. The relationship between stress and stress rate can be expressed as
σ ˙ b = C ε ˙ b ε ˙ b p = C ε ˙ b λ ˙ a ,
where C is a stiffness matrix of the material. The bond-associated plastic strain rate multiplier λ ˙ is
λ ˙ = a T C ε ˙ b a T C a .
By considering Equations (64)–(66), the bond-associated equivalent plastic strain increment is
Δ ε b p = a T C ε ˙ b Δ T 2 a T C a S b x x S b y y 2 + S b x x S b z z 2 + S b y y S b z z 2 + 6 S b x y 2 + S b y z 2 + S b z x 2 .
If the bond-associated equivalent plastic strain increment is not currently zero, the bond-associated fracture damage factor accumulates as
D b t = D b t Δ t + Δ ε b p ε b , f p
The incremental energy loss is
Δ U = U t Δ t σ HEL σ b * 2 6 G .
In the current configuration, the bond-associated final hydrostatic pressure is recalculated as
P b t = 1 3 tr τ b t , trail + Δ P b t .
Finally, the bond-associated unrotated Cauchy stress tensor and rotated Cauchy stress at the current time step can be updated to
τ b t = P b t I + S b ,
σ b t = R b t τ b t R b t T .
Thereafter, we can calculate the stress tensor P b in Equation (13), and then update the force state T b ¯ at the current time step in Equation (12).

2.4. Other Zero-Energy Mode Control Schemes

In addition to the bond-associated deformation gradient scheme proposed in BA-NOSB PD, different zero-energy mode control schemes have been proposed to suppress or eliminate numerical oscillations. In this section, we briefly introduce four zero-energy mode control schemes based on the supplemental force state. In the NOSB model based on the supplemental force state, the total force state can be written as
T ¯ total X , t ξ = T ¯ X , t ξ + T ¯ s X , t ξ ,
in which the original force state T ¯ is defined in Equation (5), and T ¯ s is the corresponding supplemental force state.
  • Control Scheme I:
In Littlewood’s work [26], a penalty scheme is used to control the zero-energy mode. This is similar to the method of controlling deformable hourglass modes in finite element analysis by adding the hourglass force to the original force state. The added hourglass force is expressed as
T ¯ s X , t ξ = c hg c h proj ξ ξ + η ξ + η Δ V X Δ V X ,
in which c hg is a user parameter controlling the level of the hourglass force, c is the bond constant of the PMB model, and c = 18 K / π δ 4 .
The projection of the hourglass vector h on the relative position vector ξ + η under the current configuration can be expressed as
h proj = h ξ + η ,
where the hourglass vector h is
h = y y ,
in which the position of the neighbor X in the current configuration is calculated by the approximative deformation gradient F as:
y = y + F ξ .
  • Control Scheme II:
In addition, Silling [27] derived the stability conditions of the NOSB PD model and proposed a new supplemental force state that can be written as
T ¯ s X , t ξ = ω ¯ ξ G c 1 ω 0 z ¯ ξ ,
where G is a user parameter controlling the level of the supplemental force state, and c 1 = c / δ is the micro-modulus in BB PD. The sum of the influence function values of all bonds at the material point X is
ω 0 = H X ω ¯ ξ d V X .
Furthermore, the non-uniform deformation state can be expressed as
z ¯ ξ = Y ¯ ξ F ξ .
  • Control Scheme III:
Li et al. [28] proposed a novel supplemental force state characterizing the non-uniform deformation state based on the linearized BB PD theory, which can be written as
T ¯ s X , t ξ = ω ¯ ξ C ¯ ξ z ¯ ξ ,
where C ¯ ξ is the elasticity coefficient tensor, defined as
C ¯ ξ = c ξ ξ / ξ 3 .
Compared with the control schemes I and II, this scheme does not require the user parameter, and avoids the complicated parameter adjustment process.
  • Control Scheme IV:
Wan et al. [29] proposed a new supplemental force state directly within the NOSB PD model, which can be expressed as
T ¯ s X , t ξ = ω ¯ ξ C e K 1 z ¯ ξ ,
where C e is the elastic modulus tensor. Similar to the Control Scheme III, this supplemental force state also does not contain the user parameter.

2.5. Numerical Discretization and Implementation

2.5.1. Numerical Discretization

The solution to the PD governing equation of motion usually requires discrete spatial domain and time domain for numerical integration. In the PD theory, the geometric model of the problem domain is presented by discretized material points. Accordingly, in the BA-NOSB PD model, the discrete form of the equation of motion for the material point X i at the instant of t can be written as
ρ 0 ( X i ) u ¨ X i , t = j N T b ¯ X i , t ξ i j T b ¯ X j , t ξ j i V X j + b X i , t ,
where N is the number of neighbors for the material point X i , and ξ i j = X j X i .
The discrete forms of the bond-associated shape tensor and deformation gradient tensor about the bond ξ i j at the material point X i is
K b ( X i , ξ i j ) = k = 1 N b ω ¯ ξ i k ξ i k ξ i k Δ V X k ,
F b ( X i , ξ i j ) = k = 1 N b ω ¯ ξ i k ξ i k + η i k ξ i k Δ V X k K b 1 ( X i , ξ i j ) ,
where N b is the number of bond-associated neighbors with respect to the bond ξ i j at the material point X i .
Finally, we use the explicit Velocity–Verlet integration algorithm to implement the time-integral scheme, which can be expressed as
u ˙ t + Δ t = u ˙ t + 1 2 u ¨ t + u ¨ t + Δ t Δ t ,
u t + Δ t = u t + u ˙ t Δ t + 1 2 u ¨ t Δ t 2 ,
where Δ t is the time step size of the time-integral cycle.

2.5.2. Artificial Viscosity

In numerical simulations of the impact or penetration problems, jump oscillations often occur in the impact domain, which may lead to unphysical deformation modes. Artificial viscosity [50,51] is widely used in penetration or impact simulations to improve the stability of numerical algorithms and to prevent unphysical overlapping or interpenetration between particles. In this work, the artificial viscosity implemented in PD theory by Lai et al. [20] is used for impact problems, which can be expressed as
i j = α c ¯ i j ϕ i j + β ϕ i j 2 ρ ¯ i j , v i j x i j < 0 0 , v i j x i j 0 ,
where
ϕ i j = δ i j v i j x i j x i j 2 + φ 2 ,
c ¯ i j = 1 2 c i + c j ,
ρ ¯ i j = 1 2 ρ i + ρ j ,
δ i j = 1 2 δ i + δ j ,
v i j = v i v j   and   x i j = x i x j .
In the above Equations, α and β are constants with values of about 1.0. The factor φ = 0.1 δ i j is taken into account to avoid numerical divergence. The velocity vector of the particle is v . The sound speed, mass density, and horizon of the particle are c , ρ , and δ , respectively.
The artificial viscosity state introduced can be written as
¯ X i , t ξ i j = ω ξ i j i i j ξ i j ,
in which is the vector differential operator. When v i j x i j < 0 ,
i i j = α c ¯ i j + 2 ϕ i j β ρ ¯ i j ϕ i j x i ,
where
ϕ i j x i = δ i j v i x i Y ¯ ξ i j Y ¯ ˙ ξ i j 1 Y ¯ ξ i j 2 + φ 2 + δ i j Y ¯ ˙ ξ i j Y ¯ ξ i j 2 Y ¯ ξ i j Y ¯ ξ i j 2 + φ 2 2 ,
where
v i x i = j = 1 N ω ¯ ξ i j Y ¯ ˙ ξ i j ξ i j Δ V X j K X i 1 F X i 1
Y ¯ ˙ ξ i j = v X j , t v X i , t
It is worth noting that the introduced artificial viscosity is not physical viscosity, but only numerical correction. This correction prevents the unphysical deformation behavior of particles in dynamic fracture problems such as impact or penetration.

2.5.3. Contact Algorithm

The impact or penetration simulation is a critical problem often encountered in fields such as industrial manufacturing, military equipment, and scientific research. Such problems often involve contact and collision analysis between two or more objects. Among the interface contact algorithms [52,53,54] that have been developed and used, the penalty method is the most commonly used [55]. In this work, the penalty method implemented in PD theory by Lai et al. [20] will be used to model impact problems. Figure 4 shows the normal impact problem for two objects A and B . The material points X i and X j are on the surfaces of the objects A and B , respectively. At the initial moment, the material point X j remains stationary, and the velocity of the material point X i is v i . The relative position vector at the current time step is g n = x j x i , and l = g n .
When l < l c , contacts occur between the material points of the master and slave, and a contact bond ξ i j = X j X i is formed. At the current moment t , the contact force of the material points X j on X i can be expressed as
f c X i , t ξ i j = g n / l C ξ V X i ln l / l c , if   l < l c 0 , otherwise ,
where l c is a standard value and is usually chosen to be slightly smaller than the grid size of the discretization. The contact stiffness C ξ is defined on the contact bond, and can be calculated as
C ξ = min ρ i , ρ j × 2 × 10 n ,
where ρ i and ρ j are the mass densities, and n is the penalty factor. It is worth noting that contact forces can only be repulsive forces.
Therefore, the total contact force at the material point X i is expressed as
F c X i , t = i = 1 M f c X i , t ξ i j , ξ i j < δ c ,
where M is the number of material points forming contact bonds, and δ c is a contact horizon.

3. Numerical Examples

In this section, four numerical examples are established to verify and demonstrate the modified BA-NOSB PD model. In the first part, two numerical examples are used to numerically model the uniaxial tension experiments of elastic bars under different loading forms to verify the effectiveness and robustness of the modified model. In the second part, two numerical experiments of the edge-on impact and normal impact of ceramics are carried out to demonstrate the ability of the modified model for impact problems of quasi-brittle materials. Good agreement is shown by comparing PD simulation results with references.

3.1. Uniaxial Tension of an Elastic Bar under Stress Loading

In this section, the numerical modeling of the uniaxial tension test of a three-dimensional elastic bar under stress loading is carried out, and the computational model is seen in Figure 5a. The configuration of the model is 80 × 10 × 1 0   mm 3 . A fixed constraint is applied at one end of the bar, and a dynamic tensile stress σ x = 5.0   kPa is applied to the other end. The specimen material used in this work is an Al2O3 ceramic with elastic modulus E = 220   GPa . The parameters of the JH2 constitutive model for the Al2O3 ceramic are listed in Table 1. As shown in Figure 5b, the loading end consists of one layer of red material points, and the fixed end consists of m layers of light green material points. In the PD simulation, the boundary conditions are as follows: the dynamic tensile stress is applied to red material points, and the fixed constraint is applied to light green material points, as seen in Figure 5b. The ceramic bar is discretized into cubic particles with the grid size Δ x , where Δ x = 1.0 × 10 3   m . The horizon size of the discretization model is δ = m · Δ x , and the time step size is Δ t = 1.0 × 10 8   s . The model is discretized into 9801 material points.
This section will verify the effectiveness of controlling zero-energy modes of the modified BA-NOSB PD model and compare the numerical stability of different control schemes. The axial displacement simulation results from different control schemes are compared with those of the original NOSB PD model. In comparison and analysis, simulation results at t = 100   μ s are selected. In this example, the grid size Δ x of the discretization model remains fixed, and five cases with horizon factors m of 1.5, 2.0, 3.0, 4.0, and 5.0 are selected for comparison and analysis of the results.
Figure 6 shows the axial displacement contours of the 3D ceramic bar calculated by the original NOSB PD model, four other zero-energy mode control schemes, and the modified BA-NOSB PD model when m = 1.5 . The displacement solution obtained from the original NOSB PD model in Figure 6a shows very violent oscillations. The result based on the Control Scheme I in Figure 6b is similar to that of the original NOSB PD model, and there are also severe oscillations. As shown in Figure 6d,e, the displacement solutions obtained by the Control Schemes III and IV are obviously incorrect, and there are obvious oscillations near the fixed end. The displacement solutions of the Control Scheme II and the modified BA-NOSB PD model show discontinuity near the fixed end, and the results of the Control Scheme II are smoother and more continuous than that of the modified BA-NOSB PD model, as shown in Figure 6c,f.
Figure 7 shows the axial displacement contours calculated by the different PD models when m = 2.0 . The displacement solution obtained from the original NOSB PD model in Figure 7a shows slight oscillations, which are more significant near the fixed end, and the displacement contour in the middle of the bar is not smooth and continuous. The results of the Control Schemes I, III, and IV are similar to that of the original NOSB PD model, all with slight oscillations, as seen in Figure 7b,d,e. In addition, the Control Scheme II and the modified BA-NOSB PD model both obtain smooth displacement solutions, but the results of the Control Scheme II have slight oscillations near the fixed end, as seen in Figure 7c,f.
Figure 8 shows the axial displacement contours calculated by the different PD models when m = 3.0 . The displacement solution obtained from the original NOSB PD model in Figure 8a shows obvious oscillations near the fixed end and slight oscillations near the loading end. In addition, the displacement contour in the middle of the bar is not smooth and continuous. The displacement solutions obtained by the Control Schemes I, III, and IV are similar. These results have obvious oscillations near the fixed end, and the displacement contours in the middle of the bar are not smooth and continuous, as shown in Figure 8b,d,e. The result obtained by the Control Scheme II in Figure 8c has slight oscillations near the fixed end. In addition, the modified BA-NOSB PD model provides a smooth displacement solution that can effectively eliminate numerical oscillations, as shown in Figure 8f.
Figure 9 shows the axial displacement contours calculated by the different PD models when m = 4.0 . There are slight oscillations in the displacement solution obtained from the original NOSB PD model in Figure 9a, and the displacement contour in the middle of the bar is not smooth and continuous. The results of the Control Schemes I, III, and IV are similar. Moreover, these results have slight oscillations near the fixed end, and the displacement contours in the middle of the bar are not smooth and continuous, as shown in Figure 9b,d,e. The result obtained by the Control Scheme II in Figure 9c shows slight oscillations near the fixed end. The modified BA-NOSB PD model provides smooth and continuous displacement contours, as shown in Figure 9f.
Figure 10 shows the axial displacement contours calculated by the different PD models when m = 5.0 . The displacement solution obtained from the original NOSB PD model in Figure 10a shows obvious oscillations near both the fixed end and the loading end. The results of the Control Schemes I, III, and IV all show obvious oscillations at the fixed end and in the middle of the bar, as shown in Figure 10b,d,e. In addition, the result of the Control Scheme II has slight oscillations near the fixed end, and the displacement contour in the middle of the bar is not smooth and continuous, as shown in Figure 10c. The displacement solution obtained from the modified BA-NOSB PD model in Figure 10f has no obvious oscillations, but it does show discontinuous features near the fixed end.
This section simulates a uniaxial tension experiment of a 3D ceramic bar under stress loading, and compares the effectiveness of different zero-energy mode control schemes to suppress or eliminate numerical oscillations. The following conclusions can be drawn from the above comparison and analysis. First, the simulation results of the original NOSB PD model show different degrees of oscillations for different values of m . Similarly, the results of the Control Schemes I, III, and IV all show different degrees of oscillations for different values of m , which fail effectively to suppress or eliminate numerical oscillations. In addition, the Control Scheme II effectively eliminates oscillations at smaller values of m ( 1.5 m 3.0 ), while slight oscillations exist at larger values of m ( 4.0 m 5.0 ). Finally, the modified BA-NOSB PD model can effectively eliminate oscillations for different values of m , showing slight discontinuous characteristics in the case of m = 1.5 or m = 5.0 .
On the other hand, by analyzing the above simulation results with numerical oscillations, two conclusions can be drawn. First, the oscillations of the result are more obvious when a large or small horizon factor m ( m 2.0 or m 4.0 ) is selected, while the oscillations of the result are lighter when a moderate horizon factor m ( m = 3.0 ) is selected. Secondly, the oscillations of the simulation results mainly exist in the region near the boundaries, such as the fixed end or the loading end.
In summary, the modified BA-NOSB PD model can effectively eliminate the numerical oscillations and unphysical deformation modes. Compared to four other zero-energy mode control schemes, the modified model has a more robust zero-energy mode control capability.

3.2. Uniaxial Tension of an Elastic Bar under Velocity Loading

In order to more accurately implement the JH2 constitutive model in the framework of BA-NOSB PD theory, the bond-associated volumetric strain μ b is defined in Section 2.2.1. In addition, for the SB PD theory, a new treatment strategy for broken bonds is proposed in Section 2.2.3. In this section, for the volumetric strain in the equation of state and the broken bonds, the following three solutions are presented to modify the original BA-NOSB PD model, as shown in Table 2.
To evaluate the above three solutions, the uniaxial tension experiment of the 3D elastic bar under velocity loading will be numerically modeled in this section, and the computational model is shown in Figure 11a. The configuration of the model is 80 × 10 × 10   mm 3 . At the initial moment, an instantaneous tensile velocity load of the same magnitude is exerted on both ends of the bar, and the magnitude of the velocity load is v X = 0.1   m / s . The sample material is the Al2O3 ceramic with the same material parameters as in Section 3.1. The JH2 constitutive model parameters of Al2O3 ceramic are shown in Table 1. In the PD simulation, the tensile velocity loads are applied to red and light green material points, respectively, as seen in Figure 11b. The grid size of the discretization is Δ x = 1.0 × 10 3   m , and the horizon factor is chosen as m = 3.0 . The time step size is Δ t = 1.0 × 10 8   s . The model is discretized into 9801 material points.
From the perspective of kinetic energy convergence, the change in the system kinetic energy over time will be compared to discuss and evaluate the above three solutions. At the initial moment, both ends of the bar are subjected to a tensile velocity load, and the initial system kinetic energy is 1.85 × 10 6   J . This initial kinetic energy is also the total energy of the whole system. Figure 12 shows the change in the total kinetic energy of the PD simulation system based on the three solutions over time. For Solution A, the ceramic bar begins to show damage at t = 65   μ s , and the system kinetic energy increases sharply to 9.902 × 10 4   J from the damage occurrence to complete brittle fracture. For Solution B, the bar begins to show damage at t = 82.5   μ s , and the system kinetic energy increases significantly to 2.701 × 10 4   J from the damage occurrence to complete brittle fracture. The final system kinetic energy for Solutions A and B is much higher than the initial value, and the system kinetic energy will increase significantly after the damage occurs, as shown in Figure 12a. Completely different from Solutions A and B, the system kinetic energy of Solution C does not change significantly from the damage appearance to brittle fracture, as shown in Figure 12b,c. Further, as shown in Figure 12c, the system kinetic energy of Solution B increases rapidly after the damage occurs, while the system kinetic energy of Solution C remains stable. Compared to the initial value, it can be seen that the system kinetic energy of Solution C is always lower than the initial value and does not increase significantly or far exceed the initial value, as shown in Figure 12d.
By comparing and analyzing the convergence curves of the system kinetic energy with the three solutions, the results are summarized as follows. First, the system kinetic energy of Solutions A and B both increase dramatically after the damage occurs, and the final values are much larger than the initial value. Second, the system kinetic energy of Solution C changes steadily over time and is always below the initial value. The results show that only Solution C satisfies the requirement of kinetic energy convergence, while both Solutions A and B do not. In addition, the system kinetic energy of Solutions A and B in the final stage is much greater than the total energy, which means that neither meets the requirement of energy conservation. From the perspective of kinetic energy convergence, Solution C provides better accuracy and robustness for the BA-NOSB PD model and is chosen as the final modified solution.

3.3. Edge-On Impact Simulation of Ceramics

In this section, the numerical experiments on the edge-on impact of ceramic plates will be carried out. The contact algorithm in Section 2.5.3 will be used to describe the interaction between the cylindrical projectile and the target plate. Further, the artificial viscosity introduced in Section 2.5.2 will be used to prevent unphysical oscillations or interpenetration in the impact zone. Figure 13 shows the computational model of a ceramic plate under edge-on impact. The configuration of the target plate is 70 × 70 × 10   mm 3 . In addition, the projectile is a cylinder with a diameter of 16   mm and a height of 21   mm . The target plate is the Al2O3 ceramic, and its material parameters and JH2 constitutive model parameters are the same as those in Section 3.1. For simplicity, the projectile is regarded as a rigid body with density ρ = 8060   kg / m 3 . With the same set-up as in the reference, free boundary conditions are applied to the ceramic target plate, and the impact velocity of the projectile is v = 85   m / s . In the PD simulation, the grid size of the discretization is Δ x = 1.0 × 10 3   m , and the horizon factor is chosen as m = 3.0 . The time step size used in the edge-on impact is Δ t = 2.0 × 10 9   s . The computational model is discretized into 55,510 material points.
Figure 14 shows the fracture patterns of the Al2O3 ceramic plate under edge-on impact in the experiment and numerical simulation. Comparing the PD simulation result with the fracture image provided by Strassburger [57], good agreement can be found. Specifically, both the delta-like fracture pattern and the radial crack propagation paths in Figure 14a can be observed from the simulation result in Figure 14b. In addition, the simulation result shows the characteristics of the fragmentation of the target edge and the semi-elliptical damage zone, which is very similar to the fracture behavior observed in the experiment. The results demonstrate that the modified BA-NOSB PD model can accurately reproduce the fracture pattern and the crack propagation paths of the Al2O3 ceramic plate under the edge-on impact experiment.
The damage evolution and crack propagation of the Al2O3 ceramic plate during edge-on impact are shown in Figure 15. After the cylindrical projectile hit the ceramic target plate, it began to penetrate the ceramic plate. At t = 2.0   μ s , the left edge of the target plate began to display gear-shaped damage, as shown in Figure 15a. At t = 4.0   μ s , the damage of the target plate further increases and extends along the impact direction, as shown in Figure 15b. With the continuous action of the projectile, five radial cracks form on the edge of the target plate, as shown at the instant of t = 6.0   μ s in Figure 15c. At t = 8.0   μ s , a triangular crack extension zone and a semicircular damage zone are formed on the target plate, as shown in Figure 15d. Subsequently, at t = 8.8   μ s , a corrugated damage zone and radial cracks are formed on the edge of the target plate, as shown in Figure 15e. Finally, with the weakening of the projectile action, the damage tended to be stable, and a fracture pattern similar to an alluvial delta is formed on the front surface of the ceramic plate, as shown at the instant of t = 9.7   μ s in Figure 15f.

3.4. Normal Impact Simulation of Ceramics

In this section, the normal impact experiment on a ceramic plate will be numerically modeled. The contact algorithm in Section 2.5.3 will be used to describe the interaction of the spherical projectile with the target plate. The artificial viscosity in Section 2.5.2 is incorporated to avoid unphysical oscillations or interpenetration in the impact area. The computational model of the normal impact simulation for the ceramic plate is shown in Figure 16. The configuration of the ceramic plate is 70 × 70 × 9   mm 3 . Moreover, the projectile is a sphere with a diameter of 6.4   mm . The material of the target plate is the Al2O3 ceramic with the same material parameters as in Section 3.1. The JH2 constitutive model parameters of the Al2O3 ceramic are shown in Table 1. For simplicity, the projectile material is set to be rigid, and its density is ρ = 7800   kg / m 3 . With the same set-up as in the reference, the ceramic plate is subjected to free boundary conditions, and the impact velocity of the projectile is v = 200   m / s . In the PD simulation, the grid size of the discretization is Δ x = 1.0 × 10 3   m , and the horizon factor is chosen as m = 3.0 . The time step size in normal impact is Δ t = 2.0 × 10 9   s . The model is discretized into 47,256 material points.
Figure 17 shows the fracture patterns of the upper surface of the Al2O3 ceramic plate under normal impact obtained from the experiment and PD simulation. Comparing the experimental result [58] and the simulation result, good agreement can be found on the fracture patterns. The PD simulation result in Figure 17b accurately reproduces the eight radial cracks along the 0°, 45°, 90°, 135°, 180°, 225°, 270°, and 315° directions of the ceramic plate in the experiment. In addition, the PD simulation result can accurately reproduce the characteristics of fragmentation and collapse in the center of the plate caused by projectile penetration. These fracture characteristics agree well with the experimental result in Figure 17b. The results show that the modified BA-NOSB PD model can accurately capture the crack propagation paths of the Al2O3 ceramic plate under normal impact, demonstrating the ability to model impact problems and the reliability of damage prediction. In particular, the established model does not take into account the asymmetric factors that may exist in the experiment such as the projectile deflection and the microscopic defects of the sample. Therefore, the PD simulation result does not show the asymmetry of the crack path in the experiment.
The damage evolution and crack growth of the Al2O3 ceramic plate during normal impact are shown in Figure 18. With the impact of the steel ball on the ceramic plate, the compression stress waves generated in the center of the target constantly spread to the boundary surfaces. The compression stress waves are reflected as the tensile stress waves once they have propagated to a free surface. When t = 2.0   μ s , the upper surface of the target forms a square damage zone. Two circumferential cracks appear in the central damage zone, and cracks initiate outward at the boundary of the damage zone, as shown in Figure 18a. With the continuous propagation, reflection, and overlapping of stress waves, the square damage zone on the target surface expands further, and a collapse is observed in the central area of the target. In addition, four radial cracks along the directions of 45°, 135°, 225°, and 315° begin to form, as shown at the instant of t = 4.0   μ s in Figure 18b. When t = 6.0   μ s , an obvious circular crack forms in the center of the target, and four radial cracks along the directions of 0°, 90°, 180°, and 270° initiate and grow outward, as seen in Figure 18c. Under the continuous penetration of the steel projectile, the crack tips along the directions of 0°, 90°, 180°, and 270° grow outward, and then obvious radial cracks form, as shown at the instant of t = 8.0   μ s in Figure 18d. When t = 8.8   μ s , four radial cracks form on the upper surface of the target along the directions of 0°, 90°, 180°, and 270°, as shown in Figure 18e. Finally, the damage tends to be stable at t = 9.5   μ s , while the center of the target produces fracture behaviors such as fragmentation and collapse, and the top surface produces eight radial cracks along the 0°, 45°, 90°, 135°, 180°, 225°, 270°, and 315° directions, as shown in Figure 18f.

4. Conclusions

In this work, a novel bond-associated non-ordinary state-based peridynamic model for quasi-brittle materials is proposed. The JH2 constitutive model is implemented in the framework of the bond-associated non-ordinary state-based peridynamic theory in the first place, to describe the mechanical behavior of quasi-brittle materials. A new calculation scheme of the volumetric strain in the constitutive model is proposed with the introduction of the bond-associated deformation gradient to improve stability and accuracy in describing the bulking effect. A new treatment strategy for broken bonds is proposed as well. In addition, considering the effect of tensile-shear deformation on material damage, a universal tensile-shear coupling bond-breaking criterion is proposed and implemented in the bond-associated non-ordinary state-based peridynamics.
The effectiveness and robustness of the modified bond-associated non-ordinary state-based peridynamic model are verified to some extent by two benchmark examples. The edge-on impact and normal impact experiments on ceramics are numerically modeled and simulated. The results reveal that the proposed model can effectively control zero-energy mode and eliminate numerical oscillations. Compared with the other four control schemes listed in this work, the modified bond-associated model has a more stable zero-energy mode control capability and convenience. In addition, from the energy conservation analysis of the system, one can find better accuracy and robustness from the modified bond-associated non-ordinary state-based peridynamic model. Good agreement can be found with the experimental data in the two demonstration benchmark tests, which implies that the modified bond-associated non-ordinary state-based peridynamic model can accurately reproduce the fracture patterns and crack growth paths of ceramics in the experiment.

Author Contributions

Conceptualization, J.Z., Y.L., X.L. (Xin Lai) and L.L.; Methodology, J.Z., Y.L., X.L. (Xin Lai) and L.L.; Software, J.Z.; Validation, J.Z.; Formal analysis, J.Z., H.M. and X.L. (Xiang Liu); Investigation, J.Z., Y.L. and H.M.; Resources, X.L. (Xin Lai), L.L. and H.M.; Writing—original draft, J.Z.; Writing—review & editing, X.L. (Xin Lai); Visualization, H.M. and X.L. (Xiang Liu); Supervision, X.L. (Xin Lai), L.L., H.M. and X.L. (Xiang Liu); Project administration, X.L. (Xin Lai) and L.L.; Funding acquisition, X.L. (Xin Lai) and L.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (Nos. 11802214, 11972267, and 51932006) and the Fundamental Research Funds for the Central Universities (WUT: 223114010).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used to support the findings of this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Walley, S.M. Historical review of high strain rate and shock properties of ceramics relevant to their application in armour. Adv. Appl. Ceram. 2010, 109, 446–466. [Google Scholar] [CrossRef]
  2. Chen, W.W.; Rajendran, A.M.; Song, B.; Nie, X. Dynamic fracture of ceramics in armor applications. J. Am. Ceram. Soc. 2007, 90, 1005–1018. [Google Scholar] [CrossRef]
  3. Song, J.H.; Wang, H.; Belytschko, T. A comparative study on finite element methods for dynamic fracture. Comput. Mech. 2008, 42, 239–250. [Google Scholar] [CrossRef]
  4. Pezzotta, M.; Zhang, Z.L.; Jensen, M.; Grande, T.; Einarsrud, M.A. Cohesive zone modeling of grain boundary microcracking induced by thermal anisotropy in titanium diboride ceramics. Comput. Mater. Sci. 2008, 43, 440–449. [Google Scholar] [CrossRef]
  5. Krishnan, K.; Sockalingam, S.; Bansal, S.; Rajan, S.D. Numerical simulation of ceramic composite armor subjected to ballistic impact. Compos. B Eng. 2010, 41, 583–593. [Google Scholar] [CrossRef]
  6. Silling, S.A. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids. 2000, 48, 175–209. [Google Scholar] [CrossRef]
  7. Silling, S.A. Dynamic Fracture Modeling with a Meshfree Peridynamic Code. In Computational Fluid and Solid Mechanics; Elsevier Science Ltd.: Amsterdam, The Netherlands, 2003; pp. 641–644. [Google Scholar]
  8. Silling, S.A.; Askari, E. A meshfree method based on the Peridynamic model of solid mechanics. Comput. Struct. 2005, 83, 1526–1535. [Google Scholar] [CrossRef]
  9. Silling, S.A.; Epton, M.; Weckner, O.; Xu, J.; Askari, E. Peridynamic states and constitutive modeling. J. Elast. 2007, 88, 151–184. [Google Scholar] [CrossRef]
  10. Ha, Y.D.; Bobaru, F. Studies of dynamic crack propagation and crack branching with Peridynamics. Int. J. Fract. 2010, 162, 229–244. [Google Scholar] [CrossRef]
  11. Ha, Y.D.; Bobaru, F. Characteristics of dynamic brittle fracture captured with Peridynamics. Eng. Fract. Mech. 2011, 78, 1156–1168. [Google Scholar] [CrossRef]
  12. Zhu, Q.Z.; Ni, T. Peridynamic formulations enriched with bond rotation effects. Int. J. Eng. Sci. 2017, 121, 118–129. [Google Scholar] [CrossRef]
  13. Chu, B.F.; Liu, Q.W.; Liu, L.S.; Lai, X.; Mei, H. A rate-dependent Peridynamic model for the dynamic behavior of ceramic materials. CMES-Comp. Model. Eng. Sci. 2020, 124, 151–178. [Google Scholar] [CrossRef]
  14. Liu, Y.X.; Liu, L.S.; Mei, H.; Liu, Q.W.; Lai, X. A modified rate-dependent Peridynamic model with rotation effect for dynamic mechanical behavior of ceramic materials. Comput. Meth. Appl. Mech. Eng. 2022, 388, 114246. [Google Scholar] [CrossRef]
  15. Kudryavtsev, O.A.; Sapozhnikov, S.B. Numerical simulations of ceramic target subjected to ballistic impact using combined DEM/FEM approach. Int. J. Mech. Sci. 2016, 114, 60–70. [Google Scholar] [CrossRef]
  16. Zhang, Y.; Qiao, P.Z. A new bond failure criterion for ordinary state-based Peridynamic mode II fracture analysis. Int. J. Fract. 2019, 215, 105–128. [Google Scholar] [CrossRef]
  17. Foster, J.T.; Silling, S.A.; Chen, W.W. Viscoplasticity using Peridynamics. Int. J. Numer. Methods Eng. 2010, 81, 1242–1258. [Google Scholar] [CrossRef]
  18. O’Grady, J.; Foster, J. Peridynamic beams: A non-ordinary, state-based model. Int. J. Solids Struct. 2014, 51, 3177–3183. [Google Scholar] [CrossRef]
  19. Lai, X.; Ren, B.; Fan, H.F.; Li, S.F.; Wu, C.T.; Regueiro, R.A.; Liu, L.S. Peridynamics simulations of geomaterial fragmentation by impulse loads. Int. J. Numer. Anal. Methods Geomech. 2015, 39, 1304–1330. [Google Scholar] [CrossRef]
  20. Lai, X.; Liu, L.S.; Li, S.F.; Zeleke, M.; Liu, Q.W.; Wang, Z. A non-ordinary state-based Peridynamics modeling of fractures in quasi-brittle materials. Int. J. Impact Eng. 2018, 111, 130–146. [Google Scholar] [CrossRef]
  21. Wu, L.W.; Huang, D.; Xu, Y.P.; Wang, L. A non-ordinary state-based Peridynamic formulation for failure of concrete subjected to impacting loads. ES-Comp. Model. Eng. Sci. 2019, 118, 561–581. [Google Scholar] [CrossRef]
  22. Wang, H.; Xu, Y.P.; Huang, D. A non-ordinary state-based Peridynamic formulation for thermo-visco-plastic deformation and impact fracture. Int. J. Mech. Sci. 2019, 159, 336–344. [Google Scholar] [CrossRef]
  23. Zhu, F.; Zhao, J.D. Peridynamic modelling of blasting induced rock fractures. J. Mech. Phys. Solids. 2021, 153, 104469. [Google Scholar] [CrossRef]
  24. Yang, S.Y.; Gu, X.; Zhang, Q.; Xia, X.Z. Bond-associated non-ordinary state-based Peridynamic model for multiple spalling simulation of concrete. Acta Mech. Sin. 2021, 37, 1104–1135. [Google Scholar] [CrossRef]
  25. Li, S.; Lai, X.; Liu, L.S. Peridynamic Modeling of Brittle Fracture in Mindlin-Reissner Shell Theory. CMES-Comp. Model. Eng. Sci. 2022, 131, 715–746. [Google Scholar] [CrossRef]
  26. Littlewood, D.J. A Nonlocal Approach to Modeling Crack Nucleation in AA 7075-T651. In Proceedings of the ASME International Mechanical Engineering Congress and Exposition, Sandia National Laboratory, Albuquerque, NM, USA, 10–18 November 2011. [Google Scholar]
  27. Silling, S.A. Stability of Peridynamic correspondence material models and their particle discretizations. Comput. Meth. Appl. Mech. Eng. 2017, 322, 42–57. [Google Scholar] [CrossRef]
  28. Li, P.; Hao, Z.M.; Zhen, W.Q. A stabilized non-ordinary state-based Peridynamic model. Comput. Meth. Appl. Mech. Eng. 2018, 339, 262–280. [Google Scholar] [CrossRef]
  29. Wan, J.; Chen, Z.; Chu, X.H.; Liu, H. Improved method for zero-energy mode suppression in Peridynamic correspondence model. Acta Mech. Sin. 2019, 35, 1021–1032. [Google Scholar] [CrossRef]
  30. Yaghoobi, A.; Chorzepa, M.G. Higher-order approximation to suppress the zero-energy mode in non-ordinary state-based Peridynamics. Comput. Struct. 2017, 188, 63–79. [Google Scholar] [CrossRef]
  31. Luo, J.; Sundararaghavan, V. Stress-point method for stabilizing zero-energy modes in non-ordinary state-based Peridynamics. Int. J. Solids Struct. 2018, 150, 197–207. [Google Scholar] [CrossRef]
  32. Cui, H.; Li, C.G.; Zheng, H. A higher-order stress point method for non-ordinary state-based Peridynamics. Eng. Anal. Bound. Elem. 2020, 117, 104–118. [Google Scholar] [CrossRef]
  33. Chen, H.L. Bond-associated deformation gradients for peridynamic correspondence model. Mech. Res. Commun. 2018, 90, 34–41. [Google Scholar] [CrossRef]
  34. Chen, H.L.; Spencer, B.W. Peridynamic bond-associated correspondence model: Stability and convergence properties. Int. J. Numer. Methods Eng. 2019, 117, 713–727. [Google Scholar] [CrossRef]
  35. Gu, X.; Zhang, Q.; Madenci, E.; Xia, X.Z. Possible causes of numerical oscillations in non-ordinary state-based Peridynamics and a bond-associated higher-order stabilized model. Comput. Meth. Appl. Mech. Eng. 2019, 357, 112592. [Google Scholar] [CrossRef]
  36. Li, P.; Hao, Z.M.; Yu, S.R.; Zhen, W.Q. Implicit implementation of the stabilized non-ordinary state-based Peridynamic model. Int. J. Numer. Methods Eng. 2020, 121, 571–587. [Google Scholar] [CrossRef]
  37. May, S.; Vignollet, J.; De Borst, R. A numerical assessment of phase-field models for brittle and cohesive fracture: Γ-convergence and stress oscillations. Eur. J. Mech. A-Solids 2015, 52, 72–84. [Google Scholar] [CrossRef]
  38. Wu, J.Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. J. Mech. Phys. Solids 2017, 103, 72–99. [Google Scholar] [CrossRef]
  39. Narayan, S.; Anand, L. A gradient-damage theory for fracture of quasi-brittle materials. J. Mech. Phys. Solids 2019, 129, 119–146. [Google Scholar] [CrossRef]
  40. Costin, L.S. A microcrack model for the deformation and failure of brittle rock. J. Geophys. Res. 1983, 88, 9485–9492. [Google Scholar] [CrossRef]
  41. Addessio, F.L.; Johnson, J.N. A constitutive model for the dynamic response of brittle materials. J. Appl. Phys. 1990, 67, 3275–3286. [Google Scholar] [CrossRef]
  42. Rajendran, A.M.; Grove, D.J. Modeling the shock response of silicon carbide, boron carbide and titanium diboride. Int. J. Impact Eng. 1996, 18, 611–631. [Google Scholar] [CrossRef]
  43. Zuo, Q.H.; Addessio, F.L.; Dienes, J.K.; Lewis, M.W. A rate-dependent damage model for brittle materials based on the dominant crack. Int. J. Solids Struct. 2006, 43, 3350–3380. [Google Scholar] [CrossRef]
  44. Zuo, Q.H.; Disilvestro, D.; Richter, J.D. A crack-mechanics based model for damage and plasticity of brittle materials under dynamic loading. Int. J. Solids Struct. 2010, 47, 2790–2798. [Google Scholar] [CrossRef]
  45. Johnson, G.R.; Holmquist, T.J. A Computational Constitutive Model for Brittle Materials Subjected to Large Strains, High Strain Rates and High Pressures. Shock Wave and High-Strain-Rate Phenomena in Materials; Marcel Dekker Inc.: New York, NY, USA, 1992; pp. 1075–1081. [Google Scholar]
  46. Johnson, G.R.; Holmquist, T.J. An Improved Computational Constitutive Model for Brittle Materials. In AIP Conference Proceedings; American Institute of Physics: College Park, MD, USA, 1994. [Google Scholar]
  47. Johnson, G.R.; Cook, W.H. A Constitutive Model and Data for Materials Subjected to Large Strains, High Strain Rates, and High Temperatures. In Proceedings of the Seventh International Symposium on Ballistics, The Hague, The Netherlands, 19–21 April 1983. [Google Scholar]
  48. Johnson, G.R.; Holmquist, T.J.; Beissel, S.R. Response of aluminum nitride (including a phase change) to large strains, high strain rates, and high pressures. J. Appl. Phys. 2003, 94, 1639–1646. [Google Scholar] [CrossRef]
  49. Rubinstein, R.; Atluri, S.N. Objectivity of incremental constitutive relations over finite time steps in computational finite deformation analyses. Comput. Meth. Appl. Mech. Eng. 1983, 36, 277–290. [Google Scholar] [CrossRef]
  50. Monaghan, J.J.; Gingold, R.A. Shock simulation by the particle method SPH. J. Comput. Phys. 1983, 52, 374–389. [Google Scholar] [CrossRef]
  51. Li, S. Smoothed particle hydrodynamics-a meshfree method, by G.R. Liu and M.B. Liu. Comput. Mech. 2004, 33, 491. [Google Scholar] [CrossRef]
  52. Belytschko, T.; Lin, J.I. A three-dimensional impact-penetration algorithm with erosion. Int. J. Impact Eng. 1987, 5, 111–127. [Google Scholar] [CrossRef]
  53. Heinstein, M.W.; Attaway, S.W.; Swegle, J.W.; Mello, F.J. A General-Purpose Contact Detection Algorithm for Nonlinear Structural Analysis Codes; Sandia National Labs: Albuquerque, NM, USA, 1993.
  54. Whirley, R.G.; Engelmann, B.E. Automatic contact algorithm in DYNA3D for crashworthiness and impact problems. Nucl. Eng. Des. 1994, 150, 225–233. [Google Scholar] [CrossRef]
  55. Bourago, N.G.; Kukudzhanov, V.N. A review of contact algorithms. Mech. Solids 2005, 40, 35–71. [Google Scholar]
  56. Cronin, D.S.; Bui, K.; Kaufmann, C.; McIntosh, G.; Berstad, T.; Cronin, D. Implementation and Validation of the Johnson-Holmquist Ceramic Material Model in LS-Dyna. In Proceedings of the 4th European LS-DYNA User’s Conference, Ulm, Germany, 22–23 May 2003. [Google Scholar]
  57. Strassburger, E. Visualization of impact damage in ceramics using the edge-on impact technique. Int. J. Appl. Ceram. Technol. 2004, 1, 235–242. [Google Scholar] [CrossRef]
  58. Simons, E.C.; Weerheijm, J.; Toussaint, G.; Sluys, L.J. An experimental and numerical investigation of sphere impact on alumina ceramic. Int. J. Impact Eng. 2020, 145, 103670. [Google Scholar] [CrossRef]
Figure 1. Definitions of different horizons. (a) Point-associated horizons; (b) bond-associated horizons.
Figure 1. Definitions of different horizons. (a) Point-associated horizons; (b) bond-associated horizons.
Materials 16 04050 g001
Figure 2. Force vector state for PD models.
Figure 2. Force vector state for PD models.
Materials 16 04050 g002
Figure 3. Schematic diagram of the relative rotation angle [14].
Figure 3. Schematic diagram of the relative rotation angle [14].
Materials 16 04050 g003
Figure 4. Schematic diagram of the contact surfaces under normal impact [20].
Figure 4. Schematic diagram of the contact surfaces under normal impact [20].
Materials 16 04050 g004
Figure 5. Schematic diagram of uniaxial tension for an elastic bar under stress loading. (a) Computational model. (b) Discretization.
Figure 5. Schematic diagram of uniaxial tension for an elastic bar under stress loading. (a) Computational model. (b) Discretization.
Materials 16 04050 g005
Figure 6. Comparison of axial displacement contours with m = 1.5 .
Figure 6. Comparison of axial displacement contours with m = 1.5 .
Materials 16 04050 g006
Figure 7. Comparison of axial displacement contours with m = 2.0 .
Figure 7. Comparison of axial displacement contours with m = 2.0 .
Materials 16 04050 g007
Figure 8. Comparison of axial displacement contours with m = 3.0 .
Figure 8. Comparison of axial displacement contours with m = 3.0 .
Materials 16 04050 g008
Figure 9. Comparison of axial displacement contours with m = 4.0 .
Figure 9. Comparison of axial displacement contours with m = 4.0 .
Materials 16 04050 g009
Figure 10. Comparison of axial displacement contours with m = 5.0 .
Figure 10. Comparison of axial displacement contours with m = 5.0 .
Materials 16 04050 g010
Figure 11. Schematic diagram of uniaxial tension of an elastic bar under velocity loading. (a) Computational model. (b) Discretization.
Figure 11. Schematic diagram of uniaxial tension of an elastic bar under velocity loading. (a) Computational model. (b) Discretization.
Materials 16 04050 g011
Figure 12. Convergence curves for the system kinetic energy with different solutions.
Figure 12. Convergence curves for the system kinetic energy with different solutions.
Materials 16 04050 g012
Figure 13. Computational model of ceramics under edge-on impact.
Figure 13. Computational model of ceramics under edge-on impact.
Materials 16 04050 g013
Figure 14. Comparison of fracture patterns for the Al2O3 ceramic plate under edge-on impact when t = 9.7   μ s (a) Experimental [57]; (b) PD simulation.
Figure 14. Comparison of fracture patterns for the Al2O3 ceramic plate under edge-on impact when t = 9.7   μ s (a) Experimental [57]; (b) PD simulation.
Materials 16 04050 g014
Figure 15. Damage evolution patterns of the Al2O3 ceramic plate under edge-on impact.
Figure 15. Damage evolution patterns of the Al2O3 ceramic plate under edge-on impact.
Materials 16 04050 g015
Figure 16. Computational model of ceramics under normal impact.
Figure 16. Computational model of ceramics under normal impact.
Materials 16 04050 g016
Figure 17. Comparison of fracture patterns for the Al2O3 ceramic plate under normal impact. (a) Experimental [58]; (b) PD simulation.
Figure 17. Comparison of fracture patterns for the Al2O3 ceramic plate under normal impact. (a) Experimental [58]; (b) PD simulation.
Materials 16 04050 g017
Figure 18. Damage evolution patterns of the Al2O3 ceramic plate under normal impact.
Figure 18. Damage evolution patterns of the Al2O3 ceramic plate under normal impact.
Materials 16 04050 g018
Table 1. JH2 constitutive model parameters of the Al2O3 ceramic [56].
Table 1. JH2 constitutive model parameters of the Al2O3 ceramic [56].
ParameterValueParameterValue
Density ρ (kg/m3)3700T (GPa)0.2
Shear modulus G (GPa)90.16 P HEL (GPa)1.46
Poisson’s ratio ν 0.22 σ HEL (GPa)2.0
A0.93 D 1 0.005
B0.31 D 2 1.0
C0.0 K 1 (GPa)130.95
M0.6 K 2 (GPa)0.0
N0.6 K 3 (GPa)0.0
ε ˙ 0 1.0 β 1.0
Table 2. Solutions to modify the original BA-NOSB PD model.
Table 2. Solutions to modify the original BA-NOSB PD model.
Solution ASolution BSolution C
μ = μ p = ρ / ρ 0 1 μ = μ b = 1 / J b 1 μ = μ b = 1 / J b 1
Y ¯ B ξ = ξ Y ¯ B ξ = ξ Y ¯ B ξ = ξ + η B
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, J.; Liu, Y.; Lai, X.; Liu, L.; Mei, H.; Liu, X. A Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model for Impact Problems of Quasi-Brittle Materials. Materials 2023, 16, 4050. https://doi.org/10.3390/ma16114050

AMA Style

Zhang J, Liu Y, Lai X, Liu L, Mei H, Liu X. A Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model for Impact Problems of Quasi-Brittle Materials. Materials. 2023; 16(11):4050. https://doi.org/10.3390/ma16114050

Chicago/Turabian Style

Zhang, Jing, Yaxun Liu, Xin Lai, Lisheng Liu, Hai Mei, and Xiang Liu. 2023. "A Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model for Impact Problems of Quasi-Brittle Materials" Materials 16, no. 11: 4050. https://doi.org/10.3390/ma16114050

APA Style

Zhang, J., Liu, Y., Lai, X., Liu, L., Mei, H., & Liu, X. (2023). A Modified Bond-Associated Non-Ordinary State-Based Peridynamic Model for Impact Problems of Quasi-Brittle Materials. Materials, 16(11), 4050. https://doi.org/10.3390/ma16114050

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop