冼劍華 蘇成

摘要：分?jǐn)?shù)階導(dǎo)數(shù)模型是描述黏彈性材料本構(gòu)關(guān)系的理想模型。進(jìn)行了分?jǐn)?shù)階導(dǎo)數(shù)線性系統(tǒng)非平穩(wěn)隨機(jī)振動的靈敏度分析。建立分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)動力響應(yīng)的時域顯式表達(dá)式；采用靈敏度分析的直接求導(dǎo)法或伴隨變量法，推導(dǎo)系統(tǒng)動力響應(yīng)靈敏度的時域顯式表達(dá)式；提出分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)響應(yīng)統(tǒng)計矩靈敏度高效計算的時域顯式方法。所提出的基于直接求導(dǎo)法和伴隨變量法的時域顯式方法，分別適用于少設(shè)計變量和多設(shè)計變量兩種情況下的響應(yīng)統(tǒng)計矩靈敏度分析。以非平穩(wěn)地震激勵下設(shè)置分?jǐn)?shù)階導(dǎo)數(shù)黏彈性阻尼器的層剪切結(jié)構(gòu)為數(shù)值算例，驗證了所提方法的計算精度和計算效率。

關(guān)鍵詞：隨機(jī)振動靈敏度；分?jǐn)?shù)階導(dǎo)數(shù)；時域顯式方法；直接求導(dǎo)法；伴隨變量法

中圖分類號： O324；TU311.3??? 文獻(xiàn)標(biāo)志碼： A??? 文章編號：1004-4523（2022）05-1058-10

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

引言

理論和實驗研究表明，分?jǐn)?shù)階導(dǎo)數(shù)模型能夠同時模擬黏彈性材料的應(yīng)力松弛特性和蠕變特性，是描述黏彈性材料本構(gòu)關(guān)系的理想模型[1?2]。近年來，在結(jié)構(gòu)振動領(lǐng)域，分?jǐn)?shù)階導(dǎo)數(shù)模型已被廣泛用于描述黏彈性阻尼器的力學(xué)行為[3?5]。

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

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

近年來提出的一類非平穩(wěn)隨機(jī)振動時域顯式方法[22?25]，通過構(gòu)建結(jié)構(gòu)動力響應(yīng)及其靈敏度的時域顯式表達(dá)式，能夠在時域內(nèi)直接建立非平穩(wěn)響應(yīng)統(tǒng)計矩及其靈敏度的顯式列式，實現(xiàn)任意時刻和自由度的降維計算，并應(yīng)用于非平穩(wěn)隨機(jī)激勵下的結(jié)構(gòu)拓?fù)鋬?yōu)化，具有理想的計算精度和計算效率。在上述研究的基礎(chǔ)上，本文將時域顯式方法進(jìn)一步發(fā)展應(yīng)用于分?jǐn)?shù)階導(dǎo)數(shù)線性系統(tǒng)的非平穩(wěn)隨機(jī)振動靈敏度分析。首先構(gòu)建分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)動力響應(yīng)的時域顯式表達(dá)式；進(jìn)而采用靈敏度分析的直接求導(dǎo)法或伴隨變量法，推導(dǎo)系統(tǒng)動力響應(yīng)靈敏度的時域顯式表達(dá)式；最終利用統(tǒng)計矩的運算規(guī)律，建立系統(tǒng)非平穩(wěn)響應(yīng)統(tǒng)計矩靈敏度的時域顯式列式。以非平穩(wěn)地震激勵下設(shè)置分?jǐn)?shù)階導(dǎo)數(shù)黏彈性阻尼器的層剪切結(jié)構(gòu)為數(shù)值算例，驗證了所提方法的計算精度和計算效率。

1 分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)動力響應(yīng)時域顯式表達(dá)式

分?jǐn)?shù)階導(dǎo)數(shù)多自由度線性系統(tǒng)的運動方程可以表達(dá)如下：

式中? M，C，K 和 Cα分別為質(zhì)量矩陣、黏性阻尼矩陣、剛度矩陣和分?jǐn)?shù)階導(dǎo)數(shù)阻尼矩陣；U （ t ），U? （ t ）和 U? （ t ）分別為位移向量、速度向量和加速度向量； D αtU （ t ）表示對位移向量 U （ t ）關(guān)于時間 t 求α（0≤α<1）階導(dǎo)數(shù)；F （ t ）為非平穩(wěn)隨機(jī)激勵；L 為隨機(jī)激勵定位向量。

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

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

式中Γ（·）表示 Gamma 函數(shù)。利用 Gamma 函數(shù)的性質(zhì)，GL 系數(shù)可以用如下遞推關(guān)系式進(jìn)行計算：

定義分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)的狀態(tài)向量如下：

采用 Newmark?β數(shù)值積分法求解式（1），能夠推導(dǎo)得到系統(tǒng)狀態(tài)向量的遞推公式為（推導(dǎo)過程見附錄）：

式中wk =（Δ t ）-αGLk （0≤ k ≤ n ）； Vi = V （ ti ），Vi -1= V （ ti -1），Vi - k = V （ ti - k ），F(xiàn)i = F （ ti ），F(xiàn)i -1= F （ ti -1），其中ti = iΔ t，ti -1=（ i -1）Δ t，ti - k =（ i - k ）Δt；Q 1，Q 2，T 和 T1只與結(jié)構(gòu)參數(shù)有關(guān)，它們的表達(dá)式以及推導(dǎo)過程見附錄。

不失一般性，假定 V0= V （0）=0和 F0= F （0）=0，基于式（6）能夠推導(dǎo)得到系統(tǒng)狀態(tài)向量的時域顯式表達(dá)式為：

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

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

一般而言，人們只關(guān)注系統(tǒng)的某些關(guān)鍵響應(yīng)，并不需要求出系統(tǒng)中所有的響應(yīng)。假設(shè)ri為所關(guān)注的關(guān)鍵響應(yīng)，如位移響應(yīng)、速度響應(yīng)、位移響應(yīng)分?jǐn)?shù)階導(dǎo)數(shù)或構(gòu)件內(nèi)力響應(yīng)等，則由式（7）能夠得到關(guān)鍵響應(yīng)ri的時域顯式表達(dá)式為：

式中?? 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 為轉(zhuǎn)換向量，其中 qD，qV 和 q α分別為針對位移響應(yīng)、速度響應(yīng)和位移響應(yīng)分?jǐn)?shù)階導(dǎo)數(shù)的轉(zhuǎn)換向量。當(dāng)ri為 Vi 中的某一位移響應(yīng)、速度響應(yīng)或位移響應(yīng)分?jǐn)?shù)階導(dǎo)數(shù)時，q 中除與該位移響應(yīng)、速度響應(yīng)或位移響應(yīng)分?jǐn)?shù)階導(dǎo)數(shù)對應(yīng)的元素為1外，其余元素均為0。當(dāng)ri為某一構(gòu)件內(nèi)力響應(yīng)時，q 依賴于相應(yīng)的本構(gòu)關(guān)系。

2 基于直接求導(dǎo)法的動力響應(yīng)靈敏度時域顯式表達(dá)式

假設(shè)θ為分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)中的一個設(shè)計參數(shù)，則對運動方程（1）兩端同時關(guān)于參數(shù)θ求導(dǎo)能夠得到以下靈敏度方程：

式中? ?M ?θ，? C ?θ，?K ?θ，? Cα?θ和?L ?θ分別為矩陣 M，C，K，Cα 和 L 關(guān)于參數(shù)θ的靈敏度；? U （ t ）?θ，? U? （ t ）?θ，? U? （ t ）?θ和?D αtU （ t ）?θ分別為響應(yīng)向量 U（ t ），U? （ t ），U? （ t ）和 DαtU （ t ）關(guān)于參數(shù)θ的靈敏度。

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

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

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

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

假定系統(tǒng)的初始條件為 V0=0，并且初始條件與設(shè)計參數(shù)θ無關(guān)，即? V0?θ=0，則基于式（6）和（14）能夠推導(dǎo)得到系統(tǒng)狀態(tài)向量靈敏度的時域顯式表達(dá)式為：

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

式中：

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

記式（9）中關(guān)鍵響應(yīng)ri關(guān)于參數(shù)θ的靈敏度為?ri ?θ，則由式（15）可得?ri ?θ的時域顯式表達(dá)式為：

式中?? 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 ）。

假設(shè)考慮 m 個設(shè)計參數(shù)，則基于直接求導(dǎo)法構(gòu)建關(guān)鍵響應(yīng)ri關(guān)于所有參數(shù)靈敏度的時域顯式表達(dá)式，計算量相當(dāng)于對分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)進(jìn)行 m 次響應(yīng)靈敏度時程分析。當(dāng)所涉及的設(shè)計參數(shù)數(shù)目較多時，上述方法計算量較大，此時可基于伴隨變量法構(gòu)建關(guān)鍵響應(yīng)靈敏度的時域顯式表達(dá)式。

3 基于伴隨變量法的動力響應(yīng)靈敏度時域顯式表達(dá)式

設(shè) r （ t（～））為所關(guān)注t（～）時刻的位移響應(yīng)或構(gòu)件內(nèi)力響應(yīng)，則響應(yīng) r （ t（～））可以表達(dá)為以下積分形式：

式中? td 為位移響應(yīng) U（ t ）的持續(xù)時長；qD為針對位移響應(yīng)的轉(zhuǎn)換向量；δ（·）為 Dirac 函數(shù)。

為了計算響應(yīng) r （ t（～））關(guān)于參數(shù)θ的靈敏度，引入一任意伴隨向量λ（ t ），并定義一個新的變量ψ（ t（～））為：

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

對式（20）左右兩端同時關(guān)于參數(shù)θ求導(dǎo)，并進(jìn)行分部積分，整理后可得：

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

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

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

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

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

若關(guān)注時刻 t（～）= ti = iΔ t 的響應(yīng)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）可得響應(yīng) Ui 的時域顯式表達(dá)式為：

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

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

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

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

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

4 非平穩(wěn)隨機(jī)振動靈敏度分析的時域顯式方法

基于關(guān)鍵響應(yīng)ri的時域顯式表達(dá)式（9），及其靈敏度?ri ?θ的時域顯式表達(dá)式（18）或（32），利用統(tǒng)計矩的運算規(guī)律能夠直接計算得到響應(yīng)ri的均值μ ri和方差σri（2）關(guān)于參數(shù)θ的靈敏度為：

或

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

式中μF （ t ）和 RF （ t，τ）分別為非平穩(wěn)隨機(jī)激勵 F （ t ）的均值函數(shù)和自相關(guān)函數(shù)。

應(yīng)當(dāng)指出，式（35）～（38）為分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)響應(yīng)統(tǒng)計矩靈敏度的時域顯式列式，它建立了系統(tǒng)響應(yīng)統(tǒng)計量靈敏度與激勵統(tǒng)計量之間的直接聯(lián)系。式（35）和（36）可以稱作基于直接求導(dǎo)法的時域顯式方法，而式（37）和（38）可以稱作基于伴隨變量法的時域顯式方法。

5 數(shù)值算例

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

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

結(jié)構(gòu)受零均值非平穩(wěn)地震激勵 F（ t ）的作用， F （ t ）取為均勻調(diào)制非平穩(wěn)隨機(jī)過程，即 F （ t ）= g （ t ） f（ t ）。g （ t ）為均勻調(diào)制函數(shù)，取為：

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

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

式中：

相應(yīng)地，F(xiàn) （ t ）的自相關(guān)函數(shù)可以表達(dá)為：

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

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

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

兩種方法的計算時間如表2所示。從表中可以看出，當(dāng)考慮40個設(shè)計參數(shù)時，基于 AVM 的 ET ? DM 相比基于 DDM 的 ETDM 具有明顯的效率優(yōu)勢。這是因為在 ETDM 中，若基于 DDM 構(gòu)建響應(yīng)靈敏度的時域顯式表達(dá)式，計算量相當(dāng)于對分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)進(jìn)行40次響應(yīng)靈敏度時程分析；而基于 AVM 構(gòu)建響應(yīng)靈敏度的時域顯式表達(dá)式，則只需要進(jìn)行1次伴隨方程求解，計算量相當(dāng)于對分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)進(jìn)行1次響應(yīng)時程分析。

應(yīng)當(dāng)指出，本文所提方法還可以進(jìn)一步應(yīng)用于分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)非平穩(wěn)隨機(jī)振動優(yōu)化設(shè)計中。特別地，當(dāng)考慮分?jǐn)?shù)階黏彈性阻尼器性能參數(shù)優(yōu)化時，由于所涉及的設(shè)計參數(shù)數(shù)目通常較少，可以采用基于 DDM 的 ETDM 進(jìn)行靈敏度計算；當(dāng)考慮分?jǐn)?shù)階黏彈性阻尼器布局拓?fù)鋬?yōu)化時，由于所涉及的設(shè)計參數(shù)數(shù)目通常較多，可以采用基于 AVM 的 ETDM 進(jìn)行靈敏度計算。

6 結(jié)論

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

參考文獻(xiàn)：

[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，F(xiàn)ailla 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ǐn)?shù)階導(dǎo)數(shù)阻尼下非線性隨機(jī)振動結(jié)構(gòu)響應(yīng)的功率譜密度估計[ J ].應(yīng)用力學(xué)學(xué)報，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ǐn)?shù)階導(dǎo)數(shù)系統(tǒng)的非平穩(wěn)響應(yīng)[ J ].應(yīng)用數(shù)學(xué)和力學(xué)，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]蘇成，徐瑞.非平穩(wěn)隨機(jī)激勵下結(jié)構(gòu)體系動力可靠度時域解法[ J ].力學(xué)學(xué)報，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）的推導(dǎo)過程

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

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

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

式中：

式中γ和β為與積分穩(wěn)定性有關(guān)的參數(shù)，本文取γ=0.5，β=0.25。

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

由式（A6）可得：

式中：

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

U? i -1也可以類似表達(dá)為：

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

式中：

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

式中

式中：

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