基于局部時均Navier-Stokes模型的分域計算方法

羅倩， 胡常莉

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

0 引言

隨著計算機技術日新月異的發展，借助數值計算方法研究流體力學相關問題，越來越受到人們的青睞。尤其是工程實際中的復雜流動問題，比如各類流體機械、船舶、水中兵器等工程領域，往往涉及湍流、相變、可壓縮等多個物理過程，這對數值模型的發展提出了很大的挑戰。其中，湍流模型直接關乎復雜流動數值方法的預測精度。

一直以來，方柱繞流與空化流動是流體力學研究中的經典與熱點問題，備受研究者的關注，建立合適的湍流模擬方法也是其數值研究的重點[1-2]。張顯雄等[3]對比了5種湍流渦黏模型在方柱繞流的數值模擬中的求解差異，結果表明：標準k-ε(k為湍動能，ε為湍流耗散率)湍流模型的計算結果程度整體弱于其余湍流模型；剪切應力輸運(SST)k-ω(ω為湍流頻率)模型的計算結果優于其余湍流模型。Long等[4]、Ji等[5]、楊龍等[6]應用大渦模擬(LES)方法進行了空化流動的數值計算，發現LES方法可以準確模擬出非定?？栈鲃拥牧鲃犹匦?。但LES方法計算資源消耗大。隨著數值研究的深入以及對湍流計算要求的提高，許多混合湍流模型被提出和應用[7-8]。局部時均Navier-Stokes(PANS)模型是由Girimaji[9]提出的，一種可從雷諾時均Navier-Stokes(RANS)平滑過渡到直接數值模擬(DNS)的混合湍流模型，被廣泛應用于流動研究中。Razi等[10]應用PANS模型研究了周期性山型流道內的流動分離問題。王亦曉等[11]采用PANS模型研究了簸箕形進水流道的水力特性。Busco等[12]采用PANS模型研究了壓水堆燃料束間隔格柵間的復雜湍流結構。劉躍等[13]應用該模型研究了繞圓柱流動過程，結果表明PANS模型可以捕捉到豐富的湍流結構。該模型在空化流動特性的研究中，也有較好的模擬效果[14-15]。但標準PANS模型對整個流場是一致性求解，未能體現橋接模型的優勢。因此，眾多學者對標準PANS模型作了進一步的發展與修正。Luo等[16]提出一種基于k-ω修正的PANS模型，開展了繞NACA66水翼的空化流動，討論了濾波器參數fk不同取值對瞬態空化渦流動預測的準確性。Zhang等[17]采用PANS模型研究了最大密度比和控制參數fk對預測水翼空化流動的影響，發現提高最大密度比和降低fk值可以獲得更好的結果。Huang等[18]根據局部網格大小和湍流長度對fk進行調整，實現在空間和時間上改變fk分布，使得PANS模型對空化湍流流動的預測能力獲得了明顯的提高。Hu等[19]考慮密度變化特點發展了一種修正PANS模型，依據當地云狀空化狀態對控制參數進行動態調節，獲得精確結果的同時提升了計算效率。

本研究基于標準PANS模型中控制參數fk的取值特點，提出一種分域計算方法，實現了同一流場采用不同湍流模型求解的效果，即可對關注區域進行精細求解，其他區域降低求解要求，減少計算消耗。通過對商業軟件CFX的二次開發實現所提方法，并應用其分別進行了柱體繞流和繞水翼空化流動的數值研究，分析了方柱繞流的流場特性和水翼云狀空化流動的空泡形態和動力特性，并與實驗結果對比，驗證了該方法在單相湍流流動和兩相湍流流動應用中的可行性。

1 數值計算方法的控制方程

1.1 控制方程

本研究的基本控制方程為連續性方程和Navier-Stokes方程：

(1)

(2)

1.2 標準PANS模型及其分域思想

PANS模型的湍動能ku和耗散率εu的輸運方程分別為

(3)

(4)

(5)

Cμ為黏度相關系數，Cμ=0.09；σku、σεu為Prantdl數，

(6)

σk=1.0,σε=1.3；Pu為湍動能生成項；

(7)

Cε1=1.44，Cε2=1.92.

(8)

PANS模型的兩個控制參數[20]分別定義為

(9)

在高雷諾數的流動中fε值通常取1：當fk=1時，說明湍流控制方程復原到RANS模型；當fk=0時，表示數值計算過程沒有湍流模型的引入，為直接求解的方式。PANS模型可以實現任何濾波器尺度對湍流流動的求解[9]，fk的取值控制濾波尺度大小，隨著fk的減小，PANS模型可釋放更多的湍流運動尺度。根據這一特點，本研究將整個計算域化分為多個區域來計算，以實現關注流域fk取小值進行計算，其他區域fk取較大值進行計算。

1.3 無量綱參數

本研究涉及的無量綱數主要有雷諾數Re、升力系數CL、阻力系數CD、壓力系數Cp和空化數σ，以及對升力系數隨時間變化曲線進行快速傅里葉變換(FFT)所得對應的斯特勞哈數St，定義分別為

(10)

(11)

(12)

(13)

(14)

(15)

式中：u∞、p分別為來流流速和當地靜壓強；L、A對應的是本文研究對象的特征長度和有效面積；υ為運動黏度；ρ為流體密度；FL、FD為研究對象所受的升力和阻力；p∞、pv分別為環境壓強和飽和蒸氣壓；f為空穴周期變化的頻率。

2 繞二維方柱流動的數值計算

2.1 計算域設置和網格無關性驗證

圖1(a)為繞方柱流場計算域示意圖。方柱邊長D為0.05 m，計算域長為50D，下邊界與方柱距離為12D. 計算流體為25 ℃的水，采用速度入口、壓力出口，對應雷諾數Re=22 000，流域上下邊界和方柱壁面為無滑移壁面。圖1(a)中的紅色虛線將整個計算域劃分為方柱近流區和遠流區兩部分進行計算，即近流區fk取值為0，遠流區fk取值為1.

圖1 方柱繞流計算域及網格示意圖Fig.1 Computational domain and mesh around the square cylinder

圖1(b)為方柱周圍的網格分布情況，計算域采用結構化網格，并對方柱近壁面區域進行網格加密處理。圖2給出了網格1(節點數240 160)、網格2(節點數304 974)、網格3(節點數525 352)和網格4(節點數973 824)4套網格下的數值計算結果。由圖2可以看出，隨著網格數目的增加，不同網格之間對應的方柱升力系數變化幅度減小，而網格數量的增加必然會加大計算消耗。因此，綜合考慮計算的適用性和經濟性，本文選用網格2進行繞方柱流動的數值研究。

圖2 不同網格計算所得升力系數Fig.2 Lift coefficients of different meshes

2.2 方柱繞流的結果討論

圖3(a)為繞方柱流動過程中的瞬時渦量圖，流體流經方柱時，方柱上下表面交替產生旋渦，旋渦隨流動發展脫落，小渦變大渦向下游運動形成卡門渦街。圖3(b)為時均流向速度u云圖，從圖中可以看出，由于方柱邊角發生流動分離，產生分離渦，方柱近壁面區域存在較大的速度梯度，且受卡門渦街的影響，方柱正后方流場區為低速區，其速度明顯小于兩側流場的速度。

圖3 渦量和速度云圖Fig.3 Contours of vorticity and velocity-u

至今為止，已有諸多學者對較大雷諾數下的方柱繞流進行了實驗和數值研究，表1列出了實驗研究和DNS方法以及本文計算方法所得的時均阻力系數和斯特勞哈數。通過對比可以看出，本文方法得到的時均阻力系數與文獻[21-22]中的實驗值和DNS方法結果相差較小，基本一致，而斯特勞哈數比實驗值和DNS方法結果值略大些，誤差在7%～16.4%左右。這是因為湍流本身具有三維性，而二維數值計算會造成流動特征頻率的過預測。

圖4(a)和圖4(b)分別為方柱表面時均壓力系數分布和中心線上時均流向速度沿x軸方向的分布。其中圖4(a)中橫坐標表示位置，與圖中方柱各邊對應，為了便于展示，時均流向速度均以初始的流體速度u0作無量綱處理?？梢钥闯?，本研究計算所得的方柱表面壓力系數的變化規律與實驗和DNS方法結果基本一致，即方柱來流方向中心點處所受壓力最大，方柱邊角產生流動分離導致壓力顯著降低。從圖4(b)中可看出，本文方法和DNS方法中心線上的時均流向速度均為由0 m/s減小至最小值，隨后漸進增加的變化趨勢，與實驗結果較為一致。然而對比發現，本文方法獲得的中心線上時均流向速度最小值與實驗值更為接近，且在x/D<3范圍，即流動較為強烈的區域，本文方法時均流向速度的變化相比DNS方法結果與實驗結果更吻合。

表1 本文方法與前人實驗和數值結果的比較

圖4 時均壓力系數和流向速度分布情況Fig.4 Distribution of time-averaged pressure coefficient and velocity-u

圖5為不同截面處時均流向速度沿y軸方向的分布，對應的截面位置如圖3(b)中黑色虛線所示。與實驗結果對比可看出，數值結果與實驗結果較為吻合。不同截面位置的時均流向速度分布特點均是由中心位置沿y軸方向漸進增大。x/D=0處，即方柱中心位置，由于流體流過方柱邊角位置發生流動分離，產生分離渦，導致方柱近壁面速度最小且為負值，在剪切層處速度為最大值。而x/D=1，如圖3(b)中黑色虛線所示，該處對應的是方柱繞流尾流的回流區域，此處的時均流向速度從負值開始漸進增大。隨著向下游流動，逐漸遠離回流區的尾流位置，如圖3(b)中x/D=3、x/D=8位置處，其剖面處中心位置的最小流向速度均已由負值增大到正值，流向速度沿y軸方向漸進增加至遠流場基本保持平緩。

圖5 不同剖面時均流向速度u的對比Fig.5 Comparison of time-averaged velocity-u at different sections

3 繞二維水翼非定?？栈鲃訑抵涤嬎?/h2>

目前，對于非定?？栈鲃拥臄抵涤嬎惴椒?，基于均質流框架并耦合Zwart等[25]空化模型和湍流模型的方法被廣泛應用，可詳見文獻[7,18,26-28]。本文也采用上述方法耦合分域計算的PANS模型對繞Clark-y水翼的非定?？栈鲃舆M行了數值計算。

3.1 計算域設置與方法說明

圖6(a)和圖6(b)分別給出了計算域設置和局部網格示意圖，計算域的選擇與實驗段[29]相同。Clark-y翼型的弦長C為0.07 m，計算域長為10C，邊界條件采用速度入口，壓力出口，流速為u∞=10 m/s，對應的雷諾數Re=7×105，調節出口壓力設定空化數σ=0.8，流動區域上下邊界為自由滑移壁面，翼型表面采用絕熱、無滑移壁面邊界。計算域網格采用結構網格劃分，對翼型近壁面采取了網格加密處理，近壁面y+值在5左右，滿足壁面函數要求。

PANS模型在繞水翼云狀空化流動模擬中控制參數fk的分域，如圖6(a)所示。將計算域分為3個區域：紅色虛線標注的近壁面區域E1，其控制參數取fk=0；遠流場區域E3，其控制參數取fk=1；E2為過渡區域，采用前人研究使用較多的fk=0.2模型進行計算。

圖6 計算域設置及網格示意圖Fig.6 2-D computational domain and mesh around the hydrofoil

3.2 繞水翼非定?？栈鲃拥慕Y果討論

表2給出了繞Clark-y水翼云狀空化流動一個周期T內空泡形態的演變過程，其中水翼前端附著空泡剛開始產生時為t0時刻。對比發現，數值計算與實驗結果較一致地捕捉到附著空泡從產生到脫落的準周期變化過程，即水翼前端開始產生附著空泡。附著空泡隨時間沿流動方向發展至最大，在回射流作用下，附著空泡逐漸出現斷裂脫落，在t0+0.7T時刻，水翼附著空泡呈現出一種藕斷絲連的狀態，斷裂后空泡隨流動向下游運動，脫落空泡形態也由沿展向方向的扁長型發展成厚度較大的團狀型。

表2 二維水翼空泡形態隨時間變化

圖7(a)和圖7(b)分別為Clark-y水翼升力系數曲線及其功率密度譜分析。與實驗結果對比發現，數值計算可以很好地模擬出水翼升力系數曲線隨時間的波動規律，與空泡形態發展過程相對應而呈周期性變化，另外，二者的時均升力系數也基本一致。將水翼升力系數通過FFT獲得功率譜密度圖，圖7中功率密度最大值反映了特征頻率，對比發現，數值計算所得斯特勞哈數(St=0.25)是略大于實驗值(St=0.22)的，這是由于實驗中存在一定的壁面效應會導致空化在非定常演變過程中產生三維U形空泡團的斷裂脫落，從而采用二維模型的數值計算易過預測其特征頻率。

圖7 翼型升力系數曲線及功率密度譜分析Fig.7 Lift coefficient and power density spectrum analysis of hydrofoil

4 結論

本文基于標準PANS模型中控制參數fk的取值特點，發展了一種對流動區域進行分域計算的方法，應用該方法分別計算了湍流流動中較為經典的方柱繞流和繞水翼非定?？栈鲃?，并基于實驗結果對比分析了該方法的可行性。得出主要結論如下：

1)基于PANS的分域計算方法可針對流場不同區域，通過設置控制參數fk的不同取值，實現對所關注流動區域的精細模擬。

2)對于較大雷諾數下的繞方柱單相流動，基于PANS的分域計算方法所得的時均阻力系數、方柱表面壓力系數和流向速度等，與DNS和實驗結果均基本吻合；對于繞Clark-y水翼的非定?？栈鲃?，該方法可獲得與實驗較一致的空穴形態演變過程和水翼升力系數的變化規律。