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

劉十一，葛耀君

（1.土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室（同濟(jì)大學(xué)），200092上海；2.上海飛機(jī)設(shè)計(jì)研究院，201310上海）

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

劉十一1，2，葛耀君1

（1.土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室（同濟(jì)大學(xué)），200092上海；2.上海飛機(jī)設(shè)計(jì)研究院，201310上海）

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

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

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

1 自激力模型表達(dá)式

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

主梁中心處的瞬時(shí)相對(duì)風(fēng)攻角θ?和瞬時(shí)相對(duì)風(fēng)速大小u?分別為

式中：x、y分別為主梁中心的橫坐標(biāo)和縱坐標(biāo)；α為主梁扭轉(zhuǎn)角，以逆時(shí)針為正.

新模型將θ?、u?及其對(duì)時(shí)間導(dǎo)數(shù)作為輸入變量.為區(qū)分不同的運(yùn)動(dòng)產(chǎn)生的氣動(dòng)力效應(yīng)，將和分解為

為簡化氣動(dòng)力表達(dá)式，下文使用瞬時(shí)相對(duì)風(fēng)速u?、主梁寬度B和空氣密度ρ對(duì)所有物理量進(jìn)行量綱約化，約化后的物理量用下標(biāo)～表示，約化方程式為

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

將主梁局部坐標(biāo)系中的無量綱氣動(dòng)力floc～表示為

新的自激力模型可表示為

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

式（9）中的fst、fqs、fus均為θ?的非線性表達(dá)式，g同時(shí)為θ?和φ的非線性表達(dá)式，靜風(fēng)力部分fst為

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

準(zhǔn)定常自激力部分fqs為

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

非定常自激力部分fus為

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

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

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

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

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

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

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

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

包含附加微分方程組的氣動(dòng)力模型可表示為

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

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

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

圖1 一個(gè)數(shù)值積分子步中的偏導(dǎo)數(shù)計(jì)算流程

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

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

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

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

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

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

圖3 扭轉(zhuǎn)振幅±9°，折算風(fēng)速13模型擬合結(jié)果

圖4 扭轉(zhuǎn)振幅±6°，折算風(fēng)速5時(shí)模型擬合結(jié)果

圖5 豎彎振幅±0.2H，折算風(fēng)速3時(shí)，模型擬合結(jié)果

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

圖6 折算風(fēng)速、扭轉(zhuǎn)振幅對(duì)無量綱等效扭轉(zhuǎn)氣動(dòng)阻尼的影響

本文用擬合的新模型計(jì)算了兩個(gè)多頻扭轉(zhuǎn)強(qiáng)迫振動(dòng)下的自激力時(shí)程，并將其與CFD模擬結(jié)果對(duì)比，如圖8、9所示.由圖可見擬合的新模型能準(zhǔn)確再現(xiàn)多頻大振幅強(qiáng)迫振動(dòng)產(chǎn)生的自激力時(shí)程.因此使用單頻強(qiáng)迫振動(dòng)試驗(yàn)擬合非線性自激力模型是可行的.

圖7 不同折算風(fēng)速、不同振幅下等效扭轉(zhuǎn)氣動(dòng)阻尼的擬合結(jié)果比較

圖8 包含折算風(fēng)速分布為7、11，振幅為±3°的扭轉(zhuǎn)頻率成分的多頻預(yù)測結(jié)果

圖9 包含折算風(fēng)速分別為5、7、11，振幅為±3°的扭轉(zhuǎn)頻率成分的多頻預(yù)測結(jié)果

4 結(jié) 論

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

2）使用CFD強(qiáng)迫振動(dòng)模擬獲得了一個(gè)橋梁斷面在不同折算風(fēng)速和不同振幅下的氣動(dòng)力時(shí)程，并擬合了新模型的一組參數(shù).針對(duì)CFD模擬結(jié)果，新模型能較好地?cái)M合不同折算風(fēng)速和振幅下的非線性氣動(dòng)力時(shí)程，并能準(zhǔn)確反映扭轉(zhuǎn)氣動(dòng)阻尼隨折算風(fēng)速和振幅的變化.使用單頻強(qiáng)迫振動(dòng)擬合的新模型能準(zhǔn)確再現(xiàn)多頻強(qiáng)迫振動(dòng)產(chǎn)生的氣動(dòng)力時(shí)程.

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

［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.

國家自然科學(xué)基金（91215302）.

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

劉十一，liushiyi_ha＠163.com.