碰撞求解

无水1324

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

cam_1980

无水，不知道你这个问题解决没有，当时我在处理间隙的时候，发现用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=[xt,yt];
程序里的计算k2n的都改成对应你的方程即可。
如果你考虑了有碰撞恢复系数的那种碰撞的话，那就要去掉程序里的f计算，而另外加有碰撞系数的数值计算。

[ 本帖最后由 cam_1980 于 2008-12-19 15:10 编辑 ]

ykczj

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

cam_1980

方程是最普通的方程，如下：

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

[ 本帖最后由 cam_1980 于 2008-12-22 14:23 编辑 ]