GDTk

L1d Reference Manual

L1d is a flow simulation code for quasi-one-dimensional gas slugs coupled to pistons. It turns out to be good at modelling the dynamics of free-piston-driven shock tunnels and expansion tubes.

The set up of a simulation involves writing a Python input script that defines a tube of varying area that contains one or more GasSlug objects and zero or more Piston objects. These components are coupled at their ends and given initial states. Over small time steps, the L1d program then updates the flow state in each cell within the gas slugs according to the constraints of mass, momentum and energy. The details of the gas state at various locations can be used to give a prediction of the performance of the physical machine.

The following sections provide brief details on many items that might go into your input script.

1. Example

When setting up a new simulation, first define the tube as a set of (x,d) break-points and identify regions of head-loss and regions where the wall-temperature varies from the nominal value. Create the GasSlugs, Pistons, and Diaphragms that will make up the gas path. Note that places where two GasSlugs join will need a GasInterface to be defined. Once all of the components have been created, assemble the gas path and then set any of the time-stepping parameters for which you want values other than the default.

Here is an example input script for the Sod shock-tube problem.

# sod.py
config.title = 'Sods ideal shock tube, 2020-04-04'
my_gm = add_gas_model('ideal-air-gas-model.lua')

# Define the tube walls.
add_break_point(0.0, 0.01)
add_break_point(3.0, 0.01)

# Create the gas-path.
left_wall = VelocityEnd(x0=0.0, vel=0.0)
driver_gas = GasSlug(p=100.0e3, vel=0.0, T=348.4, gmodel_id=my_gm, ncells=200)
interface = GasInterface(x0=0.5)
driven_gas = GasSlug(p=10.0e3, vel=0.0, T=278.7, gmodel_id=my_gm, ncells=100)
right_wall = VelocityEnd(x0=1.0, vel=0.0)
assemble_gas_path(left_wall, driver_gas, interface, driven_gas, right_wall)

# Set some time-stepping parameters
config.dt_init = 1.0e-7
config.max_time = 0.6e-3
config.max_step = 5000
add_dt_plot(0.0, 10.0e-6, 5.0e-6)
add_history_loc(0.7)

This script should define the gas path

       |                        |                        |
       +------ driver_gas ------+------ driven_gas ------+
       |                        |                        |
   left_wall                interface               right_wall

and can be invoked with the command

$ l1d4-prep --job=sod

Upon getting to the end of the user’s script, this program should then write

  1. a complete simulation parameter file ./sod/config.json

  2. A tube-definition file ./sod/tube.data

  3. State files for pistons, diaphragms and gas slugs.

Note that Python is very picky about whitespace. If you cut and paste the example from above, make sure that the lines start in the first column and that indentation is consistent with Python’s syntax rules.

2. Configuration options

There are a large number of configuration options that can be set in the input script. The options are set in the input script by adding lines of the form:

config.option = value

Here are all of the available configuration options and the default values if left unset. Note you do not have to set all of these values in the input script.

title

string, default: ""
Short title string for embedding in the parameter and solution files.

gmodels

list of GasModel objects.
You should add GasModels via the add_gas_model function described below but this list will be accessible in case you wish to do some gas state calculations in your input script.

reacting

bool, default: False
If set to True, the finite-rate chemistry will be active.

T_frozen

float, default: 300.0
Temperature above which finite-rate reactions may happen.

dt_init

float, default: 1.0e-6
The size of the time-step that will be used for the first simulation step unless the code determines that a smaller time-step is required.. After a few steps, the cfl condition takes over the determination of a suitable time-step size.

max_time

float, default: 1.0e-3
The simulation will stop if it reaches this time. It is most usual to use this critereon to stop the simulation.

max_step

int, default: 10
The simulation will be stopped if it reaches this number of steps. This is mostly used to catch the problem of the calculation taking a very long time (measured by one’s patience), possibly because the time-step size has decreased to an extremely small value.

cfl_list

list of tuples, default: [(0.0, 0.5),]
The cfl_value is the ratio of the selected time-step size divided by the allowed time-step size. The time-step size is adjusted to ensure that this cfl_value is not exceeded in any particular cell. A typical value of 0.25 seems to work well for simulations with sudden events such as diaphragm bursting, while a value as high as 0.5 may be considered for well-behaved flows. For challenging flows, with very stiff thermochemistry, you might require quite low cfl values (say 0.1). Because the simulation may have times when there is not much happening and other times when the dynamics are interesting, you can schedule different values of cfl. Add entries to this list with the add_cfl_value function.

t_order

int, default: 2
1=Euler time-stepping. This is generally cheap and cheerful.
2=predictor-corrector time-stepping, nominally second order. It is, however, twice as CPU intensive as Euler time-stepping.

x_order

int, default: 2
1=use cell averages without high-order reconstruction. Use this only if the second-order calculation is showing problems.
2=use limited reconstruction (nominally second order).

dt_plot_list

list of tuples, default: [(0.0, max_time, max_time),]
Specifies the frequency of writing complete solutions (for later plotting, maybe) and also for the writing of data at history locations. It may be convenient to have different frequencies of writing such output at different stages of the simulation. For example, free-piston driven shock tunnels have a fairly long period during which the piston travels the length of the compression tube and then a relatively short period, following diaphragm rupture, when all the interesting things happen. It is good to have low-frequency output during most of the compression process and higher-frequency output starting just before diaphragm rupture. Arranging good values may require some trial and error. Add entries to this list via the add_dt_plot function.

hloc_list

list of floats, default []
List of x-coordinates for the history locations. Add entries via the function add_history_loc.

2.1. Configuration helper functions

The configuration options that are lists of values can be managed with the following functions. The lists start empty and these functions may be called several times to build up the list entries. If you do not explicitly add any entries, the lists will take the default entries specified above.

2.1.1. cfl_list

CFL values for the time-stepping can be scheduled with

add_cfl_value(t, cfl)
t

float, no default.
Time at which cfl value takes effect.

cfl

float, no default.
New cfl value. See the discussion above for selecting a suitable value.

2.1.2. dt_plot_list

To set the frequency for writing whole solution data and (separately) history data, shedule the periods between writes with

add_dt_plot(t, dt_plot, dt_hist)
t

float, no default.
Time at which new increments take effect.

dt_plot

float, no default.
Time interval to subsequent solution write.

dt_hist

float, no default.
Time interval to next write of history data.

2.1.3. hloc_list

To specify a location at which you want the history data recorded, use

add_history_loc(x)
x

float, no default.
Location, in metres, along the tube at which the flow data will be samples a history written.

3. Gas Models

There may be one or more gas models involved in a simulation. You have to specify one when you make each GasSlug. To initialize a gas model, call the add_gas_model function.

my_gm = add_gas_model(fileName, reaction_file_1="", reaction_file_2="")

Input:

fileName

string, no default.
Name of the detailed-gas-model file.

reaction_file_1

string, default "".
Name of the detailed chemistry file for reacting gas.

reaction_file_2

string, default "".
Name of the second thermochemistry file. This second thermochemistry file is needed for only a few of the multi-T models.

Returns:

the index of the initialized gas model object. You will need this index to specify which gas model each gas slug is to use.

4. Tube

In a simulation, there is a single tube object that contains the area, and loss specification. The user’s script does not create one of these objects directly but should specify the tube details by calling the add_xxxx functions.

The following attributes are stored in the Tube object:

n

int, default: 4000
The number of small segments that will be used to describe the tube’s area distribution internal to the simulation. To enable a fast lookup process for the area calculation, the area variation between equally-spaced x-positions is taken to be linear. The default value probably won’t need to be changed except for geometries with rapidly changing cross-sections.

xd_list

List of break-point tuples defining the tube wall. Add elements to the list via the function add_break_point.

loss_region_list

list of tuples
List of head-loss regions, usually associated with sudden changes in tube cross-section and diaphragm stations. Add regions via the function add_loss_region.

T_nominal

float, default 300.0
The nominal wall temperature (in degrees K) in the absence of a patch of differing temperature.

T_patch_list

list of tuples
Regions of the tube wall that have temperature different to the nominal value can be specified via the function add_T_patch.

viscous_factor_patch_list

list of tuples
List of regions where we wish to modulate/scale the viscous wall effects. Add regions via the function add_vf_patch.

4.1. Diameter description

The tube is described as a set of (x,d)-coordinate pairs that define break points in the profile of the tube wall. You need at least 2 break points to define the tube. Linear variation of diameter between the break points is assumed.

add_break_point(x, d)
x

float x-coordinate, in metres, of the break point.

d

float diameter, in metres, of the tube wall at the break-point.

Returns the number of break points defined so far.

4.2. Head-loss regions

There is a momentum-sink term much like the so-called minor-loss terms in the fluid mechanics text books. The effect of the loss is spread over a finite region so that the cells are gradually affected as they pass through the region

add_loss_region(xL, xR, K)
xL

float Left-end location, in metres, of the loss region.

xR

float Right-end location, in metres, of the loss region.

K

float Head-loss coefficient. A value of 0.25 seems to be good for a reasonably smooth contraction such as the T4 main diaphragm station.

Returns the number of loss regions defined so far.

4.3. Temperature patches

These define sections of the tube where the wall temperature is different from the nominal value.

add_T_patch(xL, xR, T)
xL

float Left-end location, in metres, of the temperature patch.

xR

float Right-end location, in metres, of the temperature patch.

T

float Wall temperature in degrees K.

Returns the number of temperature patches defined so far.

4.4. Viscous-factor patches

These define sections of the tube where the viscous effects are scaled from the nominal value.

add_vf_patch(xL, xR, vf)
xL

float Left-end location, in metres, of the viscous-factor patch.

xR

float Right-end location, in metres, of the viscous-factor patch.

vf

float Viscous-factor for limiting viscous effects at the wall. The nominal value is 1.0, for full viscous effects. A completely inviscid wall has a value of 0.0.

Returns the number of viscous-factor patches defined so far.

5. Dynamic components

5.1. GasSlug

The principal component is a simulation is a gas slug that move back and forth within the tube. The user may create more than one gas slug to describe the initial gas properties throughout the facility.

Note that a slug needs to have appropriate left- and right-end conditions. This is achieved by creating end-condition objects such as FreeEnd and VelocityEnd objects and then assembling the gas-path via a call to the function assemble_gas_path.

my_slug = GasSlug(gmodel_id=my_gm, p=100.0e3, T=300.0, vel=0.0, massf=[1.0,],
                  ncells=10, cluster_strength=0.0,
                  viscous_effects=0, adiabatic=False,
                  hcells=[])

Most parameters have default properties so that only the user needs to override the ones that they wish to set differently. Note that the locations of the ends of the slug are communicated through end-condition objects that are attached during assembly of the gas path.

gmodel_id

int, default: None
index of the gas-model file name. You must specify a particular gas model.

p

float, default 100.0e3
Pressure in Pa.

T

float, default: 300.0
Thermal temperature, in degrees K.

T_modes

list of float, default: [] Temperatures, in K, for the other internal energy modes, if relevant. If the gas model does include other energy modes and you do not specify values for them, the thermal temperature, T, will be used.

massf

Mass fractions supplied as a list of floats or a dictionary of species names and floats. The number of mass fraction values should match the number of species expected by the selected gas model. For a single species gas, the default [1.0,] is already set for you.

vel

float, default: 0.0
Velocity in m/s.

label

string, dafault: ""
Optional label for the gas slug.

ncells

int, default: 10
Number of cells within the gas slug.

to_end_L

bool, default: False
Flag to indicate that cells should be clustered to the left end.

to_end_R

bool, default: False
Flag to indicate that cells should be clustered to the right end.

cluster_strength

float, default: 0.0
As this value approaches 1.0 from above, the clustering gets stronger. A value of zero indicates no clustering.

viscous_effects

int, default: 0
A nonzero value activates the viscous effects.
0 = inviscid equations only;
1 = include viscous source terms F_wall, loss, q, friction factor for pipe flow.

adiabatic

bool, default: False
Flag to indicate that there should be no heat transfer at the tube wall.

hcells

Either the index (int) of a single cell or a list of indices of cells for which the data are to be written every dt_his seconds, as set by add_dt_plot. Note that cells are indexed from 0 to ncells-1.

5.2. Piston

The other dynamic component that may travel back and forth in the tube is a piston.

myp = Piston(mass, diam, xL0, xR0, vel0)
mass

float, no default
Mass of piston in kg.

diam

float, no default
Face diameter, metres.

xL0

float, no default
Initial position of left-end, metres. The initial position of the piston centroid is set midway between xL0 and xR0 while piston length is the difference (xR0 - xL0).

xR0

float, no default
Initial position of right-end, metres.

vel0

float, no default
Initial velocity (of the centroid), m/s.

front_seal_f

float, default: 0.0
friction coefficient. A typical value might be 0.2.

front_seal_area

float, default: 0.0
Seal area (in m^2) over which the front-side pressure acts. This is the effective area over which the compressed gas pressed the front-side seal against the tube wall. Friction force is this area multiplied by downstream-pressure by friction coefficient.

back_seal_f

float, default: 0.0
friction coefficient. A typical value might be 0.2.

back_seal_area

float, default: 0.0
Seal area (in m^2) over which the back-side pressure acts. Friction force is this area multiplied by downstream-pressure by friction coefficient. This is for gun tunnel pistons that have flexible skirts that are pressed onto the tube wall by the pushing gas.

p_restrain

float, default: 0.0
Pressure (in Pa) at which restraint will release. Some machines, such as two-stage light-gas guns, will hold the projectile in place with some form of mechanical restraint until the pressure behind the piston reaches a critical value. The piston is then allowed to slide.

is_restrain

int, default: 0
Status flag for restraint. 0=free-to-move, 1=restrained

with_brakes

bool, default: False
Flag to indicate the presence of brakes. Such brakes, as on the T4 shock tunnel, allow free forward motion of the piston but try to prevent backward motion by applying a large frictional force at the tube wall.

brakes_on

int, default: 0
Flag to indicate the state of the brakes. 0=off, 1=on.

brakes_friction_force

float, default: 0.0
The maximum friction force (in Newtons) that the brakes can apply when they are on. This is modelled on the sliding-shoe brakes of the T4 shock tunnel, which are activated by the piston trying to travel backwards up the compression tube. The user will need to supply an estimate of this value, possibly by considering the frontal area of the piston and the maximum pressure that the brakes are expected to hold before slipping.

x_buffer

float, default: 1.0e6
Position of the stopping buffer in metres. This is the location of the piston centroid at which the piston would strike the buffer (or brake, in HEG terminology). Note that it is different to the location of the front of the piston at strike.

hit_buffer

int, default: 0
Flag to indicate state of buffer interaction. A value of 0 indicates that the piston has not (yet) hit the buffer. A value of 1 indicates that it has.

Notes

  1. The left- and right-end positions of the piston are also used to locate the ends of adjoining GasSlugs.

  2. The basic piston model has inertia but no friction. To make accurate simulations of a particular facility, it is usually important to have some account of the friction caused by gas-seals and guide-rings that may be present on the piston.

5.3. End Conditions

The end-conditions for the gas slugs provide on where the end of the gas slug is initially located, as well as what happens as the simulation proceeds. The general procedure is to define the end conditions and later make connections to the gas slugs by assembling the gas path.

5.3.1. Diaphragm

A diaphragm conditionally connects the adjacent ends of two gas slugs.

my_d = Diaphragm(x0, p_burst)
x0

float, no default
x-position in the tube, metres. This value is used to locate the end-points of the gas slugs.

p_burst

float, default: 0.0
Pressure, in Pa, at which rupture is triggered.

dxL

float, default: 0.0
Distance to the left of the diaphragm over which to average the pressure.

dxR

float, default: 0.0
Distance to the right of the diaphragm over which to average the pressure.

state

int: default 0
Indicates the state of diaphragm:
0 indicates that the diaphragm is intact/closed (with zero-velocity end condition being effectively applied to both gas slugs).
1 indicates that the diaphragm has been triggered by over pressure.
2 indicates that the diaphragm is ruptured/open and the gas slugs are interacting.

dt_hold

float, default: 0.0
Time delay, in seconds, from overpressure-trigger to actual rupture with the gas slugs being allowed to interact.

5.3.2. GasInterface

Connects two gas slugs at specified location.

my_if = GasInterface(x0)
x0

float, no default
Initial position, in metres.

5.3.3. FreeEnd

Initially locates the end of the gas slug but otherwise does not constrain it.

my_fe = FreeEnd(x0)
x0

float, no default
Initial position, in metres.

5.3.4. VelocityEnd

Specify the initial location and velocity of the end of a gas slug. To model a fixed wall, just specify a zero velocity.

my_ve = VelocityEnd(x0, vel=0.0)
x0

float, no default
Initial position, in metres.

vel

float, no default
Velocity, in m/s, of the end-point of the gas slug.

5.3.5. PistonFace

Objects of this class connect the end of a GasSlug to a Piston face. Normally, you do not need to construct these objects explicitly because the assemble_gas_path function will create and connect them implicitly.

PistonFace()

5.4. Valves

Valve objects may be constructed to restrict the motion of gas slugs at some internal point. A valve is located at a particular x-location and has a history of open-area-fraction specified by a table of times with particular fopen values. At any time, the fraction of open area is interpolated linearly from this table. For example:

v = Valve(x=-1.0, times=[0.020, 0.030], fopen=[0.0, 1.0])

This defines a valve that remains closed up until t=0.020 and is fully open by t=0.030. For this simple example, it is easy to write the sequences for times and fopen explicitly, however, more sophisticated openings may be conveniently generated via a Python function in your input script. The details are up to you but the use of numpy arrays and their associated array functions will probably be handy.

The effect of the valve is applied to the nearest internal face of a gas slug that covers the valve location. For this closest interface, the normal (fully-open) Riemann solver result of intermediate pressures and and velocity is blended with the fully-closed Riemann-solver result (of differing stagnation pressures and zero velocity) according to the fraction of open area.

6. Assembling a gas path

Assemble a gas path by making the logical connections between adjacent components. The components are assembled left-to-right, as they are supplied to the following function.

assemble_gas_path(*components)
components

An arbitrary number of arguments representing individual components or lists of components. Each component may be a GasSlug, Piston, or any other gas-path object, however, it doesn’t always make sense to connect arbitrary components. For example, connecting a GasSlug to a Piston is reasonable but connecting a Piston to a Diaphragm without an intervening GasSlug does not make sense in the context of this simulation program. Also, valves are not part of the gas path.

If you really want to make a connection manually, there is a function available to make the logical connection between a pair of components.

connect_pair(cL, cR)
cL

component object on left

cR

component object on right

7. Command-line usage

The Lagrangian simulation tools consist of two programs: l1d4-prep and l1d4. Before preparing a simulation, you need one or more detailed gas-model file(s) and, if relevant, one or more detailed chemistry files. Once you have your gas-model file(s) and input script, as described above, you are ready to simulate.

7.1. Preparation of gas files

7.1.1. prep-gas : gas model preparation

The prep-gas program is used to take a brief description of the desired gas model used in the flow simulation and produce a detailed configuration gas model configuration file for use by Eilmer at pre-processing and simulation stages. Its usage is shown here. Generally, one uses prep-gas in the first mode shown: with two arguments. The second mode simply lists available species in the Eilmer database and exits.

Usage:
 > prep-gas input output

   input    : a gas input file with user selections
   output   : detailed gas file in format ready for Eilmer4.

 > prep-gas --list-available-species

7.1.2. prep-chem : chemistry scheme preparation

prep-chem is used to take user created description of a chemistry scheme written in Lua and generate a detailed configuration file for eilmer to use at run-time. The use of prep-chem is shown here.

Usage:
 > prep-chem [--compact] gmodelfile cheminput output

   gmodelfile  : a gas model file is required as input for context
   cheminput   : input chemistry file in Lua format.
   output      : output file in format ready for Eilmer4.

Options:
   --compact   : produce a text file called 'chem-compact-notation.inp'
                 which is used to configure a GPU chemistry kernel.

7.2. Preparation of initial state

The preprocessing program is written in Python and it accepts the name of your input script, also in Python.

$ l1d4-prep --job=<myjob>

or

$ l1d4-prep --job=<myjob>.py

If your input script, <myjob>.py, is successfully processed, a summary of the objects created and connected is printed to the console. The configuration data and the initial state of the simulation is written to a set of files in newly-created directory <myjob>.

7.3. Running a simulation

Once the configuration and initial-state data are written, the main simulation code may be run.

$ l1d4 --run-simulation --job=<myjob>

By default, the progress of the simulation is printed, in summary form, to the console. The state data for the gas slugs, pistons and diaphragms is appended periodically to the files in the <myjob> directory. In that directory, there is a times.data file that lists the time instants for each time index. Beside the state data files, there are the history files, one for eash location that was specified in the input script. The format of these files is compatible with GNUPlot and has a comment as the first line, to indicate the recorded properties.

7.4. Stopping a simulation

Because the solution files are opened and appended to periodically, you can usually just kill the l1d4 process (possibly with the keyboard interrupt signal, Control-C) and all will be fine. You might be unlucky and interrupt the program while it is writing some solution data but you can then trim the solution files to tidy up.

7.5. Postprocessing

After the simulation has run to completion, the same executable program, l1d4, can be used to select data for display.

7.5.1. Selecting a time-slice of gas-slug data. 

$ l1d4 --time-slice --job=<myjob> --tindx=<int>

The data for a particular time-instant will be written to files in the current directory. The format of these data files is compatible with GNUplot.

7.5.2. History data for a piston. 

$ l1d4 --piston-history --job=<myjob> --pindx=<int>

The history data for a particular piston will be written to a file in the current directory. The format of this data file is compatible with GNUplot.

7.5.3. Generating an xt-data plot.

The history of gas-slug data, for a particular variable, over a range of time instants can be assembled into a xt-data set, with one GNUplot-compatible file for each gas slug.

$ l1d4 --xt-data --job=<myjob> --var-name=p --log10

7.6. Restarting a simulation

Because the time step adapts to the solution, it is a bit difficult to predict how many steps a simulation will require. You may have a simulation that has stopped at the maximum number of steps but you would like to restart it and proceed a bit further.

To prepare for a restart, first decide the final value of tindx that you wish to retain from the initial run and observe the size of the time step at that point in the simulation. You may wish to scan the <myjob>/times.data file to determine a suitable value for tindx and the console output from the previous run will indicate suitable a time-step size.

Edit the <myjob>/config.json file and manually change the values for max_time, max_step and dt_init, so that you are reasonably sure that the simulation will step stably on restart.

Trim the solution files to the value of tindx, from which you want to restart the simulation.

$ l1d4 --trim-solution-files --job=<myjob> --tindx-end=<int>

If you make a mistake, all is not lost at this point because the original solution files have been retained with a .backup extension to their names.

You should now be ready to restart the simulation, specifying which value of `tindx' to use when reading the initial state.

$ l1d4 --run-simulation --job=<myjob> --tindx=<int>