何朝葵

有限元方法計算格式的向量化

何朝葵

（河海大學(xué) 理學(xué)院，江蘇 南京 211100）

針對初學(xué)者反映有限元方法難以編制程序的問題，基于Python軟件的向量運算優(yōu)勢，引入形函數(shù)向量和形函數(shù)導(dǎo)數(shù)向量，從而使有限元方法的計算格式向量化．向量化后的計算格式簡單易懂，易于編寫程序，也便于根據(jù)精度要求和邊界條件的變換進行調(diào)整.

有限元方法；高斯積分；向量化；程序

有限元方法是一種非常流行的數(shù)值計算方法，廣泛應(yīng)用于結(jié)構(gòu)力學(xué)、流體力學(xué)和熱學(xué)等眾多領(lǐng)域．掌握有限元方法已成為很多工程設(shè)計從業(yè)人員的一個基本要求．有限元方法由于其具有理論性強、基礎(chǔ)要求高、實踐應(yīng)用性特別強的特點，給有限元方法課程教學(xué)的內(nèi)容安排、教學(xué)方案設(shè)計帶來了很大的挑戰(zhàn)，學(xué)生普遍反映學(xué)習的難度大，為此，眾多學(xué)者對有限元方法課程如何選擇教材，如何安排教學(xué)內(nèi)容，采用何種有效教學(xué)手段進行了深入的探討[1-8]．有限元方法課程一個關(guān)鍵的難點在于有限元方法的計算格式比較繁雜，導(dǎo)致程序的編寫比較困難，而學(xué)生判斷“懂了”和“沒懂”的關(guān)鍵指標就是程序編制是否成功，一旦程序調(diào)試成功，他們才認為掌握了有限元這種方法．在有限元方法課程教學(xué)過程中，教師不妨把教會學(xué)生編制有限元方法程序（一次插值和二次插值）作為首要任務(wù)，再思索如何簡明地講解有限元方法計算格式，引導(dǎo)他們編制相應(yīng)的程序．

1 有限元方法計算格式

有限元方法的第一步，基于加權(quán)余量原理把問題的微分形式轉(zhuǎn)化為積分形式．如對于兩點邊值問題

其中

（7）

按照有限元方法這個計算格式的流程，把關(guān)鍵點放在式（6）（7）上，采用循環(huán)語句依次對每個單元計算單元的系數(shù)矩陣和常數(shù)項矩陣里的每一個元素，并且積分還要采用高斯積分，這樣就需要多個循環(huán)語句．另外，由于3個形函數(shù)和形函數(shù)的導(dǎo)數(shù)的結(jié)構(gòu)是不同的，程序編寫起來就會很冗長，對于初學(xué)者來說，編制程序就變得非常困難．第一次講授時，20個學(xué)生中只有3個學(xué)生編制程序成功．這嚴重地影響了教學(xué)的效果．因此，如何優(yōu)化有限元方法的計算格式，使其便于初學(xué)者理解，并易于編制程序就顯得非常重要．

2 計算格式向量化

考慮到Python軟件的向量運算功能，對有限元方法計算格式采用向量形式來表示，引入形函數(shù)向量和形函數(shù)導(dǎo)數(shù)向量

那么，單元系數(shù)矩陣和單元常數(shù)項矩陣就可分別表示成

結(jié)合高斯積分

圖1　有限元方法程序流程

3 有限元方法程序

1 import numpy as np

2 x=[0，0.45，1，1.5，2]

3 n=len（x）

4 m=2；n1=m*（n-1）+1；

5 xi=[0.2113，0.7887]

6 w=[0.5，0.5]

7 n2=len（w）

8 c1=1；c2=8

9 A=np.zeros（（n1，n1））；h=np.zeros（n-1）

10 b=np.zeros（（n1，1））

11 for i in range（n-1）：

12 h[i]=x[i+1]-x[i]

13 Ai=np.zeros（（m+1，m+1））

14 bi=np.zeros（（m+1，1））

15 for j in range（n2）：

16 p=x[i]+h[i]*xi[j]+1

17 q=1

18 f=x[i]+h[i]*xi[j]

19 N=[[（2*xi[j]-1）*（xi[j]-1），4*xi[j]*（1-xi[j]），xi[j]*（2*xi[j]-1）]]

20 dN=[[4*xi[j]-3，4-8*xi[j]，4*xi[j]-1]]

21 NT=np.transpose（N）； dNT=np.transpose（dN）

22 Ai=Ai+w[j]*（p/h[i]*np.dot（dNT，dN）+h[i]*q*np.dot（NT，N））

23 bi=bi+w[j]*f*NT*h[j]

24 A[m*i：m*（i+1）+1，m*i：m*（i+1）+1]=A[m*i：m*（i+1）+1，m*i：m*（i+1）+1]+Ai

25 b[m*i：m*（i+1）+1]=b[m*i：m*（i+1）+1]+bi[0：m+1]

26 A[0，：]=0.0；A[0，0]=1；b[0]=c1

27 A[n1-1，：]=0.0；A[n1-1，n1-1]=1；b[n1-1]=c2

28 invA=np.linalg.inv（A）

29 u=np.zeros（（n1，1））

30 u=np.dot（invA，b）

31 y = np.zeros（（n1，1））

32 for i in range（n-1）：

33 y[m*i]=x[i]

34 for k in range（m-1）：

35 y[m*i+1+k]=x[i]+（k+1）*h[i]/m

36 y[n1-1]=x[n-1]

鹽城地區(qū)年平均日太陽輻照量達13 582kJ/m2，全年大部分時間太陽光照充分，6層商業(yè)裙房屋頂面積大，除放置機電各專業(yè)機房及設(shè)備，太陽能集熱板布置條件好。

37 u2 = np.zeros（（n1， 2））

38 u2[：，0]=y[：，0]； u2[：，1]=u[：，0]

39 print（u2）

程序的第5，6行是高斯積分點以及對應(yīng)的權(quán)系數(shù)，計算單元系數(shù)矩陣和單元常數(shù)項矩陣只用了2行代碼，即第22，23行．第19，20行分別為形函數(shù)向量、形函數(shù)導(dǎo)數(shù)向量，如果要采用更高次的拉格朗日插值函數(shù)，則修改第19，20行，并相應(yīng)地提高數(shù)值積分的精度，修改程序的第5，6行，以及第4行中的m即可．程序具有良好的可讀性，根據(jù)精度需求和不同的邊界條件易于做出相應(yīng)的調(diào)整.

4 應(yīng)用實例

例1邊值問題

圖2　例1的有限元方法數(shù)值解

例2 混合邊值問題

由邊界條件有，因此，在對例1采用三次插值函數(shù)計算的基礎(chǔ)上再修改程序中第8，27行：8 c1=0；g=1；27 b[n1]= b[n1]+1．運行程序，得到數(shù)值解（見圖３）．

5 結(jié)語

有限元方法程序主要步驟在于依次對每個單元計算單元系數(shù)矩陣和單元常數(shù)項矩陣，經(jīng)典的有限元方法教材編制程序用的是Fortran語言，通過多個循環(huán)語句來計算單元系數(shù)矩陣和單元常數(shù)項矩陣的元素，這可以理解為從微觀的角度來編制程序．借助Python語言矢量運算的優(yōu)勢，引入形函數(shù)向量和形函數(shù)導(dǎo)數(shù)向量，那么單元系數(shù)矩陣和單元常數(shù)項矩陣就是通過向量運算得到的，即是從宏觀的角度來理解有限元方法的計算格式，從宏觀的角度來編制程序，程序不但具有良好的可讀性，還便于根據(jù)精度要求、邊界條件的變換進行調(diào)整，這顯然達到了優(yōu)化程序的目的．由于篇幅的限制，本文是以一維問題為例來說明向量化的優(yōu)勢，實際上對于二維、三維問題這一優(yōu)勢更加明顯．把計算格式向量化后，在教學(xué)實踐中也看到了效果，編制程序成功的學(xué)生越來越多，選這門課的學(xué)生也越來越多．

[1] 劉巖．專業(yè)認證背景下有限元技術(shù)基礎(chǔ)課程改革與實踐[J]．高師理科學(xué)刊，2019，39（4）：76-79．

[2] 周騰，吳之豪，史留勇．有限元分析及混合教學(xué)法在材料力學(xué)課程中的教學(xué)研究[J]．高教學(xué)刊，2021（35）：82-86，90．

[3] 馬竹樵．應(yīng)用型本科院校《有限元方法》教學(xué)改革分析[J]．教育現(xiàn)代化，2019，6（15）：51-53．

[4] 劉楊，趙宇來，奚方權(quán)，等．有限單元法課程教學(xué)現(xiàn)狀與改革方法探究[J]．中國教育技術(shù)裝備，2019（18）：71-73．

[5] 于亞婷，杜平安．《有限元法》課程實踐教學(xué)方法探索[J]．實驗科學(xué)與技術(shù)，2008（1）：108-110．

[6] 智晉寧．有限元方法課程教學(xué)改革與實踐[J]．安徽工業(yè)大學(xué)學(xué)報（社會科學(xué)版），2010，27（5）：117-118．

[7] 張有宏，常新龍，張青，等．基于“理論引領(lǐng)、應(yīng)用并行、案例示范”理念的有限元方法教學(xué)改革與實踐[J]．教育教學(xué)論壇，2020（6）：127-128．

[8] 向家偉．“有限元方法及程序設(shè)計”課程教學(xué)實踐[J]．重慶工學(xué)院學(xué)報（自然科學(xué)版），2007（7）：171-173．

[9] 林群．微分方程數(shù)值解法基礎(chǔ)教程[M]．3版．北京：科學(xué)出版社，2017．

[10] 章本照，印建安，張宏基．流體力學(xué)數(shù)值方法[M]．北京：機械工業(yè)出版社，2003．

Vectorization of calculation format of finite element method

HE Zhaokui

（School of Science，Hohai University，Nanjing 211100，China）

Aiming at the problem that beginners reflect that the finite element method is difficult to program，based on the vector operation advantages of the software Python，the shape function vector and the shape function derivative vector are introduced，so that the calculation format of the finite element method is vectorized．The calculation format after vectorization is simple and easy to understand，easy to program，and easy to adjust according to the accuracy requirements and the transformation of boundary conditions.

finite element method；Gaussian integral；vectorization；program

1007-9831（2022）09-0062-05

O242.21∶G642.0

A

10.3969/j.issn.1007-9831.2022.09.013

2022-01-10

高等學(xué)校大學(xué)數(shù)學(xué)教學(xué)研究與發(fā)展中心項目（CMC20220203）；江蘇省高等教育教改研究課題重點項目（2021JSJG090）；2021年江蘇省高?！盎A(chǔ)課程群”專項重點課題（2021JDKT017）；河海大學(xué)課程思政示范課程建設(shè)項目（2021A05）

何朝葵（1974-），男，湖南郴州人，講師，碩士，從事有限元方法教學(xué)研究．E-mail：hzk@hhu.edu.cn