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

袁 駟，袁 全

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

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

縱觀各類自適應(yīng)求解，幾乎都有一個共同(弱)點(diǎn)：因?yàn)榻獯鹗孪任粗荒苡煤篁?yàn)誤差方法，按照有限元求解、誤差估計、更新網(wǎng)格三步循環(huán)迭代求解。這里的關(guān)鍵問題是：缺少先驗(yàn)定量的誤差估計。這是因?yàn)?，目前幾乎所有的先?yàn)誤差估計，都包含了事先不可計算的因素在其中，難以定量，只能是定性的。

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

1 一維問題簡述

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

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

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

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

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

2 模型問題和有限元解

2.1 Poisson 方程

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

2.2 二維有限元解

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

2.3 擬線法網(wǎng)格

求解域 Ω不必是規(guī)則區(qū)域，但是為適合二維EEP 超收斂計算，網(wǎng)格劃分須是“擬線法網(wǎng)格”，即先利用FEMOL 的離散方式用一組結(jié)線對求解區(qū)域進(jìn)行半離散，然后再沿結(jié)線維度進(jìn)一步離散，得到的網(wǎng)格即是求解FEMOL 常微分方程組(Ordinary Differential Equations，簡稱 ODEs)的一維有限元網(wǎng)格，也是原問題的二維有限元網(wǎng)格，如圖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 先驗(yàn)誤差控制

4 數(shù)值算例

表 1 例1 彈性扭轉(zhuǎn)問題的結(jié)果Table 1 Results of elastic torsion problem in Example 1

圖 3 方形區(qū)域彈性扭轉(zhuǎn)問題的誤差Fig. 3 Error of elastic torsion problem on a square region

例2. 非規(guī)則四邊形求解域問題

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

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

本例的區(qū)域?yàn)榉且?guī)則區(qū)域，荷載也是一個函數(shù)，是一個很典型的例題。本例將采用兩種方法估算誤差：一是基于式(8)或式(10)的“直估法”；二是基于式(6)的“先驗(yàn)法”。

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

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

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

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

圖 5 二維有限元網(wǎng)格Fig. 5 A 2D FEM mesh

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

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

5 結(jié)論

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