基于EEP 技術的一維有限元結點位移誤差計算

袁 駟，邢沁妍，袁 全

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

有限元后處理超收斂計算是近年來研究熱點之一[1?5]。根據一維有限元的超收斂理論，無論是邊值問題的Ritz 有限元還是初值問題的Galerkin有限元的位移解，m(≥1)次多項式單元可得到具有h2m階的超收斂端結點位移(真解足夠光滑)，然而單元內部的有限元解僅為hm+1階[1?3]。袁駟等[6?9]提出了有限元后處理中超收斂計算的單元能量投影(Element Energy Projection，簡稱EEP)法，能夠計算一維有限元的逐點超收斂解。該技術雖然是基于一維有限元的，但通過一種逐維(dimensionby-dimension，簡稱D-by-D)恢復策略[10?11]，已經拓展到二維和三維的有限元分析，并應用于多種問題的自適應求解[12?13]。

一維問題的EEP 公式有兩種格式：簡約格式和凝聚格式[7?9]。二者均被證明可以提供單元內任意點位移及其導數的超收斂解。其中，簡約格式的超收斂解精度[14]：位移為O(hm+2) (m>1)，導數為O(hm+1)；凝聚格式的精度[15]：位移及導數均為O(h2m)。EEP 解并不改變單元端結點的位移值，其精度已經很高了，即O(h2m)。然而，更高精度的結點位移對于結點誤差估計及其精度提升均具有重要意義。例如，對于初值問題，結點位移精度對于穩定可靠的時程積分是至關重要的，因此，需要有能夠對結點位移誤差進行有效估計和控制的技術。

本文充分利用EEP 解超收斂特性，提出一種簡便有效的結點位移誤差計算方法。本法在傳統有限元計算后，在不改變現有網格和整體剛度矩陣的情況下，用EEP 解的殘余荷載生成新的荷載向量，只經常規的回代過程，即可得到高精度的結點位移誤差。再將所估計的誤差疊加在傳統有限元解，會得到修正的結點位移，其精度具有更高階的超收斂性：對簡約格式有O(h2m+2)精度，對凝聚格式有O(h3m+mod(m,2))精度。例如，對于線性元，修正后的結點位移具有O(h4)精度。此外，由于修正后的結點位移再次具有更高階的超收斂性，因此可以用其再次改進EEP 解，進而計算另一輪超收斂結點位移。如此反復迭代，最終可獲得精確的結點位移，其中每個迭代步的計算量只是荷載向量的生成和回代。

本文以一般的二階常微分方程為模型問題，給出本法的基本思路和算法，并給出典型算例以展示本法優越的超收斂性以及多樣化的拓展和應用。

1 模型問題和EEP 解

1.1 模型問題和有限元解

本文考慮一般的二階常微分方程邊值問題(BVP)和初值問題(IVP)：

式中： ()′=d()/dx；L為式(1)中定義的線性微分算子。不失一般性，以下將u(x)稱為位移。對于非自伴情況，其伴隨算子為：

定義雙線性型和線性型如下：

則用常規的Galerkin 位移型有限元求解上述問題歸結為求解如下的弱形式問題：

式中：uh∈Sh為有限元解；Sh為m次分段多項式所組成的有限元試探(檢驗)空間。由vh的任意性，可導出通常意義下的剛度方程：

式中：K為整體剛度矩陣；D為結點位移向量；P為等效結點荷載向量。為方便起見，本文中的“結點”在未特別提及時均指“單元兩端結點”，而相應的結點位移表示為uhi。

根據一維有限元誤差估計理論，當問題式(1)的解足夠光滑且網格足夠密(單元長度h足夠小)時，由m次元(試探函數和檢驗函數同次)得到的有限元解在單元端結點上獲得O(h2m)[2?3]的超收斂性，而在單元內部則為常規的O(hm+1)[1]的收斂性，即：

式中， ||·||0和 ||·||m+1分別為整個解域上H0和Hm+1范數。式(7)的超收斂階已經非常高，通常認為是最佳超收斂階。本文提出用同次單元、同一網格進行另一輪有限元分析而計算出結點位移誤差ui?uhi近似解的方法，從而可獲得精度更高的結點位移修正解。這一方法關鍵技術是EEP 超收斂解的應用，將在下一小節中討論。

1.2 EEP 超收斂解

袁駟等[6]基于結構力學中的矩陣位移法和有限元中的投影定理[1]，提出了一種一維有限元后處理超收斂計算的新型方法，簡稱為EEP 法[7?9]。本節簡要介紹該方法的基本思想和相關公式。

2 結點位移誤差的計算

2.1 若干討論

注意到有限元解uh與EEP 解u?在單元端結點處相同，即u?i=uhi，則如果誤差e?=u?u?可以計算，在端點處e?i=ehi。因此，用相同網格、同次單元計算e?的有限元近似解 εh，其結點值 εhi就是有限元結點位移誤差的近似解，而且由于其解為結點值，仍具有超收斂性。此為本文的基本思想，簡單有效，以下給出具體做法。

2.2 問題提法

2.3 具體算法

3 算法的延伸與拓展

第2.3 節算法是一個基本算法，應用中可以有多種變化和拓展，簡述幾點如下：

1)結點位移全量。記D(1)=D+δD，算法步驟4)可改為D(1)=K?1(P+P?)，亦即本法不僅可以計算結點位移增量，也能直接計算修正后的結點位移全量。

2)修正的EEP 解。在求解了 δD得到 εh后，還可以用EEP 法再計算全新的超收斂解 ε?，此即為u?的逐點誤差估計，從而得到更精確的有限元逐點誤差估計eh≈u?+ε??uh。

3)結點的再修正。繼得到上述 ε?后，還可以用本文第2.3 節的算法再計算 ε?的誤差δe?=e??ε?的近似解 δεh(同時也是修正的位移全量的誤差近似解)。 δe?的定解方程和條件為：

4 數值算例

本節給出若干典型算例，用以驗證本法優良的內在特性和多樣化的應用場景。本文采用符號計算軟件MAPLE 計算，其優點在于用戶可以定義數值的字長。

例1. 收斂階

下面給出數值算例來驗證以上的收斂階。

表1 結點收斂階(簡約格式)Table 1 Nodal convergence orders (Simplified form)

表2 結點收斂階(凝聚格式)Table 2 Nodal convergence orders (Condensed form)

例2. 反復修正結點

表4 是4 個線性元的結點誤差，可以看到網格加密后精度隨之提高。而且，盡管沒有展示，單元次數更高也會得到更高的精度。

表3 反復修正結點位移Table 3 Repeated nodal value updating

例3. 初值問題的附加精度

本例中，考慮簡諧荷載作用下的有阻尼的單自由度體系的運動方程：

對于初值問題，沒必要形成整體剛度矩陣，有限元解可以從左到右逐個單元積分求解，因此結點位移修正算法可以逐單元進行，而不用像邊值問題待整體求解之后再作修正。如此，對于初值問題，便出現一個增強的修正方法，其特點是在計算單元右端點有限元位移解uh2時，其左端點的初始位移，可以立即采用對有限元解uh1修正后的結點值。這個逐單元修正的方法，不僅更高效，而且結點精度得到進一步的顯著提高。

表4 4 個線性元的結點誤差Table 4 Errors of nodal values of four linear elements

圖1 例3 的精確解Fig. 1 Exact solution of Example 3

表5 結點誤差和收斂階( m=1，簡約格式)Table 5 Nodal errors and convergence orders ( m=1, simplified form)

圖2 有限元解誤差Fig. 2 Errors of FE solution

圖3 整體修正誤差Fig. 3 Errors of globally updated solution

圖4 逐單元修正誤差Fig. 4 Errors of element-wise updated solution

5 結論

雖然本文的方法和算例主要針對二階常微分方程，但是所提出的方法也適用于其它一維問題，如其它階常微分方程問題以及常微分方程組問題。此外，所提出的方法有潛力通過所謂的逐維修正方法擴展到二維乃至三維有限元分析[10?11]。