基于SPH／FEM方法的淺水沖擊噴濺及其抑制結構研究

任選其， 徐 緋， 張顯鵬， 楊 揚

（西北工業大學 航空學院，西安710072）

1 引 言

淺水沖擊問題是一類特殊的結構入水沖擊問題，在二者相互作用過程中，結構受到淺水層的沖擊載荷，同時對含有底部邊界的淺水層造成擾動和擠壓，引起自由液面的劇烈噴濺。以往對于淺水沖擊問題的研究很少［1］，而淺水沖擊及其噴濺會引起如齒輪與液膜的沖擊、輪胎滑水和濺水進入發動機等很多工程問題，因此受到學者們的日益關注。此外，如何設計合理的抑制結構以減小噴濺，也是研究的熱點問題。

淺水沖擊噴濺由于其涉及到結構、地面和水的相互作用，具有強非線性，傳統的Von Karman方法和 Wagner理論方法不再適用，Korobkin［2］最早提出基于匹配漸進分析來研究兩個剛性結構對之間液膜的沖擊。文獻［3，4］先后采用Korobkin方法，研究了剛性平底和弧形底結構在不同水深條件下的入水沖擊力。但是理論方法往往忽略或簡化了噴濺的影響，難以對噴濺細節進行研究。Bukreev等［5］開展了剛性長方體準二維淺水沖擊試驗研究，發現淺水沖擊底部有射流形成，但沒分析液面噴濺機理。

噴濺抑制結構的研究主要針對水上飛機、高速艦船和輪胎等的噴濺抑制問題，常見的抑制結構包括舭彎、擋水板和翻邊等。對于深水噴濺，如水上飛機和高速艦船，常采用拱形舭彎設計，但其在高海況時的效果不佳，日本PS-1水上飛機首次采用了抑波槽設計，對深水噴濺具有較好的抑制效果，而其對淺水噴濺的抑制效果和機理仍有待研究。對于淺水噴濺，翻邊是應用最廣泛的抑制結構之一，文獻［6，7］先后開展試驗，比較了斜角、直角和拱形翻邊構型的抑制效果，其后福克公司和固特異輪胎公司通過大量全機滑跑試驗，將拱形翻邊應用于Fokker100等飛機輪胎上，并申請了相關專利［8］，但適航試驗表明，拱形翻邊對于水層較厚時的噴濺抑制效果有限，翻邊構型有待進一步優化。

噴濺及其抑制結構的研究通常需要通過大型試驗進行，耗費巨大，且難以對噴濺機理進行研究，而數值方法作為研究手段可以方便地對各種剖面構型抑制結構的效果和機理展開分析。對于自由液面的計算，一般采用ALE方法、Level Set方法和VOF方法，這些方法基本都是基于歐拉觀點進行計算，均需通過特殊算法進行液面捕捉和重構，易產生網格畸變，導致計算收斂性變差和精度下降。SPH基于拉格朗日觀點，并且具有無網格粒子特性，可以直接獲得自由液面形態，便于處理大變形問題［9］，近年來廣泛應用于自由表面流［10，11］等相關研究領域。結構入水沖擊是一種具有自由液面噴濺的流固耦合問題，采用SPH／FEM耦合方法進行建模，通過交界面處的接觸算法，將SPH粒子法模擬流體部分和FEM算法模擬固體結構部分結合起來，可以發揮方法各自的優勢。

本文基于LS-DYNA中的SPH／FEM耦合方法，首先建立長方體結構的淺水噴濺模型，并與文獻結果進行比較驗證，對淺水噴濺機理做了研究，將淺水噴濺分為三部分。其次，針對淺水噴濺抑制結構，分析了Fokker100機型采用的拱形翻邊結構對各部分噴濺的抑制機理。最后，提出一種帶有凹槽和阻水邊條的濺水抑制結構形狀，通過數值方法說明其效果。

2 SPH／FEM耦合方法

2.1 SPH方法基本方程

SPH方法是一種基于粒子離散的無網格拉格朗日方法，該方法的本質是通過光滑核函數對守恒方程進行離散，再通過粒子間的相互作用進行插值［12，13］。對于任意未知函數，采用式（1）近似x 處的函數值，

質量守恒

式中 W 為光滑核函數，具有與狄克拉δ函數相同的峰值性和正規化條件約束，h為光滑長度。用SPH方法對Navier-Stocks方程進行離散后得到，

動量守恒

能量守恒

式中 Pi為i點的各向同性壓力，mj為第j個粒子的質量，為第i個粒子在β 方向的坐標，和分別為第i個粒子的應力張量和應變張量，＝（－）為兩個粒子之間的相對速度在β方向上的分量，N為光滑長度范圍內的粒子數，ηi為粘性系數

2.2 接觸算法

實際計算中SPH和FEM相對獨立計算，每個計算歩通過罰函數接觸算法傳遞接觸狀態和作用力，保證SPH水粒子在計算過程中不會發生穿透，同時實現兩者力的相互作用［14，15］。

其中，接觸面處的法向接觸力表達式為

式中ki為接觸剛度，fs為ki的比例因子，li為SPH水粒子相對于FEM結構網格單元的穿透深度，n為結構單元的法向單位矢量。

接觸面的切向接觸力表達式為

式中μk為水粒子在網格單元表面的摩擦系數。

3 淺水沖擊噴濺機理研究

3.1 長方體結構淺水沖擊模型

建立結構淺水沖擊模型，如圖1所示。通過FEM方法建立剛性長方體結構模型，參考文獻［5］的準二維長方體入水模型尺寸，采用LS-DYNA提供的SOLID164實體單元建立結構模型，長方體結構寬度為92mm，結構總質量為6.2kg，通過附加質量加在剛體重心，材料采用RIGID材料。采用SPH建立水池模型，初始密度為998kg／m3，水池深度為80mm，寬度為800mm，粒子間距為1 mm，總粒子數為64000。厚度方向為2mm，建立準二維模型，右側設為對稱邊界條件，地面設為剛性邊界，長方體結構和水粒子的耦合通過罰函數設置點面接觸實現，接觸剛度系數ki＝0.1，比例因子fs＝1，摩擦系數μ＝0.1。

流體部分采用NULL材料模型，NULL材料模型不考慮偏應力的計算，通常用于定義流體。由于液體受到結構底部和地面之間的擠壓作用，因此考慮了液體的可壓縮性，液體的壓力和體積的關系通過Gruneisen狀態方程進行描述：

式中 C為vs－vp曲線的截距（聲速），S1，S2和S3為vs－vp曲線的斜率系數，γ0為Gruneisen系數，α為γ0的一階體積修正系數，E為材料內能，ν＝（ρ／ρ0）－1為體積變化率。根據參考文獻［16］選取狀態方程參數，列入表1。

圖1 長方體結構淺水沖擊模型Fig.1 Schematic diagram of shallow water impact model of rectangular structure

3.2 計算結果與試驗比較

入水初速為2m／s，結構重量為6.2kg時，濺水形態計算結果與Bukree試驗結果［5］的對比，如圖2所示。可以看出，兩者吻合良好，說明本文構建的SPH／FEM數值模型可以較好地模擬淺水沖擊噴濺現象。

本文通過數值模擬發現，淺水沖擊噴濺可分為初始噴濺、二次噴濺和底部射流三部分。其中，初始噴濺為入水初期產生位于噴濺端部的小股噴濺；二次噴濺為噴濺中后部的大量水流，是淺水噴濺的主體部分；底部射流為由結構底部擠出，沿地面的高速流動部分。Bukree試驗也指出，淺水沖擊有明顯的底部射流現象。

表1 流體材料狀態方程參數Tab.1 State equation parameters of fluid material

3.3 淺水噴濺來源區域及其速度分布

根據對淺水噴濺的區域劃分，分別研究初始噴濺、二次噴濺和底部射流的噴濺來源區域、速度分布和噴濺形成機理。

首先，如圖3（a）所示，對二次噴濺區的水粒子進行標記，分析各部分濺水的初始來源區域，如圖3（b）所示，黑色部分為標記后的二次噴濺區域。可以看出，初始噴濺來源于靠近結構邊緣的表層水，二次噴濺主要來源于結構外側和底部水層，底部射流來源于結構中部區域。

其次，分析各個區域特征水粒子速度隨時間的變化，如圖4所示，在底部射流區域的上表面由內到外間隔10mm標記水粒子A1～A4，在初始噴濺區表面由內至外間隔20mm標記水粒子B1～B3，在二次噴濺區水深40mm位置，間隔40mm由內至外標記水粒子C1～C3。

圖2 濺水形態計算結果與試驗比較Fig.2 Comparison of calculation results and exprimental results of spray

圖3 噴濺來源區域示意圖Fig.3 Schematic diagram of spray source area

圖4 底部射流區域的特征點選取Fig.4 Feature points of the bottom jet area

圖5 噴濺速度隨時間的變化Fig.5 Spraying speed varies with time

圖5 為初始噴濺區和二次噴濺區標記水粒子豎直方向速度隨時間的變化。可以看出，在結構與水面接觸后，初始噴濺在短時間內速度迅速增加，約5ms時達到速度峰值8m／s，之后基本保持穩定。二次噴濺的濺水速度在結構入水初期2ms內迅速增大到約1m／s，其后略有減小，隨后再次增加，到達峰值2.6m／s后，在重力作用下濺水速度逐漸減小。

圖6（a）為底部射流區豎直方向速度的變化趨勢，可以看出，表面水粒子速度初始為負值，大小等于結構的下落速度，隨后速度在短時間內明顯增大，到達速度峰值后緩慢降低。圖6（b）為水平方向速度的變化趨勢，可以發現，射流區水粒子的水平方向速度在入水的初始階段較為穩定，到達臨界時刻后迅速增至一個較高的速度峰值，然后迅速減小，最后穩定在較低速度。

比較圖6中粒子A1～A4的速度變化趨勢，可以看出A1～A4的速度變化趨勢一致，速度峰值從A1～A4依次減小，可見越靠近結構中部，底部射流速度越高，水體受到的擠壓作用越明顯。

3.4 底部射流形成過程和能量傳遞機理

從圖5和圖6可以看出，射流速度隨時間的變化具有明顯的拐點，且各個特征點的速度變化趨勢基本一致，因此結合底部射流區的A1點在速度拐點時刻的對應位置來分析射流和二次噴濺的形成機理。

圖6 底部射流區速度隨時間的變化Fig.6 Spray velocity varies with time in the bottom jet area

圖7 為t1時刻（入水后20ms）的流場速度云圖，可以看出，此時結構距離地面高度為48.3 mm，擠壓作用使結構邊緣水粒子的水平速度明顯增大，因此t1時刻為底部射流的產生時刻，此時A1點位于結構底部內側。

以射流區域的A1粒子為例說明底部射流發展過程。從圖6可知，t2時刻（入水后43ms）為其速度拐點，追蹤A1點位置可見該時刻粒子A1運動到結構邊緣位置，如圖8（a）所示。因此A1從t1～t2時刻受到結構的擠壓作用，水平方向速度由于擠壓作用在t1時刻后迅速增大。當水粒子脫離結構邊緣后，擠壓作用消失，水平方向速度減小，豎直方向的速度隨即增大。

從圖6（a）可以看出，t3時刻（入水后60ms）后豎直方向速度緩慢減小，如圖8（b）中點A1所示。可以看出，t3時刻A1粒子到達濺水根部的圓弧中點位置，豎直方向的速度達到最大值，當粒子繼續向上濺起，豎直方向速度減小并趨于穩定。

綜上所述，t1～t3時刻為射流產生和發展階段，與圖5（b）的二次噴濺速度隨時間變化比較，可以看出，該階段二次噴濺區水粒子速度再次增大并達到峰值，說明二次噴濺的形成主要是由于底部射流水平方向速度轉化為二次噴濺區粒子豎直方向的速度，將結構與地面擠壓產生的能量由內向外傳遞。

3.5 水深對淺水噴濺的影響

為了研究水深對淺水噴濺的影響，保持其他條件不變，將水深減小到40mm和16mm，建立淺水沖擊模型。

圖7 t1時刻速度云圖Fig.7 Velocity contour at t1

圖8 水粒子A的位置示意圖Fig.8 Position of water particle A

首先，比較水深對底部射流的影響，在圖4所示底部射流區A4位置取水粒子。圖9為不同水深下，A4位置水粒子水平方向的速度變化。可以看出，隨水深減小，底部射流的水平速度增大。水深為16mm時，底部射流產生時間明顯早于其他兩種水深情況，說明水深較小時，入水初期就會產生擠壓作用。

圖10為不同水深時二次噴濺區域水粒子的豎直方向速度比較，水粒子的水平位置與圖4的C1一致，深度均為水深的1／2。可以看出，隨水深減小，二次噴濺豎直方向的速度明顯增大，說明當水深減小，濺水受到的擠壓作用更加明顯，二次噴濺更為劇烈。

4 淺水噴濺抑制結構機理分析

4.1 抑制結構模型建立

翻邊結構是最常見的淺水噴濺抑制結構之一，為了研究其噴濺抑制機理，采用fokker100機型［8］使用的拱形翻邊結構截面，建立淺水沖擊模型，如圖11所示。可以看出，截面左面幾何構型為拱形翻邊結構，右面為圓弧結構，水深為16mm，粒子間距為1mm，入水速度為7m／s。結構為彈性結構，采用mooney-rivlin橡膠材料模型，材料參數參考文獻［17］的輪胎胎體參數，流體部分仍采用NULL材料模型和Gruneisen狀態方程。

圖9 底部射流速度Fig.9 Velocity of bottom jet

圖10 二次噴濺速度Fig.10 Velocity of secondary spray

4.2 彈性結構淺水沖擊噴濺區域劃分

圖12 為相同水深時，彈性結構截面淺水沖擊的噴濺形態與剛性長方體結構的比較。可以看出，彈性結構由于入水初期為圓弧形，因此初始噴濺與水平面夾角較小，二次噴濺和底部射流部分與剛性長方體結構基本一致，噴濺也可同樣劃分為三個部分。對其各部分噴濺的來源位置進行分析，可將彈性圓弧結構淺水噴濺區域進行劃分，如圖13所示，初始噴濺區為輪胎邊緣下方寬52mm，深6mm的水層；擠壓區為輪胎底部中間位置的30mm寬水域；二次噴濺區為初始噴濺區域周邊區域。

4.3 拱形翻邊構型噴濺抑制機理

圖11 帶有拱形翻邊的截面入水模型Fig.11 Water entry model with arched chine

圖12 兩種截面結構入水噴濺形態比較Fig.12 Comparison of the water spray of two cross section structures

圖13 輪胎濺水分區和標記Fig.13 Region of tire water spray

圖14 拱形翻邊結構濺水抑制效果Fig.14 Spray suppression effect of arched chine

圖14 為拱形翻邊結構的噴濺形態結果和無翻邊結構的噴濺形態比較，可以看出，初始噴濺受到翻邊結構的作用向下發生偏折，二次噴濺角度由無翻邊側的56.4°經翻邊作用后變為31.5°，部分底部射流滯留于翻邊結構下。

按圖13的速度測點位置標記水粒子，將有無翻邊結構各部分噴濺的速度進行比較，如圖15所示。圖15（a）為初始噴濺豎直方向速度比較，可以看出，拱形翻邊有效阻擋了初始噴濺，使初始噴濺的豎直方向濺水速度由15m／s變為負值，最后穩定在3m／s左右。圖15（b）為初始噴濺水平方向速度比較，可以看出，水平方向濺水速度由8m／s減小為1.2m／s。圖16分別為二次噴濺豎直方向速度和水平方向速度比較，可以看出，二次噴濺豎直方向速度沒有明顯變化，水平方向速度增大，因此噴濺與地面的夾角減小。

傳統拱形翻邊結構對淺水噴濺的抑制機理為阻擋初始噴濺，減小底部射流，可以改變二次噴濺的角度，但對二次噴濺豎直方向速度的抑制作用較小。

5 新型淺水噴濺抑制結構

5.1 濺水抑制結構設計

為了提高翻邊結構的效果，本文借鑒了水上飛機抑波槽的設計，在傳統拱形結構的下邊緣增加導流凹槽，并在端部增加彈性的阻水邊條，如圖17所示，旨在更好地阻擋初始噴濺，并利用導出的初始噴濺和底部射流有效減小二次噴濺。

圖15 初始噴濺速度比較Fig.15 Comparison of spray velocity of initial spray

工程應用說明水上飛機抑波槽在深水時具有較好的噴濺抑制效果，而其對淺水噴濺的抑制效果和抑制機理并不清楚。因此，下文對新型翻邊對淺水噴濺的抑制效果進行驗證分析。

5.2 濺水抑制效果比較

圖18為相同時刻兩種噴濺抑制結構的效果比較，左邊為拱形翻邊構型，右邊為新構型。從噴濺形態上可以看出，新構型的二次噴濺角度和高度均明顯減小。

按照3.3節的方法，選取二次噴濺區域的水粒子，進一步分析其速度變化。圖19（a）為二次噴濺豎直方向的速度比較，可以看出，新構型翻邊二次噴濺豎直方向的速度峰值比拱形翻邊降低約42%；圖19（b）為水平方向的速度比較，可以看出，水平方向速度減小約14.8%。因此，新構型翻邊可以抑制二次噴濺速度，并減小二次噴濺角度，對于淺水噴濺具有良好的抑制效果。

圖16 二次噴濺速度比較Fig.16 Comparison of spray velocity of secondary spray

圖17 新型淺水噴濺抑制結構示意圖Fig.17 Sketch map of the new suppression structure

圖18 濺水形態比較Fig.18 Comparison of spray form

5.3 新構型翻邊噴濺抑制機理

從圖19（a）可以看出，二次噴濺水流速度在4ms～6ms之間有減小和波動，該時段的噴濺形成過程如圖20所示。可以看出，初始噴濺形成后，由拱形結構導流進入導流凹槽，在導流凹槽中分為兩部分，一部分在凹槽中形成反向漩渦，如圖21所示，反向漩渦與導出的底部射流相互作用，減小了二次噴濺和底部射流的動能；另一部分由阻水邊條向側下導出，與二次噴濺相互作用后以較小的速度側向排出，進一步減小了二次噴濺速度，同時初始噴濺速度也有所降低。

圖19 新構型翻邊二次噴濺速度比較Fig.19 Comparison of spray velocity of secondary spray

圖20 新型結構濺水形成過程Fig.20 Splash formation process of new suppression structure

圖21 速度矢量圖Fig.21 Speed vector diagram

6 結 論

（1）建立了SPH／FEM淺水沖擊模型，并與文獻的試驗結果進行比較，將濺水分為初始噴濺、二次噴濺和底部射流，通過速度分析說明二次噴濺主要是結構底部對水層的擠壓作用和底部射流能量的向外傳遞共同作用引起的。

（2）分析了水深對二次噴濺和底部射流的影響，當水深減小時，擠壓作用更為明顯，二次噴濺和底部射流速度明顯增大。

（3）研究表明，拱形翻邊結構對淺水噴濺各部分的作用為導出初始噴濺，抑制底部射流，減小二次噴濺角度，但無法有效降低二次噴濺豎直方向的速度。

（4）提出了一種具有導流凹槽和阻水邊條的新型噴濺抑制構型，可在翻邊開槽處形成反向漩渦，從而更有效地抑制二次噴濺。