何朝葵

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

何朝葵

（河海大學 理學院，江蘇 南京 211100）

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

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

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

1 有限元方法計算格式

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

其中

（7）

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

2 計算格式向量化

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

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

結合高斯積分

圖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è)機房及設備，太陽能集熱板布置條件好。

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

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

39 print（u2）

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

4 應用實例

例1邊值問題

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

例2 混合邊值問題

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

5 結語

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

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

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

[3] 馬竹樵．應用型本科院?！队邢拊椒ā方虒W改革分析[J]．教育現(xiàn)代化，2019，6（15）：51-53．

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

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

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

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

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

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

[10] 章本照，印建安，張宏基．流體力學數(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

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

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