- 积分
- 0
- 注册时间
- 2013-8-30
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 Iwin4FuN 于 2018-11-23 15:09 编辑
Feon简介
Feon是湖北工业大学土木建筑与环境学院教师裴尧尧基于Python开发的一个开源免费的有限元计算框架。
版本
当前版本 1.01
安装需求
需要Numpy,可视化的话需要matplotlib,单元矩阵推到需要mpmath。
安装
使用pip安装 pip install feon,或者源码安装 python setup.py install
有限元子包
* sa---结构分析
* ffa --- 渗流分析(目前只支持一维)
* derivation --- 刚度矩阵推导
自带的单元
* Spring1D11 -一维弹簧单元
* Spring2D11 -二维弹簧单元
* Spring3D11 -三维弹簧单元
* Link1D11 - 一维杆单元
* Link2D11 - 二维弹簧单元
* Link3D11 - 三维弹簧单元
* Beam1D11- 一维梁单元
* Beam2D11-二维弹簧单元
* Beam3D11- 三维弹簧单元
* Tri2d11S---- 平面应力三角形单元
* Tri2D11 ---- 平面应变三角形单元
* Tetra3D11 ---- 四面体单元
* Quad2D11S ---- 平面应力 四边形单元
* Quad2D11 ----平面应变 四边形单元
* Plate3D11 ---Midline 板单元
* Brick3D11 ---六面体单元
单元命名方式为 "Name" + "dimension" + 'order" + "type", type 1 means elastic .
代码实例
1 二维桁架问题(图形见附件)
# -*- coding: utf-8 -*-
# ------------------------------------
# Author: YAOYAO PEI
# E-mail: yaoyao.bae@foxmail.com
# License: Hubei University of Technology License
# -------------------------------------
from feon.sa import *
from feon.tools import pair_wise
from feon.sa.draw2d import *
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
if __name__ == "__main__":
#material property
E = 210e6 #elastic modulus
A1 = 31.2e-2 #cross-section area of hanging bars
A2 = 8.16e-2 #cross-section area of others
#创建节点和单元
nds1 = []
nds2 = []
for i in range(13):
nds1.append(Node(i,0))
for i in range(11):
nds2.append(Node(i+1,-1))
els = []
for e in pair_wise(nds1):
els.append(Link2D11((e[0],e[1]),E,A1))
for e in pair_wise(nds2):
els.append(Link2D11((e[0],e[1]),E,A1))
for i in range(6):
els.append(Link2D11((nds1,nds2),E,A2))
for i in xrange(6):
els.append(Link2D11((nds2[i+5],nds1[i+7]),E,A2))
for i in range(11):
els.append(Link2D11((nds1[i+1],nds2),E,A2))
#创建有限元系统
s = System()
#将节点和单元加入到系统中
s.add_nodes(nds1,nds2)
s.add_elements(els)
#添加边界条件
s.add_node_force(nds1[0].ID,Fy = -1000)
s.add_node_force(nds1[-1].ID,Fy = -1000)
for i in range(1,12):
s.add_node_force(nds1.ID,Fy = -1900)
s.add_fixed_sup(nds1[0].ID)s.add_rolled_sup(nds1[-1].ID,"y")#求解系统
s.solve()
#显示结果
disp = [np.sqrt(nd.disp["Ux"]**2+nd.disp["Uy"]**2) for nd in s.get_nodes()]
eforce = [el.force["N"][0][0] for el in s.get_elements()]
fig = plt.figure()
ax = fig.add_subplot(211)
ax.yaxis.get_major_formatter().set_powerlimits((0,1))
ax2 = fig.add_subplot(212)
ax2.yaxis.get_major_formatter().set_powerlimits((0,1))
ax.set_xlabel(r"$Node ID$")
ax.set_ylabel(r"$Disp/m$")
ax.set_ylim([-0.05,0.05])
ax.set_xlim([-1,27])
ax.xaxis.set_minor_locator(MultipleLocator(1))
ax.plot(range(len(disp)),disp,"r*-")
ax2.set_xlabel(r"$Element ID$")
ax2.set_xlim([-1,46])
ax2.set_ylabel(r"$N/kN$")
ax2.set_ylim(-40000,40000)
ax2.xaxis.set_minor_locator(MultipleLocator(1))
for i in range(len(eforce)):
ax2.plot([i-0.5,i+0.5],[eforce,eforce],"ks-",ms = 3)
plt.show()
draw_bar_info(els[5])
2 带铰接点的钢架问题(图形见附件)
# -*- coding: utf-8 -*-
# ------------------------------------
# Author: YAOYAO PEI
# E-mail: yaoyao.bae@foxmail.com
# License: Hubei University of Technology License
# -------------------------------------
from feon.sa import *
from feon.tools import pair_wise
#定义铰接刚接单元
class BeamLink2D11(StructElement):
def __init__(self,nodes,E,A,I):
StructElement.__init__(self,nodes)
#定义材料参数
self.E = E
self.A = A
self.I = I
#设置单元节点自由度
#左边节点为3,右边节点为2
def init_unknowns(self):
self.nodes[0].init_unknowns("Ux","Uy","Phz")
self.nodes[1].init_unknowns("Ux","Uy")
self._ndof = 3
#定义坐标转换矩阵
def calc_T(self):
TBase = _calc_Tbase_for_2d_beam(self.nodes)
self._T = np.zeros((6,6))
self._T[:3,:3] = self._T[3:,3:] = TBase
#定义局部刚度矩阵
def calc_ke(self):
self._ke = _calc_ke_for_2d_beamlink(E = self.E,A = self.A,I = self.I,L = self.volume)
def _calc_ke_for_2d_beamlink(E = 1.0,A = 1.0,I = 1.0,L = 1.0):
a00 = E*A/L
a03 = -a00
a11 = 3.*E*I/L**3
a12 = 3.*E*I/L**2
a14 = -a11
a22 = 3.*E*I/L
T = np.array([[a00, 0., 0., a03, 0.,0.],
[ 0., a11, a12, 0., a14, 0.],
[ 0., a12, a22, 0.,-a12, 0.],
[a03, 0., 0., a00, 0., 0.],
[ 0., a14, -a12, 0., a11, 0.],
[ 0., 0., 0., 0., 0., 0.]])
return T
def _calc_Tbase_for_2d_beam(nodes):
x1,y1 = nodes[0].x,nodes[0].y
x2,y2 = nodes[1].x,nodes[1].y
le = np.sqrt((x2-x1)**2+(y2-y1)**2)
lx = (x2-x1)/le
mx = (y2-y1)/le
T = np.array([[lx,mx,0.],
[-mx,lx,0.],
[0.,0.,1.]])
return T
if __name__ == "__main__":
#定义材料参数
E = 210e6
A = 0.005
I = 10e-5
#创建节点和单元
n0 = Node(0,0)
n1 = Node(0,3)
n2 = Node(4,3)
n3 = Node(4,0)
n4 = Node(4,5)
n5 = Node(8,5)
n6 = Node(8,0)
e0 = Beam2D11((n0,n1),E,A,I)
e1 = BeamLink2D11((n1,n2),E,A,I)
e2 = Beam2D11((n2,n3),E,A,I)
e3 = Beam2D11((n2,n4),E,A,I)
e4 = Beam2D11((n4,n5),E,A,I)
e5 = Beam2D11((n5,n6),E,A,I)
#创建有限元系统
s = System()
s.add_nodes([n0,n1,n2,n3,n4,n5,n6])
s.add_elements([e0,e1,e2,e3,e4,e5])
s.add_node_force(1,Fx = -10)
s.add_node_force(5,Fx = -10)
s.add_fixed_sup(0,3,6)
s.solve()
print n2.disp
print e1.force
3 地下连续墙问题
# -*- coding: utf-8 -*-
# ------------------------------------
# Author: YAOYAO PEI
# E-mail: yaoyao.bae@foxmail.com
# License: Hubei University of Technology License
# -------------------------------------
from feon.sa import *
from feon.tools import pair_wise
import matplotlib.pyplot as plt
from feon.sa.draw2d import *
if __name__ == "__main__":
#材料参数
E1 = 2.85e6 #墙体弹性模量
E2 = 200e6 #支撑弹性模量
k = 15000 #基床系数
I = 0.0427 #墙的惯性矩
A = 0.8 # 墙体截面面积
A1 = 0.003 #支撑截面面积
ka = 0.6 #主动土压力系数
#创建节点
nds1 =[Node(0,-i) for i in range(10)]
nds2 = [Node(0,-(i+20)*0.5) for i in range(11)]
nds3 = [Node(-0.5,-(i+20)*0.5) for i in range(11)]
nds4 = [Node(-1.5,-2),Node(-1.5,-6)]
#创建梁单元
els=[]
for nd in pair_wise(nds1+nds2):
els.append(Beam2D11(nd,E1,A,I))
#创建土弹簧
for i in range(11):
els.append(Spring2D11((nds2,nds3),k))
#创建支撑
els.append(Link2D11((nds4[0],nds1[2]),E2,A1))
els.append(Link2D11((nds4[1],nds1[6]),E2,A1))
#创建有限元系统并添加节点和单元
s = System()
s.add_nodes(nds1,nds2,nds3,nds4)
s.add_elements(els)
nid1 = [nd.ID for nd in nds3]
nid2 = [nd.ID for nd in nds4]
#施加边界条件
s.add_fixed_sup(nid1,nid2)
for i,el in enumerate(els[:10]):
s.add_element_load(el.ID,"tri",-18*ka)
s.add_element_load(el.ID,"q",-i*18*ka)
#施加主动土压力
for el in els[10:20]:
s.add_element_load(el.ID,"q",-180*ka)
for nd in nds1:
nd.set_disp(Uy =0)
for nd in nds2:
nd.set_disp(Uy = 0)
#求解
s.solve()
#显示结果
disp = np.array([nd.disp["Ux"] for nd in nds1]+[nd.disp["Ux"] for nd in nds2])*1000
Mz = [el.force["Mz"][0][0] for el in els[:20]]
fig1,fig2,fig3 = plt.figure(),plt.figure(),plt.figure()
ax1 = fig1.add_subplot(111)
ax2 = fig2.add_subplot(111)
ax3 = fig3.add_subplot(111)
Y1 = [-i for i in range(10)]+[-(i+20)*0.5 for i in range(11)]
Y2 = [-i-0.5 for i in range(10)]+[-(i+20)*0.5-0.5 for i in range(10)]
ax1.plot(disp,Y1,"r--")
ax1.set_xlabel("$Ux/mm$")
ax1.set_ylabel("$Height/m$")
ax2.plot(Mz,Y2,"r-+")
ax2.set_xlabel("$Mz/kN.m$")
ax2.set_ylabel("$Height/m$")
for el in els[:20]:
draw_element(ax3,el,lw = 10,color = "g")
for el in els[20:31]:
draw_spring(ax3,el,color = "k")
for el in els[31:]:
draw_element(ax3,el,lw = 1.5,color = "k",marker = "s")
for nd in nds3+nds4:
draw_fixed_sup(ax3,nd,factor = (0.4,4),color ="k")
ax3.set_xlim([-2,2])
ax3.set_ylim([-16,1])
plt.show()
4、一维渗流问题
# -*- coding: utf-8 -*-
# ------------------------------------
# Author: YAOYAO PEI
# E-mail: yaoyao.bae@foxmail.com
# License: Hubei University of Technology License
# -------------------------------------
from feon.ffa import *
from feon.tools import pair_wise
import numpy as np
if __name__ == "__main__":
#渗透系数
Kxx = -2e-5
#创建节点和单元
A = np.pi*(np.linspace(0.06,0.15,7)[:-1]+0.0075)
nds = [Node(-i*0.1,0) for i in range(7)]
els = []
for i in range(6):
els.append(E1D((nds,nds[i+1]),Kxx,A))
#创建有限元系统
s = System()
s.add_nodes(nds)
s.add_elements(els)
s.add_node_head(0,0.2)
s.add_node_head(6,0.1)
#求解系统
s.solve()
#显示结果
print [nd.head["H"] for nd in nds]
print [el.velocity["Vx"] for el in els]
更多使用例子可关注《Python与有限元》一书,或者关注开源代码,地址为https://github.com/YaoyaoBae/Feon
QQ交流群 555809224
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
|