PHOENICSer常问的问题集锦[转贴]
PHOENICS Frequently-Asked Questions; a selectionCHAM Technical Report TR/010
Document last revised 21/08/02
________________________________________
Contents
1. Introduction
2. Installation
3. Boundary Conditions and Sources
4. Physical modelling
4.1 Turbulence
4.2 Multi-phase
4.3 Combustion
4.4 Melting and Solidification
4.5 Chemical Vapour Deposition (CVD)
4.6 Fire and Smoke Modelling
5. Choice of Numerical Methods
5.1 Convergence 1. Introduction
CHAM's User-Support team provides assistance to PHOENICS users world-wide, both by responding to individual enquiries and by publishing, by way of CHAM's website, a Frequently-Asked-Questions (FAQs) page.
Partly in order to draw attention to its existence, and partly because of the direct usefulness of its contents, a selection of questions and answers is now included in the hard-copy documentation of PHOENICS.
A full set of frequently asked questions and answers is available on the Internet on a special website, for maintained customers. In addition, the answers to many questions can often be found in the PHOENICS Encyclopaedia. This can be accessed from within PHOENICS-VR via the HELP and POLIS options. Alternatively, it can be accessed directly via CHAM's website at http://www.cham.co.uk
[ 本帖最后由 linda 于 2007-4-4 05:41 编辑 ] 2. Installation
________________________________________
Q1. Question: - When I click on 'compile' from the VR Environment I encounter the error message 'Bad command'. How do I cure the problem?
Answer:
This problem is usually because the path name to the script dfvars.bat is incorrect in the file /phoenics/d_utils/phoepath.bat.The relevant part of the phoepath.bat ( or phoepath.bat) file looks like:
set DFDir="c:\Program Files\Microsoft Visual Studio\DF98"
rem
echo Adding DIGITAL PHOENICS to path
path=\phoenics\d_utils\d_windf;\phoenics\d_utils;%path%
call %DFDir%\bin\dfvars.bat
When the Compaq Visual Fortran compiler was installed in the user's computer, it may not have been placed in the default directory, but in some other directory or even on a different drive.
The 'set DFDir' line should be modified to reflect the actual location of the compiler. Note that the " " marks around the path name are vital if the directory names contain spaces.
The default location for V6.0 onwards is:
c:\Program Files\Microsoft Visual Studio\DF98.
For V5.0 of the compiler, the default location is:
c:\Program Files\Devstudio\DF\bin\dfvars.bat.
In summary, the solution to the 'Bad Command' problem is to determine the location of the dfvars.bat file on your computer and then correct the path name in the phoepath.bat file.
Q2. Question: - When I click on 'compile' from the VR Environment I encounter the error message 'Out of environment space'. How do I cure the problem?
Answer:
This problem may be cured by inserting the following line in your c:config.sys file:
shell=c:\command.com c:\ /P /E:4096.
If this doesn't work, then try the following. When you click on 'compile' and then on 'GROUND', you can try increasing the environment memory for the window in which the compilation is taking place, i.e. click on the 'properties' of the window, then click on 'memory' and then under 'Initial Environment' change 'Auto' to say '4096'. This should be equivalent to the 'config-file' suggestion made above.
Q3. Question: - When I click on 'compile' from the VR Environment I encounter the error message 'Too many parameters'. How do I cure the problem?
Answer:
This error message may be cured by modifying the file /phoenics/d_utils/phoepath.bat so that the path statement:
echo Adding DIGITAL PHOENICS to path
path = \phoenics\d_utils\d_windf;\phoenics\d_utils;%path%
is replaced by
echo Adding DIGITAL PHOENICS to path
path = \phoenics\d_utils\d_windf;\phoenics\d_utils;c:\dos;c:\windows\command
This means the path is set explicitly, rather than adding to an existing path by the use of %path%. The path statement above is the bare minimum required for Windows - it should replicate the actual path on the users system.
Q4. Question: - I get black holes in my contours in the VR-Viewer. How do I cure the problem?
Answer:
Some customers have experienced some problems with some graphics cards before, including...
• some very small black holes in contours; just the size of a pixel
• large black holes in contours.
In each case, the problems were solved by downloading the latest upgrade of the graphics card driver software from the graphics card manufacturer's website.
Previously reported problems due to graphics cards were...
• MATROX G450 AGP Card. A bug in their driver software caused big black holes in contours. The customer, in Venice, solved the problem by downloading their latest driver software from their website at http://www.matrox.com
• ATI 128 Pro GL Graphics Card, with 16MB. The customer was only able to use Phoenics on 256 colour mode, and not with other options like 'true colours', for which it gave an error message, related to the ati .dll driver software.
Whatever Graphics Card/Video Board you have, it will often pay to ensure that you have the latest upgrades/service packs for their driver software installed, especially if you are experiencing problems. Sometimes their software may have bugs. You can usually download for FREE, their latest driver software from the manufacturer's website. Note that care is needed to download the software for the particular operating system being used.
Q5. Question: - The VR-Viewer sometimes draws the contours in white.
Answer:
On some machines, probably faster machines, in VR-Viewer, when plotting contours, it sometimes plots the 'slice' all in white,i.e. no colour contours can be seen. In the meantime, users can solve the problem only by exiting and restarting VR. Then, when the user next presses the contour button, he/she should not release the button quickly, but rather he/she should hold the button down that little bit longer - say a second or two. This practice seems to avoid the problem. However,once the white contour does appear, the user must exit VR and restart to clear the problem.
Q6. Question: - I have problems selecting objects with the mouse in the VR-Editor.
Answer:
We found that the picking problem, found by one user, was cured by lowering the Hardware Acceleration of the Graphics Card, from the default of Full, via the Display advanced troubleshooting settings.
[ 本帖最后由 linda 于 2007-4-4 05:43 编辑 ] 3. Boundary Conditions and Sources
________________________________________
Q1. Question: - How do I set resistance terms, such as the Ergun Equation?
Answer:
Introduction
The ERGUN equation represents the pressure drop associated with flow through a packed bed, or more generally, a porous medium.
In PHOENICS, pressure drops are applied as sinks in the momentum equation. The units of the momentum equation are those of force, so the pressure drop must also be turned into a force.
The ERGUN equation can be written as:
dp/dx = a*U1 + b*U12
where a=µ/K and b is an empirical constant. Here, µ is the dynamic viscosity of the fluid flowing through the porous media and K is an empirical constant called the permeability.
The pressure drop can be turnedinto a force, as follows:
Fx = Area * dp
•= Area * dp/dx * dx
Now Area = dy * dz;
so dx * Area = Volume
Hence:
Fx = Volume * (a*U1 + b*U12)
Implementation
The velocities used in the equation are superficial velocities, so
This can be introduced as two patches:
•PATCH(ERGUN1, VOLUME, ........
•COVAL(ERGUN1, U1, a, 0.0)
•PATCH(ERGUN2, VOLUME, ........
•COVAL(ERGUN2, U1, -b, 0.0)
Note that the -ve COefficient in patch ERGUN2 is activating a quadratic source - as stated in POLIS.
In the VR-Editor, the above can be achieved by giving an object the type User-Defined, and setting two patches with the Attributes as described above.
Q2. Question: - How do I set a Stagnation Pressure boundary condition?
Answer:
For incompressible flows, the stagnation pressure is:
Pt = P + 0.5 . . vel**2
This can be manipulated to give the mass flux:
. vel = ( 2 . . (Pt - P) )
For P1 (and also P2), a negative coefficient flags a mass source of the form:
Sm = T . ( abs ( Cm . (Vm - Pp)))
where Cm and Vm are the coefficient and value for mass (P1).
The desired stagnation source can thus be set by putting:
coefficient = - 2 *
value = Pt
type = appropriate area
These settings must be supplemented by a COVAL for the velocity component normal to the boundary, with coefficient set to ONLYMS, and value set to SAME.
PIL example:
PATCH (PSTAG, LOW, 1,N,1,NY, 1,1, 1,LSTEP)
COVAL (PSTAG, P1, -2*RHO1, PSTAG)
COVAL (PSTAG, W1, ONLYMS, SAME)
Note that for compressible flows, it is still possible to set stagnation conditions, though not in such an easy fashion.
In the VR-Editor, the above can be achieved by giving an object the type User-Defined, and setting the Attributes as described above 4. Physical Modelling
________________________________________
4.1 Turbulence
Q1. Question: -How do I obtain the wall shear stresses from my PHOENICS simulation?
Answer:
If you put STORE(STRS,SKIN,YPLS) in Group 7 of the Q1 file, PHOENICS will write to the PHI and RESULT files the wall shear stress (tau/rho), the local skin friction coefficient (tau/(rho*urel**2)), and the dimensionless wall-to-node distance (y+=ustar*dist/enul).
If the energy equation is solved, then STORE(HTCO) will cause PHOENICS to write the heat transfer coefficients to the PHI and RESULT files.
Q2. Question: -How do I introduce a log-law inlet profile without recourse to GROUND coding?
Answer:
Library case N110 provides an example of PIL Q1 coding for use with high-Reynolds number turbulence models.
Library case T207 provides an example for low-Reynolds turbulence models, but in the form of INIVAL PATCH rather than on INLET PATCH. 4.2 Multi-phase
Q1. Question: - What are the free surface capabilities in PHOENICS?
Answer:
PHOENICS has at least six features which have been used for the simulation of flows with free surfaces, namely:
(1) Height-of-liquid (HOL) method;
(2) Scalar-equation method (SEM);
(3) IPSA, with high interphase friction and an interface sharpening device;
(4) Moving porosities, with a blocked region adapting to the shape of a liquid surface;
(5) Moving body-fitted coordinates; and
(6) Surface tracking by means of moving particles.
Which method should be chosen depends on the type of application area. CHAM can offer advice if comprehensive information is supplied by the user about the kinds of free-surface problems for which he is seeking solutions. This information should include answers to at least the following questions:
(a) Is the flow steady or unsteady?
(b) Is the flow 1D, 2D or 3D?
(c) Are the fluids gases, liquids or both?
(d) Does mass cross the interface, and, if so, in accordance with what law?
(e) Does the interface remain approximately horizontal or approximately vertical?
(f)Are there any other known limitations, for example - no overturning?
(g) Is gravity (or some other body force) active?
(h) Is one, or are both, fluids, in thermal contact with a solid obstruction to the flow?
(i)Is either fluid compressible?
(j)Is the user interested in a single example or in a class of problems?
(k) If the user is interested in a class, what is the extreme range conditions contained in the class?
Some technical comments on the various methods are:
• SEM is likely to be expensive in computer timeas it needs to be run as a transient, and requires small time steps to combat numerical smearing of the interface (free surface), and to maintain numerical stability.
• HOL is more economical than SEM, but may have to be run as a transient in order to obtain convergence. In addition, this method does not allow any overturning of the interface.
• IPSA involves solving twice as many equations as the other methods unless run in equal-velocity mode, i.e. one set for the liquid phase and one for the gas. The methods will need to be run as a transient if the interface-sharpening device DONACC is used to avoid numerical smearing of the interface. 4.3 Combustion
Q1. Question: - How can I best use PHOENICS for simulating the propagation of flame, and of the associated pressure waves, in a closed cavity, initially filled with a combustible gaseous mixture, when ignition occurs at more than one point? The dimensions and times are small enough for the process to be presumed laminar. The mixture consists of hydrogen, oxygen and nitrogen; and detailed chemical kinetics must be taken into account.
Answer:
Because the thickness of the flame front will be very small compared with the dimensions of the apparatus, an accurate simulation will require the grid to be very fine.
Since the whole phenomenon is three dimensional and transient. Computing times will be so long as (probably) to be regarded as impractical. Therefore some simplifications will need to be made, as follows:
1. The influence of the initial composition, temperature and pressure, and of the chemical-kinetic constants, can best be explored by use of a one-dimensional model; for this will allow the grid to be fine enough for the reaction-diffusion interactions within the thin flame front to be computed accurately.
2. In this case it will be advisable to make use of the stretching-coordinate system exemplified in library case 109; for this eliminates numerical diffusion.
3. The CHEMKIN add-on should be used for the correct representation of the thermodynamic and chemical-kinetic properties of the gases; and if the highest realism is required, the multi-component diffusion module of PHOENICS-CVD can be activated.
4. Such a one-dimensional study will reveal under what conditions detonation (i.e. significant chemical reaction, ahead of the flame, brought about by adiabatic compression) will occur. It will also enable a flame-speed-versus-pressure-and-temperature database to be created for each initial condition.
5. The behaviour of the gas with multiple ignition points can thereafter be computed by treating the flame as a thin discontinuity of temperature and composition which propagates through space at the speed given by the data-base created in the one-dimensional study. Such a multi-dimensional thin-flame study is exemplified by the PhD work of J.Z.Wu.
6. The onset of detonation will be detected, in this multi-dimensional thin-flame study, by accurate computation of the reaction in the ahead-of-flame gas. These reactions, being rather uniform in space, do not require a fine grid. 4.4 Melting and Solidification
Q1. Question: -Should the H1 or TEM1 formulation of the energy equation be used to simulate the melting of a metallic charge in a furnace together with conductive heat transfer through the bounding refractory walls?
Answer:
The use of H1 is not recommended for composite-media problems in Versions earlier than 3.5, because the conductive heat flux was represented by qi =-(k/Cp)*dH/dxi.
The TEM1 equation is recommended because the heat flux is represented correctly as qi =-k*dT/dxi, and it is solved as an energy equation with temperature as dependent variable. Consequently, the convection and transient terms are represented in terms of enthalpy.
The latent heat effects can be introduced into the TEM1 equation by one of two methods:
(a) the use of an effective specific heat.
(b) the use of an enthalpy source term to represent the evolution of latent heat.
Method (a) is simpler to implement than Method (b), but can suffer convergence difficulties around limits of the phase-change region due to the discontinuous nature of the effective specific-heat curve at the solidus and liquidus temperatures. This leads to over-corrections in the energy equation which could be corrected by appropriate modifications to the central-coefficient term Ap in the TEM1 equation. This balance equation, like all others in PHOENICS, is solved in correction form.
The use of method (b) with the TEM1 equation presents no formulation difficulties because, as was mentioned already, the transient and convection terms are represented in terms of the enthalpy. The treatment of the latent-heat source term requires special attention to ensure convergent behaviour, as discussed by Voller and Prakash .
Further information on how to implement either of the two melting methods is not covered by free user support, but simple Q1 examples and GROUND subroutines can be purchased from CHAM's Consultancy Team. 4.5 Chemical Vapour Deposition (CVD)
Q1. Question: - How do I store thermal diffusion coefficients within CVD PHOENICS?
Answer:
The thermal diffusion coefficients can be stored by use of:
STORE(TD01,TD02,.......,TDNN)
Q2. Question: - In CVD PHOENICS what are the units of the polynomial coefficients in the THERM.DAT data file?
Answer:
The polynomial coefficients for the calculation of the thermodynamic properties are dimensionless as is evident from the equations below:
cp/R = a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4
H0/RT = a1 + a2/2*T + a3/3*T^2 + a4/4*T^3 + a5/5*T^4 + a6/T
S0/R = a1*ln(T) + a2*T +a3/2*T^2 + a4/3*T^3 + a5/4*T^4 + a6/5*T^5 + a7
Q3. Question: - How can I use species that contain more than 4 elements, when only 4 can be included in the THERMDAT data file?
Answer:
The format of the THERMDAT file was based on that used by CHEMKIN, which only provided room for 4 elements. To overcome this limitation, PHOENICS-CVD makes use of other spaces in the THERMDAT entry to enable a further 4 elements to be included.
The first 4 elements occupy columns 25-44 on line 1 of an entry: 5 columns per element - 2 for the symbol, 3 for the number of atoms of the element in each species molecule. The additional elements should be positioned as follows:
Element 5- line 1, columns 74-78
Element 6- line 4, columns 61-65
Element 7- line 4, columns 66-70
Element 8- line 4, columns 71-75
The following fictional entry provides an example:
A1B2C3D4E5F6G7H8 70590A1B 2C 3D 4G0300.00 4000.001500.00E5 1
0.01913678E+03 0.08578044E-01-0.08882060E-05-0.03574819E-08 0.06605142E-12 2
-0.06560876E+06-0.08432507E+03-0.04662286E+02 0.06091547E+00-0.04710536E-03 3
0.01968843E-06-0.03563271E-10-0.05665403E+06 0.04525264E+03F 6G 7H 8 4 4.6 Fire and Smoke Modelling.
Q1. Question: - How do I model smoke movement?
Answer:
The PHOENICS scalar variable C1 (or alternatively SMOK in FLAIR) is solved to represent the smoke concentration.
The obvious ways to introduce a fire are:
1. to have a blockage of the fluid material, usually air, with a heat release and associated source of smoke
2. to have an inlet of fluid, carrying with it a smoke concentration and an appropriate temperature (or having a fixed flux of heat instead for any combustion release).
In either case, the units for the C1 variable are to some extent up to the user.The conserved quantity is C1*density, so a natural choice is C1 = mass fraction of smoke, with units of kg/kg, but this is not essential: 'anything'/kg will do, provided that all sources are consistently treated.
The Blockage Object
Many users choose the first modelling approach above, and use a Blockage object with 'Domain material', and the 'Fixed Value' setting to specifythe value of C1 to 1.0 in the fire: C1 values elsewhere therefore give the smoke 'density' (smoke/kg) relative to that in the fire. An alternative is to use the 'Fixed Flux' setting to specify a rate of smoke creation, either for the whole fire or per unit volume.In this way the C1 value can represent a variable for which a creation rate is known (although the absence of a mass source makes this somewhat artificial).Of course, if there are several fires they can have different sources to represent their different behaviour. The C1 value elsewhere in the domain will then represent the cumulative effect of the different contributions. (If necessary, different scalars can be used for the smoke from each fire, say C2, C3 etc: in that way, it is possible to tell how much each one contributes in different locations.)
The Inlet Object
In the same way, with the second approach the value of C1 carried into the domain by thefluid (combustion product) is often set to 1.0; in that case, C1 field values should be interpreted relative to that (possibly arbitrary) value. Alternatively, a more realistic value can be set so that 'mass inflow rate times incoming C1 value' represents the mass flowrate of smoke created by the fire. (Note that the value of C1 in the fire will be rather lower than the incoming value, because the smoke will be diluted by the air in the domain as soon as it appears.)
If C1 is defined as a mass fraction (units of kg of smoke per kg of gas mixture), it is possible to derive the local smoke density (units of kg of smoke per unit volume) by multiplying C1 by the fluid density.In general, this requires storage of the density (DEN1) and another concentration variable (say CSMO); then, in Group 19 Section 6 of GROUND, the smoke density is calculated by:
CALL FN21(LBNAME('CSMO'),C1,DEN1,0.0,1.0)
Alternatively the PLANT or INFORM options of PHOENICS can be used instead of using FORTRAN coding in GROUND. Thus, with INFORM, the following PIL statement in the Q1 file achieves the desired result:
(STORED VAR CSMO IS DEN1*C1)
Fire Object in FLAIR
The FIRE object available in FLAIR has a number of options for the inert modelling of a fire including a pool fire model for transient simulations. There are three other options as discussed below.
1. Mass related heat source
Here, the source of heat and smoke associated with the fire are modelled by an inlet which may be located on top of the fire volume or throughout the fire volume.
The mass release rate from the fire m (in kg/s) appears as a source term in the continuity equation. The total energy flux Q entering the domain from the fire is:
Q = m*(Cp*Tpre+F*Hc)
where Hc is the heat of combustion in J/kg, F is the combustion efficiency factor, Cp is the specific heat of the reactants, and Tpre is the pre-combustion temperature of the reactants, i.e. the temperature of the combustible material before combustion (usually ambient). The total energy flux Q will appear as a source term in the energy equation. The source of smoke S in kg/s that appears in the smoke continuity equation is given by:
S= m*Cin
where Cin is the inlet mass fraction of smoke. The inlet mass fraction of smoke would be specified as the mass production rate of smoke in kg/s divided by the total mass release source in kg/s by the fire (reactants + smoke).
For case I104 (fire in an underground train), one might set the attributes for the fire object volume as follows:
Mass source: Fixed rate = 1.984E-3 kg/s
Heat source: Mass related:
Efficiency = 0.7 Heat =5E7 J/kg
Pre-combustion temperature: 20.0 deg C
Setting scalar: SMOK
Inlet value = 0.1
These settings for the fire volume are essentially equivalent to those used for heat with the "blockage" method, but obviously differ as already explained in respect of the treatment of the smoke. If you examined the RESULT file for such a simulation you would find that:
Nett source of R1 at patch named: FIR13 (FIRE) = 1.984E-03 (mass loss rate from fire in kg/s = continuity source)
Nett source of SMOK at patch named: FIR13 (FIRE) = 1.984E-04 (smoke production rate in kg/s = smoke continuity source)
Nett source of TEM1 at patch named: FIR13 (FIRE) = 7.002E+04 (heat released by fire volume in W = energy source)
The results show a maximum temperature of 136 deg C and a maximum smoke mass fraction of 330 ppm. The smoke mass fraction is not 0.1 at the source because of the diluting effect of the entrained air into the volume of the fire.
2. Fixed flux heat source
This option permits the user to specify a fixed heat release rate at the fire source, i.e. it is the same as the "blockage" method in terms of introducing the heat source. However, it does not include an option to fix the smoke concentration to unity throughout the fire volume.
For case I104 (fire in an underground train), the settings would be as follows:
Mass source: No source
Heat source: Fixed Flux: 7E4 W
Scalar source: Fixed Value: 0.1
No settings are required for the pre-combustion temperature.
3. Fixed temperature
This option is self-explanatory to large extent, i.e. the fire volume is represented by a fixed temperature:
Mass source: No source
Heat source: Fixed Temp: 600 degC
Scalar source: Fixed Value: 0.1 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-09kg/s
V1 7.078E-11 1.819E+02 1.2875E-08N
W1 1.880E-08 6.684E-01 1.2566E-08N
H1 8.555E-03 4.698E-01 4.0191e-03W
MIXF 1.317E-10 5.908E-01 7.7808E-11kg/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-09kg/s0.0145
V1 2.250E-04 1.2875E-08N 0.0057
W1 2.250E-04 1.2566E-08N 0.0056
H1 1.452E+02 4.0191e-03W 0.0028
MIXF 2.003E-06 7.7808E-11kg/s0.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.
页:
[1]