馬慶鵬，魏英杰，王聰，曹偉，陳超倩

(哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱150001)

0 引言

運(yùn)動體以一定初速度從空氣中穿越氣、水交界面進(jìn)入水中并在其運(yùn)動軌跡上產(chǎn)生一個狹長空穴這一過程稱為入水過程，其后的空穴稱之為入水空泡。早在19世紀(jì)末，Worthington 等［1］便利用高速閃光照相機(jī)、玻璃水槽等設(shè)備開展了鋼制球體垂直入水的小型實(shí)驗(yàn)，對球體入水空泡發(fā)展規(guī)律開展了研究，闡述了球體入水各個階段的流動現(xiàn)象。此后的一個多世紀(jì)以來，由于入水問題在諸多領(lǐng)域的廣泛應(yīng)用，尤其是諸如空投魚雷、入水射彈等入水武器方面的應(yīng)用，入水問題的研究取得了較大的進(jìn)展。

入水過程可以分為6 個階段［2］:沖擊階段、流動形成階段、空泡敞開階段、空泡閉合階段、空泡潰滅階段以及完全沾濕階段。其中，空泡閉合階段又分為面閉合和深閉合。面閉合是由于表面噴濺在大氣壓力、大氣密度、空泡口空氣流動等多種因素的共同作用下迅速向中間收縮，從而形成一個類似于凸起的鐘罩形狀的封閉面;深閉合是指空泡排開水體徑向擴(kuò)張動能降為0 后在流體靜壓作用下向軸向收縮并最終收縮于空泡軸線上一點(diǎn)這一流動現(xiàn)象?？张蓍]合尤其是深閉合是導(dǎo)致空泡潰滅的主要原因。

對于空泡閉合的研究多通過實(shí)驗(yàn)完成，研究對象多為球體、錐體、圓盤等結(jié)構(gòu)體。美國海軍軍械實(shí)驗(yàn)室的Gilbarg 等［3］通過實(shí)驗(yàn)研究了大氣壓力對球體入水噴濺和空泡的影響規(guī)律。May 等［4-6］針對球體、錐體、圓盤等多種形狀結(jié)構(gòu)體，開展了大量的低速入水實(shí)驗(yàn)，得到了低速入水空泡的發(fā)展過程及環(huán)境條件對空泡面閉合的影響規(guī)律，并分析了入水各階段結(jié)構(gòu)受到的流體動力。

20世紀(jì)后期，學(xué)者對入水問題的研究多集中在入水拍擊階段，通過理論分析及數(shù)值計(jì)算的方法對拍擊載荷及噴濺的初始形態(tài)開展了大量研究。文獻(xiàn)［7］采用邊界元法求解了二維楔形體的入水拍擊問題，得到了噴濺射流的形態(tài)發(fā)展規(guī)律及楔體表面壓力系數(shù)。文獻(xiàn)［8］針對非對稱二維楔形體入水拍擊過程開展了理論研究，分析了拍擊過程水面噴濺的發(fā)展規(guī)律。文獻(xiàn)［9］采用邊界元法求解了三維錐體垂直、傾斜條件下的入水拍擊問題，得到了不同拍擊狀態(tài)水面噴濺的形態(tài)。

以上研究主要針對入水空泡的發(fā)展機(jī)理及入水初級階段的載荷、噴濺形態(tài)及面閉合等問題，對于空泡深閉合的研究較少。此外，目前對于入水問題的研究大多為低速條件下，對于高速條件下(大于150 m/s)的實(shí)驗(yàn)及理論研究，Lundstrom［10］在美國軍方支持下開展了12.7 ～14.5 mm 穿甲彈在805 ～1 070 m/s 速度范圍內(nèi)的入水實(shí)驗(yàn)，并提出了一種基于能量守恒的空泡形態(tài)預(yù)測方法。Lee 等［11］在此基礎(chǔ)上提出了空泡深閉合的理論分析方法，對初始速度為500 ～1 500 m/s 的球體入水空泡發(fā)展過程開展了相關(guān)研究，并與數(shù)值模擬結(jié)果進(jìn)行了比較。

21世紀(jì)以來，計(jì)算機(jī)技術(shù)的發(fā)展使得計(jì)算流體力學(xué)(CFD)方法得到了更多的應(yīng)用。王聰?shù)龋?2］、何春濤等［13］針對小型結(jié)構(gòu)體低速(10 ～50 m/s)入水問題開展了全過程數(shù)值模擬研究，分析了入水空泡發(fā)展規(guī)律，并研究了空氣域壓強(qiáng)對入水空泡面閉合的影響規(guī)律。

本文基于有限體積法和流體體積(VOF)多相流模型，求解了雷諾平均的N-S 方程，考慮入水過程的空化及湍流現(xiàn)象，針對入水過程的深閉合問題，對半錐角為45°的錐頭圓柱體以500/s 初始速度自由垂直入水過程開展了數(shù)值模擬研究。主要分析了入水空泡形態(tài)發(fā)展及深閉合現(xiàn)象。得到了深閉合發(fā)生的時間及深閉合前后流場的流動特性及壓力分布變化特性。

1 數(shù)學(xué)模型

1.1 流體控制方程

本文數(shù)值計(jì)算采用VOF 多相流模型描述水、氣、汽形成的多相流動，同時假設(shè)流體為不可壓縮，忽略入水過程中由于流體粘性產(chǎn)生的熱效應(yīng)，建立了該問題的流體控制方程。其中，水相、氣相及水蒸汽相的體積分?jǐn)?shù)分別用αl、αg、αv表示，三者滿足關(guān)系式:

混合介質(zhì)的連續(xù)性方程為

動量守恒方程為

(2)式和(3)式中:ui為速度分量;ρm和μm分別為混合介質(zhì)的密度和動力粘度;μt=ρmCμk2/ε 為湍流動力粘度。ρm和μm表達(dá)式分別為

對于入水過程的空化問題，本文采用Schnerr-Sauer 空化模型進(jìn)行求解，該模型的守恒方程基于水蒸汽相建立，水蒸汽相體積分?jǐn)?shù)的輸運(yùn)方程為

式中:RB=1×10-6m 為氣核半徑;αnuc=5×10-4為不可凝結(jié)氣體體積分?jǐn)?shù);Fvap= 50 和Fcond=0.001為經(jīng)驗(yàn)常數(shù)。

1.2 湍流模型

運(yùn)動體高速入水過程具有強(qiáng)瞬時、非定常、高載荷等特性，伴隨著湍流、相變等復(fù)雜的流動過程。對于這一強(qiáng)非線性流動，本文采用k-ω SST 湍流模型進(jìn)行計(jì)算。k-ω SST 模型由Menter［14］提出，它綜合了k-ω 模型近壁穩(wěn)定性和k-ε 模型邊界層外部獨(dú)立性的優(yōu)點(diǎn)。

2 數(shù)值計(jì)算

2.1 計(jì)算模型及網(wǎng)格劃分

本文計(jì)算的運(yùn)動體外形、尺寸以及計(jì)算域二維示意圖如圖1所示。運(yùn)動體為錐頭圓柱體，柱段直徑D=10 mm，柱段長度L =5D，錐頭半錐角為45°.計(jì)算域?yàn)閳A柱形，為排除壁面邊界對中心流場影響，取計(jì)算域直徑Dd=100D. 計(jì)算域長度為230D，其中自由液面上方空氣域長度為32D，水域長度198D，運(yùn)動體錐頭頂端在初始狀態(tài)下距自由液面2D 處。

圖1 運(yùn)動體及計(jì)算域示意圖Fig.1 Sketch of projectile and computational domain

考慮到運(yùn)動體及計(jì)算域幾何形狀都為軸對稱體，本文采用二維軸對稱模型開展數(shù)值模擬研究。流域網(wǎng)格劃分及邊界條件設(shè)置如圖2所示。為保證空氣域、汽水交界面以及空泡區(qū)域計(jì)算結(jié)果精度，對計(jì)算區(qū)域網(wǎng)格進(jìn)行了局部加密處理。運(yùn)動體周圍，軸向前后2D、徑向6D 范圍計(jì)算域采用三角形網(wǎng)格加密，其余為四邊形網(wǎng)格，運(yùn)動體表面第一層網(wǎng)格厚度為D/400，加密區(qū)與外部過渡網(wǎng)格厚度為D/20.運(yùn)動體上方沿x 軸負(fù)向網(wǎng)格高度為D/20 均布，前場沿x 軸正向漸疏，初始網(wǎng)格數(shù)量為275 768.

圖2 計(jì)算網(wǎng)格示意圖Fig.2 Sketch of the mesh

2.2 數(shù)值方法

本文采用有限體積法對流體控制方程離散，采用PISO 算法對壓力-速度場的耦合進(jìn)行求解;壓力場的空間離散采用PRESTO!格式;各相體積率離散采用CICSAM 格式;綜合考慮收斂性與計(jì)算時間，對動量方程的離散采用一階迎風(fēng)格式。

計(jì)算中引入U(xiǎn)DF 自編程實(shí)現(xiàn)運(yùn)動體沿x 軸正向運(yùn)動，并采用動網(wǎng)格技術(shù)定義網(wǎng)格運(yùn)動和更新。動網(wǎng)格更新方法采用動態(tài)層法，網(wǎng)格更新過程中，為保證空氣域計(jì)算精度，新生成網(wǎng)格高度均為D/20.

3 計(jì)算結(jié)果與分析

3.1 空泡深閉合預(yù)測模型

文獻(xiàn)［11］在文獻(xiàn)［10］的基礎(chǔ)上提出了一種預(yù)測入水空泡形態(tài)的理論分析模型，入水空泡發(fā)展過程如圖3所示。

假設(shè)流體不可壓縮，忽略入水過程熱交換產(chǎn)生的能量損失及重力影響，由能量守恒原理可知，運(yùn)動體動能損失等于排開水的動能和空泡的壓力勢能之和。在此基礎(chǔ)上，忽略運(yùn)動體尺度，文獻(xiàn)［11］給出了空泡半徑的預(yù)測公式:

式中:t≥tb，t 為運(yùn)動體當(dāng)前位置時刻，tb為x =xb時刻;A(x)和B(x)分別為

圖3 入水空泡示意圖Fig.3 Illustration of water entry cavity

(8)式、(9)式中:Ep= mv2p/2 為運(yùn)動體動能;pg=px-pc，px為該深度絕對壓強(qiáng)，pc為空泡內(nèi)部壓強(qiáng)，本文取pc為飽和蒸汽壓強(qiáng)3 540 Pa;N=ln(R/r)，R為擾動流場半徑，r 為空泡半徑，本文取R/r =30.此外，對于本文模型，需在(7)式右側(cè)加上運(yùn)動體半徑r0=D/2.

(7)式的求解依賴于運(yùn)動體各時刻的速度與位移數(shù)據(jù)，對于該運(yùn)動，由牛頓第二定律可得

式中:A0為運(yùn)動體特征面積;vp為瞬時速度;Cdx=Cd0+σ 為運(yùn)動體阻力系數(shù)(其中空化數(shù)σ=2(p∞-pc)/(ρlv2p)，p∞為遠(yuǎn)場壓力)，對于半錐角45°的圓錐頭，σ=0 時的阻力系數(shù)Cd0=0.498［2］;g 為重力加速度。因此，(10)式可改寫為

給定邊界條件vp|t=0=500 m/s，求解(11)式，便可得到任意時刻運(yùn)動體的速度及位移，代入到(7)式中便可得到給定時刻下任意xb≤xt位置的空泡半徑。

令(7)式右側(cè)等于0，可得在xb位置空泡半徑為0 的兩個時刻:

式中:t1為初始時刻，td=t2即為空泡在該點(diǎn)閉合的時刻。

由(8)式、(9)式可知，對于給定的tb，A(xb)、B(xb)皆為固定值，因此由式(13)可以求得在空泡長度范圍內(nèi)，任意深度下空泡閉合所需時間，而空泡發(fā)生深閉合位置即為xd= x(min(tb)).

3.2 入水彈道

對前文描述的半錐角45°錐頭圓柱體以500 m/s初始速度自由垂直入水問題進(jìn)行數(shù)值模擬，得到了入水過程中運(yùn)動體速度vp及入水深度x/D 隨入水時間t 的變化。另一方面，由(11)式可以得出理論分析結(jié)果，兩種結(jié)果對比如圖4所示。

圖4 數(shù)值模擬與理論分析結(jié)果對比Fig.4 Comparison of thenumerical results and analytical results

由圖4可以看出，對于運(yùn)動體速度衰減曲線，兩種方法得到的結(jié)果基本一致，在入水初期及曲線后半段吻合度極高，在中段稍有區(qū)別，數(shù)值模擬結(jié)果稍微偏大。這一區(qū)別在位移曲線上得到了體現(xiàn)，可以看出，從速度值差別稍大的時間點(diǎn)開始，位移值也開始出現(xiàn)偏差，由于積分累加的緣故，在18 ms 處，兩種位移計(jì)算結(jié)果差別最大。這一差別主要由于兩種方法基于求解不同的方程，且理論分析模型忽略了空泡面閉合對運(yùn)動體運(yùn)動的影響，因此，對于本文高速條件下入水問題，這一誤差值是可以接受的。

3.3 入水空泡

3.3.1 空泡形態(tài)

圖5給出了理論求解與數(shù)值模擬在入水深度x/D 為30D 時空泡形態(tài)結(jié)果對比。從圖5可以看出，兩種方法得到的空泡外形在自由液面以下吻合較好，由于理論分析模型未考慮液面噴濺，因此水面以上區(qū)別較大。同時，對于理論分析模型中不同的N=ln(R/r)取值，空泡半徑有一定的差異。

圖5 入水深度30D 位置空泡形態(tài)對比Fig.5 Comparison of the cavity shapes at thepenetration depth H=30D

圖6給出了入水空泡的發(fā)展過程云圖，可以看出，在整個入水過程中，空泡經(jīng)歷了形成、面閉合、擴(kuò)張、收縮、深閉合、潰滅等幾個階段。在入水后1.0 ms左右，空泡發(fā)生面閉合現(xiàn)象，入水噴濺逐漸增厚并向空泡內(nèi)部侵入;在5.0 ms 后，最大空泡直徑無明顯增大，空泡在面閉合以及流體靜壓的共同作用下開始向?qū)ΨQ軸收縮，并在大約16.5 ms 左右深閉合;隨后空泡被縱向拉斷，水體迅速向兩端壓縮，空泡迅速潰滅。

3.3.2 空泡深閉合流場分析

為確定本文計(jì)算模型的深閉合位置，對空泡收縮段的流場進(jìn)行了分析。

1)深閉合區(qū)域流線分布

圖7給出了空泡收縮段附近流場流線分布圖，其中，每一時刻采取二維軸對稱顯示，豎直向下為重力方向。由圖7可以看出，在16.30 ms 時，空泡尚未閉合，此時刻流線顯示，該段空泡外流體基本沿垂直空泡壁面方向向軸線運(yùn)動，在空泡內(nèi)部，流體向下運(yùn)動。16.30～16.38 ms 這一過程，空泡壁面逐漸靠近對稱軸，空泡附近流動逐漸復(fù)雜。隨著空泡壁面的收縮，在16.36、16.37 ms 時刻，可看到空泡內(nèi)部產(chǎn)生少量水滴懸浮在空泡內(nèi)部，阻礙了空泡內(nèi)部介質(zhì)向下運(yùn)動。隨著空泡的繼續(xù)收縮，空泡內(nèi)壁與水滴間距逐漸減小，在16.38 ms 時刻時，水滴上方流體無法再向下流動，并與水滴下方介質(zhì)在水滴位置相會，這導(dǎo)致了空泡內(nèi)部流動被阻隔，雖然此時空泡內(nèi)壁尚未與對稱軸接觸，但是實(shí)際上空泡已經(jīng)發(fā)生了深閉合現(xiàn)象?？张萆铋]合后，在內(nèi)部流體交匯處逐漸形成多處渦旋，如16.5 ms 時刻流線圖所示。這些渦旋將流場分隔為數(shù)段，從圖中水滴位置變化可明顯看出空泡分離后向上下兩端收縮的趨勢。

由文獻(xiàn)［2］可知，在低速條件下的入水問題中，空泡深閉合往往發(fā)生在一點(diǎn)，而本文算例的深閉合位置并非是清晰的一點(diǎn)，而是在一段較長區(qū)域里形成數(shù)處渦旋導(dǎo)致深閉合，這一現(xiàn)象主要由高速入水特征決定。由于本文研究的運(yùn)動體直徑較小，高速條件下入水過程較長，深閉合深度較深，空泡深閉合時運(yùn)動體動能損失較大，這就導(dǎo)致深閉合位置附近水體在軸線方向上的速度梯度較小，該區(qū)域水體在極小時間內(nèi)閉合時不易區(qū)分。另一方面，由于高速入水空化效應(yīng)，空泡內(nèi)部充滿水蒸汽、氣相及少量水相的混合介質(zhì)，空泡深閉合時水體與空泡內(nèi)部混合介質(zhì)復(fù)雜的相互作用對空泡閉合點(diǎn)有較大的影響。

2)深閉合位置結(jié)果對比分析

圖6 入水空泡發(fā)展過程云圖Fig.6 Contour of cavity evoluation

上文求得運(yùn)動體速度、位移數(shù)據(jù)后，由(8)式、(9)式、(13)式得到空泡發(fā)生深閉合的時間td及位置xd，同時，也可得到此時運(yùn)動體的入水深度xtd.對于不同的N=ln(R/r)取值，上述求解結(jié)果也有一定差異，表1給出了數(shù)值模擬求解結(jié)果以及4 種不同擾動流場半徑下的理論分析結(jié)果，其中，數(shù)值模擬的閉合取閉合段軸向中心點(diǎn)位置。

圖7 深閉合位置附近流線分布Fig.7 Streamlines around the deep closure position

表1 入水空泡深閉合時間及位置表Tab.1 Times and positions of deep closure

由表1可以看出，在N 取值不同，即擾動流場半徑R 取值不同時，理論分析結(jié)果稍有差異，N 取值越大，與數(shù)值模擬結(jié)果相對誤差越小。其中，對于深閉合位置xd，N 取值遞增時，理論分析的4個結(jié)果與數(shù)值計(jì)算結(jié)果分別相差18.52%、16.03%、14.50%、12.72%。這表明，理論分析模型在擾動流場半徑的選擇上對計(jì)算結(jié)果有一定的影響。其原因一方面是由于理論分析模型簡化較多，尤其忽略了面閉合和空泡內(nèi)部氣相的影響，另一方面在于理論分析模型忽略了模型尺寸，其閉合半徑相對于數(shù)值模擬也相差一個模型半徑。

3)深閉合區(qū)域壓力場分布

圖8給出了空泡深閉合前后閉合點(diǎn)附近壓力分布云圖及水蒸汽相云圖。從圖8可以看出，在空泡閉合接近前，壓力最大值出現(xiàn)在運(yùn)動體頭部位置，空泡內(nèi)主要由低壓的水蒸汽相構(gòu)成。隨著空泡不斷收縮，水體擠壓水蒸汽相將其排開，并在閉合點(diǎn)附近出現(xiàn)一段高壓區(qū)。當(dāng)空泡收縮到一定程度后，周向水體在對稱軸上一段區(qū)域內(nèi)閉合，在撞擊點(diǎn)處壓力達(dá)到最大。隨后，閉合水體形成兩股射流向空泡兩端運(yùn)動，進(jìn)一步擠壓空泡內(nèi)部混合介質(zhì)，并在空泡分離點(diǎn)處產(chǎn)生壓力極大值。這表明，空泡深閉合過程中，周圍水體勢能逐漸轉(zhuǎn)化為動能，在與空泡內(nèi)部混合相作用時產(chǎn)生了較大壓力，壓力極值作用點(diǎn)隨著閉合射流運(yùn)動而改變直至空泡完全潰滅。

圖8 空泡深閉合壓力分布及水蒸汽相云圖Fig.8 Contour of the pressure field and vapor fraction

4 結(jié)論

本文對半錐角為45°的錐頭圓柱體以500 m/s的初始速度自由垂直入水問題開展了數(shù)值模擬研究，對入水過程的速度、位移以及空泡形態(tài)進(jìn)行了相關(guān)分析，結(jié)果與文獻(xiàn)理論分析模型得到的結(jié)果對比吻合較好。同時，研究了入水深閉合的形成過程及流場特性，得到以下結(jié)論:

1)小型運(yùn)動體高速入水深閉合發(fā)生在面閉合之后，深閉合時運(yùn)動體行程較遠(yuǎn)，閉合位置較深。

2)在空泡內(nèi)部氣、水、水蒸汽混合介質(zhì)的影響下，高速入水深閉合多發(fā)生在一段區(qū)域內(nèi)，且歷時極短，深閉合后空泡迅速潰滅。

3)在空泡深閉合區(qū)域流場壓力值較大，壓力極值出現(xiàn)在閉合區(qū)域流體撞擊位置，隨著空泡逐漸分離，高壓區(qū)始終出現(xiàn)在水體與空泡分離點(diǎn)作用位置。

References)

［1］ Worthington A M，Cole R S. Impact with a liquid surface studied by the aid of instantaneous photography［J］. Philosophical Transactions of the Royal Society，1900，194(A):175 -200.

［2］ May A. Water entry and the cavity-running behavior of missiles，AD A020429［R］. Dayton:ASTIA，1975.

［3］ Gilbarg D，Anderson R A. Influence of atmospheric pressure on the phenomena accompanying the entry of spheres into water［J］.Journal of Applied Physics，1948，9(2):127 -139.

［4］ May A，Woodhull J C. Drag coefficients of steel spheres entering water vertically［J］. Journal of Applied Physics，1948，19(12):1109 -1121.

［5］ May A，Woodhull J C. The virtual mass of a sphere entering water vertically［J］. Journal of Applied Physics，1950，21 (12):1285 -1289.

［6］ May A. Effect of surface condition of a sphere on its water-entry cavity［J］. Journal of Applied Physics，1951，22(10):1219 -1222.

［7］ ZHAO R，F(xiàn)altinsen O. Water entry of two-dimensional bodies［J］. Journal of Fluid Mechanics，1993，246:593 -612.

［8］ Semenov A，Iafrati A. On the nonlinear water entry problem of asymmetric wedges［J］. Journal of Fluid Mechanics，2006，547:231 -256.

［9］ Sun S L，Wu G X. Oblique water entry of a cone by a fully threedimensional nonlinear method［J］. Journal of Fluids and Structures，2013，42:313 -332.

［10］ Lundstrom E A，F(xiàn)ung W K. Fluid dynamic analysis of hydraulic ram Ⅲ，AD 031644［R］. Dayton:ASTIA，1976.

［11］ Lee M，Longoria R G，D E Wilson. Cavity dynamics in highspeed water entry［J］. Physics of Fluids，1997，9(3):540 -550.

［12］ 王聰，何春濤，權(quán)曉波，等. 空氣壓強(qiáng)對垂直入水空泡影響的數(shù)值研究［J］. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2012，44(5):14 -19.WANG Cong，HE Chun-tao，QUAN Xiao-bo，et al. Numerical simulation of the influence of atmospheric pressure on water-cavity formed by cylinder with vertical water-entry［J］. Journal of Harbin Institute of Technology，2012，44(5):14 -19. (in Chinese)

［13］ 何春濤，王聰，閔景新，等. 回轉(zhuǎn)體勻速垂直入水早期空泡數(shù)值模擬研究［J］. 工程力學(xué)，2012，29(4):237 -243.HE Chun-tao，WANG Cong，MIN Jing-xin，et al. Numerical simulation of early air-cavity of cylinder cone with vertical waterentry［J］. Engineering Mechanics，2012，29(4):237 -243.(in Chinese)

［14］ Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications［J］. AIAA-Journal，1994，32 (8):1598 -1605.