本人在将黄平老师所著的润滑数值计算方法的第二章不可压稳态线接触滑块程序改编成MATLAB程序时压力的计算出了问题,但是自己找不到解决办法,请大佬解惑,谢谢。程序如下:
X = zeros(121,1); % 坐标节点数组
H = zeros(121,1); % 膜厚节点数组
P = zeros(121,1); % 压力节点数组
H1=1.0; % 量纲化最大膜厚
H2=0.5; % 量纲化最小膜厚
U=1; % 量纲化速度
EDA=0.02;% 润滑油粘度
AL=0.01; % 滑块长度
KG=1;% 等于1时为直线滑块,2时为曲线滑块
N=121;% 节点取值
if KG==2
X1=-1.0;
X2=1.0;
else
X1=0.0;
X2=1.0;
end
DX=1/(N-1);
for i=1:N
if KG==1 % 直线滑块膜厚计算
X(i)=X1-(i-1.0)*DX*(X1-X2);
H(i)=H1/H2-(H1/H2-1.0)*X(i);
else % 曲线滑块膜厚计算
X(i)=X1-(i-1.0)*DX*(X1-X2);
H(i)=1.0+(H1/H2-1.0)*X(i)^2;
end
end
for i=2:N-1
P(i)=0.5;% 压力各点赋初值
end
P(1)=0.0;% 压力起点初始化
P(N)=0.0;% 压力终点初始化
IK=0;
ERO=1;
while ERO>1e-7
C1=0.0;
ALOAD=0.0;
for i=2:N-1
A1=(0.5*(H(i+1)+H(i)))^3;
A2=(0.5*(H(i)+H(i-1)))^3;
PD=P(i);
P(i)=(-0.5*DX*(H(i+1)-H(i-1))+A1*P(i+1)+A2*P(i-1))/(A1+A2);
P2(i)=0.3*PD+0.7*P(i);
if P2(i)<0
P2(i)=0;
end
C1=C1+abs(P2(i)-PD);
ALOAD=ALOAD+P2(i);
end
ERO=C1/ALOAD;% 相对误差,用于判断收敛情况
IK=IK+1;
end
X0=0;
for i=1:N
X0=X0+P(i)*X(i);
end
X0=X0*AL;
ALOAD=ALOAD*DX*AL*6.0*U*EDA*AL/(H2^2);% 承载能力计算
if KG==1
for i=1:N
P0=-(-1.0/(H(i)*H2)+H1*H2/(H1+H2)/(H2*H(i))^2+1.0/(H1+H2))/(H1/H2-1.0)*H2;
end
end @wuming524