c************************************************************************
c A message to prospective users of the GENXIX program
c
c From: Brian Spalding (dbs@cham.co.uk)
c Date: 28.07.99
c Text: This program is being made available for down-loading from the
c CHAM web-site in response to requests from prospective users.
c
c If there is sufficient interest, I can make relevant parts of
c the text of the original book available also. Please let me
c know if you would like this.
c
c************************************************************************
c gnmx.for 06.10.97
c
C This Fortran file was recreated during 1997 from printed records
c of the 1977 coding. It appears to work correctly.
c
c How to use it is explained in the book entitled:
c
c 'GENMIX, a general computer program for two-dimensional
c parabolic phenomena', by D. B. Spalding
c
c published by Pergamon Press, London, 1977
c
c************************************************************************
c
BLOCK DATA
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C 'GENMIX, A GENERAL COMPUTER PROGRAM FOR TWO-DIMENSIONAL
C PARABOLIC PHENOMENA', BY D. B. SPALDING
C REPORT NO. HTS/77/9, FEBRUARY 1977,
C IMPERIAL COLLEGE, MECHANICAL ENGINEERING DEPARTMENT, LONDON,SW72BX
C
C APPENDIX A (BASIC PROGRAM) - COMBUSTION OF METHANE AND AIR IN
C A DIVERGENT DUCT EXHAUSTING INTO THE ATMOSPHERE.
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1 PRELIMINARIES 1 1 1 1 1 1 1 1
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
COMMON/COMB/
1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
common/dbg/debug
logical debug
C
C ------------------------------------------------------------------
save
DATA KUDIF/-1/
DATA KASE,IRUN/0,0/
DATA ITEST/1/
DATA BIG,TINY/1.E20,1.E-10/
ccc data big,tiny/1.e5,1.e-5/
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 GRID AND GEOMETRY 2 2 2 2 2 2 2 2
C ------------------------------ GRID
DATA N,OMPOW/20,1./
c ----- set idimf= dimension for i, ensure that idimf.ge.n
DATA IDIMF/20/
c ----- set krad=1 for plane, 2 for axial, 3 for point symmetry
DATA KRAD,CSALFA/2,1./
C ------------------------------ GEOMETRY
DATA HEX0,XHEX0,AHEX,BHEX,CHEX/.05,0.,.01,0.,0./
DATA HIN0,XHIN0,AHIN,BHIN,CHIN/.02,0.,0.,0.,0./
DATA HDIV/.025/
DATA XU,XEND,XOUT/0.,.5,1./
DATA LASTEP,XULAST/100,2./
C ----- FOR KIND=3, XU, XHEX0 AND XHIN0=.25 (MAIN CH.2)
C
C-----------------------------------------------------------------------
chapter 3 3 3 3 3 3 3 dependent variables 3 3 3 3 3 3 3
c ----- set nf= number of dependent variables, excluding velocity
DATA NF/3/,JH,JP,JF,JOX,JTE,JPR/1,2,3,4,5,6/
c -----set novel=1 for no velocity, novel=2 otherwise
DATA NOVEL/2/
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4 4 PROPERTY DATA 4 4 4 4 4 4 4 4
C -------------------- S.I. UNITS
DATA AGRAV,GASCON/9.81,8314./
c ----- set model=1 for laminar flow,
c ----- set model=2 for 'mixing-length' model of turbulence
DATA MODEL/2/
ccc data model/1/
DATA AK,FR,CEBU,EWALL/.435,.033,.4,9./
c ----- initial almg, and adjusted value for ch.7
DATA ALMG/.09/, ALMGD/.075,.1,.14,.075/
c --- set inert=1 for inert fluid, inert=2 for chemically reactive
DATA INERT/2/
c ------------------------------ materials
c ---------- thermodynamics
DATA CFU,COX,CPR,CMIX/4*1100./
ccc data cfu,cox,cpr,cmix/4*1./
DATA WFU,WOX,WPR,WMIX/16.,32.,28.,29./
DATA HFU/4.E7/
ccc data hfu/0.0/
C ---------- CHEMICAL
DATA STOICH,ARRCON,PREEXP/4.,18.E3,1./
C ---------- TRANSPORT
DATA VISFU,VISOX,VISPR,VISMIX/4*1.E-6/
DATA PRLAM,PRTURB/.7,.86/
DATA H,UFAC/.9,.01/
C
C-----------------------------------------------------------------------
CHAPTER 5 5 5 5 5 5 5 STARTING VALUES 5 5 5 5 5 5 5 5
DATA PRESS/1.E5/
C ----- STREAM B IS PURE FUEL
DATA UB,TB,FUB,OXB/100.,1000.,1.,0./
ccc data ub,tb,fub,oxb/100.,100.,1.,0./
C ----- STREAM C IS AIR
DATA UC,TC,FUC,OXC/50.,1000.,0.,.232/
ccc data uc,tc,fuc,oxc/100.,100.,0.,.232/
c ----- set kex and kin for initial boundary type,
c ----- 1 for wall, 2 for free boundary, 3 for symmetry axis
DATA KEX,KIN/1,1/
C
C-----------------------------------------------------------------------
CHAPTER 6 6 6 6 6 6 6 6 STEP CONTROL 6 6 6 6 6 6 6 6 6
DATA FRA,DXMAX,DXRAT/1.,1.,5./DXPSI,DXRE,DXY/3*0.0/
C -----ENTRAINMENT CONTROL
DATA ULIM,PEILIM,FACEXP/.02,.05,.5/
C -----STARTING VALUES
DATA FACE,FACI,RATE,RATI/4*1./
C
C-----------------------------------------------------------------------
CHAPTER 7 7 7 7 7 7 7 7 BOUNDARY CONDITIONS 7 7 7 7 7 7
C ---------- STREAM A, THROUGH CENTRAL PIPE
DATA UA,TA,FUA,OXA/100.,2000.,0.,0./
ccc data ua,ta,fua,oxa/100.,100.,0.,0./
C ---------- STREAM D, SURROUNDING ATMOSPHERE
DATA TD,FUD,OXD/300.,0.,.232/
ccc data td,fud,oxd/100.,0.,.232/
C ----- UD IS SUPPLIED BY WAY OF THE UEX FUNCTION
C ---------- VELOCITY ALONG OUTER BOUNDARY
DATA UEX0,XUEX0,AUEX,BUEX,CUEX/5*0./
C ---------- WALL TEMPERATURE OF OUTER TUBE
DATA TWALL/299./
ccc data twall/100./
C
C-----------------------------------------------------------------------
CHAPTER 11 11 11 11 11 11 11 11 11 11 PRINT 11 11 11 11 11 11
c --- set ilplot=2 for down-stream plot, =1 for no plot
c --- set itplot=2 for cross-stream plot, =1 for no plot
DATA ILPLOT,ITPLOT/2,2/
c --- set nstat, nprof, nplot to number of steps between output of
c --- station values, profiles and cross-stream plots respectively
DATA NSTAT,NPROF,NPLOT/12,12,10000/
ccc data nstat,nprof,nplot/1,1,10000/
c ----- after xu=xout, nstat and nprof are set =24 at main, ch.11
C
END
c-----------------------------------------------------------------------
PROGRAM MAIN
common/dbg/debug
logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C-----------------------------------------------------------------------
CHAPTER 1 1 1 1 1 1 1 1 PRELIMINARIES 1 1 1 1 1 1 1 1
C-----------------------------------------------------------------------
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
COMMON/COMB/
1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
save
C
C ------------------------------------------------------------------
C ----- FUNCTIONS FOR BOUNDARY CONDITIONS
HEX(X)=HEX0+X*(AHEX+X*(BHEX+X*CHEX))
HIN(X)=HIN0+X*(AHIN+X*(BHIN+X*CHIN))
UEX(X)=UEX0+X*(AUEX+X*(BUEX+X*CUEX))
debug=.false.
C
C-----------------------------------------------------------------------
CHAPTER 2 2 2 2 2 2 2 GRID AND GEOMETRY 2 2 2 2 2 2 2 2
C SEE DATA
c ----- kind is an index which denotes a particular geometry type
if(debug) write(6,*)'>>> main'
KIND=4
IF(KRAD.EQ.1) KIND=2
IF(KRAD.EQ.2 .AND. nint(CSALFA).EQ.1) KIND=1
IF(KRAD.EQ.2 .AND. nint(CSALFA).EQ.0) KIND=3
C ----- MODIFICATIONS TO DATA
IF(KIND.NE.3) GO TO 21
XU=.25
XHEX0=.25
XHIN0=.25
21 CONTINUE
C
SNALFA=SQRT(1.-CSALFA**2)
C ----- STARTING VALUES
c!!!! change needed to avoid values too large for conversion
cccc IEND=IFIX(XEND*1.E6)
cccc IOUT=IFIX(XOUT*1.E6)
IEND=IFIX(XEND*1.E4)
IOUT=IFIX(XOUT*1.E4)
C -------------------- SUBROUTINE COMP UTE, ENTRY INIT
CALL INIT
C ------------------------------ GRID
DO 20 I=1,N
20 OM(I)=(FLOAT(I-1)/FLOAT(NM1))**OMPOW
C
C -------------------- SUBROUTINE COMP UTE, ENTRY GRID
CALL GRID
C
C-----------------------------------------------------------------------
CHAPTER 3 3 3 3 3 3 3 DEPENDENT VARIABLES 3 3 3 3 3 3 3
C SEE DATA
C U(I)= VELOCITY
C F(I,JH)= STAGNATION ENTHALPY
C F(I,JP)= PHI= OXIDANT CONCENTRATION - F(I,JF)*STOICH
C F(I,JF)= FUEL CONCENTRATION
C F(I,JOX)= OXIDANT CONCENTRATION
C F(I,JTE)= TEMPERATURE
C F(I,JPR)= PRODUCT CONCENTRATION
C
C-----------------------------------------------------------------------
CHAPTER 4 4 4 4 4 4 4 4 PROPERTY DATA 4 4 4 4 4 4 4 4
C SEE DATA
RECWFU=1./WFU
RECWOX=1./WOX
RECWPR=1./WPR
RECWMX=1./WMIX
GAMMA=CMIX/(CMIX-GASCON*RECWMX)
gamma=1.4
DO 40 J=1,NF
PRL(J)=PRLAM
RECPRL(J)=1./PRLAM
40 RECPRT(J)=1./PRTURB
C
C-----------------------------------------------------------------------
CHAPTER 5 5 5 5 5 5 5 STARTING VALUES 5 5 5 5 5 5 5 5
C SEE DATA
WB=1./(FUB*RECWFU+OXB*RECWOX+(1.-FUB-OXB)*RECWPR)
RHOB=PRESS*WB/(TB*GASCON)
WC=1./(FUC*RECWFU+OXC*RECWOX+(1.-FUC-OXC)*RECWPR)
RHOC=PRESS*WC/(TC*GASCON)
FLOB=RHOB*UB*(HDIV-HIN0)
FLOC=RHOC*UC*(HEX0-HDIV)
IF(KRAD.EQ.1) GO TO 55
XSIN=XU*SNALFA
HCOS=.5*CSALFA
FLOB=FLOB*(XSIN+HCOS*(HDIV+HIN0))
FLOC=FLOC*(XSIN+HCOS*(HEX0+HDIV))
55 CONTINUE
OMDIV=FLOB/(FLOB+FLOC+TINY)
TMIN=.5*AMIN1(TA,TB,TC,TD,TWALL)
tmin=90.
C -------------------- SEQUENCE TO PUT CELL BOUNDARY AT OMDIV.
IF(OMDIV.LE.1.E-10.OR.OMDIV.GE.(1.-1.E-10)) GO TO 53
DO 52 I=3,NM1
IF(OMINT(I)-OMDIV) 52,53,57
57 IDIV=I+1
GO TO 58
52 CONTINUE
58 FAC=OMDIV/OMINT(IDIV-1)
DO 59 I=2,IDIV
59 OM(I)=OM(I)*FAC
C -------------------- SURBROUTINE COMP UTE, ENTRY GRID
CALL GRID
53 CONTINUE
C ---------- INSERTION INTO ARRAYS
ENTHB=TB*(CFU*FUB+COX*OXB+CPR*(1.-FUB-OXB))+.5*UB**2+HFU*FUB
ENTHC=TC*(CFU*FUC+COX*OXC+CPR*(1.-FUC-OXC))+.5*UC**2+HFU*FUC
enthb=tb*(cfu*fub+cox*oxb+cpr*(1.-fub-oxb))+hfu*fub
enthc=tc*(cfu*fuc+cox*oxc+cpr*(1.-fuc-oxc))+hfu*fuc
PHIB=OXB-FUB*STOICH
PHIC=OXC-FUC*STOICH
DO 501 I=1,N
IF(OM(I).GT.OMDIV) GO TO 503
U(I)=UB
F(I,JH)=ENTHB
F(I,JP)=PHIB
F(I,JF)=FUB
GO TO 501
503 U(I)=UC
F(I,JH)=ENTHC
F(I,JP)=PHIC
F(I,JF)=FUC
501 CONTINUE
C
c --------------------------------- enter main loop at chapter 7
GO TO 700
C
C-----------------------------------------------------------------------
CHAPTER 6 6 6 6 6 6 6 6 STEP CONTROL 6 6 6 6 6 6 6 6
C SEE DATA
600 DXY=FRA*Y(NM2)
DXRE=DXY*PEI/(.5*(R(1)+R(N))*EMU(1)+TINY)
DXINC=DXLAST*DXRAT
C
C -------------------- DETERMINATION OF BOUNDARY TYPE
C -------------------- I BOUNDARY
IF(ISTEP.GE.IEND) GO TO 610
KIN=1
GO TO 611
610 IF(PSII.LE.TINY) GO TO 612
KIN=2
GO TO 611
612 KIN=3
C -------------------- E BOUNDARY
611 IF(ISTEP.GE.IOUT) GO TO 613
KEX=1
GO TO 614
613 KEX=2
614 CONTINUE
C
C -------------------- ENTRAINMENT RATES
IF(KIN.NE.2.AND.KEX.NE.2) GO TO 602
KUDIF=ISTEP
UMAX=U(1)
UMIN=U(1)
DO 615 I=2,N
UMAX=AMAX1(UMAX,U(I))
615 UMIN=AMIN1(UMIN,U(I))
UDIF=UMAX-UMIN
C --------------------------- I BOUNDARY
IF(KIN.NE.2) GO TO 601
RATI=ABS((U(2)-U(1))/(UDIF*ULIM+TINY))
RMI=(R(2)+R(3))*(EMU(2)+EMU(3))*RECYDF(2)*RATI/(1.+RATI)
FACI=FACI*RATI**FACEXP
FACI=AMAX1(0.1,AMIN1(FACI,10.))
RMI=RMI*FACI
IF(MODEL.EQ.2) RMI=AMAX1(RMI,0.4*UDIF*RHO(1)*R(1))
C --------------------------- E BOUNDARY
601 IF(KEX.NE.2) GO TO 602
RATE=ABS((U(NM1)-U(N))/(UDIF*ULIM+TINY))
RME=-(R(NM2)+R(NM1))*(EMU(NM2)+EMU(NM1))*RECYDF(NM2)*RATE/
1 (1.+RATE)
FACE=FACE*RATE**FACEXP
FACE=AMAX1(.01,AMIN1(FACE,10.))
RME=RME*FACE
IF(MODEL.EQ.2) RME=AMAX1(RME,-0.4*UDIF*RHO(N)*R(N))
C
602 DXPSI=PEI*PEILIM/(RMI-RME+TINY)
C
C --------------- SET VALUE OF DX
DX=AMIN1(DXY,DXRE,DXINC,DXPSI,DXMAX)
C
IF(ISTEP.GE.IEND) GO TO 605
IF(DX.LT.(XEND-XU)) GO TO 605
c ----- reset dx so that xu will exactly equal xend at next step
DX=XEND-XU
IEND=ISTEP+1
JUSTIN=ISTEP+1
C
605 IF(ISTEP.GE.IOUT) GO TO 606
IF(DX.LT.(XOUT-XU)) GO TO 606
c ----- reset dx so that xu will exactly equal xout at next step
DX=XOUT-XU
IOUT=ISTEP+1
JUSTEX=ISTEP+1
C
606 IF(PSII.GT.RMI*DX) GO TO 607
IF(PSII.LE.TINY) GO TO 607
c ----- reset dx so that axis is reached at next step
DX=PSII/RMI
JUSTIN=ISTEP+1
C
C ----- RESET DX SO THAT XU WILL NOT EXCEED XULAST
607 DX=AMIN1(DX,XULAST-XU)
C
C ---------- TRAP ZERO OR NEGATIVE DX
IF(DX.GT.0.) GO TO 608
IFIN=2
GO TO 1100
C
C -----DETERMINE XD
608 XD=XU+DX
DXLAST=DX
C
C ----- IF CSALFA VARIES
C RECALCULATE IT, AND SNALFA AND HCOS, HERE, FOR X=XD
GO TO 70
C
C-----------------------------------------------------------------------
CHAPTER 7 7 7 7 7 7 7 7 BOUNDARY CONDITIONS 7 7 7 7 7 7
700 ASSIGN 751 TO ISTART
C ---------- GENERAL BOUNDARY CONDITION INFORMATION
C ---------- STREAM A,THROUGH CENTRAL PIPE
C SEE DATA
ENTHA=TA*(CFU*FUA+COX*OXA+CPR*(1.-FUA-OXA))+.5*UA**2+HFU*FUA
PHIA=OXA-FUA*STOICH
RHOA=PRESS*WPR/(TA*GASCON+TINY)
FLOA=RHOA*UA*HIN(XEND)
IF(KRAD.EQ.2) FLOA=FLOA*(XEND*SNALFA+HCOS*HIN(XEND))
PSII=FLOA
PEI=FLOB+FLOC
PSIE=PSII+PEI
C ---------- STREAM D, SURROUNDING ATMOSPHERE
C SEE DATA
XUEX0=XOUT
UD=UEX0
XD=XU
ENTHD=TD*(CFU*FUD+COX*OXD+CPR*(1.-FUD-OXD))+.5*UD**2+HFU*FUD
PHID=OXD-FUD*STOICH
C ---------- OTHER RELATED INFORMATION
HDUCID=HIN0
ADUCTD=HEX0-HDUCID
IF(KRAD.EQ.2) ADUCTD=ADUCTD*(XSIN+HCOS*(HEX0+HDUCID))
AFLOWD=ADUCTD
70 CONTINUE
C --------------- BOUNDARY CONDITION FOR FORWARD STEP
C ------------------------------------------- I BOUNDARY
IF(KIN-2) 731,732,733
C ------------------------------ WALL
731 IF(ISTEP.GT.JUSTIN) GO TO 734
U(1)=0.
TAUI=0.
RMI=0.
DO 735 J=1,NF
IBIN(J)=2
735 RJTOTI(J)=0.
C -------------------- ADJUST INNER HEIGHT
734 HIND=HIN(XD-XHIN0)
GO TO 740
C ------------------------------ FREE BOUNDARY
732 IF(ISTEP.GT.JUSTIN) GO TO 736
TAUI=0.
U(1)=UA
RHO(1)=RHOA
RECRU(1)=1./(RHO(1)*U(1)+TINY)
F(1,JH)=ENTHA
F(1,JP)=PHIA
F(1,JF)=FUA
AREA=HDUCID
IF(KRAD.EQ.2) AREA=AREA*(XU*SNALFA+HCOS*HDUCID)
AFLOWD=AFLOWD+AREA
c!!!! mistyping
cccc IF(ISTEP.GT.0) GO TO 740
IF(ISTEP.eq.0) GO TO 740
736 U(1)=U(1)+DX*AGRAV*(RHO(N)-RHO(1))*RECRU(1)
GO TO 740
C ------------------------------ SYMMETRY AXIS
733 IF(ISTEP.GT.JUSTIN) GO TO 740
TAUI=0.
RMI=0.
PSII=0.
HIND=0.
U(1)=U(2)
DO 737 J=1,NF
737 F(1,J)=F(2,J)
C ---------- NO SUBSEQUENT CHANGE NEEDED
740 CONTINUE
C
C ------------------------------------------- E BOUNDARY
IF(KEX-2) 741,742,743
C ------------------------------ WALL
741 IF(ISTEP.GT.JUSTEX) GO TO 744
C ---------- FIRST STEP ONLY
U(N)=0.
RME=0.
TAUE=0.
IBEX(JH)=1
F(N,JOX)=OXC
F(N,JPR)=1.-OXC-FUC
DO 745 J=2,NF
IBEX(J)=2
745 RJTOTE(J)=0.
C ----- ADJUST ENTHANPY FIT COMPOSITION
744 CMIX=CFU*F(N,JF)+COX*F(N,JOX)+CPR*F(N,JPR)
F(N,JTE)=TWALL
F(N,JH)=CMIX*F(N,JTE)+F(N,JF)*HFU
C -------------------- ADJUST EXTERNAL HEIGHT
HEXD=HEX(XD-XHEX0)
GO TO 750
C -------------------------------------------- FREE BOUNDARY
742 IF(ISTEP.GT.JUSTEX) GO TO 746
F(N,JH)=ENTHD
F(N,JP)=PHID
F(N,JF)=FUD
F(N,JOX)=OXD
F(N,JPR)=1.-F(N,JF)-F(N,JOX)
VMIX=F(N,JF)*RECWFU+F(N,JOX)*RECWOX+F(N,JPR)*RECWPR
F(N,JTE)=TD
RHO(N)=PRESS/(VMIX*F(N,JTE)*GASCON)
U(N)=UD
RECRU(N)=1./(RHO(N)*U(N)+TINY)
C ---------- ADJUSTMENT OF MIXING LENGTH CONSTANT
ALMG=ALMGD(KIND)
C ---------- ADJUSTMENT OF DOWNSTREAM VELOCITY
746 UD=UEX(XD-XUEX0)
GO TO 750
C ---------- NO SYMMETRY AXIS
743 CONTINUE
750 GO TO ISTART, (751,800)
751 ASSIGN 800 TO ISTART
GO TO 900
C
C-----------------------------------------------------------------------
CHAPTER 8 8 8 8 8 8 8 8 ADVANCE 8 8 8 8 8 8 8 8 8 8
C ------------------------ MOMENTUM SOURCES
C -------------------------------------- PRESSURE GRADIENT
800 continue
IF(KEX.NE.2) GO TO 821
DP=(U(N)-UD)/RECRU(N)
GO TO 823
821 AFLOWU=AFLOWD
HDUCID=0.
IF(KIN.EQ.1) HDUCID=HIND
ADUCTD=HEXD-HDUCID
IF(KRAD.EQ.2) ADUCTD=ADUCTD*(XD*SNALFA+HCOS*(HEXD+HDUCID))
DA=ADUCTD-AFLOWU
DP=DA/DADP
C -------------------- WALL SHEAR AND MASS ADDITION
UBAR=0.
DO 824 I=2,NM1
824 UBAR=UBAR+(BOM(I)*U(I))
IF(KIN.EQ.2) UBAR=(UBAR-U(1))*PEI/PSIE+U(1)
UBAR=(UBAR-U(1))*PEI/PSIE+U(1)
DP=DP+DX*(-TAUI*R(1)-TAUE*R(N)+2.*RME*UBAR)/ADUCTD
DP=AMIN1(DP,.5*DPMAX)
C
C ------------------------------ COMP
823 CALL SOLVE
C
C-----------------------------------------------------------------------
CHAPTER 9 9 9 9 9 9 9 9 COMPLETE 9 9 9 9 9 9 9 9 9 9
900 CONTINUE
C ---------------------- IGNITION SEQUENCE
IF(ISTEP.GT.5) GO TO 931
IF(INERT.EQ.1) GO TO 931
T2=.5/STOICH
DO 932 I=2,NM1
F(I,JF)=T2*(ABS(F(I,JP))-F(I,JP))
932 F(I,JOX)=F(I,JP)+STOICH*F(I,JF)
931 CONTINUE
C
C ------------------------------ THERMODYNAMIC PROPERTIES
PRESS=PRESS+DP
PDGSCN=PRESS/GASCON
DO 907 I=1,N
F(I,JOX)=AMAX1(0.,F(I,JP)+STOICH*F(I,JF))
F(I,JPR)=1.-F(I,JF)-F(I,JOX)
ENTH=F(I,JH)-.5*U(I)**2-HFU*F(I,JF)
ENTH=F(I,JH)-HFU*F(I,JF)
F(I,JTE)=ENTH/(CFU*F(I,JF)+COX*F(I,JOX)+CPR*F(I,JPR))
IF(F(I,JTE).GT.TMIN) GO TO 941
IF(I.EQ.1.OR.I.EQ.N) GO TO 941
WRITE(6,942) F(I,JTE),I,ISTEP,TMIN
942 FORMAT(27H *** TEMPERATURE, F(I,JTE)=,1PE10.3,6H AT I=,I4,7H ISTEP
1=,I5/17H *** RESET =TMIN=,E10.3,23H *** MAIN CH.9 COMPLETE)
F(I,JTE)=TMIN
c write(6,*)'CFU,F(I,JF)',
c 1 CFU,F(I,JF)
c write(6,*)'COX,F(I,JOX)',
c 1 COX,F(I,JOX)
c write(6,*)'CPR,F(I,JPR)',
c 1 CPR,F(I,JPR)
c write(6,*)'hfu,enth',
c 1 hfu,enth
c write(6,*)'F(I,JH),U(I)',F(I,JH),U(I)
stop
941 VMIX=F(I,JF)*RECWFU+F(I,JOX)*RECWOX+F(I,JPR)*RECWPR
907 RHO(I)=PDGSCN/(F(I,JTE)*VMIX)
IF(KEX.EQ.1) F(N,JTE)=TWALL
DPDX=DP/DX
C
C ---------------------------- RADII AND Y@S
IF(KRAD-2) 901,902,903
C ----- KRAD=1, PLANE
901 IF(KIN.EQ.2) HIND=ABS(PSII*RECRU(1))
GO TO 909
C ----- KRAD=2, AXIAL
902 IF(KIN.NE.2) GO TO 908
HIND=ABS(PSII*RECRU(1))
HIND=2.*HIND/
1 (XD*SNALFA+SQRT((XD*SNALFA)**2+2.*HIND*CSALFA)+TINY)
GO TO 908
C ----- KRAD=3, POINT SYMMETRY
903 R(I)=0.
C ----- CHANGE ABOVE STATEMENT IF NECESSARY FOR KRAD=3
GO TO 909
908 R(1)=XU*SNALFA+HIND*CSALFA
C ------------------------------- COMP
909 CALL DISTAN
C
C-----------------------------------------------------------------------
CHAPTER 10 10 10 10 10 10 10 ADJUST 10 10 10 10 10 10 10
C
IF(KEX.EQ.2) GO TO 1022
AFLOWD=Y(N)+HIND-HDUCID
IF(KRAD.EQ.2) AFLOWD=AFLOWD*(XU*SNALFA+HCOS*(Y(N)+HIND+HDUCID))
DA1=ADUCTD/AFLOWD-1.
C -------------------- DEPENDENCE OF AREA ON PRESSURE
RECGMP=1./(GAMMA*PRESS)
DADP=0.
IF(KIN.EQ.2) DADP=PSII*RECRU(1)*(RECRU(1)*RECRU(1)*RHO(1)-RECGMP)
SUM=0.
DPMAX=BIG
DO 1025 I=2,NM1
DPMAX=AMIN1(DPMAX,RHO(I)*U(I)**2)
1025 SUM=SUM+BOM(I)*RECRU(I)*(RECRU(I)*RECRU(I)*RHO(I)-RECGMP)
DADP=DADP+PEI*SUM
c -------------------- adjustment of ps,us, etc.
IF(ABS(DA1).LT.1.E-3) GO TO 1022
DP=DA1*AFLOWD/DADP
DP=AMIN1(DP,.0*DPMAX)
PRESS=PRESS+DP
DPDX=DPDX+DP/DX
RHOFAC=1.+DP*RECGMP
DO 1027 I=2,NM1
U(I)=U(I)-DP*RECRU(I)
1027 RHO(I)=RHO(I)*RHOFAC
IF(KIN.NE.2) GO TO 1029
U(1)=U(1)-DP*RECRU(1)
1029 RHO(1)=RHO(1)*RHOFAC
RHO(N)=RHO(N)*RHOFAC
RECRU(1)=1./(RHO(1)*U(1)+TINY)
IF(KIN.NE.2) GO TO 1026
HIND=ABS(PSII*RECRU(1))
IF(KRAD.EQ.1) GO TO 1026
HIND=2.*HIND/
1 (XD*SNALFA+SQRT((XD*CSALFA)**2+2.*HIND*CSALFA)+TINY)
R(I)=XU*SNALFA+HIND*CSALFA
1026 CALL DISTAN
AFLOWD=Y(N)+HIND-HDUCID
IF(KRAD.EQ.2) AFLOWD=AFLOWD*(XU*SNALFA+HCOS*(Y(N)+HIND+HDUCID))
DA2=ADUCTD/AFLOWD-1.
1022 CONTINUE
C-----------------------------------------------------------------------
CHAPTER 11 11 11 11 11 11 11 11 11 11 PRINT 11 11 11 11 11 11
C SEE DATA
1100 CONTINUE
IF(XU.LE.XOUT) GO TO 1101
NSTAT=24
NPROF=24
1101 CONTINUE
CALL OUTPUT
C
C-----------------------------------------------------------------------
CHAPTER 12 12 12 12 12 12 12 DECIDE 12 12 12 12 12 12 12
IF(ISTEP.EQ.LASTEP) GO TO 1203
IF(XU.LT.XULAST) GO TO 1202
1203 IFIN=2
CALL OUTPUT
1202 IF(IFIN.EQ.1) GO TO 600
STOP
END
SUBROUTINE WALL(I1,OUT1,OUT2)
C/FEB. 1977 --------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
common/dbg/debug
logical debug
DIMENSION S1(2),S2(2),S3(2),S4(2),S5(2),S6(2)
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
COMMON/COMB/
1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C
C EFFECTS OF PRESSURE GRADIENT AND MASS TRANSFER ARE INCLUDED
C EFFECTS OF RADIUS VARIATION ARE NEGLECTED
C FOR VELOCITY, OUT1=BP, OUT2=T
C FOR F 'S, OUT1=FIDIF, OUT2=T
C
CHAPTER A ------------------------------- PRELIMINARIES ----------------
save
DATA SHALF/.04/, BPLAST/1./
if(debug) write(6,*)'>>> wall'
KWALL=2-1/I1
I2=I1+3-2*KWALL
I3=I1+6-4*KWALL
IF(J.GT.0) GO TO 200
CHAPTER B ------------------------------- VELOCITY ---------------------
UREF=U(I2)
RHOREF=RHO(I2)
RUREF=RHOREF*UREF
RREF=R(I2)
VREF=EMU(I1)
YREF=YI+(YE-YI)*OM(I1)
RE=RUREF*YREF/VREF
RRUREF=RREF*RUREF
AM=(RMI-(RME+RMI)*OM(I1))/RRUREF
EF=YREF*DP/(DX*RUREF*UREF)
IF(MODEL.EQ.1) GO TO 110
IF(RE.LT.132.25) GO TO 110
C ------------------------------ TURBULENT FLOW
C ---------- EXTENDED LOG LAW
ER=RE*EWALL
ARGMIN=11.5*EWALL
NIT=0
101 SHALF1=SHALF
S=SHALF**2
SLOC=S+AM+EF
IF(SLOC.GT.0.) GO TO 104
SLOC=TINY
SHALF=SQRT(ABS(AM+EF))
104 BEE=SQRT(SLOC)/AK
ARG=ER*(SHALF+(AM/(1.+BEE)+.5*EF)/SHALF)
IF(ARG.LT.ARGMIN) GO TO 110
SHALF=AK/ALOG(ARG)
IF(ABS(SHALF-SHALF1).LT..0001) GO TO 102
NIT=NIT+1
IF(NIT.LT.11) GO TO 101
102 S=SHALF**2
SAV=.5*(S+SLOC)
BP=1./(1.+BEE)
GO TO 103
C ------------------------------ LAMINAR FLOW
110 AMRE=AM*RE
FRE=EF*RE
IF(ABS(AMRE).LT..01) GO TO 111
AMRE=AMAX1(-60.,AMIN1(60.,AMRE))
EXPMRE=EXP(AMRE)
STORE=EXPMRE-1.-AMRE
AMRESQ=AMRE*AMRE
SRE=AMRE*(1.-STORE*FRE/AMRESQ)/(EXPMRE-1.)
BP=SRE*STORE/AMRESQ+FRE*(STORE-.5*AMRESQ)/(AMRESQ*AMRE)
GO TO 112
111 SRE=(2.-FRE*(1.+AMRE*.33333))/(2.+AMRE)
BP=SRE*(.5+AMRE*.16667)+FRE*(.16667+AMRE*.041667)
112 IF(SRE.GT.TINY) GO TO 113
SRE=TINY
BP=.33333
113 S=SRE/RE
SAV=S
103 T=S*RRUREF
C -------------------- UNDER-RELAX BP
BP=BPLAST+.5*(BP-BPLAST)
BPLAST=BP
OUT1=BP
OUT2=T
S1(KWALL)=SAV
S2(KWALL)=RRUREF
S3(KWALL)=UREF
S4(KWALL)=AM
S5(KWALL)=AMRE
S6(KWALL)=RE
GO TO 900
C-----------------------------------------------------------------------
CHAPTER C ------------------------ OTHER DEPENDENT VARIABLES
200 SAV=S1(KWALL)
RRUREF=S2(KWALL)
UREF=S3(KWALL)
AM=S4(KWALL)
AMRE=S5(KWALL)
RE=S6(KWALL)
IF(MODEL.EQ.1) GO TO 210
C -------------------------------- TURBULENT FLOW
PRRAT=PRL(J)*RECPRT(J)
C ----- THE FOLLOWING P-FUNCTION IS NOT TO BE USED FOR PRRAT.LT.0.5
PJAY=9.*(PRRAT-1.)/PRRAT**.25
S=SAV*RECPRT(J)/(1.+AMAX1(-.99999,PJAY*SQRT(ABS(SAV))))
OUT2=S*RRUREF
IF(J.NE.JH) GO TO 221
OUT1=(H-1.)*.5*UREF**2
GO TO 900
221 OUT1=0.
GO TO 900
C ------------------------------ LAMINAR FLOW
210 IF(ABS(AMRE).LT..01) GO TO 211
S=AM/(EXP(PRL(J)*AMRE)-1.)
GO TO 212
211 S=RECPRL(J)/(RE+.5*RE*PRL(J)*AMRE)
212 OUT2=S*RRUREF
IF(J.NE.JH) GO TO 214
OUT1=(PRL(JH)-1.)*.5*UREF**2
GO TO 900
214 OUT1=0.
C ------------------------------------------------------------------
900 IF(ITEST.EQ.1) GO TO 901
WRITE(6,9000) J,I1,OUT1,OUT2
9000 FORMAT(12H WALL TESTS,,3H J=,I3,4H I1=,I3,6H OUT1=,1PE10.3,
1 6H OUT2=,E10.3)
901 RETURN
END
SUBROUTINE PLOTS(X,IDIME,IMAX,XAXIS,Y,JDIME,JMAX,YAXIS)
common/dbg/debug
logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
C *
c subroutine for plotting j curves of y(i,j) against x(i). *
c *
c x and y are scaled to the range 0. to 1., for plotting as *
c (y-ymin)/(ymax-ymin), the maximun and minimum values are printed *
c n.b. the x and y arrays must be redefined before each call plots. *
c idime is the variable dimension for x. *
c imax is the number of x values. *
c xaxis stores the name of the x-axis. *
c jdime is the variable dimension for y. *
c jmax is the number of curves to be plotted, (up to 30). *
c the array yaxis(j) stores the names of the curves, *
c the fist character of each curve-name is used for plotting. *
c xsize alters the x-plot size by a factor of .2 to 1., in steps of .1*
c ysize is the y-plot size factor fo .2 upwards in steps of .2 *
c xsize=1., ysize=1. gives normal size plot. *
c *
c***********************************************************************
DIMENSION X(IDIME),Y(IDIME,JDIME),YAXIS(JDIME),
1 A(101),YMAX(30),YMIN(30),DIGIT(11)
c!!!! equivalence deactivated. It causes YMIN to be over-written
cc EQUIVALENCE (YMAX(1),A(1)),(YMIN(1),A(31))
CHARACTER*1 DOT,CROSS,BLANK,DIGIT,A
character*8 YAXIS
save
DATA DOT,CROSS,BLANK/1H.,1H+,1H /
1,DIGIT/1h ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H1/
if(debug) write(6,*)'>>> plots'
C***** SET PLOT SIZE FACTOR
XSIZE=0.6
YSIZE=1.4
ysize=0.2
C***** SCALING X-ARRAY TO RANGE 0 TO 100*XSIZE
XR=100.*XSIZE
XMAX=-1.E30
XMIN=+1.E30
IM=IMAX
DO 1 I=1,IM
XMAX=AMAX1(XMAX,X(I))
1 XMIN=AMIN1(XMIN,X(I))
S=XR/(XMAX-XMIN+1.E-30)
DO 2 I=1,IM
2 X(I)=(X(I)-XMIN)*S
c***** scaling y-array to range 0 to 50*ysize
YR=50.*YSIZE
JM=JMAX
DO 4 J=1,JM
YMAX(J)=-1.E30
YMIN(J)=+1.E30
DO 4 I=1,IM
cccc write(6,*)'i,j,ymax(j),ymin(j)',i,j,ymax(j),ymin(j)
YMIN(J)=AMIN1(YMIN(J),Y(I,J))
4 YMAX(J)=AMAX1(YMAX(J),Y(I,J))
DO 3 J=1,JM
cccc write(6,*)'j,ymin(j)',j,ymin(j)
S=YR/(YMAX(J)-YMIN(J)+1.E-30)
DO 3 I=1,IM
3 Y(I,J)=(Y(I,J)-YMIN(J))*S
c***** write curve names, with actual min and max values
J=1
L=IFIX(XR/10.)
K=L
GO TO 6
5 J=J+L
K=K+L
6 K=MIN0(JM,K)
WRITE(6,101) (YAXIS(I),I=J,K)
WRITE(6,103) (YMIN(I),I=J,K)
WRITE(6,104) (YMAX(I),I=J,K)
IF(K-JM) 5,7,7
7 WRITE(6,106)
c***** main loop - each pass produces a y-constant line
IX=IFIX(XSIZE*10.)
KX=IFIX(XR)+1
IY=IFIX(YR*.1)
KY=IFIX(YR)+1
N=KY+1
DO 40 M=1,KY
L=N-M
IF(L.EQ.1.OR.L.EQ.KY) GO TO 32
GO TO 33
c***** put. or + along the x-axis
32 DO 30 K=1,KX
30 A(K)=DOT
DO 31 K=1,KX,IX
31 A(K)=CROSS
GO TO 45
C***** PUT. OR + AND Y-VALUES ALONG Y-AXIS
33 A(1)=DOT
A(KX)=DOT
K=L-1
46 K=K-IY
IF(K) 48,47,46
47 A(1)=CROSS
A(KX)=CROSS
45 YL=FLOAT(L-1)/YR
GO TO 35
48 YL=-1.
C***** SEARCH FOR POINTS ALONG Y-CONSTANT LINE
35 DO 42 J=1,JM
DO 42 I=1,IM
IF(IFIX(Y(I,J)+1.5).NE.L) GO TO 42
C***** POINT FOUND- ASSIGN SYMBOL- SEARCH FOR MORE
NX=X(I)+1.5
A(NX)=YAXIS(J)
42 CONTINUE
C***** POINT Y-CONSTANT LINE
IF(nint(yl).NE.-1) GO TO 37
WRITE(6,106) (A(K),K=1,KX)
GO TO 38
37 WRITE(6,107) YL,(A(K),K=1,KX)
C***** FILL ARRAY A WITH BLANKS
38 DO 49 K=1,KX
49 A(K)=BLANK
40 CONTINUE
C***** PRINT BLANK OR X-VALUE FOR X-AIXS
L=1
K=KX-1
A(1)=DIGIT(1)
DO 51 I=IX,K,IX
L=L+1
A(I)=DOT
51 A(I+1)=DIGIT(L)
A(K)=BLANK
WRITE(6,106) (A(K),K=I,KX)
WRITE(6,100) XAXIS,XMIN,XMAX
RETURN
100 FORMAT(1h ,12HABSCISSA IS ,A8,5H MIN=,1PE9.2,5H MAX=,E9.3)
101 FORMAT(1h ,9HORDINATE ,12(A8,2X))
103 FORMAT(1H ,7H MIN ,1P11E10.2)
104 FORMAT(1H ,7H MAX ,1P11E10.2)
106 FORMAT(6X,101A1)
107 FORMAT(2H ,F3.1,1X,101A1)
END
SUBROUTINE COMP
common/dbg/debug
logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
DIMENSION A(20),B(20),C(20),CON(20),D(20),HCON(20),OMDIF(20)
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F( 120),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
C
EQUIVALENCE (A(1),DIF(1)),(C(1),SI(1)),(D(1),SIP(1))
save
C
CHAPTER A --------------------------------------------------------------
C INIT ------------- INIT ------------- INIT ------------- INIT
if(debug) write(6,*)'>>> main'
ENTRY INIT
if(debug) write(6,*)'>>> init'
C ------------------------- INITIAL VALUES AND DEFAULT VALUES
NM1=N-1
NM2=N-2
NM3=N-3
ISTEP=0
IF(KRAD.EQ.3) NOVEL=1
JUSTIN=ISTEP
JUSTEX=ISTEP
IFIN=1
DXLAST=BIG
DX=BIG
PSII=0.
BPE=1.
BPI=1.
Y(1)=0.
DP=0.
DO 13 I=1,N
EMU(I)=0.
CON(I)=0.
13 R(I)=1.
IF(NOVEL.NE.1) RETURN
DO 16 I=1,N
16 U(I)=1.
RETURN
CCHPTER B --------------------------------------------------------------
C GRID ------------- GRID ------------- GRID ------------- GRID
ENTRY GRID
if(debug) write(6,*)'>>> grid'
OMI=OM(2)
OME=1.-OM(NM1)
BOM(2)=.5*(OM(2)+OM(3))
OMINT(1)=0.
OMINT(2)=BOM(2)
DO 202 I=3,NM2
OMINT(I)=.5*(OM(I)+OM(I+1))
BOM(I)=OMINT(I)-OMINT(I-1)
OMDIF(I)=OM(I)-OM(I-1)
202 HOMDFI=.5*OMDIF(3)
BOM(NM1)=1.-OMINT(NM2)
OMINT(NM1)=1.
OMDIF(NM1)=OM(NM1)-OM(NM2)
HOMDFE=.5*OMDIF(NM1)
RETURN
CHAPTER C --------------------------------------------------------------
C DISTAN ---------- DISTAN ---------- DISTAN ---------- DISTAN
ENTRY DISTAN
if(debug) write(6,*)'>>> distan'
IF(NOVEL.NE.1) GO TO 220
DO 224 I=1,N
224 RECRU(I)=1./(RHO(I)+TINY)
GO TO 222
220 DO 221 I=1,N
RECRU(I)=1./(RHO(I)*U(I)+TINY)
IF(RECRU(I).GT.0.) GO TO 221
IFIN=2
WRITE(6,223) RECRU(I),I,ISTEP
223 FORMAT(14H *** RECRU(I)=,1PE10.3,6H AT I=,I4,11H AND ISTEP=,I5,
1 16H *** COMP DISTAN)
stop
221 CONTINUE
c -------------------- calculation of y's and r's
222 IF(KIN.EQ.1) GO TO 308
RAT=RECRU(2)*RHO(1)*U(1)
IF(KRAD.EQ.2) GO TO 307
BPI=.33333+.66667*RAT
GO TO 308
307 BPI=(R(1)*(.83333*RAT+.16667)+.5*R(2)*(RAT+1.))/(R(1)+R(2))
308 IF(KEX.EQ.1) GO TO 230
RAT=RECRU(NM1)*RHO(N)*U(N)
IF(KRAD.EQ.2) GO TO 327
BPE=.33333+.66667*RAT
GO TO 230
327 BPE=(R(N)*(.83333*RAT+.16667)+.5*R(NM1)*(RAT+1.))/(R(N)+R(NM1))
c ------------------- y's for plane flow
230 STORE=OM(2)/BPI
ADPEI(2)=(HOMDFI+STORE)*RECRU(2)
Y(2)=PEI*RECRU(2)*STORE
HPEI=.5*PEI
DO 231 I=3,NM1
ADPEI(I)=BOM(I)*RECRU(I)
231 Y(I)=Y(I-1)+HPEI*OMDIF(I)*(RECRU(I)+RECRU(I-1))
STORE=OME/BPE
ADPEI(NM1)=(HOMDFE+STORE)*RECRU(NM1)
Y(N)=Y(NM1)+PEI*RECRU(NM1)*STORE
C ------------------------------------------------
IF(KRAD-2) 270,240,280
C ---------------- Y'S AND R'S FOR AXIAL SYMMETRY
240 IF(CSALFA.LT.TINY) GO TO 260
C ------------------------------ CSALFA NE 0.
COSD2=.5*CSALFA
TWDCOS=2./CSALFA
IF(R(1).GT.TINY) GO TO 250
C ------------------------------ R(1)=0.
DO 242 I=2,N
Y(I)=SQRT(ABS(Y(I)*TWDCOS))
242 R(I)=Y(I)*CSALFA
GO TO 270
C ------------------------------ R(1) NE 0.
250 R1D2=.5*R(1)
R1D2SQ=R1D2*R1D2
DO 251 I=2,N
Y(I)=Y(I)/(R1D2+SQRT(ABS(R1D2SQ+COSD2*Y(I))))
251 R(I)=R(1)+Y(I)*CSALFA
GO TO 270
C ------------------------------ CSALFA=0.
260 RECR1=1./R(1)
DO 261 I=2,N
R(I)=R(1)
261 Y(I)=Y(I)*RECR1
GO TO 270
C -------------------------------- POINT SYMMETRY, KRAD=3
280 R1CUB=R(1)**3
DO 281 I=2,N
R(I)=(R1CUB+Y(I))**.3333333
281 Y(I)=R(I)-R(1)
C ------------------------------ GENERAL
270 YI=Y(2)
YE=Y(N)-Y(NM1)
DO 273 I=1,NM1
273 RECYDF(I)=1./(Y(I+1)-Y(I))
IF(ITEST.EQ.1) RETURN
WRITE(6,274) (RHO(I),I=1,N)
WRITE(6,275) (RECRU(I),I=1,N)
WRITE(6,276) (ADPEI(I),I=1,N)
WRITE(6,277) (Y(I),I=1,N)
WRITE(6,278) (R(I),I=1,N)
WRITE(6,279) (RECYDF(I),I=1,N)
RETURN
274 FORMAT(18H0COMP DISTAN TESTS/8H RHO(I)=/(3X,1P6E11.3))
275 FORMAT(10H RECRU(I)=/(3X,1P6E11.3))
276 FORMAT(10H ADPEI(I)=/(3X,1P6E11.3))
277 FORMAT(6H Y(I)=/(3X,1P6E11.3))
278 FORMAT(6H R(I)=/(3X,1P6E11.3))
279 FORMAT(11H RECYDF(I)=/(3X,1P6E11.3))
CHAPTER D --------------------------------------------------------------
C SOLVE ---------- SOLVE ----------- SOLVE ---------- SOLVE
ENTRY SOLVE
if(debug) write(6,*)'>>> solve'
C ---------------------------------------- PRELIMINARIES
DXDPEI=DX/PEI
CONST1=.5*DXDPEI
ENT=ABS(RMI)+ABS(RME)
IF(ENT.LE.TINY) GO TO 310
HCONI=RMI*CONST1
HCONDF=(RME-RMI)*CONST1
DO 412 I=2,NM1
HCON(I)=HCONI+HCONDF*OMINT(I)
412 CON(I)=HCON(I)+HCON(I)
C ---------------------------------------- COEFFICIENTS FOR U
310 IF(NOVEL.EQ.1) GO TO 442
J=0
C ----- CALL SUBROUTINE PHYS AT ENTRY PHYSU
CALL PHYSU
IF(KRAD-2) 410,415,411
410 DO 413 I=2,NM2
413 DIFU(I)=CONST1*(EMU(I)+EMU(I+1))*RECYDF(I)
GO TO 414
415 CONST2=.5*CONST1
DO 416 I=2,NM2
416 DIFU(I)=CONST2*(R(I+1)+R(I))*(EMU(I)+EMU(I+1))*RECYDF(I)
GO TO 414
411 CONST3=.25*CONST1
DO 419 I=2,NM2
419 DIFU(I)=CONST3*(R(I+1)+R(I))**2*(EMU(I)+EMU(I+1))*RECYDF(I)
C ----------------------------------------- A'S AND B'S
414 IF(ENT.LE.TINY) GO TO 312
DO 417 I=2,NM2
A(I)=AMAX1(0.,DIFU(I)-HCON(I),-CON(I))
417 B(I+1)=A(I)+CON(I)
GO TO 314
312 DO 315 I=2,NM2
A(I)=DIFU(I)
315 B(I+1)=A(I)
314 TI=0.
TE=0.
IF(KIN.EQ.1) CALL WALL(1,BPI,TI)
IF(KEX.EQ.1) CALL WALL(N,BPE,TE)
B(2)=AMAX1((TI+RMI)*DXDPEI,0.)
A(NM1)=AMAX1((TE-RME)*DXDPEI,0.)
cc write(6,*)'TE,RME,DXDPEI',TE,RME,DXDPEI
C ------------------------------------------ C'S AND D'S
IF(MOMSOU.EQ.0) GO TO 431
DO 418 I=2,NM1
C(I)=U(I)*BOM(I)+SI(I)
418 D(I)=A(I)+B(I)+BOM(I)
GO TO 432
431 DO 433 I=2,NM1
C(I)=U(I)*BOM(I)
433 D(I)=A(I)+B(I)+BOM(I)
432 CONTINUE
IF(ITEST.EQ.1) GO TO 404
WRITE(6,341) (DIFU(I),I=2,NM1)
WRITE(6,342) (CON(I),I=2,NM1)
WRITE(6,405) (A(I),I=2,NM1)
WRITE(6,406) (B(I),I=2,NM1)
WRITE(6,407) (C(I),I=2,NM1)
WRITE(6,408) (D(I),I=2,NM1)
341 FORMAT(23HOCOMP SOLVE TESTS FOR U/9H DIFU(I)=/(3X,1P6E11.3))
342 FORMAT(8H CON(I)=/(3X,1P6E11.3))
405 FORMAT(6H A(I)=/(3X,1P6E11.3))
406 FORMAT(6H B(I)=/(3X,1P6E11.3))
407 FORMAT(6H C(I)=/(3X,1P6E11.3))
408 FORMAT(6H D(I)=/(3X,1P6E11.3))
404 CONTINUE
C ------------------------------- ADJUST FREE-BOUNDARY VALUES
IF(KIN.EQ.2) U(1)=U(1)-DP*RECRU(1)
IF(KEX.EQ.2) U(N)=U(N)-DP*RECRU(N)
C---------------------------- SOLVE FOR DOWNSREAM U'S
C(2)=(B(2)*U(1)+C(2))/D(2)
D(2)=A(2)/D(2)
DO 421 I=3,NM1
T=1./(D(I)-B(I)*D(I-1))
D(I)=A(I)*T
421 C(I)=(B(I)*C(I-1)+C(I))*T
DO 422 IDASH=1,NM2
I=N-IDASH
422 U(I)=D(I)*U(I+1)+C(I)
IF(KIN-2) 444,445,446
444 TAUI=TI*U(2)/R(1)
GO TO 445
446 U(1)=U(2)
445 IF(KEX-2) 447,440,448
447 TAUE=TE*U(NM1)/R(N)
GO TO 440
448 U(N)=U(NM1)
440 IF(ITEST.NE.1) WRITE(6,443) (U(I),I=1,N)
443 FORMAT(6H U(I)=/(3X,1P6E11.3))
C -------------------------------------------
C ---------------------------------------------- F-SECTION
442 IF(NF.LT.1) GO TO 481
DO 480 J=1,NF
IDJ=IDIMF*(J-1)
I1J=1+IDJ
I2J=2+IDJ
INM1J=NM1+IDJ
INJ=N+IDJ
C ----- CALL SUBROUTINE PHYS AT ENTRY PHYSF
CALL PHYSF
TIF=0.
FDIFI=0.
TEF=0.
FDIFE=0.
IF(KIN.EQ.1) CALL WALL(1,FDIFI,TIF)
IF(KEX.EQ.1) CALL WALL(N,FDIFE,TEF)
IF(ITEST.EQ.1) GO TO 450
WRITE(6,451) J,(DIF(I),I=2,NM1)
451 FORMAT(24H COMP SOLVE TESTS FOR J=,I3/8H DIF(I)=/(3X,1P6E11.3))
C ---------------------------------- COEFFICIENTS FOR F S
C -------------------------------------------- A'S AND B'S
450 IF(NEWPR.EQ.1) GO TO 337
IF(ENT.LE.TINY) GO TO 335
DO 484 I=2,NM2
A(I)=AMAX1(0.,DIF(I)-HCON(I),-CON(I))
484 B(I+1)=A(I)+CON(I)
GO TO 337
335 DO 338 I=2,NM2
A(I)=DIF(I)
338 B(I+1)=A(I)
337 CONTINUE
B(2)=AMAX1((TIF+RMI)*DXDPEI,0.)
A(NM1)=AMAX1((TEF-RME)*DXDPEI,0.)
C -------------------------------------------- C'S AND D'S
GO TO (501,502,503), KSOURC
C ------------------------------ KSOURC=1, GENERAL
501 SI2=SI(2)
SINM1=SI(NM1)
DO 485 I=2,NM1
IJ=I+IDJ
D(I)=A(I)+B(I)+BOM(I)-SIP(I)
cc write(6,*)'i,ij,f(ij)',i,ij,f(ij)
485 C(I)=F(IJ)*BOM(I)+SI(I)
GO TO 504
C ------------------------------ KSOURC=2, NO SIP
502 SI2=SI(2)
SINM1=SI(NM1)
DO 505 I=2,NM1
IJ=I+IDJ
D(I)=A(I)+B(I)+BOM(I)
cc write(6,*)'i,ij,f(ij)',i,ij,f(ij)
505 C(I)=F(IJ)*BOM(I)+SI(I)
GO TO 504
C ------------------------------ KSOURC=3, NO SIP OR SI
503 SI2=0.
SINM1=0.
DO 506 I=2,NM1
IJ=I+IDJ
D(I)=A(I)+B(I)+BOM(I)
cc write(6,*)'i,ij,f(ij)',i,ij,f(ij)
506 C(I)=F(IJ)*BOM(I)
C ------------------------------
504 C(2)=C(2)-TIF*FDIFI*DXDPEI
C(NM1)=C(NM1)-TEF*FDIFE*DXDPEI
IF(KIN.GT.1) GO TO 486
IF(IBIN(J).EQ.1) GO TO 486
B(2)=0.
ccc C(2)=F(120)*BOM(2)+SI2+RJTOTI(J)*DXDPEI
C(2)=F(I2J)*BOM(2)+SI2+RJTOTI(J)*DXDPEI
D(2)=D(2)-TIF*DXDPEI
486 IF(KEX.GT.1) GO TO 491
IF(IBEX(J).EQ.1) GO TO 491
A(NM1)=0.
C(NM1)=F(INM1J)*BOM(NM1)+SINM1-RJTOTE(J)*DXDPEI
D(NM1)=D(NM1)-TEF*DXDPEI
491 CONTINUE
IF(ITEST.EQ.1) GO TO 464
WRITE(6,405) (A(I),I=2,NM1)
WRITE(6,406) (B(I),I=2,NM1)
WRITE(6,407) (C(I),I=2,NM1)
WRITE(6,408) (C(I),I=2,NM1)
C ------------------------------ SOLVE FOR DOWNSTREAM F'S
464 C(2)=(B(2)*F(I1J)+C(2))/D(2)
D(2)=A(2)/D(2)
DO 465 I=3,NM1
T=1./(D(I)-B(I)*D(I-1))
D(I)=A(I)*T
465 C(I)=(B(I)*C(I-1)+C(I))*T
DO 466 IDASH=1,NM2
I=N-IDASH
IJ=I+IDJ
466 F(IJ)=D(I)*F(IJ+1)+C(I)
C --------------------------- ADJUST F(1,J) AND F(N,J)
IF(KIN-2) 467,460,469
467 IF(IBIN(J).EQ.1) GO TO 468
F(I1J)=FDIFI+F(I2J)+(RJTOTI(J)-F(I1J)*RMI)/TIF
GO TO 460
469 F(I1J)=F(I2J)
GO TO 460
468 RJTOTI(J)=TIF*(F(I1J)-F(I2J)-FDIFI)+RMI*F(I1J)
460 IF(KEX-2) 471,470,473
471 IF(IBEX(J).EQ.1) GO TO 472
F(INJ)=FDIFE+F(INM1J)-(RJTOTE(J)-F(INJ)*RME)/TEF
GO TO 470
473 F(INJ)=F(INM1J)
GO TO 470
472 RJTOTE(J)=TEF*(F(INM1J)+FDIFE-F(INJ))+RME*F(INJ)
470 IF(ITEST.EQ.1) GO TO 480
WRITE(6,476) J,(F(I+IDJ),I=1,N)
476 FORMAT(6H F(I,,I2,1H)/(3X,1P6E11.3))
480 CONTINUE
C ------------------------------------------
481 XU=XD
PSII=PSII-RMI*DX
PSIE=PSIE-RME*DX
PEI=PSIE-PSII
ISTEP=ISTEP+1
END
SUBROUTINE PHYS
common/dbg/debug
logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
COMMON/COMB/
1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C
DIMENSION DUDY(20),EL(20),YEDGE(6)
save
C
c#### should be in blockdata
ccccc DATA KUDIF/-1/
C ------------------------------------------------------------------
CHAPTER A ---------- PHYSU ---------- PHYSU ---------- PHYSU -----------
if(debug) write(6,*)'>>> phys'
ENTRY PHYSU
if(debug) write(6,*)'>>> physu'
C ---------------------------------------- LAMINAR VISCOSITY
C SQUARE-ROOT FORMULA, WITH WEIGHTING ACCORDING TO MASS FRACTION
DO 110 I=1,N
110 EMU(I)=(VISFU*F(I,JF)+VISOX*F(I,JOX)+VISPR*F(I,JPR))*
1 SQRT(F(I,JTE))
IF(MODEL.EQ.1) GO TO 209
C ------------------------------------------------------------------
C ----------------------------- MIXING LENGTH MODEL OF TURBULENCE
IF(KUDIF.EQ.ISTEP) GO TO 102
UMAX=U(1)
UMIN=U(1)
DO 101 I=2,N
UMAX=AMAX1(UMAX,U(I))
101 UMIN=AMIN1(UMIN,U(I))
UDIF=UMAX-UMIN
102 HUDIF=.5*UDIF
DUDYMN=FR*UDIF/Y(N)
C
DO 105 I=2,NM1
105 DUDY(I)=ABS(U(I+1)-U(I-1))/(Y(I+1)-Y(I-1))
K=1
EX=DUDY(2)-DUDYMN
IF(EX.LT.0.) GO TO 103
YEDGE(K)=0.
K=2
103 DO 104 I=3,NM1
EXL=EX
EX=DUDY(I)-DUDYMN
IF(EX*EXL.GE.0.) GO TO 104
YEDGE(K)=.5*(Y(I)+Y(I-1))
IF(K.EQ.6) GO TO 107
K=K+1
104 CONTINUE
IF(EX.LT.0.) GO TO 108
YEDGE(K)=Y(N)
IF(K.EQ.6) GO TO 107
K=K+1
108 CONTINUE
DO 106 KAY=K,6
106 YEDGE(KAY)=Y(N)
107 CONTINUE
EL12=(YEDGE(2)-YEDGE(1))*ALMG
EL34=(YEDGE(4)-YEDGE(3))*ALMG
EL56=(YEDGE(6)-YEDGE(5))*ALMG
EL23=.5*(EL12+EL34)
EL45=.5*(EL34+EL56)
ASSIGN 119 TO K
DO 130 I=2,NM1
YVALUE=Y(I)
GO TO K, (119,121,123,125,127,129)
119 IF(YVALUE.LT.YEDGE(1)) GO TO 120
ASSIGN 121 TO K
121 IF(YVALUE.LT.YEDGE(2)) GO TO 122
ASSIGN 123 TO K
123 IF(YVALUE.LT.YEDGE(3)) GO TO 124
ASSIGN 125 TO K
125 IF(YVALUE.LT.YEDGE(4)) GO TO 126
ASSIGN 127 TO K
127 IF(YVALUE.LT.YEDGE(5)) GO TO 128
ASSIGN 129 TO K
GO TO 129
120 EL(I)=0.
GO TO 130
122 EL(I)=EL12
GO TO 130
124 EL(I)=EL23
GO TO 130
126 EL(I)=EL34
GO TO 130
128 EL(I)=EL45
GO TO 130
129 EL(I)=EL56
C ------------------------------ UPPER LIMITS TO MIXING LENGTH
130 EL(I)=AMIN1(EL(I),HUDIF/(DUDY(I)+TINY))
IF(KIN.NE.1) GO TO 141
DO 142 I=2,NM1
142 EL(I)=AMIN1(EL(I),AK*Y(I))
141 IF(KEX.NE.1) GO TO 143
DO 144 I=2,NM1
144 EL(I)=AMIN1(EL(I),AK*(Y(N)-Y(I)))
143 CONTINUE
C
C ------------------------------------------- VISCOSTITES
C -------------------------------------------TURBULENT CONTRIBUTION
cccc 200 DO 201 I=2,NM1
DO 201 I=2,NM1
DUDYL=DUDY(I)*EL(I)
UDMIN=UFAC*U(I)
DUDYL=AMAX1(DUDYL,UDMIN)
EMUT=RHO(I)*EL(I)*DUDYL
C ----- SIMPLE ADDITION OF THE TURBULENT AND LAMINAR CONTRIBUTIONS
EMU(I)=EMU(I)+EMUT
201 CONTINUE
C
C ------------------------------------------------MOMENTUM SOURCE
209 AGRVDX=AGRAV*DX
RPRLST=1.
MOMSOU=1
IF(ABS(DP).GT.TINY) GO TO 204
IF(ABS(AGRAV).GT.TINY) GO TO 204
MOMSOU=0
RETURN
204 continue
DO 210 I=2,NM1
210 SI(I)=ADPEI(I)*(AGRVDX*(RHO(N)-RHO(I))-DP)
GO TO 9000
C
C ------------------------------------------------------------------
chapter b ---------- physf ---------- physf ---------- physf -----------
ENTRY PHYSF
if(debug) write(6,*)'>>> physf'
IF(MODEL.EQ.2) GO TO 312
RECPR=RECPRL(J)
GO TO 310
312 RECPR=RECPRT(J)
310 NEWPR=1
IF(ABS(RECPR-RPRLST).LT.1.E-10) GO TO 314
NEWPR=2
DO 313 I=2,NM2
313 DIF(I)=DIFU(I)*RECPR
RPRLST=RECPR
c ------------------------------------------ kinetic heating source
314 IF(J.NE.JH) GO TO 3000
IF(ABS(RECPR-1.).LT.1.E-10) GO TO 320
IF(NOVEL.NE.1) GO TO 321
320 KSOURC=3
RETURN
321 DUSQP=0.
USQP=U(2)**2
DO 322 I=2,NM2
USQ=U(I+1)**2
DUSQ=(DIFU(I)-DIF(I))*(USQ-USQP)
SI(I)=.5*(DUSQ-DUSQP)
si(i)=0.0
DUSQP=DUSQ
322 USQP=USQ
KSOURC=2
SI(NM1)=-.5*DUSQP
si(nm1)=0.0
GO TO 9000
C
C ------------------------------------------------ FUEL SOURCE
3000 IF(J.NE.JF) GO TO 4000
IF(INERT.EQ.2) GO TO 342
KSOURC=3
RETURN
342 KSOURC=1
IF(MODEL.NE.1) GO TO 352
T1=DX*PREEXP*PRESS**2
T2=.5/STOICH
DO 344 I=2,NM1
FUBRNT=T2*(ABS(F(I,JP))-F(I,JP))
IF(F(I,JF).GT.FUBRNT) GO TO 346
SIP(I)=-BIG
GO TO 344
346 EXPO=EXP(-ARRCON/F(I,JTE))
F(I,JOX)=AMAX1(0.,F(I,JP)+F(I,JF)*STOICH)
TERM=-T1*EXPO*ADPEI(I)*F(I,JOX)
SIP(I)=TERM/(1.-FUBRNT/F(I,JF))
344 SI(I)=-SIP(I)*FUBRNT
GO TO 9000
352 T2=.5/STOICH
T3=1./(PHIB-PHIC)
T4=-PHIC*T3
C ----- RATE CONTROLLED BY EDDY BREAK-UP
CEBUDX=CEBU*DX
DO 353 I=2,NM1
FUBRNT=T2*(ABS(F(I,JP))-F(I,JP))
FUUNBT=T3*F(I,JP)+T4
IF(F(I,JF).LT.FUUNBT) then
SIP(I)=-ADPEI(I)*CEBUDX*DUDY(I)*RHO(I)*(FUUNBT-F(I,JF))/
1 (FUUNBT-FUBRNT+TINY)
else
SIP(I)=-BIG
endif
SI(I)=-SIP(I)*FUBRNT
cc write(6,*)'i,fubrnt,fuunbt',i,fubrnt,fuunbt
cc write(6,*)'i,si(i),sip(i)',i,si(i),sip(i)
353 continue
GO TO 9000
C
C ---------------------------------------------- JP (PHI)
4000 IF(J.NE.JP) GO TO 9000
KSOURC=3
RETURN
C
9000 IF(ITEST.EQ.1) RETURN
WRITE(6,9001) J,(SI(I),I=2,NM1)
WRITE(6,9002) (SIP(I),I=2,NM1)
9001 FORMAT(18H PHYS TESTS FOR J=,I3/7H SI(I)=/(3X,1P6E11.3))
9002 FORMAT(8H SIP(I)=/(3X,1P6E11.3))
END
c-----------------------------------------------------------------------
SUBROUTINE OUTPUT
common/dbg/debug
logical debug
C/FEB.1977 ---------- G E N M I X ---------- COPYRIGHT, D.B.SPALDING ---
COMMON/COMA/
1 ADPEI(20),BIG,BOM(20),CSALFA,DIF(20),DIFU(20),DP,DX,DXLAST,
2 EMU(20),F(20,6),IBEX(3),IBIN(3),IDIMF,IFIN,ISTEP,ITEST,J,
3 JUSTEX,JUSTIN,KEX,KIN,KRAD,KSOURC,MOMSOU,N,NEWPR,NF,NM1,NM2,
4 NM3,NOVEL,OM(20),OMINT(20),PEI,PSIE,PSII,R(20),RECRU(20),
5 RECYDF(20),RHO(20),RJTOTE(3),RJTOTI(3),RME,RMI,SI(20),SIP(20),
6 TAUE,TAUI,TINY,U(20),XD,XU,Y(20),YE,YI
COMMON/COMB/
1 AK,AGRAV,AHEX,AHIN,ALMG,ALMGD(4),ARRCON,AUEX,BHEX,BHIN,BUEX,
2 CEBU,CFU,CHEX,CHIN,CMIX,COX,CPR,CUEX,DA1,DA2,DPDX,DXINC,
3 DXMAX,DXPSI,DXRAT,DXRE,DXY,ENTHA,ENTHB,ENTHC,ENTHD,EWALL,
4 FACE,FACEXP,FACI,FLOB,FLOC,FR,FRA,FUA,FUB,FUC,FUD,GAMMA,
5 GASCON,H,HDIV,HEX0,HFU,HIN0,ILPLOT,INERT,IRUN,ITPLOT,JF,JH,
6 JOX,JP,JPR,JTE,KASE,KIND,KUDIF,LASTEP,MODEL,NPLOT,NPROF,
7 NSTAT,OMPOW,OXA,OXB,OXC,OXD,PEILIM,PHIA,PHIB,PHIC,PHID,
8 PREEXP,PRESS,PRL(3),PRLAM,PRTURB,RATE,RATI,RECPRL(3),
9 RECPRT(3),REY,STOICH,TA,TB,TC,TD,TWALL,UA,UB,UBAR,UC,UD,
1 UDIF,UEX0,UFAC,UFLUX,ULIM,VISFU,VISMIX,VISOX,VISPR,WFU,WMIX,
2 WOX,WPR,XEND,XHEX0,XHIN0,XOUT,XUEX0,XULAST
C
DIMENSION LAB(6),OUT(6),TITLE(3,4),
1 XLPLOT(100),YLAXIS(12),YLPLOT(100,12),
2 XTPLOT(20),YTAXIS(4),YTPLOT(20,4),
3 DFE(3),DFI(3),FLUX(3),STANE(3),STANI(3)
character*10 lab
character*8 ylaxis,ytaxis,title
character*4 xtaxis,xlaxis
save
C
chapter a ----------------------- initial data for printout ------------
c -------------------- cross-stream output (profile) data ----------
c ----- assign kout= no. of variables, and output labels lab(k)
DATA KOUT/6/
DATA (LAB(K),K=1,6)/'Y','U VEL','TEMP','FUEL','OXYG','OX-FU.STOI'/
C
c -------------------- transverse (cross-stream) plot data
c ----- assign nyt= no. of variables to be plotted
c --- insert dimensions, ensure that itdim.ge.n.and.jtdim.ge.nyt
DATA NYT/4/, ITDIM,JTDIM/20,4/
c ----- assign labels for plot axes
DATA XTAXIS/'Y(I)'/
DATA (YTAXIS(K),K=1,4)/'U VEL ','TEMP ','FUEL ',
1 'OXYG '/
C
c -------------------- longitudinal (down-stream) plot data
c ----- assign nyl= no. of variables to be plotted
c -- insert dimensions, ensure that ildim.ge.lastep.and.jldim.ge.nyl
DATA NYL/12/ILDIM,JLDIM/100,12/
c!!!! two lines were missing
c ----- assign labels for plot axes
DATA XLAXIS/'XU'/
DATA (YLAXIS(K),K=1,12)/'U(1) ','T(1) ','FU(1) ',
1 'OX(1) ','N,R OR Y','1,R(1) ',
1 '2,PEI ','3,RME ','4,FLUXFU',
1 '5,DPDX ','6,RATE ','7,FACE '/
C
C -------------------- TITLE DATA
DATA TITLE/'AXI-SYMM','ETRICAL ','FLOW ',
1 'PLANE FL','UW ',' ',
2 'RADIALLT','-OUTWARD',' FLOW ',
3 'VARIABLE',' CSALFA ',' '/
C
if(debug) write(6,*)'>>> output'
CHAPTER B ----------------------------- HEADINGS -----------------------
IF(ISTEP.GT.0) GO TO 1102
WRITE(6,1103) (TITLE(I,KIND),I=1,3)
1103 FORMAT(1h ,23HGENMIX, FEBRUARY 1977, ,3A8/1H ,
1 'COMBUSTION OF METHANE AND AIR IN A DUCT.')
C
PRESS1=PRESS
TEM=.5*(R(1)+R(N))
EMU1=(VISFU*F(1,JF)+VISOX*F(1,JOX)+VISPR*F(1,JPR))*
1 SQRT(F(1,JTE))
C
REY=PEI/(EMU1*TEM)
EQRAT=0.0
IF(INERT.NE.1) EQRAT=FLOB*STOICH/(FLOC+TINY)/OXC
AMACH=SQRT(PEI*UB/(GAMMA*PRESS*TEM))
C
WRITE(6,1013) KASE,IRUN,KIND,KRAD,CSALFA,MODEL,INERT,NOVEL
1013 FORMAT(1h ,5H KASE,5H IRUN,5H KIND,5H KRAD,7H CSALFA,6H MODEL,
1 6H INERT,6H NOVEL/1X,I4,3I5,F7.3,3I6/)
C
WRITE(6,1014) OMPOW, (OM(I),I=1,N)
1014 FORMAT(1h ,18H OM(I), FOR OMPOW=,F6.3/(1X,1P6E11.3))
C
WRITE(6,1010)
1 HEX0,XHEX0,AHEX,BHEX,CHEX,
2 HIN0,XHIN0,AHIN,BHIN,CHIN,
3 UEX0,XUEX0,AUEX,BUEX,CUEX,
4 XEND,XOUT,XULAST,HDIV,AGRAV
1010 FORMAT(/1h ,
1 4X,4HHEX0,6X,5HXHEX0,7X,4HAHEX,7X,4HBHEX,7X,4HCHEX/1X,1P5E11.3/
2 5X,4HHIN0,6X,5HXHIN0,7X,4HAHIN,7X,4HBHIN,7X,4HCHIN/1X,1P5E11.3/
3 5X,4HUEX0,6X,5HXUEX0,7X,4HAUEX,7X,4HBUEX,7X,4HCUEX/1X,1P5E11.3/
4 5X,4HXEND,7X,4HXOUT,5X,6HXULAST,7X,4HHDIV,6X,5HAGRAV/1X,1P5E11.3)
C
WRITE(6,1011) UA,UB,UC,UD,TA,TB,TC,TD,
2 PRESS,PREEXP,REY,EQRAT,AMACH,ULIM,PEILIM
1011 FORMAT(1h ,4X,2HUA,7X,2HUB,7X,2HUC,7X,2HUD,
1 7X,2HTA,7X,2HTB,7X,2HTC,7X,2HTD/1X,8F9.3/
2 4X,5HPRESS,3X,6HPREEXP,6X,3HREY,4X,5HEQRAT,4X,5HAMACH,5X,4HULIM,
3 3X,6HPEILIM/1X,1P7E9.2)
C
chapter c ---------------- compute output required at each step --------
1102 CONTINUE
UBAR=0.
DO 110 I=2,NM1
110 UBAR=UBAR+BOM(I)*U(I)
UFLUX=PEI*UBAR
DO 115 J=1,NF
FLUX(J)=0.
DO 116 I=2,NM1
116 FLUX(J)=FLUX(J)+BOM(I)*F(I,J)
115 FLUX(J)=PEI*FLUX(J)
C
DO 117 J=1,NF
DFI(J)=FLUX(J)/PEI-F(I,J)
117 DFE(J)=DFI(J)+F(1,J)-F(N,J)
UFLUX=UFLUX-PSIE*U(N)+U(1)*PSII
FLUX(JH)=FLUX(JH)-PSIE*ENTHD+PSII*ENTHA
FLUX(JP)=FLUX(JP)-PSIE*PHID+PSII*PHIA
FLUX(JF)=FLUX(JF)-PSIE*FUD+PSII*FUA
PRESSD=PRESS/PRESS1-1.
C
IF(ISTEP.EQ.0.OR.ILPLOT.EQ.1) GO TO 1105
c ----- assign values for downstream plot
XLPLOT(ISTEP)=XU
YLPLOT(ISTEP,1)=U(1)
YLPLOT(ISTEP,2)=F(1,JTE)
YLPLOT(ISTEP,3)=F(1,JF)
YLPLOT(ISTEP,4)=F(1,JOX)
IF(KIND-1) 111,111,114
111 YLPLOT(ISTEP,5)=R(N)
GO TO 113
114 YLPLOT(ISTEP,5)=Y(N)
113 CONTINUE
YLPLOT(ISTEP,6)=R(1)
YLPLOT(ISTEP,7)=PEI
YLPLOT(ISTEP,8)=RME
YLPLOT(ISTEP,9)=FLUX(JF)
YLPLOT(ISTEP,10)=DPDX
YLPLOT(ISTEP,11)=RATE
YLPLOT(ISTEP,12)=FACE
1105 CONTINUE
C
c ----- tests for printout
c ----- iprint=1 gives single (station) variables,
c iprint=2 adds the array (profile) variables,
c iprint-3 adds the cross-stream plots.
IPRINT=0
IF(MOD(ISTEP,NSTAT).EQ.0) IPRINT=1
IF(MOD(ISTEP,NPROF).EQ.0) IPRINT=2
IF(ISTEP.EQ.0) GO TO 1020
IF(MOD(ISTEP,NPLOT).EQ.0
1 .OR.ISTEP.EQ.JUSTEX.OR.ISTEP.EQ.JUSTIN
2 .OR.ITEST.NE.1.OR.IFIN.NE.1) IPRINT=3
1020 IF(IPRINT.EQ.0) RETURN
C
chapter d ------------------------------ station variables -------------
call writbl
WRITE(6,1030) XU,ISTEP,
1 JUSTIN,JUSTEX,DX,PRESSD,
2 KIN,KEX,DXY,DPDX,
3 PSII,PSIE,DXRE,PEI,
4 RMI,RME,DXINC,
5 R(1),R(N),DXPSI,
5 (J,J=1,4),UFLUX,(FLUX(J),J=1,NF)
1030 FORMAT(1h ,5H*** ,3HXU=,1PE10.3,2X,6HISTEP=,I5/
1 2X,7HJUSTIN=,I10,1X,7HJUSTEX=,I10,5X,3HDX=,1PE10.3,
1 8H PRESSD=,E10.3/
2 5X,4HKIN=,I10,4X,4HKEX=,I10,4X,4HDXY=,E10.3,3X,5HDPDX=,E10.3/
3 4X,5HPSII=,E10.3,3X,5HPSIE=,E10.3,3X,5HDXRE=,E10.3,
3 4X,4HPEI=,E10.3/
4 5X,4HRMI=,E10.3,4X,4HRME=,E10.3,2X,6HDXINC=,E10.3/
5 4X,5HR(1)=,E10.3,3X,5HR(N)=,E10.3,2X,6HDXPSI=,E10.3/
5 23X,3HJ =,4(I6,5X)/7H UFLUX=,E10.3,9H FLUX(J)=,(4E11.3))
C
IF(ISTEP.EQ.0) GO TO 1042
UREF=UBAR
RUREF=PEI/((R(I)+R(N))*.5*Y(N))
URUREF=1./(UREF*RUREF)
C
IF(KIN-2) 1061,1062,1063
1061 TAUID=TAUI*URUREF
DO 1025 J=1,NF
1025 STANI(J)=(RJTOTI(J)-F(1,J)*RMI)/(R(1)*DFI(J)*RUREF+TINY)
WRITE(6,1029) TAUID,(STANI(J),J=1,NF)
1029 FORMAT(1H ,6HTAUID=,1PE10.3,10H STANI(J)=,(4E11.3))
GOTO 1063
1062 WRITE(6,1069) FACI,RATI
1069 FORMAT(1H ,6H FACI=,1PE10.3,6H RATI=,E10.3)
1063 IF(KEX-2) 1081,1082,1044
cccc 1081 TAUED=TAUED*URUREF
1081 TAUED=TAUE*URUREF
DO 1027 J=1,NF
1027 STANE(J)=(RJTOTE(J)-F(N,J)*RME)/(R(N)*DFE(J)*RUREF+TINY)
WRITE(6,1028) TAUED,(STANE(J),J=1,NF)
1028 FORMAT(1H ,6HTAUED=,1PE10.3,10H STANE(J)=,(4E11.3))
GO TO 1044
1082 WRITE(6,1089) FACE,RATE
1089 FORMAT(1H ,6H FACE=,1PE10.3,5H DA2=,E10.3)
GO TO 1042
1044 WRITE(6,1047) DA1,DA2
1047 FORMAT(5H DA1=,1PE10.3,5H DA2=,E10.3)
C
CCHAPTER E ---------------------------- CROSS-STREAM PROFILES ----------
1042 IF(IPRINT.EQ.1) GO TO 1050
C
WRITE(6,1099) (LAB(K),K=1,KOUT)
DO 1091 I=1,N
OUT(1)=Y(I)
OUT(2)=U(I)
OUT(3)=F(I,JTE)
OUT(4)=F(I,JF)
OUT(5)=F(I,JOX)
OUT(6)=F(I,JP)
c ------------------------------ print profiles
1091 WRITE(6,1098) I,(OUT(K),K=1,KOUT)
C
IF(IPRINT.LT.3.OR.ITPLOT.EQ.1) GO TO 1050
c ----- assign cross-stream plots
DO 1073 I=1,N
XTPLOT(I)=Y(I)
YTPLOT(I,1)=U(I)
YTPLOT(I,2)=F(I,JTE)
YTPLOT(I,3)=F(I,JF)
YTPLOT(I,4)=F(I,JOX)
1073 CONTINUE
c --------------------- cross-stream plot output
call writbl
WRITE(6,1096) XU,ISTEP
1096 FORMAT(19H CROSS-STREAM PLOT,,4H XU=,1PE10.3,7H ISTEP=,I4)
CALL PLOTS(XTPLOT,ITDIM,N,XTAXIS,YTPLOT,JTDIM,NYT,YTAXIS)
C
chapter f -------------------------- return or terminate ---------------
1050 IF(IFIN.EQ.1) RETURN
C
call writbl
WRITE(6,112) ISTEP,LASTEP,XU,XULAST,IFIN
112 FORMAT(14H0TERMINATED AT/7H ISTEP=,I5,8H LASTEP=,I5,
1 4H XU=,1PE11.3,8H XULAST=,E11.3,6H IFIN=,I3)
IF (ILPLOT.EQ.1) RETURN
call writbl
c --------------------- downstream plot output
WRITE(6,1054) XU,ISTEP
1054 FORMAT(18H DOWN-STREAM PLOT,,4H XU=,1PE10.3,7H ISTEP=,I4)
CALL PLOTS(XLPLOT,ILDIM,ISTEP,XLAXIS,YLPLOT,JLDIM,NYL,YLAXIS)
1098 FORMAT(1H ,I3,1P10E11.3)
1099 FORMAT(1h ,2X,2HI ,10A11)
END
subroutine writbl
write(6,*)' '
end
c