The chemistry unit test (mfix/unit-tests/chemistry) can be used to test user-defined reaction rates.
- Create an inputs, build an executable, and run.
- This is a wrapper around the ODE solver.
- Initial condition defines y(t=0)
- Specify the time step (dt)
- Time march: dy/dt = (chemistry source terms)
You may need to tweak the unit test code if you are interested in specific output.
There are a few existing inputs in the test directory (mfix/unit-tests/chemistry) for specific unit tests.
Adiabatic flame example:
Let’s look at the inputs.adiabatic-flame.txt (6.7 KB) setup as an example.
The time marching section defines the duration and time step size.
The chemistry unit test does not solve for convection or diffusion. Therefore, there is no CLF constraint so the user is responsible for specifying a reasonable value.
#_______________________________________________________________________
# Time marching
#
mfix.fixed_dt = 1.0e-3
#mfix.max_step = 25
mfix.stop_time = 25.e-3
Define the chemistry solver, and specify any other options. (chemistry solver documentation)
- Particles update options (e.g., constant volume or constant density)
- Extinguished particle criteria
- Absolute and relative tolerances
- Jacobian evaluation options
#_______________________________________________________________________
# Chemistry solver settings:
#
# Define the integrator: ForwardEuler
mfix.chemistry.integrator = ForwardEuler # default
#mfix.chemistry.integrator = StiffSolver::ForwardEuler
#mfix.chemistry.integrator = StiffSolver::BackwardEuler
#mfix.chemistry.integrator = StiffSolver::VODE
Defining the geometry and grids.
Other than changing the domain size (prob_lo and prob_hi) this section really shouldn’t be changed too much. This is setup to give you a single, fully periodic computational cell.
Remember this unit test is designed for evaluating chemistry, not a full-blown CFD model.
#_______________________________________________________________________
# Domain and geometry:
#
# Number of grid cells in each direction at the coarsest level
amr.n_cell = 1 1 1
amr.blocking_factor = 1
# Maximum allowable size of each subdomain in the problem domain;
# this is used to decompose the domain for parallel calculations.
amr.max_grid_size_x = 16
amr.max_grid_size_y = 16
amr.max_grid_size_z = 16
# Maximum tile size within each grid
fabarray.mfiter_tile_size = 16 16 16
# Maximum particle tile size
particles.tile_size = 16 16 16
# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
amr.max_level = 0
# Geometry
geometry.coord_sys = 0 # 0: Cartesian
geometry.is_periodic = 1 1 1 # Is periodic in each direction?
geometry.prob_lo = 0.0 0.0 0.0 # lo corner of physical domain
geometry.prob_hi = 0.1 0.1 0.1 # hi corner of physical domain
Set what physics to include. It is doubtful that you need to modify this section.
The close-system constraint is used because it is the most consistent with setup, but it may not be a great choice for all tests. You may want to use the incompressible constraint so the dpdt term is not included in the energy equation.
#_______________________________________________________________________
# Model options:
#
mfix.advect_density = 1
mfix.advect_enthalpy = 1
mfix.solve_species = 1
mfix.constraint = IdealGasClosedSystem
Define what species are in your model.
- This should reflect how you intend to setup the species for your full CFD model.
The definition of the NASA-7 polynomials were omitted because they take up too much space.
#_______________________________________________________________________
# Species model settings
#
species.solve = CH4 O2 CO2 H2O N2
species.diffusivity = constant
species.diffusivity.constant = 2.31e-5
species.specific_heat = NASA7-Poly
# !!! NASA7 coefficients omitted for brevity !!!
Define fluid model settings.
- This should reflect how you intend to setup the fluid for your full CFD model.
- The viscosity and thermal conductivity may not be used by the unit test but are required to satisfy checks in the code.
Remember the thermodynamic pressure is what defines the state via the ideal gas law.
#_______________________________________________________________________
# Fluid model settings
#
fluid.solve = fluid
fluid.viscosity.molecular = constant
fluid.viscosity.molecular.constant = 1.0e-5
fluid.reference_temperature = 298.15
fluid.thermodynamic_pressure = 101000.
fluid.specific_heat = mixture
fluid.thermal_conductivity = constant
fluid.thermal_conductivity.constant = 0.024
fluid.species = CH4 O2 CO2 H2O N2
Define solids model settings.
This should reflect how you intend to setup the solids for your full CFD model.
#_______________________________________________________________________
# Solids properties
#
solids.types = None
dem.solve = None
pic.solve = None
Define regions.
- You probably only need a single region that spans the entire domain.
- This region is used to set the initial conditions.
#_______________________________________________________________________
# Regions for defining ICs and BCs
#
mfix.regions = full-domain
regions.full-domain.lo = 0.0 0.0 0.0
regions.full-domain.hi = 0.1 0.1 0.1
Define the chemical reactions
- This example only has one Eulerian reaction.
#_______________________________________________________________________
# Chemical reactions
#
chemistry.solve = AdiabaticFlame
chemistry.AdiabaticFlame.reaction = CH4(g) + 2O2(g) --> CO2(g) + 2H2O(g)
Define the initial conditions
#_______________________________________________________________________
# Initial Conditions
#
ic.regions = full-domain
ic.full-domain.fluid.volfrac = 1.0
ic.full-domain.fluid.velocity = 0.0 0.0 0.0
ic.full-domain.fluid.temperature = 298.15
ic.full-domain.fluid.species.CH4 = 0.046290939
ic.full-domain.fluid.species.O2 = 0.222196505
ic.full-domain.fluid.species.N2 = 0.731512556
Define the reaction rates
Implement reaction rates into the mfix_usr_reactions_rates_K.H file and copy the file into the unit-tests/chemistry/source. For this example, we’ll use a simple arrhenius expression.
// Adiabatic flame:
// CH4_Comb: CH4 + 2O2 --> CO2 + 2H2O (mol/m^3.s)
//--------------------------------------------------------------------//
// Note: The CH4 combustion rate is artificial and used for the
// adiabatic flame test case.
amrex::Real const c_O2 = fluid.molar_concentration("O2"); // [mol/m^3]
amrex::Real const c_CH4 = fluid.molar_concentration("CH4"); // [mol/m^3]
amrex::Real k0(1.0e3); // reaction rate coefficient
reactions.rate(0) = k0 * c_O2 * c_CH4;
Try it out!
Change into the chemistry unit-test directory:
cd mfix/unit-tests/chemistry
Edit GNUmakefile
The mfix_usr_reactions_rates_K.H included in the chemistry/source directory contains reaction rates for several tests. To enable the adiabatic reaction rate, edit the GNUmakefile:
diff --git a/unit-tests/chemistry/GNUmakefile b/unit-tests/chemistry/GNUmakefile
index 7bb9a838e..72ab372d6 100644
--- a/unit-tests/chemistry/GNUmakefile
+++ b/unit-tests/chemistry/GNUmakefile
@@ -27,6 +27,6 @@ USE_CUDA = FALSE
USE_ASSERTION = TRUE
TEST = TRUE
-#CPPFLAGS += -DCHEM_TEST_EULERIAN01
+CPPFLAGS += -DCHEM_TEST_EULERIAN02
Build the unit test executable
make
Run the test adiabatic flame test case.
./mfix.unit-test.chem3d.gnu.DEBUG.ex inputs.adiabatic-flame
I got an unexpected answer! What happened?
|-------+----------+----------+----------+----------+---------|
| X_CH4 | X_O2 | X_CO2 | X_H2O | X_N2 | T_f |
|-------+----------+----------+----------+----------+---------|
| 0. | 0. | 0.392077 | 0.320991 | 0.286931 | 4513.67 |
|-------+----------+----------+----------+----------+---------|
This test is designed to be “too stiff” and fail for the default Forward Euler integrator when we set the time step to fixed_dt=1.0e-3.
Option 1: Reduce the time step to 1e-4
Option 2: Use the StiffSolver::ForwardEuler
Let’s try Option 2
Edit the inputs.adiabatic-flame file to use the Forward Euler stiff solver.
#_______________________________________________________________________
# Chemistry solver settings:
#
# Define the integrator: ForwardEuler
#mfix.chemistry.integrator = ForwardEuler # default
mfix.chemistry.integrator = StiffSolver::ForwardEuler
#mfix.chemistry.integrator = StiffSolver::BackwardEuler
#mfix.chemistry.integrator = StiffSolver::VODE
Re-run the executable.
./mfix.unit-test.chem3d.gnu.DEBUG.ex inputs.adiabatic-flame
These results look better!
+-------+----------+----------+----------+----------+---------+
| X_CH4 | X_O2 | X_CO2 | X_H2O | X_N2 | T_f |
+-------+----------+----------+----------+----------+---------+
| 0. | 0.037533 | 0.126989 | 0.103965 | 0.731513 | 2066.04 |
+-------+----------+----------+----------+----------+---------+
What output is written to the screen?
- Column 1 is the simulation
time - Column 2 is fluid density
- Columns 3-7 are the five fluid species
- Column 8 is fluid enthalpy
- Column 9 is fluid temperature