李雪平 李棟泓 蘇成 魏鵬

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

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

李雪平 李棟泓 蘇成 魏鵬?

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

分別采用精細(xì)積分格式和Newmark-β積分格式推導(dǎo)了非平穩(wěn)隨機(jī)地震響應(yīng)靈敏度的時域顯式表達(dá)式.并分別將其應(yīng)用于平面框架與平面桁架問題的靈敏度分析,研究了積分時間步長對兩種積分格式的計算精度和效率的影響.發(fā)現(xiàn):當(dāng)結(jié)構(gòu)主頻處于地震主頻范圍內(nèi)時,在同等精度條件下,基于精細(xì)積分的時域表達(dá)式效率更高;而在結(jié)構(gòu)主頻偏離地震荷載主頻范圍時,基于Newmark-β積分格式的時域表達(dá)式效率更高.該研究成果可以為在考慮非平穩(wěn)隨機(jī)振動的結(jié)構(gòu)優(yōu)化問題中的數(shù)值積分算法選擇提供有效的參考和借鑒.

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

靈敏度分析即求導(dǎo)信息,在結(jié)構(gòu)優(yōu)化、最優(yōu)控制和系統(tǒng)辨別等領(lǐng)域都有廣泛的應(yīng)用.它是評價結(jié)構(gòu)參數(shù)的改變引起結(jié)構(gòu)響應(yīng)特性變化率的方法.在隨機(jī)響應(yīng)靈敏度分析中,國內(nèi)外已有相關(guān)學(xué)者對此進(jìn)行了研究.文獻(xiàn)[1]對結(jié)構(gòu)受隨機(jī)激勵的響應(yīng)靈敏度進(jìn)行分析,建立了有色噪聲下響應(yīng)協(xié)方差并推導(dǎo)了相應(yīng)靈敏度表達(dá)式.文獻(xiàn)[2]采用虛擬激勵法推導(dǎo)平穩(wěn)激勵下的響應(yīng)及靈敏度,并應(yīng)用于求解帶有不確定參數(shù)結(jié)構(gòu)在平穩(wěn)激勵下的響應(yīng).文獻(xiàn)[3]采用虛擬激勵法,多次利用狀態(tài)空間解耦及矩陣微分理論等,推導(dǎo)出多源隨機(jī)激勵下響應(yīng)均方根的參數(shù)靈敏度表達(dá)式,并應(yīng)用于汽車平順性分析.文獻(xiàn)[4]基于虛擬激勵法,分別采用Newmark-β法推導(dǎo)了響應(yīng)功率譜的一階和二階靈敏度的計算列式.文獻(xiàn)[5]對于隨機(jī)激勵下隨機(jī)結(jié)構(gòu)動力響應(yīng)的靈敏度問題進(jìn)行分析,采用基于點(diǎn)估計法求解結(jié)構(gòu)響應(yīng)均方值的均值對結(jié)構(gòu)參數(shù)的敏度信息.文獻(xiàn)[6]在虛擬激勵法的框架下,采用泰勒展開表達(dá)式求解隨機(jī)結(jié)構(gòu)受平穩(wěn)激勵的響應(yīng),并進(jìn)行靈敏度分析.

近年來,由蘇成等[7- 8]提出的時域顯式法可顯式表達(dá)出結(jié)構(gòu)的動力響應(yīng)表達(dá)式,并可高效求解結(jié)構(gòu)在非平穩(wěn)激勵下的動力響應(yīng)及動力可靠度.陳太聰?shù)萚9]基于精細(xì)積分格式提出了非平穩(wěn)隨機(jī)響應(yīng)靈敏度的時域顯式法.無論是動力響應(yīng)還是靈敏度分析,其顯式解均可表示為與結(jié)構(gòu)參數(shù)有關(guān)的系數(shù)矩陣和與荷載有關(guān)的荷載向量乘積的表達(dá)式,而其中與結(jié)構(gòu)參數(shù)有關(guān)的系數(shù)矩陣可采用精細(xì)積分法[10]、Newmark-β法[11]等積分方法來推導(dǎo)求得.而這些數(shù)值積分方法的計算效率和精度均與積分時間步長有關(guān),但目前并未有文獻(xiàn)對此進(jìn)行細(xì)致分析.結(jié)構(gòu)動力靈敏度分析是結(jié)構(gòu)動力拓?fù)鋬?yōu)化過程中所需要求解的重要信息.由于結(jié)構(gòu)優(yōu)化是一個動態(tài)逐步迭代的過程，每一步都需要求解目標(biāo)函數(shù)和約束函數(shù)對設(shè)計參數(shù)的靈敏度,這樣靈敏度的計算效率變得非常重要.因此,有必要對不同積分格式下的時域顯式靈敏度分析方法進(jìn)行研究，以便找到更高效的算法.文中首先基于時域顯式法的思想，推導(dǎo)了兩種積分格式下的非平穩(wěn)隨機(jī)響應(yīng)靈敏度的時域顯式理論公式,接著應(yīng)用上述方法分別對一個平面框架和一個平面桁架在非平穩(wěn)隨機(jī)地震激勵下的靈敏度進(jìn)行了分析，并對計算精度和效率進(jìn)行對比.

1 隨機(jī)地震響應(yīng)靈敏度分析的時域顯式表達(dá)

1.1 基于精細(xì)積分法的時域顯式表達(dá)式

考慮某n個自由度的線性結(jié)構(gòu)系統(tǒng),在隨機(jī)地震作用下的動力學(xué)方程可以寫為

(1)

(1)

式中,

(3)

I為單位矩陣,當(dāng)系統(tǒng)的初始狀態(tài)V0=0時,由文獻(xiàn)[7]可推導(dǎo)出第i時刻的響應(yīng)Vi的時域顯式表達(dá)式為

(4)

(i=1,2…,l)

式中,

(5)

其中,

(6)

式(4)中Ai只與結(jié)構(gòu)參數(shù)有關(guān),反映結(jié)構(gòu)參數(shù)對動力響應(yīng)的影響.則響應(yīng)Vi可進(jìn)一步表達(dá)為

(7)

式中,

(8)

(9)

(10)

式中,

(11)

1.2 基于Newmark-β法的時域顯式表達(dá)式

采用Newmark-β數(shù)值積分格式推導(dǎo)出動力學(xué)方程(1)在第i時刻的響應(yīng)Vi的顯式遞推關(guān)系式[12- 13]:

(12)

(i=1,2,…,l)

式中,

H21=a4(H11-I)+a6M-1K,

H22=a4H12-a5I+a6M-1C,

R3=a4R1-a6M-1,R4=a4R2,

S1=a1M+a4C,

S2=a2M+a5C,

S3=a3M+a6C,

積分參數(shù)γ和β可分別取0.5和0.25,以保證數(shù)值積分的穩(wěn)定性.注意此時遞推矩陣T2與式(6)中T1不同.當(dāng)系統(tǒng)的初始狀態(tài)V0為零向量,并且記N=T2Q2+Q1時,由遞推關(guān)系式可以推導(dǎo)出響應(yīng)Vi的時域顯式表達(dá)式:

(13)

式中,

(14)

其中,Ai也是只與結(jié)構(gòu)參數(shù)有關(guān),響應(yīng)Vi也可進(jìn)一步表達(dá)為與式(7)相同的形式,但其中的系數(shù)Bi表示如下:

(15)

由于在結(jié)構(gòu)振動分析中一般只對結(jié)構(gòu)某些關(guān)鍵部位響應(yīng)感興趣,通常并不需要求解所有的響應(yīng)向量,從而可以節(jié)省大量的計算工作量.由式(7)可知,在時刻i某一關(guān)鍵部位的響應(yīng)vi及方差可寫為

(16)

(17)

其中,φ為關(guān)注結(jié)構(gòu)響應(yīng)的定位行向量,其元素由0和1組成.

同樣,基于此思路也可推導(dǎo)動力響應(yīng)靈敏度的時域顯式求解列式.設(shè)θ代表結(jié)構(gòu)的某設(shè)計參數(shù),則對動力學(xué)方程式(1)兩端求偏導(dǎo),并考慮荷載與結(jié)構(gòu)參數(shù)無關(guān),則可得:

(18)

式中,

(19)

對比式(1)和式(19)可見,兩方程僅荷載項(xiàng)不同，因此也可以采用Newmark-β法推導(dǎo)其靈敏度的時域表達(dá)式,形式與式(10)相同,其中的系數(shù)表示如下:

(20)

由式(7)和(10)可以看出,無論是響應(yīng)還是靈敏度都可通過與結(jié)構(gòu)參數(shù)相關(guān)的系數(shù)矩陣和與荷載相關(guān)的系數(shù)矩陣求得.而與結(jié)構(gòu)相關(guān)的系數(shù)矩陣可通過精細(xì)積分法或Newmark-β法求得.根據(jù)隨機(jī)響應(yīng)的方差靈敏度的一般計算公式[14]:

(21)

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

(22)

式中,

(23)

2 數(shù)值算例

以下將通過平面框架結(jié)構(gòu)和桁架結(jié)構(gòu)的算例對比兩種積分格式在不同時間步長下所求得的響應(yīng)及靈敏度,并分析時間步長對兩種積分格式精度和效率的影響.

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

(24)

其中,t1=3s,t2=6s,t3=15s,c=0.157 2.f(t)為一零均值高斯平穩(wěn)隨機(jī)過程,其功率譜密度函數(shù)取為Kanai-Tajimi過濾白噪聲譜,即

(25)

式中ωg=14rad/s(地基土的卓越頻率),ζg=0.6(地基土的阻尼比),S0=0.07N2/Hz(譜強(qiáng)度因子),其功率譜密度圖如圖2所示.設(shè)計參數(shù)分別取為柱邊長、梁高、材料的彈性模量.計算步長分別取Δt=0.005,0.010,0.020,0.050s,則A點(diǎn)的位移方差靈敏度時程曲線如圖3所示.

圖1 平面框架結(jié)構(gòu)模型(單位：m)

圖2 功率譜密度函數(shù)

圖3 A點(diǎn)的位移方差的靈敏度時程曲線

Fig.3 Sensitivity time history of the displacement variance of nodeA

表1給出了A點(diǎn)位移方差對彈性模量的靈敏度時程.從表1中可以看出,當(dāng)積分步長為0.005 s時,基于兩種積分格式的時域顯示法求得的位移方差及靈敏度一致,以此作為精確解.隨著積分步長的變大兩種積分格式的精度受到不同程度的影響.若以最大相對誤差在5%以內(nèi)為界限,則Newmark-β法允許的最大積分步長為0.020 s,而精細(xì)積分法允許的最大積分步長為0.050 s,兩種方法所用時間比為1.4∶0.2(即7倍).此時在相同精度下基于精細(xì)積分法的時域顯式法計算效率更高,而此時結(jié)構(gòu)的主頻正處于荷載的主頻范圍內(nèi).

表1 A點(diǎn)位移方差對彈性模量的靈敏度時程

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

某輸電塔結(jié)構(gòu)如圖4所示,該結(jié)構(gòu)有限元模型采用24個平面桿單元,共20個自由度,材料彈性模量E=2.1×1011N/m2,質(zhì)量密度ρ=7 800 kg/m3,考慮輸電線的自重,節(jié)點(diǎn)9和節(jié)點(diǎn)12具有集中質(zhì)量m=1 200 kg,其余節(jié)點(diǎn)具有集中質(zhì)量m=600 kg.阻尼采用Rayleigh阻尼,其中的阻尼系數(shù)α1=0.2,α2=0.005,結(jié)構(gòu)前三階頻率分別為35.7、128.3、136.6 rad/s.考慮非平穩(wěn)地震作用與算例2.1節(jié)相同.設(shè)計參數(shù)分別取為節(jié)點(diǎn)8的x坐標(biāo)和21號桿截面面積A1.計算步長分別取Δt=0.002,0.005,0.010,0.020 s,則節(jié)點(diǎn)12的水平位移方差靈敏度時程曲線如圖5所示.

圖4 平面桁架結(jié)構(gòu)模型(單位：m)

圖5 節(jié)點(diǎn)12的位移方差靈敏度時程曲線(單位：m)

Fig.5 Sensitivity time history of the displacement variance of node 12(Unit:m)

由以上推導(dǎo)已求得結(jié)構(gòu)的位移和速度的方差及相應(yīng)的靈敏度,若關(guān)注的是桁架結(jié)構(gòu)在隨機(jī)地震下的應(yīng)力方差及相應(yīng)的靈敏度則可進(jìn)一步推導(dǎo)第j根桿的第i時刻的應(yīng)力為

(26)

式中,

(27)

相應(yīng)的第j根桿的第i時刻的應(yīng)力靈敏度為

(28)

則第j根桿的第i時刻的應(yīng)力方差為

(29)

相應(yīng)的第j根桿的第i時刻的應(yīng)力方差靈敏度為

(30)

式中,

(31)

以下選取21號桿進(jìn)行應(yīng)力方差靈敏度的求解,設(shè)計參數(shù)同上,結(jié)果如圖6所示.

同樣,當(dāng)積分步長為0.002 s時,兩種積分格式求得的位移方差及靈敏度一致,以此作為精確解.表2列出了21號桿截面面積在某些時刻點(diǎn)處對其

圖6 21號桿件應(yīng)力方差的靈敏度時程曲線

應(yīng)力方差的靈敏度值,若以最大相對誤差在5%以內(nèi)為界限,則Newmark-β法允許最大的積分步長為0.02 s,精細(xì)積分允許最大的積分步長為0.01 s.兩種方法所需時間比為33.8∶5.1(即6.62倍).即在相同的精度下,基于Newmark-β法的時域顯示法計算效率更高,而此時的結(jié)構(gòu)主頻偏離荷載主頻范圍.

表2 21號桿的應(yīng)力方差對其截面面積的靈敏度時程

3 結(jié)論

文中用兩種不同的數(shù)值算法推導(dǎo)了結(jié)構(gòu)動力響應(yīng)靈敏度的兩種時域顯式格式,通過兩個數(shù)值算例的對比發(fā)現(xiàn):在求解非平穩(wěn)隨機(jī)地震響應(yīng)的靈敏度過程中,當(dāng)結(jié)構(gòu)的主頻位于地震荷載主頻范圍內(nèi)時,在同等精度條件下,基于精細(xì)積分法的時域顯式的靈敏度計算效率更高;而在結(jié)構(gòu)的主頻偏離地震荷載主頻范圍時,基于Newmark-β法的靈敏度計算效率更高.因此,為了提高結(jié)構(gòu)動力響應(yīng)優(yōu)化的計算效率,在保證精度的前提下,對不同結(jié)構(gòu)選用具體的數(shù)值計算方法時必須考慮結(jié)構(gòu)的主頻與地震荷載頻率的分布特性.文中結(jié)論可以廣泛應(yīng)用于結(jié)構(gòu)在隨機(jī)荷載下的動力靈敏度分析，為數(shù)值積分算法的選擇提供依據(jù).

[1] 姚昌仁,麻永平.結(jié)構(gòu)隨機(jī)激勵的響應(yīng)靈敏度分析 [J].力學(xué)學(xué)報,1990,22(1):438- 445. YAO Chang-ren，MA Yong-ping.The response sensitivity analysis for structural systems in random excitation [J].Chinese Journal of Theoretical and Applied Mechanics，1990,22(1):438- 445.

[2] 王鳳陽.基于虛擬激勵的隨機(jī)振動靈敏度分析及其應(yīng)用 [D].大連:大連理工大學(xué),2010.

[3] 唐帆,王錫平,朱文海,等.多源隨機(jī)激勵系統(tǒng)的參數(shù)靈敏度分析 [J].振動與沖擊,2012,31(1):82- 85. TANG Fan，WANG Xi-ping，ZHU Wen-hai，et al.Parameter sensitivity analysis for a system with multi-source random excitation [J].Journal of Vibration and Shock，2012,31(1):82- 85.

[4] LIU Q.Sensitivity and Hessian matrix analysis of power spectral density functions for uniformly modulate devolutionary random seismic responses [J].Finite Elements in Analysis and Design，2012,48:1370- 1375.

[5] 喬紅威,呂震宙.隨機(jī)結(jié)構(gòu)隨機(jī)激勵下的響應(yīng)靈敏度分析 [J].振動與沖擊,2008(3):60- 62,179. QIAO Hong-wei，Lü Zhen-zhou.Response sensitivity analysis of stochastic structures under random excitation [J].Journal of Vibration and Shock，2008(3):60- 62,179.[6] ZHAO Y，WANG F.Sensitivity analysis for structures subjected to stationary random excitations [J].Procedia Engineering，2011,14:2475- 2480.

[7] 蘇成,徐瑞.非平穩(wěn)激勵下結(jié)構(gòu)隨機(jī)振動時域分析法[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.

[8] 蘇成,徐瑞.非平穩(wěn)隨機(jī)激勵下結(jié)構(gòu)體系動力可靠度時域解法 [J].力學(xué)學(xué)報，2010,42(3):512- 520. SU Cheng,XU Rui.Time-domain method for dynamicreliability of structural systems subjected to non-stationary random excitations [J].Chinese Journal of Theoretical and Applied Mechanics，2010,42(3):512- 520.

[9] 陳太聰,蘇成,胡智強(qiáng),等.非平穩(wěn)隨機(jī)響應(yīng)靈敏度分析的時域顯式法 [J].振動工程學(xué)報,2015,28(4):43- 49. CHEN Tai-cong，SU Cheng，HU Zhi-qiang，et al.An explicit time-domain method in sensitivity analysis of non-stationary stochastic responses [J].Journal of Vibration Engineering，2015,28(4):43- 49.

[10] 鐘萬勰.結(jié)構(gòu)動力方程的精細(xì)時程積分法 [J].大連理工大學(xué)學(xué)報,1994,34(2):131- 136. ZHONG Wan-xie.On precise time-integration method for structural dynamics [J].Journal of Dalian University of Technology，1994,34(2):131- 136.

[11] NEWMARK N M.A method of computation for structural dynamics [J].Journal of Engineering Mechanics Division，1959,85(3):67- 94.

[12] 余超.基于時域顯式法的結(jié)構(gòu)瞬時最優(yōu)控制 [D].廣州:華南理工大學(xué)，2015.

[13] SU Cheng，HUANG Huan，MA Hai-tao，et al.Random vibration analysis of hysteretic systems by explicit iteration Monte-Carlo simulation method [J].Earthquakes and Structures，2014,7(2):119- 139.

[14] CACCIOLA P，COLAJANNI P，MUSCOLINO G.A modal 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:4344- 4361.

Sensitivity Analysis of Non-Stationary Random Seismic Response Based on Time-Domain Explicit Method

LIXue-pingLIDong-hongSUChengWEIPeng

(School of Civil Engineering and Transportation∥State Key Laboratory of Subtropical Building Science，South China University of Technology，Guangzhou 510640，Guangdong，China)

In this paper，two time-domain explicit expressions of the non-stationary random seismic response sensitivity are derived respectively based on the precise integration format and the Newmark-βintegral format，and they are respectively used to analyze the sensitivity of a plane frame and a plane truss structure.Then，the effects of the integral time step on the calculation accuracy and the efficiency of the two integral formats are investigated.It is found that，when the main frequency of the structure is within the range of the main load frequency，the time-domain explicit expression based on the precise integration has a higher efficiency with the same accuracy，but when the main frequency of the structure deviates from the main load frequency range，the time-domain explicit expression based on the Newmark-βintegral format achieves a higher efficiency.The research achievements can provide an effective reference for the choice of the numerical integration algorithms used in the structure optimization consi-dering the non-stationary random vibration.

sensitivity analysis；time-domain explicit method；non-stationary random seismic response

2015- 10- 19

國家自然科學(xué)基金資助項(xiàng)目(11372004);華南理工大學(xué)中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(2014ZZ0071) Foundation items: Supported by the National Natural Science Foundation of China(11002056,11372004)

李雪平(1978-),男,博士,副研究員,主要從事結(jié)構(gòu)振動、控制、優(yōu)化與系統(tǒng)可靠度分析研究.E-mail:xueping@scut.edu.cn

? 通信作者: 魏鵬(1978-),男,博士,副教授,主要從事結(jié)構(gòu)優(yōu)化方法研究.E-mail:ctpwei@scut.edu.cn

1000- 565X(2016)10- 0110- 07

O 324;TU311

10.3969/j.issn.1000-565X.2016.10.016