Numerical Material Model for Composite Laminates in High – Velocity Impact Simulation

A numerical material model for composite laminate, was developed and integrated into the nonlinear dynamic explicit finite element programs as a material user subroutine. This model coupling nonlinear state of equation (EOS), was a macro-mechanics model, which was used to simulate the major mechanical behaviors of composite laminate under high–velocity impact conditions. The basic theoretical framework of the developed material model was introduced. An inverse flyer plate simulation was conducted, which demonstrated the advantage of the developed model in characterizing the nonlinear shock response. The developed model and its implementation were validated through a classic ballistic impact issue, i.e. projectile impacting on Kevlar29/Phenolic laminate. The failure modes and ballistic limit velocity were analyzed, and a good agreement was achieved when comparing with the analytical and experimental results. The computational capacity of this model, for Kevlar/Epoxy laminates with different architectures, i.e. plain– woven and cross–plied laminates, was further evaluated and the residual velocity curves and damage cone were accurately predicted.


INTRODUCTION
Composite materials, due to their inherently superior mechanical properties such as high-specific strength and stiffness, have been widely used in aerospace, civil, and armor protection applications (Zhang et al., 2008;Tabiei and Ivanov, 2002).To improve the ballistic performance of composite laminate, it is essential to further study its ballistic impact response and better understand its impact/penetration behaviors.Over the past decade, considerable research work has been done by many scholars from different angles or aspects.The efforts to study the ballistic behaviors of com-Latin American Journal of Solids and Structures 14 (2017) [1912][1913][1914][1915][1916][1917][1918][1919][1920][1921][1922][1923][1924][1925][1926][1927][1928][1929][1930][1931] posites mainly focus on the ballistic limit, residual velocity, energy absorption and failure modes.Three methods, i.e. experiment, theoretical analysis and numerical simulation are selected to conduct the studies.
Theoretical model can explicitly express the relationship between the ballistic performance of composites and influential parameters such as impact velocity, nose shape of projectiles, material characteristics and boundary conditions.Analytical study is a powerful tool to completely understand the damage and energy dissipation mechanisms of composites (Naik and Doshi, 2005).Analytical investigations for composite materials under low to high impact velocities have been carried out and several theoretical models have been also proposed by many researchers (Wen, 2000;Billon and Robinson, 2001;Naik and Doshi, 2005).Experimental study is another significant method because it is not only closer to the real physical process but also is the base of the other methods.Many experimental studies are reported in the literatures (Zhu et al., 1992;Wang et al., 2005;Babu et al., 2007;Zu et al. 2015), which give insight of the anti-ballistic mechanism of composites, and provide a large amount of data for the future studies.Numerous intervening parameters as mentioned above, make the impact issue of composite materials are rather complicated.Theoretical models usually depend on many assumptions and simplified conditions, which seriously hinder their utility for general engineering problems.Moreover, experiment methods are usually subjected to the technical conditions, and the internal impact process is difficult to observe.Therefore, designing composite armors or shields solely based on the theoretical models and experimental data is uneconomical and impracticable (Silva et al., 2005;Xin and Wen, 2012).
Fortunately, recent advances of computer science and numerical algorithms, offer a possibility of using numerical simulation combined with slight experiment tests for evaluating composite materials.The material model and modeling method for composite play a key role in numerical simulation.Hoof (1999) developed a numerical model for composite laminate.In this model, composite laminates were constructed ply by ply and a tie-break contact algorithm was introduced to simulate the interlaminar failure.Tabiei and Ivanov (2002) presented a micromechanical model for flexible woven fabric.The yarns' reorientation and fabric architecture were considered in this model.It could simulate the dual behaviors of the flexible fabric, i.e. the trellis mechanism behavior (before yarns locking) and generally anisotropic elastic properties (after yarns packing).Grujicic (2012) developed a fiber-level model for Kevlar ® KM2 ballistic fabric, in which yarn was represented as a bundle of collinear discrete fibers, and each fiber was discretized into a number of three-dimensional beam elements.Using this fiber-level model, the influence, of fiber transverse properties and the friction among fibers, on the ballistic behavior of fabric have been investigated.Anderson (1994) presented a constitutive relationship for anisotropic materials, in which the coupling problem of the volumetric and deviatoric had been overcome.Hayhurst (1999), Clegg (2006) and Riedel (2006) et al. implemented Anderson' approach and developed an advanced material model for Nextel and Kevlar/epoxy cloth.This model was based on smoothed-particle hydrodynamics (SPH) algorithm and used to predict the hypervelocity impact (HVI) behaviors of Nextel and Kevlar/epoxy cloth.
Composite laminate is fabricated by several single fabrics, and fabric is assembled by a number of yarns in specific patterns such as weaving or stacking.While the yarns are usually made by hundreds of the high-strength fibers (Grujicic et al., 2012).The more meticulous of the numerical model will the more reflect the microstructural of composites, and the more accuracy/fidelity of the nu-Latin American Journal of Solids and Structures 14 (2017) 1912-1931 merical results will be obtained.However, it will make the computational scale increase rapidly, and thus result in a dramatic decline in computational efficiency.On the other hand, for the engineering issues, e.g.designing of composite armors or shields, the macroscopic parameters, such as ballistic limit, residual velocity, energy absorption and penetration/perforation behavior, are the most important design references.
Aiming at such problem, a computational macro-mechanical model was developed in this paper, which was based on Anderson' methodology and had been implemented in the nonlinear dynamic explicit finite element code LS-DYNA as a user-defined material model.The organization of the paper is as follows.The theoretical framework of this model is discussed and an inverse flyer plate simulation is conducted to demonstrate the nonlinear shock response of this developed model in Section 2. To validate this model, a high-velocity penetration problem for Kevlar29/Phenolic laminate has been investigated by numerical simulation in Section 3. The computational capacity of this model for Kevlar/Epoxy laminates with different architectures is further evaluated in Section 4. The main conclusions are summarized in Section 5.

MATERIAL MODEL
As mentioned above, composite materials are generally synthesized by two or more materials through chemical method and possess themselves hierarchical structure.To enhance the engineering practicability and improve the computational efficiency, the composite material model developed in this paper is a computational macro-mechanical model, which equates the whole composite laminate into an orthotropic homogeneous material as shown in Figure 1.This model accounting for the nonlinear shock effects, aims to capture the main mechanical behaviors and evaluate the macroscopic parameters of composite laminates under high-velocity impact loading.Comparing with isotropic materials (e.g.metallic materials), the behaviors of anisotropic materials subject to impact seem to be more complex.The phenomena observed in the impact tests for composite laminates are as follows (Hayhurst et al., 1999;Silva et al., 2005): ( where ij  is the stress tensor; ij C are 9 independent elastic constants.
Pressure is defined as a third of the trace of the stresses, i.e.
  . Therefore, expanding Eq. ( 1), the expression of pressure can be obtained: Moreover, the deviatoric stress, d ij  can be expressed as:  Solids and Structures 14 (2017) 1912-1931 From the Eq. ( 2) and Eq. ( 3), it is can be seen clearly that both the volumetric and deviatoric components of strain have the contributions to the pressure and deviatoric stress.For the composite materials, a significant conclusion can be drawn that the volumetric and deviatoric response are strongly coupled, as that deviatoric strain can result in spherical stress, while volumetric strain can lead to deviatoric stress.
The first term, on the right side of Eq. ( 2), is a liner EOS for isotropic Hookean materials.However, for orthotropic materials, it represents the volumetric (thermodynamic) response.To involve the nonlinear shock effects, the first term of Eq. ( 2), i.e. the contribution to pressure from the volumetric strain, is modified by Mie-Grüneisen EOS.Meanwhile, the contribution to pressure from the deviatoric strain is remained as a correction.Thus the expression of pressure can be rewritten as: In this study, the polynomial EOS was selected and developed using Fortran language.
    where, A T  , and the other parameters will be determined by experiments or simula- tions (Anderson et al., 1994).

Failure Model
The failure model developed in this paper is an orthotropic failure model, which includes two stages, i.e. initial failure and post-failure response.The failure initiation can be based on any combination of material stress and/or strain, and the material strain criterion is adapted in this paper.Subsequent to failure initiation, stiffness and strength properties of the failed material will be updated based on the direction or modes of the failure (Silva et al., 2005).
(1) Delamination failure (11-direction) In this study, the through-thickness direction of composite laminate is defined as the 11-direction, while the in-plane principal directions are 22-and 33-directions.Delamination will be caused by the excessive stresses (or strains) in the 11-or 12-plane direction.If failure is initiated in either of these two modes, the stress in the 11-direction and the corresponding orthotropic stiffness coefficients ij C is instantaneously set to zero (Tham et al., 2008;Bandaru et al., 2015Bandaru et al., , 2016)): Eq. ( 7) shows that this failure model is independent of the directions, i.e. failure in 11-direction does not affect 22-and 33-directions.However, the delamination will in practice be associated with a reduction in shear stiffness, and the fractional residual shear stiffness is maintained through the parameter,  (0.0 to 1.0).Usually, the value of the parameter,  is obtained by experimental tests but in the absence of appropriate material data, a nominal value of 20% is typically used (Silva et al., 2005;Tham et al., 2008).
(2) In-plane failure The 22-and 33-directions are assumed to be in the plane of the composite.If failure is initiated in 22-and 33-directions, the post-failure response is similar to 11-direction, (3) Combined failure If all the three material directions fail simultaneously, the material stiffness and strength will become isotropic with no stress deviators and tensile material stresses, indicating that the material can withstand only hydrostatic pressure.

Erosion Model
Due to the orthogonality of composite laminates, some elements in the impact region will be more easily distorted during the simulations, which may lead to numerical instabilities or even calculation failure (Deka et al., 2008).In this study, element effective strain is used as the erosion criterion, similar approach has been adopted in the references (Grujicic et al., 2009;Bandaru et al., 2015Bandaru et al., , 2016)).
During the impact process, the composite material will undergo non-failure and failure states and the corresponding cell strain will also undergo this two stages.For instance, 1 that composite laminates cannot withstand any loading in the thickness direction.Therefore, the 1  will increase dramatically until it arrives to 1 erosion  (Figure 2).At this point, the meshes will be extremely distorted and the deletion mechanism will be triggered.
Figure 2: Growth of the strain in the 11-direction of composite laminate.

Inverse Flyer Plate Simulations
An inverse flyer plate test (IFPT, Figure 3), creating uniaxial compression state under high-velocity impact condition, was conducted by EMI (Ernst Mach Institute) (Hayhurst et al., 1999).The projectile plate consists of a 6 mm thick Kevlar/Epoxy plate and a 4 mm thick C45-steel backing plate whilst a 2 mm thick C45-steel plate is stationary as the witness plate (target plate).The layered projectile, carried by a polymer sabot, is accelerated by means of a single stage gas gun.The projectile's velocity is measured via trigger pins mounted at the muzzle of the gun barrel.Once the impact event of projectile/target occurs, shock waves will be generated and propagate in the witness plate, and this signal, i.e. free surface velocity history, will be recorded by a laser velocity interferometer VISAR (velocity interferometer system for any reflector) (Lässig T et al., 2015).The corresponding simulation has been implemented in this study, which can be used to verify the shock response of the developed material model (Figure 4).For comparison, the standard material model, MAT_2 (MAT_ORTHOTROPIC_ELASTIC) in LS-DYAN has been introduce.The numerical and experimental results are compared in Figure 5. Figure 5 shows that the two simulation results present a dramatic difference.Good match between the developed model case and the experiment is achieved, while the curve obtained from MAT_2 case shows a great deviation.It is encouraging and means that the developed material model is able to simulate the nonlinear shock response, i.e. phenomena (2).Structures 14 (2017) 1912-1931 3 VALIDATION OF THE DEVELOPED MATERIAL MODEL

Numerical Simulation
The theoretical framework of the developed model was discussed and its nonlinear shock response has been also demonstrated in Section 2. In this section, the performance of this model will be further verified in high-velocity impact condition.The high-velocity impact/penetration experiment is usually conducted by one-stage powder/gas gun system, the main components of which are breech, barrel, safety vessel and target chamber (Figure 6).The projectile and its sabot are launched together by the gun and the impact velocity is controlled by the amount of gun powder.In the safety vessel, the projectile and sabot will be separated before arriving the target; the projectile's velocity is measured by photoelectric diodes and the impact process is recorded by high speed camera.The classic ballistic impact case is selected form the literature (Hoof et al., 1999), in which the composite target consists of 19-plies Kevlar29/Phenolic fabric and its total thickness is 9.5 mm.The projectile is fragment simulating projectile (FSP) based on STANAG 2920 and US MIL-P-46593 (Figure 7).This is a 1.1g FSP that is fabricated using 4340 steel and heated treated to 29 HRC (Hoof et al., 1999;Silva et al., 2005).Johnson-Cook material model, including strain/strain rate hardening and thermal softening effects, is selected to model the FSP; while the target plate is modeled by the developed material model.The parameters of this two material models are summarized in Table 1 (Ramadhan et al., 2013;Ansari and Chakrabarti, 2016) and Table 2 (Tham et al., 2008;Bandaru et al., 2015) Brick elements with single integration point, which are more robust than the fully integrated elements, are used to discretize the target and projectile.Each ply of composite laminate is modeled using a single solid element, and thus 19 elements are assigned along the thickness direction.Considering computational accuracy or robust, the cell size (about 0.5 mm) of the projectile remains consistent with the target plate.An eroding surface-to-surface contact algorithm (CON-TACT_ERODING_SURFACE_TO_SURFACE) is used to simulate the contact between the projectile and composite laminate.Moreover, stiffness-based hourglass control algorithm has been also introduced to hinder the hourglass mode.The finite element impact model for projectile/target is shown in Figure 8.

Discussion of the Simulation Results
(1) Failure modes A series of simulations with different impact velocities were conducted.The case with the initial impact velocity of 590 m/s is selected as an example to discuss the simulation results.Figure 9 shows the plot of strain-11 and the red color denotes the failure strain-11 (seen Table 1).The failure mode, in the form of a crucifix, means that cracks occur and evolve along the warp and weft direction on the top surface of the target plate.This cross-shaped damage pattern was also reported in many references (Silva et al., 2005;Tham et al., 2008;Bandaru et al. 2016), which presents the characteristic of anisotropic strength degradation, i.e. phenomena (3) mentioned in Section 2. The damage pattern on the back surface is shown in Figure 10.Due to the effects of shock-wave reflection/composite-laminate bending or both, the stress state on the back surface is tensile.The phenomena of back-face bulging and fiber fracture were observed in the simulation.During the penetration process, three distinct armor-penetration stages are simulated which correspond to three major failure modes (Cheng, 2003).As shown in Figure 11, in the initial penetration stage, the projectile highly compress the composite plate which lead to shear/cut loading conditions formed at the edge of FSP.Punching shear failure is the dominate failure mode, and the damage of the fibers and matrix is mostly around the rim and underneath the projectile in this stage.With the decrease of the projectile velocity, fibers absorb the projectile's kinetic energy by stretching, and meanwhile the damage is fully extended to the outer zones.Fibers breakage is the major failure mode in the second stage.Due to the reflection of the compression wave and extensive bulging of the backface, the final stage is dominated by delamination and extensive stretching of the filaments.These three penetration stages are simulated and shown in Figure 11.The tensile and crimped meshes in the impact zone present the fiber stretching, shearing and crushing.The red color demonstrates the distribution of the failure strains ( 1 failure  ) indicating the delamination damage of the target plate.In the damage region near the back surface (about 1/3), meshes have been extremely stretched and some even have reached the erosion criteria ( 1 erosion  ).It just goes to show that more serious delamination failure (or a completely damaged state) has occurred in the final stage.However, the developed model is a macro-mechanical model, and thus the plot of failure strain-11, representing the distribution of the damage, only reflects the macro-mechanical properties or has the average mechanical meaning.
(2) Ballistic limit velocity Above-mentioned, the failure modes (qualitative analysis) of composite laminates were discussed.To further verify the developed model, the ballistic limit velocity (BLV, quantitative analysis) is also calculated.There are two standard methods to measure the BLV.The first one is the velocity history method and the second one is based on US MIL-STD-662E standard (Bandaru et al., 2016).In the first method, the maximum impact velocity where the projectile can be stopped, is utilized to be defined as the BLV ( bl V ).In the second method, BLV ( 50 V ) is the average of impact velocities with equal number of partial and full penetrations within a small velocity range of 38 m/s.To ensure the numerical reliability and accuracy, the two methods are available and compared in the present study.
The velocity time histories of projectiles with different initial velocities are shown in Figure 12a.The positive residual velocity means that the projectile has fully penetrated the target plate, while the negative represents partial penetration and rebounding.From Figure 12a, it can be observed that for the case of 590 m/s, its residual velocity is -8.03 m/s.It means that the most kinetic energy of the projectile has been absorbed by the target; negative residual velocity means that rebounding is emerged but not very large.Therefore, the impact velocity of 590 m/s can be considered as the BLV ( bl V ).For the case of 610 m/s, its residual velocity is +50.5 m/s which means the target plate has been fully penetrated.Comparing the above two cases, a phenomenon can be observed that slight increase of the initial velocity (that beyond the BLV) leads to a rapid increase of the residual velocity, similar to the report in the reference (Zhu et al., 1992).Figure 12b shows the computational result base on US standard.Six impact simulations have been conducted and among which three cases are full perforation and the other three are partial perforation.The average impact velocity, i.e.50 V is 586.5 m/s.The BVLs measured by two methods are very close indicating that the numerical result is valid or reliable in turn.
Considerable efforts of building analytic model to predict the BLV for composite laminates have been performed.Based on the assumption that the mean pressure provided by the target consists of two parts, namely quasi-static stress and dynamic stress; Wen (2000) proposed a theoretical model to calculate the BLV for composite laminates: where, M and D are the projectile's mass and diameter; T , t  and Y  are the thickness, density and elastic limit of the composite target, respectively, and  is the geometry parameter of the pro- jectile.In this case, the parameters used to calculate the BLV are listed in Table 3 (Hoof et al, 1999;Wen, 2000).The BLVs obtained from simulation (using two methods), analytical equation (Wen, 2000) and experiment (Hoof et al., 1999) are summarized and compared in Table 4.It can be observed that the numerical results either method 1 or method 2 are greatly close to the experiment, but the analytical result is a little conservative.In generally, the errors of the BLVs gotten from the three ways are acceptable, which shows that the developed model is correct or reasonable in terms of predicting the BLV.

LAMINATES
In this section, the computational capability of the developed model, for composite laminates with different architectures, has been evaluated.The residual velocity curves of two different Kevlar/Epoxy laminates, i.e. plain-woven and cross-plied laminate, have been predicted and compared with the experimental results (Wang et al., 2005).

Numerical Simulation
In the impact tests, composite laminates consist of 75% Kevlar and 25% epoxy (vol.%), which are bonded by chemical method.One kind of laminates is woven from yarns while the other one is stacked in a sequence of 0 90   (seen Figure 13a and Figure 13b).The thickness of those two kinds of composite laminates is 5 mm with 10 plies, each 0.5-mm thick.The diameter and length of the cylindrical shape projectile are 7.62 mm and 15.24 mm, respectively; the projectile is prepared from 4340 steel and heat-treated to a hardness of HRC 32  1 (Wang et al., 2005).Not considering the specific structures, both the plain-woven and cross-plied laminates are equated into an orthotropic homogeneous plate, which is modeled by the developed model (seen Figure 13c).Same to the last section, uniform 3D solid elements have been adapted and 20 elements are assigned along the thickness direction.Because a negligible deformation of the projectiles is observed in the test (Wang et al., 2005); the projectile in this simulation is modeled using rigid material model (MAT_20, ).This is a cost efficient model because that the rigid elements avoid allocating storage to record the history variables during the simulation (Zeng et al., 2005).Table 5 (Xin and Wen, 2012) gives the parameters of composite laminate, and the finite element model is shown in Figure 13.

Comparison of the Numerical and Experimental Results
To predict the residual velocities, 12 simulation cases with different initial velocities are conducted; the experimental and numerical results are summarized in Table 6. Figure 14 shows a damage cone in the impact zone at 10 s  for the 553-m/s case that is consistent with Gellert's model (Gellert et   al., 2000).Similar to Section 3, the red color represents the distribution of the failure strains ( 1 failure  ), which demonstrates the average damage evolution in the transverse direction.this two different models and the results are compared with the experiment.Because of the introduction of a polynomial EOS, the developed model can better simulate the shock response of composite laminates than MAT_2 at high-velocity impact condition.
The ballistic performance of the developed model is assessed from qualitative and quantitative analysis.The classic ballistic impact problem, FPS impacting on Kevlar29/Phenolic laminate, is investigated.Using this model, the top-face cross-shaped failure, back-face bulging failure and three typical penetration processes have been simulated.It should be notice that the failure/damage predicted by this model only possesses the macro-mechanical meaning.BVL is also calculated by two numerical methods which agrees well with results obtained from experiment and theoretical analysis.
The computational capacity of this model, for Kevlar/Epoxy laminates with different architectures, i.e. plain-woven and cross-plied laminates, is further evaluated.The residual velocity curves and the damage cone of this two different laminates have been predicted and the results are compared with the experiment and theoretical model.The general trend of the numerical and experimental curves shows a good correlation and the damage cone in the impact zone is consistent with the theoretical model.The study shows that the developed model is a macro-mechanical model, which is suitable for different laminate structures, and the major ballistic impact behaviors of composite laminate with different structures can be captured by this model.
In general, the macroscopic parameters, such as ballistic limit, residual velocity, energy absorption and penetration/perforation behavior, can be accurately calculated by this model.This indicates that equating the whole composite laminates into an orthotropic homogeneous plate, to capture their macroscopic mechanical properties, is reasonable and feasible.In the future work, further research will be conducted to improve this model and node separation technique will be also introduced.

Figure 1 :
Figure 1: Schematic representation of the equivalent conversion for composite laminate.


is the strain in the thickness direction, while 1 failure  and 1 erosion  are the corresponding failure and erosion strain components, respectively.In the non-failure stage, the 1  will slowly grow with other strain com- ponents before it arrives to 1 failure  .However, after arriving to 1 failure  , failure occurs, which causes Latin American Journal of Solids and Structures 14 (2017) 1912-1931

Figure 5 :
Figure 5: Comparison of numerical and experimental flyer velocity histories for two material models.

Figure 6 :
Figure 6: Schematic illustration of one stage powder gun.

Figure 7 :
Figure 7: Geometrical properties and mesh model of FSP.

Figure 9 :
Figure 9: Cross-shaped damage on the top surface of target plate.

Figure 11 :
Figure 11: Failure modes in three distinct armor-penetration stages.

Figure 14 :
Figure 14: Damage cone in the impact zone.
, a computational macro-mechanical model coupling nonlinear EOS is developed and implemented in the nonlinear dynamic explicit finite element code LS-DYNA.The model's theoretical framework including constitutive model, EOS, failure and erosion model has been introduced.The performance and computational capability of this developed model under high-velocity impact condition are verified and evaluated, and the following conclusions are drawn: Inverse flyer plate simulations are conducted using the developed model and standard material model (MAT_2).The back-surface velocity history curves of composite laminate is calculated by Latin American Journal of Solids and Structures 14 (2017) 1912-1931

Table 1 :
, respectively.Johnson-Cook material model parameters for the FSP.

Table 3 :
Material and physical properties of the FSP and composite laminate.

Table 4 :
Comparison of ballistic limit velocity for different methods.

Table 5 :
Composite material model for Kevlar/Epoxy panels.
Latin American Journal of Solids and Structures 14 (2017) 1912-1931 Figure 13: Composite laminates and finite element impact model.