简介这份资源是黄平老师润滑数值计算方法的配套程序包面向机械工程、摩擦学方向的学习者与研究人员用于通过代码实践理解流体动力润滑理论、求解雷诺方程等非线性偏微分方程并预测油膜厚度、摩擦系数与功率损失。压缩包共62个文件约594KB以31个f90源文件为核心覆盖点接触、线接触、滑动轴承、挤压膜、气体润滑、热效应及脂润滑等多种工况另含pdb、plg、dsw、dsp等工程配置与调试文件以及exe可执行程序、bmp位图、doc使用说明和hm、rc等资源文件便于直接编译运行与对照学习。目前已有1509人学习下载。程序按工况分模块组织涉及有限差分、有限元等数值求解技术并处理速度、压力边界条件与温度对粘度的耦合影响读者可调整几何、材料、载荷、速度等参数观察结果变化从而培养解决实际润滑工程问题的能力适合作为学术研究与工程设计的参考工具。1. 润滑数值计算方法从雷诺方程到可运行程序的落地路径做滑动轴承或齿轮箱油膜分析时很多人卡在同一个地方手上有工况参数也翻过几本润滑理论的书但真要把膜厚、压力分布、承载力算出来中间那套数值流程始终串不起来。润滑数值计算方法要解决的正是这件事——把雷诺方程、膜厚方程、能量方程这些偏微分方程通过有限差分或有限元离散成线性方程组再用迭代法求解最终输出压力场、膜厚场和温升。它适合机械设计工程师、轴承方向的研究生以及需要做润滑仿真但不想从零推导公式的从业者。标题里的“附带程序”意味着这套方法最终要落到能跑通的代码上而不是停在公式推导。下面按“方程怎么离散、程序怎么搭、参数怎么调、坑在哪”的顺序拆开讲。2. 润滑数值计算方法的方程体系与离散策略2.1 雷诺方程在稳态与瞬态下的形式差异润滑问题的核心是雷诺方程它描述流体膜压力与膜厚、速度、粘度之间的关系。稳态工况下忽略时间项二维雷诺方程写成∂/∂x (h³/η · ∂p/∂x) ∂/∂z (h³/η · ∂p/∂z) 6U · ∂h/∂x其中 h 是膜厚η 是动力粘度p 是压力U 是相对滑动速度。瞬态工况下右边多出 ∂h/∂t 项动载轴承、启停过程必须用瞬态形式。很多初学者直接拿稳态程序算冲击载荷结果压力峰值偏低就是因为丢了时间项。离散时通常用有限差分法x 和 z 方向分别划分网格。对 ∂/∂x (h³/η · ∂p/∂x) 这一项采用中心差分[ \frac{1}{\Delta x^2}\left[ \left(\frac{h^3}{\eta}\right){i1/2,j}(p{i1,j}-p_{i,j}) - \left(\frac{h^3}{\eta}\right){i-1/2,j}(p{i,j}-p_{i-1,j}) \right] ]这里的关键是 h³/η 在网格界面上的取值。常见做法是取相邻节点算术平均但在膜厚变化剧烈区域比如油膜破裂边界算术平均会引入较大误差改用调和平均更稳。我一般会在膜厚梯度大的地方加密网格同时用调和平均处理界面系数。2.2 膜厚方程与能量方程的耦合方式膜厚方程本身不复杂刚性表面下 h h₀ 几何间隙。但弹性变形下h 还要叠加变形量 δ而 δ 又依赖压力分布 p这就形成耦合。常见做法是先假设刚性膜厚算一遍压力再用压力算变形更新膜厚反复迭代直到收敛。弹性变形计算可以用影响系数法把变形写成压力分布的卷积δ(x,z) ∬ K(x-x, z-z) · p(x,z) dxdzK 是影响系数核对半无限体问题有解析表达式。程序实现时把积分离散成求和影响系数矩阵提前算好存起来每次迭代只做矩阵向量乘。能量方程用于热弹流问题形式是ρcₚ(u ∂T/∂x w ∂T/∂z) η[(∂u/∂z)² (∂w/∂z)²] k(∂²T/∂x² ∂²T/∂z²)右边第一项是粘性耗散生热第二项是热传导。能量方程和雷诺方程通过粘度-温度关系耦合η 随 T 变化通常用 Barus 或 Roelands 公式。耦合迭代顺序一般是雷诺方程 → 能量方程 → 更新粘度 → 回到雷诺方程。这个循环收敛慢需要松弛因子。2.3 网格划分与边界条件的设定原则网格划分直接影响计算精度和收敛性。x 方向滑动方向通常均匀划分z 方向轴向在端部加密因为压力梯度大。网格数不是越多越好200×200 以上时迭代收敛速度明显下降而精度提升有限。我一般先用 80×80 试算确认收敛后再加密到 160×160 做最终结果。边界条件分三类入口和出口用压力边界通常取 p0轴向端部也是 p0油膜破裂边界用 Reynolds 条件即 p0 且 ∂p/∂n0。Reynolds 条件实现时常用“负压置零”的近似处理迭代过程中如果某节点压力算出负值直接置零同时把该节点从后续迭代中“冻结”。这个方法简单但可能高估承载区更严格的做法是用 Elrod 算法或质量守恒算法但实现复杂度高不少。提示油膜破裂边界处理是润滑计算里最容易翻车的地方负压置零虽然粗糙但对大多数稳态工况够用。如果发现承载力比实验值高 15% 以上优先检查这里。3. 润滑数值计算程序的搭建与核心模块实现3.1 用 Python 实现雷诺方程有限差分离散下面是一个稳态二维雷诺方程的离散求解核心代码用 Python 写成依赖 NumPy。这段代码只处理刚性膜厚、等温条件目的是把离散逻辑跑通。import numpy as np def solve_reynolds(nx, nz, Lx, Lz, h_grid, eta, U, p_initNone, max_iter5000, tol1e-6, omega1.5): 稳态二维雷诺方程有限差分解 nx, nz: x和z方向网格数 Lx, Lz: 计算域尺寸 h_grid: 膜厚矩阵 (nx, nz) eta: 动力粘度 (标量或矩阵) U: 滑动速度 omega: SOR松弛因子 dx Lx / (nx - 1) dz Lz / (nz - 1) p np.zeros((nx, nz)) if p_init is None else p_init.copy() # 预计算界面系数 h^3/eta用调和平均 h3_eta h_grid**3 / eta coef_x np.zeros((nx-1, nz)) coef_z np.zeros((nx, nz-1)) for i in range(nx-1): coef_x[i, :] 2 * h3_eta[i, :] * h3_eta[i1, :] / (h3_eta[i, :] h3_eta[i1, :] 1e-30) for j in range(nz-1): coef_z[:, j] 2 * h3_eta[:, j] * h3_eta[:, j1] / (h3_eta[:, j] h3_eta[:, j1] 1e-30) # 右边项 6U * dh/dx用中心差分 rhs np.zeros((nx, nz)) for i in range(1, nx-1): rhs[i, :] 6 * U * (h_grid[i1, :] - h_grid[i-1, :]) / (2 * dx) for it in range(max_iter): p_old p.copy() for i in range(1, nx-1): for j in range(1, nz-1): aE coef_x[i, j] / dx**2 aW coef_x[i-1, j] / dx**2 aN coef_z[i, j] / dz**2 aS coef_z[i, j-1] / dz**2 aP aE aW aN aS p[i, j] (1 - omega) * p[i, j] omega * (aE*p[i1,j] aW*p[i-1,j] aN*p[i,j1] aS*p[i,j-1] - rhs[i,j]) / aP if p[i, j] 0: p[i, j] 0.0 # 负压置零近似Reynolds条件 # 收敛判断 err np.max(np.abs(p - p_old)) / (np.max(np.abs(p)) 1e-30) if err tol: print(f收敛于第 {it} 次迭代残差 {err:.2e}) break else: print(达到最大迭代次数未完全收敛) return p这段代码的逻辑说明先算界面系数 coef_x 和 coef_z用调和平均避免膜厚突变处的数值震荡右边项 rhs 是 6U·∂h/∂x 的离散迭代用 SOR逐次超松弛omega 取 1.5 左右通常收敛最快超过 1.8 容易发散。负压置零放在每次节点更新之后这是最简单的 Reynolds 条件近似。参数说明nx 和 nz 建议从 61 起步Lx 和 Lz 按实际轴承尺寸给h_grid 的最小值不能为零否则 h³/η 会出问题一般设最小膜厚 1e-6 m。eta 如果是等温计算给标量即可热弹流问题需要传矩阵。U 的单位是 m/s注意和 h、eta 的单位统一。3.2 弹性变形耦合与迭代收敛控制刚性解算完后如果要考虑弹性变形需要加一个变形计算模块。下面用影响系数法做简化实现def compute_deformation(p, K_matrix): 用影响系数矩阵计算弹性变形 p: 压力矩阵展平为一维 K_matrix: 影响系数矩阵形状 (n_nodes, n_nodes) p_flat p.flatten() delta_flat K_matrix p_flat return delta_flat.reshape(p.shape) def elastic_coupling(nx, nz, Lx, Lz, h0, eta, U, K_matrix, max_outer50, relax0.3): 弹性耦合迭代雷诺方程 变形更新膜厚 relax: 膜厚更新松弛因子防止振荡 h h0.copy() for outer in range(max_outer): p solve_reynolds(nx, nz, Lx, Lz, h, eta, U, max_iter2000, tol1e-5) delta compute_deformation(p, K_matrix) h_new h0 delta h (1 - relax) * h relax * h_new # 松弛更新 if outer % 5 0: print(f外循环 {outer}最大压力 {np.max(p):.3e} Pa最小膜厚 {np.min(h):.6e} m) return p, h逻辑说明外层循环先解雷诺方程得到压力再用影响系数矩阵算变形更新膜厚后回到雷诺方程。松弛因子 relax 取 0.3 左右太大容易振荡太小收敛慢。影响系数矩阵 K_matrix 需要提前根据几何和材料参数算好对半无限体问题K 的解析形式是 2/(πE) · 1/rE 是等效弹性模量。参数说明max_outer 一般 30 到 50 次足够如果 50 次还没稳定多半是松弛因子或网格有问题。判断收敛看最小膜厚和最大压力的变化连续两次外循环变化小于 1% 即可停止。3.3 结果后处理压力场与膜厚场的可视化验证算完之后必须做可视化检查否则很容易把错误结果当成收敛解。用 Matplotlib 画压力场和膜厚场import matplotlib.pyplot as plt def plot_results(p, h, Lx, Lz): fig, axes plt.subplots(1, 2, figsize(12, 5)) extent [0, Lx*1000, 0, Lz*1000] # 转成毫米 im1 axes[0].imshow(p.T, originlower, extentextent, aspectauto, cmapjet) axes[0].set_xlabel(x (mm)) axes[0].set_ylabel(z (mm)) axes[0].set_title(压力分布 (Pa)) plt.colorbar(im1, axaxes[0]) im2 axes[1].imshow(h.T*1e6, originlower, extentextent, aspectauto, cmapviridis) axes[1].set_xlabel(x (mm)) axes[1].set_ylabel(z (mm)) axes[1].set_title(膜厚分布 (μm)) plt.colorbar(im2, axaxes[1]) plt.tight_layout() plt.savefig(lubrication_results.png, dpi150) plt.show()检查要点压力场应该是单峰或双峰取决于膜厚形状峰值位置在最小膜厚稍下游膜厚场在承载区应该明显小于入口区压力场边缘不应该有剧烈震荡如果有说明网格不够密或界面系数处理有问题。我一般还会沿 x 方向中间切一条线看压力曲线是否光滑以及是否满足出口处 p0 且 ∂p/∂x0。4. 润滑数值计算避坑5 个让程序跑不出正确结果的典型问题4.1 迭代发散压力值指数增长到无穷现象迭代几次后压力值迅速增大几十步后变成 NaN 或 1e30 量级。原因松弛因子 omega 过大或者界面系数计算出负值。膜厚网格如果出现零或负值h³/η 会变成零或负数导致系数矩阵失去对角占优。解决先把 omega 降到 1.0 甚至 0.5 试算确认能收敛后再逐步提高。检查 h_grid 最小值确保所有节点膜厚大于 1e-8 m。如果膜厚确实接近零比如边界润滑需要改用更稳定的算法比如把 h³/η 在极小膜厚处做截断处理。4.2 收敛但结果明显偏大承载力比理论值高 30% 以上现象程序显示收敛但积分压力得到的承载力远高于 Sommerfeld 数或短轴承理论的估算值。原因负压置零后没有正确处理油膜破裂边界导致承载区被人为扩大。另一个常见原因是网格太粗压力峰值被抹平后积分面积反而偏大。解决加密网格到 160×160 以上再看。如果承载力仍然偏高检查负压置零逻辑是否在每次迭代都执行以及是否冻结了破裂区节点。更严格的做法是实现质量守恒算法但工作量较大。可以先和已知解析解对比比如无限长轴承的 Sommerfeld 解确认程序在简单工况下正确。4.3 弹性耦合振荡膜厚和压力来回跳不收敛现象外循环中最小膜厚和最大压力反复变化没有稳定趋势。原因松弛因子 relax 太大或者影响系数矩阵计算有误。影响系数矩阵如果对称性不对变形计算会引入非物理的振荡。解决把 relax 降到 0.1 到 0.2增加外循环次数。检查影响系数矩阵是否对称对半无限体问题K(i,j) 应该等于 K(j,i)。如果几何复杂影响系数需要用有限元算这时要确保有限元网格和润滑网格的插值正确。4.4 温度场计算不收敛粘度更新后压力场剧烈变化现象热弹流计算中能量方程解出的温度场代入粘度公式后压力场变化超过 50%下一轮温度又反向变化。原因粘温方程系数过大或者能量方程的对流项离散不稳定。温度对粘度的影响是指数级的温度稍微变化几度粘度可能差一倍。解决在粘温更新时加欠松弛比如 η_new 0.5·η_old 0.5·η_calculated。能量方程的对流项用迎风格式不要用中心差分。如果还是振荡先把粘温效应关掉算等温解再逐步打开粘温耦合每次只更新部分节点的粘度。4.5 单位不统一导致结果差几个数量级现象程序跑通但压力是 1e12 Pa 或者膜厚是 1e-2 m明显不在合理范围。原因长度用 mm、粘度用 Pa·s、速度用 m/s混在一起算。雷诺方程里 h³/η 的量纲是 m³/(Pa·s)如果 h 用 mm结果会差 10⁹ 倍。解决所有输入统一到 SI 单位长度 m压力 Pa粘度 Pa·s速度 m/s。膜厚通常给 μm记得乘 1e-6。粘度如果给的是 cP厘泊乘 1e-3 转成 Pa·s。养成在函数入口做单位检查的习惯比如断言 h_grid.min() 1e-8 且 h_grid.max() 1e-2。5. 润滑数值计算方法的进阶技巧从能跑到可信程序能跑通只是第一步要让结果可信还需要做验证和确认。我一般用三个手段第一和解析解对比。无限长轴承、无限短轴承都有解析解把程序设成对应工况看误差是否在 5% 以内。第二做网格无关性验证。同一工况用 80×80、160×160、320×320 各算一遍如果 160 和 320 的结果差小于 2%说明 160 够用。第三和实验数据对比。如果手头有轴承试验台的膜厚或摩擦力数据直接对比是最有说服力的。进阶方向有两个一是瞬态计算把时间项加进去用 Crank-Nicolson 或隐式 Euler 离散可以算启停、冲击载荷下的润滑响应。二是表面织构或粗糙度效应把膜厚方程里加入织构几何或随机粗糙度研究微造型对承载力的影响。这两个方向都需要在现有程序上改但核心离散逻辑不变。注意瞬态计算的时间步长不能太大一般取特征时间比如轴转一圈的时间的 1/100 到 1/1000。步长太大时间项离散误差会淹没物理结果。我自己的习惯是每改一次代码先用一个已知答案的算例跑一遍回归测试确认没改坏再上复杂工况。润滑计算里一个符号错误或者单位错误就能让结果完全失去意义后悔药没得吃。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/28 7:22:47 · 阅读 1.2k 润滑数值计算方法：从雷诺方程到可运行程序
简介这份资源是黄平老师润滑数值计算方法的配套程序包面向机械工程、摩擦学方向的学习者与研究人员用于通过代码实践理解流体动力润滑理论、求解雷诺方程等非线性偏微分方程并预测油膜厚度、摩擦系数与功率损失。压缩包共62个文件约594KB以31个f90源文件为核心覆盖点接触、线接触、滑动轴承、挤压膜、气体润滑、热效应及脂润滑等多种工况另含pdb、plg、dsw、dsp等工程配置与调试文件以及exe可执行程序、bmp位图、doc使用说明和hm、rc等资源文件便于直接编译运行与对照学习。目前已有1509人学习下载。程序按工况分模块组织涉及有限差分、有限元等数值求解技术并处理速度、压力边界条件与温度对粘度的耦合影响读者可调整几何、材料、载荷、速度等参数观察结果变化从而培养解决实际润滑工程问题的能力适合作为学术研究与工程设计的参考工具。1. 润滑数值计算方法从雷诺方程到可运行程序的落地路径做滑动轴承或齿轮箱油膜分析时很多人卡在同一个地方手上有工况参数也翻过几本润滑理论的书但真要把膜厚、压力分布、承载力算出来中间那套数值流程始终串不起来。润滑数值计算方法要解决的正是这件事——把雷诺方程、膜厚方程、能量方程这些偏微分方程通过有限差分或有限元离散成线性方程组再用迭代法求解最终输出压力场、膜厚场和温升。它适合机械设计工程师、轴承方向的研究生以及需要做润滑仿真但不想从零推导公式的从业者。标题里的“附带程序”意味着这套方法最终要落到能跑通的代码上而不是停在公式推导。下面按“方程怎么离散、程序怎么搭、参数怎么调、坑在哪”的顺序拆开讲。2. 润滑数值计算方法的方程体系与离散策略2.1 雷诺方程在稳态与瞬态下的形式差异润滑问题的核心是雷诺方程它描述流体膜压力与膜厚、速度、粘度之间的关系。稳态工况下忽略时间项二维雷诺方程写成∂/∂x (h³/η · ∂p/∂x) ∂/∂z (h³/η · ∂p/∂z) 6U · ∂h/∂x其中 h 是膜厚η 是动力粘度p 是压力U 是相对滑动速度。瞬态工况下右边多出 ∂h/∂t 项动载轴承、启停过程必须用瞬态形式。很多初学者直接拿稳态程序算冲击载荷结果压力峰值偏低就是因为丢了时间项。离散时通常用有限差分法x 和 z 方向分别划分网格。对 ∂/∂x (h³/η · ∂p/∂x) 这一项采用中心差分[ \frac{1}{\Delta x^2}\left[ \left(\frac{h^3}{\eta}\right){i1/2,j}(p{i1,j}-p_{i,j}) - \left(\frac{h^3}{\eta}\right){i-1/2,j}(p{i,j}-p_{i-1,j}) \right] ]这里的关键是 h³/η 在网格界面上的取值。常见做法是取相邻节点算术平均但在膜厚变化剧烈区域比如油膜破裂边界算术平均会引入较大误差改用调和平均更稳。我一般会在膜厚梯度大的地方加密网格同时用调和平均处理界面系数。2.2 膜厚方程与能量方程的耦合方式膜厚方程本身不复杂刚性表面下 h h₀ 几何间隙。但弹性变形下h 还要叠加变形量 δ而 δ 又依赖压力分布 p这就形成耦合。常见做法是先假设刚性膜厚算一遍压力再用压力算变形更新膜厚反复迭代直到收敛。弹性变形计算可以用影响系数法把变形写成压力分布的卷积δ(x,z) ∬ K(x-x, z-z) · p(x,z) dxdzK 是影响系数核对半无限体问题有解析表达式。程序实现时把积分离散成求和影响系数矩阵提前算好存起来每次迭代只做矩阵向量乘。能量方程用于热弹流问题形式是ρcₚ(u ∂T/∂x w ∂T/∂z) η[(∂u/∂z)² (∂w/∂z)²] k(∂²T/∂x² ∂²T/∂z²)右边第一项是粘性耗散生热第二项是热传导。能量方程和雷诺方程通过粘度-温度关系耦合η 随 T 变化通常用 Barus 或 Roelands 公式。耦合迭代顺序一般是雷诺方程 → 能量方程 → 更新粘度 → 回到雷诺方程。这个循环收敛慢需要松弛因子。2.3 网格划分与边界条件的设定原则网格划分直接影响计算精度和收敛性。x 方向滑动方向通常均匀划分z 方向轴向在端部加密因为压力梯度大。网格数不是越多越好200×200 以上时迭代收敛速度明显下降而精度提升有限。我一般先用 80×80 试算确认收敛后再加密到 160×160 做最终结果。边界条件分三类入口和出口用压力边界通常取 p0轴向端部也是 p0油膜破裂边界用 Reynolds 条件即 p0 且 ∂p/∂n0。Reynolds 条件实现时常用“负压置零”的近似处理迭代过程中如果某节点压力算出负值直接置零同时把该节点从后续迭代中“冻结”。这个方法简单但可能高估承载区更严格的做法是用 Elrod 算法或质量守恒算法但实现复杂度高不少。提示油膜破裂边界处理是润滑计算里最容易翻车的地方负压置零虽然粗糙但对大多数稳态工况够用。如果发现承载力比实验值高 15% 以上优先检查这里。3. 润滑数值计算程序的搭建与核心模块实现3.1 用 Python 实现雷诺方程有限差分离散下面是一个稳态二维雷诺方程的离散求解核心代码用 Python 写成依赖 NumPy。这段代码只处理刚性膜厚、等温条件目的是把离散逻辑跑通。import numpy as np def solve_reynolds(nx, nz, Lx, Lz, h_grid, eta, U, p_initNone, max_iter5000, tol1e-6, omega1.5): 稳态二维雷诺方程有限差分解 nx, nz: x和z方向网格数 Lx, Lz: 计算域尺寸 h_grid: 膜厚矩阵 (nx, nz) eta: 动力粘度 (标量或矩阵) U: 滑动速度 omega: SOR松弛因子 dx Lx / (nx - 1) dz Lz / (nz - 1) p np.zeros((nx, nz)) if p_init is None else p_init.copy() # 预计算界面系数 h^3/eta用调和平均 h3_eta h_grid**3 / eta coef_x np.zeros((nx-1, nz)) coef_z np.zeros((nx, nz-1)) for i in range(nx-1): coef_x[i, :] 2 * h3_eta[i, :] * h3_eta[i1, :] / (h3_eta[i, :] h3_eta[i1, :] 1e-30) for j in range(nz-1): coef_z[:, j] 2 * h3_eta[:, j] * h3_eta[:, j1] / (h3_eta[:, j] h3_eta[:, j1] 1e-30) # 右边项 6U * dh/dx用中心差分 rhs np.zeros((nx, nz)) for i in range(1, nx-1): rhs[i, :] 6 * U * (h_grid[i1, :] - h_grid[i-1, :]) / (2 * dx) for it in range(max_iter): p_old p.copy() for i in range(1, nx-1): for j in range(1, nz-1): aE coef_x[i, j] / dx**2 aW coef_x[i-1, j] / dx**2 aN coef_z[i, j] / dz**2 aS coef_z[i, j-1] / dz**2 aP aE aW aN aS p[i, j] (1 - omega) * p[i, j] omega * (aE*p[i1,j] aW*p[i-1,j] aN*p[i,j1] aS*p[i,j-1] - rhs[i,j]) / aP if p[i, j] 0: p[i, j] 0.0 # 负压置零近似Reynolds条件 # 收敛判断 err np.max(np.abs(p - p_old)) / (np.max(np.abs(p)) 1e-30) if err tol: print(f收敛于第 {it} 次迭代残差 {err:.2e}) break else: print(达到最大迭代次数未完全收敛) return p这段代码的逻辑说明先算界面系数 coef_x 和 coef_z用调和平均避免膜厚突变处的数值震荡右边项 rhs 是 6U·∂h/∂x 的离散迭代用 SOR逐次超松弛omega 取 1.5 左右通常收敛最快超过 1.8 容易发散。负压置零放在每次节点更新之后这是最简单的 Reynolds 条件近似。参数说明nx 和 nz 建议从 61 起步Lx 和 Lz 按实际轴承尺寸给h_grid 的最小值不能为零否则 h³/η 会出问题一般设最小膜厚 1e-6 m。eta 如果是等温计算给标量即可热弹流问题需要传矩阵。U 的单位是 m/s注意和 h、eta 的单位统一。3.2 弹性变形耦合与迭代收敛控制刚性解算完后如果要考虑弹性变形需要加一个变形计算模块。下面用影响系数法做简化实现def compute_deformation(p, K_matrix): 用影响系数矩阵计算弹性变形 p: 压力矩阵展平为一维 K_matrix: 影响系数矩阵形状 (n_nodes, n_nodes) p_flat p.flatten() delta_flat K_matrix p_flat return delta_flat.reshape(p.shape) def elastic_coupling(nx, nz, Lx, Lz, h0, eta, U, K_matrix, max_outer50, relax0.3): 弹性耦合迭代雷诺方程 变形更新膜厚 relax: 膜厚更新松弛因子防止振荡 h h0.copy() for outer in range(max_outer): p solve_reynolds(nx, nz, Lx, Lz, h, eta, U, max_iter2000, tol1e-5) delta compute_deformation(p, K_matrix) h_new h0 delta h (1 - relax) * h relax * h_new # 松弛更新 if outer % 5 0: print(f外循环 {outer}最大压力 {np.max(p):.3e} Pa最小膜厚 {np.min(h):.6e} m) return p, h逻辑说明外层循环先解雷诺方程得到压力再用影响系数矩阵算变形更新膜厚后回到雷诺方程。松弛因子 relax 取 0.3 左右太大容易振荡太小收敛慢。影响系数矩阵 K_matrix 需要提前根据几何和材料参数算好对半无限体问题K 的解析形式是 2/(πE) · 1/rE 是等效弹性模量。参数说明max_outer 一般 30 到 50 次足够如果 50 次还没稳定多半是松弛因子或网格有问题。判断收敛看最小膜厚和最大压力的变化连续两次外循环变化小于 1% 即可停止。3.3 结果后处理压力场与膜厚场的可视化验证算完之后必须做可视化检查否则很容易把错误结果当成收敛解。用 Matplotlib 画压力场和膜厚场import matplotlib.pyplot as plt def plot_results(p, h, Lx, Lz): fig, axes plt.subplots(1, 2, figsize(12, 5)) extent [0, Lx*1000, 0, Lz*1000] # 转成毫米 im1 axes[0].imshow(p.T, originlower, extentextent, aspectauto, cmapjet) axes[0].set_xlabel(x (mm)) axes[0].set_ylabel(z (mm)) axes[0].set_title(压力分布 (Pa)) plt.colorbar(im1, axaxes[0]) im2 axes[1].imshow(h.T*1e6, originlower, extentextent, aspectauto, cmapviridis) axes[1].set_xlabel(x (mm)) axes[1].set_ylabel(z (mm)) axes[1].set_title(膜厚分布 (μm)) plt.colorbar(im2, axaxes[1]) plt.tight_layout() plt.savefig(lubrication_results.png, dpi150) plt.show()检查要点压力场应该是单峰或双峰取决于膜厚形状峰值位置在最小膜厚稍下游膜厚场在承载区应该明显小于入口区压力场边缘不应该有剧烈震荡如果有说明网格不够密或界面系数处理有问题。我一般还会沿 x 方向中间切一条线看压力曲线是否光滑以及是否满足出口处 p0 且 ∂p/∂x0。4. 润滑数值计算避坑5 个让程序跑不出正确结果的典型问题4.1 迭代发散压力值指数增长到无穷现象迭代几次后压力值迅速增大几十步后变成 NaN 或 1e30 量级。原因松弛因子 omega 过大或者界面系数计算出负值。膜厚网格如果出现零或负值h³/η 会变成零或负数导致系数矩阵失去对角占优。解决先把 omega 降到 1.0 甚至 0.5 试算确认能收敛后再逐步提高。检查 h_grid 最小值确保所有节点膜厚大于 1e-8 m。如果膜厚确实接近零比如边界润滑需要改用更稳定的算法比如把 h³/η 在极小膜厚处做截断处理。4.2 收敛但结果明显偏大承载力比理论值高 30% 以上现象程序显示收敛但积分压力得到的承载力远高于 Sommerfeld 数或短轴承理论的估算值。原因负压置零后没有正确处理油膜破裂边界导致承载区被人为扩大。另一个常见原因是网格太粗压力峰值被抹平后积分面积反而偏大。解决加密网格到 160×160 以上再看。如果承载力仍然偏高检查负压置零逻辑是否在每次迭代都执行以及是否冻结了破裂区节点。更严格的做法是实现质量守恒算法但工作量较大。可以先和已知解析解对比比如无限长轴承的 Sommerfeld 解确认程序在简单工况下正确。4.3 弹性耦合振荡膜厚和压力来回跳不收敛现象外循环中最小膜厚和最大压力反复变化没有稳定趋势。原因松弛因子 relax 太大或者影响系数矩阵计算有误。影响系数矩阵如果对称性不对变形计算会引入非物理的振荡。解决把 relax 降到 0.1 到 0.2增加外循环次数。检查影响系数矩阵是否对称对半无限体问题K(i,j) 应该等于 K(j,i)。如果几何复杂影响系数需要用有限元算这时要确保有限元网格和润滑网格的插值正确。4.4 温度场计算不收敛粘度更新后压力场剧烈变化现象热弹流计算中能量方程解出的温度场代入粘度公式后压力场变化超过 50%下一轮温度又反向变化。原因粘温方程系数过大或者能量方程的对流项离散不稳定。温度对粘度的影响是指数级的温度稍微变化几度粘度可能差一倍。解决在粘温更新时加欠松弛比如 η_new 0.5·η_old 0.5·η_calculated。能量方程的对流项用迎风格式不要用中心差分。如果还是振荡先把粘温效应关掉算等温解再逐步打开粘温耦合每次只更新部分节点的粘度。4.5 单位不统一导致结果差几个数量级现象程序跑通但压力是 1e12 Pa 或者膜厚是 1e-2 m明显不在合理范围。原因长度用 mm、粘度用 Pa·s、速度用 m/s混在一起算。雷诺方程里 h³/η 的量纲是 m³/(Pa·s)如果 h 用 mm结果会差 10⁹ 倍。解决所有输入统一到 SI 单位长度 m压力 Pa粘度 Pa·s速度 m/s。膜厚通常给 μm记得乘 1e-6。粘度如果给的是 cP厘泊乘 1e-3 转成 Pa·s。养成在函数入口做单位检查的习惯比如断言 h_grid.min() 1e-8 且 h_grid.max() 1e-2。5. 润滑数值计算方法的进阶技巧从能跑到可信程序能跑通只是第一步要让结果可信还需要做验证和确认。我一般用三个手段第一和解析解对比。无限长轴承、无限短轴承都有解析解把程序设成对应工况看误差是否在 5% 以内。第二做网格无关性验证。同一工况用 80×80、160×160、320×320 各算一遍如果 160 和 320 的结果差小于 2%说明 160 够用。第三和实验数据对比。如果手头有轴承试验台的膜厚或摩擦力数据直接对比是最有说服力的。进阶方向有两个一是瞬态计算把时间项加进去用 Crank-Nicolson 或隐式 Euler 离散可以算启停、冲击载荷下的润滑响应。二是表面织构或粗糙度效应把膜厚方程里加入织构几何或随机粗糙度研究微造型对承载力的影响。这两个方向都需要在现有程序上改但核心离散逻辑不变。注意瞬态计算的时间步长不能太大一般取特征时间比如轴转一圈的时间的 1/100 到 1/1000。步长太大时间项离散误差会淹没物理结果。我自己的习惯是每改一次代码先用一个已知答案的算例跑一遍回归测试确认没改坏再上复杂工况。润滑计算里一个符号错误或者单位错误就能让结果完全失去意义后悔药没得吃。希望帮到你。本文还有配套的精品资源点击获取