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

弹流润滑数值计算方法

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