歐訓勇，陳美伊，鮑海琴，康小平

(海南熱帶海洋學院 海洋信息工程學院，海南 三亞 572022)

?

三維修正SPH方法實現(xiàn)的海水運動效果模擬*

歐訓勇，陳美伊，鮑海琴，康小平

(海南熱帶海洋學院 海洋信息工程學院，海南 三亞 572022)

圍繞SPH的標準控制方程，對SPH的標準核函數(shù)使用了所謂的線性再生核近似方法進行修正，這種方式所取得的梯度效果令人滿意。修正的SPH方法模型用來模擬流體自由表面運動。仿真效果展現(xiàn)了水平面上波的形成和傳播過程。最后在OpenGL環(huán)境下實現(xiàn)了的逼真模擬效果。

SPH；自由表面流體；流體運動

0 引言

光滑流體動力學方法(以下稱為SPH)是1977年由Lucy、Gingold和Monaghan率先提出的，至今該方法已經(jīng)得到了廣泛的研究及應(yīng)用，并成為通用性和有效性最強的應(yīng)用于物理、力學、工程等領(lǐng)域計算離散粒子的方法之一。SPH方法成為有重要意義的計算方法可以歸因于它幾個突出的優(yōu)勢：概念簡單、容易編程實現(xiàn)、處理大變形問題。作為一種無網(wǎng)格的純拉格朗日方法，SPH方法已成功應(yīng)用于一系列現(xiàn)象的數(shù)值模擬，如動態(tài)沖擊、金屬損傷和斷裂、海浪運動、傳熱過程、水下爆炸等。

當然，SPH方法也被用來解決各種流體力學問題，包括有重要的水文現(xiàn)象和離岸工程等。例如，文獻[1]中運用這一方法研究水上的重力流和孤立波；文獻[2]中分析了在一個垂直墻面和傾斜的海灘上的孤立波運動過程；文獻[3]研究了經(jīng)典的潰壩問題，包括水撞擊墻體的分析；文獻[4-5]使用SPH方法研究了物體入水等自由表面流動問題；文獻[6]研究利用SPH方法模擬了二維和三維波浪在傾斜的海灘上的破碎問題；文獻[7]利用SPH研究了彈性流體結(jié)構(gòu)相互作用的問題；文獻[8]利用不可壓縮流的SPH方法研究流體自由表面運動的模擬過程。

以上所提到關(guān)于對于流體力學問題的研究，SPH方法的各種變體已經(jīng)應(yīng)用到了很多領(lǐng)域之中。早期的研究成果中，有代表性的，例如文獻[1]中，傳統(tǒng)的SPH方法是以水作為一個弱可壓縮流體[1]，和其中所有的場變量利用一個單一的插值函數(shù)(稱為SPH光滑核)，在其作用域內(nèi)近似離散成粒子。在后來的研究文獻中，不斷提出各種改進的方法，以實現(xiàn)更好的真實流體建模，并加強了計算的準確性[4-5]。本文提出一種修正的SPH方法模擬海水自由表面運動。

1 控制方程

本文中，把水視為一個可壓縮的粘性液體，并假設(shè)它的運動發(fā)生在等溫條件下，完全是由重力引發(fā)的運動。流體質(zhì)量守恒的連續(xù)方程如下：

(1)

(2)

式中，p是流體壓力，b指的是外力，本文所指的外力在沒有其他作用力時其實就是重力。如果考慮其他力，那么b就是這些力和重力之和。

假設(shè)流體受正壓，在這種情況下，壓力是由流體密度唯一確定。根據(jù)傳統(tǒng)的SPH方法，壓力p和密度之間的變化規(guī)律如下：

(3)

式中，P0和δ0壓力和密度的參考值，參數(shù)τ的取值如果為水流體是7，為空氣則1.4。相應(yīng)地，聲音的速度表示為：

(4)

方程(3)的突出作用是表明了壓力和密度的直接關(guān)系，這就使得沒有必要解決任何額外的壓力方程(這是不可壓縮流的泊松方程)。因此，快速顯示的時間步長算法可以在計算中采用，避免在泊松方程的影響下使用隱式方法。

2 SPH理論和方程

在SPH方法中，任何系統(tǒng)都認為是由一套隨意離散分布的粒子構(gòu)成的，每個粒子都承載相關(guān)的物理信息，如：質(zhì)量、密度、動量、溫度等。由于粒子分布是離散的，獨立運動，因此該方法具有完全無網(wǎng)格的特征。場變量是通過特定的插值函數(shù)來近似的，這就是光滑核。這些函數(shù)表述了給定空間或質(zhì)點的場量值。通常核函數(shù)由限定域的非零函數(shù)構(gòu)成，這個限定域稱為支持域。設(shè)x為位置矢量，a和b為兩個粒子，場函數(shù)f(x)在粒子a處的取值由核函數(shù)W以粒子a為中心近似估算取得。由于核支持域的緊支性，只有位于這個域內(nèi)的粒子才有被插值累加。其方程如下：

(5)

式中，fb=f(xb)是粒子b的離散值，Vb是粒子b的體積，N是核支持域內(nèi)的粒子數(shù)。

光滑核W是帶兩個參數(shù)的函數(shù)，r是兩粒子a和b之間的距離，計算式子為：

(6)

h是核的光滑長度，取決于支持域的大小和核插值粒子數(shù)的離散量。

核函數(shù)的表示有多種形式，使用最為廣泛的核函數(shù)有高斯核函數(shù)、三次樣條函數(shù)、五次樣條函數(shù)。本文討論的SPH方法使用的是三次樣條核函數(shù)，其形式如下：

(7)

2.1 SPH標準方程

為了導出給定問題域的離散方程，需要微分算子離散近似，SPH方法中有幾種方法用來構(gòu)造函數(shù)梯度和微分算子的近似計算。文中使用MONAGHAN JJ提出的計算公式[1]，計算某粒子a的標量場函數(shù)的梯度，其近似的插值公式如下：

(8)

類似地，粒子a的矢量場函數(shù)的散度計算公式為：

(9)

(10)

式(10)是核函數(shù)W的在粒子a處的梯度計算式。由于式(7)是對稱的，因此Wab=Wba，式(10)有如下關(guān)系：

(11)

將式(1)和式(2)用式(8)和式(9)的離散形式表示，可得計算形式如下：

(12)

(13)

2.2 修正的SPH方法

(14)

(15)

式(14)中的系數(shù)cαβ(x)(α,β=1,2,3,4)是由應(yīng)用完整條件的插值函數(shù)及它們的導數(shù)確定的，這些函數(shù)由BELYTSCHKOT等人于1998年推導出[9]。由這些條件構(gòu)造的矩陣方程為：

AcT=I

(16)

式中，c為系數(shù)cαβ構(gòu)成的矩陣，I為單位矩陣，經(jīng)過變換A的表達式為：

(17)

在三維中，上述的修正過程要求對每個粒子a生成矩陣A及反演。利用式(14)～(16)可計算求得系數(shù)c，將修正的核函數(shù)及梯度代入式(12)～(13)得到修正核函數(shù)后的動量守恒方程。

3 SPH方法的時間步積分

數(shù)值求解的問題可以利用以下一組一階微分方程來取得：

(18)

式中，Ma和Fa遵從質(zhì)量守恒和動量守恒定律。根據(jù)式(12)和(13)的離散形式，利用修正的核函數(shù)及梯度，這兩個場量的新的計算式子如下：

(19)

(20)

式(18)中的第三個方程中的Ua是粒子a的速度修正量，該變量值采用Monaghan提出的公式進行求?。?/p>

(21)

式中，粒子a的速度是通過位于粒子a的支持域內(nèi)的所有臨近粒子b的光滑處理平均值。利用這一做法，可以防止粒子間相互滲透，和傳統(tǒng)的SPH方法對比更能確保粒子有序和有規(guī)律的運動。根據(jù)文獻[1]，式(21)中ε取值在0和1之間。

式(18)中方程組的數(shù)值求解是應(yīng)用顯式時間步的典型的預估校正類型，每個時間步分兩個階段計算。設(shè)tk和tk+1是兩個瞬時時間點，時間步長度為Δt=tk+1-tk。第一個計算階段在tk+1/2=(tk+tk+1)/2，計算取得粒子的密度、速度和位置，其計算規(guī)律如下：

(22)

第二個計算階段，t=tk+1，其計算規(guī)律如下：

(23)

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

在搜索粒子a的臨近粒子時，使用鏈式算法。計算范圍中空間區(qū)域按照匹配的支持域劃分為網(wǎng)格間距。顯然地，每個粒子僅與它四周相鄰的八個粒子存在相互作用的關(guān)系。只要有四個相鄰的粒子已經(jīng)處理過，那就沒必要把八個相鄰的粒子都進行搜索。這樣一來就把搜索算法的時間復雜度從O(N2)降到O(NlogN)的級別。當光滑長度為常量時，使用鏈表能取得最佳效果。本文使用文獻[6]提出的動態(tài)邊界條件。邊界粒子與流體粒子具有相同的行為，遵循連續(xù)性、動量守恒和狀態(tài)方程。而不同的是，流體粒子有位移行為，而邊界粒子的位移為零。

當有粒子到達邊界時，邊界的粒子密度增加。因此，由于運動方程中的壓力作用，流體粒子的受力增加。當邊界粒子和流體粒子的距離變小，受排斥作用，流體粒子的密度、壓力和作用力增加。

5 程序運行結(jié)果

根據(jù)以上各節(jié)介紹的方程及數(shù)值求解的過程，利用OpenGL三維圖形庫[10]，在MS Visual Studio 2010環(huán)境下，使用C++面向?qū)ο蟮木幊谭椒▽崿F(xiàn)了對廣闊海面上劃浪航行的浪涌效果模擬。此程序算法對實現(xiàn)海上航行模擬海浪沖擊船舷及輪船破浪行駛都有相當大的價值?？紤]到計算量，SPH程序使用粒子數(shù)總量為12 000個，即為12k。粒子質(zhì)量為0.000 205 43，密度為400，光滑核長度為0.012，粒子半徑為0.04，粒子影響間距為0.005 9，流體黏度系數(shù)為0.8，時間步長度為0.004程序運行效果截圖如圖1所示。

圖1 程序運行截圖

6 結(jié)論

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

[1] MONAGHAN J J. Gravity currents and solitary waves[J]. Physica D: Nonlinear Phenomena, 1996,98(2-4):523-533.

[2] LO E Y M, Shao Songdong. Simulation of near-shore solitary wave mechanics by an incompressible SPH method[J]. Applied Ocean Research, 2002,24(5): 275-286.

[3] COLAGROSSI A, LANDRINI M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics[J]. Journal of Computational Physics, 2003,191(2): 448-475.

[4] 龔凱. 基于光滑質(zhì)點水動力學(SPH)方法的自由表面流動數(shù)值模擬研究[D].上海:上海交通大學，2009.

[5] 鄭坤. 基于SPH方法的波浪對水平板沖擊作用研究[D]. 大連:大連理工大學，2009.

[6] DALRYMPLE R A, ROGERS B D. Numerical modeling of water waves with the SPH method[J]. Coastal Engineering, 2006,53(2-3): 141-147.

[7] ANTOCI C, GALLATI M, SIBILLA S. Numerical simulation of fluid-structure interaction by SPH[J]. Computers and Structures, 2007,85(11-14): 879-890.

[8] ATAIE-ASHTIANI B, SHOBEYRI G, FARHADI L. Modified incompressible SPH method for simulating free surface problems[J]. Fluid Dynamics Research, 2008,40(9):637-661.

[9] BELYTSCHKO T,KRONGAUZ Y,DOLBOW J,et al.On the completness of meshfree particle methods[J]. International Journal for Numerical Methods in Engineering,1998, 43(5): 785-819.

[10] 張立成，張鴿.一種OpenGL局部縮放算法及應(yīng)用[J].微型機與應(yīng)用，2013，32(19)：44-47.

Simulation of seawater motion based on a 3D modified SPH method

Ou Xunyong, Chen Meiyi, Bao Haiqin, Kang Xiaoping

(School of Marine Information Engineering, Hainan Tropical Ocean University, Sanya 572022, China)

In this paper, around the standard control equation of SPH, the so-called standard linear reproducing kernel approximation method is used to modify the standard kernel function of SPH. And the gradient effect obtained in this way is also satisfactory. The model of the modified SPH method is used to simulate the fluid motion. The simulation results show the formation and propagation of the waves on the horizontal plane. Finally in the OpenGL environment it achieves the realistic simulation.

SPH; free-surface fluid; fluid motion

海南省自然科學基金項目(20166226)

TP311.1

A

10.19358/j.issn.1674- 7720.2017.12.021

歐訓勇，陳美伊，鮑海琴，等.三維修正SPH方法實現(xiàn)的海水運動效果模擬[J].微型機與應(yīng)用，2017,36(12)：71-74.

2016-12-20)

歐訓勇(1976-)，男，碩士，副教授，主要研究方向：OpenGL三維圖形技術(shù)、虛擬現(xiàn)實技術(shù)。

陳美伊(1980-)，女，碩士，實驗師，主要研究方向：數(shù)字媒體及應(yīng)用。

鮑海琴(1982-)，通信作者，女，碩士，實驗師，主要研究方向：環(huán)境工程。E-mail：573835681@qq.com。