TALK=T;RUN( 1, 1)
TALK=T;RUN( 1, 1)
** LOAD(x220) from the x Input Library
PHOTON USE
p ; ; ; ; ;
msg Computational Domain:
gr z 1
use patgeo
msg Press Any Key to Continue...
pause
cl;gr ou z 1
msg Contours of Radiation Temperature T3:
con t3 k 1 fi;0.005
pause
cl;gr ou z 1
msg Contours of x-direction Radiation flux QRx:
con qrx k 1 fi;0.005
pause
cl;gr ou z 1
msg Contours of y-direction Radiation flux QRy:
con qry k 1 fi;0.005
pause
cl;gr ou z 1
msg Contours of Wgap:
con Wgap k 1 fi;0.005
pause
cl;gr ou z 1
msg Contours of Wdis:
con Wdis k 1 fi;0.005
pause
cl;gr ou z 1
msg Contours of LTLS:
con LTLS k 1 fi;0.005
msg Press E to exit PHOTON ...
ENDUSE
**************************************************************
#cls
TEXT(IMMERSOL: empty box; fixed wall temps
TITLE
REAL(TEM_N,TEM_S,TEM_W,TEM_E)
REAL(SIGMA,ABSORB,SCAT,MAT_W,WSL,SIZX,SIZY)
REAL(QRAD_N,QRAD_S,QRAD_W,QRAD_E)
REAL(EMIS_W,EMIS_E,EMIS_S,EMIS_N)
REAL(FEW,FNW)
** data settings
MAT_W=100;ABSORB= 1.0E-5; SCAT= 1.0E-5
MESG( A square enclosure contains a non-participating medium.
MESG(
MESG( The temperature of the west wall is 310 K.
MESG( The temperatures of the other walls are 300 K.
MESG(
MESG( All walls are assumed black, i.e.to have emissivity=1.0 .
MESG(
MESG( The task is to calculate net heat fluxes on all walls
MESG(
MESG( The problem was studied by:
MESG( A. Sanches and T.F. Smith, Surface radiation exchange for
MESG( 2-D rectangular enclosures using discrete ordinates method,
MESG( Journal of Heat Transfer, Vol. 144, pp.465-472, 1992.
MESG(
#pause
WSL=0.01;SIZX=1.0;SIZY=1.0
** opportunity to change size of box
MESG(change original size of box ? y/N
readvdu(ans,char,N)
IF(:ANS:.EQ.Y) THEN
MESGA(x_size of box = :sizx:. OK? If not insert new value..
READVDU(sizx,real,:sizx:)
MESGA(y_size of box = :sizy:. OK? If not insert new value..
READVDU(sizy,real,:sizy:)
ENDIF
#pause
** opportunity to change aspect ratio of box
MESGA(width of box = :sizx:. OK? If not insert new value..
READVDU(sizx,real,:sizx:)
TEM_S=300.;TEM_N=300.;TEM_W=310.;TEM_E=300.
EMIS_W= 1.0;EMIS_E= 1.0;EMIS_N= 1.0;EMIS_S= 1.0
SIGMA = 5.6697E-8
#cls
MESGA( Change wall temperature ? y/N
readvdu(ans,char,N)
IF(:ANS:.EQ.Y) THEN
mesg(South wall Tempereture = :TEM_S:. OK? If not insert new value..
readvdu(TEM_S,real,:TEM_S:)
mesg(North wall Tempereture = :TEM_N:. OK? If not insert new value..
readvdu(TEM_N,real,:TEM_N:)
mesg(East wall Tempereture = :TEM_E:. OK? If not insert new value..
readvdu(TEM_E,real,:TEM_E:)
mesg(West wall Tempereture = :TEM_W:. OK? If not insert new value..
readvdu(TEM_W,real,:TEM_W:)
ENDIF
#cls
MESGA( Change wall emissivities ? y/N
readvdu(ans,char,N)
IF(:ANS:.EQ.Y) THEN
mesg(South wall emissivity = :emis_s:. OK? If not insert new value..
readvdu(emis_s,real,:emis_s:)
mesg(North wall emissivity = :emis_n:. OK? If not insert new value..
readvdu(emis_n,real,:emis_n:)
mesg(East wall emissivity = :emis_e:. OK? If not insert new value..
readvdu(emis_e,real,:emis_e:)
mesg(West wall emissivity = :emis_w:. OK? If not insert new value..
readvdu(emis_w,real,:emis_w:)
ENDIF
IF ((EMIS_W.EQ.1.0).AND.(EMIS_E.EQ.1.0).AND.(EMIS_N.EQ.1.0).AND.(EM$
IS_S.EQ.1.0)) THEN
#pause
MESG(The analytical solution is as follows:
MESG(--------------------------------------
CALL ANALYTIC
ENDIF
#pause
MESGA(The nett radiation fluxes computed by IMMERSOL are best shown
MESG(by way of the Cellav values in the PROFIL print-out of:
MESG(QRX near the west and east walls, and of:
MESG(QRY near the north and south walls.
MESGA(These are to be found at the bottom of the RESULT file.
MESGA(*******************************************************
MESGA(This is a severe test for IMMERSOL which performs best for
MESG(wider opposite-facing surfaces.
MESGA(As expected, IMMERSOL predicts fluxes which are no better than
MESG(of the right order of magnitude.
MESGA(Changing the aspect ratio (SIZX/SIZY), or removal of the north
MESG(and south walls, brings better agreement.
#pause
** opportunity to remove north and south walls
MESGA(Remove north and south walls? y/N
READVDU(ans,char,n)
** opportunity change grid size
integer(ngridx,ngridy)
ngridx=40;ngridy=40
mesg(cells in x-direction = :ngridx:. OK? If not insert new value..
readvdu(ngridx,int,ngridx)
ngridx
mesg(cells in y-direction = :ngridy:. OK? If not insert new value..
readvdu(ngridy,int,ngridy)
ngridy
GROUP 3. X-direction grid specification
NREGX= 3; IREGY= 1; GRDPWR(X,1,WSL,1.0)
IREGX= 2; GRDPWR(X,ngridx,SIZX,1.0);
IREGX= 3; GRDPWR(X,1,WSL,1.0)
GROUP 4. Y-direction grid specification
NREGY= 3; IREGY= 1; GRDPWR(Y,1,WSL,1.0)
IREGY= 2; GRDPWR(Y,ngridy,SIZY,1.0);
IREGY= 3; GRDPWR(Y,1,WSL,1.0)
GROUP 7. Variables stored, solved & named
Energy transport in solids by conduction & that within the
the gap between solids by radiation
SOLVE(T3)
Heat transport coefficient within the gas by radiation
DISWAL;STORE(WGAP)
Selection of material:Material flags are stored in PRPS
STORE(PRPS)
storage for radiative heat flux & wall emsissivity
STORE(QRX,QRY,EMIS,SCAT)
GROUP 8. Terms (in differential equations) & devices
T3 Equation includes only diffision terms
TERMS(T3,N,N,Y,N,Y,N)
GROUP 9. Properties of the medium (or media)
!!!! the value of PRNDTL(T3) no longer has significance
GROUP 11. Initialization of variable or porosity fields
INIADD= F;FIINIT(PRPS)=0.0;FIINIT(T3)=0.5*(TEM_W+TEM_E)
FIINIT(EMIS)=ABSORB; FIINIT(SCAT)=SCATT
***Solid walls:
Material selection for solid boundaries
111, Steel at 27 deg c (MATH_W)
IF(:ANS:.EQ.Y) THEN
GOTO JUMP1
ENDIF
PATCH(NORTH,INIVAL, 2,NX-1, NY,NY,1,1, 1,LSTEP)
PATCH(SOUTH,INIVAL, 2,NX-1, 1,1, 1,1, 1,LSTEP)
INIT(NORTH,PRPS,0.0,MAT_W); INIT(SOUTH,PRPS,0.0,MAT_W)
INIT(NORTH,EMIS,0.0,EMIS_N); INIT(SOUTH,EMIS,0.0,EMIS_S)
LABEL JUMP1
PATCH(WEST, INIVAL, 1,1, 2, NY-1, 1,1, 1,LSTEP)
PATCH(EAST, INIVAL, NX,NX,2, NY-1, 1,1, 1,LSTEP)
INIT(WEST,PRPS,0.0,MAT_W); INIT(EAST,PRPS,0.0,MAT_W)
INIT(WEST,EMIS,0.0,EMIS_W); INIT(EAST,EMIS,0.0,EMIS_E)
GROUP 13. Boundary conditions and special sources
Wall temperatures
REAL(CONST)
CONST=1.e+5
PATCH(WL1,VOLUME, 1,1, 2,NY-1, 1,1, 1,LSTEP)
PATCH(WL2,VOLUME, NX,NX, 2,NY-1, 1,1, 1,LSTEP)
COVAL(WL1,T3,CONST,TEM_W)
COVAL(WL2,T3,0.5*CONST,TEM_E)
IF(:ANS:.EQ.Y) THEN
GOTO JUMP2
ENDIF
PATCH(WL3,VOLUME, 2,NX-1, 1,1, 1,1,1,LSTEP)
PATCH(WL4,VOLUME, 2,NX-1, NY,NY, 1,1, 1,LSTEP)
COVAL(WL3,T3,CONST,TEM_S)
COVAL(WL4,T3,CONST,TEM_N)
LABEL JUMP2
GROUP 15. Termination of sweeps
LSWEEP= 600; TSTSWP= -1
GROUP 16. Termination of iterations
SELREF=T; RESFAC= 0.01
LITER(T3)=100
RELAX(T3,LINRLX,0.5)
Group 18. Limits
VARMAX(QRY ) = 3000. ;VARMIN(QRY ) =-3000.
VARMAX(QRX ) = 3000. ;VARMIN(QRX ) =-3000.
GROUP 22. Spot-value print-out
IXMON=NX/2; IYMON= NY/2
PATCH(T3WEST,PROFIL,1,1,2,NY-2,1,1,1,LSTEP)
COVAL(T3WEST, T3,0,0)
PATCH(QRYNORTH,PROFIL,2,NX-1,NY-2,NY-2,1,1,1,LSTEP)
COVAL(QRYNORTH, QRY,0,0)
PATCH(QRYSOUTH,PROFIL,2,NX-1,2,2,1,1,1,LSTEP)
COVAL(QRYSOUTH, QRY,0,0)
PATCH(QRXWEST,PROFIL,2,2,2,NY-2,1,1,1,LSTEP)
COVAL(QRXWEST, QRX,0,0)
PATCH(QRXEAST,PROFIL,NX-2,NX-2,2,NY-2,1,1,1,LSTEP)
COVAL(QRXEAST, QRX,0,0)
IPROF=1;ORSIZ=0.2
OUTPUT(QRY ,Y,N,N,N,N,N)
OUTPUT(QRX ,Y,N,N,N,N,N)
LIBREF = 220
NXPRIN=1;NYPRIN=1
it is necessary to place the expected values here, not in the
calling q1
EX(EMIS)= 9.071E-02
EX(QRY )= 2.034E+01
EX(QRX )= 3.215E+01
EX(WGAP)= 5.800E-01
EX(LTLS)= 3.191E-02
EX(WDIS)= 1.319E-01
EX(T3 )= 3.026E+02
ENDMAIN !
SUBROUTINE ANALYTIC
*WALL EMISSIVE POWERS
REAL(EB_N,EB_S,EB_W,EB_E)
EB_N=SIGMA*TEM_N**4;EB_S=SIGMA*TEM_S**4
EB_W=SIGMA*TEM_W**4;EB_E=SIGMA*TEM_E**4
** VIEW FACTOR CALCULATIONS (crossed-string method)
REAL(FS_E,FS_W,FS_N,FN_E,FN_W,FN_S,FE_S,FE_N,FE_W,FW_S,FW_N,FW_E)
*SOUTH WALL
FS_E=((SIZX+SIZY)-(SIZX**2.+SIZY**2.)**0.5)/(2.*SIZX)
FS_W=FS_E
FS_N=1.-FS_E-FS_W
*NORTH WALL
FN_E=FS_E
FN_W=FS_W
FN_S=FS_N
*EAST WALL
FE_S=(SIZX/SIZY)*FS_E
FE_N=FE_S
FE_W=1.-FE_N-FE_S
*WEST WALL
FW_N=FE_N
FW_E=FE_W
FW_S=FE_S
MESGA(View factor calculations? y/N
READVDU(ans,char,n)
IF(:ANS:.EQ.Y) THEN
MESG( View factor between south and east, west, & north, :FS_E:$
,:FS_W:,:FS_N:
MESG( View factor between north and east, west, & south, :FN_E:$
,:FN_W:,:FN_S:
MESG( View factor between east and south, north & west, :FE_S:,$
:FE_N:,:FE_W:
MESG( View factor between west and south, north & east, :FW_S:,$
:FW_N:,:FW_E:
ENDIF
#pause
** HEAT FLUX (W/m**2)CALCULATIONS
*SOUTH WALL
QRAD_S=FS_E*(EB_S-EB_E)+FS_W*(EB_S-EB_W)+FS_N*(EB_S-EB_N)
MESG( Exact QRAD at south wall is :QRAD_S: W/m**2
QRAD_N=FN_E*(EB_N-EB_E)+FN_W*(EB_N-EB_W)+FN_S*(EB_N-EB_S)
MESG( Exact QRAD at north wall is :QRAD_N: W/m**2
QRAD_E=FE_S*(EB_E-EB_S)+FE_N*(EB_E-EB_N)+FE_W*(EB_E-EB_W)
MESG( Exact QRAD at east wall is :QRAD_E: W/m**2
QRAD_W=FW_S*(EB_W-EB_S)+FW_N*(EB_W-EB_N)+FW_E*(EB_W-EB_E)
MESG( Exact QRAD at west wall is :QRAD_W: W/m**2
MESG(
#pause
MESGA(Conservation energy check? y/N
READVDU(ans,char,n)
IF(:ANS:.EQ.Y) THEN
MESG(For rectangular enclosure having aspect ratio other
MESG(than 1 it is conventional to report heat flux per unit
MESG(length in order to check conservation of energy
MESG(South wall:
QRAD_S=QRAD_S*SIZX
MESG(Length of south wall is :sizx: m
MESG(QRAD at south wall per unit length is :QRAD_S: W/m
MESG(
MESG(North wall;
QRAD_N=QRAD_N*SIZX
MESG(Length of north wall is :sizx: m
MESG(QRAD at north wall per unit length is :QRAD_N: W/m
MESG(
MESG(East wall:
QRAD_E=QRAD_E*SIZY
MESG(Length of east wall is :sizy: m
MESG(QRAD at east wall per unit length is :QRAD_E: W/m
MESG(
MESG(East wall:
QRAD_W=QRAD_W*SIZY
MESG(Length of west wall is :sizy: m
MESG(QRAD at west wall per unit length is :QRAD_W: W/m
MESG(********Check for conservation of energy*********
MESG(The sum of all heat transfer rates must vanish)
REAL(QRADT)
QRADT=QRAD_S+QRAD_N+QRAD_W+QRAD_E
MESG(Summation of wall heat trasfer rates is :QRADT:
ENDIF
ENDSUB
STOP