NEWTON-RAPSION线接触 弹流润滑数值解法

H(i)=H(i)+DH; %膜厚 %
end %
子程序 matrix_element.m
if aa==bb
CIJ(i,j)=(-0.5/3.14)*DL(j-1)*log(0.01);
else
CIJ(i,j)=(-0.5/3.14)*DL(j-1)*log((aa-bb)^2);
A=zeros(2,2); %K1 % %
AL=zeros(2,N+1); %K2 % %
DL(j-1)=XP(j)-XP(j-1); %
if(i==j||i==j-1)
aa=XP(i);
bb=XP(j)+0.1*DL(j-1);
sum2=sum2+abs(DP(i));
sum3=sum3+P(i);
end
EP=sum2/sum3;
P=P+0.4*DP; %下山法，保证函数的绝对值稳定下降，加快收敛速度
end
H0=H0+0.3*DH0;
AL(2,i)=-AL(1,i); % %
BB=linspace(0,0,2); % %
FF=linspace(0,0,2); % %
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矩阵
H(i)=H0+XP(i)^2/2; %膜厚公式前部分 %
DH=0; %
for j=2:N+1 %
T=3.0*(pi^2)*U/(4*W^2); %雷诺方程右项系数 %
H0=0; %初始膜厚 %
XP(1)=-2.5; %
xlabel('Ｘ坐标点');
ylabel('润滑膜厚/压力值');
legend('润滑膜厚H','压力值P',-1);
hold off;
子程序 initialization.m
NIANDU0=0.08; % 初始粘度η0
U=1.0e-11; %速度参数 %
XP(i)=XP(1)+(i-1)*D; %各节点值 %
end
子程序 film_thick.m
for i=1:N+1 %
主程序 main
N=100; %节点数
P0=linspace(0,0,N+1); %初始压力P0向量
P=linspace(0,0,N+1); %压力P向量
DP=linspace(0,0,N+1); %ΔP
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
Hale Waihona Puke P=P0; DH0=1;
for i=1:N %K2 % %
AL(1,i)=T*(MIDU(NEL-1)*CIJ(NEL-1,i)+MIDU(NEL)*CIJ(NEL,i))/2; % %
end
%dangliangthick; %求HOIL
figure(1);
plot(H,'r-'); %划出膜厚形状曲线
hold on;
plot(P,'b-'); %划出压力分布曲线
hold on;
title('润滑膜形状和压力分布');
axis([1 N 0 1.5]);
end
else %
CIJ(i,j)=(-0.5/3.14)*DL(j-1)*log((XP(i)-XP(j))^2); %
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); %粘度压力关系
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向量
E=linspace(0,0,N+1); % %
for i=NEL-1:NEL % %
XP(N+1)=1; %
D=(XP(N+1)-XP(1))/N; %步长 %
for i= 2:N %
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
end %
DH=DH+CIJ(i,j)*P(j); %膜厚公式后部分δ(X) %
end %
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; %求膜厚方程
R=0.05; %当量圆柱半径 %
W=2.0e-5; %单位长度载荷（无量纲化） %
G=5000; %材料参数（无量纲化） % %
A0=2.2e-8
E0=G/A0; %当量弹性模量 %
B=(sqrt(8.0*W/pi))*R; %接触区半宽 %
PH=E0*B/(4*R); %最大接触应力 %
E(i)=MIDU(i)*H(i)^3/NIANDU(i); % %
end % %
set(gca,'Xtick',[1:5:100],'Ytick',[0:0.2:5]); %设置X轴的坐标从0到100，间距5;
%Y轴坐标从0到5，间距0.2
grid on; %画坐标分隔线