賈明曉,劉祖軍,,楊詠昕
(1.華北水利水電學(xué)院 土木與交通學(xué)院,鄭州 450011;2.同濟(jì)大學(xué) 橋梁工程系,上海 200092)
數(shù)值方法與計(jì)算機(jī)技術(shù)進(jìn)步提升了對(duì)顫振本質(zhì)的認(rèn)識(shí),已能較準(zhǔn)確預(yù)測(cè)橋梁顫振的臨界風(fēng)速、顫振頻率及顫振形態(tài)。顫振屬自激振動(dòng),其物理機(jī)理可從能量角度進(jìn)行闡釋。處于氣流中橋梁的能量反饋機(jī)制表現(xiàn)為氣流輸入結(jié)構(gòu)—?dú)饬飨到y(tǒng)中能量與結(jié)構(gòu)阻尼耗散能量之間的平衡關(guān)系。當(dāng)輸入結(jié)構(gòu)—?dú)饬飨到y(tǒng)中能量小于結(jié)構(gòu)阻尼耗能時(shí),結(jié)構(gòu)在初始擾動(dòng)下將作衰減(阻尼)振動(dòng);而當(dāng)輸入能量大于結(jié)構(gòu)阻尼耗能時(shí)結(jié)構(gòu)在初始擾動(dòng)下將作發(fā)散振動(dòng);兩者相等時(shí)結(jié)構(gòu)在初始擾動(dòng)下將作等幅簡(jiǎn)諧振動(dòng)。
對(duì)橋梁自激振動(dòng)能量分析,Scanlan[1]最早提出橋梁顫振的多模態(tài)分析方法,并從能量觀點(diǎn)對(duì)橋梁顫振穩(wěn)定性進(jìn)行研究,給出一個(gè)振動(dòng)周期內(nèi)氣流沿橋梁斷面每延米輸入的總能量與結(jié)構(gòu)耗能表達(dá)式,但此僅為理論框架,如何從能量角度對(duì)橋梁進(jìn)行顫振分析,并未給出具體方法。Larsen[2]以CFD方法為基礎(chǔ),據(jù)離散渦計(jì)算中渦旋運(yùn)動(dòng)規(guī)律提出簡(jiǎn)化分析模型。該模型描述橋梁斷面扭轉(zhuǎn)運(yùn)動(dòng)一個(gè)周期內(nèi)渦旋的運(yùn)動(dòng)情況,并通過積分估算由渦旋產(chǎn)生的氣動(dòng)力對(duì)橋梁斷面所做總功,用能量方法分析渦激力做功與結(jié)構(gòu)穩(wěn)定之關(guān)系。文獻(xiàn)[2]頗具開創(chuàng)性,但其計(jì)算模型假定在旋渦沿橫截面移動(dòng)時(shí),旋渦升力保持不變且結(jié)構(gòu)阻尼為零,此與實(shí)際情況有很大不同。劉高[3]從結(jié)構(gòu)與氣流系統(tǒng)內(nèi)部能量平衡觀點(diǎn)對(duì)系統(tǒng)顫振進(jìn)行研究,發(fā)展一種全橋多模態(tài)顫振分析方法-能量法,通過建立系統(tǒng)等效阻尼比與系統(tǒng)能量變化率之間的關(guān)系,推導(dǎo)出系統(tǒng)及各階模態(tài)等效阻尼比計(jì)算方法,并據(jù)不同風(fēng)速下系統(tǒng)能量變化率判斷系統(tǒng)顫振穩(wěn)定性。但該方法不能分析氣流與結(jié)構(gòu)間的細(xì)觀作用,尤其結(jié)構(gòu)自激振動(dòng)中系統(tǒng)能量分布及轉(zhuǎn)化情況。
在機(jī)翼氣動(dòng)彈性力學(xué)研究中,F(xiàn)razer[4]分析了維持機(jī)翼強(qiáng)迫振動(dòng)所需的能量條件;Nissim[5]提出氣動(dòng)能量概念作為對(duì)機(jī)翼主動(dòng)控制的理論基礎(chǔ);Jones[6]分析顫振發(fā)生的能量特征;Carta[7]提出用于顫振失穩(wěn)預(yù)測(cè)的能量穩(wěn)定判據(jù);Klose等[8]通過研究證明能量法用于顫振失穩(wěn)的預(yù)測(cè)被認(rèn)為特征值法在葉片高質(zhì)量比下的特殊應(yīng)用,并指出用能量法預(yù)測(cè)的失效情況。
本文通過風(fēng)洞試驗(yàn)研究平板的流線斷面顫振性能,基于流固松耦合的計(jì)算策略,利用現(xiàn)有流體軟件的用戶自定義(UDF)功能及CFD數(shù)值方法,模擬平板顫振過程,據(jù)分塊分析思想研究顫振過程中振動(dòng)模型表面不同區(qū)域吸收氣流能量特點(diǎn),分析旋渦非定常演化對(duì)模型表面不同區(qū)域壓力特性影響。
1.1 平板風(fēng)洞測(cè)振試驗(yàn)及結(jié)果分析
風(fēng)洞試驗(yàn)在同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室TJ-4邊界層風(fēng)洞中進(jìn)行,該風(fēng)洞為低速回流式風(fēng)洞,測(cè)振試驗(yàn)段尺寸為:寬 0.814 m,高 0.8 m,長(zhǎng) 2.0 m,設(shè)計(jì)最大試驗(yàn)風(fēng)速30 m/s。平板斷面風(fēng)洞測(cè)振采用模型如圖1所示,縱向長(zhǎng)度0.8 m。為保證輕質(zhì)高強(qiáng),材料采用碳纖維材質(zhì)。
圖1 平板模型(單位:mm)Fig.1 The model of plate(unit:mm)
模型的基本參數(shù):m=1.525 kg/m,Im=0.010 28 kg/m,豎向頻率fh=2.588 Hz,扭轉(zhuǎn)頻率fa=5.026 Hz。用激光位移計(jì)記錄模型試驗(yàn)風(fēng)速下位移信號(hào)。
在平板測(cè)振試驗(yàn)中,記錄的各級(jí)風(fēng)速下結(jié)構(gòu)位移響應(yīng)結(jié)果見圖2、圖3。風(fēng)速小于18 m/s時(shí)結(jié)構(gòu)位移振幅較小,基本處于靜止?fàn)顟B(tài),風(fēng)速超過18 m/s后結(jié)構(gòu)振動(dòng)位移值及方差突然增大,豎向振動(dòng)參與程度不斷增加,直之風(fēng)速達(dá)20.4 m/s時(shí)結(jié)構(gòu)出現(xiàn)顫振失穩(wěn)現(xiàn)象。
圖2 豎向位移隨風(fēng)速變化關(guān)系Fig.2 Vertical vibration displacements vs wind velocity curve
圖3 扭轉(zhuǎn)振動(dòng)位移響應(yīng)-風(fēng)速曲線Fig.3 Torsional vibration displacement vs wind velocity curve of plate section
圖4 位移幅值譜(20.4 m/s)Fig.4 Displacement amplitude spectrum(20.4 m/s)
氣動(dòng)彈性問題具有重要研究?jī)r(jià)值,當(dāng)橋梁結(jié)構(gòu)與空氣存在相對(duì)運(yùn)動(dòng)時(shí),橋梁結(jié)構(gòu)會(huì)受到氣動(dòng)力作用,氣動(dòng)力會(huì)導(dǎo)致結(jié)構(gòu)位置及形狀發(fā)生改變。結(jié)構(gòu)位形的改變引起流體占據(jù)實(shí)際空間發(fā)生變化,反過來(lái)又導(dǎo)致作用在結(jié)構(gòu)上氣動(dòng)力改變形成流體與固體的相互耦合作用。橋梁的顫振發(fā)散也是流固耦合作用的結(jié)果。
處理流固耦合問題可采用強(qiáng)耦合與松耦合兩種方法。強(qiáng)耦合求解方法以求解流體運(yùn)動(dòng)的迭代過程為主,將物體運(yùn)動(dòng)的求解耦合在流體運(yùn)動(dòng)求解中;松耦合方法為在每個(gè)時(shí)間步內(nèi)分別對(duì)流體與結(jié)構(gòu)依次求解,通過中間平臺(tái)交換數(shù)據(jù)信息,實(shí)現(xiàn)兩個(gè)場(chǎng)的耦合求解。此方法對(duì)計(jì)算機(jī)硬件要求不高,且可充分發(fā)揮CFD(Computational Fluid Dynamics)計(jì)算與CSD(Computational Structural Dynamics)計(jì)算的各自優(yōu)點(diǎn),能實(shí)現(xiàn)數(shù)值計(jì)算模塊化[9]。分析橋梁結(jié)構(gòu)風(fēng)致振動(dòng)問題時(shí)多采用松耦合的計(jì)算策略。
本文采用松耦合計(jì)算方法,用現(xiàn)有商業(yè)軟件fluent提供的基于ALE(Arbitrary Lagrange Euler)方法的動(dòng)網(wǎng)格技術(shù),實(shí)現(xiàn)顫振數(shù)值模擬,通過軟件用戶自定義函數(shù)描述模型及周圍氣流的剛體運(yùn)動(dòng),顫振數(shù)值模擬過程具體步驟為:
(1)進(jìn)行流體域內(nèi)固定斷面的非定常繞流計(jì)算,提取靜止?fàn)顟B(tài)下初始流場(chǎng)及作用在斷面上的三分力,獲得t時(shí)刻氣動(dòng)力。
(2)據(jù)計(jì)算所得升力、升力矩,通過數(shù)值方法對(duì)振動(dòng)結(jié)構(gòu)進(jìn)行結(jié)構(gòu)動(dòng)力學(xué)求解。據(jù)newmarkk-β方法求解豎向及扭轉(zhuǎn)振動(dòng)方程,獲得t+Δt時(shí)刻結(jié)構(gòu)運(yùn)動(dòng)的豎向速度與扭轉(zhuǎn)速度:
(3)將模型振動(dòng)速度賦值給包圍模型斷面并隨模型一起做剛體運(yùn)動(dòng)部分流體區(qū)域,通過fluent提供的用法自定以函數(shù)(UDF)描述模型及周圍流體運(yùn)動(dòng),采用動(dòng)網(wǎng)格方法進(jìn)行流體域求解,獲得振動(dòng)斷面氣動(dòng)升力及升力矩。
(4)重復(fù)(2)、(3)步,獲得橋梁斷面振動(dòng)響應(yīng)時(shí)程。
(5)據(jù)斷面響應(yīng)時(shí)程曲線判斷振動(dòng)是否已發(fā)散,若已發(fā)散則該風(fēng)速即為顫振臨界風(fēng)速,否則增加風(fēng)速按(1)~(4)步重新計(jì)算直到振動(dòng)發(fā)散。
具體計(jì)算流程如圖5所示。
數(shù)值模擬用商業(yè)軟件Fluent提供的RANS方法的k-ωSST兩方程模型,計(jì)算域大小參考同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)試驗(yàn)室TJ-4風(fēng)洞中段試驗(yàn)端設(shè)置,計(jì)算域沿流線長(zhǎng)度為3 m(上游1 m,下游2 m),橫向?qū)挾?.8 m。計(jì)算時(shí)壁面附近最小網(wǎng)格尺度為0.000 4 m,計(jì)算域用分塊結(jié)構(gòu)化網(wǎng)格。網(wǎng)格數(shù)量16.2萬(wàn)。
計(jì)算參數(shù)設(shè)置:動(dòng)量、湍動(dòng)能、能量耗散均采用兩階迎風(fēng)格式進(jìn)行離散,壓力速度耦合采用SIMPLE算法,求解器采用分離式,計(jì)算模式選用兩階隱式。邊界條件設(shè)定:速度入口,湍流強(qiáng)度0.5%,壓力出口,計(jì)算域上、下端設(shè)為對(duì)稱邊界條件,表面采用無(wú)滑移的壁面條件。
采用數(shù)值方法對(duì)顫振過程進(jìn)行模擬及預(yù)測(cè),數(shù)值模擬獲得該平板顫振臨界風(fēng)速為21.2 m/s與試驗(yàn)結(jié)果20.4 m/s尚有一定差距。圖6~圖8為平板運(yùn)動(dòng)時(shí)程及幅值譜,數(shù)值模擬所得顫振頻率為4.73 Hz,實(shí)測(cè)顫振頻率為 4.59 Hz。
圖5 顫振數(shù)值模擬流程圖Fig.5 The process of flutter numerical simulation
圖6 豎向振動(dòng)位移時(shí)程響應(yīng)(20.4 m/s,CFD數(shù)值模擬)Fig.6 Vertical vibration displacement vs time curve of plate section(20.5 m/s,CFD)
圖7 扭轉(zhuǎn)振動(dòng)位移時(shí)程響應(yīng)(20.4 m/s,CFD 數(shù)值模擬)Fig.7 Torsional vibration displacement vs time curve of plate section(20.4 m/s,CFD)
圖8 位移幅值譜(20.4 m/s)Fig.8 Displacement amplitude spectrum(20.4 m/s)
分塊分析主要為詳細(xì)捕捉表面壓力變化,尋找引起顫振發(fā)散的主要能量輸入特點(diǎn),考慮尾部旋渦對(duì)氣流能量輸入方式影響及影響范圍。對(duì)斷面分區(qū),分別計(jì)算不同區(qū)域的能量輸入特點(diǎn)。單位區(qū)域能量為單位面積上一個(gè)周期內(nèi)氣動(dòng)力輸入到系統(tǒng)的能量,其計(jì)算式為:
面積為S(圖9的網(wǎng)格線部分)的區(qū)域氣動(dòng)力輸入能量w(s,t)計(jì)算式為:
式中:P(x,t)為結(jié)構(gòu)表面壓力,v(x,t)為結(jié)構(gòu)運(yùn)動(dòng)速度,n(x,t)為壓力與速度夾角。
圖9 分塊分析示意圖Fig.9 The block analysis
據(jù)本文提出的分塊分析思路編制計(jì)算程序,通過提取數(shù)值計(jì)算獲得的模型表面壓力,詳細(xì)研究一個(gè)周期內(nèi)氣流向結(jié)構(gòu)輸送能量過程,對(duì)模型表面進(jìn)行分區(qū)見圖10。
圖10 模型表面分區(qū)Fig.10 The partition of the model
圖11~圖14為平板扭轉(zhuǎn)振動(dòng)時(shí)不同分區(qū)輸入到系統(tǒng)的能量。第1區(qū)域在前半個(gè)周期內(nèi)上表面消耗能量,下表面吸入能量,后半個(gè)周期則反之。一個(gè)完整周期內(nèi)輸入到系統(tǒng)的能量為正,但數(shù)值較小。從上下表面輸入到系統(tǒng)能量隨時(shí)間變化關(guān)系曲線上可看出,輸入的能量周期性明顯,如圖11所示。第4區(qū)域輸入系統(tǒng)的能量特點(diǎn)與1區(qū)基本相似,能量輸入也具有周期性特點(diǎn),二者不同之處在于第4區(qū)域消耗系統(tǒng)能量。由于第4區(qū)及第1區(qū)靠近尾渦,因此受到尾部旋渦的周期性運(yùn)動(dòng)影響,該兩區(qū)域的能量輸入具有明顯的周期性特征(圖14)。
第2區(qū)域?yàn)橄到y(tǒng)提供了較多能量,在整個(gè)振動(dòng)過程中下表面始終輸入正能量,上表面在前1/4周期消耗能量,后3/4周期內(nèi)吸入能量,總能量為正,數(shù)值較大,且增加較快,見圖12。第3區(qū)域即風(fēng)嘴部位為扭轉(zhuǎn)系統(tǒng)最主要能量來(lái)源,在一個(gè)完整周期內(nèi)上下表面均吸入系統(tǒng)能量(圖13),且該區(qū)域吸入的能量遠(yuǎn)大于其它區(qū)域?qū)ο到y(tǒng)能量的貢獻(xiàn),為引起結(jié)構(gòu)振動(dòng)發(fā)散的主要能量源。從第2、3區(qū)域能量隨時(shí)間變化關(guān)系曲線上可看出,該兩部位能量輸入周期性不明顯,受尾部旋渦運(yùn)動(dòng)影響較小。
豎向振動(dòng)主要能量輸入則主要由第2區(qū)域引起,見圖15,各分區(qū)對(duì)扭轉(zhuǎn)振動(dòng)能量的貢獻(xiàn)如圖16所示。通過上述分析知,風(fēng)嘴為扭轉(zhuǎn)振動(dòng)的主要能量輸入部位,對(duì)系統(tǒng)扭轉(zhuǎn)振動(dòng)穩(wěn)定性起重要作用,見圖17。
圖11 扭轉(zhuǎn)運(yùn)動(dòng)1區(qū)能量隨時(shí)間變化關(guān)系Fig.11 The energy of torsion in the first partition vs.time curve
圖12 扭轉(zhuǎn)運(yùn)動(dòng)2區(qū)能量隨時(shí)間變化關(guān)系Fig.12 The energy of torsion in the second partition vs.time curve
圖13 扭轉(zhuǎn)運(yùn)動(dòng)3區(qū)能量隨時(shí)間變化關(guān)系Fig.13 The energy of torsion in the third partition vs.time curve
圖14 扭轉(zhuǎn)運(yùn)動(dòng)4區(qū)能量隨時(shí)間變化關(guān)系Fig.14 The energy of torsion in the fourth partition vs.time curve
圖15 各分區(qū)扭轉(zhuǎn)振動(dòng)能量隨時(shí)間變化關(guān)系Fig.15 The energy of torsion in every partition vs.time curve
圖16 各分區(qū)豎向振動(dòng)能量隨時(shí)間變化關(guān)系Fig.16 11 The energy of vertical in every partition vs.time curve
圖17 氣流能量輸入圖Fig.17 The energy input by air
通過對(duì)CFD數(shù)值模擬計(jì)算結(jié)果研究分析發(fā)現(xiàn)平板流線型好的斷面周圍旋渦較難識(shí)別,由于振動(dòng)模型周圍流場(chǎng)千變?nèi)f化,即使模型振動(dòng)到同一位置,其周圍流場(chǎng)也完全不同。因此先用相位平均方法[10]得出模型振動(dòng)4分點(diǎn)相位的流場(chǎng)圖,后用Snapshot POD旋渦識(shí)別方法[11]重建流場(chǎng),所得流場(chǎng)特征見圖18。
由數(shù)值模擬結(jié)果看出,平板因流線型較好,振動(dòng)模型周圍流場(chǎng)中旋渦識(shí)別較困難。據(jù)POD流場(chǎng)分解結(jié)果,平板處于顫振臨界狀態(tài)時(shí),其尾部旋渦形狀為較規(guī)則圓形,旋渦渦心位于同一高度,呈直線排列,渦街結(jié)構(gòu)較穩(wěn)定,隨著模型振動(dòng)而發(fā)生擺動(dòng),但渦心連線基本保持在同一條直線上。
由4區(qū)域輸入至系統(tǒng)扭轉(zhuǎn)振動(dòng)能量隨時(shí)間的變化關(guān)系可推測(cè)出尾部旋渦對(duì)氣流輸入系統(tǒng)能量的影響??拷擦鲄^(qū)域受尾部旋渦運(yùn)動(dòng)規(guī)律作用,輸入系統(tǒng)的能量周期性特點(diǎn)明顯。為更詳細(xì)研究尾流運(yùn)動(dòng)對(duì)結(jié)構(gòu)振動(dòng)的影響,圖19給出4區(qū)域上下表面壓力分布隨時(shí)間的變化關(guān)系。為清楚說(shuō)明尾部旋渦對(duì)模型不同部位影響,定義上下表面壓力差的相對(duì)變化量為:=(CPu-CPd)/CPu,其中CPu,CPd為模型上下表面壓力最大值。4區(qū)域上下表面的相對(duì)壓力差為:(4)=2.21,(1)=1.61,(2)=1.11,(3)=0(圖20),由此看出,越靠近尾流區(qū)域,上下表面壓力差值的相對(duì)值越大,第4分區(qū)最大,而距尾流區(qū)最遠(yuǎn)的第3分區(qū)上下表面之間的壓力相對(duì)差較小,兩曲線基本重合。即越靠近尾流區(qū)域,其表面壓力變化特點(diǎn)與尾流旋渦運(yùn)動(dòng)規(guī)律越接近。平板振動(dòng)時(shí)尾部旋渦成直線排列并隨平板一起做大幅振動(dòng),會(huì)造成接近尾流區(qū)域的表面壓力相對(duì)變化較大,而遠(yuǎn)離尾流區(qū)域影響較小。
圖18 平板斷面旋渦驅(qū)動(dòng)結(jié)構(gòu)運(yùn)動(dòng)流(顫振風(fēng)速20.4 m/s,數(shù)值模擬)Fig.18 Vortex driven structural movement of plate section(20.4m/s,CFD)
圖19 平板表面4個(gè)分區(qū)上下表面壓力隨時(shí)間變化關(guān)系Fig.19 Four regional of plate girder surface pressure vs.time
圖20 平板表面4個(gè)分區(qū)上下表面相對(duì)壓力差波動(dòng)系數(shù)Fig.20 The volatility coefficient of relative pressure difference on the four model surface regional
本文通過風(fēng)洞試驗(yàn)測(cè)試了流線型較好平板斷面的顫振性能,采用CFD數(shù)值計(jì)算方法,并基于流固松耦合計(jì)算策略模擬了平板顫振過程,結(jié)合分塊分析思路定量研究氣流對(duì)振動(dòng)模型表面不同區(qū)域的能量輸入特點(diǎn)及其對(duì)模型不同部位表面壓力影響,結(jié)論如下:
(1)分塊分析計(jì)算結(jié)果表明迎風(fēng)端風(fēng)嘴是扭轉(zhuǎn)振動(dòng)能量的主要輸入部位,對(duì)系統(tǒng)扭轉(zhuǎn)振動(dòng)穩(wěn)定性起重要作用;氣流輸入到平板的主要能量在一個(gè)周期內(nèi)不斷增加,導(dǎo)致平板顫振多為突然性失穩(wěn)。
(2)由分塊分析結(jié)果知,模型尾部風(fēng)嘴處旋渦對(duì)氣流能量輸入方式影響較大。平板處于顫振狀態(tài)時(shí),尾部受旋渦的直接作用。尾端風(fēng)嘴處穩(wěn)定的渦街可控制模型尾端自由振動(dòng),使模型尾端運(yùn)動(dòng)與渦街?jǐn)[動(dòng)趨勢(shì)基本同步;尾部渦街?jǐn)[動(dòng)的規(guī)律性,造成模型靠近尾流區(qū)域能量輸入方式周期性明顯;而迎風(fēng)端的規(guī)律性較差,周期性不太明顯。
(3)分塊分析的結(jié)果表明平板尾部規(guī)律性擺動(dòng)的渦街控制結(jié)構(gòu)運(yùn)動(dòng),造成尾部風(fēng)嘴區(qū)域上下表面壓力相對(duì)波動(dòng)較大,而遠(yuǎn)離尾流區(qū)域影響較小,其上下表面壓力波動(dòng)程度相對(duì)較小。
(4)平板顫振過程中在模型尾部產(chǎn)生的旋渦尺度較小,形態(tài)基本為圓形,氣流的能量分布在該尺度相當(dāng)?shù)男郎u之間,當(dāng)模型振動(dòng)到平衡位置時(shí),渦街位于同一高度,然后隨著模型的振動(dòng)渦街發(fā)生相應(yīng)擺動(dòng),而渦心連線基本保持直線。
[1]Scanlan R H.The action of flexible bridges under wind,Ⅰ:flutter theory[J].J.of Sound and Vibration,1978,60(2):187-199.
[2] Larsen A.Aerodynamics of the tacoma narrows bridgc-60 years later[J].Joumd of Structural Engineering International,2000,10(4):243-248.
[3]劉 高.大跨度懸索橋顫振分析的能量法[J].中國(guó)公路學(xué)報(bào),2000,12(3):20-24.
LIU Gao.Flutter analysis of long-span suspension bridges by energy method[J].China Journal of Highway and Transport,2000,12(3):20-24.
[4]Frazer R A.On the power input required to maintain forced oscillations of an aeroplane wing in flight[M].ARCR&M,1939.
[5]Nissim E.Flutter suppression using active controls based on the concept of aerodynamic energy[M].NASA TND 6199,1997.
[6]Jones J G.On the energy characteristics of the aerodynamic matrix and the relationship to possible flutter[R].The Aeronautical Quarterly,Part 3,1983:212-225.
[7] Carta F O.Coupled blade-disk-shroud flutter instabilities in turbojet engine rotors[J].Journal of Engineering for Power,1967,89(3):419-426.
[8]Klose A,Heinig K.A comparison of flutter calculations based on Eigenvalue and energy method[R].AD-A2119741,1989.
[9] Lubcke H,Schmidt S,Rung T,et al.Comparison of les and rans in bluff-body flows[J].Journal of Wind Engineering and Industrial Aerodynamics,2001,89(14-15):1471-1485.
[10] Fujisawa N,Takeda G,Ike N.Phase-averaged characteristics of flow around a circular cylinder under acoustic excitation control[J].Journal of Fluids and Structures,2004,19(2):159-170.
[11] Berkooz G,Holmes P,Lumley J L.Proper orthogonal decomposition in the analysis of turbulent flows[J].Annual Review of Fluid Mechanics,1993,25:537-550.