label TOP
q1quit=t;dmpstk=f;nocomm=t;nocopy=t;nowipe=t
rset(v,0)
namsat=chkc
boolean(log1)
char(cvl1,cvl2,cvl3,cvl4)
cvl3='Correct'
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 1. Run title and other preliminaries
qcom( *
display
WELCOME TO PIL TUTORIAL NUMBER 3
********************************
This tutorial continues to teach the basic PIL commands
and specifies a simple problem of flow through the simple
labyrinth, shown below.
The Solution Geometry (Cartesian)
_____________________________________________________
|INLET |BLK2| OUTLET |
|m/a=5*rho1 Type cell P1=0.0 |
|U1=5.0;V1=0.0 Porosity=0.0 |
|KE=0.005;EP=1.536E-2 wall fric |
|____________ | | ____________|
| | |____| | |
|BLK1 | |BLK3 |
|Type=cell | |Type=cell |
|Porosity=0.0| |Porosity=0.0|
|Wall fric | |Wall fric. |
|____________|___________________________|____________|
CALL LOAD_96
You will be guided through the steps of setting up the
problem and advised about what PIL commands you should
type.
If you do not type them precisely as required, you will
be invited to try again; thereafter the expected response
will be shown on the screen.
This will happen, for example, if you type .1 where 0.1
is expected, even though both would be interpreted in the
same way be SATELLITE. So the tutorial program is stricter
than it need be.
You can quit the tutorial at any time by typing "quit".
CALL LOAD_96
* Start of Tutorial 3 *
PIL command: TEXT
-----------------
Group 1 of the Q1 file is used to specify the title
of the prediction.
The PIL command to set the title has the form:
TEXT(Required Title
The required text in this case is tutorial.3
cvl2=text(Tutorial.3
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( *
The Solution Mesh
The solution geometry and mesh now has to be defined.
There are several methods available within PHOENICS to
define the computational mesh. This Tutorial will
exemplify the method of dividing the X and Y direction
into REGIONS to accommodate the obstacles and boundary
conditions.
Each region is further divided into cells. The region and
cell distribution is tabulated on the next display.
CALL LOAD_96
Number of REGIONS in X-direction = 5
Number of REGIONS in Y-direction = 3
___________________________________________
| X region | 1 | 2 | 3 | 4 | 5 |
| No. of cells| 7 | 4 | 3 | 4 | 7 |
| Length | 0.75| 0.4 | 0.2 | 0.4 | 0.75|
| power |-1.4 | 1.0 | 1.0 | 1.0 | 1.4 |
|_____________|_____|_____|_____|_____|_____|
_______________________________
| Y region | 1 | 2 | 3 |
| No. of cells| 7 | 3 | 5 |
| Length | 0.75| 0.25| 0.5 |
| power | 1.0 | 1.0 | 1.0 |
|_____________|_____|_____|_____|
CALL LOAD_96
PIL variable: NREGX
-------------------
The first step in setting the mesh for this tutorial is
to define the number of regions in the X direction.
Then region identitifiers are 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_5'
cvl2=nregx=5
CALL LOAD_95
PIL variable: IREGX
-------------------
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
PIL command: GRDPWR
-------------------
Now that we have specified the region, the cell spacing
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, and
PWR is the power that is used to distribute the cells.
For example, GRDPWR(x,10,100,1.4) would position 10
cells in x over 100 metres, with small cells at the
beginning and large cells at the end.
A negative power will put large cells at the beginning.
Our region 1 has 7 cells, is 0.75m long, and the
power for cell distribution is -1.4.
cvl4='Please_enter_the_distribution_of_cells_in_x_region_1'
cvl2=grdpwr(x,7,0.75,-1.4)
CALL LOAD_95
PIL variable: IREGX
-------------------
Having specified the cell distribution in X region 1,
we now have to specify the cell spacing in region 2.
Again, to do this requires two operations:
1) to identify the region 2.
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, and 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 powers expand/contract cells
cvl4='Please_enter_the_index_of_the_second_x_region_to_be_modified'
cvl2=iregx=2
CALL LOAD_95
PIL command: GRDPWR(...)
------------------------
Region 2 has 4 cells of equal size, and is 0.4 m long.
cvl4='Please_enter_the_distribution_of_cells_in_x_region_2'
cvl2=grdpwr(x,4,0.4,1.0)
CALL LOAD_95
PIL variable: IREGX
-------------------
Having specified the cell distribution in X region 2
we now have to specify the cell spacing in region 3
Again, To do this requires two operations:
1) to identify region 3.
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 identity
IC is the number of cells in the region
LENGTH is the cartesian extent of the region
positive/negative power expands/contracts cells
cvl4='Please_enter_the_index_of_the_third_x_region_to_be_modified'
cvl2=iregx=3
CALL LOAD_95
PIL command GRDPWR(...)
-----------------------
Region 3 is 0.2 long has 3 cells, uniformly distributed
cvl4='Please_enter_the_distribution_of_cells_in_x_region_3'
cvl2=grdpwr(x,3,0.2,1.0)
CALL LOAD_95
PIL variable IRGEX & PIL command GRDPWR
---------------------------------------
Having specified the cell distribution in X region 3,
we now have to specify the cell spacing in region 4.
Again, To do this requires two operations:
1) to identify region 4.
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 identity
IC is the number of cells in the region
LENGTH is the cartesian extent of the region
positive/negative power expands/contracts cells
cvl4='Please_enter_the_index_of_the_x_region_to_be_modified'
cvl2=iregx=4
CALL LOAD_95
PIL command GRDPWR(...)
-----------------------
Region 4 has 4 cells uniformly distributed
and is 0.4 long
cvl4='Please_enter_the_distribution_of_cells_in_x_region_4'
cvl2=grdpwr(x,4,0.4,1.0)
CALL LOAD_95
PIL variable IREGX & PIL command GRDPWR
---------------------------------------
Having specified the cell distribution in X region 4,
we now have to specify the cell spacing in region 5.
Again, To do this requires two operations:
1) to identify region 5.
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 identity
IC is the number of cells in the region
LENGTH is the cartesian extent of the region
positive/negative power expands/contracts cells
cvl4='Please_enter_the_index_of_the_x_region_to_be_modified'
cvl2=iregx=5
CALL LOAD_95
PIL command GRDPWR(...)
-----------------------
Region 5 has 7 cells expanded with power 1.4,
and is 0.75 long.
cvl4='Please_enter_the_distribution_of_cells_in_x_region_5'
cvl2=grdpwr(x,7,0.75,1.4)
CALL LOAD_95
Congratulations,
You have now set the X direction mesh distribution.
This mesh will now be shown by use of the view
statement: (view(k,1) .
This is already in the tutorial.q1, so there is no need
to type it.
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( *
The Y-direction grid is being set for you
dcom(NREGY=3
dcom(iregy=1
dcom(grdpwr(y,7,0.75,1.0)
dcom(iregy=2
dcom(grdpwr(y,3,0.2,1.0)
dcom(iregy=3
dcom(grdpwr(y,5,0.5,1.0)
The mesh specification for the tutorial has been
completed.
This mesh will now be shown by use of the view
statement.
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 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( *
PIL command SOLVE(var)
----------------------
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 for this tutorial
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
PIL command SOLUTN(var,...)
---------------------------
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
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
PIL command STORE(var)
----------------------
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 few 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)
cvl4='Please_indicate_which_other_variables_are_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
PIL Variable PRESS0
-------------------
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.
In this case, we do NOT vary the density, but as a
measure of "Good Practice" let us set a reference
pressure, to 1.e5.
The PIL Syntax for setting a reference pressure is:
PRESS0=required real number
cvl4='Please_set_the_reference_pressure'
cvl2=press0=1.e5
CALL LOAD_95
Fluid Density
PIL variable RHO1
-----------------
PHOENICS has default values for typical variables
such as density, laminar viscosity etc.
As mentioned in the previous panel, it is possible to
activate formulae for these parameters.
For the purpose of this tutorial, the density is
constant, and its value, 1.22, must be specified.
The PIL statement :
RHO1=real value
Will specify a constant density throughout the flow
domain.
cvl4='Please_set_the_fluid_density'
cvl2=rho1=1.22
CALL LOAD_95
Laminar Viscosity
PIL variable ENUL
-----------------
For a similar reason, we must specify a value for the
Laminar viscosity. In this case it is 1.465e-5
Constant values are specified via a PIL statement:
ENUL=real number
cvl4='Please_set_the_laminar_viscosity'
cvl2=enul=1.465e-5
CALL LOAD_95
Activate Turbulence Models
PIL command TURMOD(...)
-----------------------
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 Blockages via Conpor Statements
PIL command CONPOR
------------------
In this tutorial we have three blockages.
All are type CELL, and have porosity values 0.0. The size
of each blockage is relayed via region numbers. The name
and regions for eack blockage are tabulate below. Please
make a note of this.
______________________________
| Name | X-regions | Y-regions |
| BLK1 | 1 | 1 & 2 |
| BLK2 | 3 | 2 & 3 |
| BLK3 | 5 | 1 & 2 |
|______|___________|___________|
CALL LOAD_96
GROUP 11. Initialization of variable / porosity fields
Within the structured grid environment, solids can be
represented via blocked regions of the mesh.
These blocked regions are identified by the following:
CONPOR(name,value,type,ixf,ixl,iyf,iyl,izf,izl)
Where:
Name is the blocked region identifier eg: BLK1
Value is a number from 0.0 (solid) to 1.0 (free flow).
Type is a character string denoting which cell face or
cell volume is to be affected by the conpor statement.
ixf and ixl are the start and end CELL indices of the
blockage in the I direction.
Similarly for iyf/iyl, and izf/izl.
If ixf etc. are prefixed by a #, then the blockage
indices are by REGION, and, if further prefixed by a
-ve sign, then wall friction is activated
eg. CONPOR(TB1,0.0,CELL,1,NX,-#2,-#2,#1,#NREGZ)
cvl4='Please_set_BLK1_porosity_statement_using_region_identities'
cvl2=conpor(blk1,0.0,cell,-#1,-#1,-#1,-#2,-#1,-#1)
CALL LOAD_95
We now have to set the blockage for blk2
Remember :
CONPOR(name,value,type,ixf,ixl,iyf,iyl,izf,izl)
Where:
Name is the blocked region identifier eg: BLK1
Value is a number from 0.0 (solid) to 1.0 (free flow).
Type is a character string denoting which cell face or
cell volume is to be affected by the conpor statement.
ixf and ixl are the start and end CELL indices of the
blockage in the I direction.
Similarly for iyf/iyl, and izf/izl.
If ixf etc. are prefixed by a #, then the blockage
indices are by REGION, and, if further prefixed by a
-ve sign, then wall friction is activated
eg. CONPOR(TB1,0.0,CELL,1,NX,-#2,-#2,#1,#NREGZ)
cvl4='Please_set_BLK2_porosity_statement_using_region_identities'
cvl2=conpor(blk2,0.0,cell,-#3,-#3,-#2,-#3,-#1,-#1)
CALL LOAD_95
We now have to set the blockage for blk3
Remember :
CONPOR(name,value,type,ixf,ixl,iyf,iyl,izf,izl)
Where:
Name is the blocked region identifier eg: BLK1
Value is a number from 0.0 (solid) to 1.0 (free flow).
Type is a character string denoting which cell face or
cell volume is to be affected by the conpor statement.
ixf and ixl are the start and end CELL indices of the
blockage in the I direction.
Similarly for iyf/iyl, and izf/izl.
If ixf etc. are prefixed by a #, then the blockage
indices are by REGION, and, if further prefixed by a
-ve sign, then wall friction is activated
eg. CONPOR(TB1,0.0,CELL,1,NX,-#2,-#2,#1,#NREGZ)
cvl4='Please_set_BLK3_porosity_statement_using_region_identities'
cvl2=conpor(blk3,0.0,cell,-#5,-#5,-#1,-#2,-#1,-#1)
CALL LOAD_95
Congratulations,
You have now set the blockages.
These will now be shown by use of the view statement.
please press "enter" to view the geometry thus far
and then press "enter" again to continue the tutorial
CALL LOAD_96
gclear
gview(p,0,0,1)
gdom(1,nx+1,1,ny+1,1,1,1,0)
gpatch(blk1,6,0);gpatch(blk2,6,0);gpatch(blk3,6,0)
gdraw
Set Initial Fields
PIL command FIINIT(var)
-----------------------
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_for_ke_of_0.005'
cvl2=fiinit(ke)=0.005
CALL LOAD_95
It is normally good practice to set initial fields
for both ke and ep to minimise convergence problems
cvl4='Please_set_initial_field_value_for_epsilon_to_1.536e-2'
cvl2=fiinit(ep)=1.536e-2
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 12. Patchwise adjustment of terms
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 13. Boundary conditions and special sources
qcom( *
Boundary Conditions
PIL commands INLET(...), VALUE(...)
-----------------------------------
The boundary conditions to be set are:
(i) inflow
(ii) outflow
(iii) walls
CALL LOAD_96
GROUP 13. Boundary conditions and special sources
Boundary conditions are located within the
computational volume by the pil statements
of 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. (-ve signs NOT allowed)
example:
INLET(IN,LOW,1,NX,#1,#2,1,1,#1,#1)
sets a fixed mass flow rate patch on the low cell faces
for all x cells in the first two regions in Y.
The mass flow rate boundary condition is called INLET,
and is located on WEST of region 1 of x and 3 of y)
cvl4='Please_set_name_type_and_location_of_the_inlet'
cvl2=inlet(inlet,west,#1,#1,#3,#3,#1,#1,#1,#1)
CALL LOAD_95
PIL command VALUE(...)
----------------------
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
cvl4='Please_set_mass_flow_rate_of_5*rho1_for_inlet'
cvl2=value(inlet,p1,5*rho1)
CALL LOAD_95
cvl4='Please_set_u_velocity_component_to_5.0_for_inlet'
cvl2=value(inlet,u1,5.0)
CALL LOAD_95
cvl4='Please_set_v_velocity_component_to_0.0_for_inlet'
cvl2=value(inlet,v1,0.0)
CALL LOAD_95
cvl4='Please_set_kinetic_energy_entrained_at_the_inlet_to_0.005'
cvl2=value(inlet,ke,0.005)
CALL LOAD_95
cvl4='Please_set_epsilon_being_entrained_at_the_inlet_to_1.536e-2'
cvl2=value(inlet,ep,1.536e-2)
CALL LOAD_95
Congratulations....setting all other bconds for you
you are really quite good at this now !!
dcom(outlet(out,east,#5,#5,#3,#3,#1,#1,#1,#1)
dcom(patch(top,nwall,#1,#5,#3,#3,#1,#1,#1,#1)
dcom(coval(top,u1,loglaw,0.0)
dcom(coval(top,ke,loglaw,loglaw)
dcom(coval(top,ep,loglaw,loglaw)
dcom(patch(bot,swall,#2,#4,#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)
Congratulations,
You have now set the boundary conditions.
These will now be shown by use of the view statement.
please press "enter" to view the geometry thus far
and then press "enter" again to continue the tutorial
CALL LOAD_96
gpatch(inlet,14,0);gpatch(out,4,0);gpatch(top,12,0);gpatch(bot,12,0)
gdraw
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 14. Downstream pressure for PARAB=.TRUE.
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 15. Termination of sweeps
qcom( *
GROUP 15. Termination of sweeps
PIL variable LSWEEP
-------------------
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='20_iterative_sweeps_are_required_please_set_this'
cvl2=lsweep=20
CALL LOAD_95
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 16. Termination of iterations
qcom( *
qcom( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
qcom( *
qcom( * GROUP 17. Under-relaxation devices
qcom( *
GROUP 17. Under-relaxation devices
PIL Commands REAL(vars) & RELAX(...)
------------------------------------
The iterative process can be unstable and can in some
circumstances be divergent.
To minimise the chance of instabilities and divergence
we can "slow down" or "damp" the update to any or all
of the solved or stored variables.
Two types of relaxation are available:
1) Linear relaxation for scalars such as P1,RHO1 etc
2) False time step relaxation for transport equations
such as U1,V1,W1,H1,KE,EP etc.
This value is typically:
minimum cell dimension/max velocity
The under-relaxation values for each equation are
often related to each other, and it is easier to
declare some local variables that can then be used
within the Q1 file.
PIL syntax : REAL(var1,var2,var3)
cvl4='Declare_three_real_values_and_call_them_maxv_minl_and_relx'
cvl2=real(maxv,minl,relx)
CALL LOAD_95
The three variables defined in the previous step are
to be used to define the false time step from the
following relationship
TF=minl/maxv*relx
MINL is the smallest cell size,
MAXV is the maximum velocity
RELX is a multiplier to globally alter the relaxation
Having defined the variables,they must be initialized
Maxv should be set to the estimated maximum
velocity within the flow domain
PIL syntax
variable name=value
cvl4='Set_the_value_of_maxv_to_5.0'
cvl2=maxv=5.0
CALL LOAD_95
MINL should be set to the smallest cell size
within the flow domain
PIL syntax
variable name=value
cvl4='Set_the_value_of_MINL_to_0.1'
cvl2=minl=0.1
CALL LOAD_95
RELX should be set to the overall relaxtion factor
100.0=weak relaxation,
1.0=default relaxation
0.1=strong relaxation
PIL syntax
variable name=value
cvl4='Set_the_value_of_relx_to_the_default_value_of_1.0'
cvl2=relx=1.0
CALL LOAD_95
RELAXATION on the PRESSURE correction equation
The relaxation applied to the pressure equation is
Linear relaxation.
The PIL statement to apply linear relaxation is:
RELAX(variable,LINRLX,value)
eg.
RELAX(P1,LINRLX,0.1)
would apply only one tenth of the predicted update per
iterative sweep, which is very heavily relaxed.
This tutorial should only require a factor of 0.8 to
apply 80% of the predicted update per sweep.
cvl4='Set_the_linear_under-relaxation_for_the_pressure_to_0.8'
cvl2=relax(p1,linrlx,0.8)
CALL LOAD_95
RELAXATION on the U1 velocity
The relaxation on the transport equation for velocity
is of the false time step type.
The PIL syntax:
RELAX(variable,FALSDT,value)
eg.
RELAX(U1,FALSDT,.001)
or
RELAX(U1,FALSDT,function)
In this tutorial we wish to set a function:
minl/maxv*relx
cvl4='Set_the_under-relaxation_for_the_u1_velocity_to_the_function'
cvl2=relax(u1,falsdt,minl/maxv*relx)
CALL LOAD_95
RELAXATION on the V1 velocity
V1 is also a transport equation, and will require
false time step relaxation.
It is possible to apply different values of relaxation
to each separate velocity component, but for this
tutorial we apply equal relaxation to the velocities.
Hence, we can set the function to
minl/maxv*relx
remember PIL syntax:
RELAX(var,FALSDT,function)
cvl4='Set_the_under-relaxation_for_the_v1_velocity'
cvl2=relax(v1,falsdt,minl/maxv*relx)
CALL LOAD_95
RELAXATION on the TURBULENCE model
KE and EP are also a transport equations, and will
require false time step relaxation.
Since KE and EP are functions of velocity, transported
by the velocity, and are then used to derive source
terms for the velocity, we under-relax the equations
for KE and EP in the same way as the velocity equations.
Hence, we can set the function to
minl/maxv*relx
remember PIL syntax:
RELAX(var,FALSDT,function)
cvl4='Set_the_under-relaxation_for_the_turbulent_kinetic_energy'
cvl2=relax(ke,falsdt,minl/maxv*relx)
CALL LOAD_95
well done.....setting ep relaxation
to: relax(ep,falsdt,minl/maxv*relx)
dcom(relax(ep,falsdt,minl/maxv*relx)
CALL LOAD_96
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( *
Monitoring Location
PIL variables IXMON & IYMON
---------------------------
GROUP 22. Spot-value print-out
It is possible to ask for PHOENICS to print the values
of the SOLVED variables that are present in a monitor
cell.
The location of the monitor cell is indicated by
up to three PIL statements, separated by ;
eg.
ixmon=10;iymon=30;izmon=50
For the tutorial, we only will require values
ixmon and iymon
cvl4='Please_set_the_monitor_location_to_14_in_x_and_3_in_y'
cvl2=ixmon=14;iymon=3
CALL LOAD_95
Error Print-out
PIL variable TSTSWP
-------------------
The error in the prediction can be written to the
screen at the end of every n sweeps. Where n is any
integer.
If the integer is positive, the residual data is
tabluated numbers.
If the integer is negative, both the spot value and
residual outputs are presented graphically.
eg.
TSTSWP=-5
would output graphical data every 5 sweeps.
As a general rule, the error (or residual) should
decrease as the run proceeds.
cvl4='Please_set_the_residual_output_to_be_graphical_every_sweep'
cvl2=tstswp=-1
CALL LOAD_95
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
enddis
ENDMAIN
SUBROUTINE LOAD_95
L($095)
ENDSUB
SUBROUTINE LOAD_96
L($096)
ENDSUB
SUBROUTINE LOAD_97
L($097)
ENDSUB