Modeling fragmentation with new high order finite element technology and node splitting

The modeling of fragmentation has historically been linked to the weapons industry where the main goal is to optimize a bomb or to design effective blast shields. Numerical modeling of fragmentation from dynamic loading has traditionally been modeled by legacy finite element solvers that rely on element erosion to model material failure. However this method results in the removal of too much material. This is not realistic as retaining the mass of the structure is critical to modeling the event correctly. We propose a new approach implemented in the IMPETUS AFEA SOLVER R © based on the following: New High Order Finite Elements that can easily deal with very large deformations; Stochastic distribution of initial damage that allows for a non homogeneous distribution of fragments; and a Node Splitting Algorithm that allows for material fracture without element erosion that is mesh independent. The approach is evaluated for various materials and scenarios: -Titanium ring electromagnetic compression; Hard steel Taylor bar impact, Fused silica Taylor bar impact, Steel cylinder explosion, The results obtained from the simulations are representative of the failure mechanisms observed experimentally. The main benefit of this approach is good energy conservation (no loss of mass) and numerical robustness even in complex situations.


Introduction
Random ductile and brittle fragmentation modelling is still a challenging task.Indeed, such fragmentation generally occurs in a complex multi-physics framework, as in classic tube explosion or implosion [1].
From a mechanical point of view a reliable material model is essential.At first, it is necessary to capture the strain localizations at the right time in order to deal with fragment size.Classic models for ductile materials used for this purpose are: Steinberg-Cochran-Guinan [2], Johnson-Cook [3], Zerilli-Armstrong [4], Mechanical Threshold Stress [5] and their evolutions.Most of these models take into account strain, strain rate, temperature, pressure, saturation stress and very high strain rates viscosities.They seem accurate enough to be used for the numerical simulations of the strain localizations.Moreover, random localizations shouldn't be naturally initiated in accurate numerical simulations.It is necessary to introduce a stochastic aspect in material models like random elastoplastic properties or a random initial damage distribution.
From a numerical point of view it is necessary to have a very accurate formalism in order to avoid activating fragmentation due to numerical errors.For example, Eulerian and ALE approaches must cope with interface reconstruction approximation errors.This is a very difficult task in the case of strong 3D fragmentation [6].In a pure Lagrangian approach strain localization should only occur due to initial singularities (cross section step . . . ) or wave a Corresponding author: jerome@impetus-afea.com crossing (spalling . . .).This is not the case for example with first order tetrahedron finite element approximations that can produce random fragmentation without any initial material heterogeneities [7,8].Accuracy is thus essential for a reliable random fragmentation numerical model.
Of course, random strain localization is only the first step in modelling a tube explosion.Is it necessary to apprehend 3D domain opening and multi-species interaction (high explosive and metals for example).In most cases, Eulerian or ALE approaches are chosen but here again complex 3D interface reconstruction and contacts are very challenging and leads to energy conservation problems.A pure Lagrangian approach is generally not possible.Actually, solutions to treat domain opening are often limited to element erosion (3D X-FEM or meshless approaches for multiple cracks are still not mature for industrial applications).It means that a very fine mesh is needed to limit mass loss during computation.Furthermore, very large deformations are most of the time impossible to treat with classic FE, specially the HE part.
We propose to describe and evaluate a new monolithic purely Lagrangian approach implemented in the IMPE-TUS AFEA SOLVER R .This approach is based on the following: -New High Order Finite Elements (Advanced Solid Element Technology ASET) that offers very high accuracy compared to classic first order elements -Stochastic distribution of initial damage or elastoplatic properties that allows for a realistic strain localisation EPJ Web of Conferences -Node splitting algorithm that allows for material fracture without element erosion that is mesh "independent".-New Lagrangian Particle method for HE modeling.
The approach is evaluated for various materials and scenarios: -Titanium ring electromagnetic compression; Hard steel Taylor bar impact, Fused silica Taylor bar impact, Steel cylinder explosion,

High order finite elements for transient dynamics 2.1. Description
As highlighted in the introduction numerical accuracy is primordial to capture strain localization as an initiation mechanism for fragmentation.Legacy explicit solver use reduced integration first order finite elements (1 Gauss point) + hourglass control.These elements have a very low accuracy but are very popular because of low computation time.In practice, a large amount of elements is needed to compensate for insufficient accuracy.This is consistent with the very fine mesh needed to limit mass loss for an element erosion strategy for crack modeling.Nevertheless, real 3D applications like warhead design often lead to an unrealistic number of elements.
We propose ASET, an Advanced Solid Element Technology based on fully integrated third order approximations (64 Gauss points).These elements bypass the need for hourglass control and give much superior accuracy and deformability.For sure, it is not realistic to use this type of element with erosion.Thus, we developed a node splitting algorithm (see Sect. 4.1) that allows exact mass conservation and realistic computation time.We illustrate ASET accuracy with a classic pinched cylinder test.

Wriggers's pinched cylinder test
We evaluate the ability of 3rd order hexahedron finite elements in the framework of elastoplastic bending of thin structures.This popular benchmark proposed by Wriggers et al. [9] consists of a pinched tube in the middle.This tube has an inner radius dimension r = 300 mm, a thickness t = 3 mm and a length L = 600 mm.The material is treated as an elastoplastic law, with a Young's modulus E = 3 e 3 MPa, Poisson's ratio ν = 0.3 and isotropic hardening that can be expressed as σ = 24.3+ 300ε p MPa.
In order to compare the solution described in Wriggers et al. [9] with IMPETUS we consider the plastic strain for various displacements (one-half model see Fig. 1).The observed results are very close.Unlike Wriggers who used 1600 shell elements, we used 50 3 rd order elements.We highlight the fact we only used 1 element through thickness.This leads to a very high aspect ratio, about 20.That means we ignored all classic meshing rules imposed by legacy solvers (at least 4 elements through thickness and an aspect ratio smaller than 4).Nevertheless, accuracy is preserved as proved in Fig. 2 by the Force/Displacement comparison.
It is apparent that ASET 3 rd order finite elements can deal with a strongly non-linear problem even with a very bad mesh quality.

Constitutive relations and random initialisation
As depicted in Sect. 1 the introduction of a stochastic distribution of initial damage or initial variation of elastoplastic properties is very important.The combination with a reliable constitutive relation gives the foundation for a fragmentation model.

Stochastic "Defects" distribution
A distribution function f(D) describes the number of "defects" per unit volume of matter.
Note that the maximum initial damage cannot be larger than D max .The number of defects N per unit volume of 04050-p.2 Based on the assumed distribution f(D) one can show that the probability p of having at least one initial defect larger than or equal to D 0 in a volume v is: This probability expression can be used to assign an initial defect level to each integration point in the model.The defect level is obtained by solving the expression for D 0 (given a random number p and an integration point volume v).

Ductile material model
We chose to use a slightly modified version of the popular Johnson-Cook model [3].This is a constitutive model for ductile metals with thermal softening and strain rate hardening.The yield stress is defined as: where g(D) is a damage softening and D is the damage level, defined as: That is the way the yield stress can be randomly distributed at the initial state (see Fig. 3).

Brittle material model
We chose to develop a slightly modified version of the popular Johnson-Holmquist material model [10].This is a constitutive relation for ceramic materials with different failure mechanisms in compression and tension.The material is assumed to have a pressure dependent shear resistance.At positive pressures, plastic flow is a combination of shearing and dilatation.Inelastic dilatation is interpreted as crushing that gradually reduces the shear resistance of the material.A brittle fracture criterion is used on the tensile side (see [11] for more details).

Random ductile fragmentation
In this section, we analyse 3 reference tests for random ductile fragmentation.

Petit's electromagnetic compressions test
The experiment of electromagnetic cylindrical compression was developed at the French Atomic Agency (CEA) to study the behaviour of ductile materials at large strain and high strain rate [12].This section will focus on Ti6Al4V.Clearly this particular titanium alloy is very sensitive to adiabatic shearing.One test carried out at the CEA is used as the example.
The tube inner section is observed with a frame camera.Figure 4 (axial view) shows that adiabatic shear bands appear.Experimental observations proved that he adiabatic shearing threshold does not depend only on the strain but on the strain rate history too.
The CEA numerical result is based on the Ouranos 2D Lagrangian solver and taken as reference (details can be found in [13]).We developed an approach based on third order finite elements (see Sect. 2) and a ductile material model (see Sect. 3.2).Both solvers used an equivalent loading pressure to electromagnetic load (neglecting overheating).
We obtained good correlation between the experiment and our numerical approach (see Fig. 4).One can notice that Petit's model is largely more physically advanced than ours (modified Zerilli-Armstrong).Nevertheless, our goal here is to demonstrate we can capture strain localization with relatively few high order elements (Impetus model : 2400 elements, Ouranos 2D model ∼180,000 elements).This allows extending our approach in 3D (see Fig. 5) with largely reduced computation time compared to the reference model.

Borvik's taylor bar impact
SimLab NTNU conducted very interesting Taylor tests on ARNE tool steel, hardened to nominally Rockwell 04050-p.3C values of 52 [14], see Fig. 6.An impact velocity about 250 m/s leads to quasi brittle failure.Proposed numerical model relies on cubic hexahedron elements (see Sect. 2) and a modified Johnson-Cook constitutive model (see Sect. 3.1 and 3.2).An additional numerical aspect is introduced here: Node Splitting.Indeed, the Fig. 6 model uses a very coarse mesh that prohibits the element erosion approach to generate fragments.We thus developed a node splitting technique based on the following 3 steps:

EPJ Web of Conferences
-Loop over all integration points surrounding the node that is to be split.Mark the two integration points with largest damage.-Define a vector v between these two integration points.
-Try to split the node such that the propagating crack plane has its normal in a direction as close to v as possible.
This model exhibits very good energy balance as no element is eroded during computation.Moreover, fragment mass distribution is in good agreement with experimental data.

Experiment description
This case is inspired by the experiments presented in [15] corresponding to the fragmentation of a cylinder due to HE detonation.
An AISI 1018 steel cylinder height H = 203.2mm, outside diameter De = 50.8mm and a thickness e = 3 mm is fragmented by the detonation of an explosive LX-17.The explosive LX-10 allows initiating the detonation in order to obtain a quasi-plane wave in the cylinder.In addition to the distribution of the fragments, their thickness and microstructure are studied in [15] in order to assess the corresponding fracture deformation.

Blast particle method
High explosive modeling is a challenging task especially in the case of a strong interaction with a strongly deformed/fractured domain.Most of the referenced studies used an Euler/Lagrange coupled approach.As noticed in Sect. 1, this classic approach is very complex and leads to energy conservation problems.We developed a purely Lagrangian Meshless Method applied to modeling the complete Blast Event [16].The implementation of the Blast Particle Method (BPM) in the IMPETUS Afea Solver R is described in detail in [17] and compared with experimental results in [17,18].
As a short description, the HE model is based upon the Kinetic Molecular Theory for gas.The basic assumptions are presented below, but the first two are not valid for HE so modifications to the theory are necessary.
The basic assumptions for Molecular theory are: -The average distance between particles is large compared to the particle size.-Equilibrium exists.
-Molecular collision is perfectly elastic.
Calibration for specific types of explosives is accomplished by using a traditional cylinder test, as shown in Fig. 8. Air is also modeled with the same approach as the HE. .The High explosive is modelled using the innovative Discrete Particle Method (DPM) [16].Coupling between DPM and FEM is based on contact algorithm that transfers particles impulses to the structure.DPM is independent of the structures complexity thanks to its meshless nature.This means that treatment is the same for a simple cylinder or for a complex warhead.The main advantage of this monolithic Lagrangian approach is exact mass conservation and very good energy conservation.Moreover, a close range cylinder explosion impacting a structure like a plate is trivial.This means that both fragment creation and their effect on the plate would be taken into account.A visualization of the numerical model is presented in Fig. 9.We can distinguish red dots that are related to DPM and FEM in blue.The dynamics in the simulation is in good agreement with experiments [15].The mesh size used is quite coarse and limits the analysis of fragment morphology.Indeed, experimentally obtained fragments exhibit inclined shear fractures parallel to the fragment length, whereas the fragment length is defined by flat fracture surfaces.These inclined shear fractures are not reproduced in this model due to the coarse mesh and the Node Splitting Algorithm itself.This aspect requires a specific study.Nevertheless, the obtained fragment mass distribution is acceptable as shown in Fig. 10.
Regarding computational considerations, the GPU parallelization of the solver is very efficient.In fact, the model ran in less than 1h on a standard workstation equipped with an Nvidia K20 GPU processor.This aspect is very important as it allows for real industrial applications and future developments.

Brittle random fragmentation
We now focus on evaluating the approach on a brittle material: fused silica.Fused silica is a high purity synthetic amorphous silicon dioxide widely used in the military industry as window material.It may be subjected to highenergy ballistic impacts.Under such dynamic conditions, post-yield response of the ceramic as well as the strain rate related effects become significant and should be accounted for in the constitutive modelling.In this study, we use a modified Johnson-Holmquist (J-H) model (see Sect.   an inverse calibration technique, on selected validation test configurations. Numerical simulations were performed with 3rd order hexahedron elements (see Sect.Impetus Afea Solver) [11].
Taylor impact at 89 m/s is shown in Fig. 11.Spall fracture development in the rear portion of the cylinder as as separation in the front section crushed region is exhibited.Qualitative comparison with a high speed camera selected frame is in a fairly good agreement.

Conclusions
It was demonstrated on several challenging ductile and brittle random fragmentation cases relevance of this new approach.The ASET Element Technology coupled to the Discrete Particle Method (DPM) demonstrated a very good alternative to classic approaches as it proposes a fully Lagrangian framework with exact mass conservation with reduced computation time.Industrial applications targeted are directly linked to warhead design.
Material models used in this study were quite simple as the focus was on evaluating this very innovative approach.This means that in order to be really predictive a strong effort is required to introduce more physics into the material models.
that are used to characterize the HE particle model: ρ Density E0 Internal energy D Detonation velocity ν Co-volume (Optimized in cylinder test simulations).γ Ratio heat capacities(Optimized in cylinder test simulations).

4. 3 . 3 .
Numerical model Here again the proposed numerical model for ductile steel cylinder is based on 3 rd order hexahedron elements (see Sect. 2), modified Johnson-Cook constitutive model (see Sect. 3.1 and 3.2) and Node Splitting Algorithm (see Sect. 4.1) 3.2).Parameters have been identified by Ruggiero et al.[11] by
2), a modified Johnson-Holmquist constitutive model (see Sects.3.1 and 3.3) and a Node Splitting Algorithm (see Sect. 4.1).This approach overcomes classic numerical drawbacks associated with element erosion and numerical inaccuracy of fragmentation activation.04050-p.5 EPJ Web of Conferences