This lecture explains how boundary conditions and sources are specified in PHOENICS.
The basis of the implementation is given, and some simple examples taken from library cases are given.
Full details of Boundary Conditions and Sources are available through the Satellite help-facility, the PHOENICS Encyclopaedia, and in TR/200. Many examples can be found in the Library.
NOTE: By default all domain edges are impervious to flow, frictionless and adiabatic. They represent symmetry planes or axes.
Boundary And Internal Conditions And Sources
Differential equations need to be supplemented by boundary conditions before they can be solved.
The boundary conditions which define a fluid- or heat-flow problem usually convey the necessary information about how much fluid enters the domain, where it can leave, what is its temperature on entry, what are the temperatures of the walls, etc.
Of course, fluid may be caused to enter or leave at points within the domain, not only at its external limits; and temperatures of structural elements within the domain may also be externally imposed, and exert an influence upon the flow.
PHOENICS makes no distinction between "boundary" and "internal" conditions, or between these and "source" terms, so the former term will be used here for all of them.
Treatment Of Boundary Conditions In PHOENICS
In PHOENICS, boundary conditions and sources appear on the r.h.s. of the differential equation for a variable f. Thus:
where:
Sf represents conventionally recognised source terms, such as pressure gradients or viscous heating terms. These are 'built-in' to EARTH.
Gf represents the diffusive exchange coefficient for f.
Sbc1 etc. represent various boundary conditions. These may be present only in certain regions of the domain. More terms of this kind may be also be present in other regions of the domain, and these regions may overlap.
The differential equation is integrated over a control volume to yield the finite-volume equation actually solved. The integral of the boundary source is represented in linearized form:
The finite-volume discretization of the differential equation thus yields, for each cell P in the domain, the following algebraic equation:
where:
Sf is the 'true' source
C is the coefficient
V is the value
T is the type, a geometrical multiplier
As a consequence of the integration procedure, the source is required per cell. The units of the source are (f kg/s). The type is used to convert the source from any given set of units.
Thus, if the source is defined as 'per unit volume', type supplies the cell volumes. If the source is 'per unit area', type is an appropriate cell face area.
Note that TCV is added to the numerator, and TC to the denominator. As will be shown, this allows for easy manipulation of the solution.
PIL Commands For Boundary Conditions
As set out in preceding panels, the specification of boundary conditions requires two kinds of information:
(1) and the first part of (2) is conveyed by a PATCH command:
PATCH (name, type, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL)
where:
name is a unique patch name for future reference;
type is T
IXF,IXL are the first and last IX in the patch; and similarly for y, z and t (time)
NOTE that any of the PATCH limits can be specified in terms of REGION number by using #IXF,#IXL etc.
The remainder of (2) is specified by a COVAL command, the format of which is:
COVAL (name, variable, coefficient, value),
where:
name is the patch name to which the command refers
variable is a SOLVEd-for variable
coefficient is C
value is V
How these are used can be seen by inspecting Library cases, or Menu-generated Q1 files. Some simple examples are now given.
Further information can be found in the Encyclopaedia, under the headings: Boundary Conditions, Sources, PATCH, COVAL
Examples Of Boundary Condition Settings
The types of boundary conditions which have to be provided for are:
Each of these will now be explained and illustrated with examples.
Fixed Value Boundary Condition
Practical example: we wish to fix the temperature in one corner of a cube to 0.0, and to 1.0 in the diagonally opposite corner.
Numerical practice: the value of phi can be fixed in any cell by setting C to a large number, and V to the required value.
The equation then becomes:
The PIL variable FIXVAL is provided for this purpose. A typical PATCH/COVAL would be:
PATCH (FIXED, CELL, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL) COVAL (FIXED, phi, FIXVAL, required_value)
Example From The Library
Case 100 models a solid cube of material, in which one corner is held at a low temperature, and the diagonally opposite corner is held at a high temperature:
GROUP 13. Boundary conditions and special sources **Corner at IX=IY=IZ=1 PATCH(COLD,CELL,1,1,1,1,1,1,1,1) **Fix temperature to zero COVAL(COLD,TEMP,1.E2,0.0) **Corner at IX=NX, IY=NY, IZ=NZ PATCH(HOT,CELL,NX,NX,NY,NY,NZ,NZ,1,1) **Fix temperature to 1.0 COVAL(HOT,TEMP,1.E2,1.0)
Note that the coefficient has been set to 100. This in effect sets a finite conductivity between the external value and the centre of the cell. FIXVAL would lock the cell-centre value to the external value.
Fixed Flux / Fixed Source Boundary Condition
Practical example: heat is being generated at a constant (fixed) rate.
Numerical practice: a fixed source can be put into the equation by setting C to a small number, so that the denominator is not changed, and by setting V to (source/C). T then ensures that the final source is per cell.
The equation then becomes:
The PIL variable FIXFLU is provided for this purpose. A typical PATCH/COVAL setting would be:
PATCH (SOURCE, area, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL) COVAL (SOURCE, phi, FIXFLU, source_per_unit_area)
Library case 921 concerns the prediction of the flow and temperature fields in a closed cavity with one moving wall and a heated block. The heated block appears as shown below:
Heat source in block - 10 MW/m**3
PATCH(HEATEDBL,VOLUME,NX/4+1,3*NX/4,NY/4+1,3*NY/4,1,1,1,1) COVAL(HEATEDBL,TEM1,FIXFLU,1.0E7)
The patch type VOLUME converts the source from W/m3 to W per cell, by multiplying by the appropriate cell volumes.
Practical example: One of the domain boundaries is losing heat to the surroundings. The external heat transfer coefficient, H (W/m2/K), and the external temperature, Text (K), are both known and constant.
Numerical practice: The heat source for a cell with area A is:
Q = A H (Text - Tp)
This is obviously in TC(V-phi) form if T=A, C=H and V=Text. A typical PATCH/COVAL would be:
PATCH (HEATLOSS, area, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL) COVAL (HEATLOSS, TEM1, heat_transfer_coefficient, external_temperature)
In the previous example, either or both the heat transfer coefficient and external temperature are likely to be a function of some solved-for variable or other suitable expression.
The non-linear source can always be linearized into TC(V-phi) form by arranging to update C and/or V during the course of the calculation.
The updating of C and V is signalled to EARTH by setting C or V to one of the 'ground flags' - GRND, GRND1, GRND2, ..., GRND10.
Thus the following COVALs may be seen:
COVAL (HEATLOSS, TEM1, GRNDn, external_temperature)
or
COVAL (HEATLOSS, TEM1, heat_transfer_coefficient, GRNDn)
or
COVAL (HEATLOSS, TEM1, GRNDn, GRNDn)
Many commonly-found examples have already been coded into GREX by CHAM.
Laminar Wall Boundary Condition
Laminar shear stress at stationary wall is expressed by:
where area is the cell face area, and dely is the distance from the cell face to the cell centre.
This can be put into 
form if 
.
The problem with this approach is that the density and laminar viscosity may be varying, whilst the distance to the wall will change as the grid is refined, and indeed may change from cell to cell in a BFC grid. This simple method is therefore not recommended, as all the quantities causing problems to the user are known to EARTH.
A special PATCH type is provided which automatically sets:
The coefficient is then a further multiplier, which is usually set to 1.0.
These PATCH types are NWALL, SWALL, EWALL, WWALL, HWALL and LWALL.
To further simplify matters, a WALL command generates all the PATCH and COVAL statements required.
A typical WALL command would be:
WALL (SIDE, area, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL)
Note that WALL uses the ordinary area types, not those above.
Example From The Library
The stationary and moving walls in case 921 are specified as shown:
Moving wall at South side of domain at 1 deg
WALL (MOVING,SOUTH,1,NX,1,1,1,1,1,1)
COVAL(MOVING,U1,1.0,-WALLVEL);COVAL(MOVING,TEM1,1.0,1.0)
Stationary wall at North side at 0 deg
WALL (NORTHW,NORTH,1,NX,NY,NY,1,1,1,1)
COVAL(NORTHW,U1,1.0,0.0);COVAL(NORTHW,TEM1,1.0,0.0)
Stationary wall at West side at 0 deg
WALL (WESTW,WEST,1,1,1,NY,1,1,1,1)
COVAL(WESTW,V1,1.0,0.0);COVAL(WESTW,TEM1,1.0,0.0)
Stationary wall at East side at 0 deg
WALL (EASTW,EAST,NX,NX,1,NY,1,1,1,1)
COVAL(EASTW,V1,1.0,0.0);COVAL(EASTW,TEM1,1.0,0.0)
Note that the coefficients have all been set to 1.0, and the values to the wall surface values. This suffices for laminar wall conditions.
Turbulent Wall Boundary Condition
In a turbulent flow, the near-wall grid node normally has to be in the fully-turbulent region, otherwise the assumptions in the turbulence model are invalid. The wall shear stress and heat transfer can no longer be obtained from the simple linear laminar relationships.
Unless a low-Reynolds number extension of the turbulence model is used, the normal practice is to bridge the laminar sub-layer with wall functions. These use empirical formulae for the shear stress and heat transfer coefficients.
Three types of wall function are available, selected by the COVAL settings:
Example From The Library
Library case 172 concerns the prediction of developing flow in a duct. The k-epsilon turbulence model is used. The duct surface at the north side of the duct is represented as:
GROUP 13. Boundary conditions and special sources ** North-Wall Boundary WALL (WFUN,NORTH,1,1,NY,NY,1,NZ,1,1) COVAL(WFUN,TEMP,GRND2,TWALL)
Note the use of the WALL command to locate the wall and activate wall friction effects.
The surface temperature is supplied via the COVAL for TEM1. The GRND2 in the coefficient slot activates logarithmic wall functions.
All mass flow boundary conditions are introduced as linearized sources in the continuity equation, with pressure (P1) as the variable. A mass source is thus:
where Cm and Vm are coefficient and value for P1.
At an inflow boundary, the mass flow is fixed irrespective of the internal pressure. This effect is achieved by setting Cm to FIXFLU, and Vm to the required mass flow.
The sign convention is that inflows are +ve, outflows are -ve. A fixed outflow rate can thus be fixed by setting a negative mass flow.
Example From The Library
Library case 274 concerns the flow over a simplified van geometry. The inflow boundary at the low end of the solution domain is represented as:
GROUP 13. Boundary conditions and special sources ** Upstream boundary INLET(UPSTR,LOW,#1,#NREGX,#1,#NREGY,#1,#1,1,1) VALUE(UPSTR,P1,14.0); VALUE(UPSTR,W1,14.)
Note that for an INLET, the VALUE command for P1 sets the mass flux. This is often set as RHOIN*VELIN, the (inlet density) * (inlet velocity).
The mass flux is fixed, and the in-cell pressure is allowed to float.
The VALUE command for W1 sets the velocity of the inflowing stream. In this case all other variables are taken to be 0.0 at the inlet. If they are not, then VALUE commands would have to be added.
Fixed Pressure Boundary Condition
This is the case of a mass flow boundary where the pressure is fixed irrespective of the mass flow.
As with any other variable, the pressure is fixed by putting a large number for Cm, and the required pressure for Vm.
For numerical reasons, FIXVAL tends to be too big. A Cm of about 1E3 usually suffices. The Encyclopaedia gives further guidance.
The direction of flow is then determined for each cell in the PATCH by whether Pp>Vm, or Pp< Vm. The first produces local outflow, the second local inflow.
Example From The Library
The exit boundary in case 274 is a fixed pressure boundary, set as:
Downstream boundary PATCH(DWSTR,HIGH,#1,#NREGX,#1,#NREGY,#NREGZ,#NREGZ,1,1) COVAL(DWSTR,P1,FIXP,0.) COVAL(DWSTR,U1,ONLYMS,0.0);COVAL(DWSTR,V1,ONLYMS,0.0) COVAL(DWSTR,W1,ONLYMS,0.0)
In this case, the in-cell pressure is fixed by the COVAL for P1, and the mass flux is adjusted to satisfy continuity. The direction of flow is determined by whether the in-cell pressure is > or < the fixed value.
The COVALs for U1, V1 and W1 are supplied in case part of the boundary should be an inflow - they specify velocities to be brought in. They are not used in cells where in-cell pressure > external.
PHOENICS makes no distinction between boundary conditions and source terms. Both are represented in the by now familiar TC(V - phi) form.
In the lecture so far, we have made a distinction between the 'true' source, Sf, and the boundary source, Sbc.
The true source represents the fundamental parts of an equation which are not covered by the convective, diffusive or transient terms. These are coded in EARTH, and hence are called the built-in sources.
Built-in Sources
Examples of these are:
The built-in sources can be switched off for individual variables, by setting N as the first argument of TERMS. For example:
TERMS (W1, N, P, P, P, P, P)
will deactivate the pressure gradient terms in the W1 equation.
User-set Sources
All other sources are treated as external user-set sources. Examples are:
Many commonly used examples of the second type have already been programmed by CHAM. The source code for all of these can be found in the exemplary GROUND routine, GREX, and several other files whose names start with GX.
All of these can be found in phoenics/d_earth/d_core/..., and can be viewed through POLIS. If required, they can be copied and used as templates for further modification.
Implementation
Boundary conditions (or sources) cannot always be specified through a constant coefficient and a constant value.
In many instances, the coefficient and/or the value are functions of one or several solved-for variables, the value of which cannot be foreseen at the time of data input.
As we have already seen, PHOENICS copes with these complex relationships through the insertion of FORTRAN coding in the GROUND module.
To do so, special flags (GRND, GRND1, ... GRND10) can be specified as coefficients and/or values in the COVAL command. These instruct EARTH to visit special sections of GROUND, where coefficients and/or values can be computed and set back into EARTH.
Nearly all the conditions specified in GREX are of this type.
Non-linearity
All user-defined sources are introduced in the linear 
form. PATCH
and COVAL commands are used to define the location, duration and T, C and V.
Linear sources can be inserted directly from Q1 or the Menu.
Example: A chemical species, C1, is being created throughout the whole domain at a constant rate of 10 kg/m3/s. The PIL commands are:
PATCH (CREATEC1, VOLUME, 1,NX, 1,NY, 1,NZ, 1,LSTEP) COVAL (CREATEC1, C1, FIXFLU, 10)
However, many common CFD sources cannot be expressed directly in linear form, because:
Examples Of Non-Linear Sources
Example 1: The pressure drop through a porous medium can frequently be written in the form:
This can be turned into a momentum source (force):
where vol is the cell volume.
The source is proportional to velocity squared.
Example 2: The sink of a chemical species (C1) undergoing chemical reaction can frequently be expressed by the Arrhenius law:
where a and b are empirical constants, C2 is the mass fraction of another species, and T is the temperature.
Although the expression is linear in C1, the local values of C2 and T cannot be predetermined.
Note that the Arrhenius expression is available as an option in the standard PHOENICS combustion model.
Treatment Of Non-Linear Sources
Non-linear and interconnected sources will normally allow several
representations in the linear 
form.
Practice 1: The source can be introduced as a fixed flux (FIXFLU) source, with the value of the source computed from in-store values.
This practice is not recommended, as it tends to be numerically unstable. However, in some cases it is the only way.
Practice 2: If at all possible, the source should be linearized.
Frequently, there will be more than one way to linearize a source. Only experience will show which is the most stable.
Examples Of Source Linearization
The quadratic momentum source from above can be written as:
where v* is the current in-store velocity. In the converged solution, v = v*, and the source is the required one.
This is in 
form if:
The PIL implementation for, say, the W1 equation is:
PATCH (PDROP, PHASEM, 1,NX, 1,NY, 1,NZ, 1,LSTEP) COVAL (PDROP, W1, GRND, 0.0)
Note that the type PHASEM sets T to r*Vol .
The GROUND code would set coefficient = K W1*.
The Arrhenius expression in panel 30 can be written as:
This is in C(V-phi) form if:
The PIL implementation is (assuming that the source is 'per unit volume'):
PATCH (REACT, VOLUME. 1,NX, 1,NY, 1,NZ, 1,LSTEP) COVAL (REACT, C1, GRND, 0.0)
The GROUND code would set:
So far, it has been stressed that linear sources can be specified directly from Q1, but that non-linear sources need additional coding in GROUND. Many such sources have already been provided by CHAM.
In addition to the GROUND code already supplied in GREX, a wide range of non-linear sources CAN be introduced directly from Q1 without the need for any GROUND coding. A range of these will now be described. A full list can be found in the Encyclopaedia, under the PATCH entry.
These sources are introduced by 'unconventional' coefficient settings, and/or by special PATCH names.
Normally, the coefficient in the 
source must be
positive, otherwise numerical instability may result. In view of this, negative
coefficients in the Q1 are used to flag special sources.
The following non-linear Q1-set sources will be exemplified:
These sources can also all be set in VR, using the 'User-Defined' object type.
Stagnation Pressure Boundary Condition
For incompressible flows, the stagnation pressure is:
This can be manipulated to give the mass flux:
For P1 (and also P2), a negative coefficient flags a mass source of the form:
where Cm and Vm are the coefficient and value for mass (P1).
The desired stagnation source can thus be set by putting:
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.
VRE example:
The Outlet object provides the functionality to specify the stagnation / total pressure via the 'Quadratic' option for the coefficient in the Outlet attributes dialog box. It does so in a modifed manner compared to what was described above, by already incorporating the term 2*rho1 into the formulation. This way the stagnation pressure is set to the user specified value with a coefficient of 1.0. Further details can be found at Quadtratic Option for outlet object.
If a negative coefficient is supplied for any other variable, the source introduced is:
PIL example: The quadratic pressure drop shown on panel 29 can thus be implemented without any GROUND code, as shown below:
PATCH (PDROP, PHASEM, 1,NX, 1,NY, 1,NZ, 1,LSTEP) COVAL (PDROP, W1, -K, 0.0)
Note that the negative coefficient options only apply to constant coefficients. They will not apply to GROUND-set coefficients.
Sources of the form:
can be introduced from Q1 by using 'special' PATCH names.
The format for PATCH and COVAL to be used is shown below:
PATCH (*phi2abc, type, IXF,IXL, IYF,IYL, IZF,IZL, ITF,ITL) COVAL (*phi2abc, phi1, a, b)
where:
abc - is an identifying character string
phi2 - is the name of a SOLVEd variable
phi1 - is the variable with the source
a, b - are the constants.
Note that if the name of phi2 is less than 4 characters, blanks must be left before the supplementary identifying string.
PIL example: The pressure drop from panel 19 can be generalised to give a momentum sink of:
If n = 1.98, say, this can be introduced for the W1 equation as:
PATCH (*W1 RES, PHASEM, ...... COVAL (*W1 RES, W1, K, 0.98)
Radiative Heat Loss To The Surroundings
Radiative heat transfer is based on a T4 source. Heat loss to the surroundings can thus be calculated from:
where s is the Stefan-Boltzmann constant, and E is the emissivity.
The '*name' technique is extended to cover this situation. The above source is introduced into the TEM1 equation by PATCH/COVAL commands of the form:
PATCH (*RADabcd, type, ..... COVAL (*RADabcd, TEM1, sigma * E, Text)
The variable TEMP0 is used to convert the temperature equation to Kelvin. It is thus set to 0 if TEM1 is in K, and to 273 if TEM1 is in Centigrade.
PIL example: The East side of the domain is losing heat through radiation to surroundings at 23 deg C, and the West side to surroundings at 50 deg C. The emissivity of the cold side is 0.8, and that of the hot side is 0.95.
REAL(SIGMA); SIGMA = 5.6697E-8 PATCH (*RADCOLD, EAST, NX,NX, 1,NY, 1,NZ, 1,LSTEP) COVAL (*RADCOLD, TEM1, SIGMA * 0.8, 23) PATCH (*RADHOT, WEST, 1,1, 1,NY, 1,NZ, 1,LSTEP) PATCH (*RADHOT, TEM1, SIGMA * 0.98, 50) TEMP0 = 273
Boundary conditions and sources are treated in PHOENICS as linearized sources having
the form 
.
Two PIL commands convey the information required: a PATCH command which sets "where", "when" and T, and a COVAL command which sets the values of C and V.
For a variable f, the main kinds of boundary conditions are:
Boundary conditions for mass and pressure are both treated as linearized sources in the
continuity equation, with pressure as the variable in the linear source: 
. FIXFLU is
used as coefficient for the specification of mass fluxes, and FIXVAL or FIXP for fixing
the pressure.
Boundary conditions must be supplied for all the variables when there is an inflow mass into the domain. This is done by using ONLYMS as coefficient in the COVAL command for the property, and the inflowing value as value.
Linear sources, and certain non-linear sources can be introduced directly from the Q1 or Menu, by appropriate settings of coefficient and/or value.
Non-linear sources have to be linearized, and the non-constant part programmed in GROUND. This is flagged by using GRND, GRND1 ... GRND10 for coefficient or value. Many such sources are already provided in GREX.