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

郭春雨， 郭航， 胡健

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

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

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

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

1 數值控制方程及模擬方法

1.1 控制方程

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

(1)

(2)

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

(3)

(4)

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

νt=ρΔ2Sw

(5)

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

(6)

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

(7)

(8)

1.2 數值方法

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

2 計算模型及網格和邊界條件

2.1 計算模型

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

(9)

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

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

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

2.2 網格和邊界條件

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

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

圖2 數值模擬設定Fig.2 Schematic of numerical simulation setting

3 結果與分析

3.1 網格及時間步長無關性驗證

(10)

(11)

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

(12)

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

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

(13)

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

表1 光滑圓柱網格無關性驗證Table 1 Grid independence verification on smooth cylinder

表2 光滑圓柱時間步長無關性驗證

表3 扭曲圓柱A1網格無關性驗證

表4 扭曲圓柱A1時間步長無關性驗證

3.2 力系數和斯特勞哈爾數

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

表5 扭曲圓柱A1不同旋轉角度的力系數和St

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

3.3 尾渦形成長度

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

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

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

3.4 壓力系數分布

圖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 扭曲圓柱與光滑圓柱的平均圓周壓力系數對比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區間位于圓柱近尾流區，且該區間內TKE值很大。對比圖6(b)～(e)可以發現，當θ=0°時，扭曲圓柱A2在節面和鞍面與光滑圓柱相比TKE的分布基本相同。而當θ=90°時，圓柱A2的最大TKE區間相比θ=0°時有較大的后移，同時區間內的TKE值有大幅度的減小，并且鞍面的TKE值要比節面更小。根據圖6(f)～(i)可以看出，當θ=0°時，圓柱A3的最大TKE區間相比A2也有后移，區間內的TKE值也有大幅度的減小。當θ=90°時，圓柱A3的最大TKE區間位置與圓柱A2基本相同，區間內的TKE值略有減小，并且鞍面的TKE值仍然比節面要小。

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

3.6 三維渦結構

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

圖7 三維瞬時渦結構等值面對比Fig.7 Comparison of the three-dimensional instantaneous vortex structure isosurfaces

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

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

4 結論

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

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

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

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