旋成體高速入水可壓縮性影響研究

李國良， 尤天慶， 孔德才， 李靜， 周偉江

(1.中國航天空氣動(dòng)力技術(shù)研究院 一所， 北京 100074；2.北京宇航系統(tǒng)工程研究所， 北京 100076)

0 引言

航行體入水運(yùn)動(dòng)是一個(gè)具有強(qiáng)非定常特性的復(fù)雜多相流動(dòng)過程，這一過程具有強(qiáng)瞬時(shí)及高載荷特性。航行體在入水過程中會(huì)形成一個(gè)由空氣和空化產(chǎn)生的蒸汽混合空腔，形成過程涵蓋多相之間的摻混、液體可壓縮、自由面與空泡間相互作用等復(fù)雜物理問題。現(xiàn)階段航行體高速入水問題研究更具有工程意義，它涉及武器裝備如空投魚雷、超空泡射彈、多次出入水導(dǎo)彈等的關(guān)鍵技術(shù)問題。

航行體入水問題的實(shí)驗(yàn)研究國外開展得較早[1-5]，國內(nèi)時(shí)素果等[6]采用高速攝影方法對(duì)細(xì)長體以速度60 m/s入水過程開展實(shí)驗(yàn)研究，觀察入水過程空泡形態(tài)演化規(guī)律，采用壓力傳感器獲得了物體表面壓力變化規(guī)律。蔣運(yùn)華等[7]采用實(shí)驗(yàn)方法，在實(shí)驗(yàn)水箱中開展約束和不約束運(yùn)動(dòng)體姿態(tài)的入水實(shí)驗(yàn)，分析了有無擾動(dòng)下空泡尺寸與弗勞德數(shù)和歐拉數(shù)之間的關(guān)系。何春濤等[8]通過垂直和傾斜兩種入水方式開展實(shí)驗(yàn)研究，分析了入水過程中的一系列流動(dòng)現(xiàn)象，以及空泡的生成、發(fā)展和閉合，并且開展了串聯(lián)和并聯(lián)入水時(shí)空泡的生成和相互之間的影響。

由于實(shí)驗(yàn)研究成本較高及周期偏長的原因，入水問題研究的早期更多采用數(shù)值仿真手段。夏維學(xué)等[9]采用流體體積函數(shù)(VOF)多相流模型、k-ω湍流模型(k為湍流動(dòng)能，ω為比耗散率)及動(dòng)網(wǎng)格技術(shù)，數(shù)值模擬了旋轉(zhuǎn)角速度為199 rad/s、入水速度為5.45 m/s條件下入水空泡的形成、發(fā)展、閉合和潰滅的演化過程。孫釗等[10]采用VOF多相流模型和連續(xù)表面力(CSF)模型，對(duì)親水性及疏水性球體垂直入水過程開展了數(shù)值模擬研究，分析了親水性及疏水性球體低速入水對(duì)空化初生以及流體動(dòng)力的影響。馬慶鵬等[11]基于有限體積法和VOF方法，對(duì)小型錐頭圓柱體高速垂直入水開展數(shù)值模擬，研究了入水空泡演化及深閉合過程空泡流場的流動(dòng)特性和壓力分布等問題。Marrone等[12]采用光滑粒子流體動(dòng)力學(xué)(SPH)模型對(duì)航空器和直升機(jī)水上迫降涉及的機(jī)理問題開展研究，把飛行器簡化為平板，對(duì)其在低速小攻角條件下水擊載荷進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)數(shù)據(jù)與勢流理論進(jìn)行了對(duì)比。Sean等[13]采用浸沒邊界法建立了多相流數(shù)值求解方法，克服了水、空氣、蒸汽多相之間相界面大密度比的難題，并且模擬了二維方柱和圓柱低速入水問題。Neaves等[14]忽略湍流對(duì)空泡發(fā)展的影響采用預(yù)處理方法和低耗散迎風(fēng)格式，對(duì)二維旋成體入水問題開展了數(shù)值模擬，驗(yàn)證了該方法的適用性。陳晨等[15]采用數(shù)值模擬方法研究了小型運(yùn)動(dòng)體高速傾斜入水時(shí)的空泡流動(dòng)特性，并深入研究了不同入水角度對(duì)空泡形態(tài)、力學(xué)特性的影響。

目前航行體入水問題研究通常不考慮液體的可壓縮性。本文基于壓力耦合方程組的半隱式(SIMPLE)算法，在考慮液體可壓縮性條件下，對(duì)錐柱旋成體以50 m/s、100 m/s、200 m/s、400 m/s、800 m/s速度入水的過程進(jìn)行數(shù)值仿真，研究湍流模型對(duì)入水過程的影響、液體可壓縮性對(duì)入水沖擊載荷及空泡形態(tài)的影響，以及入水速度對(duì)入水后運(yùn)動(dòng)特性的影響。

1 計(jì)算方法

1.1 空化多相流方程

混合相連續(xù)方程：

(1)

混合相動(dòng)量方程：

(2)

蒸汽相輸運(yùn)方程：

(3)

式中：ρm為混合相密度，ρm=αlρl+αvρv+αgρg，αl、αv、αg分別為液相、蒸汽相、氣相的體積分?jǐn)?shù)，ρl、ρv、ρg分別為液相、蒸汽相、氣相的密度；t為時(shí)間；u為速度，x為坐標(biāo)，下標(biāo)i和j分別為1，2，3，表示3個(gè)坐標(biāo)方向；p為壓力；μm為層流混合相黏性系數(shù)，μm=αlμl+αvμv+αgμg，μl、μv、μg分別為液相、蒸汽相、氣相的黏性系數(shù)；Re、Rc為空化源項(xiàng)；μt為混合相湍流黏性系數(shù)；g為重力加速度。

1.2 空化模型

采用Zwart-Gerber-Belamri空化模型：

(4)

式中：Fv為蒸發(fā)系數(shù)，F(xiàn)v=50；rn為成核區(qū)的體積分?jǐn)?shù)，rn=5×10-4；Rb為氣泡半徑，Rb=1 μm；pv為飽和蒸汽壓；Fc為凝結(jié)系數(shù)，F(xiàn)c=0.01.

1.3 湍流模型

湍流模型對(duì)于空泡形態(tài)的演化及運(yùn)動(dòng)過程有較大影響，本文采用剪應(yīng)力傳遞(SST)k-ω、標(biāo)準(zhǔn)(Standard)k-ε(ε為湍流動(dòng)能耗散率)、重整化群(RNG)k-ε、可實(shí)現(xiàn)(Realizable)k-ε4種湍流模型開展對(duì)比分析研究。

1.4 可壓縮狀態(tài)方程

對(duì)于高速入水，入水載荷及空泡界面主要受液相影響，因此本文中氣相和蒸汽相仍按照不可壓縮來處理。液相要考慮其可壓縮性，在實(shí)現(xiàn)液相可壓縮性時(shí)應(yīng)用Tait狀態(tài)方程。不含溫度修正的Tait狀態(tài)方程為

(5)

式中：pr和ρr為參考?jí)毫蛥⒖济芏?；K和n分別為溫度和壓力的弱函數(shù)，這里設(shè)為常數(shù)，K=3×108，n=7.

利用Tait方程定義流體密度變化的同時(shí)，還要定義流體中的聲速c方程：

(6)

1.5 數(shù)值求解方法

流場采用基于壓力的SIMPLE算法進(jìn)行求解。動(dòng)量方程中的對(duì)流項(xiàng)采用2階迎風(fēng)格式，壓力項(xiàng)離散采用PRESTO差分格式。多相流模型采用VOF+Level Set方法，時(shí)間項(xiàng)采用顯示格式。本文主要模擬高速入水問題，在入水過程中會(huì)產(chǎn)生空泡，需要多相流與動(dòng)網(wǎng)格耦合求解。動(dòng)網(wǎng)格計(jì)算采用動(dòng)態(tài)層法，為保證射彈在運(yùn)動(dòng)過程中空泡產(chǎn)生的區(qū)域網(wǎng)格狀態(tài)始終不變，采用全域網(wǎng)格隨射彈運(yùn)動(dòng)的動(dòng)網(wǎng)格方法，網(wǎng)格更新分別在流場入口及出口邊界處，在流場入口與出口新生網(wǎng)格高度始終與原網(wǎng)格高度相同。計(jì)算模型如圖1所示，模型采用錐柱旋成體，具體參數(shù)如表1所示，其中L為模型長度，D為模型直徑，θ為前部圓錐體錐角，m為模型質(zhì)量。

圖1 計(jì)算模型Fig.1 Computational model

表1 模型參數(shù)
Tab.1 Model parameters

L/mmD/mmθ/(°)m/kg50101270.01

圖2為計(jì)算域示意圖，由于計(jì)算模型為旋成體外形，采用二維軸對(duì)稱網(wǎng)格。圖2中，上部為壓力進(jìn)口邊界條件，設(shè)置為標(biāo)準(zhǔn)大氣壓力，底部為壓力出口邊界條件，按照公式p=p∞+ρlgh設(shè)置，本文中p∞為1個(gè)標(biāo)準(zhǔn)大氣壓，h表示水深，空氣域與液相域的內(nèi)側(cè)為軸對(duì)稱邊界。

圖2 計(jì)算域示意圖Fig.2 Computational domain

2 結(jié)果分析

2.1 湍流模型比較分析

湍流采用SSTk-ω、Standardk-ε、RNGk-ε、Realizablek-ε4種模型，比較湍流模型對(duì)入水速度衰減及入水深度的影響。根據(jù)文獻(xiàn)[16]提出的入水速度和空泡半徑隨時(shí)間變化的理論公式，得到瞬時(shí)加速度和空泡半徑分別為

(7)

(8)

式中：vp為瞬時(shí)速度；A0為模型橫截面；Cd為阻力系數(shù)，Cd=Cd0+σ,σ為空化數(shù)，Cd0為σ=0時(shí)的阻力系數(shù)，本文中Cd0=0.498；r(x,t)為空泡半徑；tb為前一運(yùn)動(dòng)時(shí)刻；r0為tb時(shí)刻空泡半徑；A(x)、B(x)的表達(dá)式見參考文獻(xiàn)[16]。

圖3為旋成體入水時(shí)間為0.006 s時(shí)的速度衰減曲線，可以看出4種湍流模型的數(shù)值總體趨勢上與理論解一致，但是在量值上有差異。在0.006 s時(shí)刻，理論解瞬時(shí)速度為39.08 m/s，Standardk-ε湍流模型模擬瞬時(shí)速度為33.57 m/s，SSTk-ω湍流模型模擬瞬時(shí)速度為38.45 m/s，RNGk-ε湍流模型模擬瞬時(shí)速度為37.24 m/s，Realizablek-ε湍流模型模擬瞬時(shí)速度為37.92 m/s. 從數(shù)值比較來看，使用SSTk-ω湍流模型得到的速度衰減，無論從整體趨勢上還是從瞬時(shí)點(diǎn)上都與理論解誤差最小。

圖3 旋成體入水速度隨時(shí)間衰減曲線(入水速度100 m/s)Fig.3 Velocity during water-entry vs. time(water-entry velocity of 100 m/s)

圖4 無量綱入水深度隨時(shí)間變化(入水速度100 m/s)Fig.4 Nondimensional water-entry depth vs. time (water-entry velocity of 100 m/s)

圖5 旋成體入水1 ms時(shí)數(shù)值模擬空泡形態(tài)及與理論解比較(Standard k-ε模型，入水速度100 m/s)Fig.5 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1ms of water-entry(Standard k-ε model，water-entry velocity of 100 m/s)

圖5～圖8分別為采用4種湍流模型計(jì)算入水時(shí)間為1 ms時(shí)刻空泡界面及與理論解的比較圖。由圖5～圖8可見，采用Standardk-ε、Realizablek-ε、RNGk-ε3種湍流模型數(shù)值模擬的空泡界面在水氣界面處為向內(nèi)收縮，采用SSTk-ω湍流模型數(shù)值模擬的空泡界面在自由面處為向外擴(kuò)張，造成這一現(xiàn)象的原因是因?yàn)镾STk-ω湍流模型能夠更好地模擬空泡界面剪應(yīng)力的傳遞作用，因此采用SSTk-ω湍流模型得到的空泡界面與理論值符合最好。為了驗(yàn)證本文方法的準(zhǔn)確性，進(jìn)一步在高速入水理論解及實(shí)驗(yàn)結(jié)果比對(duì)方面開展工作。

圖6 旋成體入水1 ms時(shí)數(shù)值模擬空泡形態(tài)及與理論解比較(Realizable k-ε模型，入水速度100 m/s)Fig.6 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1 ms of water-entry(Realizable k-ε model，water-entry velocity of 100 m/s)

圖7 旋成體入水1 ms時(shí)數(shù)值模擬空泡形態(tài)及與理論解比較(RNG k-ε模型，入水速度100 m/s)Fig.7 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1 ms of water-entry (RNG k-ε model，water-entry velocity at 100 m/s)

圖8 旋成體入水1 ms時(shí)數(shù)值模擬空泡形態(tài)及與理論解比較(SST k-ω模型，入水速度100 m/s)Fig.8 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1ms of water-entry (SST k-ω model，water-entry velocity of 100 m/s)

采用SSTk-ω湍流模型對(duì)入水速度200 m/s、400 m/s 2種情況開展數(shù)值模擬，圖9為入水衰減速度數(shù)值模擬與理論值的對(duì)比，圖10為無量綱入水深度數(shù)值模擬與理論值的對(duì)比。從圖9和圖10中可以看出，數(shù)值模擬與理論值具有很好的一致性，整體誤差在2%以內(nèi)。

圖9 旋成體入水速度隨時(shí)間衰減曲線Fig.9 Velocity during water-entry vs. time

參考文獻(xiàn)[15，17]的實(shí)驗(yàn)?zāi)Ｐ汀?shí)驗(yàn)狀態(tài)、網(wǎng)格總數(shù)及設(shè)置，采用本文建立的數(shù)值模擬方法與該高速入水實(shí)驗(yàn)進(jìn)行對(duì)比研究，圖11所示為入水速度衰減數(shù)值模擬與實(shí)驗(yàn)值的比較。由圖11可見，初始入水速度為141.15 m/s條件下，數(shù)值模擬與實(shí)驗(yàn)結(jié)果一致性較好，整體誤差控制在1%以內(nèi)。

圖12所示為數(shù)值模擬空泡形態(tài)發(fā)展與文獻(xiàn)[15，17]實(shí)驗(yàn)結(jié)果對(duì)比，可見數(shù)值模擬基本能夠得到實(shí)驗(yàn)中空泡發(fā)展的過程，在空泡形態(tài)上數(shù)值模擬與實(shí)驗(yàn)結(jié)果基本一致。因此，從數(shù)值模擬與理論結(jié)果的比較以及數(shù)值模擬與參考實(shí)驗(yàn)的比較來看，本文建立的數(shù)值方法是正確的。

2.2 流體可壓縮性的影響分析

高速入水會(huì)產(chǎn)生較大的沖擊載荷，液體可壓縮性對(duì)沖擊載荷有較大影響，在入水速度分別為50 m/s、 100 m/s、200 m/s、400 m/s、800 m/s 5種條件下，分析考慮液體可壓縮和不可壓縮兩種處理方式對(duì)入水力學(xué)特性的影響。表2為數(shù)值模擬結(jié)果。

根據(jù)表2的結(jié)果可知：在入水速度為50 m/s、100 m/s時(shí)，兩種處理方式最大載荷出現(xiàn)的時(shí)間相同，沖擊載荷基本相當(dāng)；在入水速度為200 m/s時(shí)，不可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水1.6×10-5s時(shí)，峰值為3 765 N，可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水1.7×10-5s時(shí)，峰值為3 474 N，二者相差7.7%，最大密度提高2%；在入水速度為400 m/s時(shí)，不可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水0.95×10-5s時(shí)，峰值為13 323 N，可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水1.1×10-5s時(shí)，峰值為12 211 N，二者相差9.1%，最大密度提高6.5%；在入水速度為800 m/s時(shí)，不可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水0.55×10-5s時(shí)，峰值為55 478 N，可壓縮流體入水沖擊載荷最大值出現(xiàn)在入水0.875×10-5s時(shí)，峰值為37 718 N，二者相差47.1%，最大密度提高32.0%.

圖13為入水速度為100 m/s時(shí)兩種處理方式?jīng)_擊載荷的變化歷程，圖中Fm為沖擊載荷；圖14為沖擊載荷最大時(shí)密度等值線圖；圖15為入水速度為800 m/s時(shí)兩種處理方式?jīng)_擊載荷的變化歷程；圖16為沖擊載荷最大時(shí)密度等值線圖。由圖13～圖16可知，入水速度在100 m/s時(shí)液體可壓縮性對(duì)入水沖擊載荷基本無影響，入水速度≥200 m/s時(shí)，液體可壓縮性對(duì)沖擊載荷有影響，沖擊載荷會(huì)變小，而且最大沖擊載荷出現(xiàn)的時(shí)間會(huì)滯后，隨著入水速度不斷增加，滯后時(shí)間會(huì)加大。

圖17、圖18分別為入水速度800 m/s兩種液體處理方式入水0.5 ms時(shí)刻的空泡形態(tài)圖，其中圖17為按照不可壓縮方式數(shù)值模擬的空泡形態(tài)圖，圖18為按照可壓縮方式數(shù)值模擬的空泡形態(tài)圖。為了比較液體可壓縮性對(duì)空泡形態(tài)的影響，對(duì)空泡輪廓線提取出來以作比較，如圖19所示，可見考慮液體可壓縮特性時(shí)空泡會(huì)有向內(nèi)收縮現(xiàn)象。從圖20可以看出，考慮液體可壓縮性時(shí)，在空泡界面附近水的密度降低，在壓差作用下周圍液體向空泡內(nèi)部擠壓，造成空泡收縮現(xiàn)象。

圖10 旋成體無量綱入水深度隨時(shí)間變化曲線Fig.10 Nondimensional depth during water-entry vs. time

圖11 旋成體入水速度數(shù)值模擬與實(shí)驗(yàn)結(jié)果[15，17]對(duì)比(入水速度141.15 m/s)Fig.11 Numerically simulated results vs. experimental results[15,17] of velocity(water-entry velocity of 141.15 m/s)

圖12 空泡形態(tài)數(shù)值模擬與實(shí)驗(yàn)結(jié)果[15，17]對(duì)比(入水速度141.15 m/s)Fig.12 Numerically simulated results vs. experimental results in Refs.[15] and [17] of cavitation configuration(water-entry velocity of 141.15 m/s)

圖13 沖擊載荷隨時(shí)間變化(入水速度100 m/s)Fig.13 Impact load vs. time(water-entry velocity of 100 m/s)

圖14 最大載荷時(shí)刻密度云圖(入水速度100 m/s)Fig.14 Density contour under maximun impact load(water-entry velocity of 100 m/s)

圖15 沖擊載荷隨時(shí)間變化(入水速度800 m/s)Fig.15 Impact load vs. time(water-entry velocity of 800 m/s)

圖16 最大載荷時(shí)刻密度云圖(入水速度800 m/s)Fig.16 Density contour under maximun impact load(water-entry velocity of 800 m/s)

表2 旋成體不同入水速度沖擊載荷數(shù)值模擬
Tab.2 Numerically simulated results of water-entry impact load at different water-entry velocities

入水速度/(m·s-1)最大載荷發(fā)生時(shí)間/10-5s最大載荷/N最大密度/(kg·m-3)不可壓縮可壓縮不可壓縮可壓縮不可壓縮可壓縮504.44.431531410001001100 2.42.4853851100010042001.61.737653474100010204000.951.101332312211100010658000.5500.875554783771810001320

圖17 旋成體入水0.5 ms時(shí)刻空泡形態(tài)圖(入水速度800 m/s，不可壓縮液體)Fig.17 Cavitation configuration of revolutionary body at 0.5 ms (water-entry velocity of 800 m/s, incompressible water)

圖18 旋成體入水0.5 ms時(shí)刻空泡形態(tài)圖(入水速度800 m/s，可壓縮液體)Fig.18 Cavitation configuration of revolutionary body at 0.5 ms(water-entry velocity of 800 m/s, compressible water)

圖19 空泡輪廓圖的比較(入水速度800 m/s)Fig.19 Comparison of cavitation outlines (water-entry velocity of 800 m/s)

圖20 旋成體入水0.5 ms時(shí)刻液體密度云圖(入水速度800 m/s，可壓縮流體)Fig.20 Density contour of water at 0.5 ms during water entry of revolutionary body (water-entry velocity of 800 m/s, compressible water)

2.3 入水速度對(duì)運(yùn)動(dòng)過程的影響

圖21 旋成體無量綱入水速度隨時(shí)間變化Fig.21 Nondimensional velocity during water-entry vs. time

圖22 旋成體入水加速度隨時(shí)間變化Fig.22 Acceleration during water-entry vs. time

3 結(jié)論

本文基于SIMPLE算法，在考慮液體可壓縮性的條件下，對(duì)錐柱旋成體在50 m/s、100 m/s、200 m/s、400 m/s、800 m/s入水速度下開展數(shù)值模擬研究。得出以下主要結(jié)論：

1)建立了高速入水問題的數(shù)值求解方法，采用SSTk-ω湍流模型進(jìn)行數(shù)值模擬的結(jié)果在入水運(yùn)動(dòng)過程及空泡形態(tài)方面與理論解一致性最好，同時(shí)與參考實(shí)驗(yàn)的結(jié)果開展了比較研究，驗(yàn)證了本文方法的準(zhǔn)確性。

2)在入水速度≤100 m/s時(shí)，液體可壓縮性對(duì)沖擊載荷基本沒有影響，在入水速度≥200 m/s時(shí)，需要考慮液體可壓縮性；隨著入水速度增加，液體可壓縮性對(duì)入水沖擊載荷影響越大，其會(huì)弱化入水沖擊載荷及延緩最大載荷出現(xiàn)的時(shí)間。

3)入水速度越大，入水過程速度衰減越快，加速度值在入水初期較大；在計(jì)算模型周圍被超空泡包裹的航行階段，隨著入水深度增加，加速度逐漸減小，而且加速度的變化逐漸平緩。