結構動力方程求解中隱式格式向顯式格式的轉換

李常青，樓夢麟，蔣麗忠

(1.中南大學 土木工程學院建筑工程系，長沙 410000;2.同濟大學 土木工程防災國家重點實驗室，上海 200092)

結構動力計算采用的逐步積分方法，有顯式算法和隱式算法之分。隱式算法需要求解耦聯方程組，當結構自由度數目很大(如上百萬個)時，求解工作量非常大。而顯式算法不需要求解方程組，且計算精度控制的要求比穩定性條件要求更高，因此顯式算法的構造一直是研究的熱點［1－6］。在關注隱式算法與顯式算法區別的同時，隱式算法與顯式算法之間是否有聯系，隱式算法是否可轉換為顯式算法，轉換以后算法的穩定性及精度會發生什么變化卻很少有討論。本文對此進行了深入的研究。

式中:［M］為質量矩陣，［C］為阻尼矩陣，［K］為剛度矩陣，{vi+1}為待求的位移向量，i+1}為廣義的荷載向量，a1，a2為常數，對不同的加速度假設有不同的值。

1 基本思路

廣義線性加速度法是最常用的隱式算法之一，其中紐馬克法與威爾遜－θ法均可看做廣義線性加速度方法，每前進一步，都需要求解耦聯線性方程組:

對式(2)矩陣求逆，得:

將式 (1)兩 邊 同乘 矩 陣

將式(3)代入式(4)，有:

級數展開式(5)中的逆矩陣，得:

式(6)級數展開的收斂條件為攝動項范數不能大于1:

將式(6)代入式(5)，得無窮級數展開的顯式算法:

式(8)中的級數展開n次，即得級數截斷后的近似公式:

式(9)無需求解耦聯線性方程組，避開了系數矩陣奇異時方程組無法求解的情況，因而提高了計算速度，節省了計算內存。

式(9)為對式(1)經求逆、級數展開、級數項截斷后得到的近似，本文稱式(9)為隱式算法(1)的伴隨顯式算法。下面分別詳細分析無條件穩定的常加速度方法的顯式化過程和條件穩定的線性加速度方法的顯式化過程。

2 常加速度方法的顯式化

2.1 公式推導

常加速度方法假定在兩相鄰時間之間，加速度為常數，其算法由下列三式構成:

對照式(1)與式(10)可知，對常加速度方法，有a1=4，a2=2，因此通過前述的顯式化過程，可得常加速度法的對應式(10)的顯式化公式為:

用式(13)求位移，式(10)求速度，式(11)求加速度，即得常加速度方法的伴隨顯式化算法。對應式(7)收斂條件為:

其中:

為任何一種矩陣范數。

考慮單自由度系統，式(14)化為

求解式(15)，得:

令 ζ=0，得:

化為h與T(結構周期)的形式為:

式(14)為級數展開的收斂條件，式(16)則是級數收斂條件在單自由度中的表現。下面的分析將表明，該逆矩陣級數展開的收斂條件對h/T的限制與按奇數項展開所得到算法的穩定性條件對 h/T的限制很接近。

導流裝置的內外壁電解著色后均獲得良好的黑色氧化膜，且肉眼看不出差異。于是采用CM-700d分光測色計分別測量了內外壁與儀器自帶白色校準板之間的總色差(ΔE)。圖3顯示了內壁和外壁隨機抽取的10個測量點的色差：內壁為91.25 ～ 92.06，平均值91.71；外壁為90.36 ～ 92.69，平均值91.41。內外壁的ΔE平均值僅相差0.3，說明內外氧化膜的顏色一致性良好。另外，內壁色差的離散度較外壁更低，說明內壁的顏色更加均勻。

2.2 精度分析

分析式(13)求解位移的精度。分別考慮n=1，2，3，4，5時的算法精度。位移計算誤差定義為σv=ν理論－v計算，通過泰勒展開及符號運算，得到如下結果:

當k=1時，位移算法精度為二階，位移誤差項為:

當k=2時，位移算法精度為三階，位移誤差項為:

當k=3時，位移算法精度為三階，位移誤差項為:

當k=4時，位移算法精度為三階，位移誤差項為:

當k=5時，位移算法精度為三階，位移誤差項為:

隨著展開項數的增加，計算精度越來越高，但不會超過相應隱式算法的最高精度。當展開2項后，已達到了三階精度。展開4項以后，計算精度增長不明顯。當n=3、4、5時，位移算法精度為三階，速度算法精度(式(11))為二階，加速度算法精度亦為二階，因此，本算法可達到二階精度。建議采用n=3。

2.3 穩定性分析

算法的穩定分析，通過計算單自由度系統得到傳遞矩陣，形式如下:

分別考慮n=1，2，3，4，5時算法穩定性的表現。由于傳遞矩陣［A］形式很復雜，尤其當5階展開時，形式更復雜。此處只列n=1時算法的傳遞矩陣［A］，其各分量為:

算法的穩定性要求:

其中:ρ(A)為矩陣A的譜半徑。

圖1為當n取不同值時，傳遞矩陣A的譜半徑隨h/T及阻尼比ξ的變化過程。

圖1 常加速度法伴隨顯式算法傳遞矩陣A譜半徑的變化Fig.1 Spectral radius of transferring matrix of constant acceleration method’s adjoint explicit method

表1 常加速度對應顯式算法穩定域大小(h時間步長、T周期)Tab.1 Stability field of constant acceleration method’s adjoint explicit method

表2 常加速度伴隨顯式算法(n=3)穩定性和級數收斂性對h/T的限制的對比Tab.2 H/T’s limited field of constant acceleration method’s adjoint explicit method

表1列出了阻尼比分別為 0，0.05，0.1，而 k 分別為1，2，3，4，5時常加速度伴隨顯式算法的穩定區間。

表1顯式:n為奇數時算法的穩定性明顯優于n為偶數時算法的穩定性，建議采用奇數次級數展開。

h/T的取值受限于算法穩定條件，見表1，同時也受限于逆矩陣展開時的收斂性條件，見式(16)。表2詳細列出了當n=3時顯式算法的穩定域大小與級數收斂條件域大小的區別。表2顯示，兩個不同要求對h/T的限制非常接近。

3 結論

廣義線性加速度法求解的耦聯線性方程組，其系數矩陣由結構的剛度陣、質量陣/時間步長的平方、阻尼矩陣/時間步長的一次方等三項線性組合而成，隨著時間步長趨于零，質量陣權重加大，而剛度陣與阻尼陣的權重減少。據此，本文構造了以質量陣項為主，剛度陣、阻尼陣組合相關項為攝動項，對系數矩陣的逆矩陣在單位陣附近進行級數展開的顯式算法。詳盡分析了常加速度法方法的顯式化過程。公式推導中，逆矩陣級數展開的收斂條件與算法本身穩定性條件均被滿足，兩者對h/T提出的限制域很接近。級數展開奇數項算法的穩定性明顯優于級數展開偶數項算法的穩定性。

分析結果表明，無條件穩定的常加速度隱式算法，通過求逆、級數展開的步驟后得到的顯式算法，不再維持原隱式算法的無條件穩定性，即使級數展開很多項，也不能明顯改善顯式算法的穩定性能。對于常加速度法的伴隨顯式算法，逆矩陣的展開項達到3時，即可使顯式算法達到相應隱式算法的精度。

本文給出了隱式算法和顯式算法之間的聯系，提供了一種構造顯式算法的新思路:即可按高精度的隱式算法構造相應精度的顯式算法。換言之，隱式算法精度有多高，對應的顯式算法精度即能達到多高。本文這種多項式加速度顯式算法的計算精度與精細積分方法［7－9］相差不多，而計算效率卻高一個數量級。

［1］李小軍，廖振鵬，杜修力.有阻尼體系動力問題的一種顯式差分解法［J］.地震工程與工程振動，1992，12(4):74－79.

［2］李小軍，廖振鵬.非線性結構動力方程求解的顯式差分格式的特性分析［J］.工程力學，1993，10(3):141－146.

［3］杜修力，王進廷.阻尼彈性結構動力計算的顯式差分法［J］.工程力學，2000，17(5):37－43.

［4］王進廷，杜修力.有阻尼體系動力分析的一種顯式差分法［J］.工程力學，2002，19(3):109－112.

［5］張曉志，程 巖，謝禮立.結構動力反應分析的三階顯式方法［J］.地震工程與工程振動，2002，22(3):1－8.

［6］陳學良，金 星，陶夏新.求解加速度反應的顯式積分格式研究［J］.地震工程與工程振動，2006，26(5):60－67.

［7］鐘萬勰.暫態歷程的精細計算方法［J］.計算力學學報，1995，1:1 －6.

［8］譚述君，鐘萬勰.非齊次動力方程Duhamel項的精細積分［J］.力學學報，2007，3:374 －380.

［9］高小科，鄧子辰，黃永安.基于三次樣條插值的精細積分法［J］.振動與沖擊，2007，26(9):75－82.