弹流润滑数值计算方法-赵江

弹流润滑数值计算方法
1.1等温线接触弹流润滑数值计算方法与程序
等温点接触弹流润滑的基本方程包括：Reynolds 方程、膜厚方程、变形方程、粘压方程、密压方程和载荷平衡方程，主要形式如下
Reynolds 方程：
x
ρh 12μ=x p ηρh x 3d )d()d d (d d s （1） 膜厚方程
)(+2R
+=)(2
0x υx h x h （2）
变形方程
∫e
x +d ln(2
-=x c s x -s x p E
πx υ))()( （3）
粘压方程
{}]
)()[.(z p p ηηη000+1+1-679+ln exp = （4）
密压方程
)..p
p
ρρ71+160+
1=0（ （5）
载荷平衡方程
∫
0=d -0e
x x x x p ω)( （6）
计算程序
计算工况参数 节点数：N=30
量纲一花起始点坐标：X0=-4.0 量纲一花终止点坐标：XE=1.5 载荷：W=1.0E5
综合弹性模量E1=2.2E11 初始粘度：EDA0=0.05 圆柱半径：R=0.05 速度:US=1.5 主程序：
FORGRAM LINEEHL
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COM4/X0,XE/COM3/E1,PH,B,R
DATA PAI,Z,P0/3.14159565,0.68,1.96E8/
DATA N,X0,XE,W,E1,EDA0,R,US/130,-4.0,1.5,1.0E5,2.2E11,0.05,0.05,1.5/
OPEN(8,FILE=’OUT.DAT’,STA TUS=’UNKNOWN’)
W1=W/(E1*R)
C W无量纲公式，
PH=E1*SQRT(0.5*W1/PAI)
C PH为最大hertz接触应力Ph，单位pa
A1=(ALOG(EDA0)+9.67)
C 粘度公式计算时的参数
A2=PH/P0
C 粘度公式计算时的参数
A3=0.59/(PH*1.E-9)
C 密度公式计算时参数
B=4.*R*PH/E1
C B为HERTZ接触区半宽，单位m
ALFA=Z*A1/P0 C粘度系数
G=ALFA*E1 C材料参数
U=EDA0*US/(2.*E1*R)
C 苏联人提出的入口区分析解中定义的
CC1=SQRT(2.*U)
AM=2.*PAI*(PH/E1)**2/CC1
ENDA=3.*(PAI/AM)**2/8
C ENDA量钢化reynolds中的栏目达
HM0=1.6*(R/B)**2*G**0.6*U**0.7*W1**(-0.13)
C hm0最小膜厚，公式为最小膜厚公式
WRITE(*,*)N,X0,XE,W,E1,EDA0,R,US
CALL SUBAK(N)
CALL EHL(N)
STOP
END
SUBROUTINE EHL(H)
C 弹流润滑子程序
DIMENSION X(1100),P(1100),H(1100),RO(1100),EPS(1100),EDA(1100),V(1100)
C X P压力、H膜厚、RO密度、EPS雷诺方程中的E、EDA粘度、V弹性变形
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XE
C ENDA量钢化reynolds中的栏目达, hm0最小膜厚
COMMON/COM3/E1,PH,B,RR
C E1当量弹性模量E，PH为最大hertz接触应力ph，B为HERTZ接触区半宽，RR
MK=1
DX=(XE-X0)/(N-1.0)
C DX等距节点间距
DO 10 I=1,N
X(I)=X0+(I-1)*DX
IF(ABS(X(I)).GE,1.0)P(I)=0.0
C 压力分布为抛物线形式,
IF(ABS(X(I)).LT,1.0)
C 压力分布为抛物线形式
10 CONTINUE
CALL HREE(N,KK,DX,X,P,H,RO,EPS,EDA,V)
CALL FZ(N,P,POLD)
C 对应于框图中,调用HREE子程序计算个节点膜厚,密度ro,和粘度EDA，
C 并调用子程序VI计算个节点的弹性变形
C 调用子程序FZ将P(I)的值赋给POLD,POLD第i个节点的前一次迭代压力值
14 KK=19
CALL ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
C 调用ITER子程序进行雷诺方程迭代，重新计算hree压力分布下的节点膜厚/粘度/密度
MK=MK+1
CALL ERROP(N,P,POLD,ERP)
C 调用ERROP子程序计算迭代前后的压力差并输出
WRITE(*,*)’ERP=’,ERP
IF(ERP.GT.1.E-5.AND.DH.GT.1.E-6)THEN
C .GT.大于
IF(MK.GE.50)THEN
MK=1
DH=0.5*DH
ENDIF
GOTO14
ENDIF
C 判断ERP和DH是否满足条件
IF(DH.LE.1.E-6)
WRITE(*,*)’pressure are not convergent’
H2=1.E3
P2=0.0
DO 106 I=1,N
IF(H(I).LT.H2)H2=H(I)
IF(P(I).GT.P2)P2=P(I)
C 找出最小膜厚H2和最大压力P2,量纲一化并输出
106 CONTINUE
H3=H2*B*B/RR
P3=P2*PH
C 这个地方还不是很懂,可能是获得最小的值然后量纲计算出有单位的值
110 FORMAT(6(1X,E12.6))
C 格式化输出1X把输出位置向右跳过n个位置
C E12.6以12个字符宽输出指数类型的浮点数，小数部分占6个字符宽
120 CONTINUE
WRITE(*,*)‘P2,H2,P3,H3=’,P2,H2,P3,H3
CALL OUTHP(N,X,P,H)
RETURN
SUBROUTINE OUTHP(N,X,P,H)
C 输出子程序
DIMENSION X(N),P(N),H(N)
DO 10 I=1,N
WRITE(8,20)X(I),P(I),H(I)
10 CONTINUE
20 FORMA T(1X,6(E12,6,1X))
RETURN
END
SUBROUTINE HREE(N,DX,X,P,H,RO,EPS,EDA,V)
DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COMAK/AK(0:1100)
DATA KK,PAI1,G0/0,0.318309886,1.570796325/
IF(KK.NE.0)GOTO3
H00=0.0
3 W1=0.0
DO 4 I=1,N
4 W1=W1+P(I)
C3=(DX*W1)/G0
DW=1.-C3
C 不是很懂
CALL VI(N,DX,P,V)
HMIN=1.E3
DO 30 I=1,N
H0=0.5*X(I)*X(I)+V(I)
C 和膜厚公式有关
IF(H0.LT.HMIN) HMIN=H0
H(I)=H0
30 CONTINUE
IF(KK.NE.0)GOTO 32
KK=1
DH=0.005*HM0
C 逐步减小法
H00=-HMIN+HM0
32 IF(DW.LT.0.0)H00=H00+DH
IF(DW.GT.0.0)H00=H00-DH
DO 60 I=1,N
H(I)=H00+H(I)
EDA(I)=EXP(A1*(-1.+(1.+A2*P(I)**Z))
C 就是粘压公式
RO(I)=(A3+1.35*P(I))/(A3+P(I))
EPS(I)=RO(I)*H(I)**3/(ENDA*EDA(I))
C 求出了雷诺方程中的一副三拉，为后面雷诺方程的迭代计算做基础
60 CONTINUE
RETURN
SUBROUTINE ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
C 雷诺方程中的迭代
DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
COMMON/COMAK/AK(0:1100)
DATA PAI1/0.318309886/
DO 100 K=1,KK
D2=0.5*(EPS(1)+EPS(2))
D3=0.5*(EPS(2)+EPS(3))
DO 70 I=2, N-1
D1=D2
D2=D3
IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2))
D8=RO(I)*AK(0)*PAI1
D9=RO(I-1)*AK(1)*PAI1
D10=1.0/(D1+D2+(D9-D8)*DX)
D11=D1*P(I-1)+D2*P(1+1)
D12=(RO(I)*H(I)-RO(I-1)*H(I-1)+(D8-D9)*P(I))*DX
P(I)=(D11-D12)*D10
IF(P(I).LT.0.0)P(I)=0.0
C 计算雷诺差分方程个系数，计算节点压力，并保证节点压力不小于0(这里还要仔细考虑)
70 CONTINUE
CALL HREE(N,DX,X,P,H,RO,EPS,EDA,V)
C 调用HREE方程重新计算各节点的膜厚，
100 CONTINUE
RETURN
END
SUBROUTINE VI(N,DX,P,V)
DIMENSION P(N),V(N)
COMMON/COMAK/AK(0:1100)
PAI1=0.318309886
C=ALOG(DX)
DO 10 I=1,N
V(I)=0.0
DO 10 J=1,N
IJ=IABS(I-J)
10 V(I)=V(I)+(AK(IJ)+C)*DX*P(J)
DO I=1,N
V(I)=-PAI1*V(I)
ENDDO
RETURN
END
SUBROUTINE SUBAK(MM)
COMMON/COMAK/AK(0:1100)
DO 10 I=1,MM
10 AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.)
END
SUBROUTINE FZ(N,P,POLD)
DIMENSION P(N),POLD(N)
DO 10 I=1,N
10 POLD(I)=P(I)
RETURN
END
SUBROUTINE ERROP(N,P,POLD,ERP)
DIMENSION P(N),POLD(N)
SD=0.0
SUM=0.0
DO 10 I=1,N
SD=SD+ABS(P(I)-POLD(I))
10 SUM=SUM+P(I)
ERP=SD/SUM
RETUREN
END
FORGRAM LINEEHL
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0
COMMON/COM4/X0,XE/COM3/E1,PH,B,R
DATA PAI,Z,P0/3.14159565,0.68,1.96E8/
DATA N,X0,XE,W,E1,EDA0,R/130,-4.0,1.5,1.0E5,2.2E11,0.05,0.05/ DATA US/1.5/
OPEN(8,FILE=’OUT.DAT’,STATUS=’UNKNOWN’)
W1=W/(E1*R)
C W无量纲公式，
PH=E1*SQRT(0.5*W1/PAI)
C PH为最大hertz接触应力Ph，单位pa
A1=(ALOG(EDA0)+9.67)
C 粘度公式计算时的参数
A2=PH/P0
C 粘度公式计算时的参数
A3=0.59/(PH*1.E-9)
C 密度公式计算时参数
B=4.*R*PH/E1
C B为HERTZ接触区半宽，单位m
ALFA=Z*A1/P0 C粘度系数
G=ALFA*E1 C材料参数
U=EDA0*US/(2.*E1*R)
C 苏联人提出的入口区分析解中定义的
CC1=SQRT(2.*U)
ENDA=3.*(PAI/AM)**2/8
C ENDA量钢化reynolds中的栏目达
HM0=1.6*(R/B)**2*G**0.6*U**0.7*W1**(-0.13)
C hm0最小膜厚，公式为最小膜厚公式
WRITE(*,*)N,X0,XE,W,E1,EDA0,R,US
CALL SUBAK(N)
CALL EHL(N)
STOP
END
SUBROUTINE EHL(H)
C 弹流润滑子程序
DIMENSION X(1100),P(1100),H(1100),RO(1100),EPS(1100),EDA(1100)
DIMENSION V(1100)
C X P压力、H膜厚、RO密度、EPS雷诺方程中的E、EDA粘度、V弹性变形
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XE
C ENDA量钢化reynolds中的栏目达, hm0最小膜厚
COMMON/COM3/E1,PH,B,RR
C E1当量弹性模量E，PH为最大hertz接触应力ph， B为HERTZ接触区半宽，RR
MK=1
DX=(XE-X0)/(N-1.0)
C DX等距节点间距
DO 10 I=1,N
X(I)=X0+(I-1)*DX
C 定位X(I)的位置,对应于框图中的计算个节点X(I),和hertz压力值
IF(ABS(X(I)).GE,1.0)P(I)=0.0
C 压力分布为抛物线形式,
IF(ABS(X(I)).LT,1.0)
C 压力分布为抛物线形式
10CONTINUE
CALL HREE(N,KK,DX,X,P,H,RO,EPS,EDA,V)
CALL FZ(N,P,POLD)
C 对应于框图中,调用HREE子程序计算个节点膜厚,密度ro,和粘度EDA，
C 并调用子程序VI计算个节点的弹性变形
C 调用子程序FZ将P(I)的值赋给POLD,POLD第i个节点的前一次迭代压力值
14 KK=19
CALL ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
C 调用ITER子程序进行雷诺方程迭代，重新计算hree压力分布下的节点膜厚/粘度/密度 MK=MK+1
CALL ERROP(N,P,POLD,ERP)
C 调用ERROP子程序计算迭代前后的压力差并输出
WRITE(*,*)’ERP=’,ERP
IF(ERP.GT.1.E-5.AND.DH.GT.1.E-6)THEN
C .GT.大于
IF(MK.GE.50)THEN
MK=1
ENDIF
GOTO14
ENDIF
C 判断ERP和DH是否满足条件
IF(DH.LE.1.E-6)
WRITE(*,*)’pressure are not convergent’
H2=1.E3
P2=0.0
DO 106 I=1,N
IF(H(I).LT.H2)H2=H(I)
IF(P(I).GT.P2)P2=P(I)
C 找出最小膜厚H2和最大压力P2,量纲一化并输出
106CONTINUE
H3=H2*B*B/RR
P3=P2*PH
C 这个地方还不是很懂,可能是获得最小的值然后量纲计算出有单位的值110FORMAT(6(1X,E12.6))
C 格式化输出1X把输出位置向右跳过n个位置
C E12.6以12个字符宽输出指数类型的浮点数，小数部分占6个字符宽120CONTINUE
WRITE(*,*)‘P2,H2,P3,H3=’,P2,H2,P3,H3
CALL OUTHP(N,X,P,H)
RETURN
END
SUBROUTINE OUTHP(N,X,P,H)
C 输出子程序
DIMENSION X(N),P(N),H(N)
DO 10 I=1,N
WRITE(8,20)X(I),P(I),H(I)
10CONTINUE
20FORMAT(1X,6(E12,6,1X))
RETURN
END
SUBROUTINE HREE(N,DX,X,P,H,RO,EPS,EDA,V)
DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COMAK/AK(0:1100)
DATA KK,PAI1,G0/0,0.318309886,1.570796325/
IF(KK.NE.0)GOTO3
H00=0.0
3 W1=0.0
DO 4 I=1,N
4 W1=W1+P(I)
C3=(DX*W1)/G0
DW=1.-C3
C 不是很懂
HMIN=1.E3
DO 30 I=1,N
H0=0.5*X(I)*X(I)+V(I)
C 和膜厚公式有关
IF(H0.LT.HMIN) HMIN=H0
H(I)=H0
30CONTINUE
IF(KK.NE.0)GOTO 32
KK=1
DH=0.005*HM0
C 逐步减小法
H00=-HMIN+HM0
32IF(DW.LT.0.0)H00=H00+DH
IF(DW.GT.0.0)H00=H00-DH
DO 60 I=1,N
H(I)=H00+H(I)
EDA(I)=EXP(A1*(-1.+(1.+A2*P(I)**Z))
C 就是粘压公式
RO(I)=(A3+1.35*P(I))/(A3+P(I))
EPS(I)=RO(I)*H(I)**3/(ENDA*EDA(I))
C 求出了雷诺方程中的一副三拉，为后面雷诺方程的迭代计算做基础
60CONTINUE
RETURN
END
SUBROUTINE ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
C 雷诺方程中的迭代
DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
COMMON/COMAK/AK(0:1100)
DATA PAI1/0.318309886/
DO 100 K=1,KK
D2=0.5*(EPS(1)+EPS(2))
D3=0.5*(EPS(2)+EPS(3))
DO 70 I=2, N-1
D1=D2
D2=D3
IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2))
D8=RO(I)*AK(0)*PAI1
D9=RO(I-1)*AK(1)*PAI1
D10=1.0/(D1+D2+(D9-D8)*DX)
D11=D1*P(I-1)+D2*P(1+1)
D12=(RO(I)*H(I)-RO(I-1)*H(I-1)+(D8-D9)*P(I))*DX
P(I)=(D11-D12)*D10
IF(P(I).LT.0.0)P(I)=0.0
C 计算雷诺差分方程个系数，计算节点压力，并保证节点压力不小于0(这里还要仔细考虑) 70CONTINUE
C 调用HREE方程重新计算各节点的膜厚，
100CONTINUE
RETURN
END
SUBROUTINE VI(N,DX,P,V)
DIMENSION P(N),V(N)
COMMON/COMAK/AK(0:1100)
PAI1=0.318309886
C=ALOG(DX)
DO 10 I=1,N
V(I)=0.0
DO 10 J=1,N
IJ=IABS(I-J)
10 V(I)=V(I)+(AK(IJ)+C)*DX*P(J)
DO I=1,N
V(I)=-PAI1*V(I)
ENDDO
RETURN
END
SUBROUTINE SUBAK(MM)
COMMON/COMAK/AK(0:1100)
DO 10 I=1,MM
10 AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.)
RETURN
END
SUBROUTINE FZ(N,P,POLD)
DIMENSION P(N),POLD(N)
DO 10 I=1,N
10 POLD(I)=P(I)
RETURN
END
SUBROUTINE ERROP(N,P,POLD,ERP)
DIMENSION P(N),POLD(N)
SD=0.0
SUM=0.0
DO 10 I=1,N
SD=SD+ABS(P(I)-POLD(I))
10 SUM=SUM+P(I)
ERP=SD/SUM
RETURN
END
PROGRAM MAIN
WRITE(*,*) 'HELLO'
WRITE(*,*)
1'HELLO'
100WRITE(*,*)'HELLO'
10STOP
END
1.2等温点接触弹流润滑数值计算方法与程序
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.6*(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)
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(*,*)’PRESSURE 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=0,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,N
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)
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,N
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=HOO-DH
DO 60 I=1,N
DO 60 J=1,N
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)GOTO2
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(I1,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 CONYINUE
RETURN
END
SUBROUTINE VI(N,DX,P,V)
DIMENSION P(N,N),V(N,N)
COMMON/COMAK/AK(0:65,0:65)
PAI1=0.0202643
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
SUBROUTINE SUBAK(MM)
COMMON/COMAK/AK(0:65,0:65)
S(X,Y)=X+SQRT(X**2+Y**2)
DO 10 I=1,N
XP=I+0.5
XM=I-0.5
DO 10 J=1,N
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)
SIMENSION 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.5,1X))
RETURN
END。