简介：本资源是一套面向机械工程与流体润滑领域初学者及进阶研究者的MATLAB数值计算实践方案，聚焦于气体静压轴承性能分析这一典型工程问题，通过有限差分法高效求解非线性雷诺方程，获得压力分布、承载力、刚度等关键特性参数。压缩包共含2个文件（18KB），包括核心求解脚本（.m）与配套技术说明文档（.docx），前者实现网格离散、边界条件处理、迭代收敛控制及结果可视化，后者详述物理模型、差分格式推导与参数设置依据，便于理解算法原理与工程适配逻辑。已有1713人学习下载，源码经实测校正，可直接运行并支持参数修改拓展，适合开展课程设计、科研建模或仿真方法入门训练。 做气体静压轴承的都知道，一上手最头疼的不是选型，而是怎么把气膜里的压力分布算出来。你查文献能看到一堆雷诺方程的变体，真到自己写代码的时候，离散格式怎么选、供气孔怎么处理、迭代怎么收敛，全是坑。这篇文章把我在MATLAB里用有限差分法求解气体静压轴承雷诺方程的过程完整过一遍，从模型简化到源码拆解，再到承载力、耗气量、刚度的计算，全部是直接能跑起来的代码。适合刚开始接触气体润滑、做气浮导轨气浮平台的工程师和研究生，也适合想快速算一个初步方案的工程人员。
1. 气体静压轴承求解思路与方案选型
1.1 气体静压轴承的润滑本质
气体静压轴承的工作原理其实不复杂：外部气源把高压气体输送到轴承面，经过供气孔（或多孔质节流器）进入气膜间隙，在间隙里形成高压气膜，靠静压把承载面托起来。和液体滑动轴承相比，气体轴承的黏度极低，摩擦力小、无污染、精度高，所以精密测量、超精密加工、气浮平台基本都离不开它。
不过“极低”这俩字是双刃剑。黏度低意味着承载能力远不如油润滑，设计的时候必须把气膜压力分布算清楚，稍不留神承载力不够，整个平台就压趴了。计算气膜压力分布的基本方程，就是气体润滑雷诺方程。
雷诺方程本质上是从Navier-Stokes方程在“薄气膜”条件下简化出来的，它假设膜厚远小于轴承平面尺寸、流动是层流、气体为等温理想气体。这几个假设对常规工况下的静压气浮轴承来说基本成立，所以用雷诺方程做工程设计是经典路线。
1.2 为什么选有限差分法而不是FEM或CFD
拿到雷诺方程，求压力分布有几种路线：解析解、有限差分、有限元、商业CFD。
解析解只存在于无限长轴承或极简单几何里，实际矩形轴承基本拿不到。CFD能算得很精细，但网格量大、计算慢，而且对一般设计需求有点杀鸡用牛刀。有限元适合复杂几何，但要自己写网格剖分，上手成本高。
有限差分法最直接：把矩形轴承区域切成规规矩矩的网格，用差分近似替代微分，得到一个代数方程组，迭代求解。规则网格下实现非常快，MATLAB里几十行代码就能跑通，对矩形推力轴承、节流孔均布的静压轴承特别合适。初版方案用它快速迭代，确认趋势后再考虑更重的工具，是性价比最高的路线。
2. 控制方程离散与无量纲化处理
2.1 雷诺方程简化与适用条件
没有相对剪切运动时（纯静压工况），二维稳态可压缩雷诺方程可以写成：
∂/∂x (p h³ ∂p/∂x) + ∂/∂y (p h³ ∂p/∂y) = 0
这里 p 是绝对压力，h 是气膜厚度。方程里为什么会有p乘在 h³ 里面？因为气体可压缩，密度随压力变化，在等温假设下 ρ ∝ p，所以流量项里会多出一个压力因子。这也是气体润滑和液体润滑最大的区别，不能直接把液体雷诺方程拿来用。
对矩形推力轴承，如果膜厚均匀（h = h0），方程表面上只剩一个变量p，但因为p自己出现在系数 ph³ 里，方程是非线性的，没法一步求解，必须迭代。
我见过不少人一开始直接把液体润滑的差分格式套到气体问题上，结果算出来的压力分布要么不收敛，要么出现负压，根因就是漏掉了这个可压缩项。
2.2 无量纲化的目的与具体步骤
数值计算里，压力量级是10⁵ Pa，膜厚是10⁻⁵ m，如果用有量纲量直接算，系数跨度十个数量级，迭代误差和浮点截断都容易出问题。最常用的办法是无量纲化。
令 X = x/L、Y = y/L（矩形轴承用同一边长做参考）、P = p/pa、H = h/h0，方程可以写成：
∂/∂X (P H³ ∂P/∂X) + (L/B)² ∂/∂Y (P H³ ∂P/∂Y) = 0
其中L、B分别是轴承长宽。无量纲化的好处：
- 网格坐标归一化到0~1，程序里不用每次乘物理尺寸
- 压力归一化到1的量级，迭代矩阵条件数更友好
- 结果可以直接对比文献里的无量纲承载力系数
实际编程时，我习惯先无量纲化，算完以后再反算回物理量。这条经验是后来对比了好几篇论文才养成的习惯，好处是收敛速度明显更快，调试代码时量级也更直观。
2.3 边界条件与供气孔建模
边界条件很直接：轴承四周暴露在大气里，压力等于环境压力 pa；供气孔处压力等于节流后压力 pd。
严格说，pd 需要联立节流器流量方程才能定，入门版本直接给固定压力，或者按供气压力的一定比例估算。我初版代码里直接设成固定压力 ps，这样程序最简单，也能反映主要物理趋势。
很多第一次接触这个问题的同学会问：供气孔那么多，为什么不是每个孔周围都围一圈高压区？因为气体是可压缩的，会从孔沿径向往四周扩散，孔之间压力会相互叠加。供气孔的位置、数量、压力直接决定压力分布形态，这是后面参数分析的重点。
3. Matlab源码实现与关键代码解析
3.1 网格生成与供气孔定位
先给出完整的主程序，我用的是100mm×100mm方形推力轴承，四角加中心共五个供气孔，孔口绝对压力0.6MPa，设计膜厚20μm：
% 气体静压矩形推力轴承 —— 有限差分法求解雷诺方程 clear; clc; close all; %% 1. 参数定义 Lx = 0.1; Ly = 0.1; % 轴承长宽 [m] h0 = 20e-6; % 设计气膜厚度 [m] pa = 101325; % 环境压力 [Pa] ps = 600000; % 供气孔绝对压力 [Pa] mu = 1.82e-5; % 空气动力黏度 [Pa.s] R = 287.1; T = 293; % 气体常数与温度 % 供气孔坐标 [x, y]，单位 m orifice = [ 0.025, 0.025; 0.025, 0.075; 0.075, 0.025; 0.075, 0.075; 0.050, 0.050 ]; %% 2. 网格划分 nx = 101; ny = 101; x = linspace(0, Lx, nx); y = linspace(0, Ly, ny); dx = x(2) - x(1); dy = y(2) - y(1); [X, Y] = meshgrid(x, y); %% 3. 初始化 h = h0 * ones(ny, nx); % 膜厚场 p = pa * ones(ny, nx); % 压力场初始化为环境压力 % 供气孔掩膜 orifice_mask = false(ny, nx); for k = 1:size(orifice, 1) [~, ix] = min(abs(x - orifice(k, 1))); [~, iy] = min(abs(y - orifice(k, 2))); orifice_mask(iy, ix) = true; end %% 4. SOR 迭代求解 omega = 1.6; % 松弛因子 tol = 1e-6; % 收敛判据 maxIter = 20000; for iter = 1:maxIter p_old = p; for j = 2:ny-1 for i = 2:nx-1 if orifice_mask(j, i) continue; end % 界面 ph^3 线性化（算术平均） pE = p(j, i+1); pW = p(j, i-1); pN = p(j+1, i); pS = p(j-1, i); F_E = (pE + p(j,i)) / 2 * ((h(j,i+1)^3 + h(j,i)^3) / 2); F_W = (pW + p(j,i)) / 2 * ((h(j,i-1)^3 + h(j,i)^3) / 2); F_N = (pN + p(j,i)) / 2 * ((h(j+1,i)^3 + h(j,i)^3) / 2); F_S = (pS + p(j,i)) / 2 * ((h(j-1,i)^3 + h(j,i)^3) / 2); aE = F_E / dx^2; aW = F_W / dx^2; aN = F_N / dy^2; aS = F_S / dy^2; aC = aE + aW + aN + aS; p_star = (aE*pE + aW*pW + aN*pN + aS*pS) / aC; p(j,i) = omega * p_star + (1 - omega) * p(j,i); end end % 边界条件与供气孔覆盖（顺序不能反） p(1,:) = pa; p(ny,:) = pa; p(:,1) = pa; p(:,nx) = pa; p(orifice_mask) = ps; err = max(abs(p - p_old), [], 'all'); if mod(iter, 200) == 0 fprintf('iter %4d, err = %.3e\n', iter, err); end if err < tol fprintf('converged at iter %4d, err = %.3e\n', iter, err); break; end end %% 5. 后处理 figure('Color','w'); surf(X*1000, Y*1000, p/1000, 'EdgeColor', 'none'); xlabel('x / mm'); ylabel('y / mm'); zlabel('p / kPa'); title('气膜压力分布'); %% 6. 特性计算 W = sum(p - pa, 'all') * dx * dy; % 承载力 N rho_a = pa / (R * T); q_top = sum(h(1,:).^3 / (12*mu) * (p(1,:) - p(2,:)) / dy * rho_a) * dx; q_bottom = sum(h(ny,:).^3 / (12*mu) * (p(ny,:) - p(ny-1,:)) / dy * rho_a) * dx; q_left = sum(h(:,1).^3 / (12*mu) * (p(:,1) - p(:,2)) / dx * rho_a) * dy; q_right = sum(h(:,nx).^3 / (12*mu) * (p(:,nx) - p(:,nx-1)) / dx * rho_a) * dy; Q = abs(q_top) + abs(q_bottom) + abs(q_left) + abs(q_right); % 总耗气量 kg/s fprintf('承载力 W = %.1f N\n', W); fprintf('平均比压 = %.2f kPa\n', W/(Lx*Ly)/1000); fprintf('总耗气量 Q = %.4e kg/s\n', Q);网格生成和供气孔定位这里有一个细节：linspace生成等距坐标，meshgrid生成网格坐标矩阵，而供气孔我用一个二维逻辑掩膜orifice_mask来标记。把距离供气孔最近的网格节点找出来，迭代时这些节点不参与差分更新，始终固定为 ps。这个“掩膜法”比在每个迭代步里用if逐点判断坐标快得多，代码也更清晰，后面如果要改孔位布局，只要改orifice矩阵就行了。
3.2 差分格式与SOR迭代核心代码
主循环是整段代码的核心，有三个点必须讲透：
第一，界面系数 ph³ 的取值。差分格式里需要界面处的 F=ph³，最简单是用两侧节点值的算术平均，即 F_{i+1/2} ≈ (p_i + p_{i+1})/2 × ((h_i³ + h_{i+1}³)/2)。这种处理相当于对非线性系数做了一次Picard线性化，在每一步迭代里系数用当前压力场计算，然后更新压力。对这个问题收敛性不错，编程也最直观。
第二，为什么用SOR。直接Gauss-Seidel迭代收敛速度太慢，尤其是网格加密到201×201以后，几千步不一定收敛。超松弛（SOR）在更新时叠加一个松弛因子ω：
p_new = (1-ω)·p_old + ω·p_star
对矩形区域，ω一般取1.5~1.8。代码里取1.6。但ω不是越大越快：一开始贪心取过1.95，直接震荡发散。如果发现迭代误差曲线不降反升，先降ω，别硬调网格。
第三，边界和供气孔的顺序。这里有一个踩过的坑：每次迭代完，必须重新强制设置边界压力=pa、供气孔压力=ps。因为内部节点差分时会把边界值也参与进来，如果不强制覆盖，孔口压力会被内部迭代“湮掉”，整个压力场最终变成均匀pa，白跑一场。这个覆盖顺序千万别弄反，我就是早期吃过这个暗亏，排查了大半天。
3.3 承载力、耗气量、刚度计算
压力分布求解之后，轴承特性就有了，这部分其实比解方程更贴近工程需求。
承载力 W = ∫∫(p - pa)dA，矩形网格上就是 sum((p - pa) * dx * dy)，也就是把每个网格单元上的表压叠加起来。注意这里用的是表压，不是绝对压力，因为环境压力已经在大气压下平衡掉了。
耗气量通过流量方程在边界上积分。质量流量 q = -h³/(12μ)·∂p/∂n·ρ，其中边界处密度 ρ = pa/(R·T)。四个边界加起来取绝对值就是总耗气量，单位kg/s。这一步对气源选型和压缩机功耗估算非常重要。
刚度 K = dW/dh。实际做法是把求解过程封装成函数solveBearing(h0, ps)，循环一组膜厚，得到W-h曲线，再用差分近似斜率：
function W = solveBearing(Lx, Ly, h0, ps, pa, mu, nx, ny, orifice, omega, tol, p_init) % 求解单个工况，返回承载力 W % 输入参数与主脚本一致，p_init 可传入上一个工况的收敛压力场 % 内部逻辑与主脚本 2~6 节相同，这里只给函数签名 end h_list = (10:5:50) * 1e-6; W_list = zeros(size(h_list)); for k = 1:length(h_list) W_list(k) = solveBearing(0.1, 0.1, h_list(k), 6e5, 101325, 1.82e-5, ... 101, 101, orifice, 1.6, 1e-6, p); end K_avg = -diff(W_list) ./ diff(h_list); % 平均刚度 N/m封装成函数以后，参数扫描非常方便，批量算膜厚、供气压力、孔径，一条for循环就搞定。这个函数化改造虽然简单，但价值很大，后面讲参数分析时全靠它。
4. 计算结果与参数分析
4.1 压力分布的几个典型特征
以100mm×100mm方形轴承、四角加中心共五个供气孔、供气绝对压力0.6MPa、膜厚20μm为例，压力分布通常呈现“火山群”形态：每个供气孔附近有一个高压尖峰，向外逐渐衰减，最终在四周降到环境压力；相邻孔之间的区域压力高于四周但低于孔中心峰值，说明孔间压力叠加确实存在。
中心孔和四角孔布局要注意一点：中心孔对承载力贡献最大，因为它周围没有边界的“泄漏”面，压力可以维持得比较高；角上的孔一半压力场被边界截断，效率相对低。这也是为什么很多静压轴承设计会把孔排布成梅花形，而不是简单四角分布。
从 surf 图上看，如果某个孔附近压力尖峰特别突兀，而其他区域压力很平，说明孔间距偏大，孔间压力叠加不足。反过来，如果整个压力面像个“平顶山”，孔间距偏小，中间区域压力分布过于均匀，承载力虽然不差，但耗气量可能偏高。
4.2 膜厚、供气压力对轴承特性的影响
计算几组膜厚对比会发现：膜厚从20μm减小到15μm，承载力明显上升，但耗气量显著下降。这个趋势跟经典理论一致：静压轴承的承载力大致与膜厚成反比，刚度随膜厚减小而增大。但要注意，膜厚太小会出现节流孔和间隙匹配失衡的问题，耗气量降低的同时，也可能出现气膜振荡，不能为了刚度无脑压缩间隙。
供气压力提高，承载力近乎线性上升，但耗气量也会上升。工程上有个指标叫“单位耗气量产生的承载力”，做设计时比单一指标更有参考意义。我曾经算过一个极端工况，供气压力从0.4MPa提到0.8MPa，承载力翻了一倍多，但耗气量涨了三倍，如果不是特别需要高刚度，这个方案性价比并不划算。
批量扫描时建议同时输出承载力、刚度、耗气量三个量，画成曲线放在一起看。很多时候单看某个指标会做出错误判断，三个量放一起，取舍关系一目了然。
4.3 网格无关性与收敛性检查
写数值代码，第一件事就是做网格无关性：分别在51×51、101×101、201×201网格下算同一个工况，对比承载力。如果101×101和201×201的结果差小于1%~2%，说明网格已经足够密，再加密只是浪费时间。
收敛性则看迭代误差曲线：理想情况是单调下降。如果曲线出现平台或周期性波动，说明松弛因子偏大，或者供气孔附近压力梯度太陡，需要局部加密网格。
我通常同时打印误差和承载力两个量，迭代过程中如果承载力已经稳定但误差还在降，说明网格精度足够，可以提前终止。这比死等误差阈值到1e-8要高效得多，尤其是批量扫描时，能省一半时间。
5. 新手最容易踩的坑与调试经验
5.1 迭代不收敛或收敛慢
最常见的原因有三个：松弛因子太大、初值给得太离谱、供气孔压力与边界压力差异过大。
初值我一般取环境压力pa，而不是0，这样迭代初期的压力梯度比较平缓，不容易震荡。供气
本文还有配套的精品资源，点击获取
