Velocity and Shape of Explosive Computation using Multi-Material and ALE Formulations

In this paper, a mathematical and numerical description of the bulk viscosity for an equation of state that is linear in density is presented. The bulk viscosity is used in many academic and industrial dynamic codes, and there is no description concerning the smearing of the shock for engineers and researchers in the manuals or in published papers. To clearly show the usefulness of the bulk viscosity, a simple one-dimensional problem is used, where a shock is developed through a pressure wave travelling inside a compressible fluid. By adding a viscous pressure to equilibrium equations, high oscillations in the front shock have been considerably attenuated, by thickening the shock over few element mesh sizes. The method is developed mathematically for one dimensional hydrodynamic problem but has been used successfully for more complex applications including high-impact problems, explosive detonation in air and underwater explosions. Application of the method to a complex problem is illustrated in calculation of the peak velocity and shape of an explosivelyformed projectile (EFP). The symmetry common to most EFPs permits their characterization using 2D axisymmetric analysis. Formation of an EFP entails volumetric expansion of the explosive and extensive plastic flow of the metal plate, both of which can be calculated using an Arbitrary Lagrangian Eulerian (ALE) method. Accordingly, a 2D axisymmetric ALE was used to calculate the velocity and shape of an EFP. The methodology was validated against EFP velocity and shape measurements published in SAND-92-1879 [Hertel 1992]. The Jones-Wilkins-Lee (JWL) equation of state (EOS) were used for the LX14 high explosive backing the copper plate. The explosive burn was initiated using a high explosive material which converts the explosive charge into a gas at high pressure and temperature. The copper plate and steel casing were included using the constitutive model developed by Johnson and Cook. An equation of state developed by Grüneisen for high-pressure simulation was used for the metals. The calculated peak velocity of the EFP was in excellent agreement with the peak velocity published by Hertel. However, the calculated shape did not agree with the experimental shadowgraph of the plate. Specifically, the calculated shape was elongated compared to the measurement and continued to elongate as long as the calculation was continued. In other words, the shape of the copper plate did not reach a dynamic equilibrium. ___________________________________ *Corresponding Author: mhamed.souli@univ-lille1.fr 32 Velocity and Shape of Explosive Computation using Multi-Material and ALE Formulations The methodology for calculating the EFP peak velocity and shape is described. The calculated results are compared to measurements from Hertel. Finally, possible sources for the inaccuracy of the calculated shape are investigated. These include the element size and formulation, initial geometry of EFP, explosive equation of state and the constitutive model for


INTRODUCTION
Simulation of supersonic flows and associated Fluid Structure Interaction (FSI) increasingly has become a focus of computational engineering for industrial and academic applications.Attempts to solve these problems encounter two major challenges: high mesh distortion and spurious oscillations in the shock front.Shocks have a narrow thickness, on the order of a few mean free path collisions in the ambient gas.Thus, to maintain accuracy, the mesh size should be scaled until the shock is resolvable by each individual element.In practice, this method is not viable because the algorithm is requested to handle a massive amount of computational time.Furthermore, the equations of conservation of mass, momentum and energy across a shock require that kinetic energy be transformed into internal energy or heat.In the absence of physical viscosity in the immediate vicinity of the shock, an artificial, unphysical one is commonly added to dissipate the excess of energy.This has the effect of thickening the shock and smearing the discontinuity into a smooth transition zone, and thus, the shock is automatically captured on the computational mesh.
Various approaches have been investigated to solve FSI problems.One of the commonly used approaches to solve these problems is the Arbitrary-Lagrangian-Eulerian (ALE) formulation [Aquelet, Souli, and Olovson (2005)], which has been used with success in the simulation of fluid with large motion, such as a sloshing fuel tank in the automotive industry and bird impact in aeronautical industry.Once the simulations are validated by experimental test results, it can be used as a design tool for the improvement of the structure involved.
Several papers have been published in the literature to discuss methods for solving mesh distortion.Meshfree methods such as SPH (Smooth Particle Hydrodynamic) [Al-Bahkali, Elkanani and Souli (2015); Al-Bahkali E., Souli  In this paper, we examine the ALE method and its description and implementation using explicit time integration method.The Lagrangian formulation, where the mesh moves with material, is mainly used to solve problems in solid mechanics but can be used for FSI problems [Moatamedi, Souli and Al-Bahkali (2014); Elkanani, Al-Bahkali and Souli (2017); Fan, Yao, and Yang et al. (2018)].For small deformation, Lagrangian formulation can solve fluid structure interface and material boundary accurately.The main limitation of this formulation for large fluid deformation or a moving structure is mesh distortion, since the fluid is solved on a moving domain due to the structure's motion.
The paper concentrates on the validation of an ALE methodology that has been implemented in an explicit finite element structural dynamics code to be able to simulate fluid structure interaction problems, where the fluid mesh can be defined by an ALE or Eulerian mesh.From an algorithmic point of view, a fluid element can contain more than one type of fluid material, for structural inflation in ambient air, an element may contain one or two different materials, inside structure pressurized gas and outside ambient air.During the simulation, state variables are computed and stored for each material in each element.An interface tracking algorithm based on Youngs' method is used to capture the interface between the two materials inside the element.This method was used successfully to model many industrial and academic applications as the sloshing tank problem.
In Section 2 of this paper we describe the ALE formulation of the equilibrium equation in a moving mesh and the advection algorithms used to solve mass, momentum and energy conservation.It is well known from previous papers that in the absence of physical viscosity, high nonphysical oscillations are generated in the immediate vicinity of the shock.These oscillations can generate high mesh distortion.To treat shock problems, and avoid high spurious oscillations, a viscous pressure is added to equilibrium equations.In Section 3 we prove in a mathematical form using a one-dimensional case that density, pressure and velocity of the new equilibrium equations vary rapidly but continued across the shock.The method has been applied to one dimensional flow but can be used for the study of more complicated flows and fluid structure interaction problems.
The aim of Section 4 is to illustrate the method for a three-dimensional problem involving an explosively-formed projectile (EFP).The symmetry common to most EFPs permits their characterization using 2D axisymmetric analysis.Formation of an EFP entails volumetric expansion of the explosive and extensive plastic flow of the metal plate, both of which can be calculated using an Arbitrary Lagrangian Eulerian (ALE) finite element method.Accordingly, the 2D axisymmetric ALE capability is applied to calculate the velocity and shape of an EFP.The methodology was validated against EFP velocity and shape measurements published SAND-92-1879 [Hertel 1992].
In the example, the calculated peak velocity of the EFP was in excellent agreement with the peak velocity published by Hertel.However, the calculated EFP shape did not agree with the experimental shadowgraph.Specifically, the calculated shape was elongated compared to the measurement and continued to elongate as long as the calculation was continued.
For the example, the calculated results are compared to measurements from Hertel.Finally, possible sources for the inaccuracy of the calculated shape are discussed.These include the element size and formulation, initial geometry of EFP, explosive equation of state and the constitutive model for the copper plate.

EQUILIBRIUM EQUATIONS IN ALE FORMULATION
Fluid problems, in which interfaces between different materials (gas and ambient air) are present, are more easily modeled by using a Lagrangian mesh.However, if for example an analysis of complex tank geometry is required, the distortion of the Lagrangian mesh makes such a method difficult to use, and many re-meshing steps are necessary for the calculation to continue.Another method to use is the Eulerian formulation.This change from a Lagrangian to an Eulerian formulation, however, introduces two problems.The first problem is the interface tracking [Young (1982)], and the second problem is the advection phase or the advection of fluid material across element boundaries.
To solve these problems, an explicit finite element method for the Lagrangian phase and a finite volume method (flux method) for the advection phase are used.We can refer to several explicit codes; see Benson [Benson (1992)] for a full description of the explicit finite element method.The advection phase has been developed for extending the range of applications that cannot use the Lagrangian formulation only.Current applications include sloshing involving a free surface and high velocity impact problems, where the target is modeled as a fluid material, thus providing a more realistic representation of the impact event by capturing large deformations.
An ALE formulation contains both pure Lagrangian and pure Eulerian formulations.The pure Lagrangian description is the approach that the mesh moves with the material, making it easy to track interfaces and to apply boundary conditions.Using an Eulerian description, the mesh remains fixed while the material passes through it.Interfaces and boundary conditions are difficult to track using this approach; however, mesh distortion is not a problem because the mesh never changes.In solid mechanics, a pure Eulerian formulation is not useful, because it can handle only a single material in an element, while an ALE formulation is capable of handling more than one material in an element.
In the ALE description, an arbitrary reference coordinate is introduced in addition to the Lagrangian and Eulerian coordinates.The material derivative with respect to the reference coordinate can be described in equation (1).Thus, substituting the relationship between the material time derivative and the reference configuration time derivative derives the ALE equations.
Where   is the Lagrangian coordinate, xi the Eulerian coordinate,   is the relative velocity.Let denote by  the velocity of the material and by  the velocity of the mesh.In order to simplify the equations we introduce the relative velocity  =  − .Thus, the governing equations for the ALE formulation are given by the following conservation equations (2) to (4): (i) Mass equation.Note that the Eulerian equations commonly used in fluid mechanics by the computational fluid dynamics (CFD) community, are derived by assuming that the velocity of the reference configuration is zero and that the relative velocity between the material and the reference configuration is therefore the material velocity.The term in the relative velocity in (3) and ( 4) is usually referred to as the advective term, and accounts for the transport of the material past the mesh.It is the additional term in the equations that makes solving the ALE equations much more difficult numerically than the Lagrangian equations, where the relative velocity is zero.
In the second phase, the advection phase, transport of mass, internal energy and momentum across cell boundaries are computed; this may be thought of as remapping the displaced mesh at the Lagrangian phase back to its original or arbitrary position.
From a discretization point of view of ( 2), ( 3) and ( 4), one-point integration is used for efficiency and to eliminate locking.The zero energy modes are controlled with an hourglass viscosity [Benson (1992)].A shock viscosity, with linear and quadratic terms, is used to resolve the shock wave [Von Neumann and Richtmyer (1950)]; a pressure term is added to the pressure in the energy equation ( 4).The resolution is advanced in time with the central difference method, which provides a second order accuracy in time using an explicit method in time.For each node, the velocity and displacement are updated using explicit time integration.
The multi-material formulation is attractive for solving a broad range of nonlinear problems in fluid and solid mechanics, because it allows arbitrary large deformations and enables free surfaces to evolve.The Lagrangian phase of the VOF method is easily implemented in an explicit ALE finite element method.Before advection, special treatment for the partially voided element is needed.For an element that is partially filled, the volume fraction satisfies   ≤ 1, and the total stress by σ is weighted by volume fraction   = .  .
In the second phase, the transport of mass, momentum and internal energy across the element boundaries is computed.This phase may be considered as a re-mapping phase.The displaced mesh from the Lagrangian phase is remapped into the initial mesh for an Eulerian formulation, or an arbitrary undistorted mesh for an ALE formulation.
In this advection phase, we solve a hyperbolic problem, or a transport problem, where the variables are density, momentum and internal energy per unit volume.Details of the numerical method used to solve the equations are described in detail in [Von Neumann and Richtmyer (1950)], where the Donor Cell algorithm, a first order advection method and the Van Leer algorithm, a second order advection method [Young (1982)] are used.As an example, the equation for mass conservation is: It is not the goal of this paper to describe the different algorithms used to solve equation ( 6); these algorithms have been described in detail in [Aquelet, Souli, and Olovson (2005); Ozdemir, Souli and Fahjan (2010)].

ARTIFICIAL SHOCK VISCOSITY
Mathematically, shocks are treated as discontinuities in the history variables, energy, pressure and velocity; but physically, shocks have a narrow thickness on the order of a few collisions mean free path in the ambient gas.For air at Standard Temperature and Pressure (STP), the mean free path is estimated to 70 nanometers, [Souli and Gabrys (2012)].Thus, to keep the accuracy of the results, the mesh size should be scaled until the shock is resolvable by each individual element.In practice, this method is not viable because the algorithm is requested to handle a massive amount of processor time, and thus numerically the shock front cannot be represented by any mesh.
To clearly describe shock phenomena from a mathematical and numerical point of view, we consider equilibrium equations described previously for a non-viscous fluid, hydrodynamic equations, called inviscid Navier-Stokes equations or Euler equations: The system of equations (6-8) is closed by using an equation of state for pressure.An important phenomenon that arises with these equations is the formation of shock; mathematically equations (6-8) develop a shock, which lead to a discontinuity, and the problem is well-posed only if the shock conditions are satisfied.These conditions, called Rankine-Hugoniot conditions, describe the relationship between the states on both sides of the shock for conservation of mass, momentum and energy across the shock and are derived by enforcing the conservation laws in integral form over a control volume that includes the shock (see Fig. 1).
In the absence of physical viscosity, high non-physical oscillations are generated in the immediate vicinity of the shock.To illustrate this phenomenon, we look at a problem of pressure wave propagating in a 10m long tube, loaded by a pressure pulse of 1 bar, applied at the left end of the tube.An initial pressure of 0.994 bar is applied to the air inside the tube as described in Fig. 2. To illustrate a one-dimensional problem, both velocities on vertical and normal directions are constrained.
For an ideal gas equation of state, to treat the shock and avoid numerical oscillations at the front shock, Von Neumann et al. [Von Neumann and Richtmyer (1950)] proposed a viscous pressure term, in the equilibrium equations (6)(7)(8), called bulk viscosity.By adding an appropriate viscous pressure  to the conservation equations, the full continuous physics is described and equations (9)(10)(11) possesses a continuous solution.
The bulk viscosity  is added to each element in compression, and no viscosity is added for element that is in tension.By adding the bulk viscosity , equation (11), we show that the equations (9-11) possess a continuous solution that satisfy conservation laws as well as Rankine-Hugoniot conditions.(11) where the viscous pressure  is given by:  = (.) 2 ..() 2  () < 0 is the fluid velocity,  a characteristic length of the shock smearing, and  a scale factor for the thickness of the smearing.
The viscous pressure irreversibly converts mechanical energy into heat; and allows numerical simulation of shock fronts by meeting the following conditions: 1. Thickness of the shock layer is of the same order as the computational grid, independently of the material and the strength of the shock.2. Effect of Q is negligible outside the shock layer and Hugoniot conditions hold outside the shock layer.3. Mass, energy and momentum are exactly conserved across the shock front.Rerunning the same problem described in Fig. 1, by using equations (9-11), with bulk viscosity, equation ( 12), Fig. 2 shows both pressures with and without the viscous pressure Q.
It is clear from these two examples that the use of the bulk viscosity in the simulation of non-viscous fluid for impact problem leads to a smooth solution that solve spurious high oscillations in front shock.Most industrial and academic software for dynamic explicit analysis refer to Von Neumann et al. [Von Neumann and Richtmyer (1950)] without giving any guidelines or details to users.For an equation of state that is linear in density, no paper has been published in the literature, showing that equations (9)(10)(11), with Q defined by equation ( 12), possess a continuous solution that satisfies conservation laws.It is the aim of this paper to show that for an equation of state that is linear in density, equations (9-11) possess a continuous solution and for the shock layer on the same order as the computational grid.
To illustrate high oscillations generated near a shock, it is useful to plot pressure along the computational domain.This can be clearly described in [Souli et al. (2018)], where oscillations can be drastically reduced in the vicinity of the shock when bulk viscosity is added to the hydrodynamic equations.Using viscous pressure term or bulk viscosity to hydrodynamic equations of equilibrium, spurious oscillations of high frequencies can be damped out without loss of accuracy.This approach has been used in many nonlinear software simulations such as LSDYNA Hallquist [Hallquist (2013)], and known as Richtmyer and Von Newman viscosity or bulk viscosity.

. Mathematical description of Artificial Shock viscosity
By adding a viscous pressure to equilibrium equations, we can mathematically prove that the solution is continuous, and the shock is smeared over a distance that is of the order of q.l, in equation (12) For clarity, we consider a one-dimensional equilibrium equation with an equation of state that is linear in density for the pressure.Let  = 1/, the specific volume, from mass equation ( 9  Using the change of variable ( 16) and ( 17), mass and momentum conservation equations ( 9) and (10)  We first combine equations ( 18) and (19) to get conservation of mass and momentum with the new variable w: Equation (18) shows that: where  is a constant with respect to the variable w that need to be defined.Using  = 1/, we get the new momentum conservation equation, that leads to an ordinary equation for  as a function of w.The equation for  will have a continuous solution that can be solved analytically.
Now, let us express the bulk viscosity Q in (15) as a function of the new variable :  = (.  . 0 ) 2 ..
To illustrate the well-functioning of the quadratic bulk viscosity, we use for the pressure an equation of state that is linear in density: using specific volume  = 1/, and  0 = 1  0 , equation of state (24) gives: Now, the pressure term P in equation ( 25) can be expressed in the momentum equation (22): using equation ( 26) away from the shock, where Q = 0: For  =  0 , we have For  =  1 , we have by replacing  defined in ( 27) into (28) we have: Now, we replace Q from (23) and C from (27) in equation ( 28): .  2 . 0 using equation ( 29) for  0 2 . 2 , we get an ordinary differential equation (30) for  as a function of : (30) which can be written as: Equation ( 31) is an ordinary differential equation that has a continuous solution for  as a function of , and can be integrated analytically to give the shape of the front shock: Figure 3: Specific volume as a continuous function of the variable  Numerically the characteristic length  represents the mesh size, using a dimensionless constant  in Fig. 2 of the order of unity, the thickness of the shock layer can be represented by 3 to 5 mesh elements.From equation (31), it is obvious that the shock layer is independent of the strength of the shock and of the material bulk modulus .
By deriving twice the specific volume  with respect to , it is straightforward to show that ( 32) is a solution of the ordinary differential equation (31).V0 V1 w V -0.5π.q.l +0.5π.q.l 4. APPLICATION AND NUMERICAL SIMULATION.
The ALE implementation detailed above was used to calculate the shape and peak velocity of an explosively formed projectile (EFP) as described in in SAND-92-1879 [Hertel 1992].
The results of the calculation were compared against experimental results detailed in the report.Measurements include the peak velocity and shape for an explosively-backed copper flyer plate.The plate is an OFHC (oxygen-free high thermal conductivity) copper explosive lens propelled by the high explosive LX-14.The casing is AISI 4340 steel.Additional details and the cross section of the assembly are provided in SAND-92-1879.

ALE Controls and Domain
The Eulerian domain was composed of area-weighted axisymmetric solids.The domain itself is shown in Figure 4.The domain element size was 0.7mm A sequential filling routine was used to define the subdomains of the air, explosive and metals.Accordingly, the domain was first filled with air.Beam elements were used to define the fill volumes for the explosive, casing and flyer plate, each of which was successively filled.

Constitutive Model
The parts of the axisymmetric model, their constitutive models, and sources of material inputs are summarized in Table 1.The Johnson-Cook damage and fracture model that is included in the constitutive model was used for characterizing damage to the flyer plate and casing material.Use of the Jones-Wilkins-Lee (JWL) EOS was the starting point for the explosive.Hertel reports JWL constants for the LX-14, and JWL is established for calculating explosive expansion.The caution in using JWL is the engineering simplifications that it employs.For example, Fuller et al. at the University of Sheffield (2016) report that the JWL EOS may overpredict peak pressures and impulses for close-in detonation.Fuller et al. indicate that the assumption of near instantaneous energy release from detonation may cause this overprediction.
Given these cautions, alternatives to the JWL EOS were considered.These included Jones-Wilkins-Lee-Baker (JWLB) and ignition and growth of reaction in a high explosive (IGRHE).JWLB was particularly promising because it is a form of the JWL developed by Baker et al. specifically for EFPs and shaped charges.The JWLB form is intended to refine the assumption of a uniformly expanding cylinder to calculate the JWL constants.
Assumption of a uniformly expanding cylinder does not account for all the work done on the cylinder by the explosive.Specifically, some of the work done by explosive products (axial flow of cylinder, bending of cylinder, etc.) is not included in the Gurney energy/volume, a parameter critical to EFP calculations.The implication is that the explosive can do more work than the JWL constants predict.Accordingly, explosives tests would be expected to produce a higher velocity than calculated using the JWL EOS.
The input constants of LX-14 are available for JWLB; in fact, these are included Hallquist 2013.LX-14 inputs for were not available at this time of this study.

Calculation Results -Peak Velocity
As noted above, the outputs of interest from the calculation were the peak velocity and shape of the EFP.The calculated peak velocity was in good agreement with the velocity measured in the experiment.Hertel reports a peak velocity of 0.228 cm/usec.The peak velocity in the calculation (at ~100 usec) using JWLB was 0.228 cm/usec, whereas it was 0.225 using JWL.
These results are compared in Table 3 and illustrated in Figure 5.The calculated velocities in Figure 5 are rigid-body velocities of the EFP along the axis of symmetry.The results summarized in Table 3 are consistent with the development of the JWLB EOS.Specifically, JWLB was intended as a refinement of the JWL EOS that applies all the explosive energy to the plate in contact.The JWLB EOS results in a slightly higher peak velocity compared to JWL.
Less clear is how these results relate to work by Fuller et al.Their finding that the JWL EOS may over-predict peak pressures and impulses for close-in detonation would suggest over-prediction of peak velocity in a case like the Hertel experiment.However, the JWL peak velocity was below the measured velocity.Of course, a close-in detonation does not act on a target identically to how explosive in contact with a plate acts on it, but it is unclear how overprediction in the close-in detonation case can coexist with a slight under-prediction in the EFP velocity.

Calculation Results -EFP Shape
Hertel's reported EFP shape is shown next to the shapes calculated with JWLB and JWL, respectively, in Figure 6, and the dimensions of the calculated and measured EFP shapes are listed in Table 4.These are the EFP shapes at the peak velocity.The basic observation is that the EFP shapes calculated are less compact than the measured shape and the shape reported by Hertel.The calculated shapes for this paper are both longer and have larger diameters.The calculated diameters are particularly exaggerated such that the calculated aspect ratio of length/diameter is 0.8 versus 1.1 in the experiment.Another difficulty with the present calculation is that the length of the EFP shape increases as long as the calculation is continued.The length, and to lesser extent the diameter, never stabilize to a constant value after reaching peak velocity.Evolution of the EFP shape is shown in Figure 7.The results for peak velocity and shape discussed above for the present calculations were performed with 0.7mm square 2D shells for the Eulerian domain.Refining to a 0.35mm changed the peak velocity 5%, and, the EFP shape was not significantly different from the 0.7mm case.It is noted that the calculation with the 0.35mm mesh was less stable than with the 0.7mm mesh and crashed before reaching the termination time of the 0.7mm calculation.

Initial Geometry
Differences in the early-time shape of the EFP suggested that possible differences in initial geometries caused the inconsistent EFP shapes.It was postulated that a difference in overpressure distribution at the explosive-casing junction caused the shape differences.The difference in overpressure distribution, in turn, was caused by differences in the geometry of the explosive-casing junction.Only a screen capture from the Hertel paper was available for developing the fill geometry for the present calculations, and imprecision in the fill geometry may have been a factor in the disagreement.
Accordingly, the sensitivity of the EFP shape to fill geometries was assessed.Calculations were performed where explosive was in the gap between the end of the copper plate and the casing.Other calculations were done where the end of the copper plate was placed in contact with the casing.Varying the fill geometry within these bounds did not significantly affect the EFP shape indicating the shape appeared insensitive to fill geometry.

Explosive Equation of State
As noted above, options for modeling the LX-14 with an EOS other than JWL were limited to JWLB.Input constants could not be obtained for IGRHE.Little difference was observed between the peak velocity and EFP shape, regardless of whether JWL or JWLB was used.

Copper Constitutive Model
The Johnson-Cook constitutive model has been observed to exhibit excessive softness.To check this possibility at a static loading rate, a coupon test was performed using the Johnson-Cook model with inputs for copper [Johnson & Cook 1985].The stress and strain values output from the coupon test were consistent with values published in the Atlas of Stress-Strain Curves.
The calculation was run with the default log-linear Johnson-Cook rate effects for copper.Input constants for the log-linear form were from Johnson and Cook [Johnson & Cook 1985].To assess sensitivity to rate effects form, Cowper-Symonds constants calculated from data published by Lindholm and Bessey (1969) were used.This change in the rate effect formulation for the copper had negligible effect on the EFP shape.
Also, whether the Johnson-Cook failure parameters (D1 through D5) were included in the calculation did not significantly affect the EFP shape.More specifically, the calculation was run two ways: (1) no entries for D1 through D5 for the copper and (2) D1 through D5 for copper as reported in [Johnson & Cook 1985], with erode=1 such that there was no erosion.Rather, for erode ≠ 0, deviatoric stresses are set to 0 upon failure.Again, difference in shape of the EFP was negligible.4.6.5.Axisymmetric Element Formulation Lastly, Benson (1990) has shown that volume-weighted axisymmetric element types underweight the nodes along the symmetry axis.The formulation is therefore not spherically symmetric.This condition causes non-physical jetting along the axis of symmetry; the jetting is a byproduct of the element formulation, not caused by the physics of the calculation.Such jetting may explain the elongation of the EFP observed in the present calculations.However, an area-weighted formulation was used for this calculation, and area-weighted formulations are known to preserve spherical symmetry [Benson 1990].As a result, it is uncertain whether the element formulation could be a contributing factor to the error in the EFP shape.

CONCLUSION
In this paper, a mathematical and numerical description of the bulk viscosity for an equation of state that is linear in density, is presented.The bulk viscosity is used in many academic and industrial dynamic codes, and there is no description concerning the smearing of the shock for engineers and researchers in the manuals or in published papers.To clearly show the usefulness of the bulk viscosity, a simple one dimensional problem is used, where a shock is developed through a pressure wave travelling inside a compressible fluid.By adding a viscous pressure to equilibrium equations, high oscillations in the front shock have been considerably attenuated, by thickening the shock over few element mesh sizes.The method is developed mathematically for one dimensional hydrodynamic problem, but has been used successfully for more complex applications including high impact problems, explosive detonation in air and underwater explosions.To solve complex three dimensional problems, a remeshing technique method is developed in the paper.The principle of an ALE method is based on the independence of the Finite Element mesh movement with respect to the material motion.In fact, the freedom of moving the mesh offered by the ALE formulation enables a combination of advantages of Lagrangian and Eulerian methods.
The 2D axisymmetric ALE finite element implementation accurately calculates the peak velocity of the copper EFP.The calculated peak velocity was within 2% of the peak velocity measured by Hertel if JWL EOS is used.Use of the JWLB EOS gave a peak velocity identical to Hertel's, within the precision of the calculation.
The calculation of the EFP shape was less accurate.The shape of the EFP was elongated compared with the measured shape.In addition, plastic flow continued as long as the calculation was continued; i.e., the EFP shape of the plate did not reach a dynamic equilibrium.
The sensitivity of the EFP shape was checked for these conditions: Eulerian domain element size and formulation; initial geometry of casing, explosive and copper plate; explosive equation of state and constitutive model for the copper plate.The shape was found to be insensitive to these conditions, when they were varied within reasonable bounds.
The excessive deformation of the copper plate was unexpected because the Johnson-Cook constitutive model is validated for multi-material ALE, though not necessarily for 2D axisymmetric ALE.It is possible that Johnson-Cook exhibits such excessive softness specifically when used with 2D axisymmetric ALE.
Another possible source of the excessive softness is the Johnson-Cook constitutive model itself.A candidate explanation is that the plasticity algorithm in Johnson-Cook fails to converge for the tolerance implemented in the LS-DYNA® solver.Accordingly, plastic flow simply continues, independently of convergence.
A limited treatment of rate effects in the Johnson-Cook constitutive model may also contribute to the error.It is known that strain rate effects in steels are function of deformation.UFC 3-340-01 (2002) Figure 4-49 has separate strength rate curves for yield and ultimate.
Copper exhibits similar behavior, as discussed in Lindholm and Bessey (1969).The copper plate undergoes large deformations, and it is likely that rate effects vary over those deformations, a property not included in the Johnson-Cook constitutive model.
It is observed in Lindholm and Bessey (1969) that strength rate effects depend on temperature, but this effect is not explicitly included in the Johnson-Cook constitutive model.
Momentum equation.The strong form of the problem governing Newtonian fluid flow in a fixed domain consists of the governing equations and suitable initial and boundary conditions.The equations governing the fluid problem are the ALE description of the Navier-Stokes equations: conditions need to be imposed for the problem to be well-posed.(iii) Energy equation.   =    , +     −      (4)

Figure 2 :
Figure 2: Shock loading and unloading in an element
. =  0 . 0 , we have:  = (. 0  0 ) are looking for solutions where the shock front propagates with a constant velocity .From a mathematical point of view, it is more convenient to use a change of variable with a constant shock speed  with respect to the undisturbed medium. =  0 − .

Table 1 .
Constitutive Model by Part Equation of StateThe parts of the axisymmetric model, their equations of state, and sources of equation of state (EOS) inputs are summarized in Table2.

Table 2 .
Equations of State by Part

Table 4 .
EFP Shape -Calculation vs. Experiment Figure 7. Shape -Fluid Density [g/cm3] vs.Time [microsec]4.6InconsistentEFP Shape -Possible CausesPossible causes of the inconsistencies between the calculated EFP shape and the experiment were investigated.Accordingly, the sensitivities of the calculated EFP shape to the following conditions were examined: