亞格子模型對三維扭曲水翼空化現象的影響

何朋朋，李子如，張孝旺，賀偉*

1 武漢理工大學 高性能艦船技術教育部重點實驗室，湖北 武漢 430063

2 武漢理工大學 船海與能源動力工程學院，湖北 武漢 430063

0 引 言

在船舶領域中，螺旋槳空化不僅會惡化推進性能，而且還會伴隨產生強烈的振動、噪聲及剝蝕等問題。由于螺旋槳為三維扭曲形式，在不同半徑位置處的空化類型各異，其中針對非定常云狀空泡初生與潰滅機制的探討，對于螺旋槳空蝕及空化噪聲的合理預測至關重要。為了更好地觀察及研究非定常云狀空化的演變過程，許多學者對展向攻角漸變的三維扭曲水翼進行了數值與試驗研究。

例如，文獻[1–3]分別通過LDV（laser doppler velocimetry）技術、高速攝影技術和PIV（particle image velocimetry）技術在空化水洞中對三維水翼的空泡形態、攻角及空化數、空泡初生及分離現象開展了試驗研究，從結果觀察到了馬蹄形云空泡，發現了空化流動特性主要與水翼攻角及空化數有關。文獻[4]以Delft Twist11N 三維扭曲水翼為對象，引入了側射流（side-entrant jet）的概念，基于試驗結果分析了側射流及主回射流（re-entrant jet）在空泡脫落中的作用。文獻[5–8] 分別結合經過湍流黏度修正的SSTk-ω 和RNGk-ε 等湍流模型和Schnerr-Sauer（S-S），Singhal 及Zwart 空化模型，對Twist11N 水翼的空化現象進行了數值模擬，并使用經過黏度修正的雷諾平均N-S （RANS）方法模擬了大尺度空泡脫落現象，但空泡脫落頻率與試驗結果的誤差較大，且云空泡形態明顯小于試驗結果，其主要原因是RANS 方法對所有尺度的湍動特征參數均進行了雷諾平均，對空泡精細結構和非定常瞬態特性的描述略顯不足。文獻[9]針對三維水翼梢渦空泡的初生問題，分別使用SSTk-ω 湍流模型、分離渦（detached-eddy simulation，DES）和大渦模擬（large-eddy simulation, LES）模型對水翼的梢渦流場進行了計算研究，結果表明，基于LES 方法的計算結果與試驗值吻合較好。文獻[10] 基于解析精度較高LES 方法，對NACA 66 水翼的非定常空化湍流進行數值研究，結果證實了壓力脈動與空泡脫落過程有關，且空泡體積加速度是造成空泡水翼周圍壓力脈動的主要原因。文獻[11] 使用兩套不同密度的網格研究了網格分辨率對三維扭曲水翼周圍云空泡流動數值模擬的影響，結果表明，對于空化流動數值模擬，細網格數值模擬結果的精度比粗網格的有顯著提高。文獻[12] 從拉格朗日觀點對三維扭曲水翼的空化流動進行分析，結果更直觀顯示出回射流和側射流的軌跡線。文獻[13]對NACA 66水翼的通氣空化流動進行數值研究，通過與試驗結果對比驗證了數值方法的有效性。總之，雖然基于LES 方法的三維扭曲水翼空化數值模擬已展開了較多的研究，但目前對于網格劃分及亞格子模型的選擇大多還依賴研究人員的經驗，缺乏客觀依據。

綜上，針對使用LES 方法的三維扭曲水翼空化數值模擬中網格密度及亞格子模型的適應性問題，本文將基于LES 的黏性多相流數值方法對圍繞Twist11N 扭翼型水翼的空化現象進行研究，重點探討不同網格密度及亞格子模型對此類水翼空化流特性的預報能力和適應性，分析該水翼云狀空化現象的非定常瞬態特性，期望研究結果可以在基于LES 方法的三維水翼及螺旋槳空化預報中為網格劃分及網格子模型的選取提供參考，并有助于剖析船舶螺旋槳中的非定常空化流場。

1 數學模型

1.1 控制方程

本文將包含汽液兩相的空化流視為密度可變的單一流體，運用基于均質平衡流理論的均質混合模型來描述汽液兩相。汽液混合空化流場的連續性方程和動量方程如下：

式中，下標m，v，l 分別指汽液混合相、汽相、液相。

1.2 大渦模擬方法

空化流動是一種復雜的湍流，本文使用具有較高解析精度的LES 方法對湍流進行處理。LES方法利用濾波器對Navier-Stokes（N-S）方程進行濾波，直接求解濾波后的N-S 方程，獲得可解尺度物理量，而濾出的亞格子尺度物理量則通過相應的亞格子模型進行求解。對于不可壓流動，經濾波處理的連續性方程及動量方程如下：

1.3 空化模型

采用基于質量輸運方程的S-S 空化模型來描述汽液兩相間的轉換，該模型是以簡化的Rayleigh-Plesset 氣泡動力學方程為基礎而建立的。

式中，nb為不氣泡數密度，一般取值為1013。

2 計算模型與設置

本文以Twist11N 扭曲型水翼為研究對象，其展向不同位置處的截面形狀均為NACA 0009 翼型，但經過了變攻角設計，從兩端向中部的展向截面攻角變化范圍為–2°～9°，具體詳見文獻[4]。圖1 所示為Twist 11N 水翼不同角度的視圖。其中：坐標原點位于水翼端部截面的弦長中點，XC表示入流方向，YC表示翼展方向，ZC表示垂直向上；LE 與TE 分別表示水翼的導邊和隨邊；弦長C=0.15 m，展長l=0.30 m。

圖1 Twist11N 水翼幾何圖Fig. 1 Geometrical views of Twist11N hydrofoil

考慮到Twist11N 水翼幾何模型具有對稱性的特點，計算時僅取水翼的一半。圖2 所示為計算域示意圖。其中，側壁處翼截面弦長的中點與上下壁面的距離分別為C，翼展方向為C，入口至導邊為2C，出口至導邊為5C。

圖2 一半Twist11N 水翼的計算域示意圖Fig. 2 Schematics of computational domain for half-span ofTwist11N hydrofoil

圖3 Twist11N 水翼周圍O 型分塊網格的劃分形式Fig. 3 Grid scheme of the O-shaped block around Twist11N hydrofoil

表1 Twist11N 水翼網格劃分參數Table 1 Grid parameters of Twist11N hydrofoil

數值研究使用ANSYS FLUENT 19.0 軟件。Twist11N 水翼在空化數 σ=1.07 時，非定常空泡流數值計算的邊界條件設置如下：入口采用速度入口條件，入口速度V∞=6.97 m/s；出口條件為壓力出口，壓力值由空化數公式σ=(pout?pv)/0.5ρlV∞2計算得到，出口壓力pout=29.0 kPa；計算域頂部、底部、側壁設置為自由滑移壁面，水翼表面為無滑移壁面，翼展中剖面為對稱面邊界條件。汽液兩相的密度、動力黏度系數及汽化壓力等參數參照試驗溫度T=24℃取值。參照文獻[17]中關于時間步長選取方法（?t=Tref/200,Tref=C/V∞），計算得到的時間步長為?t=1.076×10?4s，最終選擇時間步長為?t=1.0×10?4s。權衡計算精度和效率，每個時間步中內迭代步數設置為100。

3 結果分析

定義水翼的升力系數CL=FL/(0.5ρlV∞2A)、阻力系數CD=FD/(0.5ρlV∞2A)以及表征翼型上下表面壓力分布的壓力系數Cp=(p?pout)/(0.5ρlV∞2)，其中，FL和FD分別為水翼所受的升力和阻力，A為翼型投影面積。鑒于LES 方法計算得到的升阻力系數和壓力系數具有隨機振蕩特性，故定義了多個周期內的時均升力系數CˉL、時均阻力系數CˉD和時均壓力系數Cˉp進行對比分析。定義數值計算結果與試驗結果的誤差為：

在進行空化數值模擬前，采用第22 屆ITTC推薦的收斂比Rk對基于LES 方法的Twist11N 水翼濕流場數值模擬的網格收斂性進行判斷，使用3 套不同密度網格計算得到的時均阻力系數為單調收斂，而時均升力系數為振蕩收斂，滿足網格收斂性要求。

3.1 網格密度對空化數值模擬的影響

基于WALE 亞格子模型的LES 方法，使用G1，G2 和G3 網格對Twist11N 水翼在σ=1.07時的空化現象展開數值模擬，分別探討不同網格密度對空泡脫落頻率、云狀空泡形態和水翼升阻力系數脈動特性等的影響。

使用3 套網格計算得到空泡體積時歷曲線，經傅里葉變換得到空泡脫落頻率f，如圖4 所示。結果顯示，3 套網格計算的空泡脫落頻率分別為30.14，31.78，31.24 Hz，與試驗結果[4]（f=32.55 Hz）的誤差分別為?7.40%，?2.37%，?4.02%。由此可見， G1 粗網格計算得到的空泡脫落頻率與試驗值的誤差較大，而G2 與G3 計算得到的空泡脫落頻率誤差明顯降低，且在現有計算設置下，網格的進一步加密并沒有使空泡脫落頻率的誤差減小，反而更加偏離試驗值。

圖4 不同網格密度下計算的空泡脫落頻率Fig. 4 Shedding frequency of cavity for various grid schemes

圖5 給出的是一個周期內不同網格密度下預報的空化演變過程，采用汽相體積分數α=0.1等值面呈現的典型時刻（1～9）云狀空泡形態，并與試驗結果[4]進行對比。在所示俯視圖中，來流方向為由上至下。

圖5 不同網格模擬的空化演變過程與試驗結果對比Fig. 5 Comparison of cavitation evolution by various grid schemes and experimental results

由圖可見，與試驗結果相比，3 套網格基本都能模擬出非定常云狀空化演變過程的典型特征。具體過程包括：

1） 主脫落：附著片空泡發展到一定范圍后，發生了破裂并脫落泄出大尺度空泡云團。

2） 馬蹄形云空泡形成：破裂后的附著片空泡呈現凹型結構并繼續發展，而脫落泄出的空泡云團則發生卷曲，形成馬蹄形結構。

3） 二次脫落：馬蹄形云空泡在主流體的作用下向下游移動，而附著片空泡在“耳垂”位置發生了小尺度結構的脫落。

4） 空泡潰滅：二次脫落的小尺度空泡云團發生潰滅，附著片空泡由凹型發展成為凸型結構；主脫落的大尺度馬蹄形云空泡在下游高壓區發生潰滅，同時，發展起來的附著片空泡再次發生破裂，新的周期開始。

由圖5 的試驗結果還可以觀察到，受回射流機制的影響，在時刻1 和時刻9 即將發生主脫落的片空泡表面較不規整，其它時刻的片空泡表面較為“光滑”，3 套網格大致都模擬出了上述這種現象，且隨著網格的加密，片空泡上表面的結構形態逐漸細化，更加接近試驗中相應區域的空泡形態。然而，數值模擬結果與試驗結果相比，仍存在一些形態上的差異。例如，在試驗中觀測到的片空泡呈現出條狀或泡狀形態，但數值模擬結果并未捕捉到這些特征，主要是因為計算得到的空泡形態是由連續相分數界定得到的。

因網格密度不同，數值模擬空化演變過程中捕捉到的空泡形態存在一定差異，主要體現在：

1） 空泡長度Lc（如圖5(b)中時刻8 所示）：從圖6 給出的不同網格密度下數值計算得到的Lc/C值在一個完整周期內隨時間的變化情況可見，Lc生長速度隨時間的變化存在減慢的趨勢，這是因為片空泡發展到一定長度后，在回射流和主來流作用下翼展中部片空泡被抬升翻卷，阻礙了片空泡繼續往下游方向發展。而且，由G1 網格預報的Lc/C值明顯低于G2 和G3 網格相應時刻的預報值，且后二者數值結果與試驗結果[4]吻合較好，其中，與G2 網格數值計算結果相比，采用G3 網格在多數時刻點的Lc/C預報值與試驗結果的偏差程度更大，這從圖5 中的空泡形態上也可以明顯看出。

圖6 不同網格密度下計算的Lc/C 值對比Fig. 6 Comparison of cavity length-chord ratio for various grid schemes

2） 二次脫落：從空泡形態演變過程可見，G1網格數值計算結果在“耳垂”位置的空泡現象較簡單，G2 和G3 網格計算結果在“耳垂”位置發生多次二次脫落的現象更明顯；從圖5(c) 中的時刻1和時刻8 箭頭所指位置可見，在網格加密后可以捕捉到更精細的空泡現象。

3） 馬蹄形云空泡尺度：對于主脫落泄出而形成的馬蹄形云空泡，從圖5 中3 套網格模擬結果可見，在時刻3，脫落泄出的云空泡上部都有一個V 型的空泡結構，形態也都較為相似；在時刻4，可以明顯觀察到馬蹄形云空泡的雛形，其中G3網格數值結果中時刻4 的馬蹄形云空泡傾斜程度不同于G1，G2 網格的結果。

總體上，在時刻1～時刻4 之間，3 套網格模擬的云空泡現象較為一致。在時刻4 之后，云空泡演變過程中伴隨有收縮和潰滅現象，空泡潰滅形成的小尺度空泡與網格尺度量級更接近，網格尺度不同造成了包括流場中可捕捉到的漩渦尺度等流體湍動特征在預報精細程度上存在差異。可見，在時刻4 之后，馬蹄形云空泡的形態差異較明顯，在其展向尺度穩定后，G1 網格模擬的展向尺度小于試驗結果，而G2，G3 網格模擬的結果都與試驗結果相近。

圖7 和圖8 分別給出了3 套網格數值模擬計算得到的一個周期內升阻力系數時歷曲線。圖中，以無量綱時間t/T為橫坐標，T為相應網格密度下的周期。由圖可見，G1 網格數值計算的升阻力脈動特征相對較少，3 套網格的升阻力瞬態脈動位置較為一致，且主要集中在t/T=[0.2,0.3]及[0.8,1.0]這兩個區間范圍；3 套網格數值計算的升阻力脈動差異主要集中在后半周期，與馬蹄形云空泡的潰滅現象密切相關。此外，對比3 套網格所模擬的馬蹄形云空泡的演變過程，還可發現G2 和G3 網格所模擬的馬蹄形云空泡潰滅及回彈的現象比G1 網格的更明顯（如圖5 中時刻5～時刻8），這也解釋了G2 和G3 網格所預報的升阻力脈動更劇烈這一現象的原因。

圖7 不同網格密度下一個周期內升力系數時歷曲線Fig. 7 Time histories of the lift coefficients for various grid schemes in a shedding cycle

圖8 不同網格密度下一個周期內阻力系數時歷曲線Fig. 8 Time histories of the drag coefficients for various grid schemes in a shedding cycle

表2 不同網格密度下時均升阻力系數計算結果的對比Table 2 Numerical results of the time-averaged lift and dragcoefficient for various grid schemes

圖9 和圖10 分別給出了空泡流中不同網格密度下，展向YC=0.12 m 和YC=0.15 m 位置處水翼表面的時均壓力系數?Cˉp與兩組試驗（即Delft和EPFL 這兩個不同空泡水筒所做的試驗）數據的對比，其中Delft 試驗水翼的尺寸是EPFL的兩倍，二者以相同雷諾數進行試驗。由于部分壓力監測點的試驗結果帶有一定的不確定性[20]，在對比時，既與相同尺寸的Delft的試驗結果[18]作比較，也與EPFL的結果[18]作比較。由兩個圖可知，在水翼展向兩個截面處，不同網格密度下壓力面的時均壓力系數計算結果基本一致，而吸力面壓力的結果略有差異。主要體現在：

圖9 YC=0.12 m 展向位置處水翼表面 的分布Fig. 9 The time-averaged pressure distribution around hydrofoil surface at YC=0.12 m

圖10 YC=0.15 m 展向位置處水翼表面 的分布Fig. 10 The time-averaged pressure distribution around hydrofoil surface at YC=0.15 m

1）在0

2）在0.4 ≤x/C≤0.5區域，作為片空泡末端，該區域壓力沿弦長方向分布變化較大，對于YC=0.12 m 展向位置，G3 網格所預報的壓力轉變趨勢與試驗結果更接近，而對于YC=0.15 m 展向位置，G2 網格較好地模擬了x/C=0.4 處的壓力轉變趨勢。

3）在0.5

3.2 亞格子模型對空化數值模擬的影響

基于中密度網格G2 的LES 方法，采用WALE，WMLES 和SL 這3 種亞格子模型對Twist11N 水翼σ=1.07時的空化現象展開數值模擬，分別探討不同亞格子模型對空泡脫落頻率、云狀空泡形態以及水翼升阻力系數脈動特性等的影響。

圖11 給出了3 種亞格子模型計算所得的空泡脫落頻率f。WALE，WMLES 和SL 這3 種亞格子模型計算得到的空泡脫落頻率分別為31.78，29.19 和34.87 Hz，與試驗結果[4]（f=32.55 Hz）的誤差分別為?2.37%，?10.32% 和7.13%，其中WALE亞格子模型預報的脫落頻率誤差最小。

圖11 不同亞格子模型計算的空泡脫落頻率Fig. 11 Shedding frequency of cavity for various subgrid-scale models

圖12 給出了一個周期內不同亞格子模型所預報的空化演變過程，同樣采用汽相體積分數α=0.1等值面呈現典型時刻云狀空泡形態，并與試驗結果進行對比。

圖12 不同亞格子模型模擬的空化演變過程與試驗結果對比Fig. 12 Comparison of cavitation evolution by various subgridscale models and experimental results

由圖可見，與試驗結果相比，3 種亞格子模型基本都能夠模擬出非定常云狀空化演變過程的典型特征，但模擬得到的空泡形態差別較大，主要包括：

1） 片空泡閉合曲線輪廓形狀：對比片空泡的弦向范圍（如圖12 中所示時刻9），可以發現WALE亞格子模型計算的結果與試驗結果較吻合，片空泡范圍最大時，弦長方向可延伸至弦長的50%左右，相應地，WMLES 和SL 亞格子模型計算的結果則偏小，只延伸至弦長的40%左右，但三者導邊附近的片空泡展向延伸的長度相當。此外，由圖13 給出的一個周期內不同亞格子模型計算的Lc/C值隨時間的變化情況可知，WMLES 與SL 亞格子模型計算的Lc/C值相當，但明顯低于試驗結果[4]，而WALE 亞格子模型計算得到的Lc/C值與試驗結果吻合較好。

圖13 不同亞格子模型計算的Lc/C 值對比Fig. 13 Comparison of cavity length-chord ratios for various subgrid-scale models

2） 二次脫落：觀察圖12 中時刻3 至時刻7 的試驗結果，可見片空泡“耳垂”位置受到回射流擾動而呈現出不穩定性的區域逐漸擴大，但WMLES與SL 亞格子模型計算結果在該區域上的空泡結構與試驗相比要稀疏得多，特別是SL 亞格子模型，而WALE 亞格子模型計算結果在該區域觀察到了更豐富的小尺度空泡結構。

圖 14 YC=0.12 m 展向位置處水翼表面的分布Fig. 14 The time-averaged pressure distribution around hydrofoil surface at YC=0.12 m

3） 馬蹄形云空泡：從俯視角度看，WALE 和WMLES 亞格子模型得到的馬蹄形云空泡的尺度范圍與試驗結果吻合較好，而SL 亞格子模型計算得到的馬蹄形云空泡明顯小于試驗結果。對于翼展中部馬蹄形云空泡的發展而言，WMLES 亞格子模型在時刻1 能夠模擬出上一個周期主脫落形成的馬蹄形云空泡即將在翼面上通過隨邊的現象，而另兩個亞格子模型得到的馬蹄形云空泡則提前發生了潰滅。

表3 不同亞格子模型時均升阻力系數計算結果的對比Table 3 Numerical results of the time-averaged lift and drag coefficients for various subgrid-scale models

圖15 YC=0.15 m 展向位置處水翼表面的分布Fig. 15 The time-averaged pressure distribution around hydrofoil surface at YC=0.15 m

4 結 論

本文基于LES 方法和S-S 空化模型，通過3 套網格密度及3 種亞格子模型，研究了Delft Twist11N 三維扭曲水翼的空化演變過程、空化脫落頻率及升阻力系數等，得到如下結論：

1） 在最大壁面y+<1 的前提下，當弦長方向平均網格無量綱間距滿足 ?x+≤94 條件時，所用網格密度下的數值模擬既能夠捕捉到較多的細小空泡脫落和馬蹄形云空泡的初生及潰滅等非定常空化演變現象， 又能夠獲得具有一定精度的空泡脫落頻率、時均升阻力系數和時均壓力分布。

2） 3 種亞格子模型均能夠模擬非定常云狀空化演變過程中的典型特征。相較于WMLES 與SL 亞格子模型，WALE 亞格子模型捕捉到的片空泡及云空泡的演變與試驗結果更吻合，且預報的空泡脫落頻率、時均升阻力系數和時均壓力系數精度更高。

因此，推薦選擇基于WALE 亞格子模型的LES方法進行非定常云狀空化現象的數值模擬。