This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.
Explosive blasts have hazardous effects on the society as well as the environment. Many recent events involving blasts have displayed the need for building protective structures to reduce the amount of damage caused. A large scale physical testing is not possible for building resistive structures; hence numerical methods serve as a good alternative by providing a virtual environment. LS-Dyna code has been implemented for the analysis since it has capabilities to simulate high impact as well as Fluid Structure Interaction (FSI) problems. The blast phenomena can be simulated using ConWep airblast function or by ALE (Arbitrary Eulerian Lagrangian Formulation) techniques. In this project ALE techniques have been used to simulate the event of explosion and to study the effect on a structural plate.
The objective of this project is to study the effect of TNT explosion on a Lagrangian plate fixed at a certain distance from the point of explosion. This project also aims to setup a methodology for computationally simulating blast events using ALE formulation such that it can be further extended in complex problems.
In this project the Computational Analysis of Spherical Air Blast has been done using the ALE. Later the Fluid Structure Interaction between the Fluid and a Plate has been studied. The Plate has been modeled using Lagrangian formulation. The explosive used is TNT, and it has modeled using the 'High Explosive Burn' material model and 'JWL' Equation of State. The Air domain is modeled by using a 'Null' material model and the pressure relationship is defined by using the 'Linear Polynomial' Equation of State. The modeling is discussed in detail in the sections ahead.
2] NUMERICAL ANALYSIS
A] Explicit Dynamics
LS Dyna employs an algorithm based on modified Central Difference Time Integration (CDTI) method, which is an explicit scheme. Since only a numerical solution is possible for a non linear ordinary differential equation, this method is particularly suited for non linear problems. It requires the inversion of the lumped mass matrix as opposed to that of the global stiffness matrix in the implicit methods. In the CDTI, the equation of motion is evaluated at the previous time step (tn-1, where tn is the current timestep).
The Courant-Friedrichs-Lewy (CFL) condition is essential for the stability of the solution. The CFL is a condition in numerical equation, solving which states that, given a space discretization, a time step bigger than some computable quantity should not be taken. The condition can be viewed as a sort of discrete "light cone" condition, namely that the time step must be kept small enough so that information has enough time to propagate through the space discretization . In structural mechanics, it implies that the minimum timestep size should be small than the critical timestep size so as to satisfy the CFL condition. The minimum time step size for explicit time integration depends on the minimum element length and the sonic speed. The critical timestep size for the shell element is given by the equation 't = L / c', where L & c are the characteristic length and the sonic speed respectively .
B] Arbitrary Lagrangian Eulerian Formulation
The use of plain Lagrangian Formulation in large deformation problems leads to severe distortion of mesh rendering it incapable of yielding accurate results. Also the time step size per iteration becomes very small and consequently increases the computational time. To some extent high deformation problems can be modeled by Lagrangian techniques by using 1) Adaptive re-meshing algorithm or 2) Element deletion method. But the adaptive re-mesh algorithm prohibitively increases the computational time and is also not fully developed for three dimensional problems . The deletion of elements results in significant loss of mass. Hence there is a need for an element formulation which would avoid the mesh distortion and also be computationally efficient. The Arbitrary Lagrangian Eulerian formulation also referred to as ALE formulation is used in large strain deformation problems, mainly because it does not face problems regarding mesh distortion. This method is a combination of Lagrangian as well as Eulerian methods.
In ALE formulation, an arbitrary referential domain is defined for the description of motion which is different from the material (lagrangian) and the spatial (eulerian) domain. In the lagrangian method the mesh moves along with the material where as in the Eulerian method the mesh remains fixed and the material moves through it. Also more than one material can be included in an Eulerian formulation. In cases where there are multiple Eulerian materials, Multi Material Groups (known as ALE-MMG) are required to be used. Each of the Eulerian material is assigned to a group using the LS Dyna keyword *ALE_MULTI_MATERIAL_GROUP. It is to be noted that each part within a specific multi material group has identical material properties .
The individual drawbacks of pure Lagrangian and Eulerian formulations are discussed in paragraphs below.
As discussed previously, this method leads to severe mesh deformation since the mesh moves along with the material. In a typical high deformation lagrangian problem, the mesh so severely distorted that the volume of the element is calculated as negative by the software (known as negative volume error in LS Dyna). Eventually the solution terminates and the mesh is required to be manually corrected every time distortion occurs .This may even occur without the material reaching a failure criterion. There is an inherent limit to how much deformation a Lagrangian mesh can accommodate without mesh smoothing or remeshing..
The Eulerian mesh is fixed in space and allows material to pass through it. It is typically used for fluid calculations. . The Eulerian based finite element modeling use Navier-Stokes equations for solution. The solutions are progressed on a fixed mesh and hence avoid mesh distortions. Solving the Navier-Stokes equations are computationally more expensive and complicated as compared to that of Lagrangian formulation. Hence a purely Eulerian approach is not suitable for a fluid structure interaction problem involving explosion.
A hybrid ALE formulation which combines the advantages of both Eulerian as well as Lagrangian techniques is therefore advantageous .
Description of steps in ALE:
The ALE formulation can be defined as algorithms performing automatic rezoning . Manual rezone consists of the following three steps:
Stop Calculations after Mesh Distortion
Remap the solution from Distorted Mesh to Smoothened Mesh.
a] Stop Calculations after Mesh Distortion:
The calculations are stopped after the distortion of mesh so as to perform the mesh smoothing.
b] Mesh Smoothing:
There are many algorithms available for mesh smoothing. Some of them used in LS Dyna are listed below:
Simple Averaging: In this algorithm the coordinates of a node is the simple average of the its surrounding node coordinates.
Kikuchi's Algorithm (Volume Weighted Averaging): This algorithm uses a volume-weighted average of the coordinates of the centroids of the elements surrounding a node.
Equipotential Smoothing: Equipotential zoning is a method of making a structured mesh for finite difference or finite element calculations by using the solutions of Laplace equations as the mesh lines. This same method can also be used to smooth selected points in an unstructured three-dimensional mesh if it is at least locally structured.
Surface Smoothing: The surfaces are smoothed by extending the two-dimensional equipotential stencils to three dimensions.
Combining Smoothing Algorithms: In LS Dyna, the user can use a weighted average of the three algorithms (Equipotential, Simple averaging and Kikuchi's smoothing algorithm) to generate a composite algorithm.
The remap or the advection step maps the solution from a distorted Lagrangian mesh on to the new mesh. During the remap, the topology of the mesh is assumed to be fixed. Also the mesh motion during a step is assumed to be less than the characteristic lengths of the surrounding elements. In other words the Courant number should be less than or equal to one.
â€¦ (eq. 1)
Where, C is the Courant's Number, u is Velocity, t is Timestep and x is the Length Interval. Generally during the advection the following steps are performed:
Decide which nodes to move.
Move the boundary nodes.
Move the interior nodes.
Calculate the transport of the element-centered variables.
Calculate the momentum transport and update the velocity.
Each element solution variable needs to be transported. The total number of solution variables is dependent on the material models. A minimum of six variables need to be transported. The number of variables may go up to fifteen when the elements have strength, where components of the stress tensor and the plastic strain etc. must also be advected. The following schemes can be used for the purpose of advection:
Van Leer advection scheme
Donor Cell advection scheme
The Van Leer advection scheme is a Second order monotonic scheme whereas the Donor Cell is First order accurate scheme. During the transport calculations the errors may bring about the smoothing of the solution, reducing the peak solution variable values. The monotonic second order advection algorithms prevent the transport calculations from creating new spurious oscillations or the peak values (either minimum or maximum) for the variables. Hence Second order accuracy is required for transport of variables. The monotonic algorithms generally provide stability to the code. Also there are many constitutive material models which have a specific allowable range for their history variables. The constitutive models are defined only for the specific range. An example of this is an explosive material model. The burn fraction is required to be within the range zero to one. Other example may be that of Elastoplasticity models which require the plastic strain to be a non-negative number. The second order advection scheme is generally used for supersonic Eulerian flows, where the first order advection schemes may be inaccurate . The Van Leer advection scheme has been used in this problem.
C] Equations of States:
An Equation of State (EOS) is an equation relating the pressure, temperature and specific volume of a substance. Relations involving other properties of substance at equilibrium state are also defined as EOS . In this model we will be using two equations, one for the air and the other for explosive. Air will be modeled using Linear Polynomial eos wheras explosive with Jones-Wilkins-Lee eos.
This eos is used for modeling the air domain. The eos is linear in the internal energy per initial volume (E) is given by the equation:
Where C0, C1, C2, C3, C4, C5 and C6 are user defined constants, V is relative volume and
The Jones-Wilkins-Lee (JWL) equation of state model for explosive detonation products is given by:
In the equation, A, B, R1, R2 and Ï‰ are constants to be calibrated experimentally. V is the product volume relative to the initial explosive volume, E is the energy per unit volume, and P is the pressure. Values for these constants can be developed experimentally or can be found in the literature.
D] Material Modeling
Three types of material models have been used in this problem to model air, explosive and the lagrangian plate. Null, High Explosive Burn and Plastic Kinematic are the models used for modeling air, explosive and plate respectively.
Null material model:
This material model has been used to model the air domain. The Null material model is used especially for modeling fluids and hydrodynamic medium. It is used along with an equation of state which depends on the type of material being modeled. For air, the Linear Polynomial eos is used where as for water Gruneisen eos is used.
Explosive Burn material model
This material model in LS Dyna can be used for modeling the explosive. This material model is used along with the JWL equation of state. This model determines the lighting time for the explosive material.
Burn fractions (F) multiply the equation of states for high explosives to control the the release of chemical energy for simulating Detonations. At any instant, pressure (p) in a high explosive material is given by:
p = Fpeos (V, E)
peos - Pressure from EOS
V - Relative volume
E - Internal energy density per unit initial volume
In the Initialization phase, a lighting time (t1) is computed for each element. It is calculated as follows:
Lighting time for each element (t1) = Distance of Element Center from the detonation point (dn)
Detonation Velocity (D)
If there are multiple detonation points, closest detonation point determines t1.
Burn Fractions, F = max (F1, F2)
2 (t - t1)*D * Aemax / 3Ve if t > t1
0 if t â‰¤ t1
where, Vcj is the Chapman Jouguet relative volume & t is the Current time.
The burn can be initiated either by the Volumetric compression using the function F2. Calculation of burn fraction usually requires several timesteps for F to reach unity, thereby spreading the burn front over several elements. After reaching unity, F is held constant. If F exceeds 1, it is reset to 1. Burn fraction calculation is based on work by Wilkins (1964) and is also discussed by Giroux (1973). As an option, this material can behave as an Elastic Plastic Solid prior to detonation. Once the explosive material detonates, material behaves like a gas.
Shadow burn option :
This option is required while computing of the lighting time if there is presence of elements within the mesh which are not in the direct line of sight from the detonation points. The lighting time is based on the shortest distance through the material. The lighting times also automatically account for variations in the the detonation velocity if different explosives are used.
The lighting time increases if Inert obstacles are present in the explosive material. This is because; the detonation wave requires extra time to propagate around the obstacle. This option can be used for 2 dimensional and 3 dimensional elements.
For best results, following things are recommended :
Uniform mesh containing elements approximately of the same size should be used.
Inert obstacles (such as wave shapers) within the explosive must be larger than the characteristic element size. This is required for the automatic tracking function to function properly. In general,
Size of Inert Obstacle = 2 x Characteristic element dimension
Detonation point should be either within or on the boundary of explosive. Offset points may fail to initialize the explosion.
3] BLAST WAVE/LOAD - STRUCTURE INTERACTION
Simulating structures subjected to blast loads, can be analyzed by
A Purely Lagrangian approach or
ALE and Lagrangian formulations coupled with a Fluid-Structure Interaction (FSI) algorithm.
In a purely lagrangian approach the air blast pressure is computed with empirical equations and directly applied to Lagrangian elements of the structure. Whereas in the FSI approach, the explosive and the air are explicitly modeled and the blast wave propagating through the ALE air domain impinges on the Lagrangian structure through FSI .
The inputs required to be given to the ConWep airblast function include TNT equivalent mass, location of detonation, type of blast and surface identification for which pressure will be applied. Using this information ConWep calculates appropriate pressure to be applied to the lagrangian structure .
A. Purely Lagrangian Approach :
As stated earlier, in the purely lagrangian approach, air blast pressure is computed by using empirical equations (or functions) which directly applied to Lagrangian elements of the structure. The functions available for computing the blast pressure in LS Dyna are Brode function and Conwep (CONventional WEaPons) Function. The ConWep function is more recent and is preferred as compared to Brode's function. ConWep is a collection of conventional weapons effects calculations from the equations and curves of TM 5-855-1, "Design and Analysis of Hardened Structures to Conventional Weapons Effects" . ConWep performs a variety of conventional weapons effects calculations including an assortment of airblast routines, fragment and projectile penetrations, breach, cratering, and ground shock .
The Conwep Function in LS Dyna used in creating the blast, calculates the pressure at various points inside the cylinder with respect to the minimum distance from the point of the blast. The blast loading model was implemented in LS Dyna based on the report by Randers-Pearson and Bannister  to reduce the process of explicitly simulating the progress of the shock wave from the high explosive through the air and its interaction with the structure. Thu the purely Lagrangian approach avoids modeling the air between the explosive and structure. For a larger standoff distance a significant computational cost savings can be done. The shortcoming of the empirical blast equations and functions is their inability to account for confinement (focusing of blast due to geometry) or shadowing (when an object is blocking a surface from direct loading) of the blast waves due to their interaction with structures which may intervene between the explosive and primary structure of interest [14, 17].
The CONWEP function in LS-DYNA can be used in only in two cases:
1) Free air detonation of a spherical charge
2) Surface detonation of a hemispherical charge.
The surface detonation is the preferred option in the case of mines buried 5-20 cm deep. However, in reality the depth of the burial has a significant effect on the energy released from the explosive through the ground and onto the target. Other variables such as moisture content and soil type can also affect how the energy is transferred to the structure. Those below-ground effects are not included in the CONWEP blast model.
B. ALE and Lagrangian formulations coupled with a Fluid-Structure Interaction (FSI) algorithm:
The Fluid Structure Interaction approach has been used in this model. In the FSI approach, the explosives and air (ALE) are explicitly modeled and are coupled to the lagrangian structure by a coupling algorithm. The coupling is done in LS Dyna by using the keyword *CONSTRAINED_LAGRANGE_IN_SOLID. The Penalty coupling factor has been used in this model this Eulerian-Lagrangian Coupling. The Lagrangian structure is considered as the slave surface and the Eulerian fluid as the master. The penalty coupling factor tracks the relative displacement Lagrangian node and the Eulerian fluid material location [5, 18]. Every slave node is checked for penetration through the master surface. Incase there is penetration; an interface force is distributed to the Eulerian fluid nodes.
The magnitude of interface force is calculated as:
F = Ki * d
F = Interface force
Ki = Stiffness factor based on master and slave nodes master properties
d = Relative displacement
The interface force is solved for each time step and is considered as one of the external body forces. Therefore a total nodal force can be solved for each time step so as to derive the structural accelerations, velocities and displacements .
4] PROBLEM SETUP
The explosive considered in this problem is TNT. The explosive charge is spherical, has a radius of 4.1 cm and a charge weight of 453.02 grams. The air domain has also been considered as spherical with a radius of 164 cm. The explosive has been placed at a stand-off distance of 26.14 cm from the structure (plate). The edge length of the square plate is 50 cm whereas the thickness is 2 cm. See fig. 4.1, 4.2 and 4.3.
Fig. 4.1: Air Domain Fig. 4.2: Explosive
Fig. 4.3: FSI sensors placed on the Plate at a Stand-off Distance (26.14 cm) from the Explosive
The Explosive and Air domain are modeled by using solid ALE elements (*Section_Solid_ALE) where as the structural plate by shell elements (*Section_Shell). The parameters for the material models and equation of states for air and explosive are given in table 4.2, 4.3, 4.4 and 4.5. The target plate is made up of steel with the following properties as in table 4.a.1. Plastic-Kinematic material model has been used for the plate. The parameters have been referred from . A consistent unit system has been used for the entire model featuring cm, g and microsecondonds. Mbar is the unit for pressure.
Young's modulus (Mbar)
Yield Strength (Mbar)
Table 4.1: Material properties for Steel Plate
Table 4.2: Parameters for JWL Equation of State (Explosive)
Velocity of Detonation (cm/microsecond)
Chapman-Joguet Pressure (Mbar)
Table 4.3: Parameters for High Explosive Burn material model (Explosive)
Pressure Cutoff (Mbar)
Table 4.4: Parameters for Null material model (Air)
Table 4.5: Parameters for Linear Polynomial Equation of State
Boundary and Initial Conditions:
In this problem eighth symmetry has been considered for modeling the explosive and air whereas quarter symmetry for the structural plate. The air domain and the explosive domain have been included in the ALE multi material group. The multi-material ALE formulation allows the FE mesh to move independently of the material flow and where each element in the mesh can contain a mixture of two, or more, different materials such as air and water .
Unless constrained, the ALE mesh surface remains free under no pressure causing the pressure wave to be reflected from the mesh surface. This can be avoided by extending the mesh surface such that the wave does not reflect until the solution gets terminated. This can unnecessarily increase the computational time and memory usage. To avoid this non-reflecting boundary conditions can be been applied to the outer surface. The option available in LS Dyna can greatly reduce pressure wave reflection. It works by applying a viscous normal/shear pressure to the surface segments . The keyword used is *Boundary_Non_Reflecting. An Initial detonation point is required to be specified for the explosive using *Initial_Detonation keyword. The lighting time also can be specified using this keyword. The point specified is (0,0,0) in x,y,z co-ordinate system.
Database setup for Output results:
In order to measure the exact pressure on the structural plate, FSI sensors have been placed on it. The location of sensors placed is displayed in fig. 4.3. LS Dyna keyword *Database_FSI_Sensor has been used for the sensors. Also tracers have been placed to measure the pressure variation over time. The tracers (*Database_Tracer) either follow the moving mesh or can be fixed in space. Output of history variables for specific nodes has been recorded using *Database_History_Option. Variables such as velocity, acceleration, displacement etc. can be plotted using it.
5] RESULTS AND ANALYSIS
The data from the FSI sensors and history nodes has been displayed in figures below. Fig. 5.1 shows the pressure contour of the medium (air) with 0.3 (Mbar) as the maximum pressure when the blast wave reaches the plate. The series of iso-surface pressure plots at various timesteps is displayed in fig. 5.2.1 - 5.2.5. The fringe levels for pressure are in Mbar. The pressure exerted on the plate as a result of the blast wave is displayed in fig. 5.3.1 - 5.3.5. The time is in microsecondonds for all cases. As discussed earlier, FSI sensors have been placed on the plate at various locations. The pressure plot for various sensors is reflected in the graph (Fig. 5.4).
Fig. 5.1: Pressure Contour (Mbar) at t = 89 microsecondonds
Fig. 5.2.1: t = 10 microsecond
Fig. 5.2.2: t = 20 microsecond
Fig. 5.2.3: t = 40 microsecond
Fig. 5.2.4: t = 60 microsecond
Fig. 5.2.5: t = 110 microsecond
Fig 5.3.1: t=10
Fig 5.3.2: t=20
Fig 5.3.3: t=70
Fig 5.3.4: t=80
Fig 5.3.5: t=100
Fig. 5.4: Pressure output on plate FSI sensors
6] CONCLUSION AND FUTURE STUDY
Air blast analysis and its effect on a structural plate have been studied using ALE formulation in this project. It can be observed that the ALE formulation can prove helpful in predicting the damage to the structure as a result of blast wave interaction. The application of FSI sensors, Tracers and History nodes are some of the methods the processing of output data and studying the obtained results. This study can be further implemented in complex problems involving shadowing effect due to obstructions in the blast wave propagation. Computational Blast analysis has various applications and can extended to study the various blast resistant structures and for experimentation of new blast resistant materials. Also mine blast analysis can be studied using this approach. It can be concluded that the project yields satisfactory results.