RD-V: 0505 Shock Tube
The transitory response of a perfect gas in a long tube separated into two parts using a diaphragm is studied. The problem is well-known as the Riemann problem. The numerical results based on the finite element method with the Lagrangian and Eulerian formulations, are compared to the analytical solution.
- Lagrangian (mesh points coincident to material points)
- Eulerian (mesh points fixed)
The propagation of the gas in the tube can be studied in an analytical manner. The gas is separated into different parts characterizing the expansion wave, the shock front and the contact surface. The simulation results are compared with the analytical solution for velocity, density and pressure.
Options and Keywords Used
- Brick elements
- Lagrangian and Eulerian formulations
- Hydrodynamic viscous fluid law (/MAT/LAW6 (HYDRO or HYD_VISC)) and perfect gas modeling
- ALE boundary conditions (/ALE/BCS)
- ALE material formulation (/ALE/MAT)
Input Files
The input files used in this verification problem include:
<install_directory>/hwsolvers/demos/radioss/verification/Blast/0505_shock_tube/*
Model Description
The shock tube problem is one of the standard problems in gas dynamics. It is a very interesting test since the exact solution is known and can be compared with the simulation results. The Finite Element method using the Eulerian and Lagrangian formulations was used in the numerical models.
The initial state at time t = 0 consists of two constant states 1 and 4 with $${p}_{4}>{p}_{1}$$, $${\rho}_{4}>{\rho}_{1}$$, and ${v}_{4}={v}_{1}=0$ (table).
High Pressure Side (4) | Low Pressure Side (1) | |
---|---|---|
Pressure $$p$$ | 500000 $$\left[\text{Pa}\right]$$ | 20000 $$\left[\text{Pa}\right]$$ |
Velocity $v$ | 0 $$\left[\frac{\text{m}}{\text{s}}\right]$$ | 0 $$\left[\frac{\text{m}}{\text{s}}\right]$$ |
Density $\rho $ | 5.7487 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ | 0.22995 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ |
Temperature $$T$$ | 303 $\left[\text{K}\right]$ | 303 $\left[\text{K}\right]$ |
Just after the membrane is removed, a compression shock runs into the low pressure region, while a rarefaction (decompression) wave moves into the high pressure part of the tube. Furthermore, a contact discontinuity usually occurs.
Model Method
The hydrodynamic viscous fluid LAW6 is used to describe compressed gas.
With $$\mu =\frac{\rho}{{\rho}_{0}}-1$$
- $$p$$
- Pressure
- $${C}_{i}$$
- Hydrodynamic constants
- $${E}_{n}$$
- Internal energy per initial volume
- $\rho $
- Density
$${C}_{4}={C}_{5}=\gamma -1$$
Where, $\gamma $ is the gas constant.
Under the assumption $\gamma =const.=1.4$ (valid for low temperature range), the hydrodynamic constants ${C}_{4}={C}_{5}=0.4$.
Parameters of material LAW6 are provided in Table 2.
High Pressure Side (4) | Low Pressure Side (1) | |
---|---|---|
Initial volumetric energy density (E_{0}) | 1.25x10^{6}$\left[\frac{\mathrm{J}}{m{m}^{3}}\right]$ | 5x10^{4} $\left[\frac{\mathrm{J}}{m{m}^{3}}\right]$ |
$${C}_{4}$$ and $${C}_{5}$$ | 0.4 | 0.4 |
Density $\rho $ | 5.7487 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ | 0.22995$\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ |
Analytical Approach
Evolution of the flow pattern is illustrated in Figure 4. When the diaphragm bursts, discontinuity between the two initial states breaks into leftward and rightward moving waves, separated by a contact surface.
- Zone 1 is the low pressure gas which is not disturbed by the shock wave.
- Zone 2 (divided in 2 and 2' by the contact surface) contains the gas immediately behind the shock traveling at a constant speed. The contact surface across which the density and the temperature are discontinuous lies within this zone.
- The zone between the head and the tail of the expansion fan is noted as Zone 3. In this zone, the flow properties gradually change since the expansion process is isentropic.
- Zone 4 denotes the undisturbed high pressure gas.
Equations in Zone 2 are obtained using the normal shock relations. Pressure and the velocity are constant in Zones 2 and 2'.
The ratio of the specific heat constant of gas $\gamma $ is fixed at 1.4. It is assumed that the value does not change under the temperature effect, which is valid for the low temperature range.
Zone 4 | Zone 1 | |
---|---|---|
Pressure $$p$$ | $${p}_{4}$$ = 500000 $$\left[\text{Pa}\right]$$ | $${p}_{1}$$ = 20000 $$\left[\text{Pa}\right]$$ |
Velocity $v$ | ${v}_{4}=0$$$\left[\frac{\text{m}}{\text{s}}\right]$$ | ${v}_{1}=0$ $$\left[\frac{\text{m}}{\text{s}}\right]$$ |
Density $\rho $ | $${\rho}_{4}$$ = 5.7487 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ | $${\rho}_{1}$$ = 0.22995 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ |
Temperature $$T$$ | $${T}_{4}$$ = 303 $\left[\text{K}\right]$ | $${T}_{1}$$ = 303 $\left[\text{K}\right]$ |
High Pressure Side (4) | Low Pressure Side (1) | |
---|---|---|
α | $${a}_{4}$$=348.95 $$\left[\frac{\text{m}}{\text{s}}\right]$$ | $${a}_{1}$$= 348.95 $$\left[\frac{\text{m}}{\text{s}}\right]$$ |
R | 287.049 $$\left[\frac{\text{J}}{\text{kg}\cdot \text{K}}\right]$$ |
Analytical Solution | Results at t = 0.4 ms | |
---|---|---|
Pressure $$p$$ | $\frac{{p}_{4}}{{p}_{1}}=\frac{{p}_{2}}{{p}_{1}}{\left(1-\frac{\left(\gamma -1\right)\cdot \left(\raisebox{1ex}{${a}_{2}$}\!\left/ \!\raisebox{-1ex}{${a}_{1}$}\right.\right)\cdot \left(\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.-1\right)}{\sqrt{2\gamma \left[2\gamma +\left(\gamma +1\right)\cdot \left(\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.-1\right)\right]}}\right)}^{\frac{-2\gamma}{\gamma -1}}$ | $${p}_{2}$$ = 80941.1 $$\left[\text{Pa}\right]$$ |
Velocity $v$ | ${v}_{2}=\frac{{a}_{1}}{\gamma}\left(\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.-1\right){\left(\frac{\raisebox{1ex}{$2\gamma $}\!\left/ \!\raisebox{-1ex}{$\left(\gamma +1\right)$}\right.}{\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.+\left(\gamma -1\right)/\left(\gamma +1\right)}\right)}^{\frac{1}{2}}$ | ${v}_{2}=399.628$ $$\left[\frac{\text{m}}{\text{s}}\right]$$ |
Density $\rho $ | $${\rho}_{2}={\rho}_{2}R{T}_{2}$$ | $${\rho}_{2}$$ = 0.5786 $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ |
Temperature $$T$$ | $\frac{{T}_{1}}{{T}_{2}}=\frac{{p}_{2}}{{p}_{1}}\left(\frac{\raisebox{1ex}{$\left(\gamma +1\right)$}\!\left/ \!\raisebox{-1ex}{$\left(\gamma -1\right)$}\right.+\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.}{1+\left(\raisebox{1ex}{${p}_{2}$}\!\left/ \!\raisebox{-1ex}{${p}_{1}$}\right.\right)\cdot \left(\raisebox{1ex}{$\left(\gamma +1\right)$}\!\left/ \!\raisebox{-1ex}{$\left(\gamma -1\right)$}\right.\right)}\right)$ | ${T}_{2}=487.308$ $\left[\text{K}\right]$ |
Analytical Solution | Results at t = 0.4 ms | |
---|---|---|
Pressure $$p$$ | $${p}_{2}={p}_{2}\text{'}$$ | $${p}_{2}\text{'}=80941.1\text{\hspace{0.17em}}$$ $$\left[\text{Pa}\right]$$ |
Velocity $v$ | ${v}_{2}={v}_{2\text{'}}$ | ${v}_{2\text{'}}=399.628$ $$\left[\frac{\text{m}}{\text{s}}\right]$$ |
Density $\rho $ | $${\rho}_{2\text{'}}={\rho}_{3}\left({x}_{4/3}\right)$$ | $${\rho}_{2\text{'}}=1.5657\text{\hspace{0.17em}}$$ $\left[\frac{\text{k}\text{g}}{\text{m}{\text{m}}^{3}}\right]$ |
Temperature $$T$$ | $${\rho}_{2\text{'}}={r}_{2\text{'}}R{T}_{2\text{'}}$$ | $${T}_{2\text{'}}=180.096\text{\hspace{0.17em}}$$ $\left[\text{K}\right]$ |
Surface contact speed: ${v}_{c}-{v}_{2}$
Therefore, ${x}_{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$2\text{'}$}\right.}=0.4{v}_{s}+500=659.85\left[\mathrm{mm}\right]$
Zone 3
Where, $x=X+500$
At ${v}_{3}-{a}_{3}=-{a}_{4}=-348.95$ $$\left[\frac{\text{m}}{\text{s}}\right]$$ $\Rightarrow X=-348.95t$
$\Rightarrow {x}_{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(t=0.4\right)=360.42$ $\left[\text{mm}\right]$
At ${v}_{3}-{a}_{3}={v}_{2}-{a}_{2\text{'}}={v}_{2}-\sqrt{\frac{{p}_{2}\gamma}{{p}_{2\text{'}}}}=130.602$ $$\left[\frac{\text{m}}{\text{s}}\right]$$ $\Rightarrow -348.95t\le X\le 130.602t$
Analytical Solution | Results at t = 0.4 ms | |
---|---|---|
Pressure $$p$$ | $\frac{{p}_{3}}{{p}_{4}}={\left(1-\frac{\left(\gamma -1\right)}{2}\frac{{v}_{3}}{{a}_{4}}\right)}^{\frac{2\gamma}{\gamma -1}}$ | ${p}_{3}=500000{\left(1-0.2\left(\frac{{v}_{3}\left(X\right)}{348.95}\right)\right)}^{7}$ |
Velocity $v$ | ${v}_{3}=\frac{2}{\gamma +1}\left({a}_{4}+\frac{X}{t}\right)$ | ${v}_{3}=290.792+2.0833X$ |
Density $\rho $ | $\frac{{\rho}_{3}}{{\rho}_{4}}={\left(\frac{{p}_{3}}{{p}_{4}}\right)}^{-\gamma}$ | ${\rho}_{3}=5.7487{\left(\frac{{p}_{3}\left(X\right)}{500000}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$1.4$}\right.}$ |
Temperature $$T$$ | $\frac{{p}_{3}}{{p}_{4}}={\left(\frac{{T}_{3}}{{T}_{4}}\right)}^{\raisebox{1ex}{$\gamma $}\!\left/ \!\raisebox{-1ex}{$\gamma -1$}\right.}$ | ${T}_{3}=303{\left(\frac{{p}_{3}\left(X\right)}{500000}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3.5$}\right.}$ |
Finite Element Modeling with Lagrangian and Eulerian Formulations
Gas is modeled by 200 ALE bricks with solid property TYPE14 (general solid).
In the Lagrangian formulation, the mesh points remain coincident with the material points and the elements deform with the material. Since element accuracy and time step degrade with element distortion, the quality of the results decreases in large deformations.
In the Eulerian formulation, the coordinates of the element nodes are fixed. The nodes remain coincident with special points. Since elements are not changed by the deformation material, no degradation in accuracy occurs in large deformations.
The Lagrangian approach provides more accurate results than the Eulerian approach, due to taking into account the solved equations number.
- Material velocity
- Grid velocity
The nodes on extremities have material velocities fixed in X and Z directions. The other nodes have material and velocities fixed in X, Y and Z directions.
The ALE materials have to be declared Eulerian or Lagrangian with /ALE/MAT.
Results
Finite Element Results and Analytical Solution Comparison
- Lagrangian Formulation
- Scale factor
- $\text{\Delta}{t}_{sca}=0.1$
- Eulerian Formulation
- Scale factor
- $\text{\Delta}{t}_{sca}=0.5$