近岸波浪在剛性植被區域傳播的數值模型

王 磊,房克照, 尹 晶, 孫家文

(1. 大連理工大學 海岸和近海工程國家重點實驗室， 遼寧 大連 116024； 2. 國家海洋環境監測中心， 遼寧 大連 116023)

近岸波浪在剛性植被區域傳播的數值模型

王 磊1,房克照1, 尹 晶2, 孫家文2

(1. 大連理工大學 海岸和近海工程國家重點實驗室， 遼寧 大連 116024； 2. 國家海洋環境監測中心， 遼寧 大連 116023)

基于擴展型Boussinesq水波方程，建立了波浪在剛性植被覆蓋的近岸海域傳播的數值模型。通過在動量方程源項中引入拖曳阻力項考慮植被對波浪的衰減作用。控制方程采用有限差分和有限體積混合格式求解，模型穩定性強，具備間斷捕捉能力，能有效模擬近岸區域波浪的傳播變形、破碎和處理海岸動邊界問題。利用所建立模型對典型物理模型實驗進行模擬，計算結果與實驗結果吻合良好，表明模型可用于波浪在剛性植被覆蓋海域的數值計算。

Boussinesq水波方程;剛性植被;消浪;數值模擬；波浪傳播

近海岸波浪傳播變形及堤壩的防浪護坡，是目前港口和海岸建設中要考慮的重要工程技術問題。傳統的護岸往往采用混凝土、漿砌石等硬質材料，對生態系統會造成一定程度的破壞，干擾生態平衡。而采用自然的植被護岸，不僅能消減波能，起到防浪護岸的作用，還能維持生態系統和水動力系統的平衡，改善生態環境[1-3]。鑒于此，眾多學者針對植被消減波能效果以及對相關水動力的影響開展研究[2-11]。

考慮到問題的復雜性和現場實測及物理模型實驗的局限性，數學模型建立和計算仍然是開展植被消浪研究的重要手段[6]。相比于相位平均類模型(在波浪周期上進行平均)，相位識別類波浪模型能提供更加精細的水動力要素而被廣泛應用，在此基礎上考慮植被對波浪的作用即形成適用于植被覆蓋區域波浪傳播的數學模型。現有的此類數學模型大致可以分為三類：一是基于雷諾平均的納維—斯托克斯方程的數值模型[7]；二是淺水方程模型，主要包括完全非線性淺水方程數學模型[8]和Boussinesq方程數值模型[9-10]；三是基于緩坡方程的數值模型[11]。納維—斯托克斯方程模型雖然可以很詳細地描述流場，但是其計算代價昂貴，尤其是對于水平二維問題，因此適用性有限。完全非線性淺水方程模型相對較簡單，計算效率高，但是這類模型僅能模擬短距離波的傳播，不能考慮長距離傳播和繞射等短波效應，尤其在處理近岸植被區域波浪的傳播變形問題時，精確度不高。緩坡方程模型僅適用于坡度變化較緩情況，對于波浪非線性、反射等因素考慮不足，不適用于近岸地形復雜和波浪非線性強的情況。而Boussinesq方程兼具色散型和非線性，使其不但能考慮包括波浪傳播過程中的反射、折射、繞射等因素，擴展后還可以包括波浪破碎及海岸動邊界問題，在模擬近岸波浪的傳播變形問題中顯示出巨大的優勢[9-10]。在這類模型基礎上建立考慮植被影響的數學模型是不錯的選擇，但目前基于Boussinesq方程建立的計算模型很少，而且主要集中在一維問題[9-10]。

綜上所述，基于建立的水平二維、具備間斷捕捉格式的Boussinesq波浪模型，在守恒形式的源項中引入由于剛性植被作用引起的拖曳阻力項，來考慮剛性植被區域波浪的傳播變形，進行相關計算，并分別選取淹沒、非淹沒植被，以及植被覆蓋水平二維復雜地形上的波浪傳播變形算例進行驗證。

1 數學模型

1.1守恒形式的二維Boussinesq水波方程

圖1 坐標系統及變量定義Fig. 1 The sketch of cartesian coordinate and variable definition

在如圖1所示的笛卡爾坐標系下，守恒形式的二維Boussinesq 水波方程[12-13]：

式中：下標x和y分別表示變量對x和y的偏導數，t表示變量對時間的偏導數。U為變量矢量，F(U)和G(U)分別代表x方向和y方向的通量矢量，S(U)表示源項，分別定義為：

其中，η為波面升高，h為靜水深，d=η+h為當地總水深，g為重力加速度。P和Q分別為x和y方向的斷面流量，即P=du，Q=dv(u和v分別為x和y方向的水深平均速度)。式(2)中U和V定義如下：

式(5)中Sb為海床坡度項，Sf為水底摩擦項，Sd為色散項，Sveg為植被作用引起的阻力項，其中摩擦項τx和τy，可表示為[14]：

底摩擦系數f一般在0.001～0.01范圍內，文中f為可調參數，具體值由相關的實驗測量結果確定。φx和φy分別是x和y方向的色散項[12]，定義為：

式中：B=1/15為色散性參數，該取值使控制方程適用于中等水深范圍，即h0/L0<0.5(h0為特征波高，L0為特征波長)。τxveg和τyveg為由植被作用引起的阻力，具體求解公式見下節。

1.2植被阻力項

由于植被的存在，波浪經過植被區域時傳播形態會發生變化。假定植被為剛性植被，則可忽略植被本身的運動，波浪經過植被區域受到的作用力主要有慣性力和拖曳阻力[15]。而植被的直徑相對于波高來說很小，沒有繞射影響，僅有黏性作用，因此忽略慣性力，只考慮拖曳阻力。將剛性植被簡化為小圓柱體進行受力分析，植被引起的拖曳阻力可用莫里森公式表示[15]，即：

式中：Cd是拖曳力系數，同流體特征以及植被特征密切相關，對數值結果有決定性影響，但迄今為止其取值沒有統一認識，本文將其作為可調參數，通過實驗結果進行率定，具體取值在算例中給出。Nv為區域中植被的密度，即植被均勻分布時單位面積的植株數，Av表示植被的迎水面積，定義為[9]：

法比因為將就槍傷的疼痛，僵著半邊身體站在她對面。她對他講這么多，讓他有點尷尬，有點愧對不敢當，他又不是她的懺悔神甫，她也不是懺悔的教徒。對于常常獨處的法比，把過多地了解他人底細看成負擔，讓他不適。或許叫玉墨的這個女人在做某種不祥的準備。

式中：bv和hv分別為植被的迎水寬度和植被的高度，bv取單株植被的直徑Dv。

圖2 出水和淹沒植被的布置示意Fig. 2 Sketch of vegetation elements

圖3 植被平面布置示意Fig. 3 Layout of vegetation

考慮到植被可以以出水和淹沒兩種方式出現(如圖2所示)，式(11)中的水平速度up確定如下[10]：

式中：Φ=Vs/V為植被自身與植被區域的體積比[10]。由上式可知，當植被出水時，up=u；當植被淹沒時，up為等效速度。Vs的計算如圖3所示，Vs=Nvπbv2hv/4為植被自身所占的體積，V為植被區域所占的空間，即等于水槽寬度、水深與植被區域長度的乘積。同理可得y方向的速度vp。

1.3數值格式和邊界條件

在矩形網格系統上，采用具有總變差減小TVD(total variation diminishing)性質的有限體積-有限差分混合格式求解Boussinesq水波方程，利用高分辨率有限體積方法求解對流項而其他項仍然通過有限差分方法計算。由于充分利用了有限體積方法的守恒特性，模型克服了傳統的、基于有限差分方法求解的同類模型存在的諸多缺陷，在穩定性、波浪破碎以及海岸動邊界處理方面具有極大優勢[13-14]。數值通量的計算采用兼具中心格式簡單性和迎風格式精確性的新型數值格式——MUSTA格式[13]，時間步長滿足CFL穩定性條件限制。

計算域兩端設置為固壁邊界，模型采用在質量方程中增加源項的方法產生波浪，同時視需要在計算域末端設置海綿吸收層吸收波浪[14]。當波面升高同水深比達到0.8時，認為波浪發生破碎，控制方程中的Boussinesq高階非線性項和色散性不參與計算，方程退化為淺水非線性方程。對于海岸動邊界，采用薄層水體法[13]處理, 即給定一極小水深值dmin，這里取0.000 1 m，若計算單元水深di,j>dmin(下標i和j分別為網格單元x和y方向的標記)，判定單元為濕單元，反之，若di,j

2 模型驗證

為了驗證所建立計算模型的精度和適用性，針對幾個典型的植被區域波浪傳播的物理模型實驗進行數值模擬。

2.1淹沒植被區域波浪的傳播

該植被消浪實驗是在密西西比國家泥沙實驗室進行的[9]，實驗布置如圖4所示，實驗水槽尺寸為20.6 m×0.5 m×1.22 m(長×寬×高)。實驗中用固定在床面的木質小圓柱體代表剛性植被，植被區域長Lv=3.66 m，植被密度為Nv=350 m-2，植被的高度和直徑分別為hv=0.63 m和Dv=0.009 525 m。實驗時取靜水深h0為0.7 m，五個浪高儀分別布置在3 m、11 m、12.5 m、14 m和15.5 m處，水槽左側為造波板，右側為消波設備。實驗所用波浪為規則波，對應的波要素如表1所示。

圖4 波浪在植被區域傳播變形的實驗布置示意Fig. 4 Experimental sketch of the wave propagation in the vegetation area

工況波高/m周期/sh0/L0一0．0669411．20．322428二0．0658521．40．249571三0．0653271．60．204288四0．0679611．40．249571五0．0851791．60．204288六0．0809801．20．322428

數值水槽布置同物理模型實驗水槽基本一致，如圖5所示。計算域水槽長30 m，寬1.0 m。植被區域長度仍設為3.66 m，植被的密度、高度和直徑均與實驗設置一致。造波源距離左側邊界5 m，計算域的左右邊界分別設置4 m寬的海綿層進行消浪。數值模擬的空間和時間步長分別取為Δx=Δy= 0.05 m和Δt=0.02 s，總模擬時間為60 s。拖曳力系數Cd經實驗數據率定取值為3.5，底摩擦系數f取值0.005。由于不確定表1中波浪要素是否由G1浪高儀采集，計算時波浪要素的確定以波高計算結果同G2處實驗數據吻合為準。

數值計算結果在圖6中給出，其中x軸零點取在G2浪高儀位置處。如圖6所示，當波浪在植被區域傳播時，隨著傳播距離的增加，波高發生了不同程度的衰減，對于表1中6組不同的工況，本文模型的計算結果均和實驗數據吻合良好，表明模型不但有足夠的計算精度，也具有較好的適用性。對于該組物理模型實驗，Kuiry 等[9]也采用Boussinesq類模型進行過數值模擬，計算結果也一并在圖6中給出，可見其和本文計算結果基本一致，但存在一些細節差別。在G2處，Kuiry等的計算結果均偏低于實驗數據，而本文計算時以G2處波高確定入射波浪要素，計算結果和實驗數據吻合更好。對于所有工況，Kuiry等的計算結果均顯示波浪在植被區域傳播時會發生反射現象，入射波浪和反射波浪疊加導致波高產生波動，呈駐波形態。對于工況二～五，本文模型的計算結果也揭示同樣的現象，但對于周期較小的工況一和六并不明顯，這可能是由于計算時參數取值或數值水槽布置不同引起的。

圖5 波浪在植被區域傳播變形的數值模擬布置示意Fig. 5 Numerical sketch of the wave propagation in the vegetation area

圖6 規則波經過淹沒植被區域波高的變化Fig. 6 Variation of wave heights along the submerged vegetation area

2.2出水植被區域波浪的傳播

仍然利用圖5所示的數值水槽，取靜水深h0為 0.5 m，植被區域的長度和位置不變，植被的高度hv為0.63 m，植被的直徑Dv為0.009 4 m，拖曳力系數Cd經實驗數據率定取值為4.5，底摩擦系數f取值為0.005。數值模擬的時間和空間步長分別取為Δx=Δy=0.05 m和Δt=0.02 s，總的模擬時間為60 s。選取兩組不同初始波浪條件下植被密度值分別為350 m-2和623 m-2的工況[16]，模擬出水情況下植被區域波浪的傳播變形。

將本文Boussinesq模型的計算結果、Wu等[16]采用RANS類模型的模擬結果以及實驗數據進行了對比，數值結果如圖7、8所示。結果表明，與淹沒情況類似，波高呈現下降趨勢，而且相比植被淹沒情況，出水情況下阻力較大，波高曲線變化較快，衰減劇烈。由圖7和圖8的分別對比發現，當入射波的波浪要素一定時，隨著植被密度的增大，波高衰減速度加快，衰減幅度增大。本文的Boussinesq模型和Wu等的數值模擬在x=0處的數值結果和實驗數據一致，兩種數值模擬的曲線走勢差別不大。圖7、8中，整個數值模擬的結果和浪高儀處的實驗結果吻合依然較好，由此可見模型同樣適用于處理出水條件下波浪傳播問題，適用性強。

圖7 波高H=0.053 m、周期T=0.7 s的規則波經過出水植被區域波高的變化Fig. 7 Variation of wave heights along the emergent vegetation area for regular wave H=0.053 m and T=0.7 s

圖8 波高H=0.10 m、周期T=1.0 s的規則波經過出水植被區域波高的變化Fig. 8 Variation of wave heights along the emergent vegetation area for regular wave H=0.10 m and T=1.0 s

2.3水平二維植被區域波浪的傳播

Thuy等進行了一系列關于波浪在部分植被覆蓋的斜坡上傳播的實驗[8, 17]，具體的實驗布置如9所示，實驗水槽長15 m，寬4 m，靜水水深0.44 m，地形的斜坡有不同的坡度，其中主要的坡度為1∶20.5。植被區域長1 m，沿x方向分布在10.36 m到11.36 m之間，植被區域的寬度為0.33 m，沿y方向分布在 0.07 m到0.40 m之間。溝槽的中心線在y=0.035 m處，植被區域的中心線在y=0.235 m處，在G6浪高儀處的溝槽和植被區域的中心線分別放置兩個流速儀，植被區域和流速儀位置的平面布置如圖10所示。具有水平二維特征的植被消浪物理模型實驗非常少見，本文將其用于水平二維模型的驗證。

數值模擬時，取計算域長25 m，寬0.4 m，造波板距離左側邊界5 m，計算域的左右邊界分別設置4 m和5 m的海綿層進行消浪。初始狀態的靜水深為0.44 m，入射波的波高為0.032 m，周期為20 s。植被區域長度為1 m，與實驗布置位置一致，植被的高度hv為0.5 m，植被的直徑Dv為0.005 m，迎水寬度bv為0.005 m，植被的密度為Nv=2 200 m-2，數值模擬的空間步長分別為Δx=0.05 m和Δy=0.01 m，時間步長為Δt=0.02 s，總模擬時間為200 s。拖曳力系數Cd經實驗數據率定取值為1.5，底摩擦系數f取值為0.005。

圖9 Thuy等的實驗布置示意Fig. 9 Experimental diagram of Thuy et al.

圖10 植被區域和流速儀位置布置示意Fig. 10 The layout of the vegetation zone and the location of the velocity instrument

當波浪場穩定之后，計算得到不同流速儀處沿x方向的速度u的時間歷程曲線，如圖11所示。可以看出，當水流穩定后，兩個位置處的速度時間歷程曲線呈周期性變化，溝槽中心處流速的最大值約為40 cm/s，而植被區域中心處的最大值不足20 cm/s，顯然，植被的阻力作用導致了植被區域的速度減小，消浪效果明顯。數值模擬的速度歷程曲線與實驗結果的數據點吻合較好，進一步體現了模型處理水平二維復雜地形下部分植被覆蓋區域波浪傳播變形問題的優越性能。當植被區域寬度為0.4 m，即整個斷面均布置有植被時，計算得植被區域的波高、波峰和波谷的比較情況，如圖12所示，從圖中可以看出，隨著波浪的傳播，水深逐漸變淺，波峰增大，波谷減小，波浪的非線性增強，在這個過程中，數值模擬與實驗結果相比，無論是波峰、波谷還是波高，吻合依然很好，說明了模型處理相關水動力問題的能力，體現了模型的計算精度。

圖11 速度的時間歷程曲線Fig. 11 Time series of velocity

圖12 整個斷面均布置植被時的波高變化曲線Fig. 12 Variation of wave heights, crests and troughs through the cross section full of vegetation

3 結 語

基于已經建立的水平二維、具備間斷捕捉能力的Boussinesq類波浪模型，建立了剛性植被區域波浪傳播變形的計算模型，通過選取植被出水、淹沒情況下不同入射波波高和波浪周期的算例進行模擬。數值結果表明，波浪在近岸植被區域傳播時，波高發生不同程度的衰減，且波高衰減程度隨著植被高度和植被密度的增加而增加。同時，選取水平二維部分植被覆蓋的斜坡上波浪傳播變形的算例進行模擬，波高和流速的計算結果同實驗數據吻合良好，證明了所建立的數值模型對于復雜地形條件、剛性植被影響下的色散性波浪傳播變形過程具有較高的計算精度，適用于開展剛性植被對波浪衰減作用的數值計算研究。

[1] 吉紅香,黃本勝,邱秀云. 植被消波消浪研究綜述[J]. 水利水運工程學報, 2005(1): 75-78.(JI Hongxiang, HUANG Bensheng, QIU Xiuyun. Summary of wave dissipation by the vegetation[J]. Hydro-Science and Engineering, 2005(1): 75-78.(in Chinese))

[2] 白玉川,楊建民,胡嵋,等. 植被消浪護岸模型實驗研究[J]. 海洋工程, 2005, 23(3): 65-69.(BAI Yuchuan,YANG Jian min, HU Mei,et al. Model test of vegetation on the bank to attenuate waves and protect embankments[J]. The Ocean Engineering, 2005,23(3): 65-69.(in Chinese))

[3] 楊建民. 海岸帶邊坡防浪林消浪理論與實驗研究[J]. 海洋通報, 2008, 27(2):16-21. (YANG Jianmin. Analytical solution and physical model test of wave protection forest on coastal zones[J].Marine Science Bulletin, 2008, 27(2):16-21.(in Chinese))

[4] 王瑞雪.非淹沒剛性植物對波浪傳播變形影響實驗研究[D]．長沙：長沙理工大學, 2012. (WANG Ruixue. Laboratory investigation on wave transformation through the emergent, rigid vegetation[D]. Changsha: Changsha University of Science & Technology, 2012.(in Chinese))

[5] 耿寶磊,高峰,張宇亭. 港口導堤植被消浪的試驗研究[J]. 水運工程, 2014(7):51-57.(GENG Baolei, GAO Feng, ZHANG Yuhao. Experimental study on wave absorbing effect with vegetation on port jetty [J]. Port & Waterway Engineering, 2014(7):51-57.(in Chinese))

[6] MEI C C ,CHAN I C, LIU P L F,et al. Long waves through emergent costal vegetation[J]. Journal of Fluid Mechanics, 2011, 687:461-491.

[7] NAOT D, NEZU I, NAKAGAWA H. Hydrodynamic behavior of partly vegetated open-channels[J]. Journal of Hydraulic Engineering, 1996, 122(11): 625-633.

[8] THUY N B, TANIMOTO K, TANAKA N. Flow and potential force due to runup tsunami around a coastal forest with a gap-experiments and numerical simulations[J]. Science of Tsunami Hazards, 2010, 29(2): 43-69.

[9] KUIRY S N, WU W, DING Y. A hybrid finite-volume/finite-difference scheme for one-dimensional Boussinesq equations to simulate wave attenuation due to vegetation[C]//Proc. of the World Environment &Water Resource Congress.2011: 2114-2124.

[10] HUANG Z, YAO Y, SIM S Y ,et al. Interaction of solitary waves with emergent, rigid vegetation[J]. Ocean Engineering, 2011, 38(10): 1080-1088.

[11] 唐軍,沈永明,崔雷. 基于拋物型緩坡方程模擬近岸植被區波浪傳播[J]. 海洋學報, 2011, 33(1):7-11.(TANG Jun, SHEN Yongming, CUI Lei. Modeling costal water wave propagation in vegetation field based on parabolic mild slope equation[J]. Acta Oceanologica Sinica, 2011, 33 (1): 7-11. (in Chinese))

[12] MADSEN P A , SORENSEN S R. A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[J]. Coastal Engineering, 1992, 18(3-4): 183-204.

[13] FANG K Z, ZOU Z L, DONG P,et al. An efficient shock capturing algorithm to the extended Boussinesq wave equations[J]. Applied Ocean Research, 2013, 43:11-20.

[14] 房克照,鄒志利,孫家文,等. 擴展型Boussinesq水波方程的混合求解格式[J]. 水科學進展, 2013,24(5): 699-705. (FANG Kezhao, ZOU Zhili, SUN Jiawen ,et al. A hybrid numerical scheme for the extended Boussinesq equations[J]. Advances in Water Science, 2013, 24(5): 699-705. (in Chinese))

[15] MENDEZ F J, LOSADA I J. An empirical model to estimate the propagation of random breaking and nonbreaking waves over vegetation fields [J]. Coastal Engineering, 2004, 51(2): 103-118.

[16] WU W, ZHANG M, OZEREN Y, et al. Analysis of vegetation effect on waves using a vertical 2D RANS model[J]. Journal of Coastal Research, 2013, 29(2): 383-397.

[17] THUY N B, TANIMOTO K, TANAKA N,et al. Effect of open gap in coastal forest on tsunami run-up- investigations by experiment and numerical simulation[J]. Ocean Engineering, 2009, 36(15-16): 1258-1269.

A numerical model for coastal wave propagation in the rigid vegetation area

WANG Lei1, FANG Kezhao1, YIN Jing2, SUN Jiawen2

(1. The State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China; 2. National Marine Environmental Monitoring Center, State Oceanic Administration, Dalian 116023, China)

A numerical model for coastal wave propagation through the rigid vegetation area is developed based on the extended Boussinesq equations. Drag resistance force is introduced into the source term of momentum equations to consider rigid vegetation effect on wave attenuation. A hybrid finite-difference and finite-volume method is used to solve the governing equations. The resulting model has the advantages of strong-stability-preserving and shock-capturing and can effectively simulate wave transformation, breaking and moving shoreline in the coastal area. Typical numerical tests are conducted for model validation and the computed results are in good agreement with the experimental data, demonstrating the proposed model is applicable of predicting wave transformation in the coastal region covered with rigid vegetation.

Boussinesq equations; rigid vegetation; wave attenuation; numerical simulation;wave propagation

TV139.2; O353.2

A

10.16483/j.issn.1005-9865.2015.06.008

1005-9865(2015)06-062-08

2015-02-03

國家創新研究團隊資助項目(51221961)；國家海洋局海域管理技術重點實驗室資助項目(201206)

王 磊(1988-)，女，山東煙臺人，博士生，從事海岸動力學方面的研究。E-mail: wangleiDUT@163.com