歐訓(xùn)勇

　?。êＤ蠠釒ШＱ髮W院 電子通信工程學院,， 海南 三亞 572022）

       摘要：首先介紹孤立波的KdV方程，繼而討論了孤立波SPH方法的數(shù)值求解過程,，選擇SPH光滑核函數(shù)作為正則化高斯核函數(shù),。分析了數(shù)值求解過程的時間積分方法，給出了具體計算公式,，最后給出相應(yīng)程序中的具體參數(shù)下孤立波運動模擬效果,。

　　關(guān)鍵詞：孤立波；自由表面流體,；SPH

　　中圖分類號：TP311.1文獻標識碼：ADOI： 10.19358/j.issn.1674 7720.2016.20.019

　　引用格式：歐訓(xùn)勇. 應(yīng)用SPH方法實現(xiàn)海面孤立波運動的模擬［J］.微型機與應(yīng)用,，2016,35（20）：69 71，74.

0引言

　　孤立波是1834年由RUSSELL S觀察發(fā)現(xiàn)到的,，而Korteweg和de Vies于1895第一次通過理論模型對孤立波進行了解釋,，被稱為KdV理論。該理論基于弱色散淺水波理論,，淺水波色散由非線性效果平衡,，使得波在傳播過程中在任意距離都保持著一定的振幅和形狀,。KdV方程的精確解描述了孤立子的形狀和傳播速度,。盡管KdV理論被認為是一階近似解，很好地描述了真正的孤立波,，而高階近似解也一樣可以構(gòu)成,。參考文獻［1］中介紹了任意高階迭代，逐次逼近模型,，在第一次迭代步驟中再現(xiàn)了KdV理論,，然而，更高階求解需要數(shù)值方法,。

　　光滑粒子流體動力學（以下簡稱SPH）是一種無網(wǎng)格拉格朗日型數(shù)值方法,，由LUCY L于1977年提出。該方法最初應(yīng)用于天體物理學領(lǐng)域,，然而受海岸工程問題所驅(qū)動,，由MONAGHAN J J在1994第一次提出用于模擬流體流動。

　　在過去的二十多年中，由于在模擬流體自由表面流動有突出的性能,，SPH方法成為工程應(yīng)用領(lǐng)域最流行的粒子數(shù)值方法,，例如近海波模擬、海嘯,。本文討論利用SPH方法實現(xiàn)海面孤立波運動的三維模擬效果,。

1孤立波的方程

　　孤立波傳播的零階近似解可以用線性波方程描述，而淺水中的波速由以下式子給出：

　　式中,，g為重力加速度,，H為水域的深度。式（1）給出了孤立波傳播的一個粗略近似值,，而忽略一些特殊的性能,，如波的實際振幅和寬度，孤立波傳播過程如圖1所示,。圖1中A為波幅,，η(x)是自由表面波形狀函數(shù)。式（1）僅當A<<H時才成立,。

　　線性波的傳播方程是沒有孤立波解的,。KdV方程如下：

　　

　　式（2）適用于以不同幾何結(jié)構(gòu)構(gòu)造自由表面孤立子的形狀。函數(shù)η(x,t)表示給定位置x的表面高度,，是一個與時間有關(guān)的函數(shù),。單個自由表面孤立波的KdV方程的精確解由波的形狀函數(shù)給出，如下：

　　式中,，A是振幅,，a是孤立子水平偏移量，k為有效的波數(shù),，其取值如下：

　　一階孤立波的傳播速度計算公式為：

　　二階孤立波的傳播速度由HAL′ASZ G B做了修正［2］,，其公式如下：

2流體控制方程

　　在流體力學中，廣泛應(yīng)用歐拉方程和連續(xù)方程來描述流體運動,。描述流體運動的偏微分方程中,，局部和對流通量都包含在拉格朗日微分方程中，如下式：

　　其中,，Φ表示一個任意的標量或矢量場,。利用上式的微分算子可以得到無粘性流體力學方程如下：

　　式（8）中的v、ρ,、P,、v、g分別指流體的速度,、密度,、壓力,、運動粘度和重力加速度。在弱可壓縮流體中式(9)被定義為流體密度和壓力之間的關(guān)系：

　　盡管上述各類方程包括偏微分方程,、波傳播方程等有解析解,，但是通常情況下仍不能得證通過恰當?shù)臄?shù)值方法得到精確解，只能得到某種程度上的近似解,。然而,，這些近似方法經(jīng)常取得不利的數(shù)值特性，因此它們的一般性,、魯棒性,、實用性都受限制?？紤]層流的無粘性,，現(xiàn)今的動力流體建模都盡量避免對復(fù)雜湍流模型建模。

3孤立波的SPH數(shù)值模擬

　　3.1SPH方程

　　無網(wǎng)格拉格朗日數(shù)值方法稱為光滑粒子流體動力學是一種適用于解決式（8）方程系統(tǒng)的工具,。SPH近似數(shù)值方法是基于流體節(jié)點,，該節(jié)點稱為粒子，它們在空間中運動時,，每個粒子都攜帶有相應(yīng)的物理量,，如質(zhì)量、密度,、壓力,、速度等。這種離散方法是基于受域加權(quán)插值的,，在一個給定的點,，使用臨近區(qū)域內(nèi)的粒子，這些粒子受一個稱為光滑核函數(shù)W(ri-rj,h)的控制,，形成一個離散卷積,。其函數(shù)如下：

　　式中，i指當前粒子,，j是鄰域內(nèi)的一個粒子,，fi=f(ri)是任意流場中粒子i的位置ri，核函數(shù)Wij=W(ri-rj,，h)帶有緊支性或無限影響半徑，h稱為光滑長度,，Vj是指定粒子j的物理量值,，N是核函數(shù)Wij影響范圍內(nèi)的粒子數(shù)。公式（10）以離散化構(gòu)造一個任意流場,，流場以粒子散布在空間里,。文中討論的計算過程即正則化高斯核函數(shù)［3］，公式如下：

　　式（11）為改進的高斯核函數(shù)。式中r=|ri-rj|,。常量C0和C1由下式計算：

　　本文討論中,，影響域δ取值為3h。相應(yīng)地,，兩個一階微分算子梯度和旋度如下：

　　在SPH數(shù)值方法中通過插入數(shù)值擴散項到連續(xù)運動的方程中以保持數(shù)值穩(wěn)定性是一種較流行的實踐方法,。由于自由表面孤立子受慣性力驅(qū)動，表現(xiàn)出粘性行為,，動量擴散（物理數(shù)值）在目前的工作中可以忽略,。相反，在連續(xù)方程中密度的數(shù)值擴散項是由文獻［3］算出來的,，在參考文獻［4］中得以改善,。參考文獻［4］中ANTUONO M提出的基于線性穩(wěn)定分析，密度擴散項成為一個高效的數(shù)值阻尼振蕩工具,。

　　可壓縮性作為另一種特殊的SPH數(shù)值屬性,，可受控于近似弱可壓縮流體方程，這種狀態(tài)假設(shè)一個正壓流體的壓力和密度之間的線性關(guān)系,。本文使用的SPH數(shù)值方法的離散動力學方程如下：

　　式（14）中,，ρ0是參考密度值，f是外力之和(包括重力),，cs是聲波傳播速度,。第一公式右邊第二項是人工密度擴散項，在建模中常稱為δ SPH,，經(jīng)驗系數(shù)ξ=0.1,。Ψij項的計算公式如下：

　　式（15）中右邊第二項帶有正則化密度梯度確保流體質(zhì)量守恒,包括自由表面邊界，它通過下式計算：

　　式（16）中表示張量積,。正則張量L在流體邊界影響著離散拉普拉斯收斂,，它以修正由內(nèi)核截斷引起的離散梯度人工數(shù)值達到收斂目的。

　　要降低計算量,，弱可壓縮流體模型通常運行在聲速范圍,，但是其量要足夠大，以保持在預(yù)定義范圍,、獨立慣性和聲波內(nèi)有足夠大的絕對密度偏差,。通常取值較現(xiàn)在經(jīng)典流體速度幅值的10倍大，其計算公式和式（1）相差10倍,，具體如下：

　　式(17)中M取值為10,。

　　利用SPH方法實現(xiàn)孤立波運動的模擬效果，整個實現(xiàn)過程涉及求解孤立波SPH數(shù)值解,，確定SPH的邊界及初始條件,，以及孤立波運動的時間積分過程,。

　　3.2SPH方法的邊界和初始條件

　　SPH方法的顯著效益（至少在流體建模中）是把任意自由曲面當作自然邊界處理而不需要任何額外的計算負擔。另外,，如果流體簡單地連接著空氣,，完全可以從計算域中忽略，因為是恒壓和水密度量,。注意,，在復(fù)雜流動情況下，如破浪,，空氣可能起著重要作用,，因此它不應(yīng)該被忽視。本文中應(yīng)用兩種SPH不同的邊界條件：一種是邊界墻體和底部,，另一種是周期性的邊界,，允許在無限領(lǐng)域執(zhí)行更一般的計算［5-6］。

　　這里的周期邊界的基本量是在翼展方向形成的2D寬域展向近似平面,，帶有三維求解,。在SPH的固體邊界的模型中有幾個基礎(chǔ)性質(zhì)的變種。

　　3.3SPH方法的時間步積分

　　式（14）可以通過任意形式求解,，是穩(wěn)定的數(shù)值解析方法,。本文研究應(yīng)用二階預(yù)修正方法，第一步粒子按△t/2,，具體計算公式如下：

　　在中間狀態(tài),，密度導(dǎo)數(shù)、壓力,、外力,、粒子內(nèi)力（加速度）進行近似估算。使用新的數(shù)值以全步時間長度推進粒子運動,，其計算公式如下：

　　為了降低計算性能的要求,，保持數(shù)值穩(wěn)定，時間步長可以在每幀中自適應(yīng)選擇,。在當前的SPH模型使用CFL條件執(zhí)行計算,。

　　上式中，CFL=0.2,，Vij=Vi-Vj

4程序運行結(jié)果

　　根據(jù)以上各節(jié)介紹的方程及數(shù)值求解的過程,，利用OpenGL三維圖形庫，在C-free 5環(huán)境下,，使用C++面向?qū)ο蟮木幊谭椒▽崿F(xiàn)了對在淺水灘運動的孤立波模擬,。由于考慮計算量，SPH程序使用粒子總數(shù)為4 096個,，即為4K,。粒子質(zhì)量為0.000 205 43，密度為600,，光滑長度為0.01,，粒子半徑0.004，粒子影響間距為0.005 9,，流體黏度系數(shù)為0.2,。程序運行效果截圖如圖2所示。

5結(jié)論

　　程序運行環(huán)境為：Intel(R) Core(TM) i32348M @ 2.30 GHz CPU,，4 GB內(nèi)存,，模擬的海浪孤立波運行效果非常流暢。以上方法實現(xiàn)的算法中SPH粒子數(shù)量可達8 000個以上,。超過8 000個粒子后,，可看到模擬動畫效果出現(xiàn)卡頓了。要想達到更大的粒子數(shù)量,，采用GPU方法能取得更佳效果,。本研究在應(yīng)用于構(gòu)建模擬大水域海浪運動，將繼續(xù)深入探索GPU方法及網(wǎng)絡(luò)分塊協(xié)同渲染描繪更大水域波浪運動,。研究成果將應(yīng)用于構(gòu)建光滑粒子流體動力學方法下的船舶運動的虛擬仿真系統(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.

　?。?］ ANTUONO M, COLAGROSSI A, MARRONE S. Numerical diffusive terms in weakly compressible SPH schemes［J］. Computer Physics Communications， 2012,，183(12):2570-2580.

　?。?］ 劉瑛琦. 基于SPH方法的數(shù)值波浪水槽研究［D］.南京：河海大學，2006.

    ［6］ 高睿,，任冰,，王國玉，等. 孤立波淺化過程的SPH數(shù)值模擬［J］. 水動力學研究與進展（A輯）,2010,25(5):620-629.