CHAM FAQ
From CFD-Wiki
(→Does Phoenics have built-in turnkey models for complex flow/stress applications?) |
Phoenicsfaq (Talk | contribs) (→How can I set an anisotropic thermal conductivity?) |
||
(91 intermediate revisions not shown) | |||
Line 1: | Line 1: | ||
- | This FAQ is | + | This article contains FAQs for the PHOENICS CFD code. The structure of this FAQ page is still evolving according to the number and type of questions that are frequently asked. |
- | == | + | == PHOENICS == |
- | === What | + | === General === |
+ | ==== What is PHOENICS? ==== | ||
- | + | PHOENICS is a general-purpose commercial CFD code which is applicable to steady or unsteady, one- , two- or three-dimensional turbulent or laminar, multi-phase, compressible or incompressible flows using Cartesian, cylindrical-polar or curvilinear coordinates. The code also has a spatial marching integration option to handle parabolic and hyperbolic flows, as well as transonic free jets in the absence of recirculation zones. | |
- | === | + | The numerical procedure is of the finite-volume type in which the original partial differential equations are converted into algebraic finite-volume equations with the aid of discretisation assumptions for the transient, convection, diffusion and source terms. For this purpose, the solution domain is subdivided into a number of control volumes on a mono-block mesh using a conventional staggered-grid approach. All field variables except velocities are stored at the grid nodes, while the velocities themselves are stored at staggered cell-face locations which lie between the nodes. |
+ | |||
+ | For curvilinear coordinates, the default option is to solve the momentum equations in terms of the covariant projections of the velocity vector into the local grid directions on a structured, mono-block, staggered mesh. An alternative exists to employ the GCV method for curvilinear coordinates, whereby cell-centred Cartesian velocities or covariant velocity projections are solved on a structured mono-block or multi-block mesh. The link description for the latter can handle arbitrary block rotations, and therefore can accommodate an unstructured multi-block grid made from relatively simple, structured blocks. | ||
+ | |||
+ | For complex geometries, PHOENICS by default actually uses a Cartesian cut-cell method named PARSOL which provides an automatic, efficient, and flexible alternative to traditional boundary-fitted grid methods using curvilinear coordinates. The Cartesian cut-cell approach uses a background Cartesian or cylindrical-polar grid for the majority of the flow domain with special treatments being applied to cells which are cut by solid bodies, thus retaining a boundary-conforming grid. Specifically, the method computes the fractional areas and volumes, and employs a collection of special algorithms for computing interfacial areas, evaluating wall shear stresses, and for computing advection and diffusion near solid boundaries, etc. | ||
+ | |||
+ | The finite-volume equations for each variable are derived by integrating the partial differential equations over each control volume. Fully implicit backward differencing is employed for the transient terms, and central differencing is used for the diffusion terms. The convection terms are discretised using hybrid differencing in which the convective terms are approximated by central differences if the cell face Peclet number is less than 2 and otherwise by upwind differencing. At faces where the upwind scheme is used, physical diffusion is omitted altogether. In addition to the upwind and hybrid differencing schemes, PHOENICS is furnished with an extensive set of higher-order convection schemes, which comprise five linear schemes and twelve non-linear schemes. The linear schemes include CDS, QUICK, linear upwind and cubic upwind. The non-linear schemes employ a flux limiter to secure boundedness. These schemes include SMART, H-QUICK, UMIST, SUPERBEE, MINMOD, OSPRE, MUSCL and van-Leer harmonic. | ||
+ | |||
+ | The integration procedure results in a coupled set of algebraic finite-volume equations which express the value of a variable at a grid node in terms of the values at neighbouring grid points and the nodal value at the old time level. For unsteady flows, PHOENICS by default solves each of these equations by an implicit method, but the option exists to revert to an explicit method which is stable only when the Courant number is less than or equal to unity. The explicit method uses old-time neighbour values, whereas the implicit method uses values at the new time level. Although implicit methods allow a much larger Courant number than the explicit methods, it is not unconditionally stable since the non-linearities in the equations often limit numerical stability. | ||
+ | |||
+ | The finite-volume equations are solved iteratively using the SIMPLEST and IPSA algorithms of Spalding, which are embodied in PHOENICS for the solution of single-phase and two-phase flows, respectively. These algorithms are segregated solution methods which employ pressure-velocity coupling to enforce mass conservation by solving a pressure-correction equation and making corrections to the pressure and velocity fields. Multi-phase flows are accommodated using either an Eulerian-Lagrangian method using particle tracking, or an algebraic-slip model. The latter solves mixture continuity and momentum equations, and a volume-fraction equation is solved for each dispersed phase. Algebraic relationships are used for slip velocities relative to the mixture. | ||
+ | |||
+ | The default calculation procedure is organised in a slab-by-slab manner in which all dependent variables are solved in turn at the current slab before attention moves to the next higher slab. The slabs are thus visited in turn, from the lowermost to the uppermost, and a complete series of slab visits is referred to as a sweep through the solution domain. For parabolic and hyperbolic calculations, only one such sweep is required, with many iteration cycles at each slab for parabolic cases, and no outflow boundary condition is required because this is an outcome of the solution. For elliptic calculations, many such sweeps are conducted until convergence is attained at the current time level; in addition, the pressure-correction equation is solved in a simultaneous whole-field manner at the end of each sweep. Thereafter the solution proceeds to the next time level where the iterative process is repeated. The option exists to solve each finite-volume equation in a whole-field manner, and this is actually the default when using the automatic convergence control. The default linear equation solver for each finite-volume equation is a modified form of Stone's strongly implicit solver, but the option exists to use a Conjugate Gradient Residuals solver for each equation. | ||
+ | |||
+ | The GCV method solves the system of algebraic finite-volume equations by using a conjugate-residual linear solver with LU preconditioning, and it uses a segregated pressure-based solver strategy with an additional correction of the cell-centred momentum velocities. This provides faster convergence than conventional one-step face-velocity corrections. | ||
+ | |||
+ | The numerical solution procedure requires appropriate relaxation of the flow variables in order to procure convergence. Two types of relaxation are employed, namely inertial and linear. The former is normally applied to the velocity variables, whereas the latter is applied to all other flow variables, as and when necessary. | ||
+ | |||
+ | The PHOENICS Parallel solver subdivides the solution domain into sub-domains (spatial domain decomposition) and assigns each sub-domain to one processor on a multi-core computer. Each sub-domain has storage overlap ("halo" cells) for data exchange between the various sub-domains. The CFD code then runs simultaneously on all the processors, on its own set of data, and data is exchanged between the various sub-domains ( via the "halo" cells ) using MPI (message-passing interface). A point solver is used for whole-field solution of each of the finite-volume equations. | ||
+ | |||
+ | The convergence requirement is that for each set of finite-volume equations the sum of the absolute residual sources over the whole solution domain is less than one percent of references quantities based on the total inflow of the variable in question. An additional requirement is that the values of monitored dependent variables at a selected location do not change by more than 0.1 percent between successive iteration cycles. It is also possible to monitor the absolute values of the largest corrections to each variable anywhere in the domain. Once the largest correction falls to zero, or at least a negligible fraction of the value being corrected, then it can be assumed that convergence has been achieved, even if the sum of the residuals has not fallen below the cut-off. | ||
+ | |||
+ | === Problem Set Up === | ||
+ | ==== Meshing ==== | ||
+ | ==== Geometry ==== | ||
+ | ==== Physical properties ==== | ||
+ | ===== How can I set variable physical properties for the fluid ===== | ||
+ | |||
+ | It is best to set variable physical properties for a fluid directly in the Q1 input file by means of InForm, as in the following example which sets the properties of liquid sodium as a function of temperature: | ||
+ | |||
+ | In Group 7 of the Q1 input file: | ||
+ | |||
+ | :''''' save7begin''''' | ||
+ | :''REAL(TCRIT); TCRIT=2503.7'' | ||
+ | :''REAL(ACP,BCP,CCP,DCP)'' | ||
+ | :''ACP=1.6582; BCP=-8.479E-4; CCP=4.4541E-7; DCP=-2992.6'' | ||
+ | :''(stored of TDEK is TEM1+:TEMP0:!ZSLSTR)'' | ||
+ | :''(stored of TRAT is TDEK/:TCRIT:!ZSLSTR)'' | ||
+ | |||
+ | :'' ** density'' | ||
+ | :''(stored of SDEN is 219.0+275.32*(1.0-TRAT)+511.58*(1-TRAT)^0.5!ZSLSTR)'' | ||
+ | :'' ** thermal conductivity'' | ||
+ | :''(stored of SOTK is 124.67-0.11381*TDEK+5.5226E-5*TDEK^2-1.1842E-8*TDEK^3!ZSLSTR)'' | ||
+ | :'' ** thermodynamic specific heat'' | ||
+ | :''(stored of SOCP is 1.E3*(ACP+BCP*TDEK+CCP*TDEK^2+DCP/(TDEK)^2)!ZSLSTR)'' | ||
+ | :'' ** mean specific heat defined by integral cp*dt/T between 0 and T'' | ||
+ | :''(stored of SCPM is 1.E3*(ACP+0.5*BCP*TDEK+0.3333*CCP*TDEK^2-2.0*DCP/(TDEK)^2)!ZSLSTR)'' | ||
+ | :'' ** dynamic viscosity'' | ||
+ | :''(stored of SODV is EXP(-6.4406-0.3958*LOGE(TDEK)+556.835/TDEK)!ZSLSTR)'' | ||
+ | :'' ** kinematic viscosity'' | ||
+ | :''(stored of SOKV is SODV/SDEN!ZSLSTR)'' | ||
+ | :''''' save7end''''' | ||
+ | |||
+ | In Group 9 of the Q1 input file: | ||
+ | |||
+ | :''''' save9begin''''' | ||
+ | :'' ** use in-form set set variable physical properties'' | ||
+ | :''(property RHO1 is SDEN)'' | ||
+ | :''(property ENUL is SOKV)'' | ||
+ | :'' ** - sign denotes tehrmal conductivity is set'' | ||
+ | :''(property PRNDTL(TEM1) is -SOTK)'' | ||
+ | :'' ** set mean specific heat for use in TEM1 energy equation'' | ||
+ | :''(property CP1 is SCPM)'' | ||
+ | :''''' save9end''''' | ||
+ | ===== How can I set an anisotropic thermal conductivity? ===== | ||
+ | |||
+ | If you want to set an anisotropic thermal conductivity, this can be done using the In-Form "MODDIF" statement, which re-specifies the diffusion terms. For a case where the default thermal conductivity applies to the x-direction, the required Inform statements have the form | ||
+ | |||
+ | ''(MODDIF of <variable> at <object name> is <factor>*DIFF with NORTH)'' | ||
+ | |||
+ | ''(MODDIF of <variable> at <object name> is <factor>*DIFF with HIGH)'' | ||
+ | |||
+ | Here <variable> is the PHOENICS variable for which you want to modify the diffusion coefficients; <factor> is the numerical multiplicative factor; "DIFF" represents the unmodified diffusion term; and "with NORTH" and "with HIGH means that you are multiplying the y and the z diffusion terms, respectively. | ||
+ | |||
+ | So all you have to do is to specify the variable, which in this case is TEM1, and replace <factor> by the appropriate numerical amount, which will be the directional conductivity divided by the default conductivity (in this case the x-direction conductivity). | ||
+ | |||
+ | You can in fact generalise these statements to | ||
+ | |||
+ | ''(MODDIF of <variable> at <object name> is <any_formula>*DIFF with NORTH)'' | ||
+ | |||
+ | ''(MODDIF of <variable> at <object name> is <any_formula>*DIFF with HIGH)'' | ||
+ | |||
+ | where <any_formula> is any formula satisfying the In-Form rules. | ||
+ | |||
+ | An alternative method is to use "Group 12 patches" to modify the diffusion coefficients, as described here: [http://www.cham.co.uk/phoenics/d_polis/d_enc/patch.htm#2b] | ||
+ | |||
+ | In the PHOENICS VR Pre-processor, a non-isotropic thermal conductivity can be set by means of the PCB object, as explained here: [http://www.cham.co.uk/phoenics/d_polis/d_docs/tr326/obj-type.htm#PCB]. The PCB object instructs PHOENICS EARTH to employ a Group 12 PATCH to allow you to set the ratio between the thermal conductivity in the thickness of the PCB (shortest dimension), and in the in-plane directions (the other two dimensions). This creates a non-isotropic thermal conductivity. The VR-Editor will assign the default material conductivity to the in-plane directions of the PCB object, and the minimum-size direction of the PCB object will be flagged as the through-plane direction. If the PCB object is rotated in the Cartesian mesh to be a inclined PCB, then the (through thickness)/in-plane conductivity ratio will be incorrect because no account is taken of the inclination. | ||
+ | |||
+ | === Physical Modelling === | ||
+ | ==== General ==== | ||
+ | ==== Boundary Conditions ==== | ||
+ | ==== Turbulence Modelling ==== | ||
+ | |||
+ | === Numerical === | ||
+ | |||
+ | ===== How do I know if the solution is converged? ===== | ||
+ | There are four main tools that are used to determine whether reasonable convergence has been achieved: | ||
+ | |||
+ | * source imbalances | ||
+ | * sums of absolute residual errors | ||
+ | * spot value behaviour | ||
+ | * maximum absolute corrections | ||
+ | |||
+ | |||
+ | ''Source imbalances'' | ||
+ | |||
+ | 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 - the behaviour of the residuals, spot-value and variable maximum corrections should also be considered. | ||
+ | |||
+ | ''Sums of absolute residual errors'' | ||
+ | |||
+ | Residuals are the imbalances (or errors) in the equations for each solved variable, summed over the cells within the computational domain. These can either be absolute (SELREF=F) or relative (SELREF=T); in the latter case the imbalances are scaled with respect to a value that is intended | ||
+ | to represent the flux of the variable throughout the domain.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 levelled 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 absolute corrections'' | ||
+ | |||
+ | Apart from residual errors, 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, it would be reasonable to assume that convergence has been achieved, even if the sum of the residuals has not fallen below any specific level. You can create a plot of the maximum absolute corrections from VR, by selecting: Menu / Output / Monitor graph style / Max abs cor. This can be a very useful measure of convergence. | ||
+ | |||
+ | ===== My case is not converging. How can I secure convergence? ===== | ||
+ | ==== Relaxation ==== | ||
+ | ==== Numerical Schemes==== | ||
+ | ===== What numerical convection schemes are available in PHOENICS? ===== | ||
+ | ''For the default staggered BFC, Cylindrical-Polar and Cartesian-mesh options:'' | ||
+ | |||
+ | The default numerical convection scheme is HYBRID, and there are 19 numerical convection schemes available: | ||
+ | |||
+ | * UPWIND scheme | ||
+ | * HYBRID scheme, which is a variant of the UPWIND scheme. | ||
+ | * LUS: Linear-upwind scheme. | ||
+ | * FROMM: FROMM's scheme | ||
+ | * CUS: Cubic-upwind scheme | ||
+ | * QUICK: Quadratic upwind scheme. | ||
+ | * CDS: Central-differencing scheme. | ||
+ | * SMART: Bounded QUICK. | ||
+ | * KOREN: Bounded Cubic, TVD (Total-Variation Diminishing ) | ||
+ | * VANL1: van Leer MUSCL, bounded FROMM, TVD | ||
+ | * HQUICK: Harmonic QUICK, smooth | ||
+ | * OSPRE: Based on FROMM, smooth/continuous | ||
+ | * VANL2: Van Leer Harmonic, smooth, TVD | ||
+ | * VANALB: Based on FROMM, smooth/continuous | ||
+ | * MINMOD: Zhu-Rodi SOUCUP scheme, TVD | ||
+ | * SUPBEE: Superbee scheme, compressive limiter for sharp gradients, TVD | ||
+ | * UMIST: Bounded QUICK, TVD | ||
+ | * HCUS: Harmonic CUS, smooth | ||
+ | * CHARM: Based on QUICK, polynomial ratio. | ||
+ | |||
+ | |||
+ | The Q1 input commands to set the convection scheme are: | ||
+ | |||
+ | * UPWIND scheme for all solved variables: DIFCUT = 0 | ||
+ | * HYBRID scheme for all solved variables: DIFCUT = 0.5 or nothing, since HYBRID is the default scheme. | ||
+ | * Remaining schemes: | ||
+ | : DIFCUT = 0 | ||
+ | : SCHEME(<scheme_name>, ALL) | ||
+ | where <scheme_name> is one of the 17 remaining schemes and the keyword ALL means that <scheme_name> will be applied to all SOLVEd-for variables. | ||
+ | |||
+ | : The command to apply <scheme_name> to just the selected SOLVEd-for variables is: SCHEME(<scheme_name>, <variable_name1>, <variable_name2>,…etc.). The UPWIND scheme will be used for the remaining solved variables. | ||
+ | |||
+ | For example, the following commands set the QUICK scheme for U1 and V1, the SMART scheme for C1 and the UPWIND scheme for the remaining SOLVEd-for variables. | ||
+ | |||
+ | :DIFCUT=0 | ||
+ | :SCHEME(QUICK,U1,V1) | ||
+ | :SCHEME(SMART,C1) | ||
+ | |||
+ | Note that the HYBRID scheme can't be applied to selected variables, it can only be applied to all of them or to none. | ||
+ | |||
+ | You can find more information about all these schemes in [http://www.cham.co.uk/phoenics/d_polis/d_enc/enc_schm.htm ] and [[File:CHAM_Numerical_Schemes_Report_1995.pdf]] | ||
+ | |||
+ | ''For the BFC-GCV option:'' | ||
+ | |||
+ | The default numerical convection scheme is UPWIND. The other schemes available for GCV are: | ||
+ | |||
+ | *CDS: Central-differencing scheme. | ||
+ | *SOUP: Linear-upwind scheme. | ||
+ | *QUICK: Quadratic upwind scheme. | ||
+ | *SMART: Bounded QUICK. | ||
+ | *MINMOD: Zhu-Rodi SOUCUP scheme, TVD | ||
+ | *MUSCL: van Leer MUSCL, bounded FROMM, TVD | ||
+ | *OSHER: | ||
+ | *SUPERB: Superbee scheme, compressive limiter for sharp gradients, TVD | ||
+ | |||
+ | The Q1 commands to set the convection scheme for the GCV option are: | ||
+ | |||
+ | :SPEDAT(SET,GCVSCH,<variable_name>,C,<scheme_name>) | ||
+ | |||
+ | where <variable_name> is one of the SOLVEd-for variables and <scheme_name> is one of the 8 schemes above.For example, the following commands set the MINMOD scheme for KE and EP and leave the default UPWIND scheme for the remaining SOLVEd-for variables. | ||
+ | |||
+ | :SPEDAT(SET,GCVSCH,KE,C,MINMOD) | ||
+ | :SPEDAT(SET,GCVSCH,EP,C,MINMOD) | ||
+ | |||
+ | Please note that the SCHEME command is not compatible with GCV. | ||
+ | |||
+ | === User Programming === | ||
+ | ==== PIL ==== | ||
+ | ==== In-Form ==== | ||
+ | ===== What is In-Form? ===== | ||
+ | In-Form is a supplement to the PHOENICS Input Language (PIL) which facilitates the input | ||
+ | of problem-defining data. Specifically, it allows users to express their requirements through | ||
+ | algebraic formulae, whether for: | ||
+ | |||
+ | * space and time discretization, | ||
+ | * material properties, | ||
+ | * initial values, | ||
+ | * sources, | ||
+ | * boundary conditions, | ||
+ | * body shapes and motions, or | ||
+ | * special print-out features. | ||
+ | |||
+ | The formulae are placed in the data-input file, Q1 by way of a text editor, or via the graphical | ||
+ | user interface of PHOENICS, the Virtual-Reality Editor. | ||
+ | |||
+ | The PHOENICS input module, Satellite, thereafter transmits the implications of the formulae | ||
+ | to the solver module, EARTH. There they are recognised as instructions to perform the | ||
+ | appropriate computations. The user has nothing more to do except await and inspect the | ||
+ | computed results. | ||
+ | |||
+ | Extensive documentation of the In-Form facility is provided here: [http://www.cham.co.uk/phoenics/d_polis/d_enc/in-form.htm], and here: [http://www.cham.co.uk/documentation/tr003.pdf]. | ||
+ | |||
+ | ===== How can I obtain plots of a field variable vs time for a transient simulation?===== | ||
+ | The In-Form facility can be used to output the values of particular variables at multiple locations as a function of time, as exemplified by the following INFORM command in Group 21 of the Q1 input file: | ||
+ | |||
+ | :'' save21begin'' | ||
+ | ''(table in monplt.csv is get(TEM1[5,17,6],TEM1[15,17,6],TEM1[25,17,6],TEM1[35,17,6]) with head(T5,T15,T25,T35)!time)'' | ||
+ | :'' save21end'' | ||
+ | |||
+ | This command writes to a file named "monplt.csv" (MS Excel format) at each time step, the field values of the temperature TEM1 at the grid locations (5,17,6), (15,17,6),(25,17,6) and (35,17,6). Alternatively, you can specify the (x,y,z) cartesian coordinates of the location where you want to tabulate the computed field values. For this, you use curly brackets {} instead of straight brackets []. So, for example, to tabulate TEM1 values at (2.1,0.8,9.6), you would use the following INFORM commands: | ||
+ | |||
+ | :'' save21begin'' | ||
+ | ''(table in montem.csv is get(TEM1{2.1,0.8,9.6}) with head(TEM1)!time)'' | ||
+ | : ''save21end'' | ||
+ | |||
+ | There is an In-Form tutorial explaining the method here: [http://www.cham.co.uk/phoenics/d_polis/d_wkshp/inform/wsinf4.htm]. | ||
+ | |||
+ | ===== How can I specify the thermal properties of a material?===== | ||
+ | |||
+ | You can use In-Form to specify the thermal properties of a material, say soil designated by Material 117 in the Q1 input file, by putting the following commands in Group 9 of the Q1 file: | ||
+ | |||
+ | :'' save9begin'' | ||
+ | : ** Mean Specific heat for use in TEM1 energy equation" | ||
+ | ''(property CP1 is {value} with IMAT=117)'' | ||
+ | :'' ** Density'' | ||
+ | ''(property RHO1 is {value} with IMAT=117)'' | ||
+ | :'' ** Thermal conductivity'' | ||
+ | ''(property PRNDTL(TEM1) is -{value} with IMAT=117)'' | ||
+ | :'' save9end'' | ||
+ | |||
+ | |||
+ | The "save..." commands should start in column 3, and they are delimiters that instruct the pre-processor to read the In-Form commands. The setting of a negative value for PRNDTL(TEM1) denotes that the thermal conductivity is set, whereas as positive value sets the Prandtl number. | ||
+ | |||
+ | ==== GROUND ==== | ||
+ | |||
+ | === Solver Methodology === | ||
+ | ==== PARSOL ==== | ||
+ | ===== What is PARSOL? ===== | ||
+ | PARSOL stands for PARtial SOLids. It is a Cartesian cut-cell CFD method for solving flow in and around complex geometries. The method provides an automatic, efficient, and flexible alternative to traditional boundary-fitted grid methods using curvilinear coordinates for the representation of complex geometries. The Cartesian cut-cell approach uses a background Cartesian or cylindrical-polar grid for the majority of the flow domain with special treatments being applied to cells which are cut by solid bodies, thus retaining a boundary-conforming grid. Specifically, the method computes the fractional areas and volumes, and employs a collection of special algorithms for computing interfacial areas, evaluating wall shear stresses, and for computing advection and diffusion near solid boundaries, etc. | ||
+ | ==== S-PARSOL ==== | ||
+ | ===== What is S-PARSOL? ===== | ||
+ | S-PARSOL stands for Structured PARSOL, and it considers cut links rather than the cut edges of cells like PARSOL, a link being the line joining two cell-center nodes ( these links can be seen by clicking on the "S-PARSOL grid-line dialogue" button in the VR Environment ). S-PARSOL therefore creates no cut-cell control volumes or corresponding equations; but rather it makes changes to the finite-volume coefficients to represent the interactions between the fluid and the solid materials. | ||
+ | |||
+ | S-PARSOL has the substantial advantage over PARSOL of being able to detect thin solids without the need for high mesh densities. Also, S-PARSOL more accurately represents conjugate heat transfer through thin solids and composite media. S-PARSOL's practice of setting to zero the velocities in cut links has the advantage over PARSOL of recognising the presence of thin bodies, but this practice sometimes creates flow patterns of unpleasing appearance, which is why an improved version called X-PARSOL is being developed which is intended to combine the benefits of the earlier techniques without their deficiencies. A fairly detailed technical description of S-PARSOL is given in Section 3.3 of the following paper: | ||
+ | |||
+ | {{reference-paper|author=Spalding, D.B.|year=2013|title=Trends, Tricks & Try-ons in CFD/CHT"|rest=In Advances in Heat Transfer, Vol. 45, pp. 1-78, Editors: Ephraim M. Sparrow, Young I. Cho, John P. Abraham & John M. Gorman, Academic Press, Elsevier Inc}} A pdf copy of this paper can be accessed here: [http://www.cham.co.uk/DOCS/papers/Advances_in_Heat_Transfer_Vol_45_%20Article_00001_Electronic_Offprint_2013.pdf S-PARSOL paper] | ||
+ | |||
+ | ==== Parallel ==== | ||
+ | ==== Double Precision ==== | ||
+ | ===== How can I run in Double Precision? ===== | ||
+ | |||
+ | The default solver ("earexe.exe") file is single precision apart from the evaluation of the time fluxes, and so in this sense the default CFD solver is an hybrid module. | ||
+ | |||
+ | You can run the double-precision solver ("eardble.exe") in "c:\phoenics\d_earth\d_windf" by clicking on "Options" / "Select private solver", and then browsing to this file. The double precision version will then be run whenever you click "Run">"Solver" from the VR Environment. | ||
+ | |||
+ | ===== When do I need to run in Double Precision? ===== | ||
+ | The Double-Precision (DP) solver was first introduced to resolve convergence issues experienced by a number of users for the case of a transient, large-scale fire growing with time under conditions of pure natural convection. Specifically, incomplete convergence was observed at each time step, in the sense that regardless of the number of sweeps performed, the monitored pressure value continued to increase steadily with sweep. The investigation carried out by CHAM technical support indicated that significant rounding error was present when computing the difference between the old and new time fluxes in the transient terms of the continuity, smoke continuity and energy equations. These terms are computed by subtracting two large numbers of similar magnitude, producing a value much smaller than either number, and comparable with the round-off error of single-precision arithmetic, i.e: typically 7 significant digits. The use of the DP solver, which has an accuracy of about 15 significant digits, cured this convergence issue. In later releases of PHOENICS, the default single-precision solver executable (earexe.exe) always uses double precision to compute the time fluxes. | ||
+ | |||
+ | No systematic evaluation of the DP solver has been carried out, and so experience of its use remains rather limited. It does seem that using DP does not significantly affect execution times, and the solutions produced by a sample of cases from the PHOENICS test battery & flow library show no significant differences when run in single and double precision. This is to be expected, since round-off error is generally insignificant in single-precision arithmetic, although one can envisage cases where round-off errors might accumulate after hundreds or thousands of time steps and thereby affect the numerical solution. | ||
+ | |||
+ | Examples where one might expect double precision to make a difference are: | ||
+ | |||
+ | * transient natural-convection cases, especially if using a fully-compressible form of the equation of state; | ||
+ | * cases with a large number of fixed-pressure boundaries ("openings"), with quadratic resistance; | ||
+ | * flows which are inherently unstable, such as for example natural-convection flows, transitional flows, etc; | ||
+ | * meshes with a large difference between the largest and smallest mesh sizes, as for example encountered in applications using low-Reynolds-number turbulence models where the grid can be highly concentrated near walls, and so the differences between neighbouring grid nodes may be very small; | ||
+ | * large geometries with small but significant features; | ||
+ | * flows with large pressure and/or enthalpy values. | ||
+ | * when processing geometrical input data in absolute coordinates. | ||
+ | |||
+ | There is however a downside to using double precision; the memory requirements are double those of the single-precision (SP) solver. | ||
+ | |||
+ | The advantage of using the SP solver as the default is that larger simulations can be performed within a given amount of memory, with perhaps twice as many mesh cells as in double precision. The advantage of the DP solver is that one does not have to worry about the potentially adverse effect of round-off errors. It will however incur double the memory usage. | ||
+ | |||
+ | === Pre-Processing === | ||
+ | ==== VR Editor ==== | ||
+ | |||
+ | === Post-Processing === | ||
+ | ==== VR Viewer ==== | ||
+ | ==== AUTOPLOT ==== | ||
+ | ==== PHOTON ==== | ||
+ | |||
+ | === Installation === | ||
+ | ==== How can I download the latest version of PHOENICS? ==== | ||
+ | You can download the latest delivery version of PHOENICS 2015 (dated 31 March 2016) from CHAM’s FTP site, in either 32bit or 64bit form, by following the download instructions, a copy of which can be provided on request by emailing support@cham.co.uk with details of your name and organisation. | ||
+ | === Application Examples === | ||
+ | ==== Combustion ==== | ||
+ | |||
+ | == PHOENICS-FLAIR == | ||
+ | === General === | ||
+ | ==== What is PHOENICS-FLAIR? ==== | ||
+ | FLAIR is a special-purpose program for Heating, Ventilation and Air Conditioning | ||
+ | (HVAC) systems that are required to deliver thermal comfort, health and safety, air | ||
+ | quality, and contamination control. FLAIR provides designers with a powerful and | ||
+ | easy-to-use tool which can be used for the prediction of airflow patterns, temperature | ||
+ | distributions, and smoke movement in buildings and other enclosed spaces. | ||
+ | |||
+ | |||
+ | Extensive documentation of PHOENICS-FLAIR is provided here: [http://www.cham.co.uk/phoenics/d_polis/d_docs/tr313/tr313.htm], and here: [http://www.cham.co.uk/documentation/tr313.pdf]. | ||
+ | === Comfort Indices === | ||
+ | ==== How can I evaluate outdoor thermal comfort? ==== | ||
+ | There are a range of indices in the literature which can used for evaluating outdoor thermal comfort. Each of these indices accounts for most if not all of the following effects: air temperature and humidity, local wind speed, radiation, activity and clothing. | ||
+ | |||
+ | FLAIR provides two thermal comfort indices that have been used to assess outdoor comfort, namely the Predicted Mean Vote and the Apparent Temperature. These can be activated from the VR Environment by clicking on " Menu">"Model">"Comfort Indices - Settings", and then toggling to "ON" the required comfort indices. | ||
+ | |||
+ | The Thermal Sensation Index is gaining popularity, but it has yet to be implemented in FLAIR. However, you can use the In-Form to compute this index from the Q1 input file, as explained below, but first a discussion will be provided for the comfort indices available in FLAIR. | ||
+ | |||
+ | '''1. Predicted Mean Vote (PMV)''' | ||
+ | |||
+ | The PMV is a comfort index defined in ISO 7730 that predicts the mean value of the votes of a large group of people on a 7-point thermal sensation scale. This index was developed to provide a measure of indoor comfort, but it has been frequently applied in both original and modified forms for outdoor applications (see for example Walls et al [2015] at | ||
+ | |||
+ | http://anzasca.net/wp-content/uploads/2015/12/107_Walls_Parker_Walliss_ASA2015.pdf). The FLAIR option for PMV is described here: http://www.cham.co.uk/phoenics/d_polis/d_chron/ph361/chap5.htm#chap5.2. | ||
+ | |||
+ | '''2. Apparent Temperature''' | ||
+ | |||
+ | The Apparent Temperature is a general term for the perceived outdoor temperature, caused by the combined effects of air temperature, relative humidity, radiation and wind speed. The formulae for the Apparent Temperature used in FLAIR are those specified by the Australian Bureau of Meteorology. They are an approximation of the value provided by a mathematical model of heat balance in the human body. They can include the effects of temperature, humidity, wind-speed and radiation. | ||
+ | |||
+ | Two forms of Apparent Temperature given by the Australian Bureau have been added to PHOENICS-FLAIR based on a mathematical model of an adult walking outdoors in the shade, one with radiation and the other without. The expression without radiation is suitable for the working condition of people walking in the shade; the expression with radiation is suitable for the working condition of people in direct sunlight. | ||
+ | |||
+ | The non-radiative apparent temperature (TAPP) is computed from: | ||
+ | |||
+ | : AT = Ta + 0.33*e − 0.70*ws − 4.0 (1) | ||
+ | |||
+ | whereas the apparent temperature with radiation (TAPR) is calculated from: | ||
+ | |||
+ | : AT = Ta + 0.348*e − 0.70*ws + 0.70*Q/(ws + 10) − 4.25 (2) | ||
+ | |||
+ | where Ta is the Dry bulb temperature (°C) (TEM1 in FLAIR), e is the water-vapour partial pressure (hPa) (PVAP in FLAIR), ws is the wind speed (m/s) at 10m elevation (taken to be the local VABS in FLAIR), and Q is the net radiation absorbed per body surface unit area (W/m2), which can be expected to be in the range of -20 to 130 W/m2. | ||
+ | |||
+ | The vapour pressure e is calculated from the temperature and relative humidity using the equation: | ||
+ | |||
+ | : e = (rh/100)*6.105*exp(17.27*Ta/(237.7 + Ta ) ) (3) | ||
+ | |||
+ | where rh is the Relative Humidity [%] (RELH in FLAIR). | ||
+ | |||
+ | Note that in order to calculate either form of the Apparent Temperature, the solution of Specific Humidity and the derivation of Relative Humidity must both be turned on in FLAIR. | ||
+ | |||
+ | '''3. Thermal Sensation Index (TSI)''' | ||
+ | |||
+ | The TSI determines a measure between 0 and 7, with 4 as the most comfortable condition. The original formula was developed in Japan and included the effects of relative humidity, but was analysed further by Givoni et al [2003] to provide the simplified equation quoted in your email: | ||
+ | |||
+ | : TSI = 1.2 + 0.1115*Ta + 0.0019*SR – 0.3185*ws (4) | ||
+ | |||
+ | where Ta is the air temperature (°C), SR is the horizontal solar radiation flux (W/m2) and ws is the wind speed (m/s). Further information is provided by Givoni et al (2003) in https://www.researchgate.net/publication/223626731_Outdoor_comfort_research_issues. | ||
+ | |||
+ | SR can be interpreted as the total amount of solar radiation received from above by a surface horizontal to the ground. This will include both direct and diffuse solar radiation, where the former is solar radiation that comes in a straight line from the direction of the sun at its current position in the sky. Diffuse radiation is solar radiation that has been scattered by particles in the atmosphere, and comes equally from all directions. | ||
+ | |||
+ | Library Case 674 provides a Wind and Solar Heating example which can be used as a vehicle for computing the TSI by means of the following PIL and In-Form statements. | ||
+ | |||
+ | :'' save7begin'' | ||
+ | :''REAL(SRH,DIRSRF,DIFSRF,SOLALT,PI);PI=3.141593'' | ||
+ | :''SOLALT = 19.85976 ! Solar altitude in degrees'' | ||
+ | :''DIRSRF = 136.0 ! Direct solar radiation flux in W/m^2'' | ||
+ | :''DIFSRF = 446.0 ! Diffuse solar radiation flux in W/m^2'' | ||
+ | :''SRH = DIFSRF + DIRSRF*SIN(SOLALT*PI/180.) ! Horizontal solar radiation flux'' | ||
+ | :''(STORED of TSI is (1.2 + 0.1115*TEM1 + 0.0019*:SRH: - 0.3185*VABS) with IMAT<100)'' | ||
+ | :'' save7end'''' | ||
+ | |||
+ | |||
+ | This Library Case uses a Weather File for Gatwick in the UK ( see http://www.cham.co.uk/phoenics/d_polis/d_wkshp/windsun/windsun2.htm ), and so the foregoing solar data are taken from the PHOENICS RESULT file. | ||
+ | |||
+ | === Application Examples === | ||
+ | ==== Pollutant-gas dispersion in the atmosphere ==== | ||
+ | This Q1 input file simulates steady, 3d pollutant-gas dispersion from a roof-top vent on an isolated building in a neutral wind environment [[File:GasDispersion.pdf]] | ||
- | |||
- | |||
- | |||
[[Category: FAQ's]] | [[Category: FAQ's]] | ||
{{Stub}} | {{Stub}} |
Latest revision as of 16:21, 18 October 2016
This article contains FAQs for the PHOENICS CFD code. The structure of this FAQ page is still evolving according to the number and type of questions that are frequently asked.
Contents |
PHOENICS
General
What is PHOENICS?
PHOENICS is a general-purpose commercial CFD code which is applicable to steady or unsteady, one- , two- or three-dimensional turbulent or laminar, multi-phase, compressible or incompressible flows using Cartesian, cylindrical-polar or curvilinear coordinates. The code also has a spatial marching integration option to handle parabolic and hyperbolic flows, as well as transonic free jets in the absence of recirculation zones.
The numerical procedure is of the finite-volume type in which the original partial differential equations are converted into algebraic finite-volume equations with the aid of discretisation assumptions for the transient, convection, diffusion and source terms. For this purpose, the solution domain is subdivided into a number of control volumes on a mono-block mesh using a conventional staggered-grid approach. All field variables except velocities are stored at the grid nodes, while the velocities themselves are stored at staggered cell-face locations which lie between the nodes.
For curvilinear coordinates, the default option is to solve the momentum equations in terms of the covariant projections of the velocity vector into the local grid directions on a structured, mono-block, staggered mesh. An alternative exists to employ the GCV method for curvilinear coordinates, whereby cell-centred Cartesian velocities or covariant velocity projections are solved on a structured mono-block or multi-block mesh. The link description for the latter can handle arbitrary block rotations, and therefore can accommodate an unstructured multi-block grid made from relatively simple, structured blocks.
For complex geometries, PHOENICS by default actually uses a Cartesian cut-cell method named PARSOL which provides an automatic, efficient, and flexible alternative to traditional boundary-fitted grid methods using curvilinear coordinates. The Cartesian cut-cell approach uses a background Cartesian or cylindrical-polar grid for the majority of the flow domain with special treatments being applied to cells which are cut by solid bodies, thus retaining a boundary-conforming grid. Specifically, the method computes the fractional areas and volumes, and employs a collection of special algorithms for computing interfacial areas, evaluating wall shear stresses, and for computing advection and diffusion near solid boundaries, etc.
The finite-volume equations for each variable are derived by integrating the partial differential equations over each control volume. Fully implicit backward differencing is employed for the transient terms, and central differencing is used for the diffusion terms. The convection terms are discretised using hybrid differencing in which the convective terms are approximated by central differences if the cell face Peclet number is less than 2 and otherwise by upwind differencing. At faces where the upwind scheme is used, physical diffusion is omitted altogether. In addition to the upwind and hybrid differencing schemes, PHOENICS is furnished with an extensive set of higher-order convection schemes, which comprise five linear schemes and twelve non-linear schemes. The linear schemes include CDS, QUICK, linear upwind and cubic upwind. The non-linear schemes employ a flux limiter to secure boundedness. These schemes include SMART, H-QUICK, UMIST, SUPERBEE, MINMOD, OSPRE, MUSCL and van-Leer harmonic.
The integration procedure results in a coupled set of algebraic finite-volume equations which express the value of a variable at a grid node in terms of the values at neighbouring grid points and the nodal value at the old time level. For unsteady flows, PHOENICS by default solves each of these equations by an implicit method, but the option exists to revert to an explicit method which is stable only when the Courant number is less than or equal to unity. The explicit method uses old-time neighbour values, whereas the implicit method uses values at the new time level. Although implicit methods allow a much larger Courant number than the explicit methods, it is not unconditionally stable since the non-linearities in the equations often limit numerical stability.
The finite-volume equations are solved iteratively using the SIMPLEST and IPSA algorithms of Spalding, which are embodied in PHOENICS for the solution of single-phase and two-phase flows, respectively. These algorithms are segregated solution methods which employ pressure-velocity coupling to enforce mass conservation by solving a pressure-correction equation and making corrections to the pressure and velocity fields. Multi-phase flows are accommodated using either an Eulerian-Lagrangian method using particle tracking, or an algebraic-slip model. The latter solves mixture continuity and momentum equations, and a volume-fraction equation is solved for each dispersed phase. Algebraic relationships are used for slip velocities relative to the mixture.
The default calculation procedure is organised in a slab-by-slab manner in which all dependent variables are solved in turn at the current slab before attention moves to the next higher slab. The slabs are thus visited in turn, from the lowermost to the uppermost, and a complete series of slab visits is referred to as a sweep through the solution domain. For parabolic and hyperbolic calculations, only one such sweep is required, with many iteration cycles at each slab for parabolic cases, and no outflow boundary condition is required because this is an outcome of the solution. For elliptic calculations, many such sweeps are conducted until convergence is attained at the current time level; in addition, the pressure-correction equation is solved in a simultaneous whole-field manner at the end of each sweep. Thereafter the solution proceeds to the next time level where the iterative process is repeated. The option exists to solve each finite-volume equation in a whole-field manner, and this is actually the default when using the automatic convergence control. The default linear equation solver for each finite-volume equation is a modified form of Stone's strongly implicit solver, but the option exists to use a Conjugate Gradient Residuals solver for each equation.
The GCV method solves the system of algebraic finite-volume equations by using a conjugate-residual linear solver with LU preconditioning, and it uses a segregated pressure-based solver strategy with an additional correction of the cell-centred momentum velocities. This provides faster convergence than conventional one-step face-velocity corrections.
The numerical solution procedure requires appropriate relaxation of the flow variables in order to procure convergence. Two types of relaxation are employed, namely inertial and linear. The former is normally applied to the velocity variables, whereas the latter is applied to all other flow variables, as and when necessary.
The PHOENICS Parallel solver subdivides the solution domain into sub-domains (spatial domain decomposition) and assigns each sub-domain to one processor on a multi-core computer. Each sub-domain has storage overlap ("halo" cells) for data exchange between the various sub-domains. The CFD code then runs simultaneously on all the processors, on its own set of data, and data is exchanged between the various sub-domains ( via the "halo" cells ) using MPI (message-passing interface). A point solver is used for whole-field solution of each of the finite-volume equations.
The convergence requirement is that for each set of finite-volume equations the sum of the absolute residual sources over the whole solution domain is less than one percent of references quantities based on the total inflow of the variable in question. An additional requirement is that the values of monitored dependent variables at a selected location do not change by more than 0.1 percent between successive iteration cycles. It is also possible to monitor the absolute values of the largest corrections to each variable anywhere in the domain. Once the largest correction falls to zero, or at least a negligible fraction of the value being corrected, then it can be assumed that convergence has been achieved, even if the sum of the residuals has not fallen below the cut-off.
Problem Set Up
Meshing
Geometry
Physical properties
How can I set variable physical properties for the fluid
It is best to set variable physical properties for a fluid directly in the Q1 input file by means of InForm, as in the following example which sets the properties of liquid sodium as a function of temperature:
In Group 7 of the Q1 input file:
- save7begin
- REAL(TCRIT); TCRIT=2503.7
- REAL(ACP,BCP,CCP,DCP)
- ACP=1.6582; BCP=-8.479E-4; CCP=4.4541E-7; DCP=-2992.6
- (stored of TDEK is TEM1+:TEMP0:!ZSLSTR)
- (stored of TRAT is TDEK/:TCRIT:!ZSLSTR)
- ** density
- (stored of SDEN is 219.0+275.32*(1.0-TRAT)+511.58*(1-TRAT)^0.5!ZSLSTR)
- ** thermal conductivity
- (stored of SOTK is 124.67-0.11381*TDEK+5.5226E-5*TDEK^2-1.1842E-8*TDEK^3!ZSLSTR)
- ** thermodynamic specific heat
- (stored of SOCP is 1.E3*(ACP+BCP*TDEK+CCP*TDEK^2+DCP/(TDEK)^2)!ZSLSTR)
- ** mean specific heat defined by integral cp*dt/T between 0 and T
- (stored of SCPM is 1.E3*(ACP+0.5*BCP*TDEK+0.3333*CCP*TDEK^2-2.0*DCP/(TDEK)^2)!ZSLSTR)
- ** dynamic viscosity
- (stored of SODV is EXP(-6.4406-0.3958*LOGE(TDEK)+556.835/TDEK)!ZSLSTR)
- ** kinematic viscosity
- (stored of SOKV is SODV/SDEN!ZSLSTR)
- save7end
In Group 9 of the Q1 input file:
- save9begin
- ** use in-form set set variable physical properties
- (property RHO1 is SDEN)
- (property ENUL is SOKV)
- ** - sign denotes tehrmal conductivity is set
- (property PRNDTL(TEM1) is -SOTK)
- ** set mean specific heat for use in TEM1 energy equation
- (property CP1 is SCPM)
- save9end
How can I set an anisotropic thermal conductivity?
If you want to set an anisotropic thermal conductivity, this can be done using the In-Form "MODDIF" statement, which re-specifies the diffusion terms. For a case where the default thermal conductivity applies to the x-direction, the required Inform statements have the form
(MODDIF of <variable> at <object name> is <factor>*DIFF with NORTH)
(MODDIF of <variable> at <object name> is <factor>*DIFF with HIGH)
Here <variable> is the PHOENICS variable for which you want to modify the diffusion coefficients; <factor> is the numerical multiplicative factor; "DIFF" represents the unmodified diffusion term; and "with NORTH" and "with HIGH means that you are multiplying the y and the z diffusion terms, respectively.
So all you have to do is to specify the variable, which in this case is TEM1, and replace <factor> by the appropriate numerical amount, which will be the directional conductivity divided by the default conductivity (in this case the x-direction conductivity).
You can in fact generalise these statements to
(MODDIF of <variable> at <object name> is <any_formula>*DIFF with NORTH)
(MODDIF of <variable> at <object name> is <any_formula>*DIFF with HIGH)
where <any_formula> is any formula satisfying the In-Form rules.
An alternative method is to use "Group 12 patches" to modify the diffusion coefficients, as described here: [1]
In the PHOENICS VR Pre-processor, a non-isotropic thermal conductivity can be set by means of the PCB object, as explained here: [2]. The PCB object instructs PHOENICS EARTH to employ a Group 12 PATCH to allow you to set the ratio between the thermal conductivity in the thickness of the PCB (shortest dimension), and in the in-plane directions (the other two dimensions). This creates a non-isotropic thermal conductivity. The VR-Editor will assign the default material conductivity to the in-plane directions of the PCB object, and the minimum-size direction of the PCB object will be flagged as the through-plane direction. If the PCB object is rotated in the Cartesian mesh to be a inclined PCB, then the (through thickness)/in-plane conductivity ratio will be incorrect because no account is taken of the inclination.
Physical Modelling
General
Boundary Conditions
Turbulence Modelling
Numerical
How do I know if the solution is converged?
There are four main tools that are used to determine whether reasonable convergence has been achieved:
- source imbalances
- sums of absolute residual errors
- spot value behaviour
- maximum absolute corrections
Source imbalances
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 - the behaviour of the residuals, spot-value and variable maximum corrections should also be considered.
Sums of absolute residual errors
Residuals are the imbalances (or errors) in the equations for each solved variable, summed over the cells within the computational domain. These can either be absolute (SELREF=F) or relative (SELREF=T); in the latter case the imbalances are scaled with respect to a value that is intended to represent the flux of the variable throughout the domain.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 levelled 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 absolute corrections
Apart from residual errors, 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, it would be reasonable to assume that convergence has been achieved, even if the sum of the residuals has not fallen below any specific level. You can create a plot of the maximum absolute corrections from VR, by selecting: Menu / Output / Monitor graph style / Max abs cor. This can be a very useful measure of convergence.
My case is not converging. How can I secure convergence?
Relaxation
Numerical Schemes
What numerical convection schemes are available in PHOENICS?
For the default staggered BFC, Cylindrical-Polar and Cartesian-mesh options:
The default numerical convection scheme is HYBRID, and there are 19 numerical convection schemes available:
- UPWIND scheme
- HYBRID scheme, which is a variant of the UPWIND scheme.
- LUS: Linear-upwind scheme.
- FROMM: FROMM's scheme
- CUS: Cubic-upwind scheme
- QUICK: Quadratic upwind scheme.
- CDS: Central-differencing scheme.
- SMART: Bounded QUICK.
- KOREN: Bounded Cubic, TVD (Total-Variation Diminishing )
- VANL1: van Leer MUSCL, bounded FROMM, TVD
- HQUICK: Harmonic QUICK, smooth
- OSPRE: Based on FROMM, smooth/continuous
- VANL2: Van Leer Harmonic, smooth, TVD
- VANALB: Based on FROMM, smooth/continuous
- MINMOD: Zhu-Rodi SOUCUP scheme, TVD
- SUPBEE: Superbee scheme, compressive limiter for sharp gradients, TVD
- UMIST: Bounded QUICK, TVD
- HCUS: Harmonic CUS, smooth
- CHARM: Based on QUICK, polynomial ratio.
The Q1 input commands to set the convection scheme are:
- UPWIND scheme for all solved variables: DIFCUT = 0
- HYBRID scheme for all solved variables: DIFCUT = 0.5 or nothing, since HYBRID is the default scheme.
- Remaining schemes:
- DIFCUT = 0
- SCHEME(<scheme_name>, ALL)
where <scheme_name> is one of the 17 remaining schemes and the keyword ALL means that <scheme_name> will be applied to all SOLVEd-for variables.
- The command to apply <scheme_name> to just the selected SOLVEd-for variables is: SCHEME(<scheme_name>, <variable_name1>, <variable_name2>,…etc.). The UPWIND scheme will be used for the remaining solved variables.
For example, the following commands set the QUICK scheme for U1 and V1, the SMART scheme for C1 and the UPWIND scheme for the remaining SOLVEd-for variables.
- DIFCUT=0
- SCHEME(QUICK,U1,V1)
- SCHEME(SMART,C1)
Note that the HYBRID scheme can't be applied to selected variables, it can only be applied to all of them or to none.
You can find more information about all these schemes in [3] and File:CHAM Numerical Schemes Report 1995.pdf
For the BFC-GCV option:
The default numerical convection scheme is UPWIND. The other schemes available for GCV are:
- CDS: Central-differencing scheme.
- SOUP: Linear-upwind scheme.
- QUICK: Quadratic upwind scheme.
- SMART: Bounded QUICK.
- MINMOD: Zhu-Rodi SOUCUP scheme, TVD
- MUSCL: van Leer MUSCL, bounded FROMM, TVD
- OSHER:
- SUPERB: Superbee scheme, compressive limiter for sharp gradients, TVD
The Q1 commands to set the convection scheme for the GCV option are:
- SPEDAT(SET,GCVSCH,<variable_name>,C,<scheme_name>)
where <variable_name> is one of the SOLVEd-for variables and <scheme_name> is one of the 8 schemes above.For example, the following commands set the MINMOD scheme for KE and EP and leave the default UPWIND scheme for the remaining SOLVEd-for variables.
- SPEDAT(SET,GCVSCH,KE,C,MINMOD)
- SPEDAT(SET,GCVSCH,EP,C,MINMOD)
Please note that the SCHEME command is not compatible with GCV.
User Programming
PIL
In-Form
What is In-Form?
In-Form is a supplement to the PHOENICS Input Language (PIL) which facilitates the input of problem-defining data. Specifically, it allows users to express their requirements through algebraic formulae, whether for:
- space and time discretization,
- material properties,
- initial values,
- sources,
- boundary conditions,
- body shapes and motions, or
- special print-out features.
The formulae are placed in the data-input file, Q1 by way of a text editor, or via the graphical user interface of PHOENICS, the Virtual-Reality Editor.
The PHOENICS input module, Satellite, thereafter transmits the implications of the formulae to the solver module, EARTH. There they are recognised as instructions to perform the appropriate computations. The user has nothing more to do except await and inspect the computed results.
Extensive documentation of the In-Form facility is provided here: [4], and here: [5].
How can I obtain plots of a field variable vs time for a transient simulation?
The In-Form facility can be used to output the values of particular variables at multiple locations as a function of time, as exemplified by the following INFORM command in Group 21 of the Q1 input file:
- save21begin
(table in monplt.csv is get(TEM1[5,17,6],TEM1[15,17,6],TEM1[25,17,6],TEM1[35,17,6]) with head(T5,T15,T25,T35)!time)
- save21end
This command writes to a file named "monplt.csv" (MS Excel format) at each time step, the field values of the temperature TEM1 at the grid locations (5,17,6), (15,17,6),(25,17,6) and (35,17,6). Alternatively, you can specify the (x,y,z) cartesian coordinates of the location where you want to tabulate the computed field values. For this, you use curly brackets {} instead of straight brackets []. So, for example, to tabulate TEM1 values at (2.1,0.8,9.6), you would use the following INFORM commands:
- save21begin
(table in montem.csv is get(TEM1{2.1,0.8,9.6}) with head(TEM1)!time)
- save21end
There is an In-Form tutorial explaining the method here: [6].
How can I specify the thermal properties of a material?
You can use In-Form to specify the thermal properties of a material, say soil designated by Material 117 in the Q1 input file, by putting the following commands in Group 9 of the Q1 file:
- save9begin
- ** Mean Specific heat for use in TEM1 energy equation"
(property CP1 is {value} with IMAT=117)
- ** Density
(property RHO1 is {value} with IMAT=117)
- ** Thermal conductivity
(property PRNDTL(TEM1) is -{value} with IMAT=117)
- save9end
The "save..." commands should start in column 3, and they are delimiters that instruct the pre-processor to read the In-Form commands. The setting of a negative value for PRNDTL(TEM1) denotes that the thermal conductivity is set, whereas as positive value sets the Prandtl number.
GROUND
Solver Methodology
PARSOL
What is PARSOL?
PARSOL stands for PARtial SOLids. It is a Cartesian cut-cell CFD method for solving flow in and around complex geometries. The method provides an automatic, efficient, and flexible alternative to traditional boundary-fitted grid methods using curvilinear coordinates for the representation of complex geometries. The Cartesian cut-cell approach uses a background Cartesian or cylindrical-polar grid for the majority of the flow domain with special treatments being applied to cells which are cut by solid bodies, thus retaining a boundary-conforming grid. Specifically, the method computes the fractional areas and volumes, and employs a collection of special algorithms for computing interfacial areas, evaluating wall shear stresses, and for computing advection and diffusion near solid boundaries, etc.
S-PARSOL
What is S-PARSOL?
S-PARSOL stands for Structured PARSOL, and it considers cut links rather than the cut edges of cells like PARSOL, a link being the line joining two cell-center nodes ( these links can be seen by clicking on the "S-PARSOL grid-line dialogue" button in the VR Environment ). S-PARSOL therefore creates no cut-cell control volumes or corresponding equations; but rather it makes changes to the finite-volume coefficients to represent the interactions between the fluid and the solid materials.
S-PARSOL has the substantial advantage over PARSOL of being able to detect thin solids without the need for high mesh densities. Also, S-PARSOL more accurately represents conjugate heat transfer through thin solids and composite media. S-PARSOL's practice of setting to zero the velocities in cut links has the advantage over PARSOL of recognising the presence of thin bodies, but this practice sometimes creates flow patterns of unpleasing appearance, which is why an improved version called X-PARSOL is being developed which is intended to combine the benefits of the earlier techniques without their deficiencies. A fairly detailed technical description of S-PARSOL is given in Section 3.3 of the following paper:
Spalding, D.B. (2013), "Trends, Tricks & Try-ons in CFD/CHT"", In Advances in Heat Transfer, Vol. 45, pp. 1-78, Editors: Ephraim M. Sparrow, Young I. Cho, John P. Abraham & John M. Gorman, Academic Press, Elsevier Inc. A pdf copy of this paper can be accessed here: S-PARSOL paper
Parallel
Double Precision
How can I run in Double Precision?
The default solver ("earexe.exe") file is single precision apart from the evaluation of the time fluxes, and so in this sense the default CFD solver is an hybrid module.
You can run the double-precision solver ("eardble.exe") in "c:\phoenics\d_earth\d_windf" by clicking on "Options" / "Select private solver", and then browsing to this file. The double precision version will then be run whenever you click "Run">"Solver" from the VR Environment.
When do I need to run in Double Precision?
The Double-Precision (DP) solver was first introduced to resolve convergence issues experienced by a number of users for the case of a transient, large-scale fire growing with time under conditions of pure natural convection. Specifically, incomplete convergence was observed at each time step, in the sense that regardless of the number of sweeps performed, the monitored pressure value continued to increase steadily with sweep. The investigation carried out by CHAM technical support indicated that significant rounding error was present when computing the difference between the old and new time fluxes in the transient terms of the continuity, smoke continuity and energy equations. These terms are computed by subtracting two large numbers of similar magnitude, producing a value much smaller than either number, and comparable with the round-off error of single-precision arithmetic, i.e: typically 7 significant digits. The use of the DP solver, which has an accuracy of about 15 significant digits, cured this convergence issue. In later releases of PHOENICS, the default single-precision solver executable (earexe.exe) always uses double precision to compute the time fluxes.
No systematic evaluation of the DP solver has been carried out, and so experience of its use remains rather limited. It does seem that using DP does not significantly affect execution times, and the solutions produced by a sample of cases from the PHOENICS test battery & flow library show no significant differences when run in single and double precision. This is to be expected, since round-off error is generally insignificant in single-precision arithmetic, although one can envisage cases where round-off errors might accumulate after hundreds or thousands of time steps and thereby affect the numerical solution.
Examples where one might expect double precision to make a difference are:
- transient natural-convection cases, especially if using a fully-compressible form of the equation of state;
- cases with a large number of fixed-pressure boundaries ("openings"), with quadratic resistance;
- flows which are inherently unstable, such as for example natural-convection flows, transitional flows, etc;
- meshes with a large difference between the largest and smallest mesh sizes, as for example encountered in applications using low-Reynolds-number turbulence models where the grid can be highly concentrated near walls, and so the differences between neighbouring grid nodes may be very small;
- large geometries with small but significant features;
- flows with large pressure and/or enthalpy values.
- when processing geometrical input data in absolute coordinates.
There is however a downside to using double precision; the memory requirements are double those of the single-precision (SP) solver.
The advantage of using the SP solver as the default is that larger simulations can be performed within a given amount of memory, with perhaps twice as many mesh cells as in double precision. The advantage of the DP solver is that one does not have to worry about the potentially adverse effect of round-off errors. It will however incur double the memory usage.
Pre-Processing
VR Editor
Post-Processing
VR Viewer
AUTOPLOT
PHOTON
Installation
How can I download the latest version of PHOENICS?
You can download the latest delivery version of PHOENICS 2015 (dated 31 March 2016) from CHAM’s FTP site, in either 32bit or 64bit form, by following the download instructions, a copy of which can be provided on request by emailing support@cham.co.uk with details of your name and organisation.
Application Examples
Combustion
PHOENICS-FLAIR
General
What is PHOENICS-FLAIR?
FLAIR is a special-purpose program for Heating, Ventilation and Air Conditioning (HVAC) systems that are required to deliver thermal comfort, health and safety, air quality, and contamination control. FLAIR provides designers with a powerful and easy-to-use tool which can be used for the prediction of airflow patterns, temperature distributions, and smoke movement in buildings and other enclosed spaces.
Extensive documentation of PHOENICS-FLAIR is provided here: [7], and here: [8].
Comfort Indices
How can I evaluate outdoor thermal comfort?
There are a range of indices in the literature which can used for evaluating outdoor thermal comfort. Each of these indices accounts for most if not all of the following effects: air temperature and humidity, local wind speed, radiation, activity and clothing.
FLAIR provides two thermal comfort indices that have been used to assess outdoor comfort, namely the Predicted Mean Vote and the Apparent Temperature. These can be activated from the VR Environment by clicking on " Menu">"Model">"Comfort Indices - Settings", and then toggling to "ON" the required comfort indices.
The Thermal Sensation Index is gaining popularity, but it has yet to be implemented in FLAIR. However, you can use the In-Form to compute this index from the Q1 input file, as explained below, but first a discussion will be provided for the comfort indices available in FLAIR.
1. Predicted Mean Vote (PMV)
The PMV is a comfort index defined in ISO 7730 that predicts the mean value of the votes of a large group of people on a 7-point thermal sensation scale. This index was developed to provide a measure of indoor comfort, but it has been frequently applied in both original and modified forms for outdoor applications (see for example Walls et al [2015] at
http://anzasca.net/wp-content/uploads/2015/12/107_Walls_Parker_Walliss_ASA2015.pdf). The FLAIR option for PMV is described here: http://www.cham.co.uk/phoenics/d_polis/d_chron/ph361/chap5.htm#chap5.2.
2. Apparent Temperature
The Apparent Temperature is a general term for the perceived outdoor temperature, caused by the combined effects of air temperature, relative humidity, radiation and wind speed. The formulae for the Apparent Temperature used in FLAIR are those specified by the Australian Bureau of Meteorology. They are an approximation of the value provided by a mathematical model of heat balance in the human body. They can include the effects of temperature, humidity, wind-speed and radiation.
Two forms of Apparent Temperature given by the Australian Bureau have been added to PHOENICS-FLAIR based on a mathematical model of an adult walking outdoors in the shade, one with radiation and the other without. The expression without radiation is suitable for the working condition of people walking in the shade; the expression with radiation is suitable for the working condition of people in direct sunlight.
The non-radiative apparent temperature (TAPP) is computed from:
- AT = Ta + 0.33*e − 0.70*ws − 4.0 (1)
whereas the apparent temperature with radiation (TAPR) is calculated from:
- AT = Ta + 0.348*e − 0.70*ws + 0.70*Q/(ws + 10) − 4.25 (2)
where Ta is the Dry bulb temperature (°C) (TEM1 in FLAIR), e is the water-vapour partial pressure (hPa) (PVAP in FLAIR), ws is the wind speed (m/s) at 10m elevation (taken to be the local VABS in FLAIR), and Q is the net radiation absorbed per body surface unit area (W/m2), which can be expected to be in the range of -20 to 130 W/m2.
The vapour pressure e is calculated from the temperature and relative humidity using the equation:
- e = (rh/100)*6.105*exp(17.27*Ta/(237.7 + Ta ) ) (3)
where rh is the Relative Humidity [%] (RELH in FLAIR).
Note that in order to calculate either form of the Apparent Temperature, the solution of Specific Humidity and the derivation of Relative Humidity must both be turned on in FLAIR.
3. Thermal Sensation Index (TSI)
The TSI determines a measure between 0 and 7, with 4 as the most comfortable condition. The original formula was developed in Japan and included the effects of relative humidity, but was analysed further by Givoni et al [2003] to provide the simplified equation quoted in your email:
- TSI = 1.2 + 0.1115*Ta + 0.0019*SR – 0.3185*ws (4)
where Ta is the air temperature (°C), SR is the horizontal solar radiation flux (W/m2) and ws is the wind speed (m/s). Further information is provided by Givoni et al (2003) in https://www.researchgate.net/publication/223626731_Outdoor_comfort_research_issues.
SR can be interpreted as the total amount of solar radiation received from above by a surface horizontal to the ground. This will include both direct and diffuse solar radiation, where the former is solar radiation that comes in a straight line from the direction of the sun at its current position in the sky. Diffuse radiation is solar radiation that has been scattered by particles in the atmosphere, and comes equally from all directions.
Library Case 674 provides a Wind and Solar Heating example which can be used as a vehicle for computing the TSI by means of the following PIL and In-Form statements.
- save7begin
- REAL(SRH,DIRSRF,DIFSRF,SOLALT,PI);PI=3.141593
- SOLALT = 19.85976 ! Solar altitude in degrees
- DIRSRF = 136.0 ! Direct solar radiation flux in W/m^2
- DIFSRF = 446.0 ! Diffuse solar radiation flux in W/m^2
- SRH = DIFSRF + DIRSRF*SIN(SOLALT*PI/180.) ! Horizontal solar radiation flux
- (STORED of TSI is (1.2 + 0.1115*TEM1 + 0.0019*:SRH: - 0.3185*VABS) with IMAT<100)
- save7end''
This Library Case uses a Weather File for Gatwick in the UK ( see http://www.cham.co.uk/phoenics/d_polis/d_wkshp/windsun/windsun2.htm ), and so the foregoing solar data are taken from the PHOENICS RESULT file.
Application Examples
Pollutant-gas dispersion in the atmosphere
This Q1 input file simulates steady, 3d pollutant-gas dispersion from a roof-top vent on an isolated building in a neutral wind environment File:GasDispersion.pdf