裴志博，楊秀峰，王曉亮

（北京理工大學(xué)宇航學(xué)院，北京 100081）

1 引 言

?？諈f(xié)同作業(yè)是21 世紀(jì)發(fā)展中的重要作業(yè)形式[1]，跨域無人航行器可能是將來執(zhí)行跨域多目標(biāo)任務(wù)的重要裝備和多種不同類型傳感器的搭載平臺(tái)?？缬驘o人航行器一般是指既能在水下航行，又能在空中飛行，且能跨越水-空氣兩種流體交界面正常作業(yè)的無人航行器。早在20 世紀(jì)30年代，蘇聯(lián)就開始了跨域航行器的探索，后來美國也開展了相關(guān)的研究工作[2]。近年來，美、英等國研制了多款不同類型的跨域無人航行器，其中美國的“鷹鰩”固定翼跨域航行器已經(jīng)進(jìn)行了各個(gè)航行過程的樣機(jī)試驗(yàn)[3]。我國也開展了跨域航行器的研究，主要有航行器模型的仿真驗(yàn)證、原理樣機(jī)的研制和水空過渡等關(guān)鍵技術(shù)的驗(yàn)證等[4]。

跨域無人航行器的跨水空交界面包括“出水入空”（即出水）和“由空入水”（即入水）兩種形式。航行器入水過程非常復(fù)雜，涉及航行器與自由液面的相互作用、液面變形和破碎、航行器的運(yùn)動(dòng)穩(wěn)定性等力學(xué)問題。由于水和空氣兩種流體的性質(zhì)差別較大，比如密度相差約800 倍，這導(dǎo)致航行器在入水時(shí)受到非常大的壓力變化。除跨域航行器入水外，入水相關(guān)問題也廣泛存在于日常生活、科學(xué)研究和國防工程等領(lǐng)域，比如跳水、飛機(jī)水面降落、魚雷以及各種彈體入水等。由于對(duì)工程應(yīng)用具有重要意義，入水問題一直受到廣泛關(guān)注。

在物體入水問題的研究中，試驗(yàn)是非常重要的一種研究手段。在相關(guān)的試驗(yàn)研究中，低速入水試驗(yàn)涉及的物體類型較多，王瑞琦等[5]在試驗(yàn)中用高速攝像機(jī)記錄了不同頭型彈丸低速垂直入水過程；陸響暉等[6]對(duì)流線型和平頭型低速回轉(zhuǎn)體在不同初速度下垂直入水過程進(jìn)行了試驗(yàn)研究；余德磊等[7]開展了回轉(zhuǎn)體低速串聯(lián)入水試驗(yàn)。相比之下，高速入水的試驗(yàn)大多是關(guān)于尺寸比較小的彈體入水試驗(yàn)。其中，張偉等[8]進(jìn)行了高速彈體入水的試驗(yàn)，研究了不同工況下的彈道軌跡和空泡變化，Jiang 等[9]對(duì)彈體在不同入射角度和通氣量情況下的入水過程進(jìn)行了試驗(yàn)研究，分析了入射角和通氣量與空腔演變和閉合的關(guān)系。

隨著高性能計(jì)算機(jī)和數(shù)值方法的快速發(fā)展，數(shù)值模擬已經(jīng)成為研究科學(xué)和工程問題的一種重要方法，并且在物體入水的研究中得到越來越多的應(yīng)用。劉雙等[10]基于STAR-CCM+和Fluent，結(jié)合網(wǎng)格重構(gòu)、重疊網(wǎng)格、界面捕捉等技術(shù)，建立了相應(yīng)的數(shù)值計(jì)算模型，分析研究了二維楔形體的自由入水過程。Sun 等[11]發(fā)展了一種二維邊界元的方法，用來模擬水平圓柱的入水問題。Zhu等[12]應(yīng)用CIP（Constrained Interpolation Profile）方法模擬研究了水平圓柱垂直入水的問題。王占瑩等[13]使用VOF（Volume of Fluid）多相流模型模擬了回轉(zhuǎn)體垂直入水的過程，研究分析了回轉(zhuǎn)體在不同速度入水時(shí)不同流體動(dòng)力特性以及超空泡不同的演化過程。Sun 等[14]基于CFD、FEM和FSI 等方法，使用Star-CCM+和ABAQUS 聯(lián)合仿真研究了圓柱體高速入水過程中圓柱體的運(yùn)動(dòng)、變形和入水沖擊荷載的特征以及入水產(chǎn)生的空化氣泡的演變過程。大多傳統(tǒng)的數(shù)值模擬主要是基于網(wǎng)格類方法，比較難以處理入水過程中涉及的交界面大變形和破碎以及物體的運(yùn)動(dòng)等問題。

光滑粒子法（Smoothed Particle Hydrodynamics，SPH）是一種拉格朗日型無網(wǎng)格粒子方法[15-17]，具有拉格朗日方法準(zhǔn)確描述物質(zhì)界面的優(yōu)勢，而且避免了網(wǎng)格類方法中因網(wǎng)格變形導(dǎo)致的畸變等問題，適用于模擬自由表面、交界面、運(yùn)動(dòng)邊界和大變形等問題[18]。近年來，SPH 方法被廣泛應(yīng)用于自由表面流[19-20]、顆粒流[21-23]、液體晃動(dòng)[24-25]以及液滴撞擊[26-27]等問題。本文應(yīng)用SPH 方法模擬研究細(xì)長體入水問題。

2 數(shù)值方法

2.1 控制方程

流體控制方程采用拉格朗日形式的納維-斯托克斯（Navier-Stokes, N-S）方程：

式中，ρ是密度，t是時(shí)間，u是速度，p是壓力，μ是動(dòng)力粘度，g是重力加速度。

在SPH 方法中，流體一般作為弱可壓縮處理，本文使用的狀態(tài)方程：

其中，ρr是參考密度，c是數(shù)值聲速。為了保持流體的可壓縮性小于1%，數(shù)值聲速的取值為流體最大速度的10 倍以上。

2.2 SPH 方法

在SPH 方法中，對(duì)關(guān)于空間位置r的函數(shù)A（r） ，在點(diǎn)a處的值可以通過下面的積分求出：

式中，W(ra-r,h)是光滑函數(shù)，也稱為核函數(shù)；h是控制核函數(shù)寬度的光滑長度；dV是體積微元；Ω 是積分域。

當(dāng)流體離散成SPH 粒子時(shí)，式(4)中的積分可用粒子求和代替：

式中，下角標(biāo)a和b表示粒子，粒子a對(duì)應(yīng)的物理量有位置ra、質(zhì)量ma、密度ρa(bǔ)、體積ma/ρa(bǔ)，Aa是坐標(biāo)為ra的粒子a處的函數(shù)值，Wab=W(ra-rb,h)。根據(jù)式(5)可以得出函數(shù)A的梯度的離散形式：

式中，?aWab=Wa′b rab/rab表示粒子a的核函數(shù)在粒子b處的梯度；rab=ra-rb是粒子a到b的位移矢量；rab=|rab|是位移矢量rab的模，即粒子a與b之間的距離。式(5)和式(6)中的求和僅限于光滑函數(shù)范圍內(nèi)的粒子。

根據(jù)式(4)和式(5)可將方程(1)和(2)離散成SPH 形式：

式中，uab=ua-ub，η=0.01h2，Π 是由Monaghan[16]給出的人工粘性項(xiàng)：

式中，α取值為0.05，ρa(bǔ)b=(ρa(bǔ)+ρb)/2，μab的表達(dá)式為：

在SPH數(shù)值計(jì)算的過程中，隨著時(shí)間的推進(jìn)，粒子的密度誤差會(huì)不斷增加，導(dǎo)致壓強(qiáng)產(chǎn)生較大的波動(dòng)和震蕩。為了減少誤差，采用密度重值法每隔20 個(gè)時(shí)間步進(jìn)行一次密度修正：

式中，是修正后的粒子密度，Vb=mb/ρb是粒子的體積。

2.3 光滑函數(shù)

光滑函數(shù)決定了SPH方法中相鄰粒子間相互作用的關(guān)系，因此光滑函數(shù)的選取對(duì)計(jì)算結(jié)果的精準(zhǔn)度有著至關(guān)重要的作用，有多種核函數(shù)可供選擇[28]。三次樣條核函數(shù)是目前使用較為廣泛的鐘形核函數(shù)：

式中，s=r/h，αd稱為歸一化因子，角標(biāo)d表示問題的空間維度，在一維、二維、三維空間中歸一化因子的值分別取為1/6h、5/(14πh2)、1/(4πh3)。

采用傳統(tǒng)的鐘形核函數(shù)容易導(dǎo)致應(yīng)力不穩(wěn)定問題，文獻(xiàn)[29]提出了可以有效解決應(yīng)力不穩(wěn)定問題的雙曲型核函數(shù)：

考慮到不同光滑函數(shù)的優(yōu)勢，在本文的數(shù)值計(jì)算過程中，涉及核函數(shù)本身求和時(shí)使用三次樣條核函數(shù)；當(dāng)計(jì)算核函數(shù)導(dǎo)數(shù)時(shí)使用雙曲型核函數(shù)。

2.4 粒子遷移

在數(shù)值模擬過程中，SPH 粒子通常是復(fù)雜無序的，導(dǎo)致粒子分布的均勻性較差，降低計(jì)算的精度和穩(wěn)定性。為了使粒子分布更加均勻，Chaniotis 等[30]提出了粒子重分技術(shù)，Xu 等[31]提出了粒子遷移技術(shù)。Lind 等[32]基于菲克定律對(duì)粒子遷移技術(shù)進(jìn)行了推廣。本文采用粒子遷移技術(shù)對(duì)粒子位置進(jìn)行修正：

式中，δ ra是粒子修正后的位置相對(duì)于原位置的位移矢量：

這里，D為粒子的修正系數(shù)，取值為0.5h2。C是粒子濃度：

粒子的濃度梯度為：

由式(15)～(17)可以得出粒子位移的修正值[32]為：

Lind 等[32]在粒子濃度梯度公式中加入了類似人工應(yīng)力的修正項(xiàng)，目的是避免粒子配對(duì)不穩(wěn)定。由于本文在計(jì)算核函數(shù)梯度時(shí)使用的是雙曲核函數(shù)，已經(jīng)有效地避免了粒子配對(duì)不穩(wěn)定的問題，因此無須引入人工應(yīng)力。另外，粒子遷移技術(shù)在復(fù)雜的自由表面附近會(huì)導(dǎo)致非物理的粒子擴(kuò)散，因此本文僅對(duì)內(nèi)部粒子使用粒子遷移技術(shù)，而自由表面附近的粒子位置則保持不變。

在上述粒子位移修正過程中，需要識(shí)別自由表面的粒子信息，但是在SPH 方法中自由表面粒子的信息是未知的。目前，已經(jīng)發(fā)展了很多自由表面粒子的搜索和識(shí)別方法[33-35]，本文采用Abbas 等[35]提出的自由表面粒子判據(jù)，即?·r＜γ（γ取值為1.8）。

2.5 自由表面

在細(xì)長體入水問題中，主要涉及三類邊界條件：自由液面、固體壁面以及細(xì)長體表面。在光滑粒子法中，自由液面邊界不需要做特殊處理，但是由于自由液面的外部區(qū)域沒有粒子，為了避免自由面附近因密度梯度較大導(dǎo)致的壓力大于零，將自由表面附近的密度設(shè)置為參考密度ρr，根據(jù)狀態(tài)方程(3)保證自由面附近的壓力為零，避免數(shù)值波動(dòng)。

2.6 固體壁面

對(duì)固體壁面邊界的處理方法有很多，比如斥力邊界、反射邊界等，本文使用的是動(dòng)力邊界條件，即在壁面附近設(shè)置壁面虛粒子，在計(jì)算壁面粒子與流體粒子相互作用時(shí)，其物理性質(zhì)與流體粒子一致。固體壁面粒子與流體粒子的相互作用為[33]：

2.7 物體運(yùn)動(dòng)方程

在物體入水過程中，物體除受重力作用外，還有與水之間的相互作用，是典型的流固耦合問題。對(duì)于流體與物體之間的相互作用，先通過式(19)計(jì)算固體邊界粒子與流體粒子之間的相互作用，再通過求和得到物體與流體之間的相互作用，通過牛頓定律以及剛體運(yùn)動(dòng)方程求解物體運(yùn)動(dòng)。

物體在自身重力和流體外力的作用下，不考慮物體變形時(shí)，滿足牛頓運(yùn)動(dòng)定理：

其中，流體作用力和力矩可通過下式求出：

式中，a、b分別代表流體粒子和物體粒子，fab是流體粒子受到物體粒子的作用，F(xiàn)ba是單個(gè)物體粒子受到的流體粒子作用力，F(xiàn)是細(xì)長體整體受到的總的流體作用力，M是細(xì)長體的總質(zhì)量，Lc是物體所受力對(duì)質(zhì)心的力矩，Nc是物體相對(duì)于質(zhì)心的轉(zhuǎn)動(dòng)慣量。

3 數(shù)值模型驗(yàn)證

本節(jié)通過對(duì)比水平圓柱入水模擬結(jié)果與試驗(yàn)結(jié)果[36]，對(duì)前文所述數(shù)值方法進(jìn)行驗(yàn)證。模擬中各物理參數(shù)與試驗(yàn)一致。對(duì)比試驗(yàn)來自Greenhow等[12,36]：直徑為0.11 m 的圓柱水平放置在距水面高0.5 m 的位置自由下落并入水。試驗(yàn)分別采用了零浮力圓柱（密度與水的密度相同）和半浮力圓柱（密度為水的密度的一半）兩種圓柱。忽略空氣阻力，根據(jù)自由落體公式，可得圓柱體入水時(shí)的速度為v=2.955 m/s。

根據(jù)試驗(yàn)設(shè)置數(shù)值模型，各參數(shù)與上述試驗(yàn)保持一致，模擬從圓柱體剛接觸水面開始。數(shù)值水槽寬為1.2 m，水深為1 m，如圖1所示。水的密度為1000 kg/m3，零浮力圓柱密度為1000 kg/m3，半浮力圓柱密度為500 kg/m3。初始粒子間距為5 mm，時(shí)間步長為5.0×10-6s，整個(gè)模型共使用50354 個(gè)粒子，其中水粒子48000 個(gè)，圓柱粒子176 個(gè)，外部壁面粒子2178 個(gè)。粒子間距越小，模擬所需的粒子數(shù)就越多，數(shù)值精度越高，能夠捕捉的流場細(xì)節(jié)也越多；同時(shí)，粒子間距越小，時(shí)間步長就越小，模擬所需的計(jì)算時(shí)間就越長。圖2和圖3分別是半浮力圓柱和零浮力圓柱體入水過程中液面形態(tài)的演化過程以及與試驗(yàn)結(jié)果的對(duì)比。從圖中可以觀察到，兩種圓柱入水過程中液面形態(tài)的演變與試驗(yàn)結(jié)果一致。無論是在圓柱入水的初始階段水面高度的變化和水面飛濺方向，還是后續(xù)階段液面運(yùn)動(dòng)形態(tài)的演變都與試驗(yàn)結(jié)果吻合。另外，從圖2和圖3定性來看，圓柱體在入水過程中入水深度隨時(shí)間變化的模擬結(jié)果與試驗(yàn)結(jié)果一致。

圖1 水平圓柱入水?dāng)?shù)值模型Fig.1 Initial condition of a circular cylinder water-entry flow

圖2 半浮力圓柱入水時(shí)液面形態(tài)的演化過程（左列為數(shù)值模擬結(jié)果，右列是試驗(yàn)結(jié)果[12,36]）Fig.2 Free surface deformation of a half-buoyant cylinder(left: SPH simulation, right: experiment[12,36])

圖3 零浮力圓柱入水時(shí)液面形態(tài)的演化過程（左列為數(shù)值模擬結(jié)果，右列是試驗(yàn)結(jié)果[12,36]）Fig.3 Free surface deformation of a neutrally buoyant cylinder (left: SPH simulation, right: experiment[12,36])

圖4是半浮力圓柱入水深度隨時(shí)間變化的模擬結(jié)果與試驗(yàn)結(jié)果的對(duì)比。模擬結(jié)果與試驗(yàn)值符合的非常好，雖然最后一個(gè)試驗(yàn)值與模擬結(jié)果相差稍大。圖5展示的是零浮力圓柱體入水深度隨時(shí)間變化的模擬值與試驗(yàn)值的對(duì)比。模擬結(jié)果與試驗(yàn)結(jié)果基本一致。

圖4 半浮力圓柱入水深度變化[12,36]Fig.4 Penetration depth of a half-buoyant cylinder[12,36]

圖5 零浮力圓柱入水深度變化[12,36]Fig.5 Penetration depth of a neutrally buoyant cylinder[12,36]

與試驗(yàn)的對(duì)比結(jié)果表明，SPH 方法在模擬物體入水問題上具有可行性和可靠性。

4 模擬結(jié)果及分析

本節(jié)基于跨域無人航行器的入水問題，采用SPH 方法模擬細(xì)長體的垂直入水過程。由于細(xì)長體尺寸較小且入水速度較低，所以忽略空氣的影響。計(jì)算區(qū)域的初始設(shè)置見圖6。細(xì)長體的寬和高分別為0.05 m 和0.2 m，數(shù)值水槽的寬為2.2 m，水深為2 m。初始粒子間距取為5 mm，時(shí)間步長選取為5.0×10-6s，整個(gè)模擬系統(tǒng)共使用180242個(gè)粒子進(jìn)行離散和描述，其中水粒子包括176000個(gè)，細(xì)長體邊界粒子有264 個(gè)以及3978 個(gè)動(dòng)力壁面邊界粒子。本節(jié)分別模擬了不同密度的細(xì)長體以不同速度的入水過程，觀察細(xì)長體入水時(shí)水面形態(tài)的演變以及細(xì)長體入水深度與時(shí)間的關(guān)系。

圖6 細(xì)長體入水?dāng)?shù)值模型Fig.6 Initial condition of a slender body water-entry flow

4.1 結(jié)構(gòu)物密度的影響

首先，在初始入水速度1 m/s 和3 m/s 的情況下模擬了密度分別為1000 kg/m3、2700 kg/m3和7500 kg/m3的細(xì)長體入水過程，在細(xì)長體入水過程中自由液面的變化如圖7和圖8所示。我們可以觀察到，當(dāng)細(xì)長體剛開始接觸水面時(shí)，自由水面被打破原有的穩(wěn)定狀態(tài)，水面破碎成水滴并發(fā)生一定的飛濺，由于細(xì)長體的擠壓作用，水從細(xì)長體處開始向相對(duì)遠(yuǎn)處運(yùn)動(dòng)，細(xì)長體附近的水面有一定程度的上升，隨后水體在重力和細(xì)長體的共同作用下開始演變?yōu)椴煌男螒B(tài)。圖7和圖8分別描述了密度為1000 kg/m3、2700 kg/m3和7500 kg/m3的細(xì)長體以1 m/s 和3 m/s 的初速度入水時(shí)水體形態(tài)的變化過程，從圖中我們可以看出相同速度下結(jié)構(gòu)體從開始入水到完全進(jìn)入水中的這個(gè)過程中（對(duì)應(yīng)圖7、圖8中0.1 s 前）無論是水面高度的變化還是水面開口的形態(tài)都不隨細(xì)長體密度的改變而產(chǎn)生明顯的變化，這說明在細(xì)長體以一定的速度入水的初始階段，自身密度對(duì)水體的運(yùn)動(dòng)和自由液面的形態(tài)變化幾乎沒有影響。但是當(dāng)細(xì)長體完全進(jìn)入水面以下后，自身密度的增大會(huì)對(duì)水體向周圍的運(yùn)動(dòng)有一定的增強(qiáng)效應(yīng)，且入水速度越大，密度的影響越明顯。另外，在細(xì)長體入水深度的變化中可以發(fā)現(xiàn)相似的規(guī)律，如圖10和圖11所示。在結(jié)構(gòu)體以同一初速入水的最初始階段，相同時(shí)刻不同密度結(jié)構(gòu)物的入水深度是比較相近的，說明在細(xì)長體開始撞擊入水的過程中物體密度對(duì)細(xì)長體與水體自由液面的相互作用影響較小，但隨著入水過程的進(jìn)行，密度的對(duì)細(xì)長體入水深度的影響越來越明顯，比較兩圖，可以發(fā)現(xiàn)入水速度的增加會(huì)放大密度的影響，故在細(xì)長體開始撞擊入水的這個(gè)階段，速度對(duì)結(jié)構(gòu)物和水體之間的相互作用具有比較大的影響。

圖7 不同密度的細(xì)長體以速度v=1 m/s 入水過程中水形態(tài)的變化（從上到下結(jié)構(gòu)體密度依次為1000 kg/m3、2700 kg/m3 和7500 kg/m3）Fig.7 Free surface deformation of slender body with different densities at the speed v=1 m/s (The density of the structure from top to bottom is 1000 kg/m3,2700 kg/m3 and 7500 kg/m3)

圖8 不同密度的細(xì)長體以速度v=3 m/s 入水過程中水形態(tài)的變化（從上到下結(jié)構(gòu)體密度依次為1000 kg/m3、2700 kg/m3 和7500 kg/m3）Fig.8 Free surface deformation of slender body with different densities at the speed v=3 m/s (The density of the structure from top to bottom is 1000 kg/m3,2700 kg/m3 and 7500 kg/m3)

圖9 密度為2700 kg/m3 的細(xì)長體分別以1 m/s（上）、3 m/s（中）和5 m/s（下）的初速度入水過程中水形態(tài)的變化Fig.9 Free surface deformation of slender body with densities of 2700 kg/m3 at the different initial speeds(The speed from top to bottom is 1 m/s, 3 m/s and 5 m/s)

圖10 不同密度細(xì)長體以1 m/s 的初速入水時(shí)入水深度的變化Fig.10 Penetration depth of slender body with different densities enter water at an initial speed of 1 m/s

圖11 不同密度的細(xì)長體以3 m/s 的初速入水時(shí)入水深度的變化Fig.11 Penetration depth of slender body with different densities enter water at an initial speed of 3 m/s

4.2 入水速度的影響

接下來研究分析了在入水過程中入水速度對(duì)細(xì)長體自身運(yùn)動(dòng)和水運(yùn)動(dòng)形態(tài)演變的影響。圖9是密度為 2700 kg/m3的細(xì)長體在初速分別為1 m/s、3 m/s 和5 m/s 時(shí)入水過程的模擬結(jié)果。從上到下來看，相同時(shí)間對(duì)應(yīng)的狀態(tài)圖，水面的缺口逐漸增大，細(xì)長體相對(duì)于水面的位置也逐漸下降。說明細(xì)長體入水的初始速度越大，水面被撞開的缺口越大，水面高度上升得越多，之后細(xì)長體下降的越快。從圖12和圖13中可以觀察到，不同曲線之間的差距從開始逐漸加大到一定時(shí)間后保持不變，說明入水速度在細(xì)長體剛開始進(jìn)入水中時(shí)就對(duì)細(xì)長體的運(yùn)動(dòng)有較大的影響。初始入水速度越大，細(xì)長體進(jìn)入水中的速度越快，相同時(shí)刻細(xì)長體的入水就越深，但是不同初始速度的細(xì)長體入水深度的差距不會(huì)持續(xù)增加，而是到達(dá)一定時(shí)間后細(xì)長體入水深度在各情況之間的差距保持在一個(gè)相對(duì)穩(wěn)定的范圍內(nèi)。

圖12 密度為1000 kg/m3 的細(xì)長體以不同初速入水時(shí)入水深度的變化Fig.12 Penetration depth of slender body with a density of 1000 kg/m3 enters water at different initial speeds

圖13 密度為2700 kg/m3 的細(xì)長體以不同初速入水時(shí)入水深度的變化Fig.13 Penetration depth of slender body with a density of 2700 kg/m3 enters water at different initial speeds

5 結(jié) 論

針對(duì)跨域無人航行器入水問題，本文應(yīng)用SPH 方法模擬研究了細(xì)長體垂直入水過程。模擬結(jié)果表明：在入水過程的初始階段，入水速度對(duì)細(xì)長體與液面的相互作用有較大影響：入水速度越大，液面形態(tài)的改變越大，細(xì)長體下降的速度越快。初始入水速度對(duì)細(xì)長體入水過程的影響會(huì)逐漸減弱。在入水過程的后期階段，細(xì)長體密度的影響逐漸增大：密度越大，細(xì)長體在入水后期下落的速度越快，而且入水速度的增加對(duì)這種影響有增強(qiáng)的效果，這主要是因?yàn)殡S著細(xì)長體密度的增大，細(xì)長體自身重力的影響也增大。后續(xù)工作將研究復(fù)雜結(jié)構(gòu)的航行器以不同角度入水等問題。