基于ABAQUS顯式CEL方法的球體入水數值研究

0 引 言

入水過程是高速射彈、空投魚雷等武器從空中彈道進入水下彈道的一個重要的過渡環節。結構體從接觸自由水面瞬間，到完全進入水域內達到穩定狀態完成入水過程、形成入水空泡流動，是一個非常短暫的瞬態過程。這個短暫的入水過程因同時涉及結構體、空氣和水三相的相互作用，導致了流動介質突變引起的載荷突變、自由液面變化、空泡發展和結構體水下姿態變化、彈道發展等眾多的復雜問題。

入水問題的研究，對解決入水彈道穩定性、入水沖擊及緩沖，都有重要工程意義及理論研究價值。

1 入水及入水空泡發展

結構體以一定的初始速度穿越自由液面入水，其周圍將形成一個空氣層，伴隨著彈體下降至自由液面以下，空氣層進一步發展成為一個空腔，形成入水空泡。入水空泡從生成到潰滅呈現出復雜而有規律的發展過程，可以劃分為流動形成、空泡敞開、空泡閉合和空泡潰滅4個階段[1]。

以球體垂直入水過程為例，入水及入水空泡流動發展過程如圖1所示。

流動形成階段：從球體接觸水面開始，到形成穩定、有規律的流場這一過程。在流動形成階段，自由液面上方產生噴濺，同時球體尾部形成一個敞開的空泡，流動由此進入空泡敞開階段。

在空泡敞開階段，隨著球體的向下運動，球體入水深度不斷增大，球體前的流體質點被排擠到球體的外側流域，使入水空泡表現為：在被拉長的同時還伴隨著徑向的擴張。

球體進一步向下運動，隨著空泡繼續拉長，空泡內的壓強逐漸降低，周向擴張的流體在動能完全轉化為周向流體的壓力勢能后開始反向運動，空泡收縮并逐漸在水下軸線上的某一點閉合，這種閉合形式稱為深閉合。深閉合位置會伴有局部高壓，形成兩股反向的高速射流射入泡內，整個空泡被分為2個封閉區域并迅速分離，分別向自由液面和球面收縮并潰滅，至此完成球體入水。

圖1 球體垂直入水空泡流動發展Fig.1 Evolution of water-entry cavity during the process of sphere’s water-entry

2 有限元顯式動力學方法在水動力學中的應用

隨著計算機技術的發展，如Ansys LS-DYNA、Abaqus/Explicit等，有限元顯式動力學數值求解技術的仿真程序被開始應用于出、入水等水動力學問題的求解。

有限元的顯式動態求解方法，作為隱式求解方法的補充，最初是為了模擬高速沖擊等高速動力學問題而逐步發展而來的，特別適用于短暫、瞬時的，包含接觸、碰撞行為的非線性動力學問題的模擬。通過引入用于描述水動力學行為的狀態方程，使出、入水沖擊過程，也可通過顯式動力學方法進行求解。

2.1 AbaqusExplicit水動力學求解方法

大型通用有限元軟件Abaqus中，AbaqusExplicit為其顯式動態分析模塊，包含用于水動力學模擬的狀態方程；應用其中的耦合歐拉-拉格朗日分析（Coupled Eulerian-Lagrangian analysis，CEL）功能，即可實現相關的水動力學問題求解。

傳統有限元的拉格朗日分析（Lagrangian analysis），材料屬性與網格節點相關聯，材料伴隨著網格變形而發生形變，當解決極端變形的情況時，會由于單元的變形扭曲而失去原有的精度。而歐拉分析（Eulerian analysis）方法，其網格節點空間固定，材料不隨單元變形而是在單元間流動，可有效地解決極端變形以及包含流體流動的問題，所以諸如液體晃動、氣體流動、穿透問題等均可通過歐拉分析有效處理。

但在如本文球體入水等流固耦合（fluid-structure interaction，FSI）分析中，歐拉分析雖然可有效處理流體流動分析，但在捕捉結構流固交界面上存在一定困難。此時，可應用耦合歐拉-拉格朗日分析（CEL）功能進行求解[2]。

Abaqus CEL方法中，流固交界面的捕捉，是通過歐拉-拉格朗日接觸的定義與求解實現的。CEL方法中歐拉-拉格朗日接觸，是Abaqus/Explicit顯式通用接觸的功能擴展，它基于浸沒邊界法（immersed boundary method）原理，建立了拉格朗日體在歐拉體中運動及流固耦合接觸分析的功能，從而使歐拉單元與拉格朗日單元相關聯，并進而進行耦合求解。

2.2 AbaqusExplicit水動力學模型

本文入水過程模擬中，水介質流動視為近似不可壓縮的、粘性層流流動。采用線性Us-UpHugoniot形式的Mie-Grüneisen狀態方程描述水介質的體積響應。Mie-Grüneisen狀態方程的常用形式為：

式中：pH為雨貢紐壓力，EH為雨貢紐比能，均僅為密度的函數；Γ=Γ0ρ0/ρ，Γ0為材料常數，ρ0為參考密度。并且有如下關系為名義體積壓縮應變。

進而得到：

常用的對材料雨貢紐曲線擬合關系為：

式中：c0、s為定義線性沖擊波波速Us、粒子速度Up關系的系數，其關系如下：

進而得到線性Us-UpHugoniot形式的Mie-Grüneisen狀態方程，表述為：

通過以上狀態方程，即可定義水介質的體積強度及靜水狀態體積響應特性。同時，假設水介質的剪切響應與體積響應是非耦合的，采用牛頓粘性剪切模型，定義其粘性剪切響應。

上式中：S 為應力偏量，e˙為應變率偏量，μ 為粘性，同時，γ˙=2e˙為工程偏應變率。

綜上所述，通過（5）、（6）式，水動力狀態方程及粘性剪切模型，即完成水介質材料屬性的定義。真實水介質線性Us-UpMie-Grüneisen狀態方程的基本參數為：密度 ρ=9.832×102kg/m3，c0=1.450 6×103m/s，s=0，Γ0=0，μ=1×10-3Pa·s。

3 球體入水的數值模擬與實驗對比

Abaqus顯式CEL方法，是基于有限元平臺的水動力及流固耦合問題的數值求解方法。相對于CFD的流場求解功能，CEL方法由于對N-S方程進行了簡化，不能考慮流動邊界層效應，因此，CEL方法側重于入水空泡的基本形態及演化，以及水彈道、水動力載荷方面的應用。

本文通過數值模擬與球體垂直入水實驗[3]的比對，驗證數值模型的有效性。球體垂直入水實驗中，采用的標準臺球直徑為57.15 mm，質量m=0.17 kg，球體以5.1 m/s的觸水速度垂直入水。

采用Abaqus動態顯式CEL方法，進行入水的流固耦合求解，球體采用Lagrange離散剛體（Discrete Rigid）單元，空氣域、水域采用Euler單元，采用歐拉-拉格朗日接觸算法直接耦合結構網格和流體網格，進行入水模擬。

值得注意的是，本文低速入水問題，水、氣交界自由面為自由邊界。不同于致密固體和液體的體積壓縮，相對無約束狀態下水介質由于自由面變形，可壓縮性增加，體積性能會較實際偏軟。此時，在無約束狀態下，作為不可壓縮約束的罰參數的體積模量，通常取值比實際體積模量小2-3個量級，以表征體積模量的這種軟化特性[4]。

因此，水介質線性 Us-Up狀態方程參數為：密度 ρ=9.832×102kg/m3，c0=45.85 m/s，s=0，Γ0=0，μ=1×Pa·s，由體積模量公式 K=對應聲速下的體積模量為2.07 MPa，比實際體積模量2.07 GPa小3個量級。同時，對水介質施加重力載荷，并定義初始的壓應力場來模擬初始的水深壓力分布。

基于以上原則進行數值計算得到入水過程的入水空泡形態演化，并與實驗比對如圖2、3所示。

圖2 入水實驗入水空泡的演化Fig.2 Evolution of water-entry cavity in experiments

圖3 數值模擬入水空泡的演化Fig.3 Evolution of water-entry cavity in numerical simulation

比較表明，CEL方法捕捉到球體入水空泡演化的基本特征，入水空泡的演化過程與實驗現象基本一致，定性地表明了CEL方法進行入水模擬的有效性。同時，數值模擬也較清晰地捕捉到球體入水空泡深閉合階段，入水空泡收縮、拉斷以及伴隨形成的分別指向自由面與球體的反向回射流等現象，如圖4所示。

以自由液面為零點，豎直向下為位移的正方向，球體相對于自由液面入水的位移曲線，與實驗數據比對，如圖5所示；球體入水的速度衰減曲線，與實驗數據比對，如圖6所示。

通過實驗數據與數值模擬入水彈道的一致性校準結果，表明CEL方法對定量求解入水問題水彈道的有效性。

圖4 數值模擬入水空泡的深閉合及回射流形成Fig.4 Water-entry cavity’s deep closure and jets after closure

圖5 球體垂直入水位移曲線Fig.5 Displacement-time curve of sphere’s water-entry

圖6 球體垂直入水速度衰減曲線Fig.6 Velocity-time curve of sphere’s water-entry

值得注意的是，球體垂直入水速度衰減的實驗數據，是由光學觀測的位移數據經判讀得到的，因而實驗的速度衰減數據更多地表現為位移數據點間的平均特性。

數值模擬則更加清晰地獲得了球體入水彈道的非線性特性，主要表現為：入水初期速度非線性快速衰減；此后直至空泡敞開階段結束，速度近似呈線性衰減；至深閉合階段入水空泡頸縮，球體速度再次加速衰減；而入水空泡拉斷后，形成的指向球體的回射流作用下，又使球體瞬間加速；最后球體后端空泡不斷潰滅，球體速度呈現出在波動中衰減的特征。

4 球體入水沖擊載荷

顯式動力學方法通過小時間增量以獲得高精度解，特別適用于高速動力學問題，因而對于捕捉入水沖擊脈沖載荷也是合適的。

4.1 入水沖擊載荷的采樣頻率

為精確地捕捉入水沖擊脈沖載荷峰值，以顯式動力學方法時間增量的穩定性極限Δtstable為參考基準，作為輸出時間間隔，捕捉入水沖擊載荷峰值。

顯式方法的穩定極限用單元長度Le、材料波速cd定義為：

作為參考基準的穩定性極限Δtstable與網格尺寸直接相關，即入水沖擊脈沖載荷峰值的求解精度，直接與空間網格離散精度相關。

網格離散精度驗證中，隨著網格尺寸減小，穩定性極限Δtstable也隨著減小，載荷輸出頻率相應地增加。選用Δtstable量級的輸出時間間隔，便使得入水沖擊峰值的捕捉，僅依賴于數值求解精度。

4.2 歐拉-拉格朗日接觸的網格離散

Abaqus CEL方法中歐拉-拉格朗日顯式通用接觸模擬中，同樣需要嚴格滿足接觸定義及離散的相關要求。

Abaqus/Explicit中，剛體單元不進行光滑處理，它們精確地保持由用戶定義的表面形狀，因而本文中剛體單元定義的球體表面，需要較細致的網格密度進行幾何離散。同時，球體入水模擬中，剛體球表面是作為接觸的主控面，與整個歐拉域建立自接觸來模擬自由液面變形、入水空泡演化的。因此，為防止節點穿透，提高模擬精度，歐拉域的網格離散密度應大于接觸的主控面-剛體球表面的離散密度。

此前的研究表明，權衡計算精度與計算資源消耗，球體離散剛體網格尺寸與歐拉域網格單元尺寸之比 SR，可取為 SR=1.5～2[5]。

4.3 CEL方法入水沖擊載荷預報

基于以上分析，本文給出了CEL方法入水沖擊載荷精確性校核的基本原則：在對入水模擬進行入水彈道、特別是入水初期彈道校準的基礎上，以入水結構體剛體離散單元為基準進行網格劃分，入水剛體離散單元與歐拉域水介質網格尺寸比例始終保持SR=1.5～2；載荷輸出時間間隔采用穩定極限Δtstable相同的量級；最后，針對入水沖擊載荷峰值，進行網格離散精度的驗證。

圖7 不同球體表面離散密度對比Fig.7 Comparison of different mesh density

圖8 球體入水沖擊載荷對比Fig.8 Comparison of water-entry impact force

不同網格離散密度下，入水沖擊載荷曲線峰值對比如圖7、圖8所示。

如圖8所示，在不同的網格離散密度下，入水沖擊載荷曲線的峰值相對誤差控制在10%以內，即可認為達到入水沖擊載荷的求解精度。

5 結 語

本文針對球體入水過程，采用Abaqus顯式動力學CEL方法，完成了球體入水空泡演化、發展過程的數值模擬。通過與實驗入水彈道的對比，表明CEL方法在數值求解流固耦合問題中的有效性。在入水初期彈道校準的基礎上，提出了在數值模擬入水沖擊載荷中，求解所需的采樣頻率、歐拉-拉格朗日接觸網格離散的要求，并完成球體入水沖擊載荷的預報。

參 考 文 獻：

[1]陳九錫,張開榮譯.水彈道學模擬[M].北京：國防工業出版社,1979.Chen Jiuxi,Zhang Kairong.Hydroballistics modeling[M].Beijing：Natoinal Defense Industry Press,1979.

[2]張偉偉,金先龍.球體撞擊自由液面相關效應的數值模擬方法研究[J].船舶力學,2014,18(1-2)：28-36.Zhang Weiwei,Jin Xianlong.Numerical methods for simulating the related effects of sphere impact onto free surface[J].Journal of Ship Mechanics,2014,18(1-2)：28-36.

[3]馬慶鵬,何春濤,等.球體垂直入水空泡實驗研究[J].爆炸與沖擊,2014,3(34)：174-177.Ma Qingpeng,He Chuntao.Experimental investigation on vertical water-entry cavity of sphere[J].Explosion and Shock Waves,2014,3(34)：174-177.

[4]Abaqus 6.14/Abaqus benchmarks Guide/Analysis Tests/Adaptivity/Water sloshing in a piching tank[CP].

[5]王永虎,石秀華.空投雷彈斜入水初始彈道數值分析[J].彈道學報,2012,6(24)：93-94.Wang Yonghu,Shi Xiuhua.Numerical analysis for initial hydroballistics of airborne missile during oblique water-entry impact[J].Journal of Ballistics,2012,6(24)：93-94.