基于SPH／FEM方法的淺水沖擊噴濺及其抑制結(jié)構(gòu)研究

任選其， 徐 緋， 張顯鵬， 楊 揚(yáng)

（西北工業(yè)大學(xué) 航空學(xué)院，西安710072）

1 引 言

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

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

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

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

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

2 SPH／FEM耦合方法

2.1 SPH方法基本方程

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

質(zhì)量守恒

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

動量守恒

能量守恒

式中 Pi為i點(diǎn)的各向同性壓力，mj為第j個粒子的質(zhì)量，為第i個粒子在β 方向的坐標(biāo)，和分別為第i個粒子的應(yīng)力張量和應(yīng)變張量，＝（－）為兩個粒子之間的相對速度在β方向上的分量，N為光滑長度范圍內(nèi)的粒子數(shù)，ηi為粘性系數(shù)

2.2 接觸算法

實(shí)際計(jì)算中SPH和FEM相對獨(dú)立計(jì)算，每個計(jì)算歩通過罰函數(shù)接觸算法傳遞接觸狀態(tài)和作用力，保證SPH水粒子在計(jì)算過程中不會發(fā)生穿透，同時實(shí)現(xiàn)兩者力的相互作用［14，15］。

其中，接觸面處的法向接觸力表達(dá)式為

式中ki為接觸剛度，fs為ki的比例因子，li為SPH水粒子相對于FEM結(jié)構(gòu)網(wǎng)格單元的穿透深度，n為結(jié)構(gòu)單元的法向單位矢量。

接觸面的切向接觸力表達(dá)式為

式中μk為水粒子在網(wǎng)格單元表面的摩擦系數(shù)。

3 淺水沖擊噴濺機(jī)理研究

3.1 長方體結(jié)構(gòu)淺水沖擊模型

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

流體部分采用NULL材料模型，NULL材料模型不考慮偏應(yīng)力的計(jì)算，通常用于定義流體。由于液體受到結(jié)構(gòu)底部和地面之間的擠壓作用，因此考慮了液體的可壓縮性，液體的壓力和體積的關(guān)系通過Gruneisen狀態(tài)方程進(jìn)行描述：

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

圖1 長方體結(jié)構(gòu)淺水沖擊模型Fig.1 Schematic diagram of shallow water impact model of rectangular structure

3.2 計(jì)算結(jié)果與試驗(yàn)比較

入水初速為2m／s，結(jié)構(gòu)重量為6.2kg時，濺水形態(tài)計(jì)算結(jié)果與Bukree試驗(yàn)結(jié)果［5］的對比，如圖2所示?？梢钥闯?，兩者吻合良好，說明本文構(gòu)建的SPH／FEM數(shù)值模型可以較好地模擬淺水沖擊噴濺現(xiàn)象。

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

表1 流體材料狀態(tài)方程參數(shù)Tab.1 State equation parameters of fluid material

3.3 淺水噴濺來源區(qū)域及其速度分布

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

首先，如圖3（a）所示，對二次噴濺區(qū)的水粒子進(jìn)行標(biāo)記，分析各部分濺水的初始來源區(qū)域，如圖3（b）所示，黑色部分為標(biāo)記后的二次噴濺區(qū)域。可以看出，初始噴濺來源于靠近結(jié)構(gòu)邊緣的表層水，二次噴濺主要來源于結(jié)構(gòu)外側(cè)和底部水層，底部射流來源于結(jié)構(gòu)中部區(qū)域。

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

圖2 濺水形態(tài)計(jì)算結(jié)果與試驗(yàn)比較Fig.2 Comparison of calculation results and exprimental results of spray

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

圖4 底部射流區(qū)域的特征點(diǎn)選取Fig.4 Feature points of the bottom jet area

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

圖5 為初始噴濺區(qū)和二次噴濺區(qū)標(biāo)記水粒子豎直方向速度隨時間的變化?？梢钥闯?，在結(jié)構(gòu)與水面接觸后，初始噴濺在短時間內(nèi)速度迅速增加，約5ms時達(dá)到速度峰值8m／s，之后基本保持穩(wěn)定。二次噴濺的濺水速度在結(jié)構(gòu)入水初期2ms內(nèi)迅速增大到約1m／s，其后略有減小，隨后再次增加，到達(dá)峰值2.6m／s后，在重力作用下濺水速度逐漸減小。

圖6（a）為底部射流區(qū)豎直方向速度的變化趨勢，可以看出，表面水粒子速度初始為負(fù)值，大小等于結(jié)構(gòu)的下落速度，隨后速度在短時間內(nèi)明顯增大，到達(dá)速度峰值后緩慢降低。圖6（b）為水平方向速度的變化趨勢，可以發(fā)現(xiàn)，射流區(qū)水粒子的水平方向速度在入水的初始階段較為穩(wěn)定，到達(dá)臨界時刻后迅速增至一個較高的速度峰值，然后迅速減小，最后穩(wěn)定在較低速度。

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

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

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

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

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

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

從圖6（a）可以看出，t3時刻（入水后60ms）后豎直方向速度緩慢減小，如圖8（b）中點(diǎn)A1所示?？梢钥闯?，t3時刻A1粒子到達(dá)濺水根部的圓弧中點(diǎn)位置，豎直方向的速度達(dá)到最大值，當(dāng)粒子繼續(xù)向上濺起，豎直方向速度減小并趨于穩(wěn)定。

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

3.5 水深對淺水噴濺的影響

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

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

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

首先，比較水深對底部射流的影響，在圖4所示底部射流區(qū)A4位置取水粒子。圖9為不同水深下，A4位置水粒子水平方向的速度變化?？梢钥闯觯S水深減小，底部射流的水平速度增大。水深為16mm時，底部射流產(chǎn)生時間明顯早于其他兩種水深情況，說明水深較小時，入水初期就會產(chǎn)生擠壓作用。

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

4 淺水噴濺抑制結(jié)構(gòu)機(jī)理分析

4.1 抑制結(jié)構(gòu)模型建立

翻邊結(jié)構(gòu)是最常見的淺水噴濺抑制結(jié)構(gòu)之一，為了研究其噴濺抑制機(jī)理，采用fokker100機(jī)型［8］使用的拱形翻邊結(jié)構(gòu)截面，建立淺水沖擊模型，如圖11所示?？梢钥闯?，截面左面幾何構(gòu)型為拱形翻邊結(jié)構(gòu)，右面為圓弧結(jié)構(gòu)，水深為16mm，粒子間距為1mm，入水速度為7m／s。結(jié)構(gòu)為彈性結(jié)構(gòu)，采用mooney-rivlin橡膠材料模型，材料參數(shù)參考文獻(xiàn)［17］的輪胎胎體參數(shù)，流體部分仍采用NULL材料模型和Gruneisen狀態(tài)方程。

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

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

4.2 彈性結(jié)構(gòu)淺水沖擊噴濺區(qū)域劃分

圖12 為相同水深時，彈性結(jié)構(gòu)截面淺水沖擊的噴濺形態(tài)與剛性長方體結(jié)構(gòu)的比較?？梢钥闯觯瑥椥越Y(jié)構(gòu)由于入水初期為圓弧形，因此初始噴濺與水平面夾角較小，二次噴濺和底部射流部分與剛性長方體結(jié)構(gòu)基本一致，噴濺也可同樣劃分為三個部分。對其各部分噴濺的來源位置進(jìn)行分析，可將彈性圓弧結(jié)構(gòu)淺水噴濺區(qū)域進(jìn)行劃分，如圖13所示，初始噴濺區(qū)為輪胎邊緣下方寬52mm，深6mm的水層；擠壓區(qū)為輪胎底部中間位置的30mm寬水域；二次噴濺區(qū)為初始噴濺區(qū)域周邊區(qū)域。

4.3 拱形翻邊構(gòu)型噴濺抑制機(jī)理

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

圖12 兩種截面結(jié)構(gòu)入水噴濺形態(tài)比較Fig.12 Comparison of the water spray of two cross section structures

圖13 輪胎濺水分區(qū)和標(biāo)記Fig.13 Region of tire water spray

圖14 拱形翻邊結(jié)構(gòu)濺水抑制效果Fig.14 Spray suppression effect of arched chine

圖14 為拱形翻邊結(jié)構(gòu)的噴濺形態(tài)結(jié)果和無翻邊結(jié)構(gòu)的噴濺形態(tài)比較，可以看出，初始噴濺受到翻邊結(jié)構(gòu)的作用向下發(fā)生偏折，二次噴濺角度由無翻邊側(cè)的56.4°經(jīng)翻邊作用后變?yōu)?1.5°，部分底部射流滯留于翻邊結(jié)構(gòu)下。

按圖13的速度測點(diǎn)位置標(biāo)記水粒子，將有無翻邊結(jié)構(gòu)各部分噴濺的速度進(jìn)行比較，如圖15所示。圖15（a）為初始噴濺豎直方向速度比較，可以看出，拱形翻邊有效阻擋了初始噴濺，使初始噴濺的豎直方向?yàn)R水速度由15m／s變?yōu)樨?fù)值，最后穩(wěn)定在3m／s左右。圖15（b）為初始噴濺水平方向速度比較，可以看出，水平方向?yàn)R水速度由8m／s減小為1.2m／s。圖16分別為二次噴濺豎直方向速度和水平方向速度比較，可以看出，二次噴濺豎直方向速度沒有明顯變化，水平方向速度增大，因此噴濺與地面的夾角減小。

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

5 新型淺水噴濺抑制結(jié)構(gòu)

5.1 濺水抑制結(jié)構(gòu)設(shè)計(jì)

為了提高翻邊結(jié)構(gòu)的效果，本文借鑒了水上飛機(jī)抑波槽的設(shè)計(jì)，在傳統(tǒng)拱形結(jié)構(gòu)的下邊緣增加導(dǎo)流凹槽，并在端部增加彈性的阻水邊條，如圖17所示，旨在更好地阻擋初始噴濺，并利用導(dǎo)出的初始噴濺和底部射流有效減小二次噴濺。

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

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

5.2 濺水抑制效果比較

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

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

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

圖17 新型淺水噴濺抑制結(jié)構(gòu)示意圖Fig.17 Sketch map of the new suppression structure

圖18 濺水形態(tài)比較Fig.18 Comparison of spray form

5.3 新構(gòu)型翻邊噴濺抑制機(jī)理

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

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

圖20 新型結(jié)構(gòu)濺水形成過程Fig.20 Splash formation process of new suppression structure

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

6 結(jié) 論

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

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

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

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