TALK=T;RUN(1,1) 492.htm

  PHOTON USE
  ext;;;;
 
  gr ou z 1; use patgeo
  msg boundary condition patches. Press RETURN
  pause
  msg temperature contours. Press RETURN
  con tmp1 y m fi;0.001;con tmp1 x m fi;0.001;con tmp1 z m fi;0.001
  gr ou y m;gr ou z m
  msg press RETURN for view x
  pause;  con off;  view x;  con tmp1 x m fi;0.001
  con tmp1 y m fi;0.001;  gr ou y m;  gr ou x m;  gr ou z 1
  msg This is view x. Press RETURN for view y
  pause;  con off;  view y;  con tmp1 y m fi;0.001
  msg This is view y. Press RETURN for view z
  pause;  con off; view z;  con tmp1 z m fi;0.001;  gr ou z m
  msg This is view z. Press RETURN for view x and velocity vectors
  pause;  con off;  view x
  msg velocity vectors. Type menu for menu and further possibilities
  msg Press e to end. Otherwise enter photon-readable commands
  vec x m sh;  gr ou y m;  gr ou x m;  gr ou z 1
  ENDUSE
  Group1. Run title and other preliminaries
  
#cls
TEXT(Idealised Gas-Turbine Combustion Chamber
LIBREF=492   
readq1=t
  DISPLAY
  READQ1_BEGIN
  The following text is provided as an example of what, by use of 
  the readq1 command, can be transmitted to EARTH for writing near 
  the top of the RESULT file. Such text should leave columns 1 and 2
  blank; and it should not extend beyond the 68th column.
  
    ****************************************************************
    * This library case dates from the earliest days of PHOENICS,  *
    * when Professor WU Chung-Hua ('turbomachinery Wu', who had    *
    * returned to China from the USA) visited CHAM in 1982.        *
    *                                                              *
    * The configuration of the combustion chamber, and its being   *
    * supplied with premixed fuel vapour and air, was proposed by  *
    * Professor Wu's accompanying assistant.                       *
    ****************************************************************
  READQ1_END
TITLE  
  The shape of the combustion chamber is as shown.
  Pre-mixed fuel vapour and air enter near the axis on left.
  Secondary and dilution air enter through holes in outer wall.
  A 36-degree sector is simulated.
 
                 secondary ox. inlets   dilution inlet
           ____________ 1,2 _______________     ________
  blocked /
  region /
        /
       /                                             outlet
       |
       |_____   _
  fuel-ox. inlet|            Symmetry axis
       -- - -- -| -- - -- -  -- - -- - -- - -- - -- - -- -
 
  The flow is turbulent; the Simple Chemical Reaction Scheme is
  used; and the reaction-rate is physically controlled by means of
  the Eddy-Breakup Model.
#pause
  ENDDIS
  inform1begin
   ** Declarations of non-standard variables
CHAR(RADE,INFORM) ! for radiation and inform options
REAL(TWALL) 
INTEGER(NYB,NZB,NYIN,NZIN,NZJ1,NZJ2) ! limits for blockages etc
REAL(PI,TKEIN,EPSIN)
REAL(CHAMBR)
REAL(VINF,VIN,GVIST,RHOIN,FINF,HAIRIN,SECIN,DILIN,TAIR,TFUEL)
REAL(HFU,STOIC,FSTOI,CEBU,CPFU,CPPR,HINF,WAIR,WFU,WPR,CPAIR)

   ** Settings of non-standard variables
CHAMBR=0.065 ! Radius of combustion chamber 
PI=3.1415927 ! Pythagoras's constant
   ** Geometric specifications
NYB=2;NZB=3;NYIN=3;NZIN=3;NZJ1=7;NZJ2=10

    
WAIR=29.0; WFU=16.0; WPR=28.0 ! molecular weights
STOIC=17.24  ! stoichiometric ratio, kg air per 1 kg fuel
FSTOI=1./(1.+STOIC) ! fuel fraction in unburned stoich. mixture
CPAIR=1.5E3; CPPR=1.5E3; CPFU=1.5E3 ! specific heats
HFU=4.9E7    ! heat of reaction
PRESS0=8.0E5 ! reference pressure
TWALL=773    ! temperature of the wall, set to 500K
CEBU=1.0     ! eddy-break-up constant

  data relating to input streams
TAIR=773.0  ! temperature of in-flowing secondary and dilution air
SECIN=0.008157 ! secondary air mass flow rate
DILIN=0.0145   ! diluent-stream air mass flow rate
HAIRIN=CPAIR*TAIR ! enthalpy
VIN=-40.0   ! velocity 
RHOIN=3.606 ! density
FINF=FSTOI  ! mixture fraction of in-flowing fuel-air mixture
TFUEL=773.0 ! temperature 
VINF=100.0  ! velocity 
HINF=CPFU*TFUEL+FINF*HFU ! enthalpy
if(:inform:xx.eq.xx) then
MESG( Use In-Form in place of GRND settings Y/n ?
READVDU(INFORM,CHAR,Y) ! allow user to decide
endif
  inform1end
MESG( Include Radiation Equations y/N ?
READVDU(RADE,CHAR,N) ! allow user to decide
   
   ** End of settings of non-standard variables
   
    
    GROUP 3. X-direction grid specification
NX=6; XULAST=2.*PI/10.0; CARTES=F
XFRAC(1)=-6.0;XFRAC(2)=0.166666
    GROUP 4. Y-direction grid specification
NY=10;YVLAST=1.0
YFRAC(1)=-5.0;YFRAC(2)=0.0065;YFRAC(3)=1.0;YFRAC(4)=0.005
YFRAC(5)=1.0;YFRAC(6)=0.008;YFRAC(7)=3.0;YFRAC(8)=0.0065
    GROUP 5. Z-direction grid specification
NZ=13;ZWLAST=1.0
   ** Full length of chamber = 0.245
ZFRAC(1)=-2.0;ZFRAC(2)=0.0125;ZFRAC(3)=1.0;ZFRAC(4)=0.004
ZFRAC(5)=1.0;ZFRAC(6)=0.011;ZFRAC(7)=1.0;ZFRAC(8)=0.01
ZFRAC(9)=1.0;ZFRAC(10)=0.0235;ZFRAC(11)=7.0;ZFRAC(12)=0.0245
    GROUP 7. Variables stored, solved & named
IF(:RADE:.EQ.Y) THEN
  This set of PIL commands is useful for setting up radiation
  simulations based upon the six-flux model.
  It can be used as a separate file, eg with name: radidata .
  Then it can be "included" in the Q1 in group 7, by insertion of
  the line:
 
  INTRPT(R,RADIDATA)
 
  Declarations and settings
REAL(RADW,ABSORB,SCAT,SIGMA,EMIW,ED2ME)
SIGMA=5.6697E-8;EMIW=0.3
RADW=SIGMA*TWALL**4;ABSORB=0.2;SCAT=0.01
RADIAT(FLUX,ABSORB,SCAT,H1)
FIINIT(RADX)=RADW;FIINIT(RADY)=RADW;FIINIT(RADZ)=RADW
ED2ME=EMIW/(2.0-EMIW)
RELAX(RADX,FALSDT,0.1);RELAX(RADY,FALSDT,1.0)
RELAX(RADZ,FALSDT,1.0)
  ** end of radiation-related settings
ENDIF
 
 
 
  Group 7: variables to be solved and stored
SOLVE(H1,MIXF,FUEL); STORE(OXID,PROD,TMP1,RHO1)
SOLVE(P1,U1,V1,W1);SOLUTN(P1,Y,Y,Y,P,P,N);TURMOD(KEMODL);STORE(ENUT)
 
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,Y,Y,N,Y,N) ! switch off frictional heating
CSG3=CNGR             ! use conjugate-gradient solver

    GROUP 9. Properties of the medium  
RHO1=3GASES ! i.e. GRND6 ;  sets density appropriate to 
            ! mixture of ideal gases
RHO1A=WFU; RHO1B=WAIR; RHO1C=WPR
TMP1=SCRSNONEQ ! i.e. GRND8 ; sets temperature as 
               ! (h1 - fuel*hfu)/ mixture specific heat
TMP2A=FSTOI; TMP2B=HFU
CP1=GRND10 ! sets specific heat as weighted average of the following
CP1A=CPFU; CP1B=CPAIR; CP1C=CPPR
ENUL=3.6E-5/RHOIN
 
    Group 11  Initial values
FIINIT(W1)=10.0
FIINIT(KE)=0.1125*FIINIT(W1)*FIINIT(W1);FIINIT(RHO1)=RHOIN
FIINIT(EP)=0.1643*FIINIT(KE)**1.5/(0.27*CHAMBR)
FIINIT(RHO1)=RHOIN;FIINIT(H1)=HAIRIN;FIINIT(MIXF)=FINF
FIINIT(FUEL)=FINF
 
   ** Porosity fields
 
CONPOR(0.0,VOLUME,1,NX,1,NYB,1,NZB)
CONPOR(0.709,EAST ,1,NX,7,7,1,1); CONPOR(0.709,CELL,1,NX,7,7,1,1)
CONPOR(0.418,NORTH,1,NX,7,7,1,1); CONPOR(1.0,HIGH,1,NX,7,7,1,1)
CONPOR(0.185,EAST, 1,NX,8,8,1,1); CONPOR(0.185,CELL,1,NX,8,8,1,1)
CONPOR(0.0  ,NORTH,1,NX,8,8,1,1); CONPOR(0.877,HIGH,1,NX,8,8,1,1)
CONPOR(0.0  ,EAST ,1,NX,9,10,1,1);CONPOR(0.0,CELL,1,NX,9,10,1,1)
CONPOR(0.0  ,NORTH,1,NX,9,10,1,1);CONPOR(0.0,HIGH,1,NX,9,10,1,1)
CONPOR(0.997,EAST ,1,NX,8,8,2,2); CONPOR(0.997,CELL,1,NX,8,8,2,2)
CONPOR(0.946,NORTH,1,NX,8,8,2,2); CONPOR(1.0,HIGH,1,NX,8,8,2,2)
CONPOR(0.706,EAST ,1,NX,9,9,2,2); CONPOR(.706,CELL,1,NX,9,9,2,2)
CONPOR(0.473,NORTH,1,NX,9,9,2,2); CONPOR(1.0,HIGH,1,NX,9,9,2,2)
CONPOR(0.236,EAST ,1,NX,10,10,2,2);CONPOR(0.236,CELL,1,NX,10,10,2,2)
CONPOR(0.0  ,NORTH,1,NX,10,10,2,2);CONPOR(1.0,HIGH,1,NX,10,10,2,2)
    GROUP 13. Boundary conditions and special sources
PATCH(CHSO,PHASEM,1,NX,1,NY,1,NZ,1,1) ! reaction-rate patch
  COVAL(CHSO,FUEL,EBU,EBU) ! uses eddy-break-up model (EBU=GRND9)
COVAL(CHSO,FUEL,EBU,EBU) 
CHSOA=FSTOI; CHSOB=CEBU  ! chsoe =0; so IEBU in gxchso = 0
 
   ** Fuel/Oxidant inlet
INLET(FOIN,SOUTH,1,NX,NYIN,NYIN,NZIN,NZIN,1,1)
VALUE(FOIN,P1,RHOIN*VINF);VALUE(FOIN,V1,VINF)
TKEIN=2.5E-5*VINF*VINF;EPSIN=0.09*TKEIN**2/(9.59*ENUL)
VALUE(FOIN,KE,TKEIN);VALUE(FOIN,EP,EPSIN);VALUE(FOIN,H1,HINF)
VALUE(FOIN,MIXF,FINF);VALUE(FOIN,FUEL,FINF) 
   ** Secondary oxidant inlet 1
INLET(SOIN1,CELL,1,1,NY,NY,NZJ1,NZJ1,1,1)
VALUE(SOIN1,P1,SECIN);VALUE(SOIN1,V1,VIN)
TKEIN=7.6875E-5*VIN*VIN;EPSIN=0.09*TKEIN**2/(16.18*ENUL)
VALUE(SOIN1,KE,TKEIN);VALUE(SOIN1,EP,EPSIN);VALUE(SOIN1,H1,HAIRIN)
 
   ** Secondary oxidant inlet 2
INLET(SOIN2,CELL,1+NX/2,1+NX/2,NY,NY,NZJ1,NZJ1,1,1)
VALUE(SOIN2,P1,SECIN);VALUE(SOIN2,V1,VIN)
VALUE(SOIN2,KE,TKEIN);VALUE(SOIN2,EP,EPSIN);VALUE(SOIN2,H1,HAIRIN)
 
   ** Dilution inlet
INLET(DILUIN,CELL,1,1,NY,NY,NZJ2,NZJ2,1,1)
VALUE(DILUIN,P1,DILIN);VALUE(DILUIN,V1,VIN)
EPSIN=0.09*TKEIN**2/(21.65*ENUL);VALUE(DILUIN,KE,TKEIN)
VALUE(DILUIN,EP,EPSIN);VALUE(DILUIN,H1,HAIRIN)
   ** Outlet Boundary
OUTLET(OUTLET,HIGH,1,NX,1,NY,NZ,NZ,1,1)
   ** North walls
WALL(WALL1,NORTH,1,NX,NY,NY,3,NZJ1-1,1,1)
WALL(WALL2,NORTH,1,NX,NY,NY,NZJ1+1,NZJ2-1,1,1)
WALL(WALL3,NORTH,1,NX,NY,NY,NZJ2+1,NZ,1,1)
   ** South wall
WALL(WALL4,SOUTH,1,NX,NYB+1,NYB+1,1,NZB,1,1)
   ** Low walls
WALL(WALL5,LOW,1,NX,NYB+1,6,1,1,1,1)
WALL(WALL6,LOW,1,NX,1,NYB,NZB+1,NZB+1,1,1)
IF(:RADE:.EQ.Y) THEN
 PATCH(WALLR1,NORTH,1,NX,NY,NY,3,NZJ1-1,1,1)
 PATCH(WALLR2,NORTH,1,NX,NY,NY,NZJ1+1,NZJ2-1,1,1)
 PATCH(WALLR3,NORTH,1,NX,NY,NY,NZJ2+1,NZ,1,1)
 PATCH(WALLR4,SOUTH,1,NX,NYB+1,NYB+1,1,NZB,1,1)
 PATCH(WALLR5,LOW,1,NX,NYB+1,6,1,1,1,1)
 PATCH(WALLR6,LOW,1,NX,1,NYB,NZB+1,NZB+1,1,1)
 COVAL(WALLR1,RADY,ED2ME,RADW);COVAL(WALLR2,RADY,ED2ME,RADW)
 COVAL(WALLR3,RADY,ED2ME,RADW);COVAL(WALLR4,RADY,ED2ME,RADW)
 COVAL(WALLR5,RADZ,ED2ME,RADW);COVAL(WALLR6,RADZ,ED2ME,RADW)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=100;DENPCO=T
    GROUP 16. Termination of iterations
LITER(P1)=20;LITER(U1)=20;LITER(V1)=20;LITER(W1)=20
    GROUP 17. Under-relaxation devices
REAL(DTFV);DTFV=20.*YVLAST/VINF/NY
RELAX(U1,FALSDT,DTFV);RELAX(V1,FALSDT,DTFV);RELAX(W1,FALSDT,DTFV)
KELIN=3;RELAX(KE,LINRLX,0.5);RELAX(H1,FALSDT,10.0)
SARAH=1.0
    GROUP 18. Limits on variables or increments to them
VARMIN(TMP1)=TWALL
    GROUP 19. Data communicated by satellite to GROUND
CHSOA=FSTOI;CHSOB=CEBU;TMP2A=FSTOI;TMP2B=HFU
    GROUP 21. Print-out of variables
NPRINT=LSWEEP
    GROUP 22. Spot-value print-out
IXMON=3;IYMON=5;IZMON=5;TSTSWP=-5
    GROUP 23. Field print-out and plot control
NPLT=5;IPLTF=2;IPLTL=LSWEEP;ITABL=3
YZPR=T;IXPRF=NX/2;IXPRL=IXPRF;NYPRIN=NY/5;NZPRIN=2
NROWCO=40;ORSIZ=0.4                   ! line-printer plots
PATCH(IZ4,CONTUR,1,NX,1,NY,4,4,1,1)
PLOT(IZ4,TMP1,1.,10.)
PATCH(XSECIN1,CONTUR,1,1,1,NY,1,NZ,1,1)
PLOT(XSECIN1,MIXF,0.,10.);PLOT(XSECIN1,W1,0.,10.)
PLOT(XSECIN1,TMP1 ,0.,10.);PLOT(XSECIN1,FUEL,0.,10.)
PATCH(XSECIN2,CONTUR,1+NX/2,1+NX/2,1,NY,1,NZ,1,1)
PLOT(XSECIN2,MIXF,0.,10.);PLOT(XSECIN2,W1,0.,10.)
PLOT(XSECIN2,TMP1 ,0.,10.);PLOT(XSECIN2,FUEL,0.,10.)
LSWEEP=100
  
      The In-Form alternative
  inform9begin    
IF(:INFORM:.EQ.Y) THEN
      -----------------------
  1. statement defining temperature
(property TMP1 is MAX(:TWALL:,(H1-:HFU:*FUEL) $
/(:CPFU:*FUEL + :CPAIR:*oxid + :CPPR:*PROD)) with imat<100)
FIINIT(TMP1)= TAIR                   ! initial field settings         
  
  2. statement defining the density
CHAR(MOLWT)
MOLWT=1.0/(FUEL/:WFU: + oxid/:WAIR: + PROD/:WPR:)  
(property RHO1 is (P1+:PRESS0:)*:MOLWT:/(TMP1*8313.4))  
   
  3. statement defining the product concentration
(stored var PROD is (MIXF-FUEL)*(1+:STOIC:))  
   
  4. statement defining the oxidant (air) concentration
(stored var OXID is 1 - FUEL - PROD)  
  
  5. de-activate the settings to GRNDx for (minor) economy
  RHO1=RHOIN
  CP1=CPAIR
  TMP1=TAIR
VARMIN(FUEL)= 0.0  ! prevent temporary 
VARMIN(PROD)= 0.0  ! occurrence of negative    
VARMIN(OXID)= 0.0  ! values
ENDIF
  inform9end    
  
  inform13begin    
IF(:INFORM:.EQ.Y) THEN
  6. replace eddy-break-up patch 
char(aaa,bbb); integer(iaa,ibb)
gtparg(chso,n,aaa,iaa,bbb,ibb)
if(:aaa:.eq.y) then
+ CHSO=SKIP
endif
PATCH(WHOLE,PHASEM,1,NX,1,NY,1,NZ,1,1)
  The following statement are equivalent alternatives
  (source of FUEL at WHOLE is -CEBU*EPKE*FUEL with line)  
(source of FUEL at WHOLE is coval(CEBU*EPKE,0.0))  
  
  7. print out the reaction-rate distribution
(stored var RRAT is CEBU*EPKE*FUEL)         


  8. use longname feature for print-out clarification                               
(longname of TMP1 print as absolute_temperature_of_the_gas_Kelvin)
  
(stored var DEGF is (9./5.)*(TMP1-273))
(longname of DEGF print as degrees_Fahrenheit)
    
(longname of RRAT print as rate_of_consumption_of_fuel_kg/m^3sec)
ELSE
+ char(aaa,bbb); integer(iaa,ibb)
+ gtparg(chso,n,aaa,iaa,bbb,ibb)
+ if(:aaa:.ne.y) then
+   PATCH(CHSO,PHASEM,1,NX,1,NY,1,NZ,1,1) ! reaction-rate patch
+   COVAL(CHSO,FUEL,EBU,EBU)
+ endif
ENDIF
  inform13end