袁 駟，楊 帥，袁 全，王亦平

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

自適應(yīng)有限元分析可以提升求解質(zhì)量、精度和效率，是近年來數(shù)值計(jì)算方法研究的熱點(diǎn)[1-10]。自適應(yīng)求解的前提是，能夠?qū)τ邢拊膺M(jìn)行可靠的誤差估計(jì)，以指導(dǎo)自適應(yīng)網(wǎng)格劃分。一維有限元的誤差估計(jì)相對簡單，對于二維和三維有限元，縱觀當(dāng)前大多數(shù)后驗(yàn)誤差估計(jì)方法，都存在或部分存在如下缺點(diǎn)和弊端：① 不能在單個(gè)單元上估計(jì)誤差，需要結(jié)構(gòu)化網(wǎng)格，需要若干單元聯(lián)片；② 不能逐點(diǎn)估計(jì)誤差，難以按最大?？刂普`差；③ 線性元由于整體上缺失超收斂性，難以構(gòu)造可靠的誤差估計(jì)器；④ 誤差估計(jì)缺乏理論證明，可靠性缺乏理論保障。

對于m(≥1)次常規(guī)多項(xiàng)式單元，有限元數(shù)學(xué)理論有一個(gè)最基本的誤差估計(jì)，即若問題足夠光滑，則一維單元到三維單元在單元內(nèi)部任一點(diǎn)的位移無例外地具有O(hm+1)的收斂性[11-14]。這一估計(jì)意味著二分加密網(wǎng)格(h→h/2)或提高單元次數(shù)(m→m+1)都能獲得精度更高的解，可用來估計(jì)原有限元解的誤差。如此便很自然地派生出兩種略顯原始且笨拙的誤差估計(jì)策略：一是用二分加密網(wǎng)格上的解估計(jì)原網(wǎng)格上解的誤差，本文稱為“雙元法”；二是用高一次單元的解估計(jì)原單元解的誤差，本文稱為“雙階法”。顯然，為了得到兩套不同精度的解答，這兩種方法都需要作兩次有限元求解，計(jì)算代價(jià)過大成為其難以讓人接受的最大弊端；此外，各自還有其他不利因素，將在后文討論。

本文提出降階單元自適應(yīng)分析方法，對其初步的研究表明，該法可以克服上述所有缺點(diǎn)和弊端。降階單元的概念和作法首先由文獻(xiàn)[15]在初值問題中提出，是對初值問題中凝聚單元的逆構(gòu)思和逆運(yùn)用[16-17]，現(xiàn)已成為結(jié)構(gòu)動力計(jì)算中的一種新型高效、可按最大模自適應(yīng)步長的時(shí)程單元。本文在此基礎(chǔ)上，進(jìn)一步將其推廣到邊值問題，其基本作法是：采用m+1次常規(guī)單元做常規(guī)求解，得到常規(guī)有限元解后，略掉其m+1次項(xiàng)，得到m次的降階解，作為降階單元的解答，將略掉的m+1次項(xiàng)轉(zhuǎn)而作為降階解的誤差估計(jì)器?？梢姡惦A單元只經(jīng)一次有限元求解便得到兩套不同精度的解答，而兩套解在精度上的“階差”很自然地提供了一個(gè)可作逐點(diǎn)估計(jì)的誤差估計(jì)器，也自然地發(fā)展出以降階單元作為最終解答的自適應(yīng)算法。理論和算例均表明，該算法思路精巧、理論清晰、算法簡單、實(shí)施靈活，既可靠又通用。

本文暫只限于介紹降階線性元(m=1)及其自適應(yīng)算法，文中以常規(guī)的二階非自伴兩點(diǎn)邊值問題為模型問題進(jìn)行推導(dǎo)和說明，也給出二維Poisson 方程和彈性力學(xué)平面問題的數(shù)值算例，以展示本文方法的廣泛適用性和有效性。

1 問題描述及有限元解

1.1 模型問題及有限元解

考慮如下二階非自伴兩點(diǎn)邊值問題：

式中，p、r、q均為x的函數(shù)，且p≥p0＞0，p0為常數(shù)。定義與問題(1)相應(yīng)的雙線性型和線性型為：

采用Galerkin 法求解式(1)所示問題的近似解，可歸結(jié)為求解u∈HE1，使得：

式中：u、v分別為試探函數(shù)和檢驗(yàn)函數(shù)；HE1為滿足本質(zhì)(位移)邊界條件直到一階導(dǎo)數(shù)均平方可積的函數(shù)空間。

1.2 常規(guī)二次元

不失一般性，考慮典型區(qū)間x∈[0,h]上的單元，h為單元長度。按照常規(guī)做法，將試探函數(shù)uh和檢驗(yàn)函數(shù)vh均取為二次多項(xiàng)式；又為了后期方便，將形函數(shù)取為升階譜形式：

式中：

注意：N?3是一個(gè)兩端點(diǎn)為0 的二次“泡狀”函數(shù)；相應(yīng)的，u?h3為廣義結(jié)點(diǎn)位移，并不代單元中點(diǎn)位移。則Galerkin 有限元解歸結(jié)為求解，使得：

2 降階線性元及其特性

2.1 降階線性元的結(jié)構(gòu)

在得到上述二次單元的解答后，可將各個(gè)單元的解答寫為：

式中：

圖1 一維降階單元示意圖Fig.1 Schematic diagram of one-dimensional reduced element

圖2 二維降階單元示意圖Fig.2 Schematic diagram of two-dimensional reduced element

2.2 降階線性元的特性

降階線性元解并不等價(jià)于常規(guī)線性元解。根據(jù)一維有限元的數(shù)學(xué)理論[11-14]，常規(guī)線性元解、二次元解及降階線性元解在單元內(nèi)部和端結(jié)點(diǎn)處的收斂階，歸納于表1，由表1 可見并可推斷出如下特性：

表1 一維常規(guī)單元和降階單元位移解的收斂階Table 1 Convergence orders of displacement solutions for one-dimensional conventional and reduced elements

1) 降階線性元的解是線性元和二次元混合解：內(nèi)部是線性解，結(jié)點(diǎn)處是二次單元的解。

2) 降階線性元的解，在單元內(nèi)部的精度為h2階，在結(jié)點(diǎn)處的精度為h4階，可稱之為超線性解。

3) 結(jié)點(diǎn)處h4階的超收斂精度，可以非常有效地優(yōu)先大幅減小結(jié)點(diǎn)位移的誤差，將誤差集中在單元內(nèi)部，使得誤差模式得到統(tǒng)一。故在做誤差估計(jì)時(shí)，只估計(jì)單元內(nèi)部誤差即可，基本上可以排除結(jié)點(diǎn)誤差的影響。

4) 二次單元內(nèi)部解uh具有h3階精度，比降階線性解高一階，用來逐點(diǎn)估計(jì)uhR的誤差具有理論上的合理性和有效性，其誤差恰為，堪稱為預(yù)先內(nèi)置了最大模估計(jì)器。

綜上，降階單元的核心思想是：常規(guī)二次單元的解，包含了一個(gè)降階線性元解和一個(gè)誤差估計(jì)項(xiàng)。用降階單元的解作為最終解，用內(nèi)置的誤差估計(jì)項(xiàng)估計(jì)誤差，指導(dǎo)網(wǎng)格細(xì)分，即得到一個(gè)十分簡單方便的自適應(yīng)有限元求解策略。

3 自適應(yīng)算法

3.1 自適應(yīng)求解目標(biāo)及算法

本文的自適應(yīng)求解的終極目標(biāo)為：由用戶事先給定誤差限tol，尋找一個(gè)優(yōu)化的有限元網(wǎng)格，使得該網(wǎng)格上的最終有限元解uh按照最大模度量滿足誤差限，即逐單元滿足：

注意，按最大模自適應(yīng)得到的是一個(gè)逐點(diǎn)滿足誤差限的解答，可看作在誤差限tol意義下的數(shù)值精確解。

則，實(shí)際計(jì)算時(shí)，停機(jī)準(zhǔn)則為逐單元滿足：

可見，降階單元的誤差估計(jì)既不用逐點(diǎn)搜索，也不用逐點(diǎn)計(jì)算，簡單、精確、高效。

至此，降階單元的自適應(yīng)求解算法可歸納為：

1) 給定初始網(wǎng)格和誤差限tol。

2) 求二次有限元解uh，提取降階線性元解。

3) 逐單元檢驗(yàn)式(11)是否成立?

4) 對不滿足式(11)的單元，將其二等分(一維單元)或四等分(二維單元)，返回到2)；

5) 直至所有單元都滿足式(11)，則得到滿足tol的降階線性元解。

本小節(jié)將對本文方法與引言中提及的雙階法和雙元法分別作一簡要分析和比較。

雙階法：在同一網(wǎng)格下，分別用m次和m+1次的常規(guī)單元進(jìn)行求解，取m次單元的解作為最終解，用m+1次單元的解來估計(jì)其誤差，指導(dǎo)自適應(yīng)過程。與降階法相比，此法存在兩點(diǎn)主要弊端：一是需要在同一網(wǎng)格上進(jìn)行兩次求解，計(jì)算代價(jià)過大(降階法只需一次求解)；二是最終有限元解的端結(jié)點(diǎn)精度沒有提高(降階單元至少提高一階)，導(dǎo)致誤差模式紛雜，對某些問題自適應(yīng)過程不穩(wěn)定、最終單元數(shù)會無序增加(見第4 節(jié)的例1)。

雙元法：給定一個(gè)初始網(wǎng)格，對其所有單元進(jìn)行二分加密后得到第二個(gè)網(wǎng)格，對兩重網(wǎng)格都采用m次常規(guī)單元進(jìn)行求解，將稀疏網(wǎng)格的解作為最終解，用密集網(wǎng)格的解來估計(jì)其誤差，指導(dǎo)自適應(yīng)過程。與降階法相比，此法存在3 點(diǎn)主要弊端：一是需要進(jìn)行兩次求解，計(jì)算代價(jià)過大；二是高次元二分后增加的自由度幾乎翻倍(一個(gè)三次元二分后增加3 個(gè)自由度)，進(jìn)一步增加了計(jì)算代價(jià)；三是檢驗(yàn)解和最終解是同階收斂，簡單直接的誤差估計(jì)并不可靠。

相比之下，無論是計(jì)算效率還是計(jì)算精度，降階法比雙階法和雙元法都具有明顯優(yōu)勢，是更佳選擇。降階法巧妙地將兩次求解合二為一，在提升了降階解結(jié)點(diǎn)精度的同時(shí)，為各個(gè)單元內(nèi)置了各自的誤差估計(jì)器，自給自足、靈活方便。

4 數(shù)值算例

本文采用Fortran90 所編寫的程序代碼計(jì)算一維和二維例題，以驗(yàn)證并展示本法的有效性和可靠性。為方便起見，本文引入“誤差比”，即誤差與誤差限之比，記作又記Ne為單元數(shù)，Nadp為自適應(yīng)步數(shù)。例1.二階非自伴兩點(diǎn)邊值問題

問題描述如下：

其精確解如圖3 所示，表達(dá)式較復(fù)雜，不再給出。為了展示降階法優(yōu)于雙階法，本例中也給出雙階法的部分結(jié)果，其誤差比記作。給定誤差限tol=1×10-4，以均勻分布的16 個(gè)單元作為初始網(wǎng)格，計(jì)算過程和結(jié)果如圖4～圖9 所示。

圖3 例1 的位移精確解Fig.3 Exact displacement solution for example 1

圖4 降階法誤差比( Ne=16,Nadp=0)Fig.4 Error ratio of reduced order method(Ne=16,Nadp=0)

由所示結(jié)果可見，雙階法出現(xiàn)單元端結(jié)點(diǎn)誤差大于內(nèi)部誤差的情況(如圖5、圖7 所示)，誤差并非集中在單元內(nèi)部，使得網(wǎng)格的調(diào)整次數(shù)增多，最終經(jīng)過8 步自適應(yīng)調(diào)整，共劃分349 個(gè)單元作為最終網(wǎng)格(圖9)。而降階法在自適應(yīng)過程中始終保持著如圖4、圖6 所示的穩(wěn)定且統(tǒng)一的誤差模式，即單元端結(jié)點(diǎn)誤差遠(yuǎn)小于內(nèi)部誤差，最終僅經(jīng)過5 步自適應(yīng)調(diào)整，共劃分123 個(gè)單元作為最終網(wǎng)格(圖8)?？梢妰煞m然都給出了滿足誤差限要求的解答，但降階法相比于雙階法更加簡單、高效。例2.Poisson 方程——四邊形區(qū)域問題

圖5 雙階法誤差比( Ne=16,Nadp=0)Fig.5 Error ratio of double order method ( Ne=16,Nadp=0)

圖6 降階法誤差比( Ne=32,Nadp=1)Fig.6 Error ratio of reduced order method(Ne=32,Nadp=1)

圖7 雙階法誤差比( Ne=32,Nadp=1)Fig.7 Error ratio of double order method ( Ne=32,Nadp=1)

圖8 降階法誤差比( Ne=123,Nadp=5)Fig.8 Error ratio of reduced order method(Ne=123,Nadp=5)

圖9 雙階法誤差比( Ne=349,Nadp=8)Fig.9 Error ratio of double order method(Ne=349,Nadp=8)

問題描述如下：

給定問題精確解為：

式中：

荷載f由原方程導(dǎo)出。四邊形區(qū)域及初始網(wǎng)格如圖10 所示，給定誤差限tol=1×10-3。

圖10 四邊形區(qū)域示意圖及初始網(wǎng)格Fig.10 Schematic diagram of quadrilateral area and initial mesh

本例有意給定圖10 所示的 4×4非結(jié)構(gòu)化網(wǎng)格作為初始網(wǎng)格，用以檢驗(yàn)本法對非結(jié)構(gòu)化網(wǎng)格的適用性。計(jì)算中經(jīng)過5 步自適應(yīng)調(diào)整，得到如圖11所示的共2383 個(gè)單元的最終網(wǎng)格。圖12 給出初始網(wǎng)格下的誤差比分布，可見在單元角結(jié)點(diǎn)處誤差比最小，即結(jié)點(diǎn)位移具有更高精度，符合降階單元的預(yù)期誤差模式，也是優(yōu)勢之一。圖13 給出最終網(wǎng)格下降階線性元的誤差比分布，可以看出在 ±1以內(nèi)。本例顯示本法適用于非規(guī)則區(qū)域、非結(jié)構(gòu)化網(wǎng)格，這是本法的一個(gè)突出特色，為實(shí)際應(yīng)用提供了極大的方便靈活性。例3.Poisson 方程——SPR 算例

圖11 例2 最終網(wǎng)格劃分Fig.11 Final mesh of example 2

圖12 例2 初始網(wǎng)格域內(nèi)誤差比的分布Fig.12 Error ratio on the initial mesh of example 2

圖13 例2 最終網(wǎng)格域內(nèi)誤差比的分布Fig.13 Error ratio on the final mesh of example 2

此問題是SPR 自適應(yīng)方法用于檢驗(yàn)算法的算例[18]，問題描述如下：

給定問題精確解為

荷載f由原方程導(dǎo)出。此問題沿求解域?qū)蔷€存在應(yīng)力集中現(xiàn)象。初始網(wǎng)格取1 個(gè)單元，給定誤差限tol=1×10-3。

計(jì)算中經(jīng)過6 步自適應(yīng)調(diào)整，得到如圖14 所示的共763 個(gè)單元的自適應(yīng)最終網(wǎng)格，從中可見，網(wǎng)格逐漸向?qū)蔷€附近加密，顯示出問題的難點(diǎn)所在。圖15 給出最終網(wǎng)格下降階線性元解的誤差比分布，可以看出誤差比都在 ±1以內(nèi)。本法對存在一定高梯度的問題可以給出滿足誤差限要求的解答。

圖14 例3 最終網(wǎng)格劃分Fig.14 Final mesh of example 3

圖15 例3域內(nèi)誤差比?的分布Fig.15 Error ratio of example 3

例4.彈性力學(xué)平面問題——Cook 梁

Cook 梁問題是由Cook 首先提出的一個(gè)經(jīng)典考題[19-20]，它是如圖16 所示的一端受剪力作用的變截面梯形懸臂深梁，本例的特點(diǎn)是A點(diǎn)的應(yīng)力奇異。按平面應(yīng)力問題進(jìn)行純數(shù)值求解，給定彈性模量E=1 ，Poisson 比ν=1/3 ，厚度t=1，剪力P=1 ， 誤差限tol=5×10-2。此外，右端剪力P等效地轉(zhuǎn)化為沿截面為二次拋物線分布，以避免在上、下角點(diǎn)出現(xiàn)與理論不自洽。

圖16 Cook 梁示意圖及初始網(wǎng)格Fig.16 Schematic diagram of Cook beam and initial mesh

計(jì)算中以如圖16 所示的2 個(gè)單元作為初始網(wǎng)格，經(jīng)過6 步自適應(yīng)調(diào)整，得到如圖17 所示的共251 個(gè)單元的最終網(wǎng)格，從中可見，網(wǎng)格逐漸向A點(diǎn)附近加密，顯示出問題的難點(diǎn)所在。圖18、圖19 分別給出最終網(wǎng)格下兩個(gè)位移分量的誤差比分布，可以看出都在 ±1以內(nèi)。本法對彈性力學(xué)平面問題可以給出滿足誤差限要求的解答。

圖17 例4 最終網(wǎng)格劃分Fig.17 Final mesh of example 4

圖18 例4域內(nèi)誤差比(u)的分布Fig.18 Error ratio (u)ofexample 4

圖19 例4 域內(nèi)誤差比(v)的分布Fig.19 Error ratio(v) of example 4

5 結(jié)論

理論分析和大量的數(shù)值試驗(yàn)表明，降階法相較于其他自適應(yīng)方法有如下幾點(diǎn)優(yōu)勢：① 只需求解一個(gè)二次單元的解答，其包含所需的最終解(降階單元解)和一個(gè)誤差估計(jì)器；② 可以有效減少自適應(yīng)所需自由度數(shù)；③ 有限元解的誤差模式得到統(tǒng)一，即結(jié)點(diǎn)誤差相較單元內(nèi)部誤差為高階微量，且各單元誤差單調(diào)減小，極少出現(xiàn)單元反叛情況；④ 算法簡單，便于實(shí)施，可靠通用。

本文以二階非自伴兩點(diǎn)邊值問題為例，在初值問題降階單元的基礎(chǔ)上，進(jìn)一步提出邊值問題的降階線性元及其自適應(yīng)算法。降階單元具有高精度的結(jié)點(diǎn)解，同時(shí)內(nèi)置了有效可靠的誤差估計(jì)器，從單元構(gòu)造、數(shù)值實(shí)施來看，也更加簡單、便捷、靈活、通用，具有明顯的優(yōu)勢。