基于CFD和CSD耦合的渦激振和顫振氣彈模擬

李永樂，朱佳琪，2，唐浩俊

（1．西南交通大學橋梁工程系，成都 610031；2．嘉興市交通工程質量安全監督站，浙江 嘉興 314001）

基于CFD和CSD耦合的渦激振和顫振氣彈模擬

李永樂1，朱佳琪1，2，唐浩俊1

（1．西南交通大學橋梁工程系，成都 610031；2．嘉興市交通工程質量安全監督站，浙江 嘉興 314001）

以FLUENT為研究工具，利用微分方程的數值解法和動網格技術，基于松耦合方法將Newmark算法通過UDF嵌入Fluent軟件中，實現了CFD和CSD耦合的分析方法。通過建立二維方柱繞流模型，計算了豎向單自由度振動方柱在不同風速下的斯托羅哈數和最大振幅的變化情況，模擬了渦激共振鎖定現象，并與靜態繞流的結果進行了對比。建立了具有豎向振動和扭轉振動二自由度的薄平板模型，并識別了該平板的顫振導數，進一步對其彎扭耦合顫振臨界風速進行了逼近計算，本方法得到的顫振臨界風速與Scanlan理論公式和Selberg理論公式吻合較好。

CFD和CSD耦合分析；動網格；渦激振；顫振；氣彈效應

空氣與彈性體的相互作用是一個復雜的流固耦合作用過程。當彈性結構在空氣中運動或受空氣流動的影響時，由于空氣動力的作用，彈性體會發生變形或振動，同時會擾亂流場，氣動力又隨之發生變化，從而相互耦合。這種耦合求解過程除涉及到計算流體動力學（Computational Fluid Dynamics，CFD）外，還與計算結構動力學（Computational Structural Dynamics，CSD）聯系緊密。

20世紀90年代，數值風洞的研究得到迅速的發展，在圓柱、方柱繞流數值模擬方面，國際上先后發展了LES方法、基于LES的有限單元法、基于LES的有限差分法等，Schulz等［1］利用不可壓有限體積程序和可變形的混合網格技術、Chio等［2］開發了基于自適應網格的有限元程序對更高雷諾數流場進行了CFD研究。這些研究多數僅針對氣動繞流場，未能考慮結構自身的動力性能，以及由結構振動引起的流固耦合效應。

近年來部分研究在一定程度上考慮了流固耦合效應。陳文禮等［3］采用基于RANS法的剪切應力輸運模型（簡稱SST模型），將圓柱看作一個單自由度的彈簧振子系統，并考慮了圓柱和流場之間的雙向耦合作用，觀測到了鎖定現象和“拍”的現象，且鎖定區域與試驗得到的結果較為接近，還將數值模擬得到的圓柱振動幅值與經驗公式和實驗數據進行了比較，結果吻合較好。方志平等［4］采用數值模擬方法研究了圓柱的兩自由度渦激振動，驗證了在低雷諾數和低折減阻尼條件下，流向自由度對橫向自由度的渦激振動基本沒有影響。詹昊等［5］將Newmark法的代碼嵌入用戶自定義函數（User Defined Function，UDF），將南京大勝關長江大橋吊桿系統簡化為質量—彈簧—阻尼系統，進行了渦激振動仿真計算，并將仿真分析結果與風洞試驗結果進行了比較，兩者的振幅、發振風速吻合良好，但仿真計算結果的鎖定區間要小于風洞試驗值。這些研究在一定程度上考慮了結構自身的動力效應，但主要以渦振為主，部分研究結果的精度有待進一步驗證。

本文利用微分方程的數值解法和動網格技術，基于松耦合方法，將Newmark算法通過UDF嵌入Fluent軟件中，在實現CFD（計算流體動力學）和CSD（計算結構動力學）耦合的基礎上，除研究方柱渦振之外，還進行了薄平板彎扭耦合顫振性能的研究，并將分析結果與經典解進行了檢驗分析，相關成果對橋梁渦振和顫振分析具有直接的指導意義。

1 CFD和CSD耦合分析方法

基于采用松耦合方法，將Newmark法嵌入UDF，在單個時間步內依次求解流體方程和結構運動方程，研究兩者的耦合作用。結構豎向和扭轉運動方程見式（1）和式（2）：

橋梁結構風致振動的卓越模態通常頻率較低，流體與結構之間的耦合可以認為是弱耦合。計算過程中，首先進行一步CFD計算，提取整個流場作用于結構的荷載數據；然后，將流場作用于結構的荷載分解為對結構扭心的力矩和兩個方向的力，分別代入式（1）和式（2）中，通過Newmark法計算得到結構的位移、速度以及加速度；最后，根據上一步得到的位移，采用動網格技術，將結構移動到相應的位置，并進行下一步的CFD計算。如此反復計算，直到結構收斂于某種狀態。其中，對于顫振臨界風速的計算，則需要首先給予結構一個初始位移，通過觀察結構的振幅是否發散來判斷結構是否發生顫振，具體流程見圖1。

圖1 顫振判斷依據Fig．1 Judgment basis for flutter

2 方柱渦振分析

當氣流流經鈍體橋梁結構斷面時，周期性交替脫落的渦旋會引起橋梁的振動，這種振動就稱作渦激振動。渦激振動在較低風速下就很容易出現，且帶有自激性質，但振動的結構反過來會對渦脫形成某種反饋作用，使振幅受到限制，因此渦激振動是一種帶有自激性質的風致限幅振動［6－7］。

2.1 方柱模型設置及靜態繞流計算

參考相關文獻［8－10］的建議，方柱截面尺寸取為0．033 m×0．033 m，計算區域大小見圖2，其中B1＝10b＝0．33 m，B2＝30b＝0．99 m，H1＝H2＝10b＝0．33 m。方柱貼體網格厚度2 mm，方柱周圍5b見方（即0．165 m×0．165 m）范圍內采用三角形網格，以便在單自由度振動方柱繞流計算時能夠將其設置為動網格；該區域右側15 d×5 d范圍設置為尾流區域，采用以四邊形為主的加密網格，以保證能準確模擬方柱脫落的旋渦；周圍區域的網格尺寸稍大，以四邊形為主。方柱網格劃分見圖3。

圖2 方柱計算區域Fig．2 Computational region of square cylinder

圖3 方柱網格劃分Fig．3 Mesh distribution of square cylinder

方柱繞流采用k－ωSST模型，設置左邊界為速度入口，右邊界為壓力出口，上下邊界為對稱邊界，方柱表面為無滑移固壁，來流風向從左到右，風攻角0°。當雷諾數Re＝22000時，進行靜態繞流計算，得到方柱阻力系數為2．11，斯托羅哈數為0．13，這與相關文獻吻合較好。

2.2 豎向單自由度振動方柱繞流計算

分別對風速為3 m／s、5 m／s、7 m／s以及9 m／s的工況進行靜態繞流計算，待風場穩定之后，將UDF以及動網格激活，進行流固耦合計算。不同風速下，方柱的位移時程曲線見圖4。由圖4可知，當風速較小時，方柱振幅很小；但當風速增大到7 m／s時，方柱振幅迅速增大；而風速增大到9 m／s時，方柱振幅又減小了。這是方柱發生渦振的典型表現，而且說明渦激共振的風速應該在5 m／s至9 m／s之間。

圖4 方柱位移時程曲線Fig．4 Time-history curve of displacement of square cylinder

為了進一步確定渦激共振的風速范圍，加算了三個工況，風速分別為6 m／s、6．5 m／s和8 m／s。不同風速下，方柱最大振幅的計算結果見圖5。由圖5可知，在風速為6．5m／s時，方柱的豎向振幅出現了最大值。

圖5 不同風速下方柱位移最大值Fig．5 Themaximum displacement value of square cylinder at differentwind velocity

渦脫頻率可以通過對升力系數時程曲線做傅里葉變換得到。不同風速下，方柱的渦脫頻率見圖6。從圖6可知，靜態繞流時，方柱的渦脫落頻率隨著風速的增大而增大，基本與風速成正比。當方柱發生豎向振動，其運動對流場造成了影響，使渦脫落頻率發生了微小的變化，而當風速在6～7m／s時，渦脫頻率基本固定在1．9 Hz左右，即發生了“鎖定”現象。

3 平板顫振分析

圖6 方柱繞流渦脫頻率對比Fig．6 Comparison of vortex shedding frequencies of square cylinder

顫振是一種由氣動自激力引起的發散振動，這主要是由于振動結構能夠不斷地從空氣中吸收能量，而結構所吸收的能量大于結構阻尼在振動中所耗散的能量。

3.1 薄平板模型設置

薄平板繞流采用二維計算，湍流模型采用k－ω SST模型，平板尺寸取0．8 m×0．004 m。參考相關文獻［8－9，11］的設置，將計算區域大小設置為20B× 11B，見圖7，其中B1＝7B＝5．6 m，B2＝12B＝9．6 m，H1＝H2＝5．5B＝4．4 m。平板周圍的貼體網格采用1 mm高度的三角形網格，平板運動區域也采用三角形網格，其余區域使用四邊形網格，網格劃分見圖8。

圖7 平板計算區域Fig．7 Computational region of flat plate

圖8 平板網格劃分Fig．8 Mesh distribution of flat plate

3.2 顫振導數識別

為了驗證平板模型網格尺度是否合理，以及網格運動后的計算結果是否準確可靠。首先，對該模型進行強迫振動計算，對計算結果進行顫振導數的識別，并與平板顫振導數理論解進行對比，從而判斷該模型是否正確。引入8個無量綱的顫振導數H*i，A*i，i＝1，2，3，4，自激力可以用下式表示：

式中：L為升力，M為扭轉力矩；無量綱因子 K＝2π／V*，V*為無量綱風速，V*＝U／（fB）；h·、h分別為平板豎向速度和豎向位移；α·、α分別為平板扭轉運動速度和扭轉位移。

假設該平板的豎向振幅h0＝0．03 m，扭轉振幅α0＝3°，振動頻率均為2 Hz，計算步長0．005 s，風速取3．2 m／s、6．4 m／s、9．6 m／s、12．8 m／s、16 m／s（對應的無量綱風速分別為2、4、6、8、10）。數值分析所得平板模型顫振導數結果與理論值的對比見圖9和圖10（限于篇幅，僅列出了H*1和A*1）。可見，該平板模型計算得到的數值解與平板理論解吻合較好，說明該模型計算參數設置合理，網格劃分合適，可用于后續的分析計算。

圖9 H1*值對比Fig．9 Comparison of H1*

圖10 A*1值對比Fig．10 Comparison of A*1

3.3 顫振臨界風速模擬結果

通過CFD-CSD耦合分析法，采用數值模擬計算平板的顫振臨界風速，以研究CFD-CSD耦合分析法計算的可行性。模擬過程中不再計算顫振導數，而是通過輸出平板運動的時程，觀察運動是否發散來判斷結構是否發生顫振。因此，該計算法的關鍵在于預估顫振臨界風速的大致范圍，然后逐步縮小風速范圍，逼近顫振臨界風速及頻率。平板的計算參數取值見表1。

表1 平板計算參數Tab.1 Calculation parameters of flat plate

通過多次計算發現，在20m／s風速下，平板豎向振動和扭轉振動均在較短時間內衰減；而當風速增大至20．5 m／s時，振幅雖然發生衰減，但是衰減速度明顯變慢；當風速為21m／s時，豎向和扭轉振動總體上趨于緩慢發散；風速繼續增大至22 m／s之后，振動發散明顯，說明已經發生顫振。20．5 m／s和21 m／s風速下的振動曲線見圖11和圖12。由此可以推斷，該平板的顫振臨界風速在20．5～21 m／s之間。

圖11 豎向位移及扭轉角時程（U＝20．5 m／s）Fig．11 Time-history curve of vertical displacement and torsion angle（U＝20．5 m／s）

根據所設置的薄平板參數，按照我國《公路橋梁抗風設計規范》［12］、ECCS規范以及相關理論公式計算得到平板顫振臨界風速見表2。由表2對比可知，本模擬結果與Matsumoto［13］給出的類似Selberg公式的計算結果20．7m／s較為吻合，與Selberg公式［14］的計算結果20．2 m／s和Scanlan半逆解法的計算結果19．7 m／s也十分接近。抗風規范中的計算結果明顯偏大，基于規范的結果可能偏于不安全。

圖12 豎向位移及扭轉角時程（U＝21 m／s）Fig．12 Time-history curve of vertical displacement and torsion angle（U＝21 m／s）

表2 平板顫振臨界風速對比Tab.2 Com parison of flutter critical w ind speed of flat plate

4 結 論

（1）本文以計算流體動力學軟件FLUENT為研究工具，運用CFD方法在每一時間步長內求解空氣動力學方程，把橋梁結構視為剛性模型，對其氣動力進行精確數值求解，然后把所得氣動力作為結構荷載，求解結構動力學方程，得到結構位移。其中，結構每個振型的運動都視為一個變常系數的二階微分方程，使用Newmark方法進行數值求解。結構的位置的更新則運用動網格技術實現。然后對運動后的結構進行CFD計算，得到新的氣動力再作用到結構單元，如此反復計算，直到收斂到某一結構形狀，對橋梁進行風致振動模擬求解。

（2）方柱豎向振動時，流場受到方柱振動的影響，旋渦脫落頻率有所變化。豎向單自由度振動時，低風速下渦脫頻率小于靜態繞流，而高風速下渦脫頻率大于靜態繞流。當旋渦脫落頻率與方柱固有振動頻率接近時，發生渦激共振，振幅明顯增大；隨著風速的繼續增大，在某個風速區間范圍內，渦脫頻率不再改變，即發生了鎖定現象。

（3）薄平板作為理想的顫振研究模型，在給予其初始位移的條件下，其振動在到達顫振臨界風速附近有明顯的衰減和發散現象，由此可以不斷逼近顫振臨界風速所在的風速區間。對于薄平板，數值模擬得到的顫振臨界風速與Scanlan理論公式和Selberg理論公式吻合較好。

（4）本文實現了CFD和CSD耦合的分析方法，并對方柱的渦振和平板的顫振問題進行了分析。計算結果表明本方法能較好的模擬風致振動中的流固耦合效應。對于高雷諾數下的流固耦合問題，如何建立誤差消除機制，并在保證計算精度的前提下提高計算效率，有待進一步研究。

［1］Schulz K W，Kallinderis Y．Unsteady flow structure interaction for incompressible flows using deformable hybrid grids［J］．Journal of Computational Physics，1998，143：569－597．

［2］Chio CK，Yu W J．engineering and industrial Finite element techniques for wind engineering［J］．Journal of Wind Aerodynamics，1999，81：83－95．

［3］陳文禮，李惠．基于RANS的圓柱風致渦激振動的CFD數值模擬［J］．西安建筑科技大學學報：自然科學版，2006，38（4）：509－513．

CHENWen-li，LIHui．CFD numerical simulation of vortexinduced vibration of a circular cylinder based on a RANS method［J］．Journal of Xi'an University of Archite-cture and Technology，2006，38（4）：509－513．

［4］方平治，顧明．圓柱兩自由度渦激振動的數值模擬研究［J］．同濟大學學報：自然科學版，2008，36（3）：295－298．

FANG Ping-zhi，GU Ming．Numerical simulation for vortexinduced vibration of circular cylinder with two degree of feedoms［J］．Journal of Tongji University：natural science，2008，36（3）：295－298．

［5］詹昊，方秦漢，李萬平．鋼桁拱橋吊桿渦激振動仿真分析［J］．中國鐵道科學，2009，30（2）：31－37．

ZHAN Hao，FANG Qin-han，LI Wan-ping．Numerical simulation of vortex-induced vibration of steel truss arch bridge hanger［J］．China Railway Science，2009，30（2）：31－37．

［6］陳政清．橋梁風工程［M］．北京：人民交通出版社，2005．

［7］項海帆等．現代橋梁抗風理論與實踐［M］．北京：人民交通出版社，2005．

［8］黃林．列車風與自然風聯合作用下的車－橋耦合振動分析［D］．成都：西南交通大學，2007．

［9］李永樂，汪斌，黃林，等．平板氣動力的CFD模擬及參數研究［J］．工程力學，2009，26（3）：207－211．

LIYong-le，WANG Bin，HUANG Lin，etal．CFD simulation and parameter study on aerodynamic force of flat plate［J］．Engineering Mechanics，2009，26（3）：207－211．［10］曾鍇，汪叢軍，黃本才，等．計算風工程中幾個關鍵影響因素的分析與建議［J］．空氣動力學學報，2007，25（4）：504－508．

ZENG Kai，WANG Cong-jun，HUANG Ben-cai，et al．Suggestion and analysis of several key factors in computational wind engineering［J］．Acta aerodynamica sinica，2007，25（4）：504－508．

［11］安偉勝．超大跨度分離三箱主梁橋梁抗風性能及氣動優化研究［D］．成都：西南交通大學，2011．

［12］中交公路規劃設計院．JTG／TD60－01－2004公路橋梁抗風設計規范［S］．行業標準－交通，2004．

［13］Matsumoto M，Matsumiya H，Fujiwara S，et al．New consideration on flutter properties basing on SBS-fundamental fluttermode，similar selberg's formula，torsional divergence instability，and new coupled flutter phenomena affected by structural coupling［C］．BBAA VI International Colloquium on：Bluff Bodies Aerodynamics＆Applications，Milano，Italy，2008．

［14］Selberg A．Oscillation and aerodynamic stability of suspension bridges［J］．ACTA Polytechnica Scandinavica，Civil Engineering and Construction Series，1961，13．

Aeroelastic simulation of vortex-induced vibration and flutter based on CFD／CSD coupling solution

LIYong-le1，ZHU Jia-qi1，2，TANGHao-jun1
（1．Department of Bridge Engineering，Southwest Jiaotong University，Sichuan，Chengdu 610031，China；2．Jiaxing Communications Engineering Quality and Safety Supervision Station，Jiaxing 314001，China）

Taking advantage of the software FLUENT and using the numerical solution of differential equation and the dynamic mesh model，a CFD／CSD coupling solution based on loose coupling was realized by embedding the Newmark method into FLUNT with the help of UDF function．A 2D-square cylindermodelwas established to investigate the change of Strouhal number and themaximum vertical vortex-excited amplitude of the square cylinder under differentwind speed．The lock-in phenomenon of vortex-excited resonance was observed in the process of simulation and itwas compared with the result of static square cylinder．A 2D flat platemodelwith vertical and torsional degrees of freedom was established to identify the flutter derivatives of the flat plate and to determine the flutter critical wind speed of flutter．The simulation result agreeswellwith the criticalwind speeds of flutter calculated by using the Scanlan's formula and Selberg's formula．

CFD-CSD coupling analysis；dynamic mesh；vortex-induced vibration；flutter；aeroelastic effect

O351．2

A

10．13465／j．cnki．jvs．2015．12．015

國家科技支撐計劃課題（2012BAG05B02）；四川省杰出青年學科帶頭人計劃（2009－15－406）

2014－02－08 修改稿收到日期：2014－06－03

李永樂 男，博士，教授，博士生導師，1972年生