《有限元分析基础教程》（曾攀）笔记一-二维杆单元有限元程序（基于Python）...
曾攀老師的《有限元分析基礎(chǔ)教程》第三章有二維桿單元的推導(dǎo)，并結(jié)合一個(gè)例題進(jìn)行了解析解和基于Matlab的程序求解。但是我感覺書中的MATLAB代碼有點(diǎn)羅嗦，而且一些實(shí)現(xiàn)方法也比較麻煩，比如已經(jīng)知道了桿單元的起點(diǎn)和終點(diǎn)坐標(biāo)，仍然需要另外給出單元局部坐標(biāo)與整體坐標(biāo)的夾角，這完全沒必要。于是我就用Python重構(gòu)了這段程序，當(dāng)然并不是把書中的MATLAB代碼翻譯成python（事實(shí)上完全可以這么干，而且很快！）。比如我使用了面向?qū)ο蟮乃枷?#xff0c;把桿單元寫成了一個(gè)類，這樣思路比較清晰。
#! /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表示節(jié)點(diǎn)的位移向量，如果某個(gè)自由度的位移未知，則該處填寫‘u_Unknown’， F表示節(jié)點(diǎn)的荷載向量，如果某個(gè)自由度的荷載未知，則該處填寫‘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 單元起點(diǎn)的節(jié)點(diǎn)編號(hào)，j 單元終點(diǎn)的節(jié)點(diǎn)編號(hào)x1 y1 起點(diǎn)的坐標(biāo)，x2 y2 終點(diǎn)的坐標(biāo)，DOF 單元在總體剛度矩陣?yán)锩嫠诘奈恢?L 單元的長(zhǎng)度，cos sin 單元的方向余弦 方向正弦，k 單元?jiǎng)偠染仃?/span>"""def __init__(self, E, A, x1, y1, x2, y2, i, j):self.E = Eself.A = Aself.i = iself.j = j# 定義一個(gè)由單剛矩陣的自由度向總剛矩陣自由度轉(zhuǎn)換的數(shù)組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.Lself.sin = (y2 - y1) / self.LL, c, s = self.L, self.cos, self.sinself.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]])"定義求解單元應(yīng)力的函數(shù)"def stress(U):# 從位移矩陣U中獲取該單元兩個(gè)節(jié)點(diǎn)的1*4位移矩陣u = U(self.DOF)E, L, c, s = self.E, self.L, self.c, self.sT = np.array([-c, -s, c, s])self.bar_Stress = E / L * np.dot(T, u)return self.bar_Stress"定義總剛矩陣集成函數(shù)" 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'求解節(jié)點(diǎn)位移' 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# 獲取節(jié)點(diǎn)位移位置對(duì)應(yīng)的節(jié)點(diǎn)力矩陣f = F[np.where(U == 'u_Unknown')]u = np.linalg.solve(kk, f)U[np.where(U=='u_Unknown')] = ureturn U'求解節(jié)點(diǎn)力，必須在已經(jīng)求得節(jié)點(diǎn)位移U后才可調(diào)用本函數(shù)' def node_Force(KK, U):F = np.dot(KK, U)return F仍然以書上的例題為例
結(jié)構(gòu)的離散化同書中一致
程序中Bar2D這個(gè)類表示桿單元，實(shí)例化的時(shí)候，需要提供一下信息：
彈性模量E，面積A，起點(diǎn)i和終點(diǎn)j的編號(hào)以及各自的坐標(biāo)。
所以對(duì)于本例題，這幾個(gè)信息如下：
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)然后調(diào)用相應(yīng)的函數(shù)求解位移和荷載即可：
# 求解位移 U = node_Disaplacement(KK, U, F) # 求解節(jié)點(diǎn)力 F = node_Force(KK, U)最終求得的結(jié)果如下：
可看到最終結(jié)果與書中求得的結(jié)果是一致的，由于單位制采用m為單位，所以求得的位移單位也是m。
--------------------------------------------------------------------------
寫這段代碼的時(shí)候，讓我想起了當(dāng)年我上靳慧老師的《有限元方法》這門課。還記得當(dāng)時(shí)前面幾章的基礎(chǔ)部分是靳老師授課，后面的章節(jié)是暑期的一個(gè)小學(xué)期由加州大學(xué)爾灣分校的孫立志教授授課，孫教授采用的是全英文授課。那是我唯一的一次聽全英文授課，雖然稍微有點(diǎn)吃力，但大部分還是可以聽得懂的。國(guó)外老師的授課思路與國(guó)內(nèi)老師還是有些區(qū)別，國(guó)外老師一般從基本的屬性方法講起，我記得當(dāng)時(shí)孫老師用了不少篇幅講加權(quán)殘值法和伽遼金方法，雖然現(xiàn)在忘得差不多了，但那是一段難忘的經(jīng)歷。。
靳老師中間布置過一次作業(yè)，是一個(gè)三維桿單元結(jié)構(gòu)的例題，讓我們自己編寫程序求解該例題，然后使用通用有限元程序分析，比較二者的結(jié)果。編程語(yǔ)言不限制，但所有人都不由自主的選擇了MATLAB。本來我了逼格，我還想用FORTRAN的，但權(quán)衡了一下還是算了，學(xué)起來太麻煩。最后還是選擇了MATLAB。現(xiàn)在回想起來，當(dāng)時(shí)寫的代碼太爛了。用了一天入門matlab，然后趕鴨子上架搞出來的，為了聲明帶下標(biāo)的變量生生用了N多行eval函數(shù)，不堪回首。
轉(zhuǎn)載于:https://www.cnblogs.com/SimuLife/p/4695922.html
總結(jié)
以上是生活随笔為你收集整理的《有限元分析基础教程》（曾攀）笔记一-二维杆单元有限元程序（基于Python）...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 前端开发那点事儿
- 下一篇： Python入门——一个沙雕的表情包