本博客主要是来讲述采用Python语言,利用直接刚度法,来编写求解《有限元方法基础教程》(第五版) 的例题2.1。
(也许有的人认为商业软件已经和普及了,没有 必要再自己进行编程了,但是龙猪在这里建议大家还是去了解一下有限元的原理,这样对我们深入了解软件的原理和进一步的应用会有很大的帮助)
可参考该作者的文章:《有限元分析基础教程》(曾攀)笔记一-二维杆单元有限元程序(基于Python) - SimuLife - 博客园 https://www.cnblogs.com/SimuLife/p/4695922.html。
对于本问题的具体工作流程,同样是采用了上述作者的编码习惯。首先创建了一个弹簧单元类:
注意这里的22行,因为本问题是一维的,坐标选取的是x轴,那么相应地Y轴上的坐标始终保持为0。于是便采用了横坐标相减来计算出单元的长度。这里的start_point和end_point建议采用tuple类型,每个tuple类型里头有两个元素,分别是对应的x, y坐标。
i, j 分别是该单元的起始节点索引和终止节点索引。
在整个直接刚度法的过程中,最重要的是单元刚度矩阵的组合:
对这里的参数进行一个解释,这里的kk与SimuLife - 博客园中所述的是一致的,kk表示组合刚度矩阵,这里的ElementSet是单元集合,(这里说明一下,如果您在Abaqus中使用过Python来读取odb文件的话,会发现所有的单元其实是保存在一个单元集合中的,还有节点集合等等。)==》可参考: 《ABAQUS Python 二次开发攻略》
# 然后对总体刚度矩阵进行一个组合。
对于当前的组合,首先需要判断单元刚度矩阵在总体刚度矩阵中的位置。这个主要靠的就是该单元的起始节点索引和终止节点的索引。具体的原理可查看《有限元方法基础教程》。
# 创建一个计算节点位置的函数
# 本问题的边界条件是按照书本上的要求输入的,对于未知变量同样采用的字符"u_unknown"和“f_unknown”来代替。
这里我们注意看这个dtype函数,该dtype是表示的numpy中array的数据类型。这里一开始,我没有设置这个,采用的是默认的值。然后在通过np.where来查找u=0 和 u="u_unknown"时,会出现异常,异常如下所示:
这个异常的意思好像是在对比Python数据元素和Numpy中的元素时,Python开发者和Numpy开发者出现了分歧,Python开发者坚持采用Pythonic, 但是Numpy开发者坚持采用Numpy的数据表示方式。所以这个问题被暂时搁置了。虽然采用了dtype=object可以暂时解决这个问题,但是在程序的后面,就是通过刚度矩阵和位移数组来计算出力数组的时候,会因为数据类型的原因,抛出异常。
该异常时因为np.linalg.solve引发的,所以才会有60行那个数据类型转换的操作.
# 接下来是主函数
运行结果如下:
有什么问题可以留言。
赞助一分钱给老猪。