根据黄平老师的《润滑数值计算方法》第四章内容，推导径向滑动轴承动压润滑计算方法，并用Matlab将程序复现，最后和书中的计算结果及其论文进行对比。
对于等温过程，径向滑动轴承动压润滑的油膜压力分布计算公式为：
（1）
公式中：R为轴颈半径，单位为m；p为油膜压力，单位为Pa；h为油膜厚度，单位为m；U为轴颈速度，单位为m/s，η为润滑介质黏度，单位为Pa*s。
将公式进行无量纲化，取：
（2）
式中，c为轴承的半径间隙，单位为m；L为轴承的宽度，单位为m。
将式（2）代入式（1）得：
利用差分法对式（3）进行求解：
求解的边界条件：
油膜的终点位置必须在求解过程中确定，是浮动边界条件，在每次迭代过程中，对于P<0的各节点
令P=0，最终可以自然地确定油膜终点位置。
根据论文《Hydrodynamic lubrication analysis of journal bearing considering misalignment caused by shaft deformation》中的轴承参数：
偏心率为0.8计算的油膜压力分布为：
最大油膜压力为33.071MPa，论文中计算结果为33.06MPa。
当考虑论文中的轴颈倾斜时，需要根据倾斜角改变轴承的油膜厚度，考虑轴颈倾斜时油膜厚度为：
具体的α和β代表的含义请详见论文。
当α=0.007，β=0时，油膜压力分布为：
最大油膜压力为62.53MPa，论文中计算的最大油膜压力为63.58MPa。
当α=0.02，β=90时，油膜压力分布为：
最大油膜压力为35.06MPa，论文中计算的最大油膜压力为34.95MPa。
当轴颈倾角较大时，油膜压力的分布较为集中，需要加密网格才能计算精确。
matlab代码为：
clear;
clc;
global CO B DY Dy
%% 初始参数
B=66E-3; % 轴承宽度 m
R=30E-3; % 轴承半径 m
% CR=0.00125; % 轴承间隙比
CO=0.03E-3; % 轴承半径间隙 m
R_shaft=R-CO; % 轴颈半径 m
AN=3000; % 转速r/min
EDA=0.009; % 润滑油动力黏度
EPSON=0.8; % 偏心率
%% 量纲一化
N=101; % 周向网格点数
M=41; % 轴向网格点数
DX=2*pi/(N-1); % 周向每两个网格之间的距离
DY=1/(M-1); % 轴向每两个网格之间的距离
Dy=B/(M-1);
OMEGA=AN*2.0*pi/60; % 转速r/min转化为角速度
U=OMEGA*R; % 角速度转化为线速度
ALFA=(R/B)^2; %
beta=(DX/DY)^2;
%% 轴颈倾斜参数
titling1=0.02*2*pi/360; % gama：轴颈在轴承中的倾斜角
titling2=90*2*pi/360; % alfa：OC2和C1C3之间的夹角
%%
H=oil_film(N,M,DX,EPSON,titling1,titling2);
[P,IK]=press(N,M,DX,EPSON,beta,ALFA,H);
p=6*P*U*EDA*R/CO^2;
Wx=0;
Wy=0;
for I=1:N
SETA=(I-1)*DX;
for J=1:M
Wx=Wx+p(I,J)*sin(SETA)*DX*Dy*R;
Wy=Wy+p(I,J)*cos(SETA)*DX*Dy*R;
end
end
W=sqrt(Wx^2+Wy^2);
phi=abs(atand(Wx/Wy));
%% 计算轴瓦歪斜时的力矩
Mx=0;
My=0;
[mesh_x,mesh_y]=size(p);
DX=2*pi/mesh_x;
Dy=B/(mesh_y-1);
for I=1:mesh_x
SETA=(I-1)*DX;
for J=1:mesh_y
distance_axial=Dy*(J-1)-B/2;
Mx=Mx+p(I,J)*cos(SETA)*DX*Dy*R*distance_axial;
My=My+p(I,J)*sin(SETA)*DX*Dy*R*distance_axial;
end
end
M_angle=atand(Mx/My);
M_new=sqrt(Mx^2+My^2);
%% 油膜厚度函数
function H=oil_film(N,M,DX,EPSON,titling1,titling2)
global CO B Dy DY
H=zeros(N,M);
oil_start_ps=0;
oil_end_ps=2*pi-oil_start_ps;
for I=1:N
SETA=(I-1)*DX;
for J=1:M
H(I,J)=1+EPSON*cos(SETA)+((J-1)*Dy-B/2)*tan(titling1)*cos(SETA-titling2)/CO;
end
end
return;
end
%% 压力分布函数
function [P,IK]=press(N,M,DX,EPSON,beta,ALFA,H)
T=ones(N,M);
P=zeros(N,M);
for I=1:N
for J=2:M-1
P(I,J)=0.5;
end
end
for J=1:M
P(1,J)=0;
P(N,J)=0;
end
for I=1:N
P(I,1)=0;
P(I,M)=0;
end
IK=0;
C1=inf;
while C1>1E-12
C1=0;
ALOAD=0.0;
for I=2:N-1
I1=I-1;
I2=I+1;
for J=2:M-1
PD=P(I,J);
J1=J-1;
J2=J+1;
A1=(0.5*(H(I2,J)+H(I,J)))^3;
A2=(0.5*(H(I,J)+H(I1,J)))^3;
A3=ALFA*beta*(0.5*(H(I,J2)+H(I,J)))^3;
A4=ALFA*beta*(0.5*(H(I,J)+H(I,J1)))^3;
P(I,J)=(-DX*(H(I2,J)-H(I1,J))/2+A1*P(I2,J)+A2*P(I1,J)+A3*P(I,J2)+A4*P(I,J1))/(A1+A2+A3+A4);
P(I,J)=0.7*PD+0.3*P(I,J);
if P(I,J)<0
P(I,J)=0;
end
C1=C1+abs(P(I,J)-PD);
ALOAD=ALOAD+P(I,J);
% continue;
end
end
IK=IK+1;
C1=C1/ALOAD
end
end