基于修正光滑粒子流體動力學算法的低能量耗散數值波浪水池開發1)

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

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

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

引言

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

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

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

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

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

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

1 數值算法

1.1 δ-SPH 模型

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

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

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

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

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

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

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

首先,根據核函數導數的特性[25]有

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

2 振蕩液滴算例驗證

2.1 收斂性分析與精度驗證

振蕩液滴算例無需施加固壁邊界,初始條件簡單,被廣泛應用于SPH 等數值新算法的驗證[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) 為控制液滴初始速度場的參數.在本文模擬中,A(t0)=1 .初始壓力場分布如圖2 所示,壓力場[33]表示為

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

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

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

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

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

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

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

圖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 時,振蕩液滴在不同時刻的理論形態與修正算法δ-SPHC模擬的形態.可以發現,在兩種光滑長度條件下,δ-SPHC模擬的振蕩形態均與理論值吻合良好,且在t=8.75T周期后,液滴仍保持較為光滑的形態.由此可見,本修正算法具有較高的精度.

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

2.2 傳統δ-SPH 與修正δ-SPHC 算法的對比研究

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

圖7 δ-SPH 與δ-SPHC 動能比較:總體圖(上)與局部放大圖(下),其中黑矩形框為放大區域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呈現了兩種算法中,機械能衰減率的時歷曲線.δ-SPHC算法在兩種光滑長度系數下的計算結果中,機械能衰減率均較小,基本接近零,說明本文修正算法能夠有效減少機械能的非物理性耗散.其中,在光滑長度系數 αs=2.0 時效果更為顯著.值得一提的是,δ-SPHC算法在αs=2.0 時的機械能衰減率相比δ-SPH 算法在αs=1.3 時的機械能衰減率更小,表明本修正SPH 算法的一大優勢是能在較小光滑長度系數下實現較為準確的模擬,從而顯著減少計算時間.

圖9 不同算法和光滑長度系數下,振蕩液滴機械能衰減率時歷曲線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 內存的個人電腦上進行計算.

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

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

3 SPH 數值造波與試驗驗證

3.1 試驗與數值水池設置

3.1.1 試驗水池

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

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

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

3.1.2 SPH 數值水池設置

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

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

圖12 SPH 數值水池示意圖,紅豎線表示波高儀,選取距離波高儀一個粒子間距內自由面粒子的最大高度為波面高度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 模擬的時間步長與試驗中輸出的造波板運動數據時間間隔不一致,采用線性插值技術,根據物理試驗中造波板運動特性,更新SPH 模擬中的造波板位移xS PH、速度vS PH及加速度aS PH,即

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

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

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

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

圖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 算法中的波浪衰減現象

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

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

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

圖15 在αs=1.3 ,d/?x=60 條件下,傳統δ-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 條件下,傳統δ-SPH 模擬的波面形態Fig.16 Wave surface simulated by δ-SPH method with αs=2.0 and d/?x=60

圖17 在αs=2.0,d/?x=60 條件下,傳統δ-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 規則波模擬

4.1.1 收斂性分析與精度驗證

首先,基于規則波傳播算例,校核本文修正的δ-SPHC算法精度與收斂性.如圖18 所示,在αs=1.3時,選取xprobe=2.37 m 處數值波高結果與試驗結果進行比較,可以看出,水深方向不同數量粒子的計算結果略有差別,且隨著d/?x的增大,數值波高趨近于試驗波高.在d/?x=60 時,SPH 結果與試驗結果吻合良好.與振蕩液滴算例類似,圖18 方框中的峰值可能是由于粒子由初始的格子狀分布到各向同性均勻分布時引起的能量瞬時誤差.試驗波高曲線與SPH 數值結果曲線在初始階段第一個峰值存在相位差,可能是由于試驗波高采集與造波板運動存在一定的時間差.在αs=2.0 時,水深方向不同粒子數量的計算結果與試驗結果對比如圖19 所示.同樣可以看到,在d/?x=60 時,數值波高曲線與試驗結果基本吻合.另外,由于物理試驗水池長度僅有16 m,消波區長度設置有限,易導致試驗后期消波不徹底,造成波浪周期的變化,即隨著時間增加,周期略小于1 s;而數值水池長度L=60 m,能較為準確地重現規則波傳播,并且確保周期保持不變.因此,可以看到,在圖18、圖19 中,t=20 s 之后SPH 結果的波峰和波谷時刻與試驗值存在略微誤差,但是波高和波谷的幅值與試驗值吻合良好.基于以上討論可以得出,在水深方向粒子數d/?x=60 條件下,計算結果能夠達到預期精度要求.因此,下文將在d/?x=60 粒子設置下,進行修正δ-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 呈現了在xprobe=2.37 m 處,規則波波高數值模擬結果與試驗結果的對比.采用傳統δ-SPH 計算時,光滑長度系數 αs=1.3 條件下,由于核函數緊支域內粒子數量少,數值結果與試驗結果相差較大,而光滑長度系數 αs=2.0 時,數值波高結果接近試驗結果.但值得注意的是,采用δ-SPHC,光滑長度系數αs=1.3 時,數值結果與試驗結果能較好吻合,基本達到傳統δ-SPH 算法在光滑長度系數 αs=2.0 時的計算精度,與上文振蕩液滴算例得出的結論一致,再次說明了本文修正δ-SPHC算法的優勢,即在光滑長度系數 αs=1.3,2.0 時,計算結果與試驗結果均能較好吻合.圖21 展示了數值水池的壓力場,可見流場壓力分布均勻,在運動較為劇烈的推板附近,也能確保壓力場穩定.基于以上討論,可以認為修正的δ-SPHC算法具有較高的精度.

圖20 xprobe=2.37 m 處,不同SPH 算法下的波面高度數值結果與試驗值對比Fig.20 Comparison between different SPH simulations and experimental data for wave elevations at xprobe=2.37 m

圖21 規則波壓力云圖(δ-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的抗能量非物理性衰減結果分析

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

圖22 αs=1.3,d/?x=60 時,δ-SPHC 算法模擬的波面形態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 算法模擬的波面形態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 的波高曲線,可見相比于傳統δ-SPH 模型,修正的δ-SPHC算法在αs=1.3 時的波高結果更加準確.值得一提的是,其計算結果與傳統δ-SPH 模型在αs=2.0 時波高結果基本重合,并與試驗波高吻合良好.此外,采用δ-SPHC修正算法,在光滑長度系數 αs=2.0 時,波高計算結果與試驗波高高度吻合.

圖26 在xprobe=6.37 m 處,規則波波面高度數值模擬結果與試驗結果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為試驗測量的平均波高.可見,在兩種光滑長度系數下,δ-SPHC算法在較遠處的波高監測結果與試驗相對誤差均較小.本節討論表明,修正算法能有效解決能量非物理性衰減問題,且在較低光滑長度系數時能得到更為準確的結果.

4.1.3 不同人工黏性系數下δ-SPHC模擬結果分析

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

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

4.2 基于δ-SPHC 的不規則波模擬

4.2.1 收斂性分析與精度驗證

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

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

圖29 不規則波壓力云圖(δ-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抗能量非物理性衰減結果分析

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

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

5 結論與展望

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

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

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

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

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