空化流動隱式求解方法研究

牟 斌, 江 雄, 王建濤

(中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000)

空化流動隱式求解方法研究

牟 斌, 江 雄, 王建濤*

(中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000)

空化流動問題本質上是可壓縮流動，應用可壓縮方法開展數值模擬研究更符合物理實際。氣體和液體壓縮性的顯著差別使得低速空化流動數值模擬的剛性問題非常突出，通過引入預處理技術解決低速問題中由于特征值量值不一致導致的收斂剛性問題，提高收斂速度。同時，以預處理后的特征值構造Roe格式耗散項，提高低速流動計算精度。鑒于自然空化流動中氣液組分轉換現象和物質輸運現象并存，且氣體和液體的密度在常溫狀態下存在顯著差別，預處理后的源項雅克比矩陣的特征值與無粘通量雅克比矩陣的特征值存在量級的差異，這會使得求解過程不穩定或收斂速度極慢，即出現“源項剛性”問題。為此，本文系統推導了預處理框架下的氣液兩相流隱式求解方法，采用點隱式方法處理源項，通過直接求逆的方式增強算法的穩定性。研究中分別考察了三種不同的算子分裂方式在LU-SGS(Lower-Upper Symmetric-Gauss-Seidel)隱式迭代中的模擬效果，并提出適用于DDADI(Diagonally Dominant Alternating Direction Implicit)的算子分裂方式，在NACA0015水翼空化算例模擬中綜合比較了上述四種隱式迭代方式對低速空化問題收斂性的影響。最后，應用本方法對三維尖錐自然空化算例進行了考核模擬，捕捉到自然空泡流場中的主要特征，計算結果與試驗測量結果吻合。

空化；預處理；隱式求解；數值模擬；源項

0 引 言

隨著俄羅斯“暴風雪”超空泡魚雷面世，發展以潛射導彈、高速魚雷為代表的新型高速與超高速水中兵器，成為世界先進大國關注的重大課題。俄羅斯第二代速度達到200 m/s的魚雷正在研制過程中，美、德、英、法等國也都在進行超空泡減阻技術的基礎與應用研究，我國也于近年開展了超空泡技術的基礎研究[1-2]。

超空泡流動現象涉及到了多相流、湍流、相變、可壓縮性和非定常等復雜流動機制，迄今對于這一復雜流動的了解還十分有限。早期研究主要以理論研究為主，輔助實驗結果修正理論公式中的系數而得到一些經驗公式。隨著計算機硬件性能提升及現代CFD技術的飛速發展，超空化問題的研究已經發展出現代的基于Navier-Stokes方程的數值模擬技術。美國賓州大學學者Kunz、Lindau等，基于預處理技術，應用均質平衡流假設發展了空化流動模擬軟件—UNCLE_M，可計算自然空化、通氣空化。同時軟件耦合了六自由度方程，可以計算航行體完整的彈道和姿態，在高速水下超空泡航行體流體動力的數值模擬方面發表了不少有影響的文章，代表了目前計算研究的方向[3-10]。國內上海交大陳鑫等基于SIMPLE算法，開發了可用于模擬自然、通氣空化的軟件[1, 11-13]。SIMPLE算法中組分輸運方程與流場方程分離求解，對組分輸運方程進行松弛即可保證穩定求解，但SIMPLE方法本身由不可壓方法發展而來，在求解高速問題時需進行修正。

本文以可壓縮預處理技術為基礎，應用均質平衡流假設，發展了基于輸運方程空泡模型的空化模擬代碼。在基于輸運方程空化模型中，液體的蒸發和水蒸氣的凝結過程通過輸運方程模擬實現，源項控制相間的相互轉化。常溫條件下，水蒸汽和液態水的密度相差1000倍，源項雅克比矩陣預處理后特征值與無黏通量雅克比矩陣特征值存在量級差異，常規對角化及簡化譜半徑方法在大的源項參數下難以獲得收斂結果，出現“源項剛性”問題。在Kunz系列文章中，隱式化處理水的破壞項，水的生成項則采用顯式松弛方式。本文采用點隱式方法處理整個源項，在近似因子分裂基礎上，分別應用了三種矩陣分裂技術，考察源項矩陣處理的影響，在LUSGS、AFADI方法上均實現了穩定計算。

1 數值模擬方法

1.1 控制方程

在均勻一致、運動學、熱力學平衡假設下，水、汽、不可凝結氣體的兩相流動可以描述為一種混合流體的單相流動。加入預處理的非定常無量綱控制方程如下[14]：

上式中應用了“雙時間”方法，τ為虛擬時間變量，t為物理時間變量。守恒變量定義為：

原始變量選取：

應用原始變量有利于推導公式。E、F、G為無黏通量[14]，Ev、Fv、Gv為黏性通量，源項S為：

式中psat為水的飽和蒸汽壓。水的狀態方程選用加入溫度修正的Tait方程[15]：

氣相采用完全氣體狀態方程：

對于含任意狀態方程的流體混合物，應用Amagat相混合定律計算：

ρ=∑αiρih=∑Yihi

μ=∑αiμiCp=∑YiCpi

其中Yi為組分質量分數，αi為組分體積分數。混合物密度及組分密度通過下式相聯系：

1.2 預處理及差分格式

由于空化計算涉及的水的流動速度一般在10 m/s左右，流動馬赫數為10-3量級，傳統可壓縮求解方法會遇到雅克比矩陣特征值相差量級導致的剛性問題，必須應用預處理技術[16]。通過在控制方程時間導數項前乘以預處理矩陣改變雅克比矩陣特征值，使所有特征值處在相同量級，達到消除方程剛性的目的。考慮多相的預處理矩陣Weiss-Smith[6]矩陣為：

Γp=

(9)

變化到一般曲線坐標系推導可得基于原始變量的雅克比矩陣為：

AΓ=

(10)

雅克比矩陣其特征值及其它參數如下：

Vn=nxu+nyv+nzw

(11)

上式中，b=1時，非定常預處理轉化到定常預處理，而β=1時，方程回到無預處理情況。在求解非定常問題時，右端項需要定常雅克比矩陣及特征值[17]，方程求解時左端項和右端項特征值不匹配問題通過限制局部時間步長來解決。本文應用Roe's FDS格式空間離散，在預處理情況下，Roe格式需根據預處理后的特征向量進行重構，標準的Roe格式如下：

第一項為通量項，不作修改，第二項為Roe格式的耗散項，以新的特征值作為耗散的尺度。

修改后的Roe格式適用于全速域。當流場中出現局部超聲速，應用HCL熵修正[18]技術保持計算穩定。

1.3 隱式時間離散

首先，給出預處理方程在一般曲線坐標系下的表達式：

因此，矩陣Ap的分裂可以寫為：

為便于闡述，將式(14)在一維情形下矩陣分裂得到：

空化計算中對源項矩陣T的處理至關重要，由于液態水的密度為水蒸氣密度的1000倍以上，T的特征值往往比無黏項特征值大幾個量級，采用類似黏性譜半徑及矩陣對角化的方法會導致收斂極其緩慢，本文對其直接求逆。上式可以改寫為：

注意到：

將式(20)寫成LUSGS分裂格式：

式(20)的求解每一步迭代需要對D矩陣求逆2次。將式(19)中的源項矩陣T寫到L掃描中，可以得到如下近似因式：

式(24)與式(20)相比，L掃描完全一致，U掃描中無源項矩陣T。這樣做的好處是減少了矩陣求逆次數，可以大大節省計算時間。同時，也可將T矩陣移到U掃描，過程與移到L掃描類似。

對式(14)的離散還可采用DDADI方法，令：S=SpV/Δt，近似因子分裂得到：

這樣分裂后，三個方向的矩陣算子均可以對角化，將無法對角化的源項矩陣移到最后，與單相流動求解相比，僅需在三個方向掃描后增加一次矩陣求逆。

在含相變問題求解中，限制相變劇烈區時間步長可以獲得收斂過程平緩的結果。與文獻[19]不同，本文采用下式限制空化區時間步長：

式中下標“inv”、“vis”、“sor”分別代表無黏、黏性、源項貢獻。

1.4 湍流模型及初邊值界條件

湍流模型應用k-ωSST兩方程模型，在空化流動中，標準湍流模型過高預測空泡尾部湍流黏性，抑制空泡脫落，在流動非定常效應較強，湍流黏性系數還需要加入Coutier-Delgosha空化修正：

在不涉及底部分離算例中，初場可以選擇為來流值；在計算大分離算例時，以全濕流的計算結果為初場可以穩定計算。

邊界條件主要有遠場邊界、固壁邊界、對稱邊界、奇性軸邊界等，與單相流類似處理。

2 計算結果

2.1 NACA0015水翼空化計算

NACA0015水翼算例為國外發布的考核空化計算軟件的標準算例，計算參數：p∞=0.12×105(空化)，V∞=3.41 m/s，T∞=300 K，psat=3752 Pa，迎角已經預偏4°。

計算中Cdest=104，Cprod=100，隱式格式分別采用1.4節中的四種方法，記式(19)為LUSGS_SD，式(23)為LUSGS_SL，式(24)為LUSGS_SU，式(25)為DDADI。計算的殘差和收斂曲線及升力系數收斂曲線見圖2、圖3。水翼壁面為黏性固壁，水洞壁為無黏固壁，入口和出口指定為遠場，以來流初始化流場。

圖2的殘差收斂曲線表明，采用的四種隱式算法計算在空化情況下連續方程殘差L2模可以收斂至10-10以下。LUSGS的三種分裂方式收斂曲線在迭代4000步以后的形態完全一致；LUSGS_SD對D矩陣求逆2次，耗時最多，但在流場“暫態”期穩定性最好；LUSGS_SU、LUSGS_SL分別僅在LUSGS的U、L掃描中考慮源項作用，如果局部時間步長不采用源項限制，在計算初期殘差振蕩，甚至導致計算發散。

圖3為升力系數曲線對比，可以看到，LUSGS_SD和LUSGS_SL幾乎完全一致，而LUSGS_SU和前二者差別較大是由于局部時間步長加了限制。DDADI則由于穩定性較差，計算中加入了更多其他穩定性措施，導致收斂較慢，還有改進空間。四條曲線在迭代步數八千步后完全重合，與隱式方法不影響收斂結果的論斷相符。

從圖4的壓力系數等值線流場可以看到，隨著流動從前緣向后發展，壓力逐漸降低，當壓力低于飽和蒸汽壓時，發生自然空化；空泡隨著流動向后逐漸擴大，泡內壓力保持為常數；空泡在翼型中部閉合，閉合區引起回射流，導致出現較大的壓力梯度，與超臨界翼型翼面發生激波的圖像類似。

2.2 錐柱體算例

Rouse和McNown研究了一系列典型回轉體的空化現象，并發布了實驗結果。本文選擇了外形簡單的22.5°錐柱體算例對本文發展的方法進行進一步驗證。算例參數：V∞=4.3 m/s，T∞=300 K，psat=3589 Pa，空化數σ=0.3、0.4、0.5，全濕流狀態。

計算采用三維計算，網格取C型網格，網格維數113×65×17(流向×法向×周向)，壁面最小距離Δn=4×10-4，壁面第一層網格y+≈1-3。Cdest=104，Cprod=100。

圖5為σ=0.3計算結果，流動在過尖錐后流動發生空化，空化區壓力系數近似為常數，空化區形成大的分離區。從圖6壓力系數比較曲線看到，隨著空化數降低，空泡長度增大，回射流現象更劇烈。本文計算結果與試驗數據基本一致,較為準確地捕捉了空泡的起始以及閉合區域位置。

3 結 論

本文在可壓縮方法框架下，通過應用低速預處理技術發展了基于混合模型的多相、多組分的數值計算方法，對NACA0015水翼與錐柱體標模的驗證計算表明：

1) 源項處理方式合理，LUSGS的三種隱式分裂及DDADI方法均能獲得穩定收斂解，殘差可以降到機器零；

2) 計算結果與試驗值基本一致，所采用的算法能準確地捕捉到空泡的起始及潰滅。

下一階段將對通氣空化進行驗證計算研究，同時考慮加入六自由度計算模塊模擬出水現象。

[1]Chen Xin. An investigation of the ventilated cavitating flow[D]. Shanghai: Shanghai Jiaotong University, 2006. (in Chinese).陳鑫. 通氣空泡流研究[D]. 上海: 上海交通大學, 2006.

[2]Liu Hua, Li Jiachun, He Yousheng, et al. Suggestion on the research frame programme on hydrodynamics for the eleventh five year plan[J]. Advances in Mechanics, 2007, 37(1):142-147. (in Chinese).劉樺, 李家春, 何友聲, 等. 十一五水動力學發展規劃的建議[J]. 力學進展, 2007, 37(1): 142-147.

[3]Kunz R F, Boger D A. A preconditioned navier-stokes method for two-phase flows with application to cavitation prediction[R]. AIAA, 1999-3329, 1999.

[4]Venkateswaran S, Lindau J W, Kunz R F, et al.Computation of multiphase mixture flows with compressibility effects[J]. Journal of Computational Physics, 2002, 180: 54-77.

[5]Lindau J W, Sankaran V, Kunz R F, et al. Development of a fully-compressible multiphase Reynolds-averged Naviar-Stokes model[R]. AIAA CFD Conference, AIAA 2001-2648, Anaheim, CA, 2001.

[6]Kunz R F, Boger D A, Stinebring D R, et al. A preconditioned Navier-Stokes method for two-phase flows with application to cavitation predication[J]. Computers and Fluids, 2000, 29: 849-875.

[7]Li D, Venkateswaran S, Lindau J, et al. A unified computation formulation for multi-component and multi-phase flows[R]. AIAA-2005-1391, Reno, NV, January 10-13, 2005.

[8]Kinzel M P. Computational techniques and analyse of cavitating

fluid flows[D]. The Pennsylvania State: University the Graduate School Department of Aerospace Engineering, 2008.

[9]Kinzel M P, Willits S M, Lindau J W, et al. CFDSimulations of oscillating hydrofoils with cavitation[C]// The 44th Aerospace Sciences Meeting and Exhibit. Reno, Nevada, 2006.

[10]Kunz R F. Multi-phase CFD analysis of natural and ventilated cavitation about submerged bodies[R]. FEDSM99-7364, 1999.

[11]Chen Ying. Study of the numerical method for natural cavitating flow[D]. Shanghai: Shanghai Jiaotong University, 2009. (in Chinese).陳瑛. 自然空泡流數值模擬方法研究[D]. 上海: 上海交通大學, 2009.

[12]Guo Jianhong. Study of the numerical simulation method for complex multiphase cavitating flow[D]. Shanghai: Shanghai Jiaotong University, 2013. (in Chinese).郭建紅. 復雜多相空泡流的數值模擬方法研究[D]. 上海: 上海交通大學, 2009.

[13]Zhou Jingjun. Numerical simulation study on the ventilated supercavitating flow and hydrodynamics of vehicle[D]. Harbin: Harbin Institute of Technology, 2011. (in Chinese).周景軍. 通氣超空泡流動及航行體流體動力數值模擬研究[D]. 哈爾濱: 哈爾濱工業大學, 2013.

[14]Tian Ming. Numerical simulation of the internal two-phase flow within an aerated-liquid injector and its injection into the corresponding high-speed cross flow[D]. North Carolina State University, 2005.

[15]Koop Arjen, Hoeijmakers Harry. Numerical simulation of unsteady three-dimensional sheet cavitation[C]//The 7th International Symposium on Cavitation. Ann Arbor, Michigan, USA, 2009.

[16]Mou B, Xiao Z Y, Liu G, et al. Low speed preconditioning for Roe scheme[J]. Acta Aerodynamica Sinica, 2010, 28(2): 125-131. (in Chinese).牟斌, 肖中云, 劉剛, 等. 基于 Roe 格式的低速預處理方法研究[J]. 空氣動力學學報, 2010, 28(2): 125-131 .

[17]Buelow P E O, Schwer D A, Feng Jinzhang, et al. A preconditioned dual-time, diagonalized ADI scheme for unsteady computations[R]. AIAA-97-2101, 1997.

[18]Lin Hong-Chia. Dissipation additions to flux-difference Splitting[J]. Journal of Computational Physics, 1995, 117: 20-27.

[19]Gerlinger P. An implicit multigrid method for the simulationof chemically reacting flows[J]. Journal of Computational Physics, 1998, 146: 322-345.

Investigation on implicit numerical method for cavitation flow

Mou Bin, Jiang Xiong, Wang Jiantao*

(ComputationalAerodynamicsInstituteofChinesAerodynamicDevelopmentandResearchCenter,Mianyang621000,China)

The cavitation flow is essentially compressible, so the compressible method is more appropriate for the physical essentials. The difference of gas and liquid in compressibility intensively deteriorates the stiffness issue in low speed cavitation flow simulation. A compressible solver is adopted to numerically investigate this kind of two-phase flow with precondition technique, to shrink the convergence course by dealing with the issue of imbalance of eigenvalues in low speed flow. Meanwhile, the dissipative term of Roe’s scheme is rebuild on the basis of preconditioned system to improve the accuracy under low speed flow. Moreover, the phase transition and the material convection coexist in natural cavitation flow, and the fluids of gas and liquid vary much in density at normal temperature, these two factors are consequently followed by that the eigenvalues of source Jacobian matrix and inviscid flux Jacobian matrix differs in order of magnitude. This phenomenon is the “source stiffeness” problem and results in the unstability of the solver. The Point implicit method is applied to treat the source item, and the stability of method is enhanced by directly inversing the matrix. Three different operator splitting schemes are investigated in LU-SGS(Lower-Upper Symmetric-Gauss-Seidel) implicit iteration, and an operator splitting scheme suited for DDADI(Diagonally Dominant Alternating Direction Implicit) is developed. A comparison between the four implicit schemes is drew in NACA0015 hydrofoil low-speed cavitation case, to study the influence of different schemes on convergence. At the last part, the natural cavitation case of 3-D(three- dimension) cone is simulated by the current method, capturing the main characteristic of the natural cavity flow field, and attaining a group of simulating data that matches well with the experiment data.

cavitation; precondition; implicit method; numerical simulation; source item

0258-1825(2017)01-0027-06

2015-05-25;

2015-07-31

牟斌(1974-),男,四川什邡人,研究員,博士,研究方向:水汽兩相流. E-mail:809970229@qq.com

王建濤*(1982-),研究方向:水汽兩相流. E-mail: jtwang@ustc.edu

牟斌, 江雄, 王建濤. 空化流動隱式求解方法研究[J]. 空氣動力學學報, 2017, 35(1): 27-32.

10.7638/kqdlxxb-2015.0027 Mou B, Jiang X, Wang J T. Investigation on implicit numerical method for cavitation flow[J]. Acta Aerodynamica Sinica, 2017, 35(1): 27-32.

V211.3

A doi: 10.7638/kqdlxxb-2015.0027