python有限元例题_《有限元分析基础教程》(曾攀)笔记一-二维杆单元有限元程序(基于Python)...
曾攀老师的《有限元分析基础教程》第三章有二维杆单元的推导,并结合一个例题进行了解析解和基于Matlab的程序求解。但是我感觉书中的MATLAB代码有点罗嗦,而且一些实现方法也比较麻烦,比如已经知道了杆单元的起点和终点坐标,仍然需要另外给出单元局部坐标与整体坐标的夹角,这完全没必要。于是我就用Python重构了这段程序,当然并不是把书中的MATLAB代码翻译成python(事实上完全可以这么干,而且很快!)。比如我使用了面向对象的思想,把杆单元写成了一个类,这样思路比较清晰。
#! /usr/bin/python#-*- coding: utf-8 -*-
importmathimportnumpy as np
sqrt, cos, sin, pi=math.sqrt, math.cos, math.sin, math.pi"前处理"nodeNumber= 4KK= np.zeros((2*nodeNumber, 2*nodeNumber))"""边界条件,U表示节点的位移向量,如果某个自由度的位移未知,则该处填写‘u_Unknown’,
F表示节点的荷载向量,如果某个自由度的荷载未知,则该处填写‘f_Unknown’"""U= np.array([0, 0, ‘u_Unknown‘, 0, ‘u_Unknown‘, ‘u_Unknown‘, 0, 0], dtype=object)
F= np.array([‘f_Unknown‘, ‘f_Unknown‘, 2e4, ‘f_Unknown‘, 0, -2.5e4, ‘f_Unknown‘, ‘f_Unknown‘], dtype=object)classBar2D:"""定义二维杆单元类,该类包含杆件的基本信息:
E 弹性模量,A 杆单元面积,i 单元起点的节点编号,j 单元终点的节点编号
x1 y1 起点的坐标,x2 y2 终点的坐标,
DOF 单元在总体刚度矩阵里面所在的位置,L 单元的长度,
cos sin 单元的方向余弦 方向正弦,
k 单元刚度矩阵"""
def __init__(self, E, A, x1, y1, x2, y2, i, j):
self.E=E
self.A=A
self.i=i
self.j=j#定义一个由单刚矩阵的自由度向总刚矩阵自由度转换的数组
self.DOF = np.array([2*i-2, 2*i-1, 2*j-2, 2*j-1],dtype=np.int16)
self.L= sqrt((x1 - x2)**2 + (y1 - y2)**2)
self.cos= (x2 - x1) /self.L
self.sin= (y2 - y1) /self.L
L, c, s=self.L, self.cos, self.sin
self.k= (E*A/L)*np.array([[c*c, c*s, -c*c, -c*s],
[c*s, s*s, -c*s, -s*s],
[-c*c, -c*s, c*c, c*s],
[-c*s, -s*s, c*s, s*s]])"定义求解单元应力的函数"
defstress(U):#从位移矩阵U中获取该单元两个节点的1*4位移矩阵
u =U(self.DOF)
E, L, c, s=self.E, self.L, self.c, self.s
T= np.array([-c, -s, c, s])
self.bar_Stress= E / L *np.dot(T, u)returnself.bar_Stress"定义总刚矩阵集成函数"
defBar2D2Node_Assembly(KK, bar):for n1 in xrange(4):for n2 in xrange(4):
KK[bar.DOF[n1], bar.DOF[n2]]+=bar.k[n1, n2]returnKK‘求解节点位移‘
defnode_Disaplacement(KK, U, F):#获取缩减后的总刚矩阵
del_row_col = np.where(U ==0)
kk_delRow=np.delete(KK, del_row_col, 0)
kk_delCol= np.delete(kk_delRow, del_row_col,1)
kk=kk_delCol#获取节点位移位置对应的节点力矩阵
f = F[np.where(U == ‘u_Unknown‘)]
u=np.linalg.solve(kk, f)
U[np.where(U==‘u_Unknown‘)] =ureturnU‘求解节点力,必须在已经求得节点位移U后才可调用本函数‘
defnode_Force(KK, U):
F=np.dot(KK, U)return F
仍然以书上的例题为例
结构的离散化同书中一致
程序中Bar2D这个类表示杆单元,实例化的时候,需要提供一下信息:
弹性模量E,面积A,起点i和终点j的编号以及各自的坐标。
所以对于本例题,这几个信息如下:
E, A, x1, y1, x2, y2, x3, y3, x4, y4 = 2.95e11, 0.0001, 0, 0, 0.4, 0, 0.4, 0.3, 0, 0.3bar1= Bar2D(E, A, x1, y1, x2, y2, 1, 2)
bar2= Bar2D(E, A, x3, y3, x2, y2, 3, 2)
bar3= Bar2D(E, A, x1, y1, x3, y3, 1, 3)
bar4= Bar2D(E, A, x4, y4, x3, y3, 4, 3)
bars=[bar1, bar2, bar3, bar4]for bar inbars:
Bar2D2Node_Assembly(KK, bar)
然后调用相应的函数求解位移和荷载即可:
#求解位移
U =node_Disaplacement(KK, U, F)#求解节点力
F = node_Force(KK, U)
最终求得的结果如下:
可看到最终结果与书中求得的结果是一致的,由于单位制采用m为单位,所以求得的位移单位也是m。
--------------------------------------------------------------------------
写这段代码的时候,让我想起了当年我上靳慧老师的《有限元方法》这门课。还记得当时前面几章的基础部分是靳老师授课,后面的章节是暑期的一个小学期由加州大学尔湾分校的孙立志教授授课,孙教授采用的是全英文授课。那是我唯一的一次听全英文授课,虽然稍微有点吃力,但大部分还是可以听得懂的。国外老师的授课思路与国内老师还是有些区别,国外老师一般从基本的属性方法讲起,我记得当时孙老师用了不少篇幅讲加权残值法和伽辽金方法。虽然现在忘得差不多了,但那是一段难忘的经历,炎热的夏天,中山院,满头大汗的孙教授。
靳老师中间布置过一次作业,是一个三维杆单元结构的例题,让我们自己编写程序求解该例题,然后使用通用有限元程序分析,比较二者的结果。编程语言不限制,但所有人都不由自主的选择了MATLAB。本来我了逼格,我还想用FORTRAN的,但权衡了一下还是算了,学起来太麻烦。最后还是选择了MATLAB。现在回想起来,当时写的代码太烂了。用了一天入门matlab,然后赶鸭子上架搞出来的,为了声明带下标的变量生生用了N多行eval函数,不堪回首。
原文:http://www.cnblogs.com/SimuLife/p/4695922.html
python有限元例题_《有限元分析基础教程》(曾攀)笔记一-二维杆单元有限元程序(基于Python)...相关推荐
- matlab有限元分析教程,MATLABprogramin有限元分析基础教程曾攀.pdf
MATLABprogramin有限元分析基础教程曾攀.pdf 限元分析基础教程 曾攀 3.3.6 梁单元分析的MATLAB 程序 [MATLAB 程序]3.3.6(1) 1D 梁单元的有限元分析程序( ...
- 《有限元分析基础教程》(曾攀)笔记一-二维杆单元有限元程序(基于Python)...
曾攀老师的<有限元分析基础教程>第三章有二维杆单元的推导,并结合一个例题进行了解析解和基于Matlab的程序求解.但是我感觉书中的MATLAB代码有点罗嗦,而且一些实现方法也比较麻烦,比如 ...
- Matlab基础教程—【07】Matlab二维高层绘图操作
7.1 二维高层绘图的基本函数plot() 重要参考资料: ① 基础教程视频对应的操作纪录 ② 基础教程视频对应的PPT 辅助参考资料:Matlab基本绘图函数 1. plot()有两个参数 (1)基 ...
- 《语义网基础教程》学习笔记(二)
二.RDF概述(參考http://zh.transwiki.org/cn/rdfprimer.htm) 1.本体: 一个本体是一个概念体系(conceptualization)的显式的形式化规范. 一 ...
- python基础代码库-python爬虫基础教程:requests库(二)代码实例
get请求 简单使用 import requests ''' 想要学习Python?Python学习交流群:973783996满足你的需求,资料都已经上传群文件,可以自行下载! ''' respons ...
- python爬虫requests库_python爬虫基础教程:requests库(二)代码实例
get请求 简单使用 import requests ''' 想要学习Python?Python学习交流群:973783996满足你的需求,资料都已经上传群文件,可以自行下载! ''' respons ...
- CG基础教程-陈惟老师十二讲笔记
转自 麽洋TinyOcean:http://www.douban.com/people/Tinyocean/notes?start=50&type=note 因为看了陈惟十二讲视频没有课件,边 ...
- JavaScript基础教程速学笔记
JavaScript基础教程速学笔记 JavaScript简介 JavaScript 是 Web 的编程语言.(但是java与JavaScript的区别就是周杰与周杰伦的区别)所有现代的 HTML 页 ...
- 零基础学前端之HTML全套基础教程【学习笔记】
[前端总路线学习笔记] 文章目录 HTML全套基础教程[学习笔记] 1.系统结构 2.软件环境准备 3.HTML概述 4. 我的第一个HTML 5. HTML的基本标签 6.HTML的实体符号 7. ...
最新文章
- 任命新CFO 百度迎来首位女性高管
- ios与html数据交互,iOS iOS与html进行交互
- 测度定义_Real analysis:外测度的一个等价定义
- G6 图可视化引擎——入门教程——绘制 Tutorial 案例
- mmdetection多类目标训练查看单类准确率(AP)以及使用模型测试看结果(show)
- 在两个不同域中的WINDOWS 2003活动目录做迁移笔记
- 计算机建立excel文件,用Excel建立数据库 -电脑资料
- JQuery AJAX 的表单提交
- linux中ess33没有IP地址问题
- django实现利用邮箱进行登录
- 电脑连不上ishanghai_ishanghai用电脑肿么连网
- 安卓开发无线连接设备进行调试(adb)
- 关于学计算机趣味段子,让你开怀大笑的段子:幽默风趣,读一遍笑一遍!
- OpenWrt软件源清华大学镜像
- 2013年杰森·斯坦森动作《蜂鸟》720p.BD中英双字幕
- zuk android os 流量,原生用户最爱 Cyanogen OS版ZUK Z1固件
- MD5(单向散列算法)原理分析
- 新手选车系列之(八): 选车购车谨慎采取“一票否决制”
- 浅析城市交通现状及问题
- ECCV 2018所有论文合集
热门文章
- 鼠标连点器烦人弹窗?我直接爆破
- 如何计算心跳c语言编程,单片机心率计 电子脉搏计设计(原理图Protues仿真和C程序)...
- 新玺配资:指数回调 行情可期
- zbar下条形码和二维码的识别与解码Ⅰ|2021SC@SDUSC
- 论文阅读Generalizing A Person Retrieval Model Hetero-and Homogeneously
- IOS-navigationController切换页面失效问题最新解决办法
- Android 电量显示Widgets插件实现
- HDU - 1237简单计算器(输出问题)
- 安装配置Discuz论坛
- 各种杂项组件(1)之--Toast(提示框)