基于SPH-FEM 耦合方法的回轉體高速入水數值研究

高英杰，剛旭皓

(1.中國船級社青島分社，山東 青島 266071；2.中國船舶科學研究中心，江蘇 無錫 214082)

0 引 言

入水沖擊問題具有廣泛的工程應用背景，如船舶航行受浪砰擊，海洋工程結構物下水安裝，救生艇自由降落釋放，導彈、魚雷入水等。隨著海洋工程計算和武器裝備的發展，入水沖擊問題的研究受到持續關注。入水沖擊作為一個復雜的物理現象，伴隨著流固耦合作用，而高速入水大大加劇了流固耦合的復雜程度，同時結構物撞水瞬間受到的砰擊載荷可能導致結構運動變化和結構變形失效。因此有必要考慮流固耦合作用，分析入水結構的動態響應，為結構物高速入水設計提供參考。

對于入水沖擊問題，國內外學者做了大量研究。1929 年，Von Karman首先引入附加質量的概念，推導計算出二維楔形體入水沖擊載荷，奠定了入水沖擊理論計算的基礎。Wagner在1932 年，針對楔形體入水沖擊問題提出了小傾角平板近似理論，通過求解伯努利方程計算出二楔形體入水過程中的壓力分布。后續學者在二者的研究基礎上，通過試驗、數值計算等方式對入水沖擊問題開展了廣泛研究。Worthington采用高速攝影技術，最早記錄了小球入水的流場演變情況，分析了空泡的變化規律。MAY開展了不同形狀參數和入水速度的試驗，研究了結構物入水過程中的載荷特征。顧建農開展了彈體高速入水試驗，分析了形狀參數和入水初速度對彈道特征和運動特性的影響。張偉通過試驗分析了不同頭型彈體入水的彈道穩定性，得到了彈體入水速度衰減公式。

數值計算成本較低且不受場地、設備等現實條件限制，可較為準確地預報入水過程中結構和流場的變化情況，逐漸成為研究入水沖擊問題的關鍵方法之一。陳震開展了船舶在波浪中出入水的數值模擬，分析了入水運動參數對抨擊壓力的影響。何春濤采用有限體積法，建立了結構物垂直入水的計算模型，對比試驗結果分析了入水的運動軌跡和空泡閉合規律。馬慶鵬采用VOF 方法開展了不同入水速度下回轉體的數值計算，研究了入水沖擊載荷變化規律以及侵水深度和空泡發展的關系。李強基于Ls-dyna 軟件計算了魚雷入水過程，分析了不同頭型參數對魚雷水下運動的影響。楊力采用ALE 方法，分析了平底結構入水過程中質量、剛度對砰擊壓力的影響。路麗?；赗ANS 方程采用動網格技術，模擬了凹形體入水過程，分析了液體噴濺對角運動的影響。

目前針對結構入水的研究多集中于剛體低速入水，主要關注入水運動軌跡和空泡演變規律，考慮流固耦合的高速入水研究較少。本文進行回轉體高速入水流固耦合數值計算，分析回轉體高速入水的運動特性和結構響應，其結果可以為結構物高速入水的預報提供參考。

1 流固耦合計算模型

采用光滑粒子流體動力學（SPH）和有限元(FEM)耦合方法對高速入水過程進行模擬，其中?；剞D體模型采用有限元Lagrange 網格，水域采用SPH 粒子。FEM 方法在計算結構運動和動力響應方面,具有更高的精度和效率；SPH 方法中計算域使用粒子形式表示，無需網格劃分，在此計算介質大變形，可以有效避免有限元方法中網格畸變等情況。FEM/SPH 耦合可發揮兩者優勢，在計算入水沖擊等介質大變形的動力問題上，具有較好的準確性和較高計算效率。

1.1 SPH 方法基本理論

SPH 方法使用粒子分布表達計算域，采用核函數近似的方式將求解的偏微分方程轉換為積分形式，使用粒子近似的方式將核函數方程離散化，通過求解離散方程得到粒子場內變量。

在SPH 方法中，給定任意一個函數，其積分表達形式如下：

式中：Ω為積分域；δ (-)為狄拉克函數。

由上式可知，任意連續函數()采用狄拉克函數可以表示為以下積分形式：

式中：(-,)為核函數或稱光滑函數，核函數的選取往往會影響計算精度和效率。

根據分布積分和散度定理，可以將上式進行變換：

最后，將變換后的函數進行離散化處理，得到域內所有粒子疊加的離散形式：

式中：和為粒子編號；和為粒子的質量和密度；為粒子總數 。

所以，在SPH 方法中對于基本控制方程可以表示為如下離散形式：

1.2 界面耦合模型

在入水過程中,實現結構體和水中介質耦合作用，需要在SPH-FEM 耦合方法中，將網格和粒子之間的界面進行處理。結構體網格和水粒子之間需要采用接觸算法進行相應的計算處理。

在接觸算法中，每一時間步內計算有限元網格節點穿過SPH 粒子的距離，罰函數根據有限元網格和SPH 粒子間的穿透距離計算出兩者之間的接觸力，通過接觸力實現兩者之間的相互作用，接觸力將SPH 粒子拉出交界面，避免SPH 粒子穿越有限元網格，導致計算失效。

1.3 計算模型設置

在計算模型中，回轉體最大直徑 120 mm，頭部直徑 100 mm，回轉體長度為0.5 m，結構厚度為5 mm，其首端傾斜斜面與回轉體頭部平面所呈角度為 93.18°，入水速度為100 m/s，采用有限元網格進行劃分。水域尺寸為 1.2 m×1.2 m×3.6 m，采用SPH 粒子表示，回轉體外形及流域如圖2 所示。

圖2 計算域設置Fig.2 Computational domain settings

SPH 粒子數為1 536 000，使用 Gruneisen 狀態方程控制，回轉體采用4 節點殼單元，網格數量為520 000，材料模型考慮彈性和塑形，材質為45 號鋼材，材料性能參數見表1。

表1 材料模型參數Tab.1 Material parameters

2 計算結果分析

2.1 運動特性分析

首先進行回轉體入水速度變化分析并與理論公式對比，驗證計算方法的準確性。假設物體在入水過程中空泡內部和外部壓力差保持一致，在空化數不為定值的情況下，存在如下規律：

式中：Δ為空泡內外壓力差；C為氣流壓力降系數；σ為初始空化數；σ=0.006-0.018 ；ρ為空氣密度；ρ為水密度；v為入水初速度；v為入水后速度。

阻力系數與空化數關系：

忽略重力的情況下，入水運動方程為：

式中：m為物體的質量；A為觸水面積。

通過式(9)和式(10)，可得到速度衰減公式：

式中：衰減系數k=ρ A C/ 2m。

圖3 為衰減公式與數值計算得到的速度曲線?？梢钥闯?，觸水瞬間回轉體高速撞擊水面受到極大的瞬間沖擊載荷，速度衰減劇烈，隨著入水進程加深，速度衰減減弱，最后近似線性衰減。計算得到的速度曲線與公式曲線具有好的一致性，驗證了SPH-FEM 耦合計算模型在計算高速入水方面的有效性。

圖3 回轉體速度衰減曲線對比圖Fig.3 Curve of velocity attenuation comparison

2.2 動力特性及結構動力響應分析

王珂擬合了回轉體入水抨擊壓強峰值，入水過程中回轉體頭部中心處峰值壓強可以表示為：

式中：k=0.004 23×d+0.003 758×d+0.42，ln+15，k=+5；ρ為水的密度；為厚度；為材料彈性模量；為入水的初速度。

由式(11)可以計算出回轉體入水頭部中心錯誤壓力峰值為1.73×10Pa，采用SPH-FEM 耦合計算得到的回轉體頭部中心處壓力情況如圖4 所示，其峰值為1.92×10Pa?？梢钥闯?，數值計算得到的回轉體頭部中心處壓力峰值與文獻公式計算結果接近，數值計算結果稍大。在數值計算過程中未充分考慮高速入水過程中壓縮空氣的影響，導致結果略微偏大。高速入水是一個涉及氣-液-固三相耦合的瞬時沖擊問題，空氣高速壓縮帶來的空氣墊效應一直作為數值計算中的難點，在未來具有廣闊的研究空間。

圖4 回轉體中心單元壓力圖Fig.4 Pressure intensity curve in the central head region of the revolution body

圖5 給出了入水過程中回轉體垂直方向受力情況。從圖中可以看出，撞水瞬間回轉體頭部受到的沖擊力急劇增加，在毫秒量級時間達到峰值，峰值為1.37×10N，頭部觸水后，峰值迅速下降，隨著回轉體完全進入水面，沖擊力呈小幅度震蕩衰減，逐漸趨于穩定。

圖5 回轉體沿軸向受力圖Fig.5 Curve of the axial force of the revolution body

高速入水過程中，撞水初期會產生瞬時巨大沖擊力，可能導致結構變形甚至失效，因此有必要研究入水沖擊的結構強度問題。圖6 為回轉體頭部中心單元的等效應力曲線，圖7 給出了回轉體頭部中心單元的塑形應變曲線。本次計算考慮流固耦合效應，物體撞水后會產生變形緩沖抨擊，變形呈反復狀態。在應力曲線上表現為：出現瞬時峰值后，應力衰減呈現波動震蕩形式，隨著入水進程加深，震蕩的頻率和幅度減小，伴隨著整個入水進程。在100 m 垂直入水情況下，回轉體頭部中心處應力峰值迅速達到材料的屈服應力 355 MPa，材料發生塑形應變，塑性應變值0.146，雖然未出現破壞現象，但是在設計過程中有必要加強結構物觸水位置處的結構強度，避免結構失效。

圖6 回轉體頭部單元等效應力曲線Fig.6 Equivalent stress curve of the central element of the revolution body head

圖7 回轉體頭部單元塑形應變曲線Fig.7 Plastic strain curve of the central element of the revolution body head

2.3 入水空泡特性分析

入水空泡形態是研究空泡問題的關鍵，研究入水空泡動力學問題在基于一定的假設條件下可以推導出在入水過程中的侵徹距離與空泡半徑的函數關系。此時假定入水過程中的能量流失轉化為回轉體周圍空泡區域的流體動能和空泡介質勢能。利用分布點源理論推導出回轉體周圍流體動能為：

式中：為空泡半徑；為無量綱系數，在忽略重力作用下，對有限長度的平頭回轉體，其空泡內勢能E可以表示為：

在忽略入水初始階段入水速度的影響下，同時根據能量守恒原理有：

通過式(12)～式(16)，同時應用入水速度關系式v=(-)/(-)，最終可以求得在入水過程中的空泡半徑隨入水侵徹距離的函數關系：

式中：為時刻的位移值，為回轉體頭部的空泡半徑長度?？梢钥闯?，回轉體入水過程中的空泡尺寸與回轉體的形狀、經驗系數、阻力系數以及空化數等有關，對高速入水問題學者提出經驗系數取值一般小于3。

在得到圖8 各個時刻回轉體入水空泡體積分數的前提下，提取回轉體入水的最終狀態即0.005 s 時刻的空泡氣-液交界面輪廓，同時利用式(17)的空泡半徑與侵徹距離公式關系，取變量數值=2。對比數值模擬和理論計算得到的空泡尺寸輪廓如圖9 所示。從圖中可以看出2 條曲線表征的空泡尺寸吻合較好，證明了數值模擬中計算空泡的準確性。

圖8 回轉體入水開空泡過程Fig.8 Cavity feature during the water-entry process

圖9 回轉體泡輪廓與理論模型對比Fig.9 Comparison of the cavity model with the simulation

另一方面，理論公式為剛體入水計算結果，而本次數值計算考慮到材料的彈塑性，更加接近于實際情況。由于回轉體在入水過程中存在變形，可以看出2 條空泡輪廓曲線存在一定的差別。彈性體入水的首端空泡尺寸相比剛體結果較小，回轉體頭部為入水過程中變形最明顯的區域，在入水過程中，頭部的變形使流體進入內凹區域，導致彈性體的入水空泡形狀首部區域不再是平面，使其首部的直徑變小，隨著入水進程加深，后體區域的尺寸差異隨空泡直徑的增大逐漸減小。

2.4 入水速度對入水進程的影響

為進一步探究入水速度的影響，開展回轉體以60 m/s，80 m/s，100 m/s 速度條件下的計算。表2 給出了不同入水速度下回轉體頭部中心處壓力峰值，圖10為不同入水初速度下的歸一化速度曲線?？梢钥闯?，回轉體入水初期受到極大的沖擊力，速度衰減劇烈，之后沖擊力迅速降低，回轉體速度衰減程度放緩明顯，最后速度衰減表現為類似線性變化?；剞D體速度衰減情況受初速度影響，初速度越大速度衰減程度越大，隨著入水進程加深，不同初速度下的回轉體速度衰減程度接近，趨于穩定。

圖1 有限元網格和SPH 粒子耦合界面Fig.1 Coupling interface of FEM and SPH

表2 回轉體頭部中心壓力峰值Tab.2 Results of pressure peak of cylinders head central point

圖10 歸一化速度衰減曲線Fig.10 Curve of normalizes velocity attenuation

回轉體入水阻力系數計算公式為：

其中：為阻力值；ρ為水密度；為入水速度；為特征面積。

圖11 為初速下回轉體入水的阻力系數曲線?？梢钥闯?，在不同入水速度下回轉體阻力系數變化情況趨于一致，阻力系數峰值接近，入水速度越大，峰值持續時間越短，隨著入水進程發展，流動穩定，3 種情況的阻力系數值保持接近。由此可見，入水初速度對回轉體阻力系數的影響不明顯。

圖11 阻力系數曲線Fig.11 Curve of Cd with different velocity

3 結 語

本文采用SPH-FEM 耦合方法對回轉體高速入水過程進行數值分析，得到以下結論：1）通過與速度衰減和回轉頭部壓力峰值的理論公式對比，表明SPHFEM 耦合方法可以較為準確地預報回轉體高速入水過程中的運動及動力特征，可以較好地解決數值計算過程中由于流體介質不連續變形帶來的計算困難等問題。

2）回轉體觸水瞬間產生應力峰值，出現時間極短，受變形影響，隨后應力衰減呈現波動震蕩形式，其頭部中心處應力峰值超過了材料的屈服應力，材料發生塑形應變，在設計上可考慮對頭部結構進行加強。

3）數值計算中得到的入水空泡輪廓尺寸數據與剛性回轉體的理論空泡尺寸結果總體吻合良好，由于彈性變形的影響，在空泡輪廓頭部存在一定的差異，后體區域的尺寸差異隨空泡直徑的增大逐漸減小。

4）回轉體速度衰減受入水初速度影響，初速度越大速度衰減程度越大，隨著入水進程加深，不同初速度下的回轉體速度衰減程度接近，入水速度對阻力系數變化影響不大。