彈流潤滑數(shù)值計算方法-趙江

1、彈流潤滑數(shù)值計算方法1.1等溫線接觸彈流潤滑數(shù)值計算方法與程序等溫點接觸彈流潤滑的根本方程包括:平衡方程,主要形式如下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)計算程序計算工況參數(shù)節(jié)點數(shù)：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接觸應(yīng)力Ph,單位paA1=(ALOG(EDA0)+9.67)C粘度公式計算時的參數(shù)A2=PH/P0C粘度公式計算時的參數(shù)A3=0.59/(PH*1.E-9)C密度公式計算時參數(shù)B=4.*R*PH/E1CB為HERTZ接觸區(qū)半寬,單位mALFA=Z*A1/P0C粘度系數(shù)G=ALFA*E1C材料參數(shù)U=EDA0*US/(2.*E1*R)C蘇聯(lián)人提出的入口區(qū)分析解中定義的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接觸應(yīng)力ph,B為HERTZ接觸區(qū)半寬,RRMK=1DX=(XE-X0)/(N-1.0)CDX等距節(jié)點間距DO10I=1,NX(I)=X0+(I-1)*DXC定位X(I)的位置對應(yīng)于框圖中的計算個節(jié)點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對應(yīng)于木g圖中,調(diào)用HREE子程序計算個節(jié)點膜厚,密度ro,和粘度EDA,C并調(diào)用子程序VI計算個節(jié)點的彈性變形C調(diào)用子程序FZ將P(I)的值賦給POLD,POLD第i個節(jié)點的前一次迭代壓力值14KK=19CALLITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)/粘度/密度C調(diào)用ITER子程序進行雷諾方程迭代,重新計算hree壓力分布下的節(jié)點膜厚MK=MK+1CALLERROP(N,P,POLD,ERP)C調(diào)用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個字符寬輸出指數(shù)類型的浮點數(shù),小數(shù)局部占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、厚公式有關(guān)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計算雷諾差分方程個系數(shù),計算節(jié)點壓力,并保證節(jié)點壓力不小于0(這里還要仔細考慮)70CONTINUECALLHREE(N,DX,X,P,H,RO,EPS,EDA,V)C調(diào)用HREE方程重新計算各節(jié)點的膜厚,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接觸應(yīng)力Ph,單位paA1=(ALOGEDA0)+9.67)C粘度公式計算時的參數(shù)A2=PH/P0C粘度公式計算時的參數(shù)A3=0.59/(PH*1.E-9)C密度公式計算時參數(shù)B=4.*R*PH/E1CB為HERT菽觸區(qū)半寬,單位mALFA=Z*A1/P0C粘度系數(shù)G=ALF

16、A*E1C材料參數(shù)U=EDA0*US/(2.*E1*R)C蘇聯(lián)人提出的入口區(qū)分析解中定義的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接觸應(yīng)力ph,B為HERT接觸區(qū)半寬,RRMK=1DX=(XE-X0)/(N-1.0)CDX等距節(jié)點間距DO10I=1,NX(I)=X0+(I-1)*DXC定位X(I)的位置,對應(yīng)于框圖中的計算個節(jié)點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對應(yīng)于力I圖中,調(diào)用HRE子程序計算個節(jié)點膜厚,密度ro,和粘度EDAC并調(diào)用子程序VI計算個節(jié)點的彈性變形C調(diào)用子程序FZ等P(I)的值賦給POLD,POLDi個節(jié)點的前一次迭代壓力值14KK=19CALLITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)C調(diào)用ITER子程序進行雷諾方程迭代,重新計算hree壓力分布下的節(jié)點膜厚/粘度/

19、密度MK=MK+1CALLERROP(N,P,POLD,ERP)C調(diào)用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環(huán)口最大壓力P2,量綱一化并輸出106CONTINUEH3=H2*B*B/RRP3=P2*PHC這個地方還不是很懂,可能是獲得最小的值然后量綱計算出有單位的值110FORMAT(1X,E12.6)C格式化輸出1片巴輸出位置向右跳過n個位置CE12.6以12個字符寬輸出指數(shù)類型的浮點數(shù),小數(shù)局部占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和膜厚公式有關(guān)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計算雷諾差分方程個系數(shù),計算節(jié)點壓力,并保證節(jié)點壓力不小于0(這里還要彳？細考慮)70CONTINUECALLHREE(N,DX,X,P,H,RO,EP

25、S,EDA,V)調(diào)用HREE程重新計算各節(jié)點的膜厚,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等溫點接觸彈流潤滑數(shù)值計算方法與程序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