摩擦学原理第9章数值方法ppt课件.ppt

第九章润滑分析常用数值方法 commonnumericalmethodsforlubricationanalysis 各种流体润滑问题都涉及到在狭小间隙中的流体粘性流动 描写这种物理现象的基本方程为Reynolds方程 它的普遍形式是 这个椭圆型的偏微分方程仅仅对于特殊的间隙形状才可能求得解析解 而对于复杂的几何形状或工况条件下的润滑问题 无法用解析方法求得精确解 数值法是将偏微分方程转化为代数方程组的变换方法 它的一般原则是 首先将求解域划分成有限个数的单元 并使每一个单元充分的微小 以至于可以认为在各单元内的未知量 例如油膜压力p 相等或者依照线性变化 而不会造成很大的误差 然后 通过物理分析或数学变换方法 将求解的偏微分方程写成离散形式 即是将它转化为一组线性代数方程 该代数方程组表示了各个单元的待求未知量与周围各单元未知量的关系 最后 根据Gauss消去法或者Gauss Seidel迭代法求解代数方程组 从而求得整个求解域上的未知量 用来求解Reynolds方程的数值方法很多 最常用的是有限差分法 有限元法和边界元法 这些方法都是将求解域划分成许多个单元 但是处理方法各不相同 在有限差分法和有限元法中 代替基本方程的函数在求解域内是近似的 但完全满足边界条件 边界元方法所用的函数在求解域内完全满足基本方程 但是在边界上则近似地满足边界条件 能量方程和弹性变形方程是流体润滑问题中考虑热效应和表面弹性变形时必须求解的重要方程 在本章中也将介绍它们的数值解法 近年发展的多重网格法在润滑计算中开始得到应用 并有广泛的应用前景 本章的最后还介绍了多重网格法求解微分方程和积分方程 9 1Reynolds方程的数值解法numericalmethodsofReynoldsequation9 1 1无量纲化dimensionless 无量纲化也称为或归一化 是将有关方程的变量用合适的参考量 常数 进行转换 使变量的变化范围在1的数量级上 这个参考量称为相对单位 为了判断影响润滑问题诸因素的影响大小 以便抓住主要问题和影响的主要方面 可采用归一化的处理方法 对方程组进行适当的简化 另外 归一化还可以使所得的分析数据具有通用性和广泛性 广义雷诺方程 1 方程的无量纲化 两个因素微分方程中起作用的是变量的变化率 而不是其本身值的大小结果的广泛性 归一化 x h 1 方程的无量纲化 固定瓦径向轴承 广义雷诺方程 式中 不可压缩流体 稳定工况 流体动力粘度为常数 x y z B 不可压缩流体 稳定工况 流体动力粘度为常数 能量方程 无量纲能量方程 取相对单位 温粘方程 压粘方程 2 相似分析 相似是指动力学相似 它包括几何相似运动相似力相似 相似分析主要解决 1 确定相似条件2 导出两个相似问题之间各参数的转换关系 1 相似条件 1 物理模型与控制方程必须是相同的 单值条件 量 使问题有唯一确定解的条件 量 称为单值条件 量 纯粹由单值量组成的相似不变量 是决定问题是否相似的量 因而称之为决定性相似不变量 当相似不变量由同类量组成时称之为简单决定性相似不变量 如d B 当相似不变量由不同类量组成时称之为决定性相似判据 如 不是纯粹由单值量组成的相似不变量称为非决定性相似不变量 判据 2 单值条件 量 相似 由单值量组成的判据相等 2 相似分析的基本步骤 1 列出问题的全部控制方程和单值条件 量 2 选取合适的相对单位 3 化成无量纲形式 4 找出决定性和非决定性相似不变量 3 转换关系 算例 以360 圆柱轴承为例 定常工况 层流 等温 雷诺方程与单值条件 量 1 单值量 轴承直径D 或轴颈直径d 宽度B 半径间隙C 润滑油动力粘度 轴颈旋转角速度 和轴承所受的载荷W 方向沿y轴 方程 边界条件 式中 2 选取无量纲单位并将方程无量纲化 无量纲单位 无量纲方程 边界条件 已知函数 式中 3 找出决定性和非决定性相似不变量 模型实验 决定性 非决定性 理论计算 决定性 非决定性 例如在稳定工况下流体静压润滑的油膜厚度h为常数 若不考虑相对滑动和热效应 则粘度 也是常数 这时Reynolds方程可简化为Laplace方程 即 9 1 将式 9 1 无量纲化 令x XA y YB 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 不可压缩流体 稳定工况 流体动力粘度为常数 能量方程 无量纲能量方程 取相对单位 温粘方程 压粘方程 相似分析 相似是指动力学相似 它包括几何相似运动相似力相似 相似分析主要解决 1 确定相似条件2 导出两个相似问题之间各参数的转换关系 相似条件 1 物理模型与控制方程必须是相同的 单值条件 量 使问题有唯一确定解的条件 量 称为单值条件 量 纯粹由单值量组成的相似不变量 是决定问题是否相似的量 因而称之为决定性相似不变量 当相似不变量由同类量组成时称之为简单决定性相似不变量 如d B 当相似不变量由不同类量组成时称之为决定性相似判据 如 不是纯粹由单值量组成的相似不变量称为非决定性相似不变量 判据 2 单值条件 量 相似 由单值量组成的判据相等 2 相似分析的基本步骤 1 列出问题的全部控制方程和单值条件 量 2 选取合适的相对单位 3 化成无量纲形式 4 找出决定性和非决定性相似不变量 3 转换关系 算例 以360 圆柱轴承为例 定常工况 层流 等温 雷诺方程与单值条件 量 1 单值量 轴承直径D 或轴颈直径d 宽度B 半径间隙C 润滑油动力粘度 轴颈旋转角速度 和轴承所受的载荷W 方向沿y轴 方程 边界条件 式中 2 选取无量纲单位并将方程无量纲化 无量纲单位 无量纲方程 边界条件 已知函数 式中 3 找出决定性和非决定性相似不变量 模型实验 决定性 非决定性 理论计算 决定性 非决定性 9 1 2有限差分法finitedifferencemethod 在流体润滑计算中 有限差分方法的应用最为普遍 现将有限差分法求数值解的步骤和方法说明如下 首先 将所求解的偏微分方程无量纲化 然后 将求解域划分成等距的或不等距的网格 图9 2为等距网格 在X方向有n个节点 Y方向有m个节点 总计n m个节点 图9 2等距网格划分 网格划分的疏密程度根据计算精度要求确定 随着计算机技术的快速发展 在现阶段的润滑计算中 取m 100 500 n 100 500 即可以获得满意的精度 又能够在可以接受的时间内完成 有时 为提高计算精度 可在未知量变化剧烈的区段内细化网格 即采用两种或几种不同间距的分格 或者采用按一定比例递减的分格方法 如果用 代表所求的未知量例如油膜压力p 则变量 在整个域中的分布可以用各节点的 值来表示 根据差分原理 任意节点O i j 的一阶和二阶偏导数都可以由其周围节点的变量值来表示 图9 3差分关系 如图9 3所示 可采用中差分公式 在求解域的边界上或者根据计算要求也可采用前差分公式 或者用后差分公式 通常 中差分的精度最高 若采用下面的中差分表达式 则精度更高 例如以 表示润滑膜压力 将Reynolds方程写成二维二阶偏微分方程的标准形式其中A B C D和E均为已知量 然后 将上述方程应用到各个节点 根据中差分公式 9 5 和 9 6 用差商代替偏导数 即可求得各节点的变量 i j与相邻各节点变量的关系 这种关系可以写成 9 9 式 9 9 中各系数值随节点位置而改变 方程 9 9 是有限差分法的计算方程 对于每个节点都可以写出一个方程 而在边界上的节点变量应满足边界条件 它们的数值是已知量 这样 就可以求得一组线性代数方程 方程与未知量数目相一致 所以可以求解 采用消去法或迭代法求解代数方程组 并使计算结果满足一定的收敛精度 最终求得整个求解域上各节点的变量值 式中l m n 1 流体静压润滑hydrostaticlubrication 对无量纲化的流体静压润滑Reynolds方程式 9 2 将中差分公式 9 5 代入式 9 2 得 9 10 给出边界条件即可由方程 9 9 求得油膜压力分布的数值解 2 流体动压润滑hydrodynamiclubrication 用于不可压缩流体动压润滑轴承的Reynolds方程为 9 11 当h是x y的已知函数时 对于等粘度润滑问题而言 方程 9 11 是线性的 对于变粘度润滑问题 则需要考虑粘度随温度或压力的变化 特别是呈非牛顿性的润滑剂的粘度还受各点速度梯度的影响 则方程 9 11 变成非线性偏微分方程 求解过程较为复杂 1 准二维问题求解在润滑问题的工程计算中 往往采用准二维简化方法 其要点是对于两维变化的油膜压力 预先给定沿某一坐标方向 如轴向 的变化规律 再将这一规律代入二维的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为轴承长度 为偏心率 e c e为偏心距 c为半径间隙 然后 应用差分公式得出式 9 6 形式的计算方程 无量纲化Reynolds方程 9 1 3有限元法与边界元法finiteelementmethodandboundaryelementmethod1 有限元法finiteelementmethod 有限元法是从弹性力学计算中发展起来 继而在流体润滑计算中得到应用的一种数值计算方法 与有限差分法比较 有限元法的主要优点是 适应性强 受几何形状的限制较少 可处理各种定解条件 单元大小和节点可以任意选取 计算精度较高 但是 有限元法计算方程的构成比较复杂 如图9 5所示润滑区域 分成若干个三角形单元 在边界上存在两类边界条件 即 在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 在上述边界条件下的解 这样 有限元法是将不能直接积分求解的二维Reynolds方程转化为求泛函的驻定值 而由式 9 21 建成的计算方程可以求解 通常 有限元法的求解过程可归纳如下 1 将求解域划分成若干三角形或者四边形单元 2 按变分原理写出所求解方程的泛函 3 建立插值函数 即以单元各节点上的变量数值表出单元内任意点的数值 4 根据驻定值条件建立在单元内节点未知量的代数方程组 5 用叠加方法建立总体节点未知量的代数方程组 6 求解代数方程组 2 边界元法boundaryelementmethod 基本特点通过数学方法建立求解域内未知量与边界上未知量之间的关系 这样 只需要将边界划分成若干个单元 求解边界上未知量 进而推算求解域内未知量 所以 边界元法的主要优点是代数方程数很少 同时显著地减少了数据量 尤其是在求解二维和三维问题时更加突出 边界元法的计算流量精度要高于有限元法 并且可以方便地计算混合问题 然而 建立边界元法的计算方程在数学上十分困难 如图9 6所示的阶梯滑块 依润滑膜厚度不同可分为 1 2两部分 每一部分油膜的压力p所遵循的Reynolds方程为 9 22 图9 6阶梯滑块 根据对称性只需分析滑块的一半 即OBCE部分 其总边界s可分成s1和s2两类 s s1 s2 边界条件分别为 在s1边界上p p0 0 在s2边界上 引入一个能满足基本方程 9 20 的权函数P 根据加权余量方法可得其中 经过数学分析 求得本问题的权函数为求解域中任意点的未知量pi与边界上积分的关系为 9 23 同样 边界上任意点的未知量pi与边界上积分的关系为 9 24 以上r为i点至各点的距离 因而P和Q均为已知量 这样 由式 9 24 求得边界上各点未知量后 利用式 9 21 即可计算域内各点未知量 对于点接触弹流问题 如果假定是边界平面上单位面积内的载荷密度 则在该单位面积上有相当于作用了一集中载荷 这时载荷作用点到半无限体表面上一点的距离r为 即 弹流润滑区域中的弹性变形量为 则i点的Z方向上的位移为 式中 1 2分别为接触区表面1和表面2的泊松比 E1 E2分别为接触区表面1和表面2的弹性模量 3 FFT法FastFourierTransformmethod 式中D idx jdy 称为影响函数 而对于下图所示的网格划分情况下 对于任意一个矩形单元内施加单位均布载荷时 接触表面的弹性变形则可写为 式中a xi dx 2 b xi dx 2 c yj dy 2 d yj dy 2 E 为等效弹性模量 如果采用数值积分技术则积分式可转换为离散函数的求和式 对于点接触稳态弹流问题 通常采用有限元法等数值计算方法求出接触表面的弹性变形量 正如前面所述 当有限元法用于计算区域的网格和节点划分的比较多时 就会使计算时间变得过长 同时也占用了更多地计算资源 因此 许多研究人员在有限元法计算方法和改变计算区域的网格划分上做了许多工作 以期减少计算量并提高计算精度 例如自动网格划分技术和无网格计算技术等 但是这些方法并未能从根本上改变有限元法的基本弱点 为了适应点接触弹流润滑分析的需要 采用信号处理的方法可以快速求解这类弹流润滑问题 数字信号处理技术 信号系统分析的主要工作就是在给定系统结构与输入激励的条件下求得系统的输出响应 建立信号系统数学模型的方法可分为输入输出描述法和状态空间描述法两种 而所用的求解方法则有时间域分析法和变换域分析法 输入输出描述法着眼于系统的外部特性 一般不考虑系统的内部变量 直接建立系统输入与输出之间的函数关系 状态空间描述法着眼于系统的内部特性 建立系统内部变量之间以及内部变量与输出变量之间的函数关系 时间域分析法是以时间t或kT T为周期 为变量 直接求解定常系数的线性微分方程式 差分方程式或状态方程式 变换域分析法是应用数学的映射理论 使系统的方程式退化为代数方程式 从而简化了计算 而在各种系统中 线性非时变系统的分析具有非常主要的地位 这不仅是因为实际中的大多数系统都属于或可近似看作线性非时变系统 而且线性非时变系统的分析方法已有了完善的理论 对于随机信号的分析主要采用概率模型的统计方法 即根据输入随机信号的统计特性求得输出随机信号的统计特性 而对于确定信号的分析则主要采用数学模型的解析方法 即建立系统的动态方程式 最后求得系统输出响应的解析描述式 由于摩擦学系统的信号一般都是确定信号 因此就可以采用后一种方法建立和求解摩擦学系统的数学模型 对于任意一个函数 信号 f t 可以在区间 T 2 T 2 内用正弦函数集表示 即 式中 T为任意正整数 系数a0 an与bn为 上式就是任意函数f t 的傅立叶级数展开表达式 如果f t 满足Dirichlet条件 则可以保证存在有f t 的傅立叶级数展开表达式 并且在区间 T 2 T 2 内收敛于原函数 傅立叶级数展开表达式通常还可以写成指数形式 式中 上式为周期信号的富立叶级数指数展开表达式 当信号周期T趋于无穷时 F n 将趋近于无穷小 同时 0 2 T d n 0 nd 积分区间也扩展为 即变为 上式称为非周期信号的频谱密度函数 即频谱函数 它为非周期信号的频域表达式 同样 当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 周期卷积取其中的一个周期即可 即将其写成卷积式 i 0 1 2 N 1 j 0 1 2 M 1 可以看出 求解接触表面的弹性变形的与线性系统的随机函数的富立叶级数指数展开表达式具有相似形式 比较接触表面的弹性变形公式 和富立叶级数指数展开表达式 y n x n h n 对于离散信号则有离散卷积 这种线性卷积在弹性流体润滑计算中可以处理为圆卷积 因此 快速傅立叶变换法 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 采用数字信号处理的方法对接触问题的弹性变形和压力分布求解的一个关键问题是信号源与噪音的处理 正如前面所介绍的那样 由于接触变形和压力分布是非周期信号 因此必须对它作周期性延拓 然后进行快速傅立叶变换和快速傅立叶反变换的信号处理 最后截取有效的区域数值 在快速傅立叶反变换中将充零区域中的无关量舍去后 就可得到原来润滑区域中的实际变形量 由于FFT 快速傅立叶变换 和IFFT 快速傅立叶反变换 的运算次数仅是NLogN次的量级 N为所划分的节点总数 而用有限元法的运算次数是N2的量级 因此用FFT和IFFT求解弹流润滑可比其它数值算法 如有限元法和有限差分法快得多 9 1 4数值解法的其它问题otherproblemsofnumericalmethod1 参数变换parametertransformation 当径向轴承的偏心率较大 例如 0 8 或者楔形滑块具有较大的倾斜角时 使得最小油膜厚度hmin值很小 而在hmin附近膜厚的变化率dh dx数值很高 这就造成在hmin附近很窄的区间内油膜压力急剧地变化 此时 除非采用非常细密的网格 否则计算结果严重失真 而很密的网格又将使计算工作量增加 为了克服这一困难 可以进行参数变换 将M ph3 2代入Reynolds方程直接求解变量M 然后再根据变换关系计算出p的数值 当hmin数值很小时 在它附近的p值剧增 如果以M为变量 M的数值不至于变得很大 同时 M值在hmin附近的变化也较平缓 所以 采用M作为变量可以获得较高的计算精度 3 计算公式与多元回归法calculationformulaandMultivariantregression 数值方法的优点是对于复杂的问题能够给出较准确的解 这对于某些重要的设计和理论研究无疑是有效的手段 然而 数值方法在使用上的缺点是它求解的是个别的具体算例 缺乏一般的通用性 为了增加数值解的通用性 可以将若干组计算数据采用多元回归方法归纳成计算公式 首先 列出影响某一性能P的各个相关参数A B C D 然后 根据经验资料选择适当的函数表示各个相关参数与该性能的关系 通常采用指数函数的形式 即 9 31 最后 根据一组数量足够的 例如500个 理论计算或者实验测量的数据 采用多元回归法确定式 9 26 中的常数K和指数a b c d 等的数值 显然 这样确定的计算公式不可能十分准确地符合全部数据 而只能是具有一定的置信度 同时 还必须进行反复试算和修改才能得到满意的结果 9 2能量方程的数值解法numericalmethodforsolvingenergyequation 为了求得润滑膜中的温度分布需要求解能量方程 而在推导能量方程之前 还必须讨论润滑膜的发热和散热方式 9 2 1传导与对流散热润滑膜中热量的散失主要通过两种途径 1 沿油膜厚度方向 Z向 通过固体表面的传导散热 2 沿油膜尺寸方向 X和Y向 由润滑剂的流动而产生的对流散热 这两种散热方式所散失热量的相对比例因润滑条件而不同 现在以两块无限长平行平板间的油膜散热情况来分析它们之间的关系 如图9 8 图9 8散热分析假设静止板1的温度按线性分布 两端的温度分别为T0和T1 运动板2的表面温度保持为T0 因而温升为 T T1 T0 两平板宽度为B 其间充满润滑油 膜厚为h 现分析单位长度上通过传导和对流的散热量 1 传导散热量Hd 9 32 2 对流散热量Hv若qx为润滑油沿X方向单位长度上的容积流量 为润滑油密度 c为润滑油比热容 而油膜的平均温升为 T 2 则得 9 32 将传导散热与对流散热的比值称为Peclet数 即Peclet数 9 33 由此可见 Peclet数可用来表征润滑系统的散热情况 9 3弹性流体动压润滑数值解法numericalmethodofelastic hydrodynamiclubrication 9 3 1线接触弹流的数值解法EHLinlinecontacts1 基本方程basicequation弹流润滑计算需要同时求解以下方程组 即 1 Reynolds方程 9 41 式中 平均速度U u1 u2 2 h 和 均为x的函数 边界条件 油膜起点x x1处 p 0 油膜终点x x2处 这里 x1应根据润滑油的供应充足程度选取 通常x1 5 15 b x2为在出口区油膜自然破裂的边界 其数值在求解过程中确定 2 膜厚方程filmthicknessequation如图9 12所示 弹性圆柱接触时任意点x处的油膜厚度的表达式为 9 41 式中 hc为没有变形时的中心膜厚 R为当量曲率半径 v x 为由于压力产生的弹性变形位移 图9 12间隙形状图9 13弹性变形 3 弹性变形方程elasticdeformation对于线接触问题 接触体的长度和曲率半径远大于接触宽度 可以认为属于平面应变状态 相当于平直的弹性半无限体受分布载荷作用 如图9 13 根据弹性力学有关理论推导出表面上各点沿垂直方向的弹性位移为 9 43 其中 s是x轴上的附加坐标 它表示任意线载荷p s ds与坐标原点的距离 p s 为载荷分布函数 s1和s2为载荷p x 的起点和终点坐标 E 为当量弹性模量 c为待定常数 4 粘度 压力关系式通常采用Barus公式粘压即 9 44 5 密度 压力关系式采用根据实验曲线整理而得到的密度一压力关系式 9 45 2 Reynolds方程的求解solutionofReynoldsequation 粘压效应和弹性变形对弹流问题中的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 3 弹性变形方程的求解法elasticdeformationsolution 如果给定压力分布p x 计算方程 9 35 的积分即可得到各点的变形量v x 然而 变形方程中的积分部分是一个奇异积分 奇点为s x 此处被积函数无意义 这是弹性变形计算的困难之一 避免奇异积分的一种简便方法是采取分段积分 由于被积函数在除s x点之外都是连续的 所以可以将积分近似地处理为 这样 把奇点排除在积分区间之外 然而 这种方法的困难是如何恰当地确定 x的数值 如果选择不当将产生相当大的计算误差 另 种克服奇异积分的办法是将连续分布的压力p x 进行离散处理 其要点如下 将整个积分区间 x1 x2 划分成若干子区间 并把每一个子区间上的压力分布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 的新 旧值并使它们满足收敛精度 顺解法简单易行 通常适用于中轻载荷条件 而在重载接触时难以达到满意的精度 甚至得不到收敛解 Dowson等人提出逆解法求解线接触弹流润滑 按给定的压力p x 求得各点的dp dx值 此时Reynolds方程成为含h x 三次方的代数方程 解此方程可求得一条膜厚曲线h x 再根据压力p x 用弹性变形方程得到变形v x 从而求得另一条膜厚曲线h x 比较这两条膜厚曲线 按偏差来修正p x 以达到收敛精度 逆解法容易求得收敛结果 但是不易掌握 通常适用于重载荷弹流润滑计算 9 3 2点接触弹流的数值解法EHLinpointcontacts 点接触的一般情况是两个椭圆体相接触而形成椭圆接触区 它比线接触问题复杂得多 因此发层也较缓慢 一直到1965年Archard和Cowking对于圆接触弹流润滑提出了第一个 型的近似解 1970年 Cheng 郑绪云 对于椭圆接触弹流得出了 型解 1976年以后 Hamrock和Dowson对椭圆接触问题进行了系统的数值计算 并提出最小油膜厚度的计算公式 随后温诗铸与朱东提出了椭圆接触弹流润滑的完全数值解 这里 就求解中的主要问题作简要地介绍 1 Reynolds方程求解solutionofReynoldsequation 在一般情况下 接触点的表面速度不一定与接触区椭圆的主轴方向相吻合 此时Reynolds方程应写成 9 52 选取如图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和与线接触弹流的情况相同 为了便于求解应对Reynolds方程进行参数变换 令诱导压力q x y 为则代入方程 9 52 得 9 53 上式是考虑粘压效应的二维Reynolds方程 2 弹性变形方程求解elasticdeformationsolution 根据弹性力学可知 弹性表面上的分布压力p x y 在表面上各点产生的变形位移 x y 用下列关系表示 9 54 如图9 15所示 和 为对应于x和y的附加坐标 为求解域 式 9 54 中积分式的分母部分表示压力作用点 与要计算变形量的点 x y 之间的距离 显然 当x y 时 上述积分是奇异的 克服奇异积分我们采取的办法是 先将坐标原点平移到 式 9 54 变为 然后作极坐标变换 x rcos y rsin 则上式可得到用三角函数表示的积分结果 弹性变形计算的另一个困难是计算工作量过大 如果采用通常的数值积分方法 则在每一次迭代中由于分布压力p x y 的不同 都需要计算每一个节点的变形 而计算每一个节点的变形又必须对整个求解域计算一遍积分 这样 计算工作量太大 而且所需要占用的计算机存储单元也很多 为了克服这一困难 十分有效的办法是采用变形 柔度 矩阵 如果将求解域划分成网格 在x方向共m个节点 y方向共n个节点 设 i j 为在x方向编号为i而在y方向编号为j的节点 并且定义为在节点 i j 处有单位节点压力作用而其余节点上压力为零时 在节点 k l 上产生的变形量 于是 弹性变形方程的离散形式为 9 55 其中 k1为节点 k l 处的弹性变形量 pij为节点 i j 处的节点压力 这样 只需一次算出全部并存储起来 在迭代过程中反复计算变形时即可代入式 9 55 而不必重复计算数值积分 从而大量地减少运算工作量 然而 变形矩阵的元素共有 m n 2个 因而占用的存储单元太多 为了节约存储量 可在y方向采用等距网格 此时式中 s j l 1 故只需计算并存储共计 m n m个量 再由上式可推算全部 当然 如果在x方向也采用等距网格 只需计算和存储的元素将进一步减少至 m n 个 但会导致计算精度降低 当变形计算完成以后 代入膜厚方程 9 56 式中 Rx Ry分别为沿x y方向的当量曲率半径 再将式 9 48 代入Reynolds方程即可求解 3 Hamrock Dowson点接触膜厚公式filmthicknessformulainpointcontact Hamrock和Dowson对等温点接触弹流润滑进行了系统的数值分析 并提出了以下的油膜厚度计算公式 即Hamrock Dowson公式 9 57 9 58 式 9 57 和 9 58 中括弧内因子用以考虑端泄影响 它的大小与椭圆率k有关 椭圆率可以按下式近似计算当其它参数保持不变时 由Hamrock Dowson公式算得的油膜厚度随椭圆率的增加而迅速增大 但当k 5时 油膜厚度随k值的变化就很小 此时 由式 9 57 计算的点接触最小膜厚度和由式 9 48 求得的线接触最小膜厚度基本相同 而由式 9 58 计算得的点接触中心膜厚hc值与 公式算得的入口处的油膜厚度h0也基本相同 由此可知 对于椭圆率k 5的椭圆接触的弹流油膜厚度可以近似地采用线接触膜厚公式进行计算 9 4超薄气体润滑数值numericalmethodforsolvingultrathingaslubrication9 4 1硬盘驱动器中的稀薄气体润滑rarifiedgaslubricationinharddiskdrives 近年来气浮轴承在超精密加工中的应用愈加广泛 如计算机中用于支承高速磁头和磁盘的气膜润滑问题 并形成了润滑技术的一个新的分支 气体薄膜润滑技术 使得气体润滑技术向微观世界 向分子润滑技术迈进 图9 18磁头 磁盘的各表面层及其间的空气间隙组成 计算机硬盘驱动器的磁头 磁盘系统是将磁头连接在悬臂梁末端 依靠高速旋转磁盘产生的气浮力使磁头悬浮于磁盘表面上工作 如图9 18所示 磁头 磁盘界面由磁头 磁盘的各表面层及其间的空气间隙组成 从而形成一个膜厚很小的气体轴承 在这种结构中 由于悬臂梁上磁头的上下动态振动和磁头润滑表面的阶梯或坡度 在磁盘高速运转过程中产生楔形气体压力 支撑磁头悬浮在磁盘表面 9 4 2一般气体润滑控制方程generalgaslubricationgoverningequation 普通气体润滑方程是本问题的计算基础 当磁头飞行高度达到纳米级时 磁头超薄气膜的Knudsen数Kn 分子平均自由程 与气膜厚度 的比值 Kn 可达10左右 此时 固体与流体边界被认为将出现滑流 因此需要修正气膜间隙的气体流量 Fukui和Kanedo使用线型化玻尔兹曼方程组分析了气膜的润滑 得到了修正的雷诺方程 9 59 式中 称为轴承数 9 4 3超薄气膜润滑方程的数值解技术处理numericalsolutiontechniqueforultrathingaslubricationequation 对于流体润滑分析计算 目前工程中广泛采用有限差分法 有限元法和控制体积法 以下将采用编程简单的有限差分法对控制方程进行数值求解 但由于所求解的问题不同于传统的流体润滑问题 出现的新问题必须采用一些特殊处理方法加以解决 1 大轴承数问题superlargebearingnumberproblem 由于现在所要解决的是膜厚为十几纳米的问题 由式 9 59 可知轴承数 与最小膜厚h0的二次方成反比 所以通常的磁头飞行润滑问题时 轴承数都很大 在几十万 甚至上百万的量级上 因此 对方程 9 59 来说 含有轴承数 的剪切流项 在h0很小的时候 其数值远远大于其它的两个压力流项 如果还是按传统的润滑计算 把其作为差分的辅助项计算时 实际上是用大数修正小数 会对较小的压力带来很大的修正量 从而使得迭代过程容易出现失稳 现在我们需要求解的式 9 59 中 由于气体的可压缩性 使得剪切流项内包含了求解变量压力p 这使得该方程具有传输方程的基本形式 并不是标准的传输方程 这为求解提供了新的可能 9 63 根据以上分析 特别将方程 9 59 转换成式 9 63 的形式 在式 9 63 中将气体润滑方程的剪切流项写在等式的左端 其目的是要着重指出 可以利用气体润滑方程的剪切流中含有压力 对其进行主元求解 从而可能避免求解过程中因轴承数过大而发生计算发散 2 阶梯处膜厚突然剧烈变化 abruptvariationoffilmthickness 问题 求解时还必须处理以下两个问题 1 由于磁头的结构所决定 存在数个阶梯 在阶梯处因膜厚的突然变化 使得压力会产生剧烈变化 同时从物理意义上讲 对一个控制体来说 流入的质量流量必须与流出的质量流量相等是保证阶梯处求得正确解的根本依据 2 对于负压型磁头结构 由于磁头的中部会产生很大的负压区 绝对压力可接近0 在求解过程中由于迭代时的压力肯定会出现偏差 如果处理不好 可能会使采用Fukui流量系数时发生计算溢出 黄平等人根据热传导方程求解时 处理两种传导系数差异很大两界面时采用的一个既简单又实用的办法 提出了下面的处理磁头压力计算时 存在突变膜厚的H H算法 折算流量法 9 65 式中 式 9 65 对于膜厚没有突然变化的情况同样适用 因此可以在整个磁头的求解区域上使用 它很好地解决了膜厚突变带来的不稳定 流量不连续导致压力偏低和负压区难收敛等问题 3 计算区域的网格划分griddivisionofcalculationdomain 气膜轴承区域展开图及其求解区域网格划分见图9 20 以往方法为提高计算精度 在膜厚和压力变化剧烈区段内细化网格 但是 如果没有对阶梯部分进行流量折算的处理 即使网格划分的再细 对突变阶梯解也没有任何帮助 因为 前面介绍的H H算法能较好地解决膜厚突变 所以对网格的细分要求已经不高了 本文采用等距网格对整个计算区域进行网格划分 可以适当增加网格数来提高计算精度 对于不太复杂的磁头结构可取m和n为81左右的数目 对于复杂磁头结构一般要取161以上 本文对二维问题计算所使用的网格划分为m n 161 数值计算表明此网格数已能满足计算精度的要求 能够清晰地表达出气膜压力分布轮廓 对二维问题涉及到空气膜厚突变处为曲线或与坐标轴不平行的直线本文采用加密网格的方法来代替 虽对局部的计算压力分布有些影响但对整体计算结果影响不大 图9 20气膜轴承展开图与区域网格划分 9 4 4差分方程和迭代方法differenceequationanditerationmethod 利用迎风格式进行离散的过程与一般的差分形式没有本质上的区别 具体表达如下 如果考虑动态问题 可以把含时间项进行与剪切项一样的方法处理 在上述网格划分的基础上 以差商代替微商 整理后的得到的差分方程为 9 66 采用超松弛迭代法求解 迭代格式为 9 67 这里 和分别是在本轮迭代前和迭代后的压力数值 1是松弛系数 0 1 1 在开始迭代时 由于误差较大 应使用较小的数值 本文一般取0 1 随着迭代次数的增加 它也相应增加至1 收敛条件按下式其中 1为误差收敛精度 计算时常取0 01 0 00001 如果满足收敛条件 则终止迭代 否则 则以当前各点压力值为初始条件 改变迭代因子 1 返回步骤 5 继续迭代 直到满足收敛条件 需要指出 对于残差收敛准则本文采用强收敛条件即绝对误差而不是通常采用的相对误差条件 对于二维情况 复杂的磁头结构迭代的绝对误差会大些在0 01到0 001之间 对于简单的磁头结构绝对误差可达到0 001到0 00001 为了确保计算达到收敛精度 本文计算过程并没有设定误差范围 通过设置大的迭代次数 N 10000 来保证各点处的残差达到稳定 以此来判断是否迭代收敛 9 4 5负压型气体轴承压力计算negativepressuregasbearingcalculation1 磁头结构与压力分布magneticheadstructureandpressuredistribution 表3 3为采用的计算参数 磁头结构的示意图和三维模型如图9 21所示 图9 22分别给出在表3 3的计算参数下两个三维全域压力分布和靠近零压区域的压力分布 表3 3负压型磁头计算参数Tab 3 3Parametersofanegativeslider 图9 21负压型磁头结构示意图 unit mm a 整个区域压力分布 b 靠近零压区域压力分布图9 22负压型磁头全域压力分布图 由图9 22 a 可以看到 在滑块气体轴承靠近后端部的压力比别的部分高很多 主要因为在后端部分为气膜最薄处 在磁头的前端部分 压力随着磁头的阶梯状结构呈明显的层次分布 如果对压力分布图进行局部放大处理如图9 22 b 所示 可以看到在中间区域部分有一个明显的负压区 为了不使磁头承受太大的载荷 在负压型磁头的前端部分设计出一个盆地区域 这样可以在这些区域产生一个明显的负压用来平衡在磁头尾端最薄气膜处产生的很高正压 从而有益于磁头 磁盘系统的稳定性 2 收敛性和误差分析astringencyanderroranalysis 数值结果的收敛性用式 9 68 来判断 图9 23中给出在数值迭代过程中最大误差的收敛过程 纵坐标取全局收敛过程沿X方向的对数值 由图9 23可以看到 在初始迭代过程中 余项快速下降 第一次迭代后 最大的节点误差大概为1 0 1011 在迭代1000次后 余差项缓慢下降 通过改变迭代因子 1 余差项仍保持快速的下降趋势 直到达到所要求的收敛结果 在迭代终止时最大误差大约为6 10 4 图9 23迭代收敛过程图9 24负压型磁头全域误差图 在迭代计算过程中超松弛因子的选择对迭代收敛产生至关重要的作用 初始时因给出的迭代压力为P 1 与实际的气膜压力相差很大 此时应选用小的超松弛因子 1 迭代到一定次数后压力分布已基本趋于稳定此时选用较大的超松弛因子 1 在本文计算过程中 前1000次迭代选用 1 0 1 后期迭代过程 N 1000到N 10000 选用 1 0 9 在迭代过程中出现的尖峰主要就是由于迭代因子的变化引起的 尽管在迭代因子的转换点处 出现余差项增大 但在转换点过后 余差项迅速减小 在达到收敛状态时 余项误差处于随机状态 图9 24给出负压型磁头的全域误差图 可以看到在磁头阶梯结构的过渡处 产生较大的误差 但在迭代后在这些位置的误差已经明显降低到控制的误差范围内 3 与理论解的比较和节点数影响分析 求解控制方程时我们采用简单的有限差分法 进行求解之前要对磁头进行网格划分 从而也就引出多大的网格密度既可以满足计算要求 又可以减少计算机时问题 下面将讨论节点数对计算结果的影响 出于方便的目的 采用一维模型来进行分析 所用的计算参数如表3 5所示 表3 5计算参数Tab 3 5Parametersfortheone dimensionalproblem 图9 25是计算的厚膜图 为了尽可能地与实际相同 磁头的长度和各阶梯的高度与实际磁头基本一致 为了与理论解进行比较 我们首先对高度为10纳米的水平状态进行分析 因为在平行段几乎没有动压效应 加上轴承数很大 该状态的近似解析解为压力与膜厚的乘积等于常数 因此很容易在给出膜厚的条件下将理论解的压力计算出来 见图9 26 需要指出 理论解在尾部边界并不满足边界条件 这是因为 在轴承数很大时 Reynolds方程退化成一阶偏微分方程 无法同时满足两个边界条件 实际上 在尾部将形成边界层使压力迅速回到环境压力 在图9 26中还给出了采用H H算法计算的压力分布 由于两者的解十分接近 所以两组压力线几乎完全重合 只是数值解在尾部要按边界条件收敛到环境压力 因此有一段上升 这就是我们提到的边界层 计算结果的对比说明了 H H算法得到的解是有效的 图9 25膜厚分布图图9 26H H算法的数值解与理论解比较 另外 我们在图9 27给出节点数分别为161和481的压力比较结果 由图中可以看出 除了在气膜分布的过渡处稍有不同外 两者压力分布基本上是一致的 图9 27节点数对压力分布影响比较图 同时也给出总载荷 最大压力和最小压力的计算比较结果 如表9 2所示 由表9 2可以得到所有的误差 节点数增大为原来的三倍时 载荷的误差小于3 最大和最小压力基本不受节点数影响 载荷和压力的差别仅存在于边界层区域 在那里膜厚发生十分巨大的改变 然而这部分区域在整个磁头结构中所占的比例又十分小 因此当节点数达到一定程度时 其变化对总载荷 最大压力和最小压力等气膜特性参数没有明显的影响 表9 2节点数对特性参数的影响 4 一个算例的对比 下面我们将理论解 H H算法 平均差分法和现有国外商业软件计算的同一条件的真实磁头计算中 将几种解做一对比 并通过分析指出 当膜厚发生剧烈变化时 折算算法的有效性和不慎重处理阶梯处的可能带来计算上的错误结果 表9 3计算参数 表9 3是磁头的参数和计算条件 图9 28是采用国外商业软件计算的压力分布 图9 29则是作者采用H H算法得到的压力分布 图9 28国外商业软件结果图9 29H H算法结果 从形式上来看两者的解形状很相象 但是在数值上差异却十分明显 在给定膜厚条件下 用国外商业软件计算得到的总载荷是1 5g 最大压力是不到10个大气压 而用H H算法的到的载荷是4 5g 最大压力约30大气压 为了分析两者的结果正确与否 我们在图9 30给出了磁头中线处的膜厚和对应的线接触问题的解析解 需要指出 由于这是一个复杂膜厚的二维问题 真正意义上的解析解是不存在的 但是由于膜厚非常薄 不仅导致轴承数很大 而且在中心线膜厚很薄的侧流会很小 因此在中心线处的线接触问题的解析解有一定的参考价值 图9 30中线膜厚与线接触理论解 在图9 30中 我们把膜厚