气体润滑数值计算方法与程序

根据黄平老师的《润滑数值计算方法》第七章内容，推导气体润滑数值计算方法，并用Matlab将程序复现，最后和书中的计算结果及其论文进行对比。

对于等温过程，径向气体动压轴承中的气膜压力分布计算公式为：

（1）

公式中：R为轴颈半径，单位为m；p为油膜压力，单位为Pa，h为油膜厚度，单位为m；U为轴颈速度，单位为m/s，η为气体黏度，单位为Pa*s。

将公式进行无量纲化，取：

（2）

式中c为轴承的半径间隙，单位为m，pa为大气压，为1.013x10^5Pa。

将式（2）代入式（1）得：

（3）

利用差分法对式（3）进行求解：

求解的边界条件：

根据黄平老师书中的参数：

轴承宽度：B=60.0E-3 m；

轴颈半径：R=25.0E-3 m；

半径间隙：C0=5.0E-5 m；

工作速度：AN=6.0E4 r/min；

环境压力：PA=1.013E5 Pa；

气体动力黏度：EPA=1.79E-5 Pa*s；

偏心率：EPSON=0.7；

无量纲油膜压力计算结果如下：

% 气体轴承计算///润滑数值计算方法
clear;
clc;
%% 轴承结构参数
B=60E-3;     % 轴承宽度m
R=25E-3;    % 轴承半径m
C0=5E-5;    % 半径间隙m
%% 工作参数
AN=6E4;        % 工作速度r/min
PA=1.013E5;    % 标准大气压Pa
EDA=1.79E-5;    % 气体动力粘度Pa*s
EPSON=0.7;     % 轴颈偏心率
%% 网格参数
N=79;   % 周向网格点数量
M=22;   % 轴向网格点数量
element_area=2*pi*R/(N-1)*B/(M-1);
%% 气体轴承部分计算
[P,H]=gas_journal(B,R,C0,AN,PA,EDA,EPSON,N,M);
%%
function [P,H]=gas_journal(B,R,C0,AN,PA,EDA,EPSON,N,M)
DX=2*pi/(N-1);    % 无量纲周向两个网格点之间的距离
DY=1/(M-1);      % 无量纲轴向两个网格点之间的距离
OMEGA=AN*2*pi/60;    % 转速r/min-角速度rad/s
U=OMEGA*R;      % 角速度rad/s-线速度m/s
ALENDA=6*EDA*U*R/PA/C0^2;   % 无量纲参数
ALFA=(R/B*DX/DY)^2;
H=SUBH(N,M,DX,DY,EPSON);
[P,IK]=SUBP(N,M,DX,DY,EPSON,ALFA,ALENDA,H);
end
%% 
function H=SUBH(N,M,DX,DY,EPSON)
H=zeros(N,M);
for I=1:N
    %% I为周向
    SETA=(I-1)*DX;
    for J=1:M
        H(I,J)=1+EPSON*cos(SETA);
    end
end
end
%%
function [P,IK]=SUBP(N,M,DX,DY,EPSON,ALFA,ALENDA,H)
global C1
P=ones(N,M);
IK=0;   % 迭代次数
C1=1;
for I=2:N-1
    for J=2:M-1
        P(I,J)=1.1;
    end
end
while C1>1E-12
    C1=0;
    ALOAD=0;
    for I=2:N-1
        I1=I-1;
        I2=I+1;
        for J=2:M-1
            J1=J-1;
            J2=J+1;
            PD=P(I,J);
            A1=(0.5*(H(I2,J)+H(I,J)))^3;   % I+0.5
            A2=(0.5*(H(I,J)+H(I1,J)))^3;   % I-0.5
            A3=ALFA*(0.5*(H(I,J2)+H(I,J)))^3;   % J+0.5
            A4=ALFA*(0.5*(H(I,J)+H(I,J1)))^3;   % J-0.5
            P(I,J)=(-DX*ALENDA*(P(I+1,J)*H(I+1,J)-P(I-1,J)*H(I-1,J))+A1*P(I2,J)^2+A2*P(I1,J)^2+A3*P(I,J2)^2+A4*P(I,J1)^2)/(A1+A2+A3+A4);
            if P(I,J)<0
                P(I,J)=0;
            end
            P(I,J)=sqrt(P(I,J));
            P(I,J)=0.7*PD+0.3*P(I,J);
            if P(I,J)<0
                P(I,J)=0;
            end
            C1=C1+abs(P(I,J)-PD);
            ALOAD=ALOAD+P(I,J);
        end
    end
    IK=IK+1;
    C1=C1/ALOAD
end
end

可以看到计算程序中圆周方向上首尾没有相连，将圆周方向上第一个网格节点和最后一个网格节点相连后：

无量纲油膜压力计算结果如下：

无量纲承载力的计算公式为：

根据虞烈老师的书《可压缩气体润滑与弹性箔片气体轴承技术》中的P64页中的参数：轴承数为0.5，长径比L/D=0.5的360°圆轴承；

进行计算的结果如下：

偏心率	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
偏位角	87.23	86.83	86.02	84.54	81.79	76.47	65.78	46.47	22.73
无量纲承载力	0.0120	0.0249	0.0398	0.0583	0.0832	0.1200	0.1833	0.3313	0.8762