Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение
Подождите немного. Документ загружается.
400
wall=3 :print#9, "выпала на стык" 'выпала на стык бандажа и
корыта
goto cont
end if
if wall=0 and u*1000>=ukk(r*1000,zz*1000) then
print #9," выпала на корыто"
wall=1
end if
cont:
IF zz*1000<=hor then
print#8,using"##.### "; ZZ*1000,(R-RV)*1000,U*1000
print#9, wall;" ";
print#9,using"+###.# "; ZZ*1000,(R-
RV)*1000,U*1000,ukk(r*1000,zz*1000),WDR,WDZ,WDU,WD
' print "zz>hor:";zz;">";hor
goto 120
else
goto ext
end if
ext:
close
STOP
function cxx(Re)
IF RE>=0. AND RE<=1. then CXX=24./RE
IF RE>=1. AND RE<=390. THEN CXX=24./(RE^0.667)
IF RE>390. THEN CXX=0.45
end function
function ukk(r,z)
t=r*2*pi/nl
bk=t*sin(b1) 'ширина канала
f1=(fi1+beta1)/180*pi
f2=(fi2+fi1+beta1)/180*pi
r3=(hor+bk*(cos(b1)+cos(f1)))/(cos(b1)+cos(b2))
r1=r3-bk
r4=r1
r2=r4-bk
'разбивка на зоны по z
z1=bk*cos(b1)
z2=r3*(cos(b1)-cos(f1))
z3=r3*cos(b1)-r1*cos(f1)
401
z4=hor-bk*cos(b2)
z5=hor
'координаты центров окружностей
z8=r3*cos(b1) : u8=t-r3*sin(b1)
z9=hor-r4*cos(b2) : u9=u8+(r3-r4)*sin(f1)
select case z
case 0 to z1
ukk=sqr(r3^2-(z-z8)^2)+u8
case z1 to z2
ukk=sqr(r3^2-(z-z8)^2)+u8
case z2 to z3
ukk=sqr(r4^2-(z-z9)^2)+u9
case z3 to z4
ukk=sqr(r4^2-(z-z9)^2)+u9
case z4 to hor
ukk=sqr(r4^2-(z-z9)^2)+u9
end select
end function
sub wvap(r,u,z) 'cкорость пара
t=r*2*pi/nl
bk=t*sin(b1) 'ширина канала
f1=(fi1+beta1)/180*pi
f2=(fi2+fi1+beta1)/180*pi
r3=(hor+bk*(cos(b1)+cos(f1)))/(cos(b1)+cos(b2))
r1=r3-bk
r4=r1
r2=r4-bk
'разбивка на зоны по z
z1=bk*cos(b1)
z2=r3*(cos(b1)-cos(f1))
z3=r3*cos(b1)-r1*cos(f1)
z4=hor-bk*cos(b2)
z5=hor
'координаты центров окружностей
z8=r3*cos(b1) : u8=t-r3*sin(b1)
z9=hor-r4*cos(b2) : u9=u8+(r3-r4)*sin(f1)
'уравнения лучей b,c,d
'определение зоны (1-4)
select case z
case 0 to z1
402
ub=t-z*atn(b1)
if u<=ub then zona=1
if u>ub then zona=2
case z1+dz to z2
zona=2
case z2+dz to z3
uc=u9+r4*sin(f1)-(z-z2)*atn(f1)
if u<=uc then zona=2
if u>uc then zona=3
case z3+dz to z4
zona=3
case z4+dz to hor
ud=u9+r2*sin(b2)+(z-z4)*atn(b2)
if u>=ud then zona=3
if u<ud then zona=4
end select
select case zona
case 1:WVR=0: WVU=WV1*cos(b1) : WVZ=WV1*sin(b1)
rkor=1e5
fi=b1
case 2
a0=WV1*(r3-r1)/log(r3/r1)
rt=sqr((z-z8)^2+(u-u8)^2): WVM=a0/rt
fi=atn((u-u8)/(z8-z))
WVR=0
wvu=wvm*cos(fi) : wvz=wvm*sin(fi)
u0=u8: z0=z8: rkor=r3
case 3
a0=WV1*(r4-r2)/log(r4/r2)
rt=sqr((z-z9)^2+(u-u9)^2): WVM=a0/rt
WVR=0
select case z
case <z9: fi=atn((u-u9)/(z9-z))
case z9 : fi=1.57
case >z9: fi=1.57+atn((z-z9)/(u-u9))
end select
wvu=wvm*cos(fi) : wvz=wvm*sin(fi)
u0=u9: z0=z9 : rkor=r4
case 4: WVR=0:
WVU=-WV1*cos(b2)
WVZ=WV1*sin(b2)
rkor=1e5
403
fi=pi-b2
end select
print#10,zona,wvz,wvu,sqr(wvz^2+wvu^2),fi/3.14*180
end sub
FUNCTION RO(H)
RO=999.9-0.1475*(H-273.)-2.695E-3*(H-273.)^2 '+
END FUNCTION
function ROP(H)
ROP=6.7962-0.129299*(H-273.)+6.7314e-4*(H-273.)^2
end function
function ps(H)
PS=1.E6*EXP(82.86568+0.01028003*H-7821.541/H-
11.48776*LOG(H))
end function
function visp(H)
VISP=(8.35+3.42E-2*(H-273.)+2.E-5*(H-273.)^2)*1.E-6
end function
FUNCTION VI(H)
IV=(2461.1+2.636*(H-273.)-4.88E-3*(H-273.)^2)*1.E3 '+
END FUNCTION
FUNCTION VISL(H)
VISL=(613.-4.227*(H-273.)+9.22E-3*(H-273.)^2)*1.E-6 '+
END FUNCTION
FUNCTION YIL(H)
YIL=(12.475+3.935*(H-273.)+1.3125E-3*(H-273.)^2)*1.E3
END FUNCTION
FUNCTION CP(H)
CP=4331.-3.09*(H-273.)+1.98E-2*(H-273.)^2
END FUNCTION
FUNCTION CPL(H)
CPL=4275.5-1.2209*(H-273.)+0.59823E-2*(H-273.)^2
END FUNCTION
FUNCTION CPLS(H)
CPLS=(CP(TW0)+CP(H))/2.
END FUNCTION
FUNCTION CPDS(H)
CPDS=(CP(TD0)+CP(H))/2.
END FUNCTION
FUNCTION YIV(H)
YIV=(2461.1+2.636*(H-273.)-4.88E-3*(H-273.)^2)*1.E3
END FUNCTION
FUNCTION SIG(H)
404
SIG=(75.54-0.168*(H-273.))*1.E-3
END FUNCTION
FUNCTION VIS(H)
VIS=(613.-4.227*(H-273.)+9.22E-3*(H-273.)^2)*1.E-6
END FUNCTION
FUNCTION SQV(H)
SQV=(2447.7-1.283*(H-273.)-6.26E-3*(H-273.)^2)*1.E3
END FUNCTION
FUNCTION PP(H)
PP=1.E6*EXP(82.86568+0.01028003*H-7821.541/H-
11.48776*LOG(H))
END FUNCTION
FUNCTION YIW(H)
YIW=(12.475+3.935*(H-273.)+1.3125E-3*(H-273.)^2)*1.E3
END FUNCTION
FUNCTION PR(H)
PR=1.28-4.4E-3*(H-273.)+2.4E-5*(H-273.)^2
END FUNCTION
FUNCTION LAM(H)
LAM=(1.801+2.69E-3*(H-273.)+3.02E-5*(H-273.)^2)*1.E-2
END FUNCTION
'**********************************************************
SUB RKM(A,DX,N,Ik,B(),C())
' DIM F(15),GG(60)
DIM GG(60)
Ik=Ik+1
ON Ik GOTO 10,30,60,80
10 for I=1 to N
F(I)=B(I)
GG(4*I-3)=C(I)*DX
B(I)=F(I)+GG(4*I-3)/2.
next
A=A+DX/2
GOTO 100
30 for I=1 to N
GG(4*I-2)=C(I)*DX
B(I)=F(I)+GG(4*I-2)/2.
next
GOTO 100
60 for I=1 to N
GG(4*I-1)=C(I)*DX
B(I)=F(I)+GG(4*I-1)
405
next
A=A+DX/2
GOTO 100
80 for I=1 to N
GG(4*I)=C(I)*DX
B(I)=GG(4*I-3)+2.*(GG(4*I-2)+GG(4*I-1))
B(I)=(B(I)+GG(4*I))/6.+F(I)
next
ik=0
100 abb=0
end sub 'RKM
'**********************************************************
Файл данных DROP_3D.DAT
65, " число лопаток "
120, " D1 - внутр. радиус. колеса (мм)"
160, " D2 - наружн. радиус. колеса (мм)"
1.7, " bt - отношение хорда/шаг "
32, " beta1 - угол входа (град)"
32, " beta2 - угол выхода (град)"
42, " fi1 - угол смены радиусов лопатки"
1, " dd - шаг при вычерчивании лопатки (мм)"
0.5E-4 411.57 415.53 72 160 3800. 8314. 18. 0.5 0.167 2.E-5
DZ TV TD WD CV1 OB UGC AMV A1 Gkap DD
PROGRAM CON_STREAM
C Конденсация паровой струи в спутном потоке
C холодной воды.
DIMENSION B(10),C(10)
REAL IV,ILTL1,ID,ISTU,LAMV,ISTV1,ISTV,ISTL,MAX,KN
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='CON_ST.DAT')
OPEN (6,FILE='CON')
OPEN (8,FILE='REZCON_ST.DAT')
READ(2,*) N,DX,GV1,PV1,GL1,TL1,PL1,DNL1,
& DN1,ZNC,UGC,AMV,PI,KK
Z=ZNC
406
FN1=PI*DN1**2/4.
FNL1=PI*DNL1**2/4.-PI*DN1**2/4.
ROL1=RO(TL1)
CL1=GL1/ROL1/FNL1
ILTL1=YIL(TL1)*GL1
TV1=5.677*PV1**0.25+273.
ROV1=1.213E-5*PV1**0.9375
CV1=GV1/ROV1/FN1
FV1=GV1/ROV1/CV1
ISTV1=GV1*(YIV(TV1)+CV1**2/2.)
YISTV1=(YIV(TV1)+CV1**2/2.)
PV=PV1
CV=CV1
GV=GV1
FN=FN1
DN=DN1
TV=TV1
TVR=TV1
CL=CL1
TL=TL1
PL=PL1
ROL=RO(TL1)
VISV1=(8.35+3.42E-2*(TV1-273.)+2E-5*(TV1-273.)**2)*1E-6
REV=CV1*DN1*ROV1/VISV1
CP=CPL(TL1)
KN=SQV(TV1)/CP/(TV1-TL1)
CPV=2.E3
PRV=1.28-4.4E-3*(TV1-273.)+2.4E-5*(TV1-273.)**2
H=(0.352E-*REV**0.4*KN**0.08*PRV*CP/CPV)**1.67*DN1
COSB=1./((DN1/2./H)**2+1.)**0.5
COSG=COS(PI/2.-COSB)
DN=DN1*(H-X)/H
FN=PI*DN**2/4.
GV=CV*ROV*FN
ALFA=2.*GV1*SQV(TV1)*COSB/PI/DN1/H/(TV1-TL1)
PRGV=ALFA*(TV-TL)*PI*DN/SQV(TV)/COSB
DLGV=GV1-GV
PRFN=-PI*(DN1-DN1*X/H)*DN1/H/2.
PRINT*, ' PV=',PV,' CV=',CV, ' FN=',FN,' GV=',GV,' DN=',DN,
& ' ALFA=',ALFA,' TL1=',TL1,' PL1=',PL1,' TV=',TV,' H=',H
WRITE(6,*) ' X=',X,' PV=',PV,' CV=',CV, ' FN=',FN,' GV=',GV,
& ' DN=',DN,
407
& ' ROL=',ROL,' TL1=',TL1,' PL1=',PL1,' TV=',TV
N=6
IS=0
B(1)=PV
B(2)=CV
B(3)=TVR
B(4)=CL
B(5)=TL
B(6)=PL
K=0
X=0
140 X=X+DX
K=K+1
IF(X.GT.Z)GOTO 16
DO 150 LL=1,4
ROV=1.213E-5*PV**.9375
TV=5.677*PV**.25+273.
SQ=SQV(TV)
VISV1=(8.35+3.42E-2*(TV1-273.)+2E-5*(TV1-273.)**2)*1E-6
REV=CV1*DN1*ROV1/VISV1
CP=CPL(TL1)
KN=SQV(TV1)/CP/(TV1-TL1)
PRV=1.28-4.4E-3*(TV1-273.)+2.4E-5*(TV1-273.)**2
H=(0.352E-2*REV**0.4*KN**0.08*PRV*CP/CPV)**1.67*DN1
COSB=1./((DN1/2./H)**2+1.)**0.5
COSG=COS(PI/2.-COSB)
DN=DN1*(H-X)/H
FN=PI*DN**2/4.
FNL=PI*DNL1**2/4.-PI*DN**2/4.
GV=CV*ROV*FN
ALFA=2.*GV1*SQV(TV1)*COSB/PI/DN1/H/(TV1-TL1)
PRGV=ALFA*(TV-TL)*PI*DN/SQV(TV)/COSB
DLGV=GV1-GV
GL=GL1+DLGV
PRFN=-PI*(DN1-DN1*X/H)*DN1/H/2.
PRFNL=PI*(DN1-DN1*X/H)*DN1/H/2.
J=ALFA*(TV-TL)/SQ
LAMV=(1.801+2.69E-3*(TV-273.)+3.02E-5*(TV-273.)**2)*1.E-2
PPROV=8.805129E4*PV**62.5E-3
GCV=UGC/AMV
FNL=PI*(DNL1**2-DN**2)/4.
PVTL=((TL-273.)/5.667)**4
408
IV=YIV(TV)
ISTV=(IV+CV**2/2.)
U=ALFA*(TV-TL)/ROV/SQ
ISTU=(IV+U**2/2.)
CPV=2.E3
ADK=1.135
MAX=CV/(ADK*TV*GCV)**0.5
ISTL=YIL(TL)+CL**2/2.
C(1)=ROV*CV*(((U*COSG+CV)/GV+1./ROV/FN)*PRGV-
CV*PRFN/FN)/
& (MAX**2-1.)
C(2)=PRGV/ROV/FN-CV*PRFN/FN-MAX**2*C(1)/ROV/CV
C C(2)=-(U*COSG+CV)*PRGV/GV-C(1)/ROV/CV
C(3)=-((ISTV-ISTU)*PRGV/GV+CV*C(2))/CPV
C(4)=((CV-CL)*PRGV-FNL*C(6))/GL
C(5)=((ISTV-ISTL)*PRGV/GL-CL*C(4))/CP
C(6)=-GL*(J*PI*DN/ROL/COSB-CL*PRFNL)/FNL**2+(CV-
CL)*PRGV/FNL
CALL RKM(X,DX,N,IS,B,C)
PV=B(1)
CV=B(2)
TVR=B(3)
CL=B(4)
TL=B(5)
150 PL=B(6)
DLP=PL-PV
WRITE(6,*) ' X=',X,' PV=',PV,' CV=',CV,' GV=',GV,' MAX=',MAX
IF(K.EQ.KK)WRITE(8,*) ' X=',X, ' PV=',PV,' CV=',CV, ' FN=',FN,
& ' GV=',GV,' DN=',DN,' ROL=',ROL,' TL1=',TL1,' PL1=',PL1,'
TV=',
& TV,' IV=',IV,' U=',U,' TVR=',TVR,' H=',H,' ALFA=',ALFA,
& ' MAX=',MAX,' PRGV=',PRGV,' COSG=',COSG,'
COSB=',COSB,
& ' KN=',KN,' CL=',CL,' GL=',GL,' FNL=',FNL,' TL=',TL,' PL=',PL,
& ' DLP=',DLP
IF(K.EQ.KK)K=0
IF(MAX.LT.1.1)GOTO 140
16 STOP
END
SUBROUTINE RKM(A,DX,N,IS,B,C)
DIMENSION F(10),G(40),C(10),B(10)
IS=IS+1
409
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
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
Файл данных CON_ST.DAT
6 0.01E-5 0.005 1.E5 0.1 283. 1.E5 3.E-2 10.E-3 5.8E-3 8314. 18.
3.1416 10
N, DX, GV1, PV1, GL1,TL1, PL1, DN DN1, ZNC, UGC, AMV,
PI, KK
PROGRAM ACCEL 1
' Расчет смесительного устройства,
' состоящего из сопла Лаваля и камеры смешения
$STACK 4000
SCREEN 9 :OPTION BASE 1
VIEW:WINDOW (0,0)-(250,170)
' PROGRAM LAV27
' расчет сопла с парогенерирующей решеткой на входе
SHARED B(),C()
OPEN "I",#1,"L27.DAT"
OPEN "O",#2,"REZL27.DAT"
input #1,N,DX,GCR,GV1,TL0,PL0,TD1,PV1