基于Freefem++的有限元數值計算在計算物理中的應用

(貴州師范大學物理與電子科學學院，貴州 貴陽550001)

計算物理以物理知識和數學方法的綜合應用為主，在理工科類專業具有重要意義，特別是以科學研究為主要培養方向的物理學專業[1]。這門課程的典型特征是數值計算，涉及較多偏微分方程的數值求解，對眾多學生而言相對困難[2-3]。在眾多數值方法中，除了有限差分方法外，有限元方法也得到了廣泛的應用，在科學研究和工程實踐中具有重要意義。有限差分方法以快速計算為特征，精度一般不如有限元，特別是處理格網邊界時，這種弊端非常顯著。有限元方法基于分片插值函數，邊界更接近于真實情況，誤差相對較小。這樣的優點在實際應用中至關重要，特別是工程技術中。近十年來，隨著計算機技術的發展，有限元并行計算的應用使得費時問題得到了解決，有限元技術在工程應用中得到了快速的發展和普及。在計算物理領域，除了上述數值方法外，具有更多優點的邊界元和有限體積方法也得到了快速的發展[4-7]。當前市場上主流數值計算軟件，以有限元為主，因此，有限元在科研和工程技術中的應用至關重要。

作為綜合性高校，培養數值計算方向的學生至關重要，特別是培養有限元專業知識技能的學生更是未來教育發展方向和市場所需。目前主流有限元軟件眾多，如Ansys、MSC、Abaqus、COMSOL等，這些商業閉源軟件在業界具有良好的口碑，可靠性較好。但這些軟件由于相關模塊被封裝，學生難以理解有限元的計算原理，較難對相關模塊作進一步的科研開發，限制了教學和科研的使用。除上述軟件外，也有不少開源的有限元軟件，如Elmer、ParaFEM和Freefem++。這些軟件不但對公眾開放源代碼，而且可靠性較高，相關模塊的說明文檔也較豐富，非常適合于教學和科研。經過長時間的比較測試和應用，發現Freefem++語句非常簡潔，只需要將原方程轉換成弱解形式，使用Freefem++很容易編寫有限元程序，是最接近有限元的“計算語言”[8-9]。但該軟件也有一個弊端，沒有極坐標系的相關程序單元可供使用，僅提供直角坐標系的。為了解決這個問題，基于直角坐標系的程序單元，文章推導了極坐標系梯度與法方向的轉換公式，并應用于具體案例中，以期為Freefem++極坐標系的有限元程序設計和教學提供參考。

1 定解問題及弱解形式

問題求解區域如圖1所示。半徑R=1的1/4圓域Ω，圓心為O，邊界區域依次為Γ1、Γ2和Γ3。該區域的函數u(r,θ)滿足如下的泊松方程：

(1)

其中f表示作用力或源，r和θ分別表示徑向和橫向坐標，Δ表示拉普拉斯算符，有如下關系：

(2)

問題(1)的邊界條件，如圖1所示，滿足關系：

(3)

圖1 半徑R=1的1/4圓域Ω及其邊界，其中O表示圓心

(4)

其中Γ表示求解區域Ω的邊界，?u/?n表示u在法方向的梯度，▽表示梯度算符。

2 不同坐標系梯度和法方向轉換

(4)式在計算時，需要用到極坐標系的梯度與法方向導數。極坐標系的梯度算符有如下關系：

(5)

其中êr和êθ分別表示徑向和橫向單位矢量。(4)式表明諸如(1)式和(3)式的強解問題，可以用(4)式的弱解方程來求解，而且(4)式直接包括了強解問題的邊界條件，因此，(4)式具有很高的實用價值。由于本問題中沒有待求函數在邊界上法方向的梯度，需要對(4)式的邊界積分項進行處理。根據方向梯度的定義，可得：

(6)

依據(3)、(4)和(6)式，可以求出邊界項積分。若(x,y)和(r,θ)分別表示直角坐標和極坐標，(,)和(êr,êθ)分別表示直角坐標系和極坐標系的單位方向矢量，考慮到：

(7)

可得：

(8)

(9)

3 應用Freefem++求解定解問題

其中macroGrad(u)用于定義極坐標系中，待求函數u的梯度分量；macroLap(u,v)表示(4)式中雙線性項u·v；macroNrq用于定義法方向在極坐標系的坐標分量；BC2(u)表示(4)式中邊界積分項的法方向梯度?u/?n;problemlaplace(u,v)用于定義(4)式的弱解問題，其中on(1,u=0)和on(3,u=0)分別表示(3)式中的第一邊界條件。應用圖2所示代碼，數值計算結果如圖3(a)所示。根據數學物理方法[11-12]，可知定解問題(1)和(3)的解析解為u=rcosθsinθ，其結果如圖3(b)所示。數值解減去解析解，其平方值的分布如圖3(c)所示。很顯然，數值結果圖3(a)與解析解圖3(b)的分布一致，其誤差分布除坐標原點外，其它地方的誤差最小，誤差分布也較均勻。坐標原點處的誤差，來源于(1)式中，源f含有r，在圓點處會發散，導致數值解在圓點處不穩定所致。

圖2 基于Freefem++語言格式的有限元程序代碼

圖3 有限元數值解(a)，解析解(b)，數值解與解析解之差的平方(c)

4 結語

計算物理在理工科類專業中具有重要作用。特別是隨著科學技術的發展，有限元軟件得到普及，培養這方面的專業技術人才，既是教育發展方向，也能滿足市場需求，對提升高校競爭力至關重要。目前市場上有很多有限元軟件，大多口碑較好的是閉源商業化的，模塊不對外開放，不利于教學和科學研究的應用。為了克服這個問題，眾多科研機構開發了一系列有限元開源軟件，使得有限元方法在教學和科研中的應用成為可能。在眾多開源軟件中，Freefem++類似于C++，語句簡潔，容易上手，非常適合于教學和科研，在工程實踐中也有相關的應用。但該軟件有一個弊端，沒有極坐標系的相關程序單元供使用，軟件提供的僅有直角坐標系的相關程序單元。針對該問題，本文分析了直角坐標系與極坐標系中梯度與法方向的轉換關系，導出了極坐標系的弱解形式，設計出基于Freefem++的有限元程序。并成功地將該程序應用于1/4圓域的泊松方程中，發現數值計算結果與解析解一致。除坐標原點外，數值解與解析解兩者間較小的誤差，顯示了較好的一致性和穩定性，驗證了本文推導的極坐標系梯度和法方向計算公式的合理性。