袁 駟，袁 全

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

文獻(xiàn)[1]通過對矩陣位移法[2]和有限元法[3?4]的分析比較，得出一個結(jié)論：一維有限元的誤差主要來源于各個單元的“固端解”項。其后，基于恢復(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)有限元求解得以實現(xiàn)，其最突出的特點是可以得到按最大模逐點滿足用戶給定的誤差限的解答，可謂是數(shù)值精確解[14?15]。目前，這種自適應(yīng)有限元方法不僅有效地應(yīng)用于各種線性問題，也在特征值問題和多種非線性問題中得到了廣泛而有效的應(yīng)用[16?18]，而近期發(fā)展的網(wǎng)格局部加密技術(shù)為一類刁難奇異問題的自適應(yīng)求解提供了更高性能的求解方案[19?20]。

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

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

1 一維問題簡述

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

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

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

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

對于二維問題，半離散的有限元線法可以看作是廣義一維有限元法(結(jié)點延伸為結(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 先驗誤差控制

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ū)域為非規(guī)則區(qū)域，荷載也是一個函數(shù)，是一個很典型的例題。本例將采用兩種方法估算誤差：一是基于式(8)或式(10)的“直估法”；二是基于式(6)的“先驗法”。

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

為了實施“直估法”，即用式(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)格?！跋闰灧ā眲t用式(6)逐單元估算固端解誤差，需投入少量的計算，但比“直估法”更加精準(zhǔn)一些。

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

表 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é)果都滿足誤差限。還可看出，真實誤差均略小于預(yù)估誤差，安全而可靠。

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

5 結(jié)論

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