牟 斌, 江 雄, 王建濤

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

空化流動(dòng)隱式求解方法研究

牟 斌, 江 雄, 王建濤*

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

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

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

0 引 言

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

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

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

原始變量選取：

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

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

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

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

ρ=∑αiρih=∑Yihi

μ=∑αiμiCp=∑YiCpi

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

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

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

Γp=

(9)

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

AΓ=

(10)

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

Vn=nxu+nyv+nzw

(11)

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

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

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

1.3 隱式時(shí)間離散

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

因此，矩陣Ap的分裂可以寫(xiě)為：

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

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

注意到：

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

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

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

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

這樣分裂后，三個(gè)方向的矩陣算子均可以對(duì)角化，將無(wú)法對(duì)角化的源項(xiàng)矩陣移到最后，與單相流動(dòng)求解相比，僅需在三個(gè)方向掃描后增加一次矩陣求逆。

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

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

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

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

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

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

2 計(jì)算結(jié)果

2.1 NACA0015水翼空化計(jì)算

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

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

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

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

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

2.2 錐柱體算例

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

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

3 結(jié) 論

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

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

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

下一階段將對(duì)通氣空化進(jìn)行驗(yàn)證計(jì)算研究，同時(shí)考慮加入六自由度計(jì)算模塊模擬出水現(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).劉樺, 李家春, 何友聲, 等. 十一五水動(dòng)力學(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).周景軍. 通氣超空泡流動(dòng)及航行體流體動(dòng)力數(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]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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

牟斌, 江雄, 王建濤. 空化流動(dòng)隱式求解方法研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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