牟 斌, 江 雄, 王建濤

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

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

牟 斌, 江 雄, 王建濤*

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

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

空化；預(yù)處理；隱式求解；數(shù)值模擬；源項

0 引 言

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

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

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

原始變量選?。?/p>

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

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

氣相采用完全氣體狀態(tài)方程：

對于含任意狀態(tài)方程的流體混合物，應(yīng)用Amagat相混合定律計算：

ρ=∑αiρih=∑Yihi

μ=∑αiμiCp=∑YiCpi

其中Yi為組分質(zhì)量分?jǐn)?shù)，αi為組分體積分?jǐn)?shù)?；旌衔锩芏燃敖M分密度通過下式相聯(lián)系：

1.2 預(yù)處理及差分格式

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

Γp=

(9)

變化到一般曲線坐標(biāo)系推導(dǎo)可得基于原始變量的雅克比矩陣為：

AΓ=

(10)

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

Vn=nxu+nyv+nzw

(11)

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

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

修改后的Roe格式適用于全速域。當(dāng)流場中出現(xiàn)局部超聲速，應(yīng)用HCL熵修正[18]技術(shù)保持計算穩(wěn)定。

1.3 隱式時間離散

首先，給出預(yù)處理方程在一般曲線坐標(biāo)系下的表達(dá)式：

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

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

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

注意到：

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

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

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

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

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

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

式中下標(biāo)“inv”、“vis”、“sor”分別代表無黏、黏性、源項貢獻(xiàn)。

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

湍流模型應(yīng)用k-ωSST兩方程模型，在空化流動中，標(biāo)準(zhǔn)湍流模型過高預(yù)測空泡尾部湍流黏性，抑制空泡脫落，在流動非定常效應(yīng)較強(qiáng)，湍流黏性系數(shù)還需要加入Coutier-Delgosha空化修正：

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

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

2 計算結(jié)果

2.1 NACA0015水翼空化計算

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

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

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

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

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

2.2 錐柱體算例

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

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

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

3 結(jié) 論

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

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

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

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

[1]Chen Xin. An investigation of the ventilated cavitating flow[D]. Shanghai: Shanghai Jiaotong University, 2006. (in Chinese).陳鑫. 通氣空泡流研究[D]. 上海: 上海交通大學(xué), 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).劉樺, 李家春, 何友聲, 等. 十一五水動力學(xué)發(fā)展規(guī)劃的建議[J]. 力學(xué)進(jìn)展, 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).陳瑛. 自然空泡流數(shù)值模擬方法研究[D]. 上海: 上海交通大學(xué), 2009.

[12]Guo Jianhong. Study of the numerical simulation method for complex multiphase cavitating flow[D]. Shanghai: Shanghai Jiaotong University, 2013. (in Chinese).郭建紅. 復(fù)雜多相空泡流的數(shù)值模擬方法研究[D]. 上海: 上海交通大學(xué), 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).周景軍. 通氣超空泡流動及航行體流體動力數(shù)值模擬研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 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 格式的低速預(yù)處理方法研究[J]. 空氣動力學(xué)學(xué)報, 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]. 空氣動力學(xué)學(xué)報, 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