李鴻晶，楊筱朋，梅雨辰

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

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

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

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

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

動態(tài)荷載向量p(t)激勵下的多自由度阻尼體系運動方程寫為：

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

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

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

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

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

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

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

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

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

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

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

2.2 定積分 Dk 的計算

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

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

3 算法步驟

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

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

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

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

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

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

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

5 算例分析

通過兩個算例闡釋本文方法的特性，第一個算例用來驗證本文方法的可靠性，包括收斂性和精度；第二個算例重點闡釋本文方法的適用性。

5.1 兩自由度體系算例

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

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

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

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

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

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

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

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

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

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

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

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

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

續(xù)表 4

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

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

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

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

圖 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)定性。取截斷項數(shù)L=4、時段分段數(shù)ρ=10，采用本文IDM法計算體系每隔0.01 s 的動力響應(yīng)。選用不同的時段長度Δt，計算結(jié)果與顯式中心差分法(取Δt=0.001 s)進(jìn)行比較。圖5 和圖6 分別為最右端質(zhì)點(DOF 100)和第50 個質(zhì)點(DOF 50)的位移時程曲線。這里需要說明的是，對于該體系中心差分法的穩(wěn)定極限為Δt≤0.0018 s，因而采用Δt=0.001 s可得到穩(wěn)定的計算結(jié)果，且它是直接獲得0.01 s 時間間隔響應(yīng)所允許的最大步距。

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

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

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

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

表 6 不同方法CPU 時間對比Table 6 Comparison of CPU time with different methods

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

6 討論與結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

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

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