Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение
Подождите немного. Документ загружается.
220
C(1)=C(1)*D(1)**2/(X(4)-X(1))
C(N)=-C(N)*D(N-1)**2/(X(N)-X(N-3))
15 DO 20 I=2,N
T=D(I-1)/B(I-1)
B(I)=B(I)-T*D(I-1)
C(I)=C(I)-T*C(I-1)
20 CONTINUE
C(N)=C(N)/B(N)
DO 30 IB=1,NM1
I=N-IB
C(I)=(C(I)-D(I)*C(I+1))/B(I)
30 CONTINUE
B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0*C(N))
DO 40 I=1,NM1
B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0*C(I))
D(I)=(C(I+1)-C(I))/D(I)
C(I)=3.0*C(I)
40 CONTINUE
C(N)=3.0*C(N)
D(N)=D(N-1)
RETURN
50 B(1)=(Y(2)-Y(1))/(X(2)-X(1))
C(1)=0.0
D(1)=0.0
B(2)=B(1)
C(2)=0.0
D(2)=0.0
RETURN
END
SUBROUTINE SEVAL(N,U,X,Y,B,C,D,S,SP,SPP)
DIMENSION X(9),Y(9),B(9),C(9),D(9)
DATA I/1/
IF(I.GE.N) I=1
IF(U.LT.X(I)) GO TO 10
IF(U.LE.X(I+1)) GO TO 30
10 I=1
J=N+1
20 K=(I+J)/2
IF(U.LT.X(K)) J=K
IF(U.GE.X(K)) I=K
IF(J.GT.I+1) GO TO 20
30 DX=U-X(I)
221
S=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
SP=B(I)+DX*(2.*C(I)+3.*DX*D(I))
SPP=2.*C(I)+6.*DX*D(I)
RETURN
END
Файл исходных данных GM3.DAT
5 0.3E-3 0.15 30. 6.E5 0.1 3.E-2 1.E5 3.1416 10
100. 110. 120. 130. 140. 150. 160. 170. 180.
1.01E5 1.43E5 1.98E5 2.70E5 3.61E5 4.76E5 6.18E5 7.92E5 10.03E5
1.04E-3 1.05E-3 1.06E-3 1.07E-3 1.08E-3 1.09E-3 1.1E-3 1.11E-3
1.13E-3
1.673 1.21 0.8917 0.6683 0.5087 0.3926 0.3068 0.2426 0.1939
419.E3 461.E3 503.E3 546.E3 589.E3 632.E3 675.E3 719.E3 763.E3
2676.E3 2691.E3 2706.E3 2721.E3 2734.E3 2746.E3 2758.E3 2769.E3
2778.E3
1.307E3 1.418E3 1.528E3 1.635E3 1.739E3 1.842E3 1.943E3 2.042E3
2.1395E3
7.355E3 7.239E3 7.13E3 7.027E3 6.93E3 6.838E3 6.751E3 6.667E3
6.5858E3
2257.E3 2230.E3 2202.E3 2175.E3 2145.E3 2114.E3 2082.E3 2050.E3
2015.E3
4.22E3 4.23E3 4.24E3 4.26E3 4.28E3 4.3E3 4.33E3 4.35E3 4.38E3
71.02E3 65.91E3 61.3E3 57.11E3 53.3E3 49.82E3 46.63E3 43.66E3
40.94E3
4.22E3 4.23E3 4.25E3 4.27E3 4.29E3 4.31E3 4.35E3 4.38E3 4.417E3
(T(I),I=1,9)
(P(I),I=1,9)
(VL(I),I=1,9)
(VV(I),I=1,9)
(ENL(I),I=1,9)
(ENV(I),I=1,9)
(SL(I),I=1,9)
(SV(I),I=1,9)
(R(I),I=1,9)
(CVL(I),I=1,9)
(CVP(I),I=1,9)
(CPL(I),I=1,9)
N, DX, Z, WLV0 ,P0, XV0 D0, PSR, PI, KK
222
PROGRAM GM4
C Используются подпрограммы: spline,
C seval,rkm.
C Изоэнтропийное расширение гомогенного
C потока в сопле Лаваля - прямая задача
C ( задана геометрия канала вдоль сопла и
C начальная степень сухости-XV0 ).
DIMENSION B1(10),C1(10)
DIMENSION P(9),T(9),VL(9),VV(9),R(9),ENL(9),ENV(9),
& SL(9),SV(9),BSL(9),CSL(9),DSL(9),BSV(9),CSV(9),DSV(9),
& CPW(9),CVL(9),CVP(9),BT(9),CT(9),DT(9),BVL(9),
& DVL(9),BVV(9),CVV(9),DVV(9),BR(9),CR(9),
& DR(9),BENL(9),CENL(9),DENL(9),BENV(9),CENV(9),
& DENV(9),BCPW(9),CCPW(9),DCPW(9),BCVL(9),CCVL(9),
& DCVL(9),BCVP(9),CCVP(9),DCVP(9),BP(9),CP(9),DP(9)
REAL MAX
OPEN(2,FILE='GM4.DAT')
OPEN(4,FILE='PRN')
OPEN(6,FILE='CON')
OPEN(8,FILE='REZGM4.DAT')
READ(2,*) N,DX,Z,WLV0,P0,XV0,PSR,PI,KK,
& DC0,DCM,DCSR,ZCC,ZCM,ZCD
READ(2,*) (T(I),I=1,9)
READ(2,*) (P(I),I=1,9)
READ(2,*) (VL(I),I=1,9)
READ(2,*) (VV(I),I=1,9)
READ(2,*) (ENL(I),I=1,9)
READ(2,*) (ENV(I),I=1,9)
READ(2,*) (SL(I),I=1,9)
READ(2,*) (SV(I),I=1,9)
READ(2,*) (R(I),I=1,9)
READ(2,*) (CVL(I),I=1,9)
READ(2,*) (CVP(I),I=1,9)
READ(2,*) (CPW(I),I=1,9)
Z=ZCC+ZCM+ZCD
CALL SPLINE(9,P,T,BT,CT,DT)
CALL SPLINE(9,T,P,BP,CP,DP)
CALL SPLINE(9,T,VL,BVL,CVL,DVL)
CALL SPLINE(9,T,VV,BVV,CVV,DVV)
CALL SPLINE(9,T,ENL,BENL,CENL,DENL)
CALL SPLINE(9,T,ENV,BENV,CENV,DENV)
CALL SPLINE(9,T,SL,BSL,CSL,DSL)
223
CALL SPLINE(9,T,SV,BSV,CSV,DSV)
CALL SPLINE(9,T,R,BR,CR,DR)
CALL SPLINE(9,T,CVL,BCVL,CCVL,DCVL)
CALL SPLINE(9,T,CVP,BCVP,CCVP,DCVP)
CALL SPLINE(9,T,CPW,BCPW,CCPW,DCPW)
N=4
FC0=PI*DC0**2/4.
CALL SEVAL(9,PS,P,T,BT,CT,DT,TS,P1TS,P2TS)
CALL SEVAL(9,P0,P,T,BT,CT,DT,T0,P1T0,P2T0)
CALL SEVAL(9,TS,T,VL,BVL,CVL,DVL,VLS,P1VLS,P2VLS)
CALL SEVAL(9,T0,T,VL,BVL,CVL,DVL,VL0,P1VL0,P2VL0)
CALL SEVAL(9,TS,T,VV,BVV,CVV,DVV,VVS,P1VVS,P2VVS)
CALL SEVAL(9,T0,T,VV,BVV,CVV,DVV,VV0,P1VV0,P2VV0)
CALL
SEVAL(9,TS,T,ENL,BENL,CENL,DENL,ENLS,P1ENLS,P2ENLS)
CALL
SEVAL(9,TS,T,ENV,BENV,CENV,DENV,ENVS,P1ENVS,P2ENVS)
CALL SEVAL(9,TS,T,SL,BSL,CSL,DSL,SLS,P1SLS,P2SLS)
CALL SEVAL(9,T0,T,SL,BSL,CSL,DSL,SL0,P1SL0,P2SL0)
CALL SEVAL(9,TS,T,SV,BSV,CSV,DSV,SVS,P1SVS,P2SVS)
CALL SEVAL(9,T0,T,SV,BSV,CSV,DSV,SV0,P1SV0,P2SV0)
CALL SEVAL(9,TS,T,R,BR,CR,DR,RS,P1RS,P2RS)
CALL
SEVAL(9,TS,T,CVL,BCVL,CCVL,DCVL,CVLS,P1CVLS,P2CVLS)
CALL
SEVAL(9,TS,T,CVP,BCVP,CCVP,DCVP,CVPS,P1CVPS,P2CVPS)
CALL
SEVAL(9,T0,T,CPW,BCPW,CCPW,DCPW,CP0W,P1CP0W,P2CP0W)
CALL
SEVAL(9,TS,T,CPW,BCPW,CCPW,DCPW,CPSW,P1CPSW,P2CPSW)
FC=FC0
TS=T0
WLV=WLV0
PS=P0
DC=DC0
V0=(1.-XV0)*VL0+XV0*VV0
V=V0
GS=FC0*WLV0/V0
S0=(1.-XV0)*SL0+XV0*SV0
IS=0
B1(1)=P0
B1(2)=WLV0
224
B1(3)=V0
B1(4)=T0
K=0
X=0
WRITE(8,*)' X PS WLV SA V MAX
&TS DC XV'
120 X=X+DX
A=X
K=K+1
CALL SEVAL(9,PS,P,T,BT,CT,DT,TS,P1TS,P2TS)
CALL SEVAL(9,P0,P,T,BT,CT,DT,T0,P1T0,P2T0)
CALL SEVAL(9,TS,T,VL,BVL,CVL,DVL,VLS,P1VLS,P2VLS)
CALL SEVAL(9,T0,T,VL,BVL,CVL,DVL,VL0,P1VL0,P2VL0)
CALL SEVAL(9,TS,T,VV,BVV,CVV,DVV,VVS,P1VVS,P2VVS)
CALL SEVAL(9,T0,T,VV,BVV,CVV,DVV,VV0,P1VV0,P2VV0)
CALL
SEVAL(9,TS,T,ENL,BENL,CENL,DENL,ENLS,P1ENLS,P2ENLS)
CALL
SEVAL(9,TS,T,ENV,BENV,CENV,DENV,ENVS,P1ENVS,P2ENVS)
CALL SEVAL(9,TS,T,SL,BSL,CSL,DSL,SLS,P1SLS,P2SLS)
CALL SEVAL(9,T0,T,SL,BSL,CSL,DSL,SL0,P1SL0,P2SL0)
CALL SEVAL(9,TS,T,SV,BSV,CSV,DSV,SVS,P1SVS,P2SVS)
CALL SEVAL(9,T0,T,SV,BSV,CSV,DSV,SV0,P1SV0,P2SV0)
CALL SEVAL(9,TS,T,R,BR,CR,DR,RS,P1RS,P2RS)
CALL
SEVAL(9,TS,T,CVL,BCVL,CCVL,DCVL,CVLS,P1CVLS,P2CVLS)
CALL
SEVAL(9,TS,T,CVP,BCVP,CCVP,DCVP,CVPS,P1CVPS,P2CVPS)
CALL
SEVAL(9,T0,T,CPW,BCPW,CCPW,DCPW,CP0W,P1CP0W,P2CP0W)
CALL
SEVAL(9,TS,T,CPW,BCPW,CCPW,DCPW,CPSW,P1CPSW,P2CPSW)
DO 110 LL=1,4
TSK=TS+273.
IF(X.LT.ZCC) DC=DC0-(DC0-DCM)*X/ZCC
IF(X.GT.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.
IF(X.LT.ZCC) PRFC=-PI*(DC0-(DC0-DCM)*X/ZCC)*(DC0-
DCM)/2./ZCC
IF(X.GE.ZCC.AND.X.LE.(ZCC+ZCM)) PRFC=0.0
225
IF(X.GT.(ZCC+ZCM)) PRFC=PI*(DCM+(DCSR-DCM)*(X-
(ZCC+ZCM))/ZCD)*
& (DCSR-DCM)/(2.*ZCD)
XV=(S0-SLS)/(SVS-SLS)
V=(1.-XV)*VLS+XV*VVS
CV=(1.-XV)*CVLS+XV*CVPS
SA=((V*RS /(VVS-VLS))**2/CV/TSK)**0.5
MAX=WLV/SA
C1(1)=-WLV**2*PRFC/(MAX**2-1.)/V/FC
C1(2)=WLV*PRFC/FC/(MAX**2-1.)
C1(3)=MAX**2*V*PRFC/FC/(MAX**2-1.)
C1(4)=-TSK*(VVS-VLS)*WLV**2*PRFC/(MAX**2-1.)/RS/V/FC
110 CALL RKM(A,DX,N,IS,B1,C1)
PS=B1(1)
WLV=B1(2)
V=B1(3)
TS=B1(4)
DC=(4.*FC/PI)**0.5
IF(K.EQ.KK) THEN
WRITE(8,200) X,PS,WLV,SA,V,MAX,TS,DC,XV
J=J+1
IF(J.EQ.10) THEN
WRITE(8,*)' X PS WLV SA V MAX
&TS DC XV'
J=0
END IF
K=0
END IF
IF(X.LT.Z) GOTO 120
200 FORMAT(F6.3,F10.0,4F8.3,F8.2,2F6.3,F10.0)
IF(K.EQ.KK)WRITE(8,*) ' X=',X,' PS=',PS,' WLV=',WLV,'
SA=',SA,
& ' GS=',GS,' V=',V,' MAX=',MAX,' TS=',TS,' DC=',DC,' XV=',XV
IF(K.EQ.KK) K=0.0
IF(X.LT.Z) GOTO 120
STOP
END
Файл исходных данных GM4.DAT
4 0.3E-3 0.15 30. 6.E5 0.1 1.E5 3.1416 10
3.E-2 18.72E-3 2.8685E-2 4.E-2 1.E-3 0.105
100. 110. 120. 130. 140. 150. 160. 170. 180.
226
1.01E5 1.43E5 1.98E5 2.70E5 3.61E5 4.76E5 6.18E5 7.92E5 10.03E5
1.04E-3 1.05E-3 1.06E-3 1.07E-3 1.08E-3 1.09E-3 1.1E-3 1.11E-3
1.13E-3
1.673 1.21 0.8917 0.6683 0.5087 0.3926 0.3068 0.2426 0.1939
419.E3 461.E3 503.E3 546.E3 589.E3 632.E3 675.E3 719.E3 763.E3
2676.E3 2691.E3 2706.E3 2721.E3 2734.E3 2746.E3 2758.E3 2769.E3
2778.E3
1.307E3 1.418E3 1.528E3 1.635E3 1.739E3 1.842E3 1.943E3 2.042E3
2.1395E3
7.355E3 7.239E3 7.13E3 7.027E3 6.93E3 6.838E3 6.751E3 6.667E3
6.5858E3
2257.E3 2230.E3 2202.E3 2175.E3 2145.E3 2114.E3 2082.E3 2050.E3
2015.E3
4.22E3 4.23E3 4.24E3 4.26E3 4.28E3 4.3E3 4.33E3 4.35E3 4.38E3
71.02E3 65.91E3 61.3E3 57.11E3 53.3E3 49.82E3 46.63E3 43.66E3
40.94E3
4.22E3 4.23E3 4.25E3 4.27E3 4.29E3 4.31E3 4.35E3 4.38E3 4.417E3
(T(I),I=1,9)
(P(I),I=1,9)
(VL(I),I=1,9)
(VV(I),I=1,9)
(ENL(I),I=1,9)
(ENV(I),I=1,9)
(SL(I),I=1,9)
(SV(I),I=1,9)
(R(I),I=1,9)
(CVL(I),I=1,9)
(CVP(I),I=1,9)
(CPL(I),I=1,9)
N, DX, Z, WLV0 ,P0, XV0, PSR, PI, KK
DC0, DCM, DCSR ZCC ZCM ZCD
PROGRAM GM1QT
C Расширение гомогенного потока с трением и
C подводом теплоты в сопле Лаваля - прямая задача
C ( задана геометрия сопла и начальная
C степень сухости -XV0).
DIMENSION B(10),C(10)
REAL MAX
VISL(A)=(613.-4.227*(A-273.)+9.22E-3*(A-273.)**2)*1.E-6
VISP(A)=(8.35+3.42E-2*(A-273.)+2.E-5*(A-273.)**2)*1.E-6
RO(A)=999.9-0.1475*(A-273.)-2.695E-3*(A-273.)**2
227
SQV(A)=(2447.7-1.283*(A-273.)-6.26E-3*(A-273.)**2)*1.E3
CVL(A)=4244.1-1.583*(A-273.)+1.314E-2*(A-273.)**2
CVV(A)=135950.-799.5*(A-273.)+1.502*(A-273.)**2
SL(A)=0.087725+0.0131875*(A-273.)-9.9375E-6*(A-273.)**2
SP(A)=8.764128-0.0165849*(A-273.)+2.49062E-5*(A-273.)**2
OPEN (2,FILE='GM1QT.DAT')
OPEN (8,FILE='REZGM1QT.DAT')
READ(2,*) N,DX,Z,WLV0,P0,XV0,PSR,PI,KK,
& DC0,DCM,DCSR,ZCC,ZCM,ZCD,QW
Z=ZCC+ZCM+ZCD
N=4
FC0=PI*DC0**2/4.
T0=5.677*P0**0.25+273.
ROV0=1.213E-5*P0**0.9375
ROL0=RO(T0)
VL0=1./ROL0
VV0=1./ROV0
V0=(1.-XV0)*VL0+XV0*VV0
S0=(1.-XV0)*SL(T0)+XV0*SP(T0)
FC=FC0
T=T0
ROL=ROL0
ROV=ROV0
WLV=WLV0
V=V0
P=P0
DC=DC0
IS=0
B(1)=P
B(2)=WLV
B(3)=V
B(4)=T
J=0
K=0
X=0
WRITE(8,*)' X P WLV SA V MAX T
& DC XV PV'
120 X=X+DX
A=X
K=K+1
DO 110 LL=1,4
228
ROV=1.213E-5*P**0.9375
IF(X.LT.ZCC) DC=DC0-(DC0-DCM)*X/ZCC
IF(X.GT.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.
IF(X.LT.ZCC) PRFC=-PI*(DC0-(DC0-DCM)*X/ZCC)*(DC0-
DCM)/2./ZCC
IF(X.GE.ZCC.AND.X.LE.(ZCC+ZCM)) PRFC=0.0
IF(X.GT.(ZCC+ZCM)) PRFC=PI*(DCM+(DCSR-DCM)*(X-
(ZCC+ZCM))/ZCD)*
& (DCSR-DCM)/(2.*ZCD)
ROL=RO(T)
VV=1./ROV
VL=1./ROL
SW=SL(T)
SV=SP(T)
XV=(S0-SW)/(SV-SW)
V=(1.-XV)*VL+XV*VV
ROSM=1./V
CV=(1.-XV)*CVL(T)+XV*CVV(T)
VISW=VISL(T)
VISV=VISP(T)
VISSM=(1.-XV)*VISW+XV*VISV
R=SQV(T)
SA=((V*R /(VV-VL))**2/CV/T)**0.5
MAX=WLV/SA
PV=P0/((P0/PSR-1.)*(A/Z)**2+1.)
DIPDV=R/(VV-VL)
DIPDP=CV*T*(VV-VL)/R+V
REN=WLV*DC*ROSM/VISSM
CF=0.079/(REN**0.25)
OPR=DIPDP-DIPDV*V**2/(WLV**2)-V
C(1)=(4.*QW*V/WLV/DC+2.*CF*(V*DIPDV+WLV**2)/DC-
& V*DIPDV*PRFC/FC)/OPR
C(2)=-V*C(1)/WLV-2.*CF*WLV/DC
C(3)=V*C(2)/WLV+V*PRFC/FC
C(4)=(VV-VL)*T*C(1)/R
110 CALL RKM(A,DX,N,IS,B,C)
P=B(1)
WLV=B(2)
V=B(3)
229
T=B(4)
IF(K.EQ.KK) THEN
WRITE(8,200) X,P,WLV,SA,V,MAX,T,DC,XV,PV
J=J+1
IF(J.EQ.10) THEN
WRITE(8,*)' X P WLV SA V MAX T
& DC XV PV'
J=0
END IF
K=0
END IF
IF(X.LT.Z) GOTO 120
200 FORMAT(F6.3,F10.0,4F8.3,F8.2,2F6.3,F10.0)
STOP
END
Файл исходных данных GM1Q.DAT
4 0.3E-3 0.15 30. 6.E5 0.1 1.E5 3.1416 10
3.E-2 18.72E-3 2.8685E-2 4.E-2 1.E-3 0.105 2.E6
N, DX, Z, WLV0 ,P0, XV0, PSR, PI, KK
DC0, DCM, DCSR ZCC ZCM ZCD QW
Идентификаторы программы GM1_CILS
DC – диаметр канала
FC – площадь поперечного сечения канала
CWV – скорость гомогенного потока
V – удельный объем смеси
TS, PS – текущие значения температуры и давления потока
PROGRAM GM1_CILS
C Расширение гомогенного
C потока с трением в цилиндрическом канале
C (определяется длина канала и критическое
C давление на срезе при заданном расходе,
C начальная степень
сухости-XV0 задается ).
C Используются подпрограммы: spline, seval,rkm.
DIMENSION B1(10),C1(10)
DIMENSION P(9),T(9),VL(9),VV(9),R(9),ENL(9),ENV(9),
& SL(9),SV(9),BSL(9),CSL(9),DSL(9),BSV(9),CSV(9),DSV(9),
& CPW(9),CVW(9),CVP(9),BT(9),CT(9),DT(9),BVL(9),
& CVL(9),DVL(9),BVV(9),CVV(9),DVV(9),BR(9),CR(9),
& DR(9),BENL(9),CENL(9),DENL(9),BENV(9),CENV(9),