|
Страница 7 из 8
Приложение 1. Текст программы расчета ДУ и графики результатов расчета.
INTEGER*2 npoint/10000/,ncurv/5/,k3,i
REAL,ALLOCATABLE:: XARR(:),YARR(:,:)
REAL,ALLOCATABLE:: XARR1(:),YARR1(:,:)
REAL delenx/1.3/,deleny/1.3/,alfa/0/,alfa2/1/,w/314/,fi/0.0/,t,AMax
LOGICAL*1 poligrf/.FALSE./
DIMENSION Y(4),DY(4)
DIMENSION XOD(10),SILA(10)
DIMENSION XD(10),SSLA(10),TEMP(20),TEMP2(20),CP(20),RO(20),TEMP3(20),zh(20)
COMMON P0,P,S,S1,V,AM,AL,U,SS,ALX,alfa,alfa2,XOD,SILA,XD,SSLA,alx1,ALK,TEMP,CP,TEMP2,RO,CPVUX,ROVUX,TEMP3,zh,h
common /comA/ w,AMax NAMELIST/DATA/S,S1,V,P0,AM,AK2,AL,R,AMax,U,alx1,ALX,ALK,Y,XOD,SILA,XD,SSLA
NAMELIST/DATA2/TEMP,CP,TEMP2,RO,TEMP3,zh
OPEN(1,FILE='aa52.inp')
READ(1,NML=DATA)
READ(1,NML=DATA2)
WRITE(*,NML=DATA)
WRITE(*,NML=DATA2)
ALLOCATE ( XARR(npoint),YARR(NCURV,npoint),STAT=I)
ALLOCATE ( XARR1(npoint),YARR1(NCURV,npoint),STAT=I)
IF (I.NE.0) STOP'error'
K3=0
OPEN(3,FILE='results.txt')
X=0.
99 CONTINUE
Ht=0.0001
K3=K3+1
CALL RKYT(Y,Ht,X,DY)
WRITE(*,*) 'X=',X,'Y=',Y
XARR(K3)=X
YARR(1,K3)=-Y(1)
YARR(2,K3)=Y(2)
YARR(3,K3)=Y(3)
YARR(4,K3)=Y(4)
YARR(5,K3)=-SS
XARR1(K3)=Y(4)
YARR1(1,K3)=Y(3)/(Y(4)*h*ROVUX)
YARR1(2,K3)=CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX))
YARR1(3,K3)=0
YARR1(4,K3)=0
YARR1(5,K3)=0
WRITE(3,15)X,Y(1),Y(2),Y(3),Y(4),P,SS,CPVUX,ROVUX
15 FORMAT(2X,G12.4,8(2X,G8.3))
IF(Y(1)>=ALK)then
CALL EGRAFIC (K3,XARR,NCURV,YARR,DELENX,DELENY,poligrf)
CALL EGRAFIC (K3,XARR1,5,YARR1,DELENX,DELENY,poligrf)
STOP
endif
goto 99
END
SUBROUTINE FN(Y,X,DY)
DIMENSION Y(4),DY(4)
DIMENSION XOD(10),SILA(10)
DIMENSION XD(10),SSLA(10),TEMP(20),TEMP2(20),CP(20),RO(20),TEMP3(20),zh(20)
COMMON P0,P,S,S1,V,AM,AL,U,SS,ALX,alfa,alfa2,XOD,SILA,XD,SSLA,alx1,ALK,TEMP,CP,TEMP2,RO,CPVUX,ROVUX,TEMP3,zh,h
Z=0.
CALL LINAP(10,Y(1),XOD,SILA,P)
CALL LINAP(10,Y(1),XD,SSLA,SS)
CALL LINAP(20,Y(4),TEMP,CP,CPVUX)
CALL LINAP(20,Y(4),TEMP2,RO,ROVUX)
CALL LINAP(20,Y(4),TEMP3,zh,h)
IF(P0/Y(3)>=0.59)then
Z=SQRT(abs(2.*(CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))/(((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))-1)*(Y(3)/(Y(4)*h*ROVUX))*Y(4))*((P0/Y(3))**(2./(CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX))))-(P0/Y(3))**(((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))+1)/(CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))))))
RAS=SS*Y(3)*Z
else
Z=SQRT(abs(2.*(CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))/(((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))+1)*(Y(3)/(Y(4)*h*ROVUX))*Y(4))*(2./((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))+1))**(2./((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))-1))))
RAS=SS*Y(3)*Z
endif
IF(Y(1)<ALX)then
B1=0.
else
B1=1.
endif
DY(1)=Y(2)
DY(2)=(P-S*(Y(3)-P0))/AM
DY(3)=((CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX)))/(S*(AL-Y(1))))*((1-1/(CPVUX/(CPVUX-Y(3)/(Y(4)*h*ROVUX))))*AD(Y(1))*(B1*B2*(5000.*(Y(1)-ALX)/1.0))*abs(AI(X))+Y(3)*S*DY(1)-(Y(3)/(Y(4)*h*ROVUX))*RAS*B1*Y(4))
DY(4)=(-DY(1)*Y(4)/(AL-Y(1))+DY(3)*Y(4)/Y(3)+(Y(3)/(Y(4)*h*ROVUX))*Y(4)*Y(4)*B1*RAS/(Y(3)*S*(AL-Y(1))))
RETURN
END
SUBROUTINE RKYT(Y,H,X,DY)
DIMENSION Y(4),DY(4),AA(4),S(4),D(4)
AA(1)=Y(1)
AA(2)=Y(2)
AA(3)=Y(3)
AA(4)=Y(4)
CALL FN(Y,X,DY)
DO 3 I=1,4
S(I)=H*DY(I)
D(I)=S(I)
3 Y(I)=AA(I)+S(I)/2.
X=X+H/2.
CALL FN(Y,X,DY)
DO 4 I=1,4
S(I)=H*DY(I)
D(I)=D(I)+2.*S(I)
4 Y(I)=AA(I)+S(I)/2.
CALL FN(Y,X,DY)
DO 5 I=1,4
S(I)=H*DY(I)
D(I)=D(I)+2.*S(I)
5 Y(I)=AA(I)+S(I)
X=X+H/2.
CALL FN(Y,X,DY)
DO 6 I=1,4
Y(I)=AA(I)+(D(I)+H*DY(I))/6.
6 AA(I)=Y(I)
RETURN
END
SUBROUTINE LINAP(N,U,X,Y,ZN)
DIMENSION X(N),Y(N)
I=N-1
IF(U.GE.X(N)) GO TO 20
I=1
IF(U.LE.X(1)) GO TO 20
J=N+1
10 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 10
20 CONTINUE
UU=ABS(U)
ZNAKU=SIGN(1.0,U)
DX=UU-X(I)
ZN=Y(I)+DX*(Y(I+1)-Y(I))/(X(I+1)-X(I))
ZN=SIGN(ZN,ZNAKU)
RETURN
END
REAL FUNCTION AI(t)
REAL w,t,AMax
common /comA/ w,AMax
COMMON P0,P,S,S1
AI=AMax*sin(w*t+S1)
END
REAL FUNCTION AD(t)
COMMON P0,P,S,S1,V,AM,AL,U,SS,ALX
REAL t
if(t>ALX)then
AD=0.4*SQRT(SQRT((t-ALX)/ALX))
else
AD=0.0
endif
END
Исходные данные: &data
S=0.043
S1=0.5
P0=7.E05
AM=100.
AL=0.21
AMax=12000.
ALX=.100
ALK=0.199
Y=0.0,0.0,7.E05,293.0
XOD=.00,.02,.04,.06,.08,.10,.12,.14,.16,.20
SILA=90000.,90000.,90000.,90000.,90000.,90000.,90000.,70000.,30000.,10000.
XD=.00,.02,.04,.06,.08,.10,.11,.14,.16,.20
SSLA=.0000,.0000,.0000,0.00E-00,0.00E-00,2.12E-03,2.12E-03,4.4E-03,4.4E-03,4.4E-03
/
&data2
TEMP=300.,500.,1000.,1700.,2300.,3000.,4000.,6000.,8000.,10000.,12000.,14000.,16000.,18000.,20000.,22000.,26000.,30000.,35000.,40000.
CP=138.1,175.5,245.6,335.6,662.0,978.0,1418.0,1965.0,2227.0,2387.0,2575.,2797.,3006.,3304.,3714.,4180.,5030.,5650.,6110.,6910.
TEMP2=300.,500.,1000.,1700.,2300.,3000.,4000.,6000.,8000.,10000.,12000.,14000.,16000.,18000.,20000.,22000.,26000.,30000.,35000.,40000.
RO=93.76,56.16,28.48,15.26,4.21,1.972,1.068,0.679,0.508,0.401,0.325,0.267,0.220,0.179,0.142,0.113,0.082,0.066,0.053,0.04
TEMP3=300.,500.,1000.,1700.,2300.,3000.,4000.,6000.,8000.,10000.,12000.,14000.,16000.,18000.,20000.,22000.,26000.,30000.,35000.,40000.
zh=1.,1.,1.,1.098,2.94,4.814,6.665,6.989,7.015,7.096,7.306,7.634,8.086,8.866,10.065,11.42,13.31,14.26,15.3,17.4
/
Рис. П1.1. Результаты численного расчета пневмомеханических характеристик (12kA): 1 - ход контактов, 2 - изменение скорости, 3 - изменение давления, 4 – температура, 5 – активное сечение сопел
Рис.П1.2. Распределения газодинамических функций (12kA): 1 – Rг(Т), 2 – Кг(Т)
Рис.П1.3. Результаты численного расчета пневмомеханических характеристик (63kA): 1 - ход контактов, 2 - изменение скорости, 3 - изменение давления, 4 – температура, 5 – активное сечение сопел
Рис.1.4. Распределения газодинамических функций (100kA): 1 – Rг(Т), 2 – Кг(Т)
|
