運動方程自適應步長求解的高性能Galerkin 時程單元初探

袁 駟，袁 全

(清華大學土木工程系，土木工程安全與耐久教育部重點實驗室，北京 100084)

結構動力響應問題是工程計算的重要課題，常用的時程積分方法更是最有效的方法之一，但其計算效率、計算精度、數值穩定性等均與步長的選擇和優化相關[1?2]，因此自適應步長求解也成為近幾年研究的難點和熱點之一。就自適應求解而言，有對離散解進行誤差估計并建立自適應算法[3? 4]的研究，但在時域上逐點按最大模控制誤差的高效自適應求解更為困難，也極為少見，而這恰是本文追求的目標。

袁駟等依據有限元數學原理[5]提出的單元能量投影[6?8](element energy projection，簡稱EEP)法是一種非常有效的有限元后處理的超收斂算法，可用于進行逐點誤差估計和精度修正。該EEP 算法已應用于一維[9?10]、二維有限元[11]，乃至三維有限元[12]，并且提出一類基于EEP 技術的Galerkin有限元時程積分的自適應方法[13]。該法對常規單元構建了單步法遞推公式，用EEP 超收斂解進行誤差估計，創建了自適應步長算法公式并給出了數學證明[14]，進而再次運用EEP 技術對結點位移修正[15]來實現高效的自適應求解。

文獻[16]將二階運動方程等效地轉換為一階方程組初值問題(以下稱為一階運動方程)，構建了一類無條件穩定的Galerkin 常規單元。但該類單元采用結點位移修正技術后，蛻變為有條件穩定[16]，故需要進一步探索一類無需結點位移修正的、超高結點精度的新型單元，以期實現更加高效、更加穩定、更加可靠的按最大模控制誤差的自適應步長求解算法。

本文基于一階運動方程的非自伴隨性質，嘗試提出一種通過對檢驗函數凝聚而提高單元端結點精度的新型單元——凝聚單元。本文算法相比于文獻[16]的常規單元，在保持了無條件穩定的前提下，對任何次數單元，其端結點位移的收斂階均可提高2 階，極大提高了結點位移精度和自適應步長求解的效率。

1 問題描述與常規Galerkin 有限元

1.1 問題描述

為了表述簡明，本節以單自由度常系數運動方程為例，其可以歸結為如下二階常微分方程初值問題：

其中：

可見，將二階運動方程轉化為一階方程組后，Galerkin 法中對檢驗函數降低了要求，屬于H0空間即可，本文將充分利用這一性質。

1.2 常規Galerkin 有限元

由于可以逐單元步長積分求解，不失一般性，僅考慮兩端結點坐標為(0,h)的典型單元e。注意，對于一階運動方程的Galerkin 弱形式式(5)，若將試探函數uh取為mˉ次多項式，則檢驗函數vh通常取為低一次的mˉ?1次多項式，即：

1.3 EEP 超收斂公式

2 凝聚單元——高性能單元

本節是本文的重點，將提出一套新型高效單元——凝聚單元。研究發現，對于上述的一階運動方程的Galerkin 常規單元作少許巧妙的改進，便可在不提高單元次數且保持無條件穩定的前提下，將各次單元端結點位移的精度普遍提高2 階。本法可概括為“三部曲”，即“同次、凝聚、升階”。簡述如下：

1) 同次，是指將檢驗函數vh的次數由mˉ?1次取為和試探函數uh相同的mˉ次，亦即為vh增加了一個自由度，但并沒有增加單元的次數；

2) 凝聚，是指利用伴隨問題凝聚掉vh的一個自由度，使其還原為原本的自由度數，但更加逼近精確的檢驗函數，亦即構造出凝聚檢驗函數v?，其余則按常規有限元步驟求解；

3) 升階，是指由此得到的單元端結點解答的精度可比同次數常規單元解答的精度普遍提升2 階。

為了幫助理解該單元的構成，本節先對精確檢驗函數做一說明。

2.1 精確檢驗函數

2.2 凝聚單元的構成

2.3 凝聚單元性質

凝聚單元具有諸多優良性質，列舉幾點如下：

1) 凝聚單元，求解運動方程是無條件穩定的，求解一般初值問題是A 穩定的；

2) 凝聚單元，單元內部位移精度仍然為O(hmˉ+1)階；

4) 對于線性元，凝聚單元結點位移精度翻倍，為O(h4)階；

3 自適應步長

本文采用與文獻[13, 16]相同的自適應策略和目標，即在精確解u未知的情況下，事先給定誤差限tol，尋求一個優化的單元步長分布，使得這些單元上有限元解答uh按最大模度量，逐單元滿足：

在確定和調整步長上，同樣采用與文獻[13, 16]相同的步長公式，即：

4 數值算例

本文采用符號計算軟件Maple 計算有阻尼簡諧運動、無阻尼自由振動以及多自由度的算例，以驗證本法的有效性。為方便起見，本文引入誤差比，即誤差與誤差限之比，記作eˉh。

例1. 有阻尼簡諧運動

計算數據為：m=1，k=1，c=0.04，f=sin(0.2t)，u0=0，v0=1 ，Tˉ =256 s，初始步長h0=0.5。表1 為三次元的相關計算結果。圖1、圖2 為三次元、tol=0.001的數據結果。由表1 可知，本法的誤差得到有效控制，最大誤差比均小于1。改變單元長度次數很少，最多不到5 次，最大最小單元長度比合理，體現了本法的高效性。由圖1、圖2 可見，本法步長分布合理，誤差分布均勻，且多數單元誤差在上下限區間內。

表1 有阻尼簡諧振動數據 (256 s，三次元)Table 1 Results of damped harmonic motion (256 s, Cubic)

圖1 例1 步長分布圖Fig. 1 Step-size distribution for example 1

圖2 例1 單元位移誤差比圖Fig. 2 Displacement errors of elements for example 1

例2. 無阻尼自由振動

計算數據為：m=1，k=1，c=0,f=0，u0=0，v0=1 ，Tˉ =256 s，初始步長h=0.5。表2為三次元的相關計算結果。圖3、圖4 為三次元、tol=0.001的數據結果。由表2 和圖3 可知，本法對無阻尼自由振動只需調整步長一次，即可在全部時域上滿足誤差要求，這符合無阻尼自由振動自身的特性。由圖4 可見，其誤差分布幾乎是呈周期分布。

表2 無阻尼自由振動數據 (256 s，三次元)Table 2 Results of undamped free vibration (256 s, Cubic)

圖3 例2 步長分布圖Fig. 3 Step-size distribution for example 2

圖4 例2 單元誤差比圖Fig. 4 Displacement errors of elements for example 2

例3. 多自由度簡諧振動

圖5 例3 步長分布圖Fig. 5 Step-size distribution for example 3

圖6 例3 單元誤差比圖Fig. 6 Displacement errors of elements for example 3

物理模型來源于3 層剪切型框架結構，計算數據為：性質與單自由度相似，步長分布均勻且合理。由表3 和圖4 可見，由于其多自由度特性，盡管自適應次數和迭代次數仍然較少，但是多于單自由度問題。由圖5 可知，誤差分布多在誤差上下限區間內。

表3 多自由度簡諧運動計算數據(256 s, 三次元)Table 3 Results of multi-degree-of-freedom system motion (256 s, Cubic)

5 結論

本文以運動方程為例，提出了凝聚單元。該單元巧妙地利用一階運動方程的非自伴、低檢驗函數的性質，在保持無條件穩定的前提下，以低次單元獲得超高階收斂的結點位移精度，不失為一類高性能自適應步長時程單元。