結(jié)構(gòu)動(dòng)力方程求解中隱式格式向顯式格式的轉(zhuǎn)換

李常青，樓夢(mèng)麟，蔣麗忠

(1.中南大學(xué) 土木工程學(xué)院建筑工程系，長(zhǎng)沙 410000;2.同濟(jì)大學(xué) 土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092)

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

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

1 基本思路

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

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

將式 (1)兩 邊 同乘 矩 陣

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

級(jí)數(shù)展開(kāi)式(5)中的逆矩陣，得:

式(6)級(jí)數(shù)展開(kāi)的收斂條件為攝動(dòng)項(xiàng)范數(shù)不能大于1:

將式(6)代入式(5)，得無(wú)窮級(jí)數(shù)展開(kāi)的顯式算法:

式(8)中的級(jí)數(shù)展開(kāi)n次，即得級(jí)數(shù)截?cái)嗪蟮慕乒?

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

式(9)為對(duì)式(1)經(jīng)求逆、級(jí)數(shù)展開(kāi)、級(jí)數(shù)項(xiàng)截?cái)嗪蟮玫降慕?，本文稱式(9)為隱式算法(1)的伴隨顯式算法。下面分別詳細(xì)分析無(wú)條件穩(wěn)定的常加速度方法的顯式化過(guò)程和條件穩(wěn)定的線性加速度方法的顯式化過(guò)程。

2 常加速度方法的顯式化

2.1 公式推導(dǎo)

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

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

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

其中:

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

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

求解式(15)，得:

令 ζ=0，得:

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

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

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

2.2 精度分析

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

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

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

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

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

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

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

2.3 穩(wěn)定性分析

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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