label TOP
q1quit=t;dmpstk=f;nocomm=t;nocopy=t;nowipe=t
rset(v,0)
Preliminaries :
Satelite subroutine sat calls CHKINP for FORTRAN string comparison
namsat=chkc
Set up 4 character values, and one logical. Logical is set true
immediately prior to satrun call, and is returned false
Character strings are set to give prompts and accept data for test
boolean(log1)
char(cvl1,cvl2,cvl3,cvl4);cvl3='Correct!' !
QCOM statements are used to annotate the Q1 file, DCOM to
process command and annotate Q1 file
----- start of tutorial : -----
write group identifiers etc to Q1 file
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 1. Run title and other preliminaries
qcom( *
display
WELCOME TO PIL TUTORIAL NUMBER 1
********************************
This tutorial starts to teach you the basic PIL commands
and specifies a simple problem of flow in a channel
which can then be run by PHOENICS.
This tutorial defines the problem shown.
Wall
__________________________________________
| |
| |
| |
Inlet | | Outlet
Y | |
^ | |
| | |
| |__________________________________________|
----> X Wall
CALL LOAD_96
The tutorial is configured to expect valid PIL commands
to be entered by you at the terminal.
String comparison is then undertaken to ensure that you
have entered the correct data.
As such, it is important to note that the exact numbers
are required to be input.
* YOU MUST READ the panels presented to you
* The form of the PIL commands and other hints are given
NB: You can quit the tutorial at any time by entering "quit"
CALL LOAD_96
GRID
X-Direction Y-Direction
XLength : 0.5 YLength : 0.1
No. of regions: 1 No. of regions: 1
No. of cells : 10 No. of cells : -20
Power : 1.0 Power : 1.4
NB. Negative sign on No. of cells in Y indicates that these
are symmetrically distributed.
CALL LOAD_96
VARIABLES
Solved: P1, U1, V1 Stored : ENUT
PROPERTIES
Turbulence model : k-epsilon
Reference pressure (PRESS0): 1.E5
Density (RHO1 ): 1.22
Viscosity (ENUL ): 1.4654E-4
INITIAL FIELDS
KE=0.04 ; EP=0.13
CALL LOAD_96
BOUNDARY CONDITIONS
All PATCH LOCATIONS ARE: #1,#1,#1,#1,#1,#1,#1,#1
WALL
NAME : TOP ; BOT
TYPE : NORTH ; SOUTH
PIL var Coef Value
U1 loglaw 0.0
KE loglaw loglaw
EP loglaw loglaw
CALL LOAD_96
INLET
NAME : INLET
TYPE : WEST
PIL var Value
P1 10.0*rho1
U1 10
V1 0.0
KE 0.04
EP 0.13
CALL LOAD_96
OUTLET
NAME : OUT
TYPE : EAST
PIL var Value
P1 0.0
CALL LOAD_96
* Start of Tutorial 1 *
* * * EXPECTED title is 'tutorial.1'
Group 1 of the Q1 file is used to specify the title
of the prediction
The title should be a unique identifier. This is
presented on PHOTON output, and included in results
files.
The PIL command to set the title has the syntax:
TEXT(Required Text
The required text in this case is the title shown above
hence the input required at the prompt MUST be
text(Tutorial.1
cvl2=text(Tutorial.1
cvl4='Please_enter_the_TEXT_command_to_set_the_title_for_this_case'
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 2. Transience; time-step specification
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 3. X-direction grid specification
qcom( *
regions in x = 1
GROUP 2. Transience; time-step specification
GROUP 3. X-direction grid specification
There are several methods available within PHOENICS to
define the computational mesh. This Tutorial will
exemplify the method of REGIONS and gridpower
Hence, The first step in setting the mesh for this
tutorial is to define the number of regions in the X
direction
These region identities are then used to control the
mesh spacing according to other PIL commands
The number of regions in each direction is given
by a PIL command with the syntax
NREGdir=integer number
where :
dir is T,X,Y, or Z for time,or spatial co-ordinates
and the integer is the required number of regions
cvl4='Please_set_the_number_of_regions_in_the_x_direction_to_1'
cvl2=nregx=1
CALL LOAD_95
* * * * * iregx 1 + grdpwr
Having specified the number of regions in X,
we now have to specify the cell spacing in each region
To do this requires two operations:
1) to identify the relevant region
2) to specify the grid spacing
The identification of the relevant region is via
a PIL statement with the Syntax
IREGdir=integer number
where :
dir is T,X,Y, or Z for time,or spatial co-ordinates
the integer is the required region
cvl4='Please_enter_the_index_of_the_first_x_region_to_be_modified'
cvl2=iregx=1
CALL LOAD_95
Cell spacing in region 1 is determined via a PIL statement
with the syntax:
GRDPWR(dir,IC,LENGTH,PWR)
where :
dir is T,X,Y, or Z for time, and spatial co-ordinates
IC is the number of cells in the region
LENGTH is the extent of the region
PWR is the power that is used to distribute the cells
GRDPWR(x,10,100,1.4)
would put 10 cells in x over 100 metres, with small cells
at the beginning, large cells at the end.
Negative powers will put large cells at the beginning!
Negative numbers of cells symmetrically place cells
according to the power
reminder: xlength=0.5; no. of cells = 10; power=1.0
cvl4='Please_enter_the_distribution_of_cells_in_x_region_1'
cvl2=grdpwr(x,10,0.5,1.0)
CALL LOAD_95
do jj=1,25
mesg(
enddo
Congratulations,
You have now set the X direction mesh distribution.
This mesh will now be shown by use of the view statement.
NOTE: we have only defined the mesh in X thus far.
please press "enter" to view the mesh thus far
and then press "enter" again to continue the tutorial
CALL LOAD_96
view(k,1)
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 4. Y-direction grid specification
qcom( *
regions in y = 1
GROUP 4. Y-direction grid specification
Having completed the specification of the X Mesh,
The next task is to specify the Mesh in Y, and again
we will use the method of REGIONS and gridpower
Hence, The first step in setting the mesh in Y
is to define the number of regions in Y
These region identities are then used to control the
mesh spacing according to other PIL commands
The number of regions in each direction is given
by a PIL command with the syntax
NREGdir=integer number
where :
dir is T,X,Y, or Z for time, and spatial co-ordinates
the integer is the required number of regions
cvl4='Please_set_the_number_of_regions_in_y_to_1'
cvl2=nregy=1
CALL LOAD_95
* * * iregy 1 + grdpwr
The specification in Y is analogous to that in X.
Remember, To do this requires two operations:
1) to identify the region.
2) to specify the grid spacing
Remember, the identification of the relevant region is
via a PIL statement with the Syntax
IREGdir=integer number
And the grid spacing is specified with
GRDPWR(dir,IC,LENGTH,PWR)
where :
dir is T,X,Y, or Z for time,or spatial co-ordinates
the integer is the current region to be modified
IC is the number of cells in the region
LENGTH is the cartesian extent of the region
positive/negative power expands/contracts cells
reminder: ylength=0.1 ; No. of cells =-20; power=1.4
cvl4='Please_enter_the_index_of_the_y_region_to_be_modified'
cvl2=iregy=1
CALL LOAD_95
cvl4='Please_enter_the_distribution_of_cells_in_y_region_1'
cvl2=grdpwr(y,-20,0.1,1.4)
CALL LOAD_95
Congratulations,
You have now completed the mesh specification in the
XY plane for the tutorial
This mesh will now be shown by use of the view statement
please press "enter" to view the mesh
and then press "enter" again to continue the tutorial
CALL LOAD_96
view(k,1)
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 5. Z-direction grid specification
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 6. Body-fitted coordinates or grid distortion
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 7. Variables stored, solved & named
qcom( *
variables to be solved
GROUP 5. Z-direction grid specification
GROUP 6. Body-fitted coordinates or grid distortion
GROUP 7. Variables stored, solved & named
The PHOENICS code requires you to specify which
equations are to be solved as the solution proceeds.
PHOENICS solves the Navier-Stokes equations for fluid
flow using the SIMPLEST scheme.
As such, the relevant velocity component for each
direction and the pressure are the minimum equations
to be solved.
The PIL statement has the Syntax
SOLVE(variable list)
Variables are separated by a comma, and are
P1 pressure
U1 velocity in X direction
V1 velocity in Y direction
cvl4='Please_enter_the_names_of_the_variables_to_be_solved'
cvl2=solve(p1,u1,v1)
CALL LOAD_95
The solution process for the PRESSURE correction can
be configured to be "whole field" or slabwise
Slabwise solution is the default, and is more suited
for thin shear (parabolic), or hyperbolic (supersonic)
flows
PIL syntax: SOLUTN(var,Y,Y,N,N,N,N)
where :
var is the variable name, and the six other parameters
store,solve,elliptical,jacobi,explicit,harmonic
For more information, see the encyclopaedia entry
on SOLUTN
Elliptical (whole field) pressure fields are normal
for general engineering flows, and are essential for
flows with recirculations/separations.
This case has recirculation, and hence requires
an elliptical pressure field.
eg: SOLUTN(P1,Y,Y,Y,N,N,N)
cvl4='Change_pressure_solver_to_be_whole_field'
cvl2=solutn(p1,y,y,y,n,n,n)
CALL LOAD_95
other stored variables
The PHOENICS code can solve/store up to 50 variables
Some of which we have already activated.
Storage of auxiliary variables and properties such as
density are required according to which problem is being
examined.
The variables stored can then be under-relaxed if
required to assist convergence, and are available for
plotting within PHOTON and AUTOPLOT
This case has no problems with convergence, but will
use the k-epsilon turbulence model, and the auxiliary
store for the turbulent viscosity (ENUT) is required
The PIL statement for the storage is :
STORE(variable list)
For more information, see Encyclopaedia entry on STORE
cvl4='Please_indicate_variables_to_be_stored'
cvl2=store(enut)
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 8. Terms (in differential equations) & devices
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 9. Properties of the medium (or media)
qcom( *
reference pressure
GROUP 8. Terms (in differential equations) & devices
GROUP 9. Properties of the medium (or media)
PHOENICS can model fluids with fluid properties that
vary,for example, density as a function of temperature
and pressure
The Auxiliary equations used for relationships quite
often require reference values to avoid divisions by
zero and to minimise rounding errors.
For this case, we do NOT vary the density, but as a
measure of "Good Practice" we should set a reference
pressure.
The PIL Syntax for setting a reference pressure is:
PRESS0=required real number
( reminder ref. pressure is 1.e5 )
cvl4='Please_set_the_reference_pressure'
cvl2=press0=1.e5
CALL LOAD_95
fluid density
PHOENICS has default values for typical variables
such as density, laminar viscosity etc.
As mentioned in the previous panel, it is possible to
specify relationships for these parameters.
For the purpose of this tutorial, the density is
constant, and its value must be specified.
The PIL statement :
RHO1=real value
will specify a constant density throughout the flow
domain.
(Reminder density=1.22)
cvl4='Please_set_the_fluid_density'
cvl2=rho1=1.22
CALL LOAD_95
reference laminar viscosity
For a similar reason, we must specify a value for the
Laminar viscosity
Constant values are specified via a PIL statement:
ENUL=real number
(Reminder viscosity of fluid is 1.465e-5)
cvl4='Please_indicate_the_reference_laminar_viscosity'
cvl2=enul=1.465e-5
CALL LOAD_95
activate turbulence models
PHOENICS can use several models of turbulence.
The PIL statement to activate a turbulence model is:
TURMOD(text)
where text is:
KLMODL - The Prandtl energy model
KEMODL - The Kolmorogov k-epsilon model
(which solves KE and EP)
Other models are available, such as
fixed viscosity (PIL: ENUT=required viscosity)
Laminar (PIL: ENUT=0.0)
cvl4='Please_activate_the_k-epsilon_turbulence_model'
cvl2=turmod(kemodl)
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 10. Inter-phase-transfer processes and properties
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 11. Initialization of variable or porosity fields
qcom( *
set initial fields
PHOENICS uses an iterative solution process to reach
the predicted flow.
As such, some guesses to the initial fields should be
supplied to the iterative solver.
All initial fields are set to 1.e-10 by default
If this value is inapropriate for any variable then
we can set these to suit. with the PIL statement
FIINIT(var)=required value
eg: FIINIT(U1)=100
sets uniform U1 velocity of 100 ms-1
cvl4='Please_set_initial_field_value_of_0.04_for_ke'
cvl2=fiinit(ke)=0.04
CALL LOAD_95
do jj=1,25
mesg(
enddo
It is normally good practice to set initial fields
for both ke and ep to minimise convergence problems
cvl4='Please_set_initial_field_value_of_0.13_for_epsilon'
cvl2=fiinit(ep)=0.13
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 12. Patchwise adjustment of terms
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 13. Boundary conditions and special sources
qcom( *
Having specified the geometry and the fluid properties
we must now indicate the locations of relevant boundaries
and what exactly is happening at each boundary.
This tutorial will introduce three types of boundary
condition
* INLETs where fluid will enter the domain
* OUTLETS where fluid leaves the domain
* Wall friction
There are numerous methods and types of boundary
specification the methods demonstrated in this tutorial
are only a small sub-set of the full range.
CALL LOAD_96
GROUP 10. Inter-phase-transfer processes and properties
GROUP 11. Initialization of variable or porosity fields
GROUP 12. Patchwise variation of terms
GROUP 13. Boundary conditions and special sources
Boundary conditions are located within the
computational volume by the pil statements
with the form
PATCH(name,type,ixf,ixl,iyf,iyl,izf,izl,itf,itl)
where:
name is a unique text identifier to link patches and
the values being set on the patch.
type is the relevant cell face area, volume etc.
ixf, etc are region indices if prefixed by a #,
or else CELL indices.
example:
INLET(IN,LOW,1,NX,#1,#2,1,1,#1,#1)
sets an inlet patch on the low cell faces for all x
cells in the first two regions in Y
(reminder Inlet is called INLET, is on the west side
of the first region in X, Y and Z)
cvl4='Please_set_name,_type_and_location_of_the_inlet'
cvl2=inlet(inlet,west,#1,#1,#1,#1,#1,#1,#1,#1)
CALL LOAD_95
Having specified the patch location, we now have to
specify what mass flows etc are present in the inlet
patch.
Mass flow rates are specified by adding a source to P1
the pressure correction equation.
The source within PHOENICS is linearised to have a
coefficient and value.
The magnitude of the coefficient is assumed for the
INLET style of PATCH, and the value must be mass per
unit area.
Other variables are just set to the required value
VALUE(name,variable,value)
where: name is used to link value to the patch.
Variable is the relevant equation eg. P1
Value is the required flux
(reminder the flow rate is 10.0*rho1; inlet velocity is 10.0)
cvl4='Please_set_mass_flow_rate_for_inlet_called_inlet'
cvl2=value(inlet,p1,10.0*rho1)
CALL LOAD_95
cvl4='Please_set_u_velocity_component_for_inlet_to_10.0'
cvl2=value(inlet,u1,10.0)
CALL LOAD_95
cvl4='Please_set_v_velocity_component_for_inlet_to_0.0'
cvl2=value(inlet,v1,0.0)
CALL LOAD_95
cvl4='Please_set_kinetic_energy_entrained_at_the_inlet_to_0.04'
cvl2=value(inlet,ke,0.04)
CALL LOAD_95
cvl4='Please_set_epsilon_being_entrained_at_the_inlet_to_0.13'
cvl2=value(inlet,ep,0.13)
CALL LOAD_95
congratulations....
you are really quite good at this now !!
Having specified the the inlet patch, we now have to
specify the outlet patch.
Regions of outflow are normally specified by regions
of nomimally constant pressure in a manner similar to
inlet above, but the patch name is OUTLET.
example:
OUTLET(OUT,HIGH,1,NX,#1,#2,NZ,NZ,#1,#1)
sets an outlet over all x, first two regions in Y last Z
for the first time region
(reminder outlet is called out and is on the east of the
domain, extending over region 1 in y and z)
cvl4='Please_set_name,type_and_location_of_the_outlet'
cvl2=outlet(out,east,#1,#1,#1,#1,#1,#1,#1,#1)
CALL LOAD_95
Again, on outlet boundaries, Mass flow rates etc are
linearised to have a coefficient and value.
At an outlet, only the relative pressure need be set -
all other values are internal, and hence known to EARTH.
Required pressure values are set using :
VALUE(name,P1,value)
where:
name is used to link value to the patch.
value is the external pressure, relative to PRESS0
cvl4='Please_set_nominal_pressure_on_outlet_boundary_to_0.0'
cvl2=value(out,p1,0.0)
CALL LOAD_95
Having specified the inlet and exit patches, we now have to
specify wall friction.
Regions with wall friction can be specified using a PATCH.
example:
PATCH(TOP,NWALL,1,NX,#1,#1,1,NZ,#1,#1)
sets a north wall over all x, last cell in Y, all Z
for the first time region
the type is NWALL, SWALL EWALL,WWALL,HWALL or LWALL
as required, but ONLY for wall regions.
These are special names used to identify regions
where the distance from wall to cell centre is important
(reminder, wall is called TOP and is on the North
part of the domain, extending in all x-regions)
cvl4='Please_set_name_type_and_location_of_the_top_wall'
cvl2=patch(top,nwall,#1,#1,#1,#1,#1,#1,#1,#1)
CALL LOAD_95
The wall condition must be specified in the
linearised form having a coefficient and value.
The magnitudes of the coefficients are NOT KNOWN for the
formal PATCH statement, and we use COVAL statements to
specify the required data
COVAL(name,variable,COefficient,VALue)
where: name is used to link value to the patch.
Variable is the relevant equation eg. P1
Coeffficient can be set as follows:
Coef, influence on prediction
blasius Blasius law
loglaw Log law
genlaw Generalised log law
Value
Velocities Wall speed.
ke-ep loglaw or genlaw Bounds the Ke and ep equations
For this wall the log law method is to be activated.
cvl4='Please_set_the_coval_for_the_velocity_parallel_to_the_wall'
cvl2=coval(top,u1,loglaw,0.0)
CALL LOAD_95
cvl4='Please_set_the_coval_for_the_kinetic_energy_boundary'
cvl2=coval(top,ke,loglaw,loglaw)
CALL LOAD_95
cvl4='Please_set_the_coval_for_the_dissipation_boundary'
cvl2=coval(top,ep,loglaw,loglaw)
CALL LOAD_95
CONGRATULATIONS.
You have set the three types of Boundary to be covered
in this tutorial.
All other boundaries ( one wall in this case ) are being
set for you.
CALL LOAD_96
dcom(patch(bot,swall,#1,#1,#1,#1,#1,#1,#1,#1)
dcom(coval(bot,u1,loglaw,0.0)
dcom(coval(bot,ke,loglaw,loglaw)
dcom(coval(bot,ep,loglaw,loglaw)
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 14. Downstream pressure for PARAB=.TRUE.
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 15. Termination of sweeps
qcom( *
GROUP 14. Downstream pressure for PARAB=.TRUE.
GROUP 15. Termination of sweeps
PHOENICS is an iterative code that proceeds from your
first guess to a converged prediction.
For general engineering problems, the "whole field"
pressure is likely to be used, and the iterative
process can be described as follows:
1) A series of "sweeps" through the domain from IZ = 1
to IZ = NZ.
2) Each Z "slab" is treated in turn to derive velocity
updates to satisfy momentum, and the pressure
coefficients for the 3D solver
3) The 3D pressure field is solved, and then used to
derive velocities that satisfy continuity.
For this tutorial, we are only concerned with the
overall number of SWEEPS, with the PIL statement
LSWEEP=n
where n is the integer number of sweeps required.
cvl4='Set_the_iteration_sweeps_to_50'
cvl2=lsweep=50
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 16. Termination of iterations
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 17. Under-relaxation devices
qcom( *
dcom(real(maxv,minl,relx)
dcom(maxv=5.0
dcom(minl=0.1
dcom(relx=10.0
dcom(relax(p1,linrlx,0.8)
dcom(relax(u1,falsdt,minl/maxv*relx)
dcom(relax(v1,falsdt,minl/maxv*relx)
dcom(relax(ke,falsdt,minl/maxv*relx)
dcom(relax(ep,falsdt,minl/maxv*relx)
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 18. Limits on variables or increments to them
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 19. Data communicated by satellite to GROUND
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 20. Preliminary print-out
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 21. Print-out of variables
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 22. Spot-value print-out
qcom( *
dcom(ixmon=5
dcom(iymon=5
dcom(tstswp=-1
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 23. Field print-out and plot control
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 24. Dumps for restarts
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
You have now completed the tutorial.
You may now end satellite and write out the Q1 file, or
quit this tutorial.
'Enter' will end the satellite and save the Q1.
'Quit' will quit the tutorial.
IMPORTANT!!
If you are running from the satellite top menu, and you
want to save the Q1 you have generated, you MUST select
'Command mode', NOT 'End satellite' when you return to
the top panel, otherwise you will lose your Q1.
CALL LOAD_96
goto DONE
label QUIT
CALL LOAD_97
label DONE
ENDMAIN
SUBROUTINE LOAD_95
L($095)
ENDSUB
SUBROUTINE LOAD_96
L($096)
ENDSUB
SUBROUTINE LOAD_97
L($097)
ENDSUB