劉十一，葛耀君

（1.土木工程防災國家重點實驗室（同濟大學），200092上海；2.上海飛機設計研究院，201310上海）

非線性子系統(tǒng)的大振幅時域自激力模型

劉十一1，2，葛耀君1

（1.土木工程防災國家重點實驗室（同濟大學），200092上海；2.上海飛機設計研究院，201310上海）

針對強風作用下的大跨度橋梁或高速飛行中的機翼可能出現(xiàn)大振幅顫振響應，而現(xiàn)有時域自激力模型又無法模擬大振幅下的非線性氣動力問題.提出一種新的非定常大振幅自激氣動力模型及其參數(shù)擬合方法，新模型通過附加非線性微分方程組及附加氣動力自由度來模擬氣動力記憶效應及振幅非線性特性.對于同一主梁斷面，使用一組模型參數(shù)即可模擬不同折算風速和不同振幅下的自激力，模型參數(shù)可通過風洞試驗或CFD數(shù)值模擬結果擬合獲得.結果表明，新模型能再現(xiàn)自激力阻尼特性隨折算風速和振幅的變化，通過單頻振動擬合得到的新模型能再現(xiàn)多頻振動下的非線性自激力時程.

顫振；大振幅自激力；非線性子系統(tǒng)；記憶效應；參數(shù)擬合

模擬大跨度橋梁在強風作用下的顫振過程或顫抖振響應時，需對主梁自激力建模.平板自激力理論解表明，自激力特性與折算風速（或頻率）有關，即包含氣動力記憶效應.Scanlan等最早提出了實用的鈍體橋梁自激力經驗模型［1］，只適用于模擬小振幅情況下的自激力.由于大振幅下的非線性氣動力不滿足頻率疊加原理，必須在時域中計算.現(xiàn)有時域模型中，基于階躍響應函數(shù)的時域卷積［2］和基于有理函數(shù)近似的狀態(tài)空間方法［3-4］均可模擬非定常線性自激力.但目前尚無時域模型能完整再現(xiàn)非定常非線性自激力，一種廣泛采用的非線性氣動力模擬方法是將氣動力分為低頻和高頻分量［5-6］，低頻分量可采用非線性準定常模型計算，高頻分量則采用非定常線性模型計算，最終的總體氣動力為高頻和低頻分量的疊加；另一種改進的高、低頻疊加法［7］采用了流變模型來模擬氣動力的高頻分量.以上方法雖能模擬高頻分量的非定常特性，但無法模擬大振幅產生的非線性.Diana等［8］提出了基于非線性多項式的氣動力模型，吳騰等［9］采用人工神經網絡來模擬氣動力的非線性效應，這兩種方法均能模擬大振幅下的非線性效應，但未模擬氣動力記憶效應，屬準定常非線性模型.本文提出一種新的非定常非線性時域自激力模型.新模型中，記憶效應通過附加自由度來模擬［10］，振幅非線性效應則通過非線性的附加微分方程組（即子系統(tǒng)）及其他非線性表達式來模擬.新模型所使用的微分方程組和表達式均根據(jù)自激力特點而構造，其中包含待擬合的模型參數(shù).對于特定斷面，模型參數(shù)可通過風洞試驗或CFD數(shù)值模擬擬合獲得.本文提出的自激力模型參數(shù)與頻率、振幅無關，能反應自激力特性隨折算風速（或頻率）的變化，能模擬自激力特性隨振幅的非線性變化.

1 自激力模型表達式

本文將均勻來流下的主梁自激力分為3部分：1）靜風力部分；2）準定常部分，用于模擬當前運動對自激力的影響；3）非定常部分，用于模擬氣動力記憶效應，由非線性子系統(tǒng)生成.在新模型中，這3部分均為非線性.將來流風速用向量u表示為

主梁中心處的瞬時相對風攻角θ?和瞬時相對風速大小u?分別為

式中：x、y分別為主梁中心的橫坐標和縱坐標；α為主梁扭轉角，以逆時針為正.

新模型將θ?、u?及其對時間導數(shù)作為輸入變量.為區(qū)分不同的運動產生的氣動力效應，將和分解為

為簡化氣動力表達式，下文使用瞬時相對風速u?、主梁寬度B和空氣密度ρ對所有物理量進行量綱約化，約化后的物理量用下標～表示，約化方程式為

式中fH、fV和M分別為主梁局部坐標系中的阻力、升力和升力矩.

將主梁局部坐標系中的無量綱氣動力floc～表示為

新的自激力模型可表示為

式中：φ為附加氣動力自由度；fst為約化靜風力；fqs為約化準定常自激力部分；fus為約化非定常自激力部分，與φ有關；g=0為附加非線性微分方程組，即非線性子系統(tǒng)，它定義了φ的演化規(guī)律.由式（9）計算可得局部坐標系下的約化氣動力，經坐標變換后可轉化為全局坐標系下的氣動力，即

式（9）中的fst、fqs、fus均為θ?的非線性表達式，g同時為θ?和φ的非線性表達式，靜風力部分fst為

式中a0～a3為與頻率無關的參數(shù)向量.

準定常自激力部分fqs為

式中：rαθ?( )、rmθ?( )、rumθ?( )為瞬時相對攻角θ?的多項式向量，bi、ci、ei（i=0，..，3）為與頻率無關的參數(shù)向量.

非定常自激力部分fus為

式中：Rφθ?( )為多項式矩陣，Ei（i=0，…，3）為與頻率無關的參數(shù)矩陣.

附加非線性微分方程組g=0.

式中：φ3=；Kφ(θ?)、Kφ3(θ?)為多項式對角矩陣；kαθ?( )、kmθ?( )、kumθ?( )為多項式向量；Fi、Gi（i=0，2，4）為與頻率無關的參數(shù)對角矩陣；hi、li、pi（i=0…3）為與頻率無關的參數(shù)向量.

為保證微分方程組（16）中的氣動力自由度φ穩(wěn)定，式（17）中的Kφθ?( )和Kφ3θ?( )只包含θ?的偶數(shù)次項，且這兩個多項式對角矩陣中的元素都大于0.同樣，微分方程組（16）中只包含了φ的奇數(shù)次項.這樣就能保證非線性子系統(tǒng)的自由度φ始終是穩(wěn)定的.

將式（11）、（12）、（14）、（16）代入方程組（9）即可獲得非線性子系統(tǒng)自激力模型的完整表達式.式（11）、（13）、（15）、（17）中的所有參數(shù)均可通過風洞試驗或CFD數(shù)值模擬結果擬合獲得.模型參數(shù)數(shù)量取決于所使用的附加微分方程數(shù)量.

2 模型參數(shù)擬合方法

根據(jù)已知的主梁位移時程和氣動力時程，可擬合獲得所有模型參數(shù).擬合的基本流程為：1）假設一組模型參數(shù)的初始值.2）根據(jù)已知的主梁位移時程和當前的模型參數(shù)值，采用數(shù)值方法求解附加非線性微分方程組，然后計算出氣動力時程.3）根據(jù)計算的氣動力時程和已知的氣動力時程之差，采用一種數(shù)值優(yōu)化算法計算模型參數(shù)的修正量.4）修正模型參數(shù)，如果參數(shù)收斂，則結束，否則轉到第2步.

此方法中最關鍵的是第3步，即如何根據(jù)計算的氣動力時程和已知氣動力時程之差確定模型參數(shù)的修正量.由于模型參數(shù)數(shù)量大（大于100個），且每個迭代步計算量較大（需要求解非線性微分方程組），因此必須使用收斂較快的非線性數(shù)值優(yōu)化算法.本文使用LM算法計算每個迭代步的參數(shù)修正量.LM算法是收斂最快的數(shù)值優(yōu)化算法之一，廣泛用于神經網絡訓練等復雜的非線性優(yōu)化問題.

包含附加微分方程組的氣動力模型可表示為

式中：f t()為模型輸出的氣動力；x t()為輸入的主梁位移；向量pf、pg分別表示f、g包含的模型參數(shù).

定義模型擬合的目標函數(shù)為

指定pf、pg的初始值后，本文使用LM算法優(yōu)化pf、pg以使目標函數(shù)Π取得局部最小值.LM算法需要各離散時刻的氣動力計算值f對模型參數(shù)pf、pg的偏導數(shù)矩陣.而f與氣動力自由度φ有關，故需要各離散時刻的φ對模型參數(shù)pg的偏導數(shù)矩陣.φ為附加非線性微分方程組的解，由于φ沒有解析表達式，本文采用四階龍格庫塔法計算φ的離散時程，而φ對參數(shù)的偏導數(shù)矩陣則采用鏈式法則在數(shù)值積分過程中同時計算.計算偏導數(shù)的基本思路：因為φ的數(shù)值積分過程由基本代數(shù)操作組成，而各代數(shù)操作的偏導數(shù)可直接計算，所以積分獲得的φ的離散時程點的偏導數(shù)可通過鏈式法則間接計算獲得.四階龍格庫塔法的1個積分步由4個子步組成，子步的偏導數(shù)計算流程如圖1所示.

圖1 一個數(shù)值積分子步中的偏導數(shù)計算流程

需要指出，本文給出的模型擬合算法不僅適用于新的非線性自激力模型，還適用于傳統(tǒng)的線性狀態(tài)空間模型，是一種通用的模型擬合算法.由于新的擬合算法尋求的是時域最小二乘解，實際上對原始信號中的高頻隨機噪聲起到了時域濾波的作用，因此對原始信號的噪聲有較強的容忍能力.此算法只能獲得目標函數(shù)的局部最小值，因此擬合結果與參數(shù)初始值有關.當參數(shù)初始值在合理范圍內時，一般都能獲得理想的擬合結果，對式（11）、（13）、（15）、（17）中的模型參數(shù)，本文推薦初始值為

F0∈ 0.2，1[ ]，其他 =0.

3 參數(shù)擬合與模型驗證

本文對一個均勻來流下的箱型主梁斷面進行了CFD強迫振動繞流模擬，并通過表面壓強積分計算出三分力時程.主梁斷面形狀及尺寸如圖2所示. CFD模擬方法為LES，采用Smagorinsky靜態(tài)湍流模型，雷諾數(shù)等于10 000（相對主梁寬度）.CFD模擬了不同攻角下的靜風力時程，及不同折算風速、不同振幅下的強迫振動自激力時程.靜風力包含11個工況，攻角為-10°～10°，步長為2°.單頻扭轉強迫振動包含24個工況，折算風速為3、5、7、9、11、13，振幅為±1°、±3°、±6°、±9°.單頻豎彎強迫振動包含4個工況，折算風速為3、5、7、11，振幅為±0.2H（H為主梁高度）.此外還模擬了兩個多頻扭轉強迫振動工況，用于模型驗證.本文折算風速定義為Ur=U／（Bf）=2πU／（Bω）.本文用CFD模擬獲得的11組靜風力時程和28組單頻強迫振動時程擬合出新模型的一組模型參數(shù).擬合的模型中包含4個附加非線性微分方程組，即有4個附加氣動力自由度，模型共包含140個參數(shù).模型擬合經過約200次迭代達到收斂.圖3～5給出了3種單頻強迫振動下的阻力、升力和升力矩時程，圖中包含CFD計算及新自激力模型擬合結果.需要指出，CFD計算獲得的氣動力時程包含主梁特征紊流產生的高頻分量，無法在自激力模型中模擬.由圖3～5可見，大振幅正弦強迫振動產生的自激力呈現(xiàn)明顯的非正弦形狀，驗證了氣動力非線性的存在.為進一步分析自激力的非線性特性，本文定義等效扭轉氣動阻尼c-α為

由此可計算出任意一組扭轉強迫振動位移和氣動力時程的等效扭轉氣動阻尼.若自激力為線性，等效扭轉阻尼應與振幅無關.圖6給出了CFD模擬得到的等效扭轉氣動阻尼隨折算風速和主梁振幅的變化情況，由圖可見折算風速及振幅對氣動阻尼均有很強影響.對于當前斷面，低折算風速下的氣動阻尼隨振幅增加而增加，但高折算風速下的氣動阻尼隨振幅增加而減小.此外，扭轉振幅為1°和3°時氣動阻尼幾乎相等，由此可見，當振幅小于3°時，當前斷面自激力呈線性.

圖2 箱型主梁斷面形狀及尺寸

圖3 扭轉振幅±9°，折算風速13模型擬合結果

圖4 扭轉振幅±6°，折算風速5時模型擬合結果

圖5 豎彎振幅±0.2H，折算風速3時，模型擬合結果

圖7對比了多組等效扭轉氣動阻尼，分別由CFD模擬、擬合的新模型和擬合的線性狀態(tài)空間模型計算獲得.線性狀態(tài)空間模型使用Minimum-state格式［11］，使用CFD模擬的單頻強迫振動結果擬合.圖中包含不同振幅（3°、6°、9°）下的扭轉阻尼隨折算風速的變化曲線.由于線性狀態(tài)空間模型不能反映振幅對氣動阻尼的影響，因此只有一條曲線.擬合的新模型能較好地反映扭轉氣動阻尼隨折算風速和振幅的變化，模擬了氣動阻尼的非定常和非線性特性.

圖6 折算風速、扭轉振幅對無量綱等效扭轉氣動阻尼的影響

本文用擬合的新模型計算了兩個多頻扭轉強迫振動下的自激力時程，并將其與CFD模擬結果對比，如圖8、9所示.由圖可見擬合的新模型能準確再現(xiàn)多頻大振幅強迫振動產生的自激力時程.因此使用單頻強迫振動試驗擬合非線性自激力模型是可行的.

圖7 不同折算風速、不同振幅下等效扭轉氣動阻尼的擬合結果比較

圖8 包含折算風速分布為7、11，振幅為±3°的扭轉頻率成分的多頻預測結果

圖9 包含折算風速分別為5、7、11，振幅為±3°的扭轉頻率成分的多頻預測結果

4 結 論

1）提出一種新的非定常非線性時域自激力模型及其擬合算法，新模型中，氣動力記憶效應由附加氣動力自由度模擬，大振幅非線性效應由與瞬時相對風攻角和氣動力自由度有關的非線性表達式模擬.新模型包含一組與頻率、振幅無關的參數(shù)，這些參數(shù)可通過已知的位移和氣動力時程擬合獲得，擬合過程基于LM數(shù)值優(yōu)化算法，并使用鏈式法則計算微分方程數(shù)值解的偏導數(shù)矩陣.

2）使用CFD強迫振動模擬獲得了一個橋梁斷面在不同折算風速和不同振幅下的氣動力時程，并擬合了新模型的一組參數(shù).針對CFD模擬結果，新模型能較好地擬合不同折算風速和振幅下的非線性氣動力時程，并能準確反映扭轉氣動阻尼隨折算風速和振幅的變化.使用單頻強迫振動擬合的新模型能準確再現(xiàn)多頻強迫振動產生的氣動力時程.

3）由于CFD模擬與真實氣動力之間存在一定差別，新模型還有待通過風洞試驗進一步校核.本研究僅限于均勻來流下的自激力非線性問題，真實橋梁風振還受來流紊流產生的抖振力影響，抖振力非線性及其與自激力的相互影響還有待進一步研究.

［1］GE Y J，TANAKA H.Aerodynamic flutter analysis of cable-supported bridges by multi-modeand full-mode approaches［J］.JournalofWind Engineering and Industrial Aerodynamics，2000，86：123-153.

［2］COSTA C，BORRI C，F(xiàn)LAMAND O，et al.Time-domain buffeting simulations for wind-bridge interaction ［J］. Journal of Wind Engineering and Industrial Aerodynamics，2007，95：991-1006.

［3］CHEN X，KAREEM A.Aeroelastic analysis of bridges：effects of turbulence and aerodynamic nonlinearities［J］. Journal of Engineering Mechanics，2003，129：885-895.

［4］CHEN X Z，MATSUMOTO M，KAREEM A.Time domain flutter and buffeting response analysis of bridges［J］. Journal of Engineering Mechanics-Asce，2000，126：7-16.

［5］CHEN X Z，KAREEM A.Nonlinear response analysis of long-span bridges under turbulent winds［J］.Journal of Wind Engineering and Industrial Aerodynamics，2001，89：1335-1350.

［6］DIANA G，F(xiàn)ALCO M，BRUNI S，et al.Comparisons between wind tunnel tests on a full aeroelastic model of the proposed bridge over Stretto di Messina and numerical results［J］.Journal of Wind Engineering and Industrial Aerodynamics，1995，54-55：101-113.

［7］DIANA G，ROCCHI D，ARGENTINI T.An experimental validation of a band superposition model of the aerodynamic forces acting on multi-box deck sections［J］. Journal of Wind Engineering and Industrial Aerodynamics，2013，113：40-58.

［8］DIANA G，ROCCHI D，ARGENTINI T，et al.Aerodynamic instability of a bridge deck section model： linear and nonlinear approach to force modeling［J］.Journal of Wind Engineering and Industrial Aerodynamics，2010，98：363-374.

［9］WU T，KAREEM A.Modeling hysteretic nonlinear behavior of bridge aerodynamics via cellular automata nested neural network［J］.Journal of Wind Engineering and Industrial Aerodynamics，2011，99：378-388.

［10］LIU S，GE Y.A new unsteady aerodynamic model with unique parameters［C］／／Proceedings of BBAA7.Shanghai：China Communication Press，2012.

［11］WILDE K，OMENZETTER P，F(xiàn)UJINO Y.Suppression of bridge flutter by active deck-flaps control system ［J］. Journal of Engineering Mechanics，2001，127：80-89.

（編輯 魏希柱）

Nonlinear dynamic subsystem model for large-amplitude motion-induced aerodynamic forces of bridge decks

LIU Shiyi1，2，GE Yaojun1

（1.State Key Laboratory of Disaster Reduction in Civil Engineering（Tongji University），200092 Shanghai，China；2.Shanghai Aircraft Design and Research Institute，201310 Shanghai，China）

In simulation of long-span bridge and airfoil flutter，a large-amplitude motion-induced force model is required.However，existing time-domain motion-induced force models cannot simulate memory effects and amplitude dependency simultaneously.This paper proposes a new time-domain non-linear motion-induced force model，which employs a set of nonlinear differential equations and augmented aerodynamic degrees of freedoms to simulate memory effects and amplitude dependency.Motion-induced forces under different deck motion amplitudes，different reduced wind velocities and non-sinusoidal motions can be simulated using a single set of model parameters.A fitting algorithm for the new model is also proposed and the model parameters of a bridge deck are fitted using results obtained by CFD simulations.Numerical example shows that the new model can reproduce variation of aerodynamic damping characteristics with reduced wind velocity and amplitude.A model fitted using single-frequency forced vibration tests can reproduce forces generated by multi-frequency vibrations.

flutter；large-amplitude motion-induced force；nonlinear dynamic subsystem；memory effects；model fitting

U441.3

A

0367-6234（2015）09-0073-06

10.11918／j.issn.0367-6234.2015.09.014

2014-01-20.

國家自然科學基金（91215302）.

劉十一（1986—），男，工程師，博士研究生；葛耀君（1958—），男，教授，博士生導師.

劉十一，liushiyi_ha＠163.com.