歐訓勇
（海南熱帶海洋學院 電子通信工程學院， 海南 三亞 572022）
摘要：首先介紹孤立波的KdV方程，繼而討論了孤立波SPH方法的數(shù)值求解過程，選擇SPH光滑核函數(shù)作為正則化高斯核函數(shù)。分析了數(shù)值求解過程的時間積分方法，給出了具體計算公式，最后給出相應程序中的具體參數(shù)下孤立波運動模擬效果。
關鍵詞：孤立波；自由表面流體；SPH
中圖分類號：TP311.1文獻標識碼：ADOI： 10.19358/j.issn.1674 7720.2016.20.019
引用格式：歐訓勇. 應用SPH方法實現(xiàn)海面孤立波運動的模擬［J］.微型機與應用，2016,35（20）：69 71，74.
0引言
孤立波是1834年由RUSSELL S觀察發(fā)現(xiàn)到的，而Korteweg和de Vies于1895第一次通過理論模型對孤立波進行了解釋，被稱為KdV理論。該理論基于弱色散淺水波理論，淺水波色散由非線性效果平衡，使得波在傳播過程中在任意距離都保持著一定的振幅和形狀。KdV方程的精確解描述了孤立子的形狀和傳播速度。盡管KdV理論被認為是一階近似解，很好地描述了真正的孤立波，而高階近似解也一樣可以構成。參考文獻［1］中介紹了任意高階迭代，逐次逼近模型，在第一次迭代步驟中再現(xiàn)了KdV理論，然而，更高階求解需要數(shù)值方法。
光滑粒子流體動力學（以下簡稱SPH）是一種無網格拉格朗日型數(shù)值方法，由LUCY L于1977年提出。該方法最初應用于天體物理學領域，然而受海岸工程問題所驅動，由MONAGHAN J J在1994第一次提出用于模擬流體流動。
在過去的二十多年中，由于在模擬流體自由表面流動有突出的性能，SPH方法成為工程應用領域最流行的粒子數(shù)值方法，例如近海波模擬、海嘯。本文討論利用SPH方法實現(xiàn)海面孤立波運動的三維模擬效果。
1孤立波的方程
孤立波傳播的零階近似解可以用線性波方程描述，而淺水中的波速由以下式子給出：
式中，g為重力加速度，H為水域的深度。式（1）給出了孤立波傳播的一個粗略近似值，而忽略一些特殊的性能，如波的實際振幅和寬度，孤立波傳播過程如圖1所示。圖1中A為波幅，η(x)是自由表面波形狀函數(shù)。式（1）僅當A<<H時才成立。
線性波的傳播方程是沒有孤立波解的。KdV方程如下：
式（2）適用于以不同幾何結構構造自由表面孤立子的形狀。函數(shù)η(x,t)表示給定位置x的表面高度，是一個與時間有關的函數(shù)。單個自由表面孤立波的KdV方程的精確解由波的形狀函數(shù)給出，如下：
式中，A是振幅，a是孤立子水平偏移量，k為有效的波數(shù)，其取值如下：
一階孤立波的傳播速度計算公式為：
二階孤立波的傳播速度由HAL′ASZ G B做了修正［2］，其公式如下：
2流體控制方程
在流體力學中，廣泛應用歐拉方程和連續(xù)方程來描述流體運動。描述流體運動的偏微分方程中，局部和對流通量都包含在拉格朗日微分方程中，如下式：
其中，Φ表示一個任意的標量或矢量場。利用上式的微分算子可以得到無粘性流體力學方程如下：
式（8）中的v、ρ、P、v、g分別指流體的速度、密度、壓力、運動粘度和重力加速度。在弱可壓縮流體中式(9)被定義為流體密度和壓力之間的關系：
盡管上述各類方程包括偏微分方程、波傳播方程等有解析解，但是通常情況下仍不能得證通過恰當?shù)臄?shù)值方法得到精確解，只能得到某種程度上的近似解。然而，這些近似方法經常取得不利的數(shù)值特性，因此它們的一般性、魯棒性、實用性都受限制。考慮層流的無粘性，現(xiàn)今的動力流體建模都盡量避免對復雜湍流模型建模。
3孤立波的SPH數(shù)值模擬
3.1SPH方程
無網格拉格朗日數(shù)值方法稱為光滑粒子流體動力學是一種適用于解決式（8）方程系統(tǒng)的工具。SPH近似數(shù)值方法是基于流體節(jié)點，該節(jié)點稱為粒子，它們在空間中運動時，每個粒子都攜帶有相應的物理量，如質量、密度、壓力、速度等。這種離散方法是基于受域加權插值的，在一個給定的點，使用臨近區(qū)域內的粒子，這些粒子受一個稱為光滑核函數(shù)W(ri-rj,h)的控制，形成一個離散卷積。其函數(shù)如下：
式中，i指當前粒子，j是鄰域內的一個粒子，fi=f(ri)是任意流場中粒子i的位置ri，核函數(shù)Wij=W(ri-rj，h)帶有緊支性或無限影響半徑，h稱為光滑長度，Vj是指定粒子j的物理量值，N是核函數(shù)Wij影響范圍內的粒子數(shù)。公式（10）以離散化構造一個任意流場，流場以粒子散布在空間里。文中討論的計算過程即正則化高斯核函數(shù)［3］，公式如下：
式（11）為改進的高斯核函數(shù)。式中r=|ri-rj|。常量C0和C1由下式計算：
本文討論中，影響域δ取值為3h。相應地，兩個一階微分算子梯度和旋度如下：
在SPH數(shù)值方法中通過插入數(shù)值擴散項到連續(xù)運動的方程中以保持數(shù)值穩(wěn)定性是一種較流行的實踐方法。由于自由表面孤立子受慣性力驅動，表現(xiàn)出粘性行為，動量擴散（物理數(shù)值）在目前的工作中可以忽略。相反，在連續(xù)方程中密度的數(shù)值擴散項是由文獻［3］算出來的，在參考文獻［4］中得以改善。參考文獻［4］中ANTUONO M提出的基于線性穩(wěn)定分析，密度擴散項成為一個高效的數(shù)值阻尼振蕩工具。
可壓縮性作為另一種特殊的SPH數(shù)值屬性，可受控于近似弱可壓縮流體方程，這種狀態(tài)假設一個正壓流體的壓力和密度之間的線性關系。本文使用的SPH數(shù)值方法的離散動力學方程如下：
式（14）中，ρ0是參考密度值，f是外力之和(包括重力)，cs是聲波傳播速度。第一公式右邊第二項是人工密度擴散項，在建模中常稱為δ SPH，經驗系數(shù)ξ=0.1。Ψij項的計算公式如下：
式（15）中右邊第二項
帶有正則化密度梯度確保流體質量守恒,包括自由表面邊界，它通過下式計算：
式（16）中表示張量積。正則張量L在流體邊界影響著離散拉普拉斯收斂，它以修正由內核截斷引起的離散梯度人工數(shù)值達到收斂目的。
要降低計算量，弱可壓縮流體模型通常運行在聲速范圍，但是其量要足夠大，以保持在預定義范圍、獨立慣性和聲波內有足夠大的絕對密度偏差。通常取值較現(xiàn)在經典流體速度幅值的10倍大，其計算公式和式（1）相差10倍，具體如下：
式(17)中M取值為10。
利用SPH方法實現(xiàn)孤立波運動的模擬效果，整個實現(xiàn)過程涉及求解孤立波SPH數(shù)值解，確定SPH的邊界及初始條件，以及孤立波運動的時間積分過程。
3.2SPH方法的邊界和初始條件
SPH方法的顯著效益（至少在流體建模中）是把任意自由曲面當作自然邊界處理而不需要任何額外的計算負擔。另外，如果流體簡單地連接著空氣，完全可以從計算域中忽略，因為是恒壓和水密度量。注意，在復雜流動情況下，如破浪，空氣可能起著重要作用，因此它不應該被忽視。本文中應用兩種SPH不同的邊界條件：一種是邊界墻體和底部，另一種是周期性的邊界，允許在無限領域執(zhí)行更一般的計算［5-6］。
這里的周期邊界的基本量是在翼展方向形成的2D寬域展向近似平面，帶有三維求解。在SPH的固體邊界的模型中有幾個基礎性質的變種。
3.3SPH方法的時間步積分
式（14）可以通過任意形式求解，是穩(wěn)定的數(shù)值解析方法。本文研究應用二階預修正方法，第一步粒子按△t/2，具體計算公式如下：
在中間狀態(tài)，密度導數(shù)、壓力、外力、粒子內力（加速度）進行近似估算。使用新的數(shù)值以全步時間長度推進粒子運動，其計算公式如下：
為了降低計算性能的要求，保持數(shù)值穩(wěn)定，時間步長可以在每幀中自適應選擇。在當前的SPH模型使用CFL條件執(zhí)行計算。
上式中，CFL=0.2，Vij=Vi-Vj
4程序運行結果
根據以上各節(jié)介紹的方程及數(shù)值求解的過程，利用OpenGL三維圖形庫，在C-free 5環(huán)境下，使用C++面向對象的編程方法實現(xiàn)了對在淺水灘運動的孤立波模擬。由于考慮計算量，SPH程序使用粒子總數(shù)為4 096個，即為4K。粒子質量為0.000 205 43，密度為600，光滑長度為0.01，粒子半徑0.004，粒子影響間距為0.005 9，流體黏度系數(shù)為0.2。程序運行效果截圖如圖2所示。
5結論
程序運行環(huán)境為：Intel(R) Core(TM) i32348M @ 2.30 GHz CPU，4 GB內存，模擬的海浪孤立波運行效果非常流暢。以上方法實現(xiàn)的算法中SPH粒子數(shù)量可達8 000個以上。超過8 000個粒子后，可看到模擬動畫效果出現(xiàn)卡頓了。要想達到更大的粒子數(shù)量，采用GPU方法能取得更佳效果。本研究在應用于構建模擬大水域海浪運動，將繼續(xù)深入探索GPU方法及網絡分塊協(xié)同渲染描繪更大水域波浪運動。研究成果將應用于構建光滑粒子流體動力學方法下的船舶運動的虛擬仿真系統(tǒng)。
參考文獻
［1］ MOLTENI D, COLAGROSSI A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH［J］. Computer Physics Communications， 2009，180（6）: 861-872.
［2］ HAL′ASZ G B. Higher order corrections for shallow water solitary waves: elementary derivation and experiments［J］. European Journal of Physics， 2009，30（6）：1311-1323.
［3］ ANTUONO M, COLAGROSSI A, MARRONE S, et al. Free surface flows solved by means of SPH schemes with numerical diffusive terms［J］. Ccomputer Physics Communications,2010，181（3）:532-549.
［4］ ANTUONO M, COLAGROSSI A, MARRONE S. Numerical diffusive terms in weakly compressible SPH schemes［J］. Computer Physics Communications， 2012，183(12):2570-2580.
［5］ 劉瑛琦. 基于SPH方法的數(shù)值波浪水槽研究［D］.南京：河海大學，2006.
［6］ 高睿，任冰，王國玉，等. 孤立波淺化過程的SPH數(shù)值模擬［J］. 水動力學研究與進展（A輯）,2010,25(5):620-629.