杭州锐达数字技术有限公司
查看: 9093|回复: 52

[Python] 瞬态动力学方程组的求解程序In Python

  [复制链接]
发表于 2011-3-12 13:20 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 Rainyboy 于 2011-3-26 22:53 编辑

最近在重新阅读瞬态动力学计算方面的书籍,有一些新的体会,见:Euler法求解瞬态动力学方程的稳定性

在这个帖子中,打算进一步给出一些常见算法的实现,希望有兴趣的同志一起来测试、改进、完善。

==============瞬态动力学计算要干什么===============

工程上较复杂的振动问题,经适当的离散化后,往往可以用多自由度系统的振动理论来解决。描述其动力学行为的通常是一组自由度之间相互耦合的常微分方程,可以通过牛顿第二定理、拉格朗日方程和影响系数法等途径建立,多自由度系统的动力学方程组形如:
QQ截图未命名.jpg
其中: [M]为系统的质量矩阵,[C]为系统的阻尼矩阵, [K]为系统的刚度矩阵,{F} 为各自由度上的瞬态激振力。

瞬态动力学计算就是要用数值方法求解这个方程。通常可分为直接积分法和振型叠加法。

这个帖子给出的代码是我自己折腾这些算法时写的,至2011.3.12,搭了一个简陋的框架,完成了欧拉法和中心差分法的编码,经过了简单的测试。有兴趣想折腾折腾的同志可以拿去用,不过折腾完了希望把你的改进再发回来,大家一起折腾,呵呵!

==============程序设计:BaseClass.py的结构============

Mecha_Exception类
继承的Exception,用以表示在动力学计算中的异常,与其他异常相区分


MechanicMode类
用于储存M,K,C矩阵、表示F(t)关系的函数对象,以及系统的初始条件。其中,M、K、C为numpy.matrix对象,F为接受一个时间参数t的函数对象,初始条件也为numpy.matrix对象。在该类中,实现对动力学模型的储存和检验,目前主要做的是维数匹配。

MechanicSolutionMethod类:
用于作为动力学求解方法的基类。因此其提供一些每个求解方法都有的公用函数。
  .setMode()函数用于将模型绑定到求解方法上
  .solve()模型用于表示求解过程,求解绑定到该对象上的动力学模型
  .DrawTimeHisProc()在求解后,绘出指定自由度的的速度、加速度和位移曲线
  .ShowAllPara()是调试的时候用的,输出当前的所有参数,进行检查


==============程序设计:EularMethod.py的结构=============
EularMethod类:继承MechanicSolutionMethod,按照欧拉方法重写其solve()函数


==============程序设计:CenterDiffMethod.py的结构=============
CenterDiffMethod类:继承MechanicSolutionMethod,按照中心差分法重写其solve()函数


===============程序设计:TestCase.py的结构====================
该文件主要用于构造算例,对各种方法进行调用。
Get_A_1DOFMode()函数:顾名思义,返回一个单自由度MechanicMode对象
Get_A_2DOFsMode()函数:顾名思义,返回二个单自由度MechanicMode对象
Run()函数:基本计算流程:绑定,求解,后处理。


=================已经做的一些测试================================
A,欧拉方法在计算单自由度、二自由度模型上的正确性;B,中心差分法在计算单自由度、二自由度模型上的正确性;

例1:中心差分法计算的单自由度系统瞬态响应,步长0.3秒,结果未发散,结果的数值与理论分析吻合
mid.png
例2:欧拉法计算的单自由度瞬态响应,步长0.01,发散,与理论分析吻合
ul.png
例3:中心差分法计算二自由度系统的(某自由度的)响应,与《计算动力学》的结果一致
3.png

=================个人觉得可以继续做的地方==========================
A. 实现从文件中读入动力学模型(方便导入ANSYS等软件建立的模型)
B1. 实现更多的动力学计算方法
B2. 寻求更多的测试算例
C. 考虑如何引入非线性因素

说明:我的最终目的是通过编写程序达到对瞬态动力学计算更深的体会,上述ABC也是我个人的练习打算,要是有哪位同志也希望练习练习,不妨告知我,大家一起来。
==================2011.3.26=========================
新增Newmark方法,见附件。
测试算例1:用NEWMARK方法计算单自由度系统的响应
1.png


测试算例2:用NEWMARK方法计算二自由度系统的响应
2.png

=====================================================

Dynamics Practice.rar

3.42 KB, 阅读权限: 20, 下载次数: 25

全部代码

TestCases.rar

860 Bytes, 下载次数: 41

包含NewMark法的测试算例

NewmarkMethod.rar

876 Bytes, 下载次数: 54

NewMark法

点评

很不错  发表于 2011-3-15 08:23

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2011-3-15 08:26 | 显示全部楼层
你太有激情了,我先看看代码哈,另外,这样的方程最困难的地方在你说的瞬时上,不知道 你说的方程右端的瞬时激发力有什么特征。
看你写的这个框架 真的很不错,我很想加入。这段时间正好时间上比较充足!呵呵。
发表于 2011-3-15 08:37 | 显示全部楼层
跟踪你的贴子,还是没有看到什么叫瞬态动力学方程组,更幽默的是在百度上查这个词,出现的第一个是我现在回复的帖子。劳烦简单介绍一下,或者给点小资料。
 楼主| 发表于 2011-3-15 10:42 | 显示全部楼层
回复 2 # smtmobly 的帖子

瞬态激振力没有什么固定的,一般是简谐的,也有可能是方波冲击,也有可能是白噪声什么的,在建立方程的时候没有对它作约定……

有你python大牛&未来数学家的加入,肯定能产生很多有意思的结果啦。

我弄这个的目的是为了通过编码把现在瞬态动力学方程的数值方法中的一些考虑弄透彻(主要是参考清华大学张雄的《计算动力学》),当然也是希望能为以后的研究工作所用啦,毕竟这套方程组的求解是很多结构动力学问题的基础。

也许你的出发点与我有差别,如果有你感兴趣的东西的话,代码你拿去随便用随便改好了,只是希望改完之后还放回来,咱们一起看看能不能将各自的工作统一在一个框架下,方便交流和使用,呵呵。
发表于 2011-3-15 20:51 | 显示全部楼层
你咋这样说我呢,我都不知道自己是谁了。一起玩嘛,有兴趣就好!
发表于 2011-3-16 21:46 | 显示全部楼层
这个很有意思,
请问 Rainyboy 在什么环境下编译的?
python没用过,不知该如何下手
 楼主| 发表于 2011-3-17 11:05 | 显示全部楼层
回复 6 # wei_x 的帖子

这个是用python写的,所以要先装python解释器,python是一种解释型语言……不需要编译器……

你可以看看这个帖子,是我当时刚弄python的时候写的,里面说了在哪里下载解释器之类的:

Python:一场从零开始的奇妙旅程
http://forum.vibunion.com/forum- ... fromuid-159019.html
发表于 2011-3-18 14:34 | 显示全部楼层
恩,有意思,跟踪学习以后
 楼主| 发表于 2011-3-26 22:55 | 显示全部楼层
2011.3.26:新增Newmark方法,代码和测试用例见主楼。
发表于 2011-3-31 13:17 | 显示全部楼层
回复 9 # Rainyboy 的帖子

楼上的能不能把我带入python world

看了一下,有要好多选择
在windows下,推荐一下吧
 楼主| 发表于 2011-3-31 14:10 | 显示全部楼层
回复 10 # tenglang 的帖子

http://www.pythonxy.com/
下载pythonxy套件吧,呵呵!
发表于 2011-8-31 11:39 | 显示全部楼层
希望得到Rainyboy的指教
发表于 2011-9-19 12:42 | 显示全部楼层
回复 12 # fawcgzmg 的帖子

非常不错,这正是我最近研究齿轮的动力学方程中所遇到的问题,请问一下,如果不用PYTHON用别的编程语言能做吗,PY的程序界面做起来实在麻烦。
 楼主| 发表于 2011-9-19 13:25 | 显示全部楼层
回复 13 # fawcgzmg 的帖子

python的guidata库做起界面来很方便。见:

guidata库使用手记:快速定制科学计算程序的图形界面
http://forum.vibunion.com/forum-viewthread-tid-105433-fromuid-159019.html
发表于 2011-9-20 11:50 | 显示全部楼层
回复 1 # Rainyboy 的帖子

有QQ号码或邮箱吗?我想把日文的这篇文章发给你,我算的固有频率对不上,这个方程也不会求解
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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