洗脚的好处《弹性流体动压润滑》黄平、fortran程序注释、等温线接触

2021年1月19日发(作者：孟达)
PROGRAM GREASELINE

CHARACTER*1 S,S1,S2

CHARACTER*16 CDATE,CTIME

COMMON
/COM1/ENDA,A1,A2,A3,Z,C1,C3,CW,LMAX ,FF/COM2/EDA0/COM4/X0,XE/COM3/E1,PH,
B,U1,U2,R


DATA PAI,Z,P0/3.14159265,0.68,1.96E8 /S1,S2/1HY,1Hy/!
预定义数值


DATA
N, X0,XE,W,E1,EDA0,R,Us,CU,C1,FN/129,-4.,1.4,1.E5,2.2 6E11,0.41,0.0128,0.87,0.67,0.5,0
.846/CDATE,CTI ME/'The date is','The time is'/!
预赋值参数





OPEN(8,FILE='',STATUS='UNKNOWN' )!
建立输出结果文件

1


FORMAT(20X,A12,I2.2,':',I2.2,':',I4.4)
2


FORMAT(20X,A12,I2.2,':',I2.2,':',I2. 2,'.',I2.2)

WRITE(*,*)'
参数已经赋值？

(Y or N)?'!
询问是否在程序中更改参数


READ(*,'(A)')S

IF(.S2)THEN


GOTO 10

ENDIF

WRITE(*,*)'PH=' !
程序运行中更改数值，
PH
为
Hertz
接触压力


READ(*,*)PH

W=2.*PAI*R*PH*(PH/E1)

WRITE(*,*)'W=',W
10

CW=N+0.1




FF=1./FN

LMAX=ALOG(CW)/ALOG(2.)!
将参数量纲一化


N=2**LMAX+1

LMIN=(ALOG(CW)-ALOG(SQRT(CW)))/ALOG(2.)

LMAX=LMIN

W1=W/(E1*R)

PH=E1*SQRT(0.5*W1/PAI)

A1=(ALOG(EDA0)+9.67)

A2=PH/P0

A3=0.59/(PH*1.E-9)

B=4.*R*PH/E1

ALFA=Z*A1/P0

G=ALFA*E1

U=EDA0*US/(2.*E1*R)

C3=1.6*(R/B)**2*G**0.6*U**0.7*W1**(-0.13)

ENDA=B**(2.+FF)*(PH/2/EDA0)**FF/R**(1+FF)/US/( 2.+FF)

U1=0.5*(2.+CU)*U

U2=0.5*(2.-CU)*U




WRITE(*,*)'B,PH,G,U=',B,PH,G,U

CW=-1.13*C3

WRITE(*,*)N,X0,XE,W,E1,EDA0,R,US,PH

WRITE(8,*)N,W,E1,EDA0,R,US,B,PH,FF

WRITE(*,40)
40

FORMAT(2X,'Wait

Please',//)

CALL SUBAK(N)!
计算弹性变形系数


CALL MULTI(N)!
计算压力
P
和
H

STOP

END

SUBROUTINE MULTI(N)!
计算压力
P
和
H

REAL*8 < br>X(1100),P(1100),H(1100),RO(1100),POLD(1100),EPS (1100),EDA(1100),R(1100),K(1100),E
(1100)

COMMON
/COM1/ENDA,A1,A2,A3,Z,C1,C3,CW,LMAX ,FF/COM4/X0,XE/COM3/E1,PH,B,U1,U2,RR

DATA MK,G0/1,1.570796325/

NX=N

DX=(XE-X0)/(N-1.0)

DO 10 I=1,N

X(I)=X0+(I-1)*DX

IF(ABS(X(I)).GE.1.0)P(I)=0.0

IF(ABS(X(I)).LT.1.0)P(I)=SQRT(1.-X(I)*X(I))
10

CONTINUE

CALL
HREE(N, DX,H00,G0,X,P,H,RO,EPS,EDA)!
实现各点的膜厚、密度、粘度、弹性变形的计算


CALL FZ(N,P,POLD)
14
KK=19

CALL
ITER(N,KK,DX,H00,G0,X, P,H,RO,EPS,EDA,R)!
利用雷诺方程进行各点压力的重
新计算，并再次调用< br>HREE
子程序


MK=MK+1

CALL ERROP(N,P,POLD,ERP)!
计算迭代前后的压力差


IF(.800)THEN

GOTO 14

ENDIF

WRITE(*,*)PH,RR,B
105 IF(.800)THEN

WRITE(*,*)'Pressures are not convergent !!!'

READ(*,*)

ENDIF

FM=FRICT(N,DX,X,H,P,EDA)

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)
106 CONTINUE




DO 108 I=1,N




K(I)=P(I)*PH/1.E9

E(I)=H(I)*B*B*1.E6/RR
108 CONTINUE

H3=H2*B*B/RR

P3=P2*PH
110 FORMAT(6(1X,E12.6))
120 CONTINUE

WRITE(8,*)'P2,H2,P3,H3=',P2,H2,P3,H3

CALL OUTHP(N,X, K,E)!
实现结果的输出功能


RETURN

END

SUBROUTINE OUTHP(N,X, K,E)!
实现结果的输出功能


REAL*8

X(N), K(N),E(N)

DX=X(2)-X(1)

DO 10 I=1,N

WRITE(8,20)X(I),K(I),E(I)
10

CONTINUE
20

FORMAT(1X,6(F20.6,1X))

RETURN

END

SUBROUTINE HREE(N,DX,H00,G0,X,P,H ,RO,EPS,EDA)!
实现各点的膜厚、
密度、
粘度、
弹性变形的计算





REAL*8

X(N),P(N),H(N),RO(N),EPS(N),EDA(2200)

REAL*8

W(2200)

COMMON
/C OM1/ENDA,A1,A2,A3,Z,C1,C3,CW,K,FF/COM2/EDA0/COMAK/ AK(0:1100)




DATA KK,NW,PAI1/0,2200,0.318309886/

IF(.0)GOTO 3

HM0=C3

H00=0.0
3
W1=0.0

DO 4 I=1,N
4


W1=W1+P(I)!
此压力下的量纲一化载荷
W

C3=(DX*W1)/G0

DW=1.-C3 !
承载力判断值


CALL DISP(N,NW,K,DX,P,W)!
计算各点弹性变形量


HMIN=1.E3

DO 30 I=1,N

H0=0.5*X(I)*X(I)-PAI1*W(I)!

IF()HMIN=H0

H(I)=H0
30
CONTINUE




IF(.0)GOTO 32

KK=1

H00=-HMIN+HM0
32
H0=H00+HMIN

IF(.0.0)GOTO 48

IF(H0+0.3*CW*.0.0)HM0=H0+0.3*CW*DW

IF(H0+0.3*CW*.0.0)HM0=HM0*C3
48
H00=HM0-HMIN
50
DO 60 I=1,N
60
H(I)=H00+H(I)!
膜厚计算


DO 100 I=1,N

EDA(I)=EXP(A1*(-1.+(1.+A2*P(I))**Z) )!
各点粘度数值计算


RO(I)=1.0

EPS(I)=RO(I)*H(I)**(2+FF)*ENDA/EDA(I)**FF !
离散化
Reynolds
方程中的系数

100 CONTINUE

RETURN

END

SUBROUTINE ITER(N,KK,DX,H00,G0,X,P,H,RO,EPS,ED A,R)!
利用雷诺方程进行各点压
力的重新计算，并再次调用
HREE
子程 序


REAL*8 X(N),P(N),H(N),RO(N),EPS(N),EDA(N),R(N)

COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C3,CW,LMAX,FF/ COMAK/AK(0:1100)

DATA KG1,PAI/0,3.14159265/

IF(.0)GOTO 5

KG1=1

DX1=1./DX

DX2=DX*DX

DX3=1./DX2

DX4=DX1/PAI

DX5=DX1**(1+FF)

DXL=DX*ALOG(DX)

AK0=DX*AK(0)+DXL

AK1=DX*AK(1)+DXL
5
DO 100 K=1,KK

D2=0.5*(EPS(1)+EPS(2))

D3=0.5*(EPS(2)+EPS(3))

D5=DX1*(RO(2)*H(2)-RO(1)*H(1))

D7=DX4*(RO(2)*AK0-RO(1)*AK1)

PP=0.

DO 70 I=2,N-1

D1=D2

D2=D3