走进有限元算法：基于 Python 的有限元分析框架 Feon

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: 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
    #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: yaoyao.bae@foxmail.com
#  License: Hubei University of Technology License
# -------------------------------------
from feon.sa import *
from feon.tools import pair_wise
#define beamlink element
class 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 T
if __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 地下连续墙问题