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

葉康生，孟令寧

(清華大學(xué)土木工程系，土木工程安全與耐久教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

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

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

1 模型問題及有限元分析

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

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

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

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

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

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

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

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

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

相應(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é)點(diǎn)自由度向量。

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

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

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

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

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

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

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

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

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

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

2 超收斂求解思路及誤差估計(jì)

上述有限元求解過程可等效為如下的兩步逐維離散過程：先采用有限元線法[23](以下簡稱線法)進(jìn)行一維離散(沿x向在Ix網(wǎng)格上作m次元離散)，再對線法控制方程進(jìn)行一維離散(沿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ù)。其表達(dá)式及方程的邊界條件可參見文[23]。

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

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

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

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

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

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

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

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

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

上述p型超收斂過程等效于對二維網(wǎng)格劃分中的第j行子域采用次單元進(jìn)行有限元求解，其中在邊界取為原有限元解，其余邊界與原邊界條件一致，如圖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個(gè)行子域Ix×Iyj,j=1,2,···,ny，即可求得所有單元y向邊界的超收斂解u*(xi,y),i=1,2,···,nx+1,y∈Iy。

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

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

求解：

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

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

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

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

上述精度修復(fù)策略亦可推廣至任意四邊形網(wǎng)格。具體實(shí)施時(shí)亦是先修復(fù)單元邊界，再修復(fù)單元內(nèi)部。修復(fù)單元邊界時(shí)，同一單元的相對邊同時(shí)進(jìn)行修復(fù)，為使相鄰單元在公共邊上的修復(fù)解相同，有公共邊的單元同時(shí)參與修復(fù)計(jì)算，所有參與對邊修復(fù)計(jì)算的單元構(gòu)成一子網(wǎng)格，在該子網(wǎng)格上用m×次單元求解原控制方程，對對邊解答進(jìn)行修復(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é)點(diǎn)的位移精度有所下降，不一定具有h2m階的超收斂性，但相對于單元內(nèi)部位移，網(wǎng)格結(jié)點(diǎn)位移仍具有更高的精度。應(yīng)用本文方法，仍可有效提高解答的精度和質(zhì)量，使位移尤其是應(yīng)力的精度得到顯著提升。

3 數(shù)值算例

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

記超收斂解的誤差為：

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

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

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

圖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 有限元與超收斂過程計(jì)算時(shí)間(m=3,=4)Table 2 Computation time of FE solution and recovered solution in Example 1(m=3,=4)

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

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

圖12 例2BC網(wǎng)格線上有限元結(jié)點(diǎn)位移誤差(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)格進(jìn)行均勻二分加密得到的網(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é)點(diǎn)水平。對于(m=2,=3)情況J 剖分下規(guī)則區(qū)域位移誤差比為3.00%，K 剖分下該誤差比進(jìn)一步降低至1.74%?？梢婋S網(wǎng)格加密，規(guī)則區(qū)域修復(fù)的位移呈加速收斂之勢，即經(jīng)超收斂恢復(fù)后，規(guī)則域位移收斂階可得到提高。因非規(guī)則劃分不滿足變換的強(qiáng)正規(guī)條件，故非規(guī)則區(qū)域網(wǎng)格結(jié)點(diǎn)位移往往不具有超收斂性質(zhì)[18]。但該區(qū)域網(wǎng)格結(jié)點(diǎn)位移精度仍優(yōu)于單元內(nèi)部，本文方法可將單元域內(nèi)位移恢復(fù)至網(wǎng)格結(jié)點(diǎn)水平，故非規(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)部子域仍具有上述的超收斂效果。但在邊界層附近，導(dǎo)數(shù)精度僅能提高一階。

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

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