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

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

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

0 引言

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

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

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

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

1 計算方法

1.1 空化多相流方程

混合相連續方程：

(1)

混合相動量方程：

(2)

蒸汽相輸運方程：

(3)

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

1.2 空化模型

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

(4)

式中：Fv為蒸發系數，Fv=50；rn為成核區的體積分數，rn=5×10-4；Rb為氣泡半徑，Rb=1 μm；pv為飽和蒸汽壓；Fc為凝結系數，Fc=0.01.

1.3 湍流模型

湍流模型對于空泡形態的演化及運動過程有較大影響，本文采用剪應力傳遞(SST)k-ω、標準(Standard)k-ε(ε為湍流動能耗散率)、重整化群(RNG)k-ε、可實現(Realizable)k-ε4種湍流模型開展對比分析研究。

1.4 可壓縮狀態方程

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

(5)

式中：pr和ρr為參考壓力和參考密度；K和n分別為溫度和壓力的弱函數，這里設為常數，K=3×108，n=7.

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

(6)

1.5 數值求解方法

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

圖1 計算模型Fig.1 Computational model

表1 模型參數
Tab.1 Model parameters

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

圖2為計算域示意圖，由于計算模型為旋成體外形，采用二維軸對稱網格。圖2中，上部為壓力進口邊界條件，設置為標準大氣壓力，底部為壓力出口邊界條件，按照公式p=p∞+ρlgh設置，本文中p∞為1個標準大氣壓，h表示水深，空氣域與液相域的內側為軸對稱邊界。

圖2 計算域示意圖Fig.2 Computational domain

2 結果分析

2.1 湍流模型比較分析

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

(7)

(8)

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

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

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

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

圖5 旋成體入水1 ms時數值模擬空泡形態及與理論解比較(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種湍流模型計算入水時間為1 ms時刻空泡界面及與理論解的比較圖。由圖5～圖8可見，采用Standardk-ε、Realizablek-ε、RNGk-ε3種湍流模型數值模擬的空泡界面在水氣界面處為向內收縮，采用SSTk-ω湍流模型數值模擬的空泡界面在自由面處為向外擴張，造成這一現象的原因是因為SSTk-ω湍流模型能夠更好地模擬空泡界面剪應力的傳遞作用，因此采用SSTk-ω湍流模型得到的空泡界面與理論值符合最好。為了驗證本文方法的準確性，進一步在高速入水理論解及實驗結果比對方面開展工作。

圖6 旋成體入水1 ms時數值模擬空泡形態及與理論解比較(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時數值模擬空泡形態及與理論解比較(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時數值模擬空泡形態及與理論解比較(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-ω湍流模型對入水速度200 m/s、400 m/s 2種情況開展數值模擬，圖9為入水衰減速度數值模擬與理論值的對比，圖10為無量綱入水深度數值模擬與理論值的對比。從圖9和圖10中可以看出，數值模擬與理論值具有很好的一致性，整體誤差在2%以內。

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

參考文獻[15，17]的實驗模型、實驗狀態、網格總數及設置，采用本文建立的數值模擬方法與該高速入水實驗進行對比研究，圖11所示為入水速度衰減數值模擬與實驗值的比較。由圖11可見，初始入水速度為141.15 m/s條件下，數值模擬與實驗結果一致性較好，整體誤差控制在1%以內。

圖12所示為數值模擬空泡形態發展與文獻[15，17]實驗結果對比，可見數值模擬基本能夠得到實驗中空泡發展的過程，在空泡形態上數值模擬與實驗結果基本一致。因此，從數值模擬與理論結果的比較以及數值模擬與參考實驗的比較來看，本文建立的數值方法是正確的。

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

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

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

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

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

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

圖11 旋成體入水速度數值模擬與實驗結果[15，17]對比(入水速度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 空泡形態數值模擬與實驗結果[15，17]對比(入水速度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 沖擊載荷隨時間變化(入水速度100 m/s)Fig.13 Impact load vs. time(water-entry velocity of 100 m/s)

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

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

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

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

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

圖17 旋成體入水0.5 ms時刻空泡形態圖(入水速度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時刻空泡形態圖(入水速度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時刻液體密度云圖(入水速度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 入水速度對運動過程的影響

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

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

3 結論

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

1)建立了高速入水問題的數值求解方法，采用SSTk-ω湍流模型進行數值模擬的結果在入水運動過程及空泡形態方面與理論解一致性最好，同時與參考實驗的結果開展了比較研究，驗證了本文方法的準確性。

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

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