本文檔為有限元方法課程作業(yè)，使用python編寫，可實現(xiàn)任意單元數(shù)量的結(jié)構(gòu)計算。鑒于搜索到這篇文章的一般都是選了這門課的同學(xué)，所以基礎(chǔ)知識不多說直接開始：

1 ?程序設(shè)計

1.1 ?設(shè)計思路

有限元程序的基本功能為輸入結(jié)構(gòu)的幾何信息、材料參數(shù)、荷載和約束情況，輸出節(jié)點位移、單元應(yīng)力應(yīng)變等結(jié)果。程序采用python語言編寫，二維數(shù)據(jù)和一維數(shù)據(jù)均使用mumpy矩陣存放。

（1）生成結(jié)構(gòu)剛度矩陣

結(jié)構(gòu)剛度矩陣是有限元計算的核心?？梢允褂醚h(huán)結(jié)構(gòu)，依次生成各單元的單元剛度矩陣并添加到結(jié)構(gòu)剛度矩陣中，每次循環(huán)生成一個單元剛度矩陣并立即添加到結(jié)構(gòu)剛度矩陣中，下一次循環(huán)時不保留上一個單元的單元剛度矩陣，以節(jié)約內(nèi)存。為節(jié)省時間，程序中的單元剛度矩陣直接按式 1生成：

分塊矩陣表示為：

? ?

? ? ??

其中：

? ? ? ? ? ? ? ? ? ? ? ? ? ??

（2）生成荷載向量

荷載向量可由輸入的荷載信息矩陣直接生成，維度為2n，n為節(jié)點數(shù)。

（3）處理邊界條件

采用對角元素置1法，將剛度矩陣中與被固定節(jié)點對應(yīng)的對角元素置1，該行和列其余元素置0，同時將荷載向量對應(yīng)元素置0?？稍谧映绦蛑猩昝魅肿兞浚苯有薷膭偠染仃嚭秃奢d向量。

（4）求解大型稀疏矩陣方程

整體剛度矩陣是一個大型帶狀矩陣，對于這類矩陣產(chǎn)生的大型稀疏矩陣方程，可以調(diào)用scipy庫中的sparse.linalg.spsolve函數(shù)求解。

（5）求單元應(yīng)變

根據(jù)得到的節(jié)點位移，可由式 4~式5得到單元應(yīng)變εx,εy,γxy。用循環(huán)結(jié)構(gòu)循環(huán)計算單元應(yīng)變并存入應(yīng)變信息矩陣輸出。

? ? ? ? ? ??

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

（6）求單元應(yīng)力

根據(jù)得到的節(jié)點位移，可由式6~式7得到單力σx,σy,τxy。用循環(huán)結(jié)構(gòu)循環(huán)計算單元應(yīng)力并存入應(yīng)力信息矩陣輸出。

? ? ? ? ?

? ? ? ? ? ? ? ?式6

? ? ? ? ? ? ? ? ? ? ? ? ? ? ?式7

（6）輸出結(jié)果

輸出結(jié)果包括節(jié)點位移、單元應(yīng)變和單元應(yīng)力。

1.2 ?程序展示

（1）主要變量

EL：彈性模量

PO：泊松比

TK：厚度

NODE_SUM：節(jié)點數(shù)目

ELE_SUM：單元數(shù)目

NODE：節(jié)點坐標(biāo)矩陣，每行為節(jié)點坐標(biāo)，列按節(jié)點號排列

ELE：單元節(jié)點號矩陣，每行為單元的三個節(jié)點號，列按單元號排列，單元內(nèi)節(jié)點號逆時針排列

P_inp：荷載信息矩陣，第一列為荷載作用的節(jié)點號，第二列為x方向荷載，第三列為y方向荷載

BC：約束信息矩陣，第一列為節(jié)點號，二~三列分別為x和y向的約束情況，1為固定，0為自由

TYPE：平面問題類型，1為平面應(yīng)變問題，否則為平面應(yīng)力問題

（2）子函數(shù)

create_K( EL, PO, TK, NODE_SUM,ELE_SUM, NODE, ELE):生成整體剛度矩陣

create_P(P_inp,NODE_SUM):由輸入的節(jié)點荷載信息生成荷載向量P

do_BC(BC):處理邊界條件，置1法對角元素

solve_Uout(U,NODE_SUM):根據(jù)位移向量得到節(jié)點位移信息矩陣

solve_B(U,NODE,ELE,ELE_SUM):根據(jù)位移向量U得到結(jié)構(gòu)應(yīng)變信息矩陣

solve_S(U, NODE, ELE, ELE_SUM, EL, PO):根據(jù)位移向量U得到結(jié)構(gòu)應(yīng)力信息矩陣

有限元程序完整代碼如下：

2 ?算例驗證

建立有限元算例如圖 1所示。本文所有數(shù)值均使用國際單位制，即長度：m；時間：s，力：N，應(yīng)力：Pa。后續(xù)數(shù)值出現(xiàn)時不帶單位。

圖1 有限元算例結(jié)構(gòu)

本代碼不包含前處理部分，通過更改主程序中第一部分關(guān)鍵參數(shù)的形式輸入?yún)?shù)，按照注釋行要求輸入關(guān)鍵參數(shù)即可實現(xiàn)任意單元數(shù)量結(jié)構(gòu)的計算。彈性模量設(shè)置為20，泊松比為0.2。以2點為坐標(biāo)原點，1號單元長直角邊為x軸，短直角邊為y軸，約束1，2點x向和y向位移。荷載數(shù)值為1，作用于5號節(jié)點，方向豎直向下。主程序第一部分輸入?yún)?shù)為：

程序運行后會在當(dāng)前文件夾產(chǎn)生一個名為“result.txt”的結(jié)果文本，文本內(nèi)容如下：

——————###結(jié)果文件###——————

##節(jié)點位移##

節(jié)點號 x y

——————————————————

1.00000000 0.00000000 0.00000000

2.00000000 0.00000000 0.00000000

3.00000000 -0.21617469 -0.63980590

4.00000000 0.22453700 -0.68161747

5.00000000 0.23289931 -2.00304084

##單元應(yīng)變##

單元號 εx εy γxy

——————————————————

1.00000000 -0.10808734 0.00000000 -0.31990295

2.00000000 0.11226850 -0.04181156 0.09990295

3.00000000 0.00418116 -0.04181156 -0.22000000

##單元應(yīng)力##

單元號 σx σy τxy

——————————————————

1.00000000 -2.18358269 -0.21835827 -2.90820865

2.00000000 2.18358269 -0.61787303 0.90820865

3.00000000 -0.00000000 -0.83623130 -2.00000000

文本中數(shù)據(jù)間使用制表位隔開，可直接復(fù)制到execl中進(jìn)一步操作。

為驗證上述程序的正確性，使用ABQUS自帶的CPS3單元建模驗證，CPS3單元三結(jié)點平面應(yīng)力三角形單元，屬于常應(yīng)變單元。

位移x向和y向位移如圖 3~圖 4所示，節(jié)點位移數(shù)值和自編程序完全相符。1號單元直觀的體現(xiàn)出一階單元內(nèi)部位移線性分布的特征。

應(yīng)變?nèi)?/span>圖 5~圖 6所示，結(jié)果和自編有限元程序一致，圖 5~圖 6值觀反應(yīng)了常應(yīng)變單元的應(yīng)變分布特征。

有用的話請給個三連吧，祝各位同學(xué)高分通過