PHOTON USE
  ext
  phi
 
 
 
  msg                3D UNSTEADY HEAT CONDUCTION
  msg
  msg        The configuration:
  gr ou x 4;gr ou x 8;gr ou z 1;gr ou z max
  msg
  msg Press  to continue
  pause
  msg
  msg        Temperature distribution:
  cl;red
  con temp x max fi;.001
  con temp z max fi;.001
  msg
  msg Press  to continue
  pause
  view z;norm;gr ou z max
  msg        Temperature distribution viewed from 'Z' direction:
  msg
  msg Press e to END
  ENDUSE
 
    GROUP 1. Run title and other preliminaries
TEXT(3D Unsteady Heat Cond. Thick Pipe 
TITLE
mesg(PC/486-50 time last reported as appx.40 sec.
  DISPLAY
    This run is concerned with heat conduction. It differs from the
  preceding cases in that a transient is considered. Another
  difference is that a cylindrical-polar coordinate grid is
  employed. A schematic diagram of the geometry follows;
  note, however, that the ratio of the wall thickness to
  the inner diameter in the diagram is not to scale.
 
   ^     .-----.--------------------------.
   |   /***.-.***\                         \
   |r /**/     \**\                         \
   | |**|       |**|                         |
     |**|       |**|                         |
      \**\     /**/- - - - - - - - - - - - -/
       \***`|'***/   temp. fixed to zero   /
         `--|--'--------------------------'
            |       z---->
     Fixed heat flux
     at inner surface
  ENDDIS
 
    GROUP 2. Transience; time-step specification
   **1-second duration, in 10 equal intervals
STEADY=F;IREGT=1; GRDPWR(T,10,1.0,1.0)
 
    GROUP 3. X-direction grid specification
   **Polar coordinates; total angle is 2*pi
CARTES=F;NREGX=3
IREGX=1; GRDPWR(X,6,2.0*3.14159*3/8,1.0)
IREGX=2; GRDPWR(X,4,2.0*3.14159*2/8,1.0)
IREGX=3; GRDPWR(X,6,2.0*3.14159*3/8,1.0)
 
    GROUP 4. Y-direction grid specification
   **Pipe has 0.025 m internal radius and 0.05 wall thickness.
IREGY=1; GRDPWR(Y,5,0.05,1.0); RINNER=0.025
 
    GROUP 5. Z-direction grid specification
   **Pipe is 1.0 m long in z-direction, with 5 equal intervals.
IREGZ=1; GRDPWR(Z,5,1.0,1.0)
 
    GROUP 7. Variables stored, solved & named
   **Choose first-phase enthalpy (H1) as dependent variable
     and activate the whole-field elliptic solver.
SOLUTN(H1,Y,Y,Y,N,N,N)
     Rename it as TEMP, as mnemonic for temperature
NAME(H1)=TEMP
 
    GROUP 8. Terms (in differential equations) & devices
   **Cut out source and convection terms for pure heat conduction.
TERMS(TEMP,N,N,Y,Y,Y,Y)
     Apply 1d correction in all directions
ISOLX=1;ISOLY=1;ISOLZ=1
 
    GROUP 9. Properties of the medium (or media)
   **Thermal conductivity will be ENUL*RHO1/PRNDTL(TEMP), so
     with default RHO1=1.0 :
REAL(COND);COND=0.1234;ENUL=1.0;PRNDTL(TEMP)=1.0/COND
 
    GROUP 11. Initialization of variable or porosity fields
   ** Pipe starts at zero temperature.
FIINIT(TEMP)=0.0
 
    GROUP 13. Boundary conditions and special sources
   **Boundary conditions at x=0. and 2 pi are "cyclic".
XCYCLE=T
   **Inner surface, uniform fixed flux
PATCH(Y1S,SOUTH,#1,#NREGX,#1,#1,#1,#NREGZ,1,#NREGT);
COVAL(Y1S,TEMP,FIXFLU,0.1)
   **One quarter of the outer surface is held at zero temperature.
WALL(OUTER,NORTH,#2,#2,#NREGY,#NREGY,#1,#NREGZ,#1,#NREGT)
COVAL(OUTER,TEMP,COND,0.0)
   **One end face is also cooled
WALL(END,HIGH,#1,#NREGX,#1,#NREGY,#NREGZ,#NREGZ,#1,#NREGT)
COVAL(END,TEMP,COND,0.0)
 
    GROUP 16. Termination of iterations
   **Terminate iterations when average correction falls to
     1.0E-7 or when 200 iterations have been performed. The
     minus sign is a signal to activate progress-of-solution
LITER(TEMP)=-200;OVRRLX=1.2
 
    GROUP 21. Print-out of variables
   **Print fields of temperature every 5th time step and 4th x-plane
OUTPUT(TEMP,Y,N,N,N,Y,Y);NTPRIN=5;NXPRIN=4;YZPR=T
 
    GROUP 22. Spot-value print-out
IZMON=NZ/2+1
SPEDAT(SET,GXMONI,TRANSIENT,L,F) 
    GROUP 23. Field print-out and plot control
   **Plot variation with time of temperature at IX=4,IY=3,IZ=3
PATCH(X4Y3,PROFIL,4,4,3,3,3,3,1,LSTEP)
   **Set the scale limits to 0.0 and 0.05
PLOT(X4Y3,TEMP,0.0,0.0)
   **Plot a contour diagram for the plane IX=1
PATCH(X1P,CONTUR,1,1,1,NY,1,NZ,1,LSTEP)
   **Let the diagram have 20 temperature intervals
PLOT(X1P,TEMP,0.0,20.0)
   **Plot a contour diagram for the plane IX=4.
PATCH(X4P,CONTUR,4,4,1,NY,1,NZ,1,LSTEP);PLOT(X4P,TEMP,0.0,20.0)
   **Plot a contour diagram for the plane IZ=3; a polar plot is
     selected by making the 3rd of the PLOT command argument non
     zero.
PATCH(Z3P,CONTUR,1,NX,1,NY,3,3,1,LSTEP);PLOT(Z3P,TEMP,1.0,20.0)
PATCH(Z3P2,CONTUR,4,NX,1,NY,3,3,1,LSTEP);PLOT(Z3P2,TEMP,1.0,20.0)
PATCH(Z3P3,CONTUR,2,4,1,NY,3,3,1,LSTEP);PLOT(Z3P3,TEMP,1.0,20.0)
PATCH(Z3P1,CONTUR,1,3,1,NY,3,3,1,LSTEP);PLOT(Z3P1,TEMP,1.0,20.0)
PATCH(Z3P4,CONTUR,7,NX,2,NY-1,3,3,1,LSTEP);PLOT(Z3P4,TEMP,1.0,20.0)
   **Plot a contour diagram for part of the IZ=3 plane, and
     tabulate the values at the same time
PATCH(TABZ3P4,CONTUR,1,NX,2,NY-1,3,3,1,LSTEP)
PLOT(TABZ3P4,TEMP,1.0,20.0)
 
mesg( Press  to continue
readvdu(ans,char,n)
 
do ii=1,5
+  mesg(
enddo
mesg( Initial data that can be changed:
mesga( Pipe internal radius is set to 0.025 m
mesg( Pipe wall thickness is set to 0.05 m
mesg( Conductivity is set to 0.1234 W/m/deg
mesg( Heat flux through inner surface is 0.1 W/m**2
mesga( Do you want to change settings (y/n)? (Default n)
readvdu(ans,char,n)
 
if(:ans:.eq.y) then
+ real(rt1)
+ do ii=1,5
+   mesg(
+ enddo
+ mesga( Pipe internal radius is 0.025 m. OK?
+ mesg(If not, insert new value
+ readvdu(rt1,real,0.025)
+ RINNER=rt1
 
+ do ii=1,5
+   mesg(
+ enddo
+ mesga( Pipe wall thickness is 0.05 m. OK?
+ mesg(If not, insert new value
+ readvdu(rt1,real,0.05)
+ IREGY=1; GRDPWR(Y,5,rt1,1.0)
 
+ do ii=1,5
+   mesg(
+ enddo
+ mesga( Conductivity is 0.1234 W/m/deg. OK?
+ mesg(If not, insert new value
+ readvdu(COND,real,0.1234)
+ PRNDTL(TEMP)=1.0/COND
 
+ do ii=1,5
+   mesg(
+ enddo
+ mesga( Heat flux through inner surface is 0.1 W/m**2. OK?
+ mesg( If not, insert new value
+ readvdu(rt1,real,0.1)
+ COVAL(Y1S,TEMP,FIXFLU,rt1)
endif
liter(temp)=200
selref=t; resfac=1.e-2;TSTSWP=-1