陳軍委

(1. 上海航天控制技術(shù)研究所·上?！?01109；2.上海市空間智能控制技術(shù)重點實驗室·上?！?01109)

0 引 言

對實際動力學(xué)系統(tǒng)特性的分析要求有效的數(shù)值仿真求解算法，尤其要求在采用大步長時仍然具有高精度的算法[1]。對于不存在或可以忽略高頻的結(jié)構(gòu)動力學(xué)問題，普遍使用隱式積分方法(如Newmark、Wilson θ和Houbolt方法等)[2]或顯式方法(如Runge-Kutta法)進行瞬態(tài)分析和求解。通常，這些方法對于線性系統(tǒng)是準(zhǔn)確和穩(wěn)定的。例如，Newmark積分對于線性問題具有二階精度，并且無條件穩(wěn)定。然而，當(dāng)將上述方法應(yīng)用于非線性動力學(xué)系統(tǒng)時，為獲得精確和穩(wěn)定的數(shù)值解，選取的步長與系統(tǒng)最小周期、外部激勵周期和系統(tǒng)非線性項的變化周期相比，必須足夠小。因此，這些方法對于非線性動力學(xué)系統(tǒng)的求解效率不高，甚至容易失效[3]。

為此，研究者開發(fā)了面向非線性動力學(xué)系統(tǒng)的眾多數(shù)值方法，比如引入數(shù)值阻尼[4]，通過Lagrange乘子強制能量守恒[5]，以及采用算術(shù)能量守恒[6-7]等方法。Bathe等人提出了一種復(fù)合隱式積分算法，用于處理非線性問題[8-9]和線性剛性問題[10]。與上述富有成效的工作不同，本文作者采用了基于積分方程的方法。通常，積分格式會更具優(yōu)勢。理論上講，積分算子比微分算子更容易實現(xiàn)處理，采用積分格式能更輕易推斷出解的特性。事實上，積分格式還可以在某些情形下縮減對象的維度[11]。

具體地，采用精細積分方法(precise integration method，PIM)[12]，指數(shù)矩陣的計算可以達到計算機的精度。Zhong、Williams和Lin等人指出，基于不可分割步長，精細積分可以獲得無條件的穩(wěn)定性[12-14]。而Wang等人認為，雖然在考慮子步長的高階項截斷時，精細積分是條件穩(wěn)定的，但大多數(shù)離散模型可以輕易滿足這種穩(wěn)定條件[15]。此外，一階對角Pade近似對于精細積分也是無條件穩(wěn)定的[16]，因此，精細積分由于具備高精度和穩(wěn)定性，可被應(yīng)用于處理非線性動力學(xué)問題。

在執(zhí)行精細積分時，不可避免地需要計算關(guān)于非線性項和外部激勵的Duhamel積分。Zhong[13]和Lin等人[14]給出了外部激勵為線性、多項式、正弦和指數(shù)方程等形式或其組合形式時的積分精確計算公式。然而，這些計算公式并不適用于系統(tǒng)系數(shù)矩陣奇異或接近奇異的情況。Gu等人[17]基于多維展開的方法提出了一種修正的精細積分方法，來避免矩陣求逆。由此，原有的非齊次方程將轉(zhuǎn)化為齊次方程，但同時計算量也極大地增加了。Wang等人[15]通過顯式Gauss求積處理了Duhamel積分，但其準(zhǔn)確性極大程度地依賴于Gauss求積的點數(shù)。Li等人[18]在每一個步長內(nèi)使用了Lagrange三次插值近似非線性項，但仍然需要執(zhí)行指數(shù)矩陣的求逆操作。

本文基于精細積分方法，針對非線性動力學(xué)系統(tǒng)，提出了精細逐塊積分方法(precise block-by-block integration method，PBIM)。此方法高度精確，在使用大步長時也可保持穩(wěn)定，無須矩陣求逆便可以處理非線性項。本文余下內(nèi)容依次為：在第二部分，描述了精細逐塊積分方法，以及其精度和穩(wěn)定性分析；在第三部分，列舉了三個數(shù)值算例以驗證本方法的有效性；最后，在第四部分給出了有關(guān)結(jié)論。

1 精細逐塊積分方法

1.1 方法描述

具有時變剛度的二階動力系統(tǒng)通常以矩陣形式表示為

(1)

不失一般性，K(t)可以劃分為定常部分和時變部分，即K(t)=K1+K2(t)。式(1)由此可以變換為如下的一階形式

(2)

其中，Ip和0p分別是p×p階單位陣和零矩陣。

式(2)可以改寫為

(3)

其中，

令n=2p。顯然，A是n×n階常矩陣，g(t,z)是關(guān)于時間、狀態(tài)向量和外部激勵的n維向量方程。需要指出的是，式(3)是相比于式(1)更為通用的形式，雖然前者由后者推導(dǎo)而來。下文將以式(3)為對象進行討論。

將式(3)進行積分，得到第二類非線性Volterra積分方程，即

(4)

由于式(3)源于實際動力系統(tǒng)，所以認為其積分核滿足Lipschitz條件，即系統(tǒng)存在唯一解[11]。以步長h為離散時間間隔，記離散時間點的狀態(tài)為zi(i=0,1,2,…)，即zi=z(ti)。定義時間間隔[t2m,t2m+2](m=0,1,2,…)為一間隔塊，則在每一間隔塊內(nèi)，得

(5)

(6)

對式(5)～(6)應(yīng)用Simpson法則，得

(7)

(8)

由此，Duhamel積分可由精細積分求取[12]。在整個計算流程中，指數(shù)矩陣只需被計算一次。以eAh為例，應(yīng)用指數(shù)函數(shù)的附加定理，即2N算法，得

eAh=T2N

(9)

(10)

其中

(11)

τ非常小，故將前4項或前5項展開，便可獲得足夠高的精度。為避免出現(xiàn)運算圓整誤差，可將不變部分I和增量部分Ta分別進行存儲。將eAh進行分解，代入式(10)，可得

T2N>=>(T2)2N-1=(I+2Ta+Ta·Ta)2N-1

(12)

最終，經(jīng)N次形如Ta?2Ta+Ta·Ta的遞歸，eAh的數(shù)值結(jié)果逼近了計算機的運算精度。

需要指出的是，式(7)～(8)要同時求解未知量z2m+1和z2m+2。在此，可以應(yīng)用任何代數(shù)方程求解算法。例如，可以使用Newton迭代下山法以避免矩陣求逆，或采用同倫連續(xù)方法以實現(xiàn)全局收斂的效果[20]。

1.2 穩(wěn)定性

積分方法的數(shù)值穩(wěn)定性條件是指，對于單自由度系統(tǒng)，其迭代放大矩陣的譜半徑小于或等于1[21]。為考察精細逐塊積分方法的數(shù)值穩(wěn)定性，需考慮如下的線性齊次算例

(13)

其中，ω>0,ξ≥0。

采用Zhong和Williams提出的變換[12]，式(13)可寫為

(14)

其中，

應(yīng)用式(7)～(8)，可得

(15)

相應(yīng)的放大矩陣為

D=eAh

(16)

顯然，放大矩陣D的特征值λ滿足

|λI-eAh|=0

(17)

在不考慮Taylor展開截斷誤差時，求得特征值為

(18)

顯然，對于ξ≥0，特征值λ常滿足|λ|≤1。當(dāng)考慮Taylor展開截斷誤差時，如文獻[15, 22]所指，子步長τ趨近集中。Wang[15]指出，在這種情形下，截斷階數(shù)q=4、5和二分階數(shù)N=20對于大多數(shù)結(jié)構(gòu)而言均可保證穩(wěn)定性，只有在選取相對系統(tǒng)無阻尼固有周期的大子步長時，才會出現(xiàn)不穩(wěn)定的結(jié)果。

1.3 精度

簡便起見，考慮形如g(t,z(t))=G(t)z(t)的線性案例。對于初值問題，定常系統(tǒng)的精細積分數(shù)值解接近計算機的精度[13]。因此，精細逐塊積分的誤差滿足

(19)

(20)

其中

為進一步驗證以上結(jié)論，考慮如下的單自由度線性系統(tǒng)

(21)

其中，ω=1.0，ξ=0.05，其解析解為

x=e-ξt(c1cosηt+c2sinηt)+s1cosπt+s2sinπt

(22)

其中，

表1對比了精細逐塊積分方法、精細積分方法[12]的數(shù)值解和系統(tǒng)的解析解。當(dāng)步長為0.01s時，相比于解析解，精細逐塊積分方法可以精確至6位有效數(shù)字，高于精細積分方法。表2列出了精細逐塊積分方法在不同步長下的誤差?？梢钥闯?，當(dāng)步長減半時，誤差約縮減為1/16，這與上述收斂階數(shù)的結(jié)論非常吻合。

表1 單自由度系統(tǒng)的數(shù)值解(0s至1s)

表2 不同步長下精細逐塊積分方法的誤差

2 數(shù)值算例

為驗證精細逐塊積分方法的有效性，本部分給出了3個算例。

2.1 算例I

考慮兩自由度Fermi-Pasta-Ulam問題[23]。如圖1所示，系統(tǒng)包含一剛性線性彈簧和兩個非線性軟彈簧。本系統(tǒng)存在高頻振蕩，其Hamilton函數(shù)為

(23)

其中，qi和pi(i=1,2)分別為廣義坐標(biāo)和廣義動量矩。相應(yīng)的Hamilton正則方程為

(24)

圖1 兩自由度Fermi-Pasta-Ulam問題Fig.1 The 2-DOF Fermi-Pasta-Ulam problem

令ω=50，用精細逐塊積分方法求解式(24)。設(shè)置絕對誤差和相對誤差均為10-14，將Matlab求解器ode15s的求解結(jié)果作為參考解。表3列出了經(jīng)精細逐塊積分方法、四階Runge-Kutta法(RK4)和Newmark方法解出的q1值?？梢钥闯?，精細逐塊積分方法在步長設(shè)置為0.001s時具有9或10位有效數(shù)字的精度；當(dāng)步長增大為0.005s時，精細逐塊積分方法仍然具有7位有效數(shù)字的精度。然而，在設(shè)置步長為0.001s時，RK4方法只有6位有效數(shù)字的精度；Newmark 方法直到步長小至0.0001s時才有所收斂，并且只有1或2位有效數(shù)字的精度。

表3 經(jīng)ode15s、精細逐塊積分方法、RK4法和Newmark方法求得的q1值

圖2繪出了上述3種方法從0至80s的時間歷程響應(yīng)。精細逐塊積分方法和RK4法的步長為0.01s，Newmark方法的步長為0.0001s。顯然，Newmark方法在步長小至0.0001s時仍然不穩(wěn)定，而精細逐塊積分方法和RK4法可以獲得穩(wěn)定解。然而，如表3所示，雖然RK4法看起來長時穩(wěn)定，但不如精細逐塊積分方法精確。

2.2 算例II

考慮變系數(shù)非線性周期系統(tǒng)[24]

(25)

圖2 經(jīng)精細逐塊積分方法、RK4法和 Newmark方法求解的q1響應(yīng)Fig.2 Responses of q1 by the PBIM, the RK4 method and the Newmark method

其Galerkin解為

(26)

其精度可達10-5，可認為是參考解。分別應(yīng)用精細逐塊積分方法、RK4法和Newmark方法求解式(25)，將相應(yīng)的x1數(shù)值解的誤差繪制在圖3至圖5中。在圖3中，對精細逐塊積分方法依次選取步長h=0.1s和h=0.01s。顯然，即使步長增大了10倍，相應(yīng)的誤差仍保持在10-5。圖4展現(xiàn)了步長依次設(shè)置為h=0.1s和h=0.01s時，由RK4法獲得的x1的誤差。與精細逐塊積分方法不同，當(dāng)步長增大10倍時，誤差也自10-2至10-1放大了10倍。Newmark方法的表現(xiàn)最差。如圖5所示，Newmark方法要求非常小的步長，直到步長h=0.0001s才可保持數(shù)值穩(wěn)定。

圖3 不同步長下精細逐塊積分方法的x1的誤差Fig.3 Errors of x1 by the PBIM using different step sizes

圖4 不同步長下RK4法的x1的誤差Fig.4 Errors of x1 by the RK4 method using different step sizes

圖5 不同步長下Newmark方法的x1的誤差Fig.5 Errors of x1 by the Newmark method using different step sizes

2.3 算例III

最后，考慮如下的非線性剛性問題

(27)

圖6 精細逐塊積分方法和ode15s求解誤差比較(h=0.001s)Fig.6 Comparison of the errors by the PBIM and ode15s(h=0.001s)

3 結(jié) 論

本文提出了精細逐塊積分方法。作為隱式積分格式，此方法可以輕易滿足穩(wěn)定性條件，達到4階精度。因此，與Newmark方法和RK4法相比，對非線性動力系統(tǒng)應(yīng)用精細逐塊積分方法，可以在大步長的條件下獲得高精度和穩(wěn)定的數(shù)值解。此外，由于本方法不需要進行矩陣取逆運算，因此對具有奇異或接近奇異的系統(tǒng)矩陣的問題仍然有效。

附錄

定理A.1[11]令序列ξ0,ξ1,…滿足

其中，

得