扭曲圓柱繞流特性的大渦模擬

郭春雨， 郭航， 胡健

(哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

圓柱繞流一直是流體力學(xué)中的經(jīng)典問題，其涉及到的邊界層轉(zhuǎn)捩、流動分離、再附著及旋渦脫落等物理現(xiàn)象都是流體力學(xué)中的研究難點。圓柱繞流流動既有不固定的分離點,又有分離后的尾流和脫體渦。隨著雷諾數(shù)的增加,尾流性質(zhì),脫體渦的形態(tài)有很大的變化,具有豐富的流動現(xiàn)象。橋梁結(jié)構(gòu)物,如拱橋的吊桿、斜拉橋的拉索,具有圓柱形狀,物體表面邊界層流動在逆壓梯度下分離,在一定的雷諾數(shù)范圍產(chǎn)生規(guī)則的旋渦結(jié)構(gòu),這些交替的旋渦脫落對圓柱受到的升阻力,渦激振動以及產(chǎn)生的噪聲問題都有很大影響，因此對圓柱繞流渦脫落以及振動等問題的研究具有十分重要的意義[1]。

Schewe[2]在密封風(fēng)洞中對從亞臨界至超臨界雷諾數(shù)范圍進行了阻力、升力及斯特勞哈爾數(shù)的測量。在臨界雷諾數(shù)范圍內(nèi)觀察到了2個不連續(xù)的轉(zhuǎn)變，這可以解釋為2個臨界雷諾數(shù)的分叉點。在這2種情況下，這些轉(zhuǎn)變都伴隨著臨界值波動、對稱性破裂(出現(xiàn)穩(wěn)定升力)和遲滯現(xiàn)象。Lam等[3]利用各種實驗測量技術(shù)對雷諾數(shù)為3 000～9 000扭曲圓柱的近尾跡進行了研究，結(jié)果表明扭曲圓柱的平均渦形成長度超過了一半波長且大于光滑圓柱，并且渦形成長度也直接導(dǎo)致阻力以及相關(guān)波動升力的減小。Lam等[4]對雷諾數(shù)為3 000的不同形狀的扭曲圓柱進行了數(shù)值模擬研究，結(jié)果表明由于展長方向波浪形狀的分離線導(dǎo)致扭曲圓柱近尾流渦結(jié)構(gòu)沿展向有周期性的變化。Jung等[5]對雷諾數(shù)為3 000的3種(光滑、波浪形、扭曲形)圓柱體的繞流特性進行了數(shù)值模擬，可以得知扭曲圓柱阻力和升力與光滑圓柱相比分布減小了13%和96%，沿展長方向產(chǎn)生了周期性的橫向渦結(jié)構(gòu)，分離剪切層更加穩(wěn)定且更靠近下游。Zhou等[6]在水槽以及拖曳水池中對雷諾數(shù)為7.4×103～8×104的光滑表面、溝槽狀表面以及洼坑狀表面3種圓柱進行了研究。結(jié)果表明溝槽狀表面以及洼坑狀表面圓柱的阻力及均方根升力系數(shù)比光滑圓柱要更小，而PIV結(jié)果表明圓柱表面的粗糙度導(dǎo)致了圓柱尾流區(qū)域以及渦脫落強度相比于光滑圓柱也要更小。Jie等[7]主要研究了觸角形圓柱與光滑、橢圓以及波浪圓柱體之間的力與渦結(jié)構(gòu)等方面的不同。結(jié)果表明觸角形圓柱在升力方面比普通圓柱減小了79.2%，并在斯特勞哈爾數(shù)為0.2時只觀察到一個相對較低幅度的波峰。

本文計算采用的雷諾數(shù)范圍為亞臨界雷諾數(shù)，這是由于圓柱表面流體分離后尾流是完全湍流的[8]。與上述文獻中的方法不同，本研究采用的是相同直徑不同扭曲角度的扭曲方式，著重于改變圓柱表面形狀以達到減小圓柱受力波動以及渦激振動的目的。本研究通過三維大渦模擬方法討論了亞臨界雷諾數(shù)范圍內(nèi)扭曲圓柱的扭曲角度及旋轉(zhuǎn)角度對其水動力性能的影響。

1 數(shù)值控制方程及模擬方法

1.1 控制方程

本文使用三維大渦模擬湍流模型來對圓柱繞流進行數(shù)值模擬[9]。大渦模擬的基本思想是直接計算大尺度渦，而通過亞格子模型來對小渦進行模擬。與時間平均不同，大渦模擬方法通過空間濾波操作來分成大渦以及小渦。因此濾波后不可壓縮流體的納維-斯托克斯以及連續(xù)方程即為大渦模擬控制方程：

(1)

(2)

通常情況下亞格子模型是基于渦粘模型的，它基于人工渦流粘度的方法，其中湍流的影響被集中到湍流粘度。該方法將亞格子尺度下的動能耗散視為類似于分子擴散。因此τij的表達形式為：

(3)

(4)

本文采用的亞格子渦粘模型為WALE(wall-adapting local-eddy viscosity)模型[10]，它是一個更現(xiàn)代的亞格子尺度模型，在其公式中使用了速度梯度張量的一種新形式。WALE亞格子模型提供的混合長度形式的亞格子渦粘系數(shù)為：

νt=ρΔ2Sw

(5)

式中：ρ為流體密度；Δ為長度尺度或網(wǎng)格過濾寬度[11]；Sw為扭曲參數(shù)。長度尺度是通過網(wǎng)格體積來定義的：

(6)

式中：Cw為模型系數(shù)，取值為0.544；κ為馮卡門常數(shù)，取值為0.41。扭曲參數(shù)Sw定義為：

(7)

(8)

1.2 數(shù)值方法

本文模擬采用的是三維、不可壓縮、分離流、恒定密度以及非定常模擬，通過將有限體積法應(yīng)用于非結(jié)構(gòu)化網(wǎng)格上來對納維-斯托克斯方程進行求解。在大渦模擬中，應(yīng)用SIMPLE算法對壓力速度耦合方程進行求解，空間離散采用的是有界中心差分離散格式[12]，時間項采用的是二階隱式時間離散格式，梯度計算采用的是Hybrid Gauss-LSQ方法。

2 計算模型及網(wǎng)格和邊界條件

2.1 計算模型

本模擬選擇了4種不同的圓柱進行數(shù)值模擬，其中包括光滑圓柱以及3種不同扭曲角度的扭曲圓柱。3種圓柱的直徑和展長相同，分別為D和6D，其中扭曲圓柱是由3段波長λ=2D的相同圓柱組成。扭曲圓柱的形狀是通過扭曲角度α來定義的，而α為圓柱沿展長方向截面繞圓周上一點旋轉(zhuǎn)的角度，并規(guī)定沿逆時針方向旋轉(zhuǎn)為正，如圖1(a)所示，且α與截面展向坐標(biāo)z具有以下關(guān)系：

(9)

式中αmax為最大扭曲角度，本模擬選取的扭曲圓柱A1、A2、A3分別對應(yīng)為10°、25°、40°。

由于扭曲圓柱不具有光滑圓柱那樣的幾何對稱性，因此需要對其進行旋轉(zhuǎn)不同的角度進行計算。旋轉(zhuǎn)角度θ定義為圓柱繞其中心軸旋轉(zhuǎn)的角度，本文選取了θ=0°，45°，90°，135°，180°，225°，270°，315°共8個角度。為了更好地對圓柱繞流特性進行分析，在展向選擇了2個截面，節(jié)面和鞍面如圖1(b)所示。

圖1 扭曲圓柱Fig.1 The twisted cylinder

2.2 網(wǎng)格和邊界條件

計算域網(wǎng)格劃分采用Star CCM+生成的切割體網(wǎng)格，通過使用幾何體曲面剪切六面體模板網(wǎng)格來生成體積網(wǎng)格。圖2(a)展示的是扭曲圓柱計算域及圓柱邊界層附近的體網(wǎng)格，為了滿足大渦模擬計算要求，近壁面y+取值小于1，為了更好地捕捉圓柱尾流區(qū)的渦結(jié)構(gòu)，對下游方向的網(wǎng)格進行了加密，且邊界層處的網(wǎng)格尺寸以一個合適的比率向加密區(qū)逐漸增大。計算采用的無量綱時間步長定義為：T=U∞Δt/D，對于目前的所有模擬，計算了至少500個無量綱時間步長，對應(yīng)大約100個旋渦脫落周期以獲得更可靠的統(tǒng)計信息，這也要大于Lam[4]計算的時間步長以及脫落周期。

計算域和邊界條件如圖2(b)所示。計算域的尺寸在固定笛卡爾坐標(biāo)系下在x、y、z軸方向上分別為24D×16D×6D，其中x軸方向上游邊界距離原點為8D，下游邊界距離原點為16D，原點距離y軸方向上下2個邊界距離均為8D，展向長度則與圓柱展長相同。坐標(biāo)系原點位于圓柱底部圓心上，x軸與入口來流方向相同(順流方向)，z軸與圓柱中心軸方向平行，y軸則與x軸和z軸垂直(橫流方向)。本模擬基于圓柱直徑D和來流速度U∞的雷諾數(shù)Re=U∞D(zhuǎn)/ν=28 712。計算域的邊界條件設(shè)置如下：入口邊界設(shè)置為速度入口，出口邊界設(shè)置為壓力出口，由于考慮是無限長圓柱結(jié)構(gòu)，因此其他4個邊界均設(shè)置為對稱平面，圓柱則設(shè)置為無滑移壁面。

圖2 數(shù)值模擬設(shè)定Fig.2 Schematic of numerical simulation setting

3 結(jié)果與分析

3.1 網(wǎng)格及時間步長無關(guān)性驗證

(10)

(11)

式中：ρ為流體密度；U∞為來流速度；D為圓柱截面直徑；L為圓柱展向長度；FD和FL分別為總阻力和升力。斯特勞哈爾數(shù)為無量綱化的渦脫落頻率(fs)，其表達式為：

(12)

式中渦脫落頻率fs是通過對升力系數(shù)的時歷曲線作快速傅里葉變換得到的。

首先對光滑圓柱進行網(wǎng)格無關(guān)性驗證，表1展示了3套不同數(shù)量的網(wǎng)格，這里Si(i=1,2,3)分別代表粗、中和細(xì)網(wǎng)格對應(yīng)的計算結(jié)果。以一種類似于文獻[13]中處理非結(jié)構(gòu)化網(wǎng)格的方式，網(wǎng)格加細(xì)比是需要在進行網(wǎng)格無關(guān)性驗證之前確定的一個重要參數(shù)，網(wǎng)格加細(xì)比rG定義為：

(13)

對光滑圓柱進行時間步長無關(guān)性驗證，表2展示了3個不同無量綱化時間步長，這里Si(i=4,5,6)分別代表小、中和大時間步長對應(yīng)的計算結(jié)果。本文3個不同時間步長網(wǎng)格數(shù)量均為435萬的粗網(wǎng)格。結(jié)果表明，只有S6與實驗結(jié)果相差較大，而S4和S5則與實驗數(shù)據(jù)吻合較好。因此，在保證計算精度和效率的前提下，本文最終選擇了網(wǎng)格數(shù)量為435萬的粗網(wǎng)格，而無量綱時間步長則為T=0.02。

表1 光滑圓柱網(wǎng)格無關(guān)性驗證Table 1 Grid independence verification on smooth cylinder

表2 光滑圓柱時間步長無關(guān)性驗證

表3 扭曲圓柱A1網(wǎng)格無關(guān)性驗證

表4 扭曲圓柱A1時間步長無關(guān)性驗證

3.2 力系數(shù)和斯特勞哈爾數(shù)

為了節(jié)省大渦模擬所需的計算資源以及時間，對扭曲圓柱不同角度之間的相似性進行驗證。以α=10°的扭曲圓柱A1為例，選取了8個旋轉(zhuǎn)角度(θ=0°，45°，90°，135°，180°，225°，270°，315°)進行計算，計算結(jié)果如表5所示，從計算結(jié)果可以看出，扭曲圓柱間隔角度180°的力系數(shù)以及斯特勞哈爾數(shù)近似相同，可以初步判斷扭曲圓柱不同角度之間具有部分的相似性。

表5 扭曲圓柱A1不同旋轉(zhuǎn)角度的力系數(shù)和St

圖3 3種扭曲圓柱與光滑圓柱的對比Fig.3 Comparison of three twisted cylinders and the smooth cylinder

3.3 尾渦形成長度

圓柱尾渦形成長度是一個很重要的量，因為它影響圓柱尾部壓力以及力系數(shù)。對于尾渦形成長度的定義不是普適的，不同學(xué)者有不同的方法[15-19]。通常用尾流中心線(y=0平面)上無量綱化的平均流向速度為零的時間平均閉合點位置(PU0)以及無量綱化的流向速度波動均方根的最大值點位置(PUrms)來定義渦形成長度。因此尾渦形成長度包括平均流向速度為零(U/U∞=0)對應(yīng)的尾流閉合長度Lfc和流向速度波動均方根(u′/U∞)最大值對應(yīng)的最大湍流強度長度Lfu。

圖4為3種扭曲圓柱沿展長方向一個波長長度的尾流閉合長度Lfc和最大湍流強度長度Lfu與光滑圓柱的對比。結(jié)果表明，圓柱A1的Lfc和Lfu與光滑圓柱相比幾乎沒有增長。當(dāng)θ=0°時，扭曲圓柱A2的Lfc和Lfu與光滑圓柱相比最大分別增大了14%和24%，而在其他旋轉(zhuǎn)角度Lfc增大幅值達到23%～96%，Lfu增大幅值達到31%～105%。當(dāng)θ=0°時，扭曲圓柱A3的Lfc和Lfu與光滑圓柱相比最大分別增大了22%和27%，相比光滑圓柱有較大增幅，在其他旋轉(zhuǎn)角度Lfc最大增幅為98%，Lfu增大幅值達到39%～114%。

綜合圖4分析可知，Lfc和Lfu在θ=0°時波動幅值較小，而在其他旋轉(zhuǎn)角度時沿展長方向基本呈現(xiàn)增大的趨勢，并且隨著圓柱扭曲角度α的增大，Lfc和Lfu的波動幅值也逐漸增大，特別的扭曲圓柱A3甚至出現(xiàn)了Lfc的最小值小于光滑圓柱的情況，原因可能為隨著圓柱表面曲率逐漸增大而導(dǎo)致尾流閉合長度的波動范圍增大。當(dāng)θ=0°時，Lfc和Lfu隨著α的增大也逐漸增大；當(dāng)θ為其他值時，Lfc和Lfu隨著α的增幅變化不是很大，并且均大于θ=0°時的Lfc和Lfu的增幅。通過對比可以得出結(jié)論扭曲圓柱升阻力隨著渦形成長度Lfc和Lfu的增大而減小，當(dāng)θ=0°時,渦尾流長度相對于光滑圓柱增幅比較小，因此對圓柱尾流渦結(jié)構(gòu)的控制以及抑制幅度也比較小，使得升阻力下降也比較??；相反其他角度尾流長度得到很大的延伸表明能夠?qū)ξ矞u結(jié)構(gòu)很好地控制，從而能夠使升阻力波動最小化。

3.4 壓力系數(shù)分布

圖4 扭曲圓柱與光滑圓柱的尾流閉合長度Lfc和最大湍流強度長度Lfu對比Fig.4 Comparison of the wake closure length Lfc and maximum turbulence intensity length Lfu of twisted cylinders and the smooth cylinder

圖5 扭曲圓柱與光滑圓柱的平均圓周壓力系數(shù)對比Fig.5 Comparison of the average circumferential pressure coefficient of twisted cylinders and the smooth cylinder

3.5 湍流動能

圖6 在x-y平面的無量綱化湍流動能TKE分布對比Fig.6 Comparison of the normalized turbulent kinetic energy (TKE) distributions in the x-y plane

從圖6(a)可以看出，光滑圓柱最大TKE區(qū)間位于圓柱近尾流區(qū)，且該區(qū)間內(nèi)TKE值很大。對比圖6(b)～(e)可以發(fā)現(xiàn)，當(dāng)θ=0°時，扭曲圓柱A2在節(jié)面和鞍面與光滑圓柱相比TKE的分布基本相同。而當(dāng)θ=90°時，圓柱A2的最大TKE區(qū)間相比θ=0°時有較大的后移，同時區(qū)間內(nèi)的TKE值有大幅度的減小，并且鞍面的TKE值要比節(jié)面更小。根據(jù)圖6(f)～(i)可以看出，當(dāng)θ=0°時，圓柱A3的最大TKE區(qū)間相比A2也有后移，區(qū)間內(nèi)的TKE值也有大幅度的減小。當(dāng)θ=90°時，圓柱A3的最大TKE區(qū)間位置與圓柱A2基本相同，區(qū)間內(nèi)的TKE值略有減小，并且鞍面的TKE值仍然比節(jié)面要小。

綜合圖6分析可知，隨著α的增大，圓柱的最大值TKE區(qū)間位置逐漸后移，遠離圓柱表面，沿流向的TKE分布變化逐漸變小，且最大及整體TKE值相比光滑圓柱也在逐漸減小，尤其在α=40°時有很大程度的減幅。而隨著α的增大，鞍面的TKE值相比于節(jié)面也越來越小。當(dāng)θ=90°時整體的TKE值均小于θ=0°，這也與3.2節(jié)中升阻力系數(shù)大小以及3.3節(jié)中尾渦形成長度曲線相一致，因此可以得出結(jié)論，圓柱尾流TKE的顯著減小能夠?qū)е聢A柱的阻力以及波動升力的減小。

3.6 三維渦結(jié)構(gòu)

圖7展示了不同圓柱的三維瞬時渦結(jié)構(gòu)等值面，對比了自由剪切層發(fā)展和渦形成長度的顯著差異。

圖7 三維瞬時渦結(jié)構(gòu)等值面對比Fig.7 Comparison of the three-dimensional instantaneous vortex structure isosurfaces

圖7(a)為光滑圓柱的三維渦結(jié)構(gòu)云圖，可以觀察到在圓柱近尾流區(qū)域存在大尺度的卡門渦街，并且振蕩幅度較大。通過觀察圖7(b)～(g)能夠發(fā)現(xiàn)，圓柱A1尾流區(qū)域渦振蕩幅值與光滑圓柱相比略有減小。當(dāng)θ=0°時，圓柱A2尾流區(qū)域渦分布情況與光滑圓柱基本相同，而隨著α的增大，圓柱尾渦的脫落位置延遲得更加靠近下游，同時還能觀察到卡門渦街在遠尾流區(qū)域重新出現(xiàn)，但是其波動幅值較光滑圓柱減小了。當(dāng)θ=90°時，其他2種圓柱剪切層沿下游方向延伸得更長了，并且渦的脫落位置也更加遠離圓柱表面。同時由于初始渦結(jié)構(gòu)的形成沿下游方向更加延遲了，因此卡門渦街強度沿下游方向也更加弱化了。

綜合對圖7進行分析可知，與光滑圓柱相比，隨著α的增大，3種圓柱尾渦脫落位置以及剪切層長度沿順流方向的延伸逐漸增加，卡門渦街在更遠離圓柱表面的下游位置重新出現(xiàn)，但是其波動幅值也逐漸減小。并且可以看到，θ=90°時渦的抑制效果均比θ=0°時要更好，延伸長度也更大，這也與3.2節(jié)中力系數(shù)以及3.5節(jié)中TKE的分布情況相吻合。由于圓柱剪切層的延伸伴隨著較弱的尾渦波動，從而直接導(dǎo)致了阻力以及升力振蕩幅值的減小。

4 結(jié)論

1)扭曲圓柱A3相比光滑圓柱阻力以及升力系數(shù)減小幅度比較大，斯特勞哈爾數(shù)基本沒有變化。

2)尾流長度得到很大的延伸，同時最小壓力系數(shù)位置和大小也增大了。

3)圓柱尾流湍動能顯著減小，三維渦結(jié)構(gòu)的波動也得到了抑制。

因此扭曲圓柱A3對尾渦結(jié)構(gòu)有很好的控制以及使自由剪切層更加穩(wěn)定的發(fā)展，并且也為圓柱提供了更大的阻力減額以及更好的抑制升力波動，能夠顯著減小圓柱的渦激振動現(xiàn)象。