黃曉婷 孫鵬楠 ,2) 呂鴻冠 鐘詩蘊(yùn)

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

? (南方海洋科學(xué)與工程廣東省試驗(yàn)室(珠海),廣東珠海 519082)

引言

隨著海洋開發(fā)逐漸走向深遠(yuǎn)海,海洋環(huán)境也變得越發(fā)惡劣和復(fù)雜,準(zhǔn)確預(yù)報(bào)結(jié)構(gòu)物在極端海浪環(huán)境下的水動(dòng)力性能是確保海洋工程裝備安全性和可靠性的重要前提.由于試驗(yàn)手段存在成本高、試驗(yàn)周期長等缺點(diǎn),數(shù)值模擬方法在波浪-結(jié)構(gòu)相互作用、海岸水動(dòng)力學(xué)等[1]應(yīng)用上越發(fā)廣泛.其中,建立高精度、高效率的數(shù)值水池是預(yù)報(bào)結(jié)構(gòu)物水動(dòng)力性能的重要手段之一.國內(nèi)外學(xué)者采用了許多數(shù)值方法建立數(shù)值水池,開展波浪傳播特性的研究.基于勢流理論,Grilli 等[2]、Sung[3]、Baudic 等[4]采用高階BEM 方法(拉格朗日網(wǎng)格法)對(duì)控制方程進(jìn)行離散,構(gòu)建數(shù)值水池.但無黏無旋的理想流體假設(shè)和波浪翻卷、破碎時(shí)的網(wǎng)格畸變使得強(qiáng)非線性波浪的動(dòng)態(tài)演變過程較難模擬.對(duì)于歐拉網(wǎng)格法,如有限體積法(FVM) 等,文獻(xiàn)[5] 采用開源程序庫OpenFOAM構(gòu)建了三維數(shù)值水池,生成的波浪精度較高,但在遠(yuǎn)距離傳播時(shí),受網(wǎng)格尺寸的限制,易造成數(shù)值耗散問題,且為了模擬自由液面的演化過程,常常需要結(jié)合復(fù)雜的VOF[6-7]等自由液面捕捉方法.

相比于網(wǎng)格化算法,近年來光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics,SPH)方法在船海領(lǐng)域得到廣泛應(yīng)用[8-10],尤其在波浪與結(jié)構(gòu)物相互作用研究上優(yōu)勢突出[11-12].SPH 方法的無網(wǎng)格特性使其易于處理大變形問題,完全避免了網(wǎng)格畸變;憑借其拉格朗日粒子特性,能輕松實(shí)現(xiàn)粒子運(yùn)動(dòng)路徑的追蹤,有利于處理自由液面,方便地實(shí)現(xiàn)波浪翻卷、破碎等現(xiàn)象的模擬.但SPH 模擬中存在一定的數(shù)值耗散[13-15],常常引起能量的非物理性衰減.例如,在波浪遠(yuǎn)距離傳播模擬過程中,能量非物理性耗散表現(xiàn)為波高的非物理性減小.如果波浪衰減不能得到有效補(bǔ)償或改善,將會(huì)限制SPH 數(shù)值水池應(yīng)用的物理尺度,無法實(shí)現(xiàn)遠(yuǎn)距離、長時(shí)間海浪環(huán)境演化的模擬.

目前,提高SPH 方法的能量守恒性,減少數(shù)值耗散,成為SPH 數(shù)值水池開發(fā)中廣受關(guān)注的問題.研究人員主要從光滑長度及核函數(shù)選擇上著手解決該問題.Colagrossi 等[16]的研究表明,選擇較大的光滑長度系數(shù)可以減少模擬重力波時(shí)的數(shù)值耗散量,但此舉容易造成較高的計(jì)算成本和較長的模擬時(shí)間.Zhang 等[17]采用不同核函數(shù)進(jìn)行規(guī)則波傳播的模擬,發(fā)現(xiàn)在SPH 模擬中,Quintic 核函數(shù)的動(dòng)能衰減率更小.Macià等[18]的研究表明采用Wendland C2 核函數(shù)進(jìn)行模擬可以有效減少動(dòng)能衰減.Meng等[19]提出了一種TENO-SPH 模型與算法,提升了SPH 的數(shù)值精度.此外,加密粒子也能起到減少耗散的效果.

對(duì)于具有自由液面的流體運(yùn)動(dòng)模擬,自由液面附近粒子的支持域被截?cái)?導(dǎo)致物理量及其導(dǎo)數(shù)的計(jì)算出現(xiàn)較大誤差.即使是處于流體內(nèi)部的粒子,當(dāng)其分布不均勻時(shí),往往也無法保證理論上的二階精度.這些誤差是傳統(tǒng)SPH 格式數(shù)值耗散的主要來源之一[20].因此,在SPH 算法中,施加核函數(shù)修正技術(shù)[21]也是一種減少能量耗散的重要手段.Wen 等[22]指出,采用核函數(shù)梯度修正(kernel gradient correction,KGC)可以顯著降低波浪傳播的衰減,而無需增大光滑長度,但缺點(diǎn)是粒子間的相互作用力不再具有對(duì)稱性,從而無法保證動(dòng)量守恒.

受文獻(xiàn)[23]啟發(fā),采用對(duì)稱化的修正核函數(shù)導(dǎo)數(shù)格式,對(duì)動(dòng)量方程進(jìn)行離散,有利于減少數(shù)值耗散,但文獻(xiàn)[23]的算法需要搜索自由液面,這將帶來額外的計(jì)算負(fù)擔(dān).為此,本文基于δ-SPH 算法,開發(fā)一種新的SPH 數(shù)值波浪水池,采用另一種核函數(shù)梯度修正方式,減少SPH 數(shù)值波浪水池中能量非物理性耗散的同時(shí),確保SPH 算法格式簡單、計(jì)算穩(wěn)定、動(dòng)量守恒.

本文具體安排如下:數(shù)值方法中,首先介紹δ-SPH 模型的控制方程,隨后論述本文修正算法的推導(dǎo)過程.數(shù)值結(jié)果中,首先通過振蕩液滴基準(zhǔn)算例,驗(yàn)證和討論本修正算法對(duì)非物理性能量耗散現(xiàn)象的改善效果.隨后,在數(shù)值波浪水池的構(gòu)建上,對(duì)規(guī)則波與不規(guī)則波傳播進(jìn)行模擬.先通過數(shù)值波高與試驗(yàn)波高的對(duì)比,論證傳統(tǒng)δ-SPH 數(shù)值水池中的波浪衰減現(xiàn)象;隨后,通過傳統(tǒng)算法結(jié)果與本文修正算法結(jié)果的對(duì)比,說明本文修正算法在低能量耗散數(shù)值波浪水池開發(fā)上的有效性.

1 數(shù)值算法

1.1 δ-SPH 模型

SPH 方法是一種具有拉格朗日特性的粒子法.在計(jì)算時(shí),流場被離散成粒子,通過賦予粒子物理屬性,追蹤粒子的拉格朗日運(yùn)動(dòng)來實(shí)現(xiàn)流場物理演化過程的模擬.在δ-SPH 計(jì)算模型中,粒子密度、速度和位置變化的控制方程[24]為

式中,n=2,光滑長度h=αs?x,αs為光滑長度系數(shù),?x為粒子間距.α 為黏性系數(shù),在SPH 模擬中,通過在消波區(qū)內(nèi)增大黏性系數(shù)實(shí)現(xiàn)數(shù)值消波,具體設(shè)置為

式中,L為數(shù)值水池總長度,L?Ldamp為消波區(qū)長度.在本文模擬中,采用兩種常用的光滑長度系數(shù)αs=1.3,2.0 進(jìn)行計(jì)算.其中 αs=1.3 時(shí),由于核函數(shù)緊支域內(nèi)粒子數(shù)量更少,因此計(jì)算量小,但是離散誤差相比 αs=2.0 時(shí)更大,數(shù)值耗散也會(huì)更大.本文將開發(fā)改進(jìn)SPH 算法,力爭在αs=1.3 時(shí),也獲得較好的計(jì)算精度和較小的數(shù)值耗散.

1.2 低能量耗散的修正δ-SPHC 算法

在傳統(tǒng)SPH 方法中,梯度算子粒子近似[25]可寫為

然而,在粒子分布不均勻或核函數(shù)截?cái)嗲闆r下,該離散形式不能較好地確保計(jì)算精度.對(duì)此,許多修正方法被提出,其中常用的有修正光滑粒子法(corrective smoothed particle method,CSPM)[26-27].梯度算子采用CSPM 修正后為

經(jīng)Wen 等[22]和Zago 等[23]研究啟發(fā),采用上述CSPM 算法中修正的核函數(shù)梯度,有助于改善SPH 計(jì)算中能量的衰減問題,但是為了保證動(dòng)量守恒,也就是確保粒子對(duì)(i?j)之間相互作用力的對(duì)稱性,需要對(duì)CSPM 格式進(jìn)行改進(jìn).對(duì)此,Zago 等[23]提出了一種改進(jìn)后的修正格式,但是需要顯式搜索自由液面,這帶來了額外的計(jì)算負(fù)擔(dān),而本文基于CSPM 修正形式,參考Sun 等[28]處理梯度算子的方法,采用另一種具有對(duì)稱性的梯度算子離散形式,推導(dǎo)過程見下:

首先,根據(jù)核函數(shù)導(dǎo)數(shù)的特性[25]有

本文修正算法記為δ-SPHC.對(duì)比Zago 等[23]的算法,本修正算法避免了自由液面搜索的步驟;另外,相對(duì)于δ-SPH,δ-SPHC修正算法并不引入新的核函數(shù)修正矩陣,只是再次調(diào)用密度耗散項(xiàng)中的矩陣Mi對(duì)動(dòng)量方程中壓力梯度算子進(jìn)行修正,因此不會(huì)造成計(jì)算成本和時(shí)間的顯著增加.下文將采用振蕩液滴算例進(jìn)行修正算法抗能量非物理性衰減效果的測試和討論,并應(yīng)用于數(shù)值波浪水池規(guī)則波與不規(guī)則波的數(shù)值模擬和試驗(yàn)驗(yàn)證.本文的數(shù)值積分采用四階龍格庫塔法,具體可參考文獻(xiàn)[31].

2 振蕩液滴算例驗(yàn)證

2.1 收斂性分析與精度驗(yàn)證

振蕩液滴算例無需施加固壁邊界,初始條件簡單,被廣泛應(yīng)用于SPH 等數(shù)值新算法的驗(yàn)證[32].Antuono 等[33]采用振蕩液滴算例驗(yàn)證了δ-SPH 模型中的能量守恒問題.基于此,本文采用該算例對(duì)修正算法δ-SPHC減少能量非物理性衰減的效果進(jìn)行驗(yàn)證.

該算例中,初始速度V(u,v)、初始半徑為R的圓形液滴在向心體積力場F中振蕩,其中F(x,y)=?B2r(x,y),B=1,R=1 .圖1 展示了初始時(shí)刻t0時(shí)液滴的速度場分布[33],可表示為

圖1 振蕩液滴初始速度分布Fig.1 Initial velocity distribution

A(t0) 為控制液滴初始速度場的參數(shù).在本文模擬中,A(t0)=1 .初始?jí)毫龇植既鐖D2 所示,壓力場[33]表示為

圖2 振蕩液滴初始?jí)毫Ψ植糉ig.2 Initial pressure distribution

首先,為檢驗(yàn)SPH 修正算法的收斂性,分別在粒子分辨率為R/?x=50,R/?x=100,R/?x=200 時(shí),開展振蕩液滴運(yùn)動(dòng)特性研究.由于振蕩液滴質(zhì)心與向心體積力F的匯聚點(diǎn)重合,液滴不受其他外力作用,液滴在原位保持振蕩.因此,理論上振蕩液滴的角動(dòng)量、線動(dòng)量都守恒,即

圖3 展示在光滑長度系數(shù) αs=1.3 和 αs=2.0 條件下,不同粒子分辨率時(shí),角動(dòng)量MA隨時(shí)間的變化情況.可見,在三種粒子分辨率下,MA的值圍繞理論值MA=0 波動(dòng)的量級(jí)均在10?5以內(nèi),且隨著粒子分辨率的增加,角動(dòng)量MA的誤差范圍逐漸變小.在粒子分辨率R/?x=200 時(shí),MA趨 近于零.即在R/?x=200時(shí),角動(dòng)量誤差逐漸收斂于0,說明該修正算法在提高粒子分辨率時(shí)能實(shí)現(xiàn)良好的角動(dòng)量守恒性.

圖3 不同粒子分辨率下角動(dòng)量 MA 的時(shí)歷曲線Fig.3 Time evolution of the angular momentum MA with different particle resolutions

圖4 展示了光滑長度系數(shù) αs=1.3 和 αs=2.0 條件下,R/?x=200 分辨率時(shí),線動(dòng)量ML隨時(shí)間的變化情況.觀察可得,線動(dòng)量誤差在?10?14～ 10?14范圍內(nèi)波動(dòng),可以認(rèn)為本修正算法的線動(dòng)量實(shí)現(xiàn)了精確守恒.

圖4 線動(dòng)量時(shí)歷曲線Fig.4 Time evolution of the linear momentum

同時(shí),定義液滴振蕩過程的機(jī)械能衰減率為Eper(t)=(E(t)?Eint)/Eint,其中Eint為初始機(jī)械能,E(t) 為t時(shí)刻的機(jī)械能.液滴振蕩過程的機(jī)械能衰減率隨時(shí)間變化結(jié)果如圖5 所示.可以發(fā)現(xiàn),隨著粒子分辨率的提高,機(jī)械能衰減率Eper越來越小.在R/?x=100,R/?x=200 時(shí),兩者機(jī)械能衰減率基本接近,計(jì)算時(shí)長達(dá)t=120 時(shí),機(jī)械能衰減控制在5%以內(nèi).

圖5 αs=2.0,不同粒子分辨率時(shí)機(jī)械能衰減率 Eper 時(shí)歷曲線Fig.5 Time evolution of the decay rate of mechanical energy at different particle resolutions with αs=2.0

圖6 展示了R/?x=200 時(shí),振蕩液滴在不同時(shí)刻的理論形態(tài)與修正算法δ-SPHC模擬的形態(tài).可以發(fā)現(xiàn),在兩種光滑長度條件下,δ-SPHC模擬的振蕩形態(tài)均與理論值吻合良好,且在t=8.75T周期后,液滴仍保持較為光滑的形態(tài).由此可見,本修正算法具有較高的精度.

圖6 αs=1.3(上)和αs=2.0(下)時(shí),振蕩液滴形態(tài):理論值(虛線)和SPH 結(jié)果(壓力云圖)Fig.6 With smoothing length coefficient αs=1.3 (top) and αs=2.0 (bottom),the shapes of oscillating droplet:theoretical results (dash line) and SPH results (pressure contour)

綜上可見,修正SPH 算法在粒子分辨率R/?x=200時(shí),具有較好的數(shù)值精度和動(dòng)量與能量守恒性.因此,下文將在粒子分辨率R/?x=200 時(shí),依托該算例進(jìn)一步討論不同算法抗能量衰減效果的差異.

2.2 傳統(tǒng)δ-SPH 與修正δ-SPHC 算法的對(duì)比研究

基于上節(jié)收斂性分析與精度驗(yàn)證,本節(jié)在粒子分辨率R/?x=200 條件下,開展傳統(tǒng)δ-SPH 模型與改進(jìn)δ-SPHC模型抗能量衰減效果的對(duì)比研究.圖7展示了δ-SPH 與δ-SPHC模擬的振蕩液滴動(dòng)能隨時(shí)間的變化曲線.可見,在傳統(tǒng)算法δ-SPH 結(jié)果中,αs=1.3 時(shí),動(dòng)能隨著時(shí)間的衰減十分明顯.在αs=2.0時(shí),衰減有所減緩.但是,取h=2?x相比h=1.3?x,由于相鄰粒子數(shù)量的增多,會(huì)導(dǎo)致額外的計(jì)算成本和時(shí)間,且從圖7 的局部放大圖可以看出,在較長時(shí)間的計(jì)算下,仍會(huì)出現(xiàn)輕微動(dòng)能衰減.對(duì)比修正算法δ-SPHC,在光滑長度系數(shù)αs=1.3 時(shí),動(dòng)能幅值基本保持不變,已體現(xiàn)出較好的抗能量衰減效果;在光滑長度系數(shù)αs=2.0 時(shí),本修正算法實(shí)現(xiàn)了良好的能量守恒.此外,由局部放大圖可以看出,在修正算法αs=2.0 時(shí),動(dòng)能幅值與初始幅值始終保持一致,取得了較為精確的計(jì)算結(jié)果.而在δ-SPHC結(jié)果中,αs=1.3時(shí),動(dòng)能第一個(gè)峰值較初始值高,主要原因在第一個(gè)周期的計(jì)算中,粒子會(huì)從格子分布狀態(tài)(圖8(a))突然過渡到一種各項(xiàng)同性的均勻分布狀態(tài)(圖8(b)).在粒子分布變化過程中,存在一些粒子的瞬時(shí)分布十分不均勻,此時(shí),這不均勻性會(huì)造成.也就是說,部分粒子的運(yùn)動(dòng)并不是由于壓力梯度造成,而是粒子分布的不均勻性導(dǎo)致的數(shù)值加速度產(chǎn)生,這對(duì)粒子的總動(dòng)能造成了一定的非物理性影響.但是,從圖7 下方的放大圖可以看出,這種由于粒子分布突變引起的動(dòng)能誤差控制在1%以內(nèi).類似現(xiàn)象在Colagrossi 等[34]的研究結(jié)果也有過報(bào)道,可通過改善初始粒子分布得以解決.

圖7 δ-SPH 與δ-SPHC 動(dòng)能比較:總體圖(上)與局部放大圖(下),其中黑矩形框?yàn)榉糯髤^(qū)域Fig.7 The comparison between time evolutions of kinetic energy between δ-SPH (top) and δ-SPHC (bottom),an enlarged zone view of the results in the dark rectangle is shown in the bottom panel

圖8 振蕩液滴粒子分布Fig.8 Particle distribution of oscillating droplet

為進(jìn)一步突顯修正算法的低能量耗散特性,圖9呈現(xiàn)了兩種算法中,機(jī)械能衰減率的時(shí)歷曲線.δ-SPHC算法在兩種光滑長度系數(shù)下的計(jì)算結(jié)果中,機(jī)械能衰減率均較小,基本接近零,說明本文修正算法能夠有效減少機(jī)械能的非物理性耗散.其中,在光滑長度系數(shù) αs=2.0 時(shí)效果更為顯著.值得一提的是,δ-SPHC算法在αs=2.0 時(shí)的機(jī)械能衰減率相比δ-SPH 算法在αs=1.3 時(shí)的機(jī)械能衰減率更小,表明本修正SPH 算法的一大優(yōu)勢是能在較小光滑長度系數(shù)下實(shí)現(xiàn)較為準(zhǔn)確的模擬,從而顯著減少計(jì)算時(shí)間.

圖9 不同算法和光滑長度系數(shù)下,振蕩液滴機(jī)械能衰減率時(shí)歷曲線Fig.9 Time evolution of the decay rate of mechanical energy for the oscillating droplet with different SPH models and different smoothing length coefficient

以上算例在配有Intel(R) Core(TM)I9-10900 K CPU 和48 GB 內(nèi)存的個(gè)人電腦上進(jìn)行計(jì)算.

表1 比較了在最高粒子分辨率R/?x=200,αs=2.0 參數(shù)下,δ-SPH 與δ-SPHC模型模擬振蕩液滴的計(jì)算時(shí)長.可見,兩種算法單步計(jì)算時(shí)間相當(dāng).具體而言,雖然修正算法δ-SPHC對(duì)壓力梯度項(xiàng)進(jìn)行了核函數(shù)修正,但相對(duì)δ-SPH 算法,計(jì)算時(shí)間僅增加1.09 倍,計(jì)算量增加并不顯著.在下節(jié)中,該算法將被應(yīng)用于數(shù)值波浪水池,驗(yàn)證其改善SPH 模擬遠(yuǎn)距離波浪傳播時(shí)的波高抗衰減效果.

表1 δ-SPH 算法與δ-SPHC 算法計(jì)算效率比較Table 1 Comparison of computational efficiency between δ-SPH and δ-SPHC schemes

3 SPH 數(shù)值造波與試驗(yàn)驗(yàn)證

3.1 試驗(yàn)與數(shù)值水池設(shè)置

3.1.1 試驗(yàn)水池

作者在中山大學(xué)海洋工程與技術(shù)學(xué)院波浪水槽中開展造波試驗(yàn).如圖10 所示,該水池長度L=16 m,水深d=0.266 m,采用推板造波方式生成規(guī)則波和不規(guī)則波,推板運(yùn)動(dòng)曲線如圖11 所示.試驗(yàn)時(shí),波高儀設(shè)置在距離造波板平衡位置xprobe=2.37 m,4.37 m,5.37 m,6.37 m,15 m 處.

圖10 中山大學(xué)波浪試驗(yàn)水槽Fig.10 Experimental wave tank in Sun Yat-sen University

圖11 造波推板橫向位移隨時(shí)間變化:(a)周期為T =1 s 的規(guī)則波,(b)不規(guī)則波Fig.11 Paddle motions of the regular wave with (a) period T =1 s and(b) irregular wave

3.1.2 SPH 數(shù)值水池設(shè)置

在SPH 模擬中,固定虛粒子邊界法[35,36]能實(shí)現(xiàn)固壁邊界的精確處理,有效防止壁面穿透,且保證壁面處壓力連續(xù),因而得到了廣泛應(yīng)用.

本文SPH 數(shù)值模擬中,數(shù)值水池長度L=60 m,邊界條件采用固定虛粒子進(jìn)行施加.如圖12 所示,固定虛粒子厚度不小于核函數(shù)半徑,其壓力從流體粒子外插得到,更多細(xì)節(jié)可參考文獻(xiàn)[37].本文SPH模擬中,通過強(qiáng)制流體粒子與壁面虛粒子之間黏性力為0,實(shí)現(xiàn)壁面邊界處的自由滑移邊界條件.數(shù)值水池左側(cè)的運(yùn)動(dòng)造波板也由虛粒子組成,通過控制造波板虛粒子的總體位置發(fā)生改變,實(shí)現(xiàn)推板造波.具體實(shí)現(xiàn)方式如下.

圖12 SPH 數(shù)值水池示意圖,紅豎線表示波高儀,選取距離波高儀一個(gè)粒子間距內(nèi)自由面粒子的最大高度為波面高度Fig.12 Schematic diagram of the numerical wave tank,the red vertical line represents wave gage.Maximum height of the free-surface particle within one-particle distance from the red line is measured to represent the wave elevation

強(qiáng)制SPH 造波板按照?qǐng)D11 中試驗(yàn)所輸出的造波板位移曲線進(jìn)行運(yùn)動(dòng),從而確保SPH 模擬中,造波板運(yùn)動(dòng)路徑與試驗(yàn)一致.同時(shí),為克服SPH 模擬的時(shí)間步長與試驗(yàn)中輸出的造波板運(yùn)動(dòng)數(shù)據(jù)時(shí)間間隔不一致,采用線性插值技術(shù),根據(jù)物理試驗(yàn)中造波板運(yùn)動(dòng)特性,更新SPH 模擬中的造波板位移xS PH、速度vS PH及加速度aS PH,即

其中,上標(biāo)SPH表示數(shù)值水池模擬中的物理量,上標(biāo)k表示試驗(yàn)中造波板運(yùn)動(dòng)輸出的時(shí)間步.試驗(yàn)中造波板的速度vk與加速度ak均采用線性差分計(jì)算獲得,即

此外,值得注意的是,當(dāng)波浪傳播到數(shù)值水池另一端時(shí),往往會(huì)產(chǎn)生波浪反射,降低造波精度,這就需要在數(shù)值水池右端構(gòu)建數(shù)值消波區(qū).類似于網(wǎng)格算法中增大數(shù)值黏性來實(shí)現(xiàn)阻尼消波,本文通過提高SPH 數(shù)值水池消波區(qū)的黏性系數(shù),使波浪進(jìn)入消波區(qū)時(shí)能夠在短時(shí)間內(nèi)衰減.具體設(shè)置為 :在數(shù)值水池設(shè)置消波區(qū)Ldamp,黏性系數(shù) α 線性增大,詳見公式(4).本文計(jì)算中,Ldamp=20 m.

3.2 δ-SPH 算法在高粒子分辨率下的精度驗(yàn)證

Colagrossi 等[16]指出δ-SPH 模型能較為準(zhǔn)確地重現(xiàn)在物理黏性條件下的波浪傳播和衰減過程.因此,本文基于δ-SPH 模型開發(fā)低能量耗散的數(shù)值水池.首先,在αs=2.0 條件下,驗(yàn)證δ-SPH 計(jì)算結(jié)果的精度和收斂性.選取d/?x=15,30,60,即水深方向上粒子數(shù)分別為15,30,60 進(jìn)行規(guī)則波模擬,在三種粒子分辨率條件下監(jiān)測xprobe=2.37 處的波面高度隨時(shí)間變化,如圖13 所示.可以看出,數(shù)值波高結(jié)果變化趨勢均與試驗(yàn)波高結(jié)果吻合良好,在d/?x=60 時(shí),數(shù)值波高收斂于試驗(yàn)波高.基于以上討論,在d/?x=60 條件下,δ-SPH 模型能較為精確實(shí)現(xiàn)造波及波浪傳播的模擬.但是值得一提的是,δ-SPH 模型在αs<2.0 時(shí),容易出現(xiàn)波高的非物理性衰減,具體論述見下一節(jié).

圖13 αs=2.0 時(shí),不同粒子分辨率下波面高度的δ-SPH 模擬值與試驗(yàn)值對(duì)比Fig.13 Comparison of the wave elevations between δ-SPH results and experimental data,SPH results with different particle resolutions and αs=2.0 are provided

3.3 δ-SPH 算法中的波浪衰減現(xiàn)象

本節(jié)采用δ-SPH 模型,在光滑長度系數(shù) αs=1.3 和 αs=2.0 條件下,進(jìn)行規(guī)則波造波模擬.為了測定波高變化曲線,沿水池長度方向,距離造波板初始位置xprobe=2.37 m,4.37 m,5.37 m,6.37 m,15 m 處設(shè)置5 個(gè)數(shù)值波高儀.

在αs=1.3 時(shí),圖14 展示了計(jì)算時(shí)間為t=27 s 的波面形態(tài),可以較為明顯地看出波浪在x=10 m 處開始出現(xiàn)衰減.選取xprobe=2.37 m,5.37 m,15 m 處的波面高度時(shí)歷曲線進(jìn)行比較,如圖15 所示.在xprobe=2.37 m 處,波高與試驗(yàn)相差較大,且隨著傳播距離的增大,波高逐漸衰減,在xprobe=15 m 處,波高衰減嚴(yán)重.可以看出,在傳統(tǒng)δ-SPH 中,光滑長度系數(shù) αs=1.3 條件下無法準(zhǔn)確模擬波浪傳播.圖16 展示了δ-SPH 算法在光滑長度系數(shù) αs=2.0 時(shí)計(jì)算的波面形態(tài).與振蕩液滴類似,相比光滑長度系數(shù) αs=1.3 時(shí),衰減速度有所減緩,與Colagrossi 等[16]的結(jié)論相符.這表明,增大 αs可以有效減緩SPH 的數(shù)值衰減,但在大尺度、長時(shí)間、遠(yuǎn)距離的海浪模擬下,光滑長度系數(shù) αs取2.0 時(shí),將帶來巨大的計(jì)算量,不利于工程應(yīng)用.且比較不同距離處波高,如圖17所示,可以發(fā)現(xiàn)在xprobe=15 m 時(shí),仍存在明顯波高衰減.因此,改善SPH 方法在數(shù)值波浪模擬時(shí)能量的非物理性耗散,在較低光滑長度系數(shù)下實(shí)現(xiàn)數(shù)值水池造波的高精度高效率模擬具有重要意義.

圖14 在αs=1.3,d/?x=60 條件下,傳統(tǒng)δ-SPH 模擬的波面形態(tài)Fig.14 Wave surface simulated by δ-SPH method with αs=1.3 and d/?x=60

圖15 在αs=1.3 ,d/?x=60 條件下,傳統(tǒng)δ-SPH 模型計(jì)算的不同位置波面高度時(shí)歷曲線Fig.15 Wave elevation probed at different positions simulated by δ-SPH method with αs=1.3 and d/?x=60

圖16 在αs=2.0 ,d/?x=60 條件下,傳統(tǒng)δ-SPH 模擬的波面形態(tài)Fig.16 Wave surface simulated by δ-SPH method with αs=2.0 and d/?x=60

圖17 在αs=2.0,d/?x=60 條件下,傳統(tǒng)δ-SPH 模擬的不同位置波面高度時(shí)歷曲線Fig.17 Wave elevation probed at different positions simulated by δ-SPH method with αs=2.0 and d/?x=60

4 修正算法δ-SPHC 的波浪傳播模擬

4.1 規(guī)則波模擬

4.1.1 收斂性分析與精度驗(yàn)證

首先,基于規(guī)則波傳播算例,校核本文修正的δ-SPHC算法精度與收斂性.如圖18 所示,在αs=1.3時(shí),選取xprobe=2.37 m 處數(shù)值波高結(jié)果與試驗(yàn)結(jié)果進(jìn)行比較,可以看出,水深方向不同數(shù)量粒子的計(jì)算結(jié)果略有差別,且隨著d/?x的增大,數(shù)值波高趨近于試驗(yàn)波高.在d/?x=60 時(shí),SPH 結(jié)果與試驗(yàn)結(jié)果吻合良好.與振蕩液滴算例類似,圖18 方框中的峰值可能是由于粒子由初始的格子狀分布到各向同性均勻分布時(shí)引起的能量瞬時(shí)誤差.試驗(yàn)波高曲線與SPH 數(shù)值結(jié)果曲線在初始階段第一個(gè)峰值存在相位差,可能是由于試驗(yàn)波高采集與造波板運(yùn)動(dòng)存在一定的時(shí)間差.在αs=2.0 時(shí),水深方向不同粒子數(shù)量的計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比如圖19 所示.同樣可以看到,在d/?x=60 時(shí),數(shù)值波高曲線與試驗(yàn)結(jié)果基本吻合.另外,由于物理試驗(yàn)水池長度僅有16 m,消波區(qū)長度設(shè)置有限,易導(dǎo)致試驗(yàn)后期消波不徹底,造成波浪周期的變化,即隨著時(shí)間增加,周期略小于1 s;而數(shù)值水池長度L=60 m,能較為準(zhǔn)確地重現(xiàn)規(guī)則波傳播,并且確保周期保持不變.因此,可以看到,在圖18、圖19 中,t=20 s 之后SPH 結(jié)果的波峰和波谷時(shí)刻與試驗(yàn)值存在略微誤差,但是波高和波谷的幅值與試驗(yàn)值吻合良好.基于以上討論可以得出,在水深方向粒子數(shù)d/?x=60 條件下,計(jì)算結(jié)果能夠達(dá)到預(yù)期精度要求.因此,下文將在d/?x=60 粒子設(shè)置下,進(jìn)行修正δ-SPHC算法的抗能量非物理性衰減效果的討論.

圖18 在αs=1.3 和 d/?x=15,30,60 條件下,δ-SPHC 計(jì)算的波面高度時(shí)歷曲線與試驗(yàn)值對(duì)比Fig.18 Comparison between δ-SPHC results and experimental data for wave elevation:SPH results with αs=1.3 and d/?x=15,30,60

圖19 在αs=2.0 和 d/?x=15,30,60 條件下,δ-SPHC 計(jì)算的波面高度時(shí)歷曲線與試驗(yàn)值對(duì)比Fig.19 Comparison between δ-SPHC results and experimental data for wave elevation:SPH results with αs=2.0 and d/?x=15,30,60

圖20 呈現(xiàn)了在xprobe=2.37 m 處,規(guī)則波波高數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的對(duì)比.采用傳統(tǒng)δ-SPH 計(jì)算時(shí),光滑長度系數(shù) αs=1.3 條件下,由于核函數(shù)緊支域內(nèi)粒子數(shù)量少,數(shù)值結(jié)果與試驗(yàn)結(jié)果相差較大,而光滑長度系數(shù) αs=2.0 時(shí),數(shù)值波高結(jié)果接近試驗(yàn)結(jié)果.但值得注意的是,采用δ-SPHC,光滑長度系數(shù)αs=1.3 時(shí),數(shù)值結(jié)果與試驗(yàn)結(jié)果能較好吻合,基本達(dá)到傳統(tǒng)δ-SPH 算法在光滑長度系數(shù) αs=2.0 時(shí)的計(jì)算精度,與上文振蕩液滴算例得出的結(jié)論一致,再次說明了本文修正δ-SPHC算法的優(yōu)勢,即在光滑長度系數(shù) αs=1.3,2.0 時(shí),計(jì)算結(jié)果與試驗(yàn)結(jié)果均能較好吻合.圖21 展示了數(shù)值水池的壓力場,可見流場壓力分布均勻,在運(yùn)動(dòng)較為劇烈的推板附近,也能確保壓力場穩(wěn)定.基于以上討論,可以認(rèn)為修正的δ-SPHC算法具有較高的精度.

圖20 xprobe=2.37 m 處,不同SPH 算法下的波面高度數(shù)值結(jié)果與試驗(yàn)值對(duì)比Fig.20 Comparison between different SPH simulations and experimental data for wave elevations at xprobe=2.37 m

圖21 規(guī)則波壓力云圖(δ-SPHC,αs=2.0,d/?x=60)Fig.21 Pressure contour of regular wave (δ-SPHC,αs=2.0,d/?x=60)

4.1.2 δ-SPHC的抗能量非物理性衰減結(jié)果分析

經(jīng)過前期收斂性與精度驗(yàn)證,本節(jié)在水深方向粒子數(shù)取d/?x=60 時(shí),開展修正算法δ-SPHC在兩種光滑長度系數(shù) αs=1.3,2.0 下的抗能量非物理性衰減效果討論.在αs=1.3 時(shí),如圖22 所示,觀察在修正算法模擬下的自由面形態(tài),可以看出在波浪傳播20 m 時(shí),波面仍與推板附近波面基本持平.監(jiān)測不同距離下的波面高度如圖23 所示,與試驗(yàn)結(jié)果對(duì)比可以看出,在xprobe=15 m 處,δ-SPHC的數(shù)值波高與推板附近波高基本接近,說明能量衰減現(xiàn)象得到極大改善.圖24 展示了在αs=2.0 時(shí),規(guī)則波傳播的波面形態(tài).可見,在遠(yuǎn)距離傳播時(shí),波面仍保持在同一高度,說明了修正算法抗能量衰減的有效性.對(duì)比圖23 和圖25,可見在xprobe=15 m 位置處,αs=2.0 比 αs=1.3 計(jì)算的波高更加精確和穩(wěn)定,說明光滑長度系數(shù)的增大有利于減少離散誤差,其中在xprobe=15 m 處出現(xiàn)略微振蕩,原因主要是由于數(shù)值誤差在遠(yuǎn)距離、長時(shí)間上的累積造成,但是相比傳統(tǒng)算法,修正算法的精度得到了有效改善.

圖22 αs=1.3,d/?x=60 時(shí),δ-SPHC 算法模擬的波面形態(tài)Fig.22 The wave surface simulated by δ-SPHC with αs=1.3 and d/?x=60

圖23 αs=1.3,d/?x=60 時(shí),δ-SPHC 模擬的不同位置波面高度時(shí)歷曲線Fig.23 Wave elevations probed at different positions simulated by δ-SPHC method with αs=1.3 and d/?x=60

圖24 αs=2.0,d/?x=60 時(shí),δ-SPHC 算法模擬的波面形態(tài)Fig.24 The wave surface simulated by δ-SPHC with αs=2.0 and

圖25 αs=2.0,d/?x=60 時(shí),δ-SPHC 模擬的不同位置波面高度時(shí)歷曲線Fig.25 Wave elevations probed at different positions simulated by δ-SPHC method with αs=2.0 and d/?x=60

如圖26 所示,綜合對(duì)比δ-SPH,δ-SPHC與試驗(yàn)方法在xprobe=6.37 m 的波高曲線,可見相比于傳統(tǒng)δ-SPH 模型,修正的δ-SPHC算法在αs=1.3 時(shí)的波高結(jié)果更加準(zhǔn)確.值得一提的是,其計(jì)算結(jié)果與傳統(tǒng)δ-SPH 模型在αs=2.0 時(shí)波高結(jié)果基本重合,并與試驗(yàn)波高吻合良好.此外,采用δ-SPHC修正算法,在光滑長度系數(shù) αs=2.0 時(shí),波高計(jì)算結(jié)果與試驗(yàn)波高高度吻合.

圖26 在xprobe=6.37 m 處,規(guī)則波波面高度數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果Fig.26 Time evolution of the regular wave elevation at xprobe=6.37 m:SPH results and experimental data

為進(jìn)一步定量驗(yàn)證δ-SPHC的模擬精度,表2 給出了試驗(yàn)和SPH 模擬中,xprobe=6.37 m 測點(diǎn)處的平均波高值及相對(duì)誤差,后者定義為 ε=|HSPH?HEXP|/HEXP,HEXP為試驗(yàn)測量的平均波高.可見,在兩種光滑長度系數(shù)下,δ-SPHC算法在較遠(yuǎn)處的波高監(jiān)測結(jié)果與試驗(yàn)相對(duì)誤差均較小.本節(jié)討論表明,修正算法能有效解決能量非物理性衰減問題,且在較低光滑長度系數(shù)時(shí)能得到更為準(zhǔn)確的結(jié)果.

4.1.3 不同人工黏性系數(shù)下δ-SPHC模擬結(jié)果分析

在數(shù)值水池中,能量衰減來源主要有兩部分:一是由于流體物理黏性引起的能量耗散;二是本文所聚焦的由于數(shù)值誤差引起的能量非物理性耗散.本文雖然未在SPH 控制方程中添加物理黏性項(xiàng),但是使用了人工黏性項(xiàng)πiv,這等效于增加了流體的黏性.為進(jìn)一步討論 δ -SPHC在不同人工黏性系數(shù) α 下,對(duì)長時(shí)間、遠(yuǎn)距離波浪傳播模擬的能量耗散情況,選取四個(gè)不同黏性系數(shù),即 α=0.01,0.02,0.05,0.1 開展數(shù)值模擬.圖27 給出了在xprobe=6.37 m 處監(jiān)測的波高曲線.可見,在α=0.01 和 α=0.02 時(shí),波浪的衰減較小,與實(shí)驗(yàn)值吻合較好.在α=0.05 時(shí),由于采用了更大的人工黏性,等效增加了物理黏性,因此波高衰減有所增加.在α=0.1 時(shí),由于人工黏性引起的波高衰減更為明顯,與實(shí)驗(yàn)值相比誤差較大.綜合以上分析,在α=0.02 時(shí),能夠獲得較為穩(wěn)定和精確的波高計(jì)算結(jié)果,且與試驗(yàn)值吻合良好.因此,本文其他算例均基于 α=0.02 進(jìn)行模擬.基于此,在下文中進(jìn)一步開展不規(guī)則波生成和傳播的SPH 模擬和驗(yàn)證.

圖27 δ-SPHC 在不同黏性系數(shù)α 時(shí),xprobe=6.37 m 處波面高度時(shí)歷曲線Fig.27 Time history of wave elevation at xprobe=6.37 m simulated by δ-SPHC with different α

4.2 基于δ-SPHC 的不規(guī)則波模擬

4.2.1 收斂性分析與精度驗(yàn)證

經(jīng)過上文分析,在光滑長度系數(shù) αs=1.3 時(shí),修正算法δ-SPHC能夠準(zhǔn)確模擬規(guī)則波傳播,而本節(jié)算法,在αs=1.3 時(shí),模擬和驗(yàn)證不規(guī)則波傳播.選取d/?x=15,30,60 進(jìn)行收斂性分析,如圖28 所示,在水深方向粒子數(shù)d/?x=15,30,60 時(shí),在xprobe=2.37處探測的波高與試驗(yàn)值基本一致.此外,在不規(guī)則波模擬中,δ-SPHC實(shí)現(xiàn)了穩(wěn)定的壓力場預(yù)報(bào)(見圖29),進(jìn)一步說明了δ-SPHC方法所構(gòu)建的數(shù)值水池的精度.由于d/?x=60 時(shí),數(shù)值結(jié)果已基本收斂于試驗(yàn)值,因此下文選取d/?x=60 進(jìn)行不規(guī)則波模擬結(jié)果的討論.

圖28 αs=1.3 和不同粒子分辨率條件下,不規(guī)則波δ-SPHC 模擬結(jié)果的收斂性分析Fig.28 Convergence analysis of irregular waves simulated by δ-SPHC with αs=1.3 at different particle resolutions

圖29 不規(guī)則波壓力云圖(δ-SPHC,αs=2.0,d/?x=60)Fig.29 Pressure contour of irregular wave (δ-SPHC,αs=2.0,d/?x=60)

4.2.2 δ-SPHC抗能量非物理性衰減結(jié)果分析

如圖30 所示,監(jiān)測xprobe=6.37 處波高時(shí)歷曲線,可見在傳統(tǒng)δ-SPH 計(jì)算結(jié)果中,與試驗(yàn)結(jié)果相比,波高衰減明顯;而修正算法結(jié)果與試驗(yàn)結(jié)果基本吻合,說明不規(guī)則波模擬中波高衰減現(xiàn)象也得到較好改善.在δ-SPHC計(jì)算框架下,數(shù)值水池生成不規(guī)則波的能量非物理性耗散問題得到解決,并能夠在光滑長度系數(shù) αs=1.3 時(shí),實(shí)現(xiàn)較高的造波精度.

圖30 xprobe=6.37 處,不規(guī)則波的波面高度時(shí)歷曲線Fig.30 Time evolution of the irregular wave elevation measured at xprobe=6.37 m

5 結(jié)論與展望

本文給出一種修正的δ-SPHC算法,采用對(duì)稱化的核函數(shù)導(dǎo)數(shù)修正格式,對(duì)壓力梯度進(jìn)行離散,確保粒子對(duì)之間作用力的對(duì)稱性,實(shí)現(xiàn)了精確線動(dòng)量守恒,且無需自由面粒子搜索等額外的復(fù)雜數(shù)值處理過程.通過模擬振蕩液滴算例,說明其抗能量衰減的有效性.隨后將其應(yīng)用于數(shù)值波浪水池.研究結(jié)果表明:

(1) δ-SPH 模型存在能量非物理性耗散問題,增大光滑長度系數(shù)可以有效減緩能量衰減;

(2)修正的δ-SPHC算法具有良好的動(dòng)量守恒性,克服了傳統(tǒng)核函數(shù)梯度修正的動(dòng)量不守恒問題;

(3)修正δ-SPHC算法具有較好的低能量耗散特性.在修正算法作用下,機(jī)械能衰減得到有效控制;

(4)在修正δ-SPHC算法作用下,數(shù)值波浪水池波高衰減問題大為改善.在低光滑長度系數(shù) αs=1.3 時(shí),計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好.因此,本修正算法可以有效節(jié)約計(jì)算成本,減少計(jì)算時(shí)間.未來在模擬大尺度、長時(shí)間、遠(yuǎn)距離波浪傳播和波浪與結(jié)構(gòu)物相互作用時(shí),有望發(fā)揮重要作用.