Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение
Подождите немного. Документ загружается.
280
GV=B(2)
DK=B(3)
WGV=B(4)
WK=B(5)
TK=B(6)
110 TGV=B(7)
WRITE(4,*) ' X=',X,' PGV=',PGV,' WGV=',WGV,' WK=',WK
IF(K.EQ.KK)WRITE(6,*) ' X=',X,' PGV=',PGV,' WGV=',WGV,'
WK=',WK
IF(K.EQ.KK)WRITE(6,*) 'SK=',SK,' TGV=',TGV,' TK=',TK
IF(K.EQ.KK)WRITE(6,*) 'ROGV=',ROGV,' DSM=',DSM,'
DK=',DK
IF(K.EQ.KK)WRITE(6,*) 'DGK=',DGK,' DGKOT=',DGKOT,'
DLTGV=',DLTGV
IF(K.EQ.KK)WRITE(6,*) 'DLDK=',DLDK,'
DLDKOT=',DLDKOT,' DT=',DT
IF(K.EQ.KK)WRITE(6,*) 'MGV=',MGV,' DLD0=',DLD0
IF(K.EQ.KK)K=0
IF(X.LT.Z)GOTO 120
STOP
END
SUBROUTINE RKM(A,DX,N,IS,B,C)
DIMENSION F(10),G(40),C(10),B(10)
IS=IS+1
GOTO(10,30,60,80),IS
10 E=A
DO 20 I=1,N
F(I)=B(I)
G(4*I-3)=C(I)*DX
20 B(I)=F(I)+G(4*I-3)/2.
GOTO 50
30 DO 40 I=1,N
G(4*I-2)=C(I)*DX
40 B(I)=F(I)+G(4*I-2)/2.
50 A=E+DX/2.
GOTO 100
60 DO 70 I=1,N
G(4*I-1)=C(I)*DX
70 B(I)=F(I)+G(4*I-1)
A=E+DX
GOTO 100
80 DO 90 I=1,N
281
G(4*I)=C(I)*DX
B(I)=G(4*I-3)+2.*(G(4*I-2)+G(4*I-1))
90 B(I)=(B(I)+G(4*I))/6.+F(I)
IS=0
100 RETURN
END
Файл исходных данных PIPE.DAT
7 0.3E-3 3. 6.E5 150. 130. 3.1416
873. 373. 0.2 0.01 0.2
20.E-5 8314. 29. 18.
10 2.
N DX Z P0 WGV WK PI
TGV TK0 GG GV GK
DK0 UGC AMG AMV
KK PSD
Приложение №4
Идентификаторы программ LAV21D и LAV27A
GVB,WB,DB- расход,скорость и диаметр пузырьков
TLS- температура насыщения
DBCR- критический диаметр пузырька
DTPL - перегрев жидкости относительно температуры насыщения
при данном давлении
ROVB- плотность пара в пузырьке
ROVTS- плотность пара при температуре насыщения
PSI- коэффициент в формуле (5.13)
PSIA- коэффициент в формуле (5.16)
WBCR- работа образования пузырька критического размера
BN- число Больцмана
GIBSN- число
Гиббса
IGB- интенсивность образования пузырьков критического размера
NBS- поток пузырьков
VB- объем пузырька
TDSR- температура капель, полученная из уравнения энергии
BIO- число Био
FURC- число Фурье
SQ- удельная теплота парообразования
282
LAV21D
PROGRAM LAV21D
REAL IVTS,ILTL0,IVTV,IGB,NBS,
& LAMLS,LAMDD,MU(3),B(10),C(10)
TET(A)=2.*(SIN(A)-A*COS(A))*SIN(A)*EXP(-A**2*FURC)/(A-
SIN(A)*
& COS(A))/A
RO(A)=999.9-0.1475*(A-273.)-2.695E-3*(A-273.)**2
YILTL(A)=(12.475+3.935*(A-273.)+1.3125E-3*(A-273.)**2)*1.E3
SIG(A)=(75.54-0.168*(A-273.))*1.E-3
IV(A)=(2461.1+2.636*(A-273.)-4.88E-3*(A-273.)**2)*1.E3
CPL(A)=4275.5-1.2209*(A-273.)+0.59823E-2*(A-273.)**2
VIS(A)=(613.-4.227*(A-273.)+9.22E-3*(A-273.)**2)*1.E-6
SQV(A)=(2447.7-1.283*(A-273.)-6.26E-3*(A-273.)**2)*1.E3
OPEN (2,FILE='L21D.DAT')
OPEN (8,FILE='REZL21D.DAT')
READ(2,*) N,DX,Z,GCR,TL0,PL0,DBCR,DTPLCR,
& BN,DC0,UGC,AMV,PI,A1,PSIA,B1,LAMLS,KK,DD0,
& AA,B3,PSR,DMU0,EPS,LAMDD,A2
N=3
NP=0
ROL0=RO(TL0)
ILTL0=YILTL(TL0)*GCR
FC0=PI*DC0**2/4.
WL0=GCR/ROL0/FC0
FC=FC0
PL=PL0
TL=TL0
DB=DBCR
RBCR=DBCR/2.
TLS=5.677*PL**.25+273.
SIGL=SIG(TL)
ROVB=(PL+4.*SIGL/DB)*AMV/UGC/TLS
GVB=0.
WL=WL0
WB=WL0
IS=0
B(1)=DB
B(2)=WB
B(3)=WL
J=0
K=0
283
X=0
WRITE(8,*)' X PL WL WB TLS ROSM DTPL DC'
120 X=X+DX
K=K+1
DO 110 LL=1,4
PL=PL0/((PL0/PSR-1.)*(X/Z)**2+1.)
print*,'X,Z,PL ',X,Z,PL
TLS=5.677*PL**.25+273.
IVTS=IV(TLS)
CPLIT=CPL(TL)
GL=GCR-GVB
IF(DTPL.GE.DTPLCR.OR.NP.EQ.1)
& TL=((ILTL0-(IVTS+WB**2/2.)*GVB)/GL-
WL**2/2.)/CPLIT+273.
DTPL=TL-TLS
ROL=RO(TL)
ROLS=RO(TLS)
ROSGR=ROL*(1.-PI*2.**.5/6.)*A2
VISL=VIS(TL)
SIGL=SIG(TL)
ROVB=(PL+4.*SIGL/DB)*AMV/UGC/TLS
IF(NP.EQ.0)FC=GCR/WL/ROL
IF(DTPL.GE.DTPLCR)FC=GVB/WB/ROVB+GL/WL/ROL
IF(DTPL.LE.DTPLCR.AND.NP.EQ.1)FC=GVB/WB/ROVB+GL/WL/R
OL
DC=(4.*FC/PI)**.5
IVTS=IV(TLS)
SIGLS=SIG(TLS)
SQ=SQV(TLS)
ROVTS=1.213E-5*((TLS-273.)/5.677)**3.75
IF(DTPL.LT.DTPLCR.OR.NP.EQ.1)GOTO 130
PSI=RBCR*DTPLCR*ROVTS*SQ/2./SIGLS/TLS
WBCR=16.*PI*(SIGLS**3)*(PSI*TLS/DTPLCR/ROVTS/SQ)**2/3.
GIBSN=WBCR/TLS/BN
IGB=EXP(A1-PSIA*GIBSN)
NBS=IGB*FC*DX
GVB=NBS*PI*DB**3*ROVB/6.
NP=1
130 CPLIT=CPL(TL)
GVB=NBS*PI*DB**3*ROVB/6.
IF(NP.EQ.0)ROSM=ROL
284
IF(NP.EQ.1)ROSM=ROVB*ROL/(ROVB+GVB*(ROL-
ROVB)/GCR)
DC=(4.*FC/PI)**0.5
XV=GVB/GCR
REN=WL*DC*ROSM/VISL
CF=.0791/(REN**.25)
TW=CF*ROSM*WL**2/2.
PPL=-2.*PL0*(PL0/PSR-1.)*X/(Z**2*((PL0/PSR-1.)*(X/Z)**2
& +1.)**2)
PB=PL+4.*SIGLS/DB
CAP=B3/(PB**.56)
REB=ABS(WL-WB)*DB*ROL/VISL
IF(REB.GT.0.AND.REB.LE.1.)CX=24./REB
IF(REB.GT.1.AND.REB.LE.390.)CX=24./(REB**.667)
IF(REB.GT.390.)CX=.45
VB=PI*DB**3/6.
C(1)=8.*LAMLS*CPLIT*ROLS*ABS(TL-TLS)*(TL-
TLS)*B1*NP/
& (SQ**2*ROVTS**2*DB*WB)
IF(NP.EQ.0)C(2)=-FC*PPL/GCR-TW*PI*DC/GCR
IF(NP.EQ.1)C(2)=(-PPL/WB+.75*CX*ROL*WB*(WL/WB-1.)**2/
& DB+(WL-
WB)*(1.+.5*ROL*WL*NBS*VB/WB/GL)*3.*ROVB*C(1)/DB
& -.5*ROL*WL*(FC*PPL+TW*PI*DC)/WB/GL)/
& (ROVB+.5*ROL*(1.+WL*GVB/WB/GL))
IF(NP.EQ.0)C(3)=-FC*PPL/GCR-TW*PI*DC/GCR
IF(NP.EQ.1)C(3)=(-GVB*C(2)+(WL-
WB)*NBS*VB*3.*ROVB*C(1)/
& DB-FC*PPL-TW*PI*DC)/GL
CALL RKM(X,DX,N,IS,B,C)
DB=B(1)
WB=B(2)
110 WL=B(3)
IF(K.EQ.KK)THEN
J=J+1
IF(DTPL.LT.DTPLCR) THEN
WRITE(8,200)X,PL,WL,WB,TLS,ROSM,DTPL,DC
IF(J.EQ.10) THEN
WRITE(8,*)' X PL WL WB TLS ROSM DTPL DC'
ENDIF
ELSE
285
WRITE(8,210)X,PL,WL,WB,TLS,ROSM,DTPL,DC,ROVB,GVB,XV,N
BS,DB
& ,ROSGR,ROSM
IF(J.EQ.10) THEN
WRITE(8,*)' X PL WL WB TLS ROSM DTPL DC ROVB GVB
XV NBS DB
&ROSM'
ENDIF
PRINT *,' X=',X,' ROSM=',ROSM,' ROSGR=',ROSGR
END IF
K=0
J=0
END IF
IF(ROSM.GE.ROSGR) GOTO 120
N=4
PV=PL
GV=GVB
DD=DD0
WV=WB+AA
WD=WL-AA
TDSR=TL
GD=GL
TV=TLS
ROD=ROL
ROV=ROVB
B(1)=GV
B(2)=DD
B(3)=WV
B(4)=WD
PV=PL0/((PL0/PSR-1.)*(X/Z)**2+1.)
PPV=-2.*PL0*(PL0/PSR-1.)*X/(Z**2*((PL0/PSR-1.)
*(X/Z)**2+1.)**2)
TV=5.677*PV**.25+273.
ROV=1.2113E-5*PV**.9375
ROD=RO(TDSR)
FC=GV/ROV/WV+GD/ROD/WD
DC=(4.*FC/PI)**.5
RD=DD/2.
WRITE(8,*)' X PV DD WV WD FC DC XV FIN TV TDSR DTDV
ALFAM'
140 X=X+DX
K=K+1
286
DO 150 LL=1,4
GD=GCR-GV
RD=DD/2.
ROD=999.9-0.1475*(TDSR-273.)-2.695E-3*(TDSR-273.)**2
PV=PL0/((PL0/PSR-1.)*(X/Z)**2+1.)
PPV=-2.*PL0*(PL0/PSR-1.)*X/(Z**2*((PL0/PSR-1.)*(X/Z)**2
& +1.)**2)
ROV=1.21E-5*PV**.9375
TV=5.677*PV**.25+273.
SQ=(2447.7-1.283*(TV-273.)-6.26E-3*(TV-273.)**2)*1.E3
ILTL0=YILTL(TL0)*GCR
IVTV=(2461.1+2.636*(TV-273.)-4.88E-3*(TV-273.)**2)*1.E3
CPDIT=4275.5-1.2209*(TDSR-273.)+.59823E-2*(TDSR-273.)**2
TDSR=((ILTL0-(IVTV+WV**2/2.)*GV)/GD-WD**2/2.)
/CPDIT+273.
VISV=(8.35+3.42E-2*(TV-273.)+2.E-5*(TV-273.)**2)*1.E-6
RED=ABS(WV-WD)*DD*ROV/VISV
IF(RED.GT.0..AND.REB.LE.1.)CX=24./RED
IF(RED.GT.1..AND.RED.LE.390.)CX=24./(RED**.667)
IF(RED.GT.390.)CX=.45
DSTAR=CX*PI*DD**2*ROV*(WV-WD)**2/8.
FC=GV/ROV/WV+GD/ROD/WD
DC=(4.*FC/PI)**.5
WVD=((GV*WV**2+GD*WD**2)/(GV+GD))**.5
SWVD=(2.*CPDIT*(TL0-TV-TV*ALOG(TL0/TV)))**.5
FIN=WVD/SWVD
XV=GV/(GV+GD)
REV=WV*DC*ROV/VISV
CF=.079/REV**.25
TW=CF*ROV*WV**2/2.
CAP=B3/PV**.56
GCV=UGC/AMV
AD=LAMDD/CPDIT/ROD
ALFAM=SQ**2*CAP*ROV/TV/(2.*PI*GCV*TV)**.5
BIO=ALFAM*RD/LAMDD
TAUF=DX/WD
FURC=AD*TAUF/RD**2
MU(1)=2.
MU(2)=5.
MU(3)=8.
DO 20 I=1,3
DMU=DMU0
287
25 MU(I)=MU(I)+DMU
Y1=MU(I)/(BIO-1.)
Y2=-SIN(MU(I))/COS(MU(I))
IF(ABS(Y1-Y2).LT.EPS) GOTO 20
IF(ABS(Y1)-ABS(Y2)) 9,9,10
9 GOTO 25
10 MU(I)=MU(I)-DMU
DMU=DMU/2.
GOTO 25
20 CONTINUE
TETA1=TET(MU(1))
TETA2=TET(MU(2))
TETA3=TET(MU(3))
TETA=TETA1+TETA2+TETA3
DTDV=(TDSR-TV)*TETA
DJDV=ALFAM*DTDV/SQ
C(1)=6.*DJDV*GD/WD/ROD/DD
C(2)=-2.*DJDV/ROD/WD
C(3)=-(PPV*GV/WV/ROV+6.*DSTAR*GD/WD/PI/DD**3/ROD
& -(WV-WD)*C(1))/GV
C(4)=-PPV/ROD/WD+6.*DSTAR/PI/DD**3/ROD/WD
150 CALL RKM(X,DX,N,IS,B,C)
GV=B(1)
DD=B(2)
WV=B(3)
WD=B(4)
IF(K.EQ.KK) THEN
J=J+1
WRITE(8,220)
X,PV,DD,WV,WD,FC,DC,XV,FIN,TV,TDSR,DTDV,ALFAM
IF(J.EQ.10) THEN
WRITE(8,*)' X PV DD WV WD FC DC XV FIN TV TDSR DTDV
ALFAM'
J=0
PRINT *,' X=',X,' WV=',WV
ENDIF
K=0
ENDIF
IF(X.LT.Z) GOTO 140
200 FORMAT(F6.3,F10.0,2F8.3,F6.0,F8.3,F6.1,F6.3)
210
FORMAT(F6.3,F10.0,2F8.3,F6.0,F8.3,F6.1,F6.3,2F8.3,F6.3,F11.0,
288
& E10.3,2F8.3)
220 FORMAT(F6.3,F10.0,E10.3,2F8.3,E10.2,3F6.3,2F6.0,F11.0)
WRITE(8,220)
X,PV,DD,WV,WD,FC,DC,XV,FIN,TV,TDSR,ALFAM
STOP
END
Файл исходных данных L21D.DAT
3 0.1E-3 0.2 0.3 428. 6.E5 2.E-8 3.
1.38E-23 3.E-2 8314. 18. 3.1416 91. 15.3E-3 0.185
0.68 10 1.5E-5 0.0 70. 1.E5 1.E-1 2.E-3 0.68 1.0
N DX Z GKR TL0 PL0 DBCR DTPLCR
BN DC0 UGC AMV PI A1 PSIA B1
LAMLS KK DD0 AA B3 PSR DMU0 EPS LAMDD A2
LAV27A
PROGRAM LAV27A
C Расчет сопла Лаваля с парогенерирующей решеткой
C на входе (прямая задача).
DIMENSION B(10),C(10)
REAL JVI,JVC,LAMV
VIS(A)=(613.-4.227*(A-273.)+9.22E-3*(A-273.)**2)*1E-6
RO(A)=999.9-0.1475*(A-273.)-2.695E-3*(A-273.)**2
YIL(A)=(12.475+3.935*(A-273.)+1.3125E-3*(A-273.)**2)*1.E3
SIG(A)=(75.54-0.168*(A-273.))*1.E-3
SQV(A)=(2447.7-1.283*(A-273.)-6.26E-3*(A-273.)**2)*1.E3
CPL(A)=4275.5-1.2209*(A-273.)+0.59823E-2*(A-273.)**2
YIV(A)=(2461.1+2.636*(A-273.)-4.88E-3*(A-273.)**2)*1.E3
OPEN (2,FILE='L27.DAT')
OPEN (8,FILE='REZL27A.DAT')
READ(2,*) N,DX,GCR,GV1,TL0,PL0,TD1,PV1,
& DC1,DCM,DCSR,ZCC,ZCM,ZCD,UGC,AMV,PI,
& KK,DD1,B3,WV1
Z=ZCC+ZCM+ZCD
FC1=PI*DC1**2/4.
TD=TD1
289
PV=PV1
WV=WV1
GV=GV1
FC=FC1
DD=DD1
ROD=RO(TD)
TV=5.677*PV**0.25+273.
ROV=1.213E-5*PV**0.9375
FV=GV/ROV/WV
FD=FC-FV
GD=GCR-GV
WD=GD/ROD/FD
CPD=CPL(TD)
N=6
IS=0
B(1)=PV
B(2)=GV
B(3)=DD
B(4)=WV
B(5)=WD
B(6)=TD
J=0
K=0
X=0
WRITE(8,*)' X PV GV GD DD WV WD
& TV TD FC DC FIN'
140 X=X+DX
K=K+1
IF(X.GT.Z)GOTO 16
DO 150 LL=1,4
GD=GCR-GV
ROD=RO(TD)
ROV=1.213E-5*PV**.9375
TV=5.677*PV**.25+273.
SQ=SQV(TV)
VISV=(8.35+3.42E-2*(TV-273.)+2E-5*(TV-273.)**2)*1E-6
IF(X.LT.ZCC)DC=DC1-(DC1-DCM)*X/ZCC
IF(X.GE.ZCC.AND.X.LE.(ZCC+ZCM))DC=DCM
IF(X.GT.(ZCC+ZCM))DC=DCM+(DCSR-DCM)*(X-
(ZCC+ZCM))/ZCD
FC=PI*DC**2/4.