基于環境源項法的SST γ-Reθt轉捩模型修正

葉長亮 王福軍 湯 遠 陳 俊 鄭 源

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

0 引言

軸流泵具有大流量、低揚程的特點，在農業灌溉領域應用廣泛[1-2]。水翼作為軸流泵葉片的簡化模型，其繞流特性對于研究軸流泵葉片的繞流分布具有參考意義。文獻[3-5]通過研究水翼分離流動特性為軸流泵內部流動特性提供了借鑒，然而，目前對軸流泵葉片附近的轉捩流動仍缺乏認識。

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

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

1 數值計算方法

1.1 模型建立與網格劃分

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

1.2 網格無關性驗證

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

圖1 NACA0009鈍尾緣水翼計算域網格

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

表1 不同網格數量下翼型尾渦脫落主頻對比

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

1.3 控制方程

γ-Reθt轉捩模型由轉捩起始動量雷諾數Reθt輸運方程以及間歇因子γ輸運方程構成，其中，Reθt是判斷轉捩起始的條件，γ用于表征流動狀態，公式為

(1)

(2)

式中t——時間ρ——密度

?Uj/?xj——j方向的速度偏導

Pγ——間歇因子生成項

Eγ——間歇因子破壞項

μt——渦粘系數μ——粘度

σf——常數，取1.0

Pθt——轉捩源項

σθt——常數，取2.0

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

(3)

γeff=max(γ,γsep)

(4)

式中γsep——分離流動計算時采用的間歇因子

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

Rev——渦量雷諾數

s1——常數，取8.0

Reθc——轉捩臨界動量厚度雷諾數

Freattach、Fθt——內部關聯函數

2 計算結果分析

圖2為SSTγ-Reθt轉捩模型與SSTk-ω湍流模型計算得出的水翼邊界層切向時均速度分布的結果。圖中，y為垂直于水翼壁面的距離，h為翼型尾部厚度，Uxave為平行于翼型表面的速度分量，Ute為邊界層外緣時均速度。圖中表示了相對弦長x/L(x方向相對水翼弦長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處的速度分布曲線沿橫坐標分別移動了1、2、3、4、5個單位坐標?？梢钥闯?，在水翼前緣，即0.1倍翼型弦長處，兩個模型預測的邊界層速度分布與試驗值較為吻合，但是沿著流動方向，SSTk-ω湍流模型計算的邊界層速度分布誤差越來越大；SSTγ-Reθt轉捩模型只是在靠近水翼尾部區域，即0.8倍翼型弦長處計算的誤差較大。

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

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

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

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

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

進一步分析SSTγ-Reθt轉捩模型對不同雷諾數下邊界層轉捩的預測效果，圖5給出了邊界層完全轉捩位置xt隨雷諾數ReL的變化?？梢钥闯?，在雷諾數低于1.6×106時，邊界層轉捩位置與實測值比較一致，而在雷諾數高于1.6×106時，隨著雷諾數的增大，邊界層轉捩位置逐漸偏離實測值，說明SSTγ-Reθt轉捩模型對高雷諾數水翼邊界層轉捩預測效果不佳。

圖5 不同雷諾數下邊界層完全轉捩位置

為進一步探究SSTγ-Reθt轉捩模型在高雷諾數下預測效果不佳的原因，給出圖6所示數值計算得到的雷諾數為2.0×106時，湍流度Tu沿流向的發展，其中負值表示水翼頭部的上游位置?？梢钥闯?，沿著流動方向，湍流度Tu快速衰減，當到達水翼前緣時，湍流度未達到1%，這與試驗所測不相符。因此，高雷諾數下計算域進口必需設置較大的Tu，但過大的Tu會加速間歇因子γ預測值的增加，導致邊界層轉捩過程發展過快，文獻[18]也指出，過大的Tu會導致預測的壁面摩阻系數偏離層流值。因此高雷諾數下自由流場中湍流度衰減過快是導致轉捩位置預測不準確的主要原因。

圖6 湍流度沿流向發展

3 基于環境源項法的修正

為解決高雷諾數條件下，SSTγ-Reθt轉捩模型(后續稱為原轉捩模型)預測邊界層轉捩提前的問題，對原轉捩模型進行基于流動特征的修正，以提高在高雷諾數下邊界層轉捩預測的精度。

在原轉捩模型中，Reθt方程的計算與自由流的湍流度有關。當地自由流湍流度預測的精度直接影響了轉捩位置的預測精度。為了解決自由流動中湍流度衰減預測誤差大而導致壁面附近湍流度預測精度較低的問題，文獻[7]提出了環境源項法來控制自由流動中湍流度的衰減，即在指定區域保持環境湍流度和湍流耗散率不變，在其余區域按給定的衰減率自由衰減。文獻[19]根據這一思想，對SSTk-ω湍流模型輸運方程源項進行修改，公式為

(5)

(6)

式中ka、ωa——指定區域的環境湍動能、環境湍流比耗散率

τ——應力

β*、β3——湍流模型系數

γ3、σk3、σω3——湍流模型常數

νt——運動渦粘系數

F1——湍流模型的函數

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

ω——湍流比耗散率

σω2——常數，取0.856

根據源項均衡態假設[20]有

(7)

為此，根據NACA0009鈍尾緣水翼試驗，選取試驗中來流湍流度1%為參考對ωa取值進行校準，并定義渦粘系數比RT=ka/(νωa)對ωa進行了無量綱化，ν表示運動粘度系數，具體取值如表2所示。改變ωa，仍保證收斂精度高于10-4，從而保證計算精度。

表2 湍流度1%下環境湍流比耗散率標定設置

選取形狀因子這一參數，對不同ωa取值下的結果進行分析，由圖7可以看出，與試驗值相比，改變ωa取值，在0.2L處H12取值變化不大，而后同一位置處H12的取值隨著ωa的增加而減??；可以看出當ωa為4 000，即渦粘系數比為15時，H12整體分布與試驗較為一致。因此，對于同一湍流度下不同雷諾數流動，只需根據RT為常數即可確定ωa。

圖7 1%來流湍流度下不同環境湍流比耗散率預測的邊界層形狀因子分布

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

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

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

表3 不同雷諾數下指定區域環境湍動能ka和湍流比耗散率ωa設置

圖9為不同雷諾數下邊界層轉捩位置，通過與試驗結果的對比，雖然在高雷諾數區間預測的轉捩位置較實測值提前，而在低雷諾數區間預測的轉捩位置則稍有延后，但整體來看，修正模型在很大程度上提高了對轉捩位置預測的準確性。

圖9 不同雷諾數下邊界層轉捩位置

因此，在運用該修正轉捩模型時，首先確定來流湍流強度為1%。然后根據來流速度或者弦長雷諾數計算ka，再根據前文給定的推薦值選取ωa用于后續計算。以ReL=2.6×106為例，計算得到的ωa為5 200 s-1，渦粘系數比為19.5，因此，選擇該值進行轉捩數值計算，圖10為邊界層形狀因子分布，可以看到修正模型預測的從0.7L處迅速減小，邊界層迅速從層流狀態過渡到湍流狀態，這與試驗測量結果是比較吻合的。

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

4 修正轉捩模型算例驗證

4.1 Donaldson修型尾緣水翼

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

圖11 Donaldson修型尾緣水翼計算域網格

文獻[17]對該水翼在不同雷諾數下的流動進行了試驗研究，試驗測量結果可為本算例驗證提供參考。將基于環境源項法修正模型應用于Donaldson修型尾緣水翼繞流計算中，不同雷諾數下環境湍動能和湍流比耗散率的設置如表4所示。

表4 不同雷諾數下環境湍動能和湍流比耗散率設置

圖12為翼型弦長雷諾數為2×106時吸力面邊界層x方向時均速度分布，可以看出，修正模型預測的邊界層厚度略高于試驗值，但相對原始模型已有一定的提高。

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

圖13給出了不同雷諾數下尾跡區尾渦脫落頻率，原轉捩模型和修正轉捩模型預測結果均與試驗結果有較大偏差，其原因可能與尾緣修型對尾渦脫落的影響機理有關，試驗研究發現，尾緣修型后，原本交替排列的脫落渦會產生錯位，導致尾緣上下脫落渦產生相互作用，影響尾渦脫落頻率。但整體來看修正模型提高了對尾渦脫落頻率的預測精度，相比原轉捩模型最大提高約8%。

圖13 不同雷諾數下尾跡尾渦脫落頻率

4.2 NACA0016水翼

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

圖14 NACA0016水翼計算域網格

將基于環境源項法的修正模型應用于NACA0016水翼繞流流動，對弦長雷諾數ReL=1.7×107的流動進行計算，計算得出ka=0.001 35、ωa=2 249.7 s-1。圖15為距水翼前緣0.93L處，原轉捩模型、修正轉捩模型計算邊界層速度值與試驗值對比圖，修正轉捩模型由于準確預測了此處的層流邊界層，計算的邊界層速度分布以及邊界層厚度均與試驗結果比較一致。

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

圖16為數值計算得到的雷諾數為1.7×107時，湍流度Tu沿流向的發展?？梢钥闯?，沿著流動方向，采用修正轉捩模型，湍流度Tu快速衰減，當到達水翼前緣時，湍流度約為1%，與試驗所測相符，修正后的轉捩模型較好地保證了水翼前緣的湍流度，進而較為準確地預測出轉捩區邊界層參數。

圖16 模型修正前后計算湍流度沿流向發展

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

圖17 水翼邊界層相對厚度和壁面摩阻系數Cf分布

5 結論

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

(2)對于高雷諾數邊界層轉捩流動，壁面附近湍流度預測精度不夠，導致SSTγ-Reθt轉捩模型預測精度有所下降，通過引入環境湍動能和環境湍流比耗散率參數，建立湍流比耗散率與雷諾數的關系，修正轉捩模型的輸運方程環境源項，能夠有效提高其對高雷諾數邊界層轉捩流動的預測精度。

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