黃曉婷 孫鵬楠 ,2) 呂鴻冠 鐘詩蘊

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

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

引言

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

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

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

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

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

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

1 數(shù)值算法

1.1 δ-SPH 模型

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

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

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

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

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

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

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

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

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

2 振蕩液滴算例驗證

2.1 收斂性分析與精度驗證

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

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

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

A(t0) 為控制液滴初始速度場的參數(shù).在本文模擬中,A(t0)=1 .初始壓力場分布如圖2 所示,壓力場[33]表示為

圖2 振蕩液滴初始壓力分布Fig.2 Initial pressure distribution

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

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

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

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

圖4 線動量時歷曲線Fig.4 Time evolution of the linear momentum

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

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

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

圖6 αs=1.3(上)和αs=2.0(下)時,振蕩液滴形態(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ù)值精度和動量與能量守恒性.因此,下文將在粒子分辨率R/?x=200 時,依托該算例進一步討論不同算法抗能量衰減效果的差異.

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

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

圖7 δ-SPH 與δ-SPHC 動能比較:總體圖(上)與局部放大圖(下),其中黑矩形框為放大區(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

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

圖9 不同算法和光滑長度系數(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)存的個人電腦上進行計算.

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

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

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

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

3.1.1 試驗水池

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

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

圖11 造波推板橫向位移隨時間變化:(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]能實現(xiàn)固壁邊界的精確處理,有效防止壁面穿透,且保證壁面處壓力連續(xù),因而得到了廣泛應(yīng)用.

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

圖12 SPH 數(shù)值水池示意圖,紅豎線表示波高儀,選取距離波高儀一個粒子間距內(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

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

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

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

3.2 δ-SPH 算法在高粒子分辨率下的精度驗證

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

圖13 αs=2.0 時,不同粒子分辨率下波面高度的δ-SPH 模擬值與試驗值對比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 條件下,進行規(guī)則波造波模擬.為了測定波高變化曲線,沿水池長度方向,距離造波板初始位置xprobe=2.37 m,4.37 m,5.37 m,6.37 m,15 m 處設(shè)置5 個數(shù)值波高儀.

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

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

圖18 在αs=1.3 和 d/?x=15,30,60 條件下,δ-SPHC 計算的波面高度時歷曲線與試驗值對比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 計算的波面高度時歷曲線與試驗值對比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é)果與試驗結(jié)果的對比.采用傳統(tǒng)δ-SPH 計算時,光滑長度系數(shù) αs=1.3 條件下,由于核函數(shù)緊支域內(nèi)粒子數(shù)量少,數(shù)值結(jié)果與試驗結(jié)果相差較大,而光滑長度系數(shù) αs=2.0 時,數(shù)值波高結(jié)果接近試驗結(jié)果.但值得注意的是,采用δ-SPHC,光滑長度系數(shù)αs=1.3 時,數(shù)值結(jié)果與試驗結(jié)果能較好吻合,基本達(dá)到傳統(tǒng)δ-SPH 算法在光滑長度系數(shù) αs=2.0 時的計算精度,與上文振蕩液滴算例得出的結(jié)論一致,再次說明了本文修正δ-SPHC算法的優(yōu)勢,即在光滑長度系數(shù) αs=1.3,2.0 時,計算結(jié)果與試驗結(jié)果均能較好吻合.圖21 展示了數(shù)值水池的壓力場,可見流場壓力分布均勻,在運動較為劇烈的推板附近,也能確保壓力場穩(wěn)定.基于以上討論,可以認(rèn)為修正的δ-SPHC算法具有較高的精度.

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

圖22 αs=1.3,d/?x=60 時,δ-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 時,δ-SPHC 模擬的不同位置波面高度時歷曲線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 時,δ-SPHC 算法模擬的波面形態(tài)Fig.24 The wave surface simulated by δ-SPHC with αs=2.0 and

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

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

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

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

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

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

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

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

4.2.1 收斂性分析與精度驗證

經(jīng)過上文分析,在光滑長度系數(shù) αs=1.3 時,修正算法δ-SPHC能夠準(zhǔn)確模擬規(guī)則波傳播,而本節(jié)算法,在αs=1.3 時,模擬和驗證不規(guī)則波傳播.選取d/?x=15,30,60 進行收斂性分析,如圖28 所示,在水深方向粒子數(shù)d/?x=15,30,60 時,在xprobe=2.37處探測的波高與試驗值基本一致.此外,在不規(guī)則波模擬中,δ-SPHC實現(xiàn)了穩(wěn)定的壓力場預(yù)報(見圖29),進一步說明了δ-SPHC方法所構(gòu)建的數(shù)值水池的精度.由于d/?x=60 時,數(shù)值結(jié)果已基本收斂于試驗值,因此下文選取d/?x=60 進行不規(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 處波高時歷曲線,可見在傳統(tǒng)δ-SPH 計算結(jié)果中,與試驗結(jié)果相比,波高衰減明顯;而修正算法結(jié)果與試驗結(jié)果基本吻合,說明不規(guī)則波模擬中波高衰減現(xiàn)象也得到較好改善.在δ-SPHC計算框架下,數(shù)值水池生成不規(guī)則波的能量非物理性耗散問題得到解決,并能夠在光滑長度系數(shù) αs=1.3 時,實現(xiàn)較高的造波精度.

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

5 結(jié)論與展望

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

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

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

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

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