The finite-volume form of the CFD equations was given previously as:
aPfP |
= | aNfN + aSfS + aEfE + aWfW + aHfH + aLfL + aTfT + sources |
(1) |
Before solution, this equation is cast into correction form; the required solution is written as:
fP |
= | fP* + fP' |
(2) |
where fP* is the current value in the cell and fP' is the correction that must be applied to it. Re-writing (1) in terms of fP' leaves the term on the right-hand side equal to the imbalance in equation (1) using fP* - this is the residual for the cell:
Saneigh(fP' - fneigh') = imbalance based on fP* |
(3) |
The imbalances in the finite-volume equations are the measure of the quality of the 'solution' at each stage in the iterative process. A total residual for each solved variable is displayed on the screen while a run is in progress; this residual is defined as
quoted residual = S|cell imbalances|/RESREF(f) |
(4) |
RESREF is a normalisation value, which can be user-selected (SELREF=F) or calculated automatically (SELREF=T). More detailed information can be found in the FAQ section of the CHAM web site.
When the quoted residual falls below 1.0, solution for that variable will stop; when all quoted residuals fall below 1.0 the run itself will finish. Care is needed here: an inappropriately high RESREF will give an appearance of convergence when the solution is far from accurate. More advice is provided below on assessing whether convergence has really been achieved.
Also visible on the screen is a display of variable values at a user-specified monitor location in the grid; this location can be changed during the simulation. If convergence has been achieved spot values at all locations must be unchanging.
As an alternative to the spot value, it is possible to monitor the absolute values of the largest corrections anywhere in the domain. Once the largest correction falls to zero, or at least a negligible fraction of the value being corrected, we can assume that convergence has been achieved, even if the sum of the residuals has not fallen below the cut-off.
The display of Maximum Corrections is switched on from the Main Menu by following Output - Monitor Graph Style - Max abs corr. In the Q1 the setting is held as ISG52=2.
There are three main tools that are used to determine whether reasonable convergence has been achieved:
Source balance
Good source balance should give a discrepancy between the positive and negative sums that is a small percentage of either; values of <1% should usually be achieved, and considerably lower figures are not uncommon - but a law of diminishing returns applies, with the cost in CPU time of the extra convergence being rather greater than that required to get to a reasonable level.
Of course, not all variables can be expected to balance: the effect of solid obstructions enters the momentum equations through the pressure gradient, which is a built-in source rather than an externally imposed one - it does not therefore appear in the source balance. Typically, mass (R1/R2) and energy will balance, as will most other scalars, but care is still needed: if the value of a scalar is given a fixed value anywhere, the source required to preserve that value is NOT included in the source print-out. Source imbalance is a clear indication that convergence has not yet been achieved; source balance does not, though, necessarily indicate that convergence has been achieved - residual and spot value behaviour should also be considered.
Residual behaviour
Residuals are the imbalances (or errors) in the equations for each solved variable (see above); interpretation of residuals can be difficult and very subjective. What is certainly true is that residual values should typically go down by at least a factor of 100 from the value after the initial few sweeps (assuming that the calculation is starting from an arbitrary initial state). Eventually, the residuals are likely to level out, with small oscillations about a fairly constant value. This is usually an indication of convergence, but not always: too tight relaxation can sometimes suggest this sort of behaviour because variables are not able to change by much on each sweep, while too loose relaxation can prevent residuals falling further because the variable values are oscillating.
Spot value behaviour
Spot value behaviour is therefore useful in determining whether or not the leveled residuals can be trusted! If the spot values in a representative region of the flow have settled down to a more-or-less constant value, it is reasonable to assume (if the residual behaviour looks promising) that convergence has been achieved; if the changes are still significant, convergence has not been achieved. Some care is still needed though: apparent settling down of spot values might be caused by too-tight relaxation, resulting in a very slow drift that can be mistaken for real convergence.
Maximum Correction Behaviour
Maximum correction behaviour is probably a more reliable guide to convergence behaviour than spot value behaviour. If the largest correction is zero, the solution will not change any more and iteration can stop. If the biggest correction settles to a non-zero value, or is diminishing very slowly, this may well be a sign of too-tight relaxation.
The most widely-known technique for achieving convergence is the use of relaxation. This slows down (relaxes) the changes made to the variables from sweep to sweep. Relaxation does NOT alter the final solution, only the way in which it is achieved. Two types of relaxation are available: linear and false timestep.
Solution of the finite-volume equations generates a field of f values (fnew). Linear relaxation replaces this with
fnew |
= | fold + a(fnew - fold) | (5) |
where fold is the value from the previous sweep; a (between 0 and 1) is the relaxation coefficient. A relaxation coefficient of 0 prevents any change from the previous sweep; a value of 1 applies no relaxation.
False timestep relaxation modifies the finite-volume equations by adding an additional, pseudo-transient term:
(mass in cell)×(fold - fnew)/dtf
where dtf is the false timestep.
A large value of dtf makes the additional term small (light/loose relaxation); a small value makes the additional term large (heavy/tight relaxation).
CONWIZ is the current
default relaxation mechanism for cases set through the VR-Editor. The main features of
CONWIZ are:
It starts by making guesses about 'reference values' of:
The ideas underlying CONWIZ
No 'theory of convergence-promotion' has been invented; but the following ideas have provided guidance:
DVELDP = AREA / [(sum of CONVCOs) + MASS / (DT + DTF) ]
- In steady-flow problems the true time step, DT, equals 0;
- at the start, convection coefficients (CONVCOs) may be 0;
- so a finite false time step, DTF, is needed to keep DVELDP finite; and so to prevent small pressure differences from producing enormous velocities;
- but these are not needed once CONVCOs are finite because motion has started;
- the much-used constant-DTF relaxation is thus now recognised as intrinsically bad; it should no longer be used.
Settings for CONWIZ
The control settings for CONWIZ are set in the VR-Editor from the Main menu - Numerics - Relaxation settings panel. The default maximum increments for velocity are 100m/s per sweep. In some low-velocity cases this may be too large.
The default maximum increase in temperature is 1000° per iteration. This is often too big, especially if there are strong heat sources present. At the start of the calculation there is often no flow, so the only way the heat released can spread is by conduction, leading to impossibly high temperatures. If these are then used to derive properties, such as density, the calculation may not survive.
Under these circumstances, the maximum increment for temperature should be reset to, say, 10° per iteration. This will allow the flow field to develop sufficiently to convect the heat away before unrealistic temperatures are predicted.
The Q1 settings for maximum increment are:
SPEDAT(SET,MAXINC,variable_name,R,maximum_increment)
CONWIZ will force linear relaxation for the velocities. It will use a default value of 0.5. If the user changes these values, the new values will be used.
It will also default the relaxation for scalars to linear, but will allow the user to switch to false timestep if desired.
SARAH (Self-Adjusting Relaxation AlgoritHm) can be used to calculate a false timestep internally; user-input is still available in the form of the value assigned to the SARAH parameter. The false timestep is given by
dtf = SARAH × (internally calculated value).
Typical values of SARAH are in the range 0.1 - 0.001, but lower values may be needed. The value can be changed during run-time from the graphics monitor display.
SARAH cannot be used together with CONWIZ. It only operates on the velocities, but not the scalars.
Additional convergence controls can also be applied, although the default settings are usually adequate.
LITER and ENDIT are used, for each variable, to specify the number of iterations (within the slab or the whole-field) to be applied to the variable and the solution error at which the iteration should terminate automatically. Typical values for LITER are 20 for pressure and 10 for velocity components; 1e-3 is normal for ENDIT.
It is not generally worthwhile increasing LITER or reducing ENDIT, both of which can result in much longer run times, because the internal iterations are part of the outer loop through the domain (sweep); the non-linearities in the equation mean that there is no point trying to achieve a very precise solution on any one sweep.
Occasionally, increasing LITER for pressure can be beneficial in cases for which the pressure field is slow to reach its correct level (perhaps because two parts of the domain are connected by a link in which there is a large pressure drop) even though a realistic velocity field has been achieved.
CONWIZ will default the setting LITER for P1 to 200 to help convergence.
If PARSOL is active, increasing the LITER for the temperature equation (TEM1) to, say, 200 may help to converge the values in the cut cells.
VARMIN and VARMAX are used to specify the maximum permissible values of each variable. The default values effectively impose no constraints (VARMAX of 1e10, 1e6 or 1e4 and VARMIN = -VARMAX) but sometimes convergence can be improved by lower absolute values, usually to control excessive changes in the early stages of the simulation. In some cases these would otherwise prove fatal: extreme temperatures or pressures can have disastrous consequences because of their effect on density.
Two-phase simulations can often benefit from limits on the volume fractions (R1 and R2), to avoid the effective disappearance of one of the phases in parts of the domain. This is true as long as both phases are expected every where. If one of the phases is expected to disappear, then the VARMIN for R1 and R2 should be reduced considerably from the default 1E-6 otherwise stuck continuity errors may occur.
The most important thing to realise is that the final solution must NOT be influenced by the choice of these limits. So, if the solution shows any cells in which the variable has one of the limiting values it cannot be relied upon: the simulation must be repeated (or restarted) using less restrictive limits.
VARMIN and VARMAX are used to specify the maximum permissible values of each variable. The default values effectively impose no constraints (VARMAX of 1e10, 1e6 or 1e4 and VARMIN = -VARMAX)
The k-e turbulence model is widely used and involves significant source terms in the equations for the two turbulence properties. These source terms are linearised to aid convergence, but different linearisations can be chosen to suit the circumstances prevailing in the simulation. The parameter KELIN is used to achieve this. KELIN=0 introduces standard linearisation; KELIN=1 is designed to anticipate the effects of walls (near which turbulence creation and dissipation balance); KELIN=2 controls growth of turbulence in confined flows; KELIN=3 is a strict Newton-Raphson approach for use when the turbulence length scale changes only slowly. The default is KELIN=3, which is usually correct.
This document has provided a summary of the numerical options available within PHOENICS. It is by no means comprehensive and the reader is expected to want more detailed information on topics of particular interest.
These can be found in POLIS (PHOENICS On-Line Information System), which is part of the PHOENICS installation, or on CHAM's web site. In particular, the Encyclopaedia and Documents and Lectures entries are likely to be of value.