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

翁祥穎， 葛耀君， 陳海興

(同濟大學 土木工程防災國家重點實驗室，上海 200092)

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

1 單自由度經驗非線性渦激力模型

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

(1)

定義無量綱量：

(3)

式中：m為系統每米模態質量。

結合簡化的渦激力模型得渦振運動平衡微分方程的無量綱形式為：

(4)

其中：

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

(5)

μ2=mrY1(K)ε

(6)

(7)

式中：ξ為系統阻尼比。顯然，式(4)描述的系統振動特性與μ1，μ2密切相關：μ1，μ2符號相反時，系統會發生衰減或發散振動，不符合渦振發生時結構呈等幅振動特性。由渦振為穩定的等幅振動知μ1為正數，故定義變量：

(8)

代入以上無量綱渦振運動平衡微分方程可得標準的Van der Pol振子為：

(9)

2 Van der Pol振子近似解

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

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

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

(10)

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

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

(11)

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

(12)

(13)

(14)

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

3 擬合目標函數

結合風洞試驗實測橋梁節段模型渦振響應時程及所得二階近似解，構造表示殘差平方和的擬合目標函數為：

(15)

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

4 遺傳算法識別參數

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

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

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

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

5 橋梁節段模型渦激力氣動參數識別

5.1 風洞試驗概況

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

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

圖3 試驗系統

圖4 無量綱的豎彎渦振響應時程

5.2 μ1，μ3參數搜索

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

(17)

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

圖5 極限環與擬合限幅環

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

5.3 渦激力氣動參數識別

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

μ2=77732.02

(18)

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

表1 流線型箱梁斷面渦激力氣動參數識別結果

圖6 擬合結果與原始位移時程對比

5.4 識別結果評價

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

6 結 論

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

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

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

參 考 文 獻

[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]. 振動與沖擊, 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.