Gravitational body forces are represented in PHOENICS by way of sources of momentum in the equations for U1, U2, V1, V2, W1 and W2.
2.1 General remarks
2.2 The macros
2.3 Other means
PHOENICS provides for the gravitational force to be represented in the momentum equations in four distinct ways, namely by way of:
The the gravitational force per unit volume which appears in the vector momentum equation, F, is given by:
-
F = rho * g (1)
- -
where rho is the fluid density and g is the gravitational vector. - This formulation can be employed in PHOENICS. It is referred to below as the direct formulation.
It is often convenient to work, instead, with a perturbation density (rho - rho_ref), wherein rho_ref is a constant, often chosen as either the space-average density or as the ambient fluid density.
This practice results in the appearance of the following buoyancy-source term in the momentum equation:
F = ( rho - rho_ref ) * g (2)
- -
Thereafter it is convenient to define a reference pressure p_ref, which satisfies the hydrostatic-equilibrium equation:
grad p_ref + rho_ref * g = 0 (3)
-
Then the place of the pressure in the momentum equation can be taken by the perturbation pressure (p - p_ref). This device may also be used in PHOENICS; it is referred to as the reduced-pressure method.
Sometimes the variations of density within the field are small enough to be negligible everywhere EXCEPT in the gravitational term.
If, further, the density is supposed to be a linear function of some other solved-for variable, such as temperature, enthalpy or salinity, density variations make no explicit appearance in the momentum equations.
Formulations of this so-called (after its originator) Boussinesq type, to be found in PHOENICS, include:
rho - rho_ref = - rho_ref * dvo1dt * (T - Tref) (4) and:
rho - rho_ref = - rho_ref * (dvo1dt/cp) * (H - Href) (5)
wherein:
dvo1dt = the volumetric thermal expansion coefficient,
T = temperature, H =
enthalpy, and
cp = specific heat.
In practical circumstances arising in connexion with marine waters, the action of gravity may be influenced by both temperature and salinity (ie the dissolved-salt concentration) at the same time.
PHOENICS therefore allows the gravitational force per unit mass to depend linearly on two fluid properties, which the user is free to select.
This is allowed, however, only for Cartesian coordinate systems.
The manner of introducing the various momentum-source terms depends upon:-
All combinations can be expressed by appropriate sets of PIL statements, inserted in the Q1 file.
However, since remembering all of these is rather troublesome, certain groups of settings have been provided in the PHOENICS Input-file Libraries, to reduce the strain on memory.
These so-called "macros", or "PIL fragments", can be activated by the use of the load($....) command (or its shorter equivalents l($.... or #.... ).
These macros should be placed in the Q1 file below the statements which
determine the grid and the variables to be solved. The user is then free to
determine, after he has typed:
#gravity
so as to start the process, only whether to type:
gravdir = 1, 2, 3, 4, 5 or 6,
according to whether the gravity vector
is directed towards the north, south, east, west high or low directions, and
then:
rhoref = ... , tref = ... ,
or other statements to be described, and then:
#density or #densdiff or #bouss or #extbouss as desired.
These choices govern most situations; and what they produce can be easily inspected and then modified by hand, whenever necessary.
The statements required after gravdir has been entered depend upon which option will follow.
Whichever option is chosen, the gravitational acceleration, defaulted by PHOENICS to 9.81, can be reset by entering:
gravacc = ...... (the desired value).
If #densdiff is the macro which is about to be loaded, the user should first enter:
rhoref = ...... (the desired value of the reference density).
If #bouss is to follow, the user first enters:
tref = ...... (the reference temperature), if temperature is being solved; or
href = ...... (the reference enthalpy), if the solved-for variable is enthalpy.
If the properties of the fluid are being read from the props file, as occurs when the macro #use_props has been previously entered, that is all that is needed; for the expansion coefficient will be read from it.
Otherwise the user must enter:
dvo1dt = ..... (the appropriate coefficient)
If enthalpy is the variable, dvo1dt should be divided by cp (See eq.5).
When the following command will be #extbouss, the user is required to enter the values of:
BUOYA, BUOYB, BUOYC, IBUOYB, and IBUOYC which correspond to the
formula:
force/mass = buoya + buoyb * variable_ibuoyb + buoyc * variable_ibuoyc
The indices IBUOYB and IBUOYC could, for example be 14 and 16, if the variables in question were H1 and C1.
Thus an example of a complete set of commands for this option is:
#gravity
gravdir=2
BUOYA=0.0; BUOYB=0.000981; BUOYC= - 0.0981
IBUOYB=14; IBUOYC=16
#extbouss
Evidently, use of these macro commands minimizes the user's data-input actions.
This PIL fragment is to be found in file d_earth/d_core/inplib/074.htm
It is activated by the command #gravity, because the character variable gravity has been declared as $L074 in PIL fragment 14.
This macro declares some further variables for later use, namely:
DENSITY=$075; DENSDIFF=$076; BOUSS=$077;EXTBOUSS=$079; BUOY=$078
Then it makes some preliminary settings, and prints out some messages.
The four PIL fragments may be seen by clicking here:
density, densdiff, bouss, and extbouss.
They make the ground-number settings needed for activation of sequences in the subroutine GXBUOYSO; they then load PIL fragment buoy.
The reason for including the "start of..." and "end of..." lines is that, when these macros are used in interactive mode, and the answer "yes" is given to the "copy to Q1" question, the content of the macros is included in the new Q1 file.
Some users like subsequently to remove them by editing, replacing them by the original #... commands. The start and end markers assist this "cleaning out" process.
The functions of the buoy PIL fragment are:-
The PIL fragment buoy responds to the previously-made settings of:-
buoyopt, gravdir, rhoref, tref, href, CARTES, BFC, NX, NY, NZ
as well as to whether the temperature or the enthalpy is the energy- equation variable, and to whether fluid properties are being obtained by way the settings of the PRPS variable, and its references to the PROPS file in sub-directory d_EARTH.
The macros, as explained above, are convenient means of writing PIL statements which will then activate coding sequences to be found in subroutine GXBUOYSO.
However, it is not necessary to use them; and users who wish to have control of every move may prefer to make direct settings themselves.
The present section of this article is for such users.
The sequences in GXBUOYSO are entered when a PATCH is defined with a name beginning with the letters BUOY (either upper- or lower-case).
When direct settings are created by the user, however, any other (non-reserved) name may be used.
For cartesian coordinates, and for cylindrical-polar coordinates when the z-axis is vertical, gravity sources of the form of equation (1) above can be inserted via PATCH and COVAL statements as follows.
PATCH(GRAVITY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) ,
followed, if the z-direction is vertically upwards, by:
COVAL(GRAVITY,W1,FIXFLU, -9.81)
In two-phase flows, the following additional specification would be required:
COVAL(GRAVITY,W2,FIXFLU, -9.81)
For cartesian coordinates only, if it is the x- or y-directions which are vertically upward, the COVAL statements should be either:
COVAL(GRAVITY,U1,FIXFLU,- 9.81) or COVAL(GRAVITY,W1,FIXFLU,- 9.81)
The VAL argument, -9.81 , is of course, the gravitational acceleration in m/s**2 .
If the gravitational direction were inclined to the axes, the COVAL statements would take the form:
COVAL(GRAVITY,U1,FIXFLU, -9.81 * COS(uangle))
COVAL(GRAVITY,V1,FIXFLU, -9.81 * COS(vangle))
COVAL(GRAVITY,W1,FIXFLU, -9.81 * COS(wangle))
whereby it must be remembered than uangle, etc are not PIL variables; so the user must either insert numerical values; or declare uangle etc as PIL variables and then devise means of assigning values to them.
The PATCH type PHASEM ensures that the gravitational acceleration, set as the "VALue" argument, is multiplied by the mass of the phase in the velocity cell; this makes the momentum source equal to the force acting on the velocity cell.
For a gravity vector not aligned with any one of the coordinate directions, COVALs are required for all velocities, with the last argument of COVAL specifying the resolute of the gravitational acceleration due to gravity in the respective coordinate directions.
These can be assigned as constants only for cartesian coordinates, ie for CARTES=T and BFC=T.
Otherwise GXBUOYSO must be activated.
For cartesian, cylindrical-polar and body-fitted coordinates, options are supplied in subroutine GXBUOYSO, called from EARTH for PATCH names commencing with the characters: BUOY.
For cartesian, cylindrical-polar and body-fitted coordinates, the parameters BUOYA, BUOYB and BUOYC signify the resolutes of the gravitational-acceleration vector in the cartesian-coordinate directions XC, YC and ZC respectively.
Gravity is taken as positive in the positive cartesian-coordinate direction.
The fourth argument of COVAL (ie. the "value" ) selects one of four multipliers of the resolute of gravity in the local direction of U1, V1 and W1.
Four values are currently recognised, namely GRND1, GRND2, GRND3 and GRND4 (or DENSITY, DENSDIFF, BOUSS, LINBC).
GRND1 sets a source equal to: RHO1 * VOL * gravity resolute, which is of the form of equation (1).
GRND1 may be replaced by the word DENSITY, which is easier to remember.
GRND2 selects a multiplier equal to:
(density - BUOYD) / density
which results in a source of the form of equation (2) above.
With BUOYD set to a buoyancy reference density, this option gives the buoyancy force that arises from differences of the density field from this reference value.
The absence of the hydrostatic component in the pressure field also simplifies the specification of pressure boundary conditions; for often the 'reduced' pressure is uniform at the boundaries.
GRND2 may be replaced by the word DENSDIFF, which is easier to remember.
In the following example, the settings specify the buoyancy-force option for gravity directed along the x-direction in cartesian coordinates:
BUOYA=9.81 PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(BUOY,U1,FIXFLU,DENSDIFF)
Here BUOYB and BUOYC are left at their defaults of zero, with BUOYD being set equal to the reference density.
GRND3 selects a multiplier equal to:
BUOYE + BUOYD * ( enthalpy or temperature )
which is results in a buoyancy source of the form of equation (4) when the temperature TEM1 is solved with:
BUOYD= - dvo1dt and BUOYE = - BUOYD * Tref
where dvo1dt is the volumetric coefficient of expansion of the fluid,
This is achieved by setting BUOYD=GRND10 and BUOYE=Tref in the Q1 file.
GRND3 may be replaced by the word BOUSS, which is easier to remember.
When the enthalpy H1 is solved, the buoyancy source term is of the form of equation (5) with:
BUOYD= - dvo1dt/cp and BUOYE = - BUOYD * href
which is given by -dln{V}/dT or -{drho/dT}/rho.
For body-fitted coordinates, when gravity acts in the Cartesian-Z direction, the following settings activate the Boussinesq approximation:
BUOYC=9.81
PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
COVAL(BUOY,U1,FIXFLU,BOUSS)
COVAL(BUOY,V1,FIXFLU,BOUSS)
COVAL(BUOY,W1,FIXFLU,BOUSS)
where BUOYD and BUOYE are set as described above.
GXBUOYSO also provides for a gravitational force per unit mass or volume of fluid of the following form:
(BUOYA + BUOYB * IBUOYB + BUOYC * IBUOYC)
where IBUOYB and IBUOYC are user-defined integers and BUOYA, BUOYB and BUOYC are user-defined real variables.
The integers are the block-location indices of the dependent variables in question.
This option is selected by using GRND4 as the fourth argument of the COVAL statement, e.g.
PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(BUOY,V1,FIXFLU,LINBC)
Library case W975 makes use of this option.
In two-phase flows in horizontal ducts, in which a liquid layer lies beneath a gas layer, a horizontal force arises from the vertical variations of the depth of the liquid layer.
The PATCH-type RGRAD provides a force proportional to the horizontal gradient of the volume fraction (see RGRAD) .
This option is suitable for a duct of rectangular cross-section.
For circular-sectioned ducts, the option activated by the PATCH name LATG (which causes GREX3 to call subroutine GXLATG) is required.
For further details, the Fortran of subroutine GXLATG, in file GXSOR.F, may be consulted.
The effects of gravity on the turbulence field may be accounted for by the introduction of buoyancy source terms into the transport equations for KE and EP.
For stably-stratified flows, the user is advised that many workers recommend that such a source term should be omitted from the EP equation.
For further details on this option, the PHENC entry 'K-Epsilon turbulence model' should be consulted.