Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение
Подождите немного. Документ загружается.
380
C(2)=(Y(2)-Y(1))/D(1)
DO 10 I=2,NM1
D(I)=X(I+1)-X(I)
B(I)=2.0*(D(I-1)+D(I))
C(I+1)=(Y(I+1)-Y(I))/D(I)
C(I)=C(I+1)-C(I)
10 CONTINUE
B(1)=-D(1)
B(N)=-D(N-1)
C(1)=0.0
C(N)=0.0
IF(N.EQ.3) GO TO 15
C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
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
XXVI
XXVI
381
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)
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
Файл данных TURB8.DAT
11.97E-6 297.93E-6 42. 73.64 20. 3. 65.
40.E-3 18. 32.11 .141 1. 1. 2. 0.05E-3
8314. 18. 70. 11.6E-3 18.4E-3 4.8E-3 11.6E-3 52
1.73 3.8E3
1.1E5 9.5E4 8.5E4 7.5E4 6.5E4
380.13 375.32 371.20 368.14 364.78
1.5501 1.7779 1.9728 2.2179 2.5357
.0010495 .0010455 .0010423 .0010400 .0010375
449.19E3 428.84E3 411.49E3 398.63E3 384.45E3
2681.4E3 2680.E3 2673.5E3 2668.6E3 2663.2E3
4.218E3 4.217E3 4.216E3 4.215E3 4.214E3
375.32 371.20 368.14 364.78 361.02
1.1E5 9.5E4 8.5E4 7.5E4 6.5E4
2680.E3 2673.5E3 2668.6E3 2663.2E3 2657.E3
5E8
VISP VISK FI1 FI2 BL RBL PROP
HB ALF ZERO RM FIS DELF A1 A2
UGC AMV B1 R1 R2 R3 R4 NK
TG VOR
{SPLINE} BAR(I) TEMP(I) VPA(I) VKAP(I) ENT(I)
XXVII
382
ENTS(I) CP(I) TPA(I) DAVL(I) ENTN(I)
EPS
PROGRAM DROPLET
EXTERNAL DROP
REAL Y(5),WORK(33),CKR(6)
INTEGER IWORK(5)
COMMON CVR,CVU,CVUZ,DD,VISV,VV,VD
DATA PI/3.141592653589/,PF/180./
OPEN(2,FILE='MD.DAT')
OPEN(3,FILE='MD1.DAT')
OPEN(4,FILE='MDREZ.DAT')
OPEN(5,FILE='TUR_DR.DAT')
READ(2,*) RST,RFIN,B2G
CLOSE(2)
READ(3,*) DD,RPRINT,KI
CLOSE(3)
READ(5,*) VOR,WPL,WV,VISV,VV,VD
CLOSE(5)
WRITE(4,*)' VISV=',VISV,' VV=',VV,' VD=',VD
OM=PI*VOR/30.
WRITE(4,*)' OM=',OM,' WPL=',WPL,' DD=',DD,' RFIN=',RFIN,
& ' RPRINT=',RPRINT
B2=B2G*PI/PF
CKZ=WPL*SIN(B2)
CVR=0.
WVZ=WV*SIN(B2)
WVU=WV*COS(B2)
WRITE(4,*)' WV=',WV,' WVZ=',WVZ,' WVU=',WVU
DO 3 I=1,6
CKR(I)=.1+2*(I-1)
WRITE(4,*)' Радиальная скорость=',CKR(I)
NEQN=5
R=RST
RELERR=.6E-5
ABSERR=0.
IFLAG=1
ROUT=R
U=R*OM
CVU=WVU-U
CVUZ=WVZ/SIN(ATAN(WVZ/CVU))
NK=0
383
CKU=U-WPL*COS(B2)
CKUZ=WPL*SIN(B2)/SIN(ATAN(CKZ/CKU))
Y(1)=CKR(I)
Y(2)=CKU
Y(3)=CKUZ
Y(4)=0.
Y(5)=0.
1 CALL
RKF45(DROP,NEQN,Y,R,ROUT,RELERR,ABSERR,IFLAG,WORK,I
WORK)
CKR(I)=Y(1)
CKU=Y(2)
CKUZ=Y(3)
C=SQRT(CKR(I)**2+CKU**2+CKZ**2)
ALFA=CKZ/C*180/PI
BETA=CKU/C*180/PI
GAMMA=CKR(I)/C*180/PI
WRITE(4,1000)C,ALFA,BETA,GAMMA
1000 FORMAT(' C,ALFA,BETA,GAMMA - ',F6.2,3F7.1)
WRITE(4,1001)CKR(I),CKU,CKZ
1001 FORMAT('CKR,CKU,CKZ - ',3F6.2)
U=R*OM
CVU=WVU-U
CVUZ=WVZ/SIN(ATAN(WVZ/CVU))
PRINT*, I,R
WRITE(4,*)' U=',U,' CVU=',CVU,' CVUZ=',CVUZ
GOTO (70,2,20,30,40,50,60,70),IFLAG
2 ROUT=R+RPRINT
NK=NK+1
IF(NK.EQ.KI) WRITE(4,*) NK,R,Y
IF(NK.EQ.KI) NK=0
IF(R.LT.RFIN) GOTO 1
GOTO 3
20 PRINT *,RELERR,ABSERR
GOTO 1
30 PRINT*,' МНОГО ШАГОВ'
GOTO 1
40 ABSERR=1E-7
PRINT *,RELERR,ABSERR
GOTO 1
50 RELERR=10.*RELERR
PRINT *,RELERR,ABSERR
384
IFLAG=2
GOTO 1
60 PRINT*,' МНОГО ДАННЫХ'
IFLAG=2
GOTO 1
70 PRINT *,'НЕВЕРНЫЙ ВЫЗОВ'
STOP
3 CONTINUE
CLOSE(4)
STOP
END
SUBROUTINE DROP(R,Y,YP)
REAL Y(5),YP(5)
COMMON CVR,CVU,CVUZ,DD,VISV,VV,VD
R0(A,B)=A/B-1.
RCX=CX(Y(1),Y(2),Y(3))
YP(1)=Y(2)**2/(Y(1)*R)-RCX*R0(CVR,Y(1))
YP(2)=-Y(2)/R+RCX*R0(CVU,Y(2))
YP(3)=RCX*R0(CVR,Y(3))
YP(4)=Y(2)/(Y(1)*R)
Y(5)=Y(3)/Y(1)
RETURN
END
FUNCTION CX(RCD,UCD,ZCD)
COMMON CVR,CVU,CVUZ,DD,VISV,VV,VD
DC=SQRT((CVR-RCD)**2+(CVU-UCD)**2+(CVUZ-ZCD)**2)
RED=DC*DD/(VISV*VV)
CX=0.45
IF(RED.GE.1E-3.AND.RED.LE.6E3)
CX=24./RED+4.4/SQRT(RED)+.32
CX=.75*CX*VD/VV*DC/DD
RETURN
END
SUBROUTINE
RKF45(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,
& IWORK)
EXTERNAL F
INTEGER IWORK(5)
REAL Y(NEQN),WORK(1)
K1M=NEQN+1
K1=K1M+1
K2=K1+NEQN
385
K3=K2+NEQN
K4=K3+NEQN
K5=K4+NEQN
K6=K5+NEQN
K7=K6+1
CALL
RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK(1),
&
WORK(K1M),WORK(K1),WORK(K2),WORK(K3),WORK(K4),WOR
K(K5),
&
WORK(K6),WORK(K7),IWORK(1),IWORK(2),IWORK(3),IWORK(4),
& IWORK(5))
RETURN
END
SUBROUTINE
RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,YP,H,F1,
& F2,F3,F4,F5,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,KFLAG)
LOGICAL HFAILD,OUTPUT
REAL
Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),
& F5(NEQN)
EXTERNAL F
DATA REMIN,MAXNFE/1E-12,3000/
IF(NEQN.LT.1)GOTO 10
IF((RELERR.LT.0.).OR.(ABSERR.LT.0.))GOTO 10
MFLAG=IABS(IFLAG)
IF(MFLAG.EQ.0.OR.MFLAG.GT.8)GOTO 10
IF(MFLAG.NE.1)GOTO 20
EPS=1.
5 EPS=EPS/2.
EPSP1=EPS+1
IF(EPSP1.GT.1.)GOTO 5
U26=26*EPS
GOTO 50
10 IFGLAG=8
RETURN
20 IF(T.EQ.TOUT.AND.KFLAG.NE.3)GOTO 10
IF(MFLAG.NE.2) GOTO 25
IF(KFLAG.EQ.3.OR.INIT.EQ.0)GOTO 45
IF(KFLAG.EQ.4)GOTO 40
IF(KFLAG.EQ.5.AND.ABSERR.EQ.0.)GOTO 30
386
IF((KFLAG.EQ.6).AND.(RELERR.LE.SAVRE).AND.(ABSERR.LE.SA
VAE))
& GOTO 30
GOTO 50
25 IF(IFLAG.EQ.3)GOTO 45
IF(IFLAG.EQ.4)GOTO 40
IF(IFLAG.EQ.5.AND.ABSERR.GT.0.)GOTO 45
30 STOP
40 NFE=0
IF(MFLAG.EQ.2)GOTO 50
45 IFLAG=JFLAG
IF(KFLAG.EQ.3)MFLAG=IABS(IFLAG)
50 JFLAG=IFLAG
KFLAG=0
SAVRE=RELERR
SAVAE=ABSERR
RER=2.*EPS+REMIN
IF(RELERR.GE.RER)GOTO 55
RELERR=RER
IFLAG=3
KFLAG=3
RETURN
55 DT=TOUT-T
IF(MFLAG.EQ.1)GOTO 60
IF(INIT.EQ.0)GOTO 65
GOTO 80
60 INIT=0
KOP=0
A=T
CALL F(A,Y,YP)
NFE=1
IF(T.NE.TOUT)GOTO 65
IFLAG=2
RETURN
65 INIT=1
H=ABS(DT)
TOLN=0.
DO 70 K=1,NEQN
TOL=RELERR*ABS(Y(K))+ABSERR
IF(TOL.LE.0.)GOTO 70
TOLN=TOL
387
YPK=ABS(YP(K))
IF(YPK*H**5.GT.TOL) H=(TOL/YPK)**.2
70 CONTINUE
IF(TOLN.LE.0.)H=0.
H=AMAX1(H,U26*AMAX1(ABS(T),ABS(DT)))
JFLAG=ISIGN(2,IFLAG)
80 H=SIGN(H,DT)
IF(ABS(H).GE.2.*ABS(DT))KOP=KOP+1
IF(KOP.NE.100)GOTO 85
KOP=0
IFLAG=7
RETURN
85 IF(ABS(DT).GT.U26*ABS(T))GOTO 95
DO 90 KK=1,NEQN
90 Y(KK)=Y(KK)+DT*YP(KK)
A=TOUT
CALL F(A,Y,YP)
NFE=NFE+1
GOTO 300
95 OUTPUT=.FALSE.
SCALE=2./RELERR
AE=SCALE*ABSERR
100 HFAILD=.FALSE.
HMIN=U26*ABS(T)
DT=TOUT-T
IF(ABS(DT).GE.2.*ABS(H))GOTO 200
IF(ABS(DT).GT.ABS(H))GOTO 150
OUTPUT=.TRUE.
H=DT
GOTO 200
150 H=.5*DT
200 IF(NFE.LE.MAXNFE)GOTO 220
IFLAG=4
KFLAG=4
RETURN
220 CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1)
NFE=NFE+5
EEOET=0.
DO 250 K=1,NEQN
ET=ABS(Y(K))+ABS(F1(K))+AE
IF(ET.GT.0.)GOTO 240
IFLAG=5
388
RETURN
240 EE=ABS((-2090.*YP(K)+(21970.*F3(K)-15048.*F4(K)))+
& (22528.*F2(K)-27360.*F5(K)))
250 EEOET=AMAX1(EEOET,EE/ET)
ESTTOL=ABS(H)*EEOET*SCALE/752400.
IF(ESTTOL.LE.1.)GOTO 260
HFAILD=.TRUE.
OUTPUT=.FALSE.
S=.1
IF(ESTTOL.LT.59049.)S=.9/ESTTOL**.2
H=S*H
IF(ABS(H).GT.HMIN)GOTO 200
IFLAG=6
KFLAG=6
RETURN
260 T=T+H
DO 270 K=1,NEQN
270 Y(K)=F1(K)
A=T
CALL F(A,Y,YP)
NFE=NFE+1
S=5.
IF(ESTTOL.GT.1.889568E-4)S=.9/ESTTOL**.2
IF(HFAILD) S=AMIN1(S,1.)
H=SIGN(AMAX1(S*ABS(H),HMIN),H)
IF(OUTPUT)GOTO 300
IF(IFLAG.GT.0)GOTO 100
IFLAG=-2
RETURN
300 T=TOUT
IFLAG=2
RETURN
END
SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S)
REAL
Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),
& F5(NEQN),S(NEQN)
CH=H/4.
DO 221 K=1,NEQN
221 F5(K)=Y(K)+CH*YP(K)
CALL F(T+CH,F5,F1)
CH=3.*H/32.
389
DO 222 K=1,NEQN
222 F5(K)=Y(K)+CH*(YP(K)+3.*F1(K))
CALL F(T+3.*H/8.,F5,F2)
CH=H/2197.
DO 223 K=1,NEQN
223 F5(K)=Y(K)+CH*(1932.*YP(K)+(7296.*F2(K)-7200.*F1(K)))
CALL F(T+12.*H/13.,F5,F3)
CH=H/4104.
DO 224 K=1,NEQN
224 F5(K)=Y(K)+CH*((8341.*YP(K)-845.0*F3(K))+(29440.*F2(K)-
32832.
& *F1(K)))
CALL F(T+H,F5,F4)
CH=H/20520.
DO 225 K=1,NEQN
225 F1(K)=Y(K)+CH*((-6080.*YP(K)+(9295.*F3(K)-
5643.*F4(K)))+(41040.
& *F1(K)-28352.*F2(K)))
CALL F(T+H/2.,F1,F5)
CH=H/7618050.
DO 230 K=1,NEQN
230 S(K)=Y(K)+CH*((902880.*YP(K)+(3855735.*F3(K)-
1371249.*F4(K)))
& +(3953664.*F2(K)+277020.*F5(K)))
RETURN
END
Файл данных MD.DAT
0.12 0.16 32.
RST RFIN B2G
Файл данных MD1.DAT
2.E-5 6.666E-3 1
DD RPRINT KI
Файл данных TUR_DR.DAT
3800.0 33.628430 112.170300 1.197000E-05
1.741103 1.037506E-03
PROGRAM RK1_GPT
C Программа позволяет выполнить расчет
C движения жидкостного элемента по перу рабочей