|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
程序运行没有问题,就是轴承力fq迭代的结果一直是零,请高手帮我看看
% NewMark Function
clc,clear
%% 参数设置
global w1 w2 Ri Re Nb gama kn w2_rpm
n_One_T = 100;% 每个周期的采样点数
n_T = 100;% 采样时间占几个周期
% 滚子 球轴承参数
Ri=0.0585;% 内滚道直径
Re=0.0665;% 外滚道直径
Nb=34; % 滚子数
n_n = 0;
w2_limit1=6000;% 内圈最低转速(rpm)
w2_limit2=10000;% 内圈最高转速(rpm)
w_step = 100;% 转速变化步长(rpm)
gama = 0.00002;% 间隙(m)
kn =6.2e8;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.1)。
w1 = 0;% 外圈角速度
M=11.63*[1 0; 0 1];% 质量矩阵
C=1200*[1 0; 0 1];% 阻尼矩阵
F =1000*[0;1];% 径向力(N)
%% 响应计算循环
for w2_rpm=w2_limit1:w_step:w2_limit2 %内圈转速变化
n_n = n_n+1 % 计数变量
disp(w2_rpm)
w22 = w2_rpm*pi/30;% 转化为rad/s单位
w2 = w22;% 内圈角速度
w_cage = ( w2*Ri+w1*Re )/(Re+Ri);% 保持架
w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
T_vc = 1/w_vc;% vc周期
dt=T_vc/n_One_T;% 取点时间步长,dt随转速变化。
time=n_T*T_vc;% 总的时间
u=[0 0].';
v=[0 0]';
a=[0 0]';
q(:,1)=u; %distance
q1(:,1)=v; %speed
q2(:,1)=a; %accelerrate
%newmark kcoefficient
gama=0.5;
dt=T_vc/n_One_T;
delta=0.25;
b0=1/(delta*dt^2);
b1=gama/(delta*dt);
b2=1/(delta*dt);
b3=1/(2*delta)-1;
b4=gama/delta-1;
b5=dt*(gama/(2*delta)-1);
b6=dt*(1-gama);
b7=gama*dt;
%等效刚度矩阵
K1=b0*M+b1*C;
t_max=time; %计算时间总长
i=1;
t(1)=0; %time
QQ(:,1)=zeros(2,1);
while t(i)<t_max
fq=zeros(2,1);% 轴承力初值
% 求轴承的非线性反力
for No_ball=1:Nb
sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t(i);% 第No_ball个滚珠的位置角
Clearance(No_ball,1) = q(1,i)*sin( sita(No_ball) ) ...
+ q(2,i)*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
% 判断哪几个滚动体受到接触力
if Clearance(No_ball)<=0;
Clearance(No_ball) = 0;
end
fs = abs( (Clearance(No_ball))^1.1 );
fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball))
fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
end
Q=[fq(1,1);fq(2,1)];
Q(:,i)=[0 0]';
QQ(:,i+1)=F-Q(:,i)+M*(b0*q(:,i)+b2*q1(:,i)+b3*q2(:,i))+C*(b1*q(:,i)+b4*q1(:,i)+b5*q2(:,i));
q(:,i+1)=inv(K1)*QQ(:,i+1);
q2(:,i+1)=b0*(q(:,i+1)-q(:,i))-b2*q1(:,i)-b3*q2(:,i);
q1(:,i+1)=b1*(q(:,i+1)-q(:,i))-b4*q1(:,i)-b5*q2(:,i);
i=i+1;
t(i)=t(i-1)+dt;
end
end
[ 本帖最后由 eight 于 2007-9-5 09:37 编辑 ] |
|