[Fortran] 纯文本查看 复制代码
PROGRAM POINTEHL
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
DATA PAI,Z/3.14159265,0.68/
DATA N,PH,E1,EDA0,RX,US,X0,XE/33,0.8E9,2.21E11,0.05,0.02,1.0,-2.5,1.5/
OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
OPEN(10,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
MM=N-1
A1=ALOG(EDA0)+9.67
A2=5.1E-9*PH
A3=0.59/(PH*1.E-9)
U=EDA0*US/(2.*E1*RX)
B=PAI*PH*RX/E1
W0=2.*PAI*PH/(3.*E1)*(B/RX)**2
ALFA=Z*5.1E-9*A1
G=ALFA*E1
HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W0**(-0.073)
ENDA=12.*U*(E1/PH)*(RX/B)**3
WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
WRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,US
WRITE(*,*)' Wait please'
CALL SUBAK(MM)
CALL EHL(N,X0,XE)
STOP
END
SUBROUTINE EHL(N,X0,XE)
DIMENSION X(65),Y(65),H(4500),RO(4500),EPS(4500),EDA(4500),P(4500),POLD(4500),V(4500)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
DATA MK,G0/1,2.0943951/
CALL INITI(N,DX,X0,XE,X,Y,P,POLD)
KK=0
CALL HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
14 KK=15
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
MK=MK+1
CALL ERP(N,ER,P,POLD)
WRITE(*,*)'ER=',ER
IF(ER.GT.1.E-5.AND.DH.GT.1.E-6)THEN
IF(MK.GE.20)THEN
MK=1
DH=0.5*DH
ENDIF
GOTO 14
ENDIF
IF(DH.LE.1.E-6)WRITE(*,*)'Pressures are not convergent!!!'
CALL OUTPUT(N,DX,X,Y,H,P)
RETURN
END
SUBROUTINE ERP(N,ER,P,POLD)
DIMENSION P(N,N),POLD(N,N)
ER=0.0
SUM=0.0
DO 10 I=1,N
DO 10 J=1,N
ER=ER+ABS(P(I,J)-POLD(I,J))
POLD(I,J)=P(I,J)
SUM=SUM+P(I,J)
10 CONTINUE
ER=ER/SUM
RETURN
END
SUBROUTINE INITI(N,DX,X0,XE,X,Y,P,POLD)
DIMENSION X(N),Y(N),P(N,N),POLD(N,N)
NN=(N+1)/2
DX=(XE-X0)/(N-1.)
Y0=-0.5*(XE-X0)
DO 5 I=1,N
X(I)=X0+(I-1)*DX
Y(I)=Y0+(I-1)*DX
5 CONTINUE
DO I=1,N
D=1.-X(I)*X(I)
DO J=1,NN
C=D-Y(J)*Y(J)
IF(C.LE.0.0)P(I,J)=0.0
IF(C.GT.0.0)P(I,J)=SQRT(C)
POLD(I,J)=P(I,J)
ENDDO
ENDDO
RETURN
END
SUBROUTINE HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COMAK/AK(0:65,0:65)
DATA PAI,PAI1/3.14159265,0.2026423/
NN=(N+1)/2
CALL VI(N,DX,P,V)
HMIN=1.E3
DO 30 I=1,N
DO 30 J=1,NN
RAD=X(I)*X(I)+Y(J)*Y(J)
W1=0.5*RAD
H0=W1+V(I,J)
IF(H0.LT.HMIN)HMIN=H0
30 H(I,J)=H0
IF(KK.EQ.0)THEN
KK=1
DH=0.01*HM0
H00=-HMIN+HM0
ENDIF
W1=0.0
DO 32 I=1,N
DO 32 J=1,N
32 W1=W1+P(I,J)
W1=DX*DX*W1/G0
DW=1.-W1
IF(DW.LT.0.0)H00=H00+DH
IF(DW.GT.0.0)H00=H00-DH
DO 60 I=1,N
DO 60 J=1,NN
H(I,J)=H00+H(I,J)
EDA1=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
EDA(I,J)=EDA1
RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60 EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA1)
DO 70 J=NN+1,N
JJ=N-J+1
DO 70 I=1,N
H(I,J)=H(I,JJ)
RO(I,J)=RO(I,JJ)
EDA(I,J)=EDA(I,JJ)
70 EPS(I,J)=EPS(I,JJ)
RETURN
END
SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)
COMMON /COMAK/AK(0:65,0:65)
DATA KG1,PAI/0,3.14159265/
IF(KG1.NE.0)GOTO 2
KG1=1
AK00=AK(0,0)
AK10=AK(1,0)
2 NN=(N+1)/2
DO 100 K=1,KK
DO 70 J=2,NN
J0=J-1
J1=J+1
D2=0.5*(EPS(1,J)+EPS(2,J))
DO 70 I=2,N-1
I0=I-1
I1=I+1
D1=D2
D2=0.5*(EPS(I1,J)+EPS(I,J))
D4=0.5*(EPS(I,J0)+EPS(I,J))
D5=0.5*(EPS(I,J1)+EPS(I,J))
D8=2.0*RO(I,J)*AK00/PAI**2
D9=2.0*RO(I0,J)*AK10/PAI**2
D10=D1+D2+D4+D5+D8*DX-D9*DX
D11=D1*P(I0,J)+D2*P(I1,J)+D4*P(I,J0)+D5*P(I,J1)
D12=(RO(I,J)*H(I,J)-D8*P(I,J)-RO(I0,J)*H(I0,J)+D9*P(I,J))*DX
P(I,J)=(D11-D12)/D10
IF(P(I,J).LT.0.0)P(I,J)=0.0
70 CONTINUE
DO 80 J=1,NN
JJ=N+1-J
DO 80 I=1,N
80 P(I,JJ)=P(I,J)
CALL HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
100 CONTINUE
RETURN
END
SUBROUTINE VI(N,DX,P,V)
DIMENSION P(N,N),V(N,N)
COMMON /COMAK/AK(0:65,0:65)
PAI1=0.2026423
DO 40 I=1,N
DO 40 J=1,N
H0=0.0
DO 30 K=1,N
IK=IABS(I-K)
DO 30 L=1,N
JL=IABS(J-L)
30 H0=H0+AK(IK,JL)*P(K,L)
40 V(I,J)=H0*DX*PAI1
RETURN
END
SUBROUTINE SUBAK(MM)
COMMON /COMAK/AK(0:65,0:65)
S(X,Y)=X+SQRT(X**2+Y**2)
DO 10 I=0,MM
XP=I+0.5
XM=I-0.5
DO 10 J=0,I
YP=J+0.5
YM=J-0.5
A1=S(YP,XP)/S(YM,XP)
A2=S(XM,YM)/S(XP,YM)
A3=S(YM,XM)/S(YP,XM)
A4=S(XP,YP)/S(XM,YP)
AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
10 AK(J,I)=AK(I,J)
RETURN
END
SUBROUTINE OUTPUT(N,DX,X,Y,H,P)
DIMENSION X(N),Y(N),H(N,N),P(N,N)
A=0.0
WRITE(8,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(8,110)X(I),(H(I,J),J=1,N)
ENDDO
WRITE(10,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(10,110)X(I),(P(I,J),J=1,N)
ENDDO
110 FORMAT(66(E12.6,1X))
RETURN
END