Hybrid PD-DEM approach for modeling surface erosion by particles impact

Peridynamics (PD) theory is a promising technique for modeling solids with discontinuities. Short-range repulsive force models are commonly employed in PD impact event simulations. Despite their extensive usage, short-range force models do not take damping, friction, and tangential force components into account and hence are unable to effectively describe energy dissipation, leading to uncertainty in the calculation of contact forces. However, the accuracy of impact simulations using alternate contact models has not been extensively investigated in the context of PD impact simulations. The Discrete Element Method (DEM) has been proven to be the most reliable and effective approach to model collision processes between distinct solid objects. This work presents, a particle-based hybrid PD-DEM model to accurately predict the particle impact forces and the resulting damage to the target material. The present model brings together the unique capabilities of PD and DEM and has the potential to make use of the various DEM contact laws, which allow the development and adjustment of relevant contact forces in the normal and tangential directions. Furthermore, damping effects, friction, and intra-particle stiffness are incorporated into the simulations through DEM. The proposed method has been used for modeling material failure after being validated and verified for the contact parameters during the impact process. The predicted damage patterns and resulting material loss demonstrate good agreement with the experimental results reported in the literature.


Introduction
Solid particle erosion (SPE) of engineered materials is a dynamic process that causes material loss from a target surface due to the impingement of fast-moving solid particles [1,2].Surface erosion causes component wear, surface roughening, surface degradation, and a shortening in the functional life of the structure.Metals and polymer composite materials are increasingly being employed in a range of industries where they may be subjected to solid particle erosion.
Another important occurrence of SPE is in wind turbine blades due to solid particle impacts such as hailstones, sand, or dust particles.The blades of modern turbines are naturally subjected to high-speed winds as they travel at speeds of 80 ms −1 or higher, and are bombarded by solid (e.g.hailstone or sand) particles depending on the installation site [9].The particles can cause abrasive wear or impact erosion depending on the size and impact angle.As the coating of the blade is compromised, synergy with UV degradation speeds up the erosion rate leading to widespread damage to internal structures.Initially, the increased surface roughness increases the friction drag resulting in an earlier onset of the stall which causes a substantial aerodynamic performance penalty.Ultimately, as the severity of erosion increases, the blades' structural integrity is compromised resulting in unplanned downtime and high maintenance costs [10].Therefore, it is very crucial to investigate the solid particle erosive impact and behavior of polymeric-engineered materials under such operating conditions.
The solid particle erosion process is a complex phenomenon involving several parameters that interact simultaneously, which includes mechanical and chemical properties of the materials as well as the shape, size, velocity, and impingement angle of the impactor.Numerous experimental investigations considering particle impact tests have been conducted by researchers to understand the erosion mechanism of polymeric and brittle materials under different conditions [11][12][13][14][15][16][17][18][19][20][21][22][23].The development of numerical models to simulate erosion events has been progressively attempted in recent decades.Among the numerical methods, the finite element method (FEM) has been widely used to simulate the erosion process of ductile and brittle materials [24][25][26].Aquaro et al. [27] developed a FEM-based approach that is suitable for determining erosion rates in a variety of applications.Wang et al. [28,29] used Johnson-Cook and Johnson-Holmquist material models and established a FEM framework to simulate the erosion process.In [30][31][32] stress-strain criteria are used inside the FEM scheme to model the impact of a rigid particle on glass material.FEM solvers commonly use Lagrange multipliers or the penaltybased technique to formulate the contact constraints [33].
Even though FEM is quite accurate in predicting stress and strain fields, it experiences issues in modeling erosion, because governing partial differential equations of FEM are undefinable at geometric discontinuities of cracked surfaces.
To solve these issues in crack growth problems, Xu et al. [34] introduced the cohesive zone method (CZM) in FEM that defines kinetic relations between the elements.The fracture process, on the other hand, is a dynamic phenomenon that involves crack propagation and requires continuous geometry re-meshing, which is extremely laborious and computationally expensive.Another popular technique for simulating cohesive fractures is the boundary element method (BEM) [35,36].The BEM while using a boundary-only mesh minimizes the problem's size and compared to FEM reduces dimensionality by one.To avoid the costly re-meshing, Moës et al. [37] developed the extended finite element method (XFEM), which incorporates enrichment functions based on the partition of unity in the finite element formulation.Although XFEM can effectively predict damage patterns and allow fractures to travel through an element, it can only simulate a limited number of cracks simultaneously.Due to complex algorithms, it is still challenging for XFEM to continuously track the evolving large number of cracks and their interactions e.g.scattering of glass.
The meshfree or particle-based methods such as the discrete element method (DEM) [38], smoothed particle hydrodynamics (SPH) [39], reproducing kernel particle method (RKPM) [40], material point method (MPM) [41], meshfree Galerkin method (MFGM) or element free Galerkin method (EFGM) [42,43] provide an alternative approach for handling the inherent problem of elements and have been used for a broad range of applications.Due to their smooth and higher-order interpolation, it is quite simple to capture massive deformations in the system.Additionally, tracking algorithms or re-meshing are not required in mesh-free approaches [44].The contact methods in most of the meshfree approaches have limitations and most of them are extensions of existing methods developed for FEM [45].Meshless methods are discrete by definition, yet they are developed using local continuum procedures.The spatial partial derivatives in governing equations of meshless methods are not defined for conditions such as crack initiation and propagation and assume that the body remains continuous as it deforms.
To address this issue, Peridynamics (PD) was introduced as a non-local continuum mechanics theory by Silling [46].The advantage of this formulation is that the governing equations remain valid whether the structure is continuous or not.The PD theory is a reformulated version of classical continuum mechanics (CCM) and it is more suitable for modeling solids with discontinuities, such as cracks.In PD, the partial differential equations of CCM are replaced with integral equations.In non-local continuum models which are based on PD theory, the state of a material point is influenced by other material points that are situated within a certain region with a finite radius called a horizon [47].The PD theory incorporates damage into the material response and permits damage initiation and propagation at many spots along arbitrary internal channels without using any specific crack growth criteria.The PD theory has been effectively applied to material deformation and damage prediction applications.Silling et al. [48] successfully applied PD theory to model damage in concrete due to numerous impacts and glass plate crack fragmentation.Madenci et al. [49] developed a state-based PD model for plastic deformation.Oterkus et al. [50] applied PD to impact damage assessment and residual strength of a reinforced panel under compression.In PD, the contact model is generally implemented with the help of rigid impactors or short-range repulsive force algorithms [51].The short-range force approach simulates contact interactions between bodies using a technique similar to molecular dynamics.It applies spring-like repulsive forces between closely spaced nodes at each simulation step, yielding satisfactory outcomes in several scenarios.Another option is to use a conventional contact algorithm within a PD simulation, as demonstrated in [52].The short-range force model can be modified to include friction effects [53] if required.Additionally, to enhance its accuracy, the model may be extended to incorporate PD bonds, which it currently does not account for in its basic form.The absence of consideration for PD bonds leads to a situation where material points interact simultaneously via the material and contact models.For further details on correspondence constitutive models, particularly in relation to material interpenetration, readers may refer to the study conducted by Tupek and Radovitzky [54].Although these contact models are used for extremely complex impact problems, there is still uncertainty about the reaction forces and torques.The penalty stiffness or short-range force constants involved in these methods are arbitrarily set in most of the investigations, or their formulations are unrealistic.It is necessary to have an accurate estimation of the contact forces in order to predict the erosive impacts of particles on the material surface.Contact modeling in PD simulations is still an open area of research.
Researchers [55][56][57][58][59] have recently proposed alternatives to the aforementioned contact models and provided a generic contact modeling approach for PD impact problems.Zhang et al. [57,60,61] established a two-dimensional PD-DEM coupled with the lattice Boltzmann method (LBM) to address particle-fluid interactions and particle collisions with solid surfaces, as well as the associated damage.Jha et al. [55] also developed a two-dimensional PD coupled with a DEM [38] framework to simulate particle impacts and compressive behavior of multi-particle systems.Madenci et al. [62] and Anicode et al. [58,59] presented a three-dimensional generalized particle-based contact modeling approach in the framework of peridynamics and analyzed the particle contact parameters, and performed investigations for the progressive damage in the target material.
The present study presents a three-dimensional particlebased hybrid PD-DEM model, which combines the strengths of PD and DEM inside the LAMMPS framework.The main reason for coupling PD with DEM [38] is to take advantage of the multiple DEM [38] contact laws which are theoretically sound and are extensively validated for various materials and impact conditions.The contact model parameters can be modified to produce the appropriate DEM contact forces, damping effects, and intra-particle stiffness.The present hybrid PD-DEM model incorporates the Hertzian force-displacement law [63] for the contact force in the normal direction and the stick-slip friction model proposed by Mindlin [64] in the tangential direction.It avoids the arbitrary selection of the penalty parameter through direct contact mechanics estimations.These contact laws are extendable [65,66] for systems that involve simultaneous multi-particle interactions.This contact modeling technique in PD enables the accurate simulation of interacting entities that deform, fracture, and are composed of numerous particles.
This paper is organized as follows: The peridynamics theory of material deformation and the governing equation of motion of the impactor and target are discussed in Sect. 2. The implementation of the contact model between the impactor and the target is described in Sect. 3 of the article.Section 4 of the paper presents numerical findings for material failure in brittle and laminated targets as well as validations of impact parameters and contact forces.

Mathematical model
We have established a novel hybrid peridynamics discrete element method (PD-DEM) based meshless software library inside the LAMMPS framework to numerically model solid particle erosion.Equations arising from the conservation of momentum and constitutive models related to material deformation and stress are solved by this coupled PD-DEM.The characteristics of the target material are described by the integrodifferential equations of the PD theory, where DEM explicitly depicts the motion and interaction of discrete solid particles.

Peridynamics formulation
The peridynamics theory uses integral equations instead of the partial differential equations based classical formulation of continuum mechanics to describe the relative displacements and non-local exchange of constitutive information by applying forces between material points across finite distances [46], allowing for the natural formation of discontinuities and cracks within continuous materials.This work concentrates on the coupling of PD with DEM and a brief description of bond-based PD theory for brittle material is presented here, for more details readers are referred to Silling et al. [46,67] and Madenci et al. [47].In a bond-based PD [46], a material point x, interacts with another material point x within the interaction domain H x , as shown in Fig. 1.The interaction domain H x of the material point x is assumed to be a spherical region specified by a radius δ which is known as its horizon.Material points within the interaction domain H x of material points x are called the family members of x.
According to Silling et al. [46], the interactions of a material point x with another material point x inside the interaction domain H x are governed by the PD equation of motion as where material points are represented by spherical PD particles with diameter d m .In Eq. ( 1), ρ m is the density of the PD particle, and d is the displacement vector of a particle located at a point x at a time τ .The time derivative of the displacement vector d of each particle is correlated with the integral of an internal force field f(η, ξ ) and an external body force F b .The peridynamic force exerted on the PD particle located at the point x by all the PD particles within H x is expressed as the integral of a force density f(η, ξ ) over the volume V x , where ξ x − x and η d − d are the relative position and displacement vectors respectively.The force density f(η, ξ ), also termed as the response function representing interparticle bonds, is expressed as The bond constant c, also known as the micro-modulus function is a PD parameter obtained by equating strain energy densities from the classical theory of elasticity with peridynamics under simple loading conditions [58,59] and is given as where E is Young's modulus and ν is the Poisson's ratio of the material.In Eq. ( 2), n is a unit vector directed from x + d to x + d , and the parameter s represents bond stretch, expressed as If the bond stretch s value exceeds its critical value s c , then the bond breaks, which is an irreversible process.The critical stretch s c in 3D for bond-based peridynamics was derived by Silling et al. [51] as where G c is the fracture energy per unit area of the material.The parameter μ in Eq. ( 2) is a history-dependent function associated with material damage and its value is 0 or 1.If s ≤ 0, for live or broken bonds μ 1, otherwise, μ 0. Figure 2 shows the elastic and perfectly plastic constitutive model for the bonds.The force density function can be nonzero in both compression and tension.For elastic and plastic regions, the force density relationship can be written as The weighted ratio of the number of broken bonds to the total number of initial bonds between a material point and its family members is used to characterize local damage at a point.The local damage at a material point ranges from 0 to 1 and, according to Silling et al. [51], it can be quantified as when φ 1, it denotes a completely damaged point, all the bonds initially associated with the material point have been eliminated, and φ 0 indicates an undamaged material point i.e. all interactions are intact.

Governing equations for solid DEM particle
The DEM model is used to simulate the motion of solid particles.In DEM, both the translational and rotational motion of solid particles are governed by Newton-Euler equations of rigid body dynamics, i.e. if u j and ω j are the translational and angular velocities of the j th particle respectively, then the particle must satisfy the following equations where m j and I j are the mass and moment of inertia of j th particle, respectively.F j is the resultant force, and T j is the resultant torque acting on the particle j about the axis passing through its center.Since all of the forces and torques acting on the j th particle are added together in the vectors F j and T j , the sum can be expressed as where F c j is the contact force due to interaction of the particle j with other particles and obstacles, F ext j and T ext j are external load, F damp j and T damp j are force and torque because of damping in the system.q c j is the torque other than due to a tangential force e.g., due to rolling motion or torsion, r c j is vector connecting particle center with the contact point and n c is the total number of particles in contact with particle j.

Contact model
The primary purpose behind coupling DEM with PD is to use the DEM contact laws.The contact method used here is similar to that extensively used in the context of DEM, it is based on Hertz's theory [63] for the force-displacement relations in the normal direction, and for the tangential direction, it employs the no-slip elastic solutions for force-displacement relations proposed by Mindlin [64].It is a penalty force contact approach that prevents interpenetration between contacting particles by using a spring-damper system.The present contact model additionally includes the contact point concept and coefficient of restitution for a nonlinear spring and damper system.Here, we briefly discuss the multi-particle contact formulation; readers are referred to Anicode et al. [58,59] and Bui et al. [66] for a more detailed discussion on the formulation and evolution of the multi-particle contact forces.
Let r i be the radius of the rigid spherical particles representing the PD material points and r j be the radius of the impactor (in case of single DEM particle) or radius of subvolume particles that are combined (in case of multi-particle DEM approach) to model the solid particle e.g.sand particle.If x i and x j are position vectors of particle i and j respectively, then the contact force between the two particles based on the normal overlap D n is calculated as where, n (i) is the unit vector in the normal direction.The tangential overlap D t is obtained by using relative tangential velocity Ḋt v ( j) (i) − Ḋn over time τ as According to Bui et al. [66], the force-displacement relations in the normal direction can be extended for multiparticle interactions.The normal contact force F n on the particle i due to its interaction with neighboring particles N k becomes where K n and C n are stiffness and damping constants in the normal direction respectively and can be expressed as ϑ n is the damping ratio in the normal direction.E, R and m denotes the effective Young's modulus, radius, and mass of the two interacting particles i and j, and can be obtained as 1 m The tangential component of contact force F t on the particle i due to its interaction with neighboring particles N k is given as where K t and C t are the stiffness and damping constants in the tangential direction respectively and can be defined as ϑ t is the damping ratio in the tangential direction and G is the effective shear modulus of the two interacting particles i and j, and can be obtained as A Coulomb friction coefficient λ is used to model a stick and slip behavior [68] and F t of two interacting particles i and j is set as The external body force as in Eq. (1) acting on PD particle i becomes while the contact force on the impactor F c j in Eqs.(10) and becomes 4 Results and discussion

Validation case for normal impact
Numerous numerical studies that appropriately considered the contact forces in impact event modeling are available in the literature [58,59,[69][70][71][72].To verify and validate the proposed hybrid PD-DEM results for contact forces and impact parameters by using the previously published results [58,59,69,70], we first simulated a classical impact problem of a square steel plate subjected to a normal impact of a rigid sphere.A spherical steel particle of radius R 0.01 m  206 GPa, density ρ 7833 kg/m 3 , and Poisson's ratio v 0.28.The horizon is three times the size of PD particles and the target plate is discretized with PD particles of radius r 0.001 m as shown in Fig. 4.
The results obtained by using the present hybrid PD-DEM are plotted in Fig. 3 and the comparison of the obtained data with the results from published studies is presented in Table 1. Figure 3a, b, and c shows the time histories of normal velocity, normal displacement, and normal contact force of the impactor respectively.The predicted PD-DEM results in Table 1 are in good agreement with the literature [58,69] results, where the results obtained by using the Fig. 4 Displacement of the plate in the normal direction after impact using PD-DEM Fig. 5 Velocity of the plate in the normal direction after impact using PD-DEM original PD short-range repulsive force model (PD-SRRFM) show a clear difference from the other results.The lack of damping effects, together with the absence of friction, and intra-particle stiffness parameters are the main factors of the unrealistic PD-SRRFM findings.Figures 4 and 5 show the displacement and velocity fields of the plate in the normal Fig. 7 Force fields on the soda-lime glass wall due to the impact of spherical aluminum particle Fig. 8 Stress fields on the soda-lime glass wall due to the impact of spherical aluminum particle direction after the impacting rigid particle rebounds respectively.

Validation case for oblique impact
In this validation test case, a soda-lime glass wall is subjected to an oblique impact of a spherical solid particle made up of aluminum oxide having a radius R 5 mm.The oblique  impacting particle has a constant velocity 3.9 ms −1 and it targets the horizontal soda-lime glass wall at an angle α 85 • .The characteristics and mechanical properties of both materials are listed in Table 2.A single DEM particle represents the impactor, where the wall is geometrically defined by length l 0.04 m, width w 0.04 m, and thickness h 0.004 m and is discretized with PD particles of radius r 0.0005 m and the horizon is three times the size of PD particles.A nondimensional graphical comparison of the present PD-DEM and PD-SRRFM results with the reference literature [70] for the contact forces in the normal and tangential directions is presented in Fig. 6a and b respectively.The comparison of PD-DEM results in Fig. 6a and b shows similarity with the literature results, while the results predicted using PD-SRRFM totally differ from the reference findings.Figures 7 and 8 show the force and stress fields on the soda-lime glass wall due to the impact of spherical aluminum particle.
In validation cases of Sects.4.1 and 4.2, a good agreement with the literature has been found for the normal and tangential contact forces acting on an impactor during the contact.These comparisons of obtained results confirm that the contact model is implemented correctly in the coupled PD-DEM simulations and provide us the confidence to further model the solid particle erosion with the current hybrid method.

Erosive impact of particle on a brittle plate
In this case, damage propagation is modeled in a brittle plate due to the normal impact of a solid particle.The target plate is cylindrical in shape with diameter d 0.06 m and thickness h 0.03 m, where it is discretized with PD particles of radius r 0.00025 m resulting in a total of 688,263 particles and the horizon is three times the size of PD particles.The spherical solid particle of radius R 0.005 m described by DEM is projected normally at the center of the top flat surface of the cylinder with velocity V z −35 ms −1 and V z −70 ms −1 .The maximum time step size t 1.0 × 10 −8 s and total time t 200μs.The time step is a crucial parameter in these simulations, as its size impacts both the accuracy and stability of our results.To optimize our simulations, we implement an adaptive time stepping approach, utilizing a maximum time step of 1.0×10 −8 s and a minimum time step of 1.0×10 −10 s Fig. 9 PD discretization of the brittle plate and the solid spherical impactor for most simulations.This approach allows us to dynamically adjust the time step according to the system's complexity, ensuring both high precision and computational speed.The solid spherical impactor and the PD discretization of the brittle cylindrical plate are illustrated in Fig. 9.The material properties of the solid impactor and the brittle plate are listed in Table 3 and their values are considered similar to the ones used in Anicode et al. [59].Figures 10 and 11 show the progression of damage through the plate thickness with time for impact velocities of 35 ms −1 and 70 ms −1 , respectively.From Figs. 10 and 11, it is noticed that the damage is more localized at high velocity.Figure 12a and b show the three-dimensional view of the damaged brittle plates after the particle impact at the incident velocities of 35 ms −1 and 70 ms −1 respectively.Moreover, Figs. 13 and 14 illustrate the comparison of conical crack patterns predicted by the PD-DEM approach with the experimental results [73] at an impact velocity of 160 ms −1 and 310 ms −1 .The comparison in Figs. 13 and 14 closely resembles the experimental findings [73], confirming that the present model is a reliable technique for modeling solid particle erosion.Figure 15 shows the 3D conical crack patterns after removing the material points with damage index < 0.5 at impact velocities of 160 ms −1 and 310 ms −1 .We used the prototype micro elastic brittle (PMB) PD material model and the interested reader is encouraged to refer to articles [46,51,74] for better understand the validity criteria for the bond-based model.

Erosive impact of particle on a composite laminated plate
In this case, damage propagation is modeled in a composite laminated plate due to the normal impact of a solid particle.The target plate is cylindrical in shape with diameter d 0.06 m and thickness h 0.015 m, and is comprising of four plies.The thickness of each of the top two layers is 0.0025 m and each of the bottom two layers are 0.005 m thick.The plate is discretized with PD particles of radius r 0.00025 m resulting in a total of 372,340 PD particles and the horizon is three times the size of PD particles.A spherical solid particle of radius R 0.008 m described by DEM is projected normally at the center of the top flat surface of the cylinder with impact velocity V z −100 ms −1 .The time step size t 1.0 × 10 −8 s and total time t 200μs.Figure 16 displays the discretized multi-layer laminated plate cross-sectional view of thickness and the different colors of particles represents different layers.Here we consider two different laminated plates, first a simple laminated plate, and in the second case a hybrid laminated plate.The purpose is to test our numerical implementation and analyze the material response in both cases.In the first case, the same material properties are used for each layer of the plate.The layer structure has no difference from a single-layer plate when using the same material for each layer.Since the interaction between two layers is modeled by applying bond-based PD potential to particles at the interface which depends on material properties.Hence, with the same material, the PD potential of particles at the interface is similar to interior particles.The material parameters of the solid impactor and the laminated plate are given in Table 4. Figure 17 shows the cross-sectional images of the damage progression through the laminated plate thickness with time for an impact velocity of 100 ms −1 , when all of the plate layers are made up of the same material.Next, we consider a hybrid laminated plate, such as each layer has different material properties.The interaction between any two layers of the hybrid laminated plate is modeled by applying bond-based PD potential to particles at the interface obtained by taking average values of material properties.The material parameters of the solid impactor and the hybrid laminated plate are listed in Table 5. Figure 18 shows the progression of damage through the thickness of the hybrid laminated plate with the time when impact velocity is 100 ms −1 .The ultimate damage to the two plates along thickness is compared in Fig. 19, and it can be seen that the hybrid laminated plate has somewhat deeper and broader damage than the simple laminated plate.Figure 20 compares the damage to the laminated plates after the top layer has been removed, again it is noticed that the hybrid laminated plate has sustained more damage than the simple laminated plate.The study demonstrates how the choice of material has a significant impact on the evolution of damage profiles in laminated plates.

Conclusions
This paper presents a particle-based hybrid model that combines the peridynamics theory with the DEM to simulate   the particle impact events and predict surface erosion caused by the colliding solid particles.In this coupled framework, the particle interaction with the target material is modeled based on force-displacement relations proposed by Hertz and Mindlin in the normal and tangential directions respectively.The contact model also takes into account contact friction, damping effects, and intra-particle stiffness parameters that are normally neglected in repulsive force models adopted in PD simulations.The contact model for the hybrid PD-DEM is first tested and validated for the normal and oblique impacts on the target material without allowing any damage or failure.The time histories of the effective contact parameters of the collision process, including displacement, velocity, and reaction forces acting on the impactor , exhibit a good match with the findings of the experiments reported in the literature.Moreover, numerical simulations were carried out to simulate the erosive impact of particles on brittle and laminated targets by allowing material failure and damage.The predicted conical crack patterns in brittle material caused by the impact of spherical particles striking normally at different velocities closely resemble with the experimental findings [73].Therefore, the present model is shown to be reliable for simulating solid particle erosion under different conditions, and contact parameters can be tuned to achieve the desired effects.
, y, y´Material point (position) D Displacement vector m d Time derivative of displacement (acceleration) m/s 2

Fig. 1
Fig. 1 Bond-based Peridynamics: Interaction of material point x with material point x and material point y with material point y in undeformed and deformed states respectively

Fig. 3
Fig. 3 The evolution of impactor parameters with time using PD-DEM: a normal velocity, b normal displacement, and c normal force

Fig. 6
Fig. 6 Non-dimensional contact forces: a normal force and b tangential force

Fig. 10
Fig. 10 Cross-sectional view of damage evolution along thickness of the brittle plate with time, sequentially from a to f, when impactor velocity is 35 ms −1

Fig. 11 1 Fig. 12
Fig. 11 Cross-sectional view of damage evolution along thickness of the brittle plate with time, sequentially from a to f, when impactor velocity is 70 ms −1

Fig. 13
Fig.13 Comparison of conical crack patterns with experimental results[73] when impact velocity is 160 ms−1

Fig. 14
Fig.14 Comparison of conical crack patterns with experimental results[73] when impact velocity is 310 ms−1

Fig. 15
Fig. 15 Conical crack patterns after removing the material points with damage index < 0.5 at impact velocities a 160 ms −1 and b 310 ms −1

YoungFig. 16
Fig. 16 PD discretization of the multi-layer laminated plate thickness cross-sectional view and the solid impactor

Fig. 17
Fig.17 Cross-sectional view of damage evolution along thickness of the simple laminated plate (same material for all layers) with time, sequentially from a to f, when impact velocity is 100 ms−1

Fig. 18
Fig.18 Cross-sectional view of damage evolution along thickness of the hybrid laminated plate (layers materials are different) with time, sequentially from a to f, when impact velocity is 100 ms−1

Fig. 19
Fig.19 Comparison of damage along thickness of the laminated plates a simple, b hybrid, when impact velocity is 100 ms−1

Fig. 20
Fig.20 Comparison of damage after removing top layer of the laminated plates: a simple, b hybrid, when impact velocity is 100 ms−1

Table 1
Data comparison and thickness h 0.008 m is simply supported.The interactions between material points of the target plate are represented by bond-based PD particles.The impactor and target plate are made up of steel having the same material properties i.e.Young's modulus E

Table 2
Material properties

Table 3
Material properties

Table 4
Material properties of the composite laminated plate and solid impactor

Table 5
Material properties of the hybrid laminated plate and solid impactor