胡俊， 侯夏伊， 于勇

(北京理工大學(xué) 宇航學(xué)院，北京 100081)

空化涉及非定常、可壓縮、湍流脈動(dòng)以及質(zhì)量傳遞等多種復(fù)雜的流動(dòng)現(xiàn)象,普遍存在于流體機(jī)械、航天、船舶等領(lǐng)域[1-3]. 目前，數(shù)值模擬已經(jīng)成為空化現(xiàn)象的主要研究方法，空化模型的建立和應(yīng)用是研究空化問題的重點(diǎn)和難點(diǎn)[4]. 工程上常用的空化模型主要可分為三類：第一類空化模型基于簡化的Rayleigh-Plesset方程，如Schnerr-Sauer空化模型[5]，全空化模型[6]及Zwart-Gerbera-Belamri空化模型[7]等，該模型數(shù)值穩(wěn)定性和計(jì)算收斂性較好，廣泛應(yīng)用于空化流研究；第二類為基于經(jīng)驗(yàn)的比例空化模型，如Merkle空化模型[8]，上述模型均采用多個(gè)經(jīng)驗(yàn)系數(shù)，影響了模型的通用性；第三類為基于界面動(dòng)力學(xué)的空化模型，如Senocak提出的界面動(dòng)態(tài)模型[9]，該類模型由界面質(zhì)量和動(dòng)量守恒出發(fā)，有效消除了經(jīng)驗(yàn)系數(shù)對(duì)空化模型的影響.

目前研究空化流常用的幾種空化模型均涉及經(jīng)驗(yàn)參數(shù)，如全空化模型中的氣核質(zhì)量分?jǐn)?shù)；Schnerr-Sauer空化模型中的氣泡的數(shù)密度n；Z-G-B空化模型中的氣核體積分?jǐn)?shù)αnuc和氣泡半徑Rb，這些參數(shù)的取值會(huì)對(duì)空化流場(chǎng)的數(shù)值模擬結(jié)果產(chǎn)生很大影響. 劉艷等[10]認(rèn)為采用全空化模型對(duì)二維水翼空化流場(chǎng)計(jì)算時(shí)，空腔長度及厚度隨模型中氣核質(zhì)量分?jǐn)?shù)值增大而增大. Zhu等[11]認(rèn)為Schnerr-Sauer空化模型中的氣泡數(shù)密度n的取值對(duì)空化流場(chǎng)數(shù)值模擬結(jié)果有很大影響. Asnaghi等[12]基于大渦模擬的方法，引入瞬時(shí)剪切應(yīng)變率對(duì)Schnerr-Sauer空化模型的蒸發(fā)項(xiàng)系數(shù)進(jìn)行修正，并對(duì)二維NACA0009 MOD水翼的片狀空化流場(chǎng)進(jìn)行計(jì)算. 結(jié)果表明該方法可以有效消除考慮了氣核的Schnerr-Sauer空化模型中經(jīng)驗(yàn)參數(shù)的取值對(duì)數(shù)值模擬結(jié)果的影響.

本文在前人研究的基礎(chǔ)上，基于RANS方法，引入當(dāng)?shù)亓鲃?dòng)特征，采用平均剪切應(yīng)變率來修正Z-G-B空化模型的蒸發(fā)項(xiàng)系數(shù). 將修正后的空化模型通過UDF(user define function)嵌入FLUENT18.2商業(yè)軟件，采用Mixture多相流模型和Realizablek-ε湍流模型對(duì)攻角為2.5°，雷諾數(shù)為2×106的二維NACA0009 MOD水翼的空化流場(chǎng)以及攻角為8°，雷諾數(shù)為7×108的Clark-Y水翼的初生空化、片狀空化、云空化現(xiàn)象進(jìn)行數(shù)值模擬. 將數(shù)值模擬結(jié)果與已有文獻(xiàn)中的實(shí)驗(yàn)結(jié)果對(duì)比，分析了氣核體積分?jǐn)?shù)αnuc及氣泡半徑Rb的取值對(duì)系數(shù)修正前后的Z-G-B空化模型影響，驗(yàn)證了引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型的可行性，同時(shí)探究了引入當(dāng)?shù)亓鲃?dòng)特征修正前后的Z-G-B空化模型在模擬Clark-Y翼型的初生空化、片狀空化及云空化流場(chǎng)時(shí)的表現(xiàn).

1 數(shù)值計(jì)算方法

1.1 控制方程

采用Mixtrue多相流模型，將流場(chǎng)視為氣液兩相充分混合的單一介質(zhì)，需求解混合相的連續(xù)方程、動(dòng)量方程以及描述氣液兩相間質(zhì)量傳遞的氣相體積分?jǐn)?shù)方程為

(1)

(2)

(3)

ρm=ρvαv+ρl(1-αv)

(4)

μm=μvαv+μl(1-αv)

(5)

式中μv和μl分別為氣相動(dòng)力黏度及液相動(dòng)力黏度.

1.2 空化模型

原始Z-G-B模型基于簡化的Rayleigh-Plesset方程，并且考慮了氣核對(duì)空化初生和發(fā)展的影響，采用氣泡數(shù)密度n及氣泡半徑Rb定義氣相體積分?jǐn)?shù)

(6)

并由單個(gè)氣泡的質(zhì)量變化率直接推導(dǎo)出氣液兩相間的質(zhì)量變化率，最終改模型傳質(zhì)源項(xiàng)的表達(dá)式為

(7)

式中：Cevap和Ccond為傳質(zhì)源項(xiàng)的蒸發(fā)項(xiàng)及凝結(jié)項(xiàng)系數(shù)，分別取50和-0.01；αnuc為氣核體積分?jǐn)?shù)；Rb為氣泡半徑；Pv為空化臨界壓力. 其中，氣核體積分?jǐn)?shù)αnuc及空泡半徑Rb為該空化模型中的關(guān)鍵參數(shù)，它們的取值會(huì)對(duì)數(shù)值模擬結(jié)果產(chǎn)生一定影響.

由于空化模型中的蒸發(fā)項(xiàng)及凝結(jié)項(xiàng)系數(shù)代表了影響相變速率的相變弛豫時(shí)間，而當(dāng)?shù)亓鲃?dòng)狀態(tài)很大程度上影響相變弛豫時(shí)間，可將相變弛豫時(shí)間與局部的流動(dòng)參數(shù)建立聯(lián)系，參考文獻(xiàn)[12]引入局部瞬時(shí)剪切應(yīng)變率修正了考慮非冷凝氣體的Schnerr-Sauer空化模型，一定程度上消除了模型中經(jīng)驗(yàn)參數(shù)取值對(duì)數(shù)值模擬結(jié)果的影響，提高了空化模型的準(zhǔn)確性. 基于這一思想，本文在RANS框架下，引入時(shí)間平均剪切應(yīng)變率修正Z-G-B空化模型的系數(shù). 由于在空化現(xiàn)象中，蒸發(fā)過程對(duì)這種局部流動(dòng)狀態(tài)較為敏感，因此僅采用該方法修正Z-G-B模型的蒸發(fā)項(xiàng)系數(shù)，并用主流時(shí)間尺度對(duì)時(shí)均剪切應(yīng)變率進(jìn)行量綱—化，最終將蒸發(fā)系數(shù)修正為：

(8)

t∞=c/U∞

(9)

式中：Cmod為修正后的蒸發(fā)項(xiàng)系數(shù)；c為翼型弦長；U∞為來流速度；t∞為主流時(shí)間尺度.

1.3 計(jì)算域設(shè)置及網(wǎng)格劃分

計(jì)算域及邊界條件如圖1所示，二維NACA0009 MOD型水翼的攻角為2.5°，翼型弦長為c=100 mm. 計(jì)算區(qū)域的入口距翼型前緣2.0c，出口距翼型尾緣4.0c，上下壁面間距離為1.5c. 入口采用速度入口U∞=20 m/s，對(duì)應(yīng)雷諾數(shù)為Re=2×106，湍流強(qiáng)度設(shè)置為1%. 二維Clark-Y型水翼攻角為8°，弦長為c=70 mm，計(jì)算區(qū)域的入口距翼型前緣4.0c，出口與翼型尾緣的距離為5.0c，上下壁面間距離為2.7c. 采用速度入口U∞=10 m/s，對(duì)應(yīng)雷諾數(shù)為Re=7.0×105，湍流強(qiáng)度設(shè)置為2%. 出口均為壓力出口，上下邊界為對(duì)稱條件，翼型表面采用絕熱、無滑移固壁條件.

圖1 計(jì)算域及邊界條件

NACA009 MOD水翼的計(jì)算域及網(wǎng)格劃分與文獻(xiàn)[12]相同，將計(jì)算域劃分為9個(gè)塊，總網(wǎng)格數(shù)為34 720. 翼型附近的區(qū)域采用C型結(jié)構(gòu)化網(wǎng)格劃分，能夠較好地匹配翼型的形狀，并在翼型周圍的近壁區(qū)域進(jìn)行了網(wǎng)格加密，近壁面y+值在30～200之間，滿足壁面函數(shù)要求. 對(duì)于Clark-Y水翼的計(jì)算域，劃分了4套不同尺寸的網(wǎng)格并進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證，網(wǎng)格節(jié)點(diǎn)數(shù)分別為

N1=44 124，N2=67 596，

N3=95 834，N4=134 846，

計(jì)算對(duì)比了4套網(wǎng)格在空化數(shù)σ=3.0時(shí)Clark-Y水翼的時(shí)均升阻力系數(shù). 計(jì)算結(jié)果如表1所示：當(dāng)網(wǎng)格數(shù)大于N3時(shí)其升阻力系數(shù)差別不大，且均與實(shí)驗(yàn)值[13]吻合較好，故選擇N3作為本文的計(jì)算網(wǎng)格.

表1 Clark-Y翼型網(wǎng)格無關(guān)性驗(yàn)證

2 結(jié)果與討論

2.1 氣核體積分?jǐn)?shù)αnuc對(duì)修正前后Zwart-Gerbera-Belamr空化模型的影響

為了研究氣核體積分?jǐn)?shù)αnuc取值不同時(shí)，引入剪切應(yīng)變率修正前后的Z-G-B空化模型的表現(xiàn)，在空化數(shù)σ=0.75，0.80，0.85，0.90，氣泡半徑Rb=10-6時(shí)，改變氣核體積分?jǐn)?shù)αnuc，令αnuc=5×10-3，5×10-4，5×10-5，對(duì)二維NACA0009 MOD水翼的定?？栈鲌?chǎng)進(jìn)行數(shù)值模擬并與實(shí)驗(yàn)值[14]對(duì)比.

如圖2所示，采用原始空化模型計(jì)算時(shí)，4個(gè)空化數(shù)下水翼前緣附著的空腔長度均隨αnuc取值的增大而增大. 當(dāng)αnuc=5×10-3時(shí)，數(shù)值模擬結(jié)果與實(shí)驗(yàn)測(cè)得的空腔長度及壓力分布吻合最好；當(dāng)αnuc= 5×10-5時(shí)，由于氣核體積分?jǐn)?shù)取值過小，與實(shí)際情況相差較大，計(jì)算出的水翼附著空腔長度較短，并且難以捕捉到空腔閉合區(qū)域壓力峰值，與實(shí)驗(yàn)結(jié)果難以吻合. 隨空化數(shù)σ由0.75升高到0.90，空化流場(chǎng)趨于穩(wěn)定，氣核體積分?jǐn)?shù)αnuc的取值對(duì)流場(chǎng)中空穴長度的影響逐漸減小. 修正前的Z-G-B空化模型對(duì)氣核體積分?jǐn)?shù)αnuc的取值較為敏感，且空腔長度與氣核體積分?jǐn)?shù)αnuc成正比.

圖2 采用系數(shù)修正前后模型計(jì)算出的翼面壓力分布

引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型后，αnuc=5×10-4，5×10-5時(shí)計(jì)算出的空腔長度均增大，并且能夠捕捉到空腔閉合區(qū)域的壓力峰值，其壓力分布與實(shí)驗(yàn)值更加吻合. 當(dāng)σ=0.85，0.90,氣核體積分?jǐn)?shù)αnuc= 5×10-4，5×10-3時(shí)，計(jì)算出的空腔長度基本相等，且與實(shí)驗(yàn)值吻合較好. 引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型，可以有效減小氣核αnuc的取值對(duì)數(shù)值模擬結(jié)果的影響，在一定程度上減小了Z-G-B空化模型對(duì)該參數(shù)取值的敏感性，提高了該空化模型的可靠性及準(zhǔn)確性.

2.2 氣泡半徑Rb對(duì)修正前后Z-G-B空化模型的影響

為探究氣泡半徑Rb的取值對(duì)引入剪切應(yīng)變率修正前后Z-G-B空化模型的影響，在空化數(shù)σ=0.75，0.80，0.85及0.90，氣核體積分?jǐn)?shù)αnuc=5×10-4時(shí)，改變氣泡半徑Rb，令Rb=10-4，10-5，10-6，對(duì)二維NACA0009 MOD水翼的定?？栈鲌?chǎng)進(jìn)行數(shù)值模擬.

計(jì)算結(jié)果如圖3所示，引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型前，水翼附著的空腔長度隨Rb取值的增大而減小，Rb的取值對(duì)數(shù)值模擬結(jié)果的影響也隨空化數(shù)的升高而降低. 數(shù)值模擬結(jié)果僅在Rb=10-6時(shí)與實(shí)驗(yàn)值接近，且能夠捕捉到空腔閉合區(qū)域壓力的突增. 當(dāng)Rb=10-5，10-4時(shí)，由于其取值過大，采用原始空化模型計(jì)算出的空腔長度較短，不能捕捉空腔閉合區(qū)域壓力峰值，與實(shí)驗(yàn)結(jié)果相差較大.

圖3 采用系數(shù)修正后模型計(jì)算出的翼面壓力分布

引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型后，改變Rb取值計(jì)算得到的空腔長度均大于原始模型計(jì)算出的長度.Rb=10-5，10-6時(shí)4個(gè)空化數(shù)下計(jì)算得到的空腔長度接近，與實(shí)驗(yàn)值吻合的較好. 隨空化的升高，Rb=10-5，10-6時(shí)計(jì)算得到壓力系數(shù)分布基本一致. 但當(dāng)Rb=10-4時(shí)取值過大，引入當(dāng)?shù)亓鲃?dòng)特征修正后計(jì)算出的空穴長度雖然顯著增加，且捕捉到了空腔閉合區(qū)域壓力突增的現(xiàn)象，但仍然與實(shí)驗(yàn)值存在一定差距.

綜上所述，引入當(dāng)?shù)亓鲃?dòng)特征，采用剪切應(yīng)變率修正Z-G-B空化模型，可以在一定程度上減小模型中氣核體積分?jǐn)?shù)αnuc及氣泡半徑Rb取值對(duì)數(shù)值模擬結(jié)果的影響，有效降低了該空化模型對(duì)經(jīng)驗(yàn)參數(shù)取值的敏感性，在一定程度上提高該空化模型的適用性及準(zhǔn)確性.

2.3 修正前后Z-G-B空化模型對(duì)不同空化狀態(tài)的數(shù)值模擬

為進(jìn)一步探究引入當(dāng)?shù)亓鲃?dòng)特征對(duì)Z-G-B空化模型的影響，研究原始及修正后的Z-G-B空化模型在模擬初生空化、片狀空化及云狀空化時(shí)的表現(xiàn)，采用修正前后的Z-G-B空化模型，設(shè)定氣泡半徑Rb=10-6，氣核體積分?jǐn)?shù)αnuc=5×10-4，對(duì)二維Clark-Y水翼的空化流場(chǎng)進(jìn)行非定常計(jì)算，時(shí)間步長取為Δt=10-6s，其空化數(shù)范圍為3.0～0.55.

圖4 采用系數(shù)修正前后空化模型計(jì)算出的升阻力系數(shù)

在無空化流即σ=3.0，2.5時(shí)，采用兩空化模型得到的升阻力系數(shù)基本相等. 當(dāng)空化數(shù)σ=2.0時(shí)，修正前后的空化模型均模擬出了空化現(xiàn)象，與實(shí)驗(yàn)值相比，初生空化發(fā)生得過早. 修正后的空化模型捕捉到了空化初生時(shí)升力系數(shù)略微升高的現(xiàn)象，原始模型并未捕捉到升力系數(shù)的變化. 在空化數(shù)σ=1.8時(shí)，系數(shù)修正前后的空化模型均過早地模擬了片狀空化的發(fā)生，此時(shí)阻力系數(shù)開始迅速升高，升力系數(shù)迅速降低. 在空化數(shù)σ=0.8時(shí)，采用引入當(dāng)?shù)亓鲃?dòng)特修正后的空化模型計(jì)算出的阻力到最大值，而原始空化模型計(jì)算出的阻力系數(shù)繼續(xù)升高，沒有很好地捕捉到云空化發(fā)生時(shí)阻力系數(shù)的變化規(guī)律.

對(duì)于Clark-Y水翼空化流場(chǎng)，雖然引入當(dāng)?shù)亓鲃?dòng)特征修正后的Z-G-B空化模型計(jì)算得到的平均升阻力系數(shù)沒有原始空化模型與實(shí)驗(yàn)值吻合得好，但修正后的空化模型能夠更好地捕捉到空化流場(chǎng)發(fā)生初生空化、片狀空化及云空化時(shí)升阻力系數(shù)的變化趨勢(shì).

3 結(jié) 論

在RANS框架下，基于當(dāng)?shù)亓鲃?dòng)特征，引入平均剪切應(yīng)變率修正了Z-G-B空化模型. 采用修正前后的空化模型聯(lián)立Realizablek-ε湍流模型對(duì)二維NACA0009 MOD水翼空化流場(chǎng)及Clark-Y水翼空化流場(chǎng)進(jìn)行了數(shù)值模擬，并與已有實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比. 得到以下主要結(jié)論：

① 原始Z-G-B空化模型對(duì)氣核體積分?jǐn)?shù)αnuc及氣泡半徑Rb的取值較為敏感. 在NACA0009 MOD空化流場(chǎng)中空腔長度隨氣泡半徑Rb的取值減小，隨氣核體積分?jǐn)?shù)αnuc的增大而增大，且空化數(shù)越小時(shí)模型參數(shù)取值的改變對(duì)數(shù)值模擬結(jié)果的影響越大.

② 當(dāng)空化模型中的參數(shù)取值在一定范圍內(nèi)時(shí)，引入當(dāng)?shù)亓鲃?dòng)特征修正Z-G-B空化模型可以有效地減小氣核體積分?jǐn)?shù)αnuc及氣泡半徑Rb的取值對(duì)數(shù)值模擬結(jié)果的影響，降低了該空化模型對(duì)經(jīng)驗(yàn)參數(shù)取值的敏感性，提高了Z-G-B空化模型的準(zhǔn)確性及適用性.

③ 對(duì)于Clark-Y水翼空化流的數(shù)值模擬，雖然修正后的空化模型計(jì)算出的升阻力系數(shù)與實(shí)驗(yàn)值有一定差距，但修正后的模型能夠定性捕捉到空化初生時(shí)升力系數(shù)的略微升高及云空化發(fā)生時(shí)阻力系數(shù)達(dá)到峰值的現(xiàn)象.