Lagrangian and ALE Formulations for Soil Structure Coupling with Explosive Detonation

Simulation of Soil-Structure Interaction becomes more and more the focus of computational engineering in civil and mechanical engineering, where FEM (Finite element Methods) for soil and structural mechanics and Finite Volume for CFD (Computational Fluid Dynamics) are dominant. New advanced formulations have been developed for FSI (Fluid Structure Interaction) applications using ALE (Arbitrary Lagrangian Eulerian), mesh free and SPH (Smooth Particle Hydrodynamic) methods. In defence industry, engineers have been developing protection systems for many years to reduce the vulnerability of light armoured vehicles (LAV) against mine blast using classical Lagrangian FEM methods. To improve simulations and assist in the development of these protections, experimental tests and new numerical techniques are performed. Initial conditions such as the loading prescribed by a mine on a structure should be simulated adequately in order to conduct these numerical calculations. The effects of blast on structures often depend on how the initial conditions are estimated and applied. This article uses two methods to simulate a mine blast, namely the classical Lagrangian as well as the ALE formulations. The comparison was carried out for a simple and also a more complex target. Particle methods as SPH method can also be used for soil structure interaction.


INTRODUCTION
Theoretical and experimental analysis of explosion mine loading have been considered by several researchers over the past decades, using empirical methods as CONWEAP (Conventional Weapon) code when the explosive charge is far away from the structure, and Lagrangian description of motion for near field.
Different approaches have been explored in recent years to simulate the loading conditions generated by mine blast.This article briefly describes these approaches that include the following models: 1.
The pressure-based mine loading model 4.
The CHINOOK method 5.
The smoothed particle hydrodynamics (SPH) method.Each above techniques has limitations and therefore empirical and simplified methods using CONWEP, Westine, pressure-based mine loading models are based on analytical pressure validated for target, defense industry continues to investigate and develop additional features to predict mine blast and also to model the interaction of blast on complex structures.
The aim of this article is to compare numerical results obtained using the two different formulations, the ALE and SPH with empirical and simplified methods using CONWEP, Westine, pressure-based mine loading models.with methods experimental data.The comparison was carried out for two scenarios and the obtained numerical results were compared with the experimental data.The first scenario, which represents a typical buried mine blast, had an aluminum plate placed on four steel legs and centred over a surrogate mine filled with explosive type C4.The experimental setup, the FE models and the results are all given in Section 2. For this more complex scenario, several parametric studies were carried out for the ALE and SPH methods.For the ALE model, mesh sensitivity analysis of soil and explosive were investigated.A moving and a fixed air domain has been compared and the effect of the explosive geometry has also been studied.The effect of the properties of air material and equation of state were also studied.A mesh sensitivity analysis was conducted for the ALE model, and the results obtained with several formulations of the particle approximation theory were compared.In this article, devoted to ALE and Lagrangian formulations for soil structure interaction problems, the numerical and mathematical implementation of the ALE and SPH formulations are described.To validate the statement,, we perform a simulation of a shock wave propagation generated by explosive detonation.In Section 2, the governing equations of the ALE formulation are described.In this section, we discuss the advection algorithms used to solve mass, momentum and energy conservation in the multi-material formulation.Section 3 describes the SPH formulation, which can also been used for these applications, unlike ALE formulation which based of the Galerkin approach, SPH is a collocation method.The first target is an aluminum plate centred over a surrogate mine filled with C4 explosive which represents a typical buried mine blast scenario.The final deformation of the plate was measured for both approaches and was compared with experimental measurements.Parametric studies were conducted on both models and the best obtained results were compared to the experimental ones.Each comparison include the velocity at the center of the plate, sponson top and sponson sidewall.

ALE MULTI-MATERIAL FORMULATION
A brief description of the ALE formulation used in this paper is presented, additional details can be provided in Aquelet et al. (2005).To solve soil structure interaction problems, a Lagrangian formulation is performed for the structure and an ALE formulation for the soil and explosive materials, where soil and explosive materials can be mixed in the same element; this element is referred as mixed element, since it contains two different materials soil and explosive as described in Fig. 1.A mixture theory is used to partition the material inside the element and compute the volume weighted stress from the constitutive model of each material as described in Souli et al. (souli ,erchqui 2012).
In the ALE description, in addition to the Lagrangian and Eulerian coordinates an arbitrary referential coordinate is introduced.The material derivative with respect to the reference coordinate may be expressed in equation 1.Thus substituting the relationship between the reference configuration time derivative and material time derivative leads to the ALE equations, where i X is the Lagrangian coordinate, i x the Eulerian coordinate, i w is the relative velocity.The velocity of the material is v and the velocity of the mesh is u.In order to simplify the equations the relative velocity is introduced by w v u = − .Thus, the governing equations for the ALE formulation are given by the following conservation equations: ij σ is the stress tensor defined by .P Id σ τ = − + , whereτ is the shear stress from the constitutive model, and P the pressure.For explosive gas the pressure is computed through an equation of state defined in chapter.
For the structure, a classical elasto-plastic material model is used, where the shear strength is much higher than the volumetric strain.
Note that the Eulerian equations commonly used in fluid mechanics are derived by assuming that the velocity of the reference configuration is zero, 0 u = and that the relative velocity between the material and the reference configuration is the material velocity, w v = . The term in the relative velocity in Eq. 2.3 and Eq.2.4 is usually referred to the advective term, and accounts for the transport of the material past the mesh.An additional term, in the equations, makes solving the ALE equations much more difficult numerically than the Lagrangian equations, where the relative velocity is zero.
There are two methods to implement the ALE equations, which correspond to the two approaches taken in implementing the Eulerian viewpoint in fluid mechanics.The first method is to solve the fully coupled equations for computational fluid mechanics; this approach, used by different authors, can only deal with a single material in an element as described in an example given by Ozdemir, Souli and Fahjan (2010).The alternative method is referred to an operator split in the literature, where the calculation for each time step is divided into two phases.First, a Lagrangian phase is carried out, in which the mesh moves with the material, in this phase the changes in velocity and internal energy due to the internal and external forces, are calculated.The equilibrium equations are: In the Lagrangian phase, mass is automatically conserved, since no material flows across element boundaries.
In the second phase, the advection phase, transport of mass, energy and momentum across element boundaries are computed; this may be thought of as remapping the displaced mesh at the Lagrangian phase back to its original for Eulerian formulation or arbitrary position for ALE formulation using smoothing algorithms.From a discretization point of view of Eq. 2.5 and Eq.2.6, one point integration is used for efficiency and to eliminate locking as it is mentioned by Benson (1992).The zero energy modes are controlled with an hourglass viscosity, see Hallquist (1998).A shock viscosity with linear and quadratic terms derived by Von Neumann and Richtmeyer (1950), is used to resolve the shock wave.The resolution is advanced in time with the central difference method, which provides a second order accuracy for time integration.
For each node, the velocity and displacement are updated as follows: ) .( .
Where int F is the internal vector force and ext F the external vector force associated with body forces, coupling forces, and pressure boundary conditions, M is a diagonal lumped mass matrix.For each element of the mesh, the internal force is computed as follows: where B is the gradient matrix and Nelem is the number of elements.The time step size t ∆ , is limited by the Courant stability condition (see Benson (1992)), which may be expressed as: where l is the characteristic length of the element, and c the speed of sound through the material in the element.For a solid material, the speed of sound is defined as: where ρ is the material density, K is the module of compressibility.

SPH FORMULATION 3.1. SPH Formulation
The SPH method developed originally for solving astrophysics problem has been extended to solid mechanics by Libersky et al. (1993) to model problems involving large deformation including high velocity impact.SPH method provides many advantages in modeling severe deformation as compared to classical FEM formulation which suffers from high mesh distortion.The method was first introduced by Lucy (1977) and Gingold and Monaghan (1977) for gas dynamic problems and for problems where the main concern is a set of discrete physical particles than the continuum media.The method was extended to solve high velocity impact in solid mechanics, CFD applications governed by Navier-Stokes equations and fluid structure interaction problems.It is well known that SPH method suffers from lack of consistency that can lead to poor accuracy of motion approximation.Unlike Finite Element, interpolation in SPH method cannot reproduce constant and linear functions.
where the Dirac function satisfies: The approximation of the integral function Eq. 11 is based on the kernel approximation W, that approximates the Dirac function based on the smoothing length h.
that represents support domain of the kernel function, see Fig. 3. Taking in consideration de support domain of the kernel function, the SPH approximation of a particle   is obtained discretizing the integral into a sum over the particles that are within the kernel support domain as it is shown in Fig. 3.

(
) where the weight is the volume of the particle.
Integrating by part Eq. 14 and considering the properties of the SPH interpolation and that

( ) ( ) ( )
. 1 1. u u u ∇ = ∇ − ∇ , the SPH approximation for the gradient operator of a function is given by, ( ) �, applying the SPH interpolation on Navier-Stokes equations, one can derive a symmetric SPH formulation for Navier-Stokes equations such that the principle of action and reaction is respected and that the accuracy is improved.Finally, we have the following discretized set of equations: (ii) Momentum equation.( )

( )
For constant and linear function, The standard SPH interpolation is not exact: For It is well know from previous studies (see Villa (1999Villa ( , 2005))), that Eq. 20 and Eq.21 are exact only if the condition ∆ ℎ → 0 is fulfilled

Constitutive models and equation of states for explosive and water
In High explosive process, a rapid chemical reaction is involved, which converts the material into high-pressure gas.From a constitutive material point of view, the gas is assumed inviscid with zero shear, and the pressure is computed through JWL equation of state (Jones-Wilkins-Lee), a specific equation of state, commonly used for explosive material.There have been many equations of state proposed for gaseous products of detonation, from simple theoretically to empirically based equations of state with many adjustable parameters.
The explosive was modeled with 8-nodes elements.The equation of state determines the relation between blast pressure, change of volume and internal energy.The JWL equation of state was used in the following form: In Equation (5.1) p is the pressure, V is the relative volume: where v and vo are the current and initial element volume respectively, while A, B, C, R1, R2 and ω are material constants defined in table 1.These performance properties are based on the cylinder expansion test in controlled conditions.The first term of JWL equation, known as high-pressure term, dominates first for V close to one.The second term is influential in the JWL pressure for V close to two.Observe that in the expanded state, the relative volume is sufficiently important so that the exponential terms vanish, and JWL equation of state takes the form of an ideal gas equation of state:

Constitutive material model for soil
Geo-mechanic materials like soil and rocks behave very differently from classical elastic plastic which are ductile materials.When soil is compressed uniformly, it can withstand important loads, but during shear loading or tension, the material gets quickly into rupture.Unlike ductile materials soil acts as a brittle material and show little or no plastic deformation before fracture.Material characteristic is a critical element of numerical analysis.It can greatly influence the outcome of a numerical prediction.Many constitutive models are available to simulate soil behavior.In this paper we use a model due to which provides a simple model for soil.This model has been extensively used and validated by NASA Langley Research, where fundamental testing research for soils characterization has been performed.Soil material and many geo-materials can vary from the hard clay to the soft dry sand, and they are highly compressible.Because the material strength is pressure dependent, a volume pressure curve is necessary for constitutive modeling.
These materials typically fail in shear, where the shear failure surface is pressure dependent.When the shear failure surface is exceeded, the deviatoric stresses are limited by the failure surface and the material can then flow like a viscous fluid.If the shear strength of a material is very low as for dry sand, the model gives fluid-like behavior.Testing to determine the shear failure envelope, which depends on the confining pressure, is generally accomplished in the laboratory using a tri-axial compression test.The tri-axial test apparatus can be used for strength testing of soils, compressibility testing, and can be used to determine bulk unloading modulus by a hydrostatic compression test where the soil is loaded using a uniform pressure at the surface.In soil the shear failure surface is pressure dependent, which is a basic property of geo-materials.The pressure dependent shear surface is written in terms of pressure in the following form: where  2 = 1 2   .  is the second deviatoric invariant, p the mean pressure, and   is the deviatoric stress tensor.The coefficients  0 ,  1 ,  2 are determined from tri-axial compression tests.The material parameters for the soil material used in this paper are defined in the following table: -5.e-5 In table 1, the elastic shear modulus describes shear deformation when the soil is initially loaded.The bulk unloading modulus describes the expansion of the soil when the load is reduced.These two parameters are necessary because the loading and unloading behavior of soil is not equal due to permanent deformations.
Volumetric strain behavior is defined by the natural log of the relative volume and is negative in compression.Relative volume is the ratio of the current soil cell volume to the initial volume at the start of the calculation.The volumetric strain is represented as a 10point curve in pressure versus volume strain space.Each point on the curve is obtained from material testing at the given pressure.Figure 2 presents the scheme of the pressurevolumetric strain relationship.Starting from an initial volume (V=V0), the medium first undergoes elastic compression.An increase in pressure beyond the elastic part leads to plastic compression until full compaction is reached.An increase in pressure beyond leads to non-linear elastic compaction.Unloading process occurs according to an elastic law with the bulk unloading given in table 1.
The soil material used in this paper is a typical one used for most applications where material properties are defined in Table 2 For the structure, a classical elasto-plastic constitutive material law is used, where material properties are given in the following table : 4.2.Constitutive material for the plate A plastic kinematic material model was used to model the plate and the frame, and the corresponding material properties are given in Table 3.

COMPARISON BETWEEN ALE AND LAGRANGIAN METHODS
The problem setup is described in figure 4, where high explosive, soil and structure are meshed using 8 nodes hexahedra.The classical Lagrangian method can be used for more applications where the mesh is not highly distorted.As we can see in figure 5 after 220 microseconds, we observe high distortion in the soil mesh.If we carry the computation further in time, the soil will undergo higher mesh distortion which can lead to element negative volume and termination before physical termination time, this effect can be observed in figure 6 and 7 where pressure is plotted for both formulation.This situation does not appear in the ALE formulation, since a re-meshing procedure of the soil material is applied at each time step or after few cycles.The re-meshing procedure is described in detail Aquelet et al. (2005) To show the robustness of the ALE method, we run the same problem using ALE and Lagrangian formulations, where the Lagrangian method is stopped a soon as the soil material

Plastic compression
Full compaction starts having high mesh distortion.Fig. 8 shows good correlation for structure displacement, between the Lagrangian and ALE approaches.
In this work, we have presented the application of Lagrangian and ALE approaches for simulating blast wave propagation in water.Comparisons with experimental results from

CONCLUSIONS
In this paper, we present ALE and Lagrangian methods as well as their limitations for specific problems.Underwater explosion is commonly solved using ALE formulation, in defence industry; some of these problems are solved using SPH method.For the last decade, SPH methods are gaining in accuracy numerical stability, and the use of SPH method is becoming more common in industry for solving soil structure coupling problems.For instance, in aerospace, where bird impacts on aircraft are very common and cause significant safety threats to commercial and military aircraft.According to FAA (Federal American Aviation) regulations, aircraft should be able to land safely Souli et al (2012).For decades engineers in aerospace industry were using ALE method to simulate bird impact on aircrafts, where a viscous hydrodynamic material is used for the bird.These applications require a large ALE domain for the coupling between the bird material and the surrounding structure, mainly when the bird is spread all over the space.According to technical reports from engineers in aerospace, ALE formulation is more CPU time consuming and requires more memory allocation that SPH method.In this paper, first, we describe both ALE and Lagrangian methods, and we compare numerical results between the two methods using similar mesh size.Using a simple soil structure interaction problem, it has been observed that using same mesh size for both methods, numerical results, displacement, velocity and Von Mises stress on the structure, are under estimated with SPH method.When refining the SPH particles, where each ALE element is replaced by 4 SPH particles in two dimensional and 8 particles in three dimensions, numerical results form SPH method are in good correlation with those from ALE simulation; in terms of displacement, velocity and Von Mises stress on the structure.Since the ultimate objective is the design of structure resisting to load blast, numerical simulations from ALE and SPH methods can be included in shape design optimization with shape optimal design techniques, see Souli et al (1993), and material optimisation, see Erchiqui et al (2007).Once simulations are validated by test results, they can be used as design tool for the improvement of the system structure being involved.Once simulations are validated by test results, it can be used as design tool for the improvement of the system structure involved.

Figure 1 :
Figure 1: Lagrangian and Advection phases in one step i

Figure 3 :
Figure 3: Kernel Function and its support domain for a 2D function

Figure 8 :
Figure 8: ALE and Lagrangian structure displacement

Table 1 :
Parameters used for C-4 explosive

Table 2 :
Parameters used for Soil Material

Table 3 :
Parameters used for Structure