径向气体轴承动态特性求解

凌晨幻想

已于 2024-06-23 19:04:18 修改

阅读量444

点赞数 10

文章标签：数学建模考研

于 2024-06-23 19:03:59 首次发布

本文链接：https://blog.csdn.net/m0_53100062/article/details/139899952

版权

根据虞烈老师的《可压缩气体润滑与弹性箔片气体轴承技术》中的2.3求解动态气膜压力分布以及刚性表面气体动压轴承的刚度和阻尼，并和计算结果进行对比误差较小。

动态气膜压力Pd0的动态Reynolds方程为：

（1）

式中，P0和H0均为静态时的气膜压力和气膜厚度，而动态气膜压力Pd0为未知量，是扰动量E0和Θ0的函数。动态气膜压力对于扰动的变化率PE和PΘ在复数范围内被定义为：

（2）

将式（1）分别对扰动量E0求偏导：

（3）

对式（3）进行整理：

（4）

其中：

（5）

利用差分法将式（4）分解：

（6）

对上式进行求解得到PE，同样的方法可以求得PΘ。

求解时边界条件为：

当Ω=1时，PE的实数区域分布为：

PE的复数区域分布为：

PΘ的实数区域分布为：

PΘ的虚数区域分布为：

然后根据书中P60的公式即可求得部分刚度和阻尼参数。

（7）

式（7）中θ0为偏位角。

部分代码请移步“气体润滑数值计算方法与程序”，作者已对计算结果和书中的结果进行对比，其他人请自行对比。

function [P,H,Kxx_xy,Kyx_yy,D_xx_xy,D_yx_yy]=gas_journal(B,R,C0,AN,PA,EDA,EPSON,N,M)
DX=2*pi/(N-1);    % 无量纲周向两个网格点之间的距离
DY=B/R/(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=(DX/DY)^2;
H=SUBH(N,M,DX,DY,EPSON);
[P,IK]=SUBP(N,M,DX,DY,EPSON,ALFA,ALENDA,H);
%% 求解偏位角
Wxn=0;
Wyn=0;
for I=1:N
    SETA=(I-1)*DX;
    for J=1:M
        Wxn=Wxn+(P(I,J)-1)*sin(SETA)*DX*DY;
        Wyn=Wyn-(P(I,J)-1)*cos(SETA)*DX*DY;
    end
end
W=sqrt(Wxn^2+Wyn^2);
phi=atand(Wxn/Wyn);
%% 无量纲轴颈扰动频率
omiga=1;
[PE1,PE2]=dynamic_parameter(P,H,N,M,DX,DY,ALENDA,omiga);
PE1_real=real(PE1);
PE1_imag=imag(PE1);
PE2_real=real(PE2);
PE2_imag=imag(PE2);
%% 求解
Ky_cc=0;
Kx_cc=0;
Dx_cc=0;
Dy_cc=0;
for I=1:N
    SETA=phi/(360)*2*pi+(I-1)*DX;
    for J=1:M
        % 刚度
        Ky_cc=Ky_cc+(-PE1_real(I,J)*cos(SETA)*DX*DY);
        Kx_cc=Kx_cc+(-PE1_real(I,J)*sin(SETA)*DX*DY);
        % 阻尼
        Dy_cc=Dy_cc+(-PE1_imag(I,J)*cos(SETA)*DX*DY);
        Dx_cc=Dx_cc+(-PE1_imag(I,J)*sin(SETA)*DX*DY);
    end
end
Dy_cc=Dy_cc/omiga;
Dx_cc=Dx_cc/omiga;
%%
Ky_oo=0;
Kx_oo=0;
Dx_oo=0;
Dy_oo=0;
for I=1:N
    SETA=phi/(360)*2*pi+(I-1)*DX;
    for J=1:M
        % 刚度
        Ky_oo=Ky_oo+(-PE2_real(I,J)*cos(SETA)*DX*DY);
        Kx_oo=Kx_oo+(-PE2_real(I,J)*sin(SETA)*DX*DY);
        % 阻尼
        Dy_oo=Dy_oo+(-PE2_imag(I,J)*cos(SETA)*DX*DY);
        Dx_oo=Dx_oo+(-PE2_imag(I,J)*sin(SETA)*DX*DY);
    end
end
Dy_oo=Dy_oo/omiga;
Dx_oo=Dx_oo/omiga;
% 坐标转换
shift_matrix=[sind(phi),cosd(phi);cosd(phi),-sind(phi)];
Kxx_xy=shift_matrix*[Kx_cc;Kx_oo];
Kyx_yy=shift_matrix*[Ky_cc;Ky_oo];
D_xx_xy=shift_matrix*[Dx_cc;Dx_oo];
D_yx_yy=shift_matrix*[Dy_cc;Dy_oo];
end

function [AAA,BBB,CCC,DDD,EEE,FFF1,FFF2]=paramter_cal(P,H,N,M,DX,DY,ALENDA,omiga)
%% 初始化
AAA=zeros(N,M);
BBB=zeros(N,M);
CCC=zeros(N,M);
DDD=zeros(N,M);
EEE=zeros(N,M);
FFF1=zeros(N,M);
FFF2=zeros(N,M);
FFF3=zeros(N,M);
%%
H0=H;
P0=P;
%% 1/H0
one_H0=1./H0;
%% 
for I=1:N
    for J=2:M-1
        if I==1
            I1=N-1;
            I2=I+1;
        elseif I==N
            I1=I-1;
            I2=2;
        else
            I1=I-1;
            I2=I+1;
        end
        J1=J-1;
        J2=J+1;
        phi=(I-1)*DX;
        H0_phi=(H0(I2,J)-H0(I1,J))/(2*DX);
        H0_Y=(H0(I,J2)-H0(I,J1))/(2*DY);
        P0_phi=(P0(I2,J)-P0(I1,J))/(2*DX);
        P0_Y=(P0(I,J2)-P0(I,J1))/(2*DY);
        P0_phi_2=(P0(I2,J)+P0(I1,J)-2*P0(I,J))/(DX)^2;
        P0_Y_2=(P0(I,J2)+P0(I,J1)-2*P0(I,J))/(DY)^2;
        %%
        one_H0_phi=(one_H0(I2,J)-one_H0(I1,J))/(2*DX);
        one_H0_Y=(one_H0(I,J2)-one_H0(I,J1))/(2*DY);
        %% 参数计算
        AAA(I,J)=3*H0(I,J)^2*H0_phi*P0_phi+H0(I,J)^3*P0_phi_2+3*H0(I,J)^2*H0_Y*P0_Y+H0(I,J)^3*P0_Y_2-ALENDA*H0_phi-2i*ALENDA*H0(I,J)*omiga;
        BBB(I,J)=2*H0(I,J)^3*P0_phi+3*P0(I,J)*H0(I,J)^2*H0_phi-ALENDA*H0(I,J);
        CCC(I,J)=2*H0(I,J)^3*P0_Y+3*P0(I,J)*H0(I,J)^2*H0_Y;
        DDD(I,J)=P0(I,J)*H0(I,J)^3;
        EEE(I,J)=DDD(I,J);
        FFF1(I,J)=ALENDA*P0_phi*cos(phi)-ALENDA*sin(phi)*P0(I,J)+2i*ALENDA*omiga*P0(I,J)*cos(phi)-...
            3*H0(I,J)^2*cos(phi)*(P0_phi^2+P0_Y^2)-3*P0(I,J)*H0(I,J)^2*(-sin(phi)*P0_phi)-6*P0(I,J)*H0(I,J)*cos(phi)*(H0_phi*P0_phi+H0_Y*P0_Y)-3*P0(I,J)*H0(I,J)^2*cos(phi)*(P0_phi_2+P0_Y_2);
        FFF2(I,J)=ALENDA*P0_phi*sin(phi)+ALENDA*cos(phi)*P0(I,J)+2i*ALENDA*omiga*P0(I,J)*sin(phi)-...
            3*H0(I,J)^2*sin(phi)*(P0_phi^2+P0_Y^2)-3*P0(I,J)*H0(I,J)^2*(cos(phi)*P0_phi)-6*P0(I,J)*H0(I,J)*sin(phi)*(H0_phi*P0_phi+H0_Y*P0_Y)-3*P0(I,J)*H0(I,J)^2*sin(phi)*(P0_phi_2+P0_Y_2);
%         FFF3(I,J)=ALENDA*cos(phi)*P0_phi-ALENDA*sin(phi)*P0(I,J)+2i*ALENDA*omiga*P0(I,J)*cos(phi)-...
%             3*ALENDA*cos(phi)*(P0_phi+P0(I,J)/H0(I,J)*H0_phi)-3*P0(I,J)*H0(I,J)^3*(P0_phi*(-sin(phi)/H0(I,J)+cos(phi)*one_H0_phi)+cos(phi)*P0_Y*one_H0_Y);
    end
end
end

function [PE1,PE2]=dynamic_parameter(P,H,N,M,DX,DY,ALENDA,omiga)
[AAA,BBB,CCC,DDD,EEE,FFF1,FFF2]=paramter_cal(P,H,N,M,DX,DY,ALENDA,omiga);
%% HE、PE
C1=1;
PE1=zeros(N,M);
while C1>1E-12
    C1=0;
    ALOAD=0;
    for I=1:N
        for J=2:M-1
            if I==1
                I1=N-1;
                I2=I+1;
            elseif I==N
                I1=I-1;
                I2=2;
            else
                I1=I-1;
                I2=I+1;
            end
            J1=J-1;
            J2=J+1;
            PD=PE1(I,J);
            A1=2*DDD(I,J)/(DX)^2+2*EEE(I,J)/(DY)^2-AAA(I,J);
            A2=DDD(I,J)/(DX)^2+BBB(I,J)/(2*DX);
            A3=DDD(I,J)/(DX)^2-BBB(I,J)/(2*DX);
            A4=EEE(I,J)/(DY)^2+CCC(I,J)/(2*DY);
            A5=EEE(I,J)/(DY)^2-CCC(I,J)/(2*DY);
            PE1(I,J)=(A2*PE1(I2,J)+A3*PE1(I1,J)+A4*PE1(I,J2)+A5*PE1(I,J1)-FFF1(I,J))/A1;
%             if PE1(I,J)<0
%                 PE1(I,J)=0;
%             end
            PE1(I,J)=0.7*PD+0.3*PE1(I,J);
            C1=C1+abs(PE1(I,J)-PD);
            ALOAD=ALOAD+abs(PE1(I,J));
        end
    end
    C1=C1/ALOAD
end
%% HE、PE
C1=1;
PE2=zeros(N,M);
while C1>1E-12
    C1=0;
    ALOAD=0;
    for I=1:N
        for J=2:M-1
            if I==1
                I1=N-1;
                I2=I+1;
            elseif I==N
                I1=I-1;
                I2=2;
            else
                I1=I-1;
                I2=I+1;
            end
            J1=J-1;
            J2=J+1;
            PD=PE2(I,J);
            A1=2*DDD(I,J)/(DX)^2+2*EEE(I,J)/(DY)^2-AAA(I,J);
            A2=DDD(I,J)/(DX)^2+BBB(I,J)/(2*DX);
            A3=DDD(I,J)/(DX)^2-BBB(I,J)/(2*DX);
            A4=EEE(I,J)/(DY)^2+CCC(I,J)/(2*DY);
            A5=EEE(I,J)/(DY)^2-CCC(I,J)/(2*DY);
            PE2(I,J)=(A2*PE2(I2,J)+A3*PE2(I1,J)+A4*PE2(I,J2)+A5*PE2(I,J1)-FFF2(I,J))/A1;
%             if PE2(I,J)<0
%                 PE2(I,J)=0;
%             end
            PE2(I,J)=0.7*PD+0.3*PE2(I,J);
            C1=C1+abs(PE2(I,J)-PD);
            ALOAD=ALOAD+abs(PE2(I,J));
        end
    end
    C1=C1/ALOAD
end
end

凌晨幻想

关注

10
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
径向气体轴承动态特性求解

式中，P0和H0均为静态时的气膜压力和气膜厚度，而动态气膜压力Pd0为未知量，是扰动量E0和Θ0的函数。根据虞烈老师的《可压缩气体润滑与弹性箔片气体轴承技术》中的2.3求解动态气膜压力分布以及刚性表面气体动压轴承的刚度和阻尼，并和计算结果进行对比误差较小。部分代码请移步“气体润滑数值计算方法与程序”，作者已对计算结果和书中的结果进行对比，其他人请自行对比。然后根据书中P60的公式即可求得部分刚度和阻尼参数。对上式进行求解得到PE，同样的方法可以求得PΘ。式（7）中θ0为偏位角。
复制链接

扫一扫