孔 凡，韓仁杰，張遠(yuǎn)進(jìn)

（1.武漢理工大學(xué)土木工程與建筑學(xué)院，湖北武漢 430070；2.合肥工業(yè)大學(xué)土木與水利工程學(xué)院，安徽合肥 230009；3.武漢理工大學(xué)安全科學(xué)與應(yīng)急管理學(xué)院，湖北武漢 430070）

引言

工程動力作用具有強烈的隨機性，通常利用非平穩(wěn)隨機過程描述［1］。以完全非平穩(wěn)地震動為例：它不僅在強度上具有明顯的上升和衰減過程，而且在不同時間段上能量的頻率分布也不同［2］。因此，地震動非平穩(wěn)性體現(xiàn)為幅值非平穩(wěn)和頻率非平穩(wěn)。學(xué)者們提出了各種模型以全面反映地震動的非平穩(wěn)特性［3-6］。近年來，小波變換以其時-頻聯(lián)合分辨特性在地震動建模方面得到了廣泛發(fā)展，詳見文獻(xiàn)［7-10］。

復(fù)雜非線性系統(tǒng)在隨機激勵下的隨機動力響應(yīng)是隨機振動研究的重要課題。具體而言，利用Fokker-Planck-Kolmogorov（FPK）方程［11］、隨機平均法［12］、統(tǒng)計線性化（或非線性化）［13］、矩截斷［14］、信息理論（熵）方法［15］以及最近發(fā)展的概率密度演化方法［16］等，人們已解析或數(shù)值地得到了很多非線性隨機動力系統(tǒng)的精確解或近似解。然而，高效地計算非平穩(wěn)激勵下復(fù)雜非線性系統(tǒng)的非平穩(wěn)響應(yīng)仍然是工程隨機振動理論發(fā)展面臨的巨大挑戰(zhàn)［17］，具有重要的實際意義。統(tǒng)計線性化方法是隨機振動理論中具有廣泛適用性的非線性方法之一。它可以用于不同的無記憶非線性［18］，甚至滯回非線性系統(tǒng)［19］、平穩(wěn)或非平穩(wěn)隨機激勵［17］、多自由度［20］或連續(xù)系統(tǒng)［21］，最近的文獻(xiàn)綜述詳見文獻(xiàn)［22］。

上述研究僅考慮激勵為隨機過程的情況。然而，工程實際中，結(jié)構(gòu)往往會同時受到確定性周期和隨機激勵的聯(lián)合作用。例如，旋轉(zhuǎn)式飛機［23］經(jīng)常受到色噪聲和諧波激勵聯(lián)合作用；俘能系統(tǒng)的非線性響應(yīng)［24］；近斷層脈沖地震作用下的結(jié)構(gòu)等。為此，人們發(fā)展了聯(lián)合激勵下非線性系統(tǒng)隨機動力響應(yīng)的系列方法：Chang［25］使用高斯與非高斯線性化求得了Duffing 系統(tǒng)在諧波與高斯白噪聲激勵聯(lián)合作用下的平穩(wěn)響應(yīng)；Chen 等［26］給出了分?jǐn)?shù)階Duffing 系統(tǒng)在確定性激勵與平穩(wěn)白噪聲聯(lián)合作用下平穩(wěn)響應(yīng)的隨機平均法；Zhu 等［27］利用諧波平衡和矩截斷方法研究了Duffing 系統(tǒng)在諧波與隨機激勵聯(lián)合作用下的平穩(wěn)響應(yīng)；Zhang 等［28］給出了周期與隨機激勵聯(lián)合作用下多自由度齒輪結(jié)構(gòu)平穩(wěn)響應(yīng)的統(tǒng)計線性化方法。然而，考察聯(lián)合激勵下非平穩(wěn)隨機作用帶來的響應(yīng)非平穩(wěn)性，尚未引起人們足夠的重視。這種隨機響應(yīng)的非平穩(wěn)性主要來自于三個方面：突加激勵非平穩(wěn)、與諧波響應(yīng)耦合的非平穩(wěn)、隨機激勵自身的變慢非平穩(wěn)。

注意到統(tǒng)計線性化方法的廣泛適用性，本文提出了一種用于求解確定性周期與非平穩(wěn)隨機激勵聯(lián)合作用下、單自由度非線性系統(tǒng)非平穩(wěn)響應(yīng)的統(tǒng)計線性化方法，可視為作者發(fā)展的平穩(wěn)方法的非平穩(wěn)推廣［29-30］。該方法將系統(tǒng)響應(yīng)分解為確定性諧波和零均值隨機分量之和，將原非線性運動方程等效地化為一組耦合的、分別以確定性和隨機動力響應(yīng)為未知量的非線性微分方程。然后，利用統(tǒng)計線性化方法將非平穩(wěn)隨機激勵作用下的非線性隨機動力方程化為等效線性方程，得到關(guān)于線性隨機響應(yīng)二階矩的李雅普諾夫微分方程。最后，聯(lián)立李雅普諾夫微分方程與諧波激勵作用下的確定性微分方程，通過數(shù)值算法（如龍格-庫塔法）對這組耦合非線性常微分方程進(jìn)行求解。數(shù)值算例中，利用蒙特卡洛模擬驗證此方法的適用性和精度。

1 動力學(xué)方程

單自由度非線性系統(tǒng)在確定性周期和非平穩(wěn)隨機激勵聯(lián)合作用下的運動方程為：

式中m，c和k分別為質(zhì)量、阻尼和剛度系數(shù)；和分別為響應(yīng)的位移、速度和加速度；為和的非線性函數(shù)；Q(t)為零均值均勻調(diào)制非平穩(wěn)隨機激勵，可寫為：

式中a(t)為調(diào)制函數(shù)；Qs(t)為零均值平穩(wěn)隨機過程；F(t)為確定性周期激勵，寫為Fourier 展開形式為：

式中Fl為第l次諧波幅值。

激勵的非平穩(wěn)性導(dǎo)致響應(yīng)的非平穩(wěn)性。這里，響應(yīng)的非平穩(wěn)性主要來自于：突加激勵、諧波激勵諧和變化以及隨機激勵幅值慢變非平穩(wěn)。假設(shè)響應(yīng)x(t)可分解為均值μx(t)和方差為的零均值非平穩(wěn)隨機過程之和，即：

同時，

式中 D[·]表示方差。將式（4）代入式（1）中得：

為方便計，省略了響應(yīng)量的時間參數(shù)。對式（6）兩邊求數(shù)學(xué)期望得：

式中 E[·]表示數(shù)學(xué)期望，用式（6）減式（7）得：

式中

2 統(tǒng)計線性化方法

利用統(tǒng)計線性化方法將式（8）化為線性方程：

式中ce(t)，ke(t)為等效線性參數(shù)。當(dāng)隨機激勵Q(t)為高斯過程時，響應(yīng)x^ (t)也為高斯過程；即使Q(t)為非高斯過程，當(dāng)系統(tǒng)阻尼較小時，響應(yīng)也可近似為高斯過程。因此，ce(t)，ke(t)均為時間t的函數(shù)，且有：

以下考慮Qs(t)為白噪聲和色噪聲兩種情況。

2.1 Qs(t)為白噪聲

此時，Q(t)為均勻調(diào)制白噪聲，即：

式中w(t)為零均值高斯白噪聲，功率譜密度為S0。將式（10）化為狀態(tài)方程的形式：

與式（13）對應(yīng)的李雅普諾夫方程為：

式中v為響應(yīng)的協(xié)方差矩陣，即：

vij(t)為響應(yīng)(qi，qj)，i，j=1，2 的方差/協(xié)方差，且vij(t)=vji(t)；Θ(t)為激勵Q(t)的協(xié)方差矩陣，即：

結(jié)合響應(yīng)的高斯性可知，期望：

可得

所有未知量的初始值均為零，可利用龍格-庫塔法求解式（16）～（20），得響應(yīng)的確定性分量時程與隨機分量協(xié)方差矩陣。

2.2 Qs(t)為色噪聲

此時，Q(t)為均勻調(diào)制色噪聲，平穩(wěn)隨機過程Qs(t)的功率譜密度為S(ω)。將式（2）中的Qs(t)表示為白噪聲經(jīng)過成型濾波器的形式，即：

式中Φ為濾波器的前置輸出，n為成型濾波器維度；Φ(i)表示Φ的第i階導(dǎo)數(shù)；λ0，λ1，…，λn-1和ν0，ν1，…，νn-1為濾波器參數(shù)；w(t)為零均值高斯白噪聲，功率譜密度為S0。令：

則式（22）可以寫為：

式中

令：

則式（21）可寫為：

聯(lián)立式（8）、式（24）和式（26），得：

與式（27）對應(yīng)的李雅普諾夫方程為：

式中

且Vqq=E{qqT}，VqΦ=VΦq=E{qΦT}VΦΦ=E{ΦΦT}。

同樣地，將式（28）與式（7）聯(lián)立，可利用數(shù)值方法求解得到所有未知量需要注意的是，VΦΦ的值僅與色噪聲本身有關(guān)，可直接確定。具體的計算步驟將在算例中進(jìn)一步說明。

3 數(shù)值算例

作為算例，不失一般性地考慮具有立方非線性剛度的Duffing 振子，即：

式中ε為表示非線性強度的小量。式（7）中非線性項的數(shù)學(xué)期望為：

從而，式（7）可寫為：

同樣地，式（8）可化為：

依據(jù)式（11），等效線性參數(shù)為：

式（13）中，系數(shù)矩陣：

如果確定性激勵為單諧波的形式，即：

式中a0為確定性激勵幅值，ω0為確定性激勵頻率。

假定調(diào)制函數(shù)：

式中A為調(diào)制函數(shù)峰值，μ標(biāo)識調(diào)制函數(shù)下降速度。

3.1 Qs(t)為白噪聲

此時，Q(t)為均勻調(diào)制白噪聲，即：

將式（33）～（37）代入式（16）～（20）中，利用數(shù)值方法求得y，即所有確定性與隨機響應(yīng)分量。

3.1.1 典型參數(shù)設(shè)置的情況

選取系統(tǒng)參數(shù)m=1，c=0.4，k=1，ε=0.5；確定性激勵參數(shù)a0=1，ω0=1 rad/s；隨機激勵參數(shù)本文利用蒙特卡洛模擬（Monte Carlo simulation，MCS）驗證所建議方法的適用性。其中，利用譜表現(xiàn)方法生成白噪聲樣本10000 條，乘以調(diào)制函數(shù)得均勻調(diào)制白噪聲。利用本文所建議方法與MCS 得到了位移均值及其標(biāo)準(zhǔn)差，其對比如圖1（a），（b）所示。

圖1 Duffing 系統(tǒng)在諧波與調(diào)制白噪聲聯(lián)合作用下的位移Fig.1 Displacement of the considered Duffing system subjected to combined harmonic and modulated white noise

從圖1（a），（b）可以看出，在考慮的參數(shù)設(shè)置情況下，本文建議方法（Proposed Method，PM）與蒙特卡洛模擬的均值幾乎完全吻合。在隨機響應(yīng)分量標(biāo)準(zhǔn)差方面，所建議方法能很好地捕捉由于突加激勵（由0 上升的過程）、確定性諧和激勵分量（下降過程中呈諧和變化）和隨機激勵分量調(diào)制函數(shù)帶來的非平穩(wěn)性（慢變的上升-下降過程）。

3.1.2 非線性強度的影響

下面分析此方法在不同非線性強度下的適用性。其他參數(shù)與3.1.1 節(jié)選取相同，非線性強度系數(shù)ε取0～1。定義確定性諧和分量平均功率為：

式中T為響應(yīng)持時。定義隨機響應(yīng)分量標(biāo)準(zhǔn)差的時間平均為：采用所建議方法與蒙特卡洛模擬計算得到的確定性響應(yīng)分量平均功率對比如圖2所示，標(biāo)準(zhǔn)差時間平均對比如圖3所示。

圖2 非線性強度對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.2 Averaged power of the deterministic response component of the Duffing system with different degrees of nonlinearity

圖3 非線性強度對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.3 Time-averaged standard deviation of the stochastic response component of the Duffing system with different degrees of nonlinearity

從圖2，3 可見，在其他參數(shù)不變的情況下，確定性響應(yīng)平均功率隨著非線性強度增加而降低；位移標(biāo)準(zhǔn)差的時間平均隨非線性強度增加而降低，速度響應(yīng)標(biāo)準(zhǔn)差的時間平均隨非線性強度增加而略有減小。所有非線性強度情況下，兩種方法所得結(jié)果均吻合較好。

3.1.3 激勵幅值的影響

為研究確定性激勵幅值對此方法適用性的影響，保持其他參數(shù)不變，確定性激勵幅值a0取0～2 rad/s。利用本文所建議方法與蒙特卡洛模擬計算得到的確定性響應(yīng)分量平均功率對比如圖4所示；時變標(biāo)準(zhǔn)差的時間平均對比如圖5所示。

圖4 諧波激勵幅值對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.4 Averaged power of the deterministic response component versus the harmonic excitation amplitude

圖5 諧波激勵幅值對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.5 Influence of the harmonic excitation amplitude on the time-averaged standard deviation of the response

由圖4，5 可見，在其他參數(shù)不變的情況下，確定性響應(yīng)平均功率隨諧波激勵幅值增大而增大；位移標(biāo)準(zhǔn)差的時間平均隨諧波激勵幅值增大而減小，速度標(biāo)準(zhǔn)差的時間平均隨諧波激勵幅值增大而增大。不同激勵幅值下，兩種方法所得的結(jié)果均吻合較好。

3.1.4 確定性諧和激勵頻率的影響

為研究確定性激勵頻率對此方法適用性的影響，保持其他參數(shù)不變，確定性激勵頻率ω0取0.1～2 rad/s。采用本文所建議方法與蒙特卡洛模擬得到的確定性響應(yīng)分量的平均功率對比如圖6所示；隨機響應(yīng)非平穩(wěn)標(biāo)準(zhǔn)差的時間平均如圖7所示。

圖6 諧波激勵頻率對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.6 Influence of the harmonic excitation frequency on the averaged power of the deterministic response component

圖7 諧波激勵頻率對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.7 Influence of the harmonic excitation frequency on the time-averaged standard deviation of the response

可見，諧波激勵頻率對確定性響應(yīng)平均功率和隨機響應(yīng)標(biāo)準(zhǔn)差平均有較大影響，存在使他們達(dá)到極值的諧波激勵頻率。其中，使確定性響應(yīng)平均功率達(dá)到最大的諧波激勵頻率約為1.4 rad/s，而使隨機響應(yīng)標(biāo)準(zhǔn)差時間平均值達(dá)到最大的諧波激勵頻率約為1.5 rad/s。此外，不同激勵頻率下本文所建議的方法與蒙特卡洛模擬所得結(jié)果均吻合較好。

3.1.5 隨機激勵強度的影響

為研究隨機激勵強度對此方法適用性的影響，保持其他參數(shù)不變，功率譜強度S0取0～1。采用本文所建議方法與蒙特卡洛模擬得到的確定性響應(yīng)分量平均功率對比如圖8所示；隨機響應(yīng)時變標(biāo)準(zhǔn)差的時間平均對比如圖9所示。

圖8 隨機激勵功率譜強度對Duffing 系統(tǒng)確定性響應(yīng)功率的影響Fig.8 Influence of the stochastic excitation power spectral density on the averaged power of the deterministic response component

圖9 隨機激勵功率譜強度對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.9 Influence of the stochastic excitation power spectral density on the time-averaged standard deviation of the response

可見，隨機激勵功率譜強度對確定性響應(yīng)平均功率和隨機響應(yīng)標(biāo)準(zhǔn)差平均有較大影響。由圖8，9 可知隨著隨機激勵功率譜強度的增大，確定性響應(yīng)平均功率逐漸減小，隨機響應(yīng)標(biāo)準(zhǔn)差的時間平均逐漸增大。所有隨機激勵強度下兩種方法所得結(jié)果均吻合較好。

綜上，對于Duffing 系統(tǒng)在確定性諧波與調(diào)制白噪聲聯(lián)合作用下的響應(yīng)，本文所建議的方法在不同參數(shù)設(shè)置情況下均有良好的適用性。

3.2 Qs(t)為色噪聲

令Q(t)為均勻調(diào)制色噪聲，即：

式中Qs(t)為平穩(wěn)色噪聲。設(shè)功率譜密度為金井清譜［31］的形式：

式中ζg為場地阻尼比，ωg為特征周期。則成型濾波器的維度n=2，式（21）和（22）退化為：

式中

則式（24）中

式（26）中

且式（28）中

將式（35），（49）和式（50）代入式（28）中，并與式（31）聯(lián)立，可得一組微分方程，具體步驟可見文獻(xiàn)［18］。

用龍格-庫塔法解此非線性常微分方程組可得響應(yīng)確定性與隨機分量。

3.2.1 典型參數(shù)設(shè)置的情況

選取系統(tǒng)參數(shù)m=1，c=0.4，k=1，ε=0.5；確定性激勵參數(shù)a0=0.5，ω0=1.2 rad/s；隨機激勵參數(shù)利用譜表現(xiàn)方法生成色噪聲樣本10000條。

利用本文所建議方法與蒙特卡洛模擬得到了位移的非平穩(wěn)均值及標(biāo)準(zhǔn)差，其對比如圖10（a），（b）所示。

從圖10（a），（b）可以看出，在考慮的參數(shù)設(shè)置情況下，本文建議方法與蒙特卡洛模擬的均值較吻合。在隨機響應(yīng)分量標(biāo)準(zhǔn)差方面，所建議方法能很好地捕捉突加激勵、確定性激勵分量和隨機激勵分量調(diào)制函數(shù)帶來的非平穩(wěn)性。

圖10 Duffing系統(tǒng)在諧波與調(diào)制色噪聲激勵聯(lián)合作用下的位移Fig.10 Displacement of the considered Duffing system subjected to combined harmonic and modulated colored noise

3.2.2 非線性強度的影響

同樣地，保持其他參數(shù)不變，非線性強度系數(shù)ε取0～1。采用所建議方法與蒙特卡洛模擬計算得到的確定性響應(yīng)分量的平均功率對比如圖11所示；標(biāo)準(zhǔn)差時間平均對比如圖12所示。

圖12 非線性強度對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.12 Time-averaged standard deviation of the stochastic response component of the Duffing system with different degrees of nonlinearity

從圖11，12可以看出，在其他參數(shù)不變的情況下，確定性響應(yīng)平均功率隨著非線性強度增加先增大后降低；位移標(biāo)準(zhǔn)差的時間平均隨非線性強度增加而降低；速度標(biāo)準(zhǔn)差的時間平均在非線性程度較小時略有升高，之后隨非線性程度增加而減小。在所有非線性強度情況下，兩種方法所得結(jié)果均吻合較好。

圖11 非線性強度對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.11 Averaged power of the deterministic response component of the Duffing system with different degrees of nonlinearity

3.2.3 激勵幅值的影響

為研究確定性激勵幅值對此方法適用性的影響，保持其他參數(shù)不變，確定性激勵幅值a0取0～2。本文所建議方法與蒙特卡洛模擬計算得到的確定性響應(yīng)分量的平均功率對比如圖13所示，時變標(biāo)準(zhǔn)差的時間平均值對比如圖14所示。

圖13 諧波激勵幅值對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.13 Averaged power of the deterministic response component versus the harmonic excitation amplitude of the Duffing system

圖14 諧波激勵幅值對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.14 Influence of the harmonic excitation amplitude on the time-averaged standard deviation of the Duffing system

由圖13，14 可見，在其他參數(shù)不變的情況下，確定性響應(yīng)平均功率隨諧波激勵幅值增大而增大；位移標(biāo)準(zhǔn)差的時間平均隨諧波激勵幅值增大而減小，速度標(biāo)準(zhǔn)差的時間平均隨諧波激勵幅值增大先增大后減小，這一點與圖5所示的隨機激勵為白噪聲的情況有所不同。在所有確定性激勵幅值的情況下，兩種方法所得的結(jié)果均吻合較好。

3.2.4 確定性諧和激勵頻率的影響

為研究確定性激勵頻率對此方法適用性的影響，保持其他參數(shù)不變，確定性激勵頻率ω0取0.1～2 rad/s。本文所建議方法與蒙特卡洛模擬得到的確定性響應(yīng)分量的平均功率對比如圖15所示；隨機響應(yīng)非平穩(wěn)標(biāo)準(zhǔn)差的時間平均如圖16所示。

圖15 諧波激勵頻率對Duffing 系統(tǒng)確定性響應(yīng)平均功率的影響Fig.15 Influence of the harmonic excitation frequency on the averaged power of the deterministic response component

圖16 諧波激勵頻率對Duffing 系統(tǒng)隨機響應(yīng)標(biāo)準(zhǔn)差時間平均的影響Fig.16 Influence of the harmonic excitation frequency on the time-averaged standard deviation of the response

可見，諧波激勵頻率對確定性響應(yīng)平均功率和隨機響應(yīng)標(biāo)準(zhǔn)差的時間平均有較大影響，存在使它們達(dá)到極值的諧波激勵頻率。對于前者，達(dá)到極值的頻率為1.2 rad/s；對于后者，達(dá)到極值的頻率為1.3 rad/s。此外，不同激勵頻率下本文所建議方法與蒙特卡洛模擬所得結(jié)果均吻合較好。

綜上所述，對于Duffing 系統(tǒng)在確定性諧波與隨機激勵聯(lián)合作用下的響應(yīng)，本文所建議方法在不同參數(shù)設(shè)置的情況下，均有良好的適用性。所建議方法的計算效率較蒙特卡洛模擬有顯著優(yōu)勢。以調(diào)制白噪聲與確定性激勵聯(lián)合作用下的響應(yīng)計算為例，所建議方法僅需0.01 s，而10000 個樣本的蒙特卡洛模擬需要1.62 s，效率提高162 倍。

4 結(jié)論與展望

提出了一種用于求解確定性諧波與非平穩(wěn)隨機激勵聯(lián)合作用下，單自由度非線性隨機動力系統(tǒng)非平穩(wěn)響應(yīng)的統(tǒng)計線性化方法。首先，將系統(tǒng)響應(yīng)分解為確定性周期和零均值隨機分量之和，得到了與原非線性運動方程等效的、兩個耦合的且分別以確定性和隨機動力響應(yīng)為未知量的非線性運動微分方程。隨后，利用非平穩(wěn)統(tǒng)計線性化方法，將得到的非線性隨機運動方程化為了等效線性方程。最后，同時考慮與等效線性隨機微分方程對應(yīng)的李雅普諾夫方程和諧波激勵下的確定性運動方程，利用數(shù)值方法聯(lián)立求解了所有確定性未知量與隨機未知量。蒙特卡洛模擬驗證了此方法的適用性和精度。

數(shù)值算例表明，無論隨機噪聲為調(diào)制白噪聲或調(diào)制色噪聲，該方法均能準(zhǔn)確地求解系統(tǒng)的確定性周期響應(yīng)；求解的隨機動力響應(yīng)分量的非平穩(wěn)二階矩也具有較好的精度。此外，該方法適用于不同的參數(shù)設(shè)置情況。可進(jìn)一步將該方法拓展應(yīng)用于（多自由度）滯回非線性系統(tǒng)和分?jǐn)?shù)階非線性系統(tǒng)、（不可分）完全非平穩(wěn)激勵等情況。