Two-species relaxation towards equilibrium: "equilibrating"

Reaction model

The reactant is consumed at constant density. At any instant of time, the products will not necessarily be at the same density as the reactant because of their different compressibility. The mass reacted, its energy and the change in volume fraction are taken into account in calculating the change in product density and specific internal energy.

Equilibration method

The equilibration processes for pressure and temperature are applied in operator-split fashion (i.e. at constant total volume) over a finite time interval. A number of substeps is specified; the time interval is divided by this and each equilibration process applied in turn over the substep interval.

Pressure

The pressure difference between the phases is calculated. The bulk sound speed (dp/drho at constant entropy) is used to estimate the change in volume fraction necessary to give pressure equilibrium. The actual change in volume fraction applied is obtained assuming an exponential approach to equilibrium, with some time constant tau. The fraction of the estimated change for equilibrium is given by
f = 1 - exp(-dt / tau)

where dt is the time interval. The specific internal energy of each phase is also adjusted in a forward-time way by applying the p.dv work. This adjustment is made conservative by using the (volume-averaged) mean pressure and density in the cell.

If complete equilibration is desired, the equilibration time scale should be set very small (not zero though). A safety factor is applied to the change in volume fraction for numerical stability. If this is fs say, the number of substeps should be at least some small integer (like 5 or 10) times 1/fs.

Temperature

The temperature difference between the phases is calculated. The specific heat capacity constant volume is used to estimate the heat transfer necessary to give thermal equilibrium. The actual change in specific internal energy applied is obtained assuming an exponential approach to equilibrium, with some time constant tau. The fraction of the estimated change for equilibrium is given by
f = 1 - exp(-dt / tau)

where dt is the time interval.

If complete equilibration is desired, the equilibration time scale should be set very small (not zero though). A safety factor is applied to the change in volume fraction for numerical stability. If this is fs say, the number of substeps should be at least some small integer (like 5 or 10) times 1/fs.

Input

EOS 0
EOS 1
tau_mech tau_therm	# time scales for mechanical and thermal equilibration
p_tol t_tol		# "negligible" difference in p and T between phases
substeps		# number of times to cycle equilibration process
sfac			# relaxation factor
dfvmax			# maximum change in volume fraction per iteration
maxits			# maximum iterations over volume fraction per cycle

Notes