SPH 理論和方法在高速水動力學(xué)中的研究進(jìn)展

鐘詩蘊(yùn)，孫鵬楠*,2，呂鴻冠,2，彭玉祥，張阿漫

1 中山大學(xué) 海洋工程與技術(shù)學(xué)院，廣東 珠海 519082

2 南方海洋科學(xué)與工程廣東省實驗室，廣東 珠海 519082

3 哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001

0 引 言

面向海洋強(qiáng)國的戰(zhàn)略需要，船舶與海洋工程得到了大力發(fā)展，與之相關(guān)的結(jié)構(gòu)力學(xué)和水動力學(xué)等領(lǐng)域均產(chǎn)出了豐碩成果[1-3]。其中，高速艦船興波[4]、航行體（器）出水與入水[5]、水下爆炸與結(jié)構(gòu)毀傷[6-7]等典型高速水動力學(xué)問題與現(xiàn)代艦船及其武器裝備的綜合性能密切相關(guān)，但這些問題都具有大變形、動邊界、強(qiáng)對流、多介質(zhì)等復(fù)雜特征，而運用傳統(tǒng)理論、試驗或數(shù)值模擬方法較難得到徹底解決，因此，近年來成為了研究熱點和難點[8-9]。

高速水面航行器的興波、飛機(jī)水上滑行等過程涉及了復(fù)雜曲面結(jié)構(gòu)物與自由液面或水氣界面的復(fù)雜耦合作用，劇烈的液固耦合運動易產(chǎn)生非線性翻卷和破碎波浪[10]，甚至?xí)鹨旱螐?qiáng)烈飛濺，對航行器的結(jié)構(gòu)及內(nèi)部儀器零件的安全性構(gòu)成威脅，且誘導(dǎo)的噪聲[11]和尾跡[12-13]對航行器隱蔽性極為不利，精準(zhǔn)捕捉興波水花飛濺和預(yù)報砰擊載荷關(guān)系到航行器的隱身性、安全性設(shè)計。

在軍民用領(lǐng)域中，拋落式救生艇入水、飛行器海上回收，魚雷空投入水、潛射導(dǎo)彈出水等航行器或結(jié)構(gòu)物高速出、入水的問題[14-15]涉及瞬態(tài)強(qiáng)沖擊載荷、結(jié)構(gòu)物大幅運動、氣泡生長與潰滅等復(fù)雜過程[16]，是典型的多相流?固耦合問題[17]。這種高速跨介質(zhì)過程所產(chǎn)生的巨大砰擊載荷，一方面可能造成裝備結(jié)構(gòu)損壞、內(nèi)部零件失靈，另一方面也可能引起航行器彈道失穩(wěn)、投放失效等。因此，研究航行器的高速出、入水問題對于跨介質(zhì)武器裝備等的研發(fā)有著重要意義。

研究表明，水下爆炸沖擊波峰值及沖量相比空中爆炸有明顯提高[18]。艦船遭受近場水下爆炸時，沖擊波、高速射流和破片等首先對艦船造成局部毀傷[19]，隨后的氣泡脈動和射流沖擊可能造成艦船整體失效[7]。近場水下爆炸導(dǎo)致船體出現(xiàn)破口以后，劇烈涌流還可能會進(jìn)一步造成船體結(jié)構(gòu)損壞，最終導(dǎo)致船體沉沒。考慮到水下爆炸過程的復(fù)雜性和試驗研究的高成本，研發(fā)精準(zhǔn)的數(shù)值計算方法，預(yù)報水下爆炸載荷、結(jié)構(gòu)毀傷和艦船破損沉沒過程成為了各國海軍研究的重點[19]。

對于上述高速航行器的興波、航行體出水和入水、水下爆炸與結(jié)構(gòu)毀傷等工程問題，其背后的基礎(chǔ)力學(xué)問題離不開對復(fù)雜運動結(jié)構(gòu)物與多相流體之間的流?固耦合過程的求解。研究人員針對不同場景下的力學(xué)過程進(jìn)行簡化和凝練，提出了多種理論模型[20-21]，解釋了氣泡動力學(xué)或簡單構(gòu)型結(jié)構(gòu)物與流體作用過程背后的力學(xué)機(jī)理。然而，理論模型對于復(fù)雜的工程實際問題適用性有限。而模型試驗雖然是獲得可靠結(jié)果的重要手段[4,22]，但開展試驗通常需要較大的場地和昂貴的設(shè)備，尤其是在高速水動力問題的模型試驗中，尺度效應(yīng)明顯[23]，相關(guān)物理參數(shù)測量困難，試驗數(shù)據(jù)的可重復(fù)性較差，試錯成本也高。

近年來，隨著計算機(jī)技術(shù)的發(fā)展，涌現(xiàn)出了許多數(shù)值計算方法。目前，基于不同的離散和求解形式，數(shù)值模擬方法可分為有網(wǎng)格法和無網(wǎng)格法[24]。前者在解決大變形流?固耦合問題時容易受到網(wǎng)格畸變的困擾；后者得益于天然的拉格朗日特性及逐漸完備的粒子近似理論優(yōu)勢，在模擬大變形問題時不受邊界變形量的限制[25]，故在高速水動力問題中得到了廣泛應(yīng)用[24]。典型的無網(wǎng)格法包括離散元法（discrete element method, DEM）[26]、光滑粒子流體動力學(xué)（smoothed particle hydrodynamics, SPH）方法[27-29]、移動粒子半隱式（moving particle semi-implicit）方法[30]等。MPS 方法和SPH 方法在求解水動力問題方面難分伯仲[31-34]，但MPS 方法基于完全不可壓縮假設(shè)[35-36]，難以應(yīng)用于壓縮效應(yīng)明顯的問題，例如水下爆炸的氣泡脈動過程等；而DEM 是離散元算法[37]，主要用于處理離散粒子（例如沙礫、冰塊等顆粒流）大變形問題[38-40]，并不適用于高速水力學(xué)研究。相較而言，SPH 理論和方法在處理流體大變形和可壓縮問題時的優(yōu)勢明顯，便于開展并行計算，工程適用性強(qiáng)，更適合于求解大型高速水動力問題。

海洋結(jié)構(gòu)物在高速航行或跨介質(zhì)運動時，對自由液面以下的單相流動區(qū)域，準(zhǔn)確模擬湍流邊界層是關(guān)鍵[41]。針對這方面的應(yīng)用，有網(wǎng)格法因其湍流模型較為成熟、便于采用多分辨率網(wǎng)格技術(shù)而受到研究人員的青睞[42]。但對于自由液面或多相界面附近區(qū)域，因存在流體高速砰擊引起的液滴飛濺和水氣摻混等過程，有網(wǎng)格法的精度容易受到網(wǎng)格尺度的極大限制[24]。例如，一旦液滴飛濺到網(wǎng)格稀疏區(qū)域，易引起流體質(zhì)量不守恒造成的液滴飛濺消失。相比而言，SPH 無網(wǎng)格粒子法用于模擬液滴飛濺問題更具優(yōu)勢[4,22]。該方法基于核近似、粒子近似和虛粒子邊界等基本原理，輔之以日益完善的高精度修正理論和技術(shù)，在模擬運動邊界和大變形界面等問題時具有較強(qiáng)計算能力和廣闊發(fā)展前景。

盡管SPH 方法在計算的全過程中能夠精確滿足質(zhì)量守恒要求并擅長捕捉自由液面翻卷、破碎和噴濺等復(fù)雜現(xiàn)象，但運用傳統(tǒng)SPH 方法計算的壓力載荷一般伴隨嚴(yán)重的數(shù)值噪聲[43]，涉及復(fù)雜區(qū)域時還容易引起張力不穩(wěn)定性和數(shù)值空洞等[44]。為解決這些問題，經(jīng)過近半個世紀(jì)的不斷發(fā)展，SPH 方法在理論和數(shù)值技術(shù)方面得到了不斷優(yōu)化，計算精度和穩(wěn)定性獲得巨大提升，在新一代GPU 并行技術(shù)的加持下，顯著提高了計算效率[45]，具備了求解工程問題的能力。并且，隨著嵌入SPH的多算法耦合技術(shù)的出現(xiàn)和發(fā)展[6,46-47]，SPH 理論和方法的優(yōu)勢得以充分發(fā)揮，使其為復(fù)雜多元的高速水動力學(xué)問題提供了新的高精度、高效率求解方案。

綜上所述，為進(jìn)一步闡述SPH 理論和方法在船海高速水動力學(xué)相關(guān)問題中的研究進(jìn)展，本文將從SPH 的基本原理及數(shù)值技術(shù)研究進(jìn)展出發(fā)，圍繞高速水面航行器的水動力以及高速跨介質(zhì)航行器水動力問題、水下爆炸與結(jié)構(gòu)毀傷等方面展開綜述，并對SPH 理論和方法在求解高速水動力問題中的未來發(fā)展趨勢進(jìn)行展望，以期為涉及上述問題的相關(guān)研究提供參考。

1 SPH 理論及技術(shù)進(jìn)展

1.1 SPH 基本理論與算法

1.1.1 SPH 基本理論

SPH 方法最初是為研究天體動力學(xué)而提出的[48]，上世紀(jì)90 年代開始應(yīng)用于水動力問題的模擬[49]。以物體入水問題為例，圖1 給出了基于流體體積（volume of fluid，VOF）的歐拉網(wǎng)格法[50]及拉格朗日描述的SPH 方法[51-52]模擬圓柱體及楔形體入水以及液面飛濺的過程。在歐拉網(wǎng)格法的計算結(jié)果中，可見噴濺射流尖端存在非物理性流體消失，水氣界面較模糊的情況，而SPH 方法中粒子可以精確地捕捉到流體飛濺特征。上述結(jié)果體現(xiàn)了SPH 方法在處理大變形問題方面的優(yōu)勢，因此被廣泛應(yīng)用于船海工程領(lǐng)域[5,53]。

從原理上看，SPH 方法基于核函數(shù)近似理論計算空間導(dǎo)數(shù)，可對特定物理場的梯度、散度等采用積分形式進(jìn)行計算[54]；在數(shù)值實現(xiàn)方面，SPH方法一般將流體域離散成緊湊的拉格朗日粒子，每個粒子均具有相應(yīng)的質(zhì)量和體積，通過粒子求和代替核近似中的積分，也稱“粒子近似”。如圖2所示，基于核近似理論和粒子近似理論，SPH 方法可對控制方程進(jìn)行離散求解。圖中：f和f分別代表標(biāo)量和矢量； ρ，u，p，r，V分別表示密度、速度、壓力、位置矢量和體積，下標(biāo)i和j表示粒子編號；g，T， μ分別表示重力加速度、應(yīng)力張量和動力黏性系數(shù)；W為核函數(shù)；h為光滑長度；c0， ρ0分別為流體聲速及其參考密度。為獲得更高的壓力場、速度場計算精度，可在離散控制方程中進(jìn)一步添加耗散項[24]。隨后，在特定的初始條件和邊界條件（自由液面邊界[55]、壁面邊界[56]、流入流出邊界[57]、無反射邊界[58]等）下，將離散控制方程在時域內(nèi)進(jìn)行積分，可獲得SPH 控制方程的數(shù)值解。

圖2 SPH 方法求解流程[24]Fig. 2 The process of SPH simulation[24]

SPH 粒子攜帶各自的速度、壓力、密度等物理量，按照流體質(zhì)點本身的速度進(jìn)行拉格朗日運動。這些粒子的運動代表了流體的流動，粒子所到之處也代表了所處位置的流場速度、壓力、密度等物理量的變化。鑒于SPH 理論和方法的拉格朗日特性及無網(wǎng)格粒子特性，在模擬自由液面流動問題時，自由面的運動學(xué)條件不僅能夠自動滿足[55]，而且也能夠自動追蹤水氣和流?固界面，所以十分適合處理船海工程中的大變形、動邊界、多介質(zhì)和強(qiáng)對流的流體動力學(xué)問題[8]。

1.1.2 SPH 算法的分類

經(jīng)過近半個世紀(jì)的發(fā)展，基于不同的理論框架，SPH 方法已發(fā)展出了多種形式，主要包含兩種類型：弱可壓SPH （weakly compressible SPH,WCSPH）方法[43]和不可壓SPH （incompressible SPH, ISPH）方法[59]，并已在各自領(lǐng)域得到廣泛的應(yīng)用。對于水介質(zhì)，一般認(rèn)為其近似于不可壓縮，但為便于求解壓力，Monaghan 等[49]最早引入了狀態(tài)方程，如圖2 所示控制方程中的第4 個等式。基于狀態(tài)方程求解壓力時，積分的時間步長受到克朗條件（Courant–Friedrichs–Lewy condition），也即流體聲速的極大限制。此外，Monaghan 等[49]還提出了弱可壓假設(shè)，以馬赫數(shù)小于0.1 為標(biāo)準(zhǔn)確定人工聲速，確保流體密度變化不大于1%，以此維持流體的近似不可壓特性，所以該方法也被稱為WCSPH 方法。相反，在ISPH 方法中，流場的密度保持不變，亦即速度散度處處為0。Cummins 和Rudman[59]在SPH 方法框架內(nèi)，建立了每個時間步的泊松壓力方程，用于求解壓力，并提出了ISPH 方法來模擬不可壓縮流動問題。但是，ISPH 方法涉及大型稀疏矩陣運算和自由液面粒子搜索，計算量相對較大[60-61]，在大規(guī)模粒子數(shù)量的三維計算中挑戰(zhàn)較大。WCSPH 相關(guān)理論和方法經(jīng)過近20 年的發(fā)展，計算精度已得到明顯提高，相比于ISPH 方法，WCSPH 方法涉及的GPU 存儲量更小，計算效率更高，所以更適合于開源程序的設(shè)計。例如DualSPHysics[62]、SPHinXsys[63]等開源程序包，它們均是基于WCSPH 理論框架，目前已被廣泛應(yīng)用于模擬三維復(fù)雜自由液面流動等問題。

1.2 SPH 數(shù)值技術(shù)的進(jìn)展

1.2.1 SPH 壓力計算精度

在船海流體力學(xué)中，準(zhǔn)確計算壓力對流體砰擊載荷預(yù)報至關(guān)重要。傳統(tǒng)WCSPH 模型計算壓力場時，一般會出現(xiàn)難以避免的壓力場噪聲，其主要原因是：

1） SPH 核函數(shù)的近似誤差和數(shù)值離散誤差會導(dǎo)致密度場出現(xiàn)波動，微小的密度波動被狀態(tài)方程放大，變成劇烈的壓力場噪聲；

2） 因流體的弱可壓假設(shè)，聲速一般會被降低。在出現(xiàn)流體砰擊時，砰擊位置輻射的壓力波以較小的聲速在流場中傳播，遇到邊界時會發(fā)生多次反射和疊加，同樣也會在流場中表現(xiàn)為壓力場噪聲。

針對壓力場噪聲問題，Colagrossi 等[43]率先提出SPH 密度過濾技術(shù)，成功緩解了SPH 壓力場脈動。Vila[64]通過Riemann 求解器來計算粒子間的相互作用，提出了Riemann SPH 方法，此方法能夠?qū)PH 數(shù)值誤差和弱可壓假設(shè)造成的壓力脈動進(jìn)行快速耗散，從而提高壓力場的光滑性。Marrone 等[56]提出δ-SPH 方法并應(yīng)用于自由液面流動的模擬，該方法引入了密度耗散項Dρ，該耗散項起到與Riemann 求解器相似的耗散作用，也同樣能夠提高壓力場的計算精度。圖3 所示的δ-SPH控制方程中，Tv表示黏性力張量。相比Riemann SPH 方法，δ-SPH 方法的優(yōu)勢在于數(shù)值耗散量可通過連續(xù)方程中的密度耗散項[65]進(jìn)行定量控制，從而實時監(jiān)測流體能量的耗散量，這對于SPH 數(shù)值波浪水池的建立十分有益。

1.2.2 SPH 粒子均勻化技術(shù)

傳統(tǒng)SPH 理論基礎(chǔ)和算法精度很大程度上依賴于粒子分布的均勻性[66]。為此，Nestor 等[67]最早提出了粒子位移修正方法（particle shifting techniques，PST），而Lind 等[68]將PST 方法引入到自由液面計算中，初步解決了ISPH 模擬中的粒子不均勻及流場撕裂的問題。如圖3 中第3 個框圖所示，Sun 等[69-70]將PST 技術(shù)與δ-SPH 進(jìn)行結(jié)合，提出了δ+-SPH 方法（其中， δu為修正速度項，m表示粒子質(zhì)量），隨后又發(fā)展出了張力不穩(wěn)定性控制（tensile instability control, TIC）技術(shù)[44]，克服了高負(fù)壓區(qū)域的流場數(shù)值空洞問題，并成功應(yīng)用到中、高雷諾數(shù)的黏性邊界層流動的模擬中。其后，又進(jìn)一步建立了滿足一致性的δ+-SPH 方法，將PST 技術(shù)中位移修正以修正速度的方式表示，并在控制方程予以考慮，從而解決了早前PST 技術(shù)造成的自由液面流動中流體體積的非物理性膨脹問題，極大提高了WCSPH 理論和方法的精度及穩(wěn)定性[70]。最近，Antuono 和Sun 等[66]又建立了任意拉格朗日?歐拉(arbitrary Lagrangian-Eulerian,ALE)形式的δ-ALE-SPH 理論與方法，將質(zhì)量方程和質(zhì)量數(shù)值耗散項Dm引入到控制方程中，提高了ALE 算法的穩(wěn)定性，進(jìn)一步完善了含粒子位移修正的該SPH 理論框架。針對海洋工程應(yīng)用問題，Lü和Sun 等[5,17]針對液艙的晃蕩、物體出水和入水、黏性繞流等問題進(jìn)行了深入分析，從理論角度指出了此類問題中張力不穩(wěn)定性產(chǎn)生的原因，有針對性地提出了完備的控制技術(shù)，以克服張力不穩(wěn)定性問題。

圖3 SPH 理論和方法的演化Fig. 3 The evolution of SPH theory and method

1.2.3 高精度SPH 理論和離散格式

制約SPH 方法計算精度的另一個重要原因是，傳統(tǒng)SPH 理論和方法離散格式本身的精度不夠，誤差主要來源于核函數(shù)近似誤差和粒子近似誤差，其中后者又被稱為數(shù)值離散誤差。因此，為提高SPH 方法的計算精度，研究人員先后提出了修正SPH 算法[71]、移動最小二乘算法[72]、有限粒子法[73]以及近兩年發(fā)展的高階SPH 方法，例如WENO-SPH 方法[74]、TENO-SPH 方法[75]等。高階SPH 理論是在Riemann SPH 的基礎(chǔ)上發(fā)展而來，其通過引入網(wǎng)格法的高階計算框架來提高求解精度。如圖4 所示，在Riemann SPH 控制方程中，e表示比內(nèi)能[75]，p?和u?分別表示粒子對之間的壓力及速度[76]，是Riemann 方程求解的中間變量，也是重構(gòu)框架的作用對象。Meng 和Zhang 等[75]提出的TENO-SPH 理論和方法，有效提高了SPH 方法的計算精度和魯棒性，驗證了高階SPH 理論和應(yīng)用的可行性。此外，為減少離散誤差，可通過增加光滑長度與粒子間距的比例系數(shù)，來增加核函數(shù)半徑范圍內(nèi)的粒子數(shù)量，達(dá)到減小離散誤差的目的，使之滿足理論要求，提高SPH 方法的計算精度，具體詳見文獻(xiàn)[77]。盡管高精度SPH 方法在解決特定問題方面明顯提升了計算精度，但在物理量守恒性及算法穩(wěn)定性上仍有很大的改善空間。對于工程化應(yīng)用，因涉及了巨大的粒子數(shù)量和較長模擬時間，目前仍普遍采用傳統(tǒng)的SPH數(shù)值離散理論，輔之以δ-SPH 或Riemann SPH 等數(shù)值方法，能夠嚴(yán)格保障質(zhì)量、動量的守恒特性，實現(xiàn)對自由液面流動的長時間的穩(wěn)定數(shù)值模擬。

圖4 高階SPH 理論和方法示意圖Fig. 4 Schematics of the high-order SPH theory and method

1.2.4 SPH 方法的工程應(yīng)用

隨著SPH 理論和方法的不斷發(fā)展，其已被廣泛應(yīng)用于海洋工程領(lǐng)域。事實上，SPH 工程化應(yīng)用的重要前提是其具備了復(fù)雜結(jié)構(gòu)物的粒子化建模和離散能力，以及復(fù)雜形狀邊界條件的處理能力。經(jīng)過數(shù)十年的發(fā)展，慕尼黑工業(yè)大學(xué)Zhu 等[78]開發(fā)出了基于.STL 格式文件的粒子化填充算法，并已成功用于復(fù)雜構(gòu)型結(jié)構(gòu)物的粒子化離散。圖5 所示為飛機(jī)模型從CAD 到粒子化的離散過程，這為后續(xù)的工程化應(yīng)用研究奠定了良好的基礎(chǔ)。基于離散后的粒子，將其設(shè)置為虛擬粒子（ghost particle），采用基于固定虛粒子（fixed ghost particle）的固壁邊界理論和技術(shù)[56,79]，即可實現(xiàn)固定、運動甚至變形的壁面邊界的高精度施加處理[51, 80]。

圖5 復(fù)雜構(gòu)型結(jié)構(gòu)物的CAD 建模和粒子化離散Fig. 5 CAD modeling and particle discretization of a complex structure

1.2.5 SPH 硬件加速技術(shù)

在SPH 數(shù)值計算中，精確反映三維復(fù)雜結(jié)構(gòu)物構(gòu)型通常需要設(shè)置精細(xì)的粒子分辨率，故也必然導(dǎo)致高昂的計算成本。隨著計算機(jī)技術(shù)的快速發(fā)展，硬件并行計算技術(shù)的出現(xiàn)極大地提高了SPH 方法的計算效率，推動了SPH 方法的工程化應(yīng)用[45,60]。硬件加速技術(shù)一般可基于CPU 或GPU兩種并行計算框架。絕大部分網(wǎng)格法（例如有限體積法（FVM）等）均采用CPU 并行計算，原因是該方法需要求解大型稀疏矩陣和隱式時間積分，并不適合或較難采用GPU 并行計算。而弱可壓SPH 方法作為一種純拉格朗日粒子方法，其計算原理基于點對點的方式，且時間積分為顯式格式，這兩個特性決定了SPH 方法非常適合用于GPU并行計算。目前，大量文獻(xiàn)已表明，在同等模擬規(guī)模下SPH 方法的GPU 并行可比CPU 并行在效率上提高1～2 個數(shù)量級[45]，極大提升了SPH 方法的工程實用性，尤其是基于該方法開發(fā)的計算軟件能夠在裝有GPU 的個人電腦上運行[62]，并不完全依賴于大型超級計算機(jī)或服務(wù)器，因此在未來工程領(lǐng)域有較大的應(yīng)用潛力。

當(dāng)前，實現(xiàn)SPH 并行計算首先可基于并行計算框架，自主編制SPH 程序，例如Marrone 等[22]搭建的MPI-OpenMP 混合編程模型，設(shè)計了一種新型的三維并行SPH 求解器，成功實現(xiàn)了對三維船體興波的全局仿真，為自主構(gòu)造SPH 并行計算模型提供了參考。此外，還可以采用目前主流的并行SPH 開源程序包完成計算，例如SPHinXsys[63]、DualSPHysics[62]等開源程序包，它們均在三維復(fù)雜問題應(yīng)用方面有出色表現(xiàn)[81]，為高精度SPH 理論和方法在船海工程中的實際應(yīng)用提供了強(qiáng)大的計算能力保障。

2 水面航行器高速水動力問題

隨著上述數(shù)值技術(shù)的發(fā)展和應(yīng)用，運用SPH理論和方法求解實際工程問題的能力得到了不斷加強(qiáng)。SPH 方法利用其在處理自由液面大變形和噴濺問題方面的優(yōu)勢，可以十分方便地對高速水動力問題開展相關(guān)研究。例如，圖6 展示的是三維SPH 方法模擬的“中山大學(xué)”號（SYSU）科考船高速航行興波和飛機(jī)水上迫降過程。由圖可見，SPH 方法準(zhǔn)確再現(xiàn)了自由液面的翻卷、破碎、噴濺過程。

圖6 三維SPH 方法模擬的“中山大學(xué)”號科考船高速航行興波和飛機(jī)水上迫降Fig. 6 3D SPH simulations of the breaking wave of SYSU research vessel in high-speed and the aircraft ditching process

以下將圍繞航行器的興波與非線性波浪砰擊等問題，對SPH 理論和方法在水面航行器高速水動力問題中的應(yīng)用研究進(jìn)展進(jìn)行綜述。

2.1 高速水面航行器興波

目前，計算流?固耦合問題通常使用有限元法（FEM）、FVM 等工業(yè)化程度高、應(yīng)用發(fā)展較為成熟的求解器，例如，常使用VOF[82]或Level-Set[83]等方法模擬自由面流動。在單相、小變形、低弗勞德數(shù)的情況下，使用以上的網(wǎng)格法模擬航行器的興波和尾流可以得到精準(zhǔn)結(jié)果，但對高速航行過程中產(chǎn)生的強(qiáng)非線性自由液面流動、流體噴濺等情況，網(wǎng)格法較難精準(zhǔn)捕捉到波浪破碎、浪花飛濺等細(xì)節(jié)。圖7 展示了FVM 和SPH 方法模擬高速艦船興波的結(jié)果與Tagliafierro 等[84]的模型試驗結(jié)果的對比。由圖可見，F(xiàn)VM 只能夠模擬出微幅的液面變形，而SPH 方法能夠很好地捕捉液面的飛濺細(xì)節(jié)。可見，SPH 充分發(fā)揮了拉格朗日粒子方法的理論優(yōu)勢，能更準(zhǔn)確模擬興波特征。

圖7 高速滑行艇興波各方法模擬結(jié)果對比Fig. 7 Comparison of wave-making of high-speed planning boat

早期的計算機(jī)內(nèi)存和效率低，三維仿真計算成本高，故一般采用切片法、二維半[85-86]理論研究航行器的興波問題。Marrone 等[4]通過2D+tSPH方法對細(xì)長船體興波的問題進(jìn)行研究，對不同弗勞德數(shù)下的船艏興波演化和尾跡進(jìn)行模擬，較成功地模擬了細(xì)長船舶兩側(cè)興波特征，并通過試驗對比驗證了結(jié)果的準(zhǔn)確度。但是，2D+tSPH 方法的研究對象局限于縱向速度變化小的細(xì)長船體，對于橫剖面變化大的復(fù)雜構(gòu)型，只有進(jìn)行三維仿真才可以盡可能全面地保留興波特征。最近，Mintu等[87]成功實現(xiàn)了實尺度高速艦船船艏破波的三維模擬，驗證了三維SPH 計算方法的數(shù)值精度和技術(shù)可行性。

與真實試驗結(jié)果[84]相比，SPH 方法已可以初步模擬船體兩側(cè)噴射的水花，如圖7 所示，但目前仍未能對微小尺度的水氣混合多相流實現(xiàn)精準(zhǔn)捕捉。高密度比混合多相流模擬是未來SPH 方法的發(fā)展方向，需選用更精細(xì)的粒子分辨率、開發(fā)新的理論和數(shù)值模型進(jìn)行求解。Monaghan 和Kocharyan[88]首次提出了基于體積分?jǐn)?shù)理論的SPH混合物模型，經(jīng)過十余年的發(fā)展，F(xiàn)onty等[89-90]提出的SPH 混合物模型在模擬振蕩水柱波浪能轉(zhuǎn)換器的水?氣高速耦合作用方面表現(xiàn)優(yōu)越，可進(jìn)一步拓展應(yīng)用于航行器（體）的興波模擬中。SPH混合物模型在高速水動力學(xué)應(yīng)用領(lǐng)域發(fā)展前景廣闊，相信未來會應(yīng)用于模擬興波水花及其他大密度比多相流體擴(kuò)散、混合過程。

2.2 非線性波浪砰擊

2.2.1 SPH 數(shù)值造波與消波

水面航行器高速航行時會受到反復(fù)交替的非線性波浪砰擊載荷作用，尤其是在惡劣海況下，波浪力所引起的砰擊壓力和力矩可能會造成海洋結(jié)構(gòu)物的局部毀傷甚至整體傾覆，因此研究波浪與浮體耦合作用意義重大。近年來，數(shù)值波浪水池以及波浪與結(jié)構(gòu)物的耦合作用模型迅速發(fā)展，其中數(shù)值造波與消波技術(shù)成為研究重點。得益于SPH 方法模擬動邊界的天然優(yōu)勢，SPH 方法可通過固定虛粒子理論和技術(shù)構(gòu)建數(shù)值造波板，對其運動予以控制，模擬物理水池中的推板式[91]、活塞式[92]、搖板式[93]等造波機(jī)，生成多種形式的波浪[94]。

SPH 數(shù)值水池消波技術(shù)可分為被動消波和主動消波。被動消波是指通過斜坡和多孔結(jié)構(gòu)[95]或在數(shù)值水池中設(shè)置人工阻尼區(qū)進(jìn)行消波[92,96]，而主動消波方法是通過改變數(shù)值造波板的運動機(jī)制達(dá)到消波目的。虛粒子組成的造波板會根據(jù)反射波的變化重新調(diào)整運動方程，使反射波與新生成的波疊加形成與原波形一致的波浪[91]，從而保證結(jié)構(gòu)物受恒定的數(shù)值波浪作用。研究表明，SPH 方法可以產(chǎn)生包括不規(guī)則波[91]、畸形波[94]、孤立波[92]等各類波形，造波質(zhì)量較高，為進(jìn)一步研究浪固耦合問題提供了技術(shù)支撐和保障。

2.2.2 甲板上浪與畸形波砰擊

若巨幅波浪砰擊船體、海洋平臺等結(jié)構(gòu)物，一般會造成嚴(yán)重的波浪爬升，甚至可能翻越甲板，發(fā)生甲板上浪現(xiàn)象，導(dǎo)致甲板設(shè)備嚴(yán)重毀壞或人員傷亡等，因此，受到了研究人員的高度重視[97-98]。早期的SPH 方法在研究艦船甲板上浪問題時，為降低計算難度，通常是通過二維仿真波浪沖擊翻越固定水平平板[10]過程來進(jìn)行研究。而隨著計算機(jī)計算能力的提高，現(xiàn)已能夠?qū)ε灤叫械纳侠饲闆r進(jìn)行三維模擬，進(jìn)一步研究甲板上浪對上層建筑的破壞力[99]。為更精確預(yù)報砰擊載荷，Areu-Rangel 等[100]通過對比SPH 方法的仿真結(jié)果和精密實驗儀器的測量結(jié)果，優(yōu)化了甲板上浪的SPH 數(shù)值模型，該方法用于預(yù)報甲板上浪壓力載荷的精度得以驗證，有效性和可行性也得到證明。

傳統(tǒng)的規(guī)則波和不規(guī)則波理論難以解釋深海中時常出現(xiàn)的極大波高，因此，近年來興起了對極端海況下畸形波和畸形波砰擊海洋結(jié)構(gòu)物的研究。Sun 等[94]基于規(guī)則波和畸形波理論，采用推板和搖板式造波數(shù)值技術(shù)聯(lián)合人工阻尼消波技術(shù)，建立了數(shù)值波浪水池的SPH 數(shù)值模型，以研究規(guī)則波、聚焦波與結(jié)構(gòu)物的浪固耦合作用，結(jié)果如圖8（a）所示。針對Zhao 等[101-102]提出的畸形波?結(jié)構(gòu)物耦合運動算例，SPH 方法可預(yù)報浮體結(jié)構(gòu)物的縱蕩、升沉、縱搖運動和波浪砰擊、甲板上浪、回流等現(xiàn)象，計算結(jié)果與試驗結(jié)果吻合，如圖8（b）所示。

圖8 SPH 數(shù)值波浪水槽模擬的非線性波浪砰擊固定或浮動結(jié)構(gòu)物的結(jié)果對比Fig. 8 The results comparison of nonlinear wave slamming on a fixed or floating structure simulated in a SPH numerical wave flume

畸形波砰擊過程中會卷入氣體，波浪翻卷上涌甲板后，受氣相的影響出現(xiàn)負(fù)壓區(qū)，如圖8（a）所示。Sun 等[94]通過研究砰擊過程中的負(fù)壓現(xiàn)象，指出一次砰擊過程相當(dāng)于進(jìn)水與出水過程的耦合，非線性波浪砰擊退去后，會在結(jié)構(gòu)下方（例如外漂底部）產(chǎn)生負(fù)壓，而負(fù)壓在SPH 方法中會帶來張力不穩(wěn)定性，隨之造成非物理性的流場空洞并導(dǎo)致流場演化的錯誤。文獻(xiàn)[94] 模擬結(jié)果還表明，當(dāng)波浪砰擊甲板邊緣時，由于砰擊速度高，卷入的氣穴會形成負(fù)壓（相對于背景壓力）而導(dǎo)致流體吸附。此時，流體并不會立即在甲板上分離，而是在甲板表面形成高速吸附流，這就要求SPH 模擬能精確地預(yù)報流場中負(fù)壓形成的吸附力。對于此問題的研究，研究人員建議使用SPH 張力不穩(wěn)定性控制技術(shù)[44]。此外，波浪翻越甲板的過程還需考慮空氣的氣墊效應(yīng)[94]，否則容易導(dǎo)致壓力峰值預(yù)報失真。總之，涉及甲板上浪的浪固耦合問題還需進(jìn)一步完善SPH 的理論和方法，建立多相流SPH 模型[17]，在考慮空氣影響的情況下進(jìn)行數(shù)值模擬研究。

3 航行體跨介質(zhì)水動力問題

航行體跨介質(zhì)過程是指物體出水和入水，實質(zhì)上也是一個復(fù)雜的多相流?固耦合過程。流體與結(jié)構(gòu)物的相互強(qiáng)耦合可能導(dǎo)致相變、湍流和空化，力學(xué)機(jī)理極具復(fù)雜性和挑戰(zhàn)性。目前，對出水和入水力學(xué)機(jī)制的研究方法主要采用試驗方法[103]、數(shù)值模擬[51]或兩者結(jié)合的方法。研究表明，SPH 方法能很好地模擬此類大變形、多邊界和多相流問題[15]。

3.1 細(xì)長體入水

細(xì)長航行體入水是典型的瞬態(tài)強(qiáng)非線性動力學(xué)過程，入水瞬間的砰擊作用會造成強(qiáng)烈的液面變形，形成液滴噴濺。面對此類問題，網(wǎng)格法受網(wǎng)格尺度限制，對噴濺的模擬略顯遜色。圖9 給出了采用SPH 方法、VOF 方法計算模擬細(xì)長圓柱體（細(xì)長體）以100 m/s 高速入水的結(jié)果對比。由圖9可見，網(wǎng)格法（VOF）難以捕捉到噴濺細(xì)節(jié)，多相界面相對模糊，而SPH 方法在流場噴濺特性的捕捉方面表現(xiàn)更佳。圖10 所示為該算例中運用不同算法預(yù)報的細(xì)長體端部砰擊壓力和運動速度曲線。由圖10 可見，SPH 方法能夠精確地捕捉到瞬態(tài)沖擊載荷的峰值，獲得的細(xì)長體運動速度與CEL 方法更吻合。可見，SPH 方法能夠精確地對航行體高速入水沖擊進(jìn)行仿真，尤其是液面噴濺模擬方面的優(yōu)勢更明顯。

圖9 不同方法模擬的細(xì)長體高速入水結(jié)果對比 [50]Fig. 9 Comparison of high-speed water entry of a slender body simulated by SPH and VOF [50]

圖10 細(xì)長體高速入水砰擊壓力與運動速度變化曲線Fig. 10 Time evolution of impact pressure and velocity of the slender body during high-speed water entry

高速入水問題主要關(guān)注細(xì)長體結(jié)構(gòu)局部砰擊區(qū)域所受載荷及其附近流場的演化過程，而遠(yuǎn)場流動變化隨著距離的增大而減小。因此，為了在保證局部高精度的同時提高計算效率，Sun 等[44]采用自適應(yīng)粒子細(xì)化技術(shù)（adaptive particle refinement, APR）對大范圍的計算域分設(shè)多級粒子分辨率。該技術(shù)也能根據(jù)需求定義高分辨率區(qū)域，通過此區(qū)域動態(tài)跟蹤細(xì)長體結(jié)構(gòu)運動的過程，并可對整個彈體周圍流場區(qū)域進(jìn)行加密，精準(zhǔn)預(yù)報流體對彈體和彈道的影響，如圖11（a）所示。此外，APR 技術(shù)也可僅對彈體頭部進(jìn)行局部加密，捕捉彈頭區(qū)域的流動特征，如圖11（b）所示（ ?x表示粒子間距）。總之，APR 技術(shù)為SPH 方法實現(xiàn)高效精準(zhǔn)的仿真提供了技術(shù)支撐。

圖11 多級粒子分辨率SPH 技術(shù)模擬的細(xì)長體高速入水過程Fig. 11 High-speed water entry process of a slender body simulated by SPH with APR technique

前人基于SPH 方法的許多入水問題的研究，例如旋轉(zhuǎn)[104]、傾斜[105]入水的仿真以及空腔演化特征等，大部分忽略了空腔閉合后的內(nèi)部氣體對結(jié)構(gòu)運動的影響。首先，這是因為空腔閉合后會產(chǎn)生高速射流[17]，沖擊結(jié)構(gòu)物表面并造成運動姿態(tài)變化，所以空腔的產(chǎn)生及氣相的影響對結(jié)構(gòu)載荷與運動預(yù)報有著重要意義[103,106]，需要采用多相流SPH 模型開展計算，詳見文獻(xiàn)[17]。其次，砰擊入水的高速運動會導(dǎo)致細(xì)長體周邊產(chǎn)生空化，且目前SPH 理論框架內(nèi)尚未建立完備的空化理論和數(shù)值模型，如何做到空化、湍流[107]、多相SPH 模型的結(jié)合以實現(xiàn)精準(zhǔn)數(shù)值預(yù)報仍有待深入研究。

3.2 細(xì)長體出水

SPH 方法為細(xì)長體高速入水過程的模擬提供了可靠數(shù)值方案，對于細(xì)長體出水過程其同樣具有較精確的預(yù)報能力。此外，高速出水過程還涉及復(fù)雜氣?液?固三相的耦合，包括瞬態(tài)強(qiáng)沖擊載荷、空泡產(chǎn)生及演化和流?固耦合效應(yīng)等，這也是工程前沿所要研究的重要課題。為此，前人開展了豐富的數(shù)值和試驗研究，包括球體出水[109]和圓柱出水仿真并結(jié)合試驗校核[110]等。

在細(xì)長體高速出水時，其肩部或底部會出現(xiàn)相對背景壓力的負(fù)壓現(xiàn)象，負(fù)壓在SPH 數(shù)值計算中會導(dǎo)致強(qiáng)烈的張力不穩(wěn)定。通過建立張力不穩(wěn)定性控制理論（tensile instability control, TIC）可避免細(xì)長體出水過程的非物理性數(shù)值空洞[5]，而且，還可通過多級粒子分辨率技術(shù)對細(xì)長體附近區(qū)域的粒子進(jìn)行加密，提高局部區(qū)域的精度。 圖12 給出了細(xì)長體高速出水過程的初步模擬結(jié)果。由圖可見，對于高速出水問題研究，SPH 方法已具備較強(qiáng)的預(yù)報能力。

圖12 SPH 方法模擬的細(xì)長體高速出水過程Fig. 12 Numerical simulation of high-speed water exit process of a slender body by SPH method

細(xì)長體從空中入水的瞬態(tài)過程中砰擊載荷占主導(dǎo)地位，計算過程一般忽略黏性力，但因出水過程涉及了空化、湍流、強(qiáng)可壓過程，通常受黏性嚴(yán)重影響，所以出水與入水并非是完全逆向過程[8]。Zhang 和Sun 等[8]采用δ+-SPH 方法研究二維圓柱體出水過程，使用TIC 技術(shù)成功模擬了圓柱體底部產(chǎn)生的渦旋分離，所得結(jié)果比邊界元法（BEM）更精準(zhǔn)。

高速出水問題還需要考慮細(xì)長體出水過程肩部產(chǎn)生的空泡，這就需要開發(fā)SPH 理論框架下的空化模型，精確預(yù)報空泡的產(chǎn)生、演化和潰滅過程。空化問題具有強(qiáng)可壓特性，Sun 等[111]提出的體積自適應(yīng)粒子撕裂與融合理論和技術(shù)在應(yīng)對強(qiáng)可壓的問題時具有較強(qiáng)的適用性，能夠為開發(fā)SPH 空化模型提供基礎(chǔ)支撐。然而，目前傳統(tǒng)的計算流體力學(xué)（CFD）算法模擬空化相變的相關(guān)技術(shù)[112]暫時難以直接拓展應(yīng)用到SPH 算法中，亟需在未來研究中重點開展此項工作。

4 水下爆炸與結(jié)構(gòu)毀傷

水下爆炸是一個涉及高溫高壓爆炸產(chǎn)物劇烈膨脹的強(qiáng)非線性過程[19]，爆炸過程分別與船體結(jié)構(gòu)和海況相關(guān)，研究極具挑戰(zhàn)性，難以通過單純的理論和試驗方法開展研究，而數(shù)值模擬是實現(xiàn)實尺度預(yù)報的最佳策略之一。

從目前來看，遠(yuǎn)場水下爆炸載荷的預(yù)報理論相對完善，但近場和超近場水下爆炸載荷預(yù)報仍十分困難。SPH 方法模擬水下爆炸問題時，可自動捕捉多相界面，模擬從起爆、沖擊波傳播到氣泡膨脹的全過程。當(dāng)SPH 算法與結(jié)構(gòu)求解器進(jìn)行耦合時，可計算船體從損毀到沉沒的全過程[6,113]。本節(jié)將對SPH 理論和方法在研究水下爆炸沖擊波和氣泡脈動方面取得的進(jìn)展進(jìn)行綜述。

4.1 水下爆炸沖擊波模擬

炸藥在水下引爆后會產(chǎn)生高溫高壓的爆炸產(chǎn)物，輻射的高熱量使附近的水瞬時汽化，形成高溫高壓氣泡，并向外輻射形成沖擊波。水下爆炸中沖擊波的威力占據(jù)了總能量的一半左右。爆炸后產(chǎn)生的瞬態(tài)沖擊波在多相中的傳播，具有強(qiáng)可壓、強(qiáng)非線性等特征，這對SPH 算法的精確性和穩(wěn)定性要求極高。而以往的SPH 理論和方法基于不可壓縮[59]或弱可壓假設(shè)[49]，難以直接應(yīng)用于水下爆炸這類強(qiáng)壓縮性問題，尤其是對水下爆炸氣泡的周期性膨脹和收縮過程難以精準(zhǔn)模擬。

Riemann SPH 方法的出現(xiàn)推動了粒子方法處理可壓縮問題的應(yīng)用發(fā)展，但早期在捕捉激波時的精度不如成熟的網(wǎng)格法，在理論方面還需進(jìn)一步提高。因此，基于傳統(tǒng)SPH 理論，借鑒和引入網(wǎng)格法中處理可壓縮流動的先進(jìn)數(shù)值技術(shù)及高精度計算框架，Meng 和Zhang 等[75]提出的TENOSPH 模型與方法，在沖擊波捕捉和分辨方面具備了較高的數(shù)值精度和計算能力，增強(qiáng)了Riemann SPH 方法處理此類問題的魯棒性。Sun 等[111]針對強(qiáng)可壓問題的數(shù)值特征，提出一種新型的粒子體積自適應(yīng)理論，實現(xiàn)了粒子在膨脹作用下撕裂和收縮時的融合模擬，較準(zhǔn)確地預(yù)報了沖擊波載荷。但是，對于氣體膨脹壓縮導(dǎo)致的核函數(shù)影響域內(nèi)粒子數(shù)的變化問題，通常還需采用變光滑長度技術(shù)[114]，來提高SPH 方法的計算精度和穩(wěn)定性。

在水下爆炸的具體計算中，SPH 相關(guān)處理技術(shù)和方法取得了可觀的進(jìn)展。Wang 和Zhang 等[58]提出基于特征線原理的無反射邊界條件，有效消除了水下爆炸沖擊波計算中的邊界反射，為SPH方法用于高精度模擬水下爆炸提供了依據(jù)。彭玉祥[115]采用SPH 方法模擬文獻(xiàn)[116]中相同工況下的二維近平面固壁柱狀水下爆炸算例，研究了高壓氣泡產(chǎn)生的沖擊波和稀疏波的傳播，求解的壓力結(jié)果與參考結(jié)果高度吻合，為SPH 方法用于模擬強(qiáng)可壓瞬態(tài)流動的可行性和有效性提供了參考。同時，文獻(xiàn)[115]針對三維自由場的水下爆炸問題進(jìn)行模擬，將計算得到的沖擊波壓力時歷曲線與經(jīng)驗公式進(jìn)行了對比，兩者吻合較好，體現(xiàn)了SPH 方法預(yù)報沖擊波壓力的良好精度。

在計算三維軸對稱水下爆炸問題時，還可以通過軸對稱SPH 數(shù)值模型[117-119]進(jìn)行模擬，以減小計算量。圖13 給出了采用軸對稱SPH 方法計算得到的水下爆炸幾個典型時刻的沖擊波壓力云圖。軸對稱SPH 方法的計算量小，相對于直接三維模擬，可采用更加精細(xì)的粒子分辨率，因而計算得到的沖擊波波陣面十分銳利。由圖可見，SPH 方法模擬水下爆炸沖擊波階段已發(fā)展得相當(dāng)成熟可靠，具備了水下爆炸沖擊波階段載荷的精確預(yù)報能力。

圖13 軸對稱SPH 方法模擬的水下爆炸幾個典型時刻的沖擊波壓力云圖Fig. 13 Pressure contours of the shock wave induced by underwater explosion at typical time instant by axisymmetric SPH model

4.2 水下爆炸氣泡模擬

May 和Monaghan[120]最早采用SPH 方法研究了氣泡引起的艦船沉沒過程。在此研究領(lǐng)域，文獻(xiàn)[10]主要考慮的是巨大的上浮氣泡，而事實上更不可忽視的是水下爆炸氣泡及其在艦船附近脈動激起的鞭狀運動和射流對結(jié)構(gòu)的沖擊效應(yīng)。研究水下爆炸產(chǎn)生的大尺度脈動氣泡，需考慮環(huán)境中的靜水壓力、大氣壓力、重力及邊界等的耦合作用，以及氣泡脈動產(chǎn)生的熱交換和能量耗散，此機(jī)理十分復(fù)雜。在上述方面，前人提出了理論、試驗和計算等多種方法進(jìn)行研究，取得了豐碩成果[7,121]。

SPH 方法盡管在處理大變形和動邊界問題時具有明顯優(yōu)勢，但用于模擬水下爆炸氣泡脈動較為少見，其主要原因是，爆炸氣泡涉及了爆炸氣體的強(qiáng)壓縮和膨脹過程，給粒子法模擬帶來了挑戰(zhàn)。Joshi 等[122]對固體邊界的空泡崩塌進(jìn)行了模擬，然而忽略了氣體粒子以及氣泡的膨脹過程；Pineda 等[123]也進(jìn)行了多相模型模擬空泡崩塌，但因缺乏粒子體積自適應(yīng)技術(shù)，氣體粒子體積被過度壓縮，導(dǎo)致水?氣界面出現(xiàn)了明顯的鋸齒狀。總之，水下爆炸氣泡SPH 模擬的難點主要體現(xiàn)在氣泡脈動過程中。SPH 理論的拉格朗日粒子特性決定了氣體粒子分布會隨氣泡體積大幅度變化。為簡化問題的復(fù)雜性，目前，SPH 方法暫未計入氣泡脈動過程中的熱力學(xué)效應(yīng)和能量損失，首要任務(wù)是需要提出、驗證和優(yōu)化適用于強(qiáng)可壓問題的SPH 理論，以增強(qiáng)SPH 方法模擬氣泡脈動的可行性，提高對氣泡脈動和射流的仿真精度。

采用SPH 方法計算爆炸氣泡時，粒子不僅會因氣泡膨脹變得十分稀疏，也會因氣泡收縮而變得十分稠密，使得計算的載荷結(jié)果不準(zhǔn)確[111]。如圖14 所示，當(dāng)氣泡體積收縮至最小時（見圖14 左側(cè)），氣體粒子嚴(yán)重集聚，核函數(shù)半徑范圍內(nèi)粒子數(shù)量急劇增多；相反，當(dāng)氣泡體積膨脹到最大時，氣泡粒子的間距隨之增大，核函數(shù)半徑范圍內(nèi)粒子數(shù)量急劇減少（見圖14 右側(cè)）。以上氣泡收縮和膨脹過程中，核函數(shù)內(nèi)部粒子數(shù)量大幅變化，但氣泡外圍的水粒子間距幾乎不變（水的壓縮性較小所導(dǎo)致）。由此，造成水?氣界面的兩側(cè)粒子分布極不均衡，使得水?氣界面不光順，計算精度劇烈下降甚至計算中斷。為解決SPH 方法模擬空泡脈動時氣體粒子的體積和間距劇烈變化所帶來的計算精度及穩(wěn)定性問題，Sun 等[111]建立了基于SPH 體積自適應(yīng)理論（volume adaptive scheme,VAS）的粒子撕裂?融合技術(shù)，成功模擬了水下爆炸氣泡在自由場、自由液面、壁面等邊界下的脈動過程。

圖14 氣泡脈動時氣體粒子體積和間距變化示意圖Fig. 14 Variation of the volume and spacing of gas particles during the pulsation of a bubble

盡管SPH 理論和方法擁有精確捕捉大變形水?氣界面的優(yōu)勢，但在處理三維氣泡脈動問題時，因需要較大的計算域以避免邊界影響而使得計算量較大、效率較低。為此，Sun 等[117]構(gòu)建了軸對稱強(qiáng)可壓多相流SPH 模型，實現(xiàn)了SPH 方法對水下爆炸氣泡脈動的模擬，可精確捕捉高壓氣泡的膨脹、收縮、崩塌、射流等過程。圖15 給出了自由場潰滅氣泡半徑演化和自由液面附近氣泡膨脹及收縮時最大寬度演化的SPH 模擬結(jié)果與參考數(shù)據(jù)[124-126]的對比驗證[117]。由圖可見，軸對稱SPH方法能準(zhǔn)確預(yù)報不同條件下氣泡的膨脹和收縮動態(tài)過程。

圖15 水下爆炸氣泡脈動的模擬SPH 模擬與驗證[117]Fig. 15 SPH simulation and validation of bubble pulsation in underwater explosion [117]

本文采用Sun 等[117]的軸對稱SPH 方法，模擬了Tian 等[127]給出的自由液面以下10 m 處500 kg裝藥所產(chǎn)生的水下爆炸氣泡脈動和射流過程，如圖16 所示。結(jié)果表明，計算得到的氣泡膨脹和射流過程與Tian 等[127]提供的數(shù)值結(jié)果十分吻合。從數(shù)值模擬結(jié)果可以看出，基于Sun 等[111]提出的VAS 理論，采用粒子撕裂?融合技術(shù)，SPH 算法可用于模擬水下爆炸沖擊波和氣泡脈動的全過程。最近，F(xiàn)ang 等[119]將VAS 理論與Riemann SPH 方法結(jié)合應(yīng)用于激波及入水問題，也取得了較好的計算精度，驗證了SPH 理論和方法應(yīng)用于水下爆炸及氣泡動力學(xué)的可行性及有效性，為進(jìn)一步深入研究水下爆炸力學(xué)和載荷機(jī)理提供了更多元的求解方案。

圖16 自由液面附近水下爆炸氣泡脈動和射流的SPH 模擬Fig. 16 SPH simulations of bubble pulsation and jetting near a free surface

4.3 水下爆炸結(jié)構(gòu)毀傷模擬

如上所述，SPH 方法因其不依賴網(wǎng)格，極為適合用于模擬大變形流?固耦合問題。早在1991 年，Libersky 等[128]就將SPH 方法應(yīng)用到了固體結(jié)構(gòu)的沖擊碰撞斷裂的研究中，將SPH 方法的應(yīng)用領(lǐng)域拓寬到了固體力學(xué)領(lǐng)域。早期運用SPH 方法研究固體斷裂時都采用三維實體對結(jié)構(gòu)進(jìn)行建模，但考慮到汽車、飛機(jī)和艦船等工程結(jié)構(gòu)中廣泛采用的是薄殼結(jié)構(gòu)，采用三維實體建模會導(dǎo)致殼體厚度方向的粒子間距極小，從而限制時間積分的增量，使得計算效率降低。在2008 年，Maurel 等[129]提出了SPH 殼單元理論，解決了時間增量受殼體厚度限制的缺點，極大提高了SPH 方法在薄殼結(jié)構(gòu)動響應(yīng)方面的求解能力。此后，Caleyron 等[130]又將該理論應(yīng)用到了對薄殼結(jié)構(gòu)沖擊斷裂響應(yīng)的模擬中。同時，Ming 等[131]應(yīng)用并發(fā)展了該理論，將該理論命名為光滑粒子殼單元（smoothed particle shell，SPS），并應(yīng)用到了水下爆炸結(jié)構(gòu)毀傷的模擬中[132]。

然而，因SPH 理論和方法采用的是強(qiáng)形式的控制方程，經(jīng)過使用移動最小二乘法（moving least squares， MLS）形函數(shù)修正其核函數(shù)后，使之不能自動滿足結(jié)構(gòu)自由邊界條件，需要人工施加自由邊界條件[133-134]，增加了斷裂自由邊模擬的復(fù)雜度。針對此問題，Peng 等[135-136]采用重構(gòu)核粒子方法（reproducing kernel particle method，RKPM）建立了任意曲殼曲梁三維無網(wǎng)格法的計算模型，實現(xiàn)了加筋殼結(jié)構(gòu)的無網(wǎng)格模擬。而且，在上述基礎(chǔ)上，建立了水下爆炸結(jié)構(gòu)毀傷的流?固耦合計算模型[113]，分別采用SPH 和RKPM 方法求解水下爆炸載荷和艦船結(jié)構(gòu)的斷裂響應(yīng)，并與試驗結(jié)果進(jìn)行對比，驗證了所提計算模型的有效性和計算精度。圖17 給出了水下爆炸時艦船結(jié)構(gòu)斷裂毀傷的SPH-RKPM 數(shù)值仿真結(jié)果。由圖可見，該數(shù)值計算模型具備了預(yù)報水下爆炸時的艦船毀傷能力。目前，張阿漫教授團(tuán)隊建立了SPH-RKPM 流?固耦合動力學(xué)模型及計算方法，開發(fā)了具有自主知識產(chǎn)權(quán)的水下爆炸艦船毀傷數(shù)值計算FSLAB軟件[6]，該軟件被應(yīng)用到了艦船結(jié)構(gòu)強(qiáng)沖擊毀傷的模擬之中，廣泛論證了SPH-RKPM 耦合計算方法在水下爆炸領(lǐng)域的應(yīng)用前景，能夠為艦船防護(hù)結(jié)構(gòu)設(shè)計提供理論依據(jù)和基礎(chǔ)性技術(shù)支撐。

圖17 SPH-RKPM 算法模擬水下爆炸引起的艦船結(jié)構(gòu)毀傷過程Fig. 17 Numerical simulations of structural damage resulted from underwater explosion simulated by SPH-RKPM method

4.4 破艙涌流與結(jié)構(gòu)動態(tài)沉沒模擬

SPH 理論和方法除了能夠精確預(yù)報船體在沖擊載荷作用下的結(jié)構(gòu)破損外，對破艙涌流和結(jié)構(gòu)沉沒的模擬也具有強(qiáng)適用性。在船體結(jié)構(gòu)破損情況下，涌流灌入破口后與波浪的相互耦合作用，會使船體產(chǎn)生劇烈的橫搖、升沉等多自由度的復(fù)雜運動。Guo 等[137]采用三維SPH 方法對靜水中不同破口位置的破艙動態(tài)沉沒過程進(jìn)行了仿真，通過相關(guān)試驗驗證了SPH 方法的精度和穩(wěn)定性，為破艙動態(tài)沉沒問題提供了基準(zhǔn)算例，其舷側(cè)破口涌流的仿真結(jié)果與試驗結(jié)果對比如圖18（a）所示。 Cheng 等[138]、Ming 等[139]和Cao 等[140]構(gòu)造了SPH 數(shù)值水池，對多水密艙體在波浪耦合作用下的沉沒過程進(jìn)行研究，準(zhǔn)確模擬了艙內(nèi)涌流晃動和外部波浪的耦合過程，進(jìn)一步完善了破損結(jié)構(gòu)在波浪環(huán)境下的動態(tài)沉沒力學(xué)機(jī)理，如圖18（b）所示的沉沒過程[138]。綜上，建議未來考慮艙內(nèi)氣體的真實壓縮性，基于VAS 框架的多級粒子分辨率技術(shù)和多相流模型[141]，進(jìn)一步開發(fā)適用于破艙涌流的高精度SPH 理論和方法，以提高模擬精度和穩(wěn)定性。SPH 方法用于水下爆炸的各環(huán)節(jié)及船體破口與動態(tài)沉沒過程均具有較強(qiáng)的數(shù)值模擬能力，預(yù)期未來能夠?qū)崿F(xiàn)將SPH 方法應(yīng)用到水下爆炸全過程的精準(zhǔn)仿真。

圖18 三維SPH 方法模擬破損艙段涌流和沉沒過程的結(jié)果與試驗結(jié)果對比Fig. 18 Comparison between experimental and 3D SPH simulation results of flowing and sinking process due to a damaged ship cabin

5 高速水動力學(xué)耦合算法

考慮到SPH 數(shù)值模擬在實際工程應(yīng)用中難以規(guī)避計算成本的問題，為充分發(fā)揮SPH 模擬大變形問題的優(yōu)勢，提高其在船海力學(xué)領(lǐng)域的應(yīng)用范圍，建議SPH 方法在未來發(fā)展上應(yīng)重視與其他算法的耦合。例如，開發(fā)有/無網(wǎng)格耦合算法、黏流?勢流耦合算法等，構(gòu)建如圖19 所示的耦合計算框架，以解決船海高速水動力領(lǐng)域常見的出水和入水以及水下爆炸與沖擊問題。圖19 所示的耦合算法可在流?固耦合作用劇烈的區(qū)域采用SPH 粒子法進(jìn)行局部范圍的精準(zhǔn)計算，在較小變形的外圍區(qū)域使用黏流網(wǎng)格法進(jìn)行高效求解，在遠(yuǎn)場區(qū)域可采用勢流方法進(jìn)行大范圍海域的模擬。事實上，文獻(xiàn)[142] 中的部分耦合算法已為該框架的搭建和實現(xiàn)奠定了堅實的技術(shù)基礎(chǔ)，例如，SPH-FVM 算法[143]在大變形流動仿真上擁有較好的穩(wěn)定性，已拓展應(yīng)用至航行體的興波、高速入水等問題。此外，高階譜（higher-order spectral, HOS）方法與網(wǎng)格法求解器的耦合也得到了驗證[144]，其適合用于對大范圍海域內(nèi)的船海水動力問題進(jìn)行黏流?勢流耦合計算。圖19 所示的多算法耦合計算框架通過結(jié)合有網(wǎng)格、無網(wǎng)格和黏流、勢流數(shù)值算法的優(yōu)點進(jìn)行耦合求解，有望為SPH 理論和方法應(yīng)用于實際船海工程問題提供新的技術(shù)方案，為實現(xiàn)高精度、高效率的船海工程應(yīng)用仿真奠定基礎(chǔ)。

圖19 有/無網(wǎng)格耦合算法和黏流?勢流耦合算法求解典型高速水動力問題示意圖Fig. 19 Schematics of coupled algorithm between mesh-based and particle-based methods and between viscous-flow and potential-flow solvers for typical high-speed hydrodynamic problems

6 結(jié) 語

本文綜述了SPH 理論與方法在船海高速水動力學(xué)研究中的最新進(jìn)展。可見，SPH 理論和方法日趨得到完善，相關(guān)數(shù)值方法在高速艦船興波、細(xì)長體出水和入水、水下爆炸與沖擊等高速水動力問題中已得到廣泛應(yīng)用，取得了較高的計算精度，在精確捕捉大變形的多相界面和流?固界面等方面具備明顯優(yōu)勢。

在未來研究中，SPH 理論與方法在以下幾個方面仍有提高的空間，主要包括：

1） 在計算精度和效率方面，進(jìn)一步提高SPH 理論和方法的完備性及數(shù)值格式的收斂階數(shù)，實現(xiàn)以較少粒子數(shù)量達(dá)到預(yù)期的計算精度，提高計算效率。

2） 在邊界處理方面，仍需進(jìn)一步發(fā)展三維復(fù)雜構(gòu)型壁面邊界的SPH 數(shù)值處理技術(shù)，減弱邊界效應(yīng)對物理問題模擬結(jié)果的影響。

3） 在應(yīng)用范圍和領(lǐng)域方面，開發(fā)SPH 理論框架下的湍流模型和空化模型，加強(qiáng)高雷諾數(shù)流動和相變問題的模擬能力。同時，發(fā)展基于體積分?jǐn)?shù)的SPH 水氣混合物的理論和數(shù)值模型，進(jìn)一步拓展應(yīng)用于高速航行體的興波和水汽泡混合流模擬中。

4） 在解決實際工程問題方面，進(jìn)一步提高SPH 方法和計算程序的高性能并行計算能力，與國產(chǎn)大型超級計算機(jī)的能力結(jié)合，實現(xiàn)億級粒子數(shù)的大規(guī)模并行化模擬，增強(qiáng)SPH 理論和方法的工程化應(yīng)用能力。

7 致 謝

感謝哈爾濱工程大學(xué)明付仁副教授、王平平博士和中國海洋大學(xué)程晗博士對本文水下爆炸沖擊波和細(xì)長體出入水計算提供的幫助。