袁駟 袁全

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

文獻(xiàn)[1]通過比較桿件結(jié)構(gòu)的矩陣位移法[2]和有限元法[3-4]，得出一個(gè)結(jié)論：即有限元的誤差主要來自于其丟失的單元“固端解”項(xiàng)。其后，基于恢復(fù)單元固端解這一思想，超收斂計(jì)算的單元能量投影(element energy projection，EEP)法得以創(chuàng)立并取得長足發(fā)展，不僅對一維有限元法[5-8]，對二維有限元線法[9]、二維乃至三維有限元法[10-12]都建立了EEP超收斂算法，也得到了數(shù)學(xué)理論上的證明[12-13]。更有意義的是，采用EEP超收斂解替代精確解來估計(jì)常規(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)求解，幾乎都有一個(gè)共同點(diǎn)：因?yàn)榻獯鹗孪任粗荒苡煤篁?yàn)誤差方法，按照有限元求解、誤差估計(jì)、更新網(wǎng)格三步循環(huán)迭代求解。這里的關(guān)鍵問題是：缺少先驗(yàn)定量的誤差估計(jì)。這是因?yàn)?，目前幾乎所有的先?yàn)誤差估計(jì)，都包含了事先不可計(jì)算的因素在其中，難以定量，只能是定性的。

本文作者經(jīng)過對文獻(xiàn)[1]的反思和進(jìn)一步研究，發(fā)現(xiàn)其中的思想精華可得到更大的發(fā)揚(yáng)，使得對常見的結(jié)構(gòu)分析問題，有可能實(shí)現(xiàn)先驗(yàn)的定量的誤差估計(jì)，從而可以不經(jīng)有限元計(jì)算，便可以一舉給出滿足精度要求的網(wǎng)格劃分。

本文對這一最新進(jìn)展做一簡要介紹，并給出初步的數(shù)值結(jié)果。作為初始探索，本文限于一維正定自伴問題的Ritz有限元法，也僅限于無內(nèi)部結(jié)點(diǎn)的低次多項(xiàng)式單元。

1 方法思路

有限元的誤差主要來自于其丟失的單元“固端解”項(xiàng)。這一結(jié)論是否準(zhǔn)確？

下面分兩種情況回答。

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

圖1 二階問題物理模型

(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ù)先對有限元解的誤差做出定量估計(jì)成為可能，此即為本文之核心要義所在。為方便，本文將所提出的方法簡稱為“固端法”，以下對其做進(jìn)一步的介紹。

2 模型問題和有限元解

矩陣位移法的單元涉及軸向(拉壓)和橫向(彎曲)變形的兩種類型單元；雖然在單元上相互沒有耦合，但在整體結(jié)構(gòu)上有相互作用。本文仍然采用文獻(xiàn)[1]中的兩個(gè)常微分方程問題作為本文的模型問題：

(1)二階常微分方程(軸向變形問題，圖1)

(2)四階常微分方程(彎曲變形問題，圖2)

其中，L代表微分算子，u和w代表位移和撓度，k為彈性地基的剛度，f和q為均布載荷，均設(shè)為非負(fù)的常數(shù)。

圖2 四階問題物理模型

用有限元求解時(shí)，和矩陣位移法一致，軸向問題(1)采用線性單元，其解記為uh；彎曲問題(2)采用3次Hermite插值單元，其解記為wh，在單元上分別用結(jié)點(diǎn)位移表示為

熟知，若沒有彈性地基(k=0)，則兩類單元均為精確單元(即矩陣位移法的單元)；而有了彈性地基(k＞0)，則兩類單元都是近似單元了。

圖3 和圖4給出了典型的精確單元和近似單元解答的誤差圖，從中可以看出，精確單元結(jié)點(diǎn)處沒有誤差，近似單元結(jié)點(diǎn)處誤差也比單元內(nèi)的微小，誤差主要來自單元內(nèi)部，來自固端情況。

圖3 例1位移誤差圖

圖4 例2位移誤差圖

3 EEP解和自適應(yīng)

以上兩個(gè)問題的有限元解，都已建立了EEP超收斂解的公式，可以計(jì)算單元(ˉx1，ˉx2)上任意一點(diǎn)x∈(ˉx1，ˉx2)的EEP超收斂解，根據(jù)本文問題直接引用如下。

(1)二階問題[5]

(2)四階問題[7]

其中

特別地，當(dāng)k=0時(shí)，由于是精確單元，其結(jié)點(diǎn)位移是精確的，還滿足Luh=0和Lwh=0。此時(shí)，EEP解亦為精確解，可以直接用來估計(jì)精確單元的誤差(或計(jì)算其精確解)，計(jì)算上也更加簡潔。

二階問題

四階問題

其中?J的分量簡化為

式中，eh為有限元解在單元上的誤差。注意，此時(shí)的誤差計(jì)算并不包含有限元解。

目前，基于EEP解的自適應(yīng)有限元分析已得到長足發(fā)展。由于精確解未知，實(shí)際計(jì)算時(shí)，用EEP超收斂解代替精確解進(jìn)行后驗(yàn)誤差估計(jì)，其求解目標(biāo)是：尋求一個(gè)網(wǎng)格使得有限元解按最大模滿足給定的誤差限tol，即要求滿足

式中，e*h=u*-uh或e*h=w*-wh。這樣，自適應(yīng)分析就需要在初始網(wǎng)格上首先求有限元解，然后計(jì)算EEP解并估計(jì)各個(gè)單元的誤差e*h，對那些不滿足式(10)的單元進(jìn)行二分，形成新網(wǎng)格后再次求有限元解，如此循環(huán)迭代，直至所有單元滿足式(10)為止。

4 固端法

自適應(yīng)過程中，最耗時(shí)的是各級網(wǎng)格上的有限元求解及相應(yīng)的EEP解的計(jì)算，能否盡量減少、甚至避免有限元計(jì)算？并非沒有可能。

為此，稍加仔細(xì)地考察一下精確單元的自適應(yīng)。如前文所述，精確單元的誤差完全來自固端解。參見精確單元的EEP式(7)和式(8)，不僅所計(jì)算的誤差eh都是精確的，而且并沒有出現(xiàn)有限元解uh或wh。更利好的是，對于本文的常系數(shù)和常載荷情況，精確單元的誤差eh及其最大值ehmax≡max■■eh■■都可以很簡單地事先算出來：

二階問題精確單元

四階問題精確單元

其實(shí)，以上最大誤差就是結(jié)構(gòu)力學(xué)中熟知的，均布載荷作用下，單元兩端固定時(shí)其中點(diǎn)的最大位移(撓度)。這樣，令ehmax≤tol，便可以一舉確定出所允許的單元長度：二階問題

四階問題

式(13)和式(14)便可以用來直接確定單元的允許長度，而無需反復(fù)地自適應(yīng)迭代計(jì)算。用單元固端解預(yù)估單元的誤差并確定其允許長度，本文將其稱為固端法。

由于近似單元的誤差主要來源于固端解，因此固端法可以有效地應(yīng)用于近似單元。定性角度看，固端法要求單元長度h足夠小，以使有限元解對固端情況達(dá)到滿意的精度，否則對整體結(jié)構(gòu)恐難達(dá)到滿意的精度。對于精確單元，固端法是精確的先驗(yàn)定量誤差估計(jì)；對于近似單元，則是略掉了結(jié)點(diǎn)位移誤差(高階微量)意義下的先驗(yàn)定量估計(jì)。一句話，精確單元的自適應(yīng)有多簡單，近似單元的自適應(yīng)也可以同樣簡單。以下用算例說明。

5 數(shù)值算例

選取與文獻(xiàn)[1]相同的兩個(gè)算例。僅采用均分網(wǎng)格?？紤]兩種問題：(1)給定單元數(shù)，用固端法預(yù)估誤差，并與真實(shí)誤差比較；(2)給定誤差限tol，用固端法預(yù)估誤差、確定網(wǎng)格，并用真實(shí)誤差檢驗(yàn)。算例中，統(tǒng)一取k=f=q=1。雖然沒有單獨(dú)給出精確單元(k=0)的真實(shí)最大誤差，但此時(shí)預(yù)估最大誤差即為其真實(shí)最大誤差，可以參照，也因此沒必要給出。

例1二階問題

問題(1)的求解結(jié)果列于表1左半部，可以看出，用固端法預(yù)估的誤差(h2/8)均稍大于真實(shí)誤差，用來做誤差估計(jì)是安全可靠的。問題(2)的求解結(jié)果列于表1右半部，可以看出，預(yù)估網(wǎng)格的結(jié)果都滿足誤差限，且真實(shí)誤差均略小于預(yù)估誤差。以tol=0.005為例，由式(13)可有h=0.2，恰好5個(gè)單元，其位移誤差圖示于圖3，可直觀看到逐點(diǎn)均滿足誤差限。從表1中還可以看出，取4個(gè)單元?jiǎng)t不能滿足誤差限。

表1 例1二階問題的結(jié)果

例2四階問題

問題(1)的求解結(jié)果列于表2左半部，可以看出，用固端法預(yù)估的誤差(h4/384)均稍大于真實(shí)誤差，用來做誤差估計(jì)是安全可靠的。問題(2)的求解結(jié)果列于表2右半部，可以看出，預(yù)估網(wǎng)格的結(jié)果都滿足誤差限，且真實(shí)誤差均略小于預(yù)估誤差。以tol=0.000 05為例，由式(13)可有h=0.66，理論上需要1.5個(gè)單元，實(shí)際取大于該值的最小整數(shù)，即2個(gè)單元，其位移誤差圖示于圖4，可直觀看到逐點(diǎn)均滿足誤差限。

表2 例2四階問題的結(jié)果

6 結(jié)語

本文提出的固端法，作為先驗(yàn)定量的誤差估計(jì)方法，可以根據(jù)給定的誤差限，直接確定允許的單元長度，一舉得到允許的網(wǎng)格劃分，極大地簡化了有限元誤差估計(jì)的計(jì)算，大量減少自適應(yīng)有限元求解的迭代步驟，并大幅提升自適應(yīng)有限元求解的效率。該法在其他問題中亦獲得有效的應(yīng)用和推廣，甚至在二維有限元中也初步顯示成功，這些將另文介紹。