利用有限差分法求解雷诺方程_求解气体静压轴承的特性_matlab源码

以下用matlab编写的五点差分法离散雷诺方程用SOR(超松弛迭代法)求解压力分布的源程序： clear all r1=31e-3;%推力瓦内径 r2=71e-3;%推力瓦外径 h1=15e-6;%假设的初始最小水膜厚度 h2=30e-6;%假设的初始最大水膜厚度 theta1=50*pi/180;%推力瓦的夹角 beita=0.001*pi/180;%推力瓦周向倾角 k=6;%推力瓦块数 F1=3000;%轴向推力 n1=1000;%轴承转速 omiga=2*pi*n1/60;%推力轴承角速度 miu=0.001;%水的粘度 l=((r1+r2)/2)*theta1;%轴瓦长度 P0=1.0064e+06;%压力量纲系数 nx=40;%周向划分的网格数 ny=40;%径向划分的网格数 deltax=(50/nx)*(pi/180);%周向网格的单位长度 deltay=(r2-r1)/(r1*ny);%径向网格的单位长度 r=zeros(nx+1,ny+1);%初始化无量纲半径矩阵 theta=zeros(nx+1,ny+1);%初始化无量纲角度矩阵 H=zeros(nx+1,ny+1);%初始化无量纲水膜厚度矩阵 while 1 for i=1:nx+1 for j=1:ny+1 r(i,j)=1+(j-1)*deltay;%无量纲半径矩阵 theta(i,j)=(i-1)*deltax;%无量纲角度矩阵 H(i,j)=h2/h1-((h2-h1)/h1)*(theta(i,j)/theta1);%可倾瓦水膜厚度无量纲表达式 end end h=h1*H;%有量纲润滑膜厚度矩阵 P=zeros(nx+1,ny+1);%初始化无量纲节点压力 A=zeros(nx+1,ny+1);%初始化系数矩阵A B=zeros(nx+1,ny+1);%初始化系数矩阵B C=zeros(nx+1,ny+1);%初始化系数矩阵C D=zeros(nx+1,ny+1);%初始化系数矩阵D E=zeros(nx+1,ny+1);%初始化系数矩阵E F=zeros(nx+1,ny+1);%初始化系数矩阵F P2=zeros(nx+1,ny+1);%初始化无量纲压力矩阵P2 P3=zeros(nx+1,ny+1);%初始化无量纲压力矩阵P3 S=0;%初始化误差表达式的分子 T=0;%初始化误差表达式的分母 wucha=1;%初始化误差 error=1e-3;%压力迭代终止精度 error2=0.01;%承载力迭代终止精度 count=0;%初始化迭代次数 w=1.75;%松弛因子 while(wucha>=error) PP=P; for i=2:nx for j=2:ny P2(i,j)=P(i,j); A(i,j)=((H(i+1,j)+H(i,j))/2)/r(i,j); B(i,j)=((H(i-1,j)+H(i,j))/2)/r(i,j); C(i,j)=((deltax/deltay)^2)*((r(i,j+1)+r(i,j))/2)*((H(i,j+1)+H(i,j))/2)^3; D(i,j)=((deltax/deltay)^2)*((r(i,j-1)+r(i,j))/2)*((H(i,j-1)+H(i,j))/2)^3; E(i,j)=A(i,j)+B(i,j)+C(i,j)+D(i,j); F(i,j)=6*deltax*r(i,j)*((H(i+1,j)-H(i-1,j))/2); P3(i,j)=(A(i,j)*P(i+1,j)+B(i,j)*P(i-1,j)+C(i,j)*P(i,j+1)+D(i,j)*P(i,j-1)-F(i,j))/E(i,j); P(i,j)=(1-w)*P2(i,j)+w*P3(i,j); if P(i,j)<0 P(i,j)=0; end end end for i=2:nx for j=2:ny S=S+abs(P(i,j)-PP(i,j)); T=T+abs(P(i,j)); end end wucha=S/T; count=count+1; end load=0;%承载力初始值 for i=1:nx for j=1:ny ds=pi*(r(i,j+1)^2-r(i,j)^2)*(deltax/(2*pi)); load=load+k*P0*P(i,j)*ds*(r1^2); end end tt=abs(load-F1)/F1; if tt<error2 hmin=min(min(h)); hmax=max(max(h)); break; elseif load>F1 h1=(h1+h2)/2; h2=h1+l*tan(beita); else h2=(h1+h2)/2; h1=h2-l*tan(beita); end P0=omiga*miu*r1^2/(h1^2); end figure(1); =meshgrid((0:deltax/(pi/180):50),(0:deltay*r1:r2-r1)); mesh(x,y,h); xlabel('周向/°'); ylabel('径向r/m'); zlabel('h/m'); figure(2); =meshgrid((0:deltay*r1:r2-r1),(0:deltax/(pi/180):50)); mesh(x,y,P); xlabel('径向r/m'); ylabel('周向/°'); zlabel('P');</error2