呂亞飛，張孟杰，劉濤濤，*，陳 杰，趙 欣

（1.北京理工大學(xué) 機(jī)械與車輛學(xué)院，北京 100081；2.中國(guó)北方車輛研究所，北京 100072）

0 引言

水下垂直發(fā)射航行體具有隱蔽性、機(jī)動(dòng)性和突擊性等優(yōu)點(diǎn)，已成為各軍事大國(guó)戰(zhàn)略發(fā)展的重要方向。航行體出水載荷、彈道及姿態(tài)是水下發(fā)射技術(shù)的關(guān)鍵參數(shù)，由于穿越自由液面過程中存在復(fù)雜因素的干擾，如何確保航行體出水參數(shù)滿足要求，已經(jīng)成為亟待解決的水下發(fā)射技術(shù)難點(diǎn)之一[1]。近年來，國(guó)內(nèi)外均發(fā)展了通氣空泡技術(shù)來改善航行體在出水過程中的力學(xué)環(huán)境。作為一種特殊的通氣空泡技術(shù)，多孔排氣形成的覆蓋航行體局部表面的等壓氣膜，可以有效改善航行體出水載荷、彈道及姿態(tài)[2]，該技術(shù)相關(guān)的科學(xué)問題近年來被提出并受到廣泛關(guān)注。

多孔排氣流場(chǎng)調(diào)控技術(shù)是利用發(fā)射前預(yù)置的高壓氣腔作為氣源，在環(huán)境降壓及貼體水流剪切下依靠向航行體邊界層注入不可凝氣體而形成，伴隨著航行體的出水運(yùn)動(dòng)排氣氣泡持續(xù)空間擴(kuò)張，逐漸由條帶狀氣束發(fā)展成覆蓋航行體較大部分外表面的周向融合的層狀氣膜[3]。出于研究問題的復(fù)雜性及研究對(duì)象的敏感性，國(guó)內(nèi)外有關(guān)該技術(shù)的文獻(xiàn)較少，主要集中在簡(jiǎn)單機(jī)理試驗(yàn)的初步探索，聚焦于氣泡形態(tài)宏觀發(fā)展演化過程的描述以及氣膜形成后的整體流體動(dòng)力特性的分析[4-7]，從一定程度上提高了人們對(duì)多孔排氣流場(chǎng)結(jié)構(gòu)演化過程的認(rèn)知。

然而，從流動(dòng)特征來看，多孔排氣過程本質(zhì)上為近壁面邊界層同源多束射流流動(dòng)問題，表現(xiàn)為典型的橫射流流動(dòng)過程。對(duì)于單相介質(zhì)橫射流流動(dòng)，相關(guān)的研究表明[8-11]，流場(chǎng)發(fā)展過程中會(huì)呈現(xiàn)出反向渦對(duì)、馬蹄渦、尾跡渦和環(huán)狀剪切層渦等復(fù)雜的多尺度旋渦結(jié)構(gòu)，這些旋渦結(jié)構(gòu)在湍流形成、卷吸和流體摻混過程中起到了主導(dǎo)作用。相比于傳統(tǒng)的單相橫射流流動(dòng)，水下發(fā)射多孔排氣形成的流場(chǎng)會(huì)更加復(fù)雜，其主要原因在于氣液兩相介質(zhì)物理屬性的巨大差異使得旋渦結(jié)構(gòu)在向下游發(fā)展的過程中呈現(xiàn)出更強(qiáng)的各向異性和非定常特性。目前針對(duì)水下多孔排氣多相湍流流動(dòng)的 CFD數(shù)值模擬，主要采用標(biāo)準(zhǔn)或修正 RANS模型進(jìn)行求解[12-17]，加深了人們對(duì)大尺度空泡“回射流”推進(jìn)、非定常空泡的流動(dòng)分離、旋渦脫瀉等宏觀尺度流動(dòng)的理解，但受限于商業(yè)軟件在數(shù)學(xué)模型和求解器的靈活性，仍然無法做到精確捕捉氣體射流與水體主流相互干擾形成的渦系結(jié)構(gòu)呈現(xiàn)出的旋渦對(duì)并、合攏等復(fù)雜非定?，F(xiàn)象，導(dǎo)致預(yù)測(cè)的流場(chǎng)結(jié)構(gòu)與實(shí)驗(yàn)結(jié)果存在較大誤差。

為了進(jìn)一步探究水下航行體多孔排氣流場(chǎng)演化機(jī)制，本文以平板繞流為研究對(duì)象，基于發(fā)展的近壁面邊界層通氣氣液兩相流動(dòng)高精度算法，分析了近壁面雙孔通氣氣泡的非定常發(fā)展演化過程，探討了流場(chǎng)的旋渦結(jié)構(gòu)分布與氣泡發(fā)展之間的關(guān)系。

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

1.1 基本控制方程

本文的研究對(duì)象為氣液兩相流動(dòng)，考慮到均相流模型具有較好的收斂性和計(jì)算經(jīng)濟(jì)性，研究中將氣液兩相混合物當(dāng)作一種密度可變的均勻介質(zhì)，從而共同求解一套控制方程。同時(shí)，為了簡(jiǎn)化計(jì)算，將涉及到的氣相和液相都視為不可壓縮流體，即把流體的密度視為常量。對(duì)于不可壓縮流體，笛卡爾坐標(biāo)系下的連續(xù)性方程和動(dòng)量方程可表述為

式中，α為液相體積分?jǐn)?shù)。

1.2 VOF方法

本文采用VOF方法對(duì)氣液界面進(jìn)行追蹤捕捉。整個(gè)計(jì)算在笛卡爾坐標(biāo)系下進(jìn)行，并通過液相體積分?jǐn)?shù)α進(jìn)行界面追蹤和重構(gòu)，即滿足以下定義：

液相體積分?jǐn)?shù)α滿足連續(xù)性方程

混相流體密度ρ和ν動(dòng)力粘性系數(shù)為

但是，任何事物都存在矛盾體，即“福兮禍所伏”，要想得之先于與之，好東西的獲得通常會(huì)要求一定的付出。經(jīng)濟(jì)的急速發(fā)展，伴隨而至的是環(huán)境被破壞、河水被污染、生態(tài)失衡、氣候變暖等等，人類終于為了現(xiàn)代化生活而付出了慘痛的代價(jià)。中國(guó)人民的幸福生活一半是痛苦的，生態(tài)破壞、環(huán)境惡化給人們帶來的痛心的教訓(xùn)。

式中，ρwater和ρa(bǔ)ir分別為液相和氣相密度，νwater和νair分別為液相和氣相動(dòng)力粘性系數(shù)。

1.3 求解方法及流程

在建立了控制方程的基礎(chǔ)上，應(yīng)用 VOF方法對(duì)兩相界面進(jìn)行追蹤捕捉，并獲得相界面法向量與曲率，建立的數(shù)值計(jì)算方法的求解計(jì)算流程如圖 1所示。

圖1 數(shù)值模擬計(jì)算順序流程圖Fig.1 Flow chart of numerical simulation

第3步根據(jù)不可壓縮連續(xù)性方程可以得到公式（10）左端的速度微分項(xiàng)等于0，即

可以看出，本文采用的兩步投影法將原來不可壓縮的N-S方程進(jìn)行了分步求解，通過計(jì)算過程流體速度 u*并且利用不可壓縮連續(xù)性方程構(gòu)造了一般形式的泊松方程，解除了對(duì)傳統(tǒng)RANS湍流模型的依賴，實(shí)現(xiàn)了對(duì)湍流轉(zhuǎn)捩過程的精確捕捉。

1.4 計(jì)算域及邊界條件的設(shè)置

圖2 計(jì)算域及邊界條件設(shè)置示意圖Fig.2 Schematic diagram of computational domain and boundary conditions

2 結(jié)果與討論

2.1 數(shù)值結(jié)果驗(yàn)證

圖3 給出了不同通氣速度下數(shù)值模擬計(jì)算得到的氣泡形態(tài)與實(shí)驗(yàn)結(jié)果的對(duì)比。

圖3 不同通氣速度下實(shí)驗(yàn)觀測(cè)的氣泡形態(tài)與數(shù)值預(yù)測(cè)結(jié)果的對(duì)比Fig.3 Comparison of experimental and numerical cavity shape under different ventilation speeds

從圖中可以看出，本文建立的數(shù)值模擬方法計(jì)算得到的氣泡發(fā)展過程與實(shí)驗(yàn)結(jié)果具有較好的一致性，即氣體從通氣孔進(jìn)入流場(chǎng)后緊貼平板壁面形成連續(xù)條狀氣泡，氣泡前端為連續(xù)透明氣相區(qū)，氣液界面較為光滑；隨著流動(dòng)向下游發(fā)展，氣泡逐漸向兩側(cè)膨脹，并在 K-H不穩(wěn)定性的影響下氣液界面出現(xiàn)明顯的脈動(dòng)，氣、液兩相發(fā)生強(qiáng)烈摻混，導(dǎo)致氣泡呈現(xiàn)出不透明的混相區(qū)。同時(shí)，隨著通氣速度的逐漸增大，氣泡界面處氣、液兩相速度差增大，氣液界面失穩(wěn)位置逐漸向通氣孔靠近，連續(xù)透明區(qū)逐漸縮小。由于實(shí)驗(yàn)中無法完全確保氣體持續(xù)均勻通入，實(shí)驗(yàn)觀測(cè)到的連續(xù)透明區(qū)的氣泡界面與數(shù)值計(jì)算結(jié)果存在一定的誤差。為了進(jìn)一步定量說明數(shù)值計(jì)算方法的準(zhǔn)確性，圖4給出了實(shí)驗(yàn)測(cè)量得到的不同通氣速度下氣泡擴(kuò)張角與數(shù)值計(jì)算的對(duì)比，其中實(shí)驗(yàn)結(jié)果為多次測(cè)量的平均值，可以看出，數(shù)值計(jì)算得到的氣泡擴(kuò)張趨勢(shì)與實(shí)驗(yàn)結(jié)果也吻合較好。

圖4 不同通氣速度下實(shí)驗(yàn)測(cè)量的氣泡擴(kuò)張角與數(shù)值結(jié)果的對(duì)比Fig.4 Comparison of experimental and numerical cavity expansion angle sunder different ventilation speeds

2.2 通氣氣泡的發(fā)展演化過程

圖5給出了來流速度 U∞=5 m/s、通氣速度Ug=10 m/s下平板近壁面雙孔通氣氣泡形態(tài)隨時(shí)間的演變過程，其中氣相體積分?jǐn)?shù)取0.9。從圖中可以看出，雙孔通氣氣泡的發(fā)展呈現(xiàn)出明顯的3個(gè)階段：第1個(gè)階段為相對(duì)穩(wěn)定階段，在此發(fā)展階段，氣體在水流的沖擊下迅速沿流向彎折，由與水流正交變?yōu)榕c水流流動(dòng)方向平行并逐漸向下游運(yùn)動(dòng)，形成貼附于平板表面且具有一定厚度和展向?qū)挾鹊哪顨馀?，此時(shí)氣液界面清晰且光滑，氣泡內(nèi)部主要由氣相介質(zhì)構(gòu)成，氣泡形態(tài)較為穩(wěn)定，兩股氣泡之間呈現(xiàn)出獨(dú)立發(fā)展趨勢(shì)，沒有出現(xiàn)明顯的混摻融合現(xiàn)象。第2階段為氣液界面失穩(wěn)初步融合階段，此階段，在K-H不穩(wěn)定性的影響下，氣液界面開始出現(xiàn)失穩(wěn)，且越靠近氣泡尾部氣液界面變形越大，同時(shí)氣泡不斷卷吸環(huán)境水體進(jìn)入其內(nèi)部，造成兩股氣泡間開始接觸并發(fā)生混摻融合。由于氣液兩相間的摻混過程，氣泡末端開始出現(xiàn)湍流結(jié)構(gòu)，在湍流的作用下氣泡末端呈現(xiàn)出大量離散的多種尺度氣泡雜亂分布。第3階段為充分接觸融合階段，在此階段，氣液界面失穩(wěn)加劇，氣泡末端湍流區(qū)不斷擴(kuò)大，湍流流場(chǎng)紊亂的結(jié)構(gòu)則進(jìn)一步加劇氣液之間的摻混，引起兩股氣泡束充分接觸并融合為單一氣膜向下游發(fā)展，此時(shí)氣泡界面模糊，呈現(xiàn)出小尺度的離散空泡。另外，值得注意的是，隨著氣泡間的摻混融合過程的不斷發(fā)展，氣液界面失穩(wěn)位置和氣泡間初步接觸融合位置都發(fā)生了改變，氣液界面失穩(wěn)位置由初期的x/l=10提前至x/l=17，而氣泡間初步接觸融合位置則由初期的x/l=15提前至x/l=10，這說明氣泡間的摻混融合過程對(duì)上游氣泡的發(fā)展產(chǎn)生了明顯的影響。

圖5 平板近壁面雙孔通氣氣泡形態(tài)隨時(shí)間的演變過程Fig.5 Evolution of cavity shape for double-hole ventilated cavity near wall of the flat plate

圖6給出了平板壁面（y/l=0）處的壓力隨時(shí)間演變過程，其中黑色輪廓線為氣體體積分?jǐn)?shù)為0.5時(shí)所表征的通氣氣泡形態(tài)。從圖中可以看到，當(dāng)氣體從孔口通入流場(chǎng)后，來流發(fā)生繞流作用，在孔口附近的兩股通氣氣泡束前端及其周圍形成高壓區(qū)。隨著通氣氣泡的發(fā)展，兩股通氣氣泡束的相對(duì)穩(wěn)定區(qū)內(nèi)的壓力始終保持較低均勻分布，且比較穩(wěn)定；而x/l=6位置后的水氣摻混區(qū)域的壓力分布隨著通氣氣泡束界面的波動(dòng)呈現(xiàn)波動(dòng)性分布。隨著氣泡的進(jìn)一步發(fā)展，兩股通氣氣泡束發(fā)生摻混融合作用，兩股氣泡束間出現(xiàn)相對(duì)較強(qiáng)的高壓、低壓間隔分布，直至充分融合后，氣泡內(nèi)部的壓力表現(xiàn)為較低壓、穩(wěn)定分布。

圖6 平板表面壓力隨時(shí)間的演變過程Fig.6 Evolution of pressure of flat plate surface with time

圖7進(jìn)一步給出了y/l=0平板壁面處的周向速度隨時(shí)間演變過程，其中黑色輪廓線為氣體體積分?jǐn)?shù)0.9所表征的通氣氣泡形態(tài)。速度的正負(fù)代表著其周向運(yùn)動(dòng)方向，從圖中可以看到，氣體通入流場(chǎng)后，再水流來流的擾流作用下，孔口位置的通氣氣泡內(nèi)部產(chǎn)生均朝向孔中心位置的周向速度；同時(shí)，由于通氣氣泡束對(duì)周圍環(huán)境水體的卷吸作用，在孔口附近x/l=0～5.5范圍內(nèi)所生成的穩(wěn)定通氣氣泡外部產(chǎn)生均朝向通氣空泡外側(cè)的周向速度。隨著氣泡的發(fā)展，兩股通氣氣泡束之間不斷發(fā)生著摻混融合作用，兩股通氣氣泡束內(nèi)側(cè)的方向相反的周向速度不斷增強(qiáng)；而隨著通氣氣泡束向周圍環(huán)境水體的擴(kuò)散作用，兩者之間不斷發(fā)生著動(dòng)量交換，從而導(dǎo)致兩股通氣氣泡束外側(cè)的周向速度逐漸減弱直至趨于穩(wěn)定。此外，當(dāng)兩股通氣氣泡束充分融合形成單一空泡后，通氣氣泡尾部的斷裂脫落位置存在較強(qiáng)的周向速度分布。

圖7 平板表面周向速度隨時(shí)間的演變過程Fig.7 Evolution of circumferential velocity of flat plate surface with time

2.3 旋渦結(jié)構(gòu)特性

旋渦結(jié)構(gòu)的生成和演化對(duì)平板近壁面通氣氣泡流場(chǎng)的發(fā)展演變過程中起著關(guān)鍵性的作用。為了進(jìn)一步分析流場(chǎng)的結(jié)構(gòu)演化過程，圖8給出了t5=42.5 ms時(shí)刻下流向速度u染色的Q=0.1等值面所表征的單孔通氣氣泡流場(chǎng)中的旋渦結(jié)構(gòu)分布。從圖中可以看出在近壁面通氣氣液流場(chǎng)中，分布著多種不同尺度的旋渦結(jié)構(gòu)，主要包括馬蹄渦、反向旋轉(zhuǎn)渦對(duì)（CVP）、剪切層渦、發(fā)夾渦、Λ渦等典型旋渦結(jié)構(gòu)。隨著氣泡的形成，孔口位置處的水流繞流作用相對(duì)較強(qiáng)，從而在通氣孔前緣處生成明顯的馬蹄渦結(jié)構(gòu)。馬蹄渦結(jié)構(gòu)分布在相對(duì)穩(wěn)定區(qū)域，其與通氣氣泡在前緣分離開并表現(xiàn)為獨(dú)立發(fā)展?fàn)顟B(tài)，分布在通氣氣泡兩側(cè)的渦腿比較長(zhǎng)，達(dá)到x/l=4.5的位置處，并呈現(xiàn)為對(duì)稱分布。同時(shí)在水流的沖擊作用下，通入氣體所形成的氣泡迅速彎曲，其斷面形狀呈腎形，而在腎形截面內(nèi)出現(xiàn)空間尺度、強(qiáng)度較大的反向旋轉(zhuǎn)渦對(duì)，即CVP渦。剪切層渦卷則發(fā)生在空泡界面和平板近壁面區(qū)域，主要由流體間的速度差異所導(dǎo)致，其中氣泡界面處剪切層渦主要分布在相對(duì)穩(wěn)定區(qū)和轉(zhuǎn)捩區(qū)相接區(qū)域，且隨著氣液界面的波動(dòng)呈規(guī)則排列，而壁面處剪切層渦隨著氣泡的發(fā)展沿流向不斷增大，分布相對(duì)比較規(guī)則。另外，值得注意的是，隨著氣泡的失穩(wěn)破碎，其在流場(chǎng)下游開始斷裂脫落，流動(dòng)比較紊亂，從而在湍流區(qū)域生成大量的發(fā)夾渦，如Λ渦、Ω渦等。

圖8 Q等值面表征的近壁面單孔通氣氣泡流場(chǎng)中的旋渦結(jié)構(gòu)分布Fig.8 Vortex structure distribution presented by Q for single-hole ventilated cavity near wall of the flat plate

圖9則進(jìn)一步給出了t5=42.5 ms時(shí)刻下流向速度u染色的Q=0.1等值面所表征的雙孔通氣氣泡流場(chǎng)中的旋渦結(jié)構(gòu)的整體分布。從圖中可以看出，隨著氣泡的發(fā)展，氣泡兩側(cè)生成的壁面剪切層渦不斷卷吸著周圍環(huán)境水體向展向擴(kuò)張，并在兩孔中間一定位置發(fā)生接觸生成新的剪切層渦，新的剪切層渦卷的高度要明顯低于周圍兩束通氣氣泡界面處的剪切層渦卷。這說明氣泡束之間的初步摻混融合過程是由近壁面處的剪切層渦卷主導(dǎo)的，而非氣泡界面處的剪切層渦卷。同時(shí)，隨著兩股氣泡間摻混融合的進(jìn)一步發(fā)展，融合后產(chǎn)生的剪切層渦卷發(fā)生一定程度的抬升，并逐漸發(fā)展成為氣泡界面處的剪切層渦卷。

圖9 Q等值面表征的近壁面雙孔通氣氣泡流場(chǎng)中的旋渦結(jié)構(gòu)分布Fig.9 Vortex structure distribution presented by Q for double-hole ventilated cavity near wall of the flat plate

為了進(jìn)一步探究反向旋轉(zhuǎn)渦對(duì)在氣泡束的摻混融合過程中的作用，圖10給出了x/l=1.5、9、15.8、27.4等4個(gè)不同橫截面內(nèi)的瞬時(shí)流向渦量ωx和橫截面內(nèi)的流線分布。從圖中可以看到，在孔中心x/l=1.5橫截面內(nèi)，2個(gè)通氣孔周圍的旋渦呈現(xiàn)完全規(guī)則的對(duì)稱分布，反向旋轉(zhuǎn)渦對(duì)也分別自平板壁面孔中心處產(chǎn)生并向內(nèi)側(cè)卷吸。在x/l=9橫截面內(nèi)，2股氣泡束內(nèi)部的旋渦分布保持相似性，根據(jù)流線分布可以觀察到2束氣泡束對(duì)應(yīng)的反向旋轉(zhuǎn)渦對(duì)在內(nèi)側(cè)開始發(fā)生相互作用，而空泡輪廓顯示此時(shí)的兩股氣泡束的界面還沒有發(fā)生接觸。這說明反向旋轉(zhuǎn)渦對(duì)也是誘導(dǎo)氣泡間發(fā)生摻混融合的原因之一。在x/l=15.8橫截面內(nèi)，氣泡內(nèi)側(cè)的壁面剪切層發(fā)生相互作用并出現(xiàn)接觸交混現(xiàn)象，分布在內(nèi)側(cè)中間位置處的 CVP旋渦也相應(yīng)發(fā)生接觸并與形成的剪切層渦交匯融合在一起。隨著2股氣泡束的進(jìn)一步融合，在流場(chǎng)下游的x/l=27.4橫截面內(nèi)，2股氣泡束經(jīng)過充分融合后表現(xiàn)為單一空泡形態(tài)，此時(shí)流場(chǎng)中的旋渦分布在空泡內(nèi)部，而不再分布在空泡界面處，進(jìn)一步說明所形成的單一空泡是相對(duì)穩(wěn)定的。此外，單一空泡內(nèi)部的反向旋轉(zhuǎn)渦對(duì)變得雜亂無序，旋渦分布表現(xiàn)為紊亂、無序、小尺度渦聚集。

3 結(jié)束語

本文基于自編程序提出了一種平板近壁面通氣氣液兩相湍流流動(dòng)數(shù)值計(jì)算方法，分析了繞平板近壁面雙孔通氣氣泡的非定常發(fā)展演化過程，探討了流場(chǎng)內(nèi)部的旋渦結(jié)構(gòu)分布與氣泡發(fā)展的關(guān)系，主要結(jié)論如下：

1）該方法采用VOF方法對(duì)通氣氣泡的界面進(jìn)行追蹤，并將其結(jié)合于N-S方程中，使用兩步投影法進(jìn)行求解，從而求得流場(chǎng)速度和壓力等參數(shù)，解除了對(duì)傳統(tǒng)RANS湍流模型的依賴，實(shí)現(xiàn)了對(duì)K-H不穩(wěn)定性誘導(dǎo)的氣液界面失穩(wěn)現(xiàn)象進(jìn)行準(zhǔn)確的刻畫和對(duì)湍流轉(zhuǎn)捩過程的精確捕捉，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。

2）近壁面雙孔通氣氣泡的發(fā)展演變過程可以劃分為 3個(gè)特征流動(dòng)階段：獨(dú)立穩(wěn)定發(fā)展階段、初步摻混融合階段以及充分融合階段。兩通氣空泡束由獨(dú)立向下游發(fā)展到開始接觸發(fā)生摻混融合再到劇烈交混融合形成單一空泡，相應(yīng)的空泡界面逐漸發(fā)生失穩(wěn)。

3）壁面剪切層渦和反向旋轉(zhuǎn)渦對(duì)在近壁面雙孔通氣氣泡摻混融合過程中起著主要作用。大尺度旋渦結(jié)構(gòu)對(duì)周圍水體的卷吸作用，會(huì)導(dǎo)致通氣氣泡束的壁面剪切層渦發(fā)生接觸并促進(jìn)摻混融合，隨著單一空泡的形成，產(chǎn)生新的剪切層渦卷并逐漸與氣泡界面的剪切層渦卷融合。