一種多自由度阻尼體系動(dòng)力問(wèn)題的高階分析方法

李鴻晶，楊筱朋，梅雨辰

(南京工業(yè)大學(xué)工程力學(xué)研究所，南京 211816)

結(jié)構(gòu)動(dòng)力響應(yīng)分析是認(rèn)識(shí)動(dòng)態(tài)系統(tǒng)性態(tài)和行為的基礎(chǔ)性工作。傳統(tǒng)的動(dòng)力響應(yīng)分析方法可以大致分為兩種類(lèi)型：一類(lèi)是基于坐標(biāo)變換[1?2]的方法，如振型疊加法、頻域法等，一般適用于線(xiàn)彈性體系和經(jīng)典阻尼體系的動(dòng)力問(wèn)題；另一類(lèi)為直接積分法[3?6]，其直接對(duì)體系運(yùn)動(dòng)微分方程進(jìn)行求解，以Newmark-β 法等逐步積分法為代表，可以解決非線(xiàn)性體系和非經(jīng)典阻尼體系的動(dòng)力分析問(wèn)題。這些傳統(tǒng)動(dòng)力分析方法中有些算法具備無(wú)條件穩(wěn)定及可控算法阻尼等特性[7]，多采用低階格式，精度一般不超過(guò)二階。也有一些算法采用了高階格式，精度能達(dá)到三階以上，但所需的計(jì)算量及存儲(chǔ)量要增大不少[8]。有別于傳統(tǒng)動(dòng)力方法，精細(xì)積分法(precise time-step integration method,PTIM)[9?18]則利用體系動(dòng)力響應(yīng)積分解，通過(guò)計(jì)算矩陣指數(shù)獲得動(dòng)力響應(yīng)。由于結(jié)構(gòu)大型化、復(fù)雜化以及精細(xì)化分析等方面的需要，發(fā)展高效率的動(dòng)力分析方法始終是一種非?，F(xiàn)實(shí)的需求。

微分求積 (differential quadrature, DQ)法[19?22]是近年發(fā)展起來(lái)的一種求解微分方程(組)的數(shù)值方法，F(xiàn)ung[23?25]首先將其應(yīng)用于初值問(wèn)題，其后一些學(xué)者也開(kāi)展了相關(guān)研究[26?32]，逐漸形成了動(dòng)力DQ 分析方法。該方法數(shù)值穩(wěn)定性好，計(jì)算簡(jiǎn)便且精度高，然而由于其以位移響應(yīng)為基本未知量來(lái)離散動(dòng)力方程，所得數(shù)值格式為隱式格式，因而需要求解方程。另外，目前已提出的這類(lèi)動(dòng)力方法因穩(wěn)定性和精度方面的要求多采用不等步距的時(shí)間離散方案，求解均勻布置的采樣時(shí)刻處的響應(yīng)值時(shí)需要進(jìn)行插值，對(duì)其計(jì)算效率造成了一定的限制。與傳統(tǒng)動(dòng)力響應(yīng)DQ 分析方法不同，在文獻(xiàn)[33 ? 34]中，以單自由度體系為對(duì)象提出了一種計(jì)算Duhamel 積分的高精度算法，即積分求微方法(integral differentiation method, IDM)。該算法為高階方法，動(dòng)力響應(yīng)計(jì)算僅需要完成一次方陣與向量的乘法運(yùn)算，一次計(jì)算即可同時(shí)獲得多個(gè)時(shí)刻的響應(yīng)值。整個(gè)過(guò)程無(wú)需求解方程，且可采用等時(shí)間步距的計(jì)算方案，不必通過(guò)插值即可直接獲得各等距分布時(shí)點(diǎn)的響應(yīng)。積分求微法可視為非齊次動(dòng)力方程的一種數(shù)值解法，如果將其向多自由度體系推廣，則將引入矩陣指數(shù)。由于積分求微法無(wú)需進(jìn)行數(shù)值積分計(jì)算，與單自由度體系動(dòng)力響應(yīng)分析過(guò)程相比較，多自由度體系動(dòng)力分析僅需要處理矩陣指數(shù)的計(jì)算。

本文將逆用DQ 原理，建立多自由度阻尼體系動(dòng)力時(shí)程分析的積分求微算法，以豐富動(dòng)力響應(yīng)分析的技術(shù)手段。

1 多自由度體系狀態(tài)響應(yīng)

動(dòng)態(tài)荷載向量p(t)激勵(lì)下的多自由度阻尼體系運(yùn)動(dòng)方程寫(xiě)為：

狀態(tài)方程式(3)中的系數(shù)矩陣H 僅包含體系剛度、質(zhì)量和阻尼的特性，因而它是反映體系動(dòng)態(tài)特性的。對(duì)于線(xiàn)彈性體系動(dòng)力分析，體系特性矩陣H 將是時(shí)不變的常數(shù)矩陣。而方程的非齊次項(xiàng)中含有荷載向量p(t)，它描述了環(huán)境對(duì)體系施加的動(dòng)態(tài)作用的特性。

式(3)的解可采用 Duhamel 積分的形式表示為：

式中， eHt為矩陣指數(shù)。

形如式(5)的方程解中包含了矩陣指數(shù)及其積分式，對(duì)于實(shí)際動(dòng)力問(wèn)題，試圖利用式(5)獲得體系動(dòng)力響應(yīng)y(t)的解析解答幾乎是不可能的，但可以通過(guò)數(shù)值方法得到y(tǒng)(t)的近似結(jié)果。

2 動(dòng)力響應(yīng)積分求微分析

利用式(5)計(jì)算y(t)的難點(diǎn)在于含有矩陣指數(shù)項(xiàng)的卷積的計(jì)算。將該卷積做如下的分離：

顯見(jiàn)，求得動(dòng)力響應(yīng)y(t)的關(guān)鍵是處理好兩個(gè)問(wèn)題：一是矩陣指數(shù)T(t)的計(jì)算；二是計(jì)算出變上限定積分函數(shù)向量D(t)。兩者都要求建立起適合的數(shù)值算法。

考慮動(dòng)力持時(shí)中的任意時(shí)段[ti, tj]，該時(shí)段包含有ρ 個(gè)時(shí)步，其中ti、tj均選擇位于采樣點(diǎn)處。如圖1 所示，假定采取等距采樣，采樣周期為δ，則時(shí)段[ti, tj]的長(zhǎng)度應(yīng)為：

圖 1 動(dòng)態(tài)響應(yīng)待求解時(shí)段Fig. 1 The solved time interval for dynamic response

2.1 矩陣指數(shù)Tk 的計(jì)算

根據(jù)矩陣指數(shù)的Taylor 級(jí)數(shù)定義，矩陣指數(shù)Tk可表示為：

2.2 定積分 Dk 的計(jì)算

由式(8)可見(jiàn)，D(t)對(duì)時(shí)間t 的導(dǎo)數(shù)可表示為：

觀(guān)察上述求解動(dòng)力響應(yīng)的過(guò)程，可以發(fā)現(xiàn)除了在確定系統(tǒng)特性矩陣H 時(shí)需要對(duì)質(zhì)量矩陣m 求逆，以及在確定轉(zhuǎn)換矩陣B 時(shí)需要對(duì)Vandermonde矩陣進(jìn)行求逆計(jì)算外，其余過(guò)程都不需要矩陣求逆或者求解線(xiàn)性方程組。而對(duì)于結(jié)構(gòu)體系可采用集中質(zhì)量模型，即質(zhì)量矩陣m 為對(duì)角陣，其逆矩陣可以方便地求出。至于Vandermonde 矩陣求逆計(jì)算，可以事先算出逆矩陣結(jié)果，考慮到實(shí)際動(dòng)力分析時(shí)待求解時(shí)段 [0, ?tij]內(nèi)包含的時(shí)步數(shù)ρ 不可能很大 (一般取 ρ=10～15)，故可將轉(zhuǎn)換矩陣B 制成表格，動(dòng)力響應(yīng)分析時(shí)直接取用即可。鑒于此，按照式(32)計(jì)算積分矩陣Φ 時(shí)僅需做矩陣乘法運(yùn)算，這相當(dāng)于是一種顯式算法。

3 算法步驟

依據(jù)上述推導(dǎo)，可將多自由度阻尼體系動(dòng)力時(shí)程積分求微法總結(jié)為下述算法：

4 數(shù)值穩(wěn)定性分析

對(duì)于線(xiàn)性體系動(dòng)力方法的數(shù)值穩(wěn)定性，只需研究如下的單自由度體系方程：

阻尼一般對(duì)穩(wěn)定性起到有利作用，因此研究穩(wěn)定性時(shí)多不考慮阻尼的影響。以下分析都是在消除阻尼項(xiàng)后開(kāi)展的。觀(guān)察式(12)的動(dòng)力響應(yīng)數(shù)值格式，可得時(shí)段初始點(diǎn)與末點(diǎn)的遞推關(guān)系式：

表 1 L 不同取值時(shí)的穩(wěn)定極限值Table 1 The limit of stability with different L

由此可得出，本文方法一般為有條件穩(wěn)定方法，當(dāng)L=3、4、7、8 時(shí)，穩(wěn)定極限較大，穩(wěn)定性相對(duì)較好，再綜合考慮計(jì)算量的因素，取L=4 比較合適。由于在精細(xì)積分法中λ 一般選20，即時(shí)步小段的長(zhǎng)度 η ≈ 9.54×10?7δ，因而只需要各節(jié)點(diǎn)間距滿(mǎn)足：

即可穩(wěn)定地進(jìn)行計(jì)算。顯然，通常情況下該穩(wěn)定極限條件是極易得到滿(mǎn)足的，當(dāng)取L=4 時(shí)本文方法時(shí)步大小的取值幾乎不受穩(wěn)定性的限制。

5 算例分析

通過(guò)兩個(gè)算例闡釋本文方法的特性，第一個(gè)算例用來(lái)驗(yàn)證本文方法的可靠性，包括收斂性和精度；第二個(gè)算例重點(diǎn)闡釋本文方法的適用性。

5.1 兩自由度體系算例

選擇文獻(xiàn)[12]中的兩自由度體系算例，其特性矩陣分別為：

采用本文的IDM 法計(jì)算該體系的位移響應(yīng)和速度響應(yīng)，并以精確解為標(biāo)準(zhǔn)驗(yàn)證計(jì)算結(jié)果。為了清晰地展現(xiàn)本文方法的收斂性及精確性，將時(shí)段內(nèi)時(shí)步數(shù)ρ 取為10，時(shí)段長(zhǎng)度Δt 從1.0 s 開(kāi)始依次減為原來(lái)的1/2 時(shí)，即分別取Δt=1 s、1/2 s、1/22s、1/23s 時(shí)，計(jì)算體系在 t=0.2 s、0.4 s、0.6 s、0.8 s、1.0 s 處的響應(yīng)值。計(jì)算結(jié)果同精確解的最大絕對(duì)誤差Err 與Δt 的關(guān)系曲線(xiàn)如圖2 和圖3 所示，其中橫坐標(biāo)和縱坐標(biāo)皆采用以2 為底的對(duì)數(shù)坐標(biāo)。

觀(guān)察圖2 和圖3 可以發(fā)現(xiàn)，隨著時(shí)段長(zhǎng)度Δt 的減小，本文方法計(jì)算結(jié)果的誤差亦迅速減小。當(dāng)Δt 縮小到1/23s 時(shí)，無(wú)論是位移響應(yīng)還是速度響應(yīng)，其最大絕對(duì)誤差皆不超過(guò)10?14數(shù)量級(jí)，計(jì)算結(jié)果已非常精確，具備較完美的收斂效果。從整個(gè)時(shí)步的取值過(guò)程來(lái)看，本文IDM 法的收斂速度非?？?，幾乎是時(shí)步長(zhǎng)度每縮小一倍計(jì)算結(jié)果的絕對(duì)誤差就會(huì)減小約2?10量級(jí)。進(jìn)一步觀(guān)察可以直觀(guān)地發(fā)現(xiàn)在對(duì)數(shù)坐標(biāo)系下最大絕對(duì)誤差與時(shí)步長(zhǎng)度近似呈線(xiàn)性關(guān)系，這符合數(shù)值計(jì)算原理中誤差與時(shí)步長(zhǎng)度的關(guān)系：

式中：絕對(duì)誤差Err 對(duì)于函數(shù)可以用其范數(shù)來(lái)近似替代；q 為數(shù)值格式的收斂精度階數(shù)。

圖 2 不同Δt 下體系位移響應(yīng)最大絕對(duì)誤差Fig. 2 The maximum absolute error of the displacement response of the system with different Δt

圖 3 不同Δt 下體系速度響應(yīng)最大絕對(duì)誤差Fig. 3 The maximum absolute error of the velocity response of the system with different Δt

為合理有效地估計(jì)本文IDM 方法的收斂階數(shù)，分別選取整體誤差指標(biāo)——最大絕對(duì)誤差以及某個(gè)時(shí)刻(如t=0.4 s)時(shí)的局部絕對(duì)誤差兩個(gè)指標(biāo)作為式(41)中的誤差項(xiàng)Err，采用最小二乘法估算出體系位移響應(yīng)和速度響應(yīng)的收斂階數(shù)q，結(jié)果列于表2 和表3。

表 2 響應(yīng)收斂階數(shù)q 估計(jì)(整體誤差指標(biāo))Table 2 Evaluation of convergence order q for responses based on whole error criterion

表 3 響應(yīng)收斂階數(shù)q 估計(jì)(局部誤差指標(biāo))Table 3 Evaluation of convergence order q for responses based on local error criterion

由表2 和表3 的統(tǒng)計(jì)結(jié)果，無(wú)論是采用整體誤差指標(biāo)還是局部誤差指標(biāo)，各質(zhì)點(diǎn)的位移響應(yīng)和速度響應(yīng)的收斂階數(shù)皆在10～11，即該數(shù)值格式對(duì)于位移響應(yīng)和速度響應(yīng)的收斂精度至少為10 階。而加速度響應(yīng)是依據(jù)位移響應(yīng)和速度響應(yīng)由運(yùn)動(dòng)方程求得的，因而其具備與位移響應(yīng)、速度響應(yīng)相同的收斂精度。另一方面，從理論上來(lái)看，本文IDM 數(shù)值格式的收斂精度主要由矩陣指數(shù)的計(jì)算精度和DQ 法的逼近精度所共同決定。而矩陣指數(shù)的計(jì)算由精細(xì)積分實(shí)現(xiàn)，其數(shù)值逼近的時(shí)間節(jié)點(diǎn)間距非常小，計(jì)算精度相較于非齊次項(xiàng)DQ 逼近要高許多。因此可以認(rèn)為，本文IDM法的收斂精度主要由DQ 法的計(jì)算精度決定，它直接取決于時(shí)段內(nèi)的時(shí)步數(shù)ρ 的取值。由DQ 原理可知，在各時(shí)段Δt 內(nèi)對(duì)其積分函數(shù)D(τ)逼近的絕對(duì)誤差函數(shù)為：

從式(42)來(lái)看，其誤差函數(shù)是Δt 的ρ+1 階無(wú)窮小，即本文格式具有ρ+1 階代數(shù)精度。當(dāng)取ρ=10 時(shí)，本文方法收斂精度階數(shù)的理論值為11。考慮到計(jì)算機(jī)舍入誤差及矩陣指數(shù)的影響，實(shí)際計(jì)算時(shí)收斂精度會(huì)略微小于理論值。即實(shí)際計(jì)算時(shí)若取ρ=10，本文方法至少能達(dá)到10 階收斂精度，這同表2 和表3 中估計(jì)的結(jié)果相符合。一般地，當(dāng)計(jì)算時(shí)段內(nèi)包含的時(shí)步數(shù)為ρ 時(shí)，理論上本文方法有ρ+1 階代數(shù)精度，實(shí)際計(jì)算時(shí)至少能達(dá)到ρ 階精度。

下面將本文IDM 法與傳統(tǒng)動(dòng)力方法進(jìn)行比較，選擇的傳統(tǒng)方法包括Newmark 常平均加速度法 (取 γ=0.5、β=0.25)、Wilson-θ 法 (取 θ=1.4)和精細(xì)積分法(PTIM)。將所有這些方法的計(jì)算結(jié)果列于表4，同時(shí)給出了精確解作為比較的標(biāo)準(zhǔn)。

表 4 不同動(dòng)力方法計(jì)算的位移響應(yīng)結(jié)果Table 4 Computational results of displacement response by different dynamic methods

續(xù)表 4

觀(guān)察表4 可以發(fā)現(xiàn)，當(dāng)時(shí)段長(zhǎng)度取為Δt=1.0 s時(shí)，劃分10 個(gè)時(shí)步，本文積分求微法的計(jì)算結(jié)果與精確解已非常接近，最大相對(duì)誤差在2%以?xún)?nèi)。

注意到本文IDM 法為高階方法，此時(shí)待求解時(shí)段[ti, tj]內(nèi)含有的10 個(gè)時(shí)步的步距均為δ=0.1 s，而傳統(tǒng)的 Newmark-β 法和 Wilson-θ 法在步距取為0.05 s 時(shí)位移響應(yīng)的最大相對(duì)誤差已超過(guò)了10%，在取為0.02 s 時(shí)，也超過(guò)了2%。也就是說(shuō)，在步距是這兩種傳統(tǒng)逐步積分方法5 倍的條件下，本文IDM 法的計(jì)算結(jié)果仍然比其要精確些，這直觀(guān)顯示了本文方法高精度的特點(diǎn)。

5.2 彈簧-質(zhì)量體系算例

以圖4 所示的串聯(lián)彈簧質(zhì)量體系為例，取自由度數(shù)N=100。為了使體系固有頻率覆蓋范圍廣，彈簧剛度一半設(shè)置為較剛，另一半設(shè)置相對(duì)較柔，其特性參數(shù)見(jiàn)表5。每個(gè)質(zhì)點(diǎn)都受到相同的水平簡(jiǎn)諧荷載激勵(lì)：

圖 4 100 自由度彈簧質(zhì)量體系Fig. 4 Schematic plot of 100-degree of freedom spring-mass system

表 5 彈簧-質(zhì)量體系的基本特性Table 5 Properties of the spring-mass system

該100 自由度體系的基本頻率為0.98 rad/s，最高階頻率為1094.92 rad/s，從低頻覆蓋到極高頻，便于全面地考察本文方法的數(shù)值穩(wěn)定性。取截?cái)囗?xiàng)數(shù)L=4、時(shí)段分段數(shù)ρ=10，采用本文IDM法計(jì)算體系每隔0.01 s 的動(dòng)力響應(yīng)。選用不同的時(shí)段長(zhǎng)度Δt，計(jì)算結(jié)果與顯式中心差分法(取Δt=0.001 s)進(jìn)行比較。圖5 和圖6 分別為最右端質(zhì)點(diǎn)(DOF 100)和第50 個(gè)質(zhì)點(diǎn)(DOF 50)的位移時(shí)程曲線(xiàn)。這里需要說(shuō)明的是，對(duì)于該體系中心差分法的穩(wěn)定極限為Δt≤0.0018 s，因而采用Δt=0.001 s可得到穩(wěn)定的計(jì)算結(jié)果，且它是直接獲得0.01 s 時(shí)間間隔響應(yīng)所允許的最大步距。

圖 5 體系第100 自由度位移響應(yīng)對(duì)比Fig. 5 Comparison of the displacement responses of the DOF 100 in this system

圖 6 體系第50 自由度位移響應(yīng)對(duì)比Fig. 6 Comparison of the displacement responses of the DOF 50 in this system

觀(guān)察圖5 和圖6 的動(dòng)力響應(yīng)計(jì)算結(jié)果，IDM法在選擇大得多的步距情形下依然能夠獲得穩(wěn)定的計(jì)算結(jié)果。例如，將時(shí)步長(zhǎng)度放大至0.1 s 時(shí)，利用IDM 法計(jì)算最右端質(zhì)點(diǎn)的位移響應(yīng)，仍然沒(méi)有任何失穩(wěn)放大現(xiàn)象產(chǎn)生。此情形下其時(shí)段長(zhǎng)度已是中心差分法穩(wěn)定極限的50 倍以上，在較高自振頻率的參與下依然能得到穩(wěn)定可靠的結(jié)果，這充分顯示了本文方法良好的數(shù)值穩(wěn)定性。

下面進(jìn)一步考察本文IDM 法的計(jì)算效率。為了更直觀(guān)地體現(xiàn)算法的執(zhí)行效率，表6 給出了一些情況下IDM 法與中心差分法(CDM)耗費(fèi)的計(jì)算時(shí)間及兩種方法計(jì)算時(shí)間的比值。其中所有的計(jì)算皆是在同一臺(tái)計(jì)算機(jī)上完成的，該計(jì)算機(jī)的架構(gòu)為英特爾Xeon(R)W3565 中央處理器，主頻3.20 GHz。各類(lèi)情況的計(jì)算時(shí)間均用符號(hào)CPU 來(lái)表示，且該時(shí)間是純計(jì)算時(shí)間，不包括數(shù)據(jù)輸出及單位換算的時(shí)間。

表 6 不同方法CPU 時(shí)間對(duì)比Table 6 Comparison of CPU time with different methods

觀(guān)察表6，對(duì)于IDM 法，當(dāng)時(shí)段長(zhǎng)度取為0.05 s 時(shí)其計(jì)算時(shí)間還不到中心差分法的40%，當(dāng)時(shí)步縮小到0.025 s 時(shí)也只有CDM 法的70%左右。通過(guò)對(duì)比各響應(yīng)時(shí)程計(jì)算結(jié)果，可以發(fā)現(xiàn)對(duì)于DOF 100 質(zhì)點(diǎn)的位移響應(yīng)，IDM (Δt=0.05 s)和IDM (Δt=0.025 s)的計(jì)算結(jié)果與CDM 的結(jié)果基本一致，而對(duì)于DOF 50 質(zhì)點(diǎn)，IDM (Δt=0.025 s)的計(jì)算結(jié)果與中心差分法基本一致，而IDM (Δt=0.05 s)的精度稍差一點(diǎn)。但即便如此，在IDM 法與CDM法精度大致接近的情況下，其計(jì)算時(shí)間依然比中心差分法要少30%左右，從而在時(shí)間復(fù)雜度上IDM 法要低于傳統(tǒng)的顯式中心差分法。另外，IDM(Δt=0.025 s)在整個(gè)持時(shí)內(nèi)的時(shí)間采樣點(diǎn)數(shù)量?jī)H為CDM(Δt=0.001 s)的 40%，雖然 IDM 法結(jié)構(gòu)特性矩陣的元素?cái)?shù)量是中心差分法的4 倍，但對(duì)于線(xiàn)性時(shí)不變體系而言?xún)H需存儲(chǔ)一次，且有近1/2 是零元素。相較于中心差分法，IDM 法憑借更少的時(shí)間采樣點(diǎn)數(shù)，隨著荷載持時(shí)的延長(zhǎng)其空間復(fù)雜度的優(yōu)勢(shì)必將愈加明顯。綜上所述，本文IDM 法在計(jì)算復(fù)雜度(包括時(shí)間復(fù)雜度和空間復(fù)雜度)上要低于傳統(tǒng)的顯式CDM 法，由于CDM 已是目前效率較高的動(dòng)力分析方法，表明本文IDM 法在計(jì)算效率方面具有一定優(yōu)勢(shì)。

6 討論與結(jié)論

根據(jù)DQ 原理，函數(shù)在指定點(diǎn)處的導(dǎo)數(shù)值可以用域內(nèi)一系列離散點(diǎn)處函數(shù)值的加權(quán)和來(lái)近似表示。將DQ 原理用于偏微分方程求解時(shí)，傳統(tǒng)上都是將域內(nèi)離散點(diǎn)處的函數(shù)值視為待求解的未知量，然后利用DQ 原理將各離散點(diǎn)處的導(dǎo)數(shù)值分別以函數(shù)值線(xiàn)性加權(quán)和替代并代入微分方程，從而將微分方程轉(zhuǎn)化為代數(shù)方程組，進(jìn)而獲得微分方程的近似解答。由于DQ 法高精度、低計(jì)算消耗的特點(diǎn)，這種基于DQ 原理的微分方程數(shù)值解法僅需付出不大的計(jì)算工作量即可獲得較高精度的數(shù)值解。對(duì)于時(shí)間相關(guān)問(wèn)題，傳統(tǒng)的基于DQ原理的動(dòng)力響應(yīng)分析方法都是沿用了這種求解思路，即將各離散時(shí)刻的速度響應(yīng)和加速度響應(yīng)分別表示為時(shí)段(步)內(nèi)全部離散時(shí)刻位移響應(yīng)的加權(quán)和。但是像動(dòng)力響應(yīng)Duhamel 積分解這樣的在被積函數(shù)中含有指數(shù)函數(shù)乘積項(xiàng)的定積分式，在積分變量分離后其導(dǎo)數(shù)即為被積函數(shù)自身，所以這種類(lèi)型卷積的導(dǎo)數(shù)是非常容易求得的。同傳統(tǒng)的DQ 法應(yīng)用不同，本文的積分求微法中離散點(diǎn)處的導(dǎo)數(shù)值是在應(yīng)用DQ 原理前就已經(jīng)求得的，這樣只需通過(guò)權(quán)系數(shù)逆矩陣進(jìn)行變換即可求得位移響應(yīng)向量，即動(dòng)力響應(yīng)計(jì)算僅需執(zhí)行矩陣乘法運(yùn)算即可一次性地求出時(shí)段內(nèi)的全部位移響應(yīng)值。而權(quán)系數(shù)矩陣的求逆計(jì)算只涉及Vandermonde矩陣求逆，這可以事先求出并在程序中存儲(chǔ)。因此，本文的動(dòng)力響應(yīng)積分求微方法相當(dāng)于顯式方法，其關(guān)鍵在于對(duì)DQ 原理的逆用。

經(jīng)過(guò)逆用DQ 原理后，動(dòng)力響應(yīng)分析只需解決矩陣指數(shù)計(jì)算即可完成。PTIM 是實(shí)現(xiàn)矩陣指數(shù)計(jì)算的一種優(yōu)秀算法，具體應(yīng)用到本文動(dòng)力響應(yīng)積分求微方法中，需要注意幾個(gè)問(wèn)題：

首先，本文動(dòng)力響應(yīng)積分求微法是按照時(shí)段進(jìn)行求解的，每個(gè)待求時(shí)段內(nèi)都包含有若干個(gè)時(shí)步。根據(jù)矩陣指數(shù)加法定理，在等距時(shí)步采樣即待求時(shí)段內(nèi)各時(shí)步的步距相同的前提下(事實(shí)上動(dòng)力分析一般都會(huì)采取這種等時(shí)步方案)，各個(gè)時(shí)步端點(diǎn)的響應(yīng)值才能表示成像式(19)那樣的矩陣指數(shù)遞推格式。這樣，當(dāng)采取等距時(shí)步方案時(shí)只需要計(jì)算出待求時(shí)段內(nèi)第一個(gè)時(shí)步對(duì)應(yīng)的矩陣指數(shù)，然后通過(guò)遞推公式就可以計(jì)算出時(shí)段內(nèi)所有時(shí)步的矩陣指數(shù)。顯然，矩陣指數(shù)精細(xì)積分計(jì)算僅僅被限制在第一個(gè)時(shí)步，即整個(gè)時(shí)段動(dòng)力響應(yīng)分析只要求一次精細(xì)積分計(jì)算，而不必采用PTIM計(jì)算多個(gè)矩陣指數(shù)。這是實(shí)施本文動(dòng)力響應(yīng)積分求微法的一個(gè)要點(diǎn)。也就是說(shuō)，在本文積分求微動(dòng)力分析方法中應(yīng)用PTIM 計(jì)算矩陣指數(shù)可以獲得較高的計(jì)算效率，關(guān)鍵在于動(dòng)力時(shí)程等距分布的采樣點(diǎn)同時(shí)也被選擇為DQ 離散節(jié)點(diǎn)，這樣就只需要在第一個(gè)時(shí)步內(nèi)應(yīng)用PTIM，即在每個(gè)求解時(shí)段內(nèi)只需使用一次PTIM 完成矩陣指數(shù)計(jì)算即可。

但是，使用DQ 法時(shí)一般會(huì)選擇不等距節(jié)點(diǎn)離散方案，采取等距節(jié)點(diǎn)方案會(huì)降低DQ 分析的精度。如果采取不等距時(shí)步方案，則勢(shì)必要求多次執(zhí)行矩陣指數(shù)的精細(xì)積分計(jì)算，而且還需要進(jìn)行插值才能獲得動(dòng)力分析結(jié)果，這將極大地增加計(jì)算工作量。所以，具體進(jìn)行動(dòng)力分析時(shí)就需要權(quán)衡計(jì)算的精度和效率。本文建議采取等時(shí)步方案，不僅因?yàn)檫@是動(dòng)力分析的通常選擇，而且可以將計(jì)算工作量控制在最低范圍。由于DQ 法的高精度屬性，即便采取等時(shí)步方案時(shí)也會(huì)取得相當(dāng)高的計(jì)算精度。這從本文給出的算例中可以看出。

其次，由于PTIM 計(jì)算矩陣指數(shù)時(shí)分為2 個(gè)步驟，需要計(jì)算經(jīng)過(guò)細(xì)分的很小步長(zhǎng)內(nèi)的矩陣指數(shù)，即首先需要計(jì)算出矩陣指數(shù) eHη，然后才能通過(guò)遞推得到需要求得的矩陣指數(shù) eHδ(δ 為時(shí)步的步距，本文方法中只要求對(duì) eHδ進(jìn)行精細(xì)積分計(jì)算)。在計(jì)算 eHη時(shí)，本文引入了秦九韶算法，將冪函數(shù)求和(被截?cái)嗟腡aylor 級(jí)數(shù))計(jì)算轉(zhuǎn)化為遞推格式，每一步遞推計(jì)算都只做一次矩陣乘法運(yùn)算，有效地降低了計(jì)算工作量和數(shù)據(jù)存儲(chǔ)空間。截?cái)嚯A數(shù)L 越高，采用秦九韶算法的收益就會(huì)越大。

最后，由于本文方法逆用DQ 原理，所以不要求進(jìn)行數(shù)值積分計(jì)算，但需要計(jì)算矩陣指數(shù)e?Hδ。在使用精細(xì)積分法時(shí)，矩陣指數(shù) e?Hδ的計(jì)算可以連同 eHδ一起計(jì)算出來(lái)。只要將細(xì)分后的矩陣指數(shù) eHη的被截?cái)嗉?jí)數(shù)表達(dá)式中的奇數(shù)項(xiàng)和偶數(shù)項(xiàng)分列為2 組，然后分別使用秦九韶算法進(jìn)行計(jì)算，最后通過(guò)兩組結(jié)果相加得到 eHη，相減得到e?Hη。這樣算法可同時(shí)計(jì)算出矩陣指數(shù) eHη和e?Hη，而沒(méi)有增加額外的計(jì)算工作量。其后遞推計(jì)算 e?Hδ的過(guò)程與 eHδ完全一致，可以方便地計(jì)算出來(lái)。

另外，在文獻(xiàn)[33]中針對(duì)單自由度體系實(shí)施積分求微分析時(shí)，待求時(shí)段內(nèi)各時(shí)刻的響應(yīng)值構(gòu)成一個(gè)向量，因而可以只通過(guò)執(zhí)行一次矩陣與向量的乘法運(yùn)算就能計(jì)算出時(shí)段內(nèi)的全部響應(yīng)。而多自由度體系有所不同，因?yàn)榇髸r(shí)段內(nèi)每個(gè)時(shí)刻的響應(yīng)都是一個(gè)向量，故時(shí)段內(nèi)所有時(shí)刻的響應(yīng)實(shí)際上構(gòu)成了一個(gè)矩陣，所以時(shí)段內(nèi)全部響應(yīng)需要通過(guò)矩陣與矩陣的乘法運(yùn)算才能實(shí)現(xiàn)。

關(guān)于本文積分求微法的穩(wěn)定性問(wèn)題，在文獻(xiàn)[33]中我們已經(jīng)證明其為無(wú)條件穩(wěn)定算法。但是對(duì)于多自由度體系，由于采用了PTIM 計(jì)算矩陣指數(shù)，所以該算法的穩(wěn)定性不僅取決于積分求微計(jì)算即逆用DQ 原理本身，還與PTIM 的穩(wěn)定性有關(guān)。關(guān)于PTIM 穩(wěn)定性問(wèn)題已有相關(guān)研究[35]，結(jié)論為PTIM 是有條件穩(wěn)定算法，但穩(wěn)定性條件極易滿(mǎn)足，對(duì)于常規(guī)工程結(jié)構(gòu)幾乎不會(huì)出現(xiàn)失穩(wěn)問(wèn)題。事實(shí)上，本文之所以只計(jì)算第一個(gè)時(shí)步的矩陣指數(shù)，待求時(shí)段內(nèi)的其他時(shí)步的矩陣指數(shù)都是通過(guò)遞推的方式獲得，也從另一方面說(shuō)明了PTIM 具有非常好的穩(wěn)定性。

綜上所述，執(zhí)行本文積分求微法的主要環(huán)節(jié)在于DQ 原理逆用和矩陣指數(shù)計(jì)算。通過(guò)引入矩陣指數(shù)及相應(yīng)算法，本文將結(jié)構(gòu)動(dòng)力時(shí)程分析的積分求微法推廣至多自由度情形，通過(guò)理論分析和數(shù)值試驗(yàn)，可得到如下結(jié)論：

(1) 結(jié)構(gòu)動(dòng)力時(shí)程分析的積分求微方法能夠適用于多自由度體系的直接積分，計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)是精確收斂的。該方法是一種高階時(shí)程分析方法，包含ρ 個(gè)時(shí)步的積分求微法理論上具有(ρ+1)階代數(shù)精度，與傳統(tǒng)的低階逐步積分方法相比有較高的計(jì)算精度。

(2) 本文的積分求微法一般是有條件穩(wěn)定算法，其穩(wěn)定性與截?cái)囗?xiàng)數(shù)L 有關(guān)。當(dāng)L=4 時(shí)，穩(wěn)定極限條件極易得到滿(mǎn)足，在較大的步距和較高的體系自振頻率參與下依然能得到穩(wěn)定的計(jì)算結(jié)果，因而具備較好的數(shù)值穩(wěn)定性。

(3) 本文方法可視為顯式方法，無(wú)需進(jìn)行方程求解，動(dòng)力響應(yīng)分析表現(xiàn)為矩陣與向量的乘法運(yùn)算，一次計(jì)算能得到多個(gè)時(shí)刻的響應(yīng)值，且無(wú)需對(duì)計(jì)算結(jié)果進(jìn)行插值，因而計(jì)算效率較高，計(jì)算工作量要少于傳統(tǒng)的顯式中心差分法。

(4) 建立本文方法的關(guān)鍵環(huán)節(jié)在于對(duì)DQ 原理的逆用，對(duì)于多自由度阻尼體系動(dòng)力響應(yīng)分析，計(jì)算動(dòng)力方程非齊次項(xiàng)時(shí)不再要求對(duì)定積分式進(jìn)行特殊處理，從而可以直接應(yīng)用精細(xì)積分法計(jì)算矩陣指數(shù)，有利于充分發(fā)揮精細(xì)積分法的優(yōu)勢(shì)。

(5) 本文方法屬于高階方法，在一個(gè)相對(duì)較長(zhǎng)的時(shí)段內(nèi)再劃分出多個(gè)時(shí)步，然后同時(shí)求解。這種模式為進(jìn)一步發(fā)展動(dòng)力響應(yīng)分析并行算法提供了有利條件。