基于matlab模拟轴承混合润滑模型.zip

%%动态特性分析 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