Tutorial pq1T1; Part 2
Extending Library case 163

Summary

This second part of the tutorial shows how the editing of an existing Q1 file can extend its simulation capabilities, in this case in respect of temperature-dependence of material properties.

Instructions will be provided concerning:

The relevant features of the PHOENICS Input Language will be explained as they are introduced, one by one.

Contents

  1. Introduction
  2. Preliminary ground-clearing operations
  3. Introducing temperature-dependent material properties
  4. Extending the menu-type input capabilities
  5. Concluding remarks

1. Introduction

1.1 How this tutorial is structured

The Q1 file represented by library case 163.htm was discussed in detail in Part 1; and it forms the starting point of Part2 which proceeds by making a succession of additions to it.

These additions will always (with one exception) be placed at the bottom of the file immediately above the (to be supplied) save25end marker, thus exploiting advantageously the fact that the PHOENICS Satellite module executes sequentially the instructions which it finds in the file, and so allows earlier-made (i.e. higher-up) instructions to be over-written by later-made (i.e. lower-down) ones.

To be more precise, the lines to be added will all appear between the markers

  save25begin
  save25end

so as to ensure that the protected mode of Satellite operation will be in use.
They will be introduced in Step 1 of this tutorial.

Why 25? Because all added-at-the-end lines in a Q1 file are regarded as belonging to Group 25.

The advantage is that the Q1 file which has been created during the course of the tutorial will contain a complete record of what has been done.

The tutorial consists of a sequence of small steps, of which you are advised to omit none; for otherwise the results which PHOENICS produces will differ from those which the text of the tutorial describes.

You are of course free to use any method of editing the Q1 file and running the Satellite and EARTH that you prefer; but the tutorial is written as though you are using the facilities of the PHOENICS Commander, as follows:

  1. Activate the Commander.
  2. Click on the Input-File Library icon.
  3. Click the 'set work dir.' button in the top row, then click 'select...' and browse for the working directory you have chosen or created for this run.
  4. Click on the 'choose by case number' button.
  5. Enter 163 in the white box,
  6. Click on 'Load it'. This will bring into view the following screen:

  7. Click on 'See q1'. The q1 file will be loaded in the editing window.
  8. Click on the 'Edit Page' button above; then go to 'File' and choose 'Save as' command.
  9. Type the name of your working directory in the 'Save in' window and 'q1' in the 'File name' box. Then click on the 'Save' button.
    The q1 of the case 163 is now copied to your working directory and it is this file which we shall be working on during the tutorial.
  10. Click on 'View or edit files'
  11. Click on 'q1'.
You should now be seeing a screen looking similar to this:


You are now in a position to edit the file.


2. Preliminary ground-clearing operations

Step 1. Suppressing the questions about 'caseno'

Because the Q1 we are working with, is at present the one that contains the questions about the case number already explained in Part 1, and it would be irritating to be forced to respond at each execution of the Satellite, we shall avoid the interactive procedure of choosing the proper case number here. The quickest way to suppress the question-asking is to use the Commander editor's 'find' and 'replace' features so as to replace:

readvdu

by:

  readvdu

i.e. to de-activate the command by moving it two spaces to the right; and then to perform the same action on the 'mesg' and 'mesgm' commands.

The way of doing so can be clarified by the following scheme.

Save the changes which you have made clicking on the File tab and then on the Save File command in the list. Only after that click on the S/E icon on the left  . Then the only evidence that the Satellite is running is the brief appearance on the screen of:

EARDAT file written for RUN 1, Library Case=163

whereafter the solver begins (and soon terminates) its execution.

Please now take the actions just described and assure yourself that they have the effect described.

Although the questions are not being asked, it is still necessary to tell PHOENICS what value of caseno to use. We will choose caseno=4; and our choice must be conveyed by typing that line into the q1 file just above save1end.

However, it is not of any importance exactly which caseno you will choose (you may specify any from 0 to 4), as its only purpose is to start the simulation process in the absence of the interactive caseno input via the keyboard.

Step 2. Switch the material back to steel

Before doing anything else, insert the following two lines at the bottom of the Q1 file but above the STOP, making sure that you leave the first two spaces of the lines blank:
  save25begin
  save25end

The bottom part of the Q1 file will be like this one.

All your subsequent entries in the file should lie between these; for they will prevent your work from being obliterated by the Satellite. These marker lines ensure that you will be working always in the 'protected mode'.

The caseno of the Q1 which you are using is 4, for which the material of the slab is copper. In order to revert to the original material, steel, insert above the save25end marker, the following line:


FIINIT(PRPS)=STEEL ! To require properties of steel to be used.

It is what precedes the exclamation mark which the Satellite acts upon; the words to the right of it are comments to assist anyone who reads the file to understand what you have done.

However if earlier you have set the caseno variable different from 4, you will not have to introduce this statement. The slab material has been changed from steel to copper only for caseno=4.

Then save the changes which you have made clicking on the File tab and then on Save File, as the picture shows.

After that run the case clicking on the S/E icon on the left, to confirm that KOND has the value 43 and that #RAT is still unity.

Step 3. Change the type of slab-face patches

It will be remembered that, although the steel slab in case 163 was defined as having its faces held at 0 degrees Celsius, what was actually done was to fix the temperatures at the centres of the nearest-to-face cells to 0 degrees; which is not quite the same; and this was why the numbers .005 and .995 appeared in the formula for the theoretical temperature distribution.

However, PHOENICS does possess a means of setting the true face temperatures, namely by setting the 'types' of the patches to WWALL and EWALL and by setting the third arguments of the COVAL statements to be 1.0.

Please therefore now insert the following lines at the bottom of the Q1 file, where they will cancel the effect of the higher-up setttings.


PATCH(minXface,WWALL,1,1,1,1,1,1,1,1)   ! To locate the low-x face
PATCH(maxXface,EWALL,NX,NX,1,1,1,1,1,1) ! To locate the high-x face
COVAL(minXface,TEM1,1.0,0.0)   ! To fix the values of the true faces
COVAL(maxXface,TEM1,1.0,0.0)   ! not the cell centres, to zero

[Of course, you need not type it unless you want to. Cutting and pasting is all that is needed.]

Don't forget to save the changes you make in already familiar way.Then run the case in already familiar way, so as to observe what now appears in the RESULT file. Open the RESULT file in the Commander's editing window by clicking on its tab on the bottom, enter 'w f' into the 'find' box and click on 'find'. Instantly the line containing the words 'flow field' will appear; and it will be seen there that #RAT is no longer exactly 1.0.

We must understand the reason for this.

Step 4. Change the print-out instructions

Let us therefore take a closer look at the temperatures in the cells which are close to one of the interfaces, by setting NXPRIN=1 which is a command which says: "print out values for every cell" (which would give an undesirably large amount of print-out) and IXPRL=5 which adds "but not for values of IX above 5" (and so limits the amount).

The lines you should add to the bottom of the Q1 are as follows:


nxprin=1   ! To print values for every cell
ixprl=5    ! up to the fifth, so that the near-face temperatures 
           ! are shown in the RESULT file.

Save the changes introduced and run the case going to the 'Run modules' window of the Commander, then clicking the 'S/E' icon of the Satellite operation in 'silent 'mode, so as to observe the consequential temperatures close to the interface.

It can be deduced that the values are satisfactory; for they show that the temperature is above zero even in the first cell, and that if each temperature is divided by the first their quotients are (very nearly) 1, 3, 5, 7 and 9, as they should be if the distribution is locally linear.

[Strictly speaking, the ratios are 1.0, 2.960E+00, 4.880E+00, 6.760E+00 and 8.600E+00, as a hand calculator will inform you.
But, if you do not have one, you can tell PHOENICS to calculate the values for you, by inserting the following lines in the Q1, wherein
tem1[1,1,1]
stands for the temperature at the first point of the grid:


formula=tem1/tem1[1,1,1]
(stored var #quo is :formula:)
Then the RESULT file will be found to contain the lines:
Field Values of #QUO
 IY=   1   1.0         2.960E+00   4.880E+00   6.760E+00   8.600E+00
 IX=          1           2           3           4           5
PHOENICS is as willing to perform simple arithmetic as it is to solve differential equations.

If you do introduce the above two lines, do not forget to remove them when they have served their purpose; otherwise, although the definition which you have just given to formula will be over-written in the next step of the tutorial, the variable #quo will continue to be calculated and printed.]

We conclude that it is the theoretical temperatures which are wrong; and a moment's reflection explains why: we have forgotten to up-date the formula for #RAT so as to accord with the new patch-type settings. For this case #RAT values were as follows.

Step 5. Change the formula to fit the changed PATCH and COVAL settings

This omission is quickly remedied by adding the following lines to the bottom of the Q1 file:
     
 FORMULA=TEM1/(0.5*1.E3*(XG)*(:XULAST:-XG)/KOND) ! Correct formula
              ! where
              ! 0.5 is a constant in the theoretical expression for
              ! the parabolic profile,
              ! 1.E3 is the volumetric heat input,
              ! KOND is the thermal conductivity of the material
 (STORED var #RAT is :FORMULA:) ! Calculates the temperature ratio
Now save the changes and run the case so as to observe whether #RAT is precisely 1.0. You will see that it is very nearly so; but not quite.

Why not? Is it the result of the fact that our 100-cell grid is too coarse?

Step 6. Choose an even coarser grid

We can determine whether the accuracy worsens with increasing coarseness by inserting, just above save25end the line:
nx=5
When you run the case, you will see that #RAT has taken values very much much farther from unity.

Why?

The answer is that setting NX is not enough: it is necessary also to change the sizes of the cells.

The following lines can be found in the RESULT file:

X-coordinates of the (higher) cell faces
    1.000E-03   2.000E-03   3.000E-03   4.000E-03   5.000E-03
from which it can be seen that the largest value is 0.005, not the 0.1 which was set as XULAST, the thickness of the slab.

Step 7. Setting the grid correctly

To correct the grid, copy from higher in the Q1 file to the bottom of the file the line:
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid
for the newly specified number of cells and run again. As the calculation shows, the grid has been rebuilt to fit the new setting for the cell number and the largest value is exactly what it should be, i.e. 0.1.

The #RAT values will be found again to be not far from 1.0, 1.111 being the largest, but less close to it than for NX=100.

It does indeed appear that the coarser grid gave the worse result.

Step 8. Use the coarsest possible grid. NX=1

To understand this, set NX to its smallest possible value, namely unity, by adding to the bottom of the file:
NX=1   
GRDPWR(X,NX,XULAST,1.0) 

Running this will show that #RAT =2.0. Why? It is a consequence of the fact that PHOENICS has only one place into which to insert its heat input, namely the central point of the single cell; and this heat has to be conducted from that point to the zero-temperature faces.

In the reality which we are seeking to simulate however, the heat input is DISTRIBUTED over the thickness of the cell [see the sketch below].

Thus the PHOENICS treatment implies a 'tent-shaped' temperature profile within the cell whereas the exact solution recognises it as parabolic.

That the tent and the parabola must have the same slopes at each face is explained by the equality of heat fluxes and thermal conductivity in both cases, thus the central-point temperature for the tent being 2.0 times that for the parabola.

Step 9. Studying the influence of the number of grid cells.

Let us now try larger values of NX, namely 2, 4, 8, 16 and so on. We expect NX for IX=1 to approach 1.0 as NX increases. Does it do so?

Let us type at the bottom of the q1 file the following two statements

NX=2   
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid,
save these changes and run the case. Then check what will be the #RAT value in the first cell, as well as in the second, being the last. It is 1.333 in both cells, as in this case our approximated tent-shaped temperature diagram has turned into a trapezoid one. What conclusion can be made by comparison of these two cases? We were able, quite obviously, to reduce the simulation error from 100% (NX=1) to about 33% (NX=2).

Repeat the previous actions to make runs, changing NX in the geometrical progression from 4 to 256.

Here are the results which you should find:


   NX   1   2      4      8     16    32     64    128   256
   #RAT 2.0 1.333  1.143  1.067 1.032 1.016  1.008 .7843 .3914

Step 10. Analysis of the results obtained.

These values do tend towards 1.0 as NX increases AT FIRST. But why are the last two values BELOW 1.0? It is because we have forgotten to tell PHOENICS to enlarge the patches correspondingly.

This did not matter for the NX values which were less than 100; for the third and fourth arguments of the higher-up-in-the-Q1-file PATCH statement equalled 100, by reason of the lines:

and earlier in Q1 inside save3 markers NX has been set to 100;
so the PHOENICS solver module was clever enough to reduce this to the smaller NX value which was set later, because 100 was the largest value that the argument could take.

However, when a larger NX was set, PHOENICS had no way of knowing whether or not it was your intention to increase the grid beyond the end of the patch; therefore it did nothing, with the result which we have just seen.

The COVALS, it should be noted, do not need to be up-dated; for they relate to the named patch, wherever it might be.

Let us do so for NX=256, by copying the patch commands to below the NX setting, thus:

NX=256
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid,
PATCH(minXface,WWALL,1,1,1,1,1,1,1,1) ! to locate the low-x face
PATCH(maxXface,EWALL,NX,NX,1,1,1,1,1,1) ! to locate the high-x face
PATCH(HEATER,volume,1,nx,1,1,1,1,1,1)  ! to show that the volumetric

The result file then shows to IX=1 value of #RAT to be 1.002 , which is satisfactory.

Let us now increase NX further to 10000 by adding the following lines to the bottom of the q1 file.


NX=10000
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid,
PATCH(minXface,WWALL,1,1,1,1,1,1,1,1) ! to locate the low-x face
PATCH(maxXface,EWALL,NX,NX,1,1,1,1,1,1) ! to locate the high-x face
PATCH(HEATER,volume,1,nx,1,1,1,1,1,1)  ! to show that the volumetric
You will see an error message if you try to run EARTH. Go no further. The subject lies beyond the scope of this tutorial. Do not forget to delete the last changes introduced.

On the basis of the following analysis it is possible make the following conclusion.
A reasonable setting of of the number of cells is of great importance as it influences the simulation results: too coarse a grid introduces a great error and too fine cannot be processed by the solver; so the number of cells should always be necessary and sufficient to provide satisfactory accuracy without increasing the time of simulation too much.

Usually it is not possible to say in advance what grid would be sufficient unlike in this case, when the accurate solution is available. That is the reason why in most cases it is necessary to make several runs with different cell numbers in the calculation grid.

Step 11. Setting the thermal conductivity.

The thermal conductivity of most materials varies with temperature.

It is therefore useful to investigate the effect of such variations on the temperature distribution produced by a given heat flux. Before explaining how to do this, it will be useful to describe the various ways of introducing CONSTANT conductivities. First let us revert to the earlier grid, as follows:


NX=100                       
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid,
IXPRL=NX
NXPRIN=NX/5
and run the executable. The results should be as follows.

Step 12. How not to set the conductivity.

Now we abandon the use of the PROPS file, by nullifying the above STORE(PRPS) by way of the command:
SOLUTN(PRPS,N,N,N,N,N,N)
and, supposing it to be reasonable to do so, insert the value of KOND to be used at every point by way of the initial-value statement:
FIINIT(KOND)=43.0
Both these statements should, like all others in this tutorial, be placed just above the 'save25end' statement at the bottom of your q1 file. The last seven lines of this file should now look like this:

NX=100                       
GRDPWR(X,NX,XULAST,1.0) ! To create a uniform grid,
IXPRL=NX
NXPRIN=NX/5
SOLUTN(PRPS,N,N,N,N,N,N)
FIINIT(KOND)=43.0
  save25end
Running the case then shows that the temperatures are much higher than before and that KOND is being reported as 1.0E-05, i.e. much lower that the one which we attempted to impose.

#RAT is still unity (or very close to it), because the formula has used the low value of KOND in both cases, i.e. for simulated temperature and accurate one.

What has happened is that, although PHOENICS will have accepted 43.0 as its initial value of KOND, it did not receive the instruction to use that value, which was formerly implicit in the setting of the value of PRPS to 'steel'.

Having no other instructions, PHOENICS therefore reverted to its default practice which involves deducing thermal conductivity from the Prandtl Number.

The relevant default settings are:

PRNDTL(TEM1)=1.0  ! the Prandtl number
RHO1=1.0          ! the density
CP1= 1.0           ! the specific heat
ENUL=1.0E-05        ! the kinematic viscosity
Since the Prandtl Number is defined as ENUL*RHO1*CP1/KOND, it is not surprising that KOND has been reported as 1.0E-05.

A remark to critics: We are here studying how PHOENICS does behave, not how a newcomer might reasonably believe that it should behave.

Therefore the only conclusion we can make here will be this one:
It is not possible to set the conductivity in the most obvious way but...

When at Rome, do as the Romans do!

Although most obvious things are not always realized, we shall further show what exactly you can do in the land of PHOENICS.

Step 13. A successful way of setting the conductivity.

If a conductivity of 43.0 is wanted, in the light of what has just been said, the most obvious way is to set:
PRNDTL(TEM1)=1.0E-05/43.0
If the Satellite and solver are now run, KOND will indeed be reported as 43.0 and the temperatures will have the appropriate values, exactly as these of Step 11.

Pay attention that Field Values of PRPS are missing as we have earlier set the index of the material used to zero.

Step 14. A still more successful way of setting the conductivity

However, since neither viscosity, nor specific heat nor density have any relevance to the present problem, and indeed the viscosity of solid steel must be considered as infinite, the just-used procedure, which derives from the years when PHOENICS was concerned with fluids only, may be regarded as rather circuitous. There exists another way: the conductivity can be set directly thus:
PRNDTL(TEM1)=-43.0
whereby the minus sign is simply a signal to PHOENICS to treat the number which follows it as the thermal conductivity and not as a Prandtl number at all. Running the executable now will confirm that the results are the same as before.


3. Introducing temperature-dependent material properties

All the above methods shown in section 2 concern constant conductivities; now in section 3 a method will be described for setting those which vary. It involves the use of the In-Form (Input via Formulae) feature of PHOENICS.

Step 1. Preliminary actions to set temperature-dependent conductivity

First an example will be shown, namely that in which the thermal conductivity increases by AA% for each degree of temperature rise and that AA equals, say, 1.0.

We can express this fact in a way which PHOENICS understands by writing:


REAL(AA)               ! To declare the variable aa
AA=1.0                 ! To set its value
CHAR(CONDUCT)          ! To have a character variable with which to
CONDUCT=PRNDTL(TEM1)   ! show that PRNDTL(tem1) is the PIL variable
                       ! to be used
FORMULA=(43.0*(1.0+:AA:*TEM1/100)) ! To express the desired formula &
(PROPERTY :CONDUCT: IS :ENUL:/:FORMULA:) ! to use it via Prandtl No. 
Because the temperature rise was not very large with the previous heat input, let us increase it 1000-fold, remembering to bring down the COVAL statement which uses it:

HEATINPT=HEATINPT*1000.  ! but don't forget the COVAL which uses it 
COVAL(HEATER,TEM1,FIXFLU,HEATINPT) ! fix the heat input to HEATINPT
                                     and the formula for #RAT 
FORMULA=TEM1/(0.5*HEATINPT*(XG)*(:XULAST:-XG)/KOND) ! Correct formula
(STORED var #RAT is :FORMULA:) ! Calculates the temperature ratio
LSWEEP=5   ! because the conductiviy has to be re-calculated after
           ! the latest approximation to the temperature has been 
           ! calculated, we need several iterations. 5 may be enough
ISG21=LSWEEP ! This ensures that all 5 sweeps are completed
NPRINT=1     ! and this causes print-out to occur after every sweep,
IXMON=NX/2   ! while this will provide a line-printer plot of the
             ! centre-slab temperature versus sweep
  

Step 2. Running the case

Running the case precisely as specified proves to be disappointing: the conductivity stays at 43.0! However, the reason proves to be that we left PRNDTL(TEM1) equal to -43.0; and this meant that the In-Form settings were by-passed. Re-setting its default value is needed, thus:

PRNDTL(TEM1)=1.0
Then In-Form is allowed to do its work; but not at the first sweep where the result appear, as before, thus:
 
     Flow field at ITHYD=   1, IZ=   1, ISWEEP=     1, ISTEP=    1
     
     Field Values of #RAT
     IY=   1   1.005E+00   1.0         1.0         1.0         1.0
     IX=          1          21          41          61          81
     
     Field Values of KOND
     IY=   1   4.300E+01   4.300E+01   4.300E+01   4.300E+01   4.300E+01
     IX=          1          21          41          61          81
     
     Field Values of TEM1
     IY=   1   5.814E-01   1.895E+01   2.802E+01   2.779E+01   1.826E+01
     IX=          1          21          41          61          81
The reason is that it is only AFTER the first temperatures have been calculated that In-Form knows by how much to increase the conductivity. There are therefore sweep-to-sweep changes; but by the fifth sweep the print-out has reached it final content, namely:
     
     Field Values of #RAT
     IY=   1   1.005E+00   1.080E+00   1.111E+00   1.110E+00   1.078E+00
     IX=          1          21          41          61          81
     
     Field Values of KOND
     IY=   1   4.325E+01   5.050E+01   5.371E+01   5.363E+01   5.024E+01
     IX=          1          21          41          61          81
     
     Field Values of TEM1
     IY=   1   5.781E-01   1.743E+01   2.492E+01   2.473E+01   1.684E+01
     IX=          1          21          41          61          81    
Comparison of the two outputs shows that: Inspection of the values of #RAT, however, shows that they differ from unity. Does that mean that the calculations are incorrect?

No; what it results from is the fact that the formula which we have used for calculating #RAT no longer corresponds to the exact solution of the differential equation.

Step 3. Checking the accuracy of the PHOENICS solution

The exact solution of the differential equation for the problem which PHOENICS has just solved is somewhat complex; but if our desire is solely to check whether PHOENICS is able to solve varying-conductivity problems correctly, we can do so by making the heat-input non-uniform also.

Indeed, the task of checking is made easiest by using a concentrated rather than a distributed heat source, at the central plane of the slab.

For this, we need an odd number of cells, with NX=101, say, not 100; for only then will the grid truly have a central point. Were we to leave NX at 100, we should have to provide sources at IX=50 and IX=51, for symmetry; then the centre of the slab would lie between them. Therefore its temperature would not be calculated, and so could not be printed for checking.

This can be effected by introducing into the Q1 file the following statements:

NX=NX+1                                   ! to create an odd 
GRDPWR(X,NX,XULAST,1.0)                   ! number of cells, i.e. 101 
PATCH(HEATER,CELL,NX/2+1,NX/2+1,1,1,1,1,1,1) ! to supply heat only to the 
                                          ! central 51st cell, and at a rate
COVAL(HEATER,TEM1,FIXFLU,10.*2.*43./(xulast/2.)) ! which would make the
                                          ! temperature arbitrarily equal 10 degrees
                                          ! if the conductivity were 43.
NXPRIN=1         ! to print the results for every cell 
                 ! on either side of the central cell
IXPRF=NX/2-1     ! from 49th 
IXPRL=NX/2+3     ! up to 53d
When you do this, you will have the following results
Field Values of #RAT
 IY=   1   3.454E-01   3.524E-01   3.597E-01   3.524E-01   3.454E-01
 IX=         49          50          51          52          53
 
 Field Values of KOND
 IY=   1   4.695E+01   4.703E+01   4.710E+01   4.703E+01   4.695E+01
 IX=         49          50          51          52          53
 
 Field Values of TEM1
 IY=   1   9.183E+00   9.364E+00   9.545E+00   9.364E+00   9.183E+00
 IX=         49          50          51          52          53
and will see that the central temperature is not 10.0 degrees but 9.545, as is appropriate because the conductivity has risen, on average, by nearly 5%.

You may care to experiment with other values of AA, and other choices for centre-point temperature, and make other checks to assure yourself that PHOENICS is correctly calculating temperatures which accord with the exact solution of the differential equation; this is, for the chosen linear variation of conductivity with temperature, the quadratic equation of the form:

T + AA*T**2/2 = constant * x

Step 4. Understanding COVAL, FIXFLU and FIXVAL

Before proceeding further it is necessary to learn more about how sources of heat, or, for that matter, of mass momentum or other quantity, are represented in PHOENICS.

COVAL is a command with four arguments as seen above, thus:

COVAL(HEATER,TEM1,FIXFLU,10.*2.*43./(xulast/2.))
of which the third, FIXFLU, is used as a signal meaning:
insert a fixed (heat) flux with the following magnitude,
and the fourth is the magnitude itself.

However, other arguments can be used. For example, to fix the value of temperature at 10 degrees precisely, insert at the bottom of the Q1 file the line:

COVAL(HEATER,TEM1,FIXVAL,10.)
then make the run.

You will find that the central-point (IX=51) temperature, is indeed 1.000E+01.

It is natural to suppose that FIXVAL is another signal to the solver, signifying:
insert a fixed value (of temperature) with the following magnitude.

This is indeed so; but a more meaningful interpretation is:

insert a heat source of such a magnitude as will fix the temperature to the following value.

In order to find out more about FIXFLU and FIXVAL (and incidentally to discover some of the idiosyncrasies of the PHOENICS Satellite), paste the following lines at the bottom of your Q1 file:

NX       ! To show how SATELLITE reports values
FIXFLU   
         ! To cause the value of FIXFLU to be displayed. 
FIXVAL   
         ! To do the same for FIXVAL. (But neither will work!)
         ! So now to trick SATELLITE into revealing their values
REAL(HOWBIG)     ! To declare a real variable 
mesg(The next number is FIXFLU/10
HOWBIG=FIXFLU/10 ! To set it equal to FIXFLU/10
HOWBIG           ! NOW to find out FIXFLU's value
HOWBIG=FIXVAL/10 ! To do the same for FIXVAL
mesg(The next number is is FIXVAL/10
HOWBIG
and the run the case. Some messages will appear briefly on your screen and then disappear, probably too rapidly for you to read them.

Never mind. Such messages are always written also to the file satlog.txt, which you can open by clicking on the satlog.txt icon indicated below.

You will then find them to be as follows:

>  NX      =       101
>  FIXFLU  = FIXFLU
>  FIXVAL  = FIXVAL
> The next number is FIXFLU/10
>  HOWBIG  = 2.000000E-11
> The next number is FIXVAL/10
>  HOWBIG  = 2.000000E+09
>  EARDAT file written for RUN   1, Library Case=
The explanation is this: This is why we saw the unhelpful messages:
>  FIXFLU  = FIXFLU
>  FIXVAL  = FIXVAL
and were able to overcome this by introducing the factor of 1/10.

The essential facts which you need to know now about COVAL and its four arguments are that:

In order to consolidate your understanding, now make runs in which you replace the constant, FIXVAL, by a smaller number, e.g. 1.e5, thus:

COVAL(HEATER,TEM1,1.e5,10.)
You will find, when you run the solver, that the centre-point temperature is less than 10.0, namely 9.823 . Clearly the source which was applied was not quite big enough to raise the temperature to 10.0 .

Nett source values are printed in the RESULT file as:

 Nett source of TEM1 at patch named: MINXFACE =-8.862511E+03
 Nett source of TEM1 at patch named: MAXXFACE =-8.862512E+03
 Nett source of TEM1 at patch named: HEATER   = 1.772499E+04
wherein we note with satisfaction that the positive source at the heater is balanced by two negative sources at the slab faces equal to half its absolute magnitude.

Now try, in turn, the values 1.E6, 1.E7 and 1.E8 .
The corresponding printed central temperatures will be 9.982, 9.998, and 10.0;
and the nett heat fluxes should be: 1.8026E4, 1.8053E4 and 1.8060E4.
[However, because different computers have different round-off errors, you might be unable to get precisely the last flux value on your machine.]
Obviously the 2.E10 value ascribed to FIXVAL was bigger than necessary.

If you now re-run the case with FIXVAL as the third argument, to see what value of nett source is printed, you may be surprised to see that it is not printed at all (although the face sources are printed; and they sum satisfactorily to -1.806e+4).
The reason is that experience shows that, so successful is the technique in procuring closeness of VAR to A4, that when the quantity A3*(A4-VAR) is calculated at source-print-out time, round-off error causes inaccuracy.

Therefore, if an accurate value is required, you should use a much smaller constant than 2.E10.

The name COVAL is of course an acronym derived from COefficient-and-VALue, wherein COefficient refers to A3 and VALue to A4; and the general purpose of the command is to specify sources which depend linearly on the solved-for variable.

The very large coefficient used for fixing values and the very small one used for fixing fluxes are two extreme cases; but coefficients of intermediate value are often used.

For example, as will be seen in later tutorials, the loss of heat from a surface to a surrounding atmosphere can be represented by setting A3 as the heat-transfer coefficient and A4 as the atmospheric temperature. Then, if the surface temperature is greater than that of the atmosphere, the heat source at the surface is negative.

Step 5. Fitting the heat source to a prescribed temperature profile

In the last step, the #RAT values were of course far from unity, because we had not re-defined it appropriately to the concentrated heat source.

In the present step we shall pose the question: what would the heat-source distribution have to be in order for the temperature distribution to be that expressed by the formula selected in Step 2?

Why should we do so? Partly to satisfy the curiosity of those who are interested in such matters; but also in order to learn about how sources in PHOENICS can be specified in order to meet any requirement.

The following lines, which will be explained below, should now be pasted at the bottom of the Q1 file:


NX=NX-1                                ! to restore the original 
GRDPWR(X,NX,XULAST,1.0)                ! 100-cell grid,             
PATCH(HEATER,volume,1,nx,1,1,1,1,1,1)  ! patch location,
COVAL(HEATER,TEM1,FIXFLU,HEATINPT)     ! and heat input
FORMULA=(0.5*:HEATINPT:*(XG)*(:XULAST:-XG))/KOND ! The temperature
                  ! distribution used in the #RAT definition above  
(SOURCE OF TEM1 AT HEATER IS COVAL(1.E8,:FORMULA:)) ! source which 
                     ! will produce the desired temperature distribution 
(STORED VAR #SOR IS 1.E8*(:FORMULA:-TEM1)) ! to have it printed         
PLOT(TEM1PLOT,#SOR,0.0,0.0)               ! to have it plotted
The lines printed in black merely restore the settings to those which prevailed in steps 1 and 2; it is those printed in blue which are new.

The first of these replaces the previous specification of the heat input to the patch named HEATER, namely:

 
HEATINPT=HEATINPT*1000.  
COVAL(HEATER,TEM1,FIXFLU,HEATINPT) 
And it differs from it in several ways. Thus:
  1. It introduces the In-Form keyword 'SOURCE', which, like the earlier-encountered keywords 'STORED' and 'PROPERTY', immediately follows an opening bracket placed in the first or second column.

  2. The word COVAL appears not at the start of a following line but within the statement itself.

  3. The COVAL used earlier was a PIL command, dating from the earliest days of PHOENICS; but the COVAL in the more-recently-introduced SOURCE statement is an InForm function.

  4. The COVAL function has two arguments, in contrast to the four of the COVAL command, of which however the first (the patch name) and the second (the variable in question) do also appear in the SOURCE statement.

  5. The arguments of the COVAL function therefore correspond to the third and fourth arguments of the command; however they are no longer FIXFLU and HEATINPT, but a large number, 1.E8 and :formula:.

  6. The use of the large number as the coefficient suggests that the Step-4 method of fixing values by way of sources is being used. This is true.

  7. The presence of the formula as the second argument indicates in what way the function is superior to the command: it can accept expressions, not constants only.

  8. Those who recall the first-encountered appearance of the command in this tutorial, namely:
    FORMULA=(0.5*:HEATINPT:*(XG)*(:XULAST:-XG))/KOND)) 
    (SOURCE OF TEM1 AT HEATER IS COVAL(1.E8,:FORMULA:))  
    
    may doubt the just-made statement; for is not:
    0.5*:HEATINPT:*(XG)*(:XULAST:-XG))/KOND
    an expression?

    The answer is that it is an expression which the SATELLITE evaluates immediately, converts into a constant, and transmits via the EARDAT file to the PHOENICS solver. That same constant is there used as the heat input at every cell in the grid.

    The formula in the COVAL function, by contrast, is transmitted to the solver verbatim. There it is interpreted as an instruction to introduce a different source to each cell, according to its position (indicated by xG) and to a not-yet-known value of KOND which depends on the temperature which will prevail there.

  9. It is should be mentioned that, although the constant 1.E8 was used as the first argument of the COVAL function, an expression could have been used there too.
    Full information about In-Form and its capabilities can be found by clicking here or here.

When the case is run, it will be found that:

The last two facts are shown clearly by the line-printer plot below.


 ************************************************************
 PATCH(TEM1PLOT,PROFIL,   1, 100,   1,   1,   1,   1,     1,     1)
 PLOT(TEM1PLOT,#SOR, 0.000E+00, 0.000E+00)
 PLOT(TEM1PLOT,TEM1, 0.000E+00, 0.000E+00)
   Variable    1 = #SOR   2 = TEM1
     Minval=  8.053E+05  5.671E-01
     Maxval=  1.379E+06  2.353E+01
     Cellav=  9.810E+05  1.627E+01
 1.00 +11..+....+....+...2222222222222...+....+....+..11+
 0.90 +11            22222           22222            11+
 0.80 +  1        222                     222        1  +
 0.70 +   1     222                         222     1   +
 0.60 +    1  22                               22  1    +
 0.50 +     22                                   22     +
 0.40 +    2211                                 1122    +
 0.30 +  22    111                           111    22  +
 0.20 + 22       11111                   11111       22 +
 0.10 +22            111111111111111111111            22+
 0.00 22...+....+....+....+....+....+....+....+....+...22
      0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
 the abscissa is       X  .  min= 5.00E-04 max= 9.95E-02
 
 ************************************************************

The heat-source distribution is shown by the curve marked as 1; and the table above shows it to vary from a minimum of 8.053E+05 to a maximum of 1.379E+06, its average value being 9.810E+05.


4. Extending the menu-type input capabilities

The final section of this tutorial is concerned with creating menu-type data-input facilities for users who are unable or unwilling to enter data by way of PIL or In-Form statements edited in the Q1 file.

To illustrate the facilities, we shall first convert our heated-slab problem into that of heat loss from a fin, such as is discussed in elementary text-books on heat transfer and which is sketched below.

In preparation for this, you should now paste into the bottom of your Q1 file the following line:

label sect4
This is a label to which the Satellite will jump if it finds somewhere the instruction:
goto sect4
This statement, however, should be placed not at the bottom of Q1 but higher up, as is shown below:

Its purpose is to make the Satellite do exactly what we want in Section 4, omitting the subsequent actions.
You should now paste the following beneath the label sect4:

  deactivate unneeded items
maxxface=skip                 ! No longer needed
solutn(prps,n,n,n,n,n,n)
output(#rat,n,n,n,n,n,n)

  declare new real variables
real(hi,th,wi,alfa,lmda,Tair,Twal)
integer(set)
  set their default values
hi=0.1     ! height    10 cm
th=0.001   ! thickness 1 mm
wi=0.02    ! width     2 cm
alfa=10.0 ! heat-transfer coeff. for 5mm thick boundary layer
lmda=43.  ! thermal conductivity of steel
Tair=20.  ! degrees Celsius
Twal=80.  ! degrees Celsius
These statements first deactivate some unneeded items, then introduce and set values for new variables which appeared on the above sketch of the fin.

In order to use these we shall insert them into the PIL statements with which we are already familiar, as follows:

  Interactive statements to be placed below here
  
  Interactive statements to be placed above here


  Derived settings
xulast=wi  ! To fit the domain size to the fin
yvlast=hi
zwlast=th

prndtl(tem1)=-lmda ! to set the fin thermal conductivity 
fiinit(kond)=lmda		
COVAL(minXface,TEM1,1.0,Twal) ! To set fin root temperature

PATCH(HEATER,high,1,nx,1,1,1,1,1,1) ! To represent heat loss
COVAL(HEATER,TEM1,2.*alfa,Tair)     ! from fin (2 sides) 
lsweep=2 ! Sufficient for constant alfa and lmda
If you perform a run, you will obtain the following results:
 Nett source of TEM1 at patch named: MINXFACE = 2.261673E+00
 Nett source of TEM1 at patch named: HEATER   =-2.261456E+00
 pos. sum= 2.261673E+00 neg. sum=-2.261456E+00
 nett sum= 2.171993E-04




 ************************************************************
 PATCH(TEM1PLOT,PROFIL,   1, 100,   1,   1,   1,   1,     1,     1)
 PLOT(TEM1PLOT,TEM1, 0.000E+00, 0.000E+00)
   Variable    1 = TEM1
     Minval=  7.482E+01
     Maxval=  7.995E+01
     Cellav=  7.654E+01
 1.00 11...+....+....+....+....+....+....+....+....+....+
 0.90 + 111                                             +
 0.80 +   1111                                          +
 0.70 +      1111                                       +
 0.60 +         1111                                    +
 0.50 +            1111                                 +
 0.40 +                1111                             +
 0.30 +                    11111                        +
 0.20 +                        111111                   +
 0.10 +                              111111111          +
 0.00 +....+....+....+....+....+....+....+...111111111111
      0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
 the abscissa is       X  .  min= 1.00E-04 max= 1.99E-02
 
 ************************************************************
This shows that: We shall now make it easy to determine how these results will vary with changes of input parameters by creating a suitable menu. This is done by pasting into the Q1 between two comments
  Interactive statements to be placed below here
  
  Interactive statements to be placed above here  
the following:
  label start 
do ii=1,40 ! to clear the screen
mesg(
enddo
mesg(1. fin height=:hi: 
mesg(2. fin width=:wi:
mesg(3. fin thickness=:th:
mesg(4. fin conductivity=:lmda:
mesg(5. heat transfer coefficient=:alfa:
mesg(6. wall temperature=:Twal:
mesg(7. air temperature=:Tair:
  
mesg(choose data item to change 1,2,3,4,5,6,7 or RETURN
do ii=1,10 ! to clear the screen
mesg(
enddo  
readvdu(item,int,0)
These lines evidently send messages to the screen stating what are the current values of the seven input parameters: fin height, fin width, etc.
They then invite the user to state which data item he or she wishes to change.
The readvdu(item,int,0) statement awaits the user's choice of item and treats it as an integer with the default value of zero, in accordance with the Encyclopaedia entry: readvdu.

The result will be, that when the Satellite is run, you will be presented with the following screen:

Now paste the following lines:

if(item.eq.1) then     ! do this for item=1
mesg(enter desired height or press enter
readvdu(hi,real,hi)
goto start
endif
These allow the user to insert the desired value of the fin height, hi.
The goto start causes control to return to the question-asking page where the user should see his choice displayed.
For example, if he chooses 1.1234 meters, he should see.

The other six input items can be made accessible in the same way by pasting in the following lines, about which no further explanation is needed.

if(item.eq.2) then
mesg(enter desired width or press enter
readvdu(wi,real,wi)
goto start
endif

if(item.eq.3) then
mesg(enter desired height or press enter
readvdu(th,real,th)
goto start
endif

if(item.eq.4) then     ! do this for item=4
mesg(enter desired conductivity or press enter
readvdu(lmda,real,lmda)
lmda
goto start
endif

if(item.eq.5) then     ! 
mesg(enter desired height or press enter
readvdu(alfa,real,alfa)
goto start
endif

if(item.eq.6) then     ! 
mesg(enter desired wall temperature or press enter
readvdu(Twal,real,Twal)
goto start
endif

if(item.eq.7) then     ! 
mesg(enter desired surrounding air temperature or press enter
readvdu(Tair,real,Tair)
goto start
endif
Now if the user wishes to change the wal value to 100, he or she will have to type the number 6 in the 'Enter your answer here:' box at the bottom of the screen, press 'Enter', then type the required value of the wall temperature and press 'Enter' again.

Because the previously made changes are not retained, the user should repeat the change that he has introduced to the fin height under the number 1, which is the ordinal number of this very item.
Then the screen will show the following.

In the end you will have to press 'Enter' one more time to leave the menu. Finally, the solver run will show that the user's choices have been received by the solver from the Satellite and responded to correctly.

Thus:

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NY      =         1
 YVLAST  = 1.123400E+00
 ************************************************************
and
 PATCH(HEATER  ,HIGH  ,   1, 100,   1,   1,   1,   1,     1,     1)
 COVAL(HEATER  ,TEM1, 2.000000E+01, 2.000000E+01)
Finally the temperature distribution will be seen to be:
 ************************************************************
 PATCH(TEM1PLOT,PROFIL,   1, 100,   1,   1,   1,   1,     1,     1)
 PLOT(TEM1PLOT,TEM1, 0.000E+00, 0.000E+00)
   Variable    1 = TEM1
     Minval=  9.309E+01
     Maxval=  9.993E+01
     Cellav=  9.538E+01
 1.00 11...+....+....+....+....+....+....+....+....+....+
 0.90 + 111                                             +
 0.80 +   1111                                          +
 0.70 +      1111                                       +
 0.60 +         1111                                    +
 0.50 +            1111                                 +
 0.40 +                1111                             +
 0.30 +                    11111                        +
 0.20 +                        111111                   +
 0.10 +                              111111111          +
 0.00 +....+....+....+....+....+....+....+...111111111111
      0   .1   .2   .3   .4   .5   .6   .7   .8   .9  1.0
 the abscissa is       X  .  min= 1.00E-04 max= 1.99E-02
 
 ************************************************************
These results confirm that the desired menu-type input-data facility has been created and successfully used.


5. Concluding remarks

Having reached the end of the tutorial, you may wish to review some of the earlier steps.

This can easily be achieved by:

  1. inserting in the Q1 file the command STOP immediately below the set of commands which correspond to the step which you wish to repeat; and
  2. if these commands lie below the 'GOTO SECT4' command which you have just inserted, deactivating that line by inserting blank spaces in its first and second columns.

In the end of this rather long tutorial it is relevant to remind the user of its main contents. In Part 2 of the tutorial,

You have thus advanced significantly towards becoming one of the second of the two kinds of PHOENICS user described in the first part of this tutorial, viz. one whose superior knowledge assists less knowledgable colleagues.