管祥善，孫鵬楠，2，*，李江昊，孫龍泉

(1. 中山大學(xué) 海洋工程與技術(shù)學(xué)院，珠海 519000；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000；3. 哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001)

0 引言

航行體在波浪中的入水問(wèn)題一直是海洋工程中研究的熱點(diǎn)和難點(diǎn)，航行體在波浪條件等復(fù)雜環(huán)境下入水的彈道穩(wěn)定性，也一直是國(guó)防領(lǐng)域中十分關(guān)注的問(wèn)題。航行體在波浪中的入水過(guò)程是一種非常復(fù)雜的流體力學(xué)現(xiàn)象，涉及到水動(dòng)力學(xué)、空氣動(dòng)力學(xué)、剛體動(dòng)力學(xué)等多個(gè)學(xué)科領(lǐng)域[1]。在高速航行體入水過(guò)程中，入水沖擊載荷和運(yùn)動(dòng)響應(yīng)會(huì)影響航行體的運(yùn)動(dòng)軌跡[2-3]，而波浪在自然界廣泛存在且很難避免，波浪的相位、浪高、頻率等物理參數(shù)呈現(xiàn)強(qiáng)隨機(jī)性，也會(huì)對(duì)航行體的入水狀態(tài)和運(yùn)動(dòng)軌跡產(chǎn)生明顯的影響[4]，進(jìn)而決定發(fā)射任務(wù)的成敗[5]。研究真實(shí)波浪條件下航行體的入水過(guò)程可以更好地了解航行體在波浪中的運(yùn)動(dòng)規(guī)律以及波浪對(duì)于航行體運(yùn)行軌跡的影響，可為海洋工程和國(guó)防領(lǐng)域航行體的設(shè)計(jì)和優(yōu)化提供技術(shù)支撐，滿足國(guó)家重大戰(zhàn)略需求。

針對(duì)航行體的入水問(wèn)題，主要采用理論研究[6-7]、試驗(yàn)研究[8-11]和數(shù)值模擬[12]。 隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，數(shù)值模擬逐漸成為解決工程問(wèn)題的主要方法，計(jì)算流體動(dòng)力學(xué) （Computational Fluid Dynamics,CFD）逐漸成為水動(dòng)力學(xué)研究的一種主要手段[13-15]。但是，由于航行體波浪入水問(wèn)題中包含的復(fù)雜物理現(xiàn)象，對(duì)其進(jìn)行準(zhǔn)確數(shù)值仿真始終是CFD領(lǐng)域最具挑戰(zhàn)性的課題之一。目前，大多數(shù)的CFD數(shù)值模擬方法是基于網(wǎng)格類方法[16-18]。楊曉光等[19]采用流體體積多相流模型和動(dòng)網(wǎng)格技術(shù)建立高速入水?dāng)?shù)值模擬方法，并開(kāi)展靜水工況下航行體入水試驗(yàn)，驗(yàn)證了其數(shù)值方法的精度，并基于此數(shù)值方法對(duì)航行體在波浪條件下的高速斜入水過(guò)程進(jìn)行模擬，得到了航行體在不同相位角下的入水軌跡和運(yùn)動(dòng)響應(yīng)結(jié)果。盡管網(wǎng)格法可以精確計(jì)算航行體入水過(guò)程中的水動(dòng)力載荷，但需要提前對(duì)飛濺區(qū)域進(jìn)行網(wǎng)格加密處理才能精確地捕捉大變形液面和飛濺液滴。在模擬自由液面大變形運(yùn)動(dòng)時(shí)，拉格朗日粒子算法更具優(yōu)勢(shì)[20-22]，該算法易于處理大變形和快速移動(dòng)的自由表面，以及飛濺的離散液滴，并能捕捉和生成多相流界面，無(wú)需額外的界面處理，有利于處理多介質(zhì)多物理場(chǎng)的耦合[23-26]。光滑粒子流體動(dòng)力學(xué)（smoothed particle hydrodynamics, SPH）方法作為應(yīng)用最廣泛的粒子類方法，被很多研究者用于模擬航行體的入水問(wèn)題[27-29]。Zhao等[30]采用SPH方法，對(duì)波浪條件下的結(jié)構(gòu)物入水問(wèn)題和流場(chǎng)細(xì)節(jié)進(jìn)行了模擬分析。Yan等[31]采用更新拉格朗日粒子流體動(dòng)力學(xué)方法對(duì)復(fù)雜結(jié)構(gòu)物的入水問(wèn)題進(jìn)行了模擬，并很好地捕捉了剛體入水過(guò)程中周圍形成的水冠、空腔形狀和流場(chǎng)細(xì)節(jié)。航行體在波浪中高速入水的數(shù)值模擬，相對(duì)于靜水工況，存在著波浪建模、液面劇烈變形及液滴飛濺、剛體六自由度運(yùn)動(dòng)響應(yīng)計(jì)算等難點(diǎn)，而針對(duì)波浪的SPH方法建模，也往往采用推板造波的方式，并設(shè)置消波區(qū)來(lái)消除波浪反射，這增加了計(jì)算規(guī)模和復(fù)雜度，并且推板造波需要一定的時(shí)間形成穩(wěn)定的波浪，不適合航行體高速瞬時(shí)入水的模擬。

因此，本文采用SPH方法，提出一種新型的周期性波浪邊界技術(shù)，建立波浪條件下的數(shù)值水池，之后采用四元數(shù)法計(jì)算航行體的六自由度運(yùn)動(dòng)，通過(guò)SPH方法與四元數(shù)法結(jié)合，對(duì)高速航行體在不同波浪相位角下的入水過(guò)程展開(kāi)數(shù)值模擬，得到不同波浪相位角下航行體入水的六自由度運(yùn)動(dòng)軌跡，以期為高速入水航行器的工程設(shè)計(jì)提供計(jì)算工具和理論依據(jù)。

1 數(shù)值計(jì)算方法

1.1 SPH方法基本理論

本文采用Antuono等提出的δ-SPH方法[32]，該方法易于實(shí)現(xiàn)，與標(biāo)準(zhǔn)弱可壓光滑粒子流體動(dòng)力學(xué)（weakly compressible smoothed particle hydrodynamics,WCSPH）相比，其主要優(yōu)點(diǎn)是更具穩(wěn)定性，可避免高頻振蕩下出現(xiàn)非物理能量。

δ-SPH方法的主要控制方程如下：

其中：連續(xù)性方程和動(dòng)量方程右側(cè)最后一項(xiàng)分別為密度耗散項(xiàng)和人工黏性項(xiàng)；D/Dt表示跟隨流體微團(tuán)運(yùn)動(dòng)的導(dǎo)數(shù)，稱之為對(duì)時(shí)間t的物質(zhì)導(dǎo)數(shù)或拉格朗日導(dǎo)數(shù)；ρ、u和r分別表示流體微團(tuán)的密度、速度矢量和位置矢量；g為體積力產(chǎn)生的加速度(本文中指重力加速度)；p、V分別代表壓力和體積；α和δ分別為密度耗散項(xiàng)系數(shù)和人工黏性系數(shù)，本文中皆取0.1；下標(biāo)i、j表示粒子編號(hào)；W為核函數(shù)；h為光滑長(zhǎng)度；c0、ρ0分別為流體聲速和流體參考密度。

1.2 波浪數(shù)值模型

為研究波浪環(huán)境下高速航行體入水的沖擊載荷和運(yùn)動(dòng)軌跡，本研究基于有限水深波浪理論建立波浪數(shù)值模型[33]。波面方程為：

其中：η為波浪在t時(shí)刻x位置的高度；初始化時(shí)刻t為0；H為波高；ω=2π/T為頻率，T為波浪周期；k=2π/L為波數(shù)，L為波長(zhǎng)。L可以根據(jù)色散方程得到[33]：

其中d為靜水深度。速度勢(shì)函數(shù)[19]的表達(dá)式如下：

其中x和y是橫坐標(biāo)和縱坐標(biāo)。根據(jù)速度勢(shì)函數(shù)和水深可以得到波浪某一位置的速度和壓力[33]：

其中u和v分別是x方向和y方向的速度。根據(jù)式（5）計(jì)算的速度和壓力賦予波浪粒子初始速度和壓力值。

1.3 周期性波浪邊界技術(shù)

造波一直是建立波浪數(shù)值水池的關(guān)鍵問(wèn)題，如果采用傳統(tǒng)的推板造波方法，需要在波浪演化尾部加入海綿層消波區(qū)域來(lái)降低波浪的反射，處理相對(duì)復(fù)雜，計(jì)算域也相對(duì)較大，且波浪的穩(wěn)定也需要一定的時(shí)間，因此并不適合用于航行體高速入水工況的計(jì)算。本文提出一種新型周期性波浪邊界技術(shù)，在處理造波問(wèn)題上具有很大的優(yōu)勢(shì)：只需將波浪數(shù)值水池中的SPH粒子按照波浪數(shù)值模型賦予初始速度和壓力，流體域長(zhǎng)度滿足波長(zhǎng)的整數(shù)倍，采用周期性邊界即可實(shí)現(xiàn)波浪的穩(wěn)定流動(dòng)，且波浪的穩(wěn)定不需要時(shí)間，非常適合航行體高速入水工況的計(jì)算。因此，本文在三維波浪數(shù)值水池的邊界處理中，對(duì)波浪的入口和出口區(qū)域采用了周期性邊界技術(shù)，如圖1所示。

當(dāng)流體粒子流出計(jì)算域并進(jìn)入右側(cè)的出口區(qū)域時(shí)，粒子攜帶著所有的物理信息，立刻轉(zhuǎn)移到流體區(qū)域的最左側(cè)，作為計(jì)算域的流入粒子。出口區(qū)域和入口區(qū)域的虛粒子由流體域兩側(cè)的粒子復(fù)制生成，作為計(jì)算域的邊界條件。采用周期性邊界后，可以模擬波浪的穩(wěn)定流動(dòng)，出口邊界和入口邊界的流量可以保持守恒。采用周期性邊界的波浪數(shù)值水池模擬結(jié)果如圖2所示，可以看出施加周期性邊界之后，經(jīng)過(guò)一段時(shí)間的運(yùn)動(dòng)，波形和波浪壓力分布依然保持穩(wěn)定。

圖2 波浪的壓力分布Fig. 2 Pressure distribution of waves

為了研究計(jì)算域大小對(duì)于施加周期性邊界的數(shù)值水池的影響，本文對(duì)不同計(jì)算域的波浪水池進(jìn)行了模擬。標(biāo)準(zhǔn)波浪水池長(zhǎng)度為2L= 3.12 m，寬度為D= 0.6 m，然后建立4L×D和2L× 2D的波浪水池，得到不同計(jì)算域下波浪水池的速度場(chǎng)模擬結(jié)果，如圖3所示。從圖3中可以看出，不同計(jì)算域波浪水池的速度場(chǎng)在同一時(shí)刻保持一致，證明了計(jì)算域的大小對(duì)于波浪水池的影響很小。同時(shí)速度場(chǎng)也與線性波的理論速度場(chǎng)保持一致，證明了采用周期性邊界的波浪數(shù)值水池的有效性。為了提高計(jì)算效率，本研究中的波浪數(shù)值水池主要采取2L×D的計(jì)算域。

圖3 不同計(jì)算域下波浪的速度分布Fig. 3 Velocity of waves in different computational domains

1.4 四元數(shù)方法

流體-剛體耦合運(yùn)動(dòng)問(wèn)題中，涉及到兩種坐標(biāo)系，分別為以大地為基準(zhǔn)的慣性坐標(biāo)系（也稱固地坐標(biāo)系）和以運(yùn)動(dòng)剛體為基準(zhǔn)的非慣性坐標(biāo)系（也稱隨體坐標(biāo)系）。隨體坐標(biāo)系的坐標(biāo)原點(diǎn)Ob位于物體質(zhì)心處，如圖4所示。本文中流體流動(dòng)和剛體的平動(dòng)在固地坐標(biāo)系中計(jì)算，但是剛體的轉(zhuǎn)動(dòng)角速度和角加速度在隨體坐標(biāo)系中進(jìn)行計(jì)算。

圖4 坐標(biāo)系轉(zhuǎn)化示意圖Fig. 4 Schematic of coordinate system transformation

航行體在波浪中的運(yùn)動(dòng)軌跡非常復(fù)雜，涉及到六自由度的運(yùn)動(dòng)響應(yīng)。傳統(tǒng)的六自由度運(yùn)動(dòng)計(jì)算往往是基于歐拉角法，但歐拉角法在計(jì)算隨體坐標(biāo)系和固地坐標(biāo)系之間角速度的轉(zhuǎn)化時(shí)，角速度變換矩陣中存在著奇異值，這會(huì)導(dǎo)致歐拉角法在計(jì)算剛體翻轉(zhuǎn)過(guò)程中出現(xiàn)誤差。四元數(shù)法相對(duì)于傳統(tǒng)的歐拉角法，其角速度和四元數(shù)之間的變換矩陣中不存在奇異值，可以適應(yīng)不同姿態(tài)下剛體六自由度運(yùn)動(dòng)的計(jì)算，在剛體翻轉(zhuǎn)情況下的計(jì)算精度更高[34]。因此，本研究中采用四元數(shù)法來(lái)計(jì)算剛體的六自由度運(yùn)動(dòng)。

在數(shù)值水池中，剛體在固地坐標(biāo)系中的物體合力Fe和力矩Te可根據(jù)下列公式[12]得到：

其中，do為連體坐標(biāo)系中物體的質(zhì)心坐標(biāo)。初始四元數(shù)[q0q1q2q3]T是由初始?xì)W拉角Θe=[αβγ]T得到[17]，轉(zhuǎn)換關(guān)系如下：

四元數(shù)構(gòu)建坐標(biāo)系轉(zhuǎn)換矩陣R表達(dá)式為[17]：

隨體坐標(biāo)系下的力矩Tb和合力Fb可通過(guò)轉(zhuǎn)換矩陣R得到[35]，具體公式如下：

質(zhì)心加速度和轉(zhuǎn)動(dòng)角加速度公式[35]如下：

式中：M為物體質(zhì)量，Ib為相對(duì)質(zhì)心Ob的轉(zhuǎn)動(dòng)慣量矩陣。通過(guò)時(shí)間步積分可以得到隨體坐標(biāo)系下的質(zhì)心速度Ub和轉(zhuǎn)動(dòng)角速度?b，以及每個(gè)粒子在隨體坐標(biāo)系下的速度vb, j和加速度ab, j[35]：

固地坐標(biāo)系下的質(zhì)心速度Ue和轉(zhuǎn)動(dòng)角速度?e的計(jì)算[35]可根據(jù)以下公式：

Q為隨體坐標(biāo)系下角速度和固地坐標(biāo)系下四元數(shù)變化率之間的轉(zhuǎn)換矩陣[17]，表達(dá)式為：

從公式（13）中可以看出，隨體坐標(biāo)系中角速度到固地坐標(biāo)系四元數(shù)變化率的轉(zhuǎn)化矩陣中不存在奇異值，適合任何姿態(tài)的計(jì)算。在得到四元數(shù)變化率后，通過(guò)時(shí)間步積分可以更新四元數(shù)，進(jìn)而更新轉(zhuǎn)換矩陣R。固地坐標(biāo)系下的粒子位置和速度信息也可以實(shí)時(shí)更新[35]，具體公式如下：

四元數(shù)的完整計(jì)算流程如圖5所示。

圖5 四元數(shù)法流程Fig. 5 Flow chart of the Quaternion method

1.5 六自由度運(yùn)動(dòng)計(jì)算驗(yàn)證

為了驗(yàn)證基于四元數(shù)法的六自由度運(yùn)動(dòng)計(jì)算精度，本文同時(shí)采用自主編制的SPH程序和基于有限體積法（finite volume method, FVM）的成熟商業(yè)軟件STAR-CCM+，對(duì)方塊體六自由度落水過(guò)程進(jìn)行模擬，通過(guò)SPH與FVM結(jié)果對(duì)比，驗(yàn)證方塊體六自由度運(yùn)動(dòng)的數(shù)值計(jì)算精度。

方塊體落水示意圖如圖6所示。方塊體長(zhǎng)度為128 mm、寬度為68 mm、厚度為4 mm，質(zhì)量為0.2184 kg，質(zhì)心距方塊體底部65.2 mm、距離水面67.2 mm，方塊體連續(xù)繞y軸和z軸旋轉(zhuǎn)25°，轉(zhuǎn)動(dòng)慣量Ix、Iy、Iz分別為0.00039676 、0.00009、0.0003082 kg·m2。

圖6 方塊體落水示意圖Fig. 6 Schematic of a cuboid falling into water

圖7為采用SPH和FVM方法對(duì)方塊體落水過(guò)程的模擬結(jié)果，方塊體落水的速度變化曲線如圖8所示?？梢钥吹剑谒脑獢?shù)法的SPH方法模擬結(jié)果與FVM的模擬結(jié)果吻合得很好。

圖7 方塊體落水的SPH和FVM模擬結(jié)果Fig. 7 SPH and FVM simulation results of a cuboid falling into water

圖8 方塊體落水的的速度曲線Fig. 8 Velocity curve of a cuboid falling into water

2 航行體數(shù)值模型建立與驗(yàn)證

為了更好地與文獻(xiàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證，本文以楊曉光等開(kāi)展的高速航行體斜入水研究[19]為參考，選用相同的工況條件進(jìn)行計(jì)算。航行體質(zhì)量為1.35 kg，質(zhì)心距前端面為123 mm，轉(zhuǎn)動(dòng)慣量Ix、Iy、Iz分別為0.00059、0.005186、0.005186 kg·m2。航行體幾何示意圖如圖9所示。平靜水面下的航行體入水?dāng)?shù)值模型見(jiàn)圖10。為了降低計(jì)算規(guī)模，提高計(jì)算效率，同時(shí)保證剛體周圍流場(chǎng)的模擬更為精細(xì)，本文使用了多分辨率技術(shù)[36]，加密區(qū)域粒子隨剛體的運(yùn)動(dòng)動(dòng)態(tài)加密，如圖10（b）所示。

圖9 航行體幾何示意圖Fig. 9 Geometric diagram of the projectile

在進(jìn)行數(shù)值模擬之前，首先進(jìn)行粒子間距的收斂性分析。選擇粒子間距分別為10.0、5.0、2.5 mm，計(jì)算得到3種不同粒子間距下y方向航行體的加速度ay變化曲線，見(jiàn)圖11。當(dāng)粒子間距小于5 mm時(shí)，計(jì)算結(jié)果趨近于收斂。因此，為了提高計(jì)算效率，在后續(xù)算例中采用5 mm的粒子間距。

圖11 不同粒子間距下的y方向航行體加速度Fig. 11 The y-direction accelerations of the projectile with different particle sizes

隨后，通過(guò)在靜水狀態(tài)下的航行體高速入水模擬來(lái)驗(yàn)證本文SPH方法的可行性。文獻(xiàn)[19]的航行體高速入水試驗(yàn)在靜水狀態(tài)下開(kāi)展，航行體入水速度為60 m/s，入水傾角為20°。本文針對(duì)相同工況開(kāi)展了數(shù)值模擬研究，從軸向加速度aA對(duì)比結(jié)果圖12可見(jiàn)，SPH方法的模擬結(jié)果與試驗(yàn)結(jié)果吻合得很好，驗(yàn)證了本文SPH模型對(duì)航行體入水沖擊載荷的預(yù)報(bào)精度。

圖12 SPH方法航行體軸向加速度和試驗(yàn)結(jié)果[19]的對(duì)比Fig. 12 Comparison of the axial acceleration of the projectile between SPH simulation results and experimental results[19]

同時(shí)，在文獻(xiàn)[19]中楊曉光等也采用了網(wǎng)格類CFD方法對(duì)靜水工況下的航行體入水進(jìn)行了模擬。本文針對(duì)相同工況，開(kāi)展了無(wú)網(wǎng)格數(shù)值模擬研究。SPH方法與網(wǎng)格法在不同時(shí)刻的模擬結(jié)果對(duì)比如圖13所示，可以看出，SPH方法預(yù)報(bào)的彈體運(yùn)動(dòng)和自由面變形情況均與文獻(xiàn)結(jié)果吻合良好。此外，SPH結(jié)果中捕捉到了更為顯著的液滴飛濺，凸顯了無(wú)網(wǎng)格算法捕捉自由液面大變形的優(yōu)勢(shì)。

圖13 SPH方法航行體入水模擬結(jié)果和文獻(xiàn)結(jié)果[19]的對(duì)比Fig. 13 Comparison of water-entry simulation results of the projectile between SPH method and Ref.[19]

3 波浪環(huán)境下的航行體入水模擬

在波浪環(huán)境下航行體入水的SPH模擬中，波浪參數(shù)設(shè)置與楊曉光等[19]的數(shù)值模型一致。波浪周期T為1 s，波高H為0.15 m。在0.6 m的水深下，波長(zhǎng)L為1.56 m。數(shù)值水池計(jì)算域長(zhǎng)度為3.12 m，為兩倍波長(zhǎng)，寬度D為0.6 m。

航行體以與水平面20° 的傾角、60 m/s的速度入水，選擇4個(gè)波浪相位點(diǎn)入水展開(kāi)數(shù)值模擬，分別為0°、90°、180°和270°，入水點(diǎn)如圖14所示。航行體的質(zhì)心與對(duì)應(yīng)相位點(diǎn)的連線與水平面成20° 的傾角。

圖14 航行體入水位置的波浪相位角Fig. 14 Wave phase angles at the entry position of the projectile

3.1 不同相位角的模擬結(jié)果

航行體在波浪90° 和270° 相位角入水的模擬結(jié)果如圖15所示。

圖15 航行體入水過(guò)程的SPH模擬結(jié)果Fig. 15 SPH simulation results of the projectile in the waterentry process

從圖15中可以看出，航行體在入水后姿態(tài)角發(fā)生偏轉(zhuǎn)，頭部上揚(yáng)，在90° 的相位角入水工況中，入水15 ms后航行體開(kāi)始沖出水面。圖16給出了與楊曉光等[19]采用網(wǎng)格類方法結(jié)果的對(duì)比，可以看出本研究的SPH模擬結(jié)果與網(wǎng)格類方法的模擬結(jié)果非常吻合。其中，航行體在0° 相位角入水的模擬結(jié)果與靜水狀態(tài)下入水結(jié)果類似；而航行體在180° 相位角入水過(guò)程中，發(fā)生翻轉(zhuǎn)，頭部上揚(yáng)，航行體沖出水面。兩個(gè)相位角下航行體入水的SPH模擬結(jié)果與網(wǎng)格法模擬的航行體運(yùn)動(dòng)軌跡非常吻合。除此之外，得益于SPH粒子法的無(wú)網(wǎng)格和拉格朗日特性，可以自動(dòng)構(gòu)建自由液面，因此本研究能夠精確高效捕捉航行體入水時(shí)飛濺液滴的形態(tài)和軌跡，如圖15和圖16所示。而網(wǎng)格類CFD方法在處理液面大變形和液滴飛濺時(shí)，需要采用Level Set或者流體體積（volume of fluid,VOF）技術(shù)進(jìn)行自由液面追蹤，且追蹤精度受網(wǎng)格分辨率影響極大。圖16右側(cè)展示的網(wǎng)格法計(jì)算結(jié)果中，受網(wǎng)格尺寸限制，未能精確捕捉到液滴飛濺形態(tài)。

圖16 航行體入水過(guò)程的SPH模擬結(jié)果與參考結(jié)果[19]對(duì)比Fig. 16 The comparison between SPH simulation results and Ref.[19] of the projectile in the water-entry process

3.2 不同相位角對(duì)于航行體運(yùn)動(dòng)軌跡的影響

為了定量分析不同相位角下航行體的入水狀態(tài)和軌跡，提取了不同相位角下航行體入水后的轉(zhuǎn)動(dòng)角度和運(yùn)動(dòng)速度數(shù)據(jù)，變化曲線如圖17所示。

圖17 不同相位角下航行體入水速度和角度變化Fig. 17 The velocity and angle curves of the projectile with different phase angles

從圖17可知：4個(gè)相位角下、入水15 ms時(shí)，航行體傾斜角度全都從-20°變?yōu)檎?，都呈現(xiàn)出頭部上揚(yáng)的姿態(tài)；在270° 相位角下，航行體在入水后的翻轉(zhuǎn)速度最快，而在0° 相位角下航行體翻轉(zhuǎn)速度最慢，所以270° 相位角下航行體的x和y方向速度衰減最快，0° 相位角下航行體y方向速度衰減最慢；而在180°相位角下，航行體入水后半浮于水面上，所以在x方向上速度衰減最慢。

由于本研究中航行體入水的速度較大，波浪流動(dòng)速度的影響相對(duì)較小，相對(duì)于靜水條件下的航行體入水，波浪對(duì)于航行體入水軌跡的影響主要體現(xiàn)在不同相位角入水對(duì)于航行體實(shí)際入水角度的改變以及波浪形狀對(duì)于航行體出水時(shí)刻的影響。本文模擬的4種相位角入水工況中，根據(jù)圖15和圖16的模擬結(jié)果可知，0° 相位角下的航行體與水面的夾角最大，水動(dòng)力載荷主要集中在軸向，側(cè)向載荷相對(duì)較小，因此轉(zhuǎn)動(dòng)速度最慢，在y方向上的速度衰減最慢，彈道穩(wěn)定性方面表現(xiàn)最好；而180° 相位角下的航行體與水面的角度很小，航行體幾乎是懸浮在水面上，側(cè)向載荷較大，因此轉(zhuǎn)動(dòng)角度較快，航行體也迅速出水，由于航行體半浮于水面，受到的水阻力較小，所以在x方向上的衰減更慢；而在270° 相位角下，航行體的入水角度與波面幾乎平行，所以側(cè)向載荷較大，轉(zhuǎn)動(dòng)速度最快，但入水后航行體進(jìn)入到另一個(gè)波浪的波峰區(qū)域中，距離波面較遠(yuǎn)，出水時(shí)刻更為滯后，同時(shí)入水后的迅速轉(zhuǎn)動(dòng)也導(dǎo)致水阻力增大，所以在x和y方向的速度衰減最為明顯。靜水條件下的航行體入水軌跡發(fā)展主要受入水角度影響，而波浪條件下航行體入水軌跡受入水角度、波浪相位角和波浪形狀的多重影響，更為復(fù)雜。

4 結(jié)論

基于SPH方法，本文提出一種新型周期性波浪邊界技術(shù)并建立了波浪環(huán)境下的數(shù)值水池，簡(jiǎn)化了傳統(tǒng)推板造波過(guò)程中的繁瑣處理，可以迅速準(zhǔn)確地生成具有穩(wěn)定周期和波形的波浪。同時(shí)，將四元數(shù)法植入到SPH計(jì)算框架，消除了傳統(tǒng)歐拉角法中奇異值導(dǎo)致的誤差，實(shí)現(xiàn)了對(duì)物體入水后六自由度運(yùn)動(dòng)的精確模擬。

通過(guò)對(duì)航行體高速入水過(guò)程的模擬，SPH模擬結(jié)果與參考文獻(xiàn)中的實(shí)驗(yàn)結(jié)果及網(wǎng)格類CFD方法模擬結(jié)果吻合良好，且SPH方法能精確捕捉液滴飛濺的細(xì)節(jié)。對(duì)不同波浪相位角下航行體入水后的姿態(tài)變化的分析表明，波浪相位角會(huì)明顯影響航行體的入水運(yùn)動(dòng)軌跡，0° 相位角入水在保持彈道穩(wěn)定性方面表現(xiàn)最優(yōu)。

本研究中采用單相流模型和線性波理論，未考慮空氣效應(yīng)，也未對(duì)航行體高速入水過(guò)程中的空化現(xiàn)象進(jìn)行模擬，所以本研究模擬的是線性波條件下的不考慮氣穴效應(yīng)的航行體入水問(wèn)題。在未來(lái)研究中，擬開(kāi)發(fā)SPH空化模型，采用多相流SPH模型開(kāi)展計(jì)算，對(duì)入水彈道軌跡以及空泡演化進(jìn)行更深入的力學(xué)機(jī)理分析。