帮我看看下面的计算程序出错在那里？

tianya7071

程序运行没有问题，就是轴承力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

tianya7071

查看了好几天也没有弄出是什么原因来，着急啊，请大牛指点一下

咕噜噜

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;% 滚珠与内滚道的间隙变化。
这里有问题，sita应该是sita(i),类型要一致，你后面时间t为循环

tianya7071

有没有研究轴承非线性的
帮忙看看啊

咕噜噜

^_^，你要问什么问题，可以具体说，轴承的非线性很多方面和其他系统有很多相似性

无水1324

调试了很长时间还是没有找出错误的地方，后面的高手继续

tianya7071

我用ode方法求得的fq的结果还可以，可是用NewMark求的fq却为零，不知道是什么原因，有这方面的的高手请指点
下面是ode方法的程序
function dq = BallBrg_NonL_Sub_lyl(t,q)
%子程序

% 求解外圈固定球轴承的变柔度(VC)振动
global w1 w2 Ri Re Nb gama kn  w2_rpm
% 为了方便绘制Poincare而设置的参数
n_One_T = 100;% 每个周期的采样点数
n_T = 100;% 采样时间占几个周期
%  滚子 球轴承参数
Ri=0.0585;% 内滚道直径
Re=0.0665;% 外滚道直径
Nb=34;     % 滚子数
F = 1000;% 径向力(N)
q_initial(1:4,1) = 1e-11;% 初始值
gama = 0.00002;% 间隙(m)
kn =6.2e8;
  w_cage = ( w2*Ri+w1*Re )/(Re+Ri);% ;% 保持架转速（rad/s）
  M=11.63*[1 0; 0 1];% 质量矩阵
    C=1200*[1 0; 0 1];% 阻尼矩阵
fq=zeros(2,1);% 轴承力初值
diff_1_3 = q(1,1);% 水平方向位移
diff_2_4 = q(2,1);% 垂直方向位移
% 求轴承的非线性反力
for No_ball=1:Nb
    sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
    Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
        + diff_2_4*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
Fq(1,1)= - fq(1,1) % 水平方向外力
Fq(2,1)= - fq(2,1) + F;% 垂直方向外力
% 动力学微分方程
dq(3:4,1)=inv(M)*(Fq-C*q(3:4,1));% x和y方向加速度
dq(1:2,1)=q(3:4,1);

function BallBrg_NonL_Forum
%主程序
% 求解外圈固定球轴承的变柔度(VC)振动
clear
clc
%% 参数设置
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;     % 滚子数
F = 1000;% 径向力(N)
n_n = 0;
w2_limit1=6000;% 内圈最低转速（rpm）
w2_limit2=10400;% 内圈最高转速（rpm）
w_step = 100;% 转速变化步长(rpm)
q_initial(1:4,1) = 1e-11;% 初始值
gama = 0.00004;% 间隙(m)
kn =6.2e8;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.1)。
    w1 = 0;% 外圈角速度
     M=11.63*[1 0; 0 1];% 质量矩阵
    C=1200*[1 0; 0 1];% 阻尼矩阵
   
    %% 响应计算循环
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/100;% 取点时间步长，dt随转速变化。
    time=n_T*T_vc;% 总的时间
    n = round(time/dt);% 离散点数
    t_span(1:n) = linspace(0,time,n);% 时间数组
    [t,q]= ode45('BallBrg_NonL_Sub_lyl', t_span, q_initial);
end