根据虞烈老师的《可压缩气体润滑与弹性箔片气体轴承技术》中的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