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.
Instructions will be provided concerning:
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:
You are now in a position to edit the file.
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.
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.
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.
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 5PHOENICS 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.
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?
nx=5When 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.
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.
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.
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
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.
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.
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.0Both 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 save25endRunning 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...
Although most obvious things are not always realized, we shall further show what exactly you can do in the land of PHOENICS.
PRNDTL(TEM1)=1.0E-05/43.0If 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.
PRNDTL(TEM1)=-43.0whereby 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.
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
PRNDTL(TEM1)=1.0Then 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:
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.
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
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:
COVAL(HEATER, TEM1, FIXVAL, 1.000000E+01)rather than the next line,
COVAL(HEATER, TEM1, 2.000000E+10, 1.000000E+01)which does not.
> FIXFLU = FIXFLU > FIXVAL = FIXVALand 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.
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:
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))/KONDan 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.
When the case is run, it will be found that:
************************************************************
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.
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:
Interactive statements to be placed below here Interactive statements to be placed above herethe 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.
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.
This can easily be achieved by:
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,