基于EEP 技術(shù)的一維有限元結(jié)點(diǎn)位移誤差計(jì)算

袁 駟，邢沁妍，袁 全

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

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

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

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

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

1 模型問題和EEP 解

1.1 模型問題和有限元解

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

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

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

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

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

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

根據(jù)一維有限元誤差估計(jì)理論，當(dāng)問題式(1)的解足夠光滑且網(wǎng)格足夠密(單元長度h足夠小)時(shí)，由m次元(試探函數(shù)和檢驗(yàn)函數(shù)同次)得到的有限元解在單元端結(jié)點(diǎn)上獲得O(h2m)[2?3]的超收斂性，而在單元內(nèi)部則為常規(guī)的O(hm+1)[1]的收斂性，即：

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

1.2 EEP 超收斂解

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

2 結(jié)點(diǎn)位移誤差的計(jì)算

2.1 若干討論

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

2.2 問題提法

2.3 具體算法

3 算法的延伸與拓展

第2.3 節(jié)算法是一個(gè)基本算法，應(yīng)用中可以有多種變化和拓展，簡述幾點(diǎn)如下：

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

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

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

4 數(shù)值算例

本節(jié)給出若干典型算例，用以驗(yàn)證本法優(yōu)良的內(nèi)在特性和多樣化的應(yīng)用場景。本文采用符號計(jì)算軟件MAPLE 計(jì)算，其優(yōu)點(diǎn)在于用戶可以定義數(shù)值的字長。

例1. 收斂階

下面給出數(shù)值算例來驗(yàn)證以上的收斂階。

表1 結(jié)點(diǎn)收斂階(簡約格式)Table 1 Nodal convergence orders (Simplified form)

表2 結(jié)點(diǎn)收斂階(凝聚格式)Table 2 Nodal convergence orders (Condensed form)

例2. 反復(fù)修正結(jié)點(diǎn)

表4 是4 個(gè)線性元的結(jié)點(diǎn)誤差，可以看到網(wǎng)格加密后精度隨之提高。而且，盡管沒有展示，單元次數(shù)更高也會得到更高的精度。

表3 反復(fù)修正結(jié)點(diǎn)位移Table 3 Repeated nodal value updating

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

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

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

表4 4 個(gè)線性元的結(jié)點(diǎn)誤差Table 4 Errors of nodal values of four linear elements

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

表5 結(jié)點(diǎn)誤差和收斂階( 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 結(jié)論

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