GROUP 1. Run title and other preliminaries
TEXT(3D Steady Heat Conduction In Cube
#cls
TITLE
DISPLAY
/------/------/------/------/------/
/temperature fixed at 1.0 ----> X /|
/------/------/------/------/------/ /
/ / / / / /|/|
/------/------/------/------/------/ / /
/ / / / / /|/|/|
/------/------/------/------/------/ / / /
/ / / / / /|/|/|/|
north |------|------|------|------|------| / / / /
^ | | | | | |/|/|/|/|
| |------|------|------|------|------| / / / / high
| | | | | | |/|/|/|/ ^
y |------|------|------|------|------| / / / /
| | | | | | |/|/|/ z
| |------|------|------|------|------| / / /
| | | | | | |/|/ /
| |------|------|------|------|------| / /
| | X <-- temperature fixed at zero |/ /
south |------|------|------|------|------| low
/west----------------x----->east
#pause
The following PHOTON "use" file enables the temperature
distribution to be visualized.
PHOTON USE
ext;;;;
view 5 2 -1
msg a cube with fixed temperatures at opposite corners
msg press enter
gr x m; gr y m; gr z 1
pause;gr off;gr ou x m;gr ou y m;gr ou z 1;red
msg the temperature for z=1
msg press enter
con temp z 1 fi;0.01
con temp z 1 ;int 40
gr ou x m;gr ou y m;gr ou z 1;
pause;red
msg the temperature for y=ny
msg press enter
con temp y m fi;0.01
con temp y m ;int 40
gr ou x m;gr ou y m;gr ou z 1;
pause;red
msg the temperature for x=nx
msg press enter
con temp x m fi;0.01
con temp x m ;int 40
gr ou x m;gr ou y m;gr ou z 1;
pause;red
msg Press e to END, or ? for help
enduse
ENDDIS
#pause
integer(lit)
REAL(XLENGTH,YLENGTH,ZLENGTH)
XLENGTH=1.0;YLENGTH=1.0;ZLENGTH=1.0
domain size and grid
NX=10; NY=10; NZ=10; xulast=xlength; yvlast=ylength; zwlast=zlength
#unigrid
label top
GROUP 7. Variables stored, solved & named
**Choose first-phase enthalpy (H1) as dependent variable
and activate the whole-field elliptic solver
IF(.NOT.:ANS:.EQ.Y) THEN
SOLUTN(H1,Y,Y,Y,N,N,N);NAME(H1)=TEMP
ENDIF
** Store BLOK, the marker variable needed for the activation
of the block-correction feature in the linear-equation
solver.
STORE(BLOK)
** Indicate that it is to variable 14 that the block-
correction feature will be applied
IVARBK=14
GROUP 8. Terms (in differential equations) & devices
**For pure conduction, cut out built-in source and convection
terms
TERMS(TEMP,N,N,Y,N,Y,Y)
GROUP 9. Properties of the medium (or media)
**Thermal conductivity will be ENUL*RHO1/PRNDTL(TEMP), so :
ENUL=1.0;RHO1=1.0;PRNDTL(TEMP)=1.0
GROUP 11. Initialization of variable or porosity fields
** Define eight block-correction blocks by ascribing BLOK
values to chosen locations by way of:-
FIINIT, for block 1, and PATCH and INIT for the other
blocks
Note that INIADD is = F, so that the values set by INIT
replace the value set by FIINIT. The only region which
retains that value is that in the east, north, high
corner.
INIADD=F;FIINIT(BLOK)=1.0
West, south, low corner
PATCH(BLOK2,INIVAL,1,NX/2,1,NY/2,1,NZ/2,1,LSTEP)
INIT(BLOK2,BLOK,0.0,2.0)
East, north, low corner
PATCH(BLOK3,INIVAL,NX/2+1,NX,1,NY/2,1,NZ/2,1,LSTEP)
INIT(BLOK3,BLOK,0.0,3.0)
West, south, low corner
PATCH(BLOK4,INIVAL,1,NX/2,NY/2+1,NY,1,NZ/2,1,LSTEP)
INIT(BLOK4,BLOK,0.0,4.0)
East, north, low corner
PATCH(BLOK5,INIVAL,NX/2+1,NX,NY/2+1,NY,1,NZ/2,1,LSTEP)
INIT(BLOK5,BLOK,0.0,5.0)
West, south, high corner
PATCH(BLOK6,INIVAL,1,NX/2,1,NY/2,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK6,BLOK,0.0,6.0)
East, south, high corner
PATCH(BLOK7,INIVAL,NX/2+1,NX,1,NY/2,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK7,BLOK,0.0,7.0)
West, north, high corner
PATCH(BLOK8,INIVAL,1,NX/2,NY/2+1,NY,NZ/2+1,NZ,1,LSTEP)
INIT(BLOK8,BLOK,0.0,8.0)
GROUP 13. Boundary conditions and special sources
**Corner at IX=IY=IZ=1
PATCH(COLD,CELL,1,1,1,1,1,1,1,1)
**Fix temperature to zero
COVAL(COLD,TEMP,1.E2,0.0)
**Corner at IX=NX, IY=NY, IZ=NZ
PATCH(HOT,CELL,NX,NX,NY,NY,NZ,NZ,1,1)
**Fix temperature to 1.0
COVAL(HOT,TEMP,1.E2,1.0)
GROUP 15. Termination of sweeps
RESREF(TEMP)=1.E-6;LSWEEP=2
GROUP 16. Termination of iterations
**Terminate iterations when 500 iterations have been performed.
(The minus sign is a signal to activate progress-of-solution
print-out on the screen at a sweep frequency dictated by the
value of TSTSWP.)
LITER(TEMP)=-500
** Set the over-relaxation factor for the linear-equation
solver.
OVRRLX=1.7
** Set the frequency of application of the block=correction
feature to once per iteration
ISOLBK=1
** Set the frequencies of application of the one-dimensional
correction features in the linear-equation solver to once
per iteration for each direction.
ISOLX=1;ISOLY=1;ISOLZ=1
************************************************************
Group 17. Relaxation
RELAX(TEMP,LINRLX,1.0)
GROUP 21. Print-out of variables
**Print fields of temperature
OUTPUT(TEMP,Y,N,N,N,N,N)
GROUP 22. Spot-value print-out
IXMON=NX/2+1;IYMON=NY/2+1;IZMON=NZ/2+1;UWATCH=T
GROUP 23. Field print-out and plot control
NXPRIN=NX/5;NYPRIN=NY/5;NZPRIN=NZ/5
**Plot a profile along the line IX=nx/2,IZ=nz/2
PATCH(YLINE,PROFIL,NX/2,NX/2,1,NY,NZ/2,NZ/2,1,1)
PLOT(YLINE,TEMP,0.0,0.0);PLOT(YLINE,BLOK,0.0,0.0)
**Plot a contour diagram for the plane IX=nx/2
PATCH(XPLANE,CONTUR,NX/2,NX/2,1,NY,1,NZ,1,1)
**Let the diagram have 20 temperature intervals
PLOT(XPLANE,TEMP,0.0,20.0);PLOT(XPLANE,BLOK,0.0,20.0);ICHR=2
GROUP 24. Dumps for restarts
Interesting variants include the following:-
1. By putting FIINIT(TEMP)= ... in GROUP 11, see how the
initial guess influences solution time, final result etc.
2. By changing NX, NY and NZ, see how increasing or decreasing
the fineness of spatial subdivision influences solution
time, results, etc.
3. By changing XLENGTH, YLENGTH and ZLENGTH, see how the
solution changes when the block is made long and thin,
short and fat, etc.
4. By changing PRNDTL(TEMP), see how the value of the
conductivity changes the solution.
5. By changing the third, fourth, and further arguments of
PATCH(HOT,... and PATCH(COLD,..., explore how changing the
boundary-condition locations influence the solution.
Also, introduce more PATCHes at which boundary conditions
are provided; and display the results in different ways by
introducing more PATCHes for profile and contour plots.
6. By changing the COVAL coefficient from 1.E2 to FIXFLU,
explore the effect of having fixed-flux rather than fixed-
value boundary conditions; but remember that if all the
boundary conditions are of fixed-flux type, the solution is
not uniquely defined.
7. By changing the values of OVRRLX, ISOLBK, ISOLX, ISOLY,
ISOLZ, see how these solution-control parameters
influence the solution and the computer time and number
of iterations needed to attain it.
8. By changing the number and location of the blocks used in
the block-correction feature, establish the influences of
factors on the accuracy and computer time.