基于ABAQUS顯式CEL方法的球體入水?dāng)?shù)值研究

0 引 言

入水過程是高速射彈、空投魚雷等武器從空中彈道進(jìn)入水下彈道的一個(gè)重要的過渡環(huán)節(jié)。結(jié)構(gòu)體從接觸自由水面瞬間，到完全進(jìn)入水域內(nèi)達(dá)到穩(wěn)定狀態(tài)完成入水過程、形成入水空泡流動(dòng)，是一個(gè)非常短暫的瞬態(tài)過程。這個(gè)短暫的入水過程因同時(shí)涉及結(jié)構(gòu)體、空氣和水三相的相互作用，導(dǎo)致了流動(dòng)介質(zhì)突變引起的載荷突變、自由液面變化、空泡發(fā)展和結(jié)構(gòu)體水下姿態(tài)變化、彈道發(fā)展等眾多的復(fù)雜問題。

入水問題的研究，對(duì)解決入水彈道穩(wěn)定性、入水沖擊及緩沖，都有重要工程意義及理論研究價(jià)值。

1 入水及入水空泡發(fā)展

結(jié)構(gòu)體以一定的初始速度穿越自由液面入水，其周圍將形成一個(gè)空氣層，伴隨著彈體下降至自由液面以下，空氣層進(jìn)一步發(fā)展成為一個(gè)空腔，形成入水空泡。入水空泡從生成到潰滅呈現(xiàn)出復(fù)雜而有規(guī)律的發(fā)展過程，可以劃分為流動(dòng)形成、空泡敞開、空泡閉合和空泡潰滅4個(gè)階段[1]。

以球體垂直入水過程為例，入水及入水空泡流動(dòng)發(fā)展過程如圖1所示。

流動(dòng)形成階段：從球體接觸水面開始，到形成穩(wěn)定、有規(guī)律的流場(chǎng)這一過程。在流動(dòng)形成階段，自由液面上方產(chǎn)生噴濺，同時(shí)球體尾部形成一個(gè)敞開的空泡，流動(dòng)由此進(jìn)入空泡敞開階段。

在空泡敞開階段，隨著球體的向下運(yùn)動(dòng)，球體入水深度不斷增大，球體前的流體質(zhì)點(diǎn)被排擠到球體的外側(cè)流域，使入水空泡表現(xiàn)為：在被拉長的同時(shí)還伴隨著徑向的擴(kuò)張。

球體進(jìn)一步向下運(yùn)動(dòng)，隨著空泡繼續(xù)拉長，空泡內(nèi)的壓強(qiáng)逐漸降低，周向擴(kuò)張的流體在動(dòng)能完全轉(zhuǎn)化為周向流體的壓力勢(shì)能后開始反向運(yùn)動(dòng)，空泡收縮并逐漸在水下軸線上的某一點(diǎn)閉合，這種閉合形式稱為深閉合。深閉合位置會(huì)伴有局部高壓，形成兩股反向的高速射流射入泡內(nèi)，整個(gè)空泡被分為2個(gè)封閉區(qū)域并迅速分離，分別向自由液面和球面收縮并潰滅，至此完成球體入水。

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

2 有限元顯式動(dòng)力學(xué)方法在水動(dòng)力學(xué)中的應(yīng)用

隨著計(jì)算機(jī)技術(shù)的發(fā)展，如Ansys LS-DYNA、Abaqus/Explicit等，有限元顯式動(dòng)力學(xué)數(shù)值求解技術(shù)的仿真程序被開始應(yīng)用于出、入水等水動(dòng)力學(xué)問題的求解。

有限元的顯式動(dòng)態(tài)求解方法，作為隱式求解方法的補(bǔ)充，最初是為了模擬高速?zèng)_擊等高速動(dòng)力學(xué)問題而逐步發(fā)展而來的，特別適用于短暫、瞬時(shí)的，包含接觸、碰撞行為的非線性動(dòng)力學(xué)問題的模擬。通過引入用于描述水動(dòng)力學(xué)行為的狀態(tài)方程，使出、入水沖擊過程，也可通過顯式動(dòng)力學(xué)方法進(jìn)行求解。

2.1 AbaqusExplicit水動(dòng)力學(xué)求解方法

大型通用有限元軟件Abaqus中，AbaqusExplicit為其顯式動(dòng)態(tài)分析模塊，包含用于水動(dòng)力學(xué)模擬的狀態(tài)方程；應(yīng)用其中的耦合歐拉-拉格朗日分析（Coupled Eulerian-Lagrangian analysis，CEL）功能，即可實(shí)現(xiàn)相關(guān)的水動(dòng)力學(xué)問題求解。

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

但在如本文球體入水等流固耦合（fluid-structure interaction，F(xiàn)SI）分析中，歐拉分析雖然可有效處理流體流動(dòng)分析，但在捕捉結(jié)構(gòu)流固交界面上存在一定困難。此時(shí)，可應(yīng)用耦合歐拉-拉格朗日分析（CEL）功能進(jìn)行求解[2]。

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

2.2 AbaqusExplicit水動(dòng)力學(xué)模型

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

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

進(jìn)而得到：

常用的對(duì)材料雨貢紐曲線擬合關(guān)系為：

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

進(jìn)而得到線性Us-UpHugoniot形式的Mie-Grüneisen狀態(tài)方程，表述為：

通過以上狀態(tài)方程，即可定義水介質(zhì)的體積強(qiáng)度及靜水狀態(tài)體積響應(yīng)特性。同時(shí)，假設(shè)水介質(zhì)的剪切響應(yīng)與體積響應(yīng)是非耦合的，采用牛頓粘性剪切模型，定義其粘性剪切響應(yīng)。

上式中：S 為應(yīng)力偏量，e˙為應(yīng)變率偏量，μ 為粘性，同時(shí)，γ˙=2e˙為工程偏應(yīng)變率。

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

3 球體入水的數(shù)值模擬與實(shí)驗(yàn)對(duì)比

Abaqus顯式CEL方法，是基于有限元平臺(tái)的水動(dòng)力及流固耦合問題的數(shù)值求解方法。相對(duì)于CFD的流場(chǎng)求解功能，CEL方法由于對(duì)N-S方程進(jìn)行了簡化，不能考慮流動(dòng)邊界層效應(yīng)，因此，CEL方法側(cè)重于入水空泡的基本形態(tài)及演化，以及水彈道、水動(dòng)力載荷方面的應(yīng)用。

本文通過數(shù)值模擬與球體垂直入水實(shí)驗(yàn)[3]的比對(duì)，驗(yàn)證數(shù)值模型的有效性。球體垂直入水實(shí)驗(yàn)中，采用的標(biāo)準(zhǔn)臺(tái)球直徑為57.15 mm，質(zhì)量m=0.17 kg，球體以5.1 m/s的觸水速度垂直入水。

采用Abaqus動(dòng)態(tài)顯式CEL方法，進(jìn)行入水的流固耦合求解，球體采用Lagrange離散剛體（Discrete Rigid）單元，空氣域、水域采用Euler單元，采用歐拉-拉格朗日接觸算法直接耦合結(jié)構(gòu)網(wǎng)格和流體網(wǎng)格，進(jìn)行入水模擬。

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

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

基于以上原則進(jìn)行數(shù)值計(jì)算得到入水過程的入水空泡形態(tài)演化，并與實(shí)驗(yàn)比對(duì)如圖2、3所示。

圖2 入水實(shí)驗(yàn)入水空泡的演化Fig.2 Evolution of water-entry cavity in experiments

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

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

以自由液面為零點(diǎn)，豎直向下為位移的正方向，球體相對(duì)于自由液面入水的位移曲線，與實(shí)驗(yàn)數(shù)據(jù)比對(duì)，如圖5所示；球體入水的速度衰減曲線，與實(shí)驗(yàn)數(shù)據(jù)比對(duì)，如圖6所示。

通過實(shí)驗(yàn)數(shù)據(jù)與數(shù)值模擬入水彈道的一致性校準(zhǔn)結(jié)果，表明CEL方法對(duì)定量求解入水問題水彈道的有效性。

圖4 數(shù)值模擬入水空泡的深閉合及回射流形成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

值得注意的是，球體垂直入水速度衰減的實(shí)驗(yàn)數(shù)據(jù)，是由光學(xué)觀測(cè)的位移數(shù)據(jù)經(jīng)判讀得到的，因而實(shí)驗(yàn)的速度衰減數(shù)據(jù)更多地表現(xiàn)為位移數(shù)據(jù)點(diǎn)間的平均特性。

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

4 球體入水沖擊載荷

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

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

為精確地捕捉入水沖擊脈沖載荷峰值，以顯式動(dòng)力學(xué)方法時(shí)間增量的穩(wěn)定性極限Δtstable為參考基準(zhǔn)，作為輸出時(shí)間間隔，捕捉入水沖擊載荷峰值。

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

作為參考基準(zhǔn)的穩(wěn)定性極限Δtstable與網(wǎng)格尺寸直接相關(guān)，即入水沖擊脈沖載荷峰值的求解精度，直接與空間網(wǎng)格離散精度相關(guān)。

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

4.2 歐拉-拉格朗日接觸的網(wǎng)格離散

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

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

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

4.3 CEL方法入水沖擊載荷預(yù)報(bào)

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

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

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

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

如圖8所示，在不同的網(wǎng)格離散密度下，入水沖擊載荷曲線的峰值相對(duì)誤差控制在10%以內(nèi)，即可認(rèn)為達(dá)到入水沖擊載荷的求解精度。

5 結(jié) 語

本文針對(duì)球體入水過程，采用Abaqus顯式動(dòng)力學(xué)CEL方法，完成了球體入水空泡演化、發(fā)展過程的數(shù)值模擬。通過與實(shí)驗(yàn)入水彈道的對(duì)比，表明CEL方法在數(shù)值求解流固耦合問題中的有效性。在入水初期彈道校準(zhǔn)的基礎(chǔ)上，提出了在數(shù)值模擬入水沖擊載荷中，求解所需的采樣頻率、歐拉-拉格朗日接觸網(wǎng)格離散的要求，并完成球體入水沖擊載荷的預(yù)報(bào)。

參 考 文 獻(xiàn)：

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

[2]張偉偉,金先龍.球體撞擊自由液面相關(guān)效應(yīng)的數(shù)值模擬方法研究[J].船舶力學(xué),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]馬慶鵬,何春濤,等.球體垂直入水空泡實(shí)驗(yàn)研究[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]王永虎,石秀華.空投雷彈斜入水初始彈道數(shù)值分析[J].彈道學(xué)報(bào),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.