孫攀旭，楊 紅,2

(1. 重慶大學土木工程學院，重慶 400045；2. 重慶大學山地城鎮(zhèn)建設與新技術教育部重點實驗室，重慶 400045)

合理的阻尼模型是結構動力響應計算結果正確的關鍵之一，黏性阻尼模型和復阻尼模型是最為常用的兩種阻尼模型[1-3]。黏性阻尼模型具有數(shù)學簡易性，計算過程簡便[4-5]，但結構每周期耗散能量與外激勵頻率相關，這與大部分工程材料的試驗結果不符[6]。復阻尼模型是在試驗結果的基礎上建立的，且具有每周期耗散能量與外激勵頻率無關的優(yōu)點，但其自由振動通解中包含發(fā)散項，導致時域計算結果不能穩(wěn)定收斂[7-9]。如何克服復阻尼模型的缺陷，成為亟待解決的問題。

關于復阻尼模型的時域發(fā)散問題，常見的解決途徑是結合其他阻尼模型的優(yōu)點，構造出等效復阻尼模型，以實現(xiàn)動力響應的時域計算。例如，采用復頻率法[10-11]、應變能法[12-13]等計算出等效阻尼比，可得到等效黏性阻尼模型；Boit[14]提出了核函數(shù)為Voigt 函數(shù)的等效指數(shù)阻尼模型，為保證計算精度，可將核函數(shù)進一步拓展為Maxwell-Wiechert 函數(shù)[15-17]。廖振鵬[18]、周正華等[19]、Yang 等[20]為避免指數(shù)阻尼的復雜計算過程，采用最小二乘法得到一種等效于復阻尼模型的五參數(shù)黏彈性本構模型；Wang[21]將復阻尼矩陣等效為Rayleigh 阻尼矩陣，提出了等效Rayleigh阻尼模型，劉慶林等[22]在此基礎上進一步提出了等效Caughey 阻尼模型。

上述等效復阻尼模型的主要問題是誤差估計較為困難，且等效準則的合理性難以判定。因此，一些學者直接對修正復阻尼模型進行研究，期望解決其缺陷。朱鏡清和朱敏[23]、Pan 等[24]直接舍棄自由振動通解中發(fā)散項的計算方法，保證了計算結果的穩(wěn)定性，但這種處理在數(shù)學上不夠精確[25-26]，并且也沒有克服復阻尼模型本身存在的缺陷。Clough 和Penzien[1]、Chen 等[27]對復阻尼模型進行了修正，假定阻尼力的大小與結構體系的位移大小成正比，且與速度的方向相反，提出了遲滯阻尼模型，以保證結構每周期耗散能量與外激勵頻率無關，但存在耗能不守恒和非線性的問題。筆者[28]在遲滯阻尼模型的基礎上，利用能量守恒得到了改進遲滯阻尼模型，但沒有解決非線性的問題。朱鏡清[29]和筆者[30]在復阻尼模型的基礎上，構建出頻率相關黏性阻尼模型，以處理時域發(fā)散問題，但未考慮正、負頻率共軛的影響。Inaudi 和Kelly[31]引入正、負頻率共軛的條件，對復阻尼模型的頻域運動方程進行了修正，并提出了對應的滯變阻尼模型。滯變阻尼模型不僅保證了運動方程的完備性，且有效解決了復阻尼模型自由振動通解中的發(fā)散問題[32]，但滯變阻尼模型的有阻尼自振頻率隨著損耗因子的增加而增加，與實際不符。

在上述研究的基礎上，本文在保留滯變阻尼模型每周期耗散能量與外激勵頻率無關的優(yōu)點的情況下，通過構建改進滯變阻尼模型，在復數(shù)域內建立改進滯變阻尼模型的拉氏運動方程及相應的時域運動方程。然后依據(jù)信號識別方法，結合改進滯變阻尼模型的時域運動方程特點，提出基于傅里葉變換的時域計算方法和基于希爾伯特-黃變換的時域計算方法，以實現(xiàn)結構動力響應的時域計算。

1 改進滯變阻尼模型

1.1 復阻尼模型

在諧波作用下，大部分工程結構存在穩(wěn)態(tài)反應每周期耗散能量與外激勵頻率無關的試驗現(xiàn)象，依據(jù)試驗現(xiàn)象可構建出復阻尼模型的頻域運動方程[33]：

-?2mY(i?)+iηkY(i?)+kY(i?)=-mP(i?) (1)

式中： m 為結構的質量； k 為結構的剛度； η為結構的損耗因子；P(i?) 為p(t)的傅里葉變換項；

式中： A 為諧波的振幅； θ為諧波的振動頻率。

式中：y(t)為結構位移響應的復數(shù)表達式；實部x(t)為結構的位移響應。

采用傅里葉逆變換，式(1)對應的時域運動方程為：

m¨y(t)+iηky(t)+ky(t)=-mp(t)(4)

式(4)是建立在諧波作用下結構的穩(wěn)態(tài)反應基礎上的，僅適用于分析諧波作用下結構的穩(wěn)態(tài)反應。將式(4)進一步推廣到外激勵作用下包含結構瞬態(tài)反應和穩(wěn)態(tài)反應的時域過程中，對應的運動方程為：

m¨y(t)+iηky(t)+ky(t)=-mg(t)(5)

式中，g(t)為外激勵加速度。

文獻[34]的研究結果表明式(5)不是完備方程，求解得到的結構動力響應明顯偏小。為保證方程的完備性，利用復化對偶原則[25,35]對式(5)進行修正，可得對應的時域運動方程為：

m¨y(t)+iηky(t)+ky(t)=-m[g(t)+ig′(t)](6)

式中，g′(t) 為g(t)的復化對偶項。

式(6)對應的自由振動方程為：

式(7)的自由振動響應表達式為：

式中， A1、 A2、 A3和 A4為待定系數(shù)，可由初始條件進行確定。

由式(8)和式(9)可知復阻尼模型的自由振動響應中包含指數(shù)增長項，其時域計算結果將會隨著外部作用時間的增加，逐漸出現(xiàn)發(fā)散現(xiàn)象。

對式(6)進行傅里葉變換，得到對應的頻域運動方程為：

式中：G(i?) 為g(t)的傅里葉變換項；G′(i?)為g′(t)的傅里葉變換項。

取式(10)中結構位移響應對應的頻域運動方程，可得：

式中，X(i?) 為x(t)的傅里葉變換項。

1.2 滯變阻尼模型

為考慮負頻率的影響，將式(11)進行改進，得到基于復阻尼模型的滯變阻尼模型頻域運動方程為[31-32]：

其中：

對式(12)進行傅里葉逆變換，可得到滯變阻尼模型的時域運動方程為：

式中，xH(t) 為x(t)的希爾伯特變換項。

式(14)對應的自由振動響應的表達式為：

式中， B1和 B2為待定系數(shù)，可由初始條件進行確定。

由式(15)和式(16)可知，滯變阻尼模型的自由振動響應中不包含發(fā)散項，保證了時域計算結果的穩(wěn)定收斂。

復化對偶原則為：

希爾伯特變換原則為：

式中：[sin(θt)]H為sin(θt)的希爾伯特變換；[cos(θt)]H為cos(θt)的希爾伯特變換。

由式(17)和式(18)對比可知，復阻尼模型的復化對偶原則實質上就是希爾伯特變換。因此，滯變阻尼模型在保留了復阻尼模型的完備性基礎上，克服了復阻尼模型的時域發(fā)散缺陷。但由式(16)可知，隨著損耗因子的增加，有阻尼自振頻率 ωh將逐漸增大。這種現(xiàn)象與實際不符。

1.3 改進滯變阻尼模型

式(12)在計算過程中采用的是復平面法，?為復頻率，對應的符號函數(shù)sgn(?)與式(13)的定義是矛盾的，從而導致滯變阻尼模型的時域方程中有阻尼自振頻率隨著損耗因子的增加而增加。針對該問題，考慮到復頻率的虛部為結構的振動頻率，在頻域內本文提出改進滯變阻尼模型，即將式(12)進一步改進為：

顯然， ?為復頻率時，在頻域內定義運動方程是不準確的，此時，采用復數(shù)域內的拉氏運動方程表示改進滯變阻尼模型是更符合實際情況的做法，借助于拉普拉斯變換，構建復數(shù)域內的改進滯變阻尼模型運動方程為：

式中：X(s) 為x(t)的拉普拉斯變換；G(s) 為g(t)的拉普拉斯變換。

采用拉普拉斯逆變換，可得到改進滯變阻尼模型的時域運動方程為：

式中， ?v為結構動力響應的振動頻率。

式(21)對應的自由振動方程為：

對式(22)進行拉普拉斯變換后，可得：

式(23)包含符號函數(shù)和虛部函數(shù)，求解X(s)是困難的，因此，可采用特征值法求解式(22)。

改進滯變阻尼模型的特征值方程為：

令：

將式(25)代入式(24)，可得：

求解式(26)，可得：

進一步得到2 個共軛特征值為：

依據(jù)特征值和常數(shù)變易法，可得到對應的自由振動響應表達式為：

其中：

由式(29)和式(30)可知，本文基于復數(shù)域內的拉氏運動方程創(chuàng)建的改進滯變阻尼模型的自由振動響應中不包含發(fā)散項，時域計算結果是穩(wěn)定收斂的，可克服復阻尼模型的缺陷；隨著損耗因子的增加，有阻尼自振頻率 ωg隨著損耗因子的增加而逐漸減小，可克服傳統(tǒng)滯變阻尼模型的缺陷。

2 改進滯變阻尼模型的時域計算方法

對于基于改進滯變阻尼模型的比例阻尼體系，可直接采用傳統(tǒng)的實模態(tài)疊加法，將多自由度體系運動方程分解為單自由度體系運動方程，進而采用實模態(tài)疊加實現(xiàn)動力響應的時域計算。對于基于改進滯變阻尼模型的非比例阻尼體系，可將多自由度體系運動方程轉化為復數(shù)域內拉氏運動方程，進一步采用復模態(tài)疊加法將多自由度體系運動方程分解為單自由度體系運動方程，進而采用復模態(tài)疊加實現(xiàn)動力響應的時域計算。因此，為求解基于改進滯變阻尼模型的單自由度體系時域運動方程，本文分別提出基于傅里葉變換和希爾伯特-黃變換的單自由度體系動力響應的時域計算方法。

2.1 基于傅里葉變換的時域計算方法

基于改進滯變阻尼模型的時域運動方程中包含有振動頻率的未知項，因此，無法通過常用的時程計算方法包括常平均加速度法、Newmark-β法、中心差分法等逐步積分法[1-2]計算其時域動力響應。為此，采用快速傅里葉變換法，將外激勵加速度用傅里葉級數(shù)展開[36]，可得：

將式(31)代入式(21)，可得：

令式(32)的一般解為：

式中：xc(t)為式(32)對應的齊次方程通解；xp(t)為式(32)對應的非齊次方程特解。

令：

將式(33)和式(34)代入式(32)，可得：

式中：xc(t)為式(35)的通解；xp,0(t)為式(36)的特解；xp,n(t)為式(37)的特解。

式(35)的通解與自由振動響應相同，可表示為：

式中， D1和 D2為待定系數(shù)。

式(36)的特解為：

令式(37)的特解為：

將式(40)代入式(37)，可得：

由式(38)～式(41)可得到式(32)的一般解為：

其中：

綜上，由式(42)可實現(xiàn)基于傅里葉變換的動力響應時域計算方法。該方法是一種全局過程的計算方法，可計算出任意時刻的結構動力響應，且無條件穩(wěn)定收斂。基于傅里葉變換的時域方法可在MATLAB 軟件平臺上通過編程實現(xiàn)，對應的算法流程如圖1 所示。

圖 1 基于傅里葉變換的時域方法流程圖Fig.1 Flow chart of the proposed time-domain method based on Fourier transform

2.2 基于希爾伯特-黃變換的時域計算方法

首先將時間離散化，即按照步長 Δt對時間進行離散，任意時刻可表示為tk=kΔt(k=0,1,2,···)。采用希爾伯特-黃變換法[37]對外激勵加速度進行識別分析，由經(jīng)驗模態(tài)分解(EMD)將外激勵加速度分解為s 條本征模函數(shù)(imf)和殘余分量，進一步利用希爾伯特變換識別出imf 的瞬時頻率、瞬時振幅等信息。由此， tk時刻到tk+1時刻，外激勵加速度可表示為：

式中： θj,k為第j 條imf 分量 tk時刻到tk+1時刻的瞬時頻率； φj,k為第j 條imf 分量tk時刻的瞬時相位；Ij(t)為第j 條imf 分量的瞬時振幅；r(t)為殘余分量。

第j 條分量tk+1時刻的瞬時相位可表示為：

將式(44)代入式(21)，可得：

令式(46)的一般解為：

式中：xc(t)為式(46)對應的齊次方程通解；xp(t)為式(46)對應的非齊次方程特解。

令：

將式(47)和式(48)代入式(46)，可得：

式中：xc(t)為式(49)的通解；xp,0(t)為式(50)的特解；xp,j(t)為式(51)的特解。

式(49)的通解與自由振動響應相同，可表示為：

式(50)的特解為：

令式(51)的特解為：

將式(54)代入式(51)可得：

求解式(55)，可得：

將式(52)～式(56)代入式(47)，可得tk+1時刻結構的位移為：

tk+1

由式(57)可進一步得到 時刻結構的速度為：

式中：

依據(jù)式(57)～式(59)，可由 tk時刻結構的位移和速度計算出tk+1時刻的結構響應，從而實現(xiàn)了基于希爾伯特-黃變換的動力響應時域計算方法。該方法是一種時域逐步計算方法，同時也是無條件收斂的?；谙柌?黃變換的時域方法可在MATLAB軟件平臺上通過編程實現(xiàn)，對應的算法流程如圖2所示。

圖 2 基于希爾伯特-黃變換的時域方法流程圖Fig.2 Flow chart of the proposed time-domain method based on Hilbert-Huang transform

3 算例分析

3.1 諧波作用下的時域精確計算方法

隨機激勵可近似視為諧波的疊加，因此分析諧波作用下結構的動力響應具有代表性。以兩條諧波組合得到的外激勵信號為例，對比分析本文提出的時域計算方法。

采集到的外激勵信號通常情況下是離散的，按照步長 Δt對時間進行離散，任意時刻可表示為tk=kΔt(k=0,1,2,···)。外激勵加速度可表示為：

式中：

式中：h1(t) 和h2(t)分別為第一諧波分量和第二諧波分量的振幅函數(shù)； φk和 φk分別為第一諧波分量和第二諧波分量 tk時刻的相位；φk+1和φk+1分別為第一諧波分量和第二諧波分量tk+1時刻的相位；q1(t) 和q2(t)分別為第一諧波分量和第二諧波分量的振動頻率函數(shù)。

與基于希爾伯特-黃變換的時域計算方法相同，依據(jù)外激勵加速度的表達式，采用齊次方程通解和非齊次方程特解分別求解并疊加計算的方法，可得到諧波作用下改進阻尼模型的時域精確計算結果。

3.2 算例1

以損耗因子為0.5，自振頻率為30 rad/s 的單自由度體系為例，對體系作用外激勵R-1，激勵作用時間為15 s。外激勵R-1 的具體參數(shù)為：外激勵加速度的采樣頻率為100 Hz，對應的時間步長Δt為0.01 s，第一諧波分量和第二諧波分量的參數(shù)見式(63)，外激勵R-1 的加速度時程如圖3 所示。

考慮到快速傅里葉變換法對采樣點數(shù)的要求，首先對外激勵信號進行補零處理[1](下同)。隨后，對外激勵R-1 進行快速傅里葉變換，可得到對應的傅里葉變換譜如圖4 所示，表明其可準確識別出外激勵R-1 中兩條諧波的頻率。

圖 3 外激勵R-1 的加速度時程Fig.3 Acceleration time-history of external excitation R-1

圖 4 外激勵R-1 的傅里葉變換譜Fig.4 Fourier transform spectrum of external excitation R-1

考慮到希爾伯特-黃變換的端點效應問題，采用正弦波法對[38-39]外激勵R-1 的兩端進行延拓處理，進而對外激勵R-1 進行EMD 計算，當殘余分量為單調函數(shù)或最大值小于0.05 m/s2時停止分解[37,40](下同)。由外激勵R-1 的EMD 結果可知，imf1 和imf2 是主要分量，imf1 對應的瞬時振幅和瞬時頻率，與外激勵R-1 中第一條諧波分量(hw1)的瞬時振幅和瞬時頻率近似相等(見圖5(a)和圖6(a))；imf2 的瞬時振幅和瞬時頻率則與外激勵R-1 中第二條諧波分量(hw2)的瞬時振幅和瞬時頻率近似相等(見圖5(b)和圖6(b))。因此，希爾伯特-黃變換法可準確識別出外激勵R-1 中兩條諧波的振幅和頻率。

圖 5 外激勵R-1 的瞬時振幅對比Fig.5 Comparison of instantaneous amplitudes of external excitation R-1

圖 6 外激勵R-1 的瞬時頻率對比Fig.6 Comparison of instantaneous frequencies of external excitation R-1

分別采用基于傅里葉變換的時域計算方法(IF)、基于希爾伯特-黃變換的時域計算方法(IH)和時域精確計算方法(IP)計算該算例單自由度體系的加速度時程，結果如圖7 所示。IF、IH 和IP 計算的結構加速度時程近似相等(見圖7)，IF和IH 的加速度峰值相對誤差都小于5%(見表1)，證明了IF 和IH 的正確性。IF 是一種全局過程的計算方法，可計算任意時刻的結構動力響應，不用進行時域逐步計算，結果更為直觀。IF 利用快速傅里葉變換將任意外激勵信號分解為一系列三角函數(shù)分量和常數(shù)分量，僅產生較小的近似誤差，對應的數(shù)學證明是嚴謹?shù)?。但IF 的計算耗時遠大于IH(見表1)，這是因為外激勵信號進行傅里葉變換時，需要分解為Nc/2 條諧波，IH 進行希爾伯特-黃變換時僅分解為8 條信號分量。相比IF，IH 是一種時域逐步計算方法，計算結果依賴于前一時刻的結果；同時希爾伯特-黃變換法中存在端點效應，為初始時刻的計算結果引入了誤差，并隨著時域逐步計算影響整個過程的計算精度。由于目前常用的端點效應處理方法并不能完全消除端點效應的影響，導致IH 的計算誤差不易估計。此外，EMD 分解法和端點效應消除方法是一種經(jīng)驗分解方法，導致IH 缺乏嚴謹?shù)臄?shù)學理論支持。

圖 7 外激勵R-1 作用下的結構加速度時程響應Fig.7 Structural acceleration responses time-history under external excitation R-1

表 1 外激勵R-1 作用下結構的動力響應對比Table 1 Comparison of structural dynamic responses under external excitation R-1

3.3 算例2

采用與算例1 相同的單自由度模型，對體系作用外激勵R-2，激勵作用時間為15 s。外激勵R-2的具體參數(shù)為：外激勵加速度的采樣頻率為100 Hz，對應的時間步長 Δt為0.01 s，第一諧波分量和第二諧波分量的參數(shù)見式(64)，外激勵R-2 的加速度時程如圖8 所示。

圖 8 外激勵R-2 的加速度時程Fig.8 Acceleration time-history of external excitation R-2

與算例1 不同，算例2 的外激勵R-2 的兩條諧波的瞬時頻率是變化的。計算結果表明：由于傅里葉變換法是從能量角度對外激勵信號進行分解，對于隨時間變化的瞬時頻率，不能進行準確識別，僅能確定外激勵R-2 的頻率范圍(見圖9)。相比傅里葉變換法，希爾伯特-黃變換法可識別出外激勵R-2 中諧波分量的瞬時振幅和瞬時頻率。由外激勵R-2 的EMD 結果可知，imf1 和imf2 是主要分量，可對比分析imf1 和hw1、imf2 和hw2的瞬時振幅和瞬時頻率，結果如圖10 和圖11所示。

圖 9 外激勵R-2 的傅里葉變換譜Fig.9 Fourier transform spectrum of external excitation R-2

分別采用IF、IH 和IP 計算外激勵R-2 作用下體系的加速度時程，結果如圖12 所示。表2 所示，與算例1 的計算結果相比，IF 和IH 的計算精度都有所下降，IF 的加速度峰值相對誤差由0.42%增加到25.20%，IH 的加速度峰值相對誤差由1.28%增加到9.00%。與恒定頻率的外激勵R-1 相比，希爾伯特-黃變換法對具有瞬時變化頻率的外激勵R-2 的識別精度有所下降，但可粗略識別出其變化趨勢。相比IF，IH 的計算量更少、計算精度更高。傅里葉變換方法受限于W.Heisenberg 不確定準則，無法有效識別出外激勵信號的瞬時頻率等信息，而希爾伯特-黃變換方法具有自適應性，可有效識別出任意非平穩(wěn)外激勵信號的瞬時信息。因此，與IF 相比，IH 具有更寬泛的適用范圍。

圖 10 外激勵R-2 的瞬時振幅對比Fig.10 Comparison of instantaneous amplitudes of external excitation R-2

圖 11 外激勵R-2 的瞬時頻率對比Fig.11 Comparison of instantaneous frequencies of external excitation R-2

圖 12 外激勵R-2 作用下的結構加速度時程響應Fig.12 Structural acceleration responses time-history under external excitation R-2

表 2 外激勵R-2 作用下結構的動力響應對比Table 2 Comparison of structural dynamic responses under external excitation R-2

3.4 算例3

以損耗因子為0.5 的多自由度體系為例，質量剛度分布如圖13 所示。對體系作用外激勵R-2，激勵作用時間為15 s，結合模態(tài)疊加法，可分別采用IF、IH 和IP 計算外激勵R-2 作用下體系頂層的加速度時程，結果如圖14 所示。表3 所示，IF 的加速度峰值相對誤差為25.17%，IH 的加速度峰值相對誤差為4.05%。與算例2 的計算結果規(guī)律一致，IF 的相對誤差最大，且計算時間最長。

圖 13 模型示意圖Fig.13 Schematic of model

圖 14 外激勵R-2 作用下的結構頂層加速度時程響應Fig.14 Structural top acceleration responses time-history under external excitation R-2

表 3 外激勵R-2 作用下結構頂層的動力響應對比Table 3 Comparison of structural top dynamic responses under external excitation R-2

4 討論

通過分析改進滯變阻尼模型的特點，并對比基于傅里葉變換的時域計算方法和基于希爾伯特-黃變換的時域計算方法可知：

(1)改進滯變阻尼模型的阻尼矩陣構造容易，僅依賴于材料損耗因子和結構剛度矩陣，在后續(xù)研究中可更方便地用于非比例阻尼體系的動力響應分析，如設置了耗能阻尼器的建筑結構等。

(2)基于改進滯變阻尼模型的時域計算方法的缺陷來源于信號識別方法的局限性。相比傅里葉變換法、小波變換法，希爾伯特-黃變換法可以更加準確地識別振動頻率與時間的關系，但常用的端點效應處理法和EMD 分解法存在缺陷。在后續(xù)研究中，可采用優(yōu)化極值延拓、神經(jīng)網(wǎng)絡預測等方法減小端點效應；可采用多項式曲線、B 樣條曲線、冪函數(shù)曲線等代替?zhèn)鹘y(tǒng)樣條曲線進行插值，以提高計算精度。

5 結論

本文依據(jù)基于復阻尼模型的改進滯變阻尼模型，提出了不同的時域計算方法，并通過算例分析得出如下結論：

(1)在復阻尼模型的基礎上，提出了基于復阻尼模型的改進滯變阻尼模型，創(chuàng)建了對應的復數(shù)域內拉氏運動方程和時域運動方程。改進滯變阻尼模型不僅克服了復阻尼模型的時域計算發(fā)散缺陷，還克服了滯變阻尼模型的有阻尼自振頻率隨損耗因子增加而增加的缺陷。依據(jù)外激勵加速度的識別方法，分別提出了基于傅里葉變換、基于希爾伯特-黃變換的時域計算方法。

(2)基于傅里葉變換的時域計算方法是一種全局過程的計算方法，可直接計算任意時刻的結構動力響應，但計算量較大，且無法識別瞬時頻率。求解頻率隨時間變化的外激勵作用下結構動力響應時，其計算精度較低。

(3)基于希爾伯特-黃變換的時域計算方法是一種時域逐步積分算法，計算量較小、適用范圍更廣。EMD 分解和端點效應影響外激勵加速度的識別精度，導致該方法的計算誤差不易估計。