亞臨界雷諾數(shù)下圓柱體尾流結(jié)構(gòu)的數(shù)值模擬

周 強(qiáng)，曹曙陽，周志勇

（同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)試驗(yàn)室，上海 200092）

過去幾十年來，圓柱體繞流問題因其簡單的幾何性質(zhì)和豐富的流場特性，一直是空氣動力學(xué)及風(fēng)工程領(lǐng)域中的研究熱點(diǎn)，尤其是對尾流結(jié)構(gòu)的研究.在亞臨界雷諾數(shù)范圍內(nèi)，橋梁中的吊桿、主纜以及拉索等具有圓柱形狀的構(gòu)件，在風(fēng)的作用下產(chǎn)生規(guī)則的旋渦脫落和復(fù)雜的尾流結(jié)構(gòu)，并發(fā)生風(fēng)致振動而導(dǎo)致結(jié)構(gòu)損傷或疲勞破壞.因此，研究在亞臨界雷諾數(shù)下的圓柱體繞流問題及其尾流結(jié)構(gòu)對工程結(jié)構(gòu)設(shè)計(jì)有著十分重要的意義.

當(dāng)雷諾數(shù)Re在40～2×105范圍內(nèi)，處于亞臨界區(qū)的圓柱體繞流有著豐富的特性.已有文獻(xiàn)表明當(dāng)Re≈300～3000時，剪切層開始從圓柱體表面分離，并在尾流區(qū)呈現(xiàn)非定常的特性.同時在Re＝3900的情況下，Ong等［1］采用熱線風(fēng)速儀（hot-wire anemometry，HWA）方法測量了回轉(zhuǎn)區(qū)外（3＜x／D＜10，D為圓柱體直徑）的速度及渦量分布，Lourenco等［2］采用粒子成像測速（particle image velocimetry，PIV）技術(shù)測試了非常近尾流區(qū)（x／D≤3）流場結(jié)構(gòu)，所以此后有關(guān)圓柱繞流的數(shù)值模擬及試驗(yàn)研究大多在這典型雷諾數(shù)（Re＝3900）下進(jìn)行.Beaudan等［3］采用O形貼面網(wǎng)格和高階迎風(fēng)格式，并運(yùn)用大渦模擬（LES）方法，對不可壓縮的N-S方程進(jìn)行求解，從而首先實(shí)現(xiàn)對此雷諾數(shù)下圓柱體繞流問題的數(shù)值模擬，其結(jié)果與Lourenco等［2］的試驗(yàn)結(jié)果很接近，但在回轉(zhuǎn)區(qū)域內(nèi)的平均流向速度分布形式表現(xiàn)以及回轉(zhuǎn)長度Lr上均存在一些差別.Mittal等［4］同樣運(yùn)用LES方法，同時采用二階有限中心差分格式及C形網(wǎng)格，結(jié)果顯示在下游區(qū)的雷諾應(yīng)力分布上與Lourenco等［2］的試驗(yàn)結(jié)果有一定的差別.Kravchenko等［5］采用基于B-splines的高階方法進(jìn)行數(shù)值研究，結(jié)果表明不足夠的格子數(shù)將會導(dǎo)致回轉(zhuǎn)長度減小，并且通過對比發(fā)現(xiàn)把展向長度從πD 增加到2πD 并不能影響結(jié)果，Breuer［6］得到同樣的結(jié)果.此后 Franke等［7］、Parnaudeau等［8］也采用LES方法，Ma等［9］采用直接模擬（DNS）方法對圓柱尾流結(jié)構(gòu)進(jìn)行了研究.此外，Parnaudeau等［8］采用更先進(jìn)的PIV設(shè)備對Re＝3900下的圓柱尾流結(jié)構(gòu)進(jìn)行了試驗(yàn)研究，結(jié)果與多個計(jì)算結(jié)果很近，而與Lourenco等［2］的試驗(yàn)結(jié)果有一些微小差別.國內(nèi)學(xué)者周志勇［10］、李壽英等［11］都研究過此類問題，但均未給出詳細(xì)的尾流結(jié)構(gòu).

本文采用二階有限中心差分離散格式和O形貼面網(wǎng)格，并細(xì)化圓柱體周圍5D范圍內(nèi)的網(wǎng)格，運(yùn)用基于Smagorinsky亞格子模型［12］的LES方法，對均勻來流條件下的圓柱體繞流問題進(jìn)行了三維數(shù)值模擬，驗(yàn)證了結(jié)果的準(zhǔn)確性，并重點(diǎn)分析了尾流區(qū)的流場結(jié)構(gòu)，從而得到一些有意義的結(jié)論.

1 控制方程與建模

1.1 控制方程

進(jìn)行濾波后，LES的連續(xù)性和N-S方程為

式中：ρ為流體密度；ui，uj為速度分量；p為壓力；μ為流體運(yùn)動黏度；式中帶有上劃線的量為濾波后的場變量被定義為亞格子尺度應(yīng)力，在Smagorinsky亞格子模型中

1.2 網(wǎng)格劃分和算法

計(jì)算區(qū)域：30D×20D×2D.模型上游來流區(qū)域?yàn)?0D，下游尾流區(qū)為20D，即x方向上的長度為30D；離上下邊界各為10D，即y方向上的寬度為20D；展向長度為Lz＝2D，如圖1所示.

網(wǎng)格劃分：在圓柱周圍5D范圍，采用O型貼面網(wǎng)格，格子數(shù)為200×150；展向方向（z方向）的格子數(shù)設(shè)為60；在遠(yuǎn)離圓柱體的區(qū)域使用稀疏的網(wǎng)格，而在圓柱體附近尤其是在尾流區(qū)，采用密集網(wǎng)格；為保證計(jì)算的準(zhǔn)確性，靠近圓柱表面的第一層網(wǎng)格的厚度取為d＝3×10-3D，然后網(wǎng)格以雙曲切線表示的間距分布由內(nèi)向外伸展；共1350000個單元.

圖1 計(jì)算區(qū)域示意圖及圓柱體周圍網(wǎng)格劃分情況Fig.1 Computational domain and grid example near the cylinder

邊界條件：圓柱體表面邊界條件采用無滑移固壁邊界；來流邊界條件采用均勻來流；出口邊界采用壓力出口，即壓力梯度（內(nèi)外壓差）為零；在展向方向上（即z方向）的兩個面，本文設(shè)置為周期性邊界條件；上下側(cè)面采用入流邊界條件，設(shè)置速度為來流速度，法線方向上的速度為零.

2 數(shù)值驗(yàn)證

在數(shù)值模擬過程中，時間步長直接影響到數(shù)值結(jié)果的準(zhǔn)確性.因此本文參考相關(guān)文獻(xiàn)，設(shè)置時間步長Δt＝0.0025D／U＝6.4×10-4s.此外，在對數(shù)據(jù)結(jié)果進(jìn)行統(tǒng)計(jì)處理時，所采用的數(shù)據(jù)數(shù)量將會影響到分析結(jié)果的準(zhǔn)確性，因此需要考慮數(shù)據(jù)的統(tǒng)計(jì)時長.這里引入無量綱時間參數(shù)T＊＝nU／D作為統(tǒng)計(jì)時長（averaging time），其中n為實(shí)際流動時長，U為來流速度.因?yàn)閚＝NcT，St＝fD／U，因此得到T＊＝Nc／St，其中T為旋渦脫落周期，f為旋渦脫落頻率，Nc為旋渦脫落周期數(shù).而對于Re＝3900的圓柱體繞流問題，其St數(shù)基本沒有變化，因此統(tǒng)計(jì)時長T＊只和旋渦脫落周期數(shù)Nc有關(guān).

如表1所示，Kravchenko等［5］采用7個旋渦脫落周期內(nèi)的數(shù)據(jù)結(jié)果，即無量綱統(tǒng)計(jì)時長為T＊＝35.00.Dong 等［13］采用Nc＝40～50，即T＊≈200.00～250.00；Ma等［9］則采用 Nc＝131，T＊＝624.89；Franke等［7］認(rèn)為在對數(shù)值結(jié)果處理時，需要采用40個旋渦脫落周期以上的數(shù)據(jù)，即T＊＞200.00；Parnaudeau等［8］分別統(tǒng)計(jì)了12個和52個旋渦脫落周期的數(shù)據(jù)，對比結(jié)果發(fā)現(xiàn)兩者的結(jié)果均與試驗(yàn)結(jié)果很吻合；Beaudan等［3］僅統(tǒng)計(jì)了6個旋渦脫落周期的數(shù)據(jù)，同樣取得較好的效果.因此本文在8個旋渦脫落周期（T＊＝38.00）內(nèi)，對每個時間步的計(jì)算數(shù)據(jù)進(jìn)行采樣處理，以得到準(zhǔn)確流場信息.

本文比較了St、平均阻力系數(shù)Cd、圓柱體后方平均壓力系數(shù)Cpb、剪切層平均分離角θ、平均回轉(zhuǎn)長度Lr等平均積分分量，如表1所示.圖2比較了本文得到的圓柱體表面平均壓力系數(shù)Cp、Kravchenko等［5］（Re＝3900）用 大 渦 模 擬 得 到 的 結(jié) 果 和Norberg［14］（Re＝4020）的試驗(yàn)結(jié)果.對比表明三者結(jié)果基本吻合，因此本文進(jìn)行的數(shù)值模擬是有效和可靠的.

表1 平均積分分量對比表Tab.1 Comparison of the mean integral quantities

圖2 圓柱體表面壓力系數(shù)分布Fig.2 Pressure coefficient distribution on the cylinder surface

3 結(jié)果及討論

3.1 平均積分分量

如表1所示，本文給出并比較了本文及其他文獻(xiàn)中基于時間積分的平均積分分量，包括Cd，Cpb，θ，St，Lr（在這里引入回轉(zhuǎn)長度，其定義為圓柱體中心到Umin對應(yīng)位置的距離）以及最小來流速度與來流速度的比值Umin／U.對比本文結(jié)果與其他試驗(yàn)和數(shù)值結(jié)果，可以發(fā)現(xiàn)盡管Cd，Cpb以及Lr存在微小偏差，但還是可以認(rèn)為其結(jié)果基本一致.同時還可以發(fā)現(xiàn)一個規(guī)律：當(dāng)Lr越小時，其Cd，Cpb以及Umin／U 表現(xiàn)得越大.本文認(rèn)為產(chǎn)生這個規(guī)律的原因是：Lr體現(xiàn)了在圓柱體后方形成的初始旋渦到圓柱體表面的距離，從而影響了圓柱體后方及其表面的負(fù)壓強(qiáng)度，進(jìn)而Lr的長短決定了Cd以及Cpb的大小，因此可以認(rèn)為平均Lr起到指標(biāo)性的作用.

3.2 平均流場特征

圖3 中心線上平均流向速度分布Fig.3 Mean streamwise velocity distribution on the central line

圖4 非常近尾流區(qū)不同位置上平均流向速度分布Fig.4 Mean streamwise velocity distribution at three locations in the very near wake

此外，在x／D＝2.02處平均流向速度分布上，本文結(jié)果和Parnaudeau等［8］結(jié)果吻合得很好，但較其他結(jié)果偏大.這是由于Lr的不同導(dǎo)致的.由圖3可以看出 Kravchenko等［5］以及 Lourenco等［2］的Lr值更小，從而更早跳出回轉(zhuǎn)區(qū)，因此在x／D＝2.02處u-的負(fù)值更小，甚至變正.由此可見，Lr是體現(xiàn)尾流流場結(jié)構(gòu)的重要指標(biāo).

如圖5所示，多個數(shù)值結(jié)果以及Parnaudeau等［8］的PIV試驗(yàn)結(jié)果表明，平均豎向速度分布應(yīng)該是關(guān)于y＝0軸呈現(xiàn)反對稱性，而Lourenco等［2］的結(jié)果并未表現(xiàn)出相類似的現(xiàn)象.此外由于Lourenco等［2］的試驗(yàn)采用的長徑比L／D＝21，并且前方來流長度較短.因此本文認(rèn)為導(dǎo)致數(shù)值結(jié)果與早期試驗(yàn)間偏差的主要因素是：圓柱體前方來流長度、長徑比以及早期試驗(yàn)設(shè)備的測量誤差.

圖5 非常近尾流區(qū)不同位置上平均豎向速度分布Fig.5 Mean cross-flow velocity distribution at three locations in the very near wake

圖6給出了近尾流區(qū)三個不同位置處的平均流向速度分布.由圖可見，Ong等［1］與 Kravchenko等［5］的結(jié)果均與本文一致.同時可以看出，較非常近尾流區(qū)，此區(qū)域的平均流場速度u-的分布剖面更平滑，表明此區(qū)域內(nèi)的旋渦強(qiáng)度不如上游區(qū)域.

圖6 近尾流區(qū)不同位置處平均流向速度分布Fig.6 Mean streamwise velocity distribution at three locations in the near wake

3.3 湍流特性

從圖7可以看出流向速度的脈動（u′u′／U2）在中心平面上的分布情況，其值從圓柱體表面的零值開始，然后逐漸增大，并達(dá)到兩個峰值，然后緩慢單調(diào)地減小.對比其他數(shù)值模擬和試驗(yàn)結(jié)果可以發(fā)現(xiàn)，本文的數(shù)值模擬結(jié)果和Parnaudeau等［8］的LES結(jié)果基本一致，并且均表現(xiàn)為后一個峰值比前一個峰值更大；但是和Parnaudeau等［8］的PIV結(jié)果以及Norberg［15］（Re＝5000）的試驗(yàn)結(jié)果有些偏差，其原因還是在于來流長度和長徑比的影響.此外，由于剪切層轉(zhuǎn)變出現(xiàn)位置不同，導(dǎo)致回轉(zhuǎn)長度的不同，從而使得脈動值的峰值出現(xiàn)位置以及峰值大小均表現(xiàn)不同.對于這一差別，后文中還將進(jìn)行進(jìn)一步解釋.

圖7 中心線上流向速度脈動值分布Fig.7 Streamwise velocity fluctuation value distribution of the central line

由于旋渦的交替出現(xiàn)使得尾流區(qū)的流向速度脈動值在兩側(cè)出現(xiàn)峰值，而非在中心處，如圖8所示.這表明尾流區(qū)出現(xiàn)卡門渦街現(xiàn)象.同時由圖可見，在x／D＝1.06處出現(xiàn)尖角峰值，其原因在于x／D＝1.06處于回轉(zhuǎn)區(qū)域內(nèi)，即剪切層在這區(qū)域內(nèi)發(fā)生轉(zhuǎn)變，因此初始旋渦在這區(qū)域內(nèi)交替出現(xiàn)，從而導(dǎo)致兩端出現(xiàn)尖角峰值，隨著流動向下游發(fā)展，峰值尖角消失，而變得更圓滑.此外由圖可以發(fā)現(xiàn)多個結(jié)果基本吻合，而在x／D＝1.06和x／D＝1.54兩處本文結(jié)果都稍小.如前面所述，出現(xiàn)這樣差別的原因在于：剪切層轉(zhuǎn)變出現(xiàn)位置不同，從而導(dǎo)致回轉(zhuǎn)長度不同，進(jìn)而使得脈動值的峰值出現(xiàn)位置也不同.如圖8所示，當(dāng)回轉(zhuǎn)長度越小，其峰值越早出現(xiàn)，從而導(dǎo)致在回轉(zhuǎn)區(qū)域內(nèi)其對應(yīng)的脈動值就越大，因此本文的脈動值在x／D＝1.06和x／D＝1.54較其他結(jié)果稍偏小.在尾流區(qū)其他位置上，本文結(jié)果與其他結(jié)果吻合得很好.

圖9給出了非常近尾流區(qū)三個不同位置處豎向速度脈動值的分布.由圖可見，在非常近尾流區(qū)豎向速度的脈動值均在中心線上達(dá)到峰值，并且可以發(fā)現(xiàn)數(shù)值結(jié)果比較一致，但與早期試驗(yàn)值有些偏差，其原因所前文所述.圖10給出了近尾流區(qū)三個不同位置處的雷諾應(yīng)力（u′v′／U2）分布，可見由于遠(yuǎn)離圓柱體，旋渦強(qiáng)度較小，因此在此區(qū)域內(nèi)的雷諾應(yīng)力值較小.

圖10 近尾流區(qū)不同位置處的雷諾應(yīng)力分布Fig.10 Reynolds shear stress distribution at three locations in the near wake

4 結(jié)論

（1）對St、剪切層平均分離角θ、平均阻力系數(shù)Cd以及圓柱體后方平均壓力系數(shù)Cpb等平均積分分量進(jìn)行比較分析，所得結(jié)果與試驗(yàn)值基本一致，從而驗(yàn)證了本文數(shù)值計(jì)算的準(zhǔn)確性.

（2）詳細(xì)給出了不同剖面處的平均速度、速度脈動值以及雷諾應(yīng)力的分布，從而得到了尾流區(qū)的平均流場及湍流流場結(jié)構(gòu)，為以后的試驗(yàn)和數(shù)值模擬提供詳細(xì)的參考.

（3）通過分析與早期試驗(yàn)間存在的微小偏差，發(fā)現(xiàn)導(dǎo)致這種偏差的原因可能在于圓柱體前方來流長度、長徑比的不同，并且需要進(jìn)一步的試驗(yàn)和數(shù)值驗(yàn)證.

（4）回轉(zhuǎn)長度是亞臨界狀態(tài)下圓柱體尾流結(jié)構(gòu)中最重要的特征，是檢驗(yàn)試驗(yàn)與數(shù)值計(jì)算的參考性指標(biāo).

［1］Ong L，Wallace J.The velocity field of the turbulent very near wake of a circular cylinder［J］.Experiments in Fluids，1996，20：441.

［2］Lourenco M，Shih C.Characteristics of the plane turbulent near wake of a circular cylinder.a particle image velocimetry study（data taken from ［3］），1993.

［3］Beaudan P，Moin P.Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds number，TF-62［R］.Stanford：Department of Mechanical Engineering of Stanford University，1994.

［4］Mittal R，Moin P.Suitability of upwind-biased finite difference schemes for large-eddy simulation of turbulent flow ［J］.American Institute of Aeronautics and Astronautics Journal，1997，35（8）：1415.

［5］Kravchenko A G，Moin P.Numerical studies of flow over a circular cylinder at Re＝3900［J］.Physics of Fluids，2000，12（2）：403.

［6］Breuer M.Large eddy simulation of the sub-critical flow past a circular cylinder：numerical and modeling aspects ［J］.International Journal Numerical Methods Fluids，1998，28：1281.

［7］Franke J，F(xiàn)rank W.Large eddy simulation of the flow past a circular cylinder at Re＝3900［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90：1191.

［8］Parnaudeau P，Carlier J，Heitz D，et al.Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900［J］.Physics of Fluids，2008，20（8）：1.

［9］Ma X，Karamanos S，Karniadakis E.Dynamics and lowdimensionality of a turbulent wake ［J］.Journal Fluid Mechanics，2000，410：29.

［10］周志勇.離散渦方法用于橋梁截面氣動彈性問題的數(shù)值計(jì)算［D］.上海：同濟(jì)大學(xué)土木工程學(xué)院，2001.ZHOU Zhiyong.Numerical calculation of aeroelastic problems in bridges by discrete vortex method［D］.Shanghai：College of Civil Engineering of Tongji University，2001.

［11］李壽英，顧明.斜、直圓柱繞流的CFD模擬［J］.空氣動力學(xué)學(xué)報(bào)，2005，23（2）：222.LI Shouying，GU Ming.Numerical simulation for flow around perpendicular and oblique circular cylinders ［J］. Acta Aerodynaimca Sinica，2005，23（2）：222.

［12］Lilly K.A proposed modification of the Germano subgrid-scale closure model［J］.Physics of Fluids，1992，4（4）：633.

［13］Dong S，Karniadakis E，Ekmekci A，et al.A combined direct numerical simulation particle image velocimetry study of the turbulent near wake［J］.Journal Fluid Mechanics，2006，569：185.

［14］Norberg C.LDV-measurements in the near wake of a circular cylinder［C］∥ Proceedings of the ASME Conference on Advances in the Understanding of Bluff Body Wakes and Vortex-Induced Vibration.Washington D C：ASME，1998：1-12.

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