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

袁 駟，袁 全

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

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

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

文獻[16]將二階運動方程等效地轉(zhuǎn)換為一階方程組初值問題(以下稱為一階運動方程)，構(gòu)建了一類無條件穩(wěn)定的Galerkin 常規(guī)單元。但該類單元采用結(jié)點位移修正技術(shù)后，蛻變?yōu)橛袟l件穩(wěn)定[16]，故需要進一步探索一類無需結(jié)點位移修正的、超高結(jié)點精度的新型單元，以期實現(xiàn)更加高效、更加穩(wěn)定、更加可靠的按最大?？刂普`差的自適應步長求解算法。

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

1 問題描述與常規(guī)Galerkin 有限元

1.1 問題描述

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

其中：

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

1.2 常規(guī)Galerkin 有限元

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

1.3 EEP 超收斂公式

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

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

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

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

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

為了幫助理解該單元的構(gòu)成，本節(jié)先對精確檢驗函數(shù)做一說明。

2.1 精確檢驗函數(shù)

2.2 凝聚單元的構(gòu)成

2.3 凝聚單元性質(zhì)

凝聚單元具有諸多優(yōu)良性質(zhì)，列舉幾點如下：

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

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

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

3 自適應步長

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

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

4 數(shù)值算例

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

例1. 有阻尼簡諧運動

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

表1 有阻尼簡諧振動數(shù)據(jù) (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. 無阻尼自由振動

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

表2 無阻尼自由振動數(shù)據(jù) (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 層剪切型框架結(jié)構(gòu)，計算數(shù)據(jù)為：性質(zhì)與單自由度相似，步長分布均勻且合理。由表3 和圖4 可見，由于其多自由度特性，盡管自適應次數(shù)和迭代次數(shù)仍然較少，但是多于單自由度問題。由圖5 可知，誤差分布多在誤差上下限區(qū)間內(nèi)。

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

5 結(jié)論

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