分數階導數系統非平穩隨機振動靈敏度分析的時域顯式方法

冼劍華 蘇成

摘要：分數階導數模型是描述黏彈性材料本構關系的理想模型。進行了分數階導數線性系統非平穩隨機振動的靈敏度分析。建立分數階導數系統動力響應的時域顯式表達式；采用靈敏度分析的直接求導法或伴隨變量法，推導系統動力響應靈敏度的時域顯式表達式；提出分數階導數系統響應統計矩靈敏度高效計算的時域顯式方法。所提出的基于直接求導法和伴隨變量法的時域顯式方法，分別適用于少設計變量和多設計變量兩種情況下的響應統計矩靈敏度分析。以非平穩地震激勵下設置分數階導數黏彈性阻尼器的層剪切結構為數值算例，驗證了所提方法的計算精度和計算效率。

關鍵詞：隨機振動靈敏度；分數階導數；時域顯式方法；直接求導法；伴隨變量法

中圖分類號： O324；TU311.3??? 文獻標志碼： A??? 文章編號：1004-4523（2022）05-1058-10

DOI：10.16385/j .cnki .issn .1004-4523.2022.05.003

引言

理論和實驗研究表明，分數階導數模型能夠同時模擬黏彈性材料的應力松弛特性和蠕變特性，是描述黏彈性材料本構關系的理想模型[1?2]。近年來，在結構振動領域，分數階導數模型已被廣泛用于描述黏彈性阻尼器的力學行為[3?5]。

分數階導數系統的隨機振動分析已引起了不少學者的關注。Spanos 和 Zeldin[6]提出了分數階導數線性系統平穩隨機振動的頻域分析方法。Agraw ? al[7]給出了分數階導數單自由度線性系統平穩或非平穩隨機振動的時域解析解答。Di Paola 等[8]基于 Lyapunov 矩方程法求解了分數階導數線性振子的平穩或非平穩隨機響應。上述研究屬于線性隨機振動分析范疇。在非線性隨機振動分析方面，Spanos 和Evangelatos[9]利用蒙特卡羅模擬和統計線性化法分別求解了平穩白噪聲激勵下分數階導數非線性振子的響應統計量。孫春艷和徐偉[10]借助隨機平均法和統計線性化法，研究了分數階導數單自由度非線性系統在平穩白噪聲下的響應功率譜密度估計。 Xu 和 Li[11]采用概率密度演化法對設置分數階導數黏彈性阻尼器的單自由度非線性隨機結構進行了動力可靠度分析。更多關于分數階導數系統隨機振動分析的研究可參考文獻[12?14]。

靈敏度分析是結構優化、模型修正和結構損傷識別等問題面臨的重要課題。Kobelev[15]研究了分數階導數線性非保守系統的失穩臨界荷載靈敏度。 Martinez ? Agirre和Elejabarrieta[16]求解了分數階導數懸臂梁線性結構的特征值和特征向量靈敏度。 Lewandowski 和?aseckaPlura?[17]研究了設置分數階導數黏彈性阻尼器線性結構的動力特性靈敏度。Li 等[18]和 Yun 等[19]對黏彈性阻尼線性系統進行了動力響應靈敏度分析，他們的方法同樣適用于分數階導數系統。上述研究僅考慮了分數階導數線性系統動力特性或確定性動力響應的靈敏度，而針對分數階導數系統隨機振動靈敏度問題的研究目前尚未見到有文獻報道。另一方面，黏性阻尼線性系統隨機振動靈敏度問題的研究則相對成熟，相關研究可參考文獻[20?21]。

近年來提出的一類非平穩隨機振動時域顯式方法[22?25]，通過構建結構動力響應及其靈敏度的時域顯式表達式，能夠在時域內直接建立非平穩響應統計矩及其靈敏度的顯式列式，實現任意時刻和自由度的降維計算，并應用于非平穩隨機激勵下的結構拓撲優化，具有理想的計算精度和計算效率。在上述研究的基礎上，本文將時域顯式方法進一步發展應用于分數階導數線性系統的非平穩隨機振動靈敏度分析。首先構建分數階導數系統動力響應的時域顯式表達式；進而采用靈敏度分析的直接求導法或伴隨變量法，推導系統動力響應靈敏度的時域顯式表達式；最終利用統計矩的運算規律，建立系統非平穩響應統計矩靈敏度的時域顯式列式。以非平穩地震激勵下設置分數階導數黏彈性阻尼器的層剪切結構為數值算例，驗證了所提方法的計算精度和計算效率。

1 分數階導數系統動力響應時域顯式表達式

分數階導數多自由度線性系統的運動方程可以表達如下：

式中? M，C，K 和 Cα分別為質量矩陣、黏性阻尼矩陣、剛度矩陣和分數階導數阻尼矩陣；U （ t ），U? （ t ）和 U? （ t ）分別為位移向量、速度向量和加速度向量； D αtU （ t ）表示對位移向量 U （ t ）關于時間 t 求α（0≤α<1）階導數；F （ t ）為非平穩隨機激勵；L 為隨機激勵定位向量。

分數階導數的定義有很多種，其中最常用的有 Riemann ? Liouville （ RL ）定義、Caputo （C）定義和Grunward ? Letnikov （GL）定義。 GL 定義在數值計算中是最適用的，它可以表達為[26]：

式中Δ t 為時間步長；n 為時間步數；GLk為 GL 系數，它可以表達為[9]：

式中Γ（·）表示 Gamma 函數。利用 Gamma 函數的性質，GL 系數可以用如下遞推關系式進行計算：

定義分數階導數系統的狀態向量如下：

采用 Newmark?β數值積分法求解式（1），能夠推導得到系統狀態向量的遞推公式為（推導過程見附錄）：

式中wk =（Δ t ）-αGLk （0≤ k ≤ n ）； Vi = V （ ti ），Vi -1= V （ ti -1），Vi - k = V （ ti - k ），Fi = F （ ti ），Fi -1= F （ ti -1），其中ti = iΔ t，ti -1=（ i -1）Δ t，ti - k =（ i - k ）Δt；Q 1，Q 2，T 和 T1只與結構參數有關，它們的表達式以及推導過程見附錄。

不失一般性，假定 V0= V （0）=0和 F0= F （0）=0，基于式（6）能夠推導得到系統狀態向量的時域顯式表達式為：

式中? F[ i ]=[ F 1? F2? … Fi ]T；Ai =[ Ai，1? Ai，2? … Ai，i ]，其中系數向量 Ai，j （1≤ j ≤i≤ n ）可以由以下閉合公式進行計算：

根據式（8）所揭示的系數向量之間的內在關系，僅系數向量 Ai，1（1≤i≤ n ）需要計算和存儲，其余系數向量均可以用 Ai，1（1≤i≤ n ）表示。應當指出，系數向量 Ai，1具有明確的物理意義，它表示在 t1時刻作用的單位脈沖激勵f （ t ）（如圖1所示）下系統在ti時刻的狀態向量。因此，系數向量 Ai，1（1≤i≤ n ）的計算量相當于對分數階導數系統進行1次響應時程分析。

一般而言，人們只關注系統的某些關鍵響應，并不需要求出系統中所有的響應。假設ri為所關注的關鍵響應，如位移響應、速度響應、位移響應分數階導數或構件內力響應等，則由式（7）能夠得到關鍵響應ri的時域顯式表達式為：

式中?? a i（r）=[ a i（r），1? a i（r），2? …? a i（r），i ]，其中 a i（r），j = qT Ai，j （1≤ j ≤ i ≤ n ）；q =[ q D（T）?? q V（T）?? q α（T）]T 為轉換向量，其中 qD，qV 和 q α分別為針對位移響應、速度響應和位移響應分數階導數的轉換向量。當ri為 Vi 中的某一位移響應、速度響應或位移響應分數階導數時，q 中除與該位移響應、速度響應或位移響應分數階導數對應的元素為1外，其余元素均為0。當ri為某一構件內力響應時，q 依賴于相應的本構關系。

2 基于直接求導法的動力響應靈敏度時域顯式表達式

假設θ為分數階導數系統中的一個設計參數，則對運動方程（1）兩端同時關于參數θ求導能夠得到以下靈敏度方程：

式中? ?M ?θ，? C ?θ，?K ?θ，? Cα?θ和?L ?θ分別為矩陣 M，C，K，Cα 和 L 關于參數θ的靈敏度；? U （ t ）?θ，? U? （ t ）?θ，? U? （ t ）?θ和?D αtU （ t ）?θ分別為響應向量 U（ t ），U? （ t ），U? （ t ）和 DαtU （ t ）關于參數θ的靈敏度。

由式（10）可以看出，靈敏度方程的求解依賴于運動方程（1）的求解。由式（1）可將 U? （ t ）表達為：

將式（11）代入式（10）可得：

式中? V （ t ）如式（5）所示；L 1和 L2可以表達為：

對比式（1）和式（12）可知，靈敏度方程和運動方程在形式上是一致的。因此，與式（6）類似，系統狀態向量靈敏度的遞推公式可以表達為：

假定系統的初始條件為 V0=0，并且初始條件與設計參數θ無關，即? V0?θ=0，則基于式（6）和（14）能夠推導得到系統狀態向量靈敏度的時域顯式表達式為：

式中 F[ i ]=[ F 1 F2 … Fi ]T；Bi =[ Bi，1 Bi，2 … Bi，i ]，其中系數向量 Bi，j （1≤ j ≤i≤ n ）可以由以下閉合公式進行計算：

式中：

與式（8）類似，式（16）同樣揭示了系數向量 Bi，j （1≤ j ≤i≤ n ）之間的內在關系，僅 Bi，1（1≤i≤ n ）需要進行計算和存儲，其余系數向量均可以用 Bi，1（1≤i≤ n ）表示。與 Ai，1類似，系數向量 Bi，1同樣具有明確的物理意義，它表示在 t1時刻作用的單位脈沖激勵f （ t ）（如圖1所示）下，系統在ti時刻的狀態向量靈敏度。因此，系數向量 Bi，1（1≤i≤ n ）的計算量相當于對分數階導數系統進行1次響應靈敏度時程分析。

記式（9）中關鍵響應ri關于參數θ的靈敏度為?ri ?θ，則由式（15）可得?ri ?θ的時域顯式表達式為：

式中?? b i（r）=[ b i（r），1?? b i（r），2? …? b i（r），i ]，其中 b i（r），j = qT Bi，j （1≤ j ≤ i ≤ n ）。

假設考慮 m 個設計參數，則基于直接求導法構建關鍵響應ri關于所有參數靈敏度的時域顯式表達式，計算量相當于對分數階導數系統進行 m 次響應靈敏度時程分析。當所涉及的設計參數數目較多時，上述方法計算量較大，此時可基于伴隨變量法構建關鍵響應靈敏度的時域顯式表達式。

3 基于伴隨變量法的動力響應靈敏度時域顯式表達式

設 r （ t（～））為所關注t（～）時刻的位移響應或構件內力響應，則響應 r （ t（～））可以表達為以下積分形式：

式中? td 為位移響應 U（ t ）的持續時長；qD為針對位移響應的轉換向量；δ（·）為 Dirac 函數。

為了計算響應 r （ t（～））關于參數θ的靈敏度，引入一任意伴隨向量λ（ t ），并定義一個新的變量ψ（ t（～））為：

由于運動方程（1）在任意時刻恒成立，ψ（ t（～））恒等于 r （ t（～）），所以它們關于參數θ的靈敏度也恒等，即：

對式（20）左右兩端同時關于參數θ求導，并進行分部積分，整理后可得：

盡管式（22）對任意的伴隨向量λ（ t ）都成立，但為了避免計算位移向量的靈敏度?U （ t ）?θ，可選擇一伴隨向量以消除其中含?U （ t ）?θ的積分項，得到以下伴隨方程：

假定系統的初始位移和速度與設計參數θ無關，則有? U （0）?θ=? U? （0）?θ=0。同時，令式（22）的后三項均等于0，得到以下終值條件：

此時，式（23）和（24）組成一個關于伴隨向量λ（ t ）的終值問題。為求解該問題，可采用變量代換 s = td - t 將終值問題轉換成初值問題，即：

式中? P （ s ）=- qDδ（ s - s（～））表示在時刻 s = s（～）= td - t（～）作用的脈沖激勵；伴隨向量Λ（ s ）=λ（ td - s ）； D s（α）Λ（ s ）表示對Λ（ s ）關于 s 求α（0≤α<1）階導數，可表達為：

對比伴隨方程（25）和運動方程（1）可知兩者在形式上是一致的，所以式（25）同樣可以采用 New ? mark?β數值積分法進行求解。一旦獲得伴隨向量Λ（ s ）=Λ（ td - t ）=λ（ t ）后，由式（21）和（22）即可得到響應 r （ t（～））的靈敏度為：

若關注時刻 t（～）= ti = iΔ t 的響應ri = r （ ti ）（1≤i≤ n ），為簡單起見，可令 td = ti +1=（ i +1）Δ t，此時 P （ s ）表示在時刻 s =Δt作用的脈沖激勵。假定 U （0）= U? （0）=0和 F（0）=0，由式（1）和（2）可知 U? （ 0）=0和 D αtU （0）=0，又因為Λ（0）=0，所以式（27）可以采用梯形積分公式求解如下：

式中Λ i - k +1=Λ（ ti - k +1）（1≤ k ≤ i ≤ n ）。

由式（11）可知：

將式（29）代入式（28）可得：

式中? Vk =[ UkT?? U? kT?? D αtUkT ]T （1≤ k ≤ i ≤ n ）；L 1和 L2如式（13）所示。

由式（7）可得響應 Ui 的時域顯式表達式為：

式中? A i（U）=[ A i，（U）1? A i，（U）2? …? A i，（U）i ]，其中 A i，（U）j 取自系數向量 Ai，j （1≤ j ≤ i ≤ n ）中相應于 Ui 的列向量。

將式（7）和（31）代入式（30）可得響應靈敏度?ri ?θ的時域顯式表達式為：

式中? b i =[ b i，1?? b i，2? …? b i，i ]，其中：

由式（8）所揭示的系數向量 Ai，j （1≤ j ≤i≤ n ）之間的內在關系，可將式（33）改寫成：

從式（34）可以看出，僅系數 b（～） i（r），1（1≤i≤ n ）需要進行計算和存儲，其余系數均可以用 b（～） i（r），1（1≤i≤ n ） 表示。在已獲得系數向量 Ai，j （1≤ j ≤i≤ n ）的基礎上，若求得伴隨向量Λ（ s ），即可基于式（34）計算系數b（～） i（r），j （1≤ j ≤i≤ n ）并構建響應靈敏度?ri ?θ的時域顯式表達式（32）。需要指出的是，基于伴隨變量法構建響應ri關于不同參數靈敏度的時域顯式表達式，只需要進行1次伴隨方程求解，計算量相當于對分數階導數系統進行1次響應時程分析。因此，當所涉及的設計參數數目較多時，采用伴隨變量法構建響應靈敏度的時域顯式表達式會比采用直接求導法具有更高的計算效率，而當所涉及的設計參數數目較少時，采用直接求導法會更加簡單直接。

4 非平穩隨機振動靈敏度分析的時域顯式方法

基于關鍵響應ri的時域顯式表達式（9），及其靈敏度?ri ?θ的時域顯式表達式（18）或（32），利用統計矩的運算規律能夠直接計算得到響應ri的均值μ ri和方差σri（2）關于參數θ的靈敏度為：

或

式中? F[ i ]的均值向量和協方差矩陣可以表達為：

式中μF （ t ）和 RF （ t，τ）分別為非平穩隨機激勵 F （ t ）的均值函數和自相關函數。

應當指出，式（35）～（38）為分數階導數系統響應統計矩靈敏度的時域顯式列式，它建立了系統響應統計量靈敏度與激勵統計量之間的直接聯系。式（35）和（36）可以稱作基于直接求導法的時域顯式方法，而式（37）和（38）可以稱作基于伴隨變量法的時域顯式方法。

5 數值算例

以圖2所示設置分數階導數黏彈性阻尼器的層剪切結構為數值算例，驗證本文所提出的非平穩隨機振動靈敏度分析時域顯式方法的計算精度和計算效率。該結構每一層的質量和剛度分別為 mi =1.8×104 kg 和 ki =8.9×105 kN m （1≤i≤ N ），其中 N 為結構層數，在本算例中取 N =20。采用瑞利阻尼模型，取結構第1階和第20階模態的阻尼比為ζ=0.05。結構每一層均布置1個黏彈性阻尼器，阻尼器與樓層的夾角均為η= arccos 0.8。所布置黏彈性阻尼器的阻尼力模型取為分數階 Kelvin 模型[17]，即：

式中ud，i為第i個黏彈性阻尼器兩端節點的軸向相對位移；kd，i和 cd，i分別為第i個黏彈性阻尼器的剛度系數和阻尼系數。在本算例中，取α=0.6，kd，i = kd =3×105 kN/m和cd，i = cd =2.5×105 kN? sα/m （1≤i≤20）。

結構受零均值非平穩地震激勵 F（ t ）的作用， F （ t ）取為均勻調制非平穩隨機過程，即 F （ t ）= g （ t ） f（ t ）。g （ t ）為均勻調制函數，取為：

式中δ=0.18，ta =6 s，tb =18 s，tc =30 s 。f （ t ）為零均值平穩隨機過程，其功率譜密度函數取為Kanai Tajimi譜?[27]，即：

式中ω g =15.708 rad/s，ζg =0.6，S0=0.005 m 2/s3。 f （ t ）的自相關函數可以表達為[28]：

式中：

相應地，F （ t ）的自相關函數可以表達為：

考慮各層黏彈性阻尼器的剛度系數和阻尼系數均同時變化，此時僅有2個設計參數，分別為kd和 cd。分別采用基于直接求導法（Direct Differentiation Method ，DDM）和伴隨變量法 （Adjoint? Variable Method，AVM）的時域顯式方法（Explicit Time?Do ? main Method，ETDM）對圖2所示層剪切結構進行非平穩隨機振動靈敏度分析。同時，采用基于蒙特卡羅模擬（Monte? Carlo Simulation，MCS）的有限差分法（Finite Difference Method，FDM）計算隨機振動靈敏度的參考解，其中差分步長取為設計參數的0.2%變化量，MCS 的樣本數取為104。時程分析步長取Δt =0.02 s，時程分析總長為 T =30 s 。三種計算方法下結構頂層水平位移標準差關于kd和 cd 的靈敏度結果分別如圖3和4所示。從圖中可以看出，基于 DDM 的 ETDM 和基于 AVM 的 ETDM 計算結果基本重合，且均與基于 MCS 的 FDM 計算結果吻合良好，說明所提方法具有理想的計算精度。此外，從圖3和圖4還可以看出，頂層水平位移標準差靈敏度時程的變化趨勢明顯受到式（42）所示均勻調制函數的影響，說明在非平穩隨機激勵下分數階導數系統的隨機響應靈敏度具有明顯的非平穩特征。

上述三種方法的計算時間如表1所示。從表中可以看出，在設計參數只有2個的情況下，基于 DDM 的 ETDM 和基于 AVM 的 ETDM 計算時間分別為2.3和2.1 s，均遠少于基于 MCS 的 FDM 計算時間，說明所提方法具有理想的計算效率。這是因為在基于 MCS 的 FDM 中，需要對分數階導數系統進行大量的響應時程分析。而在 ETDM 中，若基于DDM 構建響應靈敏度的時域顯式表達式，計算量僅相當于對分數階導數系統進行2次響應靈敏度時程分析；若基于 AVM 構建響應靈敏度的時域顯式表達式，只需要進行1次伴隨方程求解，計算量相當于對分數階導數系統進行1次響應時程分析。

為了進一步考察所提方法在設計參數數目較多時的計算效率，考慮各層黏彈性阻尼器的剛度系數和阻尼系數各自不同變化，此時設計參數為kd，i和 cd，i （1≤i≤20），共計40個。采用基于 DDM 的 ET ? DM 和基于 AVM 的 ETDM 分別計算頂層水平位移標準差最大值關于剛度系數kd，i和阻尼系數 cd，i （1≤i≤20）的靈敏度，計算結果分別如圖5和6所示。從圖5和6可以看出，上述兩種方法的計算結果基本重合，這與圖3和4所觀察到的現象一致，同時還可以看出頂層水平位移標準差的最大值對底層黏彈性阻尼器剛度系數和阻尼系數的變化更敏感。

兩種方法的計算時間如表2所示。從表中可以看出，當考慮40個設計參數時，基于 AVM 的 ET ? DM 相比基于 DDM 的 ETDM 具有明顯的效率優勢。這是因為在 ETDM 中，若基于 DDM 構建響應靈敏度的時域顯式表達式，計算量相當于對分數階導數系統進行40次響應靈敏度時程分析；而基于 AVM 構建響應靈敏度的時域顯式表達式，則只需要進行1次伴隨方程求解，計算量相當于對分數階導數系統進行1次響應時程分析。

應當指出，本文所提方法還可以進一步應用于分數階導數系統非平穩隨機振動優化設計中。特別地，當考慮分數階黏彈性阻尼器性能參數優化時，由于所涉及的設計參數數目通常較少，可以采用基于 DDM 的 ETDM 進行靈敏度計算；當考慮分數階黏彈性阻尼器布局拓撲優化時，由于所涉及的設計參數數目通常較多，可以采用基于 AVM 的 ETDM 進行靈敏度計算。

6 結論

基于靈敏度分析的直接求導法和伴隨變量法，分別推導了分數階導數線性系統動力響應靈敏度的時域顯式表達式，提出了系統非平穩隨機振動靈敏度分析的時域顯式方法，可快速獲取系統關鍵響應統計矩的靈敏度。數值算例結果表明，所提出的基于直接求導法或伴隨變量法的時域顯式方法具有理想的計算精度和計算效率。當設計參數數目較少時，采用基于直接求導法的時域顯式方法更加簡單直接；當設計參數數目較多時，采用基于伴隨變量法的時域顯式方法計算效率更高。本文所提方法可與非線性隨機振動等效線性化法結合，用以求解分數階導數非線性系統的非平穩隨機振動靈敏度問題，有待進一步研究。

參考文獻：

[1]? Bagley R L，Torvik P J . A theoretical basis for the ap?plication of fractional calculus to viscoelasticity[ J ]. Jour? nal of Rheology，1983，27（3）：201-210.

[2]? Di Paola? M ，Pirrotta? A，Valenza? A . Visco-elastic be?havior through fractional calculus：an easier method for best fitting experimental results[ J ]. Mechanics of Mate ? rials，2011，43（12）：799-806.

[3]? Makris? N ， Constantinou? M? C . Fractional? derivativeMaxwell? model? for? viscous? dampers [ J ]. Journal? of Structural Engineering，1991，117（9）：2708-2724.

[4]? Shen K L，Soong T T . Modeling of viscoelastic damp?ers for structural applications[ J ]. Journal of Engineering Mechanics，1995，121（6）：694-701.

[5]? Lewandowski R，Pawlak Z . Dynamic analysis of frameswith viscoelastic dampers modelled by rheological mod? els with fractional derivatives[ J ]. Journal of Sound and Vibration，2011，330（5）：923-936.

[6]? Spanos P D，Zeldin B A . Random vibration of systemswith frequency-dependent parameters or fractional deriv? atives [ J ]. Journal? of? Engineering? Mechanics ， 1997，123：290-292.

[7]? Agrawal? O? P . Stochastic? analysis? of dynamic? systemscontaining? fractional? derivatives [ J ]. Journal? of? Sound and Vibration，2001，247（5）：927-938.

[8]? Di Paola M，Failla G，Pirrotta A . Stationary and non-stationary? stochastic response? of linear? fractional visco? elastic? systems [ J ]. Probabilistic? Engineering? Mechan? ics，2012，28：85-90.

[9]? Spanos P D，Evangelatos G I . Response of a non-linearsystem? with restoring? forces governed by? fractional de? rivatives-time domain simulation and statistical lineariza? tion? solution[ J ]. Soil? Dynamics? and? Earthquake? Engi? neering，2010，30（9）：811-821.

[10]孫春艷，徐偉.分數階導數阻尼下非線性隨機振動結構響應的功率譜密度估計[ J ].應用力學學報，2013，30（3）：401-405.

Sun C Y，Xu W . Response power spectral density esti? mate? of? a? fractionally? damped? nonlinear? oscillator [ J ].Chinese Journal of Applied Mechanics ，2013，30（3）：401-405.

[11] Xu J，Li J . Stochastic dynamic response? and reliabilityassessment of controlled structures with fractional deriv? ative model of viscoelastic dampers[ J ]. Mechanical Sys? tems and Signal Processing，2016，72-73：865-896

[12] Huang Z L，Jin X L，Lim C W，et al . Statistical analy?sis for stochastic systems including fractional derivatives [ J ]. Nonlinear Dynamics，2010，59（1-2）：339-349.

[13]李偉，趙俊鋒，李瑞紅，等. Guass白噪聲激勵下分數階導數系統的非平穩響應[ J ].應用數學和力學，2014，35（1）：63-70.

Li? W ，Zhao? J? F ，Li? R? H ，et? al . Non-stationary? re? sponse? of a? stochastic? system? with? fractional derivative damping under Gaussian white-noise excitation[ J ]. Ap? plied Mathematics and Mechanics，2014，35（1）：63-70.

[14] Fragkoulis? V? C ，Kougioumtzoglou? I? A ，Pantelous? AA，et al . Non-stationary response statistics of nonlinear oscillators with fractional derivative elements under evo? lutionary stochastic excitation[ J ]. Nonlinear Dynamics，2019，97：2291-2303.

[15] KobelevV . Sensitivity analysis of the linear nonconser?vative? systems? with? fractional? damping [ J ]. Structural and? Multidisciplinary? Optimization ， 2007， 33（3）：179-188.

[16] Martinez-Agirre M，Elejabarrieta M J . Higher order ei?gensensitivities-based numerical method for the harmon? ic analysis of viscoelastically damped structures[ J ]. In ? ternational Journal for Numerical Methods in Engineer? ing，2011，88（12）：1280-1296.

[17] Lewandowski? R ，?asecka-Plura? M . Design? sensitivityanalysis? of? structures? with? viscoelastic? dampers [ J ]. Computers and Structures，2016，164（1）：95-107.

[18] Li L，Hu Y J，Wang X L . Design sensitivity analysis ofdynamic response of nonviscously damped systems[ J ]. Mechanical? Systems? and? Signal? Processing ，2013，41（1-2）：613-638.

[19] Yun? K? S ，Youn? S? K . Design? sensitivity? analysis? fortransient? response? of? non-viscously? damped? dynamic systems[ J ]. Structural? and? Multidisciplinary? Optimiza? tion，2017，55（6）：2197-2210.

[20] Zhu M，Yang Y，Guest J K，et al . Topology optimiza?tion? for? linear? stationary? stochastic? dynamics ：applica? tions? to? frame? structures [ J ]. Structural? Safety ，2017，67：116-131.

[21] Gomez F，Spencer B F . Topology optimization frame?work for structures subjected to stationary stochastic dy? namicloads[ J ]. Structural and Multidisciplinary Optimi? zation，2019，59：813-833.

[22]蘇成，徐瑞.非平穩隨機激勵下結構體系動力可靠度時域解法[ J ].力學學報，2010，42（3）：512-520.

Su C，Xu R . Time-domain method for dynamic reliabili? ty of structural systems subjected to non-stationary ran? domexcitations[ J ]. Chinese Journal of Theoretical and Applied Mechanics，2010，42（3）：512-520.

[23] Su C，Xu R . Random vibration analysis of structures bya time-domain explicit formulation method[ J ]. Structur? al Engineering and Mechanics，2014，52（2）：239-260.

[24] Hu Z Q，Su C，Chen T C，et al . An explicit time-do?main? approach? for? sensitivity? analysis? of non-stationary random vibration problems[ J ]. Journal of Sound and Vi? bration，2016，382：122-139.

[25] Hu? Z? Q ，Wang? Z? Q ，Su? C ，et? al . Reliability? basedstructural topology optimization considering non-station? ary stochastic excitations[ J ]. KSCE Journal of Civil En ? gineering，2018，22（3）：993-1001.

[26] PodlubnyI . Fractional Differential Equations：An Intro ?duction? to? Fractional? Derivatives ， Fractional? Equa? tions，to Methods of Their Solution and Some of TheirApplications[M]. New York：Academic Press，1999.

[27] Kanai K . Semi-empirical formula for the seismic charac?teristics? of? the? ground [ J ]. Bulletin? of? the? Earthquake Research Institute，1957，35（2）：309-325.

[28] Sun? G? J ，Li H? J . Stationary? models? of random? earth?quake ground motion and their statistical properties[ J ]. Earthquake? Engineering? and? Engineering? Vibration，2004，24（6）：21-26.

[29] Newmark N W . A method of computation for structuraldynamics [ J ]. Journal of the Engineering Mechanics Di? vision，1959，85（7）：67-94.

Explicit time -domain method for sensitivity analysis of nonstationary random vibration of systems with fractional derivatives

XIAN Jian?hua1，SU Cheng1，2

（1.School of Civil Engineering and Transportation，South China University of Technology，Guangzhou 510640，China；

2.State Key Laboratory of Subtropical Building Science，South China University of Technology，Guangzhou 510640，China）

Abstract： Fractional derivative models are capable of describing the constitutive behaviors of viscoelastic materials . This paper is devoted to the sensitivity analysis of nonstationary random vibration of linear systems comprising fractional derivative terms . The explicit time-domain expressions of dynamic responses are firstly established for the system with fractional derivatives . The sensitiv ? ities of dynamic responses are then derived using the direct differentiation method （DDM） or the adjoint variable method （AVM）. On the basis of the explicit expressions of dynamic responses and their sensitivities，an explicit time-domain method （ETDM） is proposed for efficient calculation of the sensitivities of statistical moments of responses . The proposed DDM- and AVM-based ET ? DM are applicable to the scenarios with less and more design variables，respectively . A numerical example involving a shear-type structure under nonstationary seismic excitations and with viscoelastic dampers modelled by fractional derivatives is presented to validate the computational accuracy and efficiency of the proposed method .

Key words ： sensitivity of random vibration；fractional derivative；explicit time-domain method；direct differentiation method；ad? joint variable method

作者簡介：冼劍華（1994—），男，博士研究生。電話：13570433950；Email：jianhua .xian@qq .com。

通訊作者：蘇成（1968—），男，博士，教授。電話：（020）87111636；Email：cvchsu@scut .edu .cn。

附錄：公式（6）的推導過程

由式（1）可得ti時刻的運動方程為：

式中 D αtUi 可由式（2）得到：

Newmark?β數值積分法采用以下兩個基本假定[29]：

式中：

式中γ和β為與積分穩定性有關的參數，本文取γ=0.5，β=0.25。

將式（A2）～（A4）代入式（A1）可得：

由式（A6）可得：

式中：

由式（A1）可得 U? i為：

U? i -1也可以類似表達為：

將式（A10）代入式（A7）可得：

式中：

將式（A10）和（A11）代入式（A4）可得：

式中

式中：

將式（A11），（A13）和（A15）合并可得式（6），式中：