再談從矩陣位移法看有限元位移精度的損失與恢復1)

袁駟 袁全

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

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

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

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

本文對這一最新進展做一簡要介紹，并給出初步的數值結果。作為初始探索，本文限于一維正定自伴問題的Ritz有限元法，也僅限于無內部結點的低次多項式單元。

1 方法思路

有限元的誤差主要來自于其丟失的單元“固端解”項。這一結論是否準確？

下面分兩種情況回答。

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

圖1 二階問題物理模型

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

這樣，精確單元和近似單元的誤差的主要來源便得到了統(tǒng)一，即單元的固端解項。然而，求固端解是局部單元的問題，并不需要作整體的有限元求解。這就使得不經有限元求解而預先對有限元解的誤差做出定量估計成為可能，此即為本文之核心要義所在。為方便，本文將所提出的方法簡稱為“固端法”，以下對其做進一步的介紹。

2 模型問題和有限元解

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

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

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

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

圖2 四階問題物理模型

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

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

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

圖3 例1位移誤差圖

圖4 例2位移誤差圖

3 EEP解和自適應

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

(1)二階問題[5]

(2)四階問題[7]

其中

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

二階問題

四階問題

其中?J的分量簡化為

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

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

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

4 固端法

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

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

二階問題精確單元

四階問題精確單元

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

四階問題

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

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

5 數值算例

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

例1二階問題

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

表1 例1二階問題的結果

例2四階問題

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

表2 例2四階問題的結果

6 結語

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