隨機(jī)Bagley-Torvik方程的非平穩(wěn)解析解

關(guān)鍵詞：隨機(jī)振動(dòng)；分?jǐn)?shù)階導(dǎo)數(shù)；Bagley-Torvik方程；完全非平穩(wěn)；Mittag-Leffler函數(shù) 中圖分類號(hào)：0324；0321 文獻(xiàn)標(biāo)志碼：A DOI：10.16385/j.cnki.issn.1004-4523.202309028

Non-stationary analytic solution of the stochastic Bagley-Torvik equation

KONG Fan1，XU Yijian1，GUO Wenjie2，HONG Xu1，CAO Hongyou3 （1.College of Civil Engineering，Hefei University of Technology，Hefei 23ooo9，China; 2.School of Transportation Engineering，East China Jiaotong University，Nanchang 33Ool3，China; 3.School of Civil Engineering and Architecture，Wuhan Universityof Technology，Wuhan 43Oo7o，China）

Abstract：TheBagley-Torvik（B-T）equationisadiferentialequationofmotionwithfractional（3/2）orderderivativeermsthat is applied todescribethemotioofarigidplateinNewtonian，viscousfluid.Inthispaper，wedevelop nonstationaryanalytic solu tionsoftheB-Tequationwhoseinhomogeneous termisastochasticprocess.TheB-Tequationistransformedintoahalforder state-space equation inmatrix formandeigen-analysis is performedtoobtaincomplex eigenvaluesand eigenvectors.Subsequently， the generalized coordinate transformation is introduced to decouple the equation intoa system of independent 1/2 -orderdifferential equations whicharesolvedbyLaplace transformtoobtainthesolutioningeneralizedcoordinates；Thegeneralizedcoordinate solu tionisconvertedintoanaturalcoordinatesolutiontoobtaintheimpulseorstepresponsefunction.Whentheinhomogeneous term oftheequationisastochasticprocess，theLaplace transformcanbeused toderivethetime-varying frequencyresponsefunction fromwhichtheanalytical solutionofthe nonstationarystochasticresponsecanbeobtainedbyrelyingontherelationship between theexcitationandtheresponse powerspectraldensity.Thecorectnessofthemethodisverifiedbynumericalcasesusingthe Spanos-Solomos fully non-statoionary stochastic excitation as an example.

Keywords：randomvibration;fractionalderivative;Bagley-Torvik equation;fullynon-stationary;Mittag-Leflerfunction

Bagley-Torvik（以下簡(jiǎn)稱為B-T）方程是一種帶3/2分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的微分方程，由BAGLEY和TORVIK提出[1]，用于描述剛性板在牛頓流體中的振動(dòng)狀態(tài)。注意到分?jǐn)?shù)階導(dǎo)數(shù)在描述材料力學(xué)性能方面有巨大的優(yōu)勢(shì)[2]，然而，由于其記憶特性，分?jǐn)?shù)階微分方程的求解過(guò)程相較整數(shù)階方程更為困難。

人們發(fā)展了若干求解確定性B-T方程的解析和數(shù)值方法。PODLUBNY3推導(dǎo)了以無(wú)窮級(jí)數(shù)表示的格林函數(shù)，得到了該方程在確定性外力作用下的響應(yīng)；SRIVASTAVA等4利用同倫分析方法（homotopyanalysismethod，HAM）得到了與POD-LUBNY一致的結(jié)果;JENA等[5]利用Sumudu變換得到幾種特殊激勵(lì)作用下方程的精確解，與已知理論解對(duì)比驗(yàn)證了該方法的正確性；黃瀟等[基于Laplace變換給出了幾類特殊分?jǐn)?shù)階微分方程的級(jí)數(shù)解；TRINKS等通過(guò)無(wú)記憶化方法（memory-freeformulation）將分?jǐn)?shù)階導(dǎo)數(shù)轉(zhuǎn)化為整數(shù)導(dǎo)數(shù)并利用Newmark算法得到了B-T方程的數(shù)值解；WANG等[8利用Mittag-Leffler表示的指數(shù)函數(shù)給出了B-T方程的通解。此外，研究者們還嘗試使用神經(jīng)網(wǎng)絡(luò)求解B-T方程。RAJA等9結(jié)合模式搜索技術(shù)與遺傳算法求解了確定性荷載作用下的數(shù)值解；胡行華等[10]嘗試用遺傳算法（genetic algorithm，GA）-Chebyshev神經(jīng)網(wǎng)絡(luò)求解該方程。此外，其他方法還包括Laplace變換[1]、積分變換法[12]、廣義泰勒配點(diǎn)法[13]等。雖然已有眾多方法可應(yīng)用于求解該方程，但多數(shù)情況均為數(shù)值解、近似解或以無(wú)窮級(jí)數(shù)表示的復(fù)雜形式精確解析解，僅少數(shù)特例情況下能獲得簡(jiǎn)單的精確解析解。

工程動(dòng)力系統(tǒng)往往處于隨機(jī)激勵(lì)作用之下，然而，目前很少有研究者關(guān)注隨機(jī)激勵(lì)作用下B-T方程的解。分?jǐn)?shù)階隨機(jī)動(dòng)力系統(tǒng)的求解方法包括數(shù)值解法[14]、隨機(jī)平均法[15]、無(wú)記憶化方法[16]、半解析方法[7等。然而，這些方法存在數(shù)值誤差或由于各種假設(shè)導(dǎo)致無(wú)法獲得精確解析解?？追驳萚19基于SUAREZ等[18]的研究，在隨機(jī)振動(dòng)頻域法的框架下利用時(shí)變頻響函數(shù)的概念和Laplace變換技術(shù)得到了均勻調(diào)制非平穩(wěn)激勵(lì)作用下含1/2階分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的隨機(jī)動(dòng)力系統(tǒng)的非平穩(wěn)響應(yīng)解析解。然而，該方法僅考慮分?jǐn)?shù)階導(dǎo)數(shù)為1/2的情況特殊情形，未能將分?jǐn)?shù)階導(dǎo)數(shù)擴(kuò)展到大于1的狀況；此外，如何得到完全非平穩(wěn)激勵(lì)作用下分?jǐn)?shù)階系統(tǒng)的隨機(jī)動(dòng)力響應(yīng)的解析解仍是一個(gè)未解問(wèn)題。

本文提出隨機(jī)B-T方程的一種非平穩(wěn)解析解。為此，須首先得到確定性B-T方程的脈沖響應(yīng)函數(shù)，通過(guò)以下步驟獲得：將B-T方程改寫為空間狀態(tài)方程并進(jìn)行特征值分析；引入特征矩陣的廣義坐標(biāo)變換將B-T方程解耦，得到關(guān)于廣義坐標(biāo)的4個(gè)獨(dú)立的半階微分方程并利用Laplace變換求解;將廣義坐標(biāo)變換為自然坐標(biāo)，即可得到齊次以及非齊次項(xiàng)為特殊函數(shù)（階躍和脈沖函數(shù)）時(shí)B-T方程的解析解。它是一種有限項(xiàng)求和的解析解，可視為已知無(wú)窮級(jí)數(shù)解3的特殊形式。在確定脈沖響應(yīng)的基礎(chǔ)上，結(jié)合Laplace變換，推導(dǎo)非齊次項(xiàng)為完全非平穩(wěn)激勵(lì)作用時(shí)方程的均方解析解，與蒙特卡羅（MonteCarlo）方法或已知精確解法所得結(jié)果的對(duì)比驗(yàn)證了所提出方法的正確性。

1數(shù)學(xué)方法

Bagley-Torvik方程： （1常用于描述剛性板在牛頓流體中的運(yùn)動(dòng)狀態(tài)。式中， m 表示質(zhì)量；表示剛度；D表示導(dǎo)數(shù)運(yùn)算符； A 為板的面積； μ 為流體黏度； ρ 為流體密度；f表示外部荷載； 分別表示響應(yīng)位移和加速度。分?jǐn)?shù)階導(dǎo)數(shù)定義為：

式中，下標(biāo)“C\"表示Caputo定義； τ 表示0到 t 之間的積分變量； Γ 為伽馬函數(shù)。令 并將式（1）改寫為狀態(tài)空間形式[20]：

BD^1/2z（t）-Az（t）=F（t）

其中，

且

與方程（3）對(duì)應(yīng)的特征值方程可寫為：

∣G-λI∣=0

GΨ_j=λ_jΨ_j;j=1，2，3，4

式中， λ_j 和 ψ_j 分別表示第 j 個(gè)特征值及對(duì)應(yīng)的特征向量； I 為 4×4 的單位方陣；

由于矩陣 A 和 B 均為對(duì)稱正定陣，經(jīng)過(guò)對(duì)角化操作可得到如下正交關(guān)系：

其中，

觀察式（8）可知，相同特征向量的不同元素之間存在如下關(guān)系：

Ψ_2j=λ_jΨ_1j，Ψ_3j=λ_jΨ_2j，Ψ_4j=λ_jΨ_3j

即

將式（14）代人式（10）可得：

為求解方程（1），引人變換：

式中， y（t） 為計(jì)算過(guò)程中的輔助變量。

將式（16）代人式（3），左乘 ψ^?T 并注意如式（10）和（11）所示的正交條件，可得4個(gè)解耦的半階微分方程：

對(duì)上式兩邊進(jìn)行Laplace變換并將所得結(jié)果代入式（16），得：

式中， ，r=0，1，2，3;L（?） 表示Laplace變換;此外，

式中， s 為時(shí)域變量 t 經(jīng)過(guò)拉普斯變換到拉氏域中的量。且 R_j=s^-1/2y_j（0） 。式（19）右邊分子為兩項(xiàng)疊加的形式，前者只與外部激勵(lì)有關(guān)，而后者只與初始條件相關(guān)。

1. 1 確定性分析

1.1. 1 自由振動(dòng)

考慮式（1）對(duì)應(yīng)的齊次方程，即 F（s）=0 的情況，此時(shí)式（19）變?yōu)椋?/p>

將該式代人式（18）可求得含 y_j（0） 的4個(gè)表達(dá)式：

其中，

式中，erf表示高斯誤差函數(shù)。

考慮 t=0 時(shí)刻處的初始條件，即

利用式（23）提供的4個(gè)線性方程可將 y_j（0） 與初始條件的關(guān)系寫為矩陣形式：

將式（24）兩邊左乘 ψ^{scriptscriptstyleT}B 并考慮式（10）可得：

由于速度和加速度的物理有界性，分?jǐn)?shù)階項(xiàng)初

始條件為零[21]。從式（25）中可知，初始速度和位移均為零時(shí)， y_j（0）=0 。將式（25）代人式（23）可得到自由振動(dòng)下的位移與速度響應(yīng)：

1. 1.2 階躍荷載

考慮零初始條件且 f（t）=f₀u（t） 的情況，其中f₀ 為階躍幅值， u（t） 為階躍函數(shù)?？紤]式（18）和（19）可得：

利用式（28）和（29）可計(jì)算TRINKS等關(guān)注的具有特定非齊次項(xiàng)的B-T方程的解，詳見(jiàn)數(shù)值算例2.1.2。

1. 1.3 脈沖激勵(lì)

考慮零初始條件且 f（t）=δ（t） 的情況，其中δ（t） 為單位脈沖。同樣地，考慮式（18）和（19）可求得脈沖響應(yīng)函數(shù)：

當(dāng) t=0 時(shí)，由于位移的物理有界性可消去上式第一項(xiàng)，故：

基于脈沖響應(yīng)函數(shù)，可利用卷積計(jì)算確定性激勵(lì)下系統(tǒng)響應(yīng)。脈沖響應(yīng)函數(shù)同時(shí)也是隨機(jī)分析的基礎(chǔ)。

1.2 隨機(jī)分析

工程動(dòng)力系統(tǒng)常受到完全非平穩(wěn)隨機(jī)過(guò)程作用。該隨機(jī)過(guò)程的功率譜密度會(huì)在時(shí)間和頻率兩個(gè)維度上同時(shí)發(fā)生變化?？紤]零均值完全非平穩(wěn)隨機(jī)過(guò)程，根據(jù)PRIESTLEY[22]的演變隨機(jī)過(guò)程理論，具有演變功率譜密度的完全非平穩(wěn)隨機(jī)過(guò)程可表示為：

式中， a（t，ω） 為隨時(shí)間和頻率變化的慢變調(diào)制函數(shù)； Z（ω） 為正交增量過(guò)程，且：

E[dZ（ω）]=0

式中，“*”為共軛運(yùn)算符； E（?） 表示期望。

結(jié)合式（34）可得式（32）所示非平穩(wěn)隨機(jī)過(guò)程的相關(guān)函數(shù)為：

利用式（35）中 δ 函數(shù)的篩選性質(zhì)可將雙重積分化為一維積分：

E[f（t₁）f（t₂）]=

根據(jù)線性系統(tǒng)的疊加原理，隨機(jī)B-T方程的解可用Duhamel積分表示，即

故均方響應(yīng)為：

將式（35）所示的激勵(lì)相關(guān)函數(shù)代入式（38）可得：

其中，

對(duì)式（4O）兩邊進(jìn)行Laplace變換，并考慮其卷積性質(zhì)可得：

式中，上標(biāo)“ ～ ”表示Laplace域，且

為脈沖響應(yīng)函數(shù)的Laplace變換。將時(shí)頻調(diào)制函數(shù)a（t，ω） 的具體形式代入式（41）后進(jìn)行Laplace逆變換，可得到時(shí)變頻響函數(shù)，然后利用式（39）即可得響應(yīng)均方值，詳見(jiàn)下文數(shù)值算例。

2 數(shù)值算例

本節(jié)進(jìn)一步給出B-T方程在非齊次項(xiàng)為確定性函數(shù)和完全非平穩(wěn)隨機(jī)過(guò)程時(shí)的數(shù)值算例以驗(yàn)證方法的準(zhǔn)確性。以下數(shù)值算例中的參數(shù)單位均采用國(guó)際單位制。其中單位為 N/m，m 單位為 kg ，力的單位為 N，c 的單位隨分?jǐn)?shù)階數(shù)發(fā)生變化，分?jǐn)?shù)階數(shù)為3/2時(shí)單位為 kg/s^1/2 ，階數(shù) 1/2 時(shí)為 kg/s^3/2 ，一階時(shí)為kg/s ，為方便起見(jiàn)數(shù)值算例中省略 c 的單位。

2.1 確定性分析

2.1.1 自由振動(dòng)

LABECCA等23給出了自由振動(dòng)情況下B-T方程的級(jí)數(shù)解；MAHMUDOV等24利用含有兩個(gè)級(jí)數(shù)項(xiàng)的Mittag-Leffler函數(shù)研究了初始條件不為零時(shí)B-T方程的齊次解。為驗(yàn)證本文所建議方法的正確性，考慮系統(tǒng)參數(shù) k=1N/m，m=1kg 且具有初始單位位移和單位速度的情況。將初始條件代入式（25），并將位移響應(yīng)式（26）和速度響應(yīng)式（27）分別繪制于圖1和2中。圖中PM表示本文建議方法。可見(jiàn)，本文方法得到的結(jié)果能完美地符合已有級(jí)數(shù)解[24]，驗(yàn)證了方法的正確性。然而，阻尼較大時(shí)，由后者得到的結(jié)果從特定時(shí)刻開始便出現(xiàn)發(fā)散情況，如圖2所示 c=1 時(shí)的速度時(shí)程尾部實(shí)心圓點(diǎn)處。此外，與整數(shù)階系統(tǒng)自由振動(dòng)響應(yīng)類似，分?jǐn)?shù)階系統(tǒng)響應(yīng)也為幅值衰減的諧波振動(dòng)：阻尼系數(shù)越大，位移和速度的振幅隨時(shí)間增長(zhǎng)下降速率越快。

考察分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)對(duì)系統(tǒng)響應(yīng)的影響。通過(guò)構(gòu)造與式（3）中不同的 A，B 矩陣，本文方法同樣可以應(yīng)用于 1/2 階以及整數(shù)階系統(tǒng)[20]。圖3給出了 c= 0.5時(shí)分?jǐn)?shù)階數(shù)對(duì)系統(tǒng)響應(yīng)的影響。由圖3可知，相同初始條件下，整數(shù)階系統(tǒng)的幅值衰減速率大于分?jǐn)?shù)階系統(tǒng)的幅值衰減速率。

2.1.2 階躍響應(yīng)

考慮初始條件 ，參數(shù) k= 0.5N/m ， m=1 kg且非齊次項(xiàng)為方波：

f（t）={8，0?t?1

的B-T方程。將式（43）寫為 f（t）=8[u（t）-u（t- 1]，將上述非齊次項(xiàng)激勵(lì)與初始條件代人式（28）和（29），并結(jié)合B-T方程所描述系統(tǒng)的線性與時(shí)不變性，得：

非齊次項(xiàng)如式（43）所示的B-T方程也可采用PODLUBNY[3]提出的格林函數(shù)與激勵(lì)的卷積積分求解，即

其中，

式中， E_λ，μ 為二參數(shù)Mittag-Leffler函數(shù)，右上標(biāo)表示求導(dǎo)次數(shù)。

繼PODLUBNY的工作后，許多研究者嘗試用其他方法對(duì)該方程求解。TRINKS等和CENE-SIZ等[13分別利用無(wú)記憶化方法和廣義泰勒配點(diǎn)法得到了該方程的數(shù)值解；ATANACKOVIC等[25]提出了一種求數(shù)值解的新方法；RAY等[26]利用Adomain分解法給出了與PODLUBNY解類似的雙重級(jí)數(shù)形式的解析解。因此，仍采用式（46）和（47）驗(yàn)證本文提出的方法。

圖4和5分別給出了本文方法與PODLUBNY方法在不同阻尼系數(shù)時(shí)得到的位移與速度時(shí)程。由圖可知，阻尼比較小時(shí)兩者得到的結(jié)果幾乎完全重合，進(jìn)一步表明本文方法的正確性。值得注意的是，張德茂等[2曾指出：雖然PODLUBNY方法在理論上收斂，但在實(shí)際計(jì)算中有可能會(huì)出現(xiàn)尾部異常發(fā)散的現(xiàn)象。從本文圖中也可發(fā)現(xiàn)這一點(diǎn)：當(dāng)阻尼系數(shù)為1且時(shí)間 tgt;15s 時(shí)，PODLUBNY方法出現(xiàn)了嚴(yán)重的偏離現(xiàn)象（已用實(shí)心標(biāo)識(shí)給出）。反之，本文提出的方法無(wú)需計(jì)算復(fù)雜的雙重級(jí)數(shù)乘積，不僅計(jì)算相對(duì)簡(jiǎn)單，還避免了上述可能出現(xiàn)的數(shù)值發(fā)散問(wèn)題。

為探究不同分?jǐn)?shù)階導(dǎo)數(shù)下響應(yīng)函數(shù)的關(guān)系，使m=1kg c=0.5 k=1N/m 且 f₀=1N ，并將參數(shù)代入式（28），在圖6中繪制不同分?jǐn)?shù)階導(dǎo)數(shù)時(shí)B-T方程的位移響應(yīng)。由圖6可知，整數(shù)階系統(tǒng)的位移響應(yīng)很快就趨近于靜位移， 3/2 階系統(tǒng)次之， 1/2 階系統(tǒng)接近靜位移的速率非常緩慢。

2.1.3 脈沖響應(yīng)

考慮B-T方程參數(shù)為 m=1kg，k=1N/m ，激勵(lì)為 f（t）=δ（t） 且初始條件為零的情況。將上述已知條件代人式（31）得到位移響應(yīng)（如圖7所示）。由圖7可知，本文方法很好地符合PODLUBNY給出的級(jí)數(shù)解（式（47））。同樣地，由于可能的數(shù)值與截?cái)嗾`差，阻尼系數(shù)為1時(shí)PODLUBNY方法給出的結(jié)果發(fā)生了偏離，如圖7中實(shí)心標(biāo)識(shí)所示。此外，隨著方程阻尼系數(shù)增大，脈沖響應(yīng)函數(shù)幅值隨時(shí)間呈快速衰減趨勢(shì)。圖8為 c=0.5 時(shí)不同階數(shù)下的系統(tǒng)脈沖響應(yīng)函數(shù)。與圖3觀察到的結(jié)論相同，整數(shù)階系統(tǒng)脈沖響應(yīng)函數(shù)的幅值衰減速率大于分?jǐn)?shù)階系統(tǒng)的幅值衰減速率。

2.2 隨機(jī)分析

考慮一類由SPANOS等[28]提出的完全非平穩(wěn)隨機(jī)過(guò)程模型：

S（t，ω）=S₀α（ω）t²e^{-[b+β（ω）]t}

式中， α（ω） 和 β（ω） 分別為與頻率相關(guān)的函數(shù)； S₀ 為一常量平穩(wěn)功率譜密度； b 為控制參數(shù)。該模型也可表示為：

式中， a（t，ω） 為慢變調(diào)制函數(shù)，且：

Φ（ω）=S₀α（ω）

將式（50）代人式（41）并進(jìn)行Laplace變換，得到Laplace域中的時(shí)變頻響函數(shù)：

其中，

為便于進(jìn)行Laplace逆變換，將式（52）寫為部分分式形式：

對(duì)式（54）進(jìn)行Laplace逆變換即可得到時(shí)變頻響函數(shù)的解析解：

其中，

將式（55）代人式（39）并使 t₁=t₂=t 可得到響應(yīng)功率譜密度，進(jìn)而通過(guò)數(shù)值積分得到均方響應(yīng)。選取參數(shù)為 的Spanos-Solomos完全非平穩(wěn)隨機(jī)過(guò)程，其演變功率譜密度如圖9所示。選取系統(tǒng)參數(shù) k=1N/m m=1kg，c=1 的隨機(jī)B-T方程作為演示算例。圖10顯示了分?jǐn)?shù)階導(dǎo)數(shù)不同時(shí)，利用本文方法與10000條樣本的MonteCarlo模擬（MCS）得到的響應(yīng)均方值對(duì)比，其中，樣本響應(yīng)采用分?jǐn)?shù)階運(yùn)動(dòng)方程的New-mark法計(jì)算得到?？梢?jiàn)，本文建議方法能很好地符合模擬結(jié)果；另外，分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)大于整數(shù)階導(dǎo)數(shù)系統(tǒng)的響應(yīng)均方值。圖11和12分別顯示了由本文方法得到和基于MonteCarlo模擬估計(jì)的響應(yīng)演變功率譜密度三維圖。其中，后者方法利用Morlet小波（MOwavelet）[29]計(jì)算樣本響應(yīng)的小波變換，進(jìn)而通過(guò)樣本小波變換估計(jì)得到響應(yīng)的演變功率譜密

3結(jié)論

本文提出了一種求解隨機(jī)B-T方程在完全非平穩(wěn)激勵(lì)作用下隨機(jī)響應(yīng)的方法。該方法的核心在于求解B-T方程的脈沖響應(yīng)函數(shù)。通過(guò)以下步驟得到：將運(yùn)動(dòng)方程改寫為狀態(tài)空間形式并進(jìn)行特征值分析；引人廣義坐標(biāo)并利用Laplace變換計(jì)算解耦后的半階線性微分方程；通過(guò)坐標(biāo)變換得到系統(tǒng)響應(yīng)，包括自由振動(dòng)和外激勵(lì)為階躍和脈沖函數(shù)時(shí)的解。得益于廣義坐標(biāo)半階微分方程的解析解，本文得到了3/2階B-T方程的解析解。通過(guò)時(shí)變頻響函數(shù)的概念，在隨機(jī)振動(dòng)頻域分析的框架內(nèi)得到了完全非平穩(wěn)激勵(lì)下隨機(jī)B-T方程的非平穩(wěn)響應(yīng)二階矩。MonteCarlo模擬驗(yàn)證了該方法的正確性。

本文針對(duì)3/2階B-T方程導(dǎo)出的確定性響應(yīng)解析解具有有限項(xiàng)和的簡(jiǎn)單形式。與Mittag-Leffler函數(shù)表示的復(fù)雜形式相比，本文給出的脈沖響應(yīng)函數(shù)解析解為進(jìn)一步的隨機(jī)振動(dòng)分析提供了基礎(chǔ)。本文所提方法不僅能應(yīng)用于3/2階系統(tǒng)，也能夠應(yīng)用于求解任意以2為分母的有理分?jǐn)?shù)階系統(tǒng)，如1/2或3/2階；本文所提方法有待進(jìn)一步拓展到非線性分?jǐn)?shù)階隨機(jī)動(dòng)力系統(tǒng)。

參考文獻(xiàn)：

[1] TORVIK PJ，BAGLEYRL.On theappearance of thefractional derivativeinthebehaviorofrealmaterials [J].Journal of Applied Mechanics，1984，51（2）： 294-298.

[2] DIPAOLAM，PIRROTTAA，VALENZAA.Visco-elasticbehaviorthrough fractionalcalculus：an easier method forbest fitting experimental results[J].MechanicsofMaterials，2011，43（12）：799-806.

[3]PODLUBNY I. Fractional Diferential Equations[M]. San Diego：Academic Press，1999.

[4]SRIVASTAVA H，JENA R M，CHAKRAVERTY S，et al.Dynamic response analysis of fractionally damped generalized Bagley-Torvik equation subject to external loads[J]. Russian Journal of Mathematical Physics，2020，27（2）：254-268.

[5]JENA R M，CHAKRAVERTY S.Analytical solution of Bagley-Torvik equations using Sumudu transformation method[J]． SN Applied Sciences，2019，1 （3）：246.

[6］黃瀟，芮偉國(guó).分?jǐn)?shù)階廣義Bagley-Torvik方程的各種 精確解及其動(dòng)力學(xué)性質(zhì)[J].云南大學(xué)學(xué)報(bào)（自然科學(xué) 版），2018，40（1）：12-21. HUANG Xiao，RUI Weiguo.Different kinds of exact solutions of fractional generalized Bagley-Torvik equation and their dynamical properties[J]. Journal of Yunnan University（Natural Sciences Edition），2018，40 （1）：12-21.

[7]TRINKS C，RUGE P. Treatment of dynamic systems with fractional derivatives without evaluating memoryintegrals[J]. Computational Mechanics，2002，29（6）： 471-476.

[8]WANG Z H，WANG X. General solution of the Bagley-Torvik equation with fractional-order derivative[J].Communications in Nonlinear Science and Numerical Simulation，2010，15（5）：1279-1285.

[9]RAJA MA Z，KHAN JA，QURESHI IM.Solution offractional order system of Bagley-Torvik equation us ingevolutionary computational intelligence[J].Mathematical Problems in Engineering，201，2011：675075.

[10］胡行華，秦艷杰.基于GA-Chebyshev神經(jīng)網(wǎng)絡(luò)的分?jǐn)?shù) 階Bagley-Torvik方程數(shù)值解法［J].計(jì)算數(shù)學(xué)，2023， 45（1）：109-129. HUXinghua，QIN Yanjie.Numerical solution of frac tional Bagley-Torvik equations based on GA-Cheby shev neural network[J]. Mathematica Numerica Sinica， 2023，45（1）：109-129.

[11]UDDIN M，AHMAD S.On the numerical solution of Bagley-Torvik equation via the Laplace transform[J]. Tbilisi Mathematical Journal，2017，10（1）：279-284.

[12] ALI ZAFAR A，KUDRA G，AWREJCEWICZ J. An investigation of fractional Bagley-Torvik equation[J]. Entropy，2020，22（1）： 28.

[13]CENESIZ Y，KESKINY，KURNAZ A. The solution of the Bagley-Torvik equation with the generalized Taylor collocation method[J]. Journal of the Franklin Institute，2010，347（2）：452-466.

[14] SINGH M P，CHANG T S，NANDAN H. Algorithms for seismic analysis of MDOF systems with frac tional derivatives[J].Engineering Structures，2O11，33 （8）：2371-2381.

[15]HUANG Z L，JIN X L. Response and stability of a SDOF strongly nonlinear stochastic system with light damping modeled by a fractional derivative[J]. Journal of Sound and Vibration，2009，319（3-5）：1121-1135.

［16］李書進(jìn)，張志聰，孔凡，等.周期與色噪聲聯(lián)合作用下 分?jǐn)?shù)階Duffing振子非平穩(wěn)響應(yīng)的無(wú)記憶方法[J].振 動(dòng)工程學(xué)報(bào)，2023，36（4）：923-933. LI Shujin，ZHANG Zhicong，KONG Fan，et al.A memory-free method for fractional-order Duffing systems subjected to combined periodic and colored excitation[J].Journal of Vibration Engineering，2023，36 （4）：923-933.

[17］孔凡，廖海君，韓仁杰，等.聯(lián)合激勵(lì)下分?jǐn)?shù)階非線性 系統(tǒng)非平穩(wěn)響應(yīng)的半解析方法[J].振動(dòng)工程學(xué)報(bào)， 2024，37（8）：1339-1348. KONG Fan，LIAO Haijun，HAN Renjie，et al.A semi-analytical method for non-stationary response de termination of nonlinear systems subjected to combined excitation[J]. Journal of Vibration Engineering，2024， 37（8）： 1339-1348.

[18] SUAREZ L，SHOKOOH A. An eigenvector expansion method for the solution of motion containing fractional derivatives[J]. Journal of Applied Mechanics， 1997，64（3）：629-635.

[19］孔凡，許伊鍵，韓仁杰，等. 1/2 分?jǐn)?shù)階線性隨機(jī)動(dòng)力 系統(tǒng)的非平穩(wěn)響應(yīng)解析解［J].動(dòng)力學(xué)與控制學(xué)報(bào)， 2023，21（4）：23-31. KONGFan，XU Yijian，HANRenjie，et al.Analytical solution for non-stationary response of 1/2 -orderfractional linear stochastic dynamical systems[J]. Journal of Dynamics and Control，2023，21（4）：23-31.

[20]DIPAOLAM，PINNOLAFP，SPANOSPD.Analysis of multi-degree-of-freedom systems with fractional derivative elements of rational order[C]//Proceedings of the ICFDA'14 International Conference on Fractional Differentiation and Its Applications.Catania，Italy， 2014：1-6.

[21]BAGLEY RL，CALICO R A.Fractional order state equations for the control of viscoelastically damped structures[J]. Journal of Guidance，Control，and Dy namics，1991，14（2）：304-311.

[22]PRIESTLEY M B. Evolutionary spectra and non-stationary processes[J]. Journal of the Royal Statistical SocietySeriesB：Statistical Methodology，1965，27（2）： 204-229.

[23]LABECCAW，GUIMARAESO，PIQUEIRAJRC. Analytical solution of general Bagley-Torvik equation [J].Mathematical Problems in Engineering，2015， 2015（1）：591715.

[24]MAHMUDOVNI，HUSEYNOVIT，ALIEVNA， et al.Analytical approach to a class of Bagley-Torvik equations[J].TWMS Journal of Pure and Applied Mathematics，2020，11（2）：238-258.

[25] ATANACKOVIC T M， ZORICA D. On the BagleyTorvik equation[J].Journal of Applied Mechanics， 2013，80（4）：0410113.

[26]RAY S S，BERA R K.Analytical solution of the Bagley Torvik equation by Adomian decomposition method[J].Applied Mathematics and Computation， 2005，168（1）：398-410.

［27］張德茂，袁曉，高小龍.基于Simulink電路模擬仿真求 解Bagley-Torvik方程[J].四川大學(xué)學(xué)報(bào)（自然科學(xué) 版），2019，56（2）：253-259. ZHANGDemao，YUANXiao，GAO Xiaolong.Solving Bagley-Torvik equation based on Simulink circuit simulation[J]. Journal of Sichuan University （Natural Science Edition），2019，56（2）：253-259.

[28] SPANOS P T，SOLOMOS G P. Markov approximation to transient vibration[J].Journal of Engineering Mechanics，1983，109（4）：1134-1150.

[29] SPANOS P D，F(xiàn)AILLA G. Evolutionary spectra estimation using wavelets[J]. Journal of Engineering Mechanics，2004，130（8）：952-960.

第一作者：孔凡（1984一），男，博士，教授。 E-mail：kongfan@hfut.edu.cn

通信作者：洪旭（1993一），男，博士，講師。 E-mail：xhong@hfut.edu.cn