Because the solution procedures built into PHOENICS are essentially iterative, it is usual that many adjustment cycles must be carried out before the residuals have been reduced sufficiently for the computation to be terminated. Sometimes however, the residuals in one or more equations start and continue to increase rather than decrease as the sweeps proceed. This phenomenon is known as "divergence"; and it is often accompanied by the appearance of unphysically large values in some of the dependent variables. This process can be watched at the VDU screen. (See PHENC entries UWATCH and TSTSWP).
In hydrodynamic problems, it is often the values of the pressure which depart farthest from physically-realistic values; because of this, PHOENICS has a built-in trap which causes the run to be terminated when excessive values of pressure are encountered. A message explaining this reason for the abandonment is then printed.
Sometimes the problem is not that residuals grow but that they fail to diminish sufficiently as far as has been expected. Often this behaviour is associated with finite-amplitude fluctuations in the values of some of the dependent variables, which continue indefinitely. Such behaviour may be quite localised, as when a pressure fluctuation close to an aperture connected with the outside causes inflows and outflows there, which alternate from sweep to sweep.
However, it must also be recognised that, since all computations are subject to round-off error, there is often a limit beyond which no reduction of the residuals is possible. Such computations can fairly be regarded as fully converged. It is therefore always wise to set SELREF = T and RESFAC = a number which is not unreasonably small (eg 0.01). (See SELREF and RESFAC)
When divergence occurs, it is necessary to establish the cause; this if it is not faulty posing of the problem, is usually to be found in the strength of the linkages between two or more sets of equations, which are being solved in turn rather than simultaneously.
For example, if a problem of free-convection heat transfer is being solved, the two-way interaction between the temperature field and the velocity field, whereby each influences the other, is a common source of divergence.
Whether it is the source in a particular case is easily established by 'freezing' the temperature field before the divergence has progressed too far, and then seeing if the divergence persists. If it does not, the velocity-temperature link can be regarded as the source of the divergence; otherwise, the cause must be sought in some other linkage.
How is the 'freezing' to be effected? The simplest means is to under-relax heavily. If this is done through the SATELLITE, an appropriate setting is:
RELAX(H1,FALSDT,1.E-10) , or alternatively
RELAX(H1,LINRLX,1.E-10) .
In this way, one can investigate the contributions to divergence of linkages between: individual velocity components; the volume fractions and the velocity field in a two-phase process; the chemical-species concentrations in a reacting flow; the pressure and the density in a compressible flow; the turbulence energy and its dissipation rate in a turbulent flow; and many other variables.
If the non-convergence is thought to originate in the linkage between the flow field and a boundary condition, as is often true of the fluctuations mentioned above, it may be beneficial to impose a LOCAL 'freezing'. This can be effected by defining a 'source' patch for the region of interest, and then setting the third argument of the associated COVAL as a large number such as 1.E10, and the fourth argument as SAME.
If freezing by very severe under-relaxation restores convergent behaviour, it is obviously possible that modest under-relaxation will have the same qualitative tendency, while still allowing the solution to proceed so that all residuals do finally diminish (the residuals of 'frozen' variables, incidentally, do not diminish).
The use of under-relaxation is by far the most common means of securing convergence in practice. If employed indiscriminately, it can lead to waste of computer time; however, when it is applied to just those equations that have been identified as potential causes of divergence, and in just the amount that is necessary to procure convergence, it is probably the best first-recourse remedy to apply.
PHOENICS is equipped with a device called EXPERT for adjusting under relaxation factors automatically. It cannot be expected to make the correct choices in extremely complex circumstances; but it is always worth trying. (See PHENC entry EXPERT).
The VARMAX and VARMIN quantities can also be useful in preventing divergence, particularly when this is associated with a poor starting field for the iterative process, which allows some variables at first to stray outside the physically meaningful range.
VARMAX and VARMIN can be employed in two distinct ways: they may set upper and lower bounds to the values of the dependent variables themselves; or they may set limits to the changes in these values in a single adjustment cycle. (See PHENC entry VARMIN)
The values of VARMAX and VARMIN can of course be chosen differently for each dependent and auxiliary variable; but there is no built-in way of applying them over restricted regions at present. However, it is easy for a user who is willing to introduce coding in GROUND to do this for himself; for the functions FN22 and FN23 are available for applying limits to any variable over regions defined locally by IXF, IXL, IYF, etc.
Source terms can be introduced in GROUND either directly or in a linearised manner. The first involves ascribing a single quantity for each cell, viz the source itself; the second involves ascription of two quantities, the "coefficient" and "value". (See PHENC entries: COVAL and BOUNDARY CONDITIONS).
Suppose that the source of variable j, S, does depend upon j, so that it is possible to express it in a linearised manner as:
S = S# + (j-j#)*(dS/dj)#
wherein the # appended to a symbol signifies "current best estimate" of its value.
This expresion can be reformulated in terms of CO and VAL (See COVAL) ie as:
S = C*(V-j)
if C = - (dS/dj)#
and
V = j# + S#/C
If C is positive (ie S decreases as j increases), convergence is promoted by using the linearised-source form for S.
If however C is negative (ie S increases as j increases), use of this form may cause divergence. It is then better to use the first-mentioned mode, ie to employ simply:
S = S#
for otherwise the denominator of the expressions for j may fall through zero, causing unrealistically large values of j to be computed.
Many computer codes for fluid-flow simulation operate wholly in the transient mode, even when it is the ultimately-approached steady state that is of practical interest; and one reason is that, if the transient development of a flow process is followed realistically, the non-physical values which characterise divergence can hardly occur.
The fully-implicit formulations built into PHOENICS permit the steady state to be simulated more directly, albeit in an iterative way which does have some of the properties of 'time-marching'. Nevertheless, if convergence proves hard to procure, there is no reason why the user should not, after setting STEADY to F, watch the flow phenomenon develop in time.
Doing so, as well as yielding indirectly the desired steady-state solution, may well reveal the cause of the difficulty encountered by the direct approach. Thus, one feature of the flow may prove to require very small time steps for its adequate resolution; in the direct approach, small 'false' time steps would correspondingly have been needed.
Of course, if convergence still proves impossible to achieve by way of the time-dependent approach, no matter how small are the time steps which are employed, the making of some fundamental in the problem set-up should be suspected. After all, if the flow that it i desired to simulate is physically possible, there must be some way o simulating it numerically.
A recalcitrant divergence problem is sometimes best solved by changing its formulation. For example, high-Mach-number steady-state compressible flows are difficult to simulate when the density is obtained from the pressure and the temperature, the latter being deduced from the stagnation enthalpy and the kinetic energy.
The reason is that the speeds of approach of the enthalpy and the velocities to their true-solution values may be very different, with the result that densities may be temporarily computed which are very far from the true-solution values; moreover, the "staggering" of the velocity fields with respect to the scalar storage locations introduces uncertainties as to how the kinetic energy to be linked with the stagnation enthalpy is to be deduced by interpolation.
When however an equation is solved for the specific enthalpy, and the density is deduced from this quantity and the pressure, this problem may disappear entirely.
wbs