三維修正SPH方法實(shí)現(xiàn)的海水運(yùn)動(dòng)效果模擬*

歐訓(xùn)勇，陳美伊，鮑海琴，康小平

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

?

三維修正SPH方法實(shí)現(xiàn)的海水運(yùn)動(dòng)效果模擬*

歐訓(xùn)勇，陳美伊，鮑海琴，康小平

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

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

SPH；自由表面流體；流體運(yùn)動(dòng)

0 引言

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

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

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

1 控制方程

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

(1)

(2)

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

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

(3)

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

(4)

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

2 SPH理論和方程

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

(5)

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

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

(6)

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

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

(7)

2.1 SPH標(biāo)準(zhǔn)方程

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

(8)

類(lèi)似地，粒子a的矢量場(chǎng)函數(shù)的散度計(jì)算公式為：

(9)

(10)

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

(11)

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

(12)

(13)

2.2 修正的SPH方法

(14)

(15)

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

AcT=I

(16)

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

(17)

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

3 SPH方法的時(shí)間步積分

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

(18)

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

(19)

(20)

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

(21)

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

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

(22)

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

(23)

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

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

當(dāng)有粒子到達(dá)邊界時(shí)，邊界的粒子密度增加。因此，由于運(yùn)動(dòng)方程中的壓力作用，流體粒子的受力增加。當(dāng)邊界粒子和流體粒子的距離變小，受排斥作用，流體粒子的密度、壓力和作用力增加。

5 程序運(yùn)行結(jié)果

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

圖1 程序運(yùn)行截圖

6 結(jié)論

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

[5] 鄭坤. 基于SPH方法的波浪對(duì)水平板沖擊作用研究[D]. 大連:大連理工大學(xué)，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].微型機(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

海南省自然科學(xué)基金項(xiàng)目(20166226)

TP311.1

A

10.19358/j.issn.1674- 7720.2017.12.021

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

2016-12-20)

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

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

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