文丘里管內空化流動大渦模擬的驗證與確認

梁蘊致，龍云，龍新平*，程懷玉

1 武漢大學 動力與機械學院，湖北 武漢 430072

2 武漢第二船舶設計研究所，湖北 武漢 430205

3 武漢大學 水利水電學院，湖北 武漢 430072

0 引 言

隨著計算機技術的發展，數值模擬的應用領域越發廣泛。為解決模擬結果的準確性問題，部分研究者從誤差的來源出發，提出了一套數值模擬結果的驗證與確認（verification and validation，V&V）方法。驗證用于說明數值方法求解結果與對應方程準確解的偏離程度，確認則用于說明數值解與真實流動情況間的偏離程度[1]。Schaefer等[2]認為數值計算過程中存在幾何簡化等模型形式、湍流模型系數等參數、離散方法等求解方法、邊界條件等輸入量、后處理插值等結果處理這5 個方面的不確定性因素，從而使得數值模擬結果與實際流動情況產生偏離。

陳鑫等[3]認為目前幾何外形、湍流模型、邊界條件等帶來的不確定度研究已較為成熟，但針對計算網格的不確定度研究仍然有待完善。目前，對于網格不確定度的研究主要基于估算離散誤差的Richardson 外推方法[4-5]，對多套系統加密網格的數值模擬結果進行處理，從而得出采用最密一套網格數值模擬結果的不確定度。另外，現有V&V 方法大多基于RANS 模型的穩態數值模擬，例如Roache[6-7]提出的網格收斂性指數（GCI）方法及其改進版本[8-9]；Wilson 等[10]提出的修正因子法（CF）；Tao 等[11]提出的安全系數法（FSI）。ASME等組織基于上述研究提出了各自的V&V 方法的實施規程[12-13]。張涵信等[14]另辟蹊徑，考慮了網格和加權平均誤差進行兩次加權平均，得出真值，從而得到不確定度。

隨著計算機性能的不斷提高，在水動力學領域，非穩態、多相流數值模擬研究越來越普遍，大渦模擬（large-eddy simulation，LES）也被越來越多地用于數值模擬研究中，但相關的V&V 方法不夠全面，仍有待進一步研究。對于LES，類比穩態RANS 模擬結果的V&V 方法，Vreman 等[15]指出，除網格不確定度外，亞格子模型不確定度也是LES 不確定度的重要組成部分。Klein[16]將這兩個準確度階數假定為2，給出了一種LES 不確定度的三方程計算方法。Xing[17]認為，LES 不確定度還應包括網格與模型耦合帶來的不確定度，提出了七方程和五方程的LES 不確定度計算方法。對于空化這一典型的包含相變的多相流動，Long 等[18]及程懷玉等[19]基于RANS 模型對5 套系統加密的網格在水翼空化流動中的不確定度進行了計算，發現非定常空化流動增加了計算誤差，從而增大了V&V 方法的實施難度。在Xing[17]對于無空化LES V&V 實施方法的研究基礎上，Long 等[20]綜合大量計算結果，針對水翼空化流動，將2 個準確度階數固定，從而將五方程模型簡化為三方程模型，并對可變螺距（CP）槳葉進行了五方程和三方程的數值模擬準確度計算[21]。

綜上，V&V 方法被廣泛應用于計算流體力學的各領域，是檢驗數值模擬準確度的有效手段，但現有V&V 方法主要基于RANS 模型的穩態外部流動模擬，LES V&V 實施方法的研究較少，尤其是對于有限空間內部空化流動LES V&V 實施方法仍有待研究。為此，本文將對3 種空化數工況下的文丘里管內流動進行數值模擬，基于五方程模型，對有限空間內部空化流運動LES 的驗證與確認展開研究。

1 數值模型與邊界條件

1.1 幾何模型與計算域

本文計算所用文丘里管截面圖與尺寸如圖1[22-23]所示，喉管直徑Dth=10 mm，喉管部分長度Lth=Dth，進口直徑Din=5Dth，同時進、出口尺寸相等，該文丘里管的收縮段收縮角為45°，擴散段擴散角為12°。

圖1 文丘里管截面圖Fig. 1 Sketch view of Venturi

本文計算域如圖2(a)所示，為保證流動的充分發展和計算的穩定性，在文丘里管的進口及出口位置均增設了5 倍進出口直徑的直管段。本文數值模擬中，以喉管進口截面中心為坐標起點，共設置4 個監測點，分別位于x/Dth=0.5，x/Dth=5，x/Dth=9，x/Dth=13 處的壁面位置，如圖2(b)所示。

圖2 計算域及監測點位置Fig. 2 Computational domain and positions of monitoring points

1.2 網格劃分

基于Richardson 外推方法發展的一系列V&V方法需要在同樣的工況下對多套系統加密的網格展開數值模擬計算[24]，網格的系統加密即基于同一套拓撲結構，采用恒定的比例增加各個方向網格數量。為了滿足五方程大渦模擬V&V 方法的實施要求，本文數值模擬共使用了5 套網格加細比r為1.2 的系統加密結構化網格，這5 套網格以最密一套網格為G1，最稀疏一套網格為G5，其網格單元數如表1 所示。

表1 網格單元數Table 1 The number of grid cells

同時，為滿足大渦模擬的計算要求，網格在喉管和壁面位置對網格進行了加密。圖3(a)為G1 喉管位置的壁面網格，從收縮管到喉管網格平滑過渡，尺寸不斷減小，從喉管到擴散管，網格尺寸逐漸擴大。圖3(b) 為G1 喉管中部截面網格，靠近壁面位置網格高度減小。

圖3 文丘里管局部網格劃分Fig. 3 Part of grid division of Venturi

1.3 控制方程

文丘里管內空化流動為典型的兩相湍流流動，本文采用Mixture 模型將氣液兩相假定為均質流動，無相間滑動，采用體積分數的變化模擬相關流動參數。湍流模型采用LES 方法，并用Zwart-Gerber-Belamri (ZGB)空化模型模擬氣液相變過程。 LES 濾波處理后的N-S 方程為：

式中：p為壓強；ui為i方向的速度，其中，在三維笛卡爾坐標系中，下標i 為x軸正方向，j為y軸正方向，k為z軸正方向； δij為Kronecker-delta 函數；ρm和μm分別為均質混合流體的密度和動力黏度，其中ρ1為水密度，μl為水動力黏度，αv為蒸汽體積分數，ρv為蒸汽密度；式(2)帶有波浪型上劃線的項表示濾波處理項，而τij為亞格子應力，代表被過濾的小尺度脈動的作用。基于Boussinesq 假設，亞格子應力項如式(5)所示，其中μSGS為 亞格子湍流黏度，S?ij為應變率張量。本文采用的亞格子應力模型為Dynamic Smagorinsky-Lilly 模型。

空化流動模擬方面，大量算例證明，基于Rayleigh-Plesset 方程得出的ZGB 空化模型[25]有較好的普適性，且模擬精度和收斂性均較好，復雜流動中也表現出了較好的收斂性，故本文也選用ZGB 模型模擬空化過程。ZGB 模型的質量輸運方程為：

式中，me和mc分別為蒸發源項和凝結源項，表達式為：

式中：pv為飽和蒸汽壓；RB為氣泡直徑；Fvap和Fcond為經驗系數；αnuc為成核點的體積分數。

1.4 邊界條件

本文仿真分別計算了空化數為2.086，0.765及0.263 這3 種工況。空化數定義為：

式中：pout為出口壓力；uth為喉管中間截面平均流速。當空化數σ 為2.086 時，試驗中未觀察到空化現象，而空化數為0.765 及0.263 時，試驗中均發現有較為劇烈的空化現象，即本文數值模擬同時涵蓋了空化和無空化流動。

本次計算涉及的流體為水和水蒸氣，均采用25 ℃時的物性參數，飽和蒸汽壓為3 540 Pa，水的密度為997 kg/m3，黏度為8.899 ×10?4kg/(m?s)，水蒸汽的密度為0.023 08 kg/m3。

表2 為本文數值模擬的邊界條件，其來源為試驗測量數據[22-23]。數值模擬時，進口設置為速度進口，出口設置為壓力出口，壁面設置為無滑移壁面。對應于本文數值模擬采用的網格加細比，當采用更細密的一套網格時，時間步長除以1.2，即G5 網格采用的時間步長為1.44×10?5s，G4 網格采用的時間步長統一為1.2×10?5s，依此類推。

表2 計算邊界條件Table 2 Boundary conditions

本文計算使用CFX 軟件，并采用192 線程服務器進行并行計算，內存128 G，計算共耗時1 440 h。

為排除時間步長的選取對非穩態數值模擬計算結果的影響，選取G4 網格的無空化工況，另外選取了1.44×10?5s 以及1×10?5s 兩個時間步長，進行了時間步長無關性驗證。相比于時間步長為1.2×10?5s，兩時間步長下進出口壓力比m的變化均在0.5%以下，4 個監測點的壓力系數Cp變化也都在3%以下，可以認為，本次數值模擬的各套網格滿足時間步長無關性。

壓力比m是文丘里管主要性能指標之一，為出口壓力pout與進口壓力pin之比，定義如式（10）。對于進出口等徑的文丘里管，該參數表現了文丘里管中的壓能損失。壓力系數表達式如式(11)。

式中，pmon為瞬態計算結果趨于穩定后最后2 個周期內的監測點壓力平均值。

2 計算結果與試驗對比

2.1 進出口壓力比對比

圖4 為3 種空化數下，文丘里管壓力比的試驗值及5 套網格的數值模擬結果[22-23]，而G1～G5分別為從最密一套網格到最稀疏一套網格共計5 套網格的模擬值。

圖4 5 套網格模擬壓力比與試驗值的對比Fig. 4 Comparison of pressure ratios between experimental results and simulations in five grids schemes

圖4(a)中，空化數σ=2.086 時，隨著網格加密壓力比上升，在網格G2 處達到峰值。圖4(b)中，空化數σ=0.765 時，模擬結果不斷降低，逐漸向試驗值靠攏。圖4(c) 中，空化數σ=0.263時，壓力比變化幅度最小。所有網格計算結果與試驗值間的誤差均在5%以內，表明本文數值模擬的3 種空化數下，壓力比這一全局量的模擬結果與試驗值均吻合較好。

2.2 監測點壓力系數對比

文丘里管幾何模型較為簡單，但其獨特的收喉擴結構使喉管和擴散管部分的流場較為復雜，對于這一典型部分內部流動結構的研究很有必要。圖5 分別為σ=2.086，0.765，0.263 時，G1～G5網格下壓力系數Cp的模擬值與試驗結果。

圖5 文丘里管壁面監測點壓力系數Fig. 5 Pressure coefficient of monitoring points on the wall of Venturi

圖5(a)和圖5(b)中，隨著網格數的增加，數值模擬結果更加接近試驗值，最密一套網格G1 計算結果與試驗值間的誤差均很小，與試驗結果吻合較好。而圖5(c)中，x/Dth=0.5，x/Dth=5 監測點壓力系數誤差很小，而監測點x/Dth=9，x/Dth=13 處誤差相對較大。空化較弱或無空化時，各監測點壓力系數模擬值與試驗值均較為吻合，而空化較劇烈時，喉部和擴散管前部監測點仍較為吻合，僅擴散管中后部2 個監測點壓力系數偏差較大。

上述對比表明，本文所計算的3 種工況下，進出口壓力比這一全局量的試驗值與模擬值均較為吻合；對比各監測點壓力系數這一局部量時，僅空化數σ=0.263 時，擴散管中后部2 個監測點的壓力系數誤差較大，其余各監測點和各工況下的模擬結果與試驗值均較為吻合。即當無空化或空化較弱時，本文所采用數值模擬方法能夠較為準確地反映實際流動；空化較劇烈時，部分局部量偏差較大，這與求解方法和空化模型都有一定關系，有較大的改善空間，但全局量仍較為符合。

3 V&V 實施方法與結果

3.1 驗 證

LES 的誤差可以分為數值誤差δSN及模型誤差δSM兩個部分[17]：

目前，廣泛認為由直接數值模擬（DNS）方法得到精度較高的離散解可以作為數值基準值，但并非所有問題都能夠采用DNS 求解，故在大部分問題中， ?c為需要求解的未知量。根據二元泰勒級數展開法及廣義Richardson 外推方法，同時參考已有規程[12-13]對RANS 模型的誤差計算方法，展開兩誤差項并略去高階量，則LES 方法的計算誤差可以表示為

式中：cN和cM分別為數值誤差系數和模型誤差系數；pN和pM分別為數值誤差準確度階數和模型誤差準確度階數；Δn為濾波尺寸，下標n表示對應數值模擬所采用的網格；h?n為綜合考慮當地網格尺寸和瞬態計算時間步長得到的時空尺度，表達式如下：

式中：hn為當地網格尺寸；Δtn為時間步長。

對于系統加密的幾套網格，當采用的加細比為r時，則計算時相鄰網格有如下關系：

如要求解式(13)，理論上應有至少5 個等式組成方程組。本文以G1 網格的當地尺寸為網格尺寸h，采用該套網格時，數值模擬的時間步長為Δt，由5 套系統加密網格得出的數值解為 ?n，可構成的方程組如下：

方程組為復雜非線性方程，借助Matlab 最優化工具求解[26]，采用信任域方法和Jacobian 矩陣迭代法，得出壓力比m和各監測點壓力系數Cp的準確度階數如下。圖6 為壓力比m的準確度階數。如圖所示，當σ=2.086，0.263 時，pN，pM比較接近，但σ=0.765 時，pN，pM明顯低于其他2 個工況。圖7 為Cp的準確度階數。如圖所示，同一工況下，不同監測點處的準確度階數有較大差別，且同一監測點處，不同工況下的準確度階數也不相同。

圖6 壓力比準確度階數Fig. 6 The observed order of accuracy of pressure ratio

圖7 壓力系數準確度階數Fig. 7 The observed order of accuracy of pressure coefficient

文獻[21]中，pN，pM在各個工況下相差較小，可以固定這2 個準確度階數并將五方程求解簡化為三方程；但對于文丘里管這種速度、壓力變化劇烈的內部流動，準確度階數會隨監測點位置、流動工況的改變而改變，且變幅較為明顯。故該方法存在不足，需要尋找準確度階數隨監測點位置和流動工況變化的規律，從而加以改進。

LES 驗證需要進一步計算驗證不確定度，計算方法為

式中：ULES為LES 的驗證不確定度，表明 LES 結果與精確解之間的偏離程度；FSLES為不確定度安全系數，根據Roache[6]的推薦，該系數可以取為1.5 或者3。考慮到文丘里管內空化流動為三維非定常相變過程，計算中會帶來較大的不確定性，故當σ=0.765，0.263 時，取安全系數為3，而σ=2.086 時，取安全系數為1.5。

進出口壓力比m這一全局量的不確定度ULES均較小（圖8），σ=0.765 時，ULES相對較大；對應于圖4 中，σ=0.765 時模擬結果與試驗值之差也最高，這說明驗證不確定度ULES的大小能夠反映壓力比m計算誤差的大小。壓力系數Cp的驗證不確定度如圖9 所示，在σ=2.086，0.765 工況下，4 個監測點位置的驗證不確定度ULES均較小；當σ=0.263 時，x/Dth=0.5 及x/Dth=5 這2 個監測點壓力系數的ULES較小，而x/Dth=9 和x/Dth=13 這2 個監測點處ULES明顯增加。對應圖5(c)，這2 個監測點處的計算值與試驗值的誤差也偏大，且誤差最大的x/Dth=9 處，ULES也最大，即驗證不確定度ULES的大小也能夠很好地反映壓力系數Cp的誤差。

圖8 壓力比驗證不確定度Fig. 8 The verification uncertainty of pressure ratio

圖9 壓力系數的驗證不確定度Fig. 9 The verification uncertainty of pressure coefficient

3.2 確 認

驗證結束后，還需要對LES 結果進行確認，即計算并對比分析確認不確定度和對比誤差。Long 等[21]指出，相較于RANS 模型的驗證過程，LES 驗證過程已經將模型誤差納入了考慮范圍，因此 LES 的確認不確定度可以表示為

式中：Uv為確認不確定度；UD為試驗的不確定度。本文中UD為試驗所用壓力測量儀表的不確定度，壓力漂移值為壓力表量程的0.5%。

當對比誤差（試驗值D與LES 模擬值之差）小于確認不確定度時，則認為在該不確定度Uv的水平下，LES 的確認得以實現。對比誤差為

壓力比m的確認不確定度及對比誤差E如圖10所示。在3 種空化數的工況下，對比誤差E均遠小于確認不確定度Uv，即3 種空化數工況下的進出口壓力比m均實現了對應不確定度下的確認。

圖10 壓力比的確認不確定度及對比誤差Fig. 10 The validation uncertainty and comparison error of pressure ratio

各監測點處壓力系數的確認不確定度與對比誤差如圖11 所示。在空化數σ=2.086，0.765 時，各監測點壓力系數的對比誤差均小于對應的確認不確定度；而當σ=0.263 時，x/Dth=9 和x/Dth=13 這2 個監測點處的對比誤差大于確認不確定度。即σ=0.263 時，x/Dth=9 和x/Dth=13 這2 個監測點處的壓力系數沒有實現對應不確定度水平下的確認，其余各工況和各監測點的壓力系數均實現了對應不確定度水平下的確認。

圖11 壓力系數的確認不確定度及對比誤差Fig. 11 The validation uncertainty and comparison error of pressure coefficient

在空化數σ=0.263 時，文丘里管處于劇烈的空化狀態，其蒸汽體積分數云圖、試驗高速攝像和截面網格如圖12 所示。圖12 (a)和圖12(b)中，空化云延伸范圍較長，x/Dth=9 監測點位于空化云的末端，在此位置空化云脫落潰滅，同時，部分脫落的空化云還會運動到x/Dth=13 位置處再潰滅，空化云的脫落潰滅在壁面位置造成了不穩定回流，進而導致這2 個位置處的不確定度明顯上升；此外，監測點x/Dth=9 和x/Dth=13 位于擴散管的中后部，說明當空化數σ=0.263 時，本文的數值方法仍有較大的改善空間，需要對網格及求解方法做進一步調整。

圖12 σ=0.263 時文丘里管空化云形態及部分G1 網格Fig. 12 Cavitation clouds and part of G1 grid scheme in Venturi at σ=0.263

4 結 論

本文采用LES 耦合ZGB 空化模型，對3 種空化數工況下的文丘里管內瞬態空化流動進行了數值模擬，運用Matlab 最優化工具求解五方程V&V 方法中的數值基準和準確度階數，進而計算出驗證不確定度ULES、確認不確定度Uv和對比誤差E，進行了文丘里管內部空化及無空化流動大渦模擬的驗證與確認，得出如下結論：

1) 最密的一套網格數值模擬結果與試驗結果的直接對比表明：3 種空化數工況下的壓力比這一全局量的模擬值與試驗數值均較為吻合。無空化或空化較弱時，監測點壓力系數與試驗結果均吻合較好；而當σ=0.263，空化較為劇烈時，x/Dth=9 和x/Dth=13 處的壓力系數誤差偏大，其余2 個監測點數據吻合較好，數值模擬方法有較大的提升空間。

2) 驗證結果表明，壓力比這一全局量的驗證不確定度均較小；當空化較弱或無空化時，監測點壓力系數驗證不確定度均較小，而空化較為劇烈時，擴散管中后部監測點壓力系數的驗證不確定度較大。總體上驗證不確定度較大時，對比誤差也較大，說明驗證不確定度能夠反映數值模擬的準確程度。

3) 確認結果表明，在大部分工況和監測點處，均實現了對應驗證不確定度水平下壓力比的確認；對比監測點壓力系數時，僅在空化數σ=0.263，監測點x/Dth=9 和x/Dth=13 處沒有實現確認，說明空化數適中或無空化時，使用細化的網格，LES 方法能較好地模擬非定常空化流動。但空化數較小和空化極為劇烈時，需要針對該工況改善網格，并調整求解方法，以精確捕捉空化的演化過程。

4) 目前，實施有限空間內部空化流動的LES驗證與確認時，壓力比等全局量的驗證與確認較易實現，但空化較為劇烈時，完全實現各監測點壓力系數的驗證與確認還需要繼續完善數值模擬求解方法，未來還將進一步修改空化模型，使之符合更大范圍的空化工況。此外，為實現五方程V&V 方法，一個工況需要起碼5 套系統加密網格的計算結果，計算資源消耗較大，而三方程方法還需尋找準確度階數隨監測點位置和流動工況變化的規律，以做進一步改進，拓展V&V 方法的應用范圍。