基于環(huán)境源項(xiàng)法的SST γ-Reθt轉(zhuǎn)捩模型修正

葉長(zhǎng)亮 王福軍 湯 遠(yuǎn) 陳 俊 鄭 源

(1.河海大學(xué)能源與電氣學(xué)院，南京 211100；2.中國(guó)農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京 100083；3.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海 200240)

0 引言

軸流泵具有大流量、低揚(yáng)程的特點(diǎn)，在農(nóng)業(yè)灌溉領(lǐng)域應(yīng)用廣泛[1-2]。水翼作為軸流泵葉片的簡(jiǎn)化模型，其繞流特性對(duì)于研究軸流泵葉片的繞流分布具有參考意義。文獻(xiàn)[3-5]通過(guò)研究水翼分離流動(dòng)特性為軸流泵內(nèi)部流動(dòng)特性提供了借鑒，然而，目前對(duì)軸流泵葉片附近的轉(zhuǎn)捩流動(dòng)仍缺乏認(rèn)識(shí)。

現(xiàn)有的針對(duì)水力機(jī)械的轉(zhuǎn)捩研究，大多采用數(shù)值模擬方法[6-7]。文獻(xiàn)[8]在研究噴水推進(jìn)軸流泵特性時(shí)，發(fā)現(xiàn)采用轉(zhuǎn)捩模型不僅提高了外特性的預(yù)測(cè)精度，在葉片表面的壓力分布和摩阻系數(shù)分布上也與試驗(yàn)值誤差較小。文獻(xiàn)[9-10]在船用導(dǎo)管螺旋槳研究中發(fā)現(xiàn)，采用轉(zhuǎn)捩模型能成功地捕捉到層流狀態(tài)到湍流狀態(tài)的轉(zhuǎn)換。文獻(xiàn)[11]也采用相同的計(jì)算策略研究了導(dǎo)管螺旋槳上的轉(zhuǎn)捩流動(dòng)，得到的結(jié)論也驗(yàn)證了考慮轉(zhuǎn)捩的重要性。由于轉(zhuǎn)捩流動(dòng)的復(fù)雜性，針對(duì)不同的問(wèn)題，往往需要結(jié)合流動(dòng)特征對(duì)轉(zhuǎn)捩預(yù)測(cè)模型進(jìn)行修正。例如文獻(xiàn)[12]采用新的關(guān)聯(lián)函數(shù)直接建立轉(zhuǎn)捩區(qū)域長(zhǎng)度參數(shù)、臨界動(dòng)量雷諾數(shù)與當(dāng)?shù)亓鲌?chǎng)的關(guān)系，修正了SSTγ-Reθt轉(zhuǎn)捩模型，提高了水翼繞流計(jì)算精度；文獻(xiàn)[13]在關(guān)聯(lián)函數(shù)中引入了雷諾數(shù)和馬赫數(shù)，以體現(xiàn)高超聲速流動(dòng)的可壓縮性；文獻(xiàn)[14]建立了橫流轉(zhuǎn)捩關(guān)聯(lián)函數(shù)，使得γ-Reθt轉(zhuǎn)捩模型具備預(yù)測(cè)橫流轉(zhuǎn)捩的能力；為了提高對(duì)大曲率水翼繞流邊界層轉(zhuǎn)捩的預(yù)測(cè)能力，文獻(xiàn)[15-16]提出的轉(zhuǎn)捩起始位置關(guān)聯(lián)函數(shù)應(yīng)用于γ-Reθt轉(zhuǎn)捩模型并取得了較好的效果。

綜合分析發(fā)現(xiàn)，由γ-Reθt轉(zhuǎn)捩模型和SSTk-ω湍流模型耦合的計(jì)算模型(本文稱作SSTγ-Reθt轉(zhuǎn)捩模型)廣泛應(yīng)用于邊界層轉(zhuǎn)捩預(yù)測(cè)中，具有計(jì)算效率高、適用轉(zhuǎn)捩類型多的特點(diǎn)，本文采用SSTγ-Reθt轉(zhuǎn)捩模型為核心數(shù)值模型預(yù)測(cè)邊界層轉(zhuǎn)捩流動(dòng)。對(duì)于軸流泵而言，其內(nèi)部流動(dòng)雷諾數(shù)高，邊界層轉(zhuǎn)捩流動(dòng)特征復(fù)雜，需結(jié)合流場(chǎng)特性關(guān)聯(lián)SSTγ-Reθt轉(zhuǎn)捩模型的預(yù)測(cè)性能。采用廣泛用于類比軸流泵葉片的NACA0009、Donaldson修型尾緣水翼和NACA0016水翼等，結(jié)合高雷諾數(shù)繞流邊界層轉(zhuǎn)捩特性開展SSTγ-Reθt轉(zhuǎn)捩模型預(yù)測(cè)性能研究，并基于流場(chǎng)特性開展轉(zhuǎn)捩模型的修正，以期為高雷諾數(shù)葉片繞流邊界層轉(zhuǎn)捩預(yù)測(cè)提供參考。

1 數(shù)值計(jì)算方法

1.1 模型建立與網(wǎng)格劃分

文獻(xiàn)[17]對(duì)NACA0009鈍尾緣水翼繞流特性進(jìn)行了試驗(yàn)研究，并獲得了邊界層流場(chǎng)特性。為避免大范圍的水流旋轉(zhuǎn)，降低來(lái)流的紊流程度，在試驗(yàn)段上游增加了特殊裝置使來(lái)流更加穩(wěn)定；試驗(yàn)段測(cè)試的水翼攻角為0°；入口距水翼前緣0.5L(L表示水翼弦長(zhǎng))，出口距水翼后緣6L。試驗(yàn)結(jié)果表明，由于水翼的影響，進(jìn)口各位置的流線方向平均速度與整個(gè)截面的流線方向平均速度存在一定的差異(在2%以內(nèi))。不同雷諾數(shù)下，試驗(yàn)段進(jìn)口湍流度約為1%，水翼邊界層具有明顯的層流邊界層區(qū)同時(shí)觀測(cè)到自然轉(zhuǎn)捩過(guò)程。試驗(yàn)測(cè)量的水翼邊界層流場(chǎng)結(jié)果可為數(shù)值計(jì)算結(jié)果提供參考，從而分析不同湍流模型對(duì)水翼邊界層轉(zhuǎn)捩預(yù)測(cè)的影響。

1.2 網(wǎng)格無(wú)關(guān)性驗(yàn)證

文獻(xiàn)[18]開展了SSTγ-Reθt在翼型流場(chǎng)應(yīng)用中網(wǎng)格的適用性研究，就水翼壁面單元而言，網(wǎng)格在流向上的分辨率為50～150，壁面法向上的分辨率為0.77～0.97，展向上的分辨率為10～40。為了更好地在靠近壁面的邊界層區(qū)域生成密集的網(wǎng)格，使用了O形網(wǎng)格，圖1為計(jì)算域離散網(wǎng)格。為了滿足SSTγ-Reθt方法的要求，提高計(jì)算精度，捕捉流場(chǎng)的基本特征，對(duì)近壁區(qū)域進(jìn)行了局部網(wǎng)格細(xì)化。X、Y和Z方向的網(wǎng)格分辨率分別為98.6、1.03和38.4。

圖1 NACA0009鈍尾緣水翼計(jì)算域網(wǎng)格

本文對(duì)SSTγ-Reθt轉(zhuǎn)捩模型進(jìn)行了網(wǎng)格無(wú)關(guān)性研究。4組網(wǎng)格中每組水翼壁面網(wǎng)格的最大y+值都接近1。網(wǎng)格M4是由2 650 156節(jié)點(diǎn)構(gòu)成的粗網(wǎng)格。網(wǎng)格M3在水翼節(jié)點(diǎn)數(shù)上與M4相同，但在垂直于水翼壁面的方向以及翼型上下游區(qū)域被細(xì)化，因此節(jié)點(diǎn)數(shù)達(dá)到3 985 230。與網(wǎng)格M3相比，網(wǎng)格M2水翼表面包括翼型前緣和尾緣的節(jié)點(diǎn)數(shù)量都增加了1倍。網(wǎng)格M1除了沿水翼壁面網(wǎng)格節(jié)點(diǎn)數(shù)是網(wǎng)格M2的兩倍外，其他區(qū)域的節(jié)點(diǎn)數(shù)與網(wǎng)格M2相同。選取水翼尾渦脫落主頻作為水翼繞流中關(guān)鍵參數(shù)，對(duì)比4組網(wǎng)格數(shù)量下尾渦脫落主頻，如表1所示。尾渦脫落頻率比較表明，網(wǎng)格M1和網(wǎng)格M2之間的差異非常小，與試驗(yàn)值的誤差均在0.5%以內(nèi)，考慮到計(jì)算效率，因此網(wǎng)格M2將用于進(jìn)一步的計(jì)算。

表1 不同網(wǎng)格數(shù)量下翼型尾渦脫落主頻對(duì)比

為了精確獲得水翼繞流的邊界層分離和尾跡脫落頻率等流動(dòng)特性，需要足夠小的網(wǎng)格單元。本文采用ICEM CFD對(duì)計(jì)算區(qū)域劃分六面體網(wǎng)格?？刂品匠袒谟邢摅w積法在同位網(wǎng)格上進(jìn)行離散，對(duì)流項(xiàng)離散后的界面值采用混合有界迎風(fēng)格式進(jìn)行插值，最高具有二階精度；擴(kuò)散項(xiàng)離散后的界面值采用中心差分格式進(jìn)行插值；對(duì)于非定常雷諾平均模擬(URANS)，瞬態(tài)項(xiàng)離散采用具有二階精度的二階向后歐拉格式；采用耦合式解法聯(lián)立求解N-S離散方程組變量，采用多重網(wǎng)格技術(shù)以加快迭代收斂速度。時(shí)間步長(zhǎng)為5×10-5s，平均庫(kù)朗數(shù)小于1。該部分計(jì)算工作在中國(guó)農(nóng)業(yè)大學(xué)的高性能計(jì)算中心運(yùn)行，使用256 GB內(nèi)存和128個(gè)核進(jìn)行并行計(jì)算。為了獲得水翼尾跡區(qū)脫落渦的非定常特性，在尾跡區(qū)設(shè)置了監(jiān)測(cè)點(diǎn)。

1.3 控制方程

γ-Reθt轉(zhuǎn)捩模型由轉(zhuǎn)捩起始動(dòng)量雷諾數(shù)Reθt輸運(yùn)方程以及間歇因子γ輸運(yùn)方程構(gòu)成，其中，Reθt是判斷轉(zhuǎn)捩起始的條件，γ用于表征流動(dòng)狀態(tài)，公式為

(1)

(2)

式中t——時(shí)間ρ——密度

?Uj/?xj——j方向的速度偏導(dǎo)

Pγ——間歇因子生成項(xiàng)

Eγ——間歇因子破壞項(xiàng)

μt——渦粘系數(shù)μ——粘度

σf——常數(shù)，取1.0

Pθt——轉(zhuǎn)捩源項(xiàng)

σθt——常數(shù)，取2.0

當(dāng)轉(zhuǎn)捩發(fā)生在層流分離點(diǎn)附近時(shí)，由于分離剪切層的湍動(dòng)能較小，預(yù)測(cè)湍動(dòng)能發(fā)展至湍流再附著的位置與試驗(yàn)獲得的實(shí)際湍流再附著位置相差較大。為了解決這一問(wèn)題，γ-Reθt轉(zhuǎn)捩模型對(duì)間歇因子γ進(jìn)行了修正，修正的主要思路是，一旦層流邊界層分離時(shí)，允許局部間歇因子γ超過(guò)1。這將導(dǎo)致k的增加，因此將會(huì)使再附著點(diǎn)的位置更早發(fā)生[18]，具體修正公式為

(3)

γeff=max(γ,γsep)

(4)

式中γsep——分離流動(dòng)計(jì)算時(shí)采用的間歇因子

γeff——修正后有效間歇因子

Rev——渦量雷諾數(shù)

s1——常數(shù)，取8.0

Reθc——轉(zhuǎn)捩臨界動(dòng)量厚度雷諾數(shù)

Freattach、Fθt——內(nèi)部關(guān)聯(lián)函數(shù)

2 計(jì)算結(jié)果分析

圖2為SSTγ-Reθt轉(zhuǎn)捩模型與SSTk-ω湍流模型計(jì)算得出的水翼邊界層切向時(shí)均速度分布的結(jié)果。圖中，y為垂直于水翼壁面的距離，h為翼型尾部厚度，Uxave為平行于翼型表面的速度分量，Ute為邊界層外緣時(shí)均速度。圖中表示了相對(duì)弦長(zhǎng)x/L(x方向相對(duì)水翼弦長(zhǎng)L的距離)分別為0.1、0.3、0.6、0.8、0.9以及0.99處邊界層速度分布情況。為了分析不同位置的速度分布，位置x/L為0.3、0.6、0.8、0.9和0.99處的速度分布曲線沿橫坐標(biāo)分別移動(dòng)了1、2、3、4、5個(gè)單位坐標(biāo)。可以看出，在水翼前緣，即0.1倍翼型弦長(zhǎng)處，兩個(gè)模型預(yù)測(cè)的邊界層速度分布與試驗(yàn)值較為吻合，但是沿著流動(dòng)方向，SSTk-ω湍流模型計(jì)算的邊界層速度分布誤差越來(lái)越大；SSTγ-Reθt轉(zhuǎn)捩模型只是在靠近水翼尾部區(qū)域，即0.8倍翼型弦長(zhǎng)處計(jì)算的誤差較大。

圖2 不同湍流模型下水翼邊界層切向時(shí)均速度分布

本文邊界層厚度為水翼邊界層內(nèi)從壁面(速度為零)開始，沿法線方向當(dāng)x方向的速度Ux至外部速度Ue的99%位置之間的距離，記δ，即δ=y|Ux=0.99Ue。不同模型計(jì)算的無(wú)量綱邊界層相對(duì)厚度δ/h沿水翼翼弦長(zhǎng)的分布如圖3所示。隨著邊界層的發(fā)展，邊界層厚度隨著與水翼前緣距離的增加而增加。兩種模型計(jì)算得到的邊界層厚度的差異較大，SSTk-ω模型計(jì)算的厚度遠(yuǎn)大于試驗(yàn)結(jié)果，SSTγ-Reθt轉(zhuǎn)捩模型計(jì)算的水翼前緣的邊界層厚度與試驗(yàn)結(jié)果吻合較好。

圖3 不同湍流模型下水翼邊界層相對(duì)厚度δ/h分布

圖4為兩種模型下水翼邊界層形狀因子H12分布，形狀因子H12定義為邊界層位移厚度δ1與其動(dòng)量厚度δ2之比。一般認(rèn)為，當(dāng)H12超過(guò)2.6時(shí)，邊界層為層流；當(dāng)H12小于1.5時(shí)，邊界層為湍流；對(duì)于H12在1.5～2.6范圍時(shí)，邊界層處于轉(zhuǎn)捩狀態(tài)?？梢钥闯?，SSTγ-Reθt轉(zhuǎn)捩模型預(yù)測(cè)的形狀因子在x=0.2L后偏離試驗(yàn)值，并在0.65L附近下降至1.5，這說(shuō)明該位置處已完成轉(zhuǎn)捩，比試驗(yàn)中觀察到的0.85L轉(zhuǎn)捩結(jié)束較為提前。

圖4 不同湍流模型下水翼邊界層形狀因子H12分布

進(jìn)一步分析SSTγ-Reθt轉(zhuǎn)捩模型對(duì)不同雷諾數(shù)下邊界層轉(zhuǎn)捩的預(yù)測(cè)效果，圖5給出了邊界層完全轉(zhuǎn)捩位置xt隨雷諾數(shù)ReL的變化。可以看出，在雷諾數(shù)低于1.6×106時(shí)，邊界層轉(zhuǎn)捩位置與實(shí)測(cè)值比較一致，而在雷諾數(shù)高于1.6×106時(shí)，隨著雷諾數(shù)的增大，邊界層轉(zhuǎn)捩位置逐漸偏離實(shí)測(cè)值，說(shuō)明SSTγ-Reθt轉(zhuǎn)捩模型對(duì)高雷諾數(shù)水翼邊界層轉(zhuǎn)捩預(yù)測(cè)效果不佳。

圖5 不同雷諾數(shù)下邊界層完全轉(zhuǎn)捩位置

為進(jìn)一步探究SSTγ-Reθt轉(zhuǎn)捩模型在高雷諾數(shù)下預(yù)測(cè)效果不佳的原因，給出圖6所示數(shù)值計(jì)算得到的雷諾數(shù)為2.0×106時(shí)，湍流度Tu沿流向的發(fā)展，其中負(fù)值表示水翼頭部的上游位置。可以看出，沿著流動(dòng)方向，湍流度Tu快速衰減，當(dāng)?shù)竭_(dá)水翼前緣時(shí)，湍流度未達(dá)到1%，這與試驗(yàn)所測(cè)不相符。因此，高雷諾數(shù)下計(jì)算域進(jìn)口必需設(shè)置較大的Tu，但過(guò)大的Tu會(huì)加速間歇因子γ預(yù)測(cè)值的增加，導(dǎo)致邊界層轉(zhuǎn)捩過(guò)程發(fā)展過(guò)快，文獻(xiàn)[18]也指出，過(guò)大的Tu會(huì)導(dǎo)致預(yù)測(cè)的壁面摩阻系數(shù)偏離層流值。因此高雷諾數(shù)下自由流場(chǎng)中湍流度衰減過(guò)快是導(dǎo)致轉(zhuǎn)捩位置預(yù)測(cè)不準(zhǔn)確的主要原因。

圖6 湍流度沿流向發(fā)展

3 基于環(huán)境源項(xiàng)法的修正

為解決高雷諾數(shù)條件下，SSTγ-Reθt轉(zhuǎn)捩模型(后續(xù)稱為原轉(zhuǎn)捩模型)預(yù)測(cè)邊界層轉(zhuǎn)捩提前的問(wèn)題，對(duì)原轉(zhuǎn)捩模型進(jìn)行基于流動(dòng)特征的修正，以提高在高雷諾數(shù)下邊界層轉(zhuǎn)捩預(yù)測(cè)的精度。

在原轉(zhuǎn)捩模型中，Reθt方程的計(jì)算與自由流的湍流度有關(guān)。當(dāng)?shù)刈杂闪魍牧鞫阮A(yù)測(cè)的精度直接影響了轉(zhuǎn)捩位置的預(yù)測(cè)精度。為了解決自由流動(dòng)中湍流度衰減預(yù)測(cè)誤差大而導(dǎo)致壁面附近湍流度預(yù)測(cè)精度較低的問(wèn)題，文獻(xiàn)[7]提出了環(huán)境源項(xiàng)法來(lái)控制自由流動(dòng)中湍流度的衰減，即在指定區(qū)域保持環(huán)境湍流度和湍流耗散率不變，在其余區(qū)域按給定的衰減率自由衰減。文獻(xiàn)[19]根據(jù)這一思想，對(duì)SSTk-ω湍流模型輸運(yùn)方程源項(xiàng)進(jìn)行修改，公式為

(5)

(6)

式中ka、ωa——指定區(qū)域的環(huán)境湍動(dòng)能、環(huán)境湍流比耗散率

τ——應(yīng)力

β*、β3——湍流模型系數(shù)

γ3、σk3、σω3——湍流模型常數(shù)

νt——運(yùn)動(dòng)渦粘系數(shù)

F1——湍流模型的函數(shù)

ui、uj——i、j方向的速度分量

ω——湍流比耗散率

σω2——常數(shù)，取0.856

根據(jù)源項(xiàng)均衡態(tài)假設(shè)[20]有

(7)

為此，根據(jù)NACA0009鈍尾緣水翼試驗(yàn)，選取試驗(yàn)中來(lái)流湍流度1%為參考對(duì)ωa取值進(jìn)行校準(zhǔn)，并定義渦粘系數(shù)比RT=ka/(νωa)對(duì)ωa進(jìn)行了無(wú)量綱化，ν表示運(yùn)動(dòng)粘度系數(shù)，具體取值如表2所示。改變?chǔ)豠，仍保證收斂精度高于10-4，從而保證計(jì)算精度。

表2 湍流度1%下環(huán)境湍流比耗散率標(biāo)定設(shè)置

選取形狀因子這一參數(shù)，對(duì)不同ωa取值下的結(jié)果進(jìn)行分析，由圖7可以看出，與試驗(yàn)值相比，改變?chǔ)豠取值，在0.2L處H12取值變化不大，而后同一位置處H12的取值隨著ωa的增加而減?。豢梢钥闯霎?dāng)ωa為4 000，即渦粘系數(shù)比為15時(shí)，H12整體分布與試驗(yàn)較為一致。因此，對(duì)于同一湍流度下不同雷諾數(shù)流動(dòng)，只需根據(jù)RT為常數(shù)即可確定ωa。

圖7 1%來(lái)流湍流度下不同環(huán)境湍流比耗散率預(yù)測(cè)的邊界層形狀因子分布

間歇因子γ是表征流動(dòng)的一個(gè)重要參數(shù)。在邊界層層流區(qū)γ=0，而在邊界層外γ=1。在轉(zhuǎn)捩區(qū)域，隨著湍流強(qiáng)度的變化，γ介于0和1之間。圖8為采用不同渦粘系數(shù)比RT計(jì)算得出的水翼邊界層間歇因子γ的分布。當(dāng)ωa大于12 000 s-1，即渦粘系數(shù)比RT小于5時(shí)，邊界層間歇因子γ在0.4L位置處開始后迅速增加，轉(zhuǎn)捩開始發(fā)生，隨后在0.6L發(fā)展為全湍流邊界層。當(dāng)渦粘系數(shù)比為15時(shí)，間歇因子在0.8L處增加至1，即在該位置處完成轉(zhuǎn)捩，與試驗(yàn)值相符。繼續(xù)增加渦粘系數(shù)比，預(yù)測(cè)的轉(zhuǎn)捩位置逐漸向后移動(dòng)。當(dāng)渦粘系數(shù)比大于20時(shí)，轉(zhuǎn)捩結(jié)束位置的預(yù)測(cè)值穩(wěn)定在0.9L，預(yù)測(cè)的轉(zhuǎn)捩區(qū)長(zhǎng)度大于試驗(yàn)觀測(cè)值。

圖8 不同渦粘系數(shù)比下間歇因子γ分布

在弦長(zhǎng)雷諾數(shù)ReL=2×106條件下，驗(yàn)證了基于環(huán)境源項(xiàng)法修正的轉(zhuǎn)捩模型能夠有效預(yù)測(cè)轉(zhuǎn)捩位置。進(jìn)一步針對(duì)不同弦長(zhǎng)雷諾數(shù)，采用上述方法，獲得不同湍流比耗散率以及渦粘系數(shù)比，如表3所示?？梢钥闯?，隨著雷諾數(shù)的增加，對(duì)應(yīng)的渦粘系數(shù)比也隨之增加。

表3 不同雷諾數(shù)下指定區(qū)域環(huán)境湍動(dòng)能ka和湍流比耗散率ωa設(shè)置

圖9為不同雷諾數(shù)下邊界層轉(zhuǎn)捩位置，通過(guò)與試驗(yàn)結(jié)果的對(duì)比，雖然在高雷諾數(shù)區(qū)間預(yù)測(cè)的轉(zhuǎn)捩位置較實(shí)測(cè)值提前，而在低雷諾數(shù)區(qū)間預(yù)測(cè)的轉(zhuǎn)捩位置則稍有延后，但整體來(lái)看，修正模型在很大程度上提高了對(duì)轉(zhuǎn)捩位置預(yù)測(cè)的準(zhǔn)確性。

圖9 不同雷諾數(shù)下邊界層轉(zhuǎn)捩位置

因此，在運(yùn)用該修正轉(zhuǎn)捩模型時(shí)，首先確定來(lái)流湍流強(qiáng)度為1%。然后根據(jù)來(lái)流速度或者弦長(zhǎng)雷諾數(shù)計(jì)算ka，再根據(jù)前文給定的推薦值選取ωa用于后續(xù)計(jì)算。以ReL=2.6×106為例，計(jì)算得到的ωa為5 200 s-1，渦粘系數(shù)比為19.5，因此，選擇該值進(jìn)行轉(zhuǎn)捩數(shù)值計(jì)算，圖10為邊界層形狀因子分布，可以看到修正模型預(yù)測(cè)的從0.7L處迅速減小，邊界層迅速?gòu)膶恿鳡顟B(tài)過(guò)渡到湍流狀態(tài)，這與試驗(yàn)測(cè)量結(jié)果是比較吻合的。

圖10 邊界層形狀因子分布(ReL=2.6×106)

4 修正轉(zhuǎn)捩模型算例驗(yàn)證

4.1 Donaldson修型尾緣水翼

Donaldson修型尾緣水翼在鈍尾緣水翼基礎(chǔ)上，采用一條45°直線和一條三次多項(xiàng)式曲線對(duì)尾緣形狀進(jìn)行了修改[17]，由于試驗(yàn)均在同一裝置下進(jìn)行，本算例的計(jì)算域及邊界條件與NACA0009鈍尾緣水翼算例基本一致，計(jì)算域網(wǎng)格如圖11所示。

圖11 Donaldson修型尾緣水翼計(jì)算域網(wǎng)格

文獻(xiàn)[17]對(duì)該水翼在不同雷諾數(shù)下的流動(dòng)進(jìn)行了試驗(yàn)研究，試驗(yàn)測(cè)量結(jié)果可為本算例驗(yàn)證提供參考。將基于環(huán)境源項(xiàng)法修正模型應(yīng)用于Donaldson修型尾緣水翼繞流計(jì)算中，不同雷諾數(shù)下環(huán)境湍動(dòng)能和湍流比耗散率的設(shè)置如表4所示。

表4 不同雷諾數(shù)下環(huán)境湍動(dòng)能和湍流比耗散率設(shè)置

圖12為翼型弦長(zhǎng)雷諾數(shù)為2×106時(shí)吸力面邊界層x方向時(shí)均速度分布，可以看出，修正模型預(yù)測(cè)的邊界層厚度略高于試驗(yàn)值，但相對(duì)原始模型已有一定的提高。

圖12 吸力面邊界層時(shí)均速度分布(ReL=2.0×106)

圖13給出了不同雷諾數(shù)下尾跡區(qū)尾渦脫落頻率，原轉(zhuǎn)捩模型和修正轉(zhuǎn)捩模型預(yù)測(cè)結(jié)果均與試驗(yàn)結(jié)果有較大偏差，其原因可能與尾緣修型對(duì)尾渦脫落的影響機(jī)理有關(guān)，試驗(yàn)研究發(fā)現(xiàn)，尾緣修型后，原本交替排列的脫落渦會(huì)產(chǎn)生錯(cuò)位，導(dǎo)致尾緣上下脫落渦產(chǎn)生相互作用，影響尾渦脫落頻率。但整體來(lái)看修正模型提高了對(duì)尾渦脫落頻率的預(yù)測(cè)精度，相比原轉(zhuǎn)捩模型最大提高約8%。

圖13 不同雷諾數(shù)下尾跡尾渦脫落頻率

4.2 NACA0016水翼

文獻(xiàn)[21]在美國(guó)海軍大型空化水洞中進(jìn)行了NACA0016的水翼試驗(yàn)。水翼弦長(zhǎng)L為2.133 6 m，頂面呈拱形，底面平行于試驗(yàn)段兩側(cè)，試驗(yàn)段寬度為1.43L。本算例采用二維計(jì)算，為了提高近壁網(wǎng)格的正交性，在水翼周圍采用O形網(wǎng)格塊來(lái)細(xì)化網(wǎng)格，壁面網(wǎng)格法向擴(kuò)展比為1.1，得到水翼表面的最大y+約為1，最終劃分的網(wǎng)格總數(shù)約1.7×105，計(jì)算域網(wǎng)格如圖14所示。

圖14 NACA0016水翼計(jì)算域網(wǎng)格

將基于環(huán)境源項(xiàng)法的修正模型應(yīng)用于NACA0016水翼繞流流動(dòng)，對(duì)弦長(zhǎng)雷諾數(shù)ReL=1.7×107的流動(dòng)進(jìn)行計(jì)算，計(jì)算得出ka=0.001 35、ωa=2 249.7 s-1。圖15為距水翼前緣0.93L處，原轉(zhuǎn)捩模型、修正轉(zhuǎn)捩模型計(jì)算邊界層速度值與試驗(yàn)值對(duì)比圖，修正轉(zhuǎn)捩模型由于準(zhǔn)確預(yù)測(cè)了此處的層流邊界層，計(jì)算的邊界層速度分布以及邊界層厚度均與試驗(yàn)結(jié)果比較一致。

圖15 x=0.93L處邊界層速度分布

圖16為數(shù)值計(jì)算得到的雷諾數(shù)為1.7×107時(shí)，湍流度Tu沿流向的發(fā)展?？梢钥闯?，沿著流動(dòng)方向，采用修正轉(zhuǎn)捩模型，湍流度Tu快速衰減，當(dāng)?shù)竭_(dá)水翼前緣時(shí)，湍流度約為1%，與試驗(yàn)所測(cè)相符，修正后的轉(zhuǎn)捩模型較好地保證了水翼前緣的湍流度，進(jìn)而較為準(zhǔn)確地預(yù)測(cè)出轉(zhuǎn)捩區(qū)邊界層參數(shù)。

圖16 模型修正前后計(jì)算湍流度沿流向發(fā)展

圖17為NACA0016水翼吸力面邊界層相對(duì)厚度δ/L以及摩阻系數(shù)Cf的試驗(yàn)值和模擬值分布圖。可以看出，在靠近水翼中部附近，原轉(zhuǎn)捩模型與修正轉(zhuǎn)捩模型計(jì)算的邊界層相對(duì)厚度分別為0.001 53和0.000 91，相對(duì)誤差分別為112.3%和27.57%；原轉(zhuǎn)捩模型與修正轉(zhuǎn)捩模型計(jì)算的摩阻系數(shù)分別為0.005 22和0.000 71，相對(duì)誤差分別為1 100%和74.0%。在水翼尾緣附近，原轉(zhuǎn)捩模型和修正轉(zhuǎn)捩模型預(yù)測(cè)的相對(duì)厚度分別為0.012 03和0.010 80，相對(duì)誤差分別為16.4%和4.85%；原轉(zhuǎn)捩模型與修正轉(zhuǎn)捩模型計(jì)算的摩阻系數(shù)分別為0.002 6和0.002 5，相對(duì)誤差分別為30%和25%?？傮w來(lái)看，修正轉(zhuǎn)捩模型在邊界層相對(duì)厚度、摩阻系數(shù)的預(yù)測(cè)精度上比原轉(zhuǎn)捩模型均提高3倍以上。

圖17 水翼邊界層相對(duì)厚度和壁面摩阻系數(shù)Cf分布

5 結(jié)論

(1)相比于SSTk-ω模型，采用SSTγ-Reθt轉(zhuǎn)捩預(yù)測(cè)模型能夠顯著提高邊界層流場(chǎng)分布計(jì)算精度。即相同網(wǎng)格條件下，在近壁區(qū)邊界層流場(chǎng)的預(yù)測(cè)上，SSTγ-Reθt轉(zhuǎn)捩模型的預(yù)測(cè)精度高于SSTk-ω湍流模型，尤其在低雷諾數(shù)條件下，SSTγ-Reθt轉(zhuǎn)捩模型的預(yù)測(cè)精度與試驗(yàn)值更為接近。

(2)對(duì)于高雷諾數(shù)邊界層轉(zhuǎn)捩流動(dòng)，壁面附近湍流度預(yù)測(cè)精度不夠，導(dǎo)致SSTγ-Reθt轉(zhuǎn)捩模型預(yù)測(cè)精度有所下降，通過(guò)引入環(huán)境湍動(dòng)能和環(huán)境湍流比耗散率參數(shù)，建立湍流比耗散率與雷諾數(shù)的關(guān)系，修正轉(zhuǎn)捩模型的輸運(yùn)方程環(huán)境源項(xiàng)，能夠有效提高其對(duì)高雷諾數(shù)邊界層轉(zhuǎn)捩流動(dòng)的預(yù)測(cè)精度。

(3)基于SSTγ-Reθt轉(zhuǎn)捩修正模型，針對(duì)Donaldson修型尾緣水翼高雷諾數(shù)條件下的尾渦脫落頻率等典型流動(dòng)特征預(yù)測(cè)精度相比原轉(zhuǎn)捩模型最大提高約8%；針對(duì)NACA0016水翼中部轉(zhuǎn)捩區(qū)域的邊界層相對(duì)厚度、摩阻系數(shù)的預(yù)測(cè)精度相比原轉(zhuǎn)捩模型均提高3倍以上，水翼尾緣局部區(qū)域的邊界層相對(duì)厚度、摩阻系數(shù)預(yù)測(cè)精度也有顯著提升。