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

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

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

1 定解問題及弱解形式

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

(1)

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

(2)

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

(3)

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

(4)

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

2 不同坐標(biāo)系梯度和法方向轉(zhuǎn)換

(4)式在計(jì)算時(shí)，需要用到極坐標(biāo)系的梯度與法方向?qū)?shù)。極坐標(biāo)系的梯度算符有如下關(guān)系：

(5)

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

(6)

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

(7)

可得：

(8)

(9)

3 應(yīng)用Freefem++求解定解問題

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

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

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

4 結(jié)語

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