A Coupled Finite Difference Material Point Method and Its Application in Explosion Simulation

X.X.CuiX.ZhangX.ZhouY.Liuand F.Zhang

1 Introduction

A high explosive(HE)explosion is characterized by a number of challenging behaviors including the highly pressurized product gas propagating into the quiescent surrounding air and the following fluid structure interaction(FSI)with the structures nearby[Zukas and Walters(1998)].So a lot of research have been carried out to study these kind of problems.

The traditional methods can be classified into Lagrangian method and Eulerian method based on the frame of reference[Benson(1992)].The Lagrangian method has been widely used for structural analyses because of its capability of modeling history-dependent material and tracking material interface.A common practice in modeling HE explosion problems is to discretize the structure by Lagrangian finite elements and then the explosion effects are taken into account by applying the pressure load on the structure surface.For examples,a plate under air blast loading was studied by Jacinto et al.[Jacinto,Ambrosini,and Danesi(2001)]and the spallation in reinforced concrete plates subjected to blast loading was studied by Xu et al.[Xu and Lu(2006)].However,the Lagrangian finite element method(FEM)suffers from mesh tangling which deteriorates its numerical accuracy and efficiency dramatically.

Recently,many meshless methods based on Lagrangian framework have been proposed as the alternatives for the traditional finite element methods,which have showed advantages for problems associated with large deformation.Among them,the smoothed particle hydrodynamics(SPH)[Lucy(1977);Liu and Liu(2003)]and material point method(MPM)[Sulsky,Chen,and Schreyer(1994);Sulsky,Zhou,and Schreyer(1995)]have been successfully applied to HE explosion problems.Detonations of HE in air and underwater were simulated by SPH[Liu,Liu,Zong,and Lam(2003);Liu,Liu,Lam,and Zong(2003b)].The numerical tests revealed the ability of SPH in modeling explosion problems with arbitrary charge shape and different orientations.Ma et al.[Ma,Zhang,Lian,and Zhou(2009)]proposed an adaptive MPM for simulating the HE explosion problems whilst Lian et al.[Lian,Zhang,Zhou,Ma,and Zhao(2011)]extended the MPM method to the explosively driven metal problems whose numerical results agreed well with the Gurney solutions.Zhang et al.[Zhang,Zou,VanderHeyden,and Ma(2008);Zhang,Ma,and Giguere(2011)]enhanced MPM to simulate a series of fl uid-structure interactions and multi-material interactions problems.Using MPM,Banerjee[Banerjee(2004)]simulated the fragmentation of cylinders due to explosively expanding gases generated by a high energy material inside the cylinders,and Hu et al.[Hu and Chen(2006)]studied the synergistic effects of blast and fragmentation on a concrete wall.Furthermore,a comparison study of MPM and SPH in modeling hyper velocity impact problems was conducted by Ma et al.[Ma,Zhang,and Qiu(2009)].These studies concluded that these two methods possesses a great potential for simulating large deformation FSI problems at high strain-rate.However,the particle methods require more computational resource and have difficulties to construct high order scheme when simulating the fluid in high explosive(HE)explosion problems.

In contrast to the Lagrangian method,the Eulerian method employes fixed meshes so it is not plagued by mesh distortion.They usually solve the fluid region in HE explosion problem together with approaches for tracking the material interfaces and the internal history variables,such as the Youngs interface reconstruction method[Youngs(1982)],level set method[Osher and Fedkiw(2001)]and fuzzy interface method[Ning and Chen(2004)].Ma et al.[Ma,Wang,and Ning(2008)]developed a multi-material Eulerian hydrodynamic code with modified Youngs’interface reconstruction algorithm for the simulations of explosion problems such as explosion in tunnel and steel shaped charge jet.Luccioni et al.[Luccioni,Ambrosini,and Danesi(2004)]employed AUTODYN to study the structural failure of a reinforced concrete building inflicted by an air blast load.The dispersion process of the HE products in air was simulated by the three-dimensional Euler FCT solver.Wu et al.[Wu and Hao(2005)]simulated the ground shock and air blast pressure generated from surface explosions using AUTODYN2D.Furthermore,newly developed methods based on Eulerian framework such as Discontinuous Galerkin Method(DGM)[Cockburn,Hou,and Shu(1990)]has been used to solve the gaseous detonation problems[Wang,Zhang,Shu,and Ning(2012)].These studies show sound ability of Eulerian method in solving wave propagation process in HE explosion problems.However,the Eulerian description limits its ability to handle the FSI problems by one single method,so it is usually coupled with a Lagrangian method to discretize the region of structures.

There are also some mixed methods which take the advantages of both Lagrangian and Eulerian descriptions.A well-known example is the arbitrary Lagrangian-Eulerian method(ALE)[Liu,Belytschko,and Chang(1986)].The major numerical difficulty of ALE is developing an effective and efficient mesh moving scheme for complicated 3D problems.Furthermore,the numerical diffusion and dissipation still exist in ALE method.A detailed review on Lagrangian,Eulerian and their mixed methods was presented by Benson[Benson(1992)].

Since the Lagrangian method possess advantages in simulating the structures with historical variables and the Eulerian method handle the fluid better,much effort has been devoted to couple these two type of methods so as to take advantage of each method to simulate the HE explosion and the relevant large deformation problems.Fairlie et al.[Fairlie and Bergeron(2002)]described a coupled methodology for simulating the surface-laid or buried charges explosions.In the methodology,the air and explosive were modeled in an Euler-FCT grid as a single ideal gas while the surrounding soil and complex targets were modeled by Lagrangian grid.Zhang et al.[Zhang and Xu(2007)]investigated a cylindrical shell loaded by blast wave from a central charge.Finite volume method(FVM)was used to model the HE in ALE framework and FEM was adopted to model the shell in Lagrangian framework.Guillkey et al.[Guilkey,Harman,and Banerjee(2007)]developed an approach for solving full-physics FSI problems using the Eulerian description(FVM)for fluids and the Lagrangian description(MPM)for solids.To simulate FSI problems with large deformations in the structure,Gilmanov et al.[Gilmanov and Acharya(2008)]developed an effective numerical method in which the hybrid immersed boundary method(HIBM)was employed to resolve complex boundaries for the fluid flow and MPM was coupled to resolve the structural stresses and deformation.The combined method was implemented in the framework of finite difference method(FDM).Flekkoy et al.[Flekkoy,Wagner,and Feder(2000)]introduced a “hybrid model”that permits a continuum description in one region to be coupled to an atomistic description in another region.The two regions were solved by FDM and molecular dynamics(MD)respectively.

Coupling between meshless methods and FEM are also carried out to simulate the problems with large deformation.Aktay et al.[Aktay and Johnson(2007)]developed a FEM/SPH coupling technique for high velocity impact(HVI)simulation of composite panels.In the technique,contact interfaces were employed to couple the discrete smoothed particles and finite elements which were employed to model the parts undergoing large and small deformation,respectively.Zhang et al.[Zhang,Sze,and Ma(2006)]developed an explicit material point finite element method for HVI.In their method,the momentum equations were solved on a prede fined regular grid in the severely deformed region and on FE mesh elsewhere.Lian et al.[Lian,Zhang,Zhou,and Ma(2011)]developed a coupled approach in which the bodies with large and mild deformation were discretized by MPM and FEM,respectively.The interaction between two bodies was handled by a contact method and the FE nodes on the contact interface were treated as special particles.To further improve the efficiency,Lian et al.[Lian,Zhang,and Liu(2012)]proposed an adaptive material point finite element method in which material domains were initially discretized into finite elements(FE).Depending on severity of the distortion or plastic strain being developed,some elements were adaptively converted into MPM particles during the solution process.

Most of the coupling between Lagrangian methods and Euler methods divide the computational domain by an interface between fluid and structure in solving FSI problems.Individual materials occupy distinct regions in space,with interactions occurring at the material interfaces.Because of the separated nature of the materials,the interface requires additional treatment and often introduce numerical error.

In this paper,a coupled finite difference material point(CFDMP)method is developed to model the 3D HE explosion and its interaction with the surrounding structures.Taking the advantage of handling the shock wave propagation,FDM is employed to simulate a large proportion of the fluid region,while MPM is employed in the FSI region which contains the structures and the fluid near the structures.Therefore,the interface between two computational regions is located in the same material region( fluid)and the interface effect could be significantly reduced.The material interface is located in the MPM region so that the fluid-structure interaction is solved in MPM region to fully take its sound ability for simulating history-dependent material and tracking the material interface.Hence,the region involved shock wave dispersion problem is simulated by FDM and the region involved history-dependent materials and FSI problems are simulated by MPM.The interaction between FDM region and MPM region are implemented by a“bridge region”which contains only one material.MPM provides the boundary condition for FDM region by mapping the value from background grid nodes to the fictitious points outside the boundary of FDM,while FDM provides the boundary condition for MPM region by mapping value from cell-centre points to MPM interface grid nodes.The transportation between the two computational regions is implemented by moving particles in the bridge region.The proposed scheme has been implemented in our 3D explicit material point method code,MPM3D,to simulate HE explosion problems.Several numerical examples are presented to validate the efficiency and accuracy of the proposed method.

The remaining part of this paper is organized as follows.Section 2 presents the governing equations and the numerical scheme in each computational region.A description of CFDMP and the numerical implementations are presented in Section 3.Then the material models employed are introduced in Section 4.Several numerical tests are given in Section 5,and the conclusions are summarized in Section 6.

2 Governing equations and schemes

The problem domain can be divided into two computational regions in space.FDM is employed to simulate the fluid region,while the FSI region is simulated by MPM.Since the primary materials and their properties are different in two regions,different governing equations and schemes are employed as follows.

2.1 Governing equations and scheme in FDM region

The dispersion process of detonation products to the surrounding air is a flow with strong discontinuity.Owing to the extremely high detonation and dispersion speeds,the explosion process is adiabatic.The detonation products and the surrounding air can be assumed to be inviscid and compressible,which can be described by the three-dimensional compressible Euler equations

The explicit three-dimensional scheme of fractional step FDM is outlined in a time step(fromnton+1)as below[Yanenko(1971)].Takexdirection as an example,adaptive artificial viscosity[Zhang(2010)]is used to avoid the non-physical oscillations near the shock wave which can be written as

whereηis a parameter to be adjusted empirically to meet the requirements for different problems or determined according to the time step?t,spatial step?xand sound speedcas

The fractional steps method[Yanenko(1971)]is introduced to split the three dimensional problem into three one-dimensional flow problems.To reduce the artificial affect introduced by the integration sequence,the splitting can be implemented as

2.2 Governing equations and scheme in MPM region

During fluid-structure interaction,the history variables are important to describe the behavior of the structure.Therefore,the updated Lagrangian description is employed for the continuum which is governed by the momentum equation

whereVis the current material domain,σijis the Cauchy stress,ρis the current density,fiis the body force density,¨uiis the acceleration.The weak form of governing equation(10)can be obtained from the weighted residual method as[Sulsky,Chen,and Schreyer(1994)]

whereEis the energy per unit initial volume,is the strain rate,is the deviatoric stress andprepresents the pressure.

In CFDMP method,these governing equations will be solved by MPM as described in existing literature[Ma,Hanan,Komanduri,and Lu(2012)].MPM is an extension of the FLIP particle in cell(PIC)method[Brackbill and Ruppel(1986)]in computational fluid dynamics to computational solid mechanics.As a preprocessing step,we define the background grid in the FSI region,and discretize the material region by a set of particles,see Fig.1.All the material variables including mass,position,velocity,strain and stress are carried by the particles.In each time step,the particles are rigidly attached to the background grid in which the momentum equation is solved in the framework of the standard finite element method.Then,the positions and velocities of all particles are updated based on the grid nodal velocities and accelerations.Afterward,the deformed grid is discarded and a new regular grid is used in next time step,and the initial grid nodal mass and momentum can be obtained from the mass and momentum of particles.Thus,complications associated with mesh distortion are avoided.In general,a fixed regular grid can be used throughout the computation.

Figure 1:Material point discretization

Since the material domain is discretized with a set of particles,the density can be approximated as

wherenpdenotes the number of particles;is the Dirac Delta function;mpis the mass andxpis the position of particlep.Since the masses are carried by the particles,the mass conservation is automatically satisfied in MPM.

Since the particles are rigidly attached to the computational grid,the displacement of particlepcan be obtained by mapping from their grid node valuesuIusing the standard finite element interpolation functions of the grid as

whereis the interpolation function of grid nodeIevaluated at the position of particlep.The 8-node hexahedron interpolation is used whose shape function is given by

Substituting(14)and(15)into the weak form(11)and using a lumped mass matrix lead to

is the grid nodal external force.In Eq.(20),hdenotes the thickness of the boundary layer used to calculate the surface integral.The grid nodal masses can be obtained by

From time stepnton+1,the momentum equation is integrated by

The velocity and position of particles are updated by mapping the increments from background grid nodes back to particles as

Before calculating the incremental strain and spin tensors,the updated velocities of the particles are mapped back to the grid nodes to update their velocities,namely

The incremental strain and spin tensors are calculated by(take three-dimensional problems for example)

Finally,the density and stress of particles are updated by

3 Coupled finite difference material point method

Fig.2 shows a typical HE air explosion problem.A HE charge is burned into gaseous products which disperse to the surrounding air and then interact with the structure.The whole region can be divided into a fluid region and a FSI region separated by the dash line.The traditional FDM is employed to simulate the dispersion process in fluid region.When the pressure of cell-centres near the region interface,i.e.the dash line in Fig.2,attains a prescribed threshold,arrival of the shock wave front is detected and the interaction region will be activated.The interaction process is simulated by MPM so that the history variables of structure can be recorded to characterize the material damage.The detailed governing equations and schemes for the two regions have been presented in Section 2.The interaction between FDM region and MPM region is provided as follow.

Figure 2:A typical HE explosion problem

3.1 Bridge region

As shown in Fig.3,the whole problem domain is discretized tomregular cells inxdirection(y,zdirections are also applicable).The size of FDM cells and MPM background cells are the same.The cells from 0 tokare FDM cells and the cells fromk?wtomare MPM background cells,wherewdenotes the number of cells in the bridge region in directionx.Cells fromk?wtokdefine the bridge region,in which the FDM cells are coincident with the MPM background cells.As shown in Fig.3,two different materials in the MPM region are marked by circles and triangles respectively.The material in the FDM region is the same as the fluid material in the MPM region,air is taken for example here,and the structures in the MPM region could have complex geometry because it is discretized by particles.The variables of fictitious points(hollow squares)outside the FDM region in cellk+1 is obtained by interpolating the background grid nodes while the variables of interface nodes(hollow circles)of the MPM region is adjusted by interpolating the centers of cellk?w?1 after being initialized by Eqs.(18)to(21).The transport between the two computational regions is implemented by moving the particles through the cell interface between the cellsk?w?1 andk?w.The detailed methods and equations will be described in the following subsections.

Figure 3:The computational region for CFDMP

3.2 Interface boundary condition for FDM region

To solve Eq.(1)in FDM region,the variables of fictitious points(k+1)outside FDM region as shown in Fig.3 should be given.The massand momentumof the FDM’s fictitious points in time stepncan be obtained by mapping the MPM grid nodal mass and momentum via the shape function,namely(take three dimensional problems for example)

The internal energyis calculated by adding all the particles’internal energy in

3.3 Interface boundary condition for MPM region

To solve the governing equations in MPM region,the variables of grid interface nodes between cellk?w?1 and cellk?wof MPM region(see Fig.3)should be adjusted to consider the effect from the FDM region.The FDM cell-centers with cell numberk?w?1 in directionxare considered as particles and take part in the mapping process from particles to background nodes as in Eqs.(18),(19)and(21).Therefore,the mass,momentum and internal force of the background grid interface nodes of MPM are adjusted by

where the first term is the same as those in Eqs.(18),(19)and(21)in the standard MPM,while the second term is the contribution from the FDM region.The subscript“c”denotes the cells in the FDM region which are connected to the MPM interface node being adjusted,ncis the number of the cells,andpncdenotes the pressure.

3.4 Transportation between FDM and MPM regions

After integrating the governing equations in FDM region,the variables have been updated for every cell,and the transportation between the two computational regions is carried out by moving the particles through the boundary of the MPM region.We take a pair of the interface cells as an example as shown in Fig.4.Both two cells are in the FDM region and the right cell,which is a bridge cell,is also in the MPM region.The interface between these two cells are the boundary of the MPM region.

Figure 4:The transportation between FDM and MPM

Assume that the time stepnhas been solved and the transportation is needed.We first get the density,velocity and pressure of the interface by interpolating the cell center values of these two cells

The fluxes of mass and momentum transfered during this time step can be calculated as

whereis the normal velocity of the interface determined by Eq.(41).

Referring to the work by Flekkoy[Flekkoy,Wagner,and Feder(2000)]on coupling FDM and molecular dynamics(MD),particles are generated in cellk?wto guarantee the conservation of mass and momentum.The newly generated particles have the same density,velocities and pressure as they are in the “donor cell”.To make sure the newly generated particles do not have extreme distinction with existing particles in mass to improve the stability,we determine their number and mass by

whereis the mass of the existing particles.We slightly adjuststo an integer multiple of 4 to distribute particles uniformly.The furthest particles’distance from the center of the interface are calculated by their velocity in this step as

Moreover,the internal energy of the newly generated particles are given by their EOS as(air taking for example)

whereis the pressure of the newly generated particles.Since the equation is linear about mass,the conservation of internal energy is protected here.The newly generated particles in same cells have single velocity,and the mass conservation is automatically protected,so kinetic energy fl ow into the MPM region is equal to that carried by the newly generated particles.Together,the energy conservation is protected in this process.

If the particles in cellk?wcross the boundary of MPM region in this step,the transportation is from MPM to FDM.The crossing particles will not take part in the computation of MPM any more,and their conservation variables are added to the cell they move into.So the conservation is also protected in this process.For example,the particlepmoves from cellk?winto cellk?w?1 in Fig.4,so the cell-center values of cellk?w?1 should be adjusted by

3.5 Numerical implementation

One explicit step using CFDMP(fromnton+1)is summarized below.

(1)Calculate the time step for MPM and FDM by CFL criterion respectively,and take the smallest one as the time step for CFDMP.

(2)Map the massmand momentumpof all MPM particles to the background grid except the interface nodes by Eqs.(21)and(18).

(3)Map the massmand momentumpof corresponding MPM particles to the interface nodes of background grid by Eqs.(37)and(38).

(4)Compute the grid nodal internal forcefintand external forcefextexcept the interface nodes by Eqs.(19)and(20).

(5)Compute the grid nodal internal forcefintof interface nodes by Eq.(39).

(6)Integrate the momentum equation by Eq.(22).

(7)Update the fictitious points’variables by Eqs.(34)to(36).

(8)Integrate the governing equations of FDM by Eq.(6),where the operators in three directions are defined by Eqs.(7)to(9).

(9)Update the velocity and position of particles by mapping their increments back to particles by Eqs.(24)and(25).

(10)Map the velocity back to the grid nodes by Eq.(26).

(11)Calculate the incremental strain and spin tensors by Eqs.(27)and(28).

(12)Update the density of particles by Eq.(29).

(13)Update the stress of particlesby Eq.(30).

(14)Carry out the transportation process between FDM region and MPM region as described in Section 3.4 and this completes the current time step.

4 Material models

Equations of state,constitutive models and reaction models complete the whole governing equations.Brief descriptions of the models used in this paper are given below.Some of the parameters are taken from the references as well as the commercial software such as AUTODYN and LS-DYNA.

4.1 High explosive EOS

In detonation process,the reactive wave propagates at very high speed inside the HE[Zukas and Walters(1998)].The exothermic reaction is completed within a few microseconds with the HE completely converted to gaseous products.Most of the earlier works use the “artificial detonation model”[Liu and Liu(2003)]which considers the explosive as a group of gaseous products with the same energy and volume of the initial explosive charge.For most simulations in this paper,we use the “real detonation model”[Liu and Liu(2003)]which lights the explosive according to the reactive wave’s propagation,the pressure jump which occurs when the shock front arrived at material interface is better captured[Cui,Zhang,Sze,and Zhou(2013)].For saving computational resources and accelerating the simulation,we refer to the remap method in AUTODYN for air explosion problem,solving an 1D TNT explosion problem first and map the result to the 3D region as the initial condition in FDM region.Finer grid can be allocated for the 1D simulation,which is in favour of describing the strong discontinuity during the detonation process and the initial stage of the dispersion process.

We simulate the 1D TNT explosion by MPM so as to conveniently model the detonation process by “real detonation model”.In the initialization phase,a lighting timetLis calculated for each particle(MPM)by dividing the distance from the detonation point by the detonation speed.After the detonation,the gaseous products are controlled by the EOS.The real pressurepof the gaseous products is determined by multiplying the pressurepEobtained from EOS with a burn fractionFthat controls the release of chemical energy[Hallquist(1998)],namely

wherehis the characteristic size of a particle andtdenotes the current time.Several time steps are often required forFto reach unity.Once it is done,Fis kept at unity.Using this method,the discontinuous detonation wave is smoothed and assumes a continuous but rapidly changing wavefront.

After detonation,the gaseous products are described by Jones-Wilkins-Lee(JWL)EOS

Moreover,TNT with a density of 1630kg/m3and a detonation speed of 6930m/s are used in the simulation.The parameters of JWL EOS are taken from[Liu,Liu,Lam,and Zong(2003a)]asA=3.712×1011N/m2,B=3.21×109N/m2,R1=4.15,R2=0.95,ω=0.3,energy per initial volumeE0=6993×106J/m3.

4.2 Air model

Air is modeled as a null material model with the following ideal gas EOS

whereρ=1.225kg/m3ande=2.0685×105J/kg.

4.3 Concrete model with tensile damage

The concrete is modeled by Holmquist Johnson Cook(HJC)model with tensile damage.The HJC model was originally presented for concrete damage problems involving hydrostatic pressure,strain rate and compressive damage.The equivalent strength is expressed as

wheredenotes the normalized equivalent stress,σis the actual equivalent stress,represents the quasi-static uniaxial compressive strength.denotes the normalized pressure,pis the real pressure.represents the dimensionless strain rate,˙εis the real strain rate andis the reference strain rate.A,B,N,CandSmaxare normalized cohesive strength,normalized pressure hardening coefficient,pressure hardening exponent,strain rate coefficient and normalized maximum strength,respectively.Dis an index describing the material damage in the range of 0～1.According to the original HJC model[Holmquist,Johnson,and Cook(1993)],an accumulated damage failure model,also known as compression-shear damage,is considered,which is written as

whereDcdenotes the compression-shear damage parameter,△εpand△μpdenote the equivalent plastic strain and plastic volumetric strain,respectively.D1andD2are the damage constants.In order to allow for a finite amount of plastic strain to fracture,a third damage constantEfminis provided.T?=T/f′cdenotes the normalized maximum tensile hydrostatic pressure.

However,the tensile damage of concrete is not considered in the original HJC model.The tensile behavior of concrete is simply considered through maximum tensile hydrostatic pressure.A new brittle tensile failure model based on micro cracking growth of concrete was presented by Jiang similarly as the metal brittle tensile failure[Jiang(2010)].According to this model,every crack can be viewed as a sphere cavity zones with the maximum diameter of crack which is called"equivalent micro holes".This model can be formulized as

whereais damage factor of micro-crack growth,or frequency of micro-crack growth,σ0is threshold stress of damage development involving micro voids’nucleation and growth,γis dependent coefficient of the ultra threshold stress.When the tensile damage reaches the limit of damage,spall of material will occur.Considering high pressures and air voids,the equation of state(EOS)in HJC model is divided into three response regions including linear elastic zone,transition zone and full dense zone.More details can be found in papers[Lian,Zhang,Zhou,and Ma(2011);Holmquist,Johnson,and Cook(1993)].

4.4 Soil model

The soil in this work is modeled by Drucker-Prager constitutive model[Itasca(2005)].It is made up by shear failure and tension failure.For judging the shear failure region and tension failure region,functionh(σm,τ)is defined as

whereσtis the tensile strength,kφandqφare material constants which can be obtained from the cohesion and frictional angle.Whenh＞0,shear failure is employed and the yield function can be described as

Ifh＜0,tension failure is employed and yield function can be described as

The parameters of soil are taken from[Luccioni and Luege(2006)]asρ=1200 kg/m3,E=100MPa andε=0.3.The cohesion is 0.11MPa and the internal friction angle is 20°.

4.5 Steel model

The Johnson-Cook material model[Johnson and Cook(1983)]is employed in the numerical example to describe the property of the steel plate.The model accounts for the strain rate effect and has widely used to model the behavior of metal during impact and explosion.The yield stress is given by

whereεis the equivalent plastic strain,is the dimensionless plastic strain rate withis the dimensionless temperature.Tis the material’s temperature,Troomis the room temperature,andTmeltis the material’s melting temperature.The material constants are taken from the reference[Neuberger,Peles,and Rittel(2007)]to beA=950MPa,B=560MPa,n=0.26,C=0.014 andm=1.

The pressure of steel is updated by the Mie-Gr??neisen EOS as

The subscript H refers to the Hugoniot curve andμ=ρ/ρ0?1 is used to represent the compression of solid withρ0being the stress-free solid density.Moreover,γ,C0andSare the material constants which are taken asγ=2.17,C0=4569m/s andS=1.49 for the numerical example in this paper.

5 Numerical examples

The three dimensional CFDMP scheme has been implemented in our 3D explicit material point method code,MPM3D[Ma,Zhang,and Huang(2010)],to solve HE air explosion problems.Five numerical examples are presented in this section to validate the scheme and demonstrate its capabilities.

5.1 One-dimensional Shock tube problem

Sod shock tube problem[Sod(1978)]is a benchmark for validating codes for compressible fluid,so it is taken to demonstrate the FDM solver in CFDMP.As shown in Fig.5,this problem consists of a shock tube with a diaphragm separating two regions whose initial states areρleft=1.0g/mm3,pleft=1.0MPa,ρright=0.125g/mm3andpright=0.1MPa.The fluids in both regions are initially at rest.At timet=0ms,the diaphragm is ruptured.Then,the shock and the contact interface travel at different speeds.The results are usually examined att=0.143ms when the shock has traveled a distance of about 0.25mm.This problem is employed to test the capability of the FDM solver in CFDMP on simulating compressible fluid and does not involve coupling between FDM and MPM.The pro fi les of density,velocity and pressure are plotted in Fig.6 for a grid with 1000 cells,which shows that the FDM solver’s results are in excellent agreement with the analytical results.The simulation results obtained by MPM using the same cells are also plotted in the figures in which obvious numerical oscillations can be noted.Generalized interpolation material point method(GIMP)[Bardenhagen and Kober(2004)]can effectively inhibit the numerical oscillations and get better results than MPM,however,more computational resources are needed.Unlike the MPM and GIMP in which both the particles and background grid are created,the FDM solver in CFDMP method updates the variables only in the cell-center points.In this regard,the CPU times consumed by GIMP and MPM are 78s and 52s,respectively,while FDM only takes 46s.Furthermore,the convergence properties of FDM and GIMP are studied by plotting the global error norms of the results against the background cell length(h),as shown in Fig.7.The convergence rate of FDM is about 50%higher than that of GIMP.What’s more,the global errors of FDM using 500 cells are almost equal to that of GIMP using 1000 cells,which demonstrates the rationality of employing FDM to simulate the fluid region which contain shock wave propagation in HE problems in CFDMP.

Figure 5:1D shock tube problem

5.2 Two-dimensional HE explosion and interaction with a concrete slab

Figure 6:Pro files of density,velocity and pressure obtained by analytical solution,MPM,GIMP and FDM

Figure 7:Pressure convergence curves

A two-dimensional HE explosion problem is simulated as shown in Fig.8.A HE charge with a radius of 50mm detonates and drives the surrounding air to interact with a concrete slab.The HE and surrounding air are simulated by the high explosive model and air model presented in Section 4 and the boundary conditions of the fluid region are all “ fl ow out”.The concrete is simulated by the concrete model with tensile damage presented in Section 4.The computational domain(?180,160mm)×(0,500mm)is divided into FDM region(Fluid region)of(?180,120mm)×(0,500mm)and MPM region(Interaction region)of(110,160mm)×(0,500mm),as shown in Fig.8.All regions are discretized by square cells of side length 2 mm.The bridge region is(110,120mm)×(0,500mm)andw=5,i.e.,the width of the bridge region is equal to 5 times of the cell’s side length.The center of the HE charge is located at(0,250).

Fig.9 shows the colored contours of the pressure at 20μs in FDM and MPM region respectively.It can be recognized that the wave propagation transport through the bridge region correctly.To quantitatively demonstrate that CFDMP does not introduce obvious interface effect,the pressure time history of two symmetric gauge points in the air,which have the same distance from the center of the HE charge(see Fig.8),are shown in Fig.10.The gauge 1(?125,250)is in the FDM region while the gauge 2(125,250)is in MPM region.The pressure time histories of the two gauge points are in good agreement,which demonstrates that the proposed coupling method between FDM and MPM works very well.

To further study the effects of the width of the bridge region,Fig.11 plots the pressure time histories forwequals to 2 and 8,respectively,which shows that the bridge region would smooth the shock wave excessively whenwis too large.Therefore,the width of the bridge region should not be too large,andw=5 performs very well as shown in Fig.10.

The pressures in the interaction region obtained by MPM and CFDMP are compared in Fig.12.Not only the time step used in MPM simulation is much smaller than that of CFDMP due to some particles with extraordinary high speed near the contact discontinuity,but also the MPM result suffers more numerical oscillations.Besides,MPM requires more memory because it places particles in the whole region.As a result,the total computational time of MPM is 181 minutes 9 seconds,while the total computational time for CFDMP is only 50 minutes 30 seconds.

5.3 Three-dimensional HE explosion and interaction with a concrete slab

To validate the capability of CFDMP in simulating HE explosion problem and the damage effect to the structure nearby,a concrete slab under blast loads experiment[Luccioni and Luege(2006)]is studied.The geometrical con figuration of the experiment and the load locations are shown in Fig.13.Spherical HE charges of 5Kg(TNT equivalent mass of 4Kg)and 12.5Kg(TNT equivalent mass of 10Kg)of Gelamon VF80 were employed placed at 0.5m height above the top surface of the slab as shown in Fig.13(a).There were three tests on the same slab successively,the first produced a fracture parallel to the short side,so for the following deto-nations,the original slab behaved as two square independent slabs so we simulate two tests independently.The average compressive strength of the concrete(25MPa)was obtained from compression tests with the same concrete as the slab.In[Luccioni and Luege(2006)],AUTODYN was used to simulate the experiments,and their material parameters are used in our simulations.In both AUTODYN and our simulations,a remap method is used to map the 1D simulation results of the detonation process to 3D region as the initial condition in the FDM region and then simulate the propagation of blast wave in air by 3D code.

Figure 12:Contours of pressure in interaction region at t=32.5μs.(a).MPM;(b).CFDMP

Figure 13:(a).Concrete slab dimensions and charge positions;(b).Placement of the explosive charge suspended above slab[Luccioni and Luege(2006)]

In the experiment,the charges of 5Kg and 12.5Kg Gelamon produced the crushing of concrete in a circular zone of about 250mm and 300mm diameter,respectively.Hereby,Luccioni et al.[Luccioni and Luege(2006)]over fitted an empirical formula for the estimation of explosive charges from crushing diameters or crushing dimensions from explosive charges as

whereDis the diameter of the crushing zone andhis the height of the charge from the concrete slab.Wdenotes the TNT equivalent mass of the charge.

In order to reproduce the crushing or disintegration of the concrete and counteract the great distortion that can cause excessive deformation of the mesh,erosion was used in[Luccioni and Luege(2006)].Pro fit from MPM’s ability to deal with the large deformation problems,we use CFDMP to simulate this test without use of the erosion model.The particles fail when their damage value reaches unit.Fig.14 shows the numerical results for the test of 5Kg charge obtained by CFDMP,in which(a)is the failure area on the front of the slab and(b)is the tensile damage on the back of the slab.The diameter of the crushing area is about 261mm and the damage on the back of the slab expands from the center to the surrounding radially which conforms to the analysis in[Luccioni and Luege(2006)].

Figure 14:Numerical results for 5Kg charge test.(a).Concrete crushing area on the front of the slab;(b).Tensile damage on the back of the slab

To compare the CFDMP results with those obtained from the experiment and the empirical formula(67),we conducted a series of numerical tests with TNT equivalent mass of 2Kg,4Kg,10Kg and 20Kg.The results are plotted in Fig.15.For the cases of 4Kg and 10Kg,the numerical results are in good agreement with those obtained by the experiment and the empirical formula.For the case of 2Kg and 20Kg,the crushing regions also agree well with the prediction of the empirical formula.

Figure 15:Relative diameter of the crushing zone D/h as a function of the inverse of the scaled distanceW1/3/h

5.4 Response of steel plates subjected to air-blast loading

Understanding the dynamic behavior of blast loaded armor steel plates is a key to design a protection structure successfully.A series of experiments and numerical calculations were carried out by Neuberger et.al.[Neuberger,Peles,and Rittel(2007)]to study the scaling characteristics of the dynamic response of circular RHA steel plates to large bare spherical air blast.As shown in Fig.16(a),the target plate was supported by two thick armor steel plates with circular holes.The spherical TNT charges were hanged in air using fisherman’s net and were ignited from the center of the charge.The numerical model is shown in Fig.16(b),in which the charge’s weight isW,distance from the plate’s surface to the center of the charge isR,plate thickness istand diameter isD,all of them are variables in the tests.The maximum deflection of the plateδis measured during the experiments.Steel,TNT and air are simulated by the material models presented in Section 4.Two series of cases are simulated,and the experimental and numerical parameters are listed in Tab.1.The normalized peak deflection,δ/t,obtained by CFDMP are compared with those obtained by the experiments and LS-DYNA in Tab.2.

Figure 16:(a).Experimental setup[Neuberger,Peles,and Rittel(2007)];(b).The numerical model

For series 1(cases 1,2 and 3),which represents a charge of 30Kg TNT for the full scale prototype,the structural response is mostly dynamic elastic.The first two cases show that the experimental and numerical results are in excellent agreement and the scaling are hardly affected by the varying the geometry scale as stated in[Neuberger,Peles,and Rittel(2007)].The third case is the numerical result of the full scale problem given by CFDMP,it proves the conclusion again.

For series 2(cases 4,5 and 6),which represents a charge of 70Kg TNT for the full scale prototype,the increasing plastic strains arise in the structural dynamicresponse.Thus the scaling is affected,as discussed in[Neuberger,Peles,and Rittel(2007)].Case 4 and case 5 show that the experimental and numerical results are in good agreement,and the effect on the scaling are also represented.As the scaling factor vary approaching to 1,the normalized peak deflectionδ/tincrease,which tells the effect on the scaling.The last case is the numerical result of the full scale problem obtained by CFDMP,whose normalized peak deflectionδ/tis the largest which conform to the conclusion.

Table 1:Experimental and numerical cases

Table 2:Normalized peak deflection δ/t

Fig.17 shows the final con figuration of case 5.The figure obtained by CFDMP is similar to the shape of experimental photograph.Fig.18 shows the colored contours of the Mises stress at time of 0.3ms,0.5ms and 1.0ms respectively for case 2 and case 5(the scaling factor are both 2).For case 2,most of the region are in elastic period,while a large portion has entered into plastic period for case 5.So the scaling of case 5 are affected as discussed in[Neuberger,Peles,and Rittel(2007)].

5.5 Damage of concrete slab with defect subjected to air-blast loading

Figure 17:Final con figuration of case 5.(a).Experimental photograph;(b).CFDMP

Figure 18:Mises stress.(a).case 2;(b).case 5

Based on the examples shown in Section 5.3,we bring a little modification to the slab to illustrate the ability of CFDMP to solve the strong FSI problem for more complex geometry structure.As shown in Fig.19(a),in a block region right against the center of the HE charge(the red region),the concrete slab is weakened by reducing its thickness from 0.15m to 0.02m so as to simulate the effect of the defect in the slab,other conditions are the same as those in Section 5.3(5Kg charge).The tensile damage in the back of the slab is presented in Fig.19(b).All the defect region is crushed and the damage region of the slab is larger than that of the slab without defect shown in Fig.14.Besides,the damage shape of the two case are also different.The damage shape of the slab without defect is expanded from the center to the surrounding radially while the damage shape of the slab with defect is annular around the center.

Figure 19:(a).Initial condition;(b).The tensile damage in the back of the slab

6 Conclusion

A coupled finite difference material point(CFDMP)method is proposed through the bridge region to combine the advantages of FDM and MPM in this paper.It uses an Eulerian frame in the fluid region and a Lagrangian frame in the FSI region.FDM is employed in the Eulerian frame,while MPM is employed in the Lagrangian frame.So the region involving shock wave dispersion problem is simulated by FDM whereas the region involving history dependent materials and FSI is simulated by MPM.In this way,the advantages of FDM and MPM are fully utilized in different regions of the problem.Both shock tube problem and 2D HE explosion simulation have verified the accuracy and efficiency of this algorithm.CFDMP is applied to simulate the dynamic responses of concrete slab and steel slab under air blast loading,whose numerical results coincide to the available experimental results and the conclusion reported in literature.Therefore the proposed CFDMP method provides a powerful numerical tool for the study of explosion problems.What’s more,the methodology can be generalized to combine other methods which are based on different frames of reference.

Acknowledgement:Supported by the National Basic Research Program of China(2010CB832701),National Natural Science Foundation of China(11390363),Specialized Research Fund for the Doctoral Program of Higher Education(SRFDP)of Chinese Ministry of Education(20130002110078)and Tsinghua University Initiative Scientific Research Program.

Aktay,L.;Johnson,A.F.(2007):FEM/SPH coupling technique for high velocity impact simulations.Advances in Meshfree Techniques,pp.147–167.

Banerjee,B.(2004):Material point method simulations of fragmenting cylinders.In17th Engeering Mechanics Conference.

Bardenhagen,S.G.;Kober,E.M.(2004):The generalized interpolation material point method.Computer Modeling in Engineering and Sciences,vol.5(6),pp.477–495.

Benson,D.J.(1992): Computational methods in Lagrangian and Eulerian hy-drocodes.Computer Methods in Applied Mechanics and Engineering,vol.99(2-3),pp.235–394.

Brackbill,J.U.;Ruppel,H.M.(1986): Flip:a method for adaptively zoned,particle-in-cell calculations in two dimensions.Journal of Computational Physics,vol.65,pp.314–343.

Cockburn,B.;Hou,S.;Shu,C.W.(1990): Tvb runge-kutta local projection discontinuous galerkin finite element method for scalar conservation laws iv:The multidimensional case.Math Comp,vol.54,pp.545–581.

Cui,X.X.;Zhang,X.;Sze,K.Y.;Zhou,X.(2013):An alternating finite difference material point method for numerical simulation of high explosive explosion problems.Computer Modeling in Engineering and Sciences,vol.92(5),pp.507–538.

Fairlie,G.;Bergeron,D.(2002):Numerical simulation of mine blast loading on structures.In17th Military Aspects of Blast Symposium.

Flekkoy,E.G.;Wagner,G.;Feder,J.(2000): Hybrid model for combined particle and continuum dynamics.Europhysics Letters,vol.52(3),pp.271–276.

Gilmanov,A.;Acharya,S.(2008): A hybrid immersed boundary and material point method for simulating 3d fl uid-structure interaction problems.International Journal for Numerical Methods in Engineering,vol.56,pp.2151–2177.

Guilkey,J.;Harman,T.;Banerjee,B.(2007):An Eulerian Lagrangian approach for simulating explosions of energetic devices.Computers and Structures,vol.85,pp.660–674.

Hallquist,J.O.(1998):LS-DYNA Theoretical Manual.Livermore Software Technology Corporation.

Holmquist,T.J.;Johnson,G.R.;Cook,W.H.(1993):A computational constitutive model for concrete subjected to large strains,high strain rates,and high pressures.In14th International Symposium on Ballistics Quebec,Canada.

Hu,W.;Chen,Z.(2006): Model-based simulation of the synergistic effects of blast and fragmentation on a concrete wall using the MPM.International Journal of Impact Engineering,vol.32,pp.2066–2096.

Itasca(2005):FLAC-Fast Lagrangian Analysis of Continua User’s Manual(Version 5.0).InItasca Consulting Group Inc.

Jacinto,A.C.;Ambrosini,R.D.;Danesi,R.F.(2001): Experimental and computational analysis of plates under air blast loading.International Journal of Impact Engineering,vol.25,pp.927–947.

Jiang,D.(2010):Study on chracterization of damage evolution and failure law of engineering materials.PhD thesis,University of Science and Technology of China,2010.

Johnson,G.R.;Cook,W.H.(1983):A constitutive model and data for metals subjected to large strains,high strain rates and high temperatures.In7th International Symposium on Ballistics,The Hague,Netherlands:American Defense Preparedness Association,pp.541–547.

Lax,P.D.;Wendroff,B.(1964): Difference scheme for hyperbolic equations with high order of accuracy.Communications on Pure and Applied Mathematics,vol.17,pp.381–398.

Lian,Y.P.;Zhang,X.;Liu,Y.(2012):An adaptive finite element material point method and its application in extreme deformation problems.Computer Methods in Applied Mechanics and Engineering,vol.241-244(1),pp.275–285.

Lian,Y.P.;Zhang,X.;Zhou,X.;Ma,S.;Zhao,Y.L.(2011): Numerical simulation of explosively driven metal by material point method.International Journal of Impact Engineering,vol.38,pp.238–246.

Lian,Y.P.;Zhang,X.;Zhou,X.;Ma,Z.T.(2011): A femp method and its application in modeling dynamic response of reinforced concrete subjected to impact loading.Comput.Methods Appl.Mech.Engrg.,vol.200,pp.1659¨C1670.

Liu,G.R.;Liu,M.B.(2003):Smoothed particle hydrodynamics:a mesh free particle method.World Scientific,Singapore.

Liu,M.B.;Liu,G.R.;Lam,K.Y.;Zong,Z.(2003):Mesh free particle simulation of the detonation process for high explosives in shaped charge unlined cavity con figurations.Shock Waves,vol.12(6),pp.509–520.

Liu,M.B.;Liu,G.R.;Lam,K.Y.;Zong,Z.(2003): Smoothed particle hydrodynamics for numerical simulation of underwater explosion.Computational Mechanics,vol.30,pp.106–118.

Liu,M.B.;Liu,G.R.;Zong,Z.;Lam,K.Y.(2003): Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology.Computers and Fluids,vol.32,pp.305–322.

Liu,W.K.;Belytschko,T.;Chang,H.(1986):An arbitrary lagrangian-eulerian finite element method for path-dependent materials.Computer Methods in Applied Mechanics and Engineering,vol.58,pp.227–245.

Luccioni,B.M.;Ambrosini,R.D.;Danesi,R.F.(2004):Analysis of building collapse under blast loads.Engineering Structures,vol.26,pp.63–71.

Luccioni,B.M.;Luege,M.(2006):Concrete pavement slab under blast loads.International Journal of Impact Engineering,vol.32,pp.1248–1266.

Lucy,L.A.(1977):A numerical approach to the testing of the fission hypothesis.Astrophysical Journal,vol.82,pp.1013–1024.

Ma,J.;Hanan,J.C.;Komanduri,R.;Lu,H.(2012): Simulation of the deformation mechanisms of bulk metallic glass(BMG)foam using the material point method.Computer Modeling in Engineering and Sciences,vol.86(4),pp.349–382.

Ma,S.;Zhang,X.;Lian,Y.P.;Zhou,X.(2009):Simulation of high explosive explosion using adaptive material point method.Computer Modeling in Engineering and Sciences,vol.39(2),pp.101–123.

Ma,S.;Zhang,X.;Qiu,X.M.(2009):Comparison study of MPM and SPH in modeling hypervelocity impact problems.International Journal of Impact Engineering,vol.36,pp.272–282.

Ma,T.B.;Wang,C.;Ning,J.G.(2008):Multi-material Eulerian formulations and hydrocode for the simulation of explosions.Computer Modeling in Engineering and Sciences,vol.33(2),pp.155–178.

Ma,Z.T.;Zhang,X.;Huang,P.(2010): An object-oriented mpm framework for simulation of large deformation and contact of numerous grains.Computer Modeling in Engineering&Sciences,vol.55(1),pp.61–88.

Neuberger,A.;Peles,S.;Rittel,D.(2007):Scaling the response of circular plates subjected to large and close-range spherical explosions.part i:Air-blast loading.International Journal of Impact Engineering,vol.34,pp.859–873.

Ning,J.G.;Chen,L.W.(2004):Fuzzy interface treatment in eulerian method.Science in China Series E-Engineering and Materials Science,vol.47(5),pp.550–568.

Osher,S.;Fedkiw,R.P.(2001):Level set methods:an overview and some recent results.Journal of Computational Physics,vol.169,pp.463–502.

Sod,G.A.(1978): A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.Journal of Computational Physics,vol.27,pp.1–31.

Sulsky,D.;Chen,Z.;Schreyer,H.L.(1994): A particle method for history-dependent materials.Computer Methods in Applied Mechanics and Engineering,vol.118(1-2),pp.179–196.

Sulsky,D.;Zhou,S.J.;Schreyer,H.L.(1995):Application of a particle-in-cell method to solid mechanics.Computer Physics Communications,vol.87(1-2),pp.236–252.

Wang,C.;Zhang,X.X.;Shu,Q.W.;Ning,J.G.(2012):Robust high order dis-contimuous Galerkin schemes for two-dimensional gaseous detonations.Journal of Computational Physics,vol.231,pp.653–665.

Wu,C.Q.;Hao,H.(2005):Modeling of simultaneous ground shock and airblast pressure on nearby structures from surface explosions.International Journal of Impact Engineering,vol.31,pp.699–717.

Xu,K.;Lu,Y.(2006): Numerical simulation study of spallation in reinforced concrete plates subjected to blast loading.Computers and Structures,vol.84,pp.431–438.

Yanenko,N.N.(1971):The method of fractional steps.Springer-Verlag.

Youngs,D.L.(1982):Time dependent multi-material flow with large fluid distortion.Numerical Methods for Fluid Dynamics,pp.273–285.

Zhang,D.L.(2010):A course in computational fluid dynamics.Higher Education Press.

Zhang,D.Z.;Ma,X.;Giguere,P.T.(2011):Material point method enhanced by modified gradient of shape function.Journal of Computational Physics,vol.230,pp.6379–6398.

Zhang,D.Z.;Zou,Q.;VanderHeyden,W.B.;Ma,X.(2008):Material point method applied to multiphase flows.Journal of Computational Physics,vol.227,pp.3159–3173.

Zhang,X.;Sze,K.Y.;Ma,S.(2006):An explicit material point finite element method for hyper-velocity impact.International Journal for Numerical Methods in Engineering,vol.66,pp.689–706.

Zhang,Y.J.;Xu,S.L.(2007):Numerical simulation on fl ow-structure interaction loaded by a blast wave from a central charge.Journal of University of Science and Technology of China,vol.37,pp.6–12.

Zukas,J.;Walters,W.(1998):Explosive effects and applications.Spinger.