基于微分求積的結構隨機振動時域分析方法

李博文，劉冬兵，王 永，張 磊，黎 慧，呂金偉

（1.三峽大學 電氣與新能源學院，湖北 宜昌443002；2.攀枝花學院 數學與計算機學院，四川 攀枝花617000；3.國網上海市電力公司 特高壓換流站分公司，上海201413）

實際工程中存在大量的隨機振動響應問題，隨機振動的分析方法一直都是科學技術人員關注的焦點。針對此問題，常用的分析方法有功率譜方法和時域分析法[1-3]，功率譜方法中比較有代表性的方法是林家浩提出的虛擬激勵法[4]，該方法已經在橋梁、水壩和車輛等眾多工程上獲得了廣泛的應用[5-8]。時域分析法在計算結構響應的均值和方差上更加簡單，本文主要介紹一種基于微分求積公式的時域分析法，在虛擬激勵法框架下，微分求積公式也可以用來分析隨機振動響應問題，這方面的內容將另外報道。

鐘萬勰提出的精細積分法[9]在結構動力響應的計算中具有精度高的特點，因此精細積分法也被用來計算隨機振動的動力響應。張森文等[10]提出了針對激勵源為平穩和非平穩白噪聲激勵時隨機振動響應計算的精細積分時域平均法，該方法采用精細積分和Simpson積分公式對隨機振動動力學方程進行時域積分，但是沒有給出任意激勵下的一般化計算公式。蘇成等[11]提出了隨機激勵線性化假設條件下隨機振動響應的精細積分顯式遞推求解公式，但是該算法存在計算精度的問題。隨后蘇成等[12]研究了非線性系統的非平穩響應，引入了等效激勵的概念，導出了非線性系統隨機振動響應的顯示迭代公式。李雪平等[13]采用等效線性化將非線性振動體系線性化后，直接用Newmark-β顯式積分格式進行迭代求解，該算法的最大相對誤差只有3.7%，計算效率相對蒙特卡羅法也比較高。陳玉震等[14]、蘇成等[15]和袁騰飛[16]都將Newmark-β法用于隨機振動的動力響應分析，并結合相應的技術手段驗證所提算法的有效性。項盼等[17]同時考慮了結構參數的不確定性和外激勵的隨機性，提出了一種混合方法，在對不確定的結構系統分析時更為可靠。劉祥等[18]提出了一種自適應降維的改進點估計法，與蒙特卡羅法相比，計算效率提高了2 至3個數量級。文獻[19]提出了一種將虛擬激勵和改進的精細積分方法結合的計算方法，極大提高了非平穩隨機振動響應的計算效率。文獻[20]中利用復化Cotes 積分公式計算隨機振動動力學方程時間離散后的Duhamel 積分項，計算精度較高，計算效率較蒙特卡羅法要高出很多。

微分求積法(Differential Quadrature Method,DQM)是由Bellman 等提出的一種無網格數值方法[21]，其在科學和工程計算中得到了廣泛的應用[22-24]。本文從時域分析法入手，利用微分求積公式，結合精細積分法和復化積分手段，提出了一種精細積分-復化微分求積法用于隨機振動的時域分析，并推導出了一套顯式的遞推計算公式。通過分析一個非平穩隨機激勵的算例，并與文獻[20]中的算法作對比，結果表明，在保證計算精度的前提條件下，本文提出的算法具有更高的計算效率。

1 任意激勵下隨機過程描述

實際工程中實測數據表明，大多數作用于工程結構的動力外荷載除了隨時間變化外，還具備動態特性和隨機特性。根據施加激勵隨時間變化的情況可以將其分為平穩隨機激勵和非平穩隨機激勵兩種類型[1]。在實際的工程分析中，任意的隨機激勵都可以通過聯合密度函數、功率譜密度函數和概率密度函數等來表征。

連續型隨機激勵記為X(t)、Y(t)。其中X(t)的一階和2階概率密度函數可以分別表示為pX(x,t)和pX(x1,t1;x2,t2)，根據數學期望和自相關系數的定義，可以求得隨機激勵X(t)的均值和自相關系數如下，Y(t)的相關量可同理求解。

隨機激勵X(t)和Y(t)的聯合概率密度函數記為ρXY(x,t1;y,t2)，從而定義二元隨機激勵X(t)、Y(t)的互相關系數如下：

文中所提的結構隨機振動時域分析法需要在時域上將時間區間離散為一系列的時間節點，導致隨機激勵也不再連續，成為在不同時間節點上的離散點，上述參數分別記為mX(ti)、RXX(ti,tj)、RXY(ti,tj) (i,j=0,1,…,N）。其中：N表示時間域上離散節點數，根據上述公式可以求得構成協方差矩陣Cov[X，Y]的相關函數。以上描述對于平穩隨機激勵和非平穩隨機激勵均適用。

2 復化微分求積法

微分求積法可以表示成等值Runge-Kutta(RK)方法的形式[25]。受此啟發，文獻[26]提出采用微分求積法來計算非齊次動力學方程的Duhamel 積分項。為了進一步提高微分求積法計算數值積分的精度，一個可行的途徑是提高微分求積法的階數，但是這會導致相應的計算量變大，而且數值積分穩定性也會變差。此外，另外一個常用的提高數值積分精度的方法是復化求積法，其中典型的代表是復化梯形公式、復化Simpson積分和復化Cotes積分。

記可積函數f(x)在積分區間[a，b]求定積分I(f)，采用復化微分求積法給出其一個近似值。首先，將積分區間等分為m份，每等份的區間長度即為積分步長h=(b-a)/m，數值積分節點xk=a+kh，k=0,1,…,m，復化微分求積公式為

式中：xk+i/s表示一個積分單元區間[xk,xk+1]上微分求積法的積分內節點，bi是相應節點處的積分系數。根據插值型積分公式的定義來看，復化微分求積公式也屬于這一類型，不難驗證采用均勻網格的s級微分求積法具有(s-1)階代數精度。

作為具體的計算公式，下面給出s=3時2階代數精度復化DQM公式的計算公式如下：

同理，s=4時，3階代數精度復化DQM公式的計算公式如下：

顯然，當m=1時，式(4)便是傳統的微分求積公式。

下文將采用精細積分-復化微分求積法(PI-complex differential quadrature formulas，PI-CDQs)求解在任意隨機激勵下結構隨機動力響應。

3 基于PI-CDQs的結構隨機響應分析

對于包含n個自由度的線性系統，在多個隨機激勵作用下的動力微分方程可表示為

式中：M、C和K分別表示結構的質量矩陣、阻尼矩陣和剛度矩陣，U表示結構的位移向量，表示結構的速度向量和加速度向量，F(t)表示多個隨機激勵構成的列向量。引入中間變量P∈Rn×1，令：

從而可將式(7)變換成下式：

式中：

在任意積分區間[tl,tl+1]內，非齊次動力方程式(9)在tl+1時刻可以表示為如下的Duhamel 積分方程，即

式中：等式右邊第1 項中指數矩陣T=eHΔt采用精細積分法計算，具體的計算公式詳見文獻[9]。等式右邊第2 項被稱作Duhamel 積分項，其近似計算誤差是式(10)的計算誤差的主要來源。

為此，采用3階代數精度復化DQM 計算式(10)中的Duhamel 積分項。將積分區間[tl,tl+1]均勻等分為m個子區間，相應的區間積分節點為tk=tl+kh,(k=0,1,…,m)，h=(tl+1-tl)/m=Δt/m。式(10)中的Duhamel積分項可以計算如下：

為了進一步地推導Z(tl+1)的顯式表達式，將式(11)改寫成矩陣的形式：

其中：

注意到，B的構成與微分求積法的級數s和復化積分的次數m有關。由于微分求積法和復化積分都采用均勻網格，因此，積分區間[tl,tl+1]的區間長度被等分為ms等份，每個子區間的積分步長為Δt/(ms)，式(14)中各個積分節點的指數矩陣可以統一地表示為

根據式(15)，只需通過精細積分法計算出T10=e(Δt/4m)H，就可以利用指數矩陣的加法運算非常方便形成式(13)中矩陣B內所有的分塊矩陣。

由此，式(10)可以進一步表示為

1.從經濟視角上看，區塊鏈是共享經濟的價值互聯網。例如，邵奇峰等學者指出，區塊鏈在互不了解的交易雙方間建立了可靠的信任，去中心化地實現了可信的價值傳輸，因此區塊鏈被稱為價值互聯網。

假設初始時刻Z(t0)=0，可得：

設

將式(17)進一步表示成下式：

式(19)即為PI-CDQs法計算結構隨機響應的顯式表達式。當然，文中只是以3階代數精度的DQM作為例子推導出了PI-CDQs法的一般計算公式，對于其他代數精度的DQM 也不難得出類似的計算公式。

根據隨機變量的1階原點距和2階中心距的計算公式，容易依據式(19)計算tl時刻結構隨機響應Z(tl)的均值向量和協方差矩陣分別如下：

式中，協方差矩陣Cov[Wl,Wl]可由隨機激勵的相關函數組成，相關系數可以由式(2)至式(3)求得。

4 算例分析

考慮一個含雙源隨機激勵的兩自由度振動結構，其動力學方程為

式中，初始條件為x1(0)=x2(0)=0,(0)=(0)=0。

雙源隨機激勵是簡諧隨機激勵源，其中f1(t)和f2(t)是互相獨立的，其解析表達式為

式中，隨機變量Θ在區間[0，2π]上服從均勻分布，其1階和2階概率密度函數分別為式(23)和式(24)。隨機變量a服從均值為10、方差為1的正態隨機分布，其1階和2階概率密度函數分別為式(25)和式(26)，聯合概率密度為式(27)。

表1 兩源隨機激勵下2自由度體系的動力響應均值

顯然，隨機相位正弦波f1(t)的均值為0，且自相關函數僅是時間間隔τ的函數，因此f1(t)是平穩隨機激勵。f2(t)是在平穩過程a的基礎上引入了一個隨時間變化的強度包絡函數，因此f2(t)為非平穩隨機激勵。

分別采用PI-CDQs(s=4)、蒙特卡羅法和文獻[20]的方法計算該兩自由度振動系統動力響應的均值和方差。時域積分區間為[0，5]s。計算結果如表1和表2所示。仿真平臺為MATLAB R2012a，計算機處理器為AMD Ryzen5 3500U。

由表1可知，PI-CDQs(s=4，m=1)在步長Δt=0.1 s時，其計算精度和文獻[20]中的計算精度基本保持一致，最大偏差也僅有1.11%，而且PI-CDQs的計算效率較文獻[20]的計算效率高。當PI-CDQs的計算步長變為Δt=0.5 s時，隨著復化次數m的增大，其計算的偏差也隨之減小，但計算時間以毫秒級增大；而且相比于小步長的結果，計算精度幾乎沒有變化，這說明復化次數對PI-CDQs的計算效率影響不大，且可以采用大步長進行復化積分。值得注意的是文獻[20]所采用的復化Cotes 積分公式具有4階代數精度，而本文所提PI-CDQs僅有3階代數精度。

從表2可以看出，對于結構位移響應的方差，PICDQs(s=4，m=1)在步長Δt=0.1s時計算偏差和文獻[20]中的計算偏差基本保持一致，最大偏差也僅有0.317%，而且PI-CDQs的計算效率也較文獻[20]的計算效率要高。

圖1至圖2給出了采用蒙特卡羅法(Δt=0.01 s)和PI-CDQs所模擬該隨機振動體系的位移變化量x1、x2的均值響應。

圖1 結構x1位移響應的均值

從中可以看出蒙特卡羅法和PI-CDQs的仿真結果很好重合，這也說明本文算法的精度高。圖3至圖4給出了蒙特卡羅法(Δt=0.01 s)和PI-CDQs所模擬該隨機振動體系的位移變化量x1、x2的方差響應曲線。從中可以看出蒙特卡羅法和PI-CDQs的仿真結果也較好重合，進一步地說明本文算法的精度高。

圖2 結構x2位移響應的均值

圖3 結構x1位移響應的方差

表2 兩源隨機激勵下兩自由度體系的動力響應方差

5 結語

文中將復化微分求積法用于計算隨機振動體系的動力學方程中的Duhamel 積分項，并與精細積分法結合，提出了計算任意隨機激勵下結構振動響應的顯式計算公式。復化微分求積法的典型特點在于積分區間的首末端點處的被積函數值不參與積分過程，這是其與Newton-Cotes公式的不同之處。

圖4 結構x2位移響應的方差

通過與采用Cotes 積分公式計算隨機振動動力響應的算法對比，分析了在平穩隨機激勵和非平穩隨機激勵下，采用本文算法和文獻[20]中方法在都采用具有遞推關系的顯式表達式條件下所計算的隨機振動響應，仿真結果表明，本文算法的計算精度在與文獻[20]中算法計算精度相當的前提下，具有更高的計算效率。