基于SPH方法的M型快艇入水載荷仿真分析

楊一凡，陳三桂，周維星，張濤*

1 華中科技大學 船舶與海洋工程學院，湖北 武漢 430074

2 中國艦船研究設計中心，湖北 武漢 430064

0 引 言

隨著船舶向高速、輕量化方向發展，具有吃水淺、速度快和靈活性等多項優點的M型快艇越來越受到國際造船界的重視。M型快艇由主船體、輔助船體和連接橋3部分組成。其中，M型快艇的中間區域為主船體，主要作用是通過排水為船體提供浮力；兩側的圍壁為輔助船體，主要起到密封作用；主船體與輔助船體之間則通過連接橋連接。M型快艇航行時，主船體、輔助船體、連接橋與水面圍成的區域中的空氣受到壓縮而形成空氣氣墊，其氣墊效應使船體抬升并快速且平穩地航行。同時，在M型快艇高速航行中，主船體與輔助船體之間的空間內會反復進出水，船底也會受到波浪反復的砰擊，嚴重影響M型快艇的安全性。因此，研究M型快艇入水過程中船體受到的砰擊載荷具有重要意義，其空氣墊效應也是不可忽視的。

Zhang等[1]通過拖船試驗研究了M型快艇在不同波浪作用下的砰擊載荷，但總體上國內外學者對M型快艇入水問題的研究仍較少。相反，對于與M型快艇十分相似的三體船結構入水問題的研究比較多。例如，Davis等[2]最早開始采用實驗和附加質量理論研究三體船入水問題，還開展了雙體船不同截面結構的入水實驗，運用數值方法分析入水過程，并針對結構入水時雙體船縫隙中殘留有空氣的問題，采用軟連接方式將水與結構物的附加質量聯系起來。此外，曹正林等[3]對高速三體船入水時影響連接橋砰擊壓力峰值的因素也進行了相關研究。

對于結構物入水問題的研究，最初基本上采用的都是試驗方法。例如，Lin等[4]通過平板入水試驗，發現平板在砰擊水面時存在空氣層，并比較了不同高度下平板自由落體入水產生的加速度和壓力。在數值計算方面，最具代表性的是源于馮·卡曼相關研究的楔形物體入水問題[5]。然而，除了研究著陸過程中的水對水上飛機的影響外，并未考慮空氣對其的影響。Zhao等[6]基于非線性邊界元理論，研究了具有任意形狀的二維結構的入水問題，在計算模型中檢查了質量、動量和能量的守恒性。之后，有學者在此基礎上也進行了廣泛和深入的研究[7-9]。不僅如此，針對結構物入水問題的研究，還發展出了不同的理論。然而，有學者發現無論是使用邊界元法（BEM）還是有限元法（FEM），在模擬結構物砰擊水面時，評估自由表面的非線性演化過程仍然具有挑戰性。

近年來，一種新的數值計算方法——光滑粒子流體動力學（smooth particle dynamics，SPH）算法逐漸受到學者們的關注。SPH算法是一種無網格方法，因沒有大變形引起的網格變形，便于跟蹤材料的軌跡，界面處理也靈活，所以SPH算法非常適合處理結構物入水過程中的液面大變形問題。近年來，該方法成為了研究結構物入水問題的重要方法。例如，閆蕊等[10]基于SPH算法研究了空氣對平板入水中的影響， Oger等[11]提出一種在SPH算法模擬過程中提取近固體邊界局部壓力的新方法，并模擬了楔形物體入水過程，所得結果與實驗數據吻合較好。此外，Yan等[12]使用SPH單相流模型模擬了三體船的部分入水情況。

基于現有研究，本文擬采用SPH液?氣兩相流算法，對M型快艇典型截面結構入水過程中的砰擊載荷進行研究。通過模擬平板和弓形結構的入水過程，將仿真結果與文獻試驗結果進行對比，以驗證基于SPH液?氣兩相流理論計算結構物入水結果的準確性。此外，通過SPH兩相流算法模擬不同斜升角和初始入水速度下M型快艇典型截面結構的入水過程，分析砰擊載荷的變化，討論斜升角對船舶入水過程的影響。

1 SPH液?氣兩相流理論

SPH方法是一種無網格的拉格朗日算法，其開源的DualSPHysics代碼是一種硬件加速的“光滑粒子流體動力學”代碼，已被證實可用于模擬當前結構，解決自由表面流動問題[13]。

1.1 光滑核函數

采用SPH算法建立的模型其性能很大程度上取決于光滑核函數的選擇，此函數在SPH算法計算中的作用至關重要。此外，偏微分方程所有解都需轉化為核函數，這在一定程度上也決定了計算的精確性。

目前，核函數種類繁多，為了提高計算效率，新的核函數不斷涌現。而常用的則是Wendland核函數[14]，具體形式如下：

式中：W(r,h）為Wendland核函數，其中r和h分別表示任意粒子a和b的距離以及光滑長度；在二維空間，αD= 7/(4πh2)，在三維空間，αD=21/(16πh3)；q為r與h的比值。

1.2 動量方程

SPH算法中的動量方程用于移動粒子，而人工黏度方案[15]則是SPH算法中進行流體仿真的常用方法，故可以寫為：

式中：Pa和ρa分別為粒子a的壓力和密度；va為粒子a的速度；Pb和ρb分別為粒子b的壓力和密度；mb為粒子b的質量；?a為粒子a的哈密爾頓算子；Wab為核函數；g為重力加速度；Πab為黏度項，并由式（3）給出。

式中：α為引入適當的消散而需要調整的系數；為 平均聲速，=0.5(ca+cb)， 其中ca，cb分別為粒子a和b的 聲 速； μab=vab·rab/(rab2+η2)， 其中vab表 示粒子a和b的速度差值，rab表示粒子a和b的距離， 且vab=va?vb，rab=ra?rb， η2=0.01h2。

1.3 狀態方程

SPH方法中將流體視為弱可壓縮，其壓力由基于粒子密度來確定，如下所示[16-17]：

式中：γ為常數項，對于水γ取 7，對于空氣γ取 1.4；ρ為密度；ρ0和c0為設置的密度和聲速。

1.4 連續性方程

在SPH模擬弱可壓縮流體的整個過程中，通過求解質量守恒方程或連續性方程來計算密度的變化，具體如下所示：

1.5 時間步

粒子a的速度、加速度和密度變化率分別表示為：

式中：ra為粒子a的距離；Fa為粒子a的加速度；Da是粒子a的密度變化率。基于Verlet方案[18]，對這些方程進行積分。

1.6 液?氣兩相流

在單相中使用一致形式的壓力梯度及速度散度來更新液相的動量方程和連續性方程。 對于氣相，因界面處的密度有較大的不連續性，故不適合采用相同的SPH公式，而使用如下形式[19]：

式中：ui，pi，ρi，?i分別為粒子i的速度、壓力、密度和哈密爾頓算子；mj，pj，ρj分別為粒子j的質量、壓力和密度；Πij和Wij分別為黏度項和光滑核函數；a為沖擊加速度。

1.7 邊界條件

邊界處理的核心思想是由被視為滿足不可滲透條件的單獨一組流體粒子來描述邊界。

2 算法驗證

2.1 平板模型

為驗證SPH兩相流算法的精確性，本文選擇文獻[4]試驗測試的平板入水結果作為依據。該試驗采用的平板參數如下：平板質量12.5 kg，尺寸0.45 m×0.2 m，平板從距自由液面0.1 m的高度自由下落，砰擊入水。

圖1所示為本文參照文獻[4]入水試驗的平板建立的SPH模型，采用此二維模型模擬平板的入水過程。設定平板寬為0.45 m，密度為13 888.89 kg/m3，厚度為0.01 m。計算域的水域和空氣域尺寸分別如圖所示。試驗時，平板初始位置距水面的高度為0.1 m，初始入水速度為0 m/s，重力加速度為9.81 m/s2。水及空氣的人工黏度分別為0.02和0.05，粒子的間距為2 mm，計算粒子的總數為1 028 905。圖1中，x軸對應水域的底部邊緣，y軸代表模型的三維方向，z軸則對應模型的中心軸線。

圖1 采用SPH兩相流算法的平板入水模型Fig.1 Simulation model of flat plate during water entry using twophase SPH method

經模擬，分別得到平板模型入水的加速度及受到的壓力分布結果。圖2為仿真結果與文獻[4]試驗結果的對比，圖3所示為模擬仿真的壓力云圖。結合兩圖可以發現，在130 ms時，盡管平板未接觸水面，但因平板與水面間存在空氣，空氣在平板和水面的壓縮作用下對平板產生壓力，導致平板的加速度開始緩慢上升。從圖3還可以明顯看出，當平板接近水面時，平板底部的空氣迅速逸出，但仍有部分空氣殘留在平板底部，形成了“空氣墊”，該空氣墊將在入水過程中起到緩沖作用。

圖2 平板入水加速度數值仿真與試驗結果對比Fig.2 Comparison of acceleration between simulations and experimental results of flat plate during water entry

圖3 平板入水過程中的壓力云圖Fig.3 Pressure contours of flat plate during water entry

文獻[4]試驗結果針對的是初始高度為0.1 m的平板入水過程，得到的最大壓力為19.796 kPa，即對應于文獻[4]給出的入水速度1.40 m/s的結果。根據文獻[4]中壓力測點的位置布置方式，本文提取了SPH仿真計算中平板中心點的壓力結果，如圖4所示。由圖可見，其最大壓力峰值為17.940 kPa。

圖4 平板試驗中心點的測試壓力[4]Fig.4 Pressures at centre point of flat plate in test[4]

從圖2和圖4可以看出，在兩相流下通過SPH模擬獲得的加速度和最大壓力值略小于文獻[4]得到的試驗值，且模擬的加速度脈沖寬度略大于試驗結果。其原因在于，在仿真中空氣無法從y軸方向溢出，導致仿真的緩沖效果比試驗的結果更明顯。總體上，兩相流下通過SPH模擬得到的結果與試驗結果擬合得較好。

2.2 弓形模型

為進一步驗證SPH兩相流算法的精確性，針對Aarsnes[20]弓形模型的入水進行了模擬。該文獻中所采用的測試設備的總重261 kg，弓形模型總長1 m，測試段長0.1 m。對于入水過程，采用壓力傳感器P1～P3測點（圖5（a）中所示位置）記錄壓力值。弓形模型截面的碰撞速度由速度傳感器測得，而相應的加速度則通過速度差獲取。對于給定的沖擊加速度，沖擊力由式（10）給出。

式中：Fv為沖擊力；m為截面質量。

圖5（b）所示為本文相應建立的SPH模型，模型質量為261 kg/m，水和空氣的人工黏度為0.02和0.05。為了保證計算的精度，SPH模型中的粒子間距為2 mm。在本節中，由于重點不是研究弓形模型的入水，而是對SPH算法的精確性進行驗證， 因此，僅模擬了斜升角為0°且沖擊速度為0.61 m/s時弓形模型的入水過程。模擬結果與文獻[20]的試驗結果進行了對比，如圖6和圖7所示。其中，圖6所示為弓形截面所受沖擊力，圖7所示為P1～P3測點的壓力。不僅如此，還將Sun[21]使用BEM方法獲得的計算結果，以及Wang[22]使用LS-DYNA有限元方法獲得的模擬結果一起進行了比較。

圖5 弓形結構截面及SPH模型Fig.5 Arched structure's cross-section and SPH model

從圖6可以看出，與BEM和LS-DYNA方法得到的結果相比，本文采用SPH兩相流算法計算的結果更接近于文獻[4]獲得的壓力峰值及其對應的時間。

圖6 弓形結構截面的砰擊載荷圖Fig.6 Slamming loads on the arched cross-section

從圖7可以看出，本文采用兩相流SPH算法所得結果與文獻[4]中3個測點的記錄結果相吻合，這證明了本文采用的兩相流SPH算法具有較高的準確度。總之，本文在兩相流下的SPH仿真結果與文獻[4]所給出的試驗結果擬合較好，驗證了兩相流下SPH模擬計算方法的可靠性。

圖7 P1～P3測點記錄的壓力Fig.7 Pressures recorded at three measuring points

3 M型快艇結構仿真模型

圖8所示為本文采用的M型快艇截面幾何形狀（單位：mm）。該截面長度為846.8 mm，高度為305.6 mm，重量為80 kg/s。為了分析典型位置的壓力，選擇布置了4個壓力測點。P1～P4測點分別位于點A、AB和CD、DE段中心處。仿真時，本文選擇了如圖9所示4種主船體斜升角α（13°，18°，23°，28°），以研究該斜升角對M型快艇船體結構入水過程的影響。研究中，不同斜升角時的模型下落高度及質量與文獻[4]相同。

圖8 M型快艇截面示意圖Fig.8 Geometry of high-speed M-boat cross-section

圖9 不同斜升角下M型快艇截面示意圖Fig.9 Diagram of high-speed M-boat cross-section with various dead-rise angles

圖10所示為相對應的二維M型快艇截面的SPH模型。空氣域和水域尺寸分別為3 000 mm×1 000 mm和3 000 mm×1 500 mm。水及空氣的人工黏度分別為0.02和0.05，快艇截面模型的質量為80 kg/s。同時，考慮到計算時間的影響，SPH數值仿真中粒子間距設為2 mm，后文將詳細討論粒子間距的影響。圖10中，x軸對應于水域底部邊緣，y軸代表模型三維方向，z軸則對應于模型中心軸線。除了考慮主船體斜升角的影響外，本文還將討論入結構入水速度的影響，具體工況如表1所示。

表1 模擬工況Table 1 The modelling conditions

圖10 M型快艇截面的SPH模型Fig.10 The SPH model of igh-speed M-boat cross-section

4 結果分析

4.1 粒子間距的影響

本文選擇粒子間距為5，3.3，2.5，2和1.5 mm。圖11給出了斜升角為23°時計算的M型快艇截面加速度， 模擬中所有初始入水速度均取為5 m/s。從圖11可以發現，隨著粒子間距的減小，加速度曲線將逐漸穩定, 且粒子間距分別為2和1.5 mm的結果吻合較好，兩者的差別較小，說明繼續縮小粒子間距對結果的影響不大。因此，考慮到計算時間的合理性，本文選擇了2 mm作為初始粒子的間距。

圖11 不同粒子間距下M型快艇截面加速度（α=23°）Fig.11 Calculated acceleration of high-speed M-boat crosssection with various particle spacings (α =23°)

4.2 加速度結果

圖12所示為在SPH兩相流模型計算中，不同斜升角和入水速度下得到的M型快艇加速度的計算結果。模擬中，船舶截面斜升角分別為13°，18°，23°和28°，對于每個斜升角，分別選擇了5，7.5和10 m/s這3種不同初始入水速度。

從圖12中可以很好地觀察到，船舶加速度曲線有2個峰值。在船舶開始運動后，加速度將迅速增加，達到第1個峰值后快速減小，然后再較緩慢地增加并達到第2個峰值，最后又慢慢降低。根據加速度結果，可以推測出船舶截面在入水過程中遭受了2次砰擊。關于這2次砰擊的具體位置及其影響在下文詳細說明。

圖12 不同斜升角下的加速度計算結果Fig.12 Calculation results of acceleration for various dead-rise angles

保持船舶截面斜升角不變，在初始入水速度快的工況下，加速度峰值都更大，這表明受到的砰擊載荷更大，且峰值出現的時間都更早。為進一步分析斜升角對船舶入水加速度的影響，提取出不同工況下加速度的2處峰值并予以對比，結果如圖13所示。

圖13 不同初始入水速度下加速度計算結果與斜升角間的關系Fig.13 Relationship between calculated acceleration and dead-rise angle at various initial water entry velocities

從圖13可以看出，在相同初始入水速度下，隨著斜升角的增大，加速度第1個峰值減小，也即第1次砰擊載荷降低，而第2個峰值基本沒有太大變化。但是，若斜升角增大到一定程度時，可見第2個峰值顯著增加，也即第2次砰擊載荷突然增大，具體原因后文將進行分析。

4.3 測點壓力結果

圖14～圖17顯示的是圖8中所描述的壓力測點P1～P4的壓力監測結果。

從圖14～圖17中可以看出，在斜升角確定的情況下，隨著模型初始入水速度的增加，各測點的壓力峰值增大，且除測點P1外，其余測點壓力到達峰值的時間明顯提前。同時，還可以明顯看出，隨著模型初始入水速度的增加，測點P1壓力的脈寬減小。這說明船舶初始入水速度越快，受到的壓力載荷會越大，入水所需時間越短。

圖14 斜升角為13°時測點記錄的壓力變化Fig.14 Pressure histories recorded at measuring points when deadrise angle is 13°

圖15 斜升角為18°時測點記錄的壓力變化Fig.15 Pressure histories recorded at measuring points when deadrise angle is 18°

圖16 斜升角為23°時測點記錄的壓力變化Fig.16 Pressure histories recorded at measuring points when deadrise angle is 23°

圖17 斜升角為28°時測點記錄的壓力變化Fig.17 Pressure histories recorded at measuring points when deadrise angle is 28°

根據上節中加速度曲線分析中2次砰擊及對應的時間可以推斷，P1～P2測點的壓力在第1次砰擊過程中達到峰值，在第2次砰擊過程中卻波動不大。而P3～P4測點的壓力則相反，在第2次砰擊過程中達到峰值。同時，結合測點位置可以得出結論，即第1次砰擊發生在船舶主船體的龍骨處，而第2次砰擊則發生在連接橋和輔助船體之間。

4.4 M型快艇入水過程

為了更詳細了解M型快艇入水過程及過程中加速度和壓力變化，本文選取初始入水速度為5 m/s、斜升角分別為13°和28°的算例，提取典型時刻的壓力云圖，如圖18所示。

從圖18（a）并結合圖19加速度曲線可以看出：

圖18 不同斜升角下自由液面壓力云圖及高度Fig.18 Pressure contours and free surface elevation at various dead-rise angle

在t=0 ms時，船舶開始向下運動, 約7 ms時，對應于圖19中的加速度曲線第1個峰值達到最大值，此時，船舶結構上的最大壓力位于截面的龍骨處。這進一步表明加速度的第1個峰值是由水面與龍骨之間的相互作用引起的，即第1次砰擊發生在船舶的龍骨位置。

在t=21.4 ms的瞬間，對應于圖19中A處，水射流砰擊到連接橋，使得加速度曲線出現了一個小尖峰。

圖19 v=5 m/s時不同斜升角的加速度結果Fig.19 Acceleration results for various dead-rise angle at v=5 m/s

在t=25.2 ms時，船的輔助船體接觸水面，之后殘留于主船體與輔助船體間的空氣因無法溢出而受到壓縮，進而對船舶施加了一個作用力，使得加速度曲線開始緩慢上升，并在31 ms時刻瞬間達到峰值，形成加速度的第2個峰值。這表明加速度的第2個峰值發生在船的輔助船體接觸水面后，即第2次砰擊發生在連接橋和輔助船體之間，也印證了壓力分析的結果。

同時，從圖18（a）中還可以看出，空氣對于船舶入水有著非常重要的影響，特別是在M型快艇的側板接觸水面之后的入水過程。

對于圖18（b），即主船體體斜升角為28°時的M型快艇截面壓力云圖，與圖18（a）所示斜升角為13°時的壓力結果相似。但是，由于斜升角的增大，主船體底部側面的長度增加，水面淹沒主船體底部的時間變長，導致主船體船底射流砰擊連接橋的過程發生在輔助船體接觸水面之后，時間對應于圖19中的A處和圖18(b)中的t=24 ms時刻。這也是主船體斜升角為28°時加速度曲線的第2個峰值突然增加的原因所在。

圖20所示為斜升角為28°、初始入水速度v=5 m/s、有/無空氣時M型快艇入水計算得到的加速度結果。從圖20可以發現，加速度曲線的第1個峰吻合較好，而對于第2個峰值有空氣作用的結果明顯小于沒有空氣作用的結果，同時有空氣作用的加速度曲線脈寬更寬。這說明空氣在第2個峰值形成過程中起到緩沖作用，在輔助船體接觸水面后，受到壓縮空氣的作用，減緩了模型的下落，使得船體受到的砰擊載荷更均勻，起到類似于“空氣墊”的作用。

圖20 斜升角為28°時有/無空氣的加速度結果Fig.20 Acceleration results with or without air influence when dead-rise angle is 28°

在第2次砰擊發生時，因部分空氣無法及時溢出，在主船體周圍形成了氣穴，如圖21所示。這類似于文獻 [2]在三體船的入水測試中觀察到的氣穴結果，如圖22所示。根據文獻[2]的測試結果，隨著開始入水，空氣流向形成氣墊的龍骨位置。該結果表明，SPH兩相流算法可以輕松模擬出結構入水過程空氣中形成的渦旋，而有限元方法則很模擬得到。

圖21 氣穴Fig.21 Air pocket

圖22 三體船的入水測試圖[2]Fig.22 View of water entry test of trimaran [2]

5 結 論

本文采用SPH兩相流算法研究了M型快艇入水過程，利用該算法模擬了平板和弓形模型的入水問題，并將仿真結果與文獻中已有試驗結果進行了對比，從而驗證了所用算法的精確性，通過仿真模擬還觀察到了“空氣墊”現象。

通過對M型快艇截面SPH模型的仿真模擬，發現船舶在入水過程中受到了二次砰擊作用，分別發生于主船體與水面的接觸過程和側體接觸水面之后。其中，第2次砰擊發生的原因主要是，船舶與水面對空氣的壓縮作用形成壓力，并與主船體底部射流砰擊連接橋共同作用。受到壓縮的空氣在第2個峰值形成過程中起到了緩沖作用，使得側體接觸水面后，主船體受到的砰擊載荷更均勻，其作用類似于一個“空氣墊”。

斜升角對M型快艇的入水過程有很大的影響。對于第1次砰擊，因龍骨幾何形狀是楔形的，故對于較大的斜升角，龍骨處受到的砰擊載荷相對較小。對于第2次砰擊，在主船體斜升角大到一定程度時，結構受到的砰擊載荷會顯著增加，這是因為M型快艇和水面對空氣的壓縮以及主船體底部射流對連接橋砰擊的共同作用所致。因此，建議設計M型快艇時選擇合適的斜升角。

此外，第2次砰擊發生時，因部分空氣無法及時溢出會在主船體附近形成氣穴，而SPH兩相流算法能夠非常好地對此現象進行模擬和分析。