《有限元分析基礎教程》（曾攀）筆記一-二維杆單元有限元程式（基於Python）

標籤：

曾攀老師的《有限元分析基礎教程》第三章有二維杆單元的推導，並結合一個例題進行瞭解析解和基於Matlab的程式求解。但是我感覺書中的MATLAB代碼有點羅嗦，而且一些實現方法也比較麻煩，比如已經知道了杆單元的起點和終點座標，仍然需要另外給出單元局部座標與整體座標的夾角，這完全沒必要。於是我就用Python重構了這段程式，當然並不是把書中的MATLAB代碼翻譯成python（事實上完全可以這麼幹，而且很快！）。比如我使用了物件導向的思想，把杆單元寫成了一個類，這樣思路比較清晰。

#! /usr/bin/python# -*- coding: utf-8 -*-import mathimport numpy as np sqrt, cos, sin, pi = math.sqrt, math.cos, math.sin, math.pi"前處理"nodeNumber = 4KK = 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.3bar1 = 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函數，不堪回首。

《有限元分析基礎教程》（曾攀）筆記一-二維杆單元有限元程式（基於Python）