pq1t2_2.htm

Tutorial pq1t2_1; Part 2
Library case 621: Simulation of a Labyrinth Flow.

Preliminary remarks

In this tutorial, which is the second part of the tutorial accessed by clicking here, you will learn how to:

express the results of calculation in non-dimensional form;
parameterize the computational-grid distribution;
enable its parameters to be read from a separate file; and
thereby study the influence of grid refinement on the accuracy of the results.

It is accordingly provided with the following hyperlinks to the corresponding sections:

Non-dimensional form of results
Grid parameterization;
Parameters read from file; and
How grid refinement affects accuracy.
Concluding remarks.

1. Non-dimensional form of results

Reason for use

It is often useful to represent the results of computations in non-dimensional form, so that:

expectations as to their values can be formulated;
departures from expectations can be investigated and their causes determined; or, alternatively
confidence in the reliability of the software can be enjoyed when agreement shows it to be justified.

In the present case it is reasonable to consider three non-dimensional quantities, namely:

non-dimensional velocity, nd_v, defined as
vabs/invel,
where
vabs stands for the local absolute velocity, varying from place to place, and
invel is the inlet velocity;
non-dimensional pressure, nd_p, defined as
(p1-pext)/(0.5*rho1*invel**2)
where
p1 stands for the local pressure, varying from place to place,
pext is the external pressure, i.e. that of the outlet,
rho1 is (the PIL-variable name for) the density of the fluid, presumed uniform,
and
0.5*rho1*invel**2 is what is commonly known as the 'velocity head' of the inlet stream;
non-dimensional temperature increase, nd_t, defined as
(tem1-intem)*cp1*inmass/heatflux,
where:
tem1 stands for the local temperature, varying from place to place,
intem is the inlet temperature, a constant,
cp1 is (the PIL-variable name for) the specific heat of the fluid, presumed uniform,
heatflux is the heat input,
inmass represents the inlet mass flow rate calculated as inarea*rho1*invel where inarea is the inlet cross-sectional area.

Here it should be recalled that, as explained at the start of part 1 of this tutorial, all material properties (i.e. density, viscosity, specific heat, etc.) are presumed to be uniform.

Of the computed distributions of these quantities, we can reasonably expect:

that nd_v will be close to unity near the inlet, will fall below it in the larger volumes of the labyrinth, but will possibly rise above it in the gaps below IN-PLATE and above IN-BLOCK;
that, because of the presence of obstacles and of friction at walls, nd_p is likely to be above zero in most places, and largest near the inlet; but low values may appear in small-cross-sectional-area regions, where nd_vel is largest;
[One might suppose (wrongly as it will turn out) that -1.0 must be the lowest possible value, arguing that this would be a consequence of the Bernoulli equation if the duct were frictionless and the outlet area very great. We shall see.]
that nd_t will everywhere exceed zero, because there is a heat source, but no heat sink;
that nd_t will have an average value of unity at the outlet; and indeed the value would be exactly unity there if the turbulent mixing caused by the obstacles were made sufficiently intense by narrowing the gaps between them; and
that higher-than-unity values may be found within the domain nearest to the heat source in the low-z plate, because of the finite value of its surface heat-transfer coefficient.

How to introduce the non-dimensional variables

Let us first of all create a separate directory for this run in a way already described in Part 1 of the tutorial 2. To refresh your memory of this, click here.

Place there the q1 you already used in your previous runs.
We are now going to create a proper data file for subsequent runs based on the gp25_1 file. So copy this file into your new directory with help of the Commander editor.
However, if you do not want to repeat all these steps, you can copy the data file gp25_3 from the folder phoenics/d_polis/d_tuts/pq1ts/pq1t2_2 into your present working directory.

The non-dimensional variables described above are computed by PHOENICS, if the following changes are introduced into the appropriate places in the gp25_1 file of your new working directory, which file we now rename to gp25_3.

First new real variables are introduced and specified.

We shall repeat here all necessary statements to enable the use of copy and paste commands, thus avoiding any possible error.
```
real(real(invel,inmass,intem,pext,inarea,heatflux)
intem=20.0
pext=0.0
heatflux=100.
invel=0.05
inarea=yvlast*0.45*zwlast
inmass=inarea*rho1*invel
inarea
inmass
```
These are inlet temperature, velocity, cross-sectional area; and mass-flow rate, heat input and external pressure.
Next comes the group of statements for storing and solving field variables expressed via formulae. To read more about these, turn to the Encyclopaedia article about In-Form by clicking here.
```
store(vabs)
(stored var nd_v is vabs/:invel:) 
(stored var nd_p is (p1-:pext:)/(.5*:rho1:*:invel:^2)) 
(stored var nd_t is (tem1-:intem:)*:cp1:*:inmass:/:heatflux:)
mesg(variables stored)
stored
```
It is in this section that the formulae for non-dimensional variables explained earlier are introduced.
The next change concerns the object L-BLOCK, for which we replace the previous numerical value of the heat flux by a corresponding parameter:
```
> OBJ,    NAME,        L-BLOCK                  
> OBJ,    HEAT_FLUX,    0.000000E+00, heatflux
```

Similar changes should be made for the INLET object

> OBJ,    TYPE,        INLET      
> OBJ,    PRESSURE,     0.000000E+00
> OBJ,    VELOCITY,     invel, 0.000000E+00, 0.000000E+00
> OBJ,    TEMPERATURE,  intem

for which the values of velocity and temperature are replaced by the corresponding parameters, and

for the OUTLET object in respect of pressure.

> OBJ,    TYPE,        OUTLET
> OBJ,    PRESSURE,     pext

Discussion of results

Runs performed with the above settings lead to results which will be displayed and discussed as follows.

Non-dimensional velocity nd_v, as expected:
- is close to unity near the inlet,
- falls below unity in the larger volumes of the labyrinth; and
- rises above unity in gaps formed below the IN-PLATE, between the H-BLOCK and IN-BLOCK; and also
- in the vicinity of the outlet, because of the reduced cross-sectional area for the flow, and because the velocity has a vertical-direction component as well as a horizontal one.
Non-dimensional pressure nd_p

is largest at the inlet, and smallest in gaps where the flow velocity is high as has been expected.
Its minimum value is about -2.9 in the gap between the H-BLOCK and the IN-BLOCK.
This is contrary to the above supposition that -1.0 is its lowest possible value. Is the supposition wrong? Or is the grid too coarse for accuracy?
The studies to be made in section 4 will throw light on these questions.
Non-dimensional temperature nd_t

is close to zero the upper part of the labyrinth, which is farther from the heated L-BLOCK.
As has been expected, it is around unity at the outlet, although its average value cannot be deduced by eye.

It might be concluded that displays of non-dimensional variables are more informative than those of their dimensional equivalents.
However, for temperature this is not quite true; for if we seek to display its distribution inside the solids, by pressing on the wireframe toggle

, we shall get the following picture:

Although its values must exist there, they are not displayed.

[The reason is that the Viewer has been programmed to display in solid regions only a few pre-specified variables, namely temperature, PRPS and stresses and strains].

The contours of dimensional temperature below:

show that within the solids L-BLOCK and IN-BLOCK the temperature is practically uniform. This is to be explained by the high conductivities of the materials of which they are made, namely copper and aluminium.

It should be remarked that, because the default "averaged" contour option has been used for creation of the above picture, low temperatures do appear to penetrate into the block.
If however the option is switched off, as is shown in the Viewer Options window accessed by clicking on the variable button ,

a more realistic display appears, although the coarseness of the grid is also thereby revealed:

Were the IN-BLOCK to be made of glass, for example, with a conductivity of 0.002 times that of copper, the temperature distribution would be quite different: the heated lower block (L-BLOCK) would remain hot, but the upper parts of IN-BLOCK would be cooler, as is shown (with 'averaged' re-selected) here:

This result can be achieved by replacing the following line in the data file

> OBJ,    MATERIAL,    103,COPPER at 27 deg C

> OBJ,    MATERIAL,    106,GLASS

and making a new run.

[The remainder of the tutorial, however, will concern the IN-BLOCK made of copper; so do not forget to restore the previous material setting in the data file gp25_3 before proceeding.]

Concluding remarks for section 2.

The results of the simulation have coincided, more or less, with our expectations; but a doubt has been raised as to whether an unexpectedly low nd_p was the effect of grid coarseness.

It is wise therefore to explore the extent which the result depends upon the fineness of the grid; this will now be done.

2. Grid parameterization

General remarks about the grid

The computational grid which has been employed up to his point has been that which the Satellite has chosen for itself, guided only by the lines in the Q1 file:

 Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
 RSET(M,20,1,15)

together with the fact, that no instruction to the contrary having been given, the grid should fit the objects within the default tolerance.

We shall now inspect that grid and examine the decisions which the Satellite made. Thereafter we shall consider how we can over-rule these decisions by setting parameters.

Inspecting the grid

Viewed in the VR-Editor if you click on the 'Mesh toggle' button

, the grid appears thus:

However, it might be easier to examine the grid if the following procedures are made.

First activate the head-on, i.e. front view on the labyrinth objects by pressing the reset button on the Movement panel.
In the 'Reset View Parameters' window click on the 'Nearest Head-on' button and then on 'OK'.
Then convert your present 3-D image to a 2-D one. To achieve this, go to 'Settings' and select the 'Depth Effect' command. In the 'Depth Effect' window either type a greater value in the box and then click anywhere on the screen, or move the slider to the right thus increasing the distance between the observer's eyes and the labyrinth objects.

You will then see a picture like this one.

The IN-PLATE being of zero thickness has disappeared and now only the IN-BLOCK and the grid are visible.

The grid is uniformly distributed in the x-direction, but not quite so in the z-direction, as is confirmed by inspecting what is printed by the Satellite in the file q1ear (and later by EARTH in RESULT), namely:

  Group 3. X-Direction Grid Spacing
 CARTES = T
 NX = 20
 XULAST  = 2.000000E+00
 XFRAC (   1) =  5.000000E-02 ;XFRAC (   2) =  1.000000E-01
 XFRAC (   3) =  1.500000E-01 ;XFRAC (   4) =  2.000000E-01
 XFRAC (   5) =  2.500000E-01 ;XFRAC (   6) =  3.000000E-01
 XFRAC (   7) =  3.500000E-01 ;XFRAC (   8) =  4.000000E-01
 XFRAC (   9) =  4.500000E-01 ;XFRAC (  10) =  5.000000E-01
 XFRAC (  11) =  5.500000E-01 ;XFRAC (  12) =  6.000000E-01
 XFRAC (  13) =  6.500000E-01 ;XFRAC (  14) =  7.000000E-01
 XFRAC (  15) =  7.500000E-01 ;XFRAC (  16) =  8.000000E-01
 XFRAC (  17) =  8.500000E-01 ;XFRAC (  18) =  9.000000E-01
 XFRAC (  19) =  9.500000E-01 ;XFRAC (  20) =  1.000000E+00
  ************************************************************
  Group 5. Z-Direction Grid Spacing
 PARAB   =    F
 NZ      =        15
 ZWLAST  = 1.000000E+00
 ZFRAC (   1) =  5.000000E-02 ;ZFRAC (   2) =  1.125000E-01
 ZFRAC (   3) =  1.750000E-01 ;ZFRAC (   4) =  2.375000E-01
 ZFRAC (   5) =  3.000000E-01 ;ZFRAC (   6) =  3.666667E-01
 ZFRAC (   7) =  4.333333E-01 ;ZFRAC (   8) =  5.000000E-01
 ZFRAC (   9) =  5.666667E-01 ;ZFRAC (  10) =  6.333333E-01
 ZFRAC (  11) =  7.000000E-01 ;ZFRAC (  12) =  7.833333E-01
 ZFRAC (  13) =  8.666667E-01 ;ZFRAC (  14) =  9.500000E-01
 ZFRAC (  15) =  1.000000E+00

The uniformity in the x-direction is a consequence of the fact that all the x-direction positions and sizes of the objects happen to be multiples of 0.05 m; therefore all the xfrac values are such multiples also. But this accidental congruence did not occur with respect to the z-direction settings.

If you desire to know what the grid would have been like had it not been adjusted to fit the objects, paste the line:

> obj,    grid, no

beneath each object-related block of settings in file gp25_3; then run the VR Editor. The grid will then appear as follows:

wherein the red lines have disappeared and the non-fitting of the IN-BLOCK and L-BLOCK objects is obvious.

Warning

If you intend to follow instructions of this tutorial, you will be constantly switching between the PHOENICS Commander and the VR Editor.

In this case you should beware of either clicking on File and Quit to exit the VR Editor without saving all the changes made, or clearing the boxes in the 'Save Current Settings' window each time when you close the VR Editor clicking on the top-right cross.

Only then will:

you preserve your Q1 unchanged,
the incl(gp25_3) statement will be retained and
you be able to make changes, when necessary, by editing the data file gp25_3.

The grid 'region' concept

Careful comparison of the above two grids reveals that the second lacks what the first exhibits, namely several red grid lines; and it may be noted that some at least of these red lines do coincide, in part, with boundaries of objects.

Indeed, if it is recognised that the INLET, OUTLET and IN-PLATE objects are not visible in this view, it may be concluded that all of the red lines do coincide with object boundaries.

In order to confirm this, remove or deactivate the 'grid, no' lines from gp25_3, re-run the VR-Editor, then click first on the 'Mesh toggle' button and then on the 'Geometry cells toggle' button . The image which you will see is as follows:

This shows all the cells which are within objects of volumetric type or adjacent to objects of cell-face type; and the red lines do run along the boundaries of each such group of cells.

Thus is revealed how the PHOENICS Satellite fits its grid to the objects: it divides each of the x, y and z dimensions into a sufficient number of so-called 'regions' to enable their boundaries to fit the objects, and then distributes the remaining grid lines as evenly as possible between them.

In the present case there are four regions in the x-direction and six in the z-direction, facts which are confirmed by the following lines in the gp25_3 file:

  > GRID,   RSET_X_1,      5, 1.000000E+00
  > GRID,   RSET_X_2,      5, 1.000000E+00
  > GRID,   RSET_X_3,      4, 1.000000E+00
  > GRID,   RSET_X_4,      6, 1.000000E+00
  > GRID,   RSET_Y_1,      1, 1.000000E+00
  > GRID,   RSET_Z_1,      1, 1.000000E+00
  > GRID,   RSET_Z_2,      4, 1.000000E+00
  > GRID,   RSET_Z_3,      3, 1.000000E+00
  > GRID,   RSET_Z_4,      3, 1.000000E+00
  > GRID,   RSET_Z_5,      3, 1.000000E+00
  > GRID,   RSET_Z_6,      1, 1.000000E+00

Here:

the first four lines correspond to the four x-direction regions, which have 5, 5, 4 and 6 intervals (cells) respectively in each, making a total of 20;
the fifth line corresponds to the one region of the y-direction, with only one interval in it; and
the remaining six lines correspond to the six z-direction regions, having respectively 1, 4, 3, 3, 3 and 1 intervals in them, making a total of 15.

The above lines, being shifted two spaces to the right of the first column, are not active in the Q1 file; but they have been inserted by the Satellite, when working in VR-Editor mode, as a record of what it has done by way of distributing the grid lines.

It is its candour in doing so which allows us to take command.

Changing the grid distribution

The following steps should now be taken so as to convince ourselves that we are indeed in command:

Activate the above-mentioned lines by deleting the first two spaces. Then re-run Satellite and compare the new q1ear with the previous one. They should be identical.
The same will be found in respect of the RESULT file, if the EARTH solver is also run.
Let us now edit the relevant lines in gp25_3 so that they are:
```
> GRID,   RSET_X_1,      10, 1.000000E+00
> GRID,   RSET_X_2,      10, 1.000000E+00
> GRID,   RSET_X_3,      8, 1.000000E+00
> GRID,   RSET_X_4,      12, 1.000000E+00
> GRID,   RSET_Y_1,      1, 1.000000E+00
> GRID,   RSET_Z_1,      2, 1.000000E+00
> GRID,   RSET_Z_2,      8, 1.000000E+00
> GRID,   RSET_Z_3,      6, 1.000000E+00
> GRID,   RSET_Z_4,      6, 1.000000E+00
> GRID,   RSET_Z_5,      6, 1.000000E+00
> GRID,   RSET_Z_6,      2, 1.000000E+00
```
which implies that we have doubled the number of cells in each x-direction and z-direction region.
When we re-run the VR-Editor, the grid will appear thus:

The pressure contours and velocity vectors, when the Earth solver has been run, will appear in the Viewer thus:
Let us see how this mere increase in the cell number has changed the non-dimensional pressure and temperature. Whether these complexes have been changed to approach our expectations expressed in the previous part of the tutorial.
The non-dimensional pressure now looks like in this picture.

Judging by the scale we could think that the results are now even less understandable than for the previous run as the minimum value is now approaching -7. However, the region where the non-dimensional pressure is negative, i.e. in the above IN-BLOCK space has become smaller. As to the outlet area, it is practically zero.
As to the non-dimensional temperature, it is as in the following picture.

It is practically zero in the major space of the labyrinth increasing with approach to the heat source (L-BLOCK) and behind the IN-BLOCK as in the previous run. In general, the non-dimensional temperature is becoming more uniform as the heat input from the L-BLOCK is rather small and the air inlet velocity is small too.
Evidently, the grid has been refined, as is reflected in the result file thus:

even though the command describing the grid in q1, namely:
```
 Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
 RSET(M,20,1,15)
```
remained unchanged.
Within each region, the grid distribution is evidently uniform. Thus is a consequence of the fact that the last argument of each '> GRID' line is 1.0.
Let those of the four x-direction regions be altered thus:
```
> GRID,   RSET_X_1,      10, 0.5 ! power law starting from low x
> GRID,   RSET_X_2,      10, -0.5 ! power law starting from high x
> GRID,   RSET_X_3,      -8, 2.0 ! symmetrical power law
> GRID,   RSET_X_4,      -12, 2.0,g ! symmetrical geometric progression
```
Displayed in the VR-Editor, with the Mesh and Wireframe toggles pressed, the grid now appears as:

about which the following remarks can be made:
- None of the non-dimensional variables are displayed within the solids. This is understandable for pressure and velocity; but about temperature further remarks will be made below.
- A full description of the significance of the 'GRID, RSET' command can be found in the document TR 326 here.
- In the first x-region (starting from the left), the 0.5 exponent clusters the grid cells towards the right, thus enabling the boundary layer on the plate to be better simulated.
- In the second region, the minus sign in front of 0.5 causes the clustering to be at the lower-x region rather than the higher;
- In the third region, clustering at both left and right is achieved by changing the exponent to 2.0 and placing a minus in front of the number-of-intervals argument, 8.
- In the fourth region, the presence of the optional g argument indicates that the geometric-progression rather than power-law furmula is to be used.
- Further information can be found in the Encyclopaedia entry, GRDPWR.
The pressure distribution and contours which result when the solver is run with this grid is shown here:

As to the non-dimensional pressure, it will be like this one:

Although the minimum non-dimensional pressure has been reduced a little bit more as compared with the previous run, the area of negative non-dimensional pressure has been also reduced and in the space between the IN-BLOCK and the OUTLET section this parameter is stably positive. We can attribute the increase of the absolute value of the minimum pressure to the round-off error.
The non-dimensional temperature looks like in the picture which follows.

It is uniform enough except for the areas behind the WALL-W and the IN-BLOCK where its maximum values are found, but these can be also explained by flow stagnation in these zones.

Parameterizing the grid distribution

Knowing now that we are indeed in control of the grid distribution, we can see how to proceed with parameterization:
The part of the file gp25_3 which sets the grid should be changed to the following, in which three parameters are declared and set for each grid region:

the number of grid intervals, n,
the power-law exponent or geometric-progression factor p,and
the character variable which chooses between power-law or progression, g (you can find more about character variables clicking here here).

Declaring grid parameters
integer(nx1,nx2,nx3,nx4,nz1,nz2,nz3,nz4,nz5,nz6)
real   (px1,px2,px3,px4,pz1,pz2,pz3,pz4,pz5,pz6)
char (gx1,gx2,gx3,gx4,gz1,gz2,gz3,gz4,gz5,gz6)
  
  Setting grid parameters
nx1=10;  nx2=10;   nx3=-8;  nx4=-12
nz1=2;   nz2=8;    nz3=6;   nz4=6;   nz5=6;   nz6=2
px1=0.5; px2=-0.5; px3=2.0; px4=2.0
pz1=1.0; pz2= 1.0; pz3=1.0; pz4=1.0; pz5=1.0; pz6=1.0
gx1=f;   gx2=f;    gx3=f;   gx4=t
gz1=f;   gz2=f;    gz3=f;   gz4=f;   gz5=f;   gz6=f

  Reading grid parameters
> GRID,   RSET_X_1,      nx1, px1, gx1
> GRID,   RSET_X_2,      nx2, px2, gx2
> GRID,   RSET_X_3,      nx3, px3, gx3
> GRID,   RSET_X_4,      nx4, px4, gx4
                                      
> GRID,   RSET_Y_1,      1,   1.0, f  
                                      
> GRID,   RSET_Z_1,      nz1, pz1, gz1
> GRID,   RSET_Z_2,      nz2, pz2, gz2
> GRID,   RSET_Z_3,      nz3, pz3, gz3
> GRID,   RSET_Z_4,      nz4, pz4, gz4
> GRID,   RSET_Z_5,      nz5, pz5, gz5
> GRID,   RSET_Z_6,      nz6, pz6, gz6

If this is done, and a further run made, the results will be precisely as before (if indeed your last run did have the grid settings implied above).

The influences of these parameters can now be explored one by one. For example, it may be interesting to see what is the effect of switching between the power-law and geometric progression options, or to experiment with the minus sign placed before that number of intervals or the p-attribute.

More important however is to discover how many intervals are needed for the numerical solution to be accurate.

That is the subject of section 4 of this tutorial, in preparation for which you should introduce the following few extra lines just above the 'reading grid parameters' section:

    Further declarations
integer(nx0,nz0)

  Further settings
nx0= 1; nz0=nx0
nx1=nx1*nx0; nx2=nx2*nx0; nx3=nx3*nx0; nx4=nx4*nx0
nz1=nz1*nz0; nz2=nz2*nz0; nz3=nz3*nz0; nz4=nz4*nz0; nz5=nz5*nz0; nz6=nz6*nz0

By introducing the two additional parameters, nx0 and nz0, making the first equal the second and linking the other nx's and nz's to them, we have thus made it possible quickly to change the fineness of the whole grid, without losing the ability subsequently to make further refinements in individual regions of the grid.

3. Parameters read from file

Before proceeeding further, it is appropriate to consolidate what has been done so far, which has involved introducing parameters, setting their values, introducing relationships between them, and observing their effects on the flow calculations.

The nature and rationale of these actions will become clearer if some of the entries made in the file gp25_3 are tranfserred to another file which we shall call params. This will be related to q1 and gp25_3 as indicated in the following sketch:
You should therefore now paste into your Q1 the last three lines of the following, so that it contains:

  Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
 RSET(M,20,1,15)
  save3begin
incl(params)
  save3end

The params file is being included at this point so that what it contains about the grid wil not be over-written, as it otherwise would be, by the RSET command which follows.

Now create the params file itself by pasting in the following lines, most of which have appeared before in the gp25_3 file:

   params file for library 
    case 621, used in parameterised q1 
  tutorial number 2

TEXT(parameterized 621 

  PIL-parameter settings ******************************************
xulast=2.
yvlast=1.
zwlast=1. 

  non-PIL parameter declarations **********************************
  object parameters
real(invel,inmass)

         
  object-existence parameter
boolean(inpxst,inbxst,wwxst,ewxst)
       

  parameters used in non-dimensional variables
real(intem,pext,inarea,heatflux)

  grid-related parameters
integer(nx1,nx2,nx3,nx4,ny1,nz1,nz2,nz3,nz4,nz5,nz6)
real   (px1,px2,px3,px4,py1,pz1,pz2,pz3,pz4,pz5,pz6)
char(gx1,gx2,gx3,gx4,gy1,gz1,gz2,gz3,gz4,gz5,gz6)

  Further declarations
integer(nx0,nz0)

  non-PIL parameter settings **************************************
   
  object-existence parameters
inpxst=t
inbxst=t
wwxst=t
ewxst=t
    

  non-dimensional variables
intem=20.0
pext=0.0
heatflux=100.
invel=0.05
inarea=yvlast*0.45*zwlast
inmass=inarea*rho1*invel

  Grid-related parameters
nx1=10;  nx2=10;   nx3=-8;  nx4=-12
ny1=1
nz1=2;   nz2=8;    nz3=6;   nz4=6;   nz5=6;   nz6=2
px1=0.5; px2=-0.5; px3=2.0; px4=2.0; py1=1.0
pz1=1.0; pz2=1.0; pz3=1.0; pz4=1.0; pz5=1.0; pz6=1.0
gx1=;   gx2=;    gx3=;   gx4=g;   gy1=
gz1=;   gz2=;    gz3=;   gz4=;   gz5=;   gz6=

  Further settings
nx0= 1; nz0=nx0
nx1=nx1*nx0; nx2=nx2*nx0; nx3=nx3*nx0; nx4=nx4*nx0
nz1=nz1*nz0; nz2=nz2*nz0; nz3=nz3*nz0; nz4=nz4*nz0; nz5=nz5*nz0; 
nz6=nz6*nz0

Those lines which have not been seen before are printed in red. They provide for the possibility of removing from the scenario four of the objects: IN-PLATE, IN-BLOCK, WALL-W and WALL-E.

Finally we suggest that you edit the existing data file gp25_3 and save it as, for example, gp25_4, by cutting and pasting your gp25_3 in the following way:

 GVIEW(P,0.000000E+00,-1.000000E+00,0.000000E+00)
 GVIEW(UP,0.000000E+00,0.000000E+00,1.000000E+00)

store(vabs)
(stored var nd_v is vabs/:invel:) 
(stored var nd_p is (p1-:pext:)/(.5*:rho1:*:invel:^2)) 
(stored var nd_t is (tem1-:intem:)*:cp1:*:inmass:/:heatflux:) 
mesg(variables stored)
stored
  > DOM,    SIZE,        2.000000E+00, 1.000000E+00, 1.000000E+00
> DOM,    SIZE,        xulast, yvlast, zwlast
  > DOM,    MONIT,       7.500000E-01, 5.000000E-01, 4.666670E-01
> DOM,    MONIT,       .75*zwlast, .5*zwlast, .466667*zwlast
> DOM,    SCALE,       1.000000E+00, 1.000000E+00, 1.000000E+00
> DOM,    SNAPSIZE,    1.000000E-02

  Reading grid parameters
> GRID,   RSET_X_1,      nx1, px1, gx1
> GRID,   RSET_X_2,      nx2, px2, gx2
> GRID,   RSET_X_3,      nx3, px3, gx3
> GRID,   RSET_X_4,      nx4, px4, :gx4:
                                      
> GRID,   RSET_Y_1,      ny1, py1, gy1  
                                      
> GRID,   RSET_Z_1,      nz1, pz1, gz1
> GRID,   RSET_Z_2,      nz2, pz2, gz2
> GRID,   RSET_Z_3,      nz3, pz3, gz3
> GRID,   RSET_Z_4,      nz4, pz4, gz4
> GRID,   RSET_Z_5,      nz5, pz5, gz5
> GRID,   RSET_Z_6,      nz6, pz6, gz6

> OBJ,    NAME,        H-BLOCK
  > OBJ,    POSITION,    0.000000E+00, 0.000000E+00, 9.500000E-01
> OBJ,    POSITION,    .0, .0, .95*zwlast
  > OBJ,    SIZE,        2.000000E+00, 1.000000E+00, 5.000000E-02
> OBJ,    SIZE,        xulast, yvlast, .05*zwlast
> OBJ,    GEOMETRY,    cube14
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        BLOCKAGE
> OBJ,    MATERIAL,    100,ALUMINIUM at 27 deg c
  > OBJ,    grid, no
 
> OBJ,    NAME,        L-BLOCK
  > OBJ,    POSITION,    0.000000E+00, 0.000000E+00, 0.000000E+00
> OBJ,    POSITION,    0., 0., 0.
  > OBJ,    SIZE,        2.000000E+00, 1.000000E+00, 5.000000E-02
> OBJ,    SIZE,        xulast, yvlast, .05*zwlast

> OBJ,    GEOMETRY,    cube4
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        BLOCKAGE
> OBJ,    MATERIAL,    100,ALUMINIUM at 27 deg c

> OBJ,    HEAT_FLUX,    0.000000E+00, heatflux
  > OBJ,    grid, no
 
  
if(ewxst) then 
> OBJ,    NAME,        WALL-E
  > OBJ,    POSITION,    2.000000E+00, 0.000000E+00, 5.000000E-01
> OBJ,    POSITION,    xulast, .0, .5*zwlast
  > OBJ,    SIZE,        0.000000E+00, 1.000000E+00, 5.000000E-01
> OBJ,    SIZE,        .0, yvlast, .5*zwlast
> OBJ,    GEOMETRY,    cube11
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        PLATE
  > OBJ,    grid, no
endif
 
if(wwxst) then
> OBJ,    NAME,        WALL-W
  > OBJ,    POSITION,    0.000000E+00, 0.000000E+00, 0.000000E+00
> OBJ,    POSITION,    .0, .0, .0
  > OBJ,    SIZE,        0.000000E+00, 1.000000E+00, 5.000000E-01
> OBJ,    SIZE,        .0, yvlast, .5*zwlast
> OBJ,    GEOMETRY,    cube11
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        PLATE
  > OBJ,    grid, no
endif 


if(inpxst) then 
> OBJ,    NAME,        IN-PLATE
  > OBJ,    POSITION,    5.000000E-01, 0.000000E+00, 3.000000E-01
> OBJ,    POSITION,    .25*xulast, 0., .3*zwlast
  > OBJ,    SIZE,        0.000000E+00, 1.000000E+00, 6.500000E-01
> OBJ,    SIZE,        .0, yvlast, .65*zwlast
> OBJ,    GEOMETRY,    cube11
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        PLATE
> OBJ,    POROSITY,     0.000000E+00
> OBJ,    SIDE,        BOTH
  > OBJ,    grid, no
endif 


if(inbxst) then 
> OBJ,    NAME,        IN-BLOCK
  > OBJ,    POSITION,    1.000000E+00, 0.000000E+00, 5.000000E-02
> OBJ,    POSITION,    .5*xulast, .0, .05*zwlast
  > OBJ,    SIZE,        4.000000E-01, 1.000000E+00, 6.500000E-01
> OBJ,    SIZE,        .2*xulast, yvlast, .65*zwlast
> OBJ,    GEOMETRY,    cube14
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        BLOCKAGE
> OBJ,    MATERIAL,    103,COPPER at 27 deg C
  > OBJ,    MATERIAL,    111, STEEL
  > OBJ,    MATERIAL,    106, GLASS
  > OBJ,    grid, no
endif 
  


> OBJ,    NAME,        INLET
  > OBJ,    POSITION,    0.000000E+00, 0.000000E+00, 5.000000E-01
> OBJ,    POSITION,    .0, .0, .5*zwlast
  > OBJ,    SIZE,        0.000000E+00, 1.000000E+00, 4.500000E-01
> OBJ,    SIZE,        .0, yvlast, .45*zwlast
> OBJ,    GEOMETRY,    cube3t
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        INLET
> OBJ,    PRESSURE,     0.000000E+00
> OBJ,    VELOCITY,     invel, 0.000000E+00, 0.000000E+00
> OBJ,    TEMPERATURE,  intem
> OBJ,    TURB-INTENS,  5.000000E+00
  > OBJ,    grid, no
 
> OBJ,    NAME,        OUTLET
  > OBJ,    POSITION,    2.000000E+00, 0.000000E+00, 5.000000E-02
> OBJ,    POSITION,    xulast, .0, .05*zwlast
  > OBJ,    SIZE,        0.000000E+00, 1.000000E+00, 4.500000E-01
> OBJ,    SIZE,        .0, yvlast, .45*zwlast
> OBJ,    GEOMETRY,    cube12t
> OBJ,    ROTATION24,        1
> OBJ,    TYPE,        OUTLET
> OBJ,    PRESSURE,     pext
> OBJ,    TEMPERATURE,  SAME
> OBJ,    COEFFICIENT,  1.000000E+03
> OBJ,    TURBULENCE,   SAME        , SAME
  > OBJ,    grid, no

In this, the only novelties introduced in the file and emphasized by red are if...endif pairs around the statements referring to the objects which may be removed from the scene.

To confirm that these work, now set to f (i.e. false) in the params file each of the logical parameters: inpxst, inbxst, wwxst, ewxst, wherein the xst stands for 'exists'; then introduce a proper file name into the include statement in the q1 incl(gp25_4) and make a run. You should then see something like this in respect of the pressure contours and velocity vectors:

The four excluded objects indeed no longer exist.

Having separated the declaration and setting of parameters from the lines which they influence, we can now confine attention to the file params, using this as the means whereby we explore the influence of grid fineness on accuracy.

4. How grid refinement affects accuracy

What is meant by accuracy

Of great practical interest is the general question: do the predictions made by a CFD code such as PHOENICS truly represent what would be found in practice if experimental measurements were made in the postulated scenario.

If the answer is affirmative, it would be reasonable to call the calculations 'accurate'; however, accuracy of that kind is not what concerns us here, for two reasons, namely:

no such experimental measurements have been made; and
our focus is on numerical not physical accuracy.

A numerical calculation is considered accurate if the quantities of practical interest which it computes are almost independent of the fineness and distribution of the computational grid.

In order to test it, one must:

decide which quantities we shall regard as being of most practical interest; and
make computations with grids differing finenesses in order to see how their values vary.

Testing grid independence

In what follows, attention will be focussed on the maximum and minimum values of the three non-dimensional variables introduced above, namely nd_p, nd_v, and nd_t.

For simplicity we shall make the grids uniform within each region; therefore the grid-related settings in the params file appear as:

pz1=1.0; pz2= 1.0; pz3=1.0; pz4=1.0; pz5=1.0; pz6=1.0
px1=1.; px2=1.; px3=1.0; px4=1.0
pz1=1.0; pz2= 1.0; pz3=1.0; pz4=1.0; pz5=1.0; pz6=1.0
gx1=f;   gx2=f;    gx3=f;   gx4=t
gz1=f;   gz2=f;    gz3=f;   gz4=f;   gz5=f;   gz6=f
  Further settings
nx0= 1; nz0=nx0

Restore the object-existence parameter to 't'

inpxst=t
inbxst=t
wwxst=t
ewxst=t

and also type at the bottom of the q1 file:

#maxmin ! to print maximum and minimum values on the screen
lsweep=1000 ! to ensure convergence
debug=t ! to ensure that highest and lowest values are
        ! printed in the RESULT file
STOP

You are invited now to perform runs with various values of nx0 of which the first will be with nx0= 1. If you do, you will probably obtain results such as are contained in the following table.

nx0	nd_pmax	nd_pmin	nd_vmax	nd_vmin	nd_tmax	nd_tmin
1	21.02	-4.39	4.18	-1.0E-10	8.104	5.74E-05
2	25.65	-2.828	4.686	-1.0E-10	6.354	-2.05E-06
3	26.71	-3.238	4.784	-1.0E-10	5.998	-2.56E-06
4	23.42	-7.467	4.743	-1.0E-10	5.871	-3.08E-06
5	21.77	-12.332	4.979	-1.0E-10	6.073	-4.61E-06
6	20.09	-12.650	5.093	-1.0E-10	6.143	-2.31E-05
8	22.85	-9.177	4.974	-1.0E-10	6.671	-1.78E-01
10	22.24	-6.371	4.221	-1.0E-10	7.327	-2.58E-01

This shows that the grid should have 100 times as many cells (nx0=nz0=10) as has the coarsest grid (nx0=nz0=1) if the results of the computations are to be(almost) independent of the fineness of the grid.

If you performed any of these runs yourself, you will have observed that the computer time increased considerably (probably by a factor of 100 squared over the whole range!) as the grid becomes finer.

This is why users of CFD codes often content themselves with grids which are much coarser than those which would procure grid independence.

5. Concluding remarks

Grid-fineness effects

In this tutorial the focus of attention was on parameterization, especially in relation to grid fineness and distribution.

The influence of the parameters was explored only to a limited extent, in that:

refinement was made equally in respect of the x- and z direction because of the setting: nz0=nx0;
neither the exponent in the power-law expression nor the multiplication factor in the geometric-expansion one were varied; and
the results were explored only in respect of maximum and minimum values of variables.

Systematic study of the unexplored avenues has however been rendered easy by the declaration and setting of the parameters in the params file, and the provision of formulae influenced by them in the data file gp25_4.

You are therefore invited to explore their influence for yourself; and you are positively advised to parameterize any other Q1 files which you create, of which the numerical accuracy of the resulting computations is important to you, so as to be able to establish what degree of grid-independence you are achieving.

What fineness of grid you should then use will depend on the purposes for which you are making the simulation. If you seek only a qualitative idea of the nature of the flow, there will be no reason to incur the expense of using a fine grid.

Alternatively, if you are designing an important process or item of equipment, you will probably be wise to use the finest grid which the size of your computer permits, regarding the maximised accuracy as justifying the expense.

In this case you may find it useful to make use of the so-called SPINTO feature, described here. This allows the finest-grid solution to be achieved more quickly by re-starting from a succession of coarser-grid ones.

Turbulence-model effects

Although no attention has been drawn to the fact, the flow which has been simulated is ostensibly a turbulent one, as may be deduced from

the line:
```
TURMOD(KECHEN)
```
in Group 7 of the original Q1 file, which has not been changed in this respect; and
the following lines consequential lines printed in satlog.txt:
```
variables stored)
>  
>  P1  (   1)   U1  (   3)   W1  (   7)   KE  (  12)   EP  (  13)
>  VABS( 147)   ENUT( 148)   EPKE( 149)   TEM1( 150)
```
which indicate that the variables:
- KE, the kinetic energy of turbulence,
- EP, the dissipation rate of turbulence energy,
- ENUT, the effective kinematic viscosity, and
- EPKE, which equals EP divided by KE.

The above lines indicate that a so-called 'turbulence model' has been used, and specifically the one with the name KECHEN.

These terms will not be explained here; however, Tutorial 3 in the parameterized-q1 series will devoted to explaining them, and to providing a systematic introduction to the important branch of Computational Fluid Dynamics which is known as Turbulence Modelling.

Tutorial pq1t2_1; Part 2Library case 621: Simulation of a Labyrinth Flow.