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

葉康生，孟令寧

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

有限元法[1]自誕生至今已經(jīng)經(jīng)歷了長足的發(fā)展并在科學和工程領(lǐng)域得到了廣泛的應(yīng)用，目前已經(jīng)成為結(jié)構(gòu)工程領(lǐng)域最主要的分析計算方法。該法的解答精度主要依賴于網(wǎng)格和單元次數(shù)，為提高精度，常規(guī)有限元需要加密網(wǎng)格或提高單元次數(shù)，致使自由度快速增加，使得計算成本急劇提高，這在多維問題中尤為顯著。為兼顧有限元法的效率和精度，國內(nèi)外學者在天然超收斂性的基礎(chǔ)上對超收斂后處理進行了大量研究，如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 法的影響最大，該法具有計算簡便，適用性廣等優(yōu)點，但其實施依賴于應(yīng)力超收斂點，對于沒有應(yīng)力超收斂點的單元則難于實施。其他各類方法亦各有局限，研究更加通用、可靠的超收斂后處理方法仍有重要的理論意義和工程應(yīng)用價值。

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

1 模型問題及有限元分析

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

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

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

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

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

圖1 有限元網(wǎng)格Fig.1 FE mesh

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

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

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

相應(yīng)地，全域解近似為：

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

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

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

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

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

其x向、y向?qū)?shù)誤差分別為：

定義有限元導數(shù)誤差為：

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

通常情況下，雙m次矩形元域內(nèi)的網(wǎng)格結(jié)點集 Zh上位移解具有如下超收斂性質(zhì)[21]：

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

網(wǎng)格結(jié)點位移解的超收斂性是本文建立p型超收斂算法的基礎(chǔ)。

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

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

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

圖2 有限元線法網(wǎng)格Fig.2 FEMOL mesh

線法將全域解近似為：

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

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

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

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

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

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

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

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

由于求解線法控制方程時，y向網(wǎng)格可任意布置，故區(qū)間[y1,yny+1]中任一點均可設(shè)為網(wǎng)格結(jié)點，故有：

即線法網(wǎng)格結(jié)線上的解答具有h2m階的超收斂性質(zhì)。

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

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

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

圖3 單元y 向邊位移恢復(fù)(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所示。可得所有單元x向邊界的超收斂解u*(x,yj),j=1,2,···ny+1,x∈Ix，且有：

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

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

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

求解：

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

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

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

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

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

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

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

3 數(shù)值算例

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

記超收斂解的誤差為：

其x向、y向?qū)?shù)誤差分別為：

定義超收斂導數(shù)誤差為：

限于篇幅，例1、例2 導數(shù)收斂階統(tǒng)計中只給出了x向?qū)?shù)的結(jié)果，y向?qū)?shù)具有與x向?qū)?shù)同樣的收斂規(guī)律。

例1.矩形域規(guī)則劃分問題

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

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

該問題的精確解為：

圖8 例1 在(4×4)網(wǎng)格下 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)網(wǎng)格下分布的比較(m=1,=2)Fig.9 Comparison of and distribution on(4×4)mesh in Example 1(m=1,=2)

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

圖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.非規(guī)則四邊形域規(guī)則劃分問題

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

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

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

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

圖12 例2BC網(wǎng)格線上有限元結(jié)點位移誤差(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)網(wǎng)格下分布(m=2,=4)Fig.14 Distribution of on(32×32) mesh in Example 2(m=2,=4)

例3.矩形域非規(guī)則劃分問題

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

圖15 例3 模型的求解域及網(wǎng)格劃分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)

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

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

圖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 結(jié)論

在非規(guī)則域規(guī)則網(wǎng)格劃分下，求解域的內(nèi)部子域仍具有上述的超收斂效果。但在邊界層附近，導數(shù)精度僅能提高一階。

對非規(guī)則網(wǎng)格劃分，規(guī)則域的有限元解答精度明顯優(yōu)于非規(guī)則域。本文方法對規(guī)則域和非規(guī)則域的精度均有明顯提升效果，規(guī)則域的精度提升明顯優(yōu)于非規(guī)則域。

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