TALK=T;RUN(1,1) 093

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