基于轉(zhuǎn)捩SST模型凸起圓柱繞流數(shù)值研究

呂代龍, 陳少松, 周 航, 徐一航

(南京理工大學能源與動力工程學院, 江蘇南京 210018)

引 言

圓柱繞流是鈍體繞流的重要研究案例, 是流體力學的經(jīng)典問題.圓柱繞流問題是研究流動分離、渦流和渦脫落的重要基礎(chǔ), 在航空航天、船舶海洋以及兵器科學等領(lǐng)域有著重要的工程應用價值, 例如彈箭表面的電纜罩會引起流動分離, 使彈箭產(chǎn)生側(cè)向力, 該問題可以簡化為帶有局部凸起的圓柱繞流問題.

圓柱繞流問題在早期的研究過程中局限性大, 研究的成果也相對簡單, 大多以模型實驗研究為主.國內(nèi)外許多學者圍繞圓柱繞流問題做了大量研究工作, 總結(jié)了不同Reynolds數(shù)下圓柱繞流的氣動特性及變化規(guī)律.Aguirre-López等[1]采用直接數(shù)值模擬(direct numerical simulation, DNS)方法對Re=1.46×105帶凸起的圓柱繞流問題進行了數(shù)值模擬, 得到了對圓柱氣動特性影響最大的凸起位置.但是對于高Reynolds數(shù)下帶凸起的圓柱繞流近壁面流動特征、渦的產(chǎn)生與脫離及流動機理等方面的研究相對較少, 對流動過程中邊界層轉(zhuǎn)捩的影響也很少關(guān)注.

邊界層轉(zhuǎn)捩現(xiàn)象是流體力學研究中一個難點問題, 是流體從層流到湍流狀態(tài)的過渡階段.盡管關(guān)于邊界層轉(zhuǎn)捩的基礎(chǔ)理論和發(fā)生機理仍在研究發(fā)展中, 但其在工程實踐應用中的重要性日益突出.而對于帶有電纜罩的彈箭, 當電纜罩處于彈箭的非對稱位置(如彈箭的側(cè)面), 彈箭氣動特性的不對稱便會使彈箭產(chǎn)生側(cè)向力, 這種現(xiàn)象在大攻角狀態(tài)下更為明顯.因此準確預測轉(zhuǎn)捩點的位置、以及壁面的流動分離可以改善彈箭氣動特性.

1 數(shù)值計算方法

1.1 轉(zhuǎn)捩動量厚度Reynolds數(shù)輸運方程

轉(zhuǎn)捩動量厚度Reynolds數(shù)是轉(zhuǎn)捩起始點的決定因素, 其輸運方程為[11]

輸運方程的生成項為

式中,cθt為常數(shù);Fθt為開關(guān)函數(shù), 該函數(shù)的值從邊界層內(nèi)部到外部由1逐漸變?yōu)?;Reθt為當?shù)剞D(zhuǎn)捩Reynolds數(shù), 該變量由來流湍流度Tu和壓力梯度參數(shù)λθ擬合的經(jīng)驗公式得到.

1.2 間歇因子輸運方程

邊界層的轉(zhuǎn)捩、湍流模型中的渦黏性系數(shù)等均由間歇因子控制, 其輸運方程為[6]

式中,Pγ為生成項,Dγ為耗散項,γ為間歇因子.

1.3 SST k-ω湍流模型

SSTk-ω湍流模型的輸運方程如下[13]

2 數(shù)值模擬條件及網(wǎng)格

2.1 數(shù)值模擬條件

本文研究的是來流速度U∞沿x軸正方向流向直徑為D的二維帶凸起的圓柱繞流問題, 如圖1所示.圖中θ為圓柱表面的方位角, 當θ=45°時, 圓柱表面上帶凸起.分別取來流Re=1.4×105(亞臨界區(qū))，2×105，5×105(臨界區(qū))，1×106(超臨界區(qū)), 根據(jù)實驗條件選擇來流的湍流強度為0.8%.

圖1 流動示意圖

時均阻力系數(shù)CD, Strouhal數(shù)St, 時均摩擦力系數(shù)Cf和時均壓力系數(shù)Cp分別可由下列各式確定

式中,FD為圓柱受到的時均阻力,f為渦脫落的頻率,τ0為壁面剪切應力.

2.2 計算網(wǎng)格及邊界條件

如圖2所示, 采用尺寸為55D(來流方向)×40D(橫流方向)的矩形區(qū)域作為計算域, 計算域上下邊界和計算域左邊界與圓柱中心的距離均為20D, 計算域右邊界與圓柱中心的距離為35D.該尺寸保證了計算域足夠大, 可有效避免邊界對流場產(chǎn)生影響.物面法向第1層網(wǎng)格高度由y+～1確定.圖3為圓柱表面附近網(wǎng)格的示意圖.

圖2 整個計算域網(wǎng)格示意圖

圖3 圓柱表面附近網(wǎng)格示意圖

計算邊界條件定義如下: 計算域的左邊界為速度入口邊界條件, 在此邊界上指定均勻流速為u=U∞且v=0, 其中u和v分別是x和y方向上的速度分量; 圓柱表面采用壁面無滑移邊界條件, 即u=v=0; 計算域的上下邊界采用對稱邊界條件; 在計算域的右邊界選用流動出口邊界條件.

3 數(shù)值模擬結(jié)果

3.1 亞臨界區(qū)數(shù)值模擬結(jié)果分析

表1給出了Re=1.4×105時, 通過數(shù)值模擬計算得到的二維帶凸起圓柱的時均阻力系數(shù)CD和Str-ouhal數(shù)St的值, 也給出了Schewe[11]和Cantwell[14]通過實驗測得的結(jié)果及苑明順[10]數(shù)值模擬計算結(jié)果.對比表1中CD和St值, 可以看到所得結(jié)果與實驗測得的結(jié)果基本符合, 說明在亞臨界區(qū)采用轉(zhuǎn)捩 SST模型能夠很好地模擬二維圓柱繞流問題.

表1 亞臨界區(qū)二維圓柱時CD和St數(shù)值模擬計算結(jié)果

圖4給出了Re=1.4×105時通過轉(zhuǎn)捩 SST與LES模擬得到的二維圓柱上下表面的時均壓力系數(shù)Cp的分布曲線, 與Aguirre-López 等[1]通過DNS 模擬得到的結(jié)果進行對比可以看出: 與LES相比, 圓柱表面Cp的數(shù)值模擬結(jié)果與DNS得到的結(jié)果吻合較好, 在θ=40°～320°范圍內(nèi), 上表面的Cp值小于下表面, 同樣從圖5壓力云圖可以看出, 圓柱上表面的壓力小于下表面的壓力, 這是由于在上表面凸起位置的前后產(chǎn)生了3個反向旋轉(zhuǎn)的漩渦, 形成了背風面低壓區(qū), 如圖6所示.同時由于渦的存在, 使渦上方流體的流通通道變窄, 造成流動速度增大、壓強降低, 而在邊界層內(nèi), 沿物體表面的法線方向壓強保持不變, 即等于外邊界處自由流的壓強致使上表面的壓強相對于下表面較低.

圖4 亞臨界區(qū)圓柱表面時均壓力系數(shù)Cp分布曲線

圖5 亞臨界圓柱表面壓力云圖

圖6 漩渦流線圖

圖7給出了Re=1.4×105時通過數(shù)值模擬得到的二維圓柱上下表面的時均摩擦力系數(shù)Cf的分布曲線, 與Aguirre-López等[1]通過DDES方法得到的結(jié)果進行對比.由圖4, 7可知, 轉(zhuǎn)捩 SST模型能夠相對準確地模擬圓柱繞流背風面壓力和摩擦力的變化情況, 與DNS得到的結(jié)果相差不大.

圖7 亞臨界圓柱表面時均摩擦力系數(shù)Cf分布曲線

從圖7可以看出, 采用轉(zhuǎn)捩 SST模型模擬得到的下表面計算結(jié)果在θ=88°時Cf由正變?yōu)樨? 這表明此位置發(fā)生了流動分離現(xiàn)象; 而計算結(jié)果在θ=90°～130°范圍內(nèi)與DDES方法得到的不同變化的原因, 是由于流動處于亞臨界區(qū)向臨界區(qū)的過渡階段, 流動分離之后有一個分離再附的趨勢, 即有一個形成分離泡的趨勢.上表面在θ=38°處Cf的正負符號發(fā)生了變化, 這表明在凸起的影響下, 流動在θ=38°處發(fā)生了流動分離; 在θ=38°～60°范圍內(nèi)Cf處在負區(qū)間, 這是由于凸起的影響, 流動在凸起前后流動產(chǎn)生了漩渦; 在θ=60°～106°范圍內(nèi)Cf隨θ先增大后減小, 在θ=73°時Cf達到最大值, 說明流動分離之后在此范圍內(nèi)再附; 而在θ=106°處Cf的正負符號又發(fā)生了改變, 流動在這個位置又產(chǎn)生流動分離現(xiàn)象; 此處之后, 與下表面相比, 上表面的Cf圍繞Cf=0上下震蕩, 是因為表面凸起的影響, 流動分離后再附產(chǎn)生了分離泡現(xiàn)象.

對于圓柱繞流問題, 尾跡區(qū)的近壁區(qū)流動特征主要是分離與局部二次分離產(chǎn)生不同尺度漩渦的合并、成對、分叉等強烈的相互作用, 漩渦脫落的過程實際上是多個渦間及與剪切層相互作用的結(jié)果.在高Reynolds數(shù)下, 漩渦只要強度足夠大或者距離壁面足夠近, 二維或三維漩渦皆可在近壁區(qū)誘導產(chǎn)生局部分離, 形成新的漩渦.

圖8給出了凸起圓柱在背風面尾跡區(qū)流動的一個周期內(nèi)的變化過程.t1時刻尾跡區(qū)圓柱上表面附近產(chǎn)生的漩渦距離壁面非常近, 由于其誘導作用, 使尾跡區(qū)圓柱下表面在t2時刻發(fā)生了流動分離, 產(chǎn)生了新的漩渦, 且在尾跡區(qū)上表面附近的漩渦隨著時間的發(fā)展逐漸變大并沿流向向后發(fā)展, 由于其與尾跡區(qū)下表面附近的剪切層相互作用, 在t3時刻尾跡區(qū)圓柱下表面附近產(chǎn)生了3個小漩渦, 并且這3個小漩渦與之前上表面附近產(chǎn)生的小漩渦由于旋向以及來流的影響開始合并, 在t6時刻完成合并, 成為一個正向旋轉(zhuǎn)的漩渦.同時, 由于渦間以及來流的影響, 在t6時刻形成的正向漩渦和之前由誘導作用產(chǎn)生局部分離形成的漩渦開始合并, 在t7時刻完成合并, 發(fā)展成為一個大的分離渦, 然后隨著時間的推移發(fā)生漩渦脫落.

圖8 臨界區(qū)圓柱尾跡區(qū)一個周期的變化過程

綜上所述, 轉(zhuǎn)捩 SST模型對流動分離和壓力梯度等因素敏感, 能很好地模擬二維凸起圓柱表面的壓力、摩擦力變化以及尾跡近壁區(qū)的流動分離等流動特性.

3.2 臨界區(qū)數(shù)值模擬結(jié)果分析

圖9為流動處于Re=2×105和Re=5×105時數(shù)值模擬計算所得的圓柱上下表面時均壓力系數(shù)Cp分布曲線.可以看出, 兩個Reynolds數(shù)下Cp的分布曲線趨勢大致相同,Cp的值都是逐漸減小, 這也符合其變化規(guī)律; 對比兩個Reynolds數(shù)下圓柱表面的分布曲線, 由于凸起的影響, 在凸起附近的壓力系數(shù)發(fā)生突變, 之后的壓力系數(shù)逐漸與光滑表面的相近.與亞臨界區(qū)流動相同, 臨界區(qū)的圓柱上表面時均壓力系數(shù)Cp與下表面相比, 也要稍小, 其原因與亞臨界區(qū)流動相同.

圖9 臨界區(qū)圓柱表面時均壓力系數(shù)Cp分布曲線

圖10為流動處于Re=2×105和Re=5×105時數(shù)值模擬計算得到的圓柱表面時均摩擦力系數(shù)Cf分布曲線, 可以看出, 數(shù)值模擬得出了圓柱表面時均摩擦力系數(shù)Cf突然增大又回落減小的過程, 也就是圓柱邊界層的轉(zhuǎn)捩過程.

圖10 臨界區(qū)圓柱表面摩擦力系數(shù)Cf分布曲線

流動處于Re=5×105時下表面時均摩擦力系數(shù)Cf在圓柱下表面θ=130°～150°范圍內(nèi)由正變負再變正, 表明邊界層內(nèi)發(fā)生了流動分離和再附的過程, 產(chǎn)生了分離泡.Cf在圓柱上表面θ=30°～60°由正變負又變正, 這就表明由于凸起的影響, 邊界層在此處發(fā)生了流動分離, 產(chǎn)生了漩渦現(xiàn)象; 之后Cf的值突然變大, 又回落減小, 這說明流動重新附于圓柱表面上; 在圓柱上表面θ=105°～146°時Cf由正變負又變正, 這就表明在邊界層流動發(fā)生了分離和再附的過程, 即產(chǎn)生了分離泡.同時, 與流動處于亞臨界區(qū)時上表面產(chǎn)生流動分離現(xiàn)象的位置大致相同.

圖11為流動處于臨界區(qū)時圓柱上表面附近流動區(qū)域的速度矢量圖, 從左至右依次為a,b,c,d,e點.從a可以看出, 在此區(qū)域內(nèi)流向圓柱上表面外法線方向的速度梯度逐漸變大, 使流向圓柱上表面剪切應力逐漸增大, 導致時均摩擦力系數(shù)Cf逐漸變大.

圖11 臨界區(qū)二維凸起圓柱速度矢量圖

從b和c可以看出, 在圓柱上表面θ=30°～60°這個區(qū)域內(nèi)沿流向圓柱上表面氣流的速度梯度在逐漸減小, 速度矢量方向與來流方向先相同后相反再相同, 這就使得沿流向上表面時均摩擦力系數(shù)Cf逐漸減小, 最終由負又變正.

從d和e可以看出, 在此區(qū)域內(nèi)圓柱上表面氣流的速度矢量方向與來流方向先相同后相反再相同, 說明在流動發(fā)生了分離后再附, 形成了一個分離泡.

通過采用轉(zhuǎn)捩 SST模型對臨界區(qū)凸起圓柱繞流問題的計算分析, 可以看出, 氣流在凸起前后形成了3個反向旋轉(zhuǎn)的漩渦, 形成了背風面低壓區(qū), 之后隨著上表面θ的增加, 氣流先是發(fā)生了流動分離現(xiàn)象, 產(chǎn)生了分離泡, 然后發(fā)生了轉(zhuǎn)捩, 并發(fā)生了流動分離.

3.3 不同Reynolds數(shù)下模擬結(jié)果對比

圖12給出了在不同Reynolds數(shù)下得到的圓柱表面時均摩擦力系數(shù)Cf的分布曲線, 圖13為不同Reynolds數(shù)下圓柱邊界層的間歇因子γ的分布曲線, 與采用DNS模擬得到的γ的分布曲線吻合較好.

圖12 不同Reynolds數(shù)下圓柱表面時均摩擦力系數(shù)Cf分布曲線

(a)Upper surface of cylinder

從圖12可以看出, 在凸起前后, 圓柱上表面時均摩擦力Cf的符號為負, 這說明流體在此發(fā)生了流動分離, 產(chǎn)生了分離渦, 結(jié)合流場結(jié)構(gòu)可以得出, 在凸起處形成了3個逆向旋轉(zhuǎn)的小渦.同時隨著Reynolds數(shù)的增加, 氣流發(fā)生分離再附的位置逐漸向前.

當流動處于Re=2×105時, 在圓柱上表面θ=50°處的時均摩擦力系數(shù)Cf先增大后減小, 且由圖13(a)可得, 在此處邊界層的間歇因子γ從0突變到1, 表明流動在θ=50°處發(fā)生了轉(zhuǎn)捩.而在圓柱下表面, 轉(zhuǎn)捩的位置在下表面θ=90°左右處, 明顯比上表面的轉(zhuǎn)捩位置滯后, 這表明凸起可以使流動轉(zhuǎn)捩過程提前.從圖13(b)可以看出, 當Re=5×105和Re=1×106時, 圓柱下表面邊界層內(nèi)發(fā)生了轉(zhuǎn)捩過程, 其轉(zhuǎn)捩點分別位于下表面θ=89°和θ=85°處, 說明當Reynolds數(shù)處于臨界區(qū)和超臨界區(qū)時, 隨著Reynolds數(shù)的增加, 圓柱邊界層轉(zhuǎn)捩點的位置逐漸向來流方向移動.

綜上所述, 對于不同Reynolds數(shù)下的圓柱表面, 氣流會在圓柱表面凸起處產(chǎn)生3個反向旋轉(zhuǎn)的小漩渦, 發(fā)生流動分離現(xiàn)象.同時圓柱上表面的凸起可以使氣流的轉(zhuǎn)捩提前發(fā)生, 并且隨著Reynolds數(shù)的增加, 圓柱邊界層轉(zhuǎn)捩點的位置逐漸向后移動.

4 結(jié)論

本文采用轉(zhuǎn)捩 SST模型對不同Reynolds數(shù)下的二維凸起圓柱繞流問題進行了數(shù)值模擬計算, 在Re=1.4×105時, 通過對圓柱表面時均壓力系數(shù)、時均摩擦力系數(shù)和間歇因子等與DNS結(jié)果進行比較, 驗證了轉(zhuǎn)捩 SST模型在模擬圓柱繞流問題方面的準確性.之后采用該模型對不同Reynolds數(shù)下的凸起圓柱繞流問題進行了計算分析, 得到了不同Reynolds數(shù)下圓柱近壁區(qū)的流動特性和變化規(guī)律, 可以得到以下結(jié)論:

(1)當來流Reynolds數(shù)處于臨界區(qū)時, 氣流在圓柱上表面凸起處形成了3個反向旋轉(zhuǎn)的漩渦, 并在之后隨著θ的增大, 發(fā)生了流動分離和流動轉(zhuǎn)捩現(xiàn)象.

(2)對于不同Reynolds數(shù)的來流, 圓柱上表面的凸起可以使氣流發(fā)生轉(zhuǎn)捩的位置提前, 并且隨著Reynolds數(shù)的增加, 其轉(zhuǎn)捩位置逐漸后移.

(3)上表面凸起位置的前后產(chǎn)生了3個反向旋轉(zhuǎn)的漩渦, 形成了背風面低壓區(qū).同時, 由于渦的存在, 使渦上方流體的流通通道變窄, 造成其流速增大、壓強降低, 從而導致圓柱產(chǎn)生升力, 隨著來流Reynolds數(shù)的增大, 其升力逐漸變大.