二維泊松方程問題Lagrange 型有限元p 型超收斂算法

葉康生，孟令寧

(清華大學土木工程系，土木工程安全與耐久教育部重點實驗室，北京 100084)

有限元法[1]自誕生至今已經經歷了長足的發展并在科學和工程領域得到了廣泛的應用，目前已經成為結構工程領域最主要的分析計算方法。該法的解答精度主要依賴于網格和單元次數，為提高精度，常規有限元需要加密網格或提高單元次數，致使自由度快速增加，使得計算成本急劇提高，這在多維問題中尤為顯著。為兼顧有限元法的效率和精度，國內外學者在天然超收斂性的基礎上對超收斂后處理進行了大量研究，如Zienkiewicz 和Zhu 的SPR 法[2-3]、Wiberg 的SPRD法[4]、Boroomand 和Zienkiewicz 的REP 法[5]、Ubertini 的RCP 法[6]、Payen 和Bathe 的NPF-based 法[7]及SIP 法[8-9]、袁駟的EEP 法[10-12]等。這些方法中，以SPR 法的影響最大，該法具有計算簡便，適用性廣等優點，但其實施依賴于應力超收斂點，對于沒有應力超收斂點的單元則難于實施。其他各類方法亦各有局限，研究更加通用、可靠的超收斂后處理方法仍有重要的理論意義和工程應用價值。

Douglas 和Dupont[13]對二階常微分方程兩點邊值問題，提出了一種基于單元結點位移超收斂性，在單元內使用更高次基函數近似求解控制方程的超收斂算法。本課題組將該思想拓展應用至桿件自由振動問題[14]、平面曲梁和旋轉殼靜力問題[15-16]、平面曲梁和旋轉殼自由振動問題[17]、一維非線性常微分方程及常微分方程組問題[18]及一維自適應分析[19]中。以上問題均為常微分方程問題，力學中還有大量二維問題，如板、殼的彎曲，彈性力學平面問題等，這些問題往往歸結為求解偏微分方程邊值問題。本文進一步將該思想拓展應用至這類問題的有限元超收斂計算中。由于偏微分方程邊值問題亦有不同種類，為簡單起見，本文以二維泊松方程問題為模型問題探討二維C0型問題有限元的超收斂求解。C0型有限元亦分若干種類，本文僅限探討雙向同次的Lagrange型單元，這類單元具有很好的抗畸變能力[20]和解答性態。對于其他類型的單元，將另文探討。

1 模型問題及有限元分析

二維泊松方程是典型的橢圓型偏微分方程邊值問題，科學和工程中的彈性薄膜問題、滲流問題、靜電場問題、自由扭轉等問題均可由此方程描述。

本文對此方程的有限元超收斂計算開展研究，并以薄膜問題為模型進行論述，方程具體為：

式中：u(x,y)為待求位移；n為邊界外法線方向；Ω為求解域，?Ω 為 Ω的邊界，?Ω=Γu∪Γn，Γu∩Γn=?，Γu和 Γn分別為固定邊界和自由邊界，分別給定Dirichlet 和Neumann 邊界條件。該問題對應的能量泛函為：

式中，()x=?()/?x,()y=?()/?y。

對二維求解域通過網格劃分進行單元離散，以矩形網格剖分為例，如圖1 所示。記x向、y向網格劃分分別為：

圖1 有限元網格Fig.1 FE mesh

則二維有限元網格劃分可簡單表示為：I=Ix×Iy。該劃分將求解域離散為ne=nx×ny個矩形單元。將網格結點集記為Zh，它包含結點(xi,yj)，i=1,2,···,nx+1,j=1,2,···,ny+1。單元尺寸he為單元(e)的最大邊長，網格尺寸為所有單元邊長的最大值，即：

在該網格下，采用雙向m次Lagrange 單元對該問題進行求解，將單元(e)上的真解近似為雙向m次Lagrange 多項式：

式中：Ne為單元形函數矩陣；Δe為單元的結點自由度向量。

相應地，全域解近似為：

式中：Nr(x),r=1,2,···,mnx+1為x向網格Ix上的整體m次形函數；Ns(y),s=1,2,···,mny+1為y向網格Iy上的整體m次形函數；Δh為整體的結點自由度向量。

將式(8)代入式(2)，將整體能量泛函離散為：

式中：Kh為整體剛度矩陣；Ph為整體的結點荷載向量。表達式如下：

對該泛函引入固定邊界Γu上Dirichlet 邊界條件后取極值，并將邊界條件處理后的整體剛度矩陣和結點荷載向量仍記為Kh和Ph，則可得有限元方程：

由該方程可求得 Δh，代入式(8)可得有限元位移解uh。記有限元解誤差為：

其x向、y向導數誤差分別為：

定義有限元導數誤差為：

該有限元解答具有如下收斂性[1]：

通常情況下，雙m次矩形元域內的網格結點集 Zh上位移解具有如下超收斂性質[21]：

對一般四邊形網格，網格結點位移不一定具有上述h2m階的超收斂性質。文獻[22]指出，參數變換為強正規的，即滿足單元接近平行四邊形，對邊彎曲程度相近以及單元過渡時邊的方向角改變量為高階小量等條件，則四邊形等參元在網格結點上仍具有超收斂性質。而在一般邊界條件下，邊界層(邊界附近的1 層～3 層單元)網格結點解答的精度往往較差甚至喪失超收斂性。

網格結點位移解的超收斂性是本文建立p型超收斂算法的基礎。

2 超收斂求解思路及誤差估計

上述有限元求解過程可等效為如下的兩步逐維離散過程：先采用有限元線法[23](以下簡稱線法)進行一維離散(沿x向在Ix網格上作m次元離散)，再對線法控制方程進行一維離散(沿y向在Iy網格上作m次元離散)。

線法網格如圖2 所示，記(el),e=1,2,···,nx為線法單元，Lm(e-1)+2,Lm(e-1)+3,···,Lm(e-1)+m為(el)內部結線；Lm(i-1)+1,i=1,2,···,nx+1為x=xi處網格結線，dr(y),r=1,2,···,mnx+1為結線位移函數。

圖2 有限元線法網格Fig.2 FEMOL mesh

線法將全域解近似為：

式中：Nl(x)為線法整體形函數矩陣；d(y)為結線位移向量。

對式(17)引入固定邊界 Γu上的位移條件，并代入式(2)中能量泛函，通過泛函取極值，可得線法控制方程：

式中：()′=d()/dy；A、B、C、F均為y的函數。其表達式及方程的邊界條件可參見文[23]。

上述線法控制方程為常微分方程邊值問題，對該問題進行一維有限元離散(沿y向在Iy網格上作m次元離散)，將結線位移離散為：

如此經兩步逐維離散過程后的位移表達式為：

易見矩形網格分兩步逐維離散與直接進行二維Lagrange 單元離散可得到相同的離散格式，兩者代入泛函取極值可得到相同的控制方程及解答，故式(20)與式(8)中的相同，未作區分。并有如下等效關系：

因第二步離散仍為一維有限元過程，則由文獻[1]可得，結線上網格結點解有如下超收斂性質：

而對于x=xi,1≤i≤nx+1處網格結線Lm(i-1)+1上的網格結點，由式(16)、式(21)、式(22)可得如下估計式：

由于求解線法控制方程時，y向網格可任意布置，故區間[y1,yny+1]中任一點均可設為網格結點，故有：

即線法網格結線上的解答具有h2m階的超收斂性質。

為提高二維有限元單元y向邊界解答的精度和收斂階，本文基于式(22)中結線上網格結點解的超收斂性，對線法控制方程采用p型超收斂進行修復，過程如下：

線法控制方程在y向網格上的真解d，應滿足如下表達：

上述p型超收斂過程等效于對二維網格劃分中的第j行子域采用次單元進行有限元求解，其中在邊界取為原有限元解，其余邊界與原邊界條件一致，如圖3 所示。Ix×Iyjm×y=yj,y=yj+1

圖3 單元y 向邊位移恢復(m=2,=3)Fig.3 Recovery of displacements on y-direction edges of elements(m=2,=3)

依次取所有ny個行子域Ix×Iyj,j=1,2,···,ny，即可求得所有單元y向邊界的超收斂解u*(xi,y),i=1,2,···,nx+1,y∈Iy。

同理，亦可將上述逐維離散過程視為先對y向進行線法離散，后對線法控制方程在x向進行一維有限元離散。類似地，對第i列子域Ixi×Iy,i=1,2,···,nx采用×m次單元對原泊松方程問題進行有限元求解，其中在x=xi,x=xi+1邊界取為原有限元解，其余邊界與原邊界條件一致，如圖4所示?？傻盟袉卧獂向邊界的超收斂解u*(x,yj),j=1,2,···ny+1,x∈Ix，且有：

圖4 單元x 向邊位移恢復(m=2,=3)Fig.4 Recovery of displacements on x-direction edges of elements (m=2,=3)

在求得全部單元邊界(網格線)上的超收斂解后，取出二維網格劃分中的單元域1,2,···,nx,j=1,2,···,ny，采用雙次單元進行二維有限元求解，其邊界解取為前述邊界超收斂解，如圖5 所示，即用：

圖5 單元域位移恢復(m=2,=3)Fig.5 Recovery of displacements on elements (m=2,=3)

求解：

由于邊界自由度均被固定，故此步只需求解單元內部自由度。逐單元進行上述求解，即可得全域超收斂解：

由二維有限元解答的收斂性估計[1]知：

結合式(38)、式(39)可知超收斂解誤差為：

數值結果表明，u*(x,y)的導數?u*(x,y)/?x、?u*(x,y)/?y亦具有如下的超收斂性：

上述精度修復策略亦可推廣至任意四邊形網格。具體實施時亦是先修復單元邊界，再修復單元內部。修復單元邊界時，同一單元的相對邊同時進行修復，為使相鄰單元在公共邊上的修復解相同，有公共邊的單元同時參與修復計算，所有參與對邊修復計算的單元構成一子網格，在該子網格上用m×次單元求解原控制方程，對對邊解答進行修復(對邊方向取為次，剩余對邊方向保留m次)，并將單元剩余對邊的解取為原有限元解，子網格在域邊界上的邊界條件取為與原問題一致，如圖6 所示。單元內部解答修復過程與矩形網格相同。

圖6 非規則網格四邊形單元邊位移恢復(m=2,=3)Fig.6 Recovery of displacements on edges of quadrilateral elements on irregular mesh(m=2,=3)

對于任意四邊形網格，由于受單元畸變[20]的影響，網格結點的位移精度有所下降，不一定具有h2m階的超收斂性，但相對于單元內部位移，網格結點位移仍具有更高的精度。應用本文方法，仍可有效提高解答的精度和質量，使位移尤其是應力的精度得到顯著提升。

3 數值算例

本文算法已編成Fortran 程序。本節給出三個數值算例以展示其在不同求解區域、不同網格劃分條件下對有限元解答質量和精度的恢復效果。對全域或子域誤差均以最大模度量：

記超收斂解的誤差為：

其x向、y向導數誤差分別為：

定義超收斂導數誤差為：

限于篇幅，例1、例2 導數收斂階統計中只給出了x向導數的結果，y向導數具有與x向導數同樣的收斂規律。

例1.矩形域規則劃分問題

本例探討矩形域規則劃分解泊松方程問題。如圖7所示，其控制微分方程及邊界條件為：

圖7 例1 問題模型Fig.7 Problem model in Example 1

該問題的精確解為：

圖8 例1 在(4×4)網格下 eh 與 e*分布的比較(m=1,=2)Fig.8 Comparison of ehand e*distribution on(4×4)mesh in Example 1(m=1,=2)

圖9 例1 在(4×4)網格下分布的比較(m=1,=2)Fig.9 Comparison of and distribution on(4×4)mesh in Example 1(m=1,=2)

將求解域均勻剖分為(2×2)初始網格，并進行均勻二分加密，考察并統計(m=3,=4)、(m=3,=5)、(m=3,=6)、(m=3,=7)情形下本算法位移及導數的收斂階規律，結果如表1 和圖10 所示。經超收斂計算后，位移解和導數解的精度均得到顯著提高，解答收斂階隨超收斂次數的提高而提高，但都不超過網格結點位移解的收斂階h2m。收斂階規律與式(40)、式(41)中的估計吻合。表2 列出了(m=3,=4)情形下有限元求解及超收斂計算的時間對比，可見隨有限元求解規模的增大，超收斂計算的耗時占比快速下降，充分體現了本文算法的修復過程計算量小、效率高。

圖10 例1 有限元與超收斂解收斂階Fig.10 Convergence order of FE solution and recovered solution in Example 1

表1 例1 有限元與超收斂解的誤差和收斂階Table 1 Error and convergence order of FE solution and recovered solution in Example 1

表2 例1 有限元與超收斂過程計算時間(m=3,=4)Table 2 Computation time of FE solution and recovered solution in Example 1(m=3,=4)

表2 例1 有限元與超收斂過程計算時間(m=3,=4)Table 2 Computation time of FE solution and recovered solution in Example 1(m=3,=4)

例2.非規則四邊形域規則劃分問題

本例探討非規則四邊形域規則劃分解泊松方程問題。其求解域 Ω如圖11 所示，其控制微分方程及邊界條件為：

圖11 例2 模型的求解域及網格劃分Fig.11 Domain and mesh of Example 2

由于在一般邊界條件下，邊界層網格結點位移解的天然超收斂性相對較差[22]，這會導致該區域修復解的精度相對差些。為檢視本文方法對域內解答的修復效果，本文取一典型內子域 Ωλ(由均勻4×4 網格劃分中內部4 個單元構成的子域)考察解答誤差的變化情況，如圖11 所示。為探討四邊形域規則網格下有限元網格結點位移解的收斂階，采用雙二次元(m=2)、雙三次元(m=3)、雙四次元(m=4)進行有限元計算，提取(4×4)網格結點A(4.375,1.1875)位移誤差值邊界層(始終靠近邊界的一層單元)網格結點最大位移誤差值，統計其收斂階具體見表3。圖12 展示了在(64×64)規則網格下(m=3)情形下BC 網格線上有限元結點位移誤差的分布情況?？梢?，域內部網格結點位移收斂階可達最佳收斂階次h2m；而受非規則區域邊界的影響，靠近邊界的一層單元角結點位移收斂階僅為hmin(2m,m+2)。

表3 例2 有限元網格結點位移誤差和收斂階Table 3 Error and convergence order of mesh node displacement of FE solution in Example 2

圖12 例2BC網格線上有限元結點位移誤差(Fig.12 ErrorofFEnoded√isplacements on BC mesh line inExample

表4 例2 有限元和超收斂解收斂階Table 4 Convergence order of FE solution and recovered solution in Example 2

圖13 例2 有限元與超收斂解收斂階Fig.13 Convergence order of FE solution and recovered solution in Example 2

圖14 例2 在(32×32)網格下分布(m=2,=4)Fig.14 Distribution of on(32×32) mesh in Example 2(m=2,=4)

例3.矩形域非規則劃分問題

本例探討矩形域非規則四邊形網格劃分解泊松方程問題。采用HyperMesh 軟件[24]作網格劃分，得圖15 所示網格，記該網格劃分為J 剖分，對該網格進行均勻二分加密得到的網格劃分記為K 剖分。本例控制微分方程及邊界條件為：

圖15 例3 模型的求解域及網格劃分Fig.15 Domain and mesh of Example 3

表5 例2 有限元與超收斂解的誤差和收斂階(m=2)Table 5 Error and convergence order of FE solution and recovered solution in Example 2(m=2)

圖16 例3 在J 剖分下 |eh| 與 |e*|分 布的比較(m=2,=3)Fig.16 Comparison of |eh| a nd |e*| distribution on J mesh in Example 3(m=2,=3)

在最大模度量下，可見對上述問題，均可將全域位移(包括非規則區域)誤差恢復至網格結點水平。對于(m=2,=3)情況J 剖分下規則區域位移誤差比為3.00%，K 剖分下該誤差比進一步降低至1.74%。可見隨網格加密，規則區域修復的位移呈加速收斂之勢，即經超收斂恢復后，規則域位移收斂階可得到提高。因非規則劃分不滿足變換的強正規條件，故非規則區域網格結點位移往往不具有超收斂性質[18]。但該區域網格結點位移精度仍優于單元內部，本文方法可將單元域內位移恢復至網格結點水平，故非規則域位移誤差仍有顯著下降，全域位移誤差比(由非規則域控制)為40%左右，部分非規則區域位移誤差比為20%左右。

圖17 例3 在J 剖分下分布的比較(m=2,=3)Fig.17 Comparison of distribution on J mesh in Example 3(m=2,=3)

表6 例3 有限元和超收斂誤差(m=1)Table 6 Error of FE solution and recovered solution in Example 3(m=1)

表7 例3 有限元和超收斂誤差(m=2)Table 7 Error of FE solution and recovered solution in Example 3(m=2)

綜上，對于非規則網格問題，雖然有限元的求解精度較差，本文算法對解答精度仍有明顯的提升效果。

圖18 例3在J 剖分下(m=2,=3) 與分布的比較Fig.18 Comparison of(m=2,=3)anddistribution on J mesh in Example 3

表8 例3 本文算法與SPR 法的誤差對比Table 8 Error comparison of this method and SPR in Example 3

4 結論

在非規則域規則網格劃分下，求解域的內部子域仍具有上述的超收斂效果。但在邊界層附近，導數精度僅能提高一階。

對非規則網格劃分，規則域的有限元解答精度明顯優于非規則域。本文方法對規則域和非規則域的精度均有明顯提升效果，規則域的精度提升明顯優于非規則域。

本文方法通過在一系列子網格上作有限元求解獲取超收斂解，計算量小，而解答的精度和質量均有顯著提高，且方法易于推廣。故本文算法具有進一步研究的價值。