无水1324 发表于 2008-12-18 22:27

数值分析、还要非线性动力学方面的都要学习一下

cam_1980 发表于 2008-12-19 15:05

回复 楼主 无水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 编辑 ]

ykczj 发表于 2008-12-22 13:10

楼上这位兄弟能否把此问题的方程或参考文献给出,这样大家一起讨论更方便些

cam_1980 发表于 2008-12-22 14:20

回复 18楼 ykczj 的帖子

方程是最普通的方程,如下:

将其先写成范式,然后就可以用rouge-kutta法编程了。这里的方程是为了举个例子随便弄的一个方程。程序中的f就是上面方程中的f(x)计算式,如果是搞齿轮动力学的一看就明白了。

[ 本帖最后由 cam_1980 于 2008-12-22 14:23 编辑 ]
页: 1 [2]
查看完整版本: 碰撞求解