陳 鑫 ， 魯傳敬 ， 曹嘉怡 ， 李 杰 ， 郭建紅

（上海交通大學a.工程力學系；b．海洋工程國家重點實驗室；c.水動力學教育部重點實驗室，上海 200240）

1 引 言

雖然高速水下航行體通常運動在廣闊的水域中，但模型試驗卻是在有邊壁或自由面的設備中進行的，鄰近的固壁邊界對于空泡形態(tài)和水動力可能產(chǎn)生顯著的影響，即所謂洞壁效應。對閉式水洞，必須考慮阻塞問題。流動的限制促使局部流速增大，并形成比較大的空泡。在某個大于0的空化數(shù)下，流動可能被阻塞，從而無法得到更低的空化數(shù)。

關(guān)于空泡流洞壁效應的問題，以往的研究多集中在具有簡單幾何外形的自然空泡流問題上[1-2]，如平板、楔形體和圓錐等。Birkhoff等[3]以無升力楔形體的阻塞空泡流為研究對象，理論分析結(jié)果表明傳統(tǒng)試驗中定義的阻力系數(shù)對洞壁距離的影響十分敏感。Wu[4]采用不同的閉合模型，對閉式水洞中二維空泡流的外形和阻力系數(shù)進行了理論求解，并加以修正以適用相應的無界流動。

上述對洞壁效應問題的研究工作都是基于傳統(tǒng)的勢流理論，對空化流場內(nèi)部流場的精細結(jié)構(gòu)以及初生、發(fā)展、潰滅和脫落下瀉等過程，這類方法缺乏有力的手段。另一方面，勢流理論框架內(nèi)無法考慮粘性，而粘性對空泡脫體點、空泡閉合區(qū)等有直接影響。

本文以一個帶圓盤空化器的水下航行體模型為研究對象，基于均質(zhì)平衡多相流理論和輸運方程類空化模型，通過求解混合介質(zhì)的RANS方程、RNG k-ε湍流輸運方程和各相的質(zhì)量輸運方程，數(shù)值模擬了圓形截面水洞中的定常通氣空泡流動，研究了閉式空泡水洞中洞壁效應對通氣空泡流動中的空化數(shù)及流場壓力分布的影響，并擬合得到了一定空化數(shù)范圍內(nèi)計算空泡尺寸和模型阻力系數(shù)的近似公式。

2 基本方程和數(shù)值方法

本文基于均質(zhì)平衡多相流理論和輸運方程類空化模型[5-7]，研究通氣空泡流動。將由空氣、水蒸汽和水組成的混合介質(zhì)看成一種變密度單流體，各相共享同一壓力、速度場，忽略相間滑移速度和重力效應。通過引入空氣、水蒸汽和水的體積分量—αg、αv、αl，得到描述氣、汽和液多相流動的控制方程。

連續(xù)性方程：

動量方程：

能量方程：

蒸汽相連續(xù)性方程：氣相連續(xù)性方程：

考慮氣體為可壓縮理想氣體，補充一個關(guān)聯(lián)p、ρg和T的狀態(tài)方程：

方程中混合介質(zhì)的各相應滿足相容性條件：

混合介質(zhì)密度ρm和粘度μm由體積分量加權(quán)平均獲得：

此外，采用兩個獨立的輸運方程描述（4）式中汽、水間的質(zhì)量傳遞過程：

式中ui和xi分別指代速度分量和坐標方向；t為時間；p表示壓力；μt表示湍流粘度；T為溫度；kt是流體傳熱系數(shù)；cp為定壓比熱；ST為流體的內(nèi)熱源及由于粘性作用流體機械能轉(zhuǎn)換為熱能的部分；R0是氣體常數(shù)；m˙-模擬的是水到水蒸汽的轉(zhuǎn)換過程，這個量與液相體積分量以及當?shù)貕毫惋柡驼羝麎篜v之間的差值成正比，該模型與Merkle等[8]使用的模型相同；對于水蒸汽到水相變過程中的質(zhì)量傳遞率m˙+，采用了 Ginzburg-Landau 模型[9]；V∞為來流速度；Cevap、Ccond為經(jīng)驗常數(shù)，t∞是平均流時間尺度；下標 m、l、g和v分別指代混合介質(zhì)、液相、氣相和蒸汽相。

另外，采用標準壁面函數(shù)[10]的兩方程RNG k-ε湍流模式[11]，使得上述粘性方程組封閉。

3 數(shù)值方法

本文應用有限體積法離散積分微分型控制方程。壓力梯度項采用PEOSTO！格式[12]離散；動量方程的差分格式選用二階逆風格式；湍流輸運方程的差分格式選用一階逆風格式；壓力—速度耦合采用SIMPLE算法[12]。離散化后的代數(shù)方程系統(tǒng)的求解采用基于分離算法的Gauss-Seidel線性方程求解器，并結(jié)合“代數(shù)”多重網(wǎng)格法加速收斂。文中數(shù)學模型的建立及求解采用通用CFD軟件—FLUENT實現(xiàn)。

水下航行體模型由圓盤空化器、雙碗通氣裝置和軸對稱后體組成，如圖1（a）所示?？傞LL0=1 000mm，最大直徑D=75mm，圓盤空化器直徑Dn=0.2D。流動具有軸對稱性，故簡化為二維軸對稱問題求解。圖1（b）顯示的是圓盤直徑與水洞直徑之比Dn:Dt為1:40情況下的計算域。取模型對稱軸為x軸，坐標原點位于圓盤空化器頂面圓心處。入口距航行體圓盤0.513L0，出口距航行體末端L0，水洞壁距對稱軸4D。網(wǎng)格劃分采用多區(qū)塊四邊形網(wǎng)格，航行體頭體部分結(jié)構(gòu)較為復雜，圖1（c）給出模型頭部附近區(qū)塊的網(wǎng)格劃分情況；其余區(qū)塊均屬十分規(guī)則的四邊形網(wǎng)格，不再一一圖示。

邊界條件設置如下：入口邊界給定軸向速度10m/s、徑向速度為0；出口邊界壓力給定24 867.8Pa；通氣孔處為壓力入口邊界，給定總壓、總溫和靜壓；水洞壁和模型表面粗糙度為0，給定無滑移壁面邊界條件。

圖1 水下航行體模型、計算域及網(wǎng)格劃分（Dn:Dt=1:40）：（a）模型對稱面；（b）計算域與邊界條件；（c）頭體附近的網(wǎng)格Fig.1 The model,computational domain and grid(Dn:Dt=1:40):(a)The model section;(b)The computational domain and boundary conditions;(c)The grid near the head body

計算中，對于Dn:Dt分別為1:25、1:30、1:40、1:50和1:80的情況，模擬了不同空化數(shù)下的定常通氣空泡流。

4 結(jié)果與討論

由于通氣量大小不同，通氣空泡形態(tài)和壓力分布情況也有所不同。小通氣量下，通氣空泡表現(xiàn)為局部空泡形態(tài)，空泡尾部存在明顯的壓力回復（圖2（a））；而大通氣量下，由于泡內(nèi)壓力增大通氣空泡能夠跨越模型錐段從而形成超空泡，且泡內(nèi)壓力大致相同，壓力回復出現(xiàn)在模型后（圖2（b））。

圖2 典型通氣空泡形態(tài)與壓力等值線圖（Dn:Dt=1:40）Fig.2 The shape of typical ventilated cavity and pressure contour(Dn:Dt=1:40)

圖4中實線表示由計算得到的臨界空化數(shù)的線性擬合曲線；而虛線為文獻[1]中三維圓盤自然空泡流中阻塞空化數(shù)的勢流近似解。當Dn:Dt較小時兩者比較接近，而隨著Dn:Dt的增大，兩者的差別逐漸增大。差別產(chǎn)生的原因在于：一方面近似解沒有考慮粘性的影響，邊界層加強了流動的阻塞，另一方面由于通氣及后體的影響，也使得本文阻塞空化數(shù)總是大于勢流近似解。

圖 3 σc與 M之間的關(guān)系Fig.3 The relationship between σcand M

試驗中壓力的測量通常是通過安置在水洞壁面上的壓力傳感器實現(xiàn)，本文采用數(shù)值模擬的方法分析了水洞壁面以及對稱軸上的壓力變化規(guī)律，有助于了解通氣空泡流動情況下整個水洞中的壓力場分布特點，為試驗測量提供技術(shù)參考。

圖5～7是圓盤與水洞的直徑比Dn:Dt分別為 1:30、1:40 和 1:50 時，無物體流動（水洞內(nèi)未放置模型）以及一組通氣空泡流動（水洞內(nèi)放置水下航行體模型）情況下，水洞主體段沿長度方向?qū)ΨQ軸及壁面上的壓力分布曲線，無物體流動算例中入口、出口和洞壁邊界條件的設置與通氣空泡流動算例相同。

圖 4 σb與 Dn:Dt之間的關(guān)系，來自文獻[1]Fig.4 The relationship between σband Dn:Dtfrom Ref.[1]

整體上看，除模型前后一定范圍內(nèi)，無論何種計算條件下，由于粘性帶來的沿程損失，軸線和壁面上的壓力沿流動方向基本呈下降的趨勢。無物體流動情況下，水洞軸線和壁面上的壓力按線性規(guī)律減??；通氣空泡流動情況下，由于模型的置入和通入氣體帶來的能量變化，壓力分布規(guī)律較為復雜。

圖5 當Dn:Dt=1:30時，水洞對稱軸和壁面上的壓力分布Fig.5 Pressure distributions on the axis and wall surface of tunnel at Dn:Dt=1:30

圖6 當Dn:Dt=1:40時，水洞對稱軸和壁面上的壓力分布Fig.6 Pressure distributions on the axis and wall surface of tunnel at Dn:Dt=1:40

圖7 當Dn:Dt=1:50時，水洞對稱軸和壁面上的壓力分布Fig.7 Pressure distributions on the axis and wall surface of tunnel at Dn:Dt=1:50

比較圖5～7中水洞主體段起點到終點之間的壓力差值，可以看出隨著水洞直徑的減小，各種流動情況下沿程的壓力損失增大。由于計算中在出口邊界給定常壓，這使得水洞入口處壓力升高。因此，按照空化數(shù)的定義，如將水洞入口測量的壓力作為來流壓力P∞，則水洞直徑的減小將導致空化數(shù)增大。此外，由于水洞中的沿程壓力損失還與水洞壁面的粗糙度（本文中計算取為0）有關(guān)。由于水洞壁面的粗糙會引起沿程損失的增大，與光滑壁面的水洞相比得到的空化數(shù)偏大，因此在水洞試驗中應盡量保持水洞壁面光滑。

同一Dn:Dt下的通氣空化流動中，圓盤空化器前由于流動滯止壓力顯著升高，而模型后端類似繞臺階流動，也出現(xiàn)明顯的壓力峰值。距圓盤約20Dn之前及距模型末端約40Dn之后的范圍內(nèi)，與水洞中未放置模型的流動相似，由于沿程損失水洞壁面和對稱軸上的壓力分布滿足線性減小的變化規(guī)律。但是，在通氣空化發(fā)生的區(qū)域內(nèi)，由于流道變窄，流速增大，使得水洞壁面上的壓力降低。其壓力也不再滿足線性規(guī)律，而是按照類似于二次曲線的規(guī)律變化。有鑒于此，在與本文計算條件相似的水洞模型試驗中，應將來流壓力的測壓點置于距圓盤空化器至少約20Dn之前，可有效地減小模型本身對測量結(jié)果的影響。

同一水洞中，隨著通氣空化數(shù)的減小，水洞內(nèi)壓力升高的幅度逐漸增大。這是因為通氣空化數(shù)的減小對應于通氣量的增加和空泡的拉長變粗，意味著阻塞效應增強；同時通入氣體的總動能增加，從能量轉(zhuǎn)換的角度來看就表現(xiàn)為壓力的升高。此外，由于空泡影響區(qū)域的擴大，水洞壁面上低壓區(qū)域隨著空化數(shù)的減小而增大。同時，模型后水洞軸線上的壓力峰值則隨著空化數(shù)的減小而減小。

根據(jù)計算結(jié)果發(fā)現(xiàn)不同圓盤水洞直徑比下軸對稱通氣空泡的長度、最大直徑以及模型阻力系數(shù)與空化數(shù)之間的關(guān)系曲線均具有相似性，表明它們與通氣空化數(shù)以及圓盤水洞直徑比這兩個特征參數(shù)之間存在一定的函數(shù)關(guān)系。因此，本文利用數(shù)據(jù)處理工具Origin，通過非線性擬合和坐標平移等方法，給出了它們之間的近似關(guān)系式，適用于一定的通氣空化數(shù)和圓盤水洞直徑比范圍。

圖8 通氣空泡長度與空化數(shù)的歸一化曲線Fig.8 Normalized correlation curve between ventilated cavity length and cavitation number

圖9 通氣空泡最大直徑與空化數(shù)的歸一化曲線Fig.9 Normalized correlation curve between maximum ventilated cavity diameter and cavitation number

圖8～10為擬合后的歸一化關(guān)系曲線（以實線表示）與計算數(shù)據(jù)（以離散的數(shù)據(jù)點表示）的對比結(jié)果。

圖10 模型阻力系數(shù)與通氣空化數(shù)的歸一化曲線Fig.10 Normalized correlation curve between drag coefficient of model and ventialted cavitation number

5 結(jié) 論

本文基于均質(zhì)平衡多相流理論和輸運方程類空化模型，通過求解混合介質(zhì)的RANS方程、RNG k-ε湍流輸運方程和各相的質(zhì)量輸運方程，數(shù)值模擬了閉式水洞中帶圓盤空化器航行體模型的定常通氣空泡流動，分析了閉式水洞中洞壁效應對通氣空化數(shù)、壓力分布規(guī)律、空泡尺寸以及模型阻力系數(shù)的影響。在本文的計算條件下，可以得到以下結(jié)論：

（1）對于帶圓盤空化器的水下航行體模型的通氣空泡流動，計算得到的阻塞空化數(shù)線性正比于圓盤水洞直徑比，且與三維圓盤自然空泡流的勢流近似解基本一致，但由于粘性邊界層、通氣及后體的影響而略大于后者；

（2）在距航行體模型圓盤約20Dn之前以及距航行體模型末端約40Dn之后的區(qū)域內(nèi)，由于沿程損失，水洞軸線和壁面上的壓力分布滿足線性規(guī)律；而在通氣空化區(qū)域內(nèi)，水洞壁面上的壓力分布按照類似于二次曲線的規(guī)律變化；

（3）同一水洞中，隨著通氣空化數(shù)的減小，水洞內(nèi)壓力升高的幅度逐漸增大，水洞壁面上低壓區(qū)域增大，而模型后水洞軸線上的壓力峰值減??；

（4）根據(jù)計算數(shù)據(jù)，擬合得到了空泡長度、最大直徑以及模型阻力系數(shù)的近似公式，適用的通氣空化數(shù)與圓盤水洞直徑比范圍為：0.02＜σc＜0.12、0.01＜Dn:Dt＜0.04。

[1]Cohen H,Difrima R C.Wall effects in cavitating flows[C]//Second Symposium on Naval Hydrodynamics.Washington D.C.,U.S.,1958:367-390.

[2]Choi J K,Kinnas S A.Numerical water tunnel in two-and three-dimensions[J].Journal of Ship Research,1998,42:86-98.

[3]Birkhoff G,Plesset M,Simmons N.Wall effects in cavity flow II[J].Quart.Appl.Math,1952,9:413.

[4]Wu Y T,Whitney A K,Brennen C.Cavity-flow wall effects and correction rules[J].J Fluid Mech.,1971,49:223-256.

[5]Chen X,Lu C J,Wu L.An multiphase model on simulating a ventilated cavitating flow[J].Journal of Hydrodynamics,Ser.A,2005,20(Suppl.):916-920.(in Chinese)

[6]Zhang L X,Zhao W G,Shao X M.A pressure-based algorithm for cavitating flow computations[J].Journal of Hydrodynamics,Ser.B,2011,23(1):42-47.

[7]Chen Y,Lu C J.A homogenous-equilibrium-model based numerical code for cavitation flows and evaluation by computation cases[J].Journal of Hydrodynamics,Ser.B,2008,20:186-194.

[8]Merkle C L,Feng J,Buelow P E O.Computational modeling of the dynamics of sheet cavitation[C]//Third International Symposium on Cavitation.Grenoble,France,1998.

[9]Hohenberg P C,Halperin B I.Theory of dynamic critical phenomena[J].Rev.Mod.Phys.,1977,49(3):435-479.

[10]Launder B E,Spalding D B.The numerical computation of turbulent flows[J].Computer Methods in Applied Mechanics and Engineering,1974,3:269-289.

[11]Yakhot V,Orszag S A.Renormalization group analysis of turbulence I:basic theory[J].Journal of Scientific Computing,1986,1(1):1-51.

[12]Patankar S V.Numerical heat transfer and fluid flow[M].Hemisphere,Washington D.C.,1980.