基于δ-SPH方法的多孔潛堤共振反射數(shù)值模擬研究

史寶凱，陳永焜，劉 勇，DOMENICO D. MERINGOLO

(中國海洋大學(xué) 山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，青島 266100)

堆石防波堤是一種常見的海岸工程結(jié)構(gòu)，由于孔隙的存在，該結(jié)構(gòu)不僅能衰減波浪，也使水體交換成為可能。在理論和數(shù)值分析中，由于塊石形狀及分布等較為復(fù)雜，通常將其簡化為具有各向同性的可滲多孔介質(zhì)[1]。為了在防御波浪侵襲、保護(hù)岸線的同時(shí)不影響近岸景觀和水質(zhì)交換，可以采用潛堤代替出水堤[2]。但是單個(gè)多孔潛堤有可能無法滿足防護(hù)需求，可布置多排潛堤，利用多排潛堤的Bragg共振反射效應(yīng)提高掩護(hù)效果。從19世紀(jì)70年代開始，部分學(xué)者發(fā)現(xiàn)天然海灘上平行于海岸的系列沙波由于Bragg共振反射作用，具有良好的掩護(hù)效果。Kirby[3]使用拓展的緩坡方程，研究了波浪在單、雙周期正弦沙波上傳播的共振反射，模擬結(jié)果與試驗(yàn)結(jié)果符合較好。Dalrymple和Kirby[4]考慮波浪入射方向，使用邊界積分法研究了斜向入射波在單周期正弦沙波上傳播的共振反射，結(jié)果表明該方法能準(zhǔn)確模擬較大擾動的沙波地形。Cho等[5]和Jeon等[6]分別通過物理模型試驗(yàn)和VOF方法探究了多排平行潛堤的反射系數(shù)與入射波長的關(guān)系，結(jié)果表明可滲透潛堤反射系數(shù)略小于不可滲透潛堤反射系數(shù)，共振反射效應(yīng)能有效提升潛堤掩護(hù)效果，并推薦使用梯形潛堤。魏欣[7]基于物理模型試驗(yàn)和Boussinesq方程數(shù)值模型，研究了雙排梯形潛堤的反射系數(shù)與水深、潛堤坡度及其間距等因素的關(guān)系。

現(xiàn)有數(shù)值模擬研究多基于網(wǎng)格方法，在處理流體大變形等情況時(shí)對網(wǎng)格要求很高，因此本文選擇SPH(光滑粒子流體動力學(xué))研究多孔潛堤的共振反射現(xiàn)象。該方法由Monaghan[8]提出，是一種無網(wǎng)格Lagerange粒子法。目前已有很多學(xué)者利用SPH方法進(jìn)行了研究，Shao等[9]通過多孔介質(zhì)內(nèi)外采用不同的控制方程，使用SPH模型研究了波浪與多孔介質(zhì)防波堤的相互作用。Akbari等[10]通過引入表征物理量的概念，將多孔介質(zhì)內(nèi)外統(tǒng)一為一套控制方程，研究了波浪在多孔介質(zhì)斜坡堤上的爬坡和越浪過程。Ren等[11]使用體平均/密度權(quán)平均Navier-Stokes方程和LES模型，模擬了波浪通過多孔潛堤和出水堤的過程。為解決WCSPH(弱可壓縮SPH)模型中的壓力震蕩問題，Antuono等[12]在控制方程中加入擴(kuò)散項(xiàng)，提出了δ-SPH模型。經(jīng)過證明，δ-SPH模型是一種準(zhǔn)確、穩(wěn)定的數(shù)值計(jì)算方法，近年來也得到了越來越多的關(guān)注，被應(yīng)用于研究自由表面流動[13]、開孔沉箱[14]、液艙晃蕩[15]等問題。

考慮人工拋石潛堤和天然沙波、珊瑚礁等的滲透性都會對波浪運(yùn)動產(chǎn)生不可忽視的影響；同時(shí)，可以利用共振反射增強(qiáng)潛堤的掩護(hù)效果。本文將建立多孔介質(zhì)的δ-SPH數(shù)值模型，研究多孔潛堤共振反射問題，分析波高、水深、潛堤參數(shù)等對潛堤共振反射的影響，為工程設(shè)計(jì)提供參考。

1 δ-SPH數(shù)值模型

1.1 控制方程

在δ-SPH[12]數(shù)值模型中將流體假定為弱可壓縮的粘性流體，滿足質(zhì)量守恒和動量守恒，控制方程為連續(xù)方程和Navier-Stokes方程。本文引入表征物理量概念，通過在控制方程中引入阻力項(xiàng)，得到方程如下

(1)

(2)

式中：ρ為水的密度，kg/m3；u為表征速度，m/s，表現(xiàn)為流場的真實(shí)速度與孔隙率的乘積；nw為孔隙率；P為壓力，Pa；g表示重力加速度，m/s2；?表示粘度耗散項(xiàng)；α1與α2分別為多孔介質(zhì)的層流和紊流阻力系數(shù)。

壓力和密度之間的關(guān)系通過以下狀態(tài)方程顯式求解

(3)

式中：c0表示人工聲速，其取值要保證密度變化率和時(shí)間增量均控制在合理范圍之內(nèi)；ρ0表示流體的初始密度，取值為1 000 kg/m3。

(4)

(5)

式中：d50表示多孔介質(zhì)的中值粒徑。

利用SPH方法對式(1)和式(2)進(jìn)行離散，得到離散后的控制方程為

(6)

(7)

由于流體的弱可壓縮性，在模擬沖擊問題時(shí)會出現(xiàn)壓力震蕩現(xiàn)象。動量方程中等號右側(cè)第三項(xiàng)為人工粘度項(xiàng)，可以使部分沖擊表面上的動能轉(zhuǎn)化為熱能從而減少震蕩。根據(jù)Monaghan等[19]建議，參數(shù)取值如下

(8)

(9)

連續(xù)性方程中等號右側(cè)第二項(xiàng)為δ-SPH方法引入的擴(kuò)散項(xiàng)，可以進(jìn)一步降低傳統(tǒng)SPH方法在模擬流體壓強(qiáng)場時(shí)可壓縮性帶來的異常高頻震蕩，ψij為密度差值項(xiàng)，可簡單表示為ψij=ρj-ρi。

為了兼顧準(zhǔn)確性和效率，采用龍格庫塔方法進(jìn)行時(shí)間積分[20]。同時(shí)，為了進(jìn)一步保證流場壓力穩(wěn)定，每運(yùn)行30個(gè)時(shí)間步，采用下式對粒子屬性進(jìn)行平滑處理

(10)

式中：WjMLS是指MLS (Moving Least Square) 核函數(shù)

(11)

1.2 數(shù)值造波與消波

為得到穩(wěn)定的規(guī)則波波列，避免數(shù)值水槽內(nèi)發(fā)生二次反射，本文采用動量源造波方法[21]進(jìn)行造波，并在水槽兩端各設(shè)置一個(gè)阻尼區(qū)域用以吸收波浪，防止波浪反射[22]。對于二維數(shù)值波浪水槽，動量源造波項(xiàng)S=(Sx，0)表達(dá)形式為

(12)

式中：ξ=20/Wsc2，Wsc為造波區(qū)寬度，一般取為目標(biāo)波波長；ω為波浪圓頻率；D為源函數(shù)量，可表示為

(13)

為吸收阻尼區(qū)域內(nèi)的波浪，在阻尼區(qū)域的動量方程中加入一阻尼耗散項(xiàng)，該項(xiàng)表達(dá)形式為

(14)

式中：Was為消波區(qū)的布置寬度，x0為消波區(qū)入口的橫坐標(biāo)。

2 數(shù)值模型驗(yàn)證

2.1 數(shù)值波浪水槽設(shè)置

圖1給出了數(shù)值波浪水槽及潛堤的示意圖。二維水槽長32 m，高0.8 m，距左壁7 m處布置一造波源，距造波源8 m處放置兩排梯形多孔防波堤，堤前3.5 m外設(shè)置兩浪高儀，其間距在滿足利用Goda方法[23]分離入反射波計(jì)算要求的前提下，根據(jù)波況不同進(jìn)行調(diào)整。水槽左右兩端為人工消波區(qū)域。試驗(yàn)布置參考魏鑫[7]的物理模型試驗(yàn)并加以拓展，水深0.4～0.6 m，波高0.04～0.06 m，潛堤中心間距W為1.8 m、2.3 m和2.8 m。單個(gè)潛堤斷面為如圖1中所示的梯形(部分工況下退化為矩形或三角形斷面)，頂寬為wt，底寬為wb，高度為wh，斜坡坡度為δk，堤頂水深ht=h-wh。組成潛堤的多孔介質(zhì)孔隙率取0.42，中值粒徑為1.7 cm。全部數(shù)值試驗(yàn)工況設(shè)計(jì)見表1。

圖1 數(shù)值水槽布置及潛堤幾何參數(shù)Fig.1 Layout of computational domain and parameters of breakwaters

表1 試驗(yàn)工況設(shè)計(jì)Tab.1 Design of experimental conditions

2.2 數(shù)值結(jié)果驗(yàn)證

本節(jié)選取典型波浪條件(h=0.5 m，H=0.04 m)，將雙排多孔潛堤反射系數(shù)的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果[7]進(jìn)行對比，驗(yàn)證數(shù)值模型的合理性。潛堤反射系數(shù)利用Goda方法[23]分析得到，并且取至少10個(gè)穩(wěn)定波浪作為輸入數(shù)據(jù)。圖2給出了四組工況下數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果[7]的對比，圖中圓點(diǎn)為數(shù)值模擬結(jié)果，方點(diǎn)為物模試驗(yàn)結(jié)果，可以看出二者符合良好。從圖2中還可以看出，隨著kh值不斷增大(對應(yīng)的波浪周期不斷減小)，潛堤反射系數(shù)先后出現(xiàn)大小兩個(gè)峰值，前者為主頻反射，后者為次頻反射。這一現(xiàn)象即為典型的Bragg共振反射現(xiàn)象，表現(xiàn)為當(dāng)波浪在周期性連續(xù)變化的地形上傳播時(shí)，二者在特定情況下會發(fā)生共振，使得波浪反射效應(yīng)達(dá)到最強(qiáng)。

為了進(jìn)一步驗(yàn)證數(shù)值模型的合理性，在工況W=1.8 m,δk=1∶1下主頻共振時(shí)如圖3所示布置浪高儀，G1～G5浪高儀與1號潛堤前趾在波浪傳播方向上的距離分別為0.4 m、1.05 m、1.55 m、2.2 m和3.2 m。圖4給出該工況下波面變化數(shù)值模擬結(jié)果和試驗(yàn)結(jié)果[7]對比。從圖4中可以看出，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果符合良好。同時(shí)可看出波浪在經(jīng)過G1和G4測點(diǎn)時(shí)，由于潛堤的存在使得波面出現(xiàn)抬升；G5測點(diǎn)處波形前傾，表明該處波面非線性增強(qiáng)，且此處波高較入射波高變小，說明發(fā)生共振反射時(shí)，雙排潛堤掩護(hù)效果較好。

圖3 浪高儀布置示意Fig.3 Arrangement of wave altimeter gages

注：h=0.5 m,H=0.04 m,W=1.8 m,wb=0.8 m,wh=0.25 m,δk=1:1,kh= 0.8圖4 不同測點(diǎn)處波面歷時(shí)曲線的數(shù)值模擬結(jié)果和試驗(yàn)結(jié)果[7]對比Fig.4 Comparisons of free surface elevations at different wave gauges between the numerical and experimental results[7]

3 算例分析與討論

3.1 Bragg共振反射發(fā)生條件

表2給出了主頻反射發(fā)生時(shí)，波長L與2倍潛堤中心間距W的比值關(guān)系。從表中可看出：各組L/2W的數(shù)值在1.0附近波動，即當(dāng)波長約等于2倍潛堤中心間距時(shí)，雙排多孔潛堤發(fā)生主頻共振反射。

表2 主頻共振反射條件Tab.2 Main resonance reflection conditions

3.2 入射波高對反射系數(shù)的影響

圖5給出了h=0.5 m，W=2.3 m，wb=0.8 m，wh=0.25 m，δk= 1:1時(shí)，不同入射波高條件下，反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)1、2、9工況)。在數(shù)值試驗(yàn)過程中未觀察到波浪破碎現(xiàn)象。從圖5可以看出：相對波高H/h= 0.08 、0.10和0.12情況下的反射系數(shù)變化曲線基本重合，說明在波浪未破碎時(shí)，波高變化對共振反射影響較小。

3.3 堤頂水深對反射系數(shù)的影響

潛堤高度固定時(shí)，總水深h的變化會導(dǎo)致堤頂水深ht的變化。圖6-a、6-b分別給出了H=0.04 m，wb= 0.8 m，wh= 0.25 m，δk=1∶1，W=2.3 m和2.8 m時(shí)，不同相對堤頂水深情況下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)1、3、4和5、6工況)。從圖6中可以看出：隨著相對堤頂水深增加，反射系數(shù)整體減小，主要因?yàn)椴ɡ说哪芰考性谒w表層，水深增加會降低潛堤對波浪的影響效果，降低反射系數(shù)；共振反射發(fā)生點(diǎn)向kh增大方向移動，這是因?yàn)闈摰贪ㄖ行拈g距在內(nèi)的整體布置不變，主頻反射時(shí)波長仍為2倍潛堤中心間距，水深的增加使得主頻反射發(fā)生時(shí)的kh值同步增加。

考慮將水深固定，改變潛堤高度，圖7給出了h=0.5 m，H= 0.04 m，W=2.3 m，wb= 0.8 m，δk=∞時(shí)，不同矩形潛堤高度情況下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)11、15、16工況)。從圖7中可以看出：潛堤高度增加導(dǎo)致堤頂水深減小，會引起反射系數(shù)整體增大，但并未影響發(fā)生主頻反射時(shí)的波浪頻率。

3.4 潛堤中心間距對反射系數(shù)的影響

圖8-a、8-b分別給出了H=0.04 m，wb= 0.8 m，wh= 0.25 m，δk=1:1，h=0.5 m和0.6 m時(shí)，不同潛堤中心間距情況下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)1、5、7和4、6工況)。從圖8中可以看出：隨著中心間距增大，主頻反射和次頻反射點(diǎn)均向kh減小方向移動，這仍然是由發(fā)生共振反射時(shí)波長與潛堤中心間距之比規(guī)律決定；主頻反射附近曲線的峰值更加突出，但是有效反射頻帶寬度變小，說明反射效果更好但作用的波長范圍減??；中心間距增加會引起主、次頻反射系數(shù)峰值增大，但前者所受影響較小，后者增大現(xiàn)象更為明顯。

3.5 潛堤寬度對反射系數(shù)的影響

為了探究潛堤寬度對共振反射的影響，本節(jié)選取矩形潛堤作為研究對象。圖9給出了h=0.5 m，H=0.04 m，W=2.8 m，wh= 0.25 m，δk=∞時(shí)，不同潛堤寬度條件下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)10、11、14工況)。從圖9可以看出：潛堤寬度的增大使得主頻反射系數(shù)增加，而相應(yīng)的次頻反射系數(shù)卻略有減??；主頻反射系數(shù)峰值點(diǎn)向kh減小方向移動，主頻共振反射發(fā)生時(shí)波長與2倍潛堤中心間距之比變大。

3.6 潛堤坡度對反射系數(shù)的影響

潛堤頂寬固定時(shí)，坡度的增加會導(dǎo)致潛堤斷面面積的減小。圖10給出了頂寬固定，h=0.5 m，H=0.04 m，W=2.8 m，wt=0.3 m，wh= 0.25 m時(shí)，不同潛堤坡度條件下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)5、8、10工況)。從圖10中可以看出：潛堤坡度的增大使得主、次頻反射點(diǎn)向kh增大方向移動，同時(shí)主頻反射系數(shù)峰值降低；由于潛堤頂寬固定，坡度的增大使得潛堤斷面面積減小，防護(hù)效果也隨之降低，主頻共振反射發(fā)生時(shí)波長與2倍中心間距之比更接近1。

圖11給出了潛堤底寬固定，h=0.5 m，H=0.04 m，W=2.8 m，wb= 0.8 m，wh= 0.25 m時(shí)，不同潛堤坡度條件下反射系數(shù)KR隨kh值的變化曲線(分別對應(yīng)5、11、12、13工況)。此時(shí)潛堤坡度的增加會導(dǎo)致潛堤斷面面積的增加。從圖11中可以看出：主頻共振反射點(diǎn)的橫坐標(biāo)位置沒有變化；當(dāng)坡度δk=5:8時(shí)(斷面形狀退化為三角形)，主頻反射系數(shù)明顯偏低，原因是波浪能量主要集中在水體表面，潛堤靠近水體表面部分面積偏小，無法對波浪形成有效影響；對于其余三種坡度工況，主頻反射系數(shù)略有增加，但是次頻反射系數(shù)減小。

3.7 流場分析

本節(jié)以h=0.5 m，H=0.06 m，W=2.3 m，wb= 0.8，wh= 0.25 m，δk=1∶1工況組為例，分析共振反射時(shí)潛堤周圍的流場特性。圖12給出了該組工況下kh=0.65時(shí)一個(gè)周期內(nèi)的流場圖。從圖12中可以看到：t=0時(shí)一個(gè)波浪到達(dá)潛堤，波峰處水質(zhì)點(diǎn)速度較大(見圖12-a)；t=0.25T時(shí)潛堤的存在使得堤前速度迅速加大，但堤內(nèi)及附近水質(zhì)點(diǎn)仍保持較低流速(見圖12-b)；t=0.5T時(shí)隨著波浪傳播，潛堤上方波面明顯抬升，但水質(zhì)點(diǎn)速度有所減小，同時(shí)觀察到堤后的水質(zhì)點(diǎn)速度已經(jīng)變小(見圖12-c)；t= 0.75T時(shí)波峰傳播至堤間，第二排潛堤的存在使得水質(zhì)點(diǎn)速度又有小幅度增加(見圖12-d)；t=T時(shí)流場進(jìn)入下一個(gè)周期循環(huán)(見圖12-e)。整個(gè)過程中，潛堤內(nèi)部的水質(zhì)點(diǎn)速度有所變化，但速度值一直較?。徊ǚ逶诮?jīng)過潛堤的過程中，水質(zhì)點(diǎn)速度有明顯的增大—減小—小幅增加—再減小的變化趨勢，潛堤后方水域水質(zhì)點(diǎn)速度較小，說明雙排多孔潛堤可以較好地反射入射波能量，提供有效掩護(hù)。

12-e t=T注：h=0.5 m, H=0.06 m, W=2.3 m,wb = 0.8 m, wh = 0.25 m, δk=1:1,kh=0.65圖12 不同時(shí)刻潛堤周圍的流場圖Fig.12 Flow field diagram around the breakwaters at different times

4 結(jié)論

本文基于δ-SPH模型建立了數(shù)值波浪水槽，模擬了波浪與多孔介質(zhì)潛堤相互作用的過程。數(shù)值模擬結(jié)果和物理模型試驗(yàn)結(jié)果符合良好，表明該模型能有效模擬多孔潛堤的共振反射問題，可作為工程設(shè)計(jì)的有效參考手段。通過改變波浪參數(shù)和多孔潛堤參數(shù)，數(shù)值分析了波浪通過多孔潛堤時(shí)的Bragg共振反射現(xiàn)象及其影響因素，研究發(fā)現(xiàn)：(1)當(dāng)入射波長約等于2倍潛堤中心間距時(shí)，雙排多孔潛堤發(fā)生主頻反射；(2)波浪未破碎情況下，入射波高變化對反射系數(shù)影響很?。?3)潛堤高度增大或水深減小會增加潛堤的反射系數(shù)；(4)潛堤中心間距的增加會使主頻反射系數(shù)增大，但有效作用的波長范圍減小；(5)矩形潛堤寬度增大時(shí)，主頻反射系數(shù)隨之增大，次頻反射系數(shù)隨之減小，主頻反射對應(yīng)的波浪頻率減?。?6)梯形潛堤頂寬固定時(shí)，坡度增大會降低主頻反射系數(shù)，而底寬固定時(shí)，坡度增大主頻反射系數(shù)會增大。