回复 楼主 无水1324 的帖子
无水，不知道你这个问题解决没有，当时我在处理间隙的时候，发现用ode计算结果有误，于是就干脆用显式rouge-kutta法自己编了一个，结果就没有问题了（与别人已做的工作进行比较知道结果正确了）。你说的这个碰撞振动问题我以后肯定也会碰到的，因为我现在还没有考虑那种有恢复系数的碰撞，只是用了大家都熟悉的间隙函数f（x）来处理碰撞。我的函数如下（希望对你还有帮助）：tt是向量，x和y可以是标量也可以是向量（主要是用伪不动点追踪算法里用向量），函数返回值是tt末了的时间的相点。
function yyy=shuzhidy(tt,x,y)
global jieta; %阻尼比
global k0; %平均刚度
global k1; %刚度振幅
global w; %激励频率
global bb; %量纲一化齿侧间隙
global P0; %平均激励力
global P1; %激励力振幅
global period; %外周期项
global h %积分时间步长
xx=x;
yy=y;
n=max(size(tt));
for i=1:1:n-1
t=tt(i);
k11=yy;
f=(xx>bb).*(xx-bb)+(xx<-bb).*(xx+bb); %直接用逻辑式参与计算，这样可以计算x和y是向量的情况
k21=-2*jieta*yy-(k0+k1*cos(w*t))*f+P0+P1*cos(w*t);
k12=yy+h*k21/2;
f=(xx+h*k11/2>bb).*(xx+h*k11/2-bb)+(xx+h*k11/2<-bb).*(xx+h*k11/2+bb);
k22=-2*jieta*(yy+h*k21/2)-(k0+k1*cos(w*(t+h/2)))*f+P0+P1*cos(w*(t+h/2));
k13=yy+h*k22/2;
f=(xx+h*k12/2>bb).*(xx+h*k12/2-bb)+(xx+h*k12/2<-bb).*(xx+h*k12/2+bb);
k23=-2*jieta*(yy+h*k22/2)-(k0+k1*cos(w*(t+h/2)))*f+P0+P1*cos(w*(t+h/2));
k14=yy+h*k23;
f=(xx+h*k13>bb).*(xx+h*k13-bb)+(xx+h*k13<-bb).*(xx+h*k13+bb);
k24=-2*jieta*(yy+h*k23)-(k0+k1*cos(w*(t+h)))*f+P0+P1*cos(w*(t+h));
xt=h/6*(k11+2*k12+2*k13+k14)+xx;
yt=h/6*(k21+2*k22+2*k23+k24)+yy;
xx=xt;
yy=yt;
end
yyy=;
程序里的计算k2n的都改成对应你的方程即可。
如果你考虑了有碰撞恢复系数的那种碰撞的话，那就要去掉程序里的f计算，而另外加有碰撞系数的数值计算。
[ 本帖最后由 cam_1980 于 2008-12-19 15:10 编辑 ] 楼上这位兄弟能否把此问题的方程或参考文献给出，这样大家一起讨论更方便些
回复 18楼 ykczj 的帖子
方程是最普通的方程，如下：将其先写成范式，然后就可以用rouge-kutta法编程了。这里的方程是为了举个例子随便弄的一个方程。程序中的f就是上面方程中的f（x）计算式，如果是搞齿轮动力学的一看就明白了。
[ 本帖最后由 cam_1980 于 2008-12-22 14:23 编辑 ]
页:
1
[2]