5. Choice of Numerical Methods
________________________________________
5.1 Convergence
________________________________________
Q1. Question: - How do I assess Convergence?
Answer:
A number of devices are provided in PHOENICS to monitor and control the convergence of the solution.
RESIDUALS
The residuals are the quantities used in PHOENICS to monitor the convergence procedure. These are the imbalances in the FVE, defined as:
• ep = ap( p) -  (ai( i- p)) - s
•
• i=W,E,S,N,H,L,T
During the computation, the residuals are printed by Earth to the VDU, or are displayed graphically.
The frequency of display update is controlled by the PIL variable TSTSWP, which is set in the VR-editor Main menu / Output / Monitor update frequency panel. This panel also determines whether the data is shown graphically or in ASCII.
The quantities printed out are, in fact,
•  | p| / RESREF ( )
where the sum is extended to the current slab (in PARABolic cases) or to the whole field (in elliptic cases).
When the sum of errors divided by RESREF falls below 1.0, solution for that variable stops. However, the residual continues to be computed for the variable, and should it rise above 1.0, the variable re-enters the solution cycle. When the residuals for all variables fall below 1.0, sweeping stops.
RESREF
EARTH can calculate values for RESREF, based on the net sum of fluxes for a variable, if the PIL variable SELREF=T. Cases set through the VR-Editor have this as a default, but it is not set in all library cases.
The supplementary variable RESFAC acts as a tolerance - a setting of 0.01 means that sweeping stops when the errors are 1% of the reference fluxes.
These two variables can be set from the VR-Editor Main menu / Output / Relaxation control panel.
When SELREF is set to T, the value of RESREF calculated by EARTH is
RESREF() = RESFAC*TOTFLO()/(NX*NY*NZ)
TOTFLO() is the summation of all sources for variable phi, including sources associated with convection and diffusion transport. It contains terms of the form:
SOURCES: SUM1[ SUM2[ C (V - p) ]]
sum1 is over all cells, sum2 is for all sources
C is the coefficient
V the value
 p is the in-cell value of 
CONVECTION: SUMi[ (max(0,(Mdot*I)) + time-flux)]
sumi is for all directions
Mdot is the mass flow rate
time-flux for steady cases is zero, for transient is
(abs(Mdot*)T - abs(Mdot*  )Told)/DT
DIFFUSION: SUMi[ (max(0,(  *A*( i -  p)/ )))]
sumi is for all directions
 is the diffusion coefficient
 the internodal distance
A is the cell-face area
For pressure TOTFLO(P1) is SUMi[ MAX(0,(Mdot))], sumi is for all directions.
RESREF is updated at the end of each sweep, when IZSTEP=NZ. For W1 and W2 resref would be updated at the end of the sweep when IZSTEP=NZ-1.
If SELREF=F, the user is free to set their own values for RESREF for each variable. The default value of 1E-8 is small enough to ensure that sweeping does not stop before LSWEEP. RESREF values can be set from the VR-Editor Main menu / Numerics / Iteration control panel (but only if SELREF=F).
A suitable value for RESREF is often the inflow source for that variable. This can either be calculated in Q1, or obtained from the 'net sources' printout in RESULT from an earlier run.
RESREF(f ) would then be set to tolerance*inflow-flux, where tolerance can be say 0.01. Solution for the variable would stop when the sum of errors is less than 0.01 of the inflow flux of that variable.
The units of RESREF are always those of the equation in question. For the velocities, they are Newtons, for H1 or TEM1 they are Watts. The only exception is P1, where the error printed is a volume error, so RESREF(P1) has units of m^3/s
Q2. Question: - How do I choose my Relaxation settings?
Answer:
The false time step relaxation works by adding the term
 * Vol /  tf * ( p* -  p)
to the equation. p* in the in-store value from the previous sweep, and  p is the value which is about to be calculated.
If the  tf is very big, then the term is negligible. If  tf is very small, then the term becomes dominant and the value cannot change. In between, there is a large range of  tf values, which make the term approximately the same order of magnitude as the other terms in the equation. The term then acts to damp down violent sweep-to-sweep changes, which may lead to divergence.
The question now arises - how to estimate values for  tf?
It is usually best to base it on some characteristic time-scale of the process. One time scale is based on convection - (minimum distance)/(maximum velocity), another is based on diffusion - (minimum distance)2/(kinematic viscosity).
Other time scales can be based on KE/EP or buoyancy. Normally the smallest of these will be used, multiplied by some factor.
A typical starting point is the domain length divided by the inlet velocity divided by the number of cells in that direction - i.e. typical cell residence time.
The values appearing in the Relaxation Control panel of the VR 2D menu, the RELAX command in Q1, and seen from the Reset option of the graphical convergence monitor, are actual  tf. Their size depends on the problem being solved. There are no limits on 'largeness' or 'smallness'.
From Version 2.2.0 onwards, Earth has been able to calculate average time-scales for each variable, much in the same way as it can calculate residual normalising factors if SELREF=T. It does this if the PIL variable SARAH > 0. The actual false time step used for each variable is then:
 tf = SARAH * Internally_calculated_value.
Typical values of SARAH seem to be in the range 0.1 - 0.001, but can be greater or smaller. The value of SARAH can be changed during run-time from the Earth graphics monitor.
False-time step relaxation can be applied to all SOLVEd variables, except P1, R1, R2 and RS. These, and all STOREd variables, can use linear relaxation. Relaxation for pressure is rarely needed, except for highly compressible flows, and most BFC cases.
Linear relaxation is also now recommended for KE and EP, in conjunction with the source-term linearisation KELIN=3. Typical relaxation factors are LINRLX 0.4 for both KE and EP.
At the end of the day, the only correct relaxation is one which allows the case to converge in a reasonable time. Once you have achieved such a value, you will probably waste more time in trying to tune it just that bit better, than you will eventually save.
Q3. Question: - How do I convert the residuals obtained with SELREF=T to residuals normalised by reference inflow quantities?
Answer:
This question is best answered by means of an example, which in this case is taken to be 2d single-phase flow in the y-z plane with solution for mass continuity (P1), momentum (V1 and W1), enthalpy (H1) and mixture fraction (MIXF). Note that in this example the units of 'resref' and 'res sum' for P1 are kg/s rather than m^3/s because DENPCO=T in the Q1 input file.
The conversion procedure is as follows:
Near the bottom of the RESULT file you will find the following normalised residuals printout for the last sweep:
Whole-field residuals before solution
with resref values determined by EARTH
& resfac= 1.000E-03
variable resref (res sum)/resref
P1 3.718E-09 1.562E+00
V1 7.078E-11 1.819E+02
W1 1.880E-08 6.684E-01
H1 8.555E-03 4.698E-01
MIXF 1.317E-10 5.908E-01
In addition, the following nett source printout is given for all boundary conditions and source terms, including the reference inflow rates of each dependent variable:
Nett source of W1 at patch named: INLETF = 1.202E-06
Nett source of W1 at patch named: INLETA = 2.238E-04
Nett source of W1 at patch named: OUTLET =-5.453E-04
Nett source of W1 at patch named: WFNN =-2.016E-04
nett sum=-5.219E-04 pos. sum= 2.250E-04 neg. sum=-7.470E-04
Nett source of R1 at patch named: INLETF = 2.003E-06
Nett source of R1 at patch named: INLETA = 3.799E-05
Nett source of R1 at patch named: OUTLET =-3.999E-05
nett sum=-3.638E-12 pos. sum= 3.999E-05 neg. sum=-3.999E-05
Nett source of H1 at patch named: INLETF = 9.902E+01
Nett source of H1 at patch named: INLETA = 4.615E+01
Nett source of H1 at patch named: OUTLET =-1.452E+02
nett sum= 2.136E-04 pos. sum= 1.452E+02 neg. sum=-1.452E+02
Nett source of MIXF at patch named: INLETF = 2.003E-06
Nett source of MIXF at patch named: OUTLET =-2.003E-06
nett sum=-6.821E-13 pos. sum= 2.003E-06 neg. sum=-2.003E-06
The actual residuals can be recovered by multiplying the (resref) column by the (res sum)/resref column, as follows:
variable resref (res sum)/resref res sum
P1 3.718E-09 1.562E+00 5.8075E-09 kg/s
V1 7.078E-11 1.819E+02 1.2875E-08 N
W1 1.880E-08 6.684E-01 1.2566E-08 N
H1 8.555E-03 4.698E-01 4.0191e-03 W
MIXF 1.317E-10 5.908E-01 7.7808E-11 kg/s
These 'res sum' values can now be normalised by dividing them by the appropriate reference inflow rates (taken from the Nett source printout), as follows:
variable ref value res sum %error
P1 3.999E-05 5.8075E-09 kg/s 0.0145
V1 2.250E-04 1.2875E-08 N 0.0057
W1 2.250E-04 1.2566E-08 N 0.0056
H1 1.452E+02 4.0191e-03 W 0.0028
MIXF 2.003E-06 7.7808E-11 kg/s 0.0039
Q4. Question: - How do I overcome convergence problems resulting from the introduction of buoyancy forces?
Answer:
If convergence proves difficult when activating buoyancy forces, then one or more of the following comments and suggestions may prove helpful.
• Check the problem specification with reference to buoyancy effects, especially in respect of specification of the density formulae, reference density and pressure-based exit conditions
• Employ whole-field solution of the pressure-correction (P1) and energy equation (TEM1 or H1) with a fairly large number of iterations, say 30.
• Apply linear relaxation to the energy equation.
• Try using the conjugate-gradient residuals solver on the P1 equation, i.e. ENDIT(P1)=GRND1, especially if the geometry is complex.
• The use of a laminar flow model is incorrect if the Reynolds or Rayleigh number is sufficiently high, and its use may lead to convergence difficulties.
• One useful trick which often helps to get convergence for very demanding turbulent natural convection simulations is to use two-stage approach as follows:
1. Make the first run with turbulence model, say KEMODL, and higher laminar viscosity of about 0.01 and then
2. Use re-start with actual viscosity to get the final solution.
• Use false-time-step for relaxation of the momentum equations based on the buoyancy time scale. The buoyancy time scale (free convection) can be computed from: Tb = [ (Tr*H)/(DT*g) ]**0.5 /NH where Tr=reference temperature in Kelvin, H=domain height, DT=max temperature difference, i.e. temperature difference which drives the buoyancy, g=gravitational acceleration, and NH=number of cells in H direction.
• For mixed convection problems, the false time step used for the momentum equations can be taken as: Tfalse= minimum(Tc,Tb) where the forced-convection timescale is: Tc = L/U/NL where L=domain length in direction of flow, and U = inlet or typical velocity. and NL=number of cells in the L direction.
• Try heavier relaxation on the momentum equation aligned with the direction of gravity.
• Try time marching to a steady-state solution, especially if the solution appears to be neither converging or diverging. It is possible that the flow is inherently unsteady, as for example when the Rayleigh number is transitional or internal gravity waves are present. Thus, perform a transient solution (STEADY=F) by restarting from a steady-state solution field (PHI file) using as large a real time step as numerically possible. This approach will not need as much momentum relaxation or as many sweeps (per time step) as a steady-state simulation.
Q5. Question: - After the introduction of all relevant model features the simulation diverges dramatically?
Answer:
The advice is that the golden rule is to build a cfd model step-by-step on a coarse mesh rather than introduce all the necessary geometrical and physical features all at once on a fine mesh. The temptation to do the latter is always present because 'the client or management wants it now, but preferably yesterday'.
However, the experience shows a step-by-step approach is usually quicker. The reason is that in the event of convergence problems, the origin of the problems can be detected relatively easily or even identified immediately as the result of the introduction of a new feature.
In the event that the model is already built, then check the mesh is reasonable and even if it is, coarsen the mesh if it is too fine to permit quick computer turnaround. Otherwise, unless divergence is very rapid, you will be waiting forever to see the effects of any changes you make.
Next, experiment with the numerous relaxation features to try and secure convergence.
If this fails, then look at the solution after a few sweeps and examine the flow field very carefully. Sometimes you can identify the growth of the instability from a particular region of the flow.
If all this gets you no where, then you have to start simplifying further, e.g. use a uniform density if it is variable, use a fixed eddy viscosity rather than the k-e model, eliminate buoyancy in mixed-convection problems, use very high interphase friction when doing a two-phase Eulerian simulation, and so on.
The re-start from the converged run is often very helpful when dealing with highly non-linear simulation features.
The main idea is to get something that converges, and then go on from there. |