非平穩(wěn)隨機響應(yīng)靈敏度分析的時域顯式法

陳太聰，蘇 成，胡智強，馬海濤

(華南理工大學(xué)土木與交通學(xué)院，亞熱帶建筑科學(xué)國家重點實驗室，廣東 廣州510640)

非平穩(wěn)隨機響應(yīng)靈敏度分析的時域顯式法

陳太聰，蘇 成，胡智強，馬海濤

(華南理工大學(xué)土木與交通學(xué)院，亞熱帶建筑科學(xué)國家重點實驗室，廣東 廣州510640)

針對非平穩(wěn)激勵下的結(jié)構(gòu)振動問題，研究動力響應(yīng)方差靈敏度的高效時域求解算法。首先推導(dǎo)確定性動力響應(yīng)靈敏度的時域顯式表達，繼而結(jié)合方差靈敏度的一般計算公式，得到非平穩(wěn)隨機響應(yīng)方差靈敏度的時域顯式計算列式。該列式同樣適用于平穩(wěn)激勵下結(jié)構(gòu)瞬態(tài)響應(yīng)方差靈敏度的求解。以框架結(jié)構(gòu)和桁架結(jié)構(gòu)分別受不同類型的非平穩(wěn)激勵為例，時域顯示解法和其他方法的對比計算結(jié)果驗證了時域顯示解法的有效性。

隨機振動；非平穩(wěn)；靈敏度；時域；顯式法

引 言

靈敏度分析可以定量地確定系統(tǒng)參數(shù)的改變對系統(tǒng)響應(yīng)的影響［1］，在結(jié)構(gòu)優(yōu)化、最優(yōu)控制和系統(tǒng)辨識等領(lǐng)域有著廣泛的應(yīng)用。由于工程結(jié)構(gòu)面對的大量激勵如地震、風(fēng)、海浪等作用屬于隨機過程激勵，因此，開展針對隨機振動響應(yīng)的靈敏度分析具有重要的現(xiàn)實意義，相關(guān)研究日益受到關(guān)注。

根據(jù)隨機過程激勵的平穩(wěn)性質(zhì)不同，相應(yīng)響應(yīng)也可分為平穩(wěn)隨機響應(yīng)和非平穩(wěn)隨機響應(yīng)兩類。其中，由于平穩(wěn)隨機響應(yīng)較易求解［2］，相關(guān)靈敏度分析的研究較為成熟，已發(fā)展了包括有色噪聲激勵下的代數(shù)Riccati方程解法［3］、相關(guān)/非相關(guān)高斯激勵下的模態(tài)分解解法［4］、高斯/非高斯激勵下的響應(yīng)多階統(tǒng)計矩靈敏度計算的時域解法［5］、高斯激勵下的虛擬激勵解法［6］，以及隨機結(jié)構(gòu)情況下的Neumann展開結(jié)合Monte Carlo模擬解法［7］和虛擬激勵結(jié)合點估計解法［8］等多種分析方法。

而對于第2類靈敏度問題，雖然地震、風(fēng)、海浪等隨機過程激勵本質(zhì)上都是非平穩(wěn)的，相關(guān)分析更具現(xiàn)實意義，但由于非平穩(wěn)隨機響應(yīng)的求解本身較為困難，因此相應(yīng)的靈敏度分析還不易進行。近年來，Chaudhuri和Chakraborty［9］在頻域內(nèi)給出響應(yīng)功率譜和各階譜矩的靈敏度計算方法，進而得到響應(yīng)方差和可靠度的靈敏度，其中需要進行雙重頻域積分；Cacciola等［10］結(jié)合動態(tài)模型修正和模態(tài)分解方法，建立了響應(yīng)方差靈敏度的微分方程，采用遞推法和逐步積分策略進行求解；Marano等［11］建立Lyapunov微分方程，求解響應(yīng)方差的靈敏度，但仍限于單自由度系統(tǒng)的應(yīng)用；徐文濤等［12］和劉齊茂［13］都從虛擬激勵法出發(fā)，分別采用精細積分法和Newmark-β法推導(dǎo)了響應(yīng)功率譜的一階和二階靈敏度的計算列式，前文繼而給出了響應(yīng)方差靈敏度的計算列式，需要說明的是，虛擬激勵法的應(yīng)用需要進行時域和頻域內(nèi)的混合積分計算。

近年來，蘇成和徐瑞［14］提出了非平穩(wěn)隨機振動分析的時域顯式法，通過建立的結(jié)構(gòu)動力響應(yīng)時域顯式表達，可在時域內(nèi)直接計算隨機響應(yīng)的均值和方差，還可以結(jié)合Monte Carlo模擬求解構(gòu)件動力可靠度和體系動力可靠度［15］，在大型復(fù)雜結(jié)構(gòu)中的應(yīng)用顯示了良好的計算效率［16］，在隨機結(jié)構(gòu)動力學(xué)問題［17］和 非線性隨 機振動 問題［18］中也得 到 了應(yīng)用，為非平穩(wěn)隨機響應(yīng)靈敏度問題的研究打下了良好基礎(chǔ)。本文將以該方法為基礎(chǔ)，推導(dǎo)結(jié)構(gòu)動力響應(yīng)靈敏度的時域顯式表達式，并以此進一步提出非平穩(wěn)激勵下隨機響應(yīng)靈敏度分析的時域顯式法。最后通過數(shù)值算例說明本文方法的準確性和可行性。

1 動力響應(yīng)分析的時域顯式法

隨機振動分析的時域顯式法是基于結(jié)構(gòu)動力響應(yīng)的顯式表達推導(dǎo)得到的［14］。以下就把該顯式表達作一簡要介紹。n自由度結(jié)構(gòu)系統(tǒng)的運動方程可表示為

式中 K，C和M分別代表結(jié)構(gòu)系統(tǒng)的剛度矩陣、阻尼矩陣和質(zhì)量矩陣；Y，.Y和‥Y分別為位移向量、速度向量和加速度向量；L為一n×m階激勵影響矩陣；F為m維激勵向量，若結(jié)構(gòu)承受地震作用，則F為地面加速度。

式(1)可表示成狀態(tài)方程的形式

若考慮線性結(jié)構(gòu)系統(tǒng)，以及等時間步長Δt的計算，則一般地，第i時刻(ti=iΔt)的系統(tǒng)狀態(tài)可遞推地表示為

式中 矩陣Q0，Q1和Q2都依賴于結(jié)構(gòu)參數(shù)和時間步長Δt。

若初始系統(tǒng)狀態(tài)V0=0，則由遞推式(4)可以推導(dǎo)得到第i時刻系統(tǒng)狀態(tài)的時域顯式表達為

則，顯式表達式(5)可重新表示為

因此，所有時刻的系統(tǒng)狀態(tài)可計算如下

值得注意的是，系數(shù)Ai，j(0≤j≤i)的計算是一個遞推過程，若采用精細積分法［19］進行求解，可得

由式(8)和(9)可見，為了得到各時刻系統(tǒng)狀態(tài)對應(yīng)的完整系數(shù)矩陣Ai，j，只需要計算系數(shù)矩陣的前兩列Ai，0和Ai，1(i=1，2，…，l)即可，其計算量相當(dāng)于兩次確定性動力時程分析［14］。由式(6)可知，Bi也可以通過其前兩列系數(shù)完全確定。

假設(shè)第i時刻關(guān)注的某個結(jié)構(gòu)響應(yīng)為vi，不失一般性，vi可由系統(tǒng)狀態(tài)V和激勵F計算得到

式中 φ為關(guān)注結(jié)構(gòu)響應(yīng)的定位行向量，其元素由0和1組成；S1和S2分別為系統(tǒng)狀態(tài)和系統(tǒng)激勵對結(jié)構(gòu)響應(yīng)的影響矩陣。根據(jù)所關(guān)注的結(jié)構(gòu)響應(yīng)的性質(zhì)不同，S1和S2的具體取值也有所不同：①若關(guān)注結(jié)構(gòu)響應(yīng)為結(jié)點位移或速度，則S1=I，S2=0；②若為結(jié)點加速度，則S1=H，S2=W；③若為單元應(yīng)力或應(yīng)變，則在結(jié)點激勵的情況下，S1可通過單元應(yīng)力矩陣或應(yīng)變矩陣計算得到，S2=0。

由式(7)和(13)，最終可得

2 動力響應(yīng)靈敏度分析的時域顯式表達

基于以上動力響應(yīng)時域顯式表達的基本思路，以下推導(dǎo)動力響應(yīng)靈敏度的時域顯式求解列式。

設(shè)θ代表結(jié)構(gòu)的某個設(shè)計參數(shù)，在狀態(tài)方程式(2)兩端對θ求偏導(dǎo)，得

若激勵F與參數(shù)θ無關(guān)，則式(15)可化簡為

對比式(6)和(18)，再結(jié)合式(8)和(9)可知，要計算所有時刻點的系統(tǒng)狀態(tài)靈敏度，只需要計算系數(shù)矩陣和的前兩列，即可完全確定。

進一步地，結(jié)合式(13)和(17)可知，第i時刻關(guān)注結(jié)構(gòu)響應(yīng)vi對設(shè)計參數(shù)θ的靈敏度即可通過下式計算

3 隨機動力響應(yīng)靈敏度分析的時域顯式法

由式(14)可知，當(dāng)結(jié)構(gòu)激勵F為隨機激勵時，第i時刻結(jié)構(gòu)響應(yīng)vi的方差可表達為

將式(14)和(20)代入，最終整理可得第i時刻結(jié)構(gòu)響應(yīng)vi的方差對設(shè)計參數(shù)θ的靈敏度的計算表達式為

在實際工程結(jié)構(gòu)中，較多關(guān)注結(jié)點位移或速度響應(yīng)，此時，S1=I以及S2=0，則式(21)和(25)可分別化簡為

在實際計算過程中，cov(Ri，Ri)的計算量和存儲量可能過大，導(dǎo)致計算效率降低。為此，將式(31)右端的第二項展開，并整理后，得到以下更方便使用的計算表達式

需要說明的是，以上推導(dǎo)中并未對隨機激勵F的分布特性預(yù)設(shè)任何前提，因此，所得均方響應(yīng)及其靈敏度的計算列式適用于各種類型的平穩(wěn)/非平穩(wěn)隨機過程激勵情況。對于平穩(wěn)隨機過程情形，通過本文列式獲得的是全過程瞬態(tài)響應(yīng)的方差及其靈敏度，相對于常規(guī)平穩(wěn)隨機振動分析得到的穩(wěn)態(tài)階段解答，信息更為豐富。

4 數(shù)值算例

采用FORTRAN語言實現(xiàn)了非平穩(wěn)隨機響應(yīng)靈敏度分析的時域顯式解法，并進行驗證。本節(jié)給出兩種結(jié)構(gòu)的數(shù)值算例，分別與直接差分法、差分Monte Carlo模擬法和已有文獻結(jié)果進行比較，說明本文方法的可靠性。

4.1 平面框架結(jié)構(gòu)

采用圖1所示平面4層框架結(jié)構(gòu)模型，含20個梁單元，共36個自由度。所有單元的彈性模量E=3.0×1010N/m2，質(zhì)量密度ρ=2.4×103kg/m3，梁、柱截面分別為0.25 m×0.40 m和0.35 m×0.35 m，考慮Rayleigh阻尼模型，阻尼矩陣由前兩階模態(tài)阻尼比(ξ=0.05)確定。

圖1 平面框架結(jié)構(gòu)模型Fig.1 Aplanar frame structure model

考慮非平穩(wěn)限帶白噪聲激勵，F(xiàn)(t)=g(t)x(t)，其中，均勻調(diào)制函數(shù)g(t)取為

其作用頻帶范圍ω=0～200 rad/s。則該隨機激勵F(t)的相關(guān)函數(shù)可由下式計算得到

RFF(t+τ，τ)=g(t+τ)g(τ)Rxx(τ)

設(shè)計參數(shù)取為所有桿件的彈性模量E，計算總時長取為15 s，計算步長取Δt=0.005 s。

采用本文方法，計算結(jié)點 A的水平方向位移、速度和加速度的方差及其對設(shè)計參數(shù)的靈敏度隨時間的變化規(guī)律，結(jié)果如圖2～4所示。作為對比，圖中還給出了采用直接差分法和差分Monte Carlo模擬法的計算結(jié)果。其中，直接差分法的計算過程是取設(shè)計參數(shù)的1‰變化量為差分步長，按式(21)分別計算設(shè)計參數(shù)變化前后的響應(yīng)方差時程，進而用差分法計算靈敏度時程；差分Monte Carlo模擬法同樣取設(shè)計參數(shù)的1‰變化量為差分步長，根據(jù)隨機激勵模型，人工生成一定數(shù)目的激勵過程樣本并進行動力時程分析，最后統(tǒng)計得到設(shè)計參數(shù)變化前后的響應(yīng)方差時程，進而用差分法計算靈敏度時程。

圖2 A點位移的方差及其靈敏度Fig.2 Covariance and its sensitivity of the node Adisplacement

實際計算中發(fā)現(xiàn)，由于Monte Carlo模擬得到的響應(yīng)方差存在偏差，通過差分Monte Carlo模擬法計算得到的靈敏度值會發(fā)生偏差波動，但隨著樣本數(shù)增多，波動幅度會逐漸減小。在本例中，當(dāng)樣本數(shù)取為105時，靈敏度值波動幅度較小，因此，圖中給出的Monte Carlo法結(jié)果均為105樣本數(shù)的計算結(jié)果。

圖3 A點速度的方差及其靈敏度Fig.3 Covariance and its sensitivity of the node Avelocity

圖4 A點加速度的方差及其靈敏度Fig.4 Covariance and its sensitivity of the node Aacceleration

由圖示結(jié)果可以觀察得到：

(1)本文方法和直接差分法的結(jié)果基本一致，并都與差分Monte Carlo模擬法的結(jié)果趨勢基本一致，驗證了本文方法的正確性；

(2)方差的Monte Carlo模擬結(jié)果與理論計算結(jié)果的符合程度高于方差靈敏度，這是因為在差分步長變化前后，模擬得到的方差結(jié)果都存在一定的偏差，因此由差分計算得到的方差靈敏度的偏差會大于方差的偏差。由文獻［1］的分析結(jié)果可知，響應(yīng)靈敏度的差分模擬偏差會比響應(yīng)的模擬偏差大1～2個量級，本算例與之相符。如前所述，計算中已驗證，隨著Monte Carlo樣本數(shù)增加，方差的模擬偏差變小，方差靈敏度的模擬偏差也隨之減小。

4.2 平面桁架結(jié)構(gòu)

本算例取自文獻［10］，采用的輸電塔結(jié)構(gòu)如圖5所示。該結(jié)構(gòu)有限元模型采用24個平面桁架單元，共20個自由度，所有單元具有相同的軸向剛度，EA=1.210 4×108N?？紤]輸電線的自重，結(jié)點9和12具有集中質(zhì)量m=1 200 kg，其余結(jié)點具有集中質(zhì)量m=600 kg?？紤]Rayleigh阻尼模型，阻尼矩陣由前兩階模態(tài)阻尼比(ξ=0.02)確定。

考慮非平穩(wěn)地震作用，采用均勻調(diào)制的平穩(wěn)隨機過程來模擬非平穩(wěn)地面水平運動加速度‥x(t)，其中平穩(wěn)隨機過程采用以下Tajimi-Kanai功率譜

S0=0.05 m2/s3，ωg=4πrad/s，ζg=0.6均勻調(diào)制函數(shù)g(t)取為

則該隨機激勵的相關(guān)函數(shù)可由下式計算得到：

圖5 平面桁架結(jié)構(gòu)模型Fig.5 Aplanar truss structure model

設(shè)計參數(shù)取為所有桿件的軸向剛度EA，計算總時長取為15 s，計算步長取Δt=0.01 s。

圖6 結(jié)點12位移的方差靈敏度時程Fig.6 Sensitivity of covariance of the node 12 displacement

采用本文方法，計算結(jié)點12的水平方向位移的方差對設(shè)計參數(shù)的靈敏度，結(jié)果如圖6所示。由圖示結(jié)果可見，本文方法與文獻［10］中的靈敏度分析結(jié)果基本一致，進一步驗證了本文方法的正確性。

5 結(jié) 論

本文基于動力響應(yīng)分析的時域顯式法，提出了非平穩(wěn)隨機響應(yīng)靈敏度分析的時域顯式法。利用靈敏度方程與狀態(tài)方程的相似性，推出動力響應(yīng)靈敏度的時域顯式表達。繼而利用時程方差靈敏度計算的一般公式，最終得到結(jié)構(gòu)響應(yīng)方差靈敏度分析的時域顯式表達式。不同類型的算例計算結(jié)果驗證了本文方法的正確性。本文的時域顯式方法不要求隨機過程激勵具有特殊分布特性，因此，可以廣泛適用于各種類型激勵情況下的平穩(wěn)/非平穩(wěn)隨機響應(yīng)靈敏度分析。

［1］陳太聰，韓大建，蘇成.參數(shù)靈敏度分析的神經(jīng)網(wǎng)絡(luò)方法及其工程應(yīng)用［J］.計算力學(xué)學(xué)報，2004，21(6)：752—756.Chen Taicong，Han Dajian，Su Cheng.Neural network method in parameter sensitivity analysis and its application in engineering［J］.Chinese Journal of Computational Mechanics，2004，21(6)：752—756.

［2］Soong T T.Random Vibration of Mechanical and Structural Systems［M］.Englewood Cliffs：PTR Prentice Hall，1993：172—198.

［3］姚昌仁，麻永平.結(jié)構(gòu)隨機激勵的響應(yīng)靈敏度分析［J］.力學(xué)學(xué)報，1990，22(1)：438—445.Yao Changren，Ma Yongping.The response sensitivity analysis for structural systems in random excitation［J］.Acta Mechanica Sinica，1990，22(1)：438—445.

［4］Tong W H，Jiang J S，Gu S N.Dynamic design of structures under random excitation［J］.Computational Mechanics，1998，22(5)：388—394.

［5］Benfratello S，Caddemi S，Muscolino G.Gaussian and non-Gaussian stochastic sensitivity analysis of discrete structural system［J］.Computers and Structures，2000，78：425—434.

［6］徐文濤，張亞輝，林家浩.基于虛擬激勵法的車輛振動靈敏度分析及優(yōu)化［J］.機械強度，2010，32(3)：347—352.Wu Wentao，Zhang Yahui，Lin Jiahao.PEMbased sensitivity analysis for vehicle ride comfort and optimization［J］.Journal of Mechanical Strength，2010，32(3)：347—352.

［7］Bhattacharyya B，Chakraborty S.Stochastic dynamic sensitivity of uncertain structures subjected to random earthquake loading［J］.Journal of Sound and Vibration，2002，249(3)：543—556.

［8］喬紅威，呂震宙.隨機結(jié)構(gòu)隨機激勵下的響應(yīng)靈敏度分析［J］.振動與沖擊，2008，27(3)：60—62.Qiao Hongwei，Lu Zhengzhou.Response sensitivity analysis of stochastic structures under non-stationary random excitation［J］.Journal of Vibration andShock，2008，27(3)：60—62.

［9］Chaudhuri A，Chakraborty S.Sensitivity evaluation in seismic reliability analysis of structures［J］.Computer Methods in Applied Mechanics and Engineering，2004，193(1/2)：59—68.

［10］Cacciola P，Colajanni P，Muscolino G.Amodal approach for the evaluation of the response sensitivity of structural systems subjected to non-stationary random Processes［J］.Computer Methods in Applied Mechanics and Engineering，2005，194：4 344—4 361.

［11］Marano G C，Trentadue F，Morrone E，et al.Sensitivity analysis of optimum stochastic nonstationary response spectra under uncertain soil parameters［J］.Soil Dynamics and Earthquake Engineering，2008，28：1 078—1 093.

［12］Xu W T，Zhang Y H，Lin J H，et al.Sensitivity analysis and optimization of vehicle-bridge systems based on combined PEM-PIMstrategy［J］.Computers and Structures，2011，89(3/4)：339—345.

［13］Liu Q M.Sensitivity and Hessian matrix analysis of power spectral density functions for uniformly modulated evolutionary random seismic responses［J］.Finite Elements in Analysis and Design，2012，48：1 370—1 375.

［14］蘇成，徐瑞.非平穩(wěn)激勵下結(jié)構(gòu)隨機振動時域分析法［J］.工程力學(xué)，2010，27(12)：77—83.Su Cheng，Xu Rui.Random vibration analysis of structures subjected to non-stationary excitations by time domain method［J］.Engineering Mechanics，2010，27(12)：77—83.

［15］蘇成，徐瑞.非平穩(wěn)隨機激勵下結(jié)構(gòu)體系動力可靠度時域解法［J］.力學(xué)學(xué)報，2010，42(3)：512—520.Su Cheng，Xu Rui.Time-domain method for dynamic reliability of structural systems subjected to non-stationary random excitations［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(3)：512—520.

［16］蘇成，徐瑞，劉小璐，等.大跨度空間結(jié)構(gòu)抗震分析的非平穩(wěn)隨機振動時域顯式法［J］.建筑結(jié)構(gòu)學(xué)報，2011，32(11)：169—176.SU Cheng，XU Rui，LIU Xiaolu，et al.Non-stationary seismic analysis of large-span spatial structures by time-domain explicit method［J］.Journal of Building Structures，2011，32(11)：169—176.

［17］Su C，Xu R.Random vibration analysis of structures by a time-domain explicit formulation method［J］.Structural Engineering and Mechanics，2014，52(2).

［18］Su C，Huang H，Ma HT，et al.Efficient MCS for random vibration of hysteretic systems by an explicit iteration approach［J］.Earthquakes and Structures，2014，7(2)：119—139

［19］鐘萬勰.結(jié)構(gòu)動力方程的精細時程積分法［J］.大連理工大學(xué)學(xué)報，1994，34(2)：131—136.Zhong Wanxie.On precise time-integration method for structural dynamics［J］.Journal of Dalian University of Technology，1994，34(2)：131—136.

An explicit time-domain method in sensitivity analysis of non-stationary random responses

CHEN Tai-cong，SU Cheng，HU Zhi-qiang，MAHai-tao
(State Key Laboratory of Subtropical Building Science，School of Civil Engineering and Transportation，South China University of Technology，Guangzhou 510640，China)

Aiming at the structural vibration problem under non-stationary excitation，time-domain method of high efficiency is explored in the present study to determine the sensitivity of covariance of random response.Firstly，a time-domain formulation is derived for computing the sensitivity of deterministic response.Then，according to a general expression of the sensitivity of covariance，an explicit time-domain formulation is deducted to calculate the sensitivity of covariance of non-stationary random response.This formulation is also suitable for the case of stationary excitation if sensitivity of covariance of the transient response is concerned.By taking a frame and a truss subjected to different types of non-stationary excitations as examples，comparisons of the numerical results obtained with different methods illustrate the effectiveness of the proposed method.

random vibration；non-stationary；sensitivity；time-domain；explicit method

O324；TU311.3

A

：1004-4523(2015)04-0503-07

10.16385/j.cnki.issn.1004-4523.2015.04.001

陳太聰(1977—)，男，副教授。電話：13903019936；E-mail：cvchentc@scut.edu.cn

蘇成(1968—)，男，教授。電話：(020)87112755；E-mail：cvchsu@scut.edu.cn

2014-03-21；修訂日期：2014-06-22

國家自然科學(xué)基金資助項目(51078150)；亞熱帶建筑科學(xué)國家重點實驗室自主研究項目(2013ZA01，2015ZC19)；中央高?；究蒲袠I(yè)務(wù)費資助項目(2013ZZ0024)和廣西科學(xué)研究與技術(shù)開發(fā)計劃資助項目(1298011-1)