RD-V: 0500 Shyue Shock Tube (JWL EOS)

A 1D shock tube filled with detonation gas products is used as a test model for the JWL equation of state.

The shock tube problem is a standard problem in gas dynamics. It can be used as a verification test as an exact solution exists. The Jones-Wilkins-Lee equation of state is used to model the behavior of detonation product of high explosives. Here a 1D shock tube case with detonation gas products governed by the Jones-Wilkins-Lee equation of state is built, and the results compared to the analytical solution. The influence of the second order integration scheme in time and space (default in Radioss) is shown.

The following criteria is used to compare the results:

Pressure
Mass density
Velocity
Specific energy

Options and Keywords Used

Explosive material, /MAT/LAW5 (JWL)
Multi-fluid law using the FVM solver, /MAT/LAW151 (MULTIFLUID)
EULER property formulation with I_ale = 2
Second order integration scheme, /ALE/MUSCL

Input Files

Before you begin, copy the file(s) used in this problem to your working directory.

RD-V-0500_Shyue_Shock_Tube.zip

Model Description

A shock tube consists of a long tube filled with gas, initially divided into two section separated by a diaphragm.

The gas in the two sections are in different physical states: there is usually a high-pressure and low-pressure section.

At the beginning of the experiment, the separation is removed, and a compression shock runs into the low-pressure region, while a rarefaction wave moves into the high-pressure part of the tube. A contact discontinuity usually occurs between the two gases when the diaphragm is removed.

rad_ex_fig_13-3 — Figure 2. Schematic Shock Tube Problem with Pressure Distribution for Pre- and Post-diaphragm Removal

The interest of the shock tube problem is that it involves the three possible fundamental waves that can develop from uniform initial conditions: shock, rarefaction and contact waves.

Considered in one dimension, the shock tube is a Cartesian geometry Riemann problem, for which an analytical solution exists. When applied to a gas governed by a JWL equation of state, the shock tube case in 1D is called the test problem of Shyue. An analytical solution of this problem is described with analytical solutions available. ^¹

A shock tube 100 cm long is filled with detonation gas products and is divided of two chambers of equal length. The gas pressure in the left section (the high-pressure section) is 10 times the pressure in the right section (low-pressure section).

The physical properties of the detonation gas products in the right and left sections are:

Table 1. Initial Conditions
	High Pressure Section	Low Pressure Section
Pressure P	10.0 Mbar	1 Mbar
Mass density $ρ$	1.7 $\frac{g}{c m^{3}}$	1.0 $\frac{g}{c m^{3}}$
Velocity $u$	0.0 $\frac{c m}{μ s}$	0.0 $\frac{c m}{μ s}$

The following system is used: cm, ms, g, daN, Mbar

Model Method

Each section is modeled by one component.

Both are meshed with 3D /BRICK elements. As the problem is one-dimensional, the relevant mesh size is the width of the mesh in the direction of interest (here, the direction of propagation of the compression and rarefaction waves). The size of the mesh in the other directions is of no importance (Figure 3).

Figure 3. Relevant Mesh Size for 1D Model

The model is meshed with 512 elements (which accounts for a mesh size in the direction of interest of roughly 2 mm).

Figure 4. Mesh Used for the 1D Shock Tube

The detonation gas products in the two sections are modeled using /MAT/JWL which is based on the JWL equation of state. The JWL materials are then referenced by /MAT/MULTIFLUID which allows the use of the Radioss finite volume solver.

The two material laws (/MAT/JWL) that represent the high- and low-pressure gas are used in a multi-fluid law (/MAT/MULTILFUID) in order to use the Finite Volume solver.

Detonation Gas Products Modeling

This shock tube case considers detonation gas products. In Radioss, explosive materials defined with /MAT/JWL are considered solid until they detonate. To get these materials in their gaseous states, a detonator is created with a /DFS/DETPOINT card. The detonation time variable TDET is set to -1E30s (meaning minus infinity), ensuring that at the beginning at the simulation all the explosive material in the model consists of detonation gas products.

Detonation gas product behavior is governed by the JWL equation of state. This EOS relates Pressure to the specific volume of the explosive.

The JWL EOS expression is:

P_{j w l} (υ, E) = A (1 - \frac{ω}{R_{1} υ}) e^{- R_{1} υ} + B (1 - \frac{ω}{R_{2} υ}) e^{- R_{2} υ} + \frac{ω E}{υ}

Where,

$P_{j w l}$: Pressure
$υ$: Relative volume $υ = \frac{ρ_{0}}{ρ}$
$A$ , $B$ , $R_{1}$ , $R_{2}$: Parameters specific to the explosives
$ω$: Gruneisen coefficient
$E$: Detonation energy per unit volume

The /MAT/JWL cards used for the detonation gas products are defined using the following parameters:

Table 2. Values of Input Parameters for the /MAT/JWL Cards
		High Pressure Section	Low Pressure Section
Density	Initial density $ρ$	1.7 $\frac{g}{c m^{3}}$	1.0 $\frac{g}{c m^{3}}$
Density	Reference density $ρ_{0}$	1.84 $\frac{g}{c m^{3}}$	1.84 $\frac{g}{c m^{3}}$
JWL parameters	A	8.545 Mbar
	B	0.205 Mbar
	R1	4.6
	R2	1.35
	$ω$	0.25
Chapman-Jouguet parameters	Detonation velocity D_CJ	0.693 $\frac{c m}{μ s}$
	Detonation pressure P_CJ	0.21 Mbar
	Detonation Energy E₀	42.88164777 Mbar	7.233944805 Mbar

The initial density corresponds to the density of the explosive in each compartment, while the reference density is the one used to compute the relative volume used in the JWL equation of state, as:

v = \frac{ρ}{ρ_{0}}

The JWL parameters comes from the reference article. As this model is a test case for computer code verification, they are not related to any existing explosive materials.

The Chapman-Jouguet parameters used for modeling the propagation of the detonation are of no importance as this model only considers gaseous materials. The Chapman-Jouguet parameters of the TNT have been used for detonation pressure P_CJ and velocity D_CJ. The precise value of these two parameters has no influence on the results, as they are only used to model the detonation process of the solid explosive. The detonation energy E₀ is computed using the JWL equation of state expression to get an initial pressure $P_{l e f t} = 10 M b a r$ for the high-pressure gas, and $P_{r i g h t} = 1 M b a r$ for the low-pressure gas.

Two instances of the multi-material law (/MAT/MULTILFUID) are created, each one corresponding to a section and integrating the gas defined using LAW5, with a volumetric fraction of 1.

/EULER/MAT should be defined for the /MAT/MULTIFLUID materials, to indicate that they are EULERIAN materials.

The /ALE/MUSCL option activates a full second order integration scheme in time and space, which increases the accuracy of the results. The results of the simulations with and without /ALE/MUSCL will be compared to the analytical solution. ^¹

Boundary Conditions

Boundary conditions around the 1D shock tube are sliding walls, which is the default when using (/MAT/MULTIFLUID).

Engine Control

Since the Radioss 2019.0 release, the time step scale factor for all ALE and EULER elements is set by default to 0.5.

It can be modified using the keyword /DT/ALE:

/DT/ALE
0.5                  0.0

Results

Results are examined 12 $μ$ s after the removal of the separation between the two chambers.

Figure 6. Pressure Contour after 12 $μ$ s

Figure 7. Density contour and density path after 12 $μ$ s ,as well as energy balance

Spatial profile of pressure, mass density, velocity and energy are compared to the analytical results given in the reference article. These profiles are obtained by drawing a node path along the tube. To draw node paths, the elemental results visible on the animation have to be averaged at the nodes, using for example the averaging method simple.

Profiles of pressure and mass density along the tube after 12

μ

s are visible on Figure 8 and Figure 9, as well as similar profiles for velocity and specific energy on Figure 10 and Figure 11, with the curve obtained from the analytical solution.

Figure 8. Pressure Profile. (with and without /ALE/MUSCL, compared to analytical solution)

Figure 9. Mass Density Profile. (with and without /ALE/MUSCL, compared to analytical solution)

Figure 10. Velocity Profile. (with and without /ALE/MUSCL, compared to analytical solution)

Figure 11. Specific Energy Profile. (with and without /ALE/MUSCL, compared to analytical solution)

The simulation results for these four variables matches closely the analytical solution, and the fit is better for the results obtained with /MUSCL enabled (MUSCL can be disabled with /ALE/MUSCL/OFF Engine card).

The difference between numerical results and the analytical solution can be quantified by computing the £₂-norm error between the curves.

The £₂-norm error (also named Eulerian norm error) between a numerical result $Y_{n u m}$ and the exact solution $Y_{e x a c t}$ is defined as:

L_{2} e r r o r = \sqrt{\frac{\int {(Y_{n u m} - Y_{e x a c t})}^{2} d x}{\int Y_{e x a c t}^{2} d x}}

The £₂-norm error to the analytical solution on pressure, mass density, velocity are:

Table 3. £₂-norm error to the analytical solution with and without /ALE/MUSCL
	/ALE/MUSCL	Pressure	Mass Density	Velocity	Specific Energy
£₂-norm relative error	enabled	2.0%	6.5%	3.9%	4.5%
£₂-norm relative error	disabled	2.5%	9.2%	4.7%	6.7%

As shown on Table 3, the use of /ALE/MUSCL allows for better precision in the numerical results.

Conclusion

The numerical results show a good correlation with the analytical solution presented in the reference article. /MAT/JWL used within a /MAT/MULTIFLUID multi-material law allows for an accurate resolution use of the JWL equation of state.

There are clear benefits of using MUSCL scheme (default in Radioss) in terms of precision.

¹ Kamm, J.R. An Exact, Compressible One-Dimensional Riemann Solver for General, Convex Equation of State. Los Alamos National Laboratory, 2015