TALK=T;RUN( 1, 1)
 ** LOAD(x302) from the x Input Library
TEXT(Owen Furnace CH4 Diffusion Flame
TITLE
  DISPLAY
  The case considered is the zero-swirl gas-fired axisymmetric
  combustion chamber used in the experiments of Owen et al (
  16th Int. Symp. on Combustion, p105-117, 1976). The geometry
  comprises a 0.122m diameter tube into which natural gas (96%CH4)
  is injected at 288K axially to mix and burn with a coaxial
  annular air stream at 750K. The inlet air/fuel equivalence ratio
  is 0.9, and the combustor operating pressure is 3.8 atm.
  The length of the solution domain is taken as 4 furnace diameters.
 
  Calculations are made with one of 4 different combustion models:
  (a) fast-chemistry model with a 1-step global irreversible
  methane reaction to form CO2 and H2O; (b) as (a) but with an
  assumed double-delta pdf to account for the influence of
  concentration fluctuations; (c) eddy-break-up finite-rate
  chemistry model with a 1-step global irreversible reaction; and
  (d) as (c) but with a 2-step global reaction mechanism involving
  the intermediate CO.
 
  Radiative heat transfer is modelled via the P1 spherical-harmonics
  approximation for a gray medium.
  ENDDIS
    GROUP 1. Run title and other preliminaries
CHAR(CMOD);BOOLEAN(THRAD)
REAL(MDOTF,MDOTA,RJETF,ROUTA,TINA,TINF,WINF,WINA,KEINT,EPINT)
REAL(KFUEL,EFUEL,KAIR,EAIR,DENF,DENA,PI,AREAF,AREAA,DIAM)
REAL(YFUIN,YOXIN,YN2IN,WFUEL,WAIR,UGASC,RTUBE,PATM)
MDOTF=7.2E-3;MDOTA=137.E-3;YFUIN=1.0;YOXIN=0.232;YN2IN=1.-YOXIN
RJETF=3.15E-2; ROUTA=4.75E-2; RTUBE=6.1E-2;DIAM=2.*RTUBE
TINA=750.;TINF=288.;PATM=1.01325E5;PI=3.1415926;WFUEL=16.04
WAIR=28.84;UGASC=8413.43
 
PRESS0=3.8*PATM
 
AREAF=PI*RJETF*RJETF;AREAA=PI*(ROUTA*ROUTA-RJETF*RJETF)
MDOTF=7.2E-3;MDOTA=137.E-3
DENF=PRESS0*WFUEL/(UGASC*TINF);DENA=PRESS0*WAIR/(UGASC*TINA)
WINF=MDOTF/(DENF*AREAF);WINA=MDOTA/(DENA*AREAA)
KFUEL=(0.2*WINF)**2;KAIR=(0.2*WINA)**2
EFUEL=0.1643*KFUEL**1.5/(0.1*RJETF)
EAIR=0.1643*KAIR**1.5/(0.1*(ROUTA-RJETF))
MESG( Thermal radiation ? (default=N)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ THRAD=T
+ MESG( Thermal Radiation activated
ELSE
+ THRAD=F
ENDIF
    GROUP 2. Transience; time-step specification
STEADY=T
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1
    GROUP 4. Y-direction grid specification
NREGY=3
IREGY=1;GRDPWR(Y,10,RJETF,1.0)
IREGY=2;GRDPWR(Y,10,ROUTA-RJETF,1.0)
IREGY=3;GRDPWR(Y,10,RTUBE-ROUTA,1.0)
    GROUP 5. Z-direction grid specification
NREGZ=1
IREGZ=1;GRDPWR(Z,40,0.5,1.4)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);STORE(VIST,DEN1,TMP1)
SOLUTN(P1,P,P,Y,P,P,P);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(W1,P,P,P,P,P,N);TURMOD(KEMODL)
    GROUP 8. Terms (in differential equations) & devices
ENUL=4.2E-5
    GROUP 9. Properties of the medium (or media)
  *** START OF EXTENDED SCRS MODEL SETTINGS
   CFUEL=CH4 COXID=O2 CFP1=CO2 CFP2=H2O NRSYS=-2
INTEGER(NSPEC,NELEM);NSPEC=7;NELEM=4
INTEGER(ICOMB,NCSTEP,NCREAC)
MESG( Enter required combustion model
MESG( Default: 1 step finite rate EBU
MESG( The options are:
MESG(  EBU1  - 1 step finite-rate EBU
MESG(  EBU2  - 2 step finite-rate EBU
MESG(  BURN  - Infinite-rate model
MESG(  DDEL  - Infinite-rate model with Double-Delta PDF
READVDU(CMOD,CHAR,EBU1)
CASE :CMOD: OF
WHEN EBU1,4
+ MESG(1 step finite-rate EBU model
+ NCSTEP=1;NCREAC=1;ICOMB=0
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE*)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
+ SCRS(EBU,P,6.0)
+ STORE(P1RS)
WHEN EBU2,4
+ MESG(2 step finite-rate EBU model
+ NCSTEP=2;NCREAC=2;ICOMB=1
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FRATE*)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,STWO)
+ SCRS(EBU,P,6.0);SCRS(EBU,S1,1.0)
+ STORE(P1RS,S1RS)
WHEN BURN,4
+ MESG(Infinite rate model
+ NCSTEP=1;NCREAC=1;ICOMB=2
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FASTC)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
WHEN DDEL,4
+ MESG(Infinite-rate model with Double-Delta PDF
+ NCSTEP=1;NCREAC=1;ICOMB=3
+ SCRS(SYSTEM,NCSTEP,NCREAC,NELEM,FASTC)
+ SCRS(SPECIES,CH4,O2,H2,CO,H2O,CO2,N2)
+ SCRS(PROP,CHEMKIN,SCH4)
+ SCRS(PDF,DDELTA)
ENDCASE
STORE(MMWT,MCO,MCO2)
   ** Define fuel & oxidiser composition & temperature
SCRS(FUIN,YFUIN,0.0,0.0,0.0,0.0,0.0,0.0,TINF)
SCRS(OXIN,0.0,YOXIN,0.0,0.0,0.0,0.0,YN2IN,TINA)
  *** END OF EXTENDED SCRS MODEL SETTINGS
IF(THRAD) THEN
+ REAL(ABSORB,SCAT,SIGMA,EMPW,EMISW,EMISG,EMPG,TWAL,TGUES)
+ ABSORB=0.25/DIAM;SCAT=0.; EMISG=0.07;TWAL=600.;TGUES=600.
+ SIGMA=5.6697E-8; EMISW=0.9
+ EMPW=SIGMA*TWAL**4; EMPG=SIGMA*EMISG*TGUES**4
+ RADIAT(RADI,ABSORB,SCAT,H1)
+ SOLUTN(SRAD,P,P,Y,P,P,P);SOLUTN(H1,P,P,N,P,P,P)
IF(ICOMB.GT.1) THEN
+ SOLUTN(SRAD,P,P,N,P,P,P);
ENDIF
ELSE
+ STORE(H1,SRAD)
ENDIF
    GROUP 11. Initialization of variable or porosity fields
INIADD=F
KEINT=(0.15*WINA)**2; EPINT=0.1643*KEINT**1.5/(0.1*RTUBE)
FIINIT(W1)=WINA; FIINIT(KE)=KEINT; FIINIT(EP)=EPINT
FIINIT(H1)=4.724E5
FIINIT(TMP1)=TINA;FIINIT(F)=0.0;FIINIT(DEN1)=DENA
IF(THRAD) THEN
+ FIINIT(SRAD)=EMPG; FIINIT(H1)=4.72E5
ENDIF
    GROUP 13. Boundary conditions and special sources
   * FUEL-STREAM INLET boundary condition
INLET(SCRSF,LOW,1,NX,#1,#1,1,1,1,1)
VALUE(SCRSF,P1,GRND1); VALUE(SCRSF,W1,WINF)
VALUE(SCRSF,EP,EFUEL); VALUE(SCRSF,KE,KFUEL)
VALUE(SCRSF,F,1.); VALUE(SCRSF,CH4,YFUIN)
IF(THRAD) THEN
+ VALUE(SCRSF,H1,GRND3)
ENDIF
   * AIR-STREAM INLET boundary condition
INLET(SCRSO,LOW,1,NX,#2,#2,1,1,1,1)
VALUE(SCRSO,P1,GRND1); VALUE(SCRSO,W1,WINA)
VALUE(SCRSO,EP,EAIR); VALUE(SCRSO,KE,KAIR)
VALUE(SCRSO,F,0.)
IF(THRAD) THEN
+ VALUE(SCRSO,H1,GRND3)
ENDIF
 
PATCH(OUT,HIGH,1,NX,1,NY,#NREGZ,#NREGZ,1,1)
COVAL(OUT,P1,1.E2,0.);COVAL(OUT,F,ONLYMS,SAME)
IF(THRAD) THEN
+ VALUE(OUT,H1,SAME)
ENDIF
   * WALL boundary condition
PATCH(NWALL3,NWALL,1,NX,#NREGY,#NREGY,1,NZ,1,1)
COVAL(NWALL3,W1,GRND2,0.0);COVAL(NWALL3,KE,GRND2,GRND2)
COVAL(NWALL3,EP,GRND2,GRND2)
 
IF(THRAD) THEN
+ PATCH(NWALL3R,NORTH,1,NX,#NREGY,#NREGY,1,NZ,#1,1)
+ COVAL(NWALL3R,SRAD,2.*EMISW/(2.0-EMISW),EMPW)
+ PATCH(LWALL3R,LOW,1,NX,#3,#3,1,1,#1,1)
+ COVAL(LWALL3R,SRAD,2.*EMISW/(2.0-EMISW),EMPW)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=600
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=ZWLAST/NZ/WINA
RELAX(P1,LINRLX,0.5); RELAX(V1,FALSDT,DTF)
RELAX(W1,FALSDT,DTF); RELAX(KE,LINRLX,0.5)
RELAX(EP,LINRLX,0.5); KELIN = 3
RELAX(F,LINRLX,0.5) ; RELAX(DEN1,LINRLX,0.01)
 
IF(ICOMB.EQ.0) THEN
    * 1 step finite-rate EBU model (ICOMB=0)
+ RELAX(CH4,FALSDT,0.5*DTF)
ENDIF
IF(ICOMB.EQ.1) THEN
    * 2 step finite-rate EBU model (ICOMB=1)
+ RELAX(CH4,FALSDT,0.5*DTF);RELAX(CO,FALSDT,0.5*DTF)
ENDIF
 
IF(THRAD) THEN
+ RELAX(H1,FALSDT,20.*DTF); RELAX(SRAD,FALSDT,20.*DTF)
IF(ICOMB.EQ.2) THEN
  *** Infinite-rate model with radiation
+ LSWEEP=1000
+ RELAX(v1,FALSDT,5.*DTF);RELAX(W1,FALSDT,5.*DTF)
+ RELAX(H1,FALSDT,DTF); RELAX(SRAD,FALSDT,DTF)
+ RELAX(DEN1,LINRLX,1.E-4)
ENDIF
IF(ICOMB.EQ.3) THEN
  *** Infinite-rate model with Double-Delta PDF & radiation
+ LSWEEP=800
+ RELAX(H1,FALSDT,2.*DTF); RELAX(SRAD,FALSDT,2.*DTF)
ENDIF
ENDIF
    GROUP 18. Limits on variables or increments to them
VARMIN(DEN1)=1.E-10
VARMAX(F)=1.0;VARMIN(F)=1.0E-10
IF(ICOMB.EQ.3) THEN
  * Infinite-rate model with Double-Delta PDF
+ RELAX(FSQ,LINRLX,0.5)
+ VARMIN(FSQ)=1.E-10
ENDIF
    GROUP 20. Preliminary print-out
ECHO=T
    GROUP 21. Print-out of variables
TSTSWP=-1
OUTPUT(TMP1,P,P,P,P,Y,Y); OUTPUT(DEN1,P,P,P,P,Y,Y)
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=2;IZMON=6
    GROUP 23. Field print-out and plot control
 
    GROUP 24. Dumps for restarts
   RESTRT(P1,V1,W1,KE,EP,H1,TMP1,DEN1,VIST,F)
   RESTRT(CH4,O2,H2,CO,H2O,CO2,N2,SPH1)
 LIBREF  =     302
  
  print-out of mass average values elicited by spedat
SPEDAT(PRINT,NUMBER,I,10)
SPEDAT(PRINT,COMMAND1,C,AVE_DEN1)
SPEDAT(PRINT,COMMAND2,C,AVE_P1)
SPEDAT(PRINT,COMMAND3,C,AVE_TMP1)
SPEDAT(PRINT,COMMAND4,C,AVE_F)
SPEDAT(PRINT,COMMAND5,C,AVE_V1)
SPEDAT(PRINT,COMMAND6,C,AVE_W1)
SPEDAT(PRINT,COMMAND7,C,AVE_EP)
SPEDAT(PRINT,COMMAND8,C,AVE_KE)
SPEDAT(PRINT,COMMAND9,C,AVE_VIST)
SPEDAT(PRINT,COMMAND10,C,AVE_H1)
STOP