降階單元的新進展：內置了最大模誤差估計器的自適應有限元法初探

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

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

自適應有限元分析可以提升求解質量、精度和效率，是近年來數值計算方法研究的熱點[1-10]。自適應求解的前提是，能夠對有限元解進行可靠的誤差估計，以指導自適應網格劃分。一維有限元的誤差估計相對簡單，對于二維和三維有限元，縱觀當前大多數后驗誤差估計方法，都存在或部分存在如下缺點和弊端：① 不能在單個單元上估計誤差，需要結構化網格，需要若干單元聯片；② 不能逐點估計誤差，難以按最大?？刂普`差；③ 線性元由于整體上缺失超收斂性，難以構造可靠的誤差估計器；④ 誤差估計缺乏理論證明，可靠性缺乏理論保障。

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

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

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

1 問題描述及有限元解

1.1 模型問題及有限元解

考慮如下二階非自伴兩點邊值問題：

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

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

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

1.2 常規二次元

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

式中：

注意：N?3是一個兩端點為0 的二次“泡狀”函數；相應的，u?h3為廣義結點位移，并不代單元中點位移。則Galerkin 有限元解歸結為求解，使得：

2 降階線性元及其特性

2.1 降階線性元的結構

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

式中：

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

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

2.2 降階線性元的特性

降階線性元解并不等價于常規線性元解。根據一維有限元的數學理論[11-14]，常規線性元解、二次元解及降階線性元解在單元內部和端結點處的收斂階，歸納于表1，由表1 可見并可推斷出如下特性：

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

1) 降階線性元的解是線性元和二次元混合解：內部是線性解，結點處是二次單元的解。

2) 降階線性元的解，在單元內部的精度為h2階，在結點處的精度為h4階，可稱之為超線性解。

3) 結點處h4階的超收斂精度，可以非常有效地優先大幅減小結點位移的誤差，將誤差集中在單元內部，使得誤差模式得到統一。故在做誤差估計時，只估計單元內部誤差即可，基本上可以排除結點誤差的影響。

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

綜上，降階單元的核心思想是：常規二次單元的解，包含了一個降階線性元解和一個誤差估計項。用降階單元的解作為最終解，用內置的誤差估計項估計誤差，指導網格細分，即得到一個十分簡單方便的自適應有限元求解策略。

3 自適應算法

3.1 自適應求解目標及算法

本文的自適應求解的終極目標為：由用戶事先給定誤差限tol，尋找一個優化的有限元網格，使得該網格上的最終有限元解uh按照最大模度量滿足誤差限，即逐單元滿足：

注意，按最大模自適應得到的是一個逐點滿足誤差限的解答，可看作在誤差限tol意義下的數值精確解。

則，實際計算時，停機準則為逐單元滿足：

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

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

1) 給定初始網格和誤差限tol。

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

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

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

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

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

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

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

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

4 數值算例

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

問題描述如下：

其精確解如圖3 所示，表達式較復雜，不再給出。為了展示降階法優于雙階法，本例中也給出雙階法的部分結果，其誤差比記作。給定誤差限tol=1×10-4，以均勻分布的16 個單元作為初始網格，計算過程和結果如圖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)

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

圖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由原方程導出。四邊形區域及初始網格如圖10 所示，給定誤差限tol=1×10-3。

圖10 四邊形區域示意圖及初始網格Fig.10 Schematic diagram of quadrilateral area and initial mesh

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

圖11 例2 最終網格劃分Fig.11 Final mesh of example 2

圖12 例2 初始網格域內誤差比的分布Fig.12 Error ratio on the initial mesh of example 2

圖13 例2 最終網格域內誤差比的分布Fig.13 Error ratio on the final mesh of example 2

此問題是SPR 自適應方法用于檢驗算法的算例[18]，問題描述如下：

給定問題精確解為

荷載f由原方程導出。此問題沿求解域對角線存在應力集中現象。初始網格取1 個單元，給定誤差限tol=1×10-3。

計算中經過6 步自適應調整，得到如圖14 所示的共763 個單元的自適應最終網格，從中可見，網格逐漸向對角線附近加密，顯示出問題的難點所在。圖15 給出最終網格下降階線性元解的誤差比分布，可以看出誤差比都在 ±1以內。本法對存在一定高梯度的問題可以給出滿足誤差限要求的解答。

圖14 例3 最終網格劃分Fig.14 Final mesh of example 3

圖15 例3域內誤差比?的分布Fig.15 Error ratio of example 3

例4.彈性力學平面問題——Cook 梁

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

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

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

圖17 例4 最終網格劃分Fig.17 Final mesh of example 4

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

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

5 結論

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

本文以二階非自伴兩點邊值問題為例，在初值問題降階單元的基礎上，進一步提出邊值問題的降階線性元及其自適應算法。降階單元具有高精度的結點解，同時內置了有效可靠的誤差估計器，從單元構造、數值實施來看，也更加簡單、便捷、靈活、通用，具有明顯的優勢。