羅倩， 胡常莉

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

0 引言

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

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

本研究基于標(biāo)準(zhǔn)PANS模型中控制參數(shù)fk的取值特點(diǎn)，提出一種分域計(jì)算方法，實(shí)現(xiàn)了同一流場采用不同湍流模型求解的效果，即可對關(guān)注區(qū)域進(jìn)行精細(xì)求解，其他區(qū)域降低求解要求，減少計(jì)算消耗。通過對商業(yè)軟件CFX的二次開發(fā)實(shí)現(xiàn)所提方法，并應(yīng)用其分別進(jìn)行了柱體繞流和繞水翼空化流動的數(shù)值研究，分析了方柱繞流的流場特性和水翼云狀空化流動的空泡形態(tài)和動力特性，并與實(shí)驗(yàn)結(jié)果對比，驗(yàn)證了該方法在單相湍流流動和兩相湍流流動應(yīng)用中的可行性。

1 數(shù)值計(jì)算方法的控制方程

1.1 控制方程

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

(1)

(2)

1.2 標(biāo)準(zhǔn)PANS模型及其分域思想

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

(3)

(4)

(5)

Cμ為黏度相關(guān)系數(shù)，Cμ=0.09；σku、σεu為Prantdl數(shù)，

(6)

σk=1.0,σε=1.3；Pu為湍動能生成項(xiàng)；

(7)

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

(8)

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

(9)

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

1.3 無量綱參數(shù)

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

(10)

(11)

(12)

(13)

(14)

(15)

式中：u∞、p分別為來流流速和當(dāng)?shù)仂o壓強(qiáng)；L、A對應(yīng)的是本文研究對象的特征長度和有效面積；υ為運(yùn)動黏度；ρ為流體密度；FL、FD為研究對象所受的升力和阻力；p∞、pv分別為環(huán)境壓強(qiáng)和飽和蒸氣壓；f為空穴周期變化的頻率。

2 繞二維方柱流動的數(shù)值計(jì)算

2.1 計(jì)算域設(shè)置和網(wǎng)格無關(guān)性驗(yàn)證

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

圖1 方柱繞流計(jì)算域及網(wǎng)格示意圖Fig.1 Computational domain and mesh around the square cylinder

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

圖2 不同網(wǎng)格計(jì)算所得升力系數(shù)Fig.2 Lift coefficients of different meshes

2.2 方柱繞流的結(jié)果討論

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

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

至今為止，已有諸多學(xué)者對較大雷諾數(shù)下的方柱繞流進(jìn)行了實(shí)驗(yàn)和數(shù)值研究，表1列出了實(shí)驗(yàn)研究和DNS方法以及本文計(jì)算方法所得的時均阻力系數(shù)和斯特勞哈數(shù)。通過對比可以看出，本文方法得到的時均阻力系數(shù)與文獻(xiàn)[21-22]中的實(shí)驗(yàn)值和DNS方法結(jié)果相差較小，基本一致，而斯特勞哈數(shù)比實(shí)驗(yàn)值和DNS方法結(jié)果值略大些，誤差在7%～16.4%左右。這是因?yàn)橥牧鞅旧砭哂腥S性，而二維數(shù)值計(jì)算會造成流動特征頻率的過預(yù)測。

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

表1 本文方法與前人實(shí)驗(yàn)和數(shù)值結(jié)果的比較

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

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

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

3 繞二維水翼非定?？栈鲃訑?shù)值計(jì)算

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

3.1 計(jì)算域設(shè)置與方法說明

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

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

圖6 計(jì)算域設(shè)置及網(wǎng)格示意圖Fig.6 2-D computational domain and mesh around the hydrofoil

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

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

表2 二維水翼空泡形態(tài)隨時間變化

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

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

4 結(jié)論

本文基于標(biāo)準(zhǔn)PANS模型中控制參數(shù)fk的取值特點(diǎn)，發(fā)展了一種對流動區(qū)域進(jìn)行分域計(jì)算的方法，應(yīng)用該方法分別計(jì)算了湍流流動中較為經(jīng)典的方柱繞流和繞水翼非定?？栈鲃?，并基于實(shí)驗(yàn)結(jié)果對比分析了該方法的可行性。得出主要結(jié)論如下：

1)基于PANS的分域計(jì)算方法可針對流場不同區(qū)域，通過設(shè)置控制參數(shù)fk的不同取值，實(shí)現(xiàn)對所關(guān)注流動區(qū)域的精細(xì)模擬。

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