基于拉格朗日方法的通氣云狀空泡三維非定常脫落特性研究

段忠平，劉濤濤，張孟杰，黃 彪，王國玉

(北京理工大學機械與車輛學院, 北京 100081)

0 引言

對于水下高速航行體而言，為了有效利用通氣超空泡的減阻作用，必須在水下航行體表面形成特定尺度、形態可控、且相對穩定的超空泡。然而在大多數情況下，物體是從全濕流過渡到局部空泡，最后才發展為超空泡，因此通氣空泡的發展是一個漸近穩定過程：由低頻振蕩云狀空泡、霧狀空泡，逐漸轉變成半透明的空泡，并在通氣量越過某一臨界值后，生成包圍航行體的穩定超空泡[1]。

通氣云狀空泡作為通氣超空泡發展過程中的中間形態，由于空泡尺度較小，空泡發展過程容易受到流動參數、通氣參數、重力效應以及相應的流體屬性等多種因素影響而出現空泡團的斷裂脫落現象，導致流場發生強烈變化，這種空泡的脈動和斷裂等非定常流動特性一直是通氣空泡研究的熱點問題。Rashidi等[2]運用高速攝影技術觀測了通氣空泡的流動特性，指出反向射流對空泡界面、空泡的斷裂以及空泡內部的流動具有較大影響。Kinzel等[3]采用數值計算對繞回轉體通氣空泡流動特性進行了研究，闡明了反向射流是引起大尺度空泡脫落的主要原因。王國玉等[4-5]結合實驗和數值計算，基于對空泡尾部閉合區域處的壓力分析發現，空泡內外壓差形成的逆壓梯度引起空泡閉合處出現流動分離，形成回射流，回射流與空泡內部氣體相互作用產生復雜的旋渦結構，此旋渦結構與主流相互作用造成空泡斷裂脫落。王一偉等[6-7]基于反向射流的運動對通氣云狀空泡脫落機制進行了詳細分析，發現通入氣體造成反向射流呈蛇形向上游運動，在運動過程中于壁面處產生高壓，隨著空泡的繼續發展和反向射流的持續運動，兩個高壓區擴大并逐漸連接，促使空泡在高壓區的連線形成斷裂。這些研究基于對回射流的分析，促進了對通氣空泡非定常特性的認識。然而對于軸對稱回轉體而言，物體往往存在很強的三維效應，因此形成的通氣空泡流動也具有三維特性，進而造成反向射流的運動特征也發生相應的改變。Lindau等[8]基于數值計算對繞三維平頭回轉體的空化流動特性進行了研究，分析了空泡產生的起始點位置，并對空泡尾部的反向射流運動過程進行了描述，說明了其運動的三維特性。Hu等[9-11]基于混合密度分域的PANS湍流模型對平頭回轉體的云狀空化流動進行了數值計算研究，發現云狀空穴形態沿回轉體周向的分布始終呈現不對稱，導致反向射流的推進速度在圓周方向上存在差異，進而造成周向上反向射流到達回轉體肩部的時間不一致，引起空泡發生不規則斷裂。

由于通氣空泡流動的復雜性，限制了人們對其非定常流動機理的有效認識。因此，為了進一步詳細分析通氣空泡的三維非定常脫落特性，本文基于CFX軟件對繞軸對稱回轉體通氣云狀空泡發展演化過程進行了數值模擬，并基于FTLE和LCS方法對非定常通氣空泡流動進行描述，采用粒子追蹤方法對特征流動結構進行分析，揭示通氣云狀空泡的三維非定常脫落特性。

1 數值計算方法

1.1 基本控制方程

研究基于Navier-Stokes方程的計算框架，采用均質平衡流模型，則Favre平均的N-S方程連續性方程和動量方程可表述為

(1)

(2)

ρm=ρlαl+ρgαg

(3)

式中，下標i和j分別代表坐標方向;ρl，ρg分別為液相和氣相密度;αl，αg分別為液相和氣相體積分數;ρm，u和p分別為混合介質的密度、速度和壓強;μ和μt分別為混合介質的層流和湍流黏性系數。

1.2 湍流模型

由Johansen等[12]提出的濾波器湍流模型(Filter based model，FBM)通過在流域中加入濾波器函數，湍流尺度小于濾波尺度時計算采用k-ε湍流模型進行求解；當濾波尺度小于RANS模型中的湍流特征尺度時采用LES進行計算，能較好地捕捉非定常流動引起的大尺度旋渦結構。

(4)

(5)

湍流黏性由式(6)確定

(6)

式中，k，ε分別為湍動能和湍流耗散率，Pt為湍動能生成項，λ為濾波尺度。相關模型常數分別為：Cε1=1.44，Cε2=1.92，σε=1.3，σk=1.0。

1.3 網格及邊界條件設置

為了與實驗結果對比，模型的幾何尺寸及空間相對位置與實驗[13]保持一致。采用的模型為40°鈍錐頭型回轉體，其長度L=126 mm，直徑D=20 mm，計算域入口距回轉體頭部約為12D，出口距回轉體尾部為17D，上下壁面距離為4.75D。計算中，基于CFX平臺，采用速度入口和壓力出口作為邊界條件，通氣孔設置為質量流量入口，如圖1所示?；剞D體表面采用絕熱、無滑移的固壁條件。根據實驗條件，入口速度設定為U∞=6.93 m/s，Fr=15.66(Fr=U∞/(gD)1/2)，出口壓力設定為p∞=101.325 kPa，空化數σ=4.52(σ=(p∞-pv)/0.5ρv2),通氣孔設置為質量流量入口，CQ=0.06(CQ=Qin/(U∞D2))。采用包裹回轉體的C型網格拓撲結構來控制回轉體周圍的網格分布，根據實驗結果中通氣空泡的區域來確定C型網格拓撲區域的范圍。為了可以準確捕捉通氣空泡區域內的流場信息以及減少整個計算域的網格數量，包裹回轉體C型網格拓撲區域的網格尺度在0.1～0.2 mm，其他區域的網格尺度逐漸增大，網格總數為1.3×106。

圖1 計算域邊界條件與網格劃分Fig.1 The boundary conditions and mesh used in the CFD computations

2 拉格朗日擬序結構

拉格朗日擬序結構(Lagrangian Coherent Structures, LCS)指的是有限時間李雅普諾夫指數(Finite-Time Lyapunov Exponent, FTLE)代表的標量場中出現的“脊”結構，它代表該區域中無限相鄰的質點，經過流場的拉扯、卷曲和交錯等，將發散至較大的距離。該方法以流體質點在流場中對流特性的時間積分為主要的考慮因素，具有拉格朗日性質的物質線顯示旋渦邊界，對流場渦的邊界和位置有整體的客觀描述，可以描述流體質點的真實運動軌跡[14]。

假設流體質點在系統空間內的初始時刻t=t0所處位置為(x0,y0,z0)，經過T時間的運動之后所處位置為(xT,yT,zT)，則有限時間李雅普諾夫指數可以定義為

(7)

(8)

(9)

本次研究中，首先基于CFX商業軟件獲得通氣云狀空泡非定常流動流場數據，然后基于Fortran程序對計算所得到的速度場進行處理，獲得FTLE結果，然后通過Tecplot軟件對計算結果進行顯示。

3 結果與分析

3.1 通氣云狀空泡的發展演化過程

圖2分別給出了數值計算得到的通氣云狀空泡

形態隨時間的演化過程，并和實驗結果的對比情況。從圖中可以看出，數值計算結果和實驗結果吻合較好，可以較為準確地描述通氣云狀空泡發展演化過程中反向射流的推進、空泡斷裂、斷裂空泡團旋轉翻滾以及大尺度空泡團脫落等非定常特性。在該工況下，通氣云狀空泡的發展周期約為T=17 ms。在t0時刻，反向射流運動到空泡中間位置，此時空泡前端為連續透明的氣相區，尾部為水氣混合的云霧狀，并伴隨著上個周期脫落的大尺度云狀空泡團。從t0到t0+0.32T時刻，反向射流持續向上游運動，空泡尾部的水氣混相區域逐漸擴大，而透明狀的氣相區域則逐漸縮小，同時空泡尾部出現小尺度空泡團脫落。在空泡前緣附近，反向射流與通入的氣流相遇，兩股方向相反的流動相互作用，在t0+0.74T時刻，空泡表面出現凹陷和漲起，發生不規則的斷裂，漲起空泡在向下游運動的過程中不斷翻滾，直至t0+T時刻在空泡尾部再次形成大尺度U型空泡團的脫落現象。值得注意的是，通氣云狀空泡的整個發展過程中，空泡尾部始終伴隨著隨機性較強的小尺度云狀空泡團的脫落現象，同時由數值結果中的三維效果圖可以看出，繞回轉體的通氣云狀空泡呈現很強的不對稱性，且周向不同位置出現凹陷裂痕的時刻存在差異。

圖2 空泡形態隨時間的演化過程
Fig.2 The time-evolution process of ventilated cavity

為了進一步說明通氣云狀空泡的發展演化過程，圖3給出了不同的軸向位置監測點的表面壓力系數隨時間的變化曲線，其中監測點的位置如圖4所示。從圖3中可以看出，監測點x/D=0.5，x/D=0.75和x/D=1的壓力保持相對穩定，與其他監測點的變化趨勢存在顯著差異，這是由于這3個監測點位于通氣孔附近，受到尾部空泡脫落影響較小，始終位于空泡內部。隨著監測點向空泡尾部移動，空泡末端周期性大尺度空泡團的斷裂、脫落影響逐漸加劇，對應的壓力也出現了低頻、高幅度的變化。基于上述的分析可知，對于繞回轉體的通氣空泡而言，反向射流存在強烈的三維特性，反向射流與主流相互作用后并未造成空泡完全斷裂，而是表現為空泡局部斷裂，當裂痕發展至空泡尾部后出現脫落，因此圖3中只有空泡尾部監測點x/D=3.5，x/D=4壓力出現高幅度的變化。

圖3 不同軸向位置壓力隨時間的變化曲線Fig.3 Predicted pressure distribution at different axial positions

圖4 回轉體表面監測線位置示意圖Fig.4 Sketch map of line-monitors at different positions of the axisymmetric body

為了研究通氣空泡發展過程的三維特性，圖5給出了不同時刻下不同軸向監測線上的表面壓力分布，其中4條監測線(A，B,C,D)的位置如圖4所示。從圖5可以看出，不同監測線上空泡內部低壓區的壓力分布基本一致。壓力的最大峰值P1主要位于流動分離的再附著區空泡末端，且其軸向位置會隨著時間的推移發生改變，這說明空泡的脫落周期內空泡末端會發生小尺度空泡團的脫落。同時不同監測線上最大峰值P1的軸向位置在相同時刻也存在著較大差異，這表明空泡末端小尺度空泡團脫落周向位置具有一定的隨機性。由于反向射流的作用，在空泡內部低壓區會出現壓力次高峰P2，其位置會隨著時間的推移逐漸向上游移動。對比同一時刻不同監測線上的壓力分布可以看出，次高峰P2的軸向位置存在顯著差異，這說明反向射流在沿回轉體壁面向上游運動的過程中周向方向上存在著差異。

(a) t0

(b) t0+0.25T

(c) t0+0.5T

(d) t0+0.75T圖5 不同時刻下監測線上回轉體表面壓力系數分布Fig.5 Instantaneous pressure distribution of different monitoring lines

圖6統計了不同時刻下P2的位置，其中縱坐標表示P2所處位置與回轉體頭部之間的無量綱距離。從圖中可以看出，曲線的斜率表示各監測線上反向射流的運動速度。在t0時刻，C監測線上P2的位置較A，B和D的靠下游。隨著時間的推移，由曲線的斜率可以看出，不同監測線上的反向射流以不同的速率向上游運動，并且同一監測線上不同時刻反向射流的推進速度也不盡相同，其中D監測線上反向射流的平均推進速度最快。結合圖2可以看出，當反向射流與通氣氣流相遇后被推向外側，反向射流開始遠離回轉體壁面朝空泡壁面運動，從而引起空泡發生漲起。由于沿回轉體表面圓周方向上反向射流的推進速度存在差異，空泡外側高壓造成空泡中部邊界出現凹陷，形成局部斷裂；隨著時間的推移，斷裂位置逐漸向兩側擴散，從而U型空泡團脫落。

圖6 次高壓P2位置隨時間的變化Fig.6 Instantaneous positions of P2

3.2 三維非定常空泡斷裂特性分析

為了進一步詳細地分析通氣云狀階段空泡的非定常斷裂特性，本節采用拉格朗日方法進行了研究。FTLE可表示流場中各質點的動態趨勢，而非限于當前的流場，高FTLE值區域會與低FTLE值區域卷曲交錯，進而形成“脊”狀結構，即拉格朗日擬序結構，可精確捕捉流場中的旋渦邊界。由于FTLE主要用于描述兩個相鄰粒子在有限時間T范圍內的平均變動率，其值與選定的參考時間尺度或周期性運動的周期有著直接聯系，在進行相應計算時，采用時間向后積分的方法，積分時間為一個周期T=17 ms。圖7給出了通氣云狀空泡在一個脫落周期內典型時刻下縱截面內的FTLE云圖以及示蹤粒子的分布情況，圖8為相應的含氣量云圖。初始時刻，如圖7(a)所示，在回轉體的上下表面空泡內部的5個特征位置：通氣孔處、通氣孔前、通氣孔后、空泡中部、空泡尾部，呈對稱狀分別布置了5組示蹤粒子，圖中以不同顏色加以區分并依次編號為A，B，C，D和E。在一段時間內，每個示蹤粒子的位置會發生相應的改變，根據粒子的分布情況可以分析得到流場結構的變化，進而可以分析空泡的脫落特性。

觀察圖7中不同時刻下的FTLE云圖可以發現，空泡覆蓋區域的LCS結構整體呈橢球狀分布，與圖8中的空泡外形基本一致；橢球狀內部為不規則的復雜LCS結構，與空泡內部形成的若干小尺度旋渦結構一一對應，這說明LCS可以描繪出通氣云狀空泡內復雜的流動結構。同時由粒子的始末位置可以看出，示蹤粒子基本沿著LCS結構運動，回轉體上下表面初始位置呈對稱分布的粒子經過相同時間的運動后呈現出非對稱性，說明繞回轉體的通氣云狀空泡內部流動結構是非對稱的。

(a) t0

(b) t0+20.6%T

(c) t0+32.4%T

(d) t0+55.9%T

(e) t0+79.4%T

(f) t0+94.1%T圖7 不同時刻縱截面內的FTLE云圖及示蹤粒子的始末位置Fig.7 The contours of FTLE and the distributions of Lagrangian tracers on the longitudinal section

A，B，C這3組粒子主要沿水氣界面向下游運動，并且隨著時間的發展發生一定的混合，同時由原來的團狀分布變化為帶狀分布，如圖7(b)和圖8(a)所示，這說明相同軸向位置不同徑向位置處粒子的運動速度存在顯著差異。另外，回轉體上下表面粒子的運動方式存在顯著差異，回轉體上表面的粒子主要聚集在通氣孔附近且位置基本不隨時間改變，下表面的粒子以斷面形式逐漸向下游移動，粒子經過后的區域表現為透明的氣相區。

D和E兩組粒子的運動軌跡較為復雜。在t0+20.6%T時刻，D組粒子在反向射流與主流的剪切拉伸作用下，靠近回轉體表面的粒子反向朝空泡前緣移動，遠離回轉體表面的粒子隨主流向下游移動，回轉體上表面中反向移動的粒子呈現出整體的反向移動，而下表面反向移動的粒子呈現出先聚集后移動的趨勢。到t0+32.4%T時刻時，如圖7(c)所示，回轉體下表面反向移動的部分D組粒子則進一步聚集并開始偏離回轉體表面，說明此處的反向射流可能遭遇前方氣流被推向外側。隨著時間的推移，下表面的D組粒子在偏離回轉體表面向空泡界面運動的過程中受到空泡界面處C組粒子的影響逐漸轉向空泡內部，進而向下游移動，說明反向射流被推向外側遭遇空泡壁面后再次向內轉向，反向射流與空泡界面的多次接觸極易引起空泡的斷裂脫落。

在空泡的整個脫落周期內，回轉體上表面部分E組粒子沿其壁面反向朝回轉體肩部移動，部分粒子貼壁向下游發展?；剞D體下表面反向移動的E組粒子與上表面運動方式基本一致，也呈現出貼壁移動，但由于D組粒子的影響，E組粒子在t0+55.9%T時刻開始出現偏離回轉體表面，表明反向射流在此處也被推向外側。而回轉體下表面向下游移動的粒子與上表面的運動方式存在很大差異。在t0+20.6%T時刻，如圖7(a)所示，回轉體下表面已經出現部分E組粒子偏離回轉體表面，說明空泡尾部的反向射流有將空泡抬起的趨勢，從而空泡發生卷起并形成小尺度脫落。隨著時間的推移，E組粒子與A，B，C，D這4組粒子發生混合并逐漸向下游移動，直至t0+94.1%T時刻，A，B和部分C，D，E組粒子組成的混合粒子團完全脫離空泡，空泡完成大尺度的脫落過程。

(a) t0+20.6%T

(b) t0+55.9%T

(c) t0+79.4%T

(d) t0+94.1%T圖8 不同時刻縱截面內的氣體含量云圖和示蹤粒子的始末位置細節Fig.8 The contours of air fraction and the distributions of Lagrangian tracers

為了進一步研究其三維特性，圖9給出了繞回轉體不同橫截面上的FTLE云圖和部分示蹤粒子的運動軌跡。其中，x表示橫截面到回轉體頭部的距離，分別用回轉體直徑進行無量綱化處理。在初始時刻t0時，分別在空泡內部沿y=0和z=0的坐標軸線兩側布置A，B，C，D這4組粒子，由于LCS的分布通常與旋渦邊界相對應，可以用來界定具有不同流動結構的區域，因此每組粒子依據所處LCS結構位置的不同又細分為3小組，即A1，B1，C1，D1組粒子位于圓環LCS結構內部，A2，B2，C2，D2組粒子位于圓環LCS結構邊界上，A3，B3，C3，D3組粒子位于圓環LCS結構外部，如圖9所示。

圖9 x/D=1截面處初始時刻FTLE云圖、示蹤粒子位置以及含氣量云圖Fig.9 The contours of FTLE and air fraction, the tracking particle initial positions at x/D=1 section

對比不同截面的LCS結構分布可以看出，與回轉體頭部的距離越近，LCS的分布越簡單清晰。如圖10(a)所示，在靠近回轉體頭部x/D=1的截面上，LCS整體呈現較為規則的圓環狀，從粒子的運動軌跡可知，此截面上的流動主要分為3個區域：1) 圓環LCS結構內部靠近回轉體壁面的區域，如A1，B1，C1，D1組粒子所示，此區域具有明顯的周向流動，方向或呈順時針或呈逆時針或呈螺旋狀；2) 圓環LCS結構外部遠離回轉體壁面的區域，粒子主要位于空泡壁面處，如A3，B3，C3，D3組粒子所示，粒子并沒有明顯的周向運動，此區域的流動方向基本與主流方向(x方向)一致；3) 圓環LCS結構邊界區域，如A2，B2，C2，D2組粒子所示，此區域的流動兼具上述兩種運動，粒子首先呈現一定的周向運動，在周向運動過程中逐漸靠近空泡壁面，到達空泡壁面后沿主流向下游運動。

隨著截面位置向下游移動，如圖10(b)所示，在空泡中部x/D=2截面上，LCS的分布區域逐漸擴大同時變得不規則，從粒子的運動軌跡可以看出，此截面上的流動與x/D=1的截面上流動基本一致，不同之處在于x/D=2截面上第一個流動區域中粒子的運動變得更加復雜且無規律性，同時周向運動的粒子越來越多，說明周向流動的區域逐漸擴大。當流動發展至空泡尾部，如圖10(c)所示，LCS的分布更加雜亂不規則，已經不再呈現為圓環狀，幾乎所有粒子都表現出周向的運動，說明空泡尾部的三維運動特性最為明顯且復雜，也進一步說明在一個空泡脫落準周期內，空泡尾部容易出現脈動誘發多次的小尺度空泡發生脫落。

(a) x/D=1

(b) x/D=2

(c) x/D=3.5圖10 不同橫截面上FTLE云圖和部分示蹤粒子軌跡Fig.10 FTLE plot and convection process of the Lagrangian tracers in different cross sections

4 結論

本文采用CFX軟件對繞軸對稱回轉體通氣云狀空泡流動進行了三維數值模擬，并基于拉格朗日體系的有限時間李雅普諾夫指數(FTLE)、拉格朗日擬序結構(LCS)和粒子追蹤方法對通氣空泡的流動特性進行分析，得到以下結論：

1) 反向射流在周向上推進的不同步性，導致反向射流與主流相互作用后并未造成空泡的完全斷裂，而是呈現出不規則斷裂，直至在空泡尾部形成大尺度U型空泡團的準周期性的脫落，同時空泡尾部始終伴隨著不規則的小尺度空泡團脫落。

2) 縱截面上空泡覆蓋區域的LCS結構整體呈橢球狀分布，與空泡外形相對應，內部為不規律的復雜LCS結構，與空泡內部形成的若干小尺度旋渦結構相對應。粒子始末位置的非對稱性分布反映了繞軸對稱回轉體的通氣云狀空泡內部非對稱流動結構。

3) 靠近回轉體頭部橫截面上LCS整體呈現較為規則的圓環狀，隨著截面位置向下游移動，LCS的分布區域逐漸擴大同時變得不規則。粒子的運動軌跡反映了空泡內部存在明顯的周向流動，徑向方向越靠近回轉體壁面，軸向方向越靠近空泡尾部，周向流動越明顯。