袁 駟，袁 全

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

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

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

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

為了有助于整體理解，這里對文中內(nèi)容做幾點統(tǒng)一說明：

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

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

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

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

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

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

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

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

方程類別單元編號單元類型無條件穩(wěn)定 EEP計算 結(jié)點修正 凝聚計算 結(jié)點精度(收斂階) 單元精度(收斂階) 精度比二階運動方程1常規(guī)單元[13]×要無?2m m+1<2 2修正單元[13]×要要?2m+2 m+12 3常規(guī)單元[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) 總體結(jié)論是，作為自適應步長時程單元，表1 中的單元可分為三類：兩種常規(guī)單元(單元1和單元3)為基本不可用單元；兩種修正單元(單元2 和單元4)為基本可用單元；而凝聚單元和降階單元(單元5 和單元6)為高性能單元。

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

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

1.1 二階運動方程

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

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

為構(gòu)造Galerkin 弱形式，定義雙線性型和線性型：

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

1.2 常規(guī)單元

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

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

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

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

1.3 修正單元

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

鑒于二階運動方程難以構(gòu)造無條件穩(wěn)定的單元，本節(jié)介紹求解一階運動方程的有限元：常規(guī)單元和修正單元。簡言之，常規(guī)單元是無條件穩(wěn)定的，但是精度比低于2；而修正單元精度比達到期望值2，可惜又變?yōu)橛袟l件穩(wěn)定的。由于這兩種單元各有不利和不足，也難屬于高性能單元。

2.1 一階運動方程

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

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

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

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

2.3 修正單元

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

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

3 凝聚單元和降階單元

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

3.1 凝聚單元

3.2 降階單元

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

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

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

4 數(shù)值算例

表2 有阻尼簡諧振動結(jié)果 (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 層剪切型框架結(jié)構(gòu)，計算數(shù)據(jù)為：

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

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

表3 多自由度簡諧運動計算結(jié)果 (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. 與凝聚單元比較算例

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

表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 結(jié)論

本文以結(jié)構(gòu)動力響應的運動方程為例，在作者新提出的凝聚單元的基礎(chǔ)上，進一步提出降階單元?？梢哉f，這是一對姊妹單元，其特色有：

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

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

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

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