《摩擦学原理(第9章数值方法)ppt课件.ppt》由会员分享，可在线阅读，更多相关《摩擦学原理(第9章数值方法)ppt课件.ppt（173页珍藏版）》请在三一办公上搜索。
1、第九章润滑分析常用数值方法(common numerical methods for lubrication analysis),各种流体润滑问题都涉及到在狭小间隙中的流体粘性流动，描写这种物理现象的基本方程为Reynolds方程，它的普遍形式是：这个椭圆型的偏微分方程仅仅对于特殊的间隙形状才可能求得解析解，而对于复杂的几何形状或工况条件下的润滑问题，无法用解析方法求得精确解。,数值法是将偏微分方程转化为代数方程组的变换方法。它的一般原则是：首先将求解域划分成有限个数的单元，并使每一个单元充分的微小，以至于可以认为在各单元内的未知量（例如油膜压力p）相等或者依照线性变化，而不会造成很大的误差。
2、然后，通过物理分析或数学变换方法，将求解的偏微分方程写成离散形式，即是将它转化为一组线性代数方程。该代数方程组表示了各个单元的待求未知量与周围各单元未知量的关系。最后，根据Gauss消去法或者Gauss-Seidel迭代法求解代数方程组，从而求得整个求解域上的未知量。,用来求解Reynolds方程的数值方法很多，最常用的是有限差分法、有限元法和边界元法，这些方法都是将求解域划分成许多个单元，但是处理方法各不相同。在有限差分法和有限元法中，代替基本方程的函数在求解域内是近似的，但完全满足边界条件。边界元方法所用的函数在求解域内完全满足基本方程，但是在边界上则近似地满足边界条件。,能量方程和弹性变
3、形方程是流体润滑问题中考虑热效应和表面弹性变形时必须求解的重要方程，在本章中也将介绍它们的数值解法。近年发展的多重网格法在润滑计算中开始得到应用，并有广泛的应用前景。本章的最后还介绍了多重网格法求解微分方程和积分方程。,9.1Reynolds方程的数值解法numerical methods of Reynolds equation 9.1.1无量纲化 dimensionless,无量纲化也称为或归一化，是将有关方程的变量用合适的参考量(常数)进行转换，使变量的变化范围在1的数量级上。这个参考量称为相对单位。为了判断影响润滑问题诸因素的影响大小，以便抓住主要问题和影响的主要方面，可采用归一化的处
4、理方法，对方程组进行适当的简化。另外，归一化还可以使所得的分析数据具有通用性和广泛性。,广义雷诺方程,1.方程的无量纲化,两个因素 微分方程中起作用的是变量的变化率，而不是其本身值的大小 结果的广泛性,归一化,x,h,1.方程的无量纲化,固定瓦径向轴承,广义雷诺方程,式中,不可压缩流体,稳定工况,流体动力粘度为常数,x,y,z,B,不可压缩流体,稳定工况,流体动力粘度为常数,能量方程,无量纲能量方程,取相对单位,温粘方程,压粘方程,2.相似分析,相似是指动力学相似，它包括几何相似运动相似力相似,相似分析主要解决：1）确定相似条件2）导出两个相似问题之间各参数的转换关系,1）相似条件,（1）物理
5、模型与控制方程必须是相同的,单值条件（量）：使问题有唯一确定解的条件（量）称为单值条件（量）。纯粹由单值量组成的相似不变量，是决定问题是否相似的量，因而称之为决定性相似不变量。当相似不变量由同类量组成时称之为简单决定性相似不变量，如d/B。当相似不变量由不同类量组成时称之为决定性相似判据，如：。不是纯粹由单值量组成的相似不变量称为非决定性相似不变量（判据）。,（2）单值条件（量）相似，由单值量组成的判据相等,2）相似分析的基本步骤,（1）列出问题的全部控制方程和单值条件（量）,（2）选取合适的相对单位,（3）化成无量纲形式,（4）找出决定性和非决定性相似不变量,3）转换关系,算例,以360圆柱
6、轴承为例，定常工况、层流、等温。雷诺方程与单值条件（量）（1）单值量：轴承直径D（或轴颈直径d）、宽度 B、半径间隙C、润滑油动力粘度、轴颈旋转角速度和轴承所受的载荷W（方向沿y轴）。,方程,边界条件：,式中,（2）选取无量纲单位并将方程无量纲化,无量纲单位,无量纲方程,边界条件：,已知函数：,式中,（3）找出决定性和非决定性相似不变量,模型实验：,决定性,非决定性,理论计算：,决定性,非决定性,例如在稳定工况下流体静压润滑的油膜厚度h为常数，若不考虑相对滑动和热效应，则粘度也是常数。这时Reynolds方程可简化为Laplace方程，即（9.1）将式（9.1）无量纲化，令x=XA，y=YB，
7、A，B为几何尺寸；p=Ppr，pr为油腔压力；=A2/B2；则无量纲化的Reynolds方程为（9.2）,求解方程（9.2）的无量纲化边界条件为（1）在油腔内 P=1（2）在四周边缘上 P=0,若令 代入雷诺方程，得无量纲化基本方程(9.3),例如图9.1所示有限长斜面滑块。,图9.1斜面滑块,这种形式的方程被称为Poisson方程。再将P=Pc（1-Y2）代入方程（9.3），则变换成只含变量Pc和X的方程，即可求解中间断面上无量纲化压力Pc随X的变化。即（9.4）,广义雷诺方程方程的无量纲化,固定瓦径向轴承,式中,不可压缩流体,稳定工况,流体动力粘度为常数,x,y,z,B,不可压缩流体,稳定
8、工况,流体动力粘度为常数,能量方程,无量纲能量方程,取相对单位,温粘方程,压粘方程,相似分析,相似是指动力学相似，它包括几何相似运动相似力相似,相似分析主要解决：1）确定相似条件2）导出两个相似问题之间各参数的转换关系,相似条件,（1）物理模型与控制方程必须是相同的,单值条件（量）：使问题有唯一确定解的条件（量）称为单值条件（量）。纯粹由单值量组成的相似不变量，是决定问题是否相似的量，因而称之为决定性相似不变量。当相似不变量由同类量组成时称之为简单决定性相似不变量，如d/B。当相似不变量由不同类量组成时称之为决定性相似判据，如：。不是纯粹由单值量组成的相似不变量称为非决定性相似不变量（判据）。
9、,（2）单值条件（量）相似，由单值量组成的判据相等,2）相似分析的基本步骤,（1）列出问题的全部控制方程和单值条件（量）,（2）选取合适的相对单位,（3）化成无量纲形式,（4）找出决定性和非决定性相似不变量,3）转换关系,算例,以360圆柱轴承为例，定常工况、层流、等温。雷诺方程与单值条件（量）（1）单值量：轴承直径D（或轴颈直径d）、宽度 B、半径间隙C、润滑油动力粘度、轴颈旋转角速度和轴承所受的载荷W（方向沿y轴）。,方程,边界条件：,式中,（2）选取无量纲单位并将方程无量纲化,无量纲单位,无量纲方程,边界条件：,已知函数：,式中,（3）找出决定性和非决定性相似不变量,模型实验：,决定性,
10、非决定性,理论计算：,决定性,非决定性,9.1.2 有限差分法finite difference method,在流体润滑计算中，有限差分方法的应用最为普遍。现将有限差分法求数值解的步骤和方法说明如下：首先，将所求解的偏微分方程无量纲化。然后，将求解域划分成等距的或不等距的网格，图9.2为等距网格，在X方向有n个节点，Y方向有m个节点，总计nm个节点。,图9.2 等距网格划分,网格划分的疏密程度根据计算精度要求确定。随着计算机技术的快速发展，在现阶段的润滑计算中，取m=100500、n=100500，即可以获得满意的精度，又能够在可以接受的时间内完成。有时，为提高计算精度，可在未知量变化剧烈的
11、区段内细化网格，即采用两种或几种不同间距的分格，或者采用按一定比例递减的分格方法。,如果用代表所求的未知量例如油膜压力p，则变量在整个域中的分布可以用各节点的值来表示。根据差分原理，任意节点O(i,j)的一阶和二阶偏导数都可以由其周围节点的变量值来表示。,图9.3 差分关系,如图9.3所示，可采用中差分公式。在求解域的边界上或者根据计算要求也可采用前差分公式，或者用后差分公式。通常，中差分的精度最高，若采用下面的中差分表达式，则精度更高，例如 以表示润滑膜压力，将Reynolds方程写成二维二阶偏微分方程的标准形式 其中A,B,C,D和E均为已知量。,然后，将上述方程应用到各个节点，根据中差分
12、公式（9.5）和（9.6)用差商代替偏导数，即可求得各节点的变量i,j与相邻各节点变量的关系。这种关系可以写成（9.9）式（9.9）中各系数值随节点位置而改变。,方程（9.9）是有限差分法的计算方程，对于每个节点都可以写出一个方程，而在边界上的节点变量应满足边界条件，它们的数值是已知量。这样，就可以求得一组线性代数方程。方程与未知量数目相一致，所以可以求解。采用消去法或迭代法求解代数方程组，并使计算结果满足一定的收敛精度，最终求得整个求解域上各节点的变量值。,式中 lmn,1流体静压润滑hydrostatic lubrication,对无量纲化的流体静压润滑Reynolds方程式（9.2）将中
13、差分公式（9.5）代入式（9.2）得（9.10）给出边界条件即可由方程（9.9）求得油膜压力分布的数值解。,2流体动压润滑 hydrodynamic lubrication,用于不可压缩流体动压润滑轴承的Reynolds方程为（9.11）当h是x、y的已知函数时，对于等粘度润滑问题而言，方程（9.11）是线性的。对于变粘度润滑问题，则需要考虑粘度随温度或压力的变化，特别是呈非牛顿性的润滑剂的粘度还受各点速度梯度的影响，则方程（9.11）变成非线性偏微分方程，求解过程较为复杂。,（1）准二维问题求解 在润滑问题的工程计算中，往往采用准二维简化方法。其要点是对于两维变化的油膜压力，预先给定沿某一坐
14、标方向（如轴向）的变化规律，再将这一规律代入二维的Reynolds方程即变换为容易求解的一维问题。根据对无限短轴承的分析，油膜压力p沿y方向（即轴向）的分布规律为抛物线 p=pc(1-Y)（9.12）式中，pc为轴向中间断面上的压力；Y为无量纲坐标；为指数，2。,（2）二维问题求解通常的径向滑动轴承设计采用等粘度润滑计算，即假定润滑膜具有相同的粘度，同时认为间隙h只是x的函数，而不考虑安装误差和轴的弯曲变形。图9.4径向轴承展开,将轴承表面沿平面展开，如图9.4，并代入x=R，dx=Rd，则Reynolds方程变为（9.14）若令,（9.15）以上各式中，R为轴承半径；L为轴承长度；为偏心率，
15、=e/c，e为偏心距；c为半径间隙。然后，应用差分公式得出式（9.6）形式的计算方程。,无量纲化Reynolds方程,9.1.3 有限元法与边界元法finite element method and boundary element method 1有限元法 finite element method,有限元法是从弹性力学计算中发展起来，继而在流体润滑计算中得到应用的一种数值计算方法。与有限差分法比较，有限元法的主要优点是：适应性强、受几何形状的限制较少、可处理各种定解条件、单元大小和节点可以任意选取、计算精度较高。但是，有限元法计算方程的构成比较复杂。,如图9.5所示润滑区域，分成若干个三角
16、形单元。在边界上存在两类边界条件。即：在sp边界上压力为已知量，p=p0；在sq边界上流量为已知量，q=q0。图9.5润滑区有限元划分,设在e单元中的压力为pe，则定义e单元的泛函Je为（9.20）这里，A为积分面积；s为积分长度。如果润滑区域共划分为n个单元，各单元泛函的总和应为根据变分原理，泛函存在极值或驻定值的必要条件是它的变分为零，即,（9.21）由欧拉一拉格朗日公式可以证明：符合上述边界条件由Reynolds方程（9.17）求得的解p(x,y)，能够满足泛函驻定值条件（9.21）。反之，由驻定值条件（9.21）求得的解p(x,y)，必然是Reynolds方程（9.19）在上述边界条件
17、下的解。这样，有限元法是将不能直接积分求解的二维Reynolds方程转化为求泛函的驻定值，而由式（9.21）建成的计算方程可以求解。,通常，有限元法的求解过程可归纳如下：（1）将求解域划分成若干三角形或者四边形单元；（2）按变分原理写出所求解方程的泛函；（3）建立插值函数，即以单元各节点上的变量数值表出单元内任意点的数值；（4）根据驻定值条件建立在单元内节点未知量的代数方程组；（5）用叠加方法建立总体节点未知量的代数方程组；（6）求解代数方程组。,2边界元法 boundary element method,基本特点 通过数学方法建立求解域内未知量与边界上未知量之间的关系，这样，只需要将边界划分
18、成若干个单元，求解边界上未知量，进而推算求解域内未知量。所以，边界元法的主要优点是代数方程数很少，同时显著地减少了数据量，尤其是在求解二维和三维问题时更加突出。边界元法的计算流量精度要高于有限元法，并且可以方便地计算混合问题。然而，建立边界元法的计算方程在数学上十分困难。,如图9.6所示的阶梯滑块，依润滑膜厚度不同可分为1、2两部分。每一部分油膜的压力p所遵循的Reynolds方程为（9.22）图9.6阶梯滑块,根据对称性只需分析滑块的一半，即OBCE部分。其总边界s可分成s1和s2两类，s=s1+s2。边界条件分别为：在s1边界上p=p0=0；在s2边界上。引入一个能满足基本方程（9.20）
19、的权函数P，根据加权余量方法可得 其中，,经过数学分析，求得本问题的权函数为求解域中任意点的未知量pi与边界上积分的关系为（9.23）同样，边界上任意点的未知量pi与边界上积分的关系为（9.24）,以上r为i点至各点的距离，因而P和Q均为已知量。这样，由式（9.24）求得边界上各点未知量后，利用式（9.21）即可计算域内各点未知量。,对于点接触弹流问题，如果假定 是边界平面上单位面积内的载荷密度，则在该单位面积上有相当于作用了一集中载荷,这时载荷作用点 到半无限体表面上一点 的距离r为,即：弹流润滑区域中的弹性变形量为：,则i点的Z方向上的位移为,式中1、2分别为接触区表面1和表面2的泊松比，
20、E1、E2分别为接触区表面1和表面2的弹性模量。,3FFT法 Fast Fourier Transform method,式中D(idx,jdy)=称为影响函数。,而对于下图所示的网格划分情况下，对于任意一个矩形单元内施加单位均布载荷时，接触表面的弹性变形则可写为：,式中a=xi-dx/2、b=xi+dx/2、c=yj-dy/2、d=yj+dy/2，E 为等效弹性模量。,如果采用数值积分技术则积分式可转换为离散函数的求和式：,对于点接触稳态弹流问题，通常采用有限元法等数值计算方法求出接触表面的弹性变形量。正如前面所述，当有限元法用于计算区域的网格和节点划分的比较多时，就会使计算时间变得过长，同
21、时也占用了更多地计算资源。因此，许多研究人员在有限元法计算方法和改变计算区域的网格划分上做了许多工作，以期减少计算量并提高计算精度，例如自动网格划分技术和无网格计算技术等。但是这些方法并未能从根本上改变有限元法的基本弱点。为了适应点接触弹流润滑分析的需要，采用信号处理的方法可以快速求解这类弹流润滑问题。,数字信号处理技术,信号系统分析的主要工作就是在给定系统结构与输入激励的条件下求得系统的输出响应。,建立信号系统数学模型的方法可分为输入输出描述法和状态空间描述法两种，而所用的求解方法则有时间域分析法和变换域分析法。输入输出描述法着眼于系统的外部特性，一般不考虑系统的内部变量，直接建立系统输入与
22、输出之间的函数关系。状态空间描述法着眼于系统的内部特性，建立系统内部变量之间以及内部变量与输出变量之间的函数关系。时间域分析法是以时间t或kT（T为周期）为变量，直接求解定常系数的线性微分方程式、差分方程式或状态方程式。变换域分析法是应用数学的映射理论，使系统的方程式退化为代数方程式，从而简化了计算。,而在各种系统中，线性非时变系统的分析具有非常主要的地位。这不仅是因为实际中的大多数系统都属于或可近似看作线性非时变系统，而且线性非时变系统的分析方法已有了完善的理论。对于随机信号的分析主要采用概率模型的统计方法，即根据输入随机信号的统计特性求得输出随机信号的统计特性。而对于确定信号的分析则主要采
23、用数学模型的解析方法，即建立系统的动态方程式，最后求得系统输出响应的解析描述式。由于摩擦学系统的信号一般都是确定信号，因此就可以采用后一种方法建立和求解摩擦学系统的数学模型。,对于任意一个函数（信号）f（t），可以在区间（T/2,T/2）内用正弦函数集表示，即,式中，T为任意正整数。系数a0,an与bn为,上式就是任意函数f（t）的傅立叶级数展开表达式。如果f（t）满足Dirichlet条件，则可以保证存在有f（t）的傅立叶级数展开表达式，并且在区间（T/2,T/2）内收敛于原函数。傅立叶级数展开表达式通常还可以写成指数形式。,式中,上式为周期信号的富立叶级数指数展开表达式，当信号周期T趋于无
24、穷时，Fn将趋近于无穷小，同时02/Td，n0nd，积分区间也扩展为（，）。即变为：,上式称为非周期信号的频谱密度函数，即频谱函数，它为非周期信号的频域表达式。同样，当T趋于无穷时，周期信号变成非周期信号，即,上述两式分别被称为富里叶变换和富里叶逆变换,因此，为了使接触表面的弹性变形公式化成为离散卷积，则需要首先将 和 分别补零，使其长度变成，然后再作周期性(2M-1)(2N-1)-MN延拓得到p(x,y)(2M-1)(2N-1)和D(x,y)(2M-1)(2N-1)，再求p(x,y)(2M-1)(2N-1)和D(x,y)(2M-1)(2N-1)周期卷积取其中的一个周期即可。即将其写成卷积式：
25、(i=0,1,2.N-1;j=0,1,2.M-1),可以看出，求解接触表面的弹性变形的与线性系统的随机函数的富立叶级数指数展开表达式具有相似形式。,比较接触表面的弹性变形公式,和富立叶级数指数展开表达式,y（n）x（n）h（n）,对于离散信号则有离散卷积,这种线性卷积在弹性流体润滑计算中可以处理为圆卷积。因此，快速傅立叶变换法（FFT）可用于计算离散后的线性系统中。这种方法在信号处理中常称为FFT的数字滤波法或快速卷积计算。通常FFT是与频率或周期相关联的，但上式中的变量并不要求是周期函数。离散变量卷积方程要首先通过滤波方法消除其边界效应，然后对式两边进行快速傅立叶变换，则可以将上式写成频域中
26、的各点乘积方程。即：,最后采用快速傅立叶逆变换，可得出时域中的弹性变形量,式中各变数均通过FFT变换得出，即,(i=0,1,2.M-1;j=0,1,2.N-1),(i=0,1,2.M-1;j=0,1,2.N-1),采用数字信号处理的方法对接触问题的弹性变形和压力分布求解的一个关键问题是信号源与噪音的处理。正如前面所介绍的那样，由于接触变形和压力分布是非周期信号，因此必须对它作周期性延拓，然后进行快速傅立叶变换和快速傅立叶反变换的信号处理，最后截取有效的区域数值。,在快速傅立叶反变换中将充零区域中的无关量舍去后，就可得到原来润滑区域中的实际变形量，由于FFT（快速傅立叶变换）和IFFT（快速傅立
27、叶反变换）的运算次数仅是NLogN次的量级，（N为所划分的节点总数）。而用有限元法的运算次数是N2的量级。因此用FFT和IFFT求解弹流润滑可比其它数值算法，如有限元法和有限差分法快得多。,9.1.4 数值解法的其它问题other problems of numerical method 1参数变换parameter transformation,当径向轴承的偏心率较大（例如0.8），或者楔形滑块具有较大的倾斜角时，使得最小油膜厚度hmin值很小，而在hmin附近膜厚的变化率dh/dx数值很高。这就造成在hmin附近很窄的区间内油膜压力急剧地变化。此时，除非采用非常细密的网格，否则计算结果严重
28、失真，而很密的网格又将使计算工作量增加。为了克服这一困难，可以进行参数变换。将M=ph3/2代入Reynolds方程直接求解变量M，然后再根据变换关系计算出p的数值。当hmin数值很小时，在它附近的p值剧增，如果以M为变量，M的数值不至于变得很大。同时，M值在hmin附近的变化也较平缓。所以，采用M作为变量可以获得较高的计算精度。,3计算公式与多元回归法calculation formula and Multi variant regression,数值方法的优点是对于复杂的问题能够给出较准确的解，这对于某些重要的设计和理论研究无疑是有效的手段。然而，数值方法在使用上的缺点是它求解的是个别的具
29、体算例，缺乏一般的通用性。为了增加数值解的通用性，可以将若干组计算数据采用多元回归方法归纳成计算公式。首先，列出影响某一性能P的各个相关参数A、B、C、D。然后，根据经验资料选择适当的函数表示各个相关参数与该性能的关系，通常采用指数函数的形式，即（9.31）,最后，根据一组数量足够的（例如500个）理论计算或者实验测量的数据，采用多元回归法确定式（9.26）中的常数K和指数a,b,c、d等的数值。显然，这样确定的计算公式不可能十分准确地符合全部数据，而只能是具有一定的置信度。同时，还必须进行反复试算和修改才能得到满意的结果。,9.2能量方程的数值解法numerical method for s
30、olving energy equation,为了求得润滑膜中的温度分布需要求解能量方程，而在推导能量方程之前，还必须讨论润滑膜的发热和散热方式。9.2.1 传导与对流散热润滑膜中热量的散失主要通过两种途径：1沿油膜厚度方向（Z向）通过固体表面的传导散热；2沿油膜尺寸方向（X和Y向）由润滑剂的流动而产生的对流散热。这两种散热方式所散失热量的相对比例因润滑条件而不同。现在以两块无限长平行平板间的油膜散热情况来分析它们之间的关系，如图9.8。,图9.8 散热分析 假设静止板1的温度按线性分布，两端的温度分别为T0和T1，运动板2的表面温度保持为T0，因而温升为T=T1-T0。两平板宽度为B，其间充
31、满润滑油，膜厚为h。现分析单位长度上通过传导和对流的散热量。（1）传导散热量Hd（9.32）,（2）对流散热量Hv 若qx为润滑油沿X方向单位长度上的容积流量；为润滑油密度；c为润滑油比热容，而油膜的平均温升为T/2，则得（9.32）将传导散热与对流散热的比值称为Peclet数，即 Peclet数（9.33）由此可见，Peclet数可用来表征润滑系统的散热情况。,9.3 弹性流体动压润滑数值解法numerical method of elastic-hydrodynamic lubrication,9.3.1 线接触弹流的数值解法 EHL in line contacts1基本方程 basic
32、 equation弹流润滑计算需要同时求解以下方程组，即（1）Reynolds方程（9.41）式中，平均速度U=(u1+u2)/2、h、和均为x的函数。边界条件：油膜起点x=x1处，p=0,油膜终点x=x2处，这里，x1应根据润滑油的供应充足程度选取，通常x1=(515)b；x2为在出口区油膜自然破裂的边界，其数值在求解过程中确定。（2）膜厚方程 film thickness equation 如图9.12所示，弹性圆柱接触时任意点x处的油膜厚度的表达式为（9.41）式中，hc为没有变形时的中心膜厚；R为当量曲率半径；v(x)为由于压力产生的弹性变形位移。,图9.12间隙形状 图9.13 弹性
33、变形,（3）弹性变形方程 elastic deformation 对于线接触问题，接触体的长度和曲率半径远大于接触宽度，可以认为属于平面应变状态，相当于平直的弹性半无限体受分布载荷作用，如图9.13。根据弹性力学有关理论推导出表面上各点沿垂直方向的弹性位移为（9.43）其中，s是x轴上的附加坐标，它表示任意线载荷p(s)ds与坐标原点的距离；p(s)为载荷分布函数；s1和s2为载荷p(x)的起点和终点坐标；E为当量弹性模量；c为待定常数。,（4）粘度-压力关系式 通常采用Barus公式粘压即（9.44）（5）密度-压力关系式采用根据实验曲线整理而得到的密度一压力关系式（9.45）,2Reyno
34、lds方程的求解 solution of Reynolds equation,粘压效应和弹性变形对弹流问题中的Reynolds方程求解影响十分显著，这是弹流润滑计算需要特别注意的问题。此外，从弹流润滑的压力分布情况可知：压力p和它的导数值dp/dx都是在很窄的区间内急剧变化的场函数。为了使求解过程稳定，通常须采用参数变换，使变量在求解域内的变化平缓。常用的参数变换是采用诱导压力q(x)，令,并考虑粘压效应，则Reynolds方程变为（9.46）将q(x)解出以后，即可求得p(x)在弹流计算中也可用Vogenpohl变换，即令M(x)=p(x)h(x)3/2，则Reynolds方程为（9.47）
35、,3弹性变形方程的求解法 elastic deformation solution,如果给定压力分布p(x)，计算方程（9.35）的积分即可得到各点的变形量v(x)。然而，变形方程中的积分部分是一个奇异积分，奇点为s=x，此处被积函数无意义。这是弹性变形计算的困难之一。避免奇异积分的一种简便方法是采取分段积分。由于被积函数在除s=x点之外都是连续的，所以可以将积分近似地处理为,这样，把奇点排除在积分区间之外。然而，这种方法的困难是如何恰当地确定x的数值，如果选择不当将产生相当大的计算误差。另种克服奇异积分的办法是将连续分布的压力p(x)进行离散处理。其要点如下：将整个积分区间x1，x2划分成若
36、干子区间，并把每一个子区间上的压力分布p(x)近似地表示为x的多项式，如系数c1、c2和c3则根据已知的节点压力值用待定系数方法确定。例如第i个子区间xi，xi+1上的分布压力为pi(x)=c1i+c2ix+c3ix2。,于是第i个子区间xi，xi+1上的压力pi(x)在各点所产生的变形位移的积分变为 上式中各项的积分可利用关于 以及、的积分公式，这样就可求得Ii的解析解。,4求解方法,求解Reynolds方程有顺解法和逆解法之分，顺解法或称直接迭代法求解弹流问题的顺序如下：采用顺解法求解Reynolds方程是根据给定的h(x)求解p(x)，然后比较p(x)的新、旧值并使它们满足收敛精度。顺解
37、法简单易行，通常适用于中轻载荷条件，而在重载接触时难以达到满意的精度，甚至得不到收敛解。,Dowson等人提出逆解法求解线接触弹流润滑。按给定的压力p(x)求得各点的dp/dx值，此时Reynolds方程成为含h(x)三次方的代数方程，解此方程可求得一条膜厚曲线h(x)。再根据压力p(x)用弹性变形方程得到变形v(x)，从而求得另一条膜厚曲线h(x)，比较这两条膜厚曲线，按偏差来修正p(x)以达到收敛精度。逆解法容易求得收敛结果，但是不易掌握，通常适用于重载荷弹流润滑计算。,9.3.2 点接触弹流的数值解法 EHL in point contacts,点接触的一般情况是两个椭圆体相接触而形成椭
38、圆接触区，它比线接触问题复杂得多，因此发层也较缓慢。一直到1965年Archard和Cowking对于圆接触弹流润滑提出了第一个型的近似解。1970年，Cheng（郑绪云）对于椭圆接触弹流得出了型解。1976年以后，Hamrock和Dowson对椭圆接触问题进行了系统的数值计算，并提出最小油膜厚度的计算公式。随后温诗铸与朱东提出了椭圆接触弹流润滑的完全数值解，这里，就求解中的主要问题作简要地介绍。,1Reynolds方程求解 solution of Reynolds equation,在一般情况下，接触点的表面速度不一定与接触区椭圆的主轴方向相吻合，此时Reynolds方程应写成(9.52)选
39、取如图9.14所表示的坐标轴和求解域，x轴与接触椭圆的短轴相一致。若接触点两表面在x和y方向的速度分量分别为u1、u2和v1、v2，则在x和y方向的平均速度为,求解是从如图所示的矩形求解域上开始进行的。其中AB为入口边，CD为出口边，而AD和BC为端泄边，则、和用来确定求解域边界的位置，通常取=2，=4，而与出口边界有关应在求解过程中确定。图9.14点接触求解域 图9.15弹性变形,求解方程（9.52）的边界条件是：在求解域的入口和端泄边界上压力为零，即当x=-b和y=a时，p=0。而在出口边界x=b上采用Reynolds边界条件，即p=0和 与线接触弹流的情况相同，为了便于求解应对Reyno
40、lds方程进行参数变换。令诱导压力q(x,y)为则代入方程（9.52），得（9.53）上式是考虑粘压效应的二维Reynolds方程。,2弹性变形方程求解 elastic deformation solution,根据弹性力学可知：弹性表面上的分布压力p(x,y)在表面上各点产生的变形位移(x,y)用下列关系表示（9.54）如图9.15所示，和为对应于x和y的附加坐标；为求解域。式（9.54）中积分式的分母部分表示压力作用点（，）与要计算变形量的点（x，y）之间的距离。显然，当x=，y=时，上述积分是奇异的。克服奇异积分我们采取的办法是：先将坐标原点平移到（，），式（9.54）变为,然后作极坐标
41、变换：x=rcos，y=rsin，则上式可得到用三角函数表示的积分结果。弹性变形计算的另一个困难是计算工作量过大。如果采用通常的数值积分方法，则在每一次迭代中由于分布压力p(x,y)的不同，都需要计算每一个节点的变形。而计算每一个节点的变形又必须对整个求解域计算一遍积分。这样，计算工作量太大，而且所需要占用的计算机存储单元也很多。为了克服这一困难，十分有效的办法是采用变形（柔度）矩阵。,如果将求解域划分成网格，在x方向共m个节点，y方向共n个节点。设（i，j）为在x方向编号为i而在y方向编号为j的节点，并且定义为在节点（i，j）处有单位节点压力作用而其余节点上压力为零时，在节点（k,l）上产生
42、的变形量。于是，弹性变形方程的离散形式为（9.55）其中，k1为节点（k,l）处的弹性变形量；pij为节点（i，j）处的节点压力。这样，只需一次算出全部 并存储起来，在迭代过程中反复计算变形时即可代入式（9.55），而不必重复计算数值积分，从而大量地减少运算工作量。,然而，变形矩阵的元素 共有(mn)2个，因而占用的存储单元太多。为了节约存储量，可在y方向采用等距网格，此时式中，s=|j-l|+1。故只需计算并存储 共计（mn）m个量，再由上式可推算全部。当然，如果在x方向也采用等距网格，只需计算和存储的元素将进一步减少至（mn）个，但会导致计算精度降低。当变形计算完成以后，代入膜厚方程（9.
43、56)式中，Rx、Ry分别为沿x、y方向的当量曲率半径。再将式（9.48）代入Reynolds方程即可求解。,3Hamrock-Dowson 点接触膜厚公式 film thickness formula in point contact,Hamrock和Dowson对等温点接触弹流润滑进行了系统的数值分析，并提出了以下的油膜厚度计算公式，即Hamrock-Dowson公式。（9.57）（9.58）,式（9.57）和（9.58）中括弧内因子用以考虑端泄影响，它的大小与椭圆率k有关。椭圆率可以按下式近似计算 当其它参数保持不变时，由Hamrock-Dowson公式算得的油膜厚度随椭圆率的增加而迅速
44、增大。但当k5时，油膜厚度随 k值的变化就很小。此时，由式（9.57）计算的点接触最小膜厚度和由式（9.48）求得的线接触最小膜厚度基本相同；而由式（9.58）计算得的点接触中心膜厚hc值与公式算得的入口处的油膜厚度h0也基本相同。由此可知：对于椭圆率k5的椭圆接触的弹流油膜厚度可以近似地采用线接触膜厚公式进行计算。,9.4 超薄气体润滑数值 numerical method for solving ultra thin gas lubrication 9.4.1 硬盘驱动器中的稀薄气体润滑 rarified gas lubrication in hard disk drives,近年来气浮轴
45、承在超精密加工中的应用愈加广泛，如计算机中用于支承高速磁头和磁盘的气膜润滑问题，并形成了润滑技术的一个新的分支气体薄膜润滑技术，使得气体润滑技术向微观世界，向分子润滑技术迈进。图9.18 磁头/磁盘的各表面层及其间的空气间隙组成,计算机硬盘驱动器的磁头/磁盘系统是将磁头连接在悬臂梁末端，依靠高速旋转磁盘产生的气浮力使磁头悬浮于磁盘表面上工作，如图9.18所示。磁头/磁盘界面由磁头/磁盘的各表面层及其间的空气间隙组成，从而形成一个膜厚很小的气体轴承。在这种结构中，由于悬臂梁上磁头的上下动态振动和磁头润滑表面的阶梯或坡度，在磁盘高速运转过程中产生楔形气体压力，支撑磁头悬浮在磁盘表面。,9.4.2
46、一般气体润滑控制方程 general gas lubrication governing equation,普通气体润滑方程是本问题的计算基础，当磁头飞行高度达到纳米级时，磁头超薄气膜的Knudsen数Kn（分子平均自由程与气膜厚度的比值：Kn=/）可达10左右。此时，固体与流体边界被认为将出现滑流，因此需要修正气膜间隙的气体流量。Fukui和Kanedo使用线型化玻尔兹曼方程组分析了气膜的润滑，得到了修正的雷诺方程：（9.59）,式中,称为轴承数,9.4.3 超薄气膜润滑方程的数值解技术处理numerical solution technique for ultra thin gas lub
47、rication equation,对于流体润滑分析计算，目前工程中广泛采用有限差分法、有限元法和控制体积法。以下将采用编程简单的有限差分法对控制方程进行数值求解。但由于所求解的问题不同于传统的流体润滑问题，出现的新问题必须采用一些特殊处理方法加以解决。,1大轴承数问题 super large bearing number problem,由于现在所要解决的是膜厚为十几纳米的问题，由式(9.59)可知轴承数与最小膜厚h0的二次方成反比，所以通常的磁头飞行润滑问题时，轴承数都很大，在几十万，甚至上百万的量级上。因此，对方程(9.59)来说，含有轴承数的剪切流项，在h0很小的时候，其数值远远大于其
48、它的两个压力流项。如果还是按传统的润滑计算，把其作为差分的辅助项计算时，实际上是用大数修正小数，会对较小的压力带来很大的修正量，从而使得迭代过程容易出现失稳。,现在我们需要求解的式(9.59)中，由于气体的可压缩性，使得剪切流项内包含了求解变量压力p，这使得该方程具有传输方程的基本形式（并不是标准的传输方程），这为求解提供了新的可能(9.63)根据以上分析，特别将方程(9.59)转换成式(9.63)的形式。在式(9.63)中将气体润滑方程的剪切流项写在等式的左端。其目的是要着重指出：可以利用气体润滑方程的剪切流中含有压力，对其进行主元求解，从而可能避免求解过程中因轴承数过大而发生计算发散.,2
49、阶梯处膜厚突然剧烈变化(abrupt variation of film thickness)问题,求解时还必须处理以下两个问题：（1）由于磁头的结构所决定，存在数个阶梯，在阶梯处因膜厚的突然变化，使得压力会产生剧烈变化，同时从物理意义上讲，对一个控制体来说，流入的质量流量必须与流出的质量流量相等是保证阶梯处求得正确解的根本依据；（2）对于负压型磁头结构，由于磁头的中部会产生很大的负压区（绝对压力可接近0），在求解过程中由于迭代时的压力肯定会出现偏差，如果处理不好，可能会使采用Fukui流量系数时发生计算溢出。,黄平等人根据热传导方程求解时，处理两种传导系数差异很大两界面时采用的一个既简单又实
50、用的办法，提出了下面的处理磁头压力计算时，存在突变膜厚的H-H算法（折算流量法）：(9.65)式中，。式（9.65）对于膜厚没有突然变化的情况同样适用，因此可以在整个磁头的求解区域上使用。它很好地解决了膜厚突变带来的不稳定、流量不连续导致压力偏低和负压区难收敛等问题。,3计算区域的网格划分 grid division of calculation domain,气膜轴承区域展开图及其求解区域网格划分见图9.20。以往方法为提高计算精度，在膜厚和压力变化剧烈区段内细化网格。但是，如果没有对阶梯部分进行流量折算的处理，即使网格划分的再细，对突变阶梯解也没有任何帮助。因为，前面介绍的H-H算法能较好
