弹流润滑数值计算方法-赵江_第1页
弹流润滑数值计算方法-赵江_第2页
弹流润滑数值计算方法-赵江_第3页
弹流润滑数值计算方法-赵江_第4页
弹流润滑数值计算方法-赵江_第5页
已阅读5页,还剩8页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、弹流润滑数值计算方法1.1等温线接触弹流润滑数值计算方法与程序等温点接触弹流润滑的根本方程包括:平衡方程,主要形式如下Reynolds方程:膜厚方程Reynolds方程、膜厚方程、变形方程、粘压方程、密压方程和载荷dpdp一()=12的dx'ddxd(p1)dx(1)(2)(3)2xh(x)=ho+-+u(x)2R变形方程2x(x)=-e/p(x)ln(s-x)ds+c粘压方程=助explnq+9.67)-1+(1+p/Po)Z(4)密压方程(5)载荷平衡方程(6)计算程序计算工况参数节点数:N=30量纲一花起始点坐标:X0=-4.0量纲一花终止点坐标:XE=1.5载荷:W=1.0E5

2、综合弹性模量E1=2.2E11初始粘度:EDA0=0.05圆柱半径:R=0.05速度:US=1.5主程序:FORGRAMLINEEHLCOMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COM4/X0,XE/COM3/E1,PH,B,RDATAPAI,Z,P0/3.14159565,0.68,1.96E8/DATAN,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,STAS=UNKNOWN)W1=W/(E1*R)CW无量纲公式,PH=

3、E1*SQRT(0.5*W1/PAI)CPH为最大hertz接触应力Ph,单位paA1=(ALOG(EDA0)+9.67)C粘度公式计算时的参数A2=PH/P0C粘度公式计算时的参数A3=0.59/(PH*1.E-9)C密度公式计算时参数B=4.*R*PH/E1CB为HERTZ接触区半宽,单位mALFA=Z*A1/P0C粘度系数G=ALFA*E1C材料参数U=EDA0*US/(2.*E1*R)C苏联人提出的入口区分析解中定义的CC1=SQRT(2.*U)AM=2.*PAI*(PH/E1)*2/CC1ENDA=3.*(PAI/AM)*2/8CENDA量钢化reynolds中的栏目达HM0=1.6

4、*(R/B)*2*G*0.6*U*0.7*W1*(-0.13)Chm0最小膜厚,公式为最小膜厚公式WRITE(*,*)N,X0,XE,W,E1,EDA0,R,USCALLSUBAK(N)CALLEHL(N)STOPENDSUBROUTINEEHL(H)C弹流润滑子程序DIMENSIONX(1100),P(1100),H(1100),RO(1100),EPS(1100),EDA(1100),V(1100)CXP压力、H膜厚、RO密度、EPS雷诺方程中的E、EDA粘度、V弹性变形COMMON/COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XECENDA量钢化reynold

5、s中的栏目达,hm0最小膜厚COMMON/COM3/E1,PH,B,RRCE1当量弹性模量E,PH为最大hertz接触应力ph,B为HERTZ接触区半宽,RRMK=1DX=(XE-X0)/(N-1.0)CDX等距节点间距DO10I=1,NX(I)=X0+(I-1)*DXC定位X(I)的位置对应于框图中的计算个节点X(I),和hertz压力值IF(ABS(X(I).GE,1.0)P(I)=0.0C压力分布为抛物线形式,IF(ABS(X(I).LT,1.0)C压力分布为抛物线形式10CONTINUECALLHREE(N,KK,DX,X,P,H,RO,EPS,EDA,V)CALLFZ(N,P,POL

6、D)C对应于木g图中,调用HREE子程序计算个节点膜厚,密度ro,和粘度EDA,C并调用子程序VI计算个节点的弹性变形C调用子程序FZ将P(I)的值赋给POLD,POLD第i个节点的前一次迭代压力值14KK=19CALLITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)/粘度/密度C调用ITER子程序进行雷诺方程迭代,重新计算hree压力分布下的节点膜厚MK=MK+1CALLERROP(N,P,POLD,ERP)C调用ERROP子程序计算迭代前后的压力差并输出WRITE(*,*)'ERP=',ERPIF(ERP.GT.1.E-5.AND.DH.GT.1.E-6)T

7、HENC.GT.大于IF(MK.GE.50)THENMK=1DH=0.5*DHENDIFGOTO14ENDIFC判断ERP和DH是否满足条件IF(DH.LE.1.E-6)WRITE(*,*)'pressurearenotconvergent!'H2=1.E3P2=0.0DO106I=1,NIF(H(I).LT.H2)H2=H(I)IF(P(I).GT.P2)P2=P(I)C找出最小膜厚H2和最大压力P2,量纲一化并输出106CONTINUEH3=H2*B*B/RRP3=P2*PHC这个地方还不是很懂,可能是获得最小的值然后量纲计算出有单位的值110FORMAT(6(1X,E12

8、.6)C格式化输出1X把输出位置向右跳过n个位置CE12.6以12个字符宽输出指数类型的浮点数,小数局部占6个字符宽120CONTINUEWRITE(*,*)'P2,H2,P3,H3=',P2,H2,P3,H3CALLOUTHP(N,X,P,H)RETURNENDC102034CC30C32CC60SUBROUTINEOUTHP(N,X,P,H)输出子程序DIMENSIONX(N),P(N),H(N)DO10I=1,NWRITE(8,20)X(I),P(I),H(I)CONTINUEFORMAT(1X,6(E12,6,1X)RETURNENDSUBROUTINEHREE(N,D

9、X,X,P,H,RO,EPS,EDA,V)DIMENSIONX(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)DATAKK,PAI1,G0/0,0.318309886,1.570796325/IF(KK.NE.0)GOTO3H00=0.0W1=0.0DO4I=1,NW1=W1+P(I)C3=(DX*W1)/G0DW=1.-C3不是很懂CALLVI(N,DX,P,V)HMIN=1.E3DO30I=1,NH0=0.5*X(I)*X(I)+V(I)和膜

10、厚公式有关IF(H0.LT.HMIN)HMIN=H0H(I)=H0CONTINUEIF(KK.NE.0)GOTO32KK=1DH=0.005*HM0逐步减小法H00=-HMIN+HM0IF(DW.LT.0.0)H00=H00+DHIF(DW.GT.0.0)H00=H00-DHDO60I=1,NH(I)=H00+H(I)EDA(I)=EXP(A1*(-1.+(1.+A2*P(I)*Z)就是粘压公式RO(I)=(A3+1.35*P(I)/(A3+P(I)EPS(I)=RO(I)*H(I)*3/(ENDA*EDA(I)求出了雷诺方程中的一副三拉,为后面雷诺方程的迭代计算做根底CONTINUERETU

11、RNENDSUBROUTINEITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)C雷诺方程中的迭代DIMENSIONX(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)COMMON/COMAK/AK(0:1100)DATAPAI1/0.318309886/DO100K=1,KKD2=0.5*(EPS(1)+EPS(2)D3=0.5*(EPS(2)+EPS(3)DO70I=2,N-1D1=D2D2=D3IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2)D8=RO(I)*AK(0)*PAI1D9=RO(I-1)*AK(1)*PAI1D

12、10=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)*DXP(I)=(D11-D12)*D10IF(P(I).LT.0.0)P(I)=0.0C计算雷诺差分方程个系数,计算节点压力,并保证节点压力不小于0(这里还要仔细考虑)70CONTINUECALLHREE(N,DX,X,P,H,RO,EPS,EDA,V)C调用HREE方程重新计算各节点的膜厚,100CONTINUERETURNENDSUBROUTINEVI(N,DX,P,V)DIMENSIONP(N),V(N)C

13、OMMON/COMAK/AK(0:1100)PAI1=0.318309886C=ALOG(DX)DO10I=1,NV(I)=0.0DO10J=1,NIJ=IABS(I-J)10V(I)=V(I)+(AK(IJ)+C)*DX*P(J)DOI=1,NV(I)=-PAI1*V(I)ENDDORETURNENDSUBROUTINESUBAK(MM)COMMON/COMAK/AK(0:1100)DO10I=1,MM10AK(I)=(I+0.5)*(ALOG(ABS(I+0.5)-1.)-(I-0.5)*(ALOG(ABS(I-0.5)-1.)RETURNENDSUBROUTINEFZ(N,P,POLD)

14、DIMENSIONP(N),POLD(N)DO10I=1,N10POLD(I)=P(I)RETURNENDSUBROUTINEERROP(N,P,POLD,ERP)DIMENSIONP(N),POLD(N)SD=0.0SUM=0.0DO10I=1,NSD=SD+ABS(P(I)-POLD(I)10SUM=SUM+P(I)ERP=SD/SUMRETURENENDFORGRAMLINEEHLCOMMOOM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0COMMOOM4/X0,XE心OM3/E1,PH,B,RDATAPAI,Z,P0/3.14159565,0.68,1.96E8/

15、DATAN,X0,XE,W,E1,EDA0,R/130,-4.0,1.5,1.0E5,2.2E11,0.05,0.05/DATAUS/1.5/OPEgFILE='OUTDAT',STATUS=UNKNOWNW1=W/(E1*R)CW无量纲公式,PH=E1*SQRT0.5*W1/PAI)CPH为最大hertz接触应力Ph,单位paA1=(ALOGEDA0)+9.67)C粘度公式计算时的参数A2=PH/P0C粘度公式计算时的参数A3=0.59/(PH*1.E-9)C密度公式计算时参数B=4.*R*PH/E1CB为HERT菽触区半宽,单位mALFA=Z*A1/P0C粘度系数G=ALF

16、A*E1C材料参数U=EDA0*US/(2.*E1*R)C苏联人提出的入口区分析解中定义的CC1=SQRQU)AM=2.*PAI*(PH/E1)*2/CC1ENDA=3.*(PAI/AM)*2/8CENDA量钢化reynolds中的栏目达HM0=1.6*(R/B)*2*G*0.6*U*0.7*W1*(-0.13)Chm0最小膜厚,公式为最小膜厚公式WRITR*)N,X0,XE,W,E1,EDA0,R,USCALLSUBAK(N)CALLEHL(N)STOPENDSUBROUTINEHL(H)C弹流润滑子程序DIMENSIOX(1100),P(1100),H(1100),RO(1100),EPS

17、(1100),EDA(1100)DIMENSIOV(1100)CXP压力、H1厚、R蹴度、EPSf诺方程中的E、EDA占度、V单性变形COMMONDM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XECENDA量钢化reynolds中的栏目达,hm0最小膜厚COMMONDM3/E1,PH,B,RRCE1当量弹性模量E,PH最大hertz接触应力ph,B为HERT接触区半宽,RRMK=1DX=(XE-X0)/(N-1.0)CDX等距节点间距DO10I=1,NX(I)=X0+(I-1)*DXC定位X(I)的位置,对应于框图中的计算个节点X(I),和hertz压力值IF(ABSX

18、(I).GE,1.0)P(I)=0.0C压力分布为抛物线形式,IF(ABSX(I).LT,1.0)C压力分布为抛物线形式10CONTINUECALLHREE(N,KK,DX,X,P,H,RO,EPS,EDA,V)CALLFZ(N,P,POLD)C对应于力I图中,调用HRE子程序计算个节点膜厚,密度ro,和粘度EDAC并调用子程序VI计算个节点的弹性变形C调用子程序FZ等P(I)的值赋给POLD,POLDi个节点的前一次迭代压力值14KK=19CALLITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)C调用ITER子程序进行雷诺方程迭代,重新计算hree压力分布下的节点膜厚/粘度/

19、密度MK=MK+1CALLERROP(N,P,POLD,ERP)C调用ERROP程序计算迭代前后的压力差并输出WRITE*,*)'ERP=,ERPIF(ERP.GT.1.E-5.AND.DH.GT.1.E-6)THENC.GT.大于IF(MK.GE.50)THENMK=1DH=0.5*DHENDIFGOTO14ENDIFC判断ER和D是否满足条件IF(DH.LE,1,E-6)WRITE*,*)'pressurearenotconvergent!'H2=1.E3P2=0.0DO106I=1,NIF(H(I).LT.H2)H2=H(I)IF(P(I).GT.P2)P2=P(

20、I)C找出最小膜厚H环口最大压力P2,量纲一化并输出106CONTINUEH3=H2*B*B/RRP3=P2*PHC这个地方还不是很懂,可能是获得最小的值然后量纲计算出有单位的值110FORMAT(1X,E12.6)C格式化输出1片巴输出位置向右跳过n个位置CE12.6以12个字符宽输出指数类型的浮点数,小数局部占6个字符宽120CONTINUEWRITE*,*)'P2,H2,P3,H3=',P2,H2,P3,H3CALLOUTHP(N,X,P,H)RETURNENDSUBROUTINOUTHP(N,X,P,H)C输出子程序DIMENSIOX(N),P(N),H(N)DO10I

21、=1,NWRITE8,20)X(I),P(I),H(I)10CONTINUE20FORMATX,6(E12,6,1X)RETURNENDSUBROUTINEREE(N,DX,X,P,H,RO,EPS,EDA,V)DIMENSIOX(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)COMMONOM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COMAK/AK(0:1100)DATAKK,PAI1,G0/0,0.318309886,1.570796325/IF(KK.NE.0)GOTO3H00=0.03W1=0.0DO4I=1,N4W1=W1+P(

22、I)C3=(DX*W1)/G0DW=1.-C3C不是很懂CALLVI(N,DX,P,V)HMIN=1,E3DO30I=1,NH0=0.5*X(I)*X(I)+V(I)C和膜厚公式有关IF(H0.LT.HMIN)HMIN=H0H(I)=H030CONTINUEIF(KK.NE.0)GOTQ2KK=1DH=0.005*HM0C逐步减小法H00=-HMIN+HM032IF(DW.LT.0.0)H00=H00+DHIF(DW.GT.0.0)H00=H00-DHDO60I=1,NH(I)=H00+H(I)EDA(I)=EXRA1*(-1.+(1.+A2*P(I)*Z)C就是粘压公式RO(I)=(A3+1

23、.35*P(I)/(A3+P(I)EPS(I)=RO(I)*H(I)*3/(ENDA*EDA(I)C求出了雷诺方程中的一副三拉,为后面雷诺方程的迭代计算做根底60CONTINUERETURNENDSUBROUTINHER(N,KK,DX,X,P,H,RO,EPS,EDA,V)C雷诺方程中的迭代DIMENSIOX(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)COMMONOMAK/AK(0:1100)DATAPAI1/0.318309886/DO100K=1,KKD2=0.5*(EPS(1)+EPS(2)D3=0.5*(EPS(2)+EPS(3)DO70I=2,N-1

24、D1=D2D2=D3IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2)D8=RO(I)*AK(0)*PAI1D9=RO(I-1)*AK(1)*PAI1D10=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)*DXP(I)=(D11-D12)*D10IF(P(I).LT.0.0)P(I)=0.0C计算雷诺差分方程个系数,计算节点压力,并保证节点压力不小于0(这里还要彳?细考虑)70CONTINUECALLHREE(N,DX,X,P,H,RO,EP

25、S,EDA,V)调用HREE程重新计算各节点的膜厚,100CONTINUERETURNENDSUBROUTINEI(N,DX,P,V)DIMENSIOP(N),V(N)COMMONOMAK/AK(0:1100)PAI1=0.318309886C=ALOGDX)DO10I=1,NV(I)=0.0DO10J=1,NIJ=IABS(I-J)10V(I)=V(I)+(AK(IJ)+C)*DX*P(J)DOI=1,NENDSUBROUTINEUBAK(MM)ENDSUBROUTINEZ(N,P,POLD)DIMENSIOfP(N),POLD(N)DO10I=1,N10POLD(I)=P(I)RETURN

26、ENDSUBROUTINERROP(N,P,POLD,ERP)DIMENSIOP(N),POLD(N)SD=0.0SUM=0.0DO10I=1,NSD=SD+ABS(P(I)-POLD(I)10SUM=SUM+P(I)ERP=SD/SUMRETURNENDCFIXEDFORMATDEMOPROGRAMAINWRITE*,*)'HELLO'WRITE*,*)1'HELLO'100WRITE*,*)'HELLO'10STOPEND1.2等温点接触弹流润滑数值计算方法与程序ProgrampointehlCOMMON/COM1/ENDA,A1,A2,A3

27、,Z,HM0,DHDATAPAI,Z/3.14159265,0.68/DATAN,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='PRESSURESDATUS=UNKNOWN)MM=N-1A1=ALOG(EDA0)+9.67A2=5.1E-9*PHA3=0.59/(PH*1.E-9)U=EDA0*US/

28、(2.*E1*RX)B=PAI*PH*RX/E1W0=2.*PAI*PH/(3.*E1)*(B/RX)*2ALFA=Z*5.1E-9*A1G=ALFA*E1HM0=3.6*(RX/B)*2*G*0.49*U*0.68*W0*(-0.073)ENDA=12.*U*(E1/PH)*(RX/B)*3WRITE(*,*)N,X0,XE,W0,PH,E1,EDA0,RX,USWRITE(4,*)N,X0,XE,W0,PH,E1,EDA0,RX,USWRITE(*,*)'WAITPLEASE'CALLSUBAK(MM)CALLEHL(N,X0,XE)STOPENDSUBROUTINEEHL

29、(N,X0,XE)DIMENSIONX(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,DHDATAMK,G0/1,2.0943951/CALLINITI(N,DX,X0,XE,XJYPOLD)KK=0CALLHREE(N,DX,KK,H00,G0,X;YRO,EPS,EDA,V)14KK=15CALLITER(N,KK,DX,H00,G0,XYRO,EPS,EDA,V)MK=MK+1CALLERP(N,ER,POLD)WRITE(*

30、,*)'ER=,ERIF(ER.GT.1.E-5.AND,DH,GT.1.E-6)THENIF(MK.GE.20)THENMK=1DH=0.5*DHENDIFGOTO14ENDIFIF(DH.LE.1.E6)WRITE(*,*)'PRESSUREARENOTCONVERGENTCALLOUTPUT(N,DX,X,Y,P)RETURNENDSUBROUTINEERP(N,ERPOLD)DIMENSIONP(N,N),POLD(N,N)ER=0.0SUM=0.0DO10I=1,NDO10J=0,NER=ER+ABS(P(I,J)-POLD(I,J)POLD(I,J)=P(I,J)S

31、UM=SUM+P(I,J)10CONTINUEER=ER/SUMRETURNENDSUBROUTINEINITI(N,DX,X0,XE,XPPOLD)DIMENSIONX(N),Y(N),P(N,N),POLD(N,N)NN=(N+1)/2DX=(XE-X0)/(N-1.)Y0=-0.5*(XE-X0)DO5I=1,N,X(I)=X0+(I-1)*DXY(I)=Y0+(I-1)*DX5CONTINUEDOI=1,ND=1.-X(I)*X(I)DOJ=1,NC=D-Y(J)*Y(J)IF(C.LE.0.0)P(I,J)=0.0IF(C.GT.0.0)P(I,J)=SQRT(C)POLD(I,J)=P(I,J)EnddoEnddoReturnEndSubroutinehREE(n,dx,kk,h00,G0,X,YH,RO,EPS,EDA,V)DIMENSIONX(N),Y(N),O(N,N),H(N,N),RO(N,N),EPS(N,N),EDA(N,N),V(N,N)DATAPAI,PAI1/3.14159265,0.2026423/NN=(N+1)/2CALLVI(N,DX,PV)HMIN=1.E3DO30I=1,NDO30J=1,NRAD=X(I)*X(I)+Y(J)*Y(J)W1=0.5*RADH0=W1+V(I,J)IF(H0.LT.HMIN)H

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论