Re＝3 900圓柱繞流的三維大渦模擬

戰慶亮，周志勇，葛耀君

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

Re＝3 900圓柱繞流的三維大渦模擬

戰慶亮，周志勇，葛耀君

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

為研究亞臨界雷諾數范圍內圓柱繞流流場特性及三維大渦模擬方法的適用性，基于C＋＋語言及有限體積法開發了三維非結構化網格的大渦模擬計算程序.采用新的高穩定性高精度二階離散格式，及Smagorinsky亞格子模型對Re＝3 900均勻來流條件下的圓柱繞流問題進行數值模擬，并統計獲得了平均流場參數及湍流流場的詳細結構特性.結果表明：采用本文的網格、計算步長和高穩定性二階離散精度大渦模擬方法計算所得的湍流場一階統計特性和二階統計特性與實驗值吻合很好.驗證了大渦模擬程序在模擬亞臨界雷諾數下圓柱繞流流場平均值及脈動值的合理性.

Smagorinsky亞格子模型；圓柱繞流；亞臨界雷諾數；非結構化網格；二階延遲修正

圓柱繞流由于其幾何形狀簡單，一直是鈍體空氣動力學、水動力學和風工程領域中的熱點研究問題，對工程實際應用有著重大意義.圓柱繞流的流場特性取決于雷諾數（Re＝ρUD／μ，其中U為來流速度，D為圓柱直徑，μ為運動粘性系數），繞流隨流動Re數的變化呈現極復雜的變化特征.而處于亞臨界情況下的圓柱繞流，有著豐富的流場特性和復雜尾流結構.

Re＝3 900情況下圓柱引起的湍流繞流問題是典型的亞臨界雷諾數問題，隨著計算機資源的發展，Re＝3 900圓柱繞流成為亞臨界計算的一個典型算例［1］.Re＝3 900的亞臨界圓柱繞流實驗研究數據比較豐富，可以作為驗證程序計算鈍體亞臨界雷諾數繞流準確性和算法適用性的標準，包括積分量（力，回流區長度等）和局部量（速度、渦量和雷諾應力等）.文獻［2］使用熱線風速儀（hot-wire anemometry，簡稱HWA）開創性的精確測量了回流區以外尾流的速度和渦矢量場，得到了湍流的統計信息和順流及橫向多處流場能譜.為了彌補HWA在流場回流區內部測量的不足，有學者采用激光粒子成像技術（簡稱PIV）和激光多普勒儀（LDV）來進行實驗測量，如文獻［1，3］采用PIV進行了近尾流區域內的流場結構更加細致和詳細的測量.直接數值模擬（DNS）結果也用來作為一種標準，主要用來研究剪切層的失穩現象［4］.數值模擬方面，文獻［5］首次對Re＝3 900的圓柱繞流進行了三維大渦模擬數值計算，采用的是O型貼體網格，近壁處的網格密度達到網格無關要求，平均速度和雷諾應力結果與實驗結果吻合較好，但是在回流區內部流向速度差別較大；文獻［6］基于B-splines高階格式同樣進行了計算，在回流區內部及外部都與實驗較好的吻合，同時文中還指出展向網格數量不足將使剪切層提前脫落，導致不準確的結果；文獻［7］基于C型曲面網格采用二階守恒形式的中心差分大渦模擬格式，進行了同樣雷諾數的計算，不同的是展向采用了譜方法離散，文中結果與實驗值吻合較好，但是下游區的雷諾應力由于離散格式的數值粘性衰減較快.另外，文獻［8－9］同樣采用大渦模擬進行了Re＝3 900圓柱繞流模擬.文獻［10］采用二維離散渦方法研究高雷諾數圓柱繞流問題，得到了較好的氣動力結果；國內其他作者大多采用商業計算流體軟件模擬圓柱繞流流場，如文獻［11］采用Fluent軟件進行了三維大渦模擬計算；文獻［12］對高Re數圓柱繞流應用Fluent軟件進行了二維的RANS分析，來討論其適用性問題.

本文以直徑D＝1.0 m的圓柱為研究對象，通過采用面向對象C＋＋語言編寫的非結構化網格有限體積計算程序結合大渦模擬方法，采用精度高穩定性好的二階離散格式，對Re＝3 900的圓柱繞流流場進行模擬.通過對計算得到的流場平均量與相關文獻研究結果對比，評價本文程序的有效性和算法穩定性，為高Re數及復雜斷面氣動力特性數值研究提供建議.

1 模擬方法

1.1 流場控制方程及有限體積離散

本文計算采用非結構化網格下的離散程序，非結構化網格具有復雜區域適應性好、局部加密靈活和便于自適應的優點，在流體數值模擬中得到越來越廣泛的應用.但在結構化網格中比較成熟的高階格式很難直接用于非結構化網格中，因此提高精度和簡化算法仍然是計算流體動力學研究的主要方向之一.對廣義輸運方程在控制單元體上進行積分，利用高斯定理將體積分轉為面積分后得到［13］

式中：φ為待求的標量，f表示單元體的面，n為圍成單元體面的個數.在上述方程中，非結構化網格離散計算中的關鍵是提高對流項和擴散項的計算精度.上游格式的數值穩定性好，但是數值粘性較大，只有一階精度.文獻［14］提出將QUICK格式用于非結構化網格的構思，但是沒有具體實現及研究.文獻［15］推導了三角形控制體非結構化網格QUICK計算方法，但其對網格適用性有限且程序實現復雜，且難于應用到三維計算中.本文提出了一種適用于三維非結構網格中二階精度格式，且數據結構簡單.

式中下標U表示上游單元體，rUf表示上游單元體中心指向面中心的矢量.例如圖1中非結構化網格局部所示，f為面中心，C和N為單元體中心.若f面的外法向朝向單元N，則公式變為

式中Ff＝Afρfφf（u·n）f為對流通量，n為面的外法向量.本格式的精度關鍵在于梯度的計算，可以采用高斯方法和最小二乘法，這兩種格式都用到了單元體其余面的信息，包含了遠場（即上游）的高階影響，故可提高精度.梯度的計算要保證其有界性，本文采用了限制子方法（Limiter）.公式中上標0代表延遲修正項，也就是采用了顯示處理，這樣能提高計算穩定性同時不損失精度.由于速度與壓力梯度的聯系體現在動量方程中，而連續性方程中不顯含壓力梯度，本文采用Rhie-Chow動量插值技術［16］克服由非交錯網格可能引起的非物理震蕩問題.在時間上采用三層全隱式格式，具有二階精度及良好的穩定性.

圖1 非結構化有限單元體離散示意

1.3 網格劃分與參數選擇

本文計算域為：流向長度30D，橫向長度20D，展向長度3D，其中D為圓柱直徑.在圓柱周圍5D范圍內采用O型網格，網格數為205×185.展向方向的網格數30.在遠離圓柱的流場區域使用較稀疏的網格，而近尾流區采用較密的網格，圓柱表面及局部網格見圖2.靠近圓柱表面的第一層網格的厚度取為d＝0.003D，底層網格Y＋值約為1，保證大渦模擬的計算準確性，平面內單層網格數為44 000.為保證數值模擬的準確性，根據克朗數（Crount number）要求選取時間步長為0.05 s.

圖2 圓柱表面及局部網格三維視圖

2 三維LES數值計算結果與分析

圖3顯示了瞬時的分離剪切層和尾流渦脫的發展情況：圓柱兩側發展出兩個剪切層，并在圓柱尾部交替斷裂脫落，形成有一定周期規律的渦脫，在漩渦脫落區內部流場較為復雜.圖中所示為低渦量等值面，其在離開圓柱一定距離之后發生失穩破碎，且具有明顯的三維結構特征和無序性，隨著流向方向逐漸減弱.在低渦量等值面中包裹著高渦量等值面，其結構尺度則更小，破碎效應更加明顯.

圖3 繞流流場瞬時渦量圖

2.1 平均積分特性及一階統計特性

平均積分特性和一階統計針對平均流場參數進行，主要分別對固壁受到流體作用力、流場順流方向速度、橫向速度進行統計均值分析，是研究繞流流場特性的重要內容.表1對比了平均積分量（平均阻力系數Cd，斯特勞哈爾數Sr，平均回流長度r，最小平均速度Umin等量），與文獻結果的差異很小.

表1 積分特性對比

圖4為漩渦穩定脫落后的阻力與升力系數時程曲線.阻力基本穩定在一定數值上，本文得到的平均阻力系數為1.03，這與實驗和數值模擬結果都很接近.觀察圖4中升力系數時程，可以發現其均值為零，對比低雷諾數下圓柱繞流計算結果可發現，在亞臨界雷諾數區間的升力系數幅值已不再平穩，隨時間在一定范圍內變化，然而仍存在明顯的周期特性.

圖5給出了中心平面（y＝0）上平均流向速度＜u＞分布情況.圓柱壁面處＜u＞為零，在流向方向由于存在回流區，因此＜u＞的值逐漸減小到Umin＝0.33.回流區中心的位置距圓柱中心位置也稱為回流長度（recirculation length），用r表示.下游區＜u＞值逐漸變大，單調的逼近遠處流場速度.圖5顯示本文計算結果與參考文獻中的實驗結果趨勢及數值上吻合較好，尤其是與文獻［1］近年的PIV實驗吻合最好，該實驗結果采用更為先進、精細的實驗方法所得，因此認為其參數較以前文獻中數據更為可信.

圖4 阻力系數與升力系數時程曲線

圖5 中心平面處平均流向速度分布

圖6為平均流場順流速度云圖，可以更加細致直觀地分析尾流區流場結構，圖7為平均橫向速度云圖.由圖6可見，平均流向速度關于y＝0對稱，這是由于本文研究的圓柱斷面形狀是對稱的所致.在圓柱尾部有一個明顯的＜u＞小于0的回流區，表明此處流場平均速度與順流方向相反，這也是“回流區”名字的由來.文獻［1，3］實驗研究中，給出了尾流回流區域內及回流區域外不同位置處流向速度剖面圖，因此本文也選取相應位置與實驗值進行比對分析.在回流區內部，本文選取3個典型剖面進行分析，分別是x／D＝1.06，x／D＝1.54，x／D＝2.02；在回流區外部，本文也選取了3處典型位置，分別是x／D＝4.0，x／D＝7.0，x／D＝10.0來進行對比驗證.

圖8為尾流區6個位置處平均流向速度計算值與實驗值的對比，其中黑色直線為本文結果，方形為文獻［3］結果，三角為文獻［1］中PIV實驗結果，圓圈為文獻［2］中實驗結果（下文亦同）.本文計算結果與實驗結果相比得到了很好的模擬效果.在回流區內部平均流向速度減速明顯，很顯然是由于鈍體對流場的阻礙作用引起.同時還可發現平均流向速度剖面隨著流場的發展由“U”字形逐漸變為“V”字形，其減速效應越來越區域平緩.在回流區內部，本文數值模擬結果與文獻［1］的PIV實驗結果更為接近.在回流區外部，平均順流速度剖面更為平滑，且數值差異不大，表明在遠離圓柱5D范圍以外流場受圓柱影響已經較小.圖9為相應尾流不同區域處平均橫向速度剖面對比圖，在回流區內部，圖形形狀為反對稱，本文結果與文獻［1］實驗結果吻合較好，但與文獻［3］試驗值有一定誤差.

圖6 平均流向速度云圖

圖7 平均橫向速度云圖

圖8 尾流區不同位置平均流向速度分布

2.2 二階統計特性

二階統計特性主要對流場物理量的脈動值進行統計分析，代表了湍流流場的脈動特性，即雷諾應力的統計情況.對于湍流流動，速度的脈動情況、雷諾應力分布是湍流流動的重要特征.圖10所示為中心平面處順流向流場雷諾正應力的分布圖.可以看出流向速度脈動沿順流向的分布情況，在圓柱表面處其值為零，隨著遠離圓柱其值逐漸增大并達到兩個峰，然后遞減.本文與實驗值及參考文件中數值模擬結果都顯示遠離圓柱的峰值較大，峰值處流場湍流雷諾應力較大.

圖9 尾流區不同位置平均橫向速度分布

圖10 中心平面雷諾應力分布

圖11為順流向速度脈動值平均云圖.順流向速度脈動值出現兩處最大值，通過圖10得該距離為2.044 7D.流向速度脈動值與漩渦脫落關系密切，最大值位置可以代表平均漩渦脫落的位置，同時兩個峰值對應著圓柱兩側分別產生的漩渦.從流向速度脈動云圖上可以得出：漩渦被拉長之后脫落，這種脫落形態與層流下的卡門漩渦脫落形態是不一樣的.對比圖13流場不同位置順流速度脈動剖面圖，可以發現本文結構很好吻合了實驗結果.在x／D＝1.06處順流向速度脈動均值出現兩個尖角，這與實驗中觀察到的現象一致.本文認為這是由于圓柱兩側漩渦被流場拉伸至下游，而此處漩渦并未斷裂，寬度也沒有很大變化，反映在流向速度脈動剖面圖上即為兩個尖角形狀.沿流線方向兩側漩渦半徑變寬，同時發生斷裂脫落現象，對應著x／D＝1.54，x／D＝2.02處的剖面曲線變得光滑，其值也變大.在尾流回流區外部處，可以看出順流向雷諾應力平緩變小，此處的流場受圓柱繞流影響已經很小.

圖11 流向速度脈動云圖

圖12為橫向速度脈動云圖，可以看出橫向速度脈動最大峰值只有一個，這與順流向速度脈動是不同的.圖中表明橫向脈動在中心線y＝0軸線上達到最大值，說明中心線處橫向速度脈動量，即橫向湍流應力最大，此處流場橫向能量交換最為明顯.圖14中尾流區不同位置橫向速度脈動分布圖中觀察到，在x／D＝1.06處橫向速度脈動量非常小，說明在非常近圓柱處幾乎沒有橫向脈動.其余各位置剖面圖也與實驗結果有很好的吻合度.

圖12 橫向速度脈動云圖

圖13 尾流區不同位置流向速度脈動分布

圖14 尾流區不同位置橫向速度脈動分布

3 結 論

1）針對一階差分數值耗散過大和二階中心差分數值穩定性差的問題，本文采用新的高穩定性二階插值格式用于NS方程離散，并將其與延遲修正技術相結合，實現了較高精度、高穩定性的大渦模擬計算.

2）使用大渦模擬方法對亞臨界雷諾數圓柱繞流流場進行計算模擬，對平均阻力系數、斯特勞哈爾數、回流區長度等平均積分量進行比較分析，驗證其與實驗值基本一致.

3）對流場尾流中心線及回流區內部和尾流區不同位置流場一階統計特性和二階統計特性進行統計，通過與相對應實驗值進行對比，驗證了本文計算的回流區內部和回流區外部尾流流場的平均值及脈動值的合理性.

4）通過對比驗證，說明三維Smagorinsky亞格子模型以及本文離散格式可以準確模擬亞臨界雷諾數圓柱繞流問題，并考證了本文采用有限體積法所開發程序的有效性.

［1］PARNAUDEAU P，CARLIER J， HEITZ D， et al.Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3 900［J］.Physics of Fluids（1994-present），2008，20（8）：85－101.

［2］ONG L，WALLACE J.The velocity field of the turbulent very near wake of a circular cylinder［J］.Experiments in Fluids，1996，20（6）：441－453.

［3］LOURENCO L M，SHIH C.Characteristics of the plane turbulent near wake of a circular cylinder：A particle image velocimetry study［Z］.Private Communication（data taken from［8］），1993.

［4］DONG S，KARNIADAKIS G，EKMEKCI A，et al.A combined direct numerical simulation-particle image velocimetry study of the turbulentnearwake［J］.Journal of Fluid Mechanics，2006，569（1）：185－207.

［5］BEAUDAN P，MOIN P.Numerical experiments on the flow past a circular cylinder at sub－critical Reynolds number［R］.Stanford：Stanford University，1994.

［6］KRAVCHENKO A G，MOIN P.Numerical studies of flow over a circular cylinder at ReD＝3 900［J］.Physics of Fluids，2000，12（2）：403－417.

［7］MITTAL R，MOIN P.Suitability of upwind-biased finite difference schemes for large-eddy simulation of turbulent flows［J］.AIAA Journal，1997，35（8）：1415－1417.

［8］FRANKE J，FRANKW.Large eddy simulation of the flow past a circular cylinder at＜i＞Re＜／i＞＜sub＞D＜／sub＞＝3 900［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90（10）：1191－1206.

［9］LYSENKO D A，ERTESV G IS，RIAN K E.Large-eddy simulation of the flow over a circular cylinder at Reynolds number 3 900 using the Open FOAM toolbox［J］.Flow，Turbulence and Combustion，2012，89（4）：491－518.

［10］周志勇.離散渦方法用于橋梁截面氣動彈性問題的數值計算［D］.上海：同濟大學，2001.

［11］周強，曹曙陽，周志勇.亞臨界雷諾數下圓柱體尾流結構的數值模擬 ［J］.同濟大學學報（自然科學版），2013，41（1）：33－38.

［12］祝志文.高Re數圓柱繞流二維數值模擬的適用性分析［J］.振動與沖擊，2013，32（7）：98－101.

［13］FERZIGER J H，PERIC′M.Computational methods for fluid dynamics［M］.Berlin：Springer，2002.

［14］DAVIDSON L.A pressure correction method for unstructured mesheswith arbitrary control volumes［J］.International Journal for Numerical Methods in Fluids，1996，22（4）：265－281.

［15］姜華，席光.對流項二次迎風插值格式在非結構化網格中的應用［J］.西安交通大學學報，2006，11（1）：1246－1249.

［16］ISSA R I.Solution of the implicitly discretised fluid flow equations by operator-splitting［J］.Journalof Computational Physics，1986，62（1）：40－65.

［17］NORBERG C.An experimental investigation of the flow around a circular cylinder：influence of aspect ratio［J］.Journal of Fluid Mechanics，1994，258：287－316.

（編輯趙麗瑩）

3-Dimensional large eddy simulation of circular cylinder at Re＝3 900

ZHAN Qingliang，ZHOU Zhiyong，GE Yaojun

（State Key Laboratory for Disaster Reduction in Civil Engineering（Tongji University），20092 Shanghai，China）

This paper developed finite volume method program with Smagorinsky sub-grid scale model using 3-dimensional unstructuredmesh to accurately determine the flow field characteristics behind a cylinder in sub critical region and the applicability of 3-dimensional large eddy simulationmethod of this kind of calculation based on C＋＋programming language.Numerical simulation was performed for the flow over a circular cylinder at a classical subcritical Reynolds number（Re＝3 900）using a new second order defer correction scheme for unstructured mesh.The results show thatmean and fluctuating flow characteristicswere compared wellwith the existing experimental results，2-order discretization scheme large eddymethod can be used in the simulation of sub-region Reynolds number outer flow simulation.The results also proved the reliability of high accuracy and stability discretization scheme.

Smagorinsky sub grid scale model；flow over a circular cylinder；sub-critical region；unstructured mesh；second order defer correction

U441

A

0367－6234（2015）12－0075－05

10.11918／j.issn.0367-6234.2015.12.013

2015－03－30.

國家自然科學基金重大研究計劃集成項目（91215302）；國家重點基礎研究發展計劃（973計劃）（2013CB036300）.

戰慶亮（1987—），男，博士研究生；周志勇（1971—），男，研究員，博士生導師；葛耀君（1958—），男，教授，博士生導師.

周志勇，z.zhou＠tongji.edu.cn.