有限元實驗課程教學探索

黃鵬展

（新疆大學數學與系統科學學院,烏魯木齊 830046）

眾所周知,大量實際中的工程問題都可以用偏微分方程來描述,但有許多偏微分方程定解問題的解不能以實用的解析形式來表出,從而無法得到這些方程的準確解,以定量地描述客觀過程。有限元課程是數值分析的后繼課程,是計算數學專業學生必修的專業課程。課程中的有限元方法在數值計算的理論基礎上,構造求解的近似逼近方法。該課程通過一些典型、常用、有效的數值方法來解釋有限元方法的基本思想,使學生了解如何在計算機上應用這些數值方法求解一個偏微分方程定解問題。同時,通過對有限元基本概念及基本理論的學習,使學生的理論分析能力得到一定的訓練,并通過學習與實踐使得學生熟練掌握及運用有限元方法來數值求解偏微分方程。

而有限元實驗課程的設立,則能更好地為學生學習有限元服務,使得學生能夠更好、更深刻地了解有限元方法,體會有限元方法在數值求解偏微分方程時的魅力。本文主要結合筆者的親身教學經歷,從有限元的Matlab工具包、有限元軟件Free?fem++以及Matlab編程等三個方面來探索有限元實驗課程的教學方法。這三個方面從簡到難、從特定偏微分方程類型求解到較一般偏微分方程求解,循序漸進地引導學生進行有限元方法的實踐和數值的實驗。

1 有限元方法

本文以文獻〔1〕中的一個偏微分方程為例,簡單地介紹基本的有限元方法。

考慮二維區域上的Poisson方程〔2〕第三邊值問題

這里的Ω∈R2是具有光滑邊界的有界凸區域,n是?Ω上的外法線方向,α和g是?Ω上的連續函數,α=α(x,y)≥0,△是拉普拉斯算子,f是作用力。為簡便起見,在后面的計算部分,令 f=g=α=1,計算區域Ω取成單位正方形。

我們知道有限元方法基于變分原理,故先給出方程（1）的變分形式。即第一步,

第二步,給出計算區域的網格剖分。設h是一個趨近于0的正參數,記網格剖分為Kh,它由一些三角單元K組成。

第三步,試探函數空間的構造。即選取有限元子空間

其中, φi（i=1,2,…,N）是基函數,滿足

求使得，D(uh,vh)=F(vh),?vh∈ Xh，

其中ui=u(xi,yi)。

第四步,有限元離散。生成單元剛度矩陣和單元荷載向量,通過單元上局部編碼與整體編碼的對應關系得到的擴張矩陣,生成總剛度矩陣和總荷載向量,然后對約束條件進行處理,最后得到有限元方程組。由方程組的求解方法解出ui,結合基函數φi得到有限元近似解uh。

2 有限元上機實驗

有限元課程是一門計算數學專業中較難的課程,故而,我們首先對問題（1）使用較為簡單Matlab工具包來實現數值模擬,鑒于它的可視化操作,學生易于上手。從而,增強他們對這門較難課程學習的信心,克服學習前的恐懼心理和由此產生的怠學狀態。通過Matlab工具包對特定方程的模擬,能夠使學生了解所求得方程的解到底是什么樣子,初步了解有限元解法的過程。接著,采用有限元軟件Freefem++來數值求解問題（1）,該有限元軟件只需要寫出有限元的變分形式即可。最后,按照有限元編程的框架圖,使用Matlab來編制有限元求解（1）的程序。

2.1 Matlab工具包 根據文獻〔3〕中的講解,筆者以（1）為例,對Matlab工具包求解偏微分方程進行介紹。

第一步,打開Matlab軟件,在命令框中輸入pdetool,便可啟動一個可視化界面。接著,完成求解區域的構造。

第二步,選取邊界,打開Boundary Condition對話框,輸入（1）中的第三類邊界條件。

第三步,選擇PDE Mode命令,設置方程的類型,這里為橢圓型方程,接著輸入（1）中偏微分方程的系數。

第四步,利用Initialize Mesh命令,對第一步構造的求解區域進行網格剖分。如果網格不夠密,可以利用Refine Mesh命令進行加密。

第五步,選擇Solve PDE命令,求解設定的偏微分方程適定問題,給出數值解uh的圖形,見圖1。

圖1 用Matlab工具pdetool得到（1）的有限元解

Matlab工具包可以對空間二維問題快速高效求解。由于它的界面可視化,學生只需要使用用戶界面,構造出求解區域,輸入方程的類型和邊界條件,自動生成網格,不需要人工設定基函數,就可以顯示出數值解的結果。故而,它是一種利用有限元方法數值求解偏微分方程的簡單有效工具。但是這個可視化界面只是對特定的一些方程,并且邊界條件的設置也很固定。最主要的是不能夠自己選擇基函數或者構造基函數,而有限元的魅力主要正是在于它所擁有的各式各樣的基函數,并可以根據自己的需要進行構造。

故而,Matlab工具包只適合有限元方法的初學者,使他們較好較快地了解求解偏微分方程的有限元方法。在這方面,自然地它有不可替代的優勢。因此,筆者在有限元實驗課程教學時,最初幾節課,主要就是給學生們講解Matlab工具包求解偏微分方程。

2.2 有限元軟件 FreeFem++ FreeFem++〔4〕是一款數值求解偏微分方程的有限元軟件。就像它名字所說的,它是一款免費的軟件,可以直接在http://www.freefem.org這個網站上下載安裝包,并自帶它的教程。FreeFem++不是一個工具包,而是一個完整的有限元產品,具有內存需求少、計算速度快、易學好用、功能強大、應用面廣等特點〔5-8〕。故而,在有限元實驗課程隨后的課時中,筆者將講授這部分的內容。

仍以（1）為例,對有限元軟件FreeFem++求解偏微分方程進行介紹。以下部分就是用該軟件實現數值求解的FreeFem++程序,附程序注釋。

real n=80;//設置網格尺度

mesh Th=square（n,n）;//生成網格

fespace Xh（Th,P2）;//定義有限元空間

Xh u,v;//聲明試探函數和檢驗函數

func f=1;//選取作用力函數

solve Poisson（u,v）=//定義偏微分方程,即變分形式

·馬改戶 曹春生 何鄂 陳云崗 黎明 郭北平楊參軍 白羽平 冷軍 朱春林等10位雕塑、油畫大家力挺學術邀請展

int2d（Th）（dx（u）*dx（v）+dy（u）*dy（v））//雙線性項形式

+int1d（Th）（u*v）//D（u,v）的第二部分,系數為1

-int2d（Th）（v）//右端項

-int1d（Th）（v）;//F（v）的第二部分,g取1

plot（u,ps=“fig2.eps”）;//計算結果的可視化這樣就完成了采用FreeFem++軟件求解（1）的全部過程,數值結果見圖2。

圖2 利用FreeFem++軟件得到（1）的有限元解

由以上的程序我們可以看到,使用該軟件進行編程數值求解偏微分方程時,只需要寫出原問題的變分問題即可,無需進行單元局部和總體的編號、單元剛度矩陣和單元荷載向量的計算、總剛度矩陣和總荷載向量的組成、有限元線性方程組的求解等具體過程。更加主要的是采用“fespace Xh（Th,P2）”這段語句就完成了有限元空間的定義,并且使用“int2d（Th）”這個命令就實現了數值積分的功能。寥寥數行即完成了對Poisson方程第三類邊值問題的有限元方法求解。對比使用C或者Matlab編程需幾百行的代碼,FreeFem++程序大大節約了編程時間。

但是,由于它是一款免費軟件,其底層語言沒有公開,繼而我們無法了解最本質的內容。故而,它有一定的局限性,無法對底層的代碼進行修改,二次編程。并且,只能選取它所給定的一些基函數和有限元空間。因此,它適合一些想使用有限元方法做一些初步研究的研究生。如果想要驗證更加原創的一些有限元算法,筆者認為應該利用C或Matlab或者Fortran語言來實現。

2.3基于Matlab的有限元編程 鑒于以上考慮,筆者在這門課程的剩余課時里,重點講解基于Matlab的有限元方法編程。主要是根據文獻〔9-10〕中的程序框圖進行程序編制,仍舊以（1）為例。

M文件1:構造基函數 function J=Base（m,n），m,n分別是x,y軸等分點數。

M文件2:構造矩陣function A=juzhen（n,m），即構造求解未知節點所滿足的矩陣。這里為重點部分,矩陣由總剛度矩陣并利用數值積分得到。

M文件3:當f=1時矩陣的右端項function b=youduanxiang（n,m）。

M文件4:將求出的矩陣變換為一般的矩陣function B=bh（n,m），即將矩陣A轉換成與理論分析時計算的矩陣一樣。

M文件5:右端項為1時,用CG方法求解矩陣function u=gonge（n,m,detla,max），這里,max為最大循環次數,detla=1.0E-6為容許誤差。

M文件6:右端項為1時,用有限元方法得到的解 U=qiujie（n,m,detla,max）。

再利用一個主文件調用以上的M文件,給定所設的參數,就可以得到（1）的有限元方法數值解,見圖3。

圖3 利用Matlab編程得到（1）的有限元解

故而,利用文獻〔9〕、〔10〕里面的程序框架圖,我們可以編制有限元的Matlab程序,從而得到想要的數值結果。直接編程的主要優點就是算法驗證的靈活性,可以任意構造希望的基函數及有限元空間。但是它所花費的時間最長,因編制程序以及調試程序等都得消耗大量的時間。

3 結論

本文結合筆者的教學經歷,對有限元實驗課程教學進行了探索。針對數值求解偏微分方程,主要關于Matlab工具包、有限元軟件Freefem++以及Matlab有限元編程這三個方面展開討論。具體介紹了這三種方法的使用和編制,給出了每種方法的優缺點。

總之,對有限元實驗課程的教學要循序漸進。Matlab工具包最為簡單適合初學者使用,但是對于求解的問題具有局限性;有限元軟件Freefem++次之,但是對基函數選取與構造等方面的靈活性欠缺,故它適合成為研究生用來做一些初步研究的有效工具;Matlab有限元編程相比之下最為困難與耗時,但它是驗證一些有限元原創算法的必備工具。

〔1〕陸金甫,關治.偏微分方程數值解〔M〕.2版.北京:清華大學出版社,2004.

〔2〕劉相國,謝如龍,郝江鋒.二維泊松方程的交替方向迭代法〔J〕.大理學院學報,2010,9（10）:1-5.

〔3〕陸君安,尚濤,謝進,等.偏微分方程的MATLAB解法〔M〕.武漢:武漢大學出版社,2001.

〔4〕Hecht F.FreeFem++〔EB/OL〕.（2012-03-16）〔2012-10-10〕.http://www.freefem.org/ff++.

〔5〕尚月強.FreeFem++與偏微分方程有限元數值計算〔J〕.貴州師范學院學報,2012,27（12）:1-5.

〔6〕尚月強.基于MPI+FreeFem++的有限元并行計算〔J〕.軟件導刊,2012,11（2）:27-29.

〔7〕杜世森.基于FreeFEM++的一階有限元解法〔D〕.長沙:湖南大學,2011.

〔8〕詹陽烈,莊春剛，熊振華，等.基于水平集方法與Free?FEM的彈性結構拓撲優化〔J〕.中國機械工程,2009,20（11）:1326-1331.

〔9〕李開泰,黃艾香,黃慶懷.有限元方法及其應用〔M〕.北京:科學出版社,2006.

〔10〕馬逸塵,梅立泉,王阿霞.偏微分方程現代數值方法〔M〕.北京:科學出版社,2006.