自適應步長時程分析的新進展：從凝聚單元到降階單元

袁 駟，袁 全

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

結構動力響應問題是工程計算的重要課題，時程積分方法是最有效的方法之一[1?2]，而自適應步長求解也成為近幾年研究的熱點之一。就自適應求解而言，多數研究都是對離散解進行誤差估計并建立自適應算法[3?4]，但在時域上逐點按最大模控制誤差的自適應求解更為困難，也極為少見，而這恰是本文的目標。

采用有限元法進行時程分析，其解是連續的，為最大模控制誤差提供了前提條件；為了逐點控制誤差，需要有更高精度的超收斂的解。依據有限元數學原理[5]，袁駟等[6?8]提出的單元能量投影(Element Energy Projection，簡稱EEP)法是一種非常有效的后處理的超收斂算法，可用于進行逐點誤差估計和精度修正。該EEP 算法已廣泛應用于一維[9?10]、二維[11]，乃至三維有限元中[12]。EEP技術在自適應步長時程積分的有限元分析中也起到了關鍵的作用[13?15]：基于EEP 技術，可對常規單元構建單步法遞推公式[13]，可對有限元解逐點進行誤差估計[13?14]，可創建自適應步長算法公式[15]，進而可再次運用EEP 技術對結點位移修正[16]來實現更高效的自適應求解。

為了能夠成功實現高效可靠的自適應步長求解，經本文作者課題組持續研究，先后構造了常規單元[13?14]、修正單元[13?14]、凝聚單元[17?18]，乃至本文提出的降階單元，不斷提升時程單元的性能。文本將對這一系列單元做一回顧、分析和評價，從而引出最新的降階單元。降階單元是一個采用了非常規算法的常規單元，簡單實用，高效可靠，是至今最高性能的單元之一。

為了有助于整體理解，這里對文中內容做幾點統一說明：

2) 文中所說結點，如無特殊說明，均指單元端結點，亦即有限元網格結點；

4) 文中以常系數單自由度體系為例展開討論，但方法完全適用于變系數多自由度問題；

5) 文中所稱自適應步長時程單元，是指該單元配備了誤差估計器從而可以自適應調整步長的單元；

6) 由于可以逐單元積分求解，不失一般性，僅考慮兩端結點坐標為(t1,t2)=(0,h)的典型單元e，其中h為單元長度，亦為時間步長；

7) 表1 給出了文中所述6 種單元的性態比較，標×的和精度比< 2 的部分是該單元的缺點或不足，正因為這些不足，才導致我們構建了高性能的凝聚單元和降階單元。

表1 各類單元性態比較(次)Table 1 Comparison of various elements of degree

表1 各類單元性態比較(次)Table 1 Comparison of various elements of degree

注： m 次降階單元在單元上的整體精度是m+2階的，但是在估計誤差和控制誤差時被降階為m+1階的。

方程類別單元編號單元類型無條件穩定 EEP計算 結點修正 凝聚計算 結點精度(收斂階) 單元精度(收斂階) 精度比二階運動方程1常規單元[13]×要無?2m m+1<2 2修正單元[13]×要要?2m+2 m+12 3常規單元[14]√要無?2m m+1<2 4修正單元[14]×要要?2m+2一階運動方程m+12 5凝聚單元[17]√要無要2m+2 m+12 6降階單元√無無無2m+2 m+1 m+2,2

8) 總體結論是，作為自適應步長時程單元，表1 中的單元可分為三類：兩種常規單元(單元1和單元3)為基本不可用單元；兩種修正單元(單元2 和單元4)為基本可用單元；而凝聚單元和降階單元(單元5 和單元6)為高性能單元。

1 二階運動方程及其有限元求解

本節介紹求解二階運動方程的兩類有限元：常規單元和修正單元。簡言之，雖然修正單元遠優于常規單元，是基本可用單元，但由于這兩種單元都是有條件穩定的，均不屬于高性能單元。

1.1 二階運動方程

運動方程最常見的表達是如下的二階常微分方程初值問題：

式中：m、c和k分別為質量、阻尼和剛度；P為外荷載；u為動位移；u0、v0分別為初始位移和初始速度；為時域的上界。

為構造Galerkin 弱形式，定義雙線性型和線性型：

則相應的Galerkin 法歸結為求解u∈使得：

1.2 常規單元

本文將以上按常規算法構造的單元稱為二階運動方程的“常規單元”[13]，該單元簡單易行，但有以下幾點缺點(參見表1 的單元1)：

1) 其解uh在結點上是 2階收斂的，與單元精度比小于2，長時間域問題難以控制住誤差；

2) 對于線性元，精度比為1，EEP 解不具有超收斂性，是失效的；

3) 是有條件穩定的，這是一個“硬傷”。

1.3 修正單元

2 一階運動方程及其有限元求解

鑒于二階運動方程難以構造無條件穩定的單元，本節介紹求解一階運動方程的有限元：常規單元和修正單元。簡言之，常規單元是無條件穩定的，但是精度比低于2；而修正單元精度比達到期望值2，可惜又變為有條件穩定的。由于這兩種單元各有不利和不足，也難屬于高性能單元。

2.1 一階運動方程

現將二階運動方程等效地轉換成一階方程組(稱為一階運動方程)[14]：

本文將本節中按常規算法構造的單元稱為一階運動方程的“常規單元”[14]，該單元除了可以同時給出位移和速度，也能同時對位移和速度進行誤差控制這些優點之外，其最大的優點是：各個次數的常規單元均是無條件穩定的[19]。盡管如此，該單元仍有以下缺憾(參見表1 的單元3)：

1) 其解uh在結點上是2 階收斂的，與單元精度比小于2，長時間域問題難以控制住誤差；

2) 對于線性元，精度比為1，EEP 解不具有超收斂性，是失效的。

2.3 修正單元

記EEP 超收斂解的誤差e?=u?u?。將其代回式(11)的原問題，可以得到e?的控制微分方程和初值條件：

本文將按以上做法對結點位移進行修正后的單元稱為一階運動方程的“修正單元”[14]。修正單元大幅提高了結點位移的精度，使其精度比對各次單元均可達到2，特別是也解救了線性元，成為基本可用單元。本以為修正單元在無條件穩定的常規單元的基礎上作精度修正，仍可保留其無條件穩定性，但令人遺憾地發現，修正單元改進了常規單元的前兩個缺點，但將原本的優點反轉為一個缺點(參見表1 的單元4)：不是無條件穩定的。

3 凝聚單元和降階單元

一階運動方程的常規單元是無條件穩定的，所欠缺的是精度比不夠高；而修正單元提升了精度比，但卻喪失了無條件穩定性。本節討論的兩類單元均是在保留常規單元性態的基礎上，進一步提高精度比而構造的。首先介紹凝聚單元，然后推出本文提出的降階單元。簡言之，這兩類單元都是無條件穩定的，無須額外的結點修正技術即可將精度比直接提升為2，而降階單元還自帶超收斂解，無須EEP 超收斂計算。這兩類單元均屬高性能單元。

3.1 凝聚單元

3.2 降階單元

降階單元不涉及新的公式，僅需用語言表述如下：

圖1 例1 的一個單元解(h=3.5)Fig. 1 One element solution of example 1 (h=3.5)

以上一系列算法和操作使得降階單元成為一個十分簡單、簡潔、簡化的“三無”單元：無凝聚操作，無EEP 計算，無結點修正。凝聚單元和降階單元的無條件穩定性、各方面的收斂階等都有嚴格的數學證明，表現也很類似，但在實施上降階單元具有很大優勢，是一個簡單樸實、功能齊全、性能卓越的單元(參見表1 的單元6)。

4 數值算例

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

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

圖3 例1 單元位移誤差比圖Fig. 3 Error ratios of element displacements for example 1

例2. 多自由度簡諧振動

物理模型來源于3 層剪切型框架結構，計算數據為：

簡諧荷載向量為F=10(sin(10t) sin(10t)sin(10t))T，初始位移和初始速度均取0，取10 s，初始步長取為h0=0.5。表3 為三次元的相關計算結果，由其可知，對多自由度本法同樣有效，其性質與單自由度相似，步長分布均勻且合理。圖4、圖5 為三次元、tol=10?3的數據結果。由表3 和圖4 可見，由于其多自由度特性，盡管自適應次數和迭代次數仍然較少，但是多于單自由度問題。由圖5 可知，誤差分布多在誤差上下限區間內。

圖5 例2 單元誤差比圖Fig. 5 Error ratios of elements for example 2

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

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

例3. 與凝聚單元比較算例

為了對降階單元和凝聚單元的特點作一比較，用凝聚單元計算了上述兩個算例，將結果比較列于表4 和表5。可見，兩種單元均能夠圓滿完成自適應求解。總體上，相同次數的降階單元所需要的單元數和變單元長度次數均比凝聚單元少。二者的差別主要是采用了不同的誤差估計器，凝聚單元采用EEP 解估計誤差u??uh，而降階單元采用(t)估計誤差。整體上看，降階單元的表現稍優于凝聚單元。

表4 單自由度簡諧運動凝聚單元對比 (256 s,三次元)Table 4 Comparison with condensed element for damped harmonic motion (256 s,Cubic)

表5 多自由度簡諧運動凝聚單元對比 (10 s,三次元)Table 5 Comparison with condensed element for MDOF system motion (10 s,Cubic)

5 結論

本文以結構動力響應的運動方程為例，在作者新提出的凝聚單元的基礎上，進一步提出降階單元。可以說，這是一對姊妹單元，其特色有：

(1)均是無條件穩定的；

(3)均具備有效可靠的誤差估計器；

(4)同屬于高性能自適應步長時程單元；

(5)從單元構造、實施和計算來看，降階單元更加簡單、簡潔、簡便，也更具優勢。