弹流润滑数值计算程序迭代失败！！求助！！

小戈 · 发表于 2019-4-18 14:59:38

如题，程序出自《弹性流体动压润滑方法与程序》，运行失败，程序如下。求助出错原因

[Fortran] 纯文本查看 复制代码

001        PROGRAM LINEEHL
002                COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COM4/X0,XE/COM3/E1,PH,B,R
003                DATA PAI,Z,P0/3.14159265,0.68,1.96E8/
004                DATA N,X0,XE,W,E1,EDA0,R,Us/130,-4.0,1.5,1.0E5,2.2E11,0.05,0.05,1.5/
005                OPEN(8,FILE='OUT.DAT',STATUS='UNKNOWN')
006                W1=W/(E1*R)
007                PH=E1*SQRT(0.5*W1/PAI)
008                A1=(ALOG(EDA0)+9.67)
009                A2=PH/P0
010                A3=0.59/(PH*1.E-9)
011                B=4.*R*PH/E1
012                ALFA=Z*A1/P0
013                G=ALFA*E1
014                U=EDA0*US/(2.*E1*R)
015                CC1=SQRT(2.*U)
016                AM=2.*PAI*(PH/E1)**2/CC1
017                ENDA=3.*(PAI/AM)**2/8.
018                HM0=1.6*(R/B)**2*G**0.6*U**0.7*W1**(-0.13)
019                WRITE(*,*)N,X0,XE,W,E1,EDA0,R,US
020                CALL SUBAK(N)
021                CALL EHL(N)
022                STOP
023                END
024                SUBROUTINE EHL(N)
025                DIMENSION X(1100),P(1100),H(1100),RO(1100),POLD(1100),EPS(1100),EDA(1100),V(1100)
026                COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XE
027                COMMON /COM3/E1,PH,B,RR
028                MK=1
029                DX=(XE-X0)/(N-1.0)
030                DO 10 I=1,N
031                X(I)=X0+(I-1)*DX
032                IF(ABS(X(I)).GE.1.0)P(I)=0.0
033                IF(ABS(X(I)).LT.1.0)P(I)=SQRT(1.-X(I)*X(I))
03410      CONTINUE
035                CALL HREE(N,DX,X,P,H,RO,EPS,EDA,V)
036                CALL FZ(N,P,POLD)
03714                KK=19
038                CALL ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
039                MK=MK+1
040                CALL ERROP(N,P,POLD,ERP)
041                WRITE(*,*)'ERP=',ERP
042                IF(ERP.GT.1.E-5.AND.DH.GT.1.E-7)THEN
043                IF(MK.GE.50)THEN
044                MK=1
045                DH=0.5*DH
046                ENDIF
047                GOTO 14
048                ENDIF
049                IF(DH.LE.1.E-7)WRITE(*,*)'Pressures are not convergent!!!'
050                 H2=1.E3
051                P2=0.0
052                DO 106 I=1,N
053                IF(H(I).LT.H2)H2=H(I)
054                IF(P(I).GT.P2)P2=P(I)
055106     CONTINUE
056                H3=H2*B*B/RR
057                P3=P2*PH
058110                FORMAT(6(1X,E12.6))
059120                CONTINUE
060                WRITE(*,*)'P2,H2,P3,H3=',P2,H2,P3,H3
061                CALL OUTHP(N,X,P,H)
062                RETURN
063                END
064                SUBROUTINE OUTHP(N,X,P,H)
065                DIMENSION X(N),P(N),H(N)
066                DO 10 I=1,N
067                WRITE(8,20)X(I),P(I),H(I)
06810      CONTINUE
06920      FORMAT(1X,6(E12.6,1X))
070                RETURN
071                END
072                SUBROUTINE HREE(N,DX,X,P,H,RO,EPS,EDA,V)
073                DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
074                COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COMAK/AK(0:1100)
075                DATA KK,PAI1,G0/0,0.318309886,1.570796325/
076                IF(KK.NE.0)GOTO 3
077                H00=0.0
0783                W1=0.0
079                DO 4 I=1,N
0804       W1=W1+P(I)
081                C3=(DX*W1)/G0
082                DW=1.-C3
083                CALL VI(N,DX,P,V)
084                HMIN=1.E3
085                DO 30 I=1,N
086                H0=0.5*X(I)*X(I)+V(I)
087                IF(H0.LT.HMIN)HMIN=H0
088                H(I)=H0
08930                CONTINUE
090                IF(KK.NE.0)GOTO 32
091                KK=1
092                DH=0.005*HM0
093                H00=-HMIN+HM0
09432                IF(DW.LT.0.0)H00=H00+DH
095                IF(DW.GT.0.0)H00=H00-DH
096                DO 60 I=1,N
097                H(I)=H00+H(I)
098                EDA(I)=EXP(A1*(-1.+(1.+A2*P(I))**Z))
099                RO(I)=(A3+1.35*P(I))/(A3+P(I))
100                EPS(I)=RO(I)*H(I)**3/(ENDA*EDA(I))
10160                CONTINUE
102                RETURN
103                END
104                SUBROUTINE ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
105                DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
106                COMMON /COMAK/AK(0:1100)
107                DATA PAI1/0.318309886/
108                DO 100 K=1,KK
109                D2=0.5*(EPS(1)+EPS(2))
110                D3=0.5*(EPS(2)+EPS(3))
111                DO 70 I=2,N-1
112                D1=D2
113                D2=D3
114                IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2))
115                D8=RO(I)*AK(0)*PAI1
116                D9=RO(I-1)*AK(1)*PAI1
117                D10=1.0/(D1+D2+(D9-D8)*DX)
118                D11=D1*P(I-1)+D2*P(I+1)
119                D12=(RO(I)*H(I)-RO(I-1)*H(I-1)+(D8-D9)*P(I))*DX
120                P(I)=(D11-D12)*D10
121                IF(P(I).LT.0.0)P(I)=0.0
12270      CONTINUE
123                CALL HREE(N,DX,X,P,H,RO,EPS,EDA,V)
124100     CONTINUE
125                RETURN
126                END
127                SUBROUTINE VI(N,DX,P,V)
128                DIMENSION P(N),V(N)
129                COMMON /COMAK/AK(0:1100)
130                PAI1=0.318309886
131                C=ALOG(DX)
132                DO 10 I=1,N
133                V(I)=0.0
134                DO 10 J=1,N
135                IJ=IABS(I-J)
13610                V(I)=V(I)+(AK(IJ)+C)*DX*P(J)
137                DO I=1,N
138                V(I)=-PAI1*V(I)
139                ENDDO
140                RETURN
141                END
142 
143                SUBROUTINE SUBAK(MM)
144                COMMON /COMAK/AK(0:1100)
145                DO 10 I=0,MM
14610      AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.)
147                RETURN
148                END
149 
150                SUBROUTINE FZ(N,P,POLD)
151                DIMENSION P(N),POLD(N)
152                DO 10 I=1,N
15310      POLD(I)=P(I)
154                RETURN
155                END
156 
157                SUBROUTINE ERROP(N,P,POLD,ERP)
158                DIMENSION P(N),POLD(N)
159                SD=0.0
160                SUM=0.0
161                DO 10 I=1,N
162                SD=SD+ABS(P(I)-POLD(I))
163                POLD(I)=P(I)
16410      SUM=SUM+P(I)
165                ERP=SD/SUM
166                RETURN
167                END

li913 · 发表于 2019-4-19 08:35:08

给出错误提示。

L25 · 发表于 2022-11-8 19:26:39

你好，请问解决了吗

[其他行业算法] 弹流润滑数值计算程序迭代失败 ！！求助！！

[其他行业算法] 弹流润滑数值计算程序迭代失败！！求助！！