|
Feon 是以有限元编程教学和算法研究为目的，基于 Python 开发的有限元分析框架，其开发者为湖南工业大学土木建筑与环境学院的 Dr. Pei Yaoyao。目前的版本为 1.0.0 版本。
1. 安装方法
安装 Feon 之前，首先需要安装拓展包：Numpy、可视化组件 Matplotlib、矩阵运算 Mpmath。
安装有两种途径：
Using pip:$pip install feon
或者
$Python setup.py install
框架工包含三个包：
sa：结构分析包 ffa：流体分析包 derivation：刚度矩阵包
2. 支持的单元类型
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 .
3. 解决算例
3.1 二维桁架问题
# -*- coding: utf-8 -*-# ------------------------------------# Author: YAOYAO PEI# E-mail: [email protected]# License: Hubei University of Technology License# -------------------------------------from feon.sa import *from feon.tools import pair_wisefrom feon.sa.draw2d import *import matplotlib.pyplot as pltfrom matplotlib.ticker import MultipleLocatorif __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 #create nodes and elements 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[i],nds2[i]),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[i]),E,A2)) #create FEA system s = System() #add nodes and elements into the system s.add_nodes(nds1,nds2) s.add_elements(els) #apply boundry condition 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[i].ID,Fy = -1900) s.add_fixed_sup(nds1[0].ID) s.add_rolled_sup(nds1[-1].ID,"y") #solve the system s.solve() #show results 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[i],eforce[i]],"ks-",ms = 3) plt.show() draw_bar_info(els[5])
3.2 带铰接点的钢架问题
# -*- coding: utf-8 -*-# ------------------------------------# Author: YAOYAO PEI# E-mail: [email protected]# License: Hubei University of Technology License# -------------------------------------from feon.sa import *from feon.tools import pair_wise#define beamlink elementclass BeamLink2D11(StructElement): def __init__(self,nodes,E,A,I): StructElement.__init__(self,nodes) self.E = E self.A = A self.I = I #define node degree of freedom, left node has three dofs #while the right node has only two def init_unknowns(self): self.nodes[0].init_unknowns("Ux","Uy","Phz") self.nodes[1].init_unknowns("Ux","Uy") self._ndof = 3 #transformative matrix 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 #stiffness matrix 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 Tif __name__ == "__main__": #materials E = 210e6 A = 0.005 I = 10e-5 #nodes and elements 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) #system 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.3 地下连续墙问题
请您登录后阅读全文， 登录 或者 注册
|