Tek-Tips is the largest IT community on the Internet today!
Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!
-
Congratulations Rhinorhino on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!
Convert Fortran 77 to Fortran 90
- Thread starter alikami
- Start date
- Thread starter
- #5
Thanks for replying...
Actually the program is running no problem in it..
I want to ,convert to matlab using f2matlab code.
But it accepts f90 ...so problem.
I tried to use convert.f90 but some problem
I am attaching the code , if someone can help
Regards
Aali
PARAMETER(NZM=100,NWM=50,NTM=200,PI=3.1415926)
EXTERNAL F1,FVC
LOGICAL FIRST
COMMON G,GE,AJ,RG,RP,TS,AS,AT
COMMON /COM/NZ,NW,M
COMMON /BDA/BA,BN,TI,AI
COMMON /FCDA/I,H,DTIME,IT
COMMON /S/GP1,GM1,SA,SB,SC
DIMENSION Y(3),V1(NZM+1),P1(NZM+1),T1(NZM+1),V2(NZM+1),
# P2(NZM+1),T2(NZM+1),V20(NZM+1),P20(NZM+1),T20(NZM+1),
# WW(NZM+1),WW1(NZM+1),W1(3,5),WM(NZM+1),P(NZM+1,NTM),
# T(NZM+1,NTM),V(NZM+1,NTM),R(NZM+1,NTM),TIME(NTM),
# BR(NZM+1),BR1(NZM+1),PMG(NTM),W0(NWM+1),
# AP(NZM+1,NWM+1),DAP(NZM+1),Z(NZM+1),APP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),FORCE(NTM),CF(NTM),
# ABH(NWM+1)
REAL LHEAD,LIMDA,ISPP(NTM),ISPF(NTM)
CHARACTER CHR*10
OPEN(1,FILE='INTER_I1.DAT',STATUS='OLD')
READ(1,*)CHR,RP
READ(1,*)CHR,PB
READ(1,*)CHR,DTWORK
READ(1,*)CHR,DTTAIL
READ(1,*)CHR,ICOR
READ(1,*)CHR,LCOR
READ(1,*)CHR,PCON1
READ(1,*)CHR,PCON2
READ(1,*)CHR,DE
READ(1,*)CHR,DT
READ(1,*)CHR,E
READ(1,*)CHR,G
READ(1,*)CHR,RG
READ(1,*)CHR,BA
READ(1,*)CHR,BN
READ(1,*)CHR,TS
READ(1,*)CHR,WMAX
READ(1,*)CHR,TI
READ(1,*)CHR,AI
READ(1,*)CHR,LHEAD
READ(1,*)CHR,NZ
READ(1,*)CHR,NW
CLOSE(1,STATUS='KEEP')
OPEN(2,FILE='INTER_I2.DAT',STATUS='OLD')
READ(2,*)(W0(J),J=1,NW+1)
READ(2,*)(WM(J),J=1,NZ+1)
READ(2,*)(Z(J),J=1,NZ+1)
READ(2,*)((AP(J,I),I=1,NW+1),J=1,NZ+1)
READ(2,*)VCH,AB0
READ(2,*)(ABH(I),I=1,NW+1)
CLOSE(2)
DO 5 I=1,NZ+1
IF(LHEAD.LE.ABS(Z(I)-Z(1))) THEN
II=I
ZLH=Z(I)
GOTO 10
ENDIF
5 CONTINUE
10 M=NZ+1-II+1
DO 15 I=1,M
Z(I)=Z(II+I-1)-ZLH
WM(I)=WM(II+I-1)
DO 15 J=1,NW+1
AP(I,J)=AP(II+I-1,J)
15 CONTINUE
AT=PI*DT**2/4.0
GP1=G+1.0
GM1=G-1.0
SC=0.5*(G+1.)/(G-1.0)
SA=(2.0/(G+1.))**SC
GE=SQRT(G)*SA
SB=0.5*(G-1.0)
AJ=AT/AP(M,1)
FIRST=.TRUE.
AS=SQRT(2.0*G*RG*TS/(G+1.0))
H=Z(2)-Z(1)
DO 20 I=1,M
V20(I)=100.0
P20(I)=5.0
T20(I)=100.0
PT(I)=100.
TT(I)=50.
VT(I)=30.
WW1(I)=0.0
WW(I)=0.
DAP(I)=0.0
APP(I)=AP(I,1)
20 CONTINUE
AEPS=1.E-4
BEPS=1.E-4
C30------CEPS is toll on mass flow rate------------------------------------C
CEPS=1.E-3
DTIME=0.0
A=0.
B=1.
C30--------------------------------------------------------------------C
C30-C35 Call the funtion FINE to calculate the velocity coefficient, C
C thus to calculate the value of first cross-section pressure C
C in the initial time. And give the number of time nodes ,the C
C----- time,the burning area an initinal value ---------C
CALL FINE(A,B,FVC,AEPS,LIMDA)
AIL=SIMP(0.,LIMDA,F1,BEPS,1E-2)
PO=(RP*BA*AB0*AS*2.*(G+1.)/(APP(M)*G*AIL))**(1./(1.-BN))
OPEN(2311,FILE='anlys1.out',STATUS='unknown')
write(2311,*)po
IT=0
TIM=0.0
ITER=0
AH=ABH(1)
25 LL=1
30 Y(2)=ABS(PO)
BB=RP*AH*BA*Y(2)**(BN-1.)*RG/APP(1)
CP=G/(G-1.0)*RG
AA=0.5*(BB**2)/CP
C35--------------------------------------------------------------------C
C35-C40 Calculate the values of following parameters:temperature, C
C---- velocity,pressure,combustion gas density ------C
IF(AA.EQ.0.) THEN
Y(3)=TS
ELSE
Y(3)=(-1.+SQRT(1.+4.*AA*TS))/(2.*AA)
ENDIF
Y(1)=BB*Y(3)
IF(Y(1).LE.1.E-3) Y(1)=0.01
RR=Y(2)/(RG*Y(3))
C40--------------------------------------------------------------------C
WRITE(*,*)'LL====',LL,' TIME===',TIM
C45--------------------------------------------------------------------C
C45-C50 Use the Merson method to solve the differential equation group,C
C then use the solved parameters
ressure,velocity,temperature, C
C call the subroutine FLOAD2 to aquire the value of the C
C following parameters:combustion gas density,total pressure C
C in the first cross-section,adiabatic burning temperature, C
C----- quadratic value of specific acoustic velocity. --------------C
X=Z(1)
DO 35 I=1,M
IF(I.EQ.1) GOTO 31
CALL MS(3,X,Y,CEPS,H,FIRST,W1,AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
31 CALL FLOAD2(Y,RR,PS,TSG,A2)
CALL BURN (Y(2),BR(I))
C------ CALL BURN (Y(2),BR(I),Y(1))---------C
V2(I)=Y(1)
P2(I)=Y(2)
T2(I)=Y(3)
35 CONTINUE
C50--------------------------------------------------------------------C
C50-C55 Middle-ring replacement----mass flow rate replacement ------C
IF(LL.EQ.LCOR) GOTO 45
GML=(Y(2)/(RG*Y(3)))*V2(M)*APP(M)
GMT=GE*AT*PS*PCON1/SQRT(RG*TSG)
AB=GMT-GML
WAB=AB/GMT
C55--------------------------------------------------------------------C
IF(ABS(WAB).LE.CEPS) GOTO 45
IF(LL.GE.2) GOTO 42
C60--------------------------------------------------------------------C
C60-C65 Use the ratio method to refine the value of pressure in the C
C---- initial time. --------C
AB1=AB
PO1=PO
PO=PO*GML/GMT
40 LL=LL+1
write(2311,2312)tim,ll,iter,po,br(1),gml,gmt,ab
2312 format(f5.2,1x,i3,1x,i3,1x,f10.1,1x,f8.5,1x,f6.2,1x,f6.2,
# 1x,f6.5,1x,f6.4)
GOTO 30
C65--------------------------------------------------------------------C
C65-C70 If the time is not the initial time,the linear interpolation C
C----- is used to refine the pressure value. -------C
42 P3=PO+AB*(PO-PO1)/(AB1-AB)
PO1=PO
AB1=AB
PO=P3
GOTO 40
C70--------------------------------------------------------------------C
45 IF(TIM.EQ.0.0.OR.ITER.EQ.ICOR) GOTO 55
C75--------------------------------------------------------------------C
DO 50 I=1,M
IF(ABS(V2(I)-V20(I))/V2(I).GT.CEPS) GOTO 62
IF(ABS(P2(I)-P20(I))/P2(I).GT.CEPS) GOTO 62
IF(ABS(T2(I)-T20(I))/T2(I).GT.CEPS) GOTO 62
50 CONTINUE
C80--------------------------------------------------------------------C
55 DO 60 I=1,M
V1(I)=V2(I)
P1(I)=P2(I)
T1(I)=T2(I)
BR1(I)=BR(I)
WW1(I)=WW(I)
60 CONTINUE
C85--------------------------------------------------------------------C
C85-C90 When the precision meets the requirements,time step length C
C is added to the time node,then in the second time the values C
C--- of those parameters are given. -------C
IT=IT+1
DO 65 I=1,M
P(I,IT)=P1(I)
T(I,IT)=T1(I)
V(I,IT)=V1(I)
R(I,IT)=P1(I)/(RG*T1(I))
65 CONTINUE
C90--------------------------------------------------------------------C
P(M,IT)=PCON1*P1(M)
TIME(IT)=TIM
C95--------------------------------------------------------------------C
IF(TIM.EQ.0.0) POO1=PO
PAV=PO/POO1
POO1=PO
PO=PO*PAV
PMG(IT)=GMT
DTIME=DTWORK
TIM=TIM+DTIME
ITER=0
GO TO 75
C100-------------------------------------------------------------------C
C100-C105 Calculate values of the partial differentiation of pressure, C
C--- temperature,velocity to the time -------C
62 DO 70 I=1,M
VT(I)=(V2(I)-V1(I))/DTIME
PT(I)=(P2(I)-P1(I))/DTIME
TT(I)=(T2(I)-T1(I))/DTIME
V20(I)=V2(I)
P20(I)=P2(I)
T20(I)=T2(I)
70 CONTINUE
C105-------------------------------------------------------------------C
75 ITER=ITER+1
WRITE(*,*)'ITER==',ITER
C110-------------------------------------------------------------------C
C110-C115 The outer-ring replacement ---- burning rate replacement ---C
DO 80 I=1,M
BR(I)=0.5*(BR1(I)+BR(I))
WW(I)=WW1(I)+BR(I)*DTIME
IF(WW(I).GT.WMAX) GOTO 85
IF(WW(I).GT.WM(I)) WW(I)=WM(I)
IF(WW(I).EQ.WM(I)) THEN
APP(I)=AP(I,NW+1)
ELSE
CALL GEOM(I,WW(I),APP(I),W0,AP,NZM,NWM)
ENDIF
80 CONTINUE
C115-------------------------------------------------------------------C
C115-C120 Use the interpolation method to calculate the head-end C
C burning area,And call the funtion FAT to calculate nozzle C
C throat area under the conditions of nozzle erosion or C
C------ unerosion ---------C
CALL AHEAD(WW(1),W0,ABH,AH,NWM)
CALL FAT(TIM,AT,E,DT)
GO TO 25
C120-------------------------------------------------------------------C
C120-C125 The subroutine TAIL is called to calculate the motor pressureC
C decreasing step when the combustion pressure 'Pc'is less thanC
C---- one tenth of Pcmax,the motor is considered to end . --------C
85 ITWORK=IT
PMAX=0.0
DO 90 I=1,ITWORK
IF(PMAX.LT.P(1,I)) PMAX=P(1,I)
90 CONTINUE
CALL TAIL(IT,P,R,T,V,AP,PMG,TIME,ITWORK,DTTAIL,VCH,NZM,
# NWM,NTM,PMAX)
ITTAIL=IT
CALL OUTPUT(Z,P,V,R,T,PMG,CF,FORCE,TIME,ISPP,ISPF,ITWORK,
# ITTAIL,PB,PCON2,DE,E,DT,NZM,NTM,PMAX)
C125-------------------------------------------------------------------C
END
C*************************SUB OUTPUT***********************************C
C This subroutine is used to calculate and output the required C
C values of inter ballistic parameters and the the stress C
C----- adjustment. ---------C
SUBROUTINE OUTPUT(Z,P,V,R,T,PMG,CF,FORCE,TIME,ISPP,ISPF,
# ITWORK,ITTAIL,PB,PCON2,DE,E,DT,NZM,NTM,PMAX)
C1---------------------------------------------------------------------C
PARAMETER(PI=3.1415926)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
COMMON /COM/NZ,NW,M
DIMENSION Z(NZM+1),P(NZM+1,NTM),V(NZM+1,NTM),R(NZM+1,NTM),
# T(NZM+1,NTM),FORCE(NTM),PMG(NTM),TIME(NTM),CF(NTM)
REAL IP,ISPP(NTM),ISPF(NTM),MTAV,ISPAV,LIMDA
C5---------------------------------------------------------------------C
OPEN(1,FILE='VERST_I2.DAT',STATUS='UNKNOWN')
C10--------------------------------------------------------------------C
AE=PI*DE**2/4.0
DO 15 I=1,ITTAIL
CALL FAT(TIME(I),AT,E,DT)
EPS=AE/AT
CALL PARA (LIMDA,EPS)
PAIP=1.0-(G-1.)/(G+1.)*LIMDA*LIMDA
CFO=GE*SQRT(2.0*G/(G-1.)*(1.0-PAIP))
CF(I)=(CFO+EPS*(PAIP**(G/(G-1.))-PB/P(M,I)))*PCON2
IF(CF(I).LT.0.0) CF(I)=0.0
FORCE(I)=CF(I)*P(M,I)*AT
15 CONTINUE
C15--------------------------------------------------------------------C
PP=0.0
FF=0.0
PMM=0.0
DO 10 I=1,ITTAIL-1
DTT=TIME(I+1)-TIME(I)
PP=PP+0.5*(P(M,I)+P(M,I+1))*DTT
FF=FF+0.5*(FORCE(I)+FORCE(I+1))*DTT
PMM=PMM+0.5*(PMG(I)+PMG(I+1))*DTT
ISPP(I+1)=PP/PMM
ISPF(I+1)=FF/PMM
10 CONTINUE
C20--------------------------------------------------------------------C
ISPP(1)=0.0
ISPF(1)=0.0
PCAV=PP/TIME(ITTAIL)
FOAV=FF/TIME(ITTAIL)
MTAV=PMM/TIME(ITTAIL)
IP=FF
ISPAV=FOAV/MTAV
C25--------------------------------------------------------------------C
WRITE(1,*)'''PMAX'',',PMAX,','
CLOSE(1)
C30--------------------------------------------------------------------C
C I_PS_Z------- the stated pressure of every node C
C versus Z-corrdinate. C
C I_PE_Z------- the pressure of every node C
C versus Z-corrdinate. C
C the following 'I_TS_Z,I_TE_Z,I_VS_Z,I_VE_Z,I_RS_Z,I_RE_Z' C
C are used like the above example . C
C----------------------------------------------------------------------C
OPEN(1,FILE='I_PS_Z.DAT',STATUS='UNKNOWN')
DO 25 I=1,M
P(I,1)=P(I,1)/1000000.
WRITE(1,*)Z(I),P(I,1)
25 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_PE_Z.DAT',STATUS='UNKNOWN')
DO 30 I=1,M
P(I,ITWORK)=P(I,ITWORK)/1000000.
WRITE(1,*)Z(I),P(I,ITWORK)
30 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TS_Z.DAT',STATUS='UNKNOWN')
DO 35 I=1,M
WRITE(1,*)Z(I),T(I,1)
35 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TE_Z.DAT',STATUS='UNKNOWN')
DO 40 I=1,M
WRITE(1,*)Z(I),T(I,ITWORK)
40 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VS_Z.DAT',STATUS='UNKNOWN')
DO 45 I=1,M
WRITE(1,*)Z(I),V(I,1)
45 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VE_Z.DAT',STATUS='UNKNOWN')
DO 50 I=1,M
WRITE(1,*)Z(I),V(I,ITWORK)
50 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RS_Z.DAT',STATUS='UNKNOWN')
DO 55 I=1,M
WRITE(1,*)Z(I),R(I,1)
55 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RE_Z.DAT',STATUS='UNKNOWN')
DO 60 I=1,M
WRITE(1,*)Z(I),R(I,ITWORK)
60 CONTINUE
CLOSE(1)
C35--------------------------------------------------------------------C
C I_PH_T----- the pressure in the head-end cross-section C
C versus time . C
C I_PE_T----- the pressure in the aft-end cross-section C
C versus time . C
C the following 'I_TH_T,I_TE_T,I_VH_T,I_VE_T,I_RH_T,I_RE_T' C
C are used like the above example. C
C----------------------------------------------------------------------C
OPEN(1,FILE='I_PH_T.DAT',STATUS='UNKNOWN')
P(1,1)=P(1,1)*1000000.
P(1,ITWORK)=P(1,ITWORK)*1000000.
DO 65 I=1,ITTAIL
C------- In file I_PH_T.DAT, pressure unit is MPa.--------------------C
P(1,I)=P(1,I)/1000000.
C----------------------------------------------------------------------C
WRITE(1,*)TIME(I),P(1,I)
65 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_PE_T.DAT',STATUS='UNKNOWN')
P(M,1)=P(M,1)*1000000.
P(M,ITWORK)=P(M,ITWORK)*1000000.
DO 70 I=1,ITTAIL
C------- In file I_PE_T.DAT, pressure unit is MPa.--------------------C
P(M,I)=P(M,I)/1000000.
C----------------------------------------------------------------------C
WRITE(1,*)TIME(I),P(M,I)
70 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TH_T.DAT',STATUS='UNKNOWN')
DO 75 I=1,ITTAIL
WRITE(1,*)TIME(I),T(1,I)
75 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TE_T.DAT',STATUS='UNKNOWN')
DO 80 I=1,ITTAIL
WRITE(1,*)TIME(I),T(M,I)
80 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VH_T.DAT',STATUS='UNKNOWN')
DO 85 I=1,ITTAIL
WRITE(1,*)TIME(I),V(1,I)
85 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VE_T.DAT',STATUS='UNKNOWN')
DO 90 I=1,ITTAIL
WRITE(1,*)TIME(I),V(M,I)
90 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RH_T.DAT',STATUS='UNKNOWN')
DO 95 I=1,ITTAIL
WRITE(1,*)TIME(I),R(1,I)
95 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RE_T.DAT',STATUS='UNKNOWN')
DO 100 I=1,ITTAIL
WRITE(1,*)TIME(I),R(M,I)
100 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_CF_T.DAT',STATUS='UNKNOWN')
DO 105 I=1,ITTAIL
WRITE(1,*)TIME(I),CF(I)
105 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_F_T.DAT',STATUS='UNKNOWN')
DO 110 I=1,ITTAIL
WRITE(1,*)TIME(I),FORCE(I)
110 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_M_T.DAT',STATUS='UNKNOWN')
DO 115 I=1,ITTAIL
WRITE(1,*)TIME(I),PMG(I)
115 CONTINUE
CLOSE(1)
C OPEN(1,FILE='I_ISP_T.DAT',STATUS='UNKNOWN')
C DO 120 I=1,ITTAIL
C WRITE(1,*)TIME(I),ISPP(I)
C120 CONTINUE
C CLOSE(1)
OPEN(1,FILE='I_ISP_T.DAT',STATUS='UNKNOWN')
DO 125 I=1,ITTAIL
WRITE(1,*)TIME(I),ISPF(I)
125 CONTINUE
CLOSE(1)
C40--------------------------------------------------------------------C
OPEN(1,FILE='INTER_O.DAT',STATUS='UNKNOWN')
WRITE(1,*)'The internal ballistic performance:'
WRITE(1,*)'Time node number of work time : ITWORK=',ITWORK
WRITE(1,*)'Work time of rocket engine : ta=',TIME(ITWORK)
WRITE(1,*)'Time node number including work time and tail time :'
WRITE(1,*)' ITTAIL=',ITTAIL
WRITE(1,*)'Sum of work time and tail time :tb=',TIME(ITTAIL)
WRITE(1,*)'Maximum chamber pressure : Pmax=',PMAX
WRITE(1,*)'Average chamber pressure : Pcav=',PCAV
WRITE(1,*)'Average force : Fav=',FOAV
WRITE(1,*)'Average mass flow rate : Mav=',MTAV
WRITE(1,*)'Total impulse : Ip=',IP
WRITE(1,*)'Average specific impulse : Isp=',ISPAV
CLOSE(1)
C45--------------------------------------------------------------------C
RETURN
END
C*********************SUB MS******************************************C
C This subroutine is MERSON METHOD .it is used to solve the C
C---- diffierential equation group. ----------C
SUBROUTINE MS(N,X,Y,EPS,H,FIRST,W,AP,APP,DAP,W0,WW,VT,PT,TT,
# NZM,NWM)
C1---------------------------------------------------------------------C
LOGICAL INC,FIRST
DIMENSION Y(N),W(N,5),AP(NZM+1,NWM+1),APP(NZM+1),DAP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),W0(NWM+1),WW(NZM+1)
C5---------------------------------------------------------------------C
IF(.NOT.FIRST) GO TO 2
LOCP=1
FIRST=.FALSE.
2 HC=H/FLOAT(LOCP)
LOC=0
4 CALL FCT(X,Y,W(1,3),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
DO 6 I=1,N
W(I,1)=W(I,3)*HC/3.0+Y(I)
6 CONTINUE
CALL FCT(HC/3.0+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 8 I=1,N
W(I,1)=(W(I,3)+W(I,4))*HC/6.0+Y(I)
8 CONTINUE
CALL FCT(HC/3.0+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 10 I=1,N
W(I,1)=(W(I,3)+3.0*W(I,4))*HC*0.125+Y(I)
10 CONTINUE
CALL FCT(HC*0.5+X,W(1,1),W(1,5),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 12 I=1,N
W(I,1)=W(I,3)*HC*0.5-W(I,4)*HC*1.5+W(I,5)*HC*2.0+Y(I)
12 CONTINUE
CALL FCT(HC+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
DO 14 I=1,N
W(I,2)=(W(I,3)+W(I,4)+W(I,5)*4.0)*HC/6.0+Y(I)
14 CONTINUE
C10--------------------------------------------------------------------C
INC=.TRUE.
DO 18 I=1,N
ER=ABS(W(I,1)-W(I,2))*0.2
IF(ABS(W(I,2)).GT.1.0) ER=ER/ABS(W(I,2))
IF(ER.LT.EPS) GO TO 18
HC=HC*0.5
LOCP=LOCP+LOCP
LOC=LOC+LOC
GO TO 4
18 IF(64.0*ER.GT.EPS) INC=.FALSE.
C15--------------------------------------------------------------------C
X=X+HC
DO 20 I=1,N
Y(I)=W(I,2)
20 CONTINUE
LOC=LOC+1
IF(LOC.EQ.LOCP) RETURN
C20--------------------------------------------------------------------C
IF(.NOT.INC.OR.LOC.NE.LOC/2*2) GO TO 4
C25--------------------------------------------------------------------C
HC=2.0*HC
LOC=LOC/2
LOCP=LOCP/2
GO TO 4
C30--------------------------------------------------------------------C
END
C*******************SUB FCT*******************************************C
C This subroutine is used to calculate the values of all the C
C---- right of the equations in the differential equation groups.---C
SUBROUTINE FCT(X,Y,YX,AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
C1--------------------------------------------------------------------C
COMMON /FCDA/I,H,DTIME,IT
COMMON G,GE,AJ,RG,RP,TS,AS,AT
DIMENSION Y(3),YX(3),AP(NZM+1,NWM+1),APP(NZM+1),DAP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),W0(NWM+1),WW(NZM+1)
C5--------------------------------------------------------------------C
K=I
IF(I.EQ.2) K=I+1
K2=K-2
K1=K-1
Q=X/H-FLOAT(K2)
H2=H
CALL FLOAD2(Y,R,PS,TSG,A2)
CALL BURN(Y(2),BR)
C------ CALL BURN(Y(2),BR,Y(1))----------C
C10-------------------------------------------------------------------C
DO 73 II=1,K
CALL GEOM(II, DAP(II)=APP1-APP(II)
73 CONTINUE
C15--------------------------------------------------------------------C
APT=FX(APP(K2),APP(K1),APP(K),Q,H,H2)
DAPP=FX(DAP(K2),DAP(K1),DAP(K),Q,H,H2)
DAPT=DFX(APP(K2),APP(K1),APP(K),Q,H,H2)
VTT=FX(VT(K2),VT(K1),VT(K),Q,H,H2)
PTT=FX(PT(K2),PT(K1),PT(K),Q,H,H2)
TTT=FX(TT(K2),TT(K1),TT(K),Q,H,H2)
A=((RP-R)*DAPP-R*Y(1)*DAPT)/APT
B=-RP*Y(1)*DAPP/(APT*R)
C=(DAPP*(RP*G*RG*TS+0.5*(G-1.)*RP*Y(1)*Y(1)-Y(2))
# -G*Y(2)*Y(1)*DAPT)/APT
DV=A2-Y(1)*Y(1)
YX(1)=(Y(1)*(VTT-B)+(C-PTT)/R)/DV
YX(2)=(G*Y(2)*(B-VTT)+Y(1)*(PTT-C))/DV
YX(3)=Y(3)*(C-RG*Y(3)*A)/(Y(1)*Y(2))-(G-1.)*Y(3)/Y(1)*YX(1)-
# TTT/Y(1)
C20-------------------------------------------------------------------C
RETURN
END
C*******************SUB FLOAD2*****************************************C
C This subroutine is used to calculate some parameters such as: C
C---- combustion gas density,specific heat ratio,etc .---------------C
SUBROUTINE FLOAD2(Y,R,PS,TSG,A2)
DIMENSION Y(3)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
REAL LIMDA
R=Y(2)/(Y(3)*RG)
CP=G*RG/(G-1.)
TSG=Y(3)+Y(1)**2/(2.*CP)
AC=SQRT(2.*G*RG*TSG/(G+1.))
LIMDA=Y(1)/AC
A=1.0-LIMDA*LIMDA*(G-1.0)/(G+1.0)
PS=Y(2)*A**(G/(1.0-G))
A2=G*RG*Y(3)
RETURN
END
C****************SUB BURN*********************************************C
C This subroutine is used to calculate the propellant burning C
C----- rate. ---------C
SUBROUTINE BURN(P,BR)
C------- SUBROUTINE BURN(P,BR,V) ---------C
COMMON /BDA/BA,BN,TI,AI
C=EXP(AI*(TI-20.0))
C------- PARAMETER(VTH=*****,VKU=****) ---C
C------- IF(V.LE.VTH) THEN ----------C
C------- BR=BA*C*P**BN ---------------C
C------- ELSE --------------------------C
C------- BR=(1.0+VKU*(V-VTH))*BA*C*P**BN
C------- ENDIF --------------------------C
BR=BA*C*P**BN
RETURN
END
C***********************FUN FX*****************************************C
C----- This function 'FX' is the Lagrange interpolation . --------C
FUNCTION FX(YO,Y1,Y2,Q,H1,H2)
FX=(Q*Q*H1*H1*(H2*YO-(H1+H2)*Y1+H1*Y2)-Q*H1*(H2*H2*YO+(H1*H1
# -H2*H2)*Y1-H1*H1*Y2)+H1*H2*(H1+H2)*Y1)/(H1*H2*(H1+H2))
RETURN
END
C**********************FUN DFX*****************************************C
C---- This function DFX is a differential function -------C
FUNCTION DFX(YO,Y1,Y2,Q,H1,H2)
DFX=(2.0*Q*H1*(H2*YO-(H1+H2)*Y1+H1*Y2)-(H2*H2*YO+
# (H1*H1-H2*H2)*Y1-H1*H1*Y2))/(H1*H2*(H1+H2))
RETURN
END
C********************SUB FAT*******************************************C
C This subroutine 'FAT' is used to calculate the nozzle throat C
C----- area under the conditions of nozzle erosion or unerosion -----C
SUBROUTINE FAT(TIME,AT,E,DT)
PARAMETER(PI=3.1415926)
AT=PI*DT**2/4.
IF(TIME.GE.15) AT=PI*(DT+2*E*(TIME-15))**2/4.
RETURN
END
C***********************SUB AHEAD**************************************C
C This subroutine is used to calculate the head-end burning area C
C----- according to the propellant grain geometry calculations. ----C
SUBROUTINE AHEAD(WW,W0,ABH,AH,NWM)
DIMENSION W0(NWM+1),ABH(NWM+1)
J=1
10 J1=J+1
IF(WW.GE.W0(J).AND.WW.LE.W0(J1)) GOTO 20
J=J1
GOTO 10
20 AH=FLINE(ABH(J),ABH(J1),W0(J),W0(J1),WW)
RETURN
END
C**********************FUN FLINE***************************************C
C----- Use the linear interpolation to interpolate the given nodes.---C
FUNCTION FLINE(Y1,Y2,X1,X2,X)
FLINE=Y1+(Y2-Y1)*(X-X1)/(X2-X1)
RETURN
END
C*********************SUB GEOM*****************************************C
C This subroutine calls the interpolation function subroutine FX C
C to interpolate the propellant grain geometry,then gets the C
C value of the area's averge ratio of change to the time,in C
C---- every time node. ----------C
SUBROUTINE GEOM(I,WW,APP,W0,AP,NZM,NWM)
COMMON /COM/NZ,NW,M
DIMENSION W0(NWM+1),AP(NZM+1,NWM+1)
IF(WW.GE.W0(NW+1)) THEN
APP=AP(I,NW+1)+(AP(I,NW+1)-AP(I,NW))/(W0(NW+1)-
# W0(NW))*(WW-W0(NW+1))
RETURN
ENDIF
J=1
10 JP1=J+1
IF(WW.GE.W0(J).AND.WW.LE.W0(JP1)) THEN
APP=FLINE(AP(I,J),AP(I,JP1),W0(J),W0(JP1),WW)
ELSE
J=JP1
GOTO 10
ENDIF
RETURN
END
C*******************SUB TAIL******************************************C
C This subroutine is used to calculate and output these para- C
C----- meters' values in the tail phase. ---------------C
SUBROUTINE TAIL(IT,P,R,T,V,AP,PMG,TIME,ITWORK,DTTAIL,VCH,
# NZM,NWM,NTM,PMAX)
C1---------------------------------------------------------------------C
COMMON /COM/NZ,NW,M
DIMENSION P(NZM+1,NTM),R(NZM+1,NTM),T(NZM+1,NTM),V(NZM+1,NTM),
# PMG(NTM),TIME(NTM),AP(NZM+1,NWM+1)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
C5---------------------------------------------------------------------C
TIM=TIME(IT)
TT=0.0
A=2.0*VCH
A1=GE*SQRT(RG*TS)*AT*(G-1.0)
CD=GE/SQRT(RG*TS)
C=2.0*G/(G-1.0)
D=0.1*PMAX
C10--------------------------------------------------------------------C
10 TIM=TIM+DTTAIL
TT=TT+DTTAIL
IT=IT+1
B=A+A1*TT
P(1,IT)=P(1,ITWORK)*(A/B)**C
R(1,IT)=R(1,ITWORK)*(P(1,IT)/P(1,ITWORK))**(1./G)
T(1,IT)=P(1,IT)/(RG*R(1,IT))
V(1,IT)=0.0
P(M,IT)=P(M,ITWORK)*(A/B)**C
R(M,IT)=R(M,ITWORK)*(P(M,IT)/P(M,ITWORK))**(1./G)
T(M,IT)=P(M,IT)/(RG*R(M,IT))
PMG(IT)=CD*P(M,IT)*AT
V(M,IT)=PMG(IT)/(R(M,IT)*AP(M,NW+1))
TIME(IT)=TIM
C15--------------------------------------------------------------------C
IF(P(M,IT).LE.D) GOTO 20
GOTO 10
C20--------------------------------------------------------------------C
20 RETURN
END
C*******************SUB PARA*******************************************C
C This subroutine is called to calculate the velocity coefficientC
C in the nozzle exit cross-section by calling the subroutine C
C----- SMACH which is used to calculate it . -------------C
SUBROUTINE PARA(LIMDA,EPS)
COMMON /S/GP1,GM1,SA,SB,SC
REAL M1,M2,LIMDA
M1=1.01
M2=4.0
CALL SMACH(EPS,EM,M1,M2)
EM2=EM*EM
LIMDA=SQRT(0.5*GP1*EM2/(1.0+GM1*0.5*EM2))
RETURN
END
C*********************SUB SMACH****************************************C
C This subroutine is used to calculate the mach number in the C
C----- nozzle exit cross-section. ---------C
SUBROUTINE SMACH(E,M,M1,M2)
COMMON /S/GP1,GM1,SA,SB,SC
REAL M,M1,M2,MO
FUN(X,SA,SB,SC)=SA*(1.0+SB*X*X)**SC/X
F2=FUN(M2,SA,SB,SC)
10 F10=FUN(M1,SA,SB,SC)
T=(F2-F10)/(M2-M1)
MO=M1+(E-F10)/T
IF(ABS(MO-M1).LE.1E-5) GO TO 20
M1=MO
GO TO 10
20 M=MO
RETURN
END
C*********************SUB FINE*****************************************C
C This subroutine calls the function FVC,and uses the Two-cuttingC
C---- Method to calculate the solution of the given function .-------C
SUBROUTINE FINE(A,B,F,AEPS,LIMDA)
LOGICAL FAG
REAL LIMDA
FAG=F(A).GT.0.0
5 C=0.5*(A+B)
IF((B-A).GE.AEPS) GO TO 20
10 LIMDA=C
RETURN
20 IF(F(C)) 25,10,35
25 IF(FAG) GO TO 40
30 A=C
GO TO 5
35 IF(FAG) GO TO 30
40 B=C
GO TO 5
END
C*********************FUN SIMP*****************************************C
C---- Simpson integration method . ---------C
FUNCTION SIMP(A,B,F,EPS,HO)
H=B-A
RP=F(A)+F(B)
N=1
2 X=A-0.5*H
RC=0.0
DO 4 I=1,N
X=X+H
RC=F(X)+RC
4 CONTINUE
R2=(4.0*RC+RP)*H/6.0
IF(N.EQ.1.OR.ABS(H).GT.HO) GOTO 6
X=R2-R1
IF(ABS(R2).GE.1.0) X=X/R2
IF(ABS(X).LT.EPS) GO TO 8
6 R1=R2
RP=2.0*RC+RP
N=N+N
H=0.5*H
GO TO 2
8 SIMP=R2
RETURN
END
C********************FUN F1********************************************C
C It is used to as an extenal function of Simpson integration . C
C----- its value is returned by F1. -------------C
FUNCTION F1(LIMDA)
COMMON /BDA/BA,BN,TI,AI
COMMON G,GE,AJ,RG,RP,TS,AS,AT
REAL LIMDA
A=(1.0+LIMDA*LIMDA)**(2.0-BN)
B=(1.0-LIMDA*LIMDA*(G-1.0)/(G+1.0))**BN
C=1.0
F1=4.0*(1.0-LIMDA*LIMDA)/(A*B*C)
RETURN
END
C******************FUN FVC*********************************************C
C This function helps subroutine FINE to solve the velocity coe- C
C---- fficient in the combustion chamber exit cross-section . -------C
FUNCTION FVC(X)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
FVC=AJ-X*(1.0-X*X*(G-1.0)/(G+1.0))**(1.0/(G-1.0))*
# ((G+1.0)/2.0)**(1.0/(G-1.0))
RETURN
END
Actually the program is running no problem in it..
I want to ,convert to matlab using f2matlab code.
But it accepts f90 ...so problem.
I tried to use convert.f90 but some problem
I am attaching the code , if someone can help
Regards
Aali
PARAMETER(NZM=100,NWM=50,NTM=200,PI=3.1415926)
EXTERNAL F1,FVC
LOGICAL FIRST
COMMON G,GE,AJ,RG,RP,TS,AS,AT
COMMON /COM/NZ,NW,M
COMMON /BDA/BA,BN,TI,AI
COMMON /FCDA/I,H,DTIME,IT
COMMON /S/GP1,GM1,SA,SB,SC
DIMENSION Y(3),V1(NZM+1),P1(NZM+1),T1(NZM+1),V2(NZM+1),
# P2(NZM+1),T2(NZM+1),V20(NZM+1),P20(NZM+1),T20(NZM+1),
# WW(NZM+1),WW1(NZM+1),W1(3,5),WM(NZM+1),P(NZM+1,NTM),
# T(NZM+1,NTM),V(NZM+1,NTM),R(NZM+1,NTM),TIME(NTM),
# BR(NZM+1),BR1(NZM+1),PMG(NTM),W0(NWM+1),
# AP(NZM+1,NWM+1),DAP(NZM+1),Z(NZM+1),APP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),FORCE(NTM),CF(NTM),
# ABH(NWM+1)
REAL LHEAD,LIMDA,ISPP(NTM),ISPF(NTM)
CHARACTER CHR*10
OPEN(1,FILE='INTER_I1.DAT',STATUS='OLD')
READ(1,*)CHR,RP
READ(1,*)CHR,PB
READ(1,*)CHR,DTWORK
READ(1,*)CHR,DTTAIL
READ(1,*)CHR,ICOR
READ(1,*)CHR,LCOR
READ(1,*)CHR,PCON1
READ(1,*)CHR,PCON2
READ(1,*)CHR,DE
READ(1,*)CHR,DT
READ(1,*)CHR,E
READ(1,*)CHR,G
READ(1,*)CHR,RG
READ(1,*)CHR,BA
READ(1,*)CHR,BN
READ(1,*)CHR,TS
READ(1,*)CHR,WMAX
READ(1,*)CHR,TI
READ(1,*)CHR,AI
READ(1,*)CHR,LHEAD
READ(1,*)CHR,NZ
READ(1,*)CHR,NW
CLOSE(1,STATUS='KEEP')
OPEN(2,FILE='INTER_I2.DAT',STATUS='OLD')
READ(2,*)(W0(J),J=1,NW+1)
READ(2,*)(WM(J),J=1,NZ+1)
READ(2,*)(Z(J),J=1,NZ+1)
READ(2,*)((AP(J,I),I=1,NW+1),J=1,NZ+1)
READ(2,*)VCH,AB0
READ(2,*)(ABH(I),I=1,NW+1)
CLOSE(2)
DO 5 I=1,NZ+1
IF(LHEAD.LE.ABS(Z(I)-Z(1))) THEN
II=I
ZLH=Z(I)
GOTO 10
ENDIF
5 CONTINUE
10 M=NZ+1-II+1
DO 15 I=1,M
Z(I)=Z(II+I-1)-ZLH
WM(I)=WM(II+I-1)
DO 15 J=1,NW+1
AP(I,J)=AP(II+I-1,J)
15 CONTINUE
AT=PI*DT**2/4.0
GP1=G+1.0
GM1=G-1.0
SC=0.5*(G+1.)/(G-1.0)
SA=(2.0/(G+1.))**SC
GE=SQRT(G)*SA
SB=0.5*(G-1.0)
AJ=AT/AP(M,1)
FIRST=.TRUE.
AS=SQRT(2.0*G*RG*TS/(G+1.0))
H=Z(2)-Z(1)
DO 20 I=1,M
V20(I)=100.0
P20(I)=5.0
T20(I)=100.0
PT(I)=100.
TT(I)=50.
VT(I)=30.
WW1(I)=0.0
WW(I)=0.
DAP(I)=0.0
APP(I)=AP(I,1)
20 CONTINUE
AEPS=1.E-4
BEPS=1.E-4
C30------CEPS is toll on mass flow rate------------------------------------C
CEPS=1.E-3
DTIME=0.0
A=0.
B=1.
C30--------------------------------------------------------------------C
C30-C35 Call the funtion FINE to calculate the velocity coefficient, C
C thus to calculate the value of first cross-section pressure C
C in the initial time. And give the number of time nodes ,the C
C----- time,the burning area an initinal value ---------C
CALL FINE(A,B,FVC,AEPS,LIMDA)
AIL=SIMP(0.,LIMDA,F1,BEPS,1E-2)
PO=(RP*BA*AB0*AS*2.*(G+1.)/(APP(M)*G*AIL))**(1./(1.-BN))
OPEN(2311,FILE='anlys1.out',STATUS='unknown')
write(2311,*)po
IT=0
TIM=0.0
ITER=0
AH=ABH(1)
25 LL=1
30 Y(2)=ABS(PO)
BB=RP*AH*BA*Y(2)**(BN-1.)*RG/APP(1)
CP=G/(G-1.0)*RG
AA=0.5*(BB**2)/CP
C35--------------------------------------------------------------------C
C35-C40 Calculate the values of following parameters:temperature, C
C---- velocity,pressure,combustion gas density ------C
IF(AA.EQ.0.) THEN
Y(3)=TS
ELSE
Y(3)=(-1.+SQRT(1.+4.*AA*TS))/(2.*AA)
ENDIF
Y(1)=BB*Y(3)
IF(Y(1).LE.1.E-3) Y(1)=0.01
RR=Y(2)/(RG*Y(3))
C40--------------------------------------------------------------------C
WRITE(*,*)'LL====',LL,' TIME===',TIM
C45--------------------------------------------------------------------C
C45-C50 Use the Merson method to solve the differential equation group,C
C then use the solved parameters
ressure,velocity,temperature, CC call the subroutine FLOAD2 to aquire the value of the C
C following parameters:combustion gas density,total pressure C
C in the first cross-section,adiabatic burning temperature, C
C----- quadratic value of specific acoustic velocity. --------------C
X=Z(1)
DO 35 I=1,M
IF(I.EQ.1) GOTO 31
CALL MS(3,X,Y,CEPS,H,FIRST,W1,AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
31 CALL FLOAD2(Y,RR,PS,TSG,A2)
CALL BURN (Y(2),BR(I))
C------ CALL BURN (Y(2),BR(I),Y(1))---------C
V2(I)=Y(1)
P2(I)=Y(2)
T2(I)=Y(3)
35 CONTINUE
C50--------------------------------------------------------------------C
C50-C55 Middle-ring replacement----mass flow rate replacement ------C
IF(LL.EQ.LCOR) GOTO 45
GML=(Y(2)/(RG*Y(3)))*V2(M)*APP(M)
GMT=GE*AT*PS*PCON1/SQRT(RG*TSG)
AB=GMT-GML
WAB=AB/GMT
C55--------------------------------------------------------------------C
IF(ABS(WAB).LE.CEPS) GOTO 45
IF(LL.GE.2) GOTO 42
C60--------------------------------------------------------------------C
C60-C65 Use the ratio method to refine the value of pressure in the C
C---- initial time. --------C
AB1=AB
PO1=PO
PO=PO*GML/GMT
40 LL=LL+1
write(2311,2312)tim,ll,iter,po,br(1),gml,gmt,ab
2312 format(f5.2,1x,i3,1x,i3,1x,f10.1,1x,f8.5,1x,f6.2,1x,f6.2,
# 1x,f6.5,1x,f6.4)
GOTO 30
C65--------------------------------------------------------------------C
C65-C70 If the time is not the initial time,the linear interpolation C
C----- is used to refine the pressure value. -------C
42 P3=PO+AB*(PO-PO1)/(AB1-AB)
PO1=PO
AB1=AB
PO=P3
GOTO 40
C70--------------------------------------------------------------------C
45 IF(TIM.EQ.0.0.OR.ITER.EQ.ICOR) GOTO 55
C75--------------------------------------------------------------------C
DO 50 I=1,M
IF(ABS(V2(I)-V20(I))/V2(I).GT.CEPS) GOTO 62
IF(ABS(P2(I)-P20(I))/P2(I).GT.CEPS) GOTO 62
IF(ABS(T2(I)-T20(I))/T2(I).GT.CEPS) GOTO 62
50 CONTINUE
C80--------------------------------------------------------------------C
55 DO 60 I=1,M
V1(I)=V2(I)
P1(I)=P2(I)
T1(I)=T2(I)
BR1(I)=BR(I)
WW1(I)=WW(I)
60 CONTINUE
C85--------------------------------------------------------------------C
C85-C90 When the precision meets the requirements,time step length C
C is added to the time node,then in the second time the values C
C--- of those parameters are given. -------C
IT=IT+1
DO 65 I=1,M
P(I,IT)=P1(I)
T(I,IT)=T1(I)
V(I,IT)=V1(I)
R(I,IT)=P1(I)/(RG*T1(I))
65 CONTINUE
C90--------------------------------------------------------------------C
P(M,IT)=PCON1*P1(M)
TIME(IT)=TIM
C95--------------------------------------------------------------------C
IF(TIM.EQ.0.0) POO1=PO
PAV=PO/POO1
POO1=PO
PO=PO*PAV
PMG(IT)=GMT
DTIME=DTWORK
TIM=TIM+DTIME
ITER=0
GO TO 75
C100-------------------------------------------------------------------C
C100-C105 Calculate values of the partial differentiation of pressure, C
C--- temperature,velocity to the time -------C
62 DO 70 I=1,M
VT(I)=(V2(I)-V1(I))/DTIME
PT(I)=(P2(I)-P1(I))/DTIME
TT(I)=(T2(I)-T1(I))/DTIME
V20(I)=V2(I)
P20(I)=P2(I)
T20(I)=T2(I)
70 CONTINUE
C105-------------------------------------------------------------------C
75 ITER=ITER+1
WRITE(*,*)'ITER==',ITER
C110-------------------------------------------------------------------C
C110-C115 The outer-ring replacement ---- burning rate replacement ---C
DO 80 I=1,M
BR(I)=0.5*(BR1(I)+BR(I))
WW(I)=WW1(I)+BR(I)*DTIME
IF(WW(I).GT.WMAX) GOTO 85
IF(WW(I).GT.WM(I)) WW(I)=WM(I)
IF(WW(I).EQ.WM(I)) THEN
APP(I)=AP(I,NW+1)
ELSE
CALL GEOM(I,WW(I),APP(I),W0,AP,NZM,NWM)
ENDIF
80 CONTINUE
C115-------------------------------------------------------------------C
C115-C120 Use the interpolation method to calculate the head-end C
C burning area,And call the funtion FAT to calculate nozzle C
C throat area under the conditions of nozzle erosion or C
C------ unerosion ---------C
CALL AHEAD(WW(1),W0,ABH,AH,NWM)
CALL FAT(TIM,AT,E,DT)
GO TO 25
C120-------------------------------------------------------------------C
C120-C125 The subroutine TAIL is called to calculate the motor pressureC
C decreasing step when the combustion pressure 'Pc'is less thanC
C---- one tenth of Pcmax,the motor is considered to end . --------C
85 ITWORK=IT
PMAX=0.0
DO 90 I=1,ITWORK
IF(PMAX.LT.P(1,I)) PMAX=P(1,I)
90 CONTINUE
CALL TAIL(IT,P,R,T,V,AP,PMG,TIME,ITWORK,DTTAIL,VCH,NZM,
# NWM,NTM,PMAX)
ITTAIL=IT
CALL OUTPUT(Z,P,V,R,T,PMG,CF,FORCE,TIME,ISPP,ISPF,ITWORK,
# ITTAIL,PB,PCON2,DE,E,DT,NZM,NTM,PMAX)
C125-------------------------------------------------------------------C
END
C*************************SUB OUTPUT***********************************C
C This subroutine is used to calculate and output the required C
C values of inter ballistic parameters and the the stress C
C----- adjustment. ---------C
SUBROUTINE OUTPUT(Z,P,V,R,T,PMG,CF,FORCE,TIME,ISPP,ISPF,
# ITWORK,ITTAIL,PB,PCON2,DE,E,DT,NZM,NTM,PMAX)
C1---------------------------------------------------------------------C
PARAMETER(PI=3.1415926)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
COMMON /COM/NZ,NW,M
DIMENSION Z(NZM+1),P(NZM+1,NTM),V(NZM+1,NTM),R(NZM+1,NTM),
# T(NZM+1,NTM),FORCE(NTM),PMG(NTM),TIME(NTM),CF(NTM)
REAL IP,ISPP(NTM),ISPF(NTM),MTAV,ISPAV,LIMDA
C5---------------------------------------------------------------------C
OPEN(1,FILE='VERST_I2.DAT',STATUS='UNKNOWN')
C10--------------------------------------------------------------------C
AE=PI*DE**2/4.0
DO 15 I=1,ITTAIL
CALL FAT(TIME(I),AT,E,DT)
EPS=AE/AT
CALL PARA (LIMDA,EPS)
PAIP=1.0-(G-1.)/(G+1.)*LIMDA*LIMDA
CFO=GE*SQRT(2.0*G/(G-1.)*(1.0-PAIP))
CF(I)=(CFO+EPS*(PAIP**(G/(G-1.))-PB/P(M,I)))*PCON2
IF(CF(I).LT.0.0) CF(I)=0.0
FORCE(I)=CF(I)*P(M,I)*AT
15 CONTINUE
C15--------------------------------------------------------------------C
PP=0.0
FF=0.0
PMM=0.0
DO 10 I=1,ITTAIL-1
DTT=TIME(I+1)-TIME(I)
PP=PP+0.5*(P(M,I)+P(M,I+1))*DTT
FF=FF+0.5*(FORCE(I)+FORCE(I+1))*DTT
PMM=PMM+0.5*(PMG(I)+PMG(I+1))*DTT
ISPP(I+1)=PP/PMM
ISPF(I+1)=FF/PMM
10 CONTINUE
C20--------------------------------------------------------------------C
ISPP(1)=0.0
ISPF(1)=0.0
PCAV=PP/TIME(ITTAIL)
FOAV=FF/TIME(ITTAIL)
MTAV=PMM/TIME(ITTAIL)
IP=FF
ISPAV=FOAV/MTAV
C25--------------------------------------------------------------------C
WRITE(1,*)'''PMAX'',',PMAX,','
CLOSE(1)
C30--------------------------------------------------------------------C
C I_PS_Z------- the stated pressure of every node C
C versus Z-corrdinate. C
C I_PE_Z------- the pressure of every node C
C versus Z-corrdinate. C
C the following 'I_TS_Z,I_TE_Z,I_VS_Z,I_VE_Z,I_RS_Z,I_RE_Z' C
C are used like the above example . C
C----------------------------------------------------------------------C
OPEN(1,FILE='I_PS_Z.DAT',STATUS='UNKNOWN')
DO 25 I=1,M
P(I,1)=P(I,1)/1000000.
WRITE(1,*)Z(I),P(I,1)
25 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_PE_Z.DAT',STATUS='UNKNOWN')
DO 30 I=1,M
P(I,ITWORK)=P(I,ITWORK)/1000000.
WRITE(1,*)Z(I),P(I,ITWORK)
30 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TS_Z.DAT',STATUS='UNKNOWN')
DO 35 I=1,M
WRITE(1,*)Z(I),T(I,1)
35 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TE_Z.DAT',STATUS='UNKNOWN')
DO 40 I=1,M
WRITE(1,*)Z(I),T(I,ITWORK)
40 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VS_Z.DAT',STATUS='UNKNOWN')
DO 45 I=1,M
WRITE(1,*)Z(I),V(I,1)
45 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VE_Z.DAT',STATUS='UNKNOWN')
DO 50 I=1,M
WRITE(1,*)Z(I),V(I,ITWORK)
50 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RS_Z.DAT',STATUS='UNKNOWN')
DO 55 I=1,M
WRITE(1,*)Z(I),R(I,1)
55 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RE_Z.DAT',STATUS='UNKNOWN')
DO 60 I=1,M
WRITE(1,*)Z(I),R(I,ITWORK)
60 CONTINUE
CLOSE(1)
C35--------------------------------------------------------------------C
C I_PH_T----- the pressure in the head-end cross-section C
C versus time . C
C I_PE_T----- the pressure in the aft-end cross-section C
C versus time . C
C the following 'I_TH_T,I_TE_T,I_VH_T,I_VE_T,I_RH_T,I_RE_T' C
C are used like the above example. C
C----------------------------------------------------------------------C
OPEN(1,FILE='I_PH_T.DAT',STATUS='UNKNOWN')
P(1,1)=P(1,1)*1000000.
P(1,ITWORK)=P(1,ITWORK)*1000000.
DO 65 I=1,ITTAIL
C------- In file I_PH_T.DAT, pressure unit is MPa.--------------------C
P(1,I)=P(1,I)/1000000.
C----------------------------------------------------------------------C
WRITE(1,*)TIME(I),P(1,I)
65 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_PE_T.DAT',STATUS='UNKNOWN')
P(M,1)=P(M,1)*1000000.
P(M,ITWORK)=P(M,ITWORK)*1000000.
DO 70 I=1,ITTAIL
C------- In file I_PE_T.DAT, pressure unit is MPa.--------------------C
P(M,I)=P(M,I)/1000000.
C----------------------------------------------------------------------C
WRITE(1,*)TIME(I),P(M,I)
70 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TH_T.DAT',STATUS='UNKNOWN')
DO 75 I=1,ITTAIL
WRITE(1,*)TIME(I),T(1,I)
75 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_TE_T.DAT',STATUS='UNKNOWN')
DO 80 I=1,ITTAIL
WRITE(1,*)TIME(I),T(M,I)
80 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VH_T.DAT',STATUS='UNKNOWN')
DO 85 I=1,ITTAIL
WRITE(1,*)TIME(I),V(1,I)
85 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_VE_T.DAT',STATUS='UNKNOWN')
DO 90 I=1,ITTAIL
WRITE(1,*)TIME(I),V(M,I)
90 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RH_T.DAT',STATUS='UNKNOWN')
DO 95 I=1,ITTAIL
WRITE(1,*)TIME(I),R(1,I)
95 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_RE_T.DAT',STATUS='UNKNOWN')
DO 100 I=1,ITTAIL
WRITE(1,*)TIME(I),R(M,I)
100 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_CF_T.DAT',STATUS='UNKNOWN')
DO 105 I=1,ITTAIL
WRITE(1,*)TIME(I),CF(I)
105 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_F_T.DAT',STATUS='UNKNOWN')
DO 110 I=1,ITTAIL
WRITE(1,*)TIME(I),FORCE(I)
110 CONTINUE
CLOSE(1)
OPEN(1,FILE='I_M_T.DAT',STATUS='UNKNOWN')
DO 115 I=1,ITTAIL
WRITE(1,*)TIME(I),PMG(I)
115 CONTINUE
CLOSE(1)
C OPEN(1,FILE='I_ISP_T.DAT',STATUS='UNKNOWN')
C DO 120 I=1,ITTAIL
C WRITE(1,*)TIME(I),ISPP(I)
C120 CONTINUE
C CLOSE(1)
OPEN(1,FILE='I_ISP_T.DAT',STATUS='UNKNOWN')
DO 125 I=1,ITTAIL
WRITE(1,*)TIME(I),ISPF(I)
125 CONTINUE
CLOSE(1)
C40--------------------------------------------------------------------C
OPEN(1,FILE='INTER_O.DAT',STATUS='UNKNOWN')
WRITE(1,*)'The internal ballistic performance:'
WRITE(1,*)'Time node number of work time : ITWORK=',ITWORK
WRITE(1,*)'Work time of rocket engine : ta=',TIME(ITWORK)
WRITE(1,*)'Time node number including work time and tail time :'
WRITE(1,*)' ITTAIL=',ITTAIL
WRITE(1,*)'Sum of work time and tail time :tb=',TIME(ITTAIL)
WRITE(1,*)'Maximum chamber pressure : Pmax=',PMAX
WRITE(1,*)'Average chamber pressure : Pcav=',PCAV
WRITE(1,*)'Average force : Fav=',FOAV
WRITE(1,*)'Average mass flow rate : Mav=',MTAV
WRITE(1,*)'Total impulse : Ip=',IP
WRITE(1,*)'Average specific impulse : Isp=',ISPAV
CLOSE(1)
C45--------------------------------------------------------------------C
RETURN
END
C*********************SUB MS******************************************C
C This subroutine is MERSON METHOD .it is used to solve the C
C---- diffierential equation group. ----------C
SUBROUTINE MS(N,X,Y,EPS,H,FIRST,W,AP,APP,DAP,W0,WW,VT,PT,TT,
# NZM,NWM)
C1---------------------------------------------------------------------C
LOGICAL INC,FIRST
DIMENSION Y(N),W(N,5),AP(NZM+1,NWM+1),APP(NZM+1),DAP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),W0(NWM+1),WW(NZM+1)
C5---------------------------------------------------------------------C
IF(.NOT.FIRST) GO TO 2
LOCP=1
FIRST=.FALSE.
2 HC=H/FLOAT(LOCP)
LOC=0
4 CALL FCT(X,Y,W(1,3),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
DO 6 I=1,N
W(I,1)=W(I,3)*HC/3.0+Y(I)
6 CONTINUE
CALL FCT(HC/3.0+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 8 I=1,N
W(I,1)=(W(I,3)+W(I,4))*HC/6.0+Y(I)
8 CONTINUE
CALL FCT(HC/3.0+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 10 I=1,N
W(I,1)=(W(I,3)+3.0*W(I,4))*HC*0.125+Y(I)
10 CONTINUE
CALL FCT(HC*0.5+X,W(1,1),W(1,5),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,
# NWM)
DO 12 I=1,N
W(I,1)=W(I,3)*HC*0.5-W(I,4)*HC*1.5+W(I,5)*HC*2.0+Y(I)
12 CONTINUE
CALL FCT(HC+X,W(1,1),W(1,4),AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
DO 14 I=1,N
W(I,2)=(W(I,3)+W(I,4)+W(I,5)*4.0)*HC/6.0+Y(I)
14 CONTINUE
C10--------------------------------------------------------------------C
INC=.TRUE.
DO 18 I=1,N
ER=ABS(W(I,1)-W(I,2))*0.2
IF(ABS(W(I,2)).GT.1.0) ER=ER/ABS(W(I,2))
IF(ER.LT.EPS) GO TO 18
HC=HC*0.5
LOCP=LOCP+LOCP
LOC=LOC+LOC
GO TO 4
18 IF(64.0*ER.GT.EPS) INC=.FALSE.
C15--------------------------------------------------------------------C
X=X+HC
DO 20 I=1,N
Y(I)=W(I,2)
20 CONTINUE
LOC=LOC+1
IF(LOC.EQ.LOCP) RETURN
C20--------------------------------------------------------------------C
IF(.NOT.INC.OR.LOC.NE.LOC/2*2) GO TO 4
C25--------------------------------------------------------------------C
HC=2.0*HC
LOC=LOC/2
LOCP=LOCP/2
GO TO 4
C30--------------------------------------------------------------------C
END
C*******************SUB FCT*******************************************C
C This subroutine is used to calculate the values of all the C
C---- right of the equations in the differential equation groups.---C
SUBROUTINE FCT(X,Y,YX,AP,APP,DAP,W0,WW,VT,PT,TT,NZM,NWM)
C1--------------------------------------------------------------------C
COMMON /FCDA/I,H,DTIME,IT
COMMON G,GE,AJ,RG,RP,TS,AS,AT
DIMENSION Y(3),YX(3),AP(NZM+1,NWM+1),APP(NZM+1),DAP(NZM+1),
# VT(NZM+1),PT(NZM+1),TT(NZM+1),W0(NWM+1),WW(NZM+1)
C5--------------------------------------------------------------------C
K=I
IF(I.EQ.2) K=I+1
K2=K-2
K1=K-1
Q=X/H-FLOAT(K2)
H2=H
CALL FLOAD2(Y,R,PS,TSG,A2)
CALL BURN(Y(2),BR)
C------ CALL BURN(Y(2),BR,Y(1))----------C
C10-------------------------------------------------------------------C
DO 73 II=1,K
CALL GEOM(II, DAP(II)=APP1-APP(II)
73 CONTINUE
C15--------------------------------------------------------------------C
APT=FX(APP(K2),APP(K1),APP(K),Q,H,H2)
DAPP=FX(DAP(K2),DAP(K1),DAP(K),Q,H,H2)
DAPT=DFX(APP(K2),APP(K1),APP(K),Q,H,H2)
VTT=FX(VT(K2),VT(K1),VT(K),Q,H,H2)
PTT=FX(PT(K2),PT(K1),PT(K),Q,H,H2)
TTT=FX(TT(K2),TT(K1),TT(K),Q,H,H2)
A=((RP-R)*DAPP-R*Y(1)*DAPT)/APT
B=-RP*Y(1)*DAPP/(APT*R)
C=(DAPP*(RP*G*RG*TS+0.5*(G-1.)*RP*Y(1)*Y(1)-Y(2))
# -G*Y(2)*Y(1)*DAPT)/APT
DV=A2-Y(1)*Y(1)
YX(1)=(Y(1)*(VTT-B)+(C-PTT)/R)/DV
YX(2)=(G*Y(2)*(B-VTT)+Y(1)*(PTT-C))/DV
YX(3)=Y(3)*(C-RG*Y(3)*A)/(Y(1)*Y(2))-(G-1.)*Y(3)/Y(1)*YX(1)-
# TTT/Y(1)
C20-------------------------------------------------------------------C
RETURN
END
C*******************SUB FLOAD2*****************************************C
C This subroutine is used to calculate some parameters such as: C
C---- combustion gas density,specific heat ratio,etc .---------------C
SUBROUTINE FLOAD2(Y,R,PS,TSG,A2)
DIMENSION Y(3)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
REAL LIMDA
R=Y(2)/(Y(3)*RG)
CP=G*RG/(G-1.)
TSG=Y(3)+Y(1)**2/(2.*CP)
AC=SQRT(2.*G*RG*TSG/(G+1.))
LIMDA=Y(1)/AC
A=1.0-LIMDA*LIMDA*(G-1.0)/(G+1.0)
PS=Y(2)*A**(G/(1.0-G))
A2=G*RG*Y(3)
RETURN
END
C****************SUB BURN*********************************************C
C This subroutine is used to calculate the propellant burning C
C----- rate. ---------C
SUBROUTINE BURN(P,BR)
C------- SUBROUTINE BURN(P,BR,V) ---------C
COMMON /BDA/BA,BN,TI,AI
C=EXP(AI*(TI-20.0))
C------- PARAMETER(VTH=*****,VKU=****) ---C
C------- IF(V.LE.VTH) THEN ----------C
C------- BR=BA*C*P**BN ---------------C
C------- ELSE --------------------------C
C------- BR=(1.0+VKU*(V-VTH))*BA*C*P**BN
C------- ENDIF --------------------------C
BR=BA*C*P**BN
RETURN
END
C***********************FUN FX*****************************************C
C----- This function 'FX' is the Lagrange interpolation . --------C
FUNCTION FX(YO,Y1,Y2,Q,H1,H2)
FX=(Q*Q*H1*H1*(H2*YO-(H1+H2)*Y1+H1*Y2)-Q*H1*(H2*H2*YO+(H1*H1
# -H2*H2)*Y1-H1*H1*Y2)+H1*H2*(H1+H2)*Y1)/(H1*H2*(H1+H2))
RETURN
END
C**********************FUN DFX*****************************************C
C---- This function DFX is a differential function -------C
FUNCTION DFX(YO,Y1,Y2,Q,H1,H2)
DFX=(2.0*Q*H1*(H2*YO-(H1+H2)*Y1+H1*Y2)-(H2*H2*YO+
# (H1*H1-H2*H2)*Y1-H1*H1*Y2))/(H1*H2*(H1+H2))
RETURN
END
C********************SUB FAT*******************************************C
C This subroutine 'FAT' is used to calculate the nozzle throat C
C----- area under the conditions of nozzle erosion or unerosion -----C
SUBROUTINE FAT(TIME,AT,E,DT)
PARAMETER(PI=3.1415926)
AT=PI*DT**2/4.
IF(TIME.GE.15) AT=PI*(DT+2*E*(TIME-15))**2/4.
RETURN
END
C***********************SUB AHEAD**************************************C
C This subroutine is used to calculate the head-end burning area C
C----- according to the propellant grain geometry calculations. ----C
SUBROUTINE AHEAD(WW,W0,ABH,AH,NWM)
DIMENSION W0(NWM+1),ABH(NWM+1)
J=1
10 J1=J+1
IF(WW.GE.W0(J).AND.WW.LE.W0(J1)) GOTO 20
J=J1
GOTO 10
20 AH=FLINE(ABH(J),ABH(J1),W0(J),W0(J1),WW)
RETURN
END
C**********************FUN FLINE***************************************C
C----- Use the linear interpolation to interpolate the given nodes.---C
FUNCTION FLINE(Y1,Y2,X1,X2,X)
FLINE=Y1+(Y2-Y1)*(X-X1)/(X2-X1)
RETURN
END
C*********************SUB GEOM*****************************************C
C This subroutine calls the interpolation function subroutine FX C
C to interpolate the propellant grain geometry,then gets the C
C value of the area's averge ratio of change to the time,in C
C---- every time node. ----------C
SUBROUTINE GEOM(I,WW,APP,W0,AP,NZM,NWM)
COMMON /COM/NZ,NW,M
DIMENSION W0(NWM+1),AP(NZM+1,NWM+1)
IF(WW.GE.W0(NW+1)) THEN
APP=AP(I,NW+1)+(AP(I,NW+1)-AP(I,NW))/(W0(NW+1)-
# W0(NW))*(WW-W0(NW+1))
RETURN
ENDIF
J=1
10 JP1=J+1
IF(WW.GE.W0(J).AND.WW.LE.W0(JP1)) THEN
APP=FLINE(AP(I,J),AP(I,JP1),W0(J),W0(JP1),WW)
ELSE
J=JP1
GOTO 10
ENDIF
RETURN
END
C*******************SUB TAIL******************************************C
C This subroutine is used to calculate and output these para- C
C----- meters' values in the tail phase. ---------------C
SUBROUTINE TAIL(IT,P,R,T,V,AP,PMG,TIME,ITWORK,DTTAIL,VCH,
# NZM,NWM,NTM,PMAX)
C1---------------------------------------------------------------------C
COMMON /COM/NZ,NW,M
DIMENSION P(NZM+1,NTM),R(NZM+1,NTM),T(NZM+1,NTM),V(NZM+1,NTM),
# PMG(NTM),TIME(NTM),AP(NZM+1,NWM+1)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
C5---------------------------------------------------------------------C
TIM=TIME(IT)
TT=0.0
A=2.0*VCH
A1=GE*SQRT(RG*TS)*AT*(G-1.0)
CD=GE/SQRT(RG*TS)
C=2.0*G/(G-1.0)
D=0.1*PMAX
C10--------------------------------------------------------------------C
10 TIM=TIM+DTTAIL
TT=TT+DTTAIL
IT=IT+1
B=A+A1*TT
P(1,IT)=P(1,ITWORK)*(A/B)**C
R(1,IT)=R(1,ITWORK)*(P(1,IT)/P(1,ITWORK))**(1./G)
T(1,IT)=P(1,IT)/(RG*R(1,IT))
V(1,IT)=0.0
P(M,IT)=P(M,ITWORK)*(A/B)**C
R(M,IT)=R(M,ITWORK)*(P(M,IT)/P(M,ITWORK))**(1./G)
T(M,IT)=P(M,IT)/(RG*R(M,IT))
PMG(IT)=CD*P(M,IT)*AT
V(M,IT)=PMG(IT)/(R(M,IT)*AP(M,NW+1))
TIME(IT)=TIM
C15--------------------------------------------------------------------C
IF(P(M,IT).LE.D) GOTO 20
GOTO 10
C20--------------------------------------------------------------------C
20 RETURN
END
C*******************SUB PARA*******************************************C
C This subroutine is called to calculate the velocity coefficientC
C in the nozzle exit cross-section by calling the subroutine C
C----- SMACH which is used to calculate it . -------------C
SUBROUTINE PARA(LIMDA,EPS)
COMMON /S/GP1,GM1,SA,SB,SC
REAL M1,M2,LIMDA
M1=1.01
M2=4.0
CALL SMACH(EPS,EM,M1,M2)
EM2=EM*EM
LIMDA=SQRT(0.5*GP1*EM2/(1.0+GM1*0.5*EM2))
RETURN
END
C*********************SUB SMACH****************************************C
C This subroutine is used to calculate the mach number in the C
C----- nozzle exit cross-section. ---------C
SUBROUTINE SMACH(E,M,M1,M2)
COMMON /S/GP1,GM1,SA,SB,SC
REAL M,M1,M2,MO
FUN(X,SA,SB,SC)=SA*(1.0+SB*X*X)**SC/X
F2=FUN(M2,SA,SB,SC)
10 F10=FUN(M1,SA,SB,SC)
T=(F2-F10)/(M2-M1)
MO=M1+(E-F10)/T
IF(ABS(MO-M1).LE.1E-5) GO TO 20
M1=MO
GO TO 10
20 M=MO
RETURN
END
C*********************SUB FINE*****************************************C
C This subroutine calls the function FVC,and uses the Two-cuttingC
C---- Method to calculate the solution of the given function .-------C
SUBROUTINE FINE(A,B,F,AEPS,LIMDA)
LOGICAL FAG
REAL LIMDA
FAG=F(A).GT.0.0
5 C=0.5*(A+B)
IF((B-A).GE.AEPS) GO TO 20
10 LIMDA=C
RETURN
20 IF(F(C)) 25,10,35
25 IF(FAG) GO TO 40
30 A=C
GO TO 5
35 IF(FAG) GO TO 30
40 B=C
GO TO 5
END
C*********************FUN SIMP*****************************************C
C---- Simpson integration method . ---------C
FUNCTION SIMP(A,B,F,EPS,HO)
H=B-A
RP=F(A)+F(B)
N=1
2 X=A-0.5*H
RC=0.0
DO 4 I=1,N
X=X+H
RC=F(X)+RC
4 CONTINUE
R2=(4.0*RC+RP)*H/6.0
IF(N.EQ.1.OR.ABS(H).GT.HO) GOTO 6
X=R2-R1
IF(ABS(R2).GE.1.0) X=X/R2
IF(ABS(X).LT.EPS) GO TO 8
6 R1=R2
RP=2.0*RC+RP
N=N+N
H=0.5*H
GO TO 2
8 SIMP=R2
RETURN
END
C********************FUN F1********************************************C
C It is used to as an extenal function of Simpson integration . C
C----- its value is returned by F1. -------------C
FUNCTION F1(LIMDA)
COMMON /BDA/BA,BN,TI,AI
COMMON G,GE,AJ,RG,RP,TS,AS,AT
REAL LIMDA
A=(1.0+LIMDA*LIMDA)**(2.0-BN)
B=(1.0-LIMDA*LIMDA*(G-1.0)/(G+1.0))**BN
C=1.0
F1=4.0*(1.0-LIMDA*LIMDA)/(A*B*C)
RETURN
END
C******************FUN FVC*********************************************C
C This function helps subroutine FINE to solve the velocity coe- C
C---- fficient in the combustion chamber exit cross-section . -------C
FUNCTION FVC(X)
COMMON G,GE,AJ,RG,RP,TS,AS,AT
FVC=AJ-X*(1.0-X*X*(G-1.0)/(G+1.0))**(1.0/(G-1.0))*
# ((G+1.0)/2.0)**(1.0/(G-1.0))
RETURN
END
- Status
Similar threads
- Locked
- Question
