孔凡 侯召旭 徐軍 李書(shū)進(jìn)

摘要： 結(jié)構(gòu)隔震和耗能減振裝置有時(shí)會(huì)同時(shí)體現(xiàn)出動(dòng)力特性的頻率依賴和滯回關(guān)系，分?jǐn)?shù)階導(dǎo)數(shù)模型能以較少的參數(shù)模擬力?位移關(guān)系的頻率依賴性。采用諧波平衡法研究簡(jiǎn)諧激勵(lì)下同時(shí)具有滯回特性和分?jǐn)?shù)階阻尼單元的系統(tǒng)穩(wěn)態(tài)響應(yīng)。利用激勵(lì)和響應(yīng)過(guò)程的Fourier級(jí)數(shù)展開(kāi)并取諧波平衡后，求得滯回響應(yīng)和位移響應(yīng)Fourier級(jí)數(shù)之間的關(guān)系;將滯回微分方程寫(xiě)成余量的形式并結(jié)合Galerkin方法和Levenberg?Marquardt算法求得響應(yīng)的Fourier級(jí)數(shù);分別對(duì)硬化和軟化Bouc?Wen滯回系統(tǒng)的數(shù)值模擬顯示所建議方法的精度。研究結(jié)果表明，分?jǐn)?shù)階數(shù)和穩(wěn)態(tài)位移幅值之間的關(guān)系依賴于系統(tǒng)和滯回參數(shù)以及激勵(lì)頻率。

關(guān)鍵詞： 多諧波平衡; Bouc?Wen滯回; Levenberg?Marquardt算法; 分?jǐn)?shù)階導(dǎo)數(shù)

中圖分類(lèi)號(hào)： TU311.3; O322? ? 文獻(xiàn)標(biāo)志碼： A? ? 文章編號(hào)： 1004-4523（2021）03-0552-07

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

引? 言

極端災(zāi)害性荷載作用下的工程結(jié)構(gòu)構(gòu)件和節(jié)點(diǎn)常表現(xiàn)出明顯的力?位移滯回關(guān)系，該關(guān)系表明恢復(fù)力不僅依賴于運(yùn)動(dòng)當(dāng)前狀態(tài)，而且依賴于其復(fù)雜的歷史狀態(tài)[1]。實(shí)際上，人們通常利用合理的結(jié)構(gòu)設(shè)計(jì)，期望在大震作用下結(jié)構(gòu)的特殊部位（如框架結(jié)構(gòu)的梁端）首先出現(xiàn)塑性鉸從而耗散地震輸入能量，從而達(dá)到結(jié)構(gòu)大震不倒的目的[1]。此外，一些用于結(jié)構(gòu)隔震和耗能減振的控制裝置，也通常表現(xiàn)出明顯滯回特性，如約束屈曲支撐、鉛芯橡膠隔震墊和黏滯阻尼器等裝置[2]。為了準(zhǔn)確合理描述結(jié)構(gòu)構(gòu)件和控制裝置中出現(xiàn)的滯回特性，人們發(fā)展了多種滯回模型，包括Ramberg?Osgood模型、多線性模型以及Bouc?Wen模型等[3]。Bouc?Wen模型具有接近實(shí)際情況的光滑曲線、僅用一個(gè)微分方程就可表現(xiàn)剛度退化和卸載拐點(diǎn)等復(fù)雜滯回特性，以及與運(yùn)動(dòng)微分方程聯(lián)合后形成易于求解的擴(kuò)階運(yùn)動(dòng)方程等優(yōu)點(diǎn)。

另外，許多應(yīng)用于工程結(jié)構(gòu)或控制裝置中的材料都存在本構(gòu)關(guān)系的頻率相關(guān)效應(yīng)，即本構(gòu)關(guān)系依賴于加載頻率[4]。傳統(tǒng)上，人們使用簡(jiǎn)單的Kelvin或Maxwell以及它們復(fù)雜的串并聯(lián)形式[5]（標(biāo)準(zhǔn)線性固體模型，Standard Linear Solid Models， SLS）部分描述這種頻率依賴特性。注意到，SLS模型需要大量參數(shù)才能精確地同時(shí)擬合試驗(yàn)得到的儲(chǔ)能和耗能模量數(shù)據(jù)。在此過(guò)程中，SLS模型的某些擬合參數(shù)易出現(xiàn)負(fù)值，從而失去物理意義[6]。反之，分?jǐn)?shù)階導(dǎo)數(shù)[7]模型以少量參數(shù)就能精確地描述本構(gòu)關(guān)系的頻率依賴性[6]，引起了越來(lái)越多的關(guān)注。該模型的控制裝置應(yīng)用包括：黏彈性阻尼[8]、隔震裝置[9]、液體黏滯阻尼器[10]、調(diào)諧液柱阻尼[11]，調(diào)諧質(zhì)量阻尼器[12]等。

可見(jiàn)，工程結(jié)構(gòu)、控制裝置和受控結(jié)構(gòu)可能會(huì)同時(shí)出現(xiàn)力?位移關(guān)系的滯回和頻率依賴特性。以一種新型的廣泛應(yīng)用于半主動(dòng)控制裝置的磁敏彈性體為例，本質(zhì)上作為一種橡膠黏彈性材料，具有明顯的本構(gòu)頻率依賴性[13]。為此，Behrooz等[14]利用三參數(shù)SLS模型和Bouc?Wen滯回模型的并聯(lián)形式、Eem等[15]采用Maxwell模型和Ramberg?Osgood滯回模型的并聯(lián)形式模擬磁敏彈性體的力?位移關(guān)系。Behrooz和Eem的研究均關(guān)注裝置的滯回力學(xué)性能，對(duì)黏彈性材料的頻率依賴性卻只以簡(jiǎn)單的SLS模型模擬。相反，在分?jǐn)?shù)階導(dǎo)數(shù)非線性動(dòng)力學(xué)研究方面，人們往往較多關(guān)注分?jǐn)?shù)階導(dǎo)數(shù)線性[16?17]或無(wú)記憶非線性[18?19]動(dòng)力系統(tǒng)，對(duì)具有滯回特性的分?jǐn)?shù)階系統(tǒng)研究較少。文獻(xiàn)調(diào)查顯示只有少量研究[20?21]關(guān)注滯回分?jǐn)?shù)動(dòng)力系統(tǒng)的隨機(jī)和穩(wěn)態(tài)動(dòng)力響應(yīng)。其中，文獻(xiàn)[20]通過(guò)等效線性化方法將滯回分?jǐn)?shù)階系統(tǒng)化為等效線性參數(shù)依賴于幅值的線性系統(tǒng)，得到了滯回分?jǐn)?shù)階系統(tǒng)頻率?響應(yīng)幅值曲線。該曲線只具有主共振峰值，而無(wú)法反應(yīng)超諧共振峰值。

考察動(dòng)力系統(tǒng)在諧波激勵(lì)下的穩(wěn)態(tài)響應(yīng)是理解系統(tǒng)動(dòng)力特性重要的一環(huán)，因此，本文與文獻(xiàn)[20]一樣，關(guān)注諧波激勵(lì)下滯回分?jǐn)?shù)階系統(tǒng)穩(wěn)態(tài)響應(yīng)。但是，本文與文獻(xiàn)[20]采用的方法不同，本文采用諧波平衡法研究具有不同分?jǐn)?shù)階導(dǎo)數(shù)的滯回系統(tǒng)的穩(wěn)態(tài)響應(yīng)，以期得到頻率?響應(yīng)幅值的完整信息。首先，發(fā)展分?jǐn)?shù)階滯回系統(tǒng)的諧波平衡法，將滯回運(yùn)動(dòng)微分轉(zhuǎn)化為一組代數(shù)方程;隨后，利用Levenberg?Marquardt算法對(duì)與該代數(shù)方程對(duì)應(yīng)的最小二乘問(wèn)題求解;最后，考察硬化和軟化的Bouc?Wen滯回系統(tǒng)，得到它們?cè)诓煌謹(jǐn)?shù)階數(shù)情況下的頻率?響應(yīng)幅值曲線。與時(shí)域逐步積分法的對(duì)比顯示了所建議方法的精度。

1 分?jǐn)?shù)階滯回系統(tǒng)的諧波平衡法

考慮具有分?jǐn)?shù)階阻尼項(xiàng)的單自由度滯回動(dòng)力系統(tǒng)

式中? 分別為質(zhì)量和線彈性剛度系數(shù);分別為系統(tǒng)位移和加速度響應(yīng);為滯回位移;為屈服前后剛度比;為分?jǐn)?shù)階阻尼系數(shù);為階Riemann?Liouville意義上的分?jǐn)?shù)階導(dǎo)數(shù)，即

下文為方便計(jì)，將Caputo分?jǐn)?shù)階導(dǎo)數(shù)簡(jiǎn)寫(xiě)為并簡(jiǎn)稱(chēng)為分?jǐn)?shù)階導(dǎo)數(shù)。

另一方面，將位移響應(yīng)和激勵(lì)在區(qū)間上進(jìn)行周期延拓后，可將它們寫(xiě)為Fourier級(jí)數(shù)的形式：

式（18）表明了滯回位移的Fourier級(jí)數(shù)向量可表示為位移的Fourier級(jí)數(shù)向量的形式。然而，這兩個(gè)級(jí)數(shù)向量都未知，需引進(jìn)它們之間的第二個(gè)關(guān)系式。為此，引入和之間的第二個(gè)關(guān)系式，即形如

利用Galerkin法將余量式（20）在Fourier基上投影，可得到第二個(gè)與之間的關(guān)系式，隨后將此關(guān)系式代入式（18）后可將向量消除，即可得到以為未知量的方程

式?中? 為余量的Fouier級(jí)數(shù)。求解式（21）可得位移響應(yīng)的Fourier級(jí)數(shù)，對(duì)響應(yīng)的Fourier級(jí)數(shù)作Fourier逆變換后可得位移響應(yīng)時(shí)程。

2 代數(shù)方程數(shù)值求解的Levenberg-Marquardt算法

式（21）的求解可通過(guò)許多數(shù)值算法實(shí)現(xiàn)，本文采用一種高效的Levenberg?Marquardt（LM）算法。LM算法最早由Marquardt[23]于1963年提出，是一種應(yīng)用廣泛的非線性最小二乘算法，可視為最速下降法和線性化方法（泰勒級(jí)數(shù)）的一種組合;可參閱相關(guān)數(shù)值計(jì)算參考書(shū)目。它將代數(shù)方程轉(zhuǎn)化為求解最小二乘的形式，即

式中? 均為控制滯回環(huán)形狀的模型參數(shù)。

3 數(shù)值算例

利用與逐步積分法得到的相關(guān)結(jié)果對(duì)比驗(yàn)證所建議方法的適用性和精度。作為演示，本算例采用具有微分形式的Bouc?Wen滯回模型，但本文發(fā)展的方法同時(shí)也適用其他滯回關(guān)系。注意到，具有Bouc?Wen滯回的整數(shù)階動(dòng)力系統(tǒng)響應(yīng)的數(shù)值解可通過(guò)增量形式的Newmark???Newton???Raphson （NNR）迭代法實(shí)現(xiàn)，其本質(zhì)是將動(dòng)力方程轉(zhuǎn)化為增量靜力方程，難點(diǎn)在于切線性變剛度的求解和滯回模型拐點(diǎn)的處理。在此基礎(chǔ)上，為求解Bouc?Wen滯回的分?jǐn)?shù)階動(dòng)力系統(tǒng)響應(yīng)，作者發(fā)展了一種基于Riemann?Liouville離散分?jǐn)?shù)階導(dǎo)數(shù)的NNR方法：其中，與Newmark?1/6法處理加速度項(xiàng)對(duì)應(yīng)，在對(duì)離散式（5）右邊積分時(shí)，須假定離散區(qū)間內(nèi)的速度呈線性變化。限于篇幅和行文重點(diǎn)，該增量形式的數(shù)值方法茲不贅述，擬另文發(fā)表。

首先，考察所建議方法在特定參數(shù)設(shè)置下計(jì)算得到的動(dòng)力響應(yīng)時(shí)程和滯回環(huán)。選取動(dòng)力系統(tǒng)參數(shù)，;滯回參數(shù)，;激勵(lì)參數(shù)，，。系統(tǒng)質(zhì)量和剛度參數(shù)選取時(shí)體現(xiàn)了振子的“歸一化”;阻尼參數(shù)選取時(shí)體現(xiàn)了工程系統(tǒng)的小阻尼特性;滯回參數(shù)選取時(shí)體現(xiàn)了明顯的滯回效應(yīng)。圖1（a）?（c）所示為在上述參數(shù)設(shè)置情況下所建議方法和逐步積分法得到的，和滯回環(huán)對(duì)比。由圖1（a）和（b）可知：除了位移時(shí)程的開(kāi)始階段，多諧波平衡法得到的平穩(wěn)響應(yīng)能很好地吻合逐步積分法得到的結(jié)果;時(shí)程開(kāi)始的誤差由多諧波平衡法忽略初始條件導(dǎo)致。

隨后，考察分?jǐn)?shù)階數(shù)對(duì)軟化Bouc?Wen系統(tǒng)響應(yīng)幅值的影響。除分?jǐn)?shù)階數(shù)外，其他模型及激勵(lì)參數(shù)與上例相同。圖2所示為多諧波平衡法和逐步積分法計(jì)算得到的響應(yīng)和幅值與分?jǐn)?shù)階數(shù)之間的關(guān)系曲線。由圖2可見(jiàn)：所建議方法與時(shí)域逐步積分法吻合較好;的幅值隨著的增大而減小，的幅值隨著的增大稍有減小但幾乎保持不變。圖3給出了時(shí)的滯回環(huán)。由圖3可知，滯回環(huán)的飽滿程度隨著的增大而減小。

接著，考察分?jǐn)?shù)階數(shù)對(duì)硬化Bouc?Wen系統(tǒng)響應(yīng)幅值的影響。此時(shí)，取滯回參數(shù)，，其他系統(tǒng)、滯回和激勵(lì)參數(shù)與上例同。圖4所示為多諧波平衡法和逐步積分法計(jì)算得到的響應(yīng)（和）幅值與分?jǐn)?shù)階數(shù)之間的關(guān)系曲線。由圖可知，在該情況下，和幅值均隨著的增大而增大，前者增大幅值更加明顯。圖5所示為時(shí)系統(tǒng)的滯回環(huán)?？梢?jiàn)，所有滯回環(huán)均呈現(xiàn)明顯的硬化特性，分?jǐn)?shù)階數(shù)對(duì)硬化滯回環(huán)的影響較軟化滯回環(huán)時(shí)小。

分?jǐn)?shù)階數(shù)對(duì)硬化和軟化系統(tǒng)響應(yīng)幅值的影響不同。為深入考察這一點(diǎn)，本文利用諧波平衡法和逐步積分法得到了分?jǐn)?shù)階導(dǎo)數(shù)的情況下，上述軟化和硬化系統(tǒng)的激勵(lì)頻率?響應(yīng)幅值曲線，如圖6和7所示。由圖6?7可見(jiàn)，在所有情況下，諧波平衡法均能很好地吻合逐步積分法給出的結(jié)果。此外，應(yīng)注意到，文[20]利用線性化方法研究滯回分?jǐn)?shù)階系統(tǒng)的頻率響應(yīng)曲線，只能得到主共振峰值，無(wú)法得到超諧共振峰值;而本文采用的諧波平衡法能同時(shí)得到主共振和超諧共振峰值。文[24]采用增量諧波平衡法和平均法研究了分?jǐn)?shù)階硬Duffing振子的頻率響應(yīng)曲線，也得到類(lèi)似的結(jié)論。由圖6?7可見(jiàn)：軟特性Bouc?Wen分?jǐn)?shù)階系統(tǒng)的激勵(lì)頻率?響應(yīng)幅值曲線左側(cè)較陡;反之，硬特性Bouc?Wen分?jǐn)?shù)階系統(tǒng)的激勵(lì)頻率?響應(yīng)幅值曲線右側(cè)較陡。本例采用的激勵(lì)頻率（）處于圖6所示軟特性峰值右側(cè)，分?jǐn)?shù)階數(shù)增大時(shí)響應(yīng)幅值呈明顯遞減趨勢(shì);同時(shí)，激勵(lì)頻率處于圖7所示硬特性峰值左側(cè)，分?jǐn)?shù)階數(shù)增大時(shí)響應(yīng)幅值呈小幅增加趨勢(shì)。由以上分析不難得出：對(duì)于特定滯回系統(tǒng)，分?jǐn)?shù)階數(shù)對(duì)響應(yīng)幅值的影響與激勵(lì)頻率大小有密切關(guān)系。例如，對(duì)于如圖6所示的軟特性分?jǐn)?shù)階系統(tǒng)，與時(shí)的情況不同，當(dāng)激勵(lì)頻率時(shí)，響應(yīng)幅值隨分?jǐn)?shù)階數(shù)的增加而增加。另外，分?jǐn)?shù)階數(shù)對(duì)響應(yīng)幅值的影響還與滯回參數(shù)有關(guān)，如本例所示軟或硬特性系統(tǒng)參數(shù)。值得注意的是，當(dāng)系統(tǒng)阻尼和滯回耗能較小時(shí)，非線性系統(tǒng)的激勵(lì)頻率?響應(yīng)幅值曲線易出現(xiàn)分岔現(xiàn)象，進(jìn)而導(dǎo)致更加復(fù)雜的分?jǐn)?shù)階數(shù)和響應(yīng)幅值之間的關(guān)系。

現(xiàn)對(duì)所提出方法在參數(shù)設(shè)置和計(jì)算效率方面作如下評(píng)述。首先，在計(jì)算式（24）所示的Jacobian矩陣時(shí)，每個(gè)迭代步均涉及時(shí)域和頻域之間的相互轉(zhuǎn)化，必須進(jìn)行Fourier變換[25]。因此，使用快速Fourier變換（Fast Fourier Transform， FFT）及其逆變換能加快速度。其次，系統(tǒng)非線性的存在，可能會(huì)產(chǎn)生高頻響應(yīng);因此，在設(shè)置響應(yīng)上限截止頻率時(shí)，須覆蓋響應(yīng)可能的頻率區(qū)間。本文采用的上限截止頻率大于10倍激勵(lì)頻率。最后，頻率泄露可能影響收斂速度并降低計(jì)算精度。為此，在設(shè)置采樣頻率時(shí)，建議在激勵(lì)頻率上取離散頻率點(diǎn)，即激勵(lì)頻率為離散頻率步長(zhǎng)的整數(shù)倍。計(jì)算每個(gè)激勵(lì)頻率點(diǎn)的響應(yīng)幅值時(shí)，均須綜合考慮上述要求，合理指定計(jì)算參數(shù)。例如，采用與圖1相同的系統(tǒng)參數(shù)設(shè)置，選取上限截止頻率，頻域內(nèi)采樣間隔，共100個(gè)采樣點(diǎn)。此時(shí)，第10個(gè)頻率離散點(diǎn)等于激勵(lì)頻率。試算表明，僅5次迭代就可達(dá)到預(yù)定收斂目標(biāo)，耗時(shí)0.07 s。

本文結(jié)合Riemman?Liouville定義的離散格式和Newmark?Newton?Raphson迭代方法在時(shí)域中求解Bouc?Wen分?jǐn)?shù)階導(dǎo)數(shù)系統(tǒng)的動(dòng)力響應(yīng)。Riemman?Liouville離散格式中，當(dāng)前積分步的運(yùn)動(dòng)狀態(tài)須表示為之前所有離散時(shí)間點(diǎn)運(yùn)動(dòng)狀態(tài)的加權(quán)和。其中，每個(gè)歷史時(shí)間點(diǎn)的權(quán)重須及時(shí)更新，大大降低了計(jì)算效率。此外，與整數(shù)階非線性動(dòng)力系統(tǒng)的Newton?Raphson迭代方法一樣，更新切線剛度與尋找滯回位移拐點(diǎn)也增加了時(shí)域方法的計(jì)算量。計(jì)算經(jīng)驗(yàn)表明：時(shí)域方法的計(jì)算時(shí)長(zhǎng)比時(shí)間離散點(diǎn)數(shù)增長(zhǎng)更快;同時(shí)，須采用非常小的時(shí)間步長(zhǎng)以保證計(jì)算精度。例如，采用與上述頻域方法相同的系統(tǒng)參數(shù)設(shè)置，使時(shí)域方法達(dá)到與頻域方法相同的計(jì)算精度時(shí)，計(jì)算耗時(shí)0.71 s。

4 結(jié)論與展望

本文發(fā)展了諧和激勵(lì)下滯回分?jǐn)?shù)階動(dòng)力系統(tǒng)平穩(wěn)解的諧波平衡法。首先，利用諧波平衡法將滯回分?jǐn)?shù)階微分方程轉(zhuǎn)化為了一組以響應(yīng)Fourier級(jí)數(shù)為未知量的非線性代數(shù)方程;利用Levenberg?Marquardt算法求解了與非線性代數(shù)方程對(duì)應(yīng)的最小二乘問(wèn)題，得到了響應(yīng)Fourier級(jí)數(shù)。研究表明，諧波平衡法得到的響應(yīng)時(shí)程、滯回環(huán)、頻率?幅值曲線均與逐步積分法得到的相關(guān)結(jié)果吻合較好。該方法能同時(shí)得到滯回分?jǐn)?shù)階系統(tǒng)頻率?幅值曲線的主共振和超諧共振峰值。

對(duì)具有不同滯回特性的系統(tǒng)，考察了分?jǐn)?shù)階數(shù)對(duì)響應(yīng)幅值大小的影響。文中給出了在特定軟特性和硬特性參數(shù)的情況下，激勵(lì)頻率等于線性系統(tǒng)自振頻率時(shí)，分?jǐn)?shù)階數(shù)對(duì)響應(yīng)幅值影響的規(guī)律。研究發(fā)現(xiàn)，分?jǐn)?shù)階數(shù)對(duì)響應(yīng)幅值影響較復(fù)雜，與系統(tǒng)特性、滯回特性、激勵(lì)頻率均有一定關(guān)系;擬另文發(fā)表。

進(jìn)一步的工作可在其他滯回模型、多自由度系統(tǒng)以及滯回分?jǐn)?shù)階系統(tǒng)的混沌和分岔等方面展開(kāi)。

參考文獻(xiàn)：

[1] 李 杰，李國(guó)強(qiáng).地震工程學(xué)導(dǎo)論[M].北京：地震出版社，1992.

[2] 歐進(jìn)萍.結(jié)構(gòu)振動(dòng)控制[M].北京：科學(xué)出版社，2003.

[3] 潘鵬，張耀庭.建筑抗震設(shè)計(jì)理論與方法[M].北京：科學(xué)出版社，2017.

[4] Xu Z-D， Wang D-X， Shi C-F. Model， tests and application design for viscoelastic dampers[J]. Journal of Vibration and Control， 2011， 17： 1359-1370.

[5] Xu Z-D. Earthquake mitigation study on viscoelastic dampers for reinforced concrete structures[J]. Journal of Vibration and Control， 2007， 13： 29-43.

[6] Di Paola M， Pirrotta A， Valenza A. Visco-elastic behavior through fractional calculus： An easier method for best fitting experimental results[J]. Mechanics of Materials， 2011， 43（12）： 799-806.

[7] Oldham K， Spanier J. The Fractional Calculus Theory and Applications of Differentiation and Integration to Arbitrary Order[M]. Elsevier， 1974.

[8] Xu Z-D， Xu C， Hu J. Equivalent fractional Kelvin model and experimental study on viscoelastic damper[J]. Journal of Vibration and Control， 2015， 21（13）： 2536-2552.

[9] Hwang J， Wang J. Seismic response prediction of HDR bearings using fractional derivative Maxwell model[J]. Engineering Structures， 1998， 20： 849-856.

[10] Makris N， Dargush G， Constantinou M. Dynamic analysis of generalized viscoelastic fluids[J]. Journal of Engineering Mechanics， 1993， 119： 1663-1679.

[11] Di Matteo A， Iacono F L， Navarra G， et al. Innovative modeling of tuned liquid column damper motion?[J]. Communications in Nonlinear Science and Numerical Simulation， 2015， 23： 229-244.

[12] Rüdinger F. Tuned mass damper with fractional derivative damping?[J]. Engineering Structures， 2006， 28： 1774-1779.

[13] Li Y， Li J. On rate-dependent mechanical model for adaptive magnetorheological elastomer base isolator [J]. Smart Materials and Structures， 2017， 26： 045001.

[14] Behrooz M， Wang X， Gordaninejad F. Modeling of a new semi-active/passive magnetorheological elastomer isolator[J]. Smart Materials and Structures， 2014， 23： 045013.

[15] Eem S H， Jung H J， Koo J H. Modeling of magneto-rheological elastomers for harmonic shear deformation[J]. IEEE Transactions on Magnetics， 2012， 48： 3080-3083.

[16] Chang T S， Singh M P. Seismic analysis of structures with a fractional derivative model of viscoelastic dampers[J]. Earthquake Engineering and Engineering Vibration， 2002， 1： 251-260.

[17] Singh M P， Chang T S， Nandan H. Algorithms for seismic analysis of MDOF systems with fractional derivatives?[J]. Engineering Structures， 2011， 33： 2371-2381.

[18] Failla G， Pirrotta A. On the stochastic response of a fractionally-damped Duffing oscillator[J]. Communications in Nonlinear Science and Numerical Simulation， 2012， 17： 5131-5142.

[19] Chen L， Wang W， Li Z， et al. Stationary response of Duffing oscillator with hardening stiffness and fractional derivative[J]. International Journal of Non-Linear Mechanics， 2013， 48： 44-50.

[20] Spanos P D， Di Matteo A， Pirrotta A. Steady-state dynamic response of various hysteretic systems endowed with fractional derivative elements[J]. Nonlinear Dynamics， 2019，98（4）：3113-3124.

[21] Di Matteo A， Spanos P D， Pirrotta A. Approximate survival probability determination of hysteretic systems with fractional derivative elements[J]. Probabilistic Engineering Mechanics， 2018， 54： 138-146.

[22] Tren?evski K， Tomovski ?. On fractional derivatives of some functions of exponential type[J]. Publikacije Elektrotehni?kog fakulteta-serija： Matematika， 2002，（13）： 77-84.

[23] Marquardt D W. An algorithm for least-squares estimation of nonlinear parameters[J]. Journal of the Society for Industrial and Applied Mathematics， 1963， 11（2）： 431-441.

[24] Shen Y J， Wen S F， Li X H， et al. Dynamical analysis of fractional-order nonlinear oscillator by incremental harmonic balance method[J]. Nonlinear Dynamics， 2016， 85： 1457-1467.

[25] Wong C W， Ni Y Q， Lau S L. Steady-state oscillation of hysteretic differential model. I： Response analysis[J]. Journal of Engineering Mechanics， 1994， 120（11）： 2271-2298.

作者簡(jiǎn)介： 孔? 凡（1984-），男，博士，副教授，碩士生導(dǎo)師。E-mail： kongfan@whut.edu.cn