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

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

(中國海洋大學 山東省海洋工程重點實驗室，青島 266100)

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

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

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

1 δ-SPH數值模型

1.1 控制方程

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

(1)

(2)

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

壓力和密度之間的關系通過以下狀態方程顯式求解

(3)

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

(4)

(5)

式中：d50表示多孔介質的中值粒徑。

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

(6)

(7)

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

(8)

(9)

連續性方程中等號右側第二項為δ-SPH方法引入的擴散項，可以進一步降低傳統SPH方法在模擬流體壓強場時可壓縮性帶來的異常高頻震蕩，ψij為密度差值項，可簡單表示為ψij=ρj-ρi。

為了兼顧準確性和效率，采用龍格庫塔方法進行時間積分[20]。同時，為了進一步保證流場壓力穩定，每運行30個時間步，采用下式對粒子屬性進行平滑處理

(10)

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

(11)

1.2 數值造波與消波

為得到穩定的規則波波列，避免數值水槽內發生二次反射，本文采用動量源造波方法[21]進行造波，并在水槽兩端各設置一個阻尼區域用以吸收波浪，防止波浪反射[22]。對于二維數值波浪水槽，動量源造波項S=(Sx，0)表達形式為

(12)

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

(13)

為吸收阻尼區域內的波浪，在阻尼區域的動量方程中加入一阻尼耗散項，該項表達形式為

(14)

式中：Was為消波區的布置寬度，x0為消波區入口的橫坐標。

2 數值模型驗證

2.1 數值波浪水槽設置

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

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

表1 試驗工況設計Tab.1 Design of experimental conditions

2.2 數值結果驗證

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

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

圖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 不同測點處波面歷時曲線的數值模擬結果和試驗結果[7]對比Fig.4 Comparisons of free surface elevations at different wave gauges between the numerical and experimental results[7]

3 算例分析與討論

3.1 Bragg共振反射發生條件

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

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

3.2 入射波高對反射系數的影響

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

3.3 堤頂水深對反射系數的影響

潛堤高度固定時，總水深h的變化會導致堤頂水深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時，不同相對堤頂水深情況下反射系數KR隨kh值的變化曲線(分別對應1、3、4和5、6工況)。從圖6中可以看出：隨著相對堤頂水深增加，反射系數整體減小，主要因為波浪的能量集中在水體表層，水深增加會降低潛堤對波浪的影響效果，降低反射系數；共振反射發生點向kh增大方向移動，這是因為潛堤包括中心間距在內的整體布置不變，主頻反射時波長仍為2倍潛堤中心間距，水深的增加使得主頻反射發生時的kh值同步增加。

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

3.4 潛堤中心間距對反射系數的影響

圖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時，不同潛堤中心間距情況下反射系數KR隨kh值的變化曲線(分別對應1、5、7和4、6工況)。從圖8中可以看出：隨著中心間距增大，主頻反射和次頻反射點均向kh減小方向移動，這仍然是由發生共振反射時波長與潛堤中心間距之比規律決定；主頻反射附近曲線的峰值更加突出，但是有效反射頻帶寬度變小，說明反射效果更好但作用的波長范圍減??；中心間距增加會引起主、次頻反射系數峰值增大，但前者所受影響較小，后者增大現象更為明顯。

3.5 潛堤寬度對反射系數的影響

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

3.6 潛堤坡度對反射系數的影響

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

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

3.7 流場分析

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

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 不同時刻潛堤周圍的流場圖Fig.12 Flow field diagram around the breakwaters at different times

4 結論

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