线性拟合与润滑求解.zip_MATLAB 滑动轴承_ontoj8b_有限差分 压力_润滑数值计算_滑动轴承

function q1(eps,deps,Bd,m,n) %有限差分法计算有限长径向滑动轴承挤压膜压力分布及承载量 %偏心率：eps %挤压速度：deps %轴承径宽比：Bd %有限差分网格划分，轴承周向方向网格数：m %有限差分网格划分，轴承轴向方向网格数：n zspan=[-1,1];fispan=[0,pi];%区间范围 delfi=(fispan(2)-fispan(1))/m; fi=fispan(1):delfi:fispan(2); delz=(zspan(2)-zspan(1))/n; z=zspan(1):delz:zspan(2); ndeps=length(deps);%deps向量长度 %%%%%%%%%%%%%%%%% pf=zeros(m+1,n+1); H=zeros((m-1)*(n-1),(m-1)*(n-1));%压力系数矩阵 g=zeros(m-1,n-1); %%%%%%%%%%%%%边值条件 pf(1,:)=0; pf(m+1,:)=0; pf(:,1)=0; pf(:,n+1)=0; %%%%%%%%%%%%%%%%% p=pf([2:m],[2:n]);%解 %%%%%%%%%%%%%%% %求有限差分法计算方程的系数 for a=1:ndeps %取ndeps组deps(dε/dt)值，分别求解对应压力分布和承载力 for j=1:n-1 for i=1:m-1 A1=1; B1=Bd^2; C1=-3*eps*sin(fi(i))/(1+eps*cos(fi(i))); E1=cos(fi(i))*deps(a)/(1+eps*cos(fi(i)))^3; K=2*(A1/delfi^2+B1/delz^2); F=E1/K; A=(A1/delfi^2+C1/(2*delfi))/K; B=(A1/delfi^2-C1/(2*delfi))/K; C=B1/delz^2/K; D=B1/delz^2/K; E=1; %%%%%%%%%%%%%%%%% %构建压力系数矩阵 g(i,j)=F; ii=(j-1)*(m-1)+i; H(ii,ii)=-E; if mod(ii,(m-1))==1 g(i,j)=g(i,j)-B*pf(1,j);%右端项组装 if i>1 H(ii,ii-1)=0;%MATLAB中矩阵行列与一般写法相反 end else H(ii,ii-1)=B;%MATLAB中矩阵行列与一般写法相反 end if mod(ii,(m-1))==0 g(i,j)=g(i,j)-A*pf(m-1,j);%右端项组装 if j<n-1 H(ii,ii+1)=0;%MATLAB中矩阵行列与一般写法相反 end else H(ii,ii+1)=A;%MATLAB中矩阵行列与一般写法相反 end if j>1 H(ii,ii-m+1)=D; else g(i,j)=g(i,j)-D*pf(i,1); end; if j<n-1 H(ii,ii+m-1)=C; else g(i,j)=g(i,j)-C*pf(i,n-1); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %用雅克比迭代求解压力系数矩阵 p(:)=in_jacobi(H,g(:),eye(length(g(:)),1)); pf([2:m],[2:n])=p; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %将负压各点置零 for j=1:n+1 for i=1:m+1 if pf(i,j)<0 pf(i,j)=0; end end end %%数值积分求承载力 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sum0=0; sum1=0; sum3=0; sum4=0; sum5=0; sum6=0; for i=1:m+1 for k=2:2:n sum0=sum0+pf(i,k); end for k=3:2:n-1 sum1=sum1+pf(i,k); end G(i)=delz*(pf(i,1)+4*sum1+2*sum0+pf(i,n+1))/3; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=2:2:m sum3=sum3+G(i)*sin(fi(i)); end for i=3:2:m-1 sum4=sum4+G(i)*sin(fi(i)); end for i=2:2:m sum5=sum5+G(i)*cos(fi(i)); end for i=3:2:m-1 sum6=sum6+G(i)*cos(fi(i)); end F1=-delfi*(G(1)*sin(fi(1))+4*sum3+2*sum4+G(m+1)*sin(fi(m+1)))/3; F2=-delfi*(G(1)*cos(fi(1))+4*sum5+2*sum6+G(m+1)*cos(fi(m+1)))/3; F(a)=sqrt(F1^2+F2^2); Fw(a)=F(a) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %绘制挤压速度dε/dt与承载力Fw之间的关系曲线 figure; a=1:ndeps; plot(deps,Fw(a)); xlabel('挤压速度dε/dt'); ylabel('承载力Fw'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %绘制压力三维分布图 X=1:1:(n+1); Y=(1:1:(m+1))/(m+1)*pi; Y=Y'; Z=pf([1:m+1],[1:n+1]); figure; surf(X,Y,Z); ylabel('圆周方向'); xlabel('轴向方向'); zlabel('量纲一的油膜压力'); axis([-inf,inf,-inf,inf,-inf,inf]); %set(gcf,'color','white'); %shading faceted grid on Colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %雅克比迭代函数 function x=in_jacobi(A,b,x1) exps=1.0e-10; count=0; n=length(x1); for row=1:n, A(row,[1:row-1,row+1:n])=-A(row,[1:row-1,row+1:n])/A(row,row); b(row,1)=b(row,1)/A(row,row); A(row,row)=0; end x=x1+2*exps; while max(abs(x-x1))>exps & count<10000, x=A*x1+b; temp=x1; x1=x; x=temp; count=count+1; if x<0 x=0; end; end; in_jacobi=x; disp('叠代次数是:'); count