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.

This famous experiment is interesting for observing the shock-wave propagation. Moreover, this case uses the representation of perfect gas and compares the different formulations: The ALE uses Lagrangian or Eulerian.

rad_ex_13_shock-tube
Figure 1.
The first part of the study deals with the modeling description of perfect gas with the hydrodynamic viscous fluid LAW6. The purpose is to test the different formulations:
  • 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

Input Files

Refer to Access the Model Files to download the required model file(s).

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.

A shock tube consists of a long tube filled with the same gas in two different physical states. The tube is divided into two parts, separated by a diaphragm. The initial state is defined by the values for density, pressure and velocity, as shown in Figure 2 and Figure 3. All the viscous effects are negligible along the tube sides; it is also assumed that there is no motion in the beginning.

rad_ex_fig_13-1
Figure 2. Shock Tube Sketch

rad_ex_fig_13-2
Figure 3. Initial States with Discontinuities

The initial state at time t = 0 consists of two constant states 1 and 4 with p4>p1p4>p1, ρ4>ρ1ρ4>ρ1, and v4=v1=0v4=v1=0 (table).

Table 1. Initial Conditions in the Shock Tube
  High Pressure Side (4) Low Pressure Side (1)
Pressure pp 500000 [Pa][Pa] 20000 [Pa][Pa]
Velocity vv 0 [ms][ms] 0 [ms][ms]
Density ρρ 5.7487 [kgmm3][kgmm3] 0.22995 [kgmm3][kgmm3]
Temperature TT 303 [K][K] 303 [K][K]

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.

The general equation describing pressure is:(1) p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ)Ep=C0+C1μ+C2μ2+C3μ3+(C4+C5μ)E

With μ=ρρ01μ=ρρ01

Where,
pp
Pressure
CiCi
Hydrodynamic constants
EnEn
Internal energy per initial volume
ρρ
Density
Perfect gas is modeled by setting all coefficients:(2) C0=C1=C2=C3=0C0=C1=C2=C3=0 And

C4=C5=γ1C4=C5=γ1

Where, γγ is the gas constant.

Then the initial internal energy, per initial volume is calculated from initial pressure:(3) E0=p0γ1E0=p0γ1

Under the assumption γ=const.=1.4γ=const.=1.4 (valid for low temperature range), the hydrodynamic constants C4=C5=0.4C4=C5=0.4.

Gas pressure is described by:(4) p=(C4+C5μ)Ep=(C4+C5μ)E (5) p=(0.4+0.4ρρ0ρ0)Ep=(0.4+0.4ρρ0ρ0)E (6) E0=p00.4E0=p00.4

Parameters of material LAW6 are provided in Table 2.

Table 2. Material Properties of Gas in LAW6
  High Pressure Side (4) Low Pressure Side (1)
Initial volumetric energy density (E0) 1.25x106[Jmm3][Jmm3] 5x104 [Jmm3][Jmm3]
C4C4 and C5C5 0.4 0.4
Density ρρ 5.7487 [kgmm3][kgmm3] 0.22995[kgmm3][kgmm3]

Analytical Approach

The shock tube problem has an analytical solution of time before the shock hits the extremity of the tube. 1

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

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.

Each wave pattern is composed of a contact discontinuity in the middle and a shock or a rarefaction wave on the left and the right sides separating the uniform state solution. The shock wave moves at a supersonic speed into the low pressure side. A one-dimensional problem is considered.

rad_ex_fig_13-4
Figure 5. Shock Diagram, Expansion Waves and Contact Surface
There are four distinct zones marked 1, 2, 3 and 4 in Figure 5.
  • 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 γγ 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.

The analytical solution to the Riemann problem is indicated at t=0.4 ms. A solution is given according to the distinct zones and continuity must be checked. Evolution in Zones 2 and 3 is dependent on the constant conditions of Zone 1 and 4. The analytical equations use pressure, velocity, density, temperature, speed of sound through gas and a specific gas constant. Equations in Zone 2 are obtained using normal shock relations and the gas velocity in Zone 2 is constant throughout. The shock wave and the surface contact speeds make it possible to define the position of the zone limits.
Zone 4 Zone 1
Pressure pp p4p4 = 500000 [Pa][Pa] p1p1 = 20000 [Pa][Pa]
Velocity vv v4=0v4=0[ms][ms] v1=0v1=0 [ms][ms]
Density ρρ ρ4ρ4 = 5.7487 [kgmm3][kgmm3] ρ1ρ1 = 0.22995 [kgmm3][kgmm3]
Temperature TT T4T4 = 303 [K][K] T1T1 = 303 [K][K]
Speed of sound through gas:(7) a=pγρa=pγρ
Specific gas constant:(8) R=a2TγR=a2Tγ
High Pressure Side (4) Low Pressure Side (1)
α a4a4=348.95 [ms][ms] a1a1= 348.95 [ms][ms]
R 287.049 [JkgK][JkgK]
Table 3. Zone 2
Analytical Solution Results at t = 0.4 ms
Pressure pp p4p1=p2p1(1(γ1)(a2a1)(p2p11)2γ[2γ+(γ+1)(p2p11)])2γγ1p4p1=p2p1(1(γ1)(a2/a1)(p2/p11)2γ[2γ+(γ+1)(p2/p11)])2γγ1 p2p2 = 80941.1 [Pa][Pa]
Velocity vv v2=a1γ(p2p11)(2γ(γ+1)p2p1+(γ1)/(γ+1))12v2=a1γ(p2/p11)(2γ/(γ+1)p2/p1+(γ1)/(γ+1))12 v2=399.628v2=399.628 [ms][ms]
Density ρρ ρ2=ρ2RT2ρ2=ρ2RT2 ρ2ρ2 = 0.5786 [kgmm3][kgmm3]
Temperature TT T1T2=p2p1((γ+1)(γ1)+p2p11+(p2p1)((γ+1)(γ1)))T1T2=p2p1((γ+1)/(γ1)+p2/p11+(p2/p1)((γ+1)/(γ1))) T2=487.308 [K]
Shock wave speed:(9) vs=a1γ+12γ(p2p11)+1=663.166
[ms]
Therefore, x21=0.4vs+500=765.266[mm]
Table 4. Zone 2'
Analytical Solution Results at t = 0.4 ms
Pressure p p2=p2' p2'=80941.1 [Pa]
Velocity v v2=v2' v2'=399.628 [ms]
Density ρ ρ2'=ρ3(x4/3) ρ2'=1.5657 [kgmm3]
Temperature T ρ2'=r2'RT2' T2'=180.096 [K]

Surface contact speed: vcv2

Therefore, x22'=0.4vs+500=659.85[mm]

Zone 3

Zone 3 is defined as:(10) a4Xtv3a3

Where, x=X+500

At v3a3=a4=348.95 [ms] X=348.95t

x43(t=0.4)=360.42 [mm]

At v3a3=v2a2'=v2p2γp2'=130.602 [ms] 348.95tX130.602t

x32'(t=0.4)=552.24 [mm]
Table 5. Zone 3
Analytical Solution Results at t = 0.4 ms
Pressure p p3p4=(1(γ1)2v3a4)2γγ1 p3=500000(10.2(v3(X)348.95))7
Velocity v v3=2γ+1(a4+Xt) v3=290.792+2.0833X
Density ρ ρ3ρ4=(p3p4)γ ρ3=5.7487(p3(X)500000)11.4
Temperature T p3p4=(T3T4)γγ1 T3=303(p3(X)500000)13.5
Continuity verifications:(11) v3(X32')=v2'(X32') (12) v3(X43)=v4(X43) (13) p3(X32')=p2'(X32') (14) p3(X43)=p4(X43) (15) ρ3(X32')=ρ2'(X32') (16) ρ3(X43)=ρ4(X43) (17) T3(X32')=T2'(X32') (18) T3(X43)=T4(X43)

Finite Element Modeling with Lagrangian and Eulerian Formulations

Gas is modeled by 200 ALE bricks with solid property TYPE14 (general solid).

The model consists of regular mesh and elements, the size of which is 5 mm x 5 mm x 5 mm.

rad_ex_fig_13-5
Figure 6. Mesh Used for Lagrangian and Eulerian Approaches

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.

For the ALE boundary conditions (/ALE/BCS), constraints are applied on:
  • 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

Simulation results along the tube axis at 0.4 ms are shown in the following diagrams.

ex_13_pressureA
Figure 7. Pressure

ex_13_densityA
Figure 8. Density

ex_13_velocityA
Figure 9. Velocity
Lagrangian Formulation
Scale factor
Δtsca=0.1
Eulerian Formulation
Scale factor
Δtsca=0.5

Pressure Distribution

Pressure wave produced in the shock-tube at t = 0.4 ms
ex_13_eulerian_pressure ex_13_pressure-lang
Figure 10. Different Approaches and Animations Regarding Pressure, Density and Velocity
1
J. D. Anderson Jr., Modern Compressible Flow with Historical Perspective, McGraw Hill Professional Publishing, 2nd ed., Oct. 1989