李常青，樓夢麟，蔣麗忠

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

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

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

1 基本思路

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

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

將式 (1)兩 邊 同乘 矩 陣

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

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

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

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

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

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

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

2 常加速度方法的顯式化

2.1 公式推導(dǎo)

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

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

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

其中:

為任何一種矩陣范數(shù)。

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

求解式(15)，得:

令 ζ=0，得:

化為h與T(結(jié)構(gòu)周期)的形式為:

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

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

2.2 精度分析

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

當(dāng)k=1時，位移算法精度為二階，位移誤差項為:

當(dāng)k=2時，位移算法精度為三階，位移誤差項為:

當(dāng)k=3時，位移算法精度為三階，位移誤差項為:

當(dāng)k=4時，位移算法精度為三階，位移誤差項為:

當(dāng)k=5時，位移算法精度為三階，位移誤差項為:

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

2.3 穩(wěn)定性分析

算法的穩(wěn)定分析，通過計算單自由度系統(tǒng)得到傳遞矩陣，形式如下:

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

算法的穩(wěn)定性要求:

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

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

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

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

表2 常加速度伴隨顯式算法(n=3)穩(wěn)定性和級數(shù)收斂性對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時常加速度伴隨顯式算法的穩(wěn)定區(qū)間。

表1顯式:n為奇數(shù)時算法的穩(wěn)定性明顯優(yōu)于n為偶數(shù)時算法的穩(wěn)定性，建議采用奇數(shù)次級數(shù)展開。

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

3 結(jié)論

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

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

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

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

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

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

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

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

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

［7］鐘萬勰.暫態(tài)歷程的精細(xì)計算方法［J］.計算力學(xué)學(xué)報，1995，1:1 －6.

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

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