非線性渦激力模型氣動參數(shù)識別新方法

翁祥穎， 葛耀君， 陳海興

(同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗室，上海 200092)

流經(jīng)鈍體斷面的流體在鈍體下游表面形成漩渦且有規(guī)律脫落引起的物體振動稱渦振。渦振為典型的限幅自激振動，由于Van der Pol 振子具有類似特性，現(xiàn)有橋梁渦振單自由度渦激力模型均基于Van der Pol振子構(gòu)造而成。此模型均含與振動系統(tǒng)折算頻率相關(guān)的氣動參數(shù)，使用該模型時需據(jù)風(fēng)洞試驗結(jié)果識別系統(tǒng)的渦激力參數(shù)。Ehsan等[1]通過賦予處于渦振狀態(tài)的節(jié)段模型大于其渦振振幅初位移并測量模型振動衰減時程以確定模型中與非線性阻尼項相關(guān)的氣動參數(shù)，識別時直接假定渦振為簡單簡諧振動。而該簡諧振動并非基于非線性渦激力模型結(jié)構(gòu)渦振運(yùn)動平衡方程的解，該假定與方程解析解間的較大差別必對渦激力模型氣動參數(shù)識別結(jié)果產(chǎn)生影響。Gupta[2]基于不變嵌入理論提出非線性最小二乘識別法。對不同初始值該方法會獲得不同結(jié)果甚至擬合不收斂，但并未就如何正確設(shè)定初始值給出建議。Marra等[3-4]基于緩變參數(shù)法亦提出渦激力參數(shù)識別方法，然而在將識別的氣動參數(shù)代入渦激力模型并反算節(jié)段模型響應(yīng)時發(fā)現(xiàn)反算響應(yīng)功率譜與實(shí)測相比除渦脫頻率外出現(xiàn)一額外卓越頻率。另外，以上識別方法其本質(zhì)均為非線性最小二乘法，當(dāng)給擬合目標(biāo)設(shè)定不同初始值時，擬合所得結(jié)果亦不同，不便于研究及應(yīng)用。對此，本文基于遺傳算法提出具有較高穩(wěn)定性的非線性單自由度渦激力氣動參數(shù)識別方法。

1 單自由度經(jīng)驗非線性渦激力模型

該模型由Scanlan[5]提出，因其具有較高精度而在橋梁渦振研究中廣泛應(yīng)用。基于Van der Pol振子，該模型由含非線性效應(yīng)的氣動阻尼、氣動剛度自激力及瞬時簡諧強(qiáng)迫力共同構(gòu)成：

(1)

定義無量綱量：

(3)

式中：m為系統(tǒng)每米模態(tài)質(zhì)量。

結(jié)合簡化的渦激力模型得渦振運(yùn)動平衡微分方程的無量綱形式為：

(4)

其中：

μ1=-[2ξKn-mrY1(K)]

(5)

μ2=mrY1(K)ε

(6)

(7)

式中：ξ為系統(tǒng)阻尼比。顯然，式(4)描述的系統(tǒng)振動特性與μ1，μ2密切相關(guān)：μ1，μ2符號相反時，系統(tǒng)會發(fā)生衰減或發(fā)散振動，不符合渦振發(fā)生時結(jié)構(gòu)呈等幅振動特性。由渦振為穩(wěn)定的等幅振動知μ1為正數(shù)，故定義變量：

(8)

代入以上無量綱渦振運(yùn)動平衡微分方程可得標(biāo)準(zhǔn)的Van der Pol振子為：

(9)

2 Van der Pol振子近似解

Van der Pol振子含非線性阻尼項，其精確的閉合解析解目前仍難以獲得。傳統(tǒng)的近似求解方法包括KBM法、諧波平衡法等均要求振子中阻尼系數(shù)μ1為接近于零的小參數(shù)。對渦振而言，由于μ1與氣動參數(shù)直接相關(guān)，無法先了解其大小，因而無法保證傳統(tǒng)解法是否適用。為此本文采用不依賴小參數(shù)的近似級數(shù)求解方法-同倫分析法(Homotopy Analysis Method)[6]。該法求解非線性微分方程核心為在方程未知量中引入控制參數(shù)p，使未知量變?yōu)楹琾的多元函數(shù)，p由0變到1的過程中，未知量則由微分方程對應(yīng)的線性方程解變到原始的非線性微分方程解。

為滿足同倫分析法要求各未知量相對于控制參數(shù)p求各階導(dǎo)數(shù)條件，需先將式(9)中兩未知量-頻率Kn及振動幅值a由隱式變?yōu)轱@式，為此定義變量為：

τ=ωs,y(s)=au(τ)

(10)

代入式(9)得適用于同倫分析法的基本方程為：

ω2u″(τ)+μ3u(τ)=μ1ω[1-a2u2(τ)]u′(τ)

(11)

由同倫分析法得Van der Pol振子二階近似解為：

(12)

(13)

(14)

重復(fù)以上過程，可得更高階近似解。但三階以上的解是關(guān)于頻率大于等于7ωn的高頻成分，考慮橋梁節(jié)段模型渦振響應(yīng)主導(dǎo)振動頻率較單一，二階近似解已能較好反映出振動特性。

3 擬合目標(biāo)函數(shù)

結(jié)合風(fēng)洞試驗實(shí)測橋梁節(jié)段模型渦振響應(yīng)時程及所得二階近似解，構(gòu)造表示殘差平方和的擬合目標(biāo)函數(shù)為：

(15)

式中：ym(s)，ye(μ1,μ3,s)為實(shí)測渦振響應(yīng)與計算所得響應(yīng)；s為無量綱時間；N為位移響應(yīng)時程樣本容量。

4 遺傳算法識別參數(shù)

對非線性擬合問題，傳統(tǒng)方法如牛頓法、擬牛頓法、梯度下降法等與梯度相關(guān)的擬合方法給出的擬合結(jié)果往往只是所給初始值附近局部極值，無法識別定義域范圍內(nèi)全局最優(yōu)值，故此方法識別結(jié)果與擬合目標(biāo)初始值設(shè)定直接相關(guān)，不便于研究及應(yīng)用。為克服該缺陷，采用旨在搜索全局最優(yōu)值的啟發(fā)式擬合法-遺傳算法，該算法主要思想為模擬自然界生物優(yōu)勝劣汰的進(jìn)化規(guī)律。擬合對象取值范圍及擬合目標(biāo)函數(shù)一旦確定，遺傳算法則通過既定的遺傳策略，如種群數(shù)量、選擇、交叉、變異等操作，產(chǎn)生代表更優(yōu)解的下一代個體直至收斂獲得最終擬合值[7]。

基于此擬合思路，編制擬合μ1，μ3的遺傳算法程序?？紤]試驗中實(shí)測響應(yīng)信號均含隨機(jī)噪聲，本文分別采用光滑及含高斯白噪聲的有噪位移時程信號驗證所提方法的可行性及程序編制的正確性。針對渦振發(fā)生時橋梁節(jié)段模型發(fā)生較大、規(guī)則的限幅自激振動，采集的振動信號較光滑，故設(shè)置有噪信號的噪信比為2.0%。數(shù)值驗證中原始光滑信號在參數(shù)μ1=0.70，μ3=1 280條件下數(shù)值積分由式(9)獲得。遺傳算法其它參數(shù)設(shè)置：種群個體數(shù)M=50，交叉概率Pc=0.7，個體變異概率Pm=0.05，最大迭代次數(shù)N=300；參數(shù)μ1，μ3的搜索范圍分別為[0,1]、[302,402]；采用線性排序結(jié)合精英選擇的混合選擇策略。兩種數(shù)值模擬試驗均進(jìn)行30次以觀察識別結(jié)果的穩(wěn)定性。對光滑位移時程信號，各次識別結(jié)果變化較小且均接近于μ1，μ3目標(biāo)值，最終識別結(jié)果為μ1=0.708 4,μ3=1 279.99。對有噪的位移時程信號，各次試驗識別的μ3很穩(wěn)定，最大誤差僅0.3%；但μ1相對較離散。觀察并統(tǒng)計識別結(jié)果發(fā)現(xiàn)，μ3微量差別引起的相位差使識別結(jié)果殘差平方和變化劇烈并導(dǎo)致μ1識別結(jié)果較離散。為此進(jìn)行兩輪識別：第一輪同時識別μ1及μ3，因μ3較穩(wěn)定故先確定μ3；第二輪僅識別μ1，將μ3視為已知數(shù)以消除識別過程中相位差再次變化對μ1識別產(chǎn)生的劇烈影響。

圖1 原始有噪位移時程與擬合結(jié)果對比

由于加入高斯白噪聲，有噪信號識別結(jié)果穩(wěn)定性不及光滑信號，達(dá)到搜索停止條件時種群多樣性有時仍較大，因而搜索到最優(yōu)解難度、時間隨之增加。統(tǒng)計30次識別結(jié)果，據(jù)殘差平方和最小原則可確定本次數(shù)值驗證進(jìn)行30次試驗中最優(yōu)識別結(jié)果為μ1=0.704,μ3=1 280.00，與目標(biāo)值基本一致。圖1為原始有噪位移時程與據(jù)識別結(jié)果計算所得位移時程。由圖1看出，兩者吻合較好。

5 橋梁節(jié)段模型渦激力氣動參數(shù)識別

5.1 風(fēng)洞試驗概況

試驗在同濟(jì)大學(xué)TJ-3風(fēng)洞進(jìn)行。圖2為風(fēng)洞試驗節(jié)段模型所用流線型閉口箱梁斷面。箱梁節(jié)段模型寬2D=0.553 m，一階豎彎模態(tài)質(zhì)量M=10.1 kg/m，頻率f=5.664 Hz，系統(tǒng)阻尼比ξ=1.25%。均勻流場中風(fēng)速U=6.0 m/s時節(jié)段模型出現(xiàn)穩(wěn)定豎彎渦振。試驗系統(tǒng)見圖3。圖4為無量綱的豎彎渦振位移時程。對應(yīng)于該鎖定風(fēng)速，系統(tǒng)折算頻率Kn=3.221。試驗中用兩個采樣頻率100 Hz的激光位移計測量節(jié)段模型位移時程。

圖2 渦振試驗主梁節(jié)段模型流線型斷面

圖3 試驗系統(tǒng)

圖4 無量綱的豎彎渦振響應(yīng)時程

5.2 μ1，μ3參數(shù)搜索

結(jié)合式(8)、(13)，得無量綱實(shí)測渦振與Van der Pol振子位移幅值間關(guān)系為：

(17)

其中：At為無量綱實(shí)測豎彎渦振位移響應(yīng)幅值。該比值可通過迭代法確定，初始值分子取2。

圖5 極限環(huán)與擬合限幅環(huán)

橋梁節(jié)段模型發(fā)生渦振時系統(tǒng)振動頻率與零風(fēng)速時固有頻率相比變化較小，因而確定μ3的搜索范圍為[33,42]。由于Van der Pol振子在相空間內(nèi)為限幅環(huán)運(yùn)動，隨μ1的改變，限幅環(huán)形狀亦發(fā)生變化。利用此對應(yīng)關(guān)系可初步確定μ1搜索范圍。圖5(a)對應(yīng)μ1=0,0.5,1三個不同值限幅環(huán)[8]。圖5(b)為由實(shí)測豎彎渦振響應(yīng)擬合所得限幅環(huán)，發(fā)現(xiàn)其長軸方向與y軸偏差很小。對比圖5(a)結(jié)果，確定μ1的搜索范圍為[0,0.5]。遺傳算法其它參數(shù)設(shè)置同上節(jié)數(shù)值試驗。

5.3 渦激力氣動參數(shù)識別

以圖4無量綱位移響應(yīng)為擬合目標(biāo)，重復(fù)50次擬合，發(fā)現(xiàn)擬合結(jié)果大都收斂于μ1=0.394，μ3=3.242。圖6為擬合結(jié)果與原始樣本時程對比，可發(fā)現(xiàn)二者幅值、相位均吻合較好。識別出μ1，μ3后，由樣本無量綱幅值均值A(chǔ)t=4.50×10-3，通過式(17)迭代最終得μ2的值：

μ2=77732.02

(18)

據(jù)μ1,μ2,μ3定義最終確定該主梁斷面渦激力氣動參數(shù)，見表1。

表1 流線型箱梁斷面渦激力氣動參數(shù)識別結(jié)果

圖6 擬合結(jié)果與原始位移時程對比

5.4 識別結(jié)果評價

分析識別所得非線性渦激力模型三個氣動參數(shù)發(fā)現(xiàn)，識別的Y1(K)與相似橋梁斷面在基本一致的折算頻率下采用基于緩變參數(shù)的廣義諧波函數(shù)KBM法[4]識別結(jié)果較接近。與氣動剛度相關(guān)的Y2(K)為負(fù)值且絕對值相對較小，表明在該節(jié)段模型渦振鎖定區(qū)，氣動剛度作用使模型振動卓越頻率略高于固有頻率。通過該鎖定風(fēng)速下模型豎彎渦振響應(yīng)頻譜圖發(fā)現(xiàn)，系統(tǒng)振動的卓越頻率由零風(fēng)速f=5.664 Hz提高至f=5.694 Hz。該微量的提高已表明Y2(K)識別結(jié)果的合理性。與文獻(xiàn)[1,3]相比，本文識別的非線性阻尼項氣動參數(shù)ε較大。分析非線性渦激力模型式(1)，其模擬渦振限幅自激特性機(jī)理通過非線性氣動阻尼隨質(zhì)點(diǎn)振動過程變化進(jìn)而導(dǎo)致系統(tǒng)總阻尼發(fā)生正負(fù)交替變化，表明非線性阻尼項中εx2/D2必為與1有相同量級的數(shù)。基于此，本文采用無量綱位移量級10-3較文獻(xiàn)[1,3]的10-2小?？紤]折算頻率K及試驗中采用的橋梁斷面差別，本文對參數(shù)ε的識別結(jié)果具有合理性。

6 結(jié) 論

(1)通過對單自由度經(jīng)驗非線性渦激力模型進(jìn)行變換，將其轉(zhuǎn)化為標(biāo)準(zhǔn)Van der Pol振子形式。在求析解時，考慮非線性振動傳統(tǒng)求解方法對運(yùn)動平衡方程中小參數(shù)的依賴性，選擇與系統(tǒng)參數(shù)值大小無關(guān)的同倫分析法獲得Van der Pol振子的近似解。

(2)擬合過程引入能在待擬合參數(shù)定義域范圍內(nèi)進(jìn)行全局最優(yōu)搜索的遺傳算法進(jìn)行參數(shù)識別。光滑及有噪位移時程數(shù)值試驗結(jié)果表明本文所提方法具有可行性及穩(wěn)定性。

(3)利用該方法擬合風(fēng)洞試驗實(shí)測試數(shù)據(jù)時，針對遺傳算法擬合所需參數(shù)范圍因在氣動參數(shù)未知情況下較難確定問題，通過比較相平面內(nèi)極限環(huán)方法予以解決。用該方法對流線型閉口箱梁節(jié)段模型的渦激氣動參數(shù)進(jìn)行識別。擬合結(jié)果表明，該方法能有效用于風(fēng)洞試驗且識別結(jié)果穩(wěn)定性較高。

參 考 文 獻(xiàn)

[1]Ehsan F,Scanlan R H,ASCE M, Vortex-induced vibrations of flexible bridges[J]. Journal of Engineering Mechanics, 1990, 116(6):1392-1411.

[2]Gupta H, Sarkar P P, Mehta K C.Identification of vortex- induced-response parameters in time domain[J]. Journal of Engineering Mechanics-Asce, 1996. 122(11):1031-1037.

[3]Marra A M, Mannini C,Bartoli G.Van der Pol-type equation for modeling vortex-induced oscillations of bridge decks[J].Journal of Wind Engineering and Industrial Aerodynamics, 2011,99(6-7):776-785.

[4]周 濤, 朱樂東, 郭震山.經(jīng)驗非線性渦激力模型參數(shù)識別[J]. 振動與沖擊, 2011,30(3):115-118,144.

ZHOU Tao, ZHU Le-dong, GUO Zhen-shan.Parameters identification on nonlinear empirical model for vortex-induced vibration(VIV)[J]. Journal of Vibration and Shock, 2011, 30(3): 115-118,144.

[5]Scanlan R H.State-of-the-art methods for calculating flutter, vortex-induced, and buffeting response of bridge structures [M].Nat.Tech. Information Service,1981.

[6]LIAO Shi-jun.An analytic approximate approach for free oscillations of self-excited systems[J]. International Journal of Non-Linear Mechanics, 2004,39(2):271-280.

[7]Holland J H.Genetic algorithms[J]. Scientific American, 1992, 267(1):66-72.

[8]López J L,López-Ruiz R.Approximating the amplitude andform of limit cycles in the weakly nonlinear regime of Liénardsystems[J]. Chaos, Solitons & Fractals, 2007,34(4):1307-1317.