4.1 Non-Stiff Reactions in Laminar Flow
(a) Species Equations
(b) Energy Equation
4.1.1 Activation
4.2 Activation
4.2.1 Algorithm
4.2.2 Implementation
4.2.3 Activation and Control
If the reactions that take place are expected to be non-stiff, or only weakly stiff, then the reaction source terms may be introduced as source terms in the finite-volume equations and then solved in a point-by-point manner by the PHOENICS solver. The manner in which these source terms are introduced in the PHOENICS species and energy equations are described below.
If the reactions that take place are expected to be non-stiff, or only weakly stiff, then the reaction source terms are included in the following manner.
For species k the reaction source term So ( in g/s) is written as:
So(k) = VOL.W(k).w(k) (4.1)
where VOL is the cell volume (cm^3), W is the molecular mass (g/mol), and w is the net molar production rate (mol/cm^3 s).
Equation (4.1) can be rewritten in terms of the molar creation rate C(k) and the molar destruction rate; the latter is expressed using the characteristic destruction time-scale, tau(k), for species k, leading to:
So(k) = VOL.W(k).(C(k) - [X(k)]/tau(k)) (4.2)
where [X(k)] is the concentration of species k in mole/cm^3, ie.,
[X(k)] = Y(k).Rho/W(k) (4.3)
So, for our CO, and VALue we have
CO = Rho/tau(k) (4.4)
VAL = W(k).C(k)/CO (4.5)
which gives us an approximate linearisation in terms of Y(k).
The heat release term in the enthalpy equation is a little more difficult. What is done is to calculate the heat-release rate when GXCHKI is visited for ther first species so that when the heat-source is required for the TEM1 equation, it is available in F-array storage.
In completely blocked cells for which the volume porosity is no greater than 1.E-10, or for cells in which the value of the PRPS variable indicates a solid, ie. it is not smaller than 100, the CO and VAL terms are not set for the reactions or for the temperature.
In order to make this method work effectively with fairly stiff systems of reactions, it is advisable to use a form of self- adaptive false-timestep under-relaxation. The false-timestep used is the lesser of the EARTH false-timestep, the characteristic creation timescale, and the characteristic destruction timescale. This additional under-relaxation is added directly into the Ap term of the species equations.
For the TEM1 equation it is best to use linear under-relaxation with a factor of around 0.1. In addition if the variable QDOT is stored (to hold the heat release-rate), linear under-relaxation may also be applied to QDOT which helps greatly to control the convergence of the temperature variable.
It should be stressed that the correct relaxation parameters are problem dependent and so little more than broad guide-lines can be given here.
The chemical reaction source terms are activated by setting PATCHes with names for which the first 6 characters are;
CHEMK1
and with a PATCH-type of VOLUME, and COVALs for the CHEMKIN species and, if solved, TEM1. For example,
PATCH(CHEMK1,VOLUME,1,NX,1,NX,1,NY,1,NZ,1,LSTEP) DO K=1,KK-1 + COVAL(CHEMK1,C:K:,GRND9,GRND9) ENDDO COVAL(CHEMK1,TEM1,GRND9,GRND9)
where KK-1 is the number of species solved.
When the reaction source terms dominate the terms in the difference equations for advection and for diffusion, and when the time-scales for the reactions are widely varying, then an implicit method is required which solves all species simultaneously. Ideally one would like to solve all species, and temperature simultaneously in all cells.
Whilst this is possible for 1-dimensional problems, the storage requirements for a 2-dimensional case, let alone a 3-dimensional case, are prohibitive. So, the approach that has been followed in the current CHEMKIN interface, is to solve all species and the enthalpy equation implicitly in a Point-By-Point manner, thus recognising the predominance of the species-to-species coupling over the cell-to-cell coupling.
The algorithm adopted is to solve the equation
Res(k) = Su(k) - Ap(k).(Y(k)-Y(k;0)) + M.w(k).W(k) = 0
for each species, and for enthalpy, to solve,
Res(h') = Su(h') - Ap(h').(h' - h'(0)) -
K
M.Sum(h0(k).W(k).w(k)) = 0
k=1
where Ap and Su include the contributions due to other source terms, and due to the convective and diffusive fluxes, but do not include any contribution from the reactions. In this formula Y(k;0) and h'(0) are the values of mass-fractions and enthalpy on entry to the solution procedure.
Solution of the equations is by means of a Newton-Raphson method, with self-adjustable false-timestep under-relaxation.
The algorithm is implemented as follows:
When cells are blocked the behaviour is modified as follows. If the value of the volume porosity, VPOR, is no greater than 1.E-10, or if the value of PRPS is greater than 197 when EGWF=T, thereby indicating that for a scalar variable there is no transport into the blockage, then no solution is performed for the cell.
If, however, the value of PRPS indicates a solid for which transport of scalars is active, then the mass fractions are not solved, but the temperature, TEM1, if solved, is solved using the standard point-by- point method.
The solution is performed by calling a routine, PHOTWO, that drives the Sandia TWO-PoiNT boundary problem solver (TWOPNT) in two stages, if required. If the energy equation is solved, then a solution is first obtained for the species at constant temperature, only after this has been achieved does the solver, if required, proceed to introduce the energy equation.
In the event that a solution cannot be obtained, TWOPNT is allowed to apply a false time-step under-relaxation and will reduce the time-step by a user-specified factor until either a solution is obtained or a minimum time-step is encountered whereupon a failure is reported.
Once a solution is obtained, TWOPNT will attempt to increase the false time-step, and will eventually eliminate the false time-step under-relaxation altogether. The false time-step under-relaxation introduced by TWOPNT is in addition to any under-relaxation introduced into Ap by PHOENICS itself. For the initial time-step, TWOPNT takes the value of false-timestep set for the first CHEMKIN variable.
The meaning of the residuals is changed to be the greatest absolute residual prior to solution instead of the mean of the absolute residuals. Sweeping, or for PARABolic cases slab iterations, may be halted when the residuals for all CHEMKIN variables fall below the appropriate RESREF values and the change in the variables in a sweep or slab iteration falls below the appropriate ENDIT.
The user is advised to set the species and temperature equation RESREF values <0.0, as PHOENICS is prone to switch off the solution of these equations when using the TWOPNT solver. For NZ> 1, the user may observe that the species and temperature residuals may eventually trap to a value of (NZ+1). This is because the residuals calculated by the CHEMKIN Interface have fallen below this limiting value computed in EARTH. This problem, which does not affect the solution, will be rectified in latter releases.
Positive RESREF values may be used with caution for the TWOPNT solver, but the user is advised to set values of the order of 1.E-6 or even smaller so as to prevent PHOENICS switching off the solver.
The implicit TWOPNT solver is activated by setting
CHSOA = GRND9
The action of the solver may be controlled by making Q1 settings for the various solver parameters by way of the PIL command SPEDAT. This command is used as follows:
SPEDAT(SET,CHEM,VNAME,TYPE,VALUE)
where VNAME is the solver control parameter (see below), VALUE is the value of VNAME and TYPE can be R for real, I for integer, L for logical and C for character. In pre PHOENICS 2.2 releases of CHEMKIN, it should be mentioned that the less-flexibile READ-Q1 facility was used instead of SPEDAT.
The following solver parameters are available:
CALL1 - if 1 then for problems for which the energy (TEM1) equation is solved a preliminary solution with fixed temperature should be performed, if 2 or 0 then only the full solution should be performed
NJAC - maximum no. of iterations of N-R solver before the Jacobian is updated without internal false- timestepping (Default: 20)
ITJAC - maximum no. of iterations of N-R solver before the Jacobian is updated with false-timestepping (Default: 50)
IRETIR - no. of timesteps to be taken before the false- timestep is increased (Default: 50)
NUMDT - the number of time-steps to be taken in the event of a failure of the Newton iteration (Default: 50)
NINIT - no. of false-timesteps taken by the N-R solver prior to attempting a solution without internal false- timestepping (Default: 0)
UFAC - factor by which the time-steps are increased when IRETIR steps have been performed (Default: 2)
DFAC - factor by which the time-steps are decreased when a time-step fails to converge (Default: 4.1)
DTMIN - the lower limit below which the false-timstep is not permitted to fall (Default: 1.E-10)
DTMAX - the upper limit above which the false-timestep is not permitted to rise (Default: 0.1)
ATOL - the absolute tolerance used by TWOPNT to determine whether the steady-state problem has converged according to the formula
dVAR >(ATOL + RTOL.VAR)
(Default: 1.E-4)
RTOL - the relative tolerance used by TWOPNT to determine whether the steady-state problem has converged according to the formula given above (Default: 1.E-8)
ATIM - the absolute tolerance used by TWOPNT to determine whether the false-transient problem has converged according to the formula given above (Default: 1.E-4)
RTIM - the relative tolerance used by TWOPNT to determine whether the false-transient problem has converged according to the formula given above (Default: 1.E-8)
An example on how to set one of the foregoing solver-control parameters is SPEDAT(SET,CHEM,NJAC,I,50).
wbs