ESTCN is a state-to-state calculation program for estimating flow conditions in reflected-shock tunnels. With a gas-model file already set up, it is intended to be used as a command-line tool for quick calculations of tunnel conditions.

1. Getting started

The estcn program is built upon the core gas models and the Python library that is wrapped around that gas-model library. This is part of a larger gas-dynamics toolkit and general getting started notes can be found at http://cfcfd.mechmining.uq.edu.au/docs/getting-started There, you will see how to get a copy of the source code, a list of what other software you will need to build and install the tool kit, and a collection of environment variables that need to be set.

To install the estcn program, move to the gas source directory and use the make utility.

cd dgd/src/gas
make install

This will also install files associated with the gas models.

2. Getting command help

When run as an application, this program takes its input as command line arguments, performs the requested calculations and outputs the gas-state results.

To see what specific inputs are required, start the program as:

estcn --help

3. Example shock-tunnel calculation

To do a calculation, you are going to need a gas model file in your working directory. There is a guide (http://cfcfd.mechmining.uq.edu.au/pdfs/gas-user-guide.pdf) to the gas-model module that is worth reading but, for the moment, let’s just copy a gas-model file from the samples in the source code repository.

cp ${DGD_REPO}/src/gas/sample-data/cea-lut-air-version-test.lua ./cea-lut-air.lua

This is a look-up-table model of air in thermochemical equilibrium that has been prepared by running the CEA2 program over a large number of gas states and tabulating the results.

Which particular input parameters you need to supply to estcn depends on the chosen task, however, a typical low-enthalpy flow condition for the T4 shock tunnel may start with a test gas (air) at room temperature (T1=300 K) and a little above atmospheric pressure (p1=125 kPa). The observed shock speed, Vs, for this particular shot was 2414 m/s and the observed nozzle-supply pressure relaxed to 34.37 MPa. With the old Mach 4 nozzle having an area ratio of 27, the flow conditions in the facility may be computed using the command:

estcn --task=stn --gas=cea-lut-air.lua \
      --T1=300 --p1=125.0e3 --Vs=2414 --pe=34.37e6 --ar=27.0

Note that this command is given on a single logical line, with the \ character indicating continuation of the line.

The gas is assumed to remain in thermochemical equilibrium and the flow processing is done in decoupled quasi-one-dimensional wave processes such as shock waves and expansion fans. For the reflected-shock tunnel, this means that the initial, quiescent test gas (state 1) is first processed by the incident shock (to get to state 2) and subsequently by the reflected shock (to get to state 5). The incident shock sets the inflow conditions for the reflected shock but there is no further interaction.

shock tube process fig from stn report

This condition has an enthalpy of 5.43 MJ/kg. (Look for the label (H5s-H1) in the console output.) The nozzle-exit condition, labelled as State 7, has a pressure of 93.6 kPa and a static temperature of 1284 degrees K, with a flow speed of 2.95 km/s (V7).

The full console output is shown below. Note that some of the lines are quite long and may be wrapped in the HTML view below and in your console.

Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-lut-air.lua, p1: 125000 Pa, T1: 300 K, massf=[1.0], Vs: 2414 m/s
Write pre-shock condition.
Start incident-shock calculation.
Start reflected-shock calculation.
Start calculation of isentropic relaxation.
Start isentropic relaxation to throat (Mach 1)
Start isentropic relaxation to nozzle exit of given area.
Done with reflected shock tube calculation.
State 1: pre-shock condition
  p: 125000 Pa, T: 300 K, rho: 1.45152 kg/m**3, u: 215959 J/kg, h: 302075 J/kg
  R: 287.055 J/(kg.K), gam: 1.39053, Cp: 1022.1 J/(kg.K), a: 345.998 m/s, s: 6796.3 J/(kg.K)
State 2: post-shock condition.
  p: 7.3156e+06 Pa, T: 2630.41 K, rho: 9.68285 kg/m**3, u: 2.39478e+06 J/kg, h: 3.1503e+06 J/kg
  R: 287.225 J/(kg.K), gam: 1.28907, Cp: 1280.85 J/(kg.K), a: 971.095 m/s, s: 8128.65 J/(kg.K)
  V2: 361.874 m/s, Vg: 2052.13 m/s
State 5: reflected-shock condition.
  p: 5.94876e+07 Pa, T: 4551.26 K, rho: 44.3175 kg/m**3, u: 5.09065e+06 J/kg, h: 6.43295e+06 J/kg
  R: 294.93 J/(kg.K), gam: 1.28602, Cp: 1326.08 J/(kg.K), a: 1277.77 m/s, s: 8446.72 J/(kg.K)
  Vr: 573.6 m/s
State 5s: equilibrium condition (relaxation to pe)
  p: 3.437e+07 Pa, T: 4160.97 K, rho: 28.207 kg/m**3, u: 4.51266e+06 J/kg, h: 5.73115e+06 J/kg
  R: 292.838 J/(kg.K), gam: 1.2852, Cp: 1319.62 J/(kg.K), a: 1215.36 m/s, s: 8446.72 J/(kg.K)
Enthalpy difference (H5s - H1): 5.42908e+06 J/kg
State 6: Nozzle-throat condition (relaxation to M=1)
  p: 1.93221e+07 Pa, T: 3787.56 K, rho: 17.5341 kg/m**3, u: 3.96172e+06 J/kg, h: 5.06369e+06 J/kg
  R: 290.946 J/(kg.K), gam: 1.28474, Cp: 1312.73 J/(kg.K), a: 1155.39 m/s, s: 8446.72 J/(kg.K)
  V6: 1155.39 m/s, M6: 0.999999, mflux6: 20258.7 kg/s/m**2
State 7: Nozzle-exit condition (relaxation to correct mass flux)
  p: 93702.4 Pa, T: 1283.58 K, rho: 0.254313 kg/m**3, u: 1.01045e+06 J/kg, h: 1.37891e+06 J/kg
  R: 287.051 J/(kg.K), gam: 1.31935, Cp: 1185.91 J/(kg.K), a: 696.505 m/s, s: 8446.72 J/(kg.K)
  V7: 2950.34 m/s, M7: 4.23591, mflux7: 20258.3 kg/s/m**2, area_ratio: 27, pitot: 2.14969e+06 Pa
  pitot7_on_p5s: 0.0625456

For this particular example, we have selected to stop the expansion at a particular nozzle area ratio. Alternatively, we may stop the expansion at a particular Pitot pressure by specifying --task=stnp and a suitable ratio for the option --pp_on_pe. If you don’t want to specify a relaxation pressure with option --pe, the reflected-shock conditions (state 5) will be used directly as the nozzle supply conditions.

3.1. Getting species mass fractions

If you are interested in the chemical species fractions within the air test gas, you can do the same state-to-state calculation with a gas model that more-directly uses the NASA Glenn CEA2 program. This time, copy the gas model file:

cp ${DGD_REPO}/src/gas/sample-data/cea-air5species-gas-model.lua .

Note that you need to have the CEA2 executable program in the installation directory, along with its database input files. Since CEA2 is not ours to give away, you need to get it from an appropriate source.

Once you have your copy of the CEA2 program in place, run the same shock tunnel calculation (for the same conditions as above) with the command:

estcn --task=stn --gas=cea-air5species-gas-model.lua \
      --T1=300 --p1=125.0e3 --Vs=2414 --pe=34.37e6 --ar=27.0

This time, the calculation takes a bit longer because our gas-model code is calling out to the CEA2 program for the gas behaviour but you will now get the mass fractions of the chemical species for air at each of the states. Look for the dictionary labelled CEA-massf for each gas state in the console output below.

Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-air5species-gas-model.lua, p1: 125000 Pa, T1: 300 K, massf=[1.0], Vs: 2414 m/s
Write pre-shock condition.
Start incident-shock calculation.
Start reflected-shock calculation.
Start calculation of isentropic relaxation.
Start isentropic relaxation to throat (Mach 1)
Start isentropic relaxation to nozzle exit of given area.
Done with reflected shock tube calculation.
State 1: pre-shock condition
  p: 125000 Pa, T: 300 K, rho: 1.4458 kg/m**3, u: -84587 J/kg, h: 1871.1 J/kg
  R: 288.198 J/(kg.K), gam: 1.3985, Cp: 1011.4 J/(kg.K), a: 347.7 m/s, s: 6830.1 J/(kg.K)
  CEA-massf: {'O': 0.0, 'NO': 0.0, 'O2': 0.23292, 'N2': 0.76708, 'N': 0.0}
State 2: post-shock condition.
  p: 7.2897e+06 Pa, T: 2615.79 K, rho: 9.66369 kg/m**3, u: 2.096e+06 J/kg, h: 2.85035e+06 J/kg
  R: 288.378 J/(kg.K), gam: 1.24383, Cp: 1471.1 J/(kg.K), a: 970.2 m/s, s: 8158.4 J/(kg.K)
  CEA-massf: {'O': 0.00071876, 'NO': 0.027856, 'O2': 0.21735, 'N2': 0.75408, 'N': 1.4384e-06}
  V2: 361.162 m/s, Vg: 2052.84 m/s
State 5: reflected-shock condition.
  p: 5.9375e+07 Pa, T: 4529.8 K, rho: 44.284 kg/m**3, u: 4.79256e+06 J/kg, h: 6.13334e+06 J/kg
  R: 295.995 J/(kg.K), gam: 1.16491, Cp: 2090.9 J/(kg.K), a: 1276.8 m/s, s: 8478.4 J/(kg.K)
  CEA-massf: {'O': 0.029813, 'N': 0.00016904, 'O2': 0.13695, 'NO': 0.12407, 'N2': 0.70899}
  Vr: 572.859 m/s
State 5s: equilibrium condition (relaxation to pe)
  p: 3.437e+07 Pa, T: 4143.27 K, rho: 28.225 kg/m**3, u: 4.21674e+06 J/kg, h: 5.43447e+06 J/kg
  R: 293.903 J/(kg.K), gam: 1.16797, Cp: 2043.6 J/(kg.K), a: 1214.8 m/s, s: 8478.4 J/(kg.K)
  CEA-massf: {'N2': 0.71741, 'O': 0.021906, 'O2': 0.15435, 'N': 6.6776e-05, 'NO': 0.10627}
Enthalpy difference (H5s - H1): 5.43258e+06 J/kg
State 6: Nozzle-throat condition (relaxation to M=1)
  p: 1.93258e+07 Pa, T: 3771.61 K, rho: 17.547 kg/m**3, u: 3.66604e+06 J/kg, h: 4.76744e+06 J/kg
  R: 292.024 J/(kg.K), gam: 1.17551, Cp: 1955.9 J/(kg.K), a: 1155 m/s, s: 8478.4 J/(kg.K)
  CEA-massf: {'N2': 0.72633, 'O': 0.014711, 'O2': 0.17168, 'N': 2.2302e-05, 'NO': 0.087255}
  V6: 1155.02 m/s, M6: 1.00001, mflux6: 20267.1 kg/s/m**2
State 7: Nozzle-exit condition (relaxation to correct mass flux)
  p: 93940.5 Pa, T: 1280.98 K, rho: 0.25446 kg/m**3, u: 714610 J/kg, h: 1.08378e+06 J/kg
  R: 288.198 J/(kg.K), gam: 1.31535, Cp: 1202.1 J/(kg.K), a: 696.8 m/s, s: 8478.4 J/(kg.K)
  CEA-massf: {'N': 0.0, 'O': 0.0, 'O2': 0.23272, 'NO': 0.00036334, 'N2': 0.76691}
  V7: 2949.81 m/s, M7: 4.23337, mflux7: 20266.4 kg/s/m**2, area_ratio: 27, pitot: 2.15104e+06 Pa
  pitot7_on_p5s: 0.0625849

4. Other examples

Subset calculations of the shock-tube flow processing can be done by selecting a different task. For example, just the incident shock processing can be computed with the --task=ishock, specifying only the gas, initial pressure, temperature and incident shock speed.

4.1. Incident shock

Here is an example from Huber’s Table IV for a speed of 37.06 ft/s at a geopotential altitude of 173500 feet. Note that we expect ionization to be a significant effect at these conditions, so we need to use a CEA2 gas model that includes the relevant chemical species.

cp ${DGD_REPO}/src/gas/sample-data/cea-air13species-gas-model.lua .
estcn --task=ishock --gas=cea-air13species-gas-model.lua --p1=59 --T1=283 --Vs=11296

The expected pressure (from Table IV) is 86.5 kPa and the temperature is 12000 K, quite close to the values computed by estcn and shown below.

Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-air13species-gas-model.lua, p1: 59 Pa, T1: 283 K, massf=[1.0], Vs: 11296 m/s
Write pre-shock condition.
Start incident-shock calculation.
State 1: pre-shock condition
  p: 59 Pa, T: 283 K, rho: 0.00072614 kg/m**3, u: -96470 J/kg, h: -15218 J/kg
  R: 287.113 J/(kg.K), gam: 1.40039, Cp: 1004.2 J/(kg.K), a: 337.3 m/s, s: 8946.9 J/(kg.K)
  CEA-massf: {'e-': 0.0, 'O+': 0.0, 'N2': 0.75566, 'N': 0.0, 'O2+': 0.0, 'N2+': 0.0, 'Ar': 0.01283, 'O': 0.0, 'O2': 0.23151, 'NO': 0.0, 'Ar+': 0.0, 'N+': 0.0, 'NO+': 0.0}
State 2: post-shock condition.
  p: 86686 Pa, T: 12034 K, rho: 0.0111608 kg/m**3, u: 5.57475e+07 J/kg, h: 6.35146e+07 J/kg
  R: 645.436 J/(kg.K), gam: 1.05916, Cp: 11555.9 J/(kg.K), a: 3017.2 m/s, s: 18077.8 J/(kg.K)
  CEA-massf: {'e-': 4.8858e-06, 'O+': 0.020462, 'N2': 0.00040306, 'N': 0.64896, 'O2+': 0.0, 'N2+': 7.0373e-05, 'Ar': 0.011137, 'O': 0.21099, 'O2': 0.23151, 'NO': 3.1054e-05, 'Ar+': 0.0016928, 'N+': 0.10617, 'NO+': 7.1533e-05}
  V2: 734.94 m/s, Vg: 10561.1 m/s

4.2. Pitot pressure

Using the test flow conditions from the exit of the Mach 4 nozzle, we can then compute the expected Pitot pressure to be 2.14 MPa.

estcn --gas=cea-lut-air.lua --task=pitot --p1=93.6e3 --T1=1284 --V1=2.95e3
Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-lut-air.lua, p1: 93600 Pa, T1: 1284 K, massf: [1.0] V1: 2950 m/s
Pitot condition:
  p: 2.1462e+06 Pa, T: 3875.76 K, rho: 1.8455 kg/m**3, u: 4.56771e+06 J/kg, h: 5.73065e+06 J/kg
  R: 300.054 J/(kg.K), gam: 1.29539, Cp: 1315.84 J/(kg.K), a: 1176.2 m/s, s: 9268.14 J/(kg.K)

4.3. Total condition

The hypothetical stagnation conditions for a specified free stream can be computed as:

estcn --gas=cea-lut-air.lua --task=total --p1=93.6e3 --T1=1284 --V1=2.95e3
Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-lut-air.lua, p1: 93600 Pa, T1: 1284 K, massf: [1.0], V1: 2950 m/s
Total condition:
  p: 3.4273e+07 Pa, T: 4160.5 K, rho: 28.1302 kg/m**3, u: 4.51229e+06 J/kg, h: 5.73066e+06 J/kg
  R: 292.842 J/(kg.K), gam: 1.28521, Cp: 1319.61 J/(kg.K), a: 1215.29 m/s, s: 8447.42 J/(kg.K)

4.4. Cone-surface pressure

The conditions on the surface of a conical pressure probe (with 15 degrees half-angle) can be computed as:

estcn --gas=cea-lut-air.lua --task=cone --sigma-deg=15 --p1=93.6e3 --T1=1284 --V1=2.95e3
Equilibrium Shock Tube Conditions, with Nozzle
  Version: 11-Jun-2020
Input parameters:
  gasFileName is cea-lut-air.lua, p1: 93600 Pa, T1: 1284 K, massf: [1.0], V1: 2950 m/s, sigma: 15 degrees
Free-stream condition:
  p: 93600 Pa, T: 1284 K, rho: 0.253952 kg/m**3, u: 1.01083e+06 J/kg, h: 1.3794e+06 J/kg
  R: 287.051 J/(kg.K), gam: 1.31933, Cp: 1185.96 J/(kg.K), a: 696.612 m/s, s: 8447.42 J/(kg.K)
Shock angle: 0.366598 (rad), 21.0045 (deg)
Cone-surface velocity: 2784.53 m/s
Cone-surface condition:
  p: 271013 Pa, T: 1680.12 K, rho: 0.561944 kg/m**3, u: 1.38319e+06 J/kg, h: 1.86546e+06 J/kg
  R: 287.049 J/(kg.K), gam: 1.30525, Cp: 1227.41 J/(kg.K), a: 790.217 m/s, s: 8471.86 J/(kg.K)