帮我看看下面关于轴承力的程序错误在那里？

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

[ 本帖最后由 eight 于 2007-9-5 09:37 编辑 ]

花如月

这么长的程序，结果不对肯定是程序的问题。运行没问题不代表程序没问题，仔细再检查检查

eight

原帖由 花如月 于 2007-9-4 20:31 发表

                               
登录/注册后可看大图

这么长的程序，结果不对肯定是程序的问题。运行没问题不代表程序没问题，仔细再检查检查

恩，同意，除非有搞相同专业的高手路过，否则还是只能靠楼主自己解决。多检查，多调试一下

tianya7071

newmark法我看了好象没有问题，可能迭代求非线性轴承力出现问题了，就是不知道问题在那里，有matlab高手吗？帮我看看

xjzuo

请将你的问题出错情形,以及你想求什么讲一下.
另: 检查"下面判断语句"是否合适:
if Clearance(No_ball)<=0;
       Clearance(No_ball) = 0;
    end

xjzuo

还有"Clearance(No_ball,1)"这种用法是否不好?
以及调试、检查一下 Clearance(No_ball)是否>0 ？