李鴻晶,楊筱朋,梅雨辰
(南京工業(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ù)手段。
動態(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é)果。
利用式(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
根據(jù)矩陣指數(shù)的Taylor 級數(shù)定義,矩陣指數(shù)Tk可表示為:
由式(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)于是一種顯式算法。
依據(jù)上述推導(dǎo),可將多自由度阻尼體系動力時程積分求微法總結(jié)為下述算法:
對于線性體系動力方法的數(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)定性的限制。
通過兩個算例闡釋本文方法的特性,第一個算例用來驗證本文方法的可靠性,包括收斂性和精度;第二個算例重點闡釋本文方法的適用性。
選擇文獻(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é)果仍然比其要精確些,這直觀顯示了本文方法高精度的特點。
以圖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)勢。
根據(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)分析并行算法提供了有利條件。