matlab_雷诺方程_有限差分法_数值计算_

% 雷诺方程求解过程----不可压缩流体密度rou=常数，v=0，膜厚h不变 clc,clear deltax=0.01*10^(-3); % deltax为板长方向的变化量0.01mm h1=5*10^(-6)*ones(1,501); % h1为初始粉末的膜厚，1um,3um,5um,10um h2=1*10^(-6)*(0+1.3*randn(1,501)); % h2为均值是0，方差为1.3的微凸体高度---可能有问题 h=h1-h2; % h为实际的膜厚 for i=1:1:501 if h(1,i)<0 h(1,i)=0; % h(1,i)为任意点处的压力，小于0时取为0 end end U=2; % U为上下板的相对滑移速度，初始令为2m/s，4m/s，6m/s B=0:0.01:5; % 求解雷诺方程过程 p=ones(1,501); % p为压力，初始赋值为1行501列的数组 p(1,1)=0; % p(1,1)为初始起始点的类流体压力，由边界条件为0 p(1,501)=0; % p(1,501)为终点的类流体压力，由边界条件为0 pp=ones(1,501); % pp为前一次的压力，初始赋值为1行500列的数组 eta0=0.014; % eta0为初始粘度，令为0.014 eta=eta0*ones(1,501); % eta为初始粘度，令为0.014 alpha=2.2*10^(-8); % alpha为流体的粘压系数，令为2.2*10^(-8) c=alpha/15; % c为流体的粘压系数，令为2.2*10^(-8)/15 z=0.68; p0=5.1*10^(-9); % 参数a1、a3、a4、a0共四个参数 pp=0; % pp为初始总压力，令为0 error=10^(-6); % error为理论的误差，初始值令为10^(-6) wucha=1; % wucha为实际的误差，初始值令为1 count1=0; while ( wucha>error) % wucha>error作为判断条件 pp=p; % pp为前一次的类流体压力值,为后面计算误差所用 pS=0; % pS为这一次和上一次求得整个板的类流体压力总和之差 pT=0; % pT为整个板的类流体压力总和 for i=2:1:500 a1=1/(2*deltax^2).*((h(1,i).^3/eta(1,i))+(h(1,i-1).^3/eta(1,i-1))); a3=-1/(2*deltax^2).*((h(1,i+1).^3/eta(1,i+1))+(h(1,i).^3/eta(1,i))+(h(1,i-1).^3/eta(1,i-1))); a4=1/(2*deltax^2).*((h(1,i+1).^3/eta(1,i+1))+(h(1,i).^3/eta(1,i))); a0=3*U/deltax.*(h(1,i+1)-h(1,i-1)); p(1,i)=(a0-(a1.*p(1,i-1)+a4.*p(1,i+1)))./a3; % p为任意点处压力 eta(1,i)=eta0*exp(alpha*p(1,i)); % Barus粘压公式 pp=p(i)+pp; % pp为总压力 if p(1,i)<0 % p(1,i)为任意点处的压力，不可能小于0, 取反或者令为0 p(1,i)=p(1,i-1); end end for i=1:1:500 pS=pS+p(i)-pp(i); % pS为这一次和上一次求得整个板的类流体压力总和之差 pT=pT+p(i); % pT为整个板的类流体压力总和 end wucha=abs(pS/pT); count1=count1+1; end h p count1 wucha figure(1) plot(B,h,'r-') % 板长方向的膜厚分布图 figure(2) plot(B,p,'r-') % 板长方向的压力分布图