[Fortran] 纯文本查看 复制代码
PROGRAM NGDE
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER (MAX=42)
PARAMETER (MAX1=42)
DIMENSION V(MAX1+1),XN(MAX1),X(MAX1,MAX1,MAX1),COLF(MAX1,MAX1), &
XN1S(MAX1),ZETA(MAX1),DP(MAX1),XM(MAX1)
PI=3.1415926
XKB=1.38D-23
R=8.314
XNA=6.02225D23
T=1773.D0
COOLRATE=1000.D0
P=1.
XMW=26.95D-3
RHO=2700.D0
A=948.D0
B=0.202
C=13.07
D=36373.D0
AMU=1.966D-6
BMU=147.47
Q=10.**(12./(MAX-2))
V1=XMW/(RHO*XNA)
V(1)=V1
DO 100 I=1,MAX-1
XN(I)=0.
XN1S(I)=0.
100 CONTINUE
DO 101 I=2,MAX
V(I)=V(1)*(Q**(I-1))
DP(I)=(6.*V(I)/PI)**0.333333
XM(I)=RHO*V(I)
101 CONTINUE
DP(1)=(6.*V(1)/PI)**0.333333
XM(1)=RHO*V(1)
DO 102 K=2,MAX-1
DO 102 I=2,MAX-1
DO 102 J=2,MAX-1
IF (V(K) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K+1)) THEN
X(I,J,K)=(V(K+1)-V(I)-V(J))/(V(K+1)-V(K))
ELSE IF (V(K-1) .LE. (V(I)+V(J)) .AND. (V(I)+V(J)) .LT. V(K)) THEN
X(I,J,K)=(V(I)+V(J)-V(K-1))/(V(K)-V(K-1))
ELSE
X(I,J,K)=0.D0
END IF
102 CONTINUE
TIME=0.D0
STEP=1.D-4
S=1.001
TEMP1=6.*V1/PI
S1=PI*(TEMP1**0.66666666667)
XMASS=RHO*V1
PS=EXP(C-D/T)*101325.*P
XNS=PS/(XKB*T)
XN(1)=S*XNS
200 FORMAT( 'Time',5X,'Temperature',5X,'XN(1)',5X,'XJK',5X,'S',5X,'Diameter',5X,'Particle volume',5X,'Total number')
WRITE(*,200)
ICOUNTER=1
ICOUNTER2=1
DO WHILE (T .GT. 1600)
SIGMA=(A-B*T)*1.D-3
PS=EXP(C-D/T)*101325.*P
XNS=PS/(XKB*T)
S=XN(1)/XNS
THETA=(S1*SIGMA)/(XKB*T)
AA=(2.*SIGMA)/(PI*XMASS)
BB=THETA-(4.*(THETA**3.))/(27.*(DLOG(DABS(S))**2.))
XJK=(XNS**2.)*S*V1*(AA**0.5)*EXP(BB)
CC=0.6666667*THETA/DLOG(DABS(S))
XKSTAR=CC**3
DPSTAR=4.*SIGMA*V1/(XKB*T*DLOG(DABS(S)))
VSTAR=PI*(DPSTAR**3.)/6.
XNTOT=0.D0
DO I=2,MAX-1
XNTOT=XNTOT+XN(I)
END DO
VTOT=0.D0
VTOTAL=0.D0
DO I=2,MAX-1
VTOT=VTOT+XN(I)*V(I)
END DO
VTOTAL=VTOT+XN(1)*V(1)
VAV=VTOT/XNTOT
DPAV=(6.*VAV/PI)**0.3333333
IF (XNTOT .GT. 100.) STEP=1.D-5
IF (TIME .GT. 0.14) STEP=5.D-8
201 FORMAT(8(1X,E14.8))
202 FORMAT(8(1X,E14.8))
203 FORMAT(/'time=',E14.8)
204 FORMAT('N',I2,2(2X,E14.8))
OPEN(11,FILE='output1.DAT',STATUS='UNKNOWN',ACCESS='APPEND')
OPEN(12,FILE='output2.DAT',STATUS='UNKNOWN',ACCESS='APPEND')
ICOUNTER=ICOUNTER+1
ICOUNTER2=ICOUNTER2+1
IF (ICOUNTER .EQ. 400) THEN
WRITE(*,201)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
WRITE(11,202)TIME,T,XN(1),XJK,S,DPAV,VTOT,XNTOT
ICOUNTER=1
END IF
IF (ICOUNTER2 .EQ. 40000) THEN
WRITE(*,203)TIME
WRITE(12,203)TIME
DO I=1,MAX-1
WRITE(*,204)I,DP(I),XN(I)
WRITE(12,204)I,DP(I),XN(I)
END DO
ICOUNTER2=1
END IF
XKMIN=1.D-9
XMU=AMU*(T**1.5)/(BMU+T)
XLAMBDA=(XMU/(P*101325.))*SQRT(PI*R*T/(2.*0.04))
DO 103 I=1,MAX-1
DO 103 J=1,MAX-1
XKN1=(2.*XLAMBDA)/DP(I)
XKN2=(2.*XLAMBDA)/DP(J)
D1=(XKB*T)/(3.*PI*XMU*DP(I))*((5.+4.*XKN1+6.*(XKN1**2)+18.*(XKN1**3))/(5.-XKN1+(8.+PI)*(XKN1**2)))
D2=(XKB*T)/(3.*PI*XMU*DP(J))*((5.+4.*XKN2+6.*(XKN2**2)+18.*(XKN2**3))/(5.-XKN2+(8.+PI)*(XKN2**2)))
C1=SQRT((8.*XKB*T)/(PI*XM(I)))
C2=SQRT((8.*XKB*T)/(PI*XM(J)))
XL1=(8.*D1)/(PI*C1)
XL2=(8.*D2)/(PI*C2)
G1=((DP(I)+XL1)**3)-((DP(I)**2+XL1**2)**1.5)/(3.*DP(I)*XL1)-DP(I)
G2=((DP(J)+XL2)**3)-((DP(J)**2+XL2**2)**1.5)/(3.*DP(J)*XL2)-DP(J)
COLF(I,J)=2.*PI*(D1+D2)*(DP(I)+DP(J))/((DP(I)+DP(J))/(DP(I)+DP(J)+2.*SQRT(G1**2+G2**2))+(8.*(D1+D2))/(SQRT(C1**2+C2**2)*(DP(I)+DP(J))))
IF (COLF(I,J) .LT. XKMIN) THEN
XKMIN=COLF(I,J)
END IF
103 CONTINUE
DO 104 K=2,MAX-1
IF (VSTAR .LT. V1) THEN
ZETA(2)=VSTAR/V(2)
ELSE IF (V(K-1) .LE. VSTAR .AND. VSTAR .LT. V(K)) THEN
ZETA(K)=VSTAR/V(K)
ELSE
ZETA(K)=0.D0
END IF
104 CONTINUE
DO 105 I=1,MAX-1
DP(I)=(6.*V(I)/PI)**0.333333
XN1S(I)=XNS*EXP(4.*SIGMA*XMW/(R*T*RHO*DP(I)))
105 CONTINUE
T1=0.D0
T2=0.D0
T3=0.D0
T4=0.D0
DO 106 K=2,MAX-1
SUM1=0.D0
SUM2=0.D0
IF (STEP .NE. 1.D-4) THEN
IF (XN(1) .GT. XN1S(K-1)) THEN
IF (K .EQ. 2) THEN
ADDTERM=0.D0
ELSE
ADDTERM=(V(1)/(V(K)-V(K-1)))*COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
T1=T1+COLF(1,K-1)*(XN(1)-XN1S(K-1))*XN(K-1)
END IF
END IF
IF (XN(1) .LT. XN1S(K+1)) THEN
ADDTERM=-(V(1)/(V(K+1)-V(K)))*COLF(1,K+1)*(XN(1)-XN1S(K+1))*XN(K+1)
T2=T2+COLF(1,K+1)*(-XN(1)+XN1S(K+1))*XN(K+1)
END IF
IF (XN(1) .LT. XN1S(K)) THEN
SUBTERM=-(V(1)/(V(K)-V(K-1)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
T3=T3+COLF(1,K)*(-XN(1)+XN1S(K))*XN(K)
END IF
IF (XN(1) .GT. XN1S(K)) THEN
SUBTERM=(V(1)/(V(K+1)-V(K)))*COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
T4=T4+COLF(1,K)*(XN(1)-XN1S(K))*XN(K)
END IF
END IF
DO I=2,MAX-1
SUM2=SUM2+COLF(K,I)*XN(I)
DO J=2,K
SUM1=SUM1+X(I,J,K)*COLF(I,J)*XN(I)*XN(J)
END DO
END DO
XN(K)=XN(K)+STEP*(0.5*SUM1-XN(K)*SUM2+XJK*ZETA(K)+ADDTERM-SUBTERM)
106 CONTINUE
XN(1)=XN(1)-XJK*XKSTAR*STEP-STEP*(T1+T4-T3-T2)
T=T-STEP*COOLRATE
TIME=TIME+STEP
CLOSE(11)
CLOSE(12)
END DO
STOP
END
目前程式已經可以跑了, 我還在繼續測試結果是否合理