木本植被覆蓋岸坡上波浪爬升過程的數(shù)值模擬研究

徐海玨，胡萍，白玉川*，楊波

(1.天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300350；2.天津大學(xué) 河流海岸工程泥沙研究所，天津 300350；3.AECOM(天津)工程顧問有限公司，天津 300072)

1 引言

近年來，海平面上升、海岸線侵蝕等問題對于沿海地區(qū)產(chǎn)生了不利影響[1]。伴隨著沿海住宅和生態(tài)景觀的建設(shè)開發(fā)，人們愈發(fā)強(qiáng)調(diào)海岸防護(hù)措施的功能性和生態(tài)友好性。近岸木本植物構(gòu)成的生態(tài)緩沖帶作為一種新型的生態(tài)防護(hù)結(jié)構(gòu)，不僅能建立起近岸的綠色生態(tài)系統(tǒng)持續(xù)改善環(huán)境，還可以顯著地消減波能，減少波浪對近岸結(jié)構(gòu)物的影響，更能與傳統(tǒng)的近海工程防護(hù)措施共同形成工程植物防護(hù)體系[2-3]?；谶@些優(yōu)勢，利用木本植物形成的植被緩沖區(qū)的生態(tài)軟防護(hù)結(jié)構(gòu)正在逐漸取代傳統(tǒng)的硬質(zhì)海岸保護(hù)結(jié)構(gòu)。因此，分析植被和波浪因素對消浪效果的影響，深刻理解波浪－植物之間相互作用的機(jī)理，如何準(zhǔn)確預(yù)測波浪在植被影響下的傳播過程，提升其對海岸的防護(hù)效果是當(dāng)前亟待解決的問題。

近年來，針對植被消浪的研究日漸深入，且學(xué)者多采用理論推導(dǎo)、試驗(yàn)分析和現(xiàn)場調(diào)研3種方式。相對于理論推導(dǎo)，物理縮比試驗(yàn)和現(xiàn)場觀測研究可給出針對具體工程的波浪變形現(xiàn)象的定性描述。在物理模型試驗(yàn)方面，白玉川等[4]選取檜柏樹枝作為防浪樹的試驗(yàn)?zāi)Ｐ?，探究了防浪林種植寬度、底坡高度、底坡高程對消浪效果的影響；吉紅香[5]用剛性樹干和帶有柔性枝葉的模型樹模擬堤外灘地防浪林，分析了灘地水深、浪高、波長等因素對植物消浪的影響，以及波高、波浪爬高等的分布規(guī)律；周悅等[6]進(jìn)行剛?cè)峤M合型植被消浪特性的物理試驗(yàn)研究，探究了防浪林的樹木大小、剛?cè)岬炔煌M合方式對消浪效果的影響；Yao等[7]通過PVC直管模擬植被冠層，研究推導(dǎo)剛性植被阻力系數(shù)的方法。在現(xiàn)場觀測方面，Mazda等[8]和Quartel等[9]通過對波浪入射條件和傳播過程中的波浪變形的分析，分別討論了越南Tong King三角洲和紅河三角洲造林對海岸保護(hù)的作用；曹大正等[10]、李志強(qiáng)等[11]和田野[12]基于現(xiàn)場的波浪和泥沙的動態(tài)觀測，分別研究了互米花草、紅樹林、無瓣海桑等近岸植物的消浪促淤效果。盡管有很多學(xué)者做出了創(chuàng)新的成果，模型試驗(yàn)和實(shí)地調(diào)研仍有其各自的局限性：實(shí)驗(yàn)室模型相似性選擇的確定的合理性有待進(jìn)一步討論，且受到很多實(shí)驗(yàn)室客觀條件（動力、側(cè)限、波浪反射）的影響；實(shí)地調(diào)研需要大量的人力物力，還會對環(huán)境產(chǎn)生一定程度的影響。由于波浪爬坡過程是一較為復(fù)雜的非線性過程，其演變過程受多方面因素的影響，因此有必要通過理論推導(dǎo)和數(shù)值模擬對其進(jìn)行討論，并給出沿程波面和流場變化的細(xì)致的描述。在理論推導(dǎo)方面，Bai等[13-14]、李紹武等[15]從N-S方程、大渦模擬等理論出發(fā)，建立了波浪淺水傳播變形模型；Dalrymple等[16]基于線性波和剛性圓柱體空間分布的摩擦模型對局部波浪耗散進(jìn)行了計(jì)算；Fernando和Losada[17]基于理論推導(dǎo)和經(jīng)驗(yàn)系數(shù)相結(jié)合的方式，用剛性直桿代替植被，給出了可變深度植被場作用下波浪傳播模型；楊建民[18]引入反映剛性植物樹干拖曳力作用的水流阻力系數(shù)λ，得到波動勢函數(shù)的一階和二階近似解；唐軍等[19]基于考慮波浪折射、繞射等綜合效應(yīng)的拋物型緩坡方程，建立了近岸植被影響下波浪傳播運(yùn)動的數(shù)學(xué)模型；Maza等[20]基于IHFOAM的三維數(shù)值方法來研究海嘯波與紅樹林的相互作用，并用拖曳系數(shù)的宏觀方法計(jì)算波浪力；Zhu和Chen[21]討論了由植被阻力引起的耗散和近共振三位一體的相互作用對單個(gè)諧波衰減的影響。但目前的理論描述多將植被簡化為剛性直桿，缺乏考慮植被枝干的分區(qū)作用和植被密度、樹干高度和樹枝角度等植被特性對波浪的影響，且對在傾斜底坡這一條件下的植被的消浪效果缺乏描述和分析。因此，本文參考了以海南島文昌海岸近岸的植被（圖1）為代表的近岸木本植物形態(tài)，建立了木本植被保護(hù)下波浪沿斜坡爬升的表面波衰減理論模型。針對植被上、下層對波浪不同的作用形態(tài)，分別引入相應(yīng)的阻力系數(shù)以反映植被樹枝和樹干對波浪的拖曳力作用，并采用有限差分法離散方程和標(biāo)記網(wǎng)絡(luò)法（Marker and Cell Method，MAC）來跟蹤自由曲面上的水顆粒軌跡。以波浪沿坡度為1/30的斜坡爬升為算例，分別討論了有無植被作用下沿程波形和流場的變化，并分析了植物枝干的直徑、高度、密度、角度等植被特性和波浪因素對消浪效果的影響，得到植被消浪的基本規(guī)律。

2 理論模型

理論模型將波浪簡化為垂直平面上的二維波，得到波浪的連續(xù)性方程和修正的Navier-Stokes（N-S）方程，并引入阻力分量來反映植被枝干對波浪的作用。

2.1 控制方程

參照圖2所示的坐標(biāo)系，水平x坐標(biāo)以向右作為其正方向，垂直z坐標(biāo)以向上作為其正方向，以i=tanθ為底坡斜率，θ為坡角?？紤]單向波通過木本植被緩沖區(qū)傳播到岸邊的情況下，波浪會同時(shí)受到植被和傾斜底床的作用，木本植被海岸示意圖如圖3所示。

在實(shí)際的海岸通常有一個(gè)從深海到海岸的斜坡，數(shù)值模型的分析區(qū)域包含了一段平坦的河床和一段傾斜的河床。以線性減小水深為例，水深公式為

模型在水平方向上可分為3個(gè)區(qū)域討論。區(qū)域I（-∞＜x≤x0，-d0≤z≤0）是在有植物作用之前的平坦河床部分，這個(gè)區(qū)域的入射波在平底上傳播且區(qū)域中的反射波主要來自植物反射。區(qū)域II（x0＜x≤x1，-d(x)≤z≤0）則被有著不同密度和直徑的樹干的木本植物所覆蓋。波浪從區(qū)域I傳播通過植物到達(dá)區(qū)域III，而有一些反射波從區(qū)域III開始傳播通過植物到達(dá)區(qū)域I，由于植物樹干與波浪之間的相互作用，波能被大量耗散，即波浪衰減。區(qū)域 III（x1＜x≤xb，-d(x)≤z≤0）是在植物之后的區(qū)域，入射波的動能已經(jīng)大大減弱。有時(shí)，流速降至遠(yuǎn)低于沉積物起動速度的值，以致于大多數(shù)泥沙顆粒無法沿著斜坡移動。

圖1 海南島文昌海岸近岸木本植被緩沖帶Fig.1 Nearshore woody vegetation at Wenchang coast in Hainan Island

圖2 無植被海岸示意圖Fig.2 Sketch of beach without vegetation

圖3 木本植被護(hù)岸示意圖Fig.3 Sketch of beach with woody vegetation

一般來說，在沒有植被的情況下，破波點(diǎn)存在于在區(qū)域II或區(qū)域III中的某個(gè)地方。但由于植被的存在，區(qū)域II的消波率發(fā)生了很大的變化，引起流場的變化，并最終導(dǎo)致了破波點(diǎn)位置變化。在這3個(gè)區(qū)域上，波浪可簡化為二維波并提出了考慮波浪傳播整體效果的連續(xù)介質(zhì)等效模型，得到對應(yīng)的連續(xù)性方程和修正的N-S方程：

式中，u、w分別為x、z方向的速度分量；t為時(shí)間；ρ為水的密度；p為壓力；ν為運(yùn)動黏滯系數(shù)；νt為動態(tài)渦流黏度；g為重力加速度；Fx和Fz分別是植物樹干對波浪作用力在x方向和z方向的分力。

Massel等[22]已經(jīng)證明了植物樹干對流體的作用力主要為阻力，因此，F(xiàn)x和Fz可以直接寫成阻力分量，而Fx和Fz在不同區(qū)域的表達(dá)是不同的，具體如下：區(qū)域I：如果流體滿足非黏性、不可壓縮假定，則表達(dá)式為

區(qū)域II：能量耗散主要是由于流體質(zhì)點(diǎn)與底部的相互作用，以及在沒有植物的傾斜底坡上攀爬的波浪的流體質(zhì)點(diǎn)之間的相互作用。如果在流場中增加了木本植被作用，那么枝干與流體之間的相互作用也就增加了，具體表現(xiàn)為枝干的阻力。由于前兩種情況已經(jīng)在許多研究中被頻繁地討論，本研究主要聚焦于最后一種情況，即增加的枝干的阻力。此時(shí)，仍假設(shè)流體滿足非黏性、不可壓縮假定，但植被的阻力不可忽略，則表達(dá)式為

首先，考慮單個(gè)小區(qū)域的木本植物作用，并將單位面積植株的數(shù)量記為N。然而，即使對同一種樹，枝干的密度和直徑在上層或下層是不同的。一般的木本植物上層有著濃密且直徑較小的樹枝，而下層為直徑較大的樹干。因此，我們在該模型中分區(qū)引入了不同的阻力表達(dá)式來區(qū)分樹干和樹枝的作用。

在植物下層 [-d(x)≤z＜-(d(x)-d1)]，單位面積上樹干的數(shù)量NL明顯等于植株的數(shù)量N，記d1為對應(yīng)植物平均樹干高度（hL）的下層水深。假定第i個(gè)植株的樹干直徑為DLi，那么拖曳力表示為

式中，F(xiàn)lx和Flz分別是下層植物樹干對波浪作用力在x方向和z方向的分力；Cdmi是第i棵植株的拖曳力系數(shù)；DLi是第i棵植株的樹干的的直徑。

在植物上層[-(d(x)-d1)≤z≤0]，單位面積上樹枝的數(shù)量NU要大于植株的數(shù)量N。此外，并不是所有的樹枝都是垂直于靜水面的，其阻力表達(dá)式定義為

式中，F(xiàn)ux和Fuz分別是上層植物樹枝對波浪作用力在x方向和z方向的分力；θi是樹枝與垂直方向的夾角大小，一般滿足正態(tài)分布；DUi是第i棵植株的樹枝的直徑。

根據(jù)Massel等[22]的假設(shè)，將這些角度均用平均角度θm來表示時(shí)，理論結(jié)果與實(shí)驗(yàn)結(jié)果相差不大。因此，為簡化模型，用平均角度θm代替各個(gè)樹枝與垂直方向的夾角大小θi，得到簡化后的植物上層阻力表達(dá)式定義，

如上文所述，Cdmi是一個(gè)修正系數(shù)，它不僅包括樹干和樹枝對流體的阻力，還包括不同圓柱形樹干和樹枝后面的渦旋之間的相互作用力。假設(shè)植被的生長相對較為稀疏或樹干和樹枝對流體的阻力起主導(dǎo)作用，則可以在計(jì)算中忽略渦旋之間的相互作用，那么對單個(gè)樹枝而言系數(shù)Cdi與雷諾數(shù)Re相關(guān)。根據(jù)SPM[23]得到的結(jié)果，Cdi的表達(dá)式可以寫成

式中，Re為雷諾數(shù)；Cdi是單棵植株的拖曳力系數(shù)。

引入以下影響因素：

式中，I(Re,N)是由于植株密度改變引起的影響因素，它的值是由另一個(gè)小范圍數(shù)值模型的結(jié)果決定的。

區(qū)域III：這一區(qū)域內(nèi)沒有樹的作用，即不存在流體和樹枝之間的相互作用，且經(jīng)過植被的消能作用后，在此區(qū)域的波浪形態(tài)表現(xiàn)為前進(jìn)波和破碎波。因此，在區(qū)域III中流體不應(yīng)被視為非黏性的，而應(yīng)視為黏性的旋轉(zhuǎn)湍流，

式中，νt是動態(tài)渦流黏度，νt=2C(x,z,t)Δ2|S|。其中，|S|為應(yīng)變率的絕對值，其公式表示為|S|=[2(Sij·Sij)]1/2，Sij=[(?ui/?xj)+(?uj/?xi)]/2，Sij可表示為Sxx，Sxz，Szx，Szz，i=x,z且j=x,z；Δ為與G相關(guān)的特征過濾值；C為與湍流結(jié)構(gòu)密切相關(guān)的參數(shù)；G為網(wǎng)格過濾值。

2.2 邊界條件

自由表面邊界介于剪切層和無剪切環(huán)流間，通常假設(shè)外部流完全不受湍流影響，因此在這個(gè)邊界上，所有的湍流應(yīng)力和通量都可以設(shè)為0。底邊界按可滑移條件處理，且滿足垂直于邊界的速度為0和邊界網(wǎng)格內(nèi)的速度散度為0兩個(gè)條件。底坡與水面相交點(diǎn)附近的網(wǎng)格既含有底邊界又含有自由表面其處理稍為復(fù)雜，但因近岸岸灘和本文算例的底坡均很緩，則含有底邊界和自由面的網(wǎng)格的壓力可近似等于靜水壓力，以簡化這些網(wǎng)格的處理。

3 數(shù)值求解方法

本文的數(shù)值求解方法采用有限差分法離散方程，并用MAC法來追蹤自由曲面上的水顆粒軌跡。每個(gè)物理變量的網(wǎng)格劃分和分布如圖4，其中dx、dy分別為x和z方向的空間步長。

圖4 MAC法單元圖Fig.4 MAC method element diagram

3.1 數(shù)值方法求解流場

MAC方法采用了網(wǎng)格交錯(cuò)技術(shù)，將壓力、湍流剪切應(yīng)力等變量定義于網(wǎng)格中心，將速度 (u,w)等變量定義于網(wǎng)格邊線中點(diǎn)，這樣便于滿足單元內(nèi)速度散度為0這一條件。在方程中，對時(shí)間采用向前差分，對空間采用中心差分。由動量方程導(dǎo)出的差分方程以及由連續(xù)性方程導(dǎo)出的差分方程即為MAC方法需求解的差分方程。

求解泊松方程采用松弛迭代法，其迭代關(guān)系式為

式中，上標(biāo)m和m+1分別表示第m和m+1次迭代；ω是松弛因子

其中，

式中，上標(biāo)m+1/2表示由m和m+1兩個(gè)迭代步導(dǎo)出的量。

3.2 施主格式的穩(wěn)定性分析

施主格式數(shù)值耗散比中心差分格式更小，為了滿足穩(wěn)定性條件以及減少差分格式的數(shù)值耗散，本研究一律采用施主格式。為得到方程的收斂解，需對NS方程導(dǎo)出的差分方程和壓力泊松方程的差分方程的數(shù)值穩(wěn)定性進(jìn)行分析。N-S方程為非線性方程，其穩(wěn)定性分析是相當(dāng)困難的。根據(jù)Von Neumann法分析，對于波浪運(yùn)動，直角坐標(biāo)系下二維問題的MAC方法的時(shí)間穩(wěn)定性條件可簡化為

泊松方程的系數(shù)矩陣是對稱正定矩陣時(shí)，取松弛因子 0 ＜ω＜2，松弛迭代是收斂的，且存在著使收斂速度最快的最優(yōu) ω*值，一般來說，取 ω ＞1的超松弛迭代比取 ω ＜1的欠松弛方程迭代收斂速度要快。但在NS方程和壓力泊松方程的耦合問題中，每一時(shí)間步泊松方程系數(shù)矩陣是不相同的，因此無法分析其收斂條件，也無法確定最優(yōu)的 ω*值。在實(shí)際計(jì)算中欠松弛迭代法較為適用，即求解壓力泊松方程時(shí)取前一時(shí)間步長的壓力值作為迭代的初始值，每一步時(shí)間增量很小，壓力變化也不大，因此可以保證泊松方程的迭代過程穩(wěn)定地趨于收斂。

4 數(shù)值求解算例及消浪效果影響規(guī)律分析

數(shù)值求解是根據(jù)給定的邊界條件，利用MAC方法對波浪的傳播和爬升過程進(jìn)行模擬。由于模型存在一定的局限性(恒定水深、常浪入射條件、植物均勻分布)，本文暫不進(jìn)行理論模型結(jié)果與實(shí)際木本植被緩沖帶觀測值的比較，這二者比較的結(jié)果將隨著研究的深入進(jìn)行后續(xù)討論。但為給相關(guān)的物理模型試驗(yàn)研究提供對照驗(yàn)證，也為后續(xù)物理試驗(yàn)的開展提供依據(jù)，本節(jié)將給出數(shù)值求解模型的算例結(jié)果，并對木本植物群的消浪效果影響規(guī)律進(jìn)行分析。

本文以波浪沿1/30的斜坡爬升為算例，通過給出入射波的基本參數(shù)（入射水深DS，入射波高HS，入射波長LS，入射波周期TS），分別討論有無植被作用下的波形變化和速度場情況。此外，在考慮木本植被作用的情況下，通過引入不同的阻力系數(shù)，討論了植物枝干的直徑、高度、密度等因素對波浪能量衰減的影響。此外，本節(jié)將得到的部分算例結(jié)果與相關(guān)文章中的試驗(yàn)結(jié)果規(guī)律進(jìn)行對照，驗(yàn)證了數(shù)值模型的有效性，可供實(shí)際參考。

4.1 有無植被作用下波浪沿傾斜海灘傳播的數(shù)值模擬研究

波浪沿斜堤或斜坡向上傳播時(shí)，不斷向上爬升必然會出現(xiàn)波形的不斷變化，本小節(jié)通過對比分析有無植被作用情況下的波浪變形和流場情況，討論近岸植被對波浪衰減的影響。

4.1.1 入射波浪條件

數(shù)值計(jì)算區(qū)域參照圖2所示，x=0 m為波浪入射邊界，給定計(jì)算水平長度xb=32 m，其中平坦河床段x0=2.0 m，傾斜河床段xb-x0=30 m，i=1/30。在初始時(shí)刻計(jì)算區(qū)域內(nèi)的水體是靜止的，即初始速度場為0 m/s，初始壓力場按靜水壓力場給出。以二階Stokes波為例，給定的各項(xiàng)波浪參數(shù)為：入射水深DS=1.0 m，入射波高HS=0.25 m，入射波周期TS=2.9 s，入射波長LS=8.5 m和波陡 δ =0.029。波浪通過入流邊界處的水質(zhì)點(diǎn)速度和波面位置按相對應(yīng)的理論公式變化而產(chǎn)生。根據(jù)時(shí)間穩(wěn)定性條件，網(wǎng)格尺度取dx=0.125 m，dy=0.025 m，時(shí)間步長dt=0.01 s，網(wǎng)格的劃分和時(shí)間步長的選取綜合考慮了差分格式的穩(wěn)定性、盡可能小的數(shù)值耗散和合理的計(jì)算工作量。

4.1.2 無植被作用情況下波浪傳播數(shù)值模擬及模型合理性驗(yàn)證

通過調(diào)整入射波波前與靜水面相交點(diǎn)到波峰的距離，進(jìn)而改變?nèi)肷溥吔绲牟嫖恢?，得?條波浪在1/30的底坡上傳播的波形曲線（圖5），并以左邊界位置的波面高度HL=0.03 m為例繪其傳播流場（圖 6）。

通過波浪傳播的波形圖可以發(fā)現(xiàn)，在沒有植被的拖曳力作用時(shí)，破浪爬坡過程波高上升趨勢明顯，隨著水深變淺，波高增大，波長變短，波前峰變陡，而波峰后逐漸出現(xiàn)波面變化較小的臺階式爬升、繼而破碎現(xiàn)象，這種水平不對稱的變化與流函數(shù)波理論解是一致的。對應(yīng)的，以左邊界位置的波面高度HL=0.03 m為例生成的流場圖中，隨著波前峰趨于變陡波峰處的水平速度不斷增大，豎向速度大小變化加劇，并在最大水平速度達(dá)到波速時(shí)發(fā)生破碎，波浪沿斜坡傳播流場（圖6）不同位置處的局部放大圖如圖7所示。

4.1.3 有植被作用情況下波浪傳播數(shù)值模擬

圖5 無植被作用下波浪沿斜坡（i=1/30）傳播波形Fig.5 The wave form of the wave propagation along slope (i=1/30) without vegetation

圖6 無植被作用下波浪沿斜坡（i=1/30）傳播流場Fig.6 The flow field of the wave propagation along slope (i=1/30) without vegetation

圖7 無植被作用下流場局部放大圖（流場所處的位置見圖6）Fig.7 The partial enlarged drawing of the flow field without vegetation (the positon of the flow field see Fig.6)

考慮植被作用情況下，參數(shù)Cdmi設(shè)定為：Cdmi=0.7，為與無植被作用情況進(jìn)行初步對比，將單位長度內(nèi)的植被參數(shù)設(shè)為：樹干數(shù)NL=1，樹枝數(shù)NU=6，樹干高度hL=0.6 m，樹干直徑DL=0.2 m，樹枝直徑DU=0.04 m，樹枝傾斜度 θ=π/6，植被水平方向覆蓋寬度lT=27 m，得到波浪在植被均勻覆蓋的1/30的底坡上傳播的波形（圖10）和流場（圖11）。其中，圖11為布置了均勻分布的植被波浪傳播流場示意圖，用以反映植物隨著底坡高程變化對流場的影響。

在木本植被沿斜坡均勻覆蓋的情況下，波浪傳播至在水平中段區(qū)域波高仍有明顯的上升趨勢，但隨著植被拖曳力的持續(xù)作用，波浪傳播至x=16.0 m后有明顯的衰減趨勢，波浪雖也產(chǎn)生波長變短，波前峰變陡的水平不對稱現(xiàn)象，但波浪趨于平緩的向岸線方向爬升而不發(fā)生破碎。對應(yīng)的，同樣以左邊界位置的波面高度HL=0.03 m為例生成的流場圖中，相對于無植被作用情況，隨著波前峰近岸趨于變陡，波峰處的水平速度和豎向速度變化不劇烈，流場狀態(tài)趨于平穩(wěn)。

有無植被作用情況下，波浪沿斜坡向上傳播的不同位置處的波峰下沿垂線上的水平速度矢量分布如圖12所示。波浪在沿斜坡爬升過程中，雖在前半段水平速度均隨著爬坡而增大，但在木本植被的消浪作用下，波浪的波面高度趨于穩(wěn)定，波峰下沿垂線上的水平速度值相對較小，尤其在中后段a、b兩種情況的波峰垂線上的水平速度分布差異明顯：隨著植被消能作用的疊加，有植被作用下的波浪水平速度持續(xù)減小，趨于平緩的向坡頂爬升；而無植被消能作用下的波浪在近坡頂處的水平速度仍持續(xù)增大，波浪波動較大，趨于破碎。

圖8 瀕臨破碎波峰下水平速度垂向分布比較Fig.8 Comparison of vertical distribution of horizontal velocity under breaking wave crest field

圖9 二階Stokes波模擬結(jié)果驗(yàn)證Fig.9 Verification of second-order Stokes wave simulation results

圖10 植被作用下波浪沿斜坡（i=1/30）傳播波形Fig.10 The wave form diagram of the wave propagation along slope (i=1/30) with vegetation

圖11 植被作用下波浪沿斜坡（i=1/30）傳播流場Fig.11 The flow field diagram of the wave propagation along slope (i=1/30) with vegetation

圖12 無植被作用（a）和有植被作用（b）波峰下水平速度垂向分布對比Fig.12 Comparison of vertical distribution of horizontal velocity under wave crest field with vegetation (a) and without vegetation (b)

圖13 無植被作用（a）和有植被作用（b）波峰位置波面高度沿程變化Fig.13 The variation of wave surface height on wave crest along distance with vegetation (a) and without vegetation (b)

為更直觀的體現(xiàn)木本植被的消浪作用，給出了有無植被兩種情況下波浪的波面高度沿程變化（圖13），且考慮波浪爬升的水面變化和二階stokes波的波形特點(diǎn)，以沿程波峰位置處的波面高度代替波高進(jìn)行分析。由圖可見，植被作用下的波面位置變化的總體趨勢為前半段穩(wěn)定增高，而后半段穩(wěn)定減小，且趨于平緩的爬升而不發(fā)生破碎。

4.2 植被特性對消浪效果的影響規(guī)律分析

考慮植被作用情況下，為綜合考慮木本植被對波能衰減的影響機(jī)理，需考慮的木本植被特性為：植被密度（單位長度內(nèi)植株的數(shù)量）、樹枝密度（每棵植株的樹枝數(shù)）、樹干高度，樹干直徑、樹枝直徑和樹枝傾斜度等因素。其中，由于植物在栽種后的生長過程中其樹干直徑和樹枝直徑短期內(nèi)變化緩慢，因此本文不對其進(jìn)行討論，而通過分析沿程波峰位置處的波面位置變化，著重討論植被密度、樹枝密度、樹干高度和樹枝傾斜度對波能衰減的影響機(jī)理。本節(jié)給定的波浪參數(shù)同上節(jié)，表1為考慮不同的植被特性條件下的數(shù)值模擬工況。

4.2.1 植被密度對消浪效果的影響

圖14是對應(yīng)工況1～4的沿程波峰位置處的波面位置變化，x軸為波浪向岸傳播的水平方向距離，z軸為波峰位置處的波面位置，圖中對應(yīng)各個(gè)工況折線沿程的最后數(shù)據(jù)點(diǎn)為波浪傳播至破碎前的波峰位置處的波面高度值。

表1 植被特性模擬工況表Table 1 The arrangement of simulation condition about vegetation characteristics

圖14 植被密度與波面位置關(guān)系Fig.14 The relationship between vegetation density and wave surface position

植被作用下，工況1～4通過控制樹干高度、樹枝傾斜度為定值，即通過同比例改變樹干、樹枝的數(shù)目，實(shí)現(xiàn)植被分布密度的改變。樹干與樹枝的比例為1:6，單位長度內(nèi)的植株的變化范圍為 0 .1≤NL≤3，其中NL=0.1（工況1）的情況為接近無植被作用的極端情況。在這一范圍內(nèi)，工況1與工況2、3、4的沿程波面高度變化對比明顯，且沿程高度差不斷增大。在x∈[22, 24]這一區(qū)段內(nèi)，波面高度工況1與工況4的高差值達(dá)到了0.13 m。而工況2～4的波面高度均先呈現(xiàn)穩(wěn)定的、較小幅度的爬升，繼而在水平傳播方向上的中后段出現(xiàn)較為明顯的下降趨勢。產(chǎn)生這種現(xiàn)象可能主要由于以下兩個(gè)原因：一方面，在植被消浪和傾斜底坡兩種條件綜合影響的情況下，波浪存在沿程消能波高減小和斜坡爬升波高增大兩種趨勢，在計(jì)算區(qū)域的中前段，傾斜底坡對波浪爬升的影響作用較明顯，而在中后段，植被的沿程拖曳力對波能衰減的作用逐漸顯現(xiàn)，并出現(xiàn)較大幅度的波面高度的減小；另一方面，當(dāng)水位較低時(shí)，即波浪傳播至近岸水深較淺的區(qū)域時(shí)，此時(shí)波能較小，且樹干對波浪施加與其傳播方向一致的拖曳力，起主要消浪作用，使波能迅速衰減。

對比工況2、3、4的情況，工況3的植被密度為工況2的2倍，工況4的植被密度為工況2的3倍。受植被密度影響，工況3的沿程波面高度略小于工況2，且波浪破碎區(qū)域相同；工況4雖然植被密度較大，但在沿程x＜26 m的區(qū)段內(nèi)，工況4相對于工況2、3波面高度較高，而在后段波面高度減弱更加明顯，波浪更加平穩(wěn)地向近岸方向傳播而不發(fā)生破碎。工況4這種沿程波面高度的變化可能是由于，在植被覆蓋密度較大的情況下，植物對水體的拖曳力疊加，植被區(qū)域內(nèi)水體波動較為劇烈，導(dǎo)致了在x＜26 m的區(qū)段內(nèi)波面位置較高的現(xiàn)象，而與此同時(shí)，也增加了能量的損耗，使透過植被區(qū)后的波能大大減弱，對應(yīng)的，波面位置則出現(xiàn)明顯的減弱現(xiàn)象。由此可見，單純增加植被的覆蓋密度來增加其對波浪的衰減作用的方式并不合理，一方面經(jīng)濟(jì)性欠佳，另一方面此方式也可能會增加沿程水體的波動，而對生態(tài)、岸灘防護(hù)等起到不利的影響。因此，在本小節(jié)的4種工況中，工況2更符合自然植被的實(shí)際分布情況，也較為經(jīng)濟(jì)可行，綜合效果更好。

4.2.2 樹枝密度對消浪效果的影響

相對于工況2而言，工況5～7對樹枝密度進(jìn)行了調(diào)整。圖15是對應(yīng)工況2、5～7的沿程波峰位置處的波面位置變化，圖中坐標(biāo)軸和折線的表示意義與圖14一致。由圖15可以看出，在x∈[0, 12]這一區(qū)段內(nèi)，隨著樹枝密度的增大，沿程波面高度的變化不明顯，而后則出現(xiàn)波面高度的較大差異：樹枝密度越大，波面位置的變化越趨于平緩。在x∈[0, 12]這一區(qū)段內(nèi)，植被覆蓋段區(qū)域?yàn)?0 m，其長度約等于區(qū)域內(nèi)的波長，而植被對波浪消能的作用是在傳播中不斷疊加的，而在沿程第一段波長的長度范圍內(nèi)，樹枝的效果尚未體現(xiàn)。但隨著波浪不斷向前傳播，樹枝密度對波面位置的影響愈發(fā)顯著，其中工況2的沿程波面高度最為穩(wěn)定，且在x∈[18, 20]這一區(qū)段內(nèi)，工況2與其他工況相比有較明顯的波浪衰減，波面高差約為0.05 m。

圖15 樹枝密度與波面位置關(guān)系Fig.15 The relationship between branch density and wave surface position

4.2.3 樹干高度對消浪效果的影響

工況8～11對樹干高度進(jìn)行調(diào)整，即改變了水體的上、下分層界限的高度（對應(yīng)圖3中的d1）。圖16是對應(yīng)工況2、8～11的沿程波峰位置處的波面位置變化，圖中坐標(biāo)軸和折線的表示意義與圖14一致。類似的，在沿程第一段波長的區(qū)段內(nèi)，沿程波面高度的區(qū)別不大，樹干高度的變化對波浪衰減的影響體現(xiàn)不明顯，而在此之后的區(qū)段，樹干高度對波面高度變化的影響效果較為明顯。

圖16 樹干高度與波面位置關(guān)系Fig.16 The relationship between vegetation density and wave surface position

樹干高度的變化改變了圖3中水體上下層分界線的位置，依據(jù)波浪傳播過程中大部分波浪能量集中在表層水體內(nèi)這一特點(diǎn)，本文討論樹干高度變化的影響主要分析其對表層水體的影響，其主要體現(xiàn)在以下兩個(gè)方面：（1）改變樹干和樹枝在水體中的作用區(qū)域，即改變了水體不同區(qū)域的拖曳力大?。唬?）樹干-樹枝組合作用、樹枝主導(dǎo)作用、樹干主導(dǎo)作用3種對水體表面紊動消能的影響。

由圖16可知，植被的樹干高度過高或過低，其工況對應(yīng)的沿程波浪的衰減效果均較差：工況8和工況11分別對應(yīng)樹干高度hL=0.2 m和hL=1.0 m的兩種情況，相對其他3種工況（工況2、9、10），其沿程波面高度值總體較大。其中，工況8對應(yīng)的波浪在沿程傳播過程中，在區(qū)域I和區(qū)域II的中前段，即水深較深的區(qū)域，表面波完全受樹枝的拖曳力作用。而由于樹枝的直徑較小、作用力方向與波浪傳播方向有傾角等原因，其波浪衰減效果較差。工況11對應(yīng)的波浪在區(qū)域I同時(shí)受到樹干和樹枝的拖曳力作用，隨著波浪沿斜坡向上爬升，水深逐漸減小，樹干也逐漸完全代替樹枝對波浪施加與其傳播方向一致的拖曳力。

對比5種工況，分析其波浪沿程波面高度的變化規(guī)律可以得到結(jié)論：在水深較深、波能較大的情況下，相比于樹枝主導(dǎo)作用和樹干主導(dǎo)作用兩種情況，樹枝與樹干組合作用于水體表面時(shí)，消浪效果更好。這可能是由于組合作用下，既有樹干較強(qiáng)的拖曳力作用，也有樹枝形態(tài)對水體的擾動作用，增加了水體紊動造成的波能損耗。而對比工況8、9與其他3種工況在中后段的波面高度變化情況，可以發(fā)現(xiàn)：在水深較淺、波能衰減程度較大的情況下，若樹干代替樹枝對波浪施加拖曳力，波能衰減較明顯，樹干起主要消浪作用。綜合以上幾種工況，工況9（樹干高度hL=0.4 m）消浪效果最好。

4.2.4 樹枝傾斜角度對消浪效果的影響

調(diào)整樹枝的傾斜角度，即改變樹枝與波面線的夾角。當(dāng)其夾角為90°時(shí)，即樹枝對波浪施加與其速度方向一致的拖曳力，消浪效果最好。由圖17，對比4種工況，可以發(fā)現(xiàn)樹枝平均角度 θ =π/6（工況9）和θ=π/3（工況14）兩種情況的消浪效果較優(yōu)。在波浪傳播過程中，波形存在峰谷不對稱和波峰的前后不對稱，且隨著波浪不斷向近岸方向傳播，波形的不對稱現(xiàn)象也在逐漸加劇。因此可以認(rèn)為工況9和工況14在波浪的傳播過程中，相對于其他兩種工況，樹枝作用于波浪的拖曳力方向與波浪速度方向總體更趨于一致，能增加入射波浪在豎直方向上勢能和水平方向動能的消耗，消浪效果較優(yōu)。

圖17 樹枝角度與波面位置關(guān)系Fig.17 The relationship between branch tilt angle and wave surface position

4.3 波浪因素對消浪效果的影響規(guī)律分析

采用無因次參數(shù)分析影響植被消浪效果的波浪因素，波浪參數(shù)一般包括：相對波高Hs/d，波陡H/Ls（ δ），相對水深d/Ls（表2），植被覆蓋區(qū)域取消浪效果相對較優(yōu)的工況9的植被特性參數(shù)。其中，植被覆蓋區(qū)域植被特性不變時(shí)，相對水深改變的同時(shí)會引起樹干、樹枝入水深度的改變，并產(chǎn)生影響，其本質(zhì)上與樹干高度對消浪效果的影響機(jī)理相同，故本節(jié)主要討論相對波高Hs/d和波陡H/Ls（δ）兩個(gè)波浪因素對消浪效果的影響規(guī)律。

由于波浪因素已發(fā)生改變，故單純分析沿程波面位置的變化不具有代表性，故取有植被作用下的沿程波高為HT，無植被作用下沿程波高為HI，采用透浪系數(shù)Kt=HT/HI研究植被作用下對波高的影響，透浪系數(shù)越小，消浪效果越好。由3.1節(jié)的分析可知，隨著水深變淺，波長沿程逐步減小，故將沿程分為[0, 7)、[7, 14)、[14, 21)、[21, 27)、[27, 32)5 個(gè)不等間距的區(qū)間段，并對區(qū)段內(nèi)的透浪系數(shù)及其對應(yīng)的沿程波高值的變化規(guī)律進(jìn)行分析。

表2 植被特性模擬工況表Table 2 The arrangement of simulation condition about vegetation characteristics

4.3.1 相對波高HS/d對消浪效果的影響

圖18是對應(yīng)工況a～d的透浪系數(shù)隨相對波高變化，x軸為波浪向岸傳播的水平方向距離，z軸為透浪系數(shù)大小。圖19是相應(yīng)的沿程波高隨相對波高變化圖，z軸為波高大小。圖中折線上最后一個(gè)數(shù)據(jù)點(diǎn)均為波浪傳播至破碎前的波峰位置處的波面高度值。相對波高的變化范圍為 0 .10≤H0/d≤0.25，透浪系數(shù)沿程變化規(guī)律可分為兩段：在區(qū)域I和區(qū)域II的前段，波浪爬升、水流的慣性等因素導(dǎo)致波浪的變形情況較復(fù)雜，透浪系數(shù)無明顯規(guī)律；而在區(qū)域II的中后段，波浪變形趨于穩(wěn)定，在這一范圍內(nèi)，隨著相對波高的增大，透浪系數(shù)減小趨勢明顯。這是由于相對波高越大，枝干對水體表層作用的范圍越大，波能損耗越大。對應(yīng)的，由圖19可見，在區(qū)域I和區(qū)域II的前段，沿程波高值大小與相對波高大小呈正比關(guān)系，而隨著沿程趨于近岸方向，波高值的差距呈縮小趨勢，雖然波高值的變化還受到爬升、破碎的影響，但總體而言，隨著相對波高的增大，沿程波高值衰減速度相對較快。

圖18 相對波高與透浪系數(shù)關(guān)系Fig.18 The relationship between relatively wave height and transmitted coefficient

圖19 相對波高與沿程波高關(guān)系Fig.19 The relationship between relatively wave height and wave height along the path

4.3.2 波陡HS/LS(δ)對消浪效果的影響

圖20是對應(yīng)工況a、e～g的透浪系數(shù)隨波陡變化，x軸為波浪向岸傳播的水平方向距離，z軸為透浪系數(shù)。圖21是相應(yīng)的沿程波高隨波陡變化，z軸為波高大小。

圖20 波陡與透浪系數(shù)關(guān)系Fig.20 The relationship between wave steepness and transmitted coefficient

圖21 波陡與沿程波高關(guān)系Fig.21 The relationship between wave steepness and wave height along the path

波陡的變化范圍為 0 .029≤δ≤0.048，透浪系數(shù)沿程變化規(guī)律可分為兩段：在區(qū)域I和區(qū)域II的前段，與圖18類似的，波浪爬升、水流的慣性等因素導(dǎo)致波浪變形變化較復(fù)雜，透浪系數(shù)無明顯規(guī)律；而在區(qū)域II的中后段，波浪變形趨于穩(wěn)定，在這一范圍內(nèi)隨著波陡的增大透浪系數(shù)呈現(xiàn)增大趨勢。由圖21可以看出，在給定的波陡范圍內(nèi)，植被的消浪效果較好，波浪沿程平緩的向岸方向傳播，波浪沿程不會發(fā)生破碎等劇烈紊動現(xiàn)象，即隨著波陡的增大，不會產(chǎn)生的較大的波能損耗。而此時(shí)，波長對消浪效果起主導(dǎo)作用的可能性較大，大波浪較大的往復(fù)運(yùn)動幅度增加了水體與植被之間的相互作用，即大波浪作用時(shí)的消浪效果要好于小波浪的消浪效果。類似的，楊建民[18]、郭曉宇等[26]根據(jù)研究植被消浪的試驗(yàn)數(shù)據(jù)及現(xiàn)象觀察也得到了相同的結(jié)論。根據(jù)線性波理論，同一波高條件下，長波的水質(zhì)點(diǎn)速度沿水深變化較短波慢，相應(yīng)的長波受到的影響可能會比較顯著[26]。

5 結(jié)論

針對近岸種植木本植物這一新型的海岸生態(tài)軟防護(hù)措施，本文考慮植被枝干的分區(qū)作用，在N-S方程中分別考慮樹枝和樹干的拖曳力影響，提出了木本植被保護(hù)下波浪沿斜坡爬升的表面波衰減理論模型，并采用有限差分法離散方程和MAC方法來跟蹤自由曲面上的水顆粒軌跡。以波浪沿1/30的斜坡爬升為算例，對比討論了有無植被作用下波浪沿傾斜海灘上傳播過程，并將算例結(jié)果與以往試驗(yàn)結(jié)果規(guī)律進(jìn)行對照，驗(yàn)證了數(shù)值模型的有效性。在木本植被均勻覆蓋岸坡這一條件下，分別討論了植被分布密度、單株植物樹枝密度、樹干高度、樹枝傾斜角度4種植被特性和相對波高、波陡兩種波浪因素對植被消浪效果的影響，得到植被消浪的基本規(guī)律。

針對植被特性對消浪效果的影響規(guī)律研究，結(jié)果表明：（1）單位長度內(nèi)的植株的變化范圍為0.1≤NL≤3時(shí)，在保證消浪效果的情況下，從經(jīng)濟(jì)性和岸灘防護(hù)的兩方面考慮，NL=1時(shí)效果最好，實(shí)際工程中，單純增加植被的覆蓋密度來增加其對波浪的衰減作用的做法有待進(jìn)一步討論，其可能導(dǎo)致植被區(qū)域內(nèi)水體波動較為劇烈，局部區(qū)段內(nèi)波面高度變化較大，不利于生態(tài)和岸灘防護(hù)；（2）單植株的樹枝數(shù)變化范圍為3≤NU≤6時(shí)，樹枝密度越大，波面位置變化越趨于平緩，在NU=6時(shí)消浪效果最好；（3）樹干高度變化范圍為0.2 m≤hL≤1.0 m時(shí)，樹干高度為hL=0.4 m消浪效果最好，在水深較深，波能較大的情況下，樹枝-樹干組合情況增加了水體紊動造成的波能損耗，消浪效果較好；在水深較淺、波能衰減程度較大的情況下，樹干起主要消浪作用；（4）樹枝傾斜角度在 π/6≤ θ≤π/3范圍內(nèi)時(shí)，θ=π/6和 θ=π/3兩種情況樹枝作用于波浪的拖曳力方向與波浪速度方向總體更趨于一致，消浪效果較優(yōu)。針對波浪因素對消浪效果的影響規(guī)律研究，結(jié)果表明：相對波高的變化范圍為 0.10≤HS/d≤0.25，透射系數(shù)總體隨著波陡與相對波高的增加呈緩慢減小趨勢；波陡變化范圍為 0 .012≤δ≤0.048，在這一范圍內(nèi)，波浪沿程傳播平緩，波陡的變化對消浪效果影響不明顯，波長對消浪效果起主導(dǎo)作用，大波浪作用時(shí)的消浪效果好于小波浪。文中的計(jì)算結(jié)果和消浪的基本規(guī)律可為實(shí)際的護(hù)岸工程和生態(tài)景觀設(shè)計(jì)提供參考。

植被護(hù)岸這一軟防護(hù)措施在實(shí)際使用中，多設(shè)置在傾斜底坡的海灘或近岸結(jié)構(gòu)上，而由于波浪與植被本身作用十分復(fù)雜，因此，提出考慮植被形態(tài)和斜坡爬升的表面波衰減理論模型、分析波浪作用下植被參數(shù)與波要素對消浪效果的影響規(guī)律，對工程實(shí)際具有重要的意義。但由于波浪在植被覆蓋區(qū)域圍內(nèi)的紊動情況復(fù)雜，且未考慮樹枝部分的柔性結(jié)構(gòu)對消浪效果的影響，對植被特性的討論還有待進(jìn)一步研究與分析。