main.m

运用New-Rapshon方法计算线接触常温弹流润滑程序

N=100;                            %节点数
P0=linspace(0,0,N+1);            %初始压力P0向量
P=linspace(0,0,N+1);             %压力P向量
DP=linspace(0,0,N+1);            %ΔP
XP=linspace(0,0,N+1);            %各个节点值X
H=linspace(0,0,N+1);             %膜厚H
F=linspace(0,0,N+1);             %F向量
B=linspace(0,0,N+1);             %B向量
DF=linspace(0,0,N+1);            %ΔF
DL=linspace(0,0,N+1);            %求CIJ时的一个中间变量
MIDU=linspace(0,0,N+1);          %密度压力关系
NIANDU=linspace(0,0,N+1);        %粘度压力关系
CIJ=zeros(N+1,N+1);              %C矩阵
KK=zeros(N,N);               %最终的K矩阵N*N
KK1=zeros(N+1,N+1);              %最终的K1矩阵
KK2=zeros(N+1,N+1);              %最终的K2矩阵
LNODS=zeros(N+1,2);              %中间矩阵

initialization;         %初始化
for i=1:N+1
    if(XP(i)^2<=1)
        P0(i)=sqrt(1-XP(i)^2);
    end
end
P=P0;
DH0=1;
MAXIT=100;
ITER=0;
while (abs(DH0)>0.01&&ITER<=MAXIT)
    ITER=ITER+1;
    IT=0;
    EP=1;
    while (EP>0.01&&IT<=MAXIT)
        IT=IT+1;
        film_thick;          %求膜厚方程
        matrix_final;        %求解最终矩阵
       DP(2:N+1)=KK\DF(2:N+1)'; %求解方程
       DH0=DP(N+1);
       DP(N+1)=0;
       sum2=0;
       sum3=0;
       for i=1:N
           sum2=sum2+abs(DP(i));
           sum3=sum3+P(i);
       end
       EP=sum2/sum3; 
       P=P+0.4*DP;    %下山法，保证函数的绝对值稳定下降，加快收敛速度
   end
    H0=H0+0.3*DH0;
end
%dangliangthick;        %求HOIL

figure(1);
plot(H,'r-');    %划出膜厚形状曲线
hold on;
plot(P,'b-');    %划出压力分布曲线
hold on;
title('润滑膜形状和压力分布');
axis([1 N 0 1.5]);
set(gca,'Xtick',[1:5:100],'Ytick',[0:0.2:5]);     %设置X轴的坐标从0到100，间距5;
                                                    %Y轴坐标从0到5，间距0.2

grid on;                                            %画坐标分隔线
xlabel('Ｘ坐标点');
ylabel('润滑膜厚/压力值');
legend('润滑膜厚H','压力值P',-1);
hold off;