《有限元分析基础教程》（曾攀）笔记一-二维杆单元有限元程序（基于Python）

曾攀老师的《有限元分析基础教程》第三章有二维杆单元的推导，并结合一个例题进行了解析解和基于Matlab的程序求解。但是我感觉书中的MATLAB代码有点罗嗦，而且一些实现方法也比较麻烦，比如已经知道了杆单元的起点和终点坐标，仍然需要另外给出单元局部坐标与整体坐标的夹角，这完全没必要。于是我就用Python重构了这段程序，当然并不是把书中的MATLAB代码翻译成python（事实上完全可以这么干，而且很快！）。比如我使用了面向对象的思想，把杆单元写成了一个类，这样思路比较清晰。

#! /usr/bin/python
# -*- coding: utf-8 -*-

import math
import numpy as np 

sqrt, cos, sin, pi = math.sqrt, math.cos, math.sin, math.pi
"前处理"
nodeNumber = 4
KK = 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)


class Bar2D:
    """定义二维杆单元类，该类包含杆件的基本信息：
    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]])
    "定义求解单元应力的函数"
    def stress(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)
        return self.bar_Stress

"定义总刚矩阵集成函数"
def Bar2D2Node_Assembly(KK, bar):
    for  n1 in xrange(4):
        for n2 in xrange(4):
            KK[bar.DOF[n1], bar.DOF[n2]] += bar.k[n1, n2]
    return KK

'求解节点位移'
def node_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')] = u
    return U
    

'求解节点力，必须在已经求得节点位移U后才可调用本函数'
def node_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.3
bar1 = 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 in bars:
    Bar2D2Node_Assembly(KK, bar)

然后调用相应的函数求解位移和荷载即可：

# 求解位移
U = node_Disaplacement(KK, U, F)
# 求解节点力
F = node_Force(KK, U)

最终求得的结果如下：

可看到最终结果与书中求得的结果是一致的，由于单位制采用m为单位，所以求得的位移单位也是m。

--------------------------------------------------------------------------

写这段代码的时候，让我想起了当年我上靳慧老师的《有限元方法》这门课。还记得当时前面几章的基础部分是靳老师授课，后面的章节是暑期的一个小学期由加州大学尔湾分校的孙立志教授授课，孙教授采用的是全英文授课。那是我唯一的一次听全英文授课，虽然稍微有点吃力，但大部分还是可以听得懂的。国外老师的授课思路与国内老师还是有些区别，国外老师一般从基本的属性方法讲起，我记得当时孙老师用了不少篇幅讲加权残值法和伽辽金方法，虽然现在忘得差不多了，但那是一段难忘的经历。。

靳老师中间布置过一次作业，是一个三维杆单元结构的例题，让我们自己编写程序求解该例题，然后使用通用有限元程序分析，比较二者的结果。编程语言不限制，但所有人都不由自主的选择了MATLAB。本来我了逼格，我还想用FORTRAN的，但权衡了一下还是算了，学起来太麻烦。最后还是选择了MATLAB。现在回想起来，当时写的代码太烂了。用了一天入门matlab，然后赶鸭子上架搞出来的，为了声明带下标的变量生生用了N多行eval函数，不堪回首。