声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2441|回复: 6

[动力学和稳定性] 微分方程编程求解(有偿求助)

[复制链接]
发表于 2010-11-23 20:10 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
有一个项目,模拟球自由落体与地面发生碰撞的过程。方程很简单,就是一个二阶常微分方程。难点在于要时刻判断什么时候与地面接触,接触之后经弹性力作用发生反弹,之后再发生碰撞直到静止。需要matlab编程,且用变步长数值解法。请各位帮忙,有偿酬谢!
impact vibration.jpg
回复
分享到:

使用道具 举报

发表于 2010-11-23 20:53 | 显示全部楼层
该不是什么课的大作业吧……

点评

赞成: 2.0
赞成: 2
  发表于 2010-11-24 23:53
 楼主| 发表于 2010-11-23 21:00 | 显示全部楼层
回复 2 # Rainyboy 的帖子

请问,你会求解这类问题吗?

点评

不怎么会,MATLAB编程不太熟悉,来学习学习,呵呵  发表于 2010-11-23 21:23
发表于 2010-11-24 12:56 | 显示全部楼层
初步分析了下,这个问题并不复杂。没有碰撞的时候是自由落体运动,可以用大步长,出现碰撞时用小步长进行数值计算,判断条件题目已给。变步长的数值计算方法应该是比较成熟的了。自己找找看
分析问题的过程需要自己进行,中间有问题可以询问。整个问题放在这里,不加自己的分析和思考那是不太负责的

点评

赞成: 5.0
赞成: 5
  发表于 2010-11-24 23:54
 楼主| 发表于 2010-11-24 17:42 | 显示全部楼层
本帖最后由 greatchina 于 2010-11-24 17:44 编辑

回复 4 # hustxyon
你好,这是我已经编好的程序,可以解决问题,但是计算的步数太多了,希望可以修改一下步长。请多多指教,谢谢
clc;
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%碰撞接触刚度计算
r2=0.02;      %球半径;
r3=0.04;      %假设地面也为一球面
v2=0.3;       %泊松比
v3=0.3;       %泊松比
E2=207000000000;%弹性模量
E3=207000000000;%弹性模量
sg2=(1-v2^2)/E2;
sg3=(1-v3^2)/E3;
k=4/(3*(sg2+sg3))*(r2*r3/(-r2+r3))^0.5;%接触刚度。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=1.5;
e=0.9;            %e恢复系数;
m=10;             %球质量kg
g=9.8;            %重力加速度
h=1;              %自由落体初始高度m
vvcc=[];
q=[];
FF=[];
VVn=[];
dd=[];
y0=[h;0];          %初值(位置/速度)
tt=0:0.00001:4;    %计算时间,注意步长选取很小,否则计算的反弹力会很大,以至球反弹的高度比自由落体前的高度还要高,不符合实际(问题所在)。
for j=1:400000     %(按时间划分为100000个计算周期)
    d=y0(1)+r2;    %球最下端位置
    Vn=y0(2);      %接触相对速度
    if d<=0
        dd=[dd,d];        %记录球最下端位置
        vvcc=[vvcc,Vn];   %同时,记录球的速度
        vcc=vvcc(1);      %取出碰撞接触开始时刻速度;
        F=k*(abs(d))^n*(1+3*(1-e^2)/4*Vn/vcc);  %计算接触反弹力,与前期计算的刚度K和接触深度abs(d)及速度有关。   
        FF=[FF,F];        %记录反弹力的变化
    else
        F=0;        %没有发生接触,所以接触力为零
        vvcc=[];    %清空速度记录
    end
    Fn=F;           %取出反弹力,用于以下微分方程求解
    t0=tt(j);tf=tt(j+1);%开始计算
    tspan=linspace(t0,tf);
    [t,y]=ode45(@fangcheng,tspan,y0,[],m,Fn,g);
    y0=y(size(y,1),:); %y的最后一行所有列付给y0。也就是将本次计算结果,作为下次计算的初值。
    q=[q;y0];%收集结果
end
plot(tt(1:400000),q(:,1));%绘制球自由落体运动

function Dydt=fangcheng(t,y,m,Fn,g);
Dydt=[y(2);
      (-m*g+Fn)/m];

%问题描述:为了保证接触反弹力计算结果的正确性,计算中步长选取很小,导致计算时间很长。
%程序修改方向:在接触之前采用大步长计算,在即将发生接触之前,步长马上调整为小步长。(问题所在)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


1.jpg
发表于 2010-11-24 23:55 | 显示全部楼层
直接用ode45求解器求解是实现不了变步长的,有两个思路可选择
1、常微分方程变步长数值算法已经有很多成熟的算法,可以解决变步长的问题,你可以去查找一下相关资料。
2、程序里用while循环其实就可以用ode45实现两个状态用两种步长的功能,再想想看
发表于 2012-12-19 16:09 | 显示全部楼层
请问这个问题解决了吗?我也在做这方面的研究,急切需要高手指点,万分感谢!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 03:40 , Processed in 0.084929 second(s), 28 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表