袁 駟，袁 全

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

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

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

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

1 一維問(wèn)題簡(jiǎn)述

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

1) 精確單元：其形函數(shù)是齊次控制微分方程的通解，結(jié)點(diǎn)位移是精確的，固端解不折不扣地就是精確單元內(nèi)部的誤差。常截面桿件矩陣位移法的單元即是精確單元，參見(jiàn)文獻(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)收斂的。所以，對(duì)于近似單元，特別是高次單元，相對(duì)于單元內(nèi)部位移的誤差，結(jié)點(diǎn)位移的誤差是高階微量，可以合理地將其略去，亦即近似單元誤差的主要來(lái)源亦為固端解。

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

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

2 模型問(wèn)題和有限元解

2.1 Poisson 方程

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

2.2 二維有限元解

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

2.3 擬線法網(wǎng)格

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

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

2.4 有限元線法

以下簡(jiǎn)單介紹FEMOL[22]的概念。圖2 為一個(gè)典型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)問(wèn)題的結(jié)果Table 1 Results of elastic torsion problem in Example 1

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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