《有限元分析基礎教程》（曾攀）筆記一-二維桿單元有限元程序（基於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函數，不堪回首。