Appropriate Coupling Solvers for the Numerical Simulation of Rolled Homogeneous Armor Plate Response Subjected to Blast Loading

Abstract

Rolled homogeneous armor (RHA) plate subjected to blast loading is a complex problem involving the nonlinear fluid-structure interaction. The numerical techniques using the spatial discretization scheme that has been provided as a solver in the AUTODYN computer code will be used in this study in order to predict the RHA response subjected to explosive (TNT) blast loading. The final deflection will be used as a reference in order to identify the suitable solver for both materials RHA and TNT; then the plastic deformation will be chosen in the simulation process. Instead of using the same solver for RHA and TNT domains, the optimization of solver can be achieved if it is only used in an appropriate domain, or in other words, a different domain will be using different solver. The solvers, which were available in AUTODYN, were used in the analysis of impact and explosion or fluid-structure interaction. Therefore, in this paper, we will determine the suitable solver for both materials (TNT and RHA plate), and the appropriate interaction coupling solver will be obtained. Defining TNT and RHA plates using the Arbitrary Lagrangian Eulerian solver has found the best coupling solver for this case study when compared with existing experimental data. This coupling solver will be used for future analysis in simulating blast-loading phenomena.

1. Introduction

In the blast phenomena, interaction between fluid and structure, also called fluid-structure interaction (FSI), normally will occur, and there is no single method that can be used for all conditions in FSI analysis [1]. The governing partial deferential equation for FSI model needs to be solved in both time and space domain with the basic physic principles involving the conservation of mass, momentum, and energy. The solution over the time domain can be achieved by an explicit method [2]. It can be obtained by utilizing different spatial discretization such as Lagrange, Euler, and Arbitrary Lagrangian Eulerian (ALE) or mesh-free method also known as Smooth Particle Hydrodynamic (SPH) methods [1]. However, the basic solvers for explicit integration numerical wave-codes (sometimes termed “hydro-codes”) can be utilized as an outline with their associated strengths and weaknesses [3]. Air, plate, and trinitrotoluene (TNT) are three different domains in this model analysis. Each domain has a solver that is suitable to be used. The numerical solver be used in AUTODYN generally fall into the following methods which are Lagrange, Euler, ALE, and SPH methods. With intelligent selection of suitable solver for various regimes, an optimal solution in terms of accuracy and efficiency can be achieved. The appropriate solvers in AUTODYN will be carried out, and the effect on a final result will be discussed further.

2. Method of Discretization

Numerical analysis needs to be subdivided into problem domain analysis to the nodes, elements, and spatial discretization. The spatial discretization is performed by representing the fields and structures of the problem using computational point in space, and connected with each other through a computational grid. Usually, the fine grid will lead to more accuracy of the result. The most popular spatial discretization has been widely used are Langrage, Euler, ALE and SPH and are provided in the AUTODYN computer code.

2.1. The Lagrange Solver

Lagrange's equation has been use to formulate the equation of motion, which represents the response of structure subjected to any external load. Simple Lagrange's governing equation can be derived by considering conservative system, where all internal forces have a potential energy and Lagrangian discretization scheme is usually a based on the finite element method.

In the Lagrangian scheme, the elements (and hence the nodes) move with the material, while in the Eulerian scheme the nodes are fixed in space, but the material is allowed to flow as described in [3]. Thus, the characteristics make the Lagrangian scheme more suitable for modeling solid materials, while the Eulerian scheme is more suitable for modeling fluid materials. The advantages and limitations of the Lagrangian solver are as summarized in Table 1.

Table 1:

Lagrangian solver summary [3].

Advantages	Limitations	Enhancement to overcome limitation

Efficiency, fewer computations per time step relative to other solvers	Element distortion can lead to small-time step	Rezoning and erosion

Clear definition of material interface and boundaries	Cell distortion can lead to grid tangling inaccuracies	Rezoning and erosion

Good time history Information	Thin sectionneeds small-time step	Use of structural solver

Good for strength modeling	Complex logic for sliding interface

Simpler code

2.2. The Euler Solver

The Euler method of space discretization as described in [3], where the numerical meshes or grids are fixed in space and the physical material flows through the mesh. It is typically well suited for the description of the material behavior of large deformation or flows situation and consequently, by definition, it does not result in grid distortion due to the fixed grid scenario.

Most of the time, the Euler is usually used for representing fluids and gases. However, to describe for solid behavior, additional calculations are required to transport the solid stress tensor and the histories of the material through the grid. The advantages and limitations of the Eulerian solver are as summarized in Table 2.

Table 2:

Eulerian solver summary [3].

Advantages	Limitations	Enhancement to overcome limitation

No grid distortion	More computations per cycle than Lagrange	Simplified Euler formulation without strength

Large deformation handle	Diffusion of material boundaries	Sophisticated interface tracking implementations

Allows mixing of different materials within cells	Need finer zoning for similar accuracy compared to Lagrange Need interactive equation of states solver for multiple material cells	Use of higher order techniques

Rezoning not required	Need void cell defined in space where material may flow to very large meshes	Use of dynamic memory allocation to store information only where needed

	Shock diffusion	Higher order implementations

	Less flexible for strength modeling than Lagrange solvers

	Thin section needs many cells	Couple to alternative solvers

2.3. The Arbitrary Lagrange-Euler (ALE) Solver

The ALE is a special solver that incorporates Lagrange and Euler solvers in a single governing equation and provides a full coupling between the blast wave and the structural response. The interaction of solid-type material in Lagrangian and fluid-type material in Eulerian generated a FSI.

In the differential form of the conservative equation of mass, momentum and energy are more readily obtained from the corresponding well-known Eulerian form than from the ALE form; the conservative equation is to be replaced in the various conservative term, such as material velocity, mass density and stress tensor of material.

The ALE solver, also allows for “automatic rezoning,” which can be quite useful for certain problems. With these features, this solver is perhaps appropriate for modeling solid, fluid, and gas. Thus, generally, this solver is suitable for a variety of fluid-structure interaction problem analyses. Baylot and Bevins [4] had performed the analysis of blast wave propagation using the ALE solver and the structure response in the Lagrangian scheme simultaneously. However, like any other discretization schemes, the ALE solver also has its advantages and limitations as summarized in Table 3.

Table 3:

ALE solver summary [3].

Advantages	Limitations	Enhancement to overcome limitation

Wide range of applicability with arbitrary mesh motions (automatic rezoning)	Need to specify mesh motion constraints

Clear definition of material interfaces and boundaries	Cell distortion can lead to grid tangling and inaccuracies	Rezoning and erosion

Good time history information	Thin sections need small-time step	Use shell solver

Good for strength modeling	Complex logic for sliding interface

2.4. The Smooth Particle Hydrodynamic (SPH) Solver

Based on the Lagrangian space discretization scheme, SPH was developed as an alternative method without any nodes or grid approaches or also known as a meshless method. The SPH solver was initially used in astrophysics [5]. It had the potential to be efficient and accurate in material deformation as well as flexible in terms of the inclusion of specific material models.

The idea behind SPH is to replace the equations of fluid dynamics by equations for particles. The SPH particle is not only interacting mass point but also interpolation point used to calculate the value of physics variables based on the data from neighboring SPH particle, scaled by the weighting function. Since there is no grid defined, the SPH method does not suffer from the normal problem of grid tangling in large deformation problem.

However, according to Quanet al. [1] the SPH method requires a sort of the particle in order to locate the current neighboring particles, which makes the computational times per cycle more expensive than the mesh based on Lagrangian technique. This can make the SPH method less efficient than the Lagrange's method.

3. Validation of the Few Solvers

The appropriate solvers discussed in Section 2 were validated against one of the series of blast test experiments performed by Neuberger et al. [6]. The experiment consisted of a circular target plate with a diameter of 1000 mm that was clamped with two thick armor steel rings, tightened together with bolts, fully clamped, and subjected to 8.75 kg of TNT located at 200 mm standoff distance. The experimental setup is as shown in Figure 1.

Figure 1:

The experimental setup [6].

The experimental setup as shown in Figure 1 will be modeled in AUTODYN 3D and the equivalent simulation model was developed as shown in Figure 2.

Figure 2:

A three-dimensional simulation model created in AUTODYN 3D.

The simulation was conducted by using a high speed computer processor Intel Core i7-2600K CPU@ 3.40 GHz (8 CPUs), 3.7 GHz. The model was created in a domain size of (500 × 500 × 500) mm with 5 mm element size defined as air and symmetrical in the “X” and “Z” axes. A quarter symmetry model was used in the simulation in order to reduce the computational time.

The same mechanical property of the RHA material in [7] was used in the simulation model as shown in Table 4 where t is plate thickness, σ_y is yield stress, σ_UTS is ultimate tensile stress, ε_L is strain, E is Modulus Young, ν is Poisson's ratio, andA, B, n, c, and m are material dependent parameters, and may be determined from empirical fit of flow stress data.

Table 4:

Material properties for RHA [7].

Properties	Value

t	3–20 mm
σ_y	950 MPa
σ_UTS	1250 MPa
ε_L	9%
E	210 GPa
ν	0.28
Hardness	380–430 HBW
A	950 MPa
B	560 MPa
n	0.26
c	0.0014
m	1

There are four solvers that will be used to calculate interaction between TNT explosive and the RHA plate with respect to the EoS for RHA and TNT material. Johnson and Cook EoS has been selected to represent RHA material as following:

Y = A (1 + \frac{B}{A} ε^{n}) (1 + C \ln {\dot{ε}}^{*}) (1 - T^{* m}),

(1)

where A, B, C, m, and n are constant, ε is equivalent plastic strain, ${\dot{ε}}^{*}$ is the strain rate nondimensionalized by reference strain rate of 1/s, and T^* is nondimensional temperature. Parameter A, the initial (ε ≈ 0) yield strength of the material at a plastic strain rate of $\dot{ε} = 1 / s$ at the room temperature (298 K), is modified by a strain-hardening factor (containing parameter, B and n), a strain-rate-hardening factor (containing parameter C), and a thermal-softening factor (containing parameter m), while T^*, defined by

T^{*} = \frac{T - T_{r}}{T_{m} - T_{r}},

(2)

where T_r is room temperature and T_m is the melting temperature of material, 1783 K for RHA. Equation (2) is the form used in USArmy Research Laboratory (ARL) [8] and is valid for T_r ≤ T ≤ T_m, the region of interest in most blast and ballistics application.

For the blast application, this paper used Jones Wilkins Lee (JWL) EoS to represent TNT blast loading. The JWL EoS has been used by AUTODYN to accurately describe the pressure-volume-energy behavior of the detonation product as follows:

P = A (1 - \frac{ω}{R_{1} V}) e^{- R_{1} V} + B (1 - \frac{ω}{R_{2} V}) e^{- R_{2} V} + \frac{ω E}{V},

(3)

where V is the volume of detonation product or, P is pressure in megabar, and E is the energy in Mbcc/cc, while A, B, R₁, R₂ and ω are constant parameters based on the type of explosive. For TNT, it was identified by Lee et al. [9], that A = 5.242, B = 0.07678, C = 0.01082, R₁ = 4.2, R₂ = 1.10 and ω = 0.30.

4. Results and Discussion

4.1. Interaction of Coupling Solver

TNT as an explosive material will produce a blast force that propagated through the air as a medium and hit the RHA plate as the target object. This phenomenon involved the interaction between fluid and structure.

The analysis started when TNT as high-explosive material detonated, it interacted with the surrounding which were defined as an air ideal gas. Air as a medium containing the explosive charges striked the RHA plate. The TNT and plate domain interacted. In the numerical simulation, interaction involved two or more type, domains using the same or different type of solvers. Three material variable, that is; air, TNT, and RHA plate, were involved in the case response of RHA plate subjected to blast loading. For air, the medium has been set according to Chung Kim Yuen et al. [10] had established the Ideal Gas Euler solver which was the most appropriate solver for air (Table 5). In this paper, the appropriate solver for TNT and RHA plate solver will be determined. There are four solvers for TNT that had been considered as follows Lagrange, ALE, Euler, and SPH while the two solvers for RHA plate; Lagrange and ALE. The interaction diagram between TNT and RHA plate is as illustrated in Figure 3. The appropriate solver interaction will be used in future analysis.

Table 5:

Properties of air [10].

Properties	Value

ρ	1.225 kg/m³
T	288.2 K
C _p	1.005 kJ/kgK
C _v	0.718 kJ/kgK
e	2.068 × 10⁵ kJ/kg

Figure 3:

The diagram interaction between solvers to represent TNT and RHA plates.

4.2. The Euler-Lagrange Solvers Interaction

In this type of interaction, it was involved in the interaction between RHA plates represented by the Lagrange solver, and TNT represented by Euler solver. This interaction was also called Euler-Lagrange coupling, and was commonly used to simulate an interaction between fluid and structure domains.

Simultaneous analysis in the numerical simulation, for both domains; Lagragian grids provide the geometry constraint in order to allow the material to flow in the Eulerian grids. At the Euler-Lagrange interface, the Lagrange grids act as structure geometry inside the flow boundary, while the Euler grids will provide a pressure or heat boundary to the Lagrangian domain. As the Lagrange's grids move or distort, the effect of interaction will produce the deflections of a structure as illustrated in Figure 4.

Figure 4:

The Euler-Lagrange interaction process.

In the finite element analysis (FEA), the Euler solver is commonly used for fluid analysis by using the finite volume method (FVM) or the finite difference method (FDM) discretization scheme, and the solution involves large relative deformation, whereas the Lagrange used the finite element method (FEM) discretization scheme for structure analysis for a small relative deformation. Both solvers provide different rates of deformation and will lead to the small-time step. This problem frequently happens in the numerical simulation. Coupling solver for this case showed the plate deflection stopped at 0.15 ms due to the small-time step as shown in Figure 5.

Figure 5:

RHA deflection captured using Euler Lagrange interaction before computation was stoped due to the small time step.

However, according to Birnbaum et al. [3], for two dimensional (2D) cases, Euler-Lagrange coupling showed very powerful and stable coupling for FSI problem, whereas in three-dimensional (3D) cases, the computational requirement time was excessive even on supercomputers. This stems from the substantial and complex intersection calculation that must be performed for each time step. In this paper, the coupling Euler-Lagrange solvers were unable to produce a good results when compared with the experimental data. Thus, other solver interactions should be considered.

4.3. The Smooth Particle Hydrodynamic (SPH) Interaction

Previous numerical simulations used SPH solver to see the trajectory particle on the impact analysis and commonly was used for the porous material such as soil, sand and water. Toussaint and Durocher [11], Barsotti et al., [12] and Quan et al., [1] used SPH solver to define for sand and Campbell and Vignjevic [13] for water.

In this paper, numerical simulation will be involved in the interaction between RHA plate and TNT with respect to the SPH, Lagrange, and ALE solvers. The case is as follows;

RHA plate with Lagrange solver and TNT with SPH solver;

RHA plate with ALE solver and TNT with SPH solver.

The limitations of Euler-Lagrange coupling in AUTODYN 3D simulation lead to explore on coupling for another solver interaction to obtain the promising result. It should also evaluate the coupling between the mesh-less technique (SPH) representative of TNT and the traditional Lagrange and ALE solvers representative of RHA plate. By using the SPH method, TNT was discretized as the continuum through a set of the nodes without the connective mesh. The node or element assumed a physical meaning; they represent material particle carrying properties such as pressure and thermal and impacted the RHA plate as shown in Figure 6 representing case (a) and Figure 8 representing case (b). The deformation results of the RHA plate for both cases (a) and (b) are as shown in Figures 7 and 9. The approximate final deformation results had large discrepancies when compared with experimental data. Thus, the SPH-Lagrange coupling was not the appropriate solution for this particular case.

Figure 6:

The RHA plate (Lagrange)-TNT (SPH) interaction process.

Figure 7:

Case (a) RHA plate response to Lagrange solver and TNT to SPH solver.

Figure 8:

The RHA plate (ALE)-TNT (SPH) interaction process.

Figure 9:

Case (b) RHA plate response to ALE solver and TNT to SPH solver.

4.4. The Lagrange and ALE Interaction (RHA Plate and TNT Represented by ALE/Lagrange Solver)

This simulation will involve the four interaction cases between RHA plate and TNT with respect to the Lagrange and ALE solvers. The cases are as follows:

RHA plate with ALE solver and TNT with ALE solver;

RHA plate with Lagrange solver and TNT with Lagrange solver;

RHA plate with Lagrange solver and TNT with ALE solver;

RHA plate with ALE solver and TNT with Lagrange solver (Figure 17).

A surface on Lagrange's domain interacted with another surface of a different domain allowing for impact and sliding, gap opening and closing, and mesh distortion between the bodies. The interaction between TNT and RHA plate using Lagrange and ALE element for all cases found was a similar pattern on the grid distortion on the RHA plate as shown in Figures 10, 12, 14, and 16. Consequently, the predictions on the final deflection of the RHA plate were in good agreement with all cases on the Lagrange and ALE interaction as shown in Figures 11, 13, and 15 despite differences in the computation time taken and cycle interaction as shown in Table 6.

Table 6:

Result consolidation.

Neuberger et al. [6] experimental data Residual deflection = 64 mm

TNT solver	Plate solver
	Lagrange	ALE

Lagrange	Residual deflection = 68 mm	Residual deflection = 60 mm
	Computation time = 66 hours	Computation time = 108 hours
	Cycle of iteration = 400000	Cycle of iteration = 500000

ALE	Residual deflection = 60 mm	Residual deflection = 63 mm
	Computation time = 102 hours	Computation time = 96 hours
	Cycle of iteration = 530000	Cycle of iteration = 430000

SPH	Residual deflection = 25 mm	Residual deflection = 32.64 mm
	Computation time = 80 hours	Computation time = 40 hours
	Cycle of iteration = 326000	Cycle of iteration = 327000

Euler	Residual deflection = 12 mm	Residual deflection = 12 mm
	Computation time = 4 hours	Computation time = 4 hours
	Cycle of iteration = 1200	Cycle of iteration = 1200
	*Iteration top due to time step to small	*Iteration top due to time step to small

Figure 10:

The RHA plate (ALE)-TNT (ALE) interaction process.

Figure 11:

Case (a) RHA plate response to ALE solver and TNT to ALE solver.

Figure 12:

The RHA plate (Lagrange)-TNT (Lagrange) interaction process.

Figure 13:

Case (b) RHA plate response to Lagrange solver and TNT to Lagrange solver.

Figure 14:

The RHA plate (Lagrange)-TNT (ALE) interaction process.

Figure 15:

Case (c) RHA plate response to Lagrange solver and TNT to ALE solver.

Figure 16:

The RHA plate (ALE)-TNT (Lagrange) interaction process.

Figure 17:

Case (d) RHA plate response with ALE solver and TNT with Lagrange solver.

However, taking into account the computation time taken, cycles of iterations had been obtained and consolidated in Table 6 and deformation results are as shown in Figure 11, interaction using ALE solver for both material TNT and RHA plate showed an appropriate interaction selection in this case. The deformation result shown in Figure 11 was in very good agreement with when compared with Neuberger et al. [6] experimental data.

5. Conclusions

Different numerical solvers have certain advantages and limitations. It is critical that analyst must understand which solver or combination of solver is appropriate to be used in a particular case or problem of interest.

Use of computer code or software package intelligence is considered as a key in order to obtain the good prediction result. No doubt, the coupling techniques are extremely powerful for any FSI analyses, although they have inherent limitations that must be recognized. The solver investigation must be done before the appropriate solver is obtained for a particular case. However, the suitable coupling solver can only be selected after the analyst performed numerical trial and validation against experiment.

This paper presents the appropriate solver coupling of a RHA plate subjected to blast loading and has been compared with Neuberger et al. [6] experimental data. For this case, it was found out that using the ALE solver to represent both the TNT and RHA plate in the numerical simulation, managed to produce good agreement with the experimental test data. Thus, the ALE solver will be chosen as the coupling solver to simulate similar cases for future analysis.

Footnotes

Acknowledgments

Financial assistance from the Royal Malaysian Air Force (RMAF) towards this research is hereby acknowledged. Opinions and views either directly or indirectly from various parties, including the staff and lecturers from the Engineering Research Center of National Defense University (UPNM) are greatly appreciated.

References

Quan

Birnbaum

N. K.

Cowler

M. S.

Gerber

B. I.

Clegg

R. A.

, and Hayhurst

C. J.

, “Numerical simulation of structural deformation under shock and impact loads using a coupled multi-solver,” in Proceedings of the 5th Asia-Pacific Conference on Shock and Impact Load on Structures, Hunan, China, November 2003.

Tannehill

J. C.

Anderson

D. A.

, and Pletcher

R. H.

, Computational Fluid Mechanic and Heat Transfer, Taylor & Francis, 1798.

Birnbaum

N. K.

Nigel Francis

, and Bence Gerber

, “Coupled techniques for the simulation of fluid-structure and impact problem,” Century Dynamic, 2003, http://hsrlab.gatech.edu/AUTODYN/papers/paper63.pdf.

Baylot

J. T.

and Bevins

T. L.

, “Effect of responding and failing structural components on the airblast pressures and loads on and inside of the structure,” Computers and Structures, vol. 85, no. 11–14, pp. 891–910, 2007.

Gigold

R. A.

and Monaghan

J. J.

, “Smoothed particle hyrodynamics: theory and application to non spherical star,” Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977.

Neuberger

Peles

, and Rittel

, “Springback of circular clamped armor steel plates subjected to spherical air-blast loading,” International Journal of Impact Engineering, vol. 36, no. 1, pp. 53–60, 2009.

Neuberger

Peles

, and Rittel

, “Scaling the response of circular plates subjected to large and close-range spherical explosion—part I: air-blast loading,” International Journal of Impact Engineering, vol. 34, pp. 859–873, 2007.

Meyer

H. W.

Jr. and Kleponis

D. S.

, An Analysis of Parameter for the Johnson-Cook Strength Model for 2-in-Thick Rolled Homogenous Armor, Weapon and Materials Research Directorate Army Research Laboratory, 2001.

Lee

Finger

, and Collins

, JWL Equation of State Coefficients for High Explosives, Lawrence Livermore Laboratory Calofornia, 1973.

10.

Chung Kim Yuen

Langdon

G. S.

Nurick

G. N.

Pickering

E. G.

, and Balden

V. H.

, “Response of V-shape plates to localised blast load: experiments and numerical simulation,” International Journal of Impact Engineering, vol. 46, pp. 97–109, 2012.

11.

Toussaint

and Durocher

, “Finite element simulation using SPH particles as loading on typical light armoured vehicles,” in Proceedings of the 10th International LS-DYNA User Conference, East Lansing, Mich, USA, June 2008.

12.

Barsotti

M. A.

Puryear

J. M. H.

Stevens

D. J.

Alberson

R. M.

, and McMahon

, “Modeling mine blast with SPH,” in Proceedings of the 12th International LS-DYNA User Conference, East Lansing, Mich, USA, June 2012.

13.

Campbell

J. C.

and Vignjevic

, “Simulation structural response to water impact,” International Journal of Impact Engineering, vol. 49, pp. 1–10, 2012.