基于SPH-ALE方法的波浪水槽數(shù)值模擬

李雪臨,路 寬,韓林生,王 靜,王花梅,石建軍,朱 銳

(國家海洋技術(shù)中心,天津 300112)

?

基于SPH-ALE方法的波浪水槽數(shù)值模擬

李雪臨,路 寬,韓林生,王 靜,王花梅,石建軍,朱 銳

(國家海洋技術(shù)中心,天津 300112)

為了精確模擬波浪傳播,基于光滑粒子流體動力學(xué)及任意拉格朗日歐拉(SPH-ALE)方法建立二維數(shù)值波浪水槽,在原SPH方法中引入近似黎曼求解器替代人工黏性項,采用排斥力邊界條件防止流體粒子穿透固邊界,海綿層內(nèi)采用指數(shù)型衰減函數(shù)來消除水槽末端的波浪反射,并對規(guī)則波的傳播進(jìn)行數(shù)值模擬。結(jié)果表明:與采用人工黏性項的原SPH方法相比,SPH-ALE方法能夠無衰減地模擬波浪傳播,并可有效減小固邊界附近的粒子壓力振蕩。

波浪傳播;光滑粒子;流體動力學(xué);任意拉格朗日歐拉方法;數(shù)值模擬；波浪水槽;排斥力邊界條件;海綿層

光滑粒子流體動力學(xué)(SPH)方法因其無網(wǎng)格、自適應(yīng)以及拉格朗日性質(zhì),不需要特殊處理即可實現(xiàn)對復(fù)雜自由表面的精確模擬,特別適用于處理含大變形自由表面的流體動力學(xué)問題。Monaghan[1]于1994年首次將SPH方法應(yīng)用于自由表面流動的數(shù)值模擬。Doring等[2]采用SPH方法研究了楔體入水和波浪與方箱的相互作用問題,模擬得到的楔體入水過程中垂向加速度和角加速度的歷時曲線與實驗結(jié)果吻合較好;Padova等[3]對不同Ursell數(shù)對應(yīng)的非線性波和JONSWAP譜對應(yīng)的隨機(jī)波在規(guī)則或不規(guī)則底坡上的傳播和崩破、卷破、破碎進(jìn)行了SPH數(shù)值模擬和物理模型試驗研究,探討了不同的人工黏性經(jīng)驗系數(shù)和粒子初始間距與光滑長度之比對計算精度和效率的影響;高睿等[4]建立了求解孤立波沿斜坡傳播、變形和破碎的數(shù)值模型,分析了孤立波淺化過程中波面、波高的變化特征以及波浪的不同破碎形態(tài);鄭興等[5]提出了二階精度近似的SPH方法(K2_SPH),并對潰壩、涌浪和大幅液艙晃蕩等強(qiáng)非線性自由表面問題進(jìn)行了研究,驗證了K2_SPH方法在計算精度上的提升。以上研究中動量方程的擴(kuò)散項均采用Monaghan提出的人工黏性項,但在利用數(shù)值波浪水槽研究波浪傳播及其與結(jié)構(gòu)物的相互作用時,人工黏性在大多數(shù)情況下往往因耗散過度而造成波浪傳播過程的能量衰減[6],限制了此類SPH方法的應(yīng)用。

為了避免人工黏性耗散性的影響,任冰等[7]基于Parshikov等[8]提出的Godunov型粒子動力學(xué)(GPH)方法,在修正的SPH方法(CSPM)中引入近似黎曼解,對動量方程的壓力梯度項進(jìn)行修正,并合理設(shè)置結(jié)構(gòu)物的沖擊邊界,利用所建立的二維數(shù)值波浪水槽模擬了規(guī)則波對浪濺區(qū)結(jié)構(gòu)物的沖擊過程,并利用物理模型試驗結(jié)果進(jìn)行了驗證。由于不需要人工黏性來保持?jǐn)?shù)值穩(wěn)定,波浪傳播過程中的能量衰減較小,同時可減少壓力場的數(shù)值噪音。Vila[9]將任意拉格朗日-歐拉方法與SPH方法相結(jié)合,提出另一種混合的SPH-ALE方法,同樣引入近似黎曼解替代人工黏性,該方法采用守恒變量描述控制方程,與采用原始變量的Parshikov方法相比具有更小的數(shù)值振蕩。

本文基于Vila提出的SPH-ALE方法建立數(shù)值波浪水槽,采用排斥力邊界條件防止流體粒子穿透固邊界,采用海綿層消波技術(shù)消除水槽末端的波浪反射,利用所建立的二維數(shù)值波浪水槽對規(guī)則波的傳播進(jìn)行了數(shù)值模擬,并與采用人工黏性的原SPH方法的計算結(jié)果進(jìn)行了對比。

1 數(shù)值模型

主要對比采用人工黏性項的SPH模型和采用守恒變量描述的SPH-ALE模型。兩種模型采用相同的邊界條件和時間積分算法,主要區(qū)別在于控制方程中的連續(xù)方程和動量方程。計算程序是在SPHysics[10]的基礎(chǔ)上進(jìn)行二次開發(fā)。

1.1 控制方程的離散

1.1.1 連續(xù)方程和動量方程

a. SPH方法

采用人工黏性項的SPH粒子近似方程為

(1)

式中:r為粒子位置;v為粒子速度;ρ為粒子密度;m為粒子質(zhì)量;P為粒子壓力;g為重力加速度;Wij為光滑核函數(shù);i、j分別代表相互作用的兩個粒子。人工黏性項Πij為

(2)

b. SPH-ALE方法

Vila考慮求解相互作用的兩個粒子間的一維黎曼問題，SPH中引入近似黎曼解,推導(dǎo)得到采用守恒變量描述的SPH-ALE方程為

(3)

本文采用HLLC(Harten Lax van Leer constant)近似黎曼求解器求解兩粒子接觸不連續(xù)的一維黎曼問題,采用二階迎風(fēng)格式MUSCL(monotone upstream centered schemes for conservation laws)[11],利用斜率限制器β重構(gòu)粒子間界面兩側(cè)的狀態(tài)并計算界面處通量,避免了Godunov一階迎風(fēng)格式的耗散性,從而具有更高的計算精度。

SPH-ALE方法的主要優(yōu)點是能夠有效減小流體粒子的壓力和速度振蕩[12]。由于不需要采用人工黏性項保持?jǐn)?shù)值穩(wěn)定性,該方法可有效地減少數(shù)值耗散,能更精確地模擬波浪傳播及壓力場,從而精確預(yù)測作用于結(jié)構(gòu)物上的波浪力。

1.1.2 狀態(tài)方程

在弱可壓縮流體SPH方法中,采用Tait狀態(tài)方程計算粒子壓力:

(4)

式中:γ為絕熱指數(shù);c0為參考密度ρ0下的聲速。由于弱可壓縮條件的限制,聲速應(yīng)為流體最大速度的10倍以上,以保證流體密度的變化在1%以內(nèi)。

1.1.3 XSPH方法

Monaghan[13]提出的XSPH方法能夠保持粒子有序,經(jīng)修正后粒子的速度由式(5)表示:

(5)

式中參數(shù)ε可取0.5。XSPH方法使粒子i的速度更接近于其附近粒子的平均速度。

1.2 邊界條件

1.2.1 排斥力邊界條件

采用由Monaghan等[14]提出、Rogers等[15]修正的排斥力邊界條件,能有效防止流體粒子穿透固壁,且易于處理復(fù)雜邊界。邊界粒子施加在流體粒子上、沿固邊界法線方向的作用力為

f=nR(ψ)P(ξ)ε(z,u⊥)

(6)

式中:n為邊界粒子處的法向量;ψ為流體粒子到固邊界的垂直距離,ξ為流體粒子與邊界粒子確定的向量在該邊界粒子處切向的投影,邊界粒子處的法向量和切向量由與其相鄰的邊界粒子坐標(biāo)來確定;u⊥為流體粒子速度在固邊界法向上的投影。

1.2.2 造波邊界條件

在SPH方法中通過給定運動固邊界粒子的位移和速度,模擬造波板生成目標(biāo)波浪。本文采用推板式造波,生成線性規(guī)則波時造波板粒子的水平位移Xp和速度Up分別如式(7)和式(8)所示:

ωt

(7)

(8)

式中:S為推板沖程;ω為波浪圓頻率。

1.2.3 海綿層消波

采用類似于網(wǎng)格法的海綿層消波技術(shù),在數(shù)值水槽末端設(shè)置一定長度的海綿層(至少一個波長)吸收反射波,選擇指數(shù)型衰減函數(shù)[16]對位于海綿層內(nèi)的流體粒子速度進(jìn)行衰減,衰減函數(shù)形式如式(9)所示:

f(x)=1-exp(-α(Ls-(x-x0)))

(9)

式中:Ls為海綿層長度;x0為海綿層起始位置坐標(biāo);α為衰減系數(shù),針對具體問題通過數(shù)值試驗合理確定。

1.3 時間積分算法

本文采用的Sympletic時間積分算法在無摩擦或無黏性情況下具備時間可逆性[17],非常適合應(yīng)用于無網(wǎng)格粒子法。Sympletic算法在半個時間步后粒子密度和位置為

(10)

(11)

2 波浪傳播的數(shù)值模擬

2.1 計算參數(shù)設(shè)置

數(shù)值計算域如圖1所示。數(shù)值水槽的長度L=17 m,推板式造波機(jī)位于水槽左端,距離水槽左邊界0.1 m,入射波為規(guī)則波,波高H=0.1 m,周期T=1.5 s,波長Lw=2.83 m。水槽右端設(shè)置長度Ls=3 m的海綿層以防止波浪反射,衰減系數(shù)取α=3.0。水槽內(nèi)水深d=0.5 m,初始狀態(tài)下流體粒子呈均勻等距分布，Δx=Δz=0.02 m,計算域內(nèi)的粒子總數(shù)為22 104個,其中邊界粒子1 004個,流體粒子21 100個。

圖1 數(shù)值計算域示意圖(單位：m)

分別采用SPH方法和SPH-ALE方法對水槽內(nèi)的波浪傳播進(jìn)行數(shù)值模擬。SPH方法中人工黏性項的自由參數(shù)υ取0.1或0.5,υ<0.1時可能會造成數(shù)值計算不穩(wěn)定。SPH-ALE方法的計算結(jié)果對斜率限制器β的取值比較敏感,取值較小時(β≈1.0)仍會導(dǎo)致波浪傳播的衰減,取值較大時(β≈1.5)則會造成波幅隨時間增長的非物理效應(yīng)[19],這里取β=1.3。

為了得到瞬時波面和波面歷時曲線,采用Gómez-Gesteira等[20]提出的自由表面確定方法,從自由表面之上某一位置開始,以0.01 h為空間步長垂直向下按照式(12)依次估算各點的密度值。當(dāng)估算密度首次超過參考密度的1/2時,將此位置標(biāo)記為自由表面所在位置。

(12)

2.2 數(shù)值結(jié)果分析

圖2為SPH方法、SPH-ALE方法計算的距造波板4 m處的波面歷時曲線與線性理論解的對比結(jié)果(圖中η為波面高程)。從圖2可以看出,兩種方法模擬的波浪相位均與線性理論解吻合較好,但SPH方法因人工黏性的耗散計算波高值明顯偏小,υ=0.1時波面穩(wěn)定后的計算平均波高比目標(biāo)波高偏低約30%;而SPH-ALE方法的計算波高與目標(biāo)波高偏差很小,波面穩(wěn)定后的計算平均波高與目標(biāo)波高的偏差僅為2.4%。此外,SPH-ALE方法計算的波面為非標(biāo)準(zhǔn)的正弦曲線,波峰處較尖陡而波谷處較平坦,能夠比較真實地反映出波浪的非線性。

圖2 波面歷時曲線的數(shù)值解與線性理論解對比

海綿層內(nèi)距造波板16 m處兩種方法計算的波面歷時曲線如圖3所示。由圖3可以看到,兩種數(shù)值方法均表現(xiàn)出很好的消波效果,在距離海綿層起始位置2 m處的波面高程差最大不超過0.003 m,小于入射波高的3%。

圖3 海綿層內(nèi)(x=16 m)波面歷時曲線

圖4 瞬時波面計算值與線性理論解對比

圖4分別為t=10 s和t=15 s時,兩種方法計算得到的瞬時波面與線性理論解的對比結(jié)果。其中,橫坐標(biāo)代表距造波板的水平距離。SPH方法人工黏性項的自由參數(shù)取υ=0.1和υ=0.5兩種情況。由圖4可看出,SPH方法計算的瞬時波面因人工黏性的耗散沿水槽長度方向逐漸衰減,且隨人工黏性項自由參數(shù)υ的增大而衰減加劇,而自由參數(shù)υ取值過小又會造成數(shù)值計算的不穩(wěn)定,因此使用SPH方法難以保證波浪傳播數(shù)值模擬的精度。SPH-ALE方法計算的瞬時波面沿水槽長度方向幾乎無任何衰減,與線性理論解的相位和波高值均吻合較好,同時波面形狀也反映出了波浪的非線性,波面相對于線性理論解整體上略有抬升。

圖5為t=15 s時SPH方法(υ=0.1)和SPH-ALE方法計算得到的流體壓力場分布圖。由圖5可以看出,兩種方法算得的自由表面都比較光滑,自由表面附近的壓力分布均勻合理;但在造波板和水槽底部邊界附近(圖中黑色封口框范圍),SPH方法的計算結(jié)果呈現(xiàn)出明顯的粒子壓力不連續(xù)現(xiàn)象,而SPH-ALE方法算得的壓力場仍分布均勻,能夠有效減小固邊界附近的粒子壓力振蕩。

圖5 計算流體壓力場分布

3 結(jié) 論

基于SPH-ALE方法建立二維數(shù)值波浪水槽,分別采用SPH方法和SPH-ALE方法對規(guī)則波在數(shù)值水槽內(nèi)的傳播進(jìn)行數(shù)值模擬,數(shù)值解和線性理論解的對比結(jié)果表明:因人工黏性的耗散作用,SPH方法計算的波面歷時曲線和瞬時波面與線性理論解相比均有明顯的衰減,而SPH-ALE方法的計算波面與線性理論解的相位和波高值均吻合很好,并可真實反映波浪的非線性;此外,在造波板和水槽底部邊界附近,SPH方法計算的壓力場呈現(xiàn)出明顯的不連續(xù)現(xiàn)象,而SPH-ALE方法能夠有效地減小固邊界附近的粒子壓力振蕩。

[ 1 ] MONAGHAN J J.Simulating free surface flows with SPH[J].Journal of Computational Physics,1994,110(2):399-406.

[ 2 ] DORING M,OGER G,ALESSANDRINI B,et al.SPH simulations of floating bodies in waves[C]//Proceedings of ASME 2004 23rd International Conference on Offshore Mechanics and Arctic Engineering.New York: American Society of Mechanical Engineers,2004: 741-747.

[ 3 ] PADOVA D,DALRYMPLE R A,MOSSA M,et al.SPH simulations of regular and irregular waves and their comparison with experimental data[EB/OL].[2015-09-08].http://arxiv.org/ftp/arxiv/papers/0911/0911.1872.pdf.[ 4 ] 高睿,任冰,王國玉,等.孤立波淺化過程的SPH數(shù)值模擬[J].水動力學(xué)研究與進(jìn)展(A輯),2010,25(5):620-629.(GAO Rui,REN Bing,WANG Guoyu,et al.SPH model of solitary waves shoaling on a mild sloping beach[J].Chinese Journal of Hydrodynamics,2010,25(5):620-629.(in Chinese))

[ 5 ] 鄭興,段文洋.K2_SPH方法及其對二維非線性水波的模擬[J].計算物理,2011,28(5):659-666.(ZHENG Xing,DUAN Wenyang.K2-SPH method and application for 2D nonlinear water wave simulation[J].Chinese Journal of Computational Physics,2011,28(5):659-666.(in Chinese))

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

[ 7 ] 任冰,高睿,金釗,等.波浪對透空式結(jié)構(gòu)物沖擊作用的光滑粒子流體動力學(xué)數(shù)值模擬[J].海洋學(xué)報,2012,34(1):163-177.(REN Bing,GAO Rui,JIN Zhao,et al.The numerical simulation of wave slamming on an open-piled structure using the SPH model[J].Acta Oceanologica Sinica,2012,34(1): 163-177.(in Chinese))

[ 8 ] PARSHIKOV A N,MEDIN S A,LOUKASHENKO I I,et al.Improvements in SPH method by means of interparticle contact algorithm and analysis of perforation tests at moderate projectile velocities[J].International Journal of Impact Engineering,2000,24(8):779-796.

[ 9 ] VILA J P.On particle weighted methods and smooth particle hydrodynamics[J].Mathematical Models and Methods in Applied Sciences,1999,9(2):161-209.

[11] TORO E F.Shock-capturing methods for free-surface shallow flows[M].San Francisco: Wiley,2001.

[12] ROGERS B D,DALRYMPLE R A,STANSBY P K.Simulation of caisson breakwater movement using 2-D SPH[J].Journal of Hydraulic Research,2010,48(Sup1):135-141.

[13] MONAGHAN J J.On the problem of penetration in particle methods[J].Journal Computational Physics,1989,82(1):1-15.

[14] MONAGHAN J J,KOS A.Solitary waves on a Creatan Beach[J].Journal of Waterway,Port,Coastal and Ocean Engineering,1999,125(3):145-154.

[15] ROGERS B D,DALRYMPLE R A.SPH modeling of tsunami waves[J].Advances in Coastal and Ocean Engineering,2008,10:75-100.

[16] XU Rui.An improved incompressible smoothed particle hydrodynamics method and its application in free-Surface simulations[D].Manchester: University of Manchester,2009.

[17] LEIMKUHLER B J,REICH S,SKEEL R D.Integration methods for molecular dynamic IMA volume in mathematics and its application[M].Berlin: Springer,1997.

[18] MONAGHAN J J.Smoothed particle hydrodynamics[J].Reports on Progress in Physics,2005,68(8):1703-1759.

[19] ROGERS B D.Refined localized modelling of coastal flow features using adaptive quadtree grids[D].Oxford: Uiniversity of Oxford,2001.

Numerical simulation of wave flume based on SPH-ALE method

//LI Xuelin, LU kuan, HAN Linsheng, WANG Jing, WANG Huamei, SHI Jianjun, ZHU Rui

(NationalOceanTechnologyCenter,Tianjin300112,China)

In order to simulate wave propagation accurately, a two-dimensional numerical wave flume was constructed based on the smoothed particle hydrodynamics and arbitrary Lagrange Euler (SPH-ALE) method. An approximate Riemann solver was introduced to replace the artificial viscosity in the standard SPH model. The repulsive force boundary condition was used at the solid wall to prevent fluid particles from penetrating the solid boundary. The exponential damping function was used in the sponge layer to eliminate wave reflection at the end of the flume. The regular wave propagation was simulated in the constructed wave flume. The results show that, in comparison with the standard SPH method using artificial viscosity, the SPH-ALE method can simulate the wave propagation without attenuation and effectively reduce the pressure fluctuations near the solid boundaries.

wave propagation; smoothed particle; hydrodynamics; arbitrary Lagrange Euler method; numerical simulation; wave flume; repulsive force boundary condition; sponge layer

10.3880/j.issn.1006-7647.2016.06.008

國家自然科學(xué)基金 (51309056)

李雪臨(1981—),男,工程師,博士,主要從事波浪與結(jié)構(gòu)物相互作用研究。E-mail:xuedut@163.com

TV139.2

A

1006-7647(2016)06-0039-05

2015-09-08 編輯：鄭孝宇)