✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
混合润滑向心滑动轴承的优越性远非滚动轴承可比拟,但由于影响边界膜强度的因素十分复杂,所以其现今还缺乏有系统理论依据的计算方法。对此,可根据向心滑动轴承的工作状况,并以维持边界润滑为计算准则,对其进行间接的、条件性的验算。
⛄ 代码
%%动态特性分析
clc;clear;clf
global ps;
global row;
global pa;
global d;
global v;
D0=0.022; %轴颈直径,m
R=D0/2;%轴颈半径，m
c=0.00002; %平均气膜厚度，m
D2=0.02203; %轴承内径,m
L=0.03; %轴承长度,m
u=17.5e-6; %空气的粘度 Pa.s
omiga=40000*2*pi/60; % 轴转速 rad/s
pa=101325;%环境压力,pa
row=1.293;%空气密度，kg/m^3
R1=287.24;%气体常数，N*m(kg*k)^-1
ps=5;%无量纲供气压力，供气压力P0=pa*Ps
lamda=(6*u*omiga*(R/c)^2)/pa;%无量纲参数,类似的还有X=Rx,Y=Ly，h=Hc
mc=24*u*R^2/(c^3*pa*row);%流量因子
mc1=-row*pa*c^3/(24*u);%计算流量用到的
mc2=row*omiga*R*c/2;%计算流量用到的
%% 设置参数
e=0.5; %轴颈偏心率
d=0.0003;%节流孔直径，m
N=10;%单排节流孔数量
l=2/10;%节流孔至轴承端面距离与轴承长度的比值
m=60; %周向等分数；
n=60; %轴向等分数；
N0=m/N;
n1=n*l;%第一排节流孔位置
n2=n*(1-l)+2;%第二排节流孔位置
%% 对x,y归一化
deltay=1/n; %轴向等分间距
deltax=2*pi/m; %周向角度等分大小
%% 计算各点处气膜厚度
for i=1:n+1 %轴向
for j=1:m+1%周向
H(i,j)=1-e*cos((j-1)*deltax);%气膜厚度,在最小气膜厚度处展开
end
end
%%
%雷诺方程迭代公式：A(i,j)*(P(i,j+1))^2+B(i,j)*(P(i,j-1))^2+C(i,j)*(P(i+1,j))^2+D(i,j)*(P(i-1,j))^2-E(i,j)*P(i,j)^2=lamda*deltax*(P(i,j+1)*H(i,j+1)-P(i,j-1)*H(i,j-1))
%% 迭代系数
for i=1:n+1
for j=1:m+1
%最小气膜处展开的迭代系数
A(i,j)=H(i,j)^3+(3/2)*H(i,j)^2*deltax*e*sin((j-1)*deltax);
B(i,j)=H(i,j)^3-(3/2)*H(i,j)^2*deltax*e*sin((j-1)*deltax);
C(i,j)=H(i,j)^3*(deltax/deltay)^2*(R/L)^2;
D(i,j)=H(i,j)^3*(deltax/deltay)^2*(R/L)^2;
E(i,j)=2*H(i,j)^3+2*H(i,j)^3*(deltax/deltay)^2*(R/L)^2;
end
end
%% 迭代过程计算
ERR=1.0e-5; %误差限
%各点赋初值
P0=ones(n+1,m+1); %最小气膜处展开
P1=ones(2,N+1);%节流器出口压力pd
Q=ones(2,N+1);%节流器出口流量
vv=ones(2,N+1);%节流器出口压力pd'
beta=ones(2,N+1);%beta=Pd/Ps，节流器出口压力与供气压力比值，也称节流压力比
%迭代计算
ks=1;%迭代计数
while ks>0 %流量收敛
k=1;%迭代计数
P1(1,1:N+1)=vv(1,1:N+1);%节流器出口压力
P1(2,1:N+1)=vv(2,1:N+1);%节流器出口压力
%P0=ones(n+1,m+1);
while k>0%最小气膜厚度处展开，计算气膜压力分布，无量纲
P=P0; %下次迭代赋值
for i=2:n%i==1和i==n+1为大气边界
for j=1:m+1%j==1和j==m+1为重合边界
P0(n1,N0*(0:N)+1)=P1(1,1:N+1);%节流孔边界
P0(n2,N0*(0:N)+1)=P1(2,1:N+1);%节流孔边界
if j==1
P0(i,j)=((A(i,j)*(P(i,j+1))^2+B(i,j)*(P0(i,m))^2+C(i,j)*(P(i+1,j))^2+D(i,j)*(P0(i-1,j))^2-lamda*deltax*(P(i,j+1)*H(i,j+1)-P0(i,m)*H(i,m)))/E(i,j))^0.5;%重合边界
else
if j==m+1
P0(i,j)=((A(i,j)*(P(i,2))^2+B(i,j)*(P0(i,j-1))^2+C(i,j)*(P(i+1,j))^2+D(i,j)*(P0(i-1,j))^2-lamda*deltax*(P(i,2)*H(i,2)-P0(i,j-1)*H(i,j-1)))/E(i,j))^0.5;%重合边界
else
P0(i,j)=-0.2*P(i,j)+1.2*((A(i,j)*(P(i,j+1))^2+B(i,j)*(P0(i,j-1))^2+C(i,j)*(P(i+1,j))^2+D(i,j)*(P0(i-1,j))^2-lamda*deltax*(P(i,j+1)*H(i,j+1)-P0(i,j-1)*H(i,j-1)))/E(i,j))^0.5;%计算各点压力值,超松弛迭代
%P0(i,j)=((A(i,j)*(P(i,j+1))^2+B(i,j)*(P0(i,j-1))^2+C(i,j)*(P(i+1,j))^2+D(i,j)*(P0(i-1,j))^2-lamda*deltax*(P(i,j+1)*H(i,j+1)-P0(i,j-1)*H(i,j-1)))/E(i,j))^0.5;%计算各点压力值，高斯-赛德尔迭代
end
end
if P0(i,j)<0
P0(i,j)=0; %气膜破裂，置零,空穴现象
break;
end
end
end
%% 误差控制
sum1=0;
sum2=0;
for s=1:n+1
for t=1:m+1
sum1=sum1+abs((P0(s,t)-P(s,t)));%相邻两次迭代间的偏差
sum2=sum2+abs(P(s,t));
end
end
sum=sum1/sum2;%相对误差
if sum<=ERR%判断收敛
break;
end
k=k+1;
%% 遍历各处 i ,j
end
%%计算节流孔流出流量
beta(1,1:N+1)=P(n1,N0*(0:N)+1)/ps;
beta(2,1:N+1)=beta(1,1:N+1);
%节流孔流量方程
if beta(1,1:N+1)>0.528
Q(1,1:N+1)=0.8*(pi*d^2/4)*ps*(2*row*pa).^0.5*(3.5*(beta(1,1:N+1).^(10/7)-beta(1,1:N+1).^(12/7))).^0.5;%节流孔出口质量流量
Q(2,1:N+1)=Q(1,1:N+1);%节流孔出口质量流量
else
Q(1,1:N+1)=0.8*(pi*d^2/4)*ps*(2*row*pa).^0.5*0.7^0.5*(5/6)^3;%节流孔出口质量流量
Q(2,1:N+1)=Q(1,1:N+1);%节流孔出口质量流量
end
%%计算每个节流孔流入间隙流量，取包围节流孔处deltax*deltay面积的区域
%赋初值
M1=ones(2,N+1);%x方向
M2=ones(2,N+1);%z方向
M3=ones(2,N+1);%-x方向
M4=ones(2,N+1);%-z方向
for j1=1:N
N1=N0*(j1-1)+1;
%节流孔出口与进入轴承间隙的流量差
%M01(1,j1)=row*c*deltax*deltay*H(n1,N1)*P(n1,N1);
%M01(2,j1)=row*c*deltax*deltay*H(n2,N1)*P(n2,N1);
%进入轴承间隙流量,M1为x正向，M3为x负向，M2为y正向，M4为y负向
if j1==1
%第一排
M1(1,j1)=(mc1/R)*(deltay/deltax)*(P(n1,N1+1)^2-P(n1,N1)^2)*((H(n1,N1+1)+H(n1,N1))/2)^3+mc2*deltay*((H(n1,N1+1)+H(n1,N1))/2)*((P(n1,N1+1)+P(n1,N1))/2);
M2(1,j1)=(mc1/L)*(deltax/deltay)*(P(n1+1,N1)^2-P(n1,N1)^2)*((H(n1+1,N1)+H(n1,N1))/2)^3;
M3(1,j1)=(mc1/R)*(deltay/deltax)*(P(n1,m)^2-P(n1,N1)^2)*((H(n1,m)+H(n1,N1))/2)^3+mc2*deltay*((H(n1,m)+H(n1,N1))/2)*((P(n1,m)+P(n1,N1))/2);
M4(1,j1)=(mc1/L)*(deltax/deltay)*(P(n1-1,N1)^2-P(n1,N1)^2)*((H(n1-1,N1)+H(n1,N1))/2)^3;
%第二排
M1(2,j1)=M1(1,j1);
M2(2,j1)=M2(1,j1);
M3(2,j1)=M3(1,j1);
M4(2,j1)=M4(1,j1);
else
M1(1,j1)=(mc1/R)*(deltay/deltax)*(P(n1,N1+1)^2-P(n1,N1)^2)*((H(n1,N1+1)+H(n1,N1))/2)^3+mc2*deltay*((H(n1,N1+1)+H(n1,N1))/2)*((P(n1,N1+1)+P(n1,N1))/2);
M2(1,j1)=(mc1/L)*(deltax/deltay)*(P(n1+1,N1)^2-P(n1,N1)^2)*((H(n1+1,N1)+H(n1,N1))/2)^3;
M3(1,j1)=(mc1/R)*(deltay/deltax)*(P(n1,N1-1)^2-P(n1,N1)^2)*((H(n1,N1-1)+H(n1,N1))/2)^3+mc2*deltay*((H(n1,N1-1)+H(n1,N1))/2)*((P(n1,N1-1)+P(n1,N1))/2);
M4(1,j1)=(mc1/L)*(deltax/deltay)*(P(n1-1,N1)^2-P(n1,N1)^2)*((H(n1-1,N1)+H(n1,N1))/2)^3;
M1(2,j1)=M1(1,j1);
M2(2,j1)=M2(1,j1);
M3(2,j1)=M3(1,j1);
M4(2,j1)=M4(1,j1);
end
end
M0=(M1+M2+M3+M4);
M0(1:2,N+1)=M0(1:2,1);
%%判断流量是否收敛
for kk=1:2
for kk1=1:N+1
Q1(kk,kk1)=abs(Q(kk,kk1)-M0(kk,kk1))/M0(kk,kk1);
end
end
Q2=max(max(Q1));
if Q2<1
break;
end
%%修正节流孔出口压力
%已知每个节流孔处进入间隙流量，用最小二乘法求出每个节流孔处的出口压力pd'
for k1=1:2
for k2=1:11
k0=1;a=1;b=5;tol=0.0001;
v(k1,k2)=M0(k1,k2);
fa=fun(a);
while k0>0
t(k1,k2)=(a+b)/2;
ft=fun(t(k1,k2));
if ft==0|(b-a)/2<tol
break;
end
if ft(k1,k2).*fa(k1,k2)>0
a=t(k1,k2);
else
b=t(k1,k2);
end
k0=k0+1;
end
vv(k1,k2)=t(k1,k2);%节流孔处的出口压力pd'
end
end
ks=ks+1;
end
%%绘图
sitas=(0:m)*deltax; %周向
Ls=(0:n)*deltay; %轴向
[X,Y]=meshgrid(Ls,sitas);
mesh(X,Y,P')%画网格
title '气膜内压力分布'
xlabel '无量纲轴向'
ylabel '无量纲周向'
zlabel '无量纲压力'
function f=fun(x)
global ps;
global row;
global pa;
global d;
global v;
f=v-0.8*(pi*d^2/4)*ps*(2*row*pa)^0.5*(3.5*((x/ps)^(10/7)-(x/ps)^(12/7)))^0.5;
end
⛄ 运行结果
⛄ 参考文献
[1] 周广武. 水润滑橡胶合金轴承混合润滑分析与动力学性能优化. 重庆大学, 2013.
[2] 魏立队, 段树林, 武起立,等. 基于柔性整机体模型的柴油机主轴承混合润滑特性研究[J]. 润滑与密封, 2012, 000(010):16-21.
[3] 王子嘉, 茆志伟, 张进杰,等. 基于混合润滑模型的连杆大端轴承接触特性研究.
[4] 王家序, 倪小康, 韩彦峰,等. 轴向往复运动下微槽轴承混合润滑数值模型[J]. 吉林大学学报(工学版), 2019, v.49;No.203(03):221-229.
