固端法：二維有限元先驗定量誤差估計與控制

袁 駟，袁 全

(清華大學土木工程系，北京 100084)

文獻[1]通過對矩陣位移法[2]和有限元法[3?4]的分析比較，得出一個結論：一維有限元的誤差主要來源于各個單元的“固端解”項。其后，基于恢復單元固端解這一思想，超收斂計算的單元能量投影(Element Energy Projection，簡稱EEP)法得以創立并取得長足發展，不僅對一維有限元法 (Finite Element Method, 簡稱 FEM)[5?8]，對二維有限元線法(Finite Element Method of Lines, 簡稱FEMOL、線法)[9]、二維乃至三維有限元法[10?11]都建立了EEP 超收斂算法，也得到了數學理論上的證明[12?13]。更有意義的是，用EEP 超收斂解替代精確解來估計常規有限元解的誤差，使得基于EEP 技術的自適應有限元求解得以實現，其最突出的特點是可以得到按最大模逐點滿足用戶給定的誤差限的解答，可謂是數值精確解[14?15]。目前，這種自適應有限元方法不僅有效地應用于各種線性問題，也在特征值問題和多種非線性問題中得到了廣泛而有效的應用[16?18]，而近期發展的網格局部加密技術為一類刁難奇異問題的自適應求解提供了更高性能的求解方案[19?20]。

縱觀各類自適應求解，幾乎都有一個共同(弱)點：因為解答事先未知，只能用后驗誤差方法，按照有限元求解、誤差估計、更新網格三步循環迭代求解。這里的關鍵問題是：缺少先驗定量的誤差估計。這是因為，目前幾乎所有的先驗誤差估計，都包含了事先不可計算的因素在其中，難以定量，只能是定性的。

本文作者經過對文獻[1]的反思和進一步研究，在文獻[21]中對于一維Ritz 有限元法提出了先驗定量誤差估計的“固端法”，從而可以不經有限元計算，便可以一舉給出滿足精度要求的網格劃分。繼“固端法”在一維有限元中的初步成功，本文嘗試將其拓展到二維有限元。作為初始探索，本文以Poisson 方程為例，采用最常規的4 結點線性元，用EEP 技術預先做出誤差估計，能直接給出滿足精度要求的網格劃分。本文對這一最新進展做一簡要介紹，并給出初步的數值結果。

1 一維問題簡述

一維有限元的誤差主要來自于其丟失的單元“固端解”項[1]，可以分兩種情況論述。

1) 精確單元：其形函數是齊次控制微分方程的通解，結點位移是精確的，固端解不折不扣地就是精確單元內部的誤差。常截面桿件矩陣位移法的單元即是精確單元，參見文獻[1]中的圖1：其狀態(Ⅱ)為有限元解，而狀態(Ⅰ)為固端解，亦是有限元解的誤差。

2) 近似單元：其形函數不是齊次控制微分方程的通解，結點位移是近似的，單元內部誤差由固端解和非固端的有限元解共同組成。但是，有限元的數學理論已有證明，有限元的端結點位移相比于單元內部位移是超收斂的，而且具有最佳超收斂性[4]；以C0類單元為例(文獻[1]的表2)， m次單元在單元內部為 O(hm+1)收斂，而在結點上是O(h2m)收斂的。所以，對于近似單元，特別是高次單元，相對于單元內部位移的誤差，結點位移的誤差是高階微量，可以合理地將其略去，亦即近似單元誤差的主要來源亦為固端解。

這樣，精確單元和近似單元的誤差的主要來源便得到了統一，即單元的固端解項。然而，更大的利好是，求固端解是局部單元的問題，并不需要作整體的有限元求解。這就使得不經有限元整體求解而預先對有限元解的誤差做出定量估計成為可能。

對于二維問題，半離散的有限元線法可以看作是廣義一維有限元法(結點延伸為結線)，其誤差的主要來源為線法單元兩端結線(單元邊結線)固定的解，所以仍可以說是來源于“固端解”。而二維有限元法可以看作用一維有限元求解線法的結線位移，而當結線位移事先固定為零時，線法和有限元法就沒有區別了，在有限元網格上也可以定義固端解了。

2 模型問題和有限元解

2.1 Poisson 方程

本文的分析以二維Poisson 方程為模型問題，其形式如下：

2.2 二維有限元解

用二維有限元求解此問題時，本文采用雙向m 次四邊形單元，單元試探函數為：

2.3 擬線法網格

求解域 Ω不必是規則區域，但是為適合二維EEP 超收斂計算，網格劃分須是“擬線法網格”，即先利用FEMOL 的離散方式用一組結線對求解區域進行半離散，然后再沿結線維度進一步離散，得到的網格即是求解FEMOL 常微分方程組(Ordinary Differential Equations，簡稱 ODEs)的一維有限元網格，也是原問題的二維有限元網格，如圖1所示。

圖 1 二維問題的逐維離散Fig. 1 Dimension-by-dimension discretization of 2D problems

2.4 有限元線法

以下簡單介紹FEMOL[22]的概念。圖2 為一個典型FEMOL 二次單元的幾何映射。

圖 2 FEMOL 單元映射Fig. 2 Mapping of FEMOL element

2.5 線法的 EEP 解

3 固端法

3.1 固端解

3.2 線性元公式

3.3 先驗誤差控制

4 數值算例

表 1 例1 彈性扭轉問題的結果Table 1 Results of elastic torsion problem in Example 1

圖 3 方形區域彈性扭轉問題的誤差Fig. 3 Error of elastic torsion problem on a square region

例2. 非規則四邊形求解域問題

本例考慮一個Poisson 方程定義在如圖4 所示的非規則四邊形求解域上，四周為齊次本質邊界條件，問題描述如下：

f 由u 反求得到，形式復雜，在此略去。

本例的區域為非規則區域，荷載也是一個函數，是一個很典型的例題。本例將采用兩種方法估算誤差：一是基于式(8)或式(10)的“直估法”；二是基于式(6)的“先驗法”。

圖 4 非規則四邊形求解域Fig. 4 An irregular quadrilateral solution domain

為了實施“直估法”，即用式(8)或式(10)直接估算單元大小，需對本例荷載和區域作一些近似處理：本例的荷載 f是函數，大約在區域中心，可找到其最大值 f0?0.9；本例的區域非矩形，將其近似為邊長為5 的方形區域，以便估算沿邊長所用的單元數。以上處理，相當于在 5×5的區域上有均布荷載 f0作用的Poisson 方程。這樣就可以用式(8)或式(10)直接估算單元大小(模仿例1)，并直接確定網格?！跋闰灧ā眲t用式(6)逐單元估算固端解誤差，需投入少量的計算，但比“直估法”更加精準一些。

問題(Ⅰ)求解結果列于表2。其中直估法的計算極為簡單，例如對于2×2 網格(參見圖5)，由式(8)可直估誤差為 f0h2/8=0.9·(5/2)2/8=0.7031?？梢钥闯?，簡簡單單的直估法雖然在網格比較稀疏時比先驗法估計得要偏高一些，但隨著網格加密，二者完全趨于一致。這是因為，本例的最大誤差發生在區域中心，而當網格很密、單元很小時，區域中心單元上的荷載趨于常數 f0，因此與直估法趨于一致。這也輔證了直估法的合理性。還可以看到，無論是直估法還是先驗法，所估算的誤差都略大于真實誤差；作為誤差估計器，是偏于安全和可靠的。

表 2 例2 非規則區域問題(I)的結果Table 2 Case (I) results of irregular domain problem in Example 2

圖 5 二維有限元網格Fig. 5 A 2D FEM mesh

問題(Ⅱ)求解結果列于表3。最后一列給出最大誤差比，可以看出，其值均小于1，亦即預估網格的結果都滿足誤差限。還可看出，真實誤差均略小于預估誤差，安全而可靠。

表 3 例2 非規則區域問題(II)的結果Table 3 Case (II) results of irregular domain problem in Example 2

5 結論

本文提出二維有限元先驗定量誤差估計的“固端法”，可以根據給定的誤差限，直接確定允許的單元大小，一舉得到允許的網格劃分，極大地簡化了二維有限元誤差估計的計算，大量減少自適應有限元求解的迭代步驟，并大幅提升自適應有限元求解的效率。該法的更為深入、廣泛、系統的研究成果將另文介紹。