摩擦学原理(第9章数值方法)课件

第九章

润滑分析常用数值方法

(commonnumericalmethodsforlubricationanalysis)

第九章

润滑分析常用数值方法

(commonn1

各种流体润滑问题都涉及到在狭小间隙中的流体粘性流动，描写这种物理现象的基本方程为Reynolds方程，它的普遍形式是：这个椭圆型的偏微分方程仅仅对于特殊的间隙形状才可能求得解析解，而对于复杂的几何形状或工况条件下的润滑问题，无法用解析方法求得精确解。摩擦学原理(第9章数值方法)ppt课件2数值法是将偏微分方程转化为代数方程组的变换方法。它的一般原则是：首先将求解域划分成有限个数的单元，并使每一个单元充分的微小，以至于可以认为在各单元内的未知量（例如油膜压力p）相等或者依照线性变化，而不会造成很大的误差。然后，通过物理分析或数学变换方法，将求解的偏微分方程写成离散形式，即是将它转化为一组线性代数方程。该代数方程组表示了各个单元的待求未知量与周围各单元未知量的关系。最后，根据Gauss消去法或者Gauss-Seidel迭代法求解代数方程组，从而求得整个求解域上的未知量。数值法是将偏微分方程转化为代数方程组的变换方法。它的一般原则3用来求解Reynolds方程的数值方法很多，最常用的是有限差分法、有限元法和边界元法，这些方法都是将求解域划分成许多个单元，但是处理方法各不相同。在有限差分法和有限元法中，代替基本方程的函数在求解域内是近似的，但完全满足边界条件。边界元方法所用的函数在求解域内完全满足基本方程，但是在边界上则近似地满足边界条件。用来求解Reynolds方程的数值方法很多，最常用的是有限差4能量方程和弹性变形方程是流体润滑问题中考虑热效应和表面弹性变形时必须求解的重要方程，在本章中也将介绍它们的数值解法。近年发展的多重网格法在润滑计算中开始得到应用，并有广泛的应用前景。本章的最后还介绍了多重网格法求解微分方程和积分方程。摩擦学原理(第9章数值方法)ppt课件59.1 Reynolds方程的数值解法

numericalmethodsofReynoldsequation

9.1.1无量纲化

dimensionless无量纲化也称为或归一化，是将有关方程的变量用合适的参考量(常数)进行转换，使变量的变化范围在1的数量级上。这个参考量称为相对单位。

为了判断影响润滑问题诸因素的影响大小，以便抓住主要问题和影响的主要方面，可采用归一化的处理方法，对方程组进行适当的简化。

另外，归一化还可以使所得的分析数据具有通用性和广泛性。9.1 Reynolds方程的数值解法

numerical6广义雷诺方程广义雷诺方程7类型名称及简图特点类型名称及简图特点单油楔固定瓦圆筒轴承(α=360°)结构简单，制造方便，有较大承载能力，但高速稳定性差多油楔固定瓦椭圆轴承流量较大、温升较低。旋转精度和高速稳定性优于单油楔圆轴承但承载能力略有降低工艺性比多油楔轴承好部分瓦轴承(α≤180°)结构简单，制造方便，有较大承载能力。功耗，温升都低于圆筒轴承。高速稳定性差，用于载荷方向基本不变的重载轴承双油楔借位轴承同上，用于单向旋转的轴承浮动环轴承环随轴颈旋转，其转速约为轴颈转速的1/2，润滑油流量大，温升低，高速稳定性好，用于小尺寸高速轻载轴承双向三油楔轴承高速稳定性好，工艺性不如圆筒轴承及椭圆轴承类型名称及简图特点类型名称及简图特点单圆筒轴承(α=360°8类型名称及简图特点类型名称及简图特点多油楔固定瓦多沟轴承结构简单，制造方便，便承载能力低，仅用轻载轴承，高速稳定性略优于圆筒轴承单向三油楔轴承同上。用于单向旋转的轴承螺旋槽轴承利用螺旋的泵入作用和槽面阶梯产生动压承载油膜，温升低，高速稳定性好阶梯轴承同上，承载能力较低，用于小型轴承类型名称及简图特点类型名称及简图特点多多沟轴承结构简单，制造9类型名称及简图特点

多油楔可倾瓦可倾瓦弹性支承轴承高速稳定性较好，特别透用于高速轻载轴承，但工艺性较差

可倾瓦摆动支承轴承同上。但工艺性较好，大、中、小型轴承均适用

多油楔联合轴承动静压联合轴承承载能力大，温升低，功耗小，定心性和稳定性好，特别适于频繁起动的场合，工艺性差，制造较困难但瓦面结构复杂

类型名称及简图特点

多可倾瓦弹性支承轴承高速稳定性较好，特别10类型名称及简图特点类型名称及简图特点

固定瓦推力轴承多油沟推力轴承同多油沟径向轴承。只能在轻载下使用固定瓦斜-平面推力轴承允许轴承有起动载荷

斜面推力轴承用于单向旋转，无起动载荷情况阶梯面推力轴承结构结构简单，用于小尺寸轴承

螺旋槽推力轴承同螺旋槽径向轴承

联合轴承动静压联合推力轴承同动静压联合径向轴承可倾瓦可倾瓦弹性支承推力轴承同可倾瓦弹性支承径向轴承

类型名称及简图特点类型名称及简图特点

固多油沟推力轴承同多油111.方程的无量纲化两个因素微分方程中起作用的是变量的变化率，而不是其本身值的大小结果的广泛性归一化xh1.方程的无量纲化两个因素归一化xh121.方程的无量纲化固定瓦径向轴承广义雷诺方程1.方程的无量纲化固定瓦径向轴承广义雷诺方程13式中

π2π0UVXY式中π2π0UVXY14摩擦学原理(第9章数值方法)ppt课件15不可压缩流体稳定工况不可压缩流体稳定工况16流体动力粘度为常数xyzB流体动力粘度为常数xyzB17不可压缩流体稳定工况流体动力粘度为常数不可压缩流体稳定工况流体动力粘度为常数18能量方程无量纲能量方程取相对单位能量方程无量纲能量方程取相对单位19温粘方程压粘方程温粘方程压粘方程202.相似分析相似是指动力学相似，它包括几何相似运动相似力相似相似分析主要解决：1）确定相似条件2）导出两个相似问题之间各参数的转换关系2.相似分析相似是指动力学相似，它包括相似分析主要解决：211）相似条件（1）物理模型与控制方程必须是相同的单值条件（量）：使问题有唯一确定解的条件（量）称为单值条件（量）。纯粹由单值量组成的相似不变量，是决定问题是否相似的量，因而称之为决定性相似不变量。当相似不变量由同类量组成时称之为简单决定性相似不变量，如d/B。当相似不变量由不同类量组成时称之为决定性相似判据，如：。不是纯粹由单值量组成的相似不变量称为非决定性相似不变量（判据）。（2）单值条件（量）相似，由单值量组成的判据相等1）相似条件（1）物理模型与控制方程必须是相同的单值条件（量222）相似分析的基本步骤（1）列出问题的全部控制方程和单值条件（量）（2）选取合适的相对单位（3）化成无量纲形式（4）找出决定性和非决定性相似不变量3）转换关系2）相似分析的基本步骤（1）列出问题的全部控制方程和单值条件23算例以360°圆柱轴承为例，定常工况、层流、等温。雷诺方程与单值条件（量）（1）单值量：轴承直径D（或轴颈直径d）、宽度B、半径间隙C、润滑油动力粘度μ、轴颈旋转角速度ω和轴承所受的载荷W（方向沿y轴）。方程边界条件：

式中

算例以360°圆柱轴承为例，定常工况、层流、等温。方程边界条24（2）选取无量纲单位并将方程无量纲化

无量纲单位无量纲方程边界条件：已知函数：式中（2）选取无量纲单位并将方程无量纲化无量纲单位无量纲方程25（3）找出决定性和非决定性相似不变量

模型实验：决定性非决定性理论计算：决定性非决定性（3）找出决定性和非决定性相似不变量

模型实验：决定性非决定26例如在稳定工况下流体静压润滑的油膜厚度h为常数，若不考虑相对滑动和热效应，则粘度也是常数。这时Reynolds方程可简化为Laplace方程，即

（9.1）将式（9.1）无量纲化，令x=XA，y=YB，A，B为几何尺寸；p=Ppr，pr为油腔压力；=A2/B2；则无量纲化的Reynolds方程为

（9.2）例如在稳定工况下流体静压润滑的油膜厚度h为常数，若不考虑27求解方程（9.2）的无量纲化边界条件为（1）在油腔内P=1（2）在四周边缘上P=0

求解方程（9.2）的无量纲化边界条件为28若令

代入雷诺方程，得无量纲化基本方程(9.3)

例如图9.1所示有限长斜面滑块。图9.1 斜面滑块若令例如图9.1所示有限长斜面滑块。图9.1 斜面滑块29这种形式的方程被称为Poisson方程。再将P=Pc（1-Y2）代入方程（9.3），则变换成只含变量Pc和X的方程，即可求解中间断面上无量纲化压力Pc随X的变化。即（9.4）这种形式的方程被称为Poisson方程。30广义雷诺方程方程的无量纲化固定瓦径向轴承广义雷诺方程方程的无量纲化固定瓦径向轴承31式中

π2π0UVXY式中π2π0UVXY32摩擦学原理(第9章数值方法)ppt课件33不可压缩流体稳定工况不可压缩流体稳定工况34流体动力粘度为常数xyzB流体动力粘度为常数xyzB35不可压缩流体稳定工况流体动力粘度为常数不可压缩流体稳定工况流体动力粘度为常数36能量方程无量纲能量方程取相对单位能量方程无量纲能量方程取相对单位37温粘方程压粘方程温粘方程压粘方程38相似分析相似是指动力学相似，它包括几何相似运动相似力相似相似分析主要解决：1）确定相似条件2）导出两个相似问题之间各参数的转换关系相似分析相似是指动力学相似，它包括相似分析主要解决：39相似条件（1）物理模型与控制方程必须是相同的单值条件（量）：使问题有唯一确定解的条件（量）称为单值条件（量）。纯粹由单值量组成的相似不变量，是决定问题是否相似的量，因而称之为决定性相似不变量。当相似不变量由同类量组成时称之为简单决定性相似不变量，如d/B。当相似不变量由不同类量组成时称之为决定性相似判据，如：。不是纯粹由单值量组成的相似不变量称为非决定性相似不变量（判据）。（2）单值条件（量）相似，由单值量组成的判据相等相似条件（1）物理模型与控制方程必须是相同的单值条件（量）402）相似分析的基本步骤（1）列出问题的全部控制方程和单值条件（量）（2）选取合适的相对单位（3）化成无量纲形式（4）找出决定性和非决定性相似不变量3）转换关系2）相似分析的基本步骤（1）列出问题的全部控制方程和单值条件41算例以360°圆柱轴承为例，定常工况、层流、等温。雷诺方程与单值条件（量）（1）单值量：轴承直径D（或轴颈直径d）、宽度B、半径间隙C、润滑油动力粘度μ、轴颈旋转角速度ω和轴承所受的载荷W（方向沿y轴）。方程边界条件：

式中

算例以360°圆柱轴承为例，定常工况、层流、等温。方程边界条42（2）选取无量纲单位并将方程无量纲化

无量纲单位无量纲方程边界条件：已知函数：式中（2）选取无量纲单位并将方程无量纲化无量纲单位无量纲方程43（3）找出决定性和非决定性相似不变量

模型实验：决定性非决定性理论计算：决定性非决定性（3）找出决定性和非决定性相似不变量

模型实验：决定性非决定449.1.2有限差分法

finitedifferencemethod

在流体润滑计算中，有限差分方法的应用最为普遍。现将有限差分法求数值解的步骤和方法说明如下：首先，将所求解的偏微分方程无量纲化。然后，将求解域划分成等距的或不等距的网格，图9.2为等距网格，在X方向有n个节点，Y方向有m个节点，总计nm个节点。

9.1.2有限差分法

finitedifference45图9.2等距网格划分图9.2等距网格划分46网格划分的疏密程度根据计算精度要求确定。随着计算机技术的快速发展，在现阶段的润滑计算中，取m=100~500、n=100~500，即可以获得满意的精度，又能够在可以接受的时间内完成。有时，为提高计算精度，可在未知量变化剧烈的区段内细化网格，即采用两种或几种不同间距的分格，或者采用按一定比例递减的分格方法。网格划分的疏密程度根据计算精度要求确定。随着计算机技术的快速47如果用代表所求的未知量例如油膜压力p，则变量在整个域中的分布可以用各节点的值来表示。根据差分原理，任意节点O(i,j)的一阶和二阶偏导数都可以由其周围节点的变量值来表示。i,ji+1,ji-1,ji,j+1i,j-1NESWxOxyy图9.3差分关系如果用代表所求的未知量例如油膜压力p，则变量在整48如图9.3所示，可采用中差分公式。在求解域的边界上或者根据计算要求也可采用前差分公式，或者用后差分公式。 通常，中差分的精度最高，若采用下面的中差分表达式，则精度更高，例如以表示润滑膜压力，将Reynolds方程写成二维二阶偏微分方程的标准形式其中A,B,C,D和E均为已知量。 如图9.3所示，可采用中差分公式。49然后，将上述方程应用到各个节点，根据中差分公式（9.5）和（9.6)用差商代替偏导数，即可求得各节点的变量i,j与相邻各节点变量的关系。这种关系可以写成（9.9）式（9.9）中各系数值随节点位置而改变。然后，将上述方程应用到各个节点，根据中差分公式（9.5）和（50方程（9.9）是有限差分法的计算方程，对于每个节点都可以写出一个方程，而在边界上的节点变量应满足边界条件，它们的数值是已知量。这样，就可以求得一组线性代数方程。方程与未知量数目相一致，所以可以求解。采用消去法或迭代法求解代数方程组，并使计算结果满足一定的收敛精度，最终求得整个求解域上各节点的变量值。

式中l＝m×n方程（9.9）是有限差分法的计算方程，对于每个节点都可以511．流体静压润滑

hydrostaticlubrication

对无量纲化的流体静压润滑Reynolds方程式（9.2）将中差分公式（9.5）代入式（9.2）得

（9.10）

给出边界条件即可由方程（9.9）求得油膜压力分布的数值解。

1．流体静压润滑

hydrostaticlubricati522．流体动压润滑

hydrodynamiclubrication

用于不可压缩流体动压润滑轴承的Reynolds方程为（9.11）当h是x、y的已知函数时，对于等粘度润滑问题而言，方程（9.11）是线性的。对于变粘度润滑问题，则需要考虑粘度随温度或压力的变化，特别是呈非牛顿性的润滑剂的粘度还受各点速度梯度的影响，则方程（9.11）变成非线性偏微分方程，求解过程较为复杂。2．流体动压润滑

hydrodynamiclubrica53（1）准二维问题求解在润滑问题的工程计算中，往往采用准二维简化方法。其要点是对于两维变化的油膜压力，预先给定沿某一坐标方向（如轴向）的变化规律，再将这一规律代入二维的Reynolds方程即变换为容易求解的一维问题。根据对无限短轴承的分析，油膜压力p沿y方向（即轴向）的分布规律为抛物线

p=pc(1-Y)（9.12）式中，pc为轴向中间断面上的压力；Y为无量纲坐标；为指数，＝2。（1）准二维问题求解54（2）二维问题求解通常的径向滑动轴承设计采用等粘度润滑计算，即假定润滑膜具有相同的粘度，同时认为间隙h只是x的函数，而不考虑安装误差和轴的弯曲变形。图9.4 径向轴承展开（2）二维问题求解55将轴承表面沿平面展开，如图9.4，并代入x=R，dx=Rd，则Reynolds方程变为（9.14）若令将轴承表面沿平面展开，如图9.4，并代入x=R，dx=Rd56（9.15）以上各式中，R为轴承半径；L为轴承长度；为偏心率，=e/c，e为偏心距；c为半径间隙。然后，应用差分公式得出式（9.6）形式的计算方程。无量纲化Reynolds方程无量纲化Reynolds方程579.1.3有限元法与边界元法

finiteelementmethodandboundaryelementmethod

1．有限元法

finiteelementmethod

有限元法是从弹性力学计算中发展起来，继而在流体润滑计算中得到应用的一种数值计算方法。与有限差分法比较，有限元法的主要优点是：适应性强、受几何形状的限制较少、可处理各种定解条件、单元大小和节点可以任意选取、计算精度较高。但是，有限元法计算方程的构成比较复杂。9.1.3有限元法与边界元法

finiteelement58如图9.5所示润滑区域，分成若干个三角形单元。在边界上存在两类边界条件。即：在sp边界上压力为已知量，p=p0；在sq边界上流量为已知量，q=q0。

图9.5润滑区有限元划分

摩擦学原理(第9章数值方法)ppt课件59设在e单元中的压力为pe，则定义e单元的泛函Je为

（9.20）

这里，A为积分面积；s为积分长度。如果润滑区域共划分为n个单元，各单元泛函的总和应为根据变分原理，泛函存在极值或驻定值的必要条件是它的变分为零，即设在e单元中的压力为pe，则定义e单元的泛函Je为60

（9.21）由欧拉一拉格朗日公式可以证明：符合上述边界条件由Reynolds方程（9.17）求得的解p(x,y)，能够满足泛函驻定值条件（9.21）。反之，由驻定值条件（9.21）求得的解p(x,y)，必然是Reynolds方程（9.19）在上述边界条件下的解。这样，有限元法是将不能直接积分求解的二维Reynolds方程转化为求泛函的驻定值，而由式（9.21）建成的计算方程可以求解。

摩擦学原理(第9章数值方法)ppt课件61通常，有限元法的求解过程可归纳如下：（1）将求解域划分成若干三角形或者四边形单元；（2）按变分原理写出所求解方程的泛函；（3）建立插值函数，即以单元各节点上的变量数值表出单元内任意点的数值；（4）根据驻定值条件建立在单元内节点未知量的代数方程组；（5）用叠加方法建立总体节点未知量的代数方程组；（6）求解代数方程组。通常，有限元法的求解过程可归纳如下：622．边界元法

boundaryelementmethod

基本特点通过数学方法建立求解域内未知量与边界上未知量之间的关系，这样，只需要将边界划分成若干个单元，求解边界上未知量，进而推算求解域内未知量。所以，边界元法的主要优点是代数方程数很少，同时显著地减少了数据量，尤其是在求解二维和三维问题时更加突出。边界元法的计算流量精度要高于有限元法，并且可以方便地计算混合问题。然而，建立边界元法的计算方程在数学上十分困难。2．边界元法

63如图9.6所示的阶梯滑块，依润滑膜厚度不同可分为1、2两部分。每一部分油膜的压力p所遵循的Reynolds方程为

（9.22）图9.6 阶梯滑块如图9.6所示的阶梯滑块，依润滑膜厚度不同可分为1、2两64根据对称性只需分析滑块的一半，即OBCE部分。其总边界s可分成s1和s2两类，s=s1+s2。边界条件分别为：在s1边界上p=p0=0；在s2边界上。引入一个能满足基本方程（9.20）的权函数P，根据加权余量方法可得

其中，根据对称性只需分析滑块的一半，即OBCE部分。其总边界s可分65经过数学分析，求得本问题的权函数为求解域中任意点的未知量pi与边界上积分的关系为（9.23）同样，边界上任意点的未知量pi与边界上积分的关系为（9.24）以上r为i点至各点的距离，因而P和Q均为已知量。这样，由式（9.24）求得边界上各点未知量后，利用式（9.21）即可计算域内各点未知量。经过数学分析，求得本问题的权函数为以上r为i点至各点的距离，66对于点接触弹流问题，如果假定是边界平面上单位面积内的载荷密度，则在该单位面积上有相当于作用了一集中载荷x,y,zr这时载荷作用点到半无限体表面上一点的距离r为即：弹流润滑区域中的弹性变形量为：则i点的Z方向上的位移为式中ν1、ν2分别为接触区表面1和表面2的泊松比，E1、E2分别为接触区表面1和表面2的弹性模量。3．FFT法

FastFourierTransformmethod对于点接触弹流问题，如果假定是边界平面67式中D(idx,jdy)=称为影响函数。xy计算网格的划分dxdyNM而对于下图所示的网格划分情况下，对于任意一个矩形单元内施加单位均布载荷时，接触表面的弹性变形则可写为：式中a=xi-dx/2、b=xi+dx/2、c=yj-dy/2、d=yj+dy/2，Eˊ＝为等效弹性模量。

如果采用数值积分技术则积分式可转换为离散函数的求和式：式中D(idx,jdy)=xy计68对于点接触稳态弹流问题，通常采用有限元法等数值计算方法求出接触表面的弹性变形量。正如前面所述，当有限元法用于计算区域的网格和节点划分的比较多时，就会使计算时间变得过长，同时也占用了更多地计算资源。因此，许多研究人员在有限元法计算方法和改变计算区域的网格划分上做了许多工作，以期减少计算量并提高计算精度，例如自动网格划分技术和无网格计算技术等。但是这些方法并未能从根本上改变有限元法的基本弱点。为了适应点接触弹流润滑分析的需要，采用信号处理的方法可以快速求解这类弹流润滑问题。数字信号处理技术信号系统分析的主要工作就是在给定系统结构与输入激励的条件下求得系统的输出响应。信号系统的分类时变信号非时变信号时间关系线性信号系统非线性信号系统确定信号随机信号系统特性信号类型信号类型对于点接触稳态弹流问题，通常采用有限元法等数值计算方法求69建立信号系统数学模型的方法可分为输入输出描述法和状态空间描述法两种，而所用的求解方法则有时间域分析法和变换域分析法。输入输出描述法着眼于系统的外部特性，一般不考虑系统的内部变量，直接建立系统输入与输出之间的函数关系。状态空间描述法着眼于系统的内部特性，建立系统内部变量之间以及内部变量与输出变量之间的函数关系。时间域分析法是以时间t或kT（T为周期）为变量，直接求解定常系数的线性微分方程式、差分方程式或状态方程式。变换域分析法是应用数学的映射理论，使系统的方程式退化为代数方程式，从而简化了计算。

而在各种系统中，线性非时变系统的分析具有非常主要的地位。这不仅是因为实际中的大多数系统都属于或可近似看作线性非时变系统，而且线性非时变系统的分析方法已有了完善的理论。对于随机信号的分析主要采用概率模型的统计方法，即根据输入随机信号的统计特性求得输出随机信号的统计特性。而对于确定信号的分析则主要采用数学模型的解析方法，即建立系统的动态方程式，最后求得系统输出响应的解析描述式。由于摩擦学系统的信号一般都是确定信号，因此就可以采用后一种方法建立和求解摩擦学系统的数学模型。

建立信号系统数学模型的方法可分为输入输出描述法和状态70对于任意一个函数（信号）f（t），可以在区间（－T/2,T/2）内用正弦函数集表示，即式中，T为任意正整数。系数a0,an与bn为上式就是任意函数f（t）的傅立叶级数展开表达式。如果f（t）满足Dirichlet条件，则可以保证存在有f（t）的傅立叶级数展开表达式，并且在区间（－T/2,T/2）内收敛于原函数。傅立叶级数展开表达式通常还可以写成指数形式。式中对于任意一个函数（信号）f（t），可以在区间（－T/2,T/71上式为周期信号的富立叶级数指数展开表达式，当信号周期T趋于无穷时，F[n]将趋近于无穷小，同时ω0＝2π/T→dω，nω0→ndω→ω，积分区间也扩展为（－∞，∞）。即变为：上式称为非周期信号的频谱密度函数，即频谱函数，它为非周期信号的频域表达式。同样，当T趋于无穷时，周期信号变成非周期信号，即上述两式分别被称为富里叶变换和富里叶逆变换上式为周期信号的富立叶级数指数展开表达式，当信号周期T趋于无72因此，为了使接触表面的弹性变形公式化成为离散卷积，则需要首先将和分别补零，使其长度变成，然后再作周期性(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)周期卷积取其中的一个周期即可。即将其写成卷积式：(i=0,1,2….N-1;j=0,1,2….M-1)可以看出，求解接触表面的弹性变形的与线性系统的随机函数的富立叶级数指数展开表达式具有相似形式。比较接触表面的弹性变形公式和富立叶级数指数展开表达式y（n）＝x（n）h（n）＝

对于离散信号则有离散卷积因此，为了使接触表面的弹性变形公式化成为离散卷积，则需要首先73这种线性卷积在弹性流体润滑计算中可以处理为圆卷积。因此，快速傅立叶变换法（FFT）可用于计算离散后的线性系统中。这种方法在信号处理中常称为FFT的数字滤波法或快速卷积计算。通常FFT是与频率或周期相关联的，但上式中的变量并不要求是周期函数。离散变量卷积方程要首先通过滤波方法消除其边界效应，然后对式两边进行快速傅立叶变换，则可以将上式写成频域中的各点乘积方程。即：最后采用快速傅立叶逆变换，可得出时域中的弹性变形量式中各变数均通过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)这种线性卷积在弹性流体润滑计算中可以处理为圆卷积。因此，快速74采用数字信号处理的方法对接触问题的弹性变形和压力分布求解的一个关键问题是信号源与噪音的处理。正如前面所介绍的那样，由于接触变形和压力分布是非周期信号，因此必须对它作周期性延拓，然后进行快速傅立叶变换和快速傅立叶反变换的信号处理，最后截取有效的区域数值。

在快速傅立叶反变换中将充零区域中的无关量舍去后，就可得到原来润滑区域中的实际变形量，由于FFT（快速傅立叶变换）和IFFT（快速傅立叶反变换）的运算次数仅是NLogN次的量级，（N为所划分的节点总数）。而用有限元法的运算次数是N2的量级。因此用FFT和IFFT求解弹流润滑可比其它数值算法，如有限元法和有限差分法快得多。

采用数字信号处理的方法对接触问题的弹性变形和压力分布求解的一759.1.4数值解法的其它问题

otherproblems

ofnumericalmethod

1．参数变换parametertransformation

当径向轴承的偏心率较大（例如>0.8），或者楔形滑块具有较大的倾斜角时，使得最小油膜厚度hmin值很小，而在hmin附近膜厚的变化率dh/dx数值很高。这就造成在hmin附近很窄的区间内油膜压力急剧地变化。此时，除非采用非常细密的网格，否则计算结果严重失真，而很密的网格又将使计算工作量增加。为了克服这一困难，可以进行参数变换。将M=ph3/2代入Reynolds方程直接求解变量M，然后再根据变换关系计算出p的数值。当hmin数值很小时，在它附近的p值剧增，如果以M为变量，M的数值不至于变得很大。同时，M值在hmin附近的变化也较平缓。所以，采用M作为变量可以获得较高的计算精度。9.1.4数值解法的其它问题

otherproblems763．计算公式与多元回归法

calculationformulaandMultivariantregression

数值方法的优点是对于复杂的问题能够给出较准确的解，这对于某些重要的设计和理论研究无疑是有效的手段。然而，数值方法在使用上的缺点是它求解的是个别的具体算例，缺乏一般的通用性。为了增加数值解的通用性，可以将若干组计算数据采用多元回归方法归纳成计算公式。首先，列出影响某一性能P的各个相关参数A、B、C、D。然后，根据经验资料选择适当的函数表示各个相关参数与该性能的关系，通常采用指数函数的形式，即（9.31）3．计算公式与多元回归法

calculationform77最后，根据一组数量足够的（例如500个）理论计算或者实验测量的数据，采用多元回归法确定式（9.26）中的常数K和指数a,b,c、d等的数值。显然，这样确定的计算公式不可能十分准确地符合全部数据，而只能是具有一定的置信度。同时，还必须进行反复试算和修改才能得到满意的结果。摩擦学原理(第9章数值方法)ppt课件789.2 能量方程的数值解法

numericalmethodforsolvingenergyequation

为了求得润滑膜中的温度分布需要求解能量方程，而在推导能量方程之前，还必须讨论润滑膜的发热和散热方式。9.2.1传导与对流散热润滑膜中热量的散失主要通过两种途径：1．沿油膜厚度方向（Z向）通过固体表面的传导散热；2．沿油膜尺寸方向（X和Y向）由润滑剂的流动而产生的对流散热。这两种散热方式所散失热量的相对比例因润滑条件而不同。现在以两块无限长平行平板间的油膜散热情况来分析它们之间的关系，如图9.8。9.2 能量方程的数值解法

numericalmethod79

图9.8散热分析假设静止板1的温度按线性分布，两端的温度分别为T0和T1，运动板2的表面温度保持为T0，因而温升为T=T1-T0。两平板宽度为B，其间充满润滑油，膜厚为h。现分析单位长度上通过传导和对流的散热量。（1）传导散热量Hd

（9.32）摩擦学原理(第9章数值方法)ppt课件80（2）对流散热量Hv

若qx为润滑油沿X方向单位长度上的容积流量；为润滑油密度；c为润滑油比热容，而油膜的平均温升为T/2，则得（9.32）将传导散热与对流散热的比值称为Peclet数，即

Peclet数（9.33）

由此可见，Peclet数可用来表征润滑系统的散热情况。

（2）对流散热量Hv819.3 弹性流体动压润滑数值解法

numericalmethodofelastic-hydrodynamiclubrication9.3.1线接触弹流的数值解法EHLinlinecontacts1．基本方程basicequation弹流润滑计算需要同时求解以下方程组，即（1）Reynolds方程（9.41）式中，平均速度U=(u1+u2)/2、h、和均为x的函数。边界条件：油膜起点x=x1处，p=09.3 弹性流体动压润滑数值解法

numericalme82油膜终点x=x2处，这里，x1应根据润滑油的供应充足程度选取，通常x1=(5~15)b；x2为在出口区油膜自然破裂的边界，其数值在求解过程中确定。（2）膜厚方程filmthicknessequation

如图9.12所示，弹性圆柱接触时任意点x处的油膜厚度的表达式为（9.41）式中，hc为没有变形时的中心膜厚；R为当量曲率半径；v(x)为由于压力产生的弹性变形位移。油膜终点x=x2处，83图9.12间隙形状图9.13弹性变形摩擦学原理(第9章数值方法)ppt课件84（3）弹性变形方程elasticdeformation

对于线接触问题，接触体的长度和曲率半径远大于接触宽度，可以认为属于平面应变状态，相当于平直的弹性半无限体受分布载荷作用，如图9.13。根据弹性力学有关理论推导出表面上各点沿垂直方向的弹性位移为（9.43）其中，s是x轴上的附加坐标，它表示任意线载荷p(s)ds与坐标原点的距离；p(s)为载荷分布函数；s1和s2为载荷p(x)的起点和终点坐标；E'为当量弹性模量；c为待定常数。（3）弹性变形方程elasticdeformati85（4）粘度-压力关系式通常采用Barus公式粘压即（9.44）（5）密度-压力关系式采用根据实验曲线整理而得到的密度一压力关系式（9.45）（4）粘度-压力关系式862．Reynolds方程的求解

solutionofReynoldsequation

粘压效应和弹性变形对弹流问题中的Reynolds方程求解影响十分显著，这是弹流润滑计算需要特别注意的问题。此外，从弹流润滑的压力分布情况可知：压力p和它的导数值dp/dx都是在很窄的区间内急剧变化的场函数。为了使求解过程稳定，通常须采用参数变换，使变量在求解域内的变化平缓。常用的参数变换是采用诱导压力q(x)，令

2．Reynolds方程的求解

87并考虑粘压效应，则Reynolds方程变为

（9.46）将q(x)解出以后，即可求得p(x)在弹流计算中也可用Vogenpohl变换，即令M(x)=p(x)[h(x)]3/2，则Reynolds方程为（9.47）并考虑粘压效应，则Reynolds方程变为883．弹性变形方程的求解法

elasticdeformationsolution

如果给定压力分布p(x)，计算方程（9.35）的积分即可得到各点的变形量v(x)。然而，变形方程中的积分部分是一个奇异积分，奇点为s=x，此处被积函数无意义。这是弹性变形计算的困难之一。避免奇异积分的一种简便方法是采取分段积分。由于被积函数在除s=x点之外都是连续的，所以可以将积分近似地处理为3．弹性变形方程的求解法

89这样，把奇点排除在积分区间之外。然而，这种方法的困难是如何恰当地确定x的数值，如果选择不当将产生相当大的计算误差。另－种克服奇异积分的办法是将连续分布的压力p(x)进行离散处理。其要点如下：将整个积分区间[x1，x2］划分成若干子区间，并把每一个子区间上的压力分布p(x)近似地表示为x的多项式，如系数c1、c2和c3则根据已知的节点压力值用待定系数方法确定。例如第i个子区间[xi，xi+1]上的分布压力为pi(x)=c1i+c2ix+c3ix2。这样，把奇点排除在积分区间之外。然而，这种方法的困难是如何恰90于是第i个子区间[xi，xi+1]上的压力pi(x)在各点所产生的变形位移的积分变为上式中各项的积分可利用关于以及、的积分公式，这样就可求得Ii的解析解。于是第i个子区间[xi，xi+1]上的压力pi(x)在各点所914．求解方法求解Reynolds方程有顺解法和逆解法之分，顺解法或称直接迭代法求解弹流问题的顺序如下：采用顺解法求解Reynolds方程是根据给定的h(x)求解p(x)，然后比较p(x)的新、旧值并使它们满足收敛精度。顺解法简单易行，通常适用于中轻载荷条件，而在重载接触时难以达到满意的精度，甚至得不到收敛解。4．求解方法求解Reynolds方程有顺解法和逆解法之分，顺92

Dowson等人提出逆解法求解线接触弹流润滑。按给定的压力p(x)求得各点的dp/dx值，此时Reynolds方程成为含h(x)三次方的代数方程，解此方程可求得一条膜厚曲线h(x)。再根据压力p(x)用弹性变形方程得到变形v(x)，从而求得另一条膜厚曲线h(x)，比较这两条膜厚曲线，按偏差来修正p(x)以达到收敛精度。逆解法容易求得收敛结果，但是不易掌握，通常适用于重载荷弹流润滑计算。Dowson等人提出逆解法求解线接触弹流润滑。按给定的压力939.3.2点接触弹流的数值解法

EHLinpointcontacts点接触的一般情况是两个椭圆体相接触而形成椭圆接触区，它比线接触问题复杂得多，因此发层也较缓慢。一直到1965年Archard和Cowking对于圆接触弹流润滑提出了第一个Грубин型的近似解。1970年，Cheng（郑绪云）对于椭圆接触弹流得出了Грубин型解。1976年以后，Hamrock和Dowson对椭圆接触问题进行了系统的数值计算，并提出最小油膜厚度的计算公式。随后温诗铸与朱东提出了椭圆接触弹流润滑的完全数值解，这里，就求解中的主要问题作简要地介绍。9.3.2点接触弹流的数值解法

941．Reynolds方程求解

solutionofReynoldsequation

在一般情况下，接触点的表面速度不一定与接触区椭圆的主轴方向相吻合，此时Reynolds方程应写成

(9.52)选取如图9.14所表示的坐标轴和求解域，x轴与接触椭圆的短轴相一致。若接触点两表面在x和y方向的速度分量分别为u1、u2和v1、v2，则在x和y方向的平均速度为1．Reynolds方程求解

95求解是从如图所示的矩形求解域上开始进行的。其中AB为入口边，CD为出口边，而AD和BC为端泄边，则、和用来确定求解域边界的位置，通常取=2，=4，而与出口边界有关应在求解过程中确定。

图9.14点接触求解域 图9.15弹性变形

求解是从如图所示的矩形求解域上开始进行的。其中AB为入口边，96求解方程（9.52）的边界条件是：在求解域的入口和端泄边界上压力为零，即当x=-b和y=a时，p=0。而在出口边界x=b上采用Reynolds边界条件，即p=0和与线接触弹流的情况相同，为了便于求解应对Reynolds方程进行参数变换。令诱导压力q(x,y)为则代入方程（9.52），得（9.53）上式是考虑粘压效应的二维Reynolds方程。求解方程（9.52）的边界条件是：在求解域的入口和端泄边界上972．弹性变形方程求解elasticdeformationsolution根据弹性力学可知：弹性表面上的分布压力p(x,y)在表面上各点产生的变形位移(x,y)用下列关系表示

（9.54）如图9.15所示，和为对应于x和y的附加坐标；为求解域。式（9.54）中积分式的分母部分表示压力作用点（，）与要计算变形量的点（x，y）之间的距离。显然，当x=，y=时，上述积分是奇异的。克服奇异积分我们采取的办法是：先将坐标原点平移到（，），式（9.54）变为2．弹性变形方程求解elasticdeformati98然后作极坐标变换：x=rcos，y=rsin，则上式可得到用三角函数表示的积分结果。弹性变形计算的另一个困难是计算工作量过大。如果采用通常的数值积分方法，则在每一次迭代中由于分布压力p(x,y)的不同，都需要计算每一个节点的变形。而计算每一个节点的变形又必须对整个求解域计算一遍积分。这样，计算工作量太大，而且所需要占用的计算机存储单元也很多。为了克服这一困难，十分有效的办法是采用变形（柔度）矩阵。然后作极坐标变换：x=rcos，y=rsin，则99如果将求解域划分成网格，在x方向共m个节点，y方向共n个节点。设（i，j）为在x方向编号为i而在y方向编号为j的节点，并且定义为在节点（i，j）处有单位节点压力作用而其余节点上压力为零时，在节点（k,l）上产生的变形量。于是，弹性变形方程的离散形式为（9.55）其中，k1为节点（k,l）处的弹性变形量；pij为节点（i，j）处的节点压力。这样，只需一次算出全部并存储起来，在迭代过程中反复计算变形时即可代入式（9.55），而不必重复计算数值积分，从而大量地减少运算工作量。

如果将求解域划分成网格，在x方向共m个节点，y方向共n个节点100

然而，变形矩阵的元素共有(mn)2个，因而占用的存储单元太多。为了节约存储量，可在y方向采用等距网格，此时式中，s=|j-l|+1。故只需计算并存储共计（mn）m个量，再由上式可推算全部。当然，如果在x方向也采用等距网格，只需计算和存储的元素将进一步减少至（mn）个，但会导致计算精度降低。当变形计算完成以后，代入膜厚方程（9.56)式中，Rx、Ry分别为沿x、y方向的当量曲率半径。再将式（9.48）代入Reynolds方程即可求解。然而，变形矩阵的元素共有(mn)2个，因而占用的存1013．Hamrock-Dowson点接触膜厚公式

filmthicknessformulainpointcontactHamrock和Dowson对等温点接触弹流润滑进行了系统的数值分析，并提出了以下的油膜厚度计算公式，即Hamrock-Dowson公式。

（9.57）（9.58）3．Hamrock-Dowson点接触膜厚公式

102式（9.57）和（9.58）中括弧内因子用以考虑端泄影响，它的大小与椭圆率k有关。椭圆率可以按下式近似计算当其它参数保持不变时，由Hamrock-Dowson公式算得的油膜厚度随椭圆率的增加而迅速增大。但当k>5时，油膜厚度随k值的变化就很小。此时，由式（9.57）计算的点接触最小膜厚度和由式（9.48）求得的线接触最小膜厚度基本相同；而由式（9.58）计算得的点接触中心膜厚hc值与Грубин公式算得的入口处的油膜厚度h0也基本相同。由此可知：对于椭圆率k>5的椭圆接触的弹流油膜厚度可以近似地采用线接触膜厚公式进行计算。式（9.57）和（9.58）中括弧内因子用以考虑端泄影响，它1039.4 超薄气体润滑数值

numericalmethodforsolvingultrathingaslubrication

9.4.1硬盘驱动器中的稀薄气体润滑

rarifiedgaslubricationinharddiskdrives近年来气浮轴承在超精密加工中的应用愈加广泛，如计算机中用于支承高速磁头和磁盘的气膜润滑问题，并形成了润滑技术的一个新的分支—气体薄膜润滑技术，使得气体润滑技术向微观世界，向分子润滑技术迈进。

图9.18磁头/磁盘的各表面层及其间的空气间隙组成9.4 超薄气体润滑数值

numericalm104计算机硬盘驱动器的磁头/磁盘系统是将磁头连接在悬臂梁末端，依靠高速旋转磁盘产生的气浮力使磁头悬浮于磁盘表面上工作，如图9.18所示。磁头/磁盘界面由磁头/磁盘的各表面层及其间的空气间隙组成，从而形成一个膜厚很小的气体轴承。在这种结构中，由于悬臂梁上磁头的上下动态振动和磁头润滑表面的阶梯或坡度，在磁盘高速运转过程中产生楔形气体压力，支撑磁头悬浮在磁盘表面。

计算机硬盘驱动器的磁头/磁盘系统是将磁头连接在悬臂梁末1059.4.2一般气体润滑控制方程

generalgaslubricationgoverningequation

普通气体润滑方程是本问题的计算基础，当磁头飞行高度达到纳米级时，磁头超薄气膜的Knudsen数Kn（分子平均自由程与气膜厚度的比值：Kn=/）可达10左右。此时，固体与流体边界被认为将出现滑流，因此需要修正气膜间隙的气体流量。Fukui和Kanedo使用线型化玻尔兹曼方程组分析了气膜的润滑，得到了修正的雷诺方程：（9.59）式中称为轴承数9.4.2一般气体润滑控制方程

1069.4.3超薄气膜润滑方程的数值解技术处理

numericalsolutiontechniqueforultrathingaslubricationequation

对于流体润滑分析计算，目前工程中广泛采用有限差分法、有限元法和控制体积法。以下将采用编程简单的有限差分法对控制方程进行数值求解。但由于所求解的问题不同于传统的流体润滑问题，出现的新问题必须采用一些特殊处理方法加以解决。9.4.3超薄气膜润滑方程的数值解技术处理

numeric1071．大轴承数问题superlargebearingnumberproblem

由于现在所要解决的是膜厚为十几纳米的问题，由式(9.59)可知轴承数与最小膜厚h0的二次方成反比，所以通常的磁头飞行润滑问题时，轴承数都很大，在几十万，甚至上百万的量级上。因此，对方程(9.59)来说，含有轴承数的剪切流项，在h0很小的时候，其数值远远大于其它的两个压力流项。如果还是按传统的润滑计算，把其作为差分的辅助项计算时，实际上是用大数修正小数，会对较小的压力带来很大的修正量，从而使得迭代过程容易出现失稳。

1．大轴承数问题superlargebearingn108现在我们需要求解的式(9.59)中，由于气体的可压缩性，使得剪切流项内包含了求解变量压力p，这使得该方程具有传输方程的基本形式（并不是标准的传输方程），这为求解提供了新的可能(9.63)根据以上分析，特别将方程(9.59)转换成式(9.63)的形式。在式(9.63)中将气体润滑方程的剪切流项写在等式的左端。其目的是要着重指出：可以利用气体润滑方程的剪切流中含有压力，对其进行主元求解，从而可能避免求解过程中因轴承数过大而发生计算发散.现在我们需要求解的式(9.59)中，由于气体的可压缩性，使得1092．阶梯处膜厚突然剧烈变化(abruptvariationoffilmthickness)问题求解时还必须处理以下两个问题：（1）由于磁头的结构所决定，存在数个阶梯，在阶梯处因膜厚的突然变化，使得压力会产生剧烈变化，同时从物理意义上讲，对一个控制体来说，流入的质量流量必须与流出的质量流量相等是保证阶梯处求得正确解的根本依据；（2）对于负压型磁头结构，由于磁头的中部会产生很大的负压区（绝对压力可接近0），在求解过程中由于迭代时的压力肯定会出现偏差，如果处理不好，可能会使采用Fukui流量系数时发生计算溢出。2．阶梯处膜厚突然剧烈变化(abruptvariation110黄平等人根据热传导方程求解时，处理两种传导系数差异很大两界面时采用的一个既简单又实用的办法，提出了下面的处理磁头压力计算时，存在突变膜厚的H-H算法（折算流量法）：

(9.65)

式中，。式（9.65）对于膜厚没有突然变化的情况同样适用，因此可以在整个磁头的求解区域上使用。它很好地解决了膜厚突变带来的不稳定、流量不连续导致压力偏低和负压区难收敛等问题。黄平等人根据热传导方程求解时，处理两种传导系数差异很大两界面1113．计算区域的网格划分

griddivisionofcalculationdomain气膜轴承区域展开图及其求解区域网格划分见图9.20。以往方法为提高计算精度，在膜厚和压力变化剧烈区段内细化网格。但是，如果没有对阶梯部分进行流量折算的处理，即使网格划分的再细，对突变阶梯解也没有任何帮助。因为，前面介绍的H-H算法能较好地解决膜厚突变，所以对网格的细分要求已经不高了，本文采用等距网格对整个计算区域进行网格划分，可以适当增加网格数来提高计算精度。对于不太复杂的磁头结构可取m和n为81左右的数目，对于复杂磁头结构一般要取161以上。本文对二维问题计算所使用的网格划分为m=n=161，数值计算表明此网格数已能满足计算精度的要求，能够清晰地表达出气膜压力分布轮廓。对二维问题涉及到空气膜厚突变处为曲线或与坐标轴不平行的直线本文采用加密网格的方法来代替，虽对局部的计算压力分布有些影响但对整体计算结果影响不大。3．计算区域的网格划分

gridd112图9.20气膜轴承展开图与区域网格划分摩擦学原理(第9章数值方法)ppt课件1139.4.4差分方程和迭代方法

differenceequationand

iterationmethod

利用迎风格式进行离散的过程与一般的差分形式没有本质上的区别，具体表达如下。如果考虑动态问题，可以把含时间项进行与剪切项一样的方法处理。在上述网格划分的基础上，以差商代替微商。整理后的得到的差分方程为

(9.66)采用超松弛迭代法求解,迭代格式为(9.67)

9.4.4差分方程和迭代方法

differenceequ114这里，和分别是在本轮迭代前和迭代后的压力数值。1是松弛系数，0<11。在开始迭代时，由于误差较大，应使用较小的数值，本文一般取0.1，随着迭代次数的增加，它也相应增加至1。

收敛条件按下式其中，这里，和分别是在本轮迭代前和迭代后的压力数值。115

1

为误差收敛精度，计算时常取0.01～0.00001。如果满足收敛条件，则终止迭代；否则，则以当前各点压力值为初始条件，改变迭代因子1，返回步骤（5）继续迭代，直到满足收敛条件。需要指出：对于残差收敛准则本文采用强收敛条件即绝对误差而不是通常采用的相对误差条件。对于二维情况，复杂的磁头结构迭代的绝对误差会大些在0.01到0.001之间，对于简单的磁头结构绝对误差可达到0.001到0.00001。为了确保计算达到收敛精度，本文计算过程并没有设定误差范围，通过设置大的迭代次数（N=10000）来保证各点处的残差达到稳定，以此来判断是否迭代收敛。1为误差收敛精度，计算时常取0.01～0.116

9.4.5负压型气体轴承压力计算

negativepressuregasbearingcalculation

1．磁头结构与压力分布

magneticheadstructureandpressuredistribution

表3-3为采用的计算参数。磁头结构的示意图和三维模型如图9.21所示。图9.22分别给出在表3-3的计算参数下两个三维全域压力分布和靠近零压区域的压力分布。参数数值(单位)参数数值(单位)特征膜厚h010×10-9(m)磁头长度L1.25×10-3(m)转动速度U15(m/s)磁头宽度B1.00×10-3(m)俯仰角1.0×10-4(rad)气体粘度1.83×10-5(Pa·s)池深hl1.114×10-7(m)一级阶梯hu1.778×10-6(m)轴承数Λ137250大气压力pa1×105(Pa)表3-3负压型磁头计算参数Tab.3-3Parametersofanegativeslider

9.4.5负压型气体轴承压力计算

negativepr117

图9.21负压型磁头结构示意图(unit:mm)摩擦学原理(第9章数值方法)ppt课件118(a)整个区域压力分布(b)靠近零压区域压力分布

图9.22负压型磁头全域压力分布图摩擦学原理(第9章数值方法)ppt课件119由图9.22(a)可以看到，在滑块气体轴承靠近后端部的压力比别的部分高很多。主要因为在后端部分为气膜最薄处。在磁头的前端部分，压力随着磁头的阶梯状结构呈明显的层次分布。如果对压力分布图进行局部放大处理如图9.22(b)所示，可以看到在中间区域部分有一个明显的负压区。为了不使磁头承受太大的载荷，在负压型磁头的前端部分设计出一个盆地区域，这样可以在这些区域产生一个明显的负压用来平衡在磁头尾端最薄气膜处产生的很高正压，从而有益于磁头/磁盘系统的稳定性。摩擦学原理(第9章数值方法)ppt课件1202．收敛性和误差分析

astringencyanderroranalysis数值结果的收敛性用式(9.68)来判断。图9.23中给出在数值迭代过程中最大误差的收敛过程，纵坐标取全局收敛过程沿X方向的对数值。由图9.23可以看到，在初始迭代过程中，余项快速下降，第一次迭代后，最大的节点误差大概为1.01011。在迭代1000次后，余差项缓慢下降，通过改变迭代因子1

，余差项仍保持快速的下降趋势，直到达到所要求的收敛结果，在迭代终止时最大误差大约为610-4。2．收敛性和误差分析

121

图9.23迭代收敛过程图9.24负压型磁头全域误差图摩擦学原理(第9章数值方法)ppt课件122在迭代计算过程中超松弛因子的选择对迭代收敛产生至关重要的作用，初始时因给出的迭代压力为P＝1，与实际的气膜压力相差很大，此时应选用小的超松弛因子1，迭代到一定次数后压力分布已基本趋于稳定此时选用较大的超松弛因子1。在本文计算过程中，前10