二维雷诺方程数值求解：离散化与SOR迭代实现

简介：面向机械润滑与流体力学领域的Matlab程序包，聚焦二维雷诺方程求解与油膜压力分布预测，适用于轴承、机械密封等润滑系统设计及性能优化。压缩包内共2个m文件，文件总大小仅2KB，包含数据处理/参数定义模块与主求解程序，结构紧凑，可快速载入几何参数并运行计算。已有560人学习下载。程序基于粘性流体动压润滑理论，采用有限差分等数值离散方法迭代求解雷诺方程，支持设定不同速度、载荷及边界条件，最终输出油膜压力分布结果。借助该程序，工程师可以直观评估润滑性能、判断承载能力与油膜厚度变化，进而改进结构设计、减少摩擦磨损并延长设备寿命。同时，文件注释清晰、代码体量小，便于初学者理解算法流程，也可作为教学演示或二次开发基础。

1. 程序.rar里的二维雷诺方程：润滑压力分布到底在算什么

从共享盘拷来的程序.rar解压后，里面多半是几个Fortran或C源文件，外加一份格式奇特的DAT输入数据。很多做滑动轴承、齿轮或密封分析的工程师，卡住的不是编译，而是二维雷诺方程本身：膜厚、粘度和速度都给了，油膜压力分布怎么稳定解出来。这个方程描述润滑油膜在x和z两个方向上的压力建立过程，左侧是压力流，右侧是剪切流与挤压效应；输出压力场之后，承载力、摩擦力和泄漏量全部从这里派生。对于要读懂老代码、准备迁移到Python的人来说，先把离散格式和边界条件摸清楚，比移植代码本身更省时间。

2. 二维雷诺方程的离散化：把偏微分方程变成可解的代数方程

2.1 油膜雷诺方程里每一项的物理来源

二维稳态雷诺方程最常见的形式是：

∂/∂x (h³/η · ∂p/∂x) + ∂/∂z (h³/η · ∂p/∂z) = 6U ∂h/∂x

左侧两项代表压力梯度引起的Poiseuille流动在膜厚方向积分后的净流出；等号右侧的6U∂h/∂x是楔形间隙在剪切流动下产生的压力源项。U是运动表面的速度，η是润滑油动力粘度，h是油膜厚度。稳态假设下挤压项12∂h/∂t被丢弃，只适用于转速稳定、不出现大幅振动冲击的工况。

二维体现在压力p同时随x和z变化，而y方向（膜厚方向）因为间隙远小于特征长度，压力梯度被忽略。这个降维假设是雷诺方程能用工程方法求解的前提。程序.rar里的老代码绝大多数也是在这个假设下用有限差分做的，如果把膜厚方向也纳入计算，问题会变成N-S方程或广义Reynolds方程，网格量和收敛难度完全不同。

2.2 中心差分格式与系数组装

对内部节点(i,j)，取x方向步长Δx、z方向步长Δz，用中心差分离散：

对x方向的扩散项，界面处需要h³的值。这里有一个关键选择：用算术平均(hm = (h_i³ + h_{i+1}³)/2)还是调和平均。对连续膜厚，两者差距不大；但遇到阶梯轴承、轴瓦边缘台阶这类膜厚突变，算术平均会在突变位置造成压力过冲，调和平均更稳定。常见的离散方式是：

∂/∂x(h³ ∂p/∂x) ≈ [h_{i+1/2}³(p_{i+1,j}-p_{i,j}) - h_{i-1/2}³(p_{i,j}-p_{i-1,j})] / Δx²

其中 h_{i+1/2}³ = 2 h_{i,j}³ h_{i+1,j}³ / (h_{i,j}³ + h_{i+1,j}³)。右端源项 6U∂h/∂x 用中心差分：6U(h_{i+1,j} - h_{i-1,j}) / (2Δx)。

下面的Python函数演示如何组装稀疏矩阵形式的线性系统。注意，这里直接用节点膜厚做调和平均，避免单独插值：

import numpy as np from scipy.sparse import lil_matrix def assemble_2d_reynolds(h, eta, U, dx, dz): nx, nz = h.shape n = nx * nz A = lil_matrix((n, n)) b = np.zeros(n) def idx(i, j): return i * nz + j for i in range(1, nx - 1): for j in range(1, nz - 1): row = idx(i, j) # 界面膜厚用调和平均，h3 代表膜厚三次方 h3_E = 2 * h[i, j]**3 * h[i, j+1]**3 / (h[i, j]**3 + h[i, j+1]**3 + 1e-30) h3_W = 2 * h[i, j]**3 * h[i, j-1]**3 / (h[i, j]**3 + h[i, j-1]**3 + 1e-30) h3_N = 2 * h[i, j]**3 * h[i-1, j]**3 / (h[i, j]**3 + h[i-1, j]**3 + 1e-30) h3_S = 2 * h[i, j]**3 * h[i+1, j]**3 / (h[i, j]**3 + h[i+1, j]**3 + 1e-30) # 扩散项系数 aE = h3_E / dx**2 aW = h3_W / dx**2 aN = h3_N / dz**2 aS = h3_S / dz**2 aC = aE + aW + aN + aS A[row, idx(i, j+1)] = aE A[row, idx(i, j-1)] = aW A[row, idx(i-1, j)] = aN A[row, idx(i+1, j)] = aS A[row, row] = -aC # 源项：楔形剪切项 b[row] = 6.0 * U * (h[i, j+1] - h[i, j-1]) / (2.0 * dx) / eta # 注意 b 里的粘度放在分母还是分子，取决于方程写法 # 左侧已经含 1/eta 的话，这里就不应再除 # 程序.rar 的老代码里最容易错的就是这个位置 return A.tocsr(), b

这段代码的组装逻辑是逐节点算系数、填入稀疏矩阵。边界节点全部保持为0，对应环境压力边界。注意代码里b[row]的粘度位置容易混淆：如果把方程写成 ∂/∂x(h³∂p/∂x)+... = 6ηU∂h/∂x，左侧就不需要再乘粘度修正，直接对压力求解，b里也就不再除以eta。老程序里常见的单位错乱，十有八九出在这个系数上。

上面的代码只是组装，真正求解还要加边界条件。稀疏矩阵用scipy.sparse.linalg.spsolve直接解即可，但对大网格内存不够，所以后面会回到SOR迭代。

2.3 边界条件：入口、出口、侧边与对称面怎么落到节点上

常见的边界类型有四种。环境压力边界直接置p=0，适合暴露在大气中的侧边和出口；给定压力边界比如供油槽，直接把对应节点赋固定值；对称面用∂p/∂z=0，需要用一阶单侧差分把对称轴上的节点用相邻行表达；出口处如果压力会降到负值，则要考虑空化边界。

在程序.rar的老程序里，边界处理通常不是写在方程里，而是通过“虚拟节点”实现的。例如对称面处布置一圈虚拟网格，令p(i,1)=p(i,0)。Python里更简单的做法是在迭代循环里固定边界行，不参与更新。固定值边界的实现很直接：

def apply_boundary_conditions(p, fixed_mask, fixed_values): # fixed_mask 为 True 的节点保持固定压力，不参与迭代更新 p[fixed_mask] = fixed_values[fixed_mask] return p

固定值边界处理中最容易犯的错误是把给定压力边界的节点也写进SOR的残差统计里，导致残差不下降。正确做法是残差只统计内部自由度，也就是mask为False的节点。

3. 用Python实现二维雷诺方程求解器：从参数表到压力分布

3.1 问题参数与无量纲化

以滑动轴承的简化模型为例，参数表如下：

直接代入国际单位制后，膜厚三次方在10的负13次方量级，压力在10的6次方量级，数值上很容易出现上溢或矩阵病态。常见做法是先做无量纲化：令x*=x/L、z*=z/B、H=h/h0，压力用P = p·h0²/(6ηUL)。无量纲化后的方程是：

∂/∂x*(H³ ∂P/∂x*) + (L/B)² ∂/∂z*(H³ ∂P/∂z*) = ∂H/∂x*

长宽比L/B以平方形式出现在交叉项里。对细长轴承，宽度方向扩散很弱，压力分布接近一维；对短轴承，z方向项占主导，这就是短轴承近似理论的基础。无量纲化后求解变量都在0~1之间，SOR迭代的收敛速度明显更快。

3.2 SOR迭代：老程序里的经典解法

虽然现在可以直接用稀疏矩阵求解器，但从二维雷诺方程的特性看，SOR迭代仍是程序.rar里最常见的老代码方案。原因有两个：一是系数矩阵虽稀疏但非对称，spsolve在网格数超过200×100时内存和耗时增长很快；二是SOR迭代天然支持边界条件的实时修正，特别是空化负压置零处理可以嵌在迭代里。

SOR迭代公式对每个内部节点是：

p_new = (1-ω)p_old + ω[(aE·pE + aW·pW + aN·pN + aS·pS - b) / aC]

ω是松弛因子，0<ω<2。ω=1就是Gauss-Seidel；ω>1是超松弛，通常取1.2~1.7。网格越密，ω需要取小一点，否则高频误差发散。

一个能直接跑的二维稳态求解器如下：

import numpy as np def solve_reynolds_2d(h, U, eta, L, B, nx=120, nz=60, omega=1.5, max_iter=50000, tol=1e-8): dx = L / (nx - 1) dz = B / (nz - 1) h3 = h ** 3 p = np.zeros((nx, nz)) rhs_factor = 6.0 * U * eta for it in range(max_iter): p_old = p.copy() max_res = 0.0 for i in range(1, nx - 1): for j in range(1, nz - 1): h3_E = 2 * h3[i, j] * h3[i, j+1] / (h3[i, j] + h3[i, j+1] + 1e-30) h3_W = 2 * h3[i, j] * h3[i, j-1] / (h3[i, j] + h3[i, j-1] + 1e-30) h3_N = 2 * h3[i, j] * h3[i-1, j] / (h3[i, j] + h3[i-1, j] + 1e-30) h3_S = 2 * h3[i, j] * h3[i+1, j] / (h3[i, j] + h3[i+1, j] + 1e-30) aE = h3_E / dx**2 aW = h3_W / dx**2 aN = h3_N / dz**2 aS = h3_S / dz**2 aC = aE + aW + aN + aS # 源项：中央差分 b = rhs_factor * (h[i, j+1] - h[i, j-1]) / (2.0 * dx) p_star = (aE * p[i, j+1] + aW * p[i, j-1] + aN * p[i-1, j] + aS * p[i+1, j] - b) / aC # SOR加权更新 p[i, j] = (1 - omega) * p[i, j] + omega * p_star # 负压置零：Gümbel空化条件 p = np.maximum(p, 0.0) # 残差统计（只统计内部节点） res = np.linalg.norm(p - p_old) / np.linalg.norm(p + 1e-30) if res < tol: break return p

这个求解器的更新逻辑是：逐个内部节点算出当前最优值，再与旧值做加权平均。空化处理直接放在每个迭代步之后，把所有负压截断为零，等效于最简单的空化模型。残差用的是相对变化量，避免绝对压力本身很大时误判。

注意事项：omega取值过大时，迭代初期残差反而上升，然后发散。如果程序.rar里的老代码用1.7，你的网格比它密2倍，不能直接照搬。一般先从1.2试跑，残差下降稳定后再慢慢加到1.5。h3在界面处用调和平均是为了防止薄膜厚区域的数值刚度。

3.3 压力分布后处理：最大压力、承载力与流量

求解完成后，第一步是看压力分布是否合理。对收敛的楔形滑块，压力峰值应当出现在收敛间隙中后部，而不是入口处；最大压力无量纲值一般在0.1~0.5之间，量纲值换算回去再与工程经验对比。

承载力是压力分布在轴承面积上的双重积分：

def capacity(p, dx, dz): # 双重梯形积分计算承载力 Wx = np.trapezoid(p, axis=0, dx=dx) W = np.trapezoid(Wx, dx=dz) return W

流量计算公式是每个截面上速度流和压力流的代数和。对x方向的体积流量：

Q_x = ∫(U·h/2 - h³/(12η)·∂p/∂x) dz

用中心差分近似∂p/∂x之后，可以在任意两个x位置算Q_x。稳态、不可压、无空化时，各截面流量应当守恒。我一般把流量检查放在承载力计算后面，它比压力极值更能暴露离散错误。如果两个截面的流量偏差超过5%，基本可以肯定是网格太粗或者膜厚不光滑。

4. 润滑压力分布求解的收敛判据、空化与网格陷阱

4.1 收敛判据：用相对残差而不是绝对残差

二维雷诺方程求解遇到“不收敛”时，反应不应该是盲目加大迭代次数。第一个要查的是判据本身。压力场的量级可能从几千到几百万帕，用绝对残差|p_new-p_old|<1e-3在某些单位下永远无法满足，在另一些单位下又过于宽松。标准化做法是统计内部节点上的相对变化量，前面代码里用的就是这种。

残差曲线还有一个参考价值：SOR迭代如果omega略大，残差曲线会先下降再反弹，这时候需要降低omega而不是加大迭代次数。如果曲线是平缓下降但极慢，多半是网格太密，omega要相应调小，或者改用多重网格/共轭梯度。老程序.rar里的固定迭代次数习惯，在现代机器上完全没有必要保留，建议改成残差收敛+最大迭代次数的双条件。

4.1.1 残差阈值怎么定

对大多数润滑分析，相对残差1e-6到1e-8已经足够。压力分布云图在1e-4就能看出形状，但承载力积分需要更高精度。注意：残差收敛不等于结果正确，只代表迭代没有发散。真实性要靠边界条件和物理校验。

4.2 Reynolds空化边界：把负压截断不是唯一方式

前面代码用的负压置零是Gümbel条件。严格来说，Reynolds空化条件要求破裂边界上同时满足p=0和∂p/∂x=0。这两种模型对压力分布的区别在出口区，承载力差异通常只有百分之几，但如果要算空化区域的流量，就必须用Reynolds条件找破裂位置。

程序.rar里很多老代码只在迭代中做p=max(p,0)，这是Gümbel，不是Reynolds。想升级成Reynolds条件，常见做法是先不做空化跑出p，找出负压区的边界，然后在边界上把∂p/∂x的差分改为单侧差分。工程上大多数轴承设计用Gümbel足够，要严谨推荐用软件里的空化模型。判断是否需要升级的简单办法：负压区面积占比超过20%，或者出口区有回流迹象时，再考虑改边界条件。

4.3 网格分辨率与膜厚梯度的离散误差

网格数选择直接影响压力峰值的精度。对楔形滑块这类膜厚连续变化的工况，50×25网格的峰值压力误差可能在5%左右，加密到200×100时误差通常小于1%。但要注意，膜厚如果来自有限元变形数据或实测粗糙度，插值后直接算∂h/∂x会出现锯齿，这种锯齿在压力场上表现为周期性振荡。

处理办法有两个：一是先对膜厚做平滑，例如Savitzky-Golay滤波；二是在膜厚突变处改用调和平均、在平缓区域维持算术平均。下面是网格收敛性检查的参考做法：

表中的偏差是相对上一级网格的差异。如果从100×50加密到200×100后，最大压力和承载力的变化仍在2%以上，说明还没进入网格收敛区间，再加密没有意义，要先回去检查离散格式或平滑膜厚。

需要注意对称面边界的离散：对称面∂p/∂z=0应该在边界节点用一阶单侧差分，直接复制相邻行也是一种近似。对于压力梯度较大的出口区，复制方式会在对称轴附近产生虚假的等值线扭曲。

5. 验证二维雷诺方程程序的三个实用技巧

5.1 用无限长滑块解析解做基准

一维稳态情况下，雷诺方程存在解析解。对固定斜平面滑块，dp/dx = 6ηU(h - h*)/h³，h*由进出口压力为零的条件确定。把二维程序的长宽比L/B设成0.1以下，宽度方向压力梯度自然趋近于零，算出来的压力分布就应该逼近这个一维解析解。验证时取最大压力的位置和峰值做对比，偏差在2%以内，程序基本没问题。

5.2 与老程序的输出数据对齐检查

把程序.rar里的输入数据导出为CSV格式，老程序算出的压力文件也转成同样网格的数组，用matplotlib把两条等值线叠在一起。如果峰值位置对不上，先查坐标方向定义：x是运动方向还是周向展开方向。如果只是峰值压力差很多，查单位，常见的坑是粘度用了cP但方程要求Pa·s。老代码注释里写的“黏度”并不一定是SI单位。

import matplotlib.pyplot as plt # old_p: 老程序读出的压力场，new_p: 当前Python结果 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) c1 = ax1.contourf(old_p, levels=20) ax1.set_title("legacy") c2 = ax2.contourf(new_p, levels=c1.levels) ax2.set_title("python") plt.show()

上面的levels用c1.levels而不是默认值，能保证两张图用同一套色标，肉眼比较时不会因为自动缩放掩盖差异。

5.3 流量守恒粗查比压力云图更敏感

承载力和压力峰值都可以通过调边界条件“凑”出来，但流量守恒是积分量，对离散误差更敏感。在x方向任意取两个截面，用前面给的流量公式分别积分，稳态工况下两者应该几乎相等。偏差超过3%，基本可以确定是网格太粗、膜厚锯齿或边界条件不完整。流量守恒检查可以写成一个独立函数，每次改完网格都跑一遍，这是最快暴露离散错误的方法。

本文还有配套的精品资源，点击获取