植被對波浪作用下床面切應(yīng)力影響的數(shù)值模擬分析

葛昭佩，唐軍*，趙楚嫣

( 1. 大連理工大學(xué) 海岸和近海工程國家重點實驗室，遼寧 大連 116023)

1 引言

波浪作用下的床面切應(yīng)力對海床泥沙起動具有關(guān)鍵性作用，是模擬分析近岸泥沙輸移、岸灘演變的重要參數(shù)。近岸植被具有防浪固灘護岸的作用，研究近岸植被對波浪作用下床面切應(yīng)力的影響有助于理解其防護海岸的作用機制。

目前關(guān)于波浪作用下床面切應(yīng)力的研究取得了一定成果。Jonsson[1]從理論公式出發(fā)，推導(dǎo)了層流邊界層的床面切應(yīng)力及邊界層厚度公式，并引入波浪底摩阻系數(shù)及邊界層流態(tài)判別方法，給出了紊流邊界層的波浪摩阻系數(shù)經(jīng)驗公式；孔令雙等[2]基于邊界層理論，導(dǎo)出了波流環(huán)境中床面切應(yīng)力，并建立了泥沙起動及水體挾沙力公式；蔡翠蘇[3]基于水槽試驗，探討了規(guī)則波和不規(guī)則波作用下底摩阻系數(shù)計算方法；齊富康等[4]基于渤海海峽觀測數(shù)據(jù)，分別應(yīng)用Soulsby和Van Rijn 的模型計算波流環(huán)境中床面切應(yīng)力；Lin和Zhang[5]建立三維數(shù)值模型，模擬計算線性波、Stokes 波、孤立波等條件下層流邊界層床面切應(yīng)力，并分別與解析解比較分析；滕涌等[6]基于ECOMSED波流耦合底邊界層模型，分析了波高、近底流速及水深對波流條件下床面切應(yīng)力的影響；Larsen 和Fuhrman[7–8]基于waves2Foam 模擬海嘯波的傳播及爬升過程，著重討論了海嘯波作用下的邊界層及床面切應(yīng)力并提出了一種預(yù)測海嘯波作用下瞬時邊界層厚度和床面切應(yīng)力的方法。相對裸床，植被影響下波浪水動力特性復(fù)雜，目前關(guān)于植被水域波浪作用下床面切應(yīng)力的計算以經(jīng)驗公式為主。Wang 等[9]基于水槽試驗分析了水流和波流條件下植被水域平均流速及湍流動能，并由湍流動能推算植被水域床面切應(yīng)力，發(fā)現(xiàn)純水流時植被水域床面切應(yīng)力沿程降低，而波流條件下植被水域床面切應(yīng)力增大并在離開植被區(qū)后顯著大幅降低；Reidenbach 和Thomas[10]觀測了弗吉尼亞海岸保護區(qū)水動力，并應(yīng)用經(jīng)驗公式分別計算波、流的床面切應(yīng)力，發(fā)現(xiàn)與純水流時相比波流條件下床面切應(yīng)力大幅增大，位于大葉藻水域測點處的床面切應(yīng)力始終小于裸床處且小于泥沙起動臨界切應(yīng)力，有效保護了海床免受侵蝕；陳家貴和沈小雄[11]基于Jonsson[1]提出的層流邊界層公式，結(jié)合水槽試驗數(shù)據(jù)，分析了入射波高及柔性植被對最大床面切應(yīng)力的影響；李勰等[12]基于Jonsson[1]提出的邊界層最大切應(yīng)力公式及Luhar 等[13]提出的植被帶中邊界層流速公式，分析了剛性植被根莖對床面切應(yīng)力沿程變化的影響。

總體來看，目前關(guān)于植被對波浪或波流作用下床面切應(yīng)力影響的研究主要基于水槽試驗或現(xiàn)場觀測數(shù)據(jù)并結(jié)合經(jīng)驗公式分析計算，但由于植被影響下的流速剖面及雷諾應(yīng)力分布復(fù)雜，經(jīng)驗公式未能很好地描述波浪條件下植被水域床面切應(yīng)力[9,14–15]。本文以O(shè)penFOAM 中的waves2Foam 求解器[16]為基礎(chǔ)，在其中引入植被源項，建立植被水域波流數(shù)值水槽，通過計算全水深速度場及邊界層內(nèi)速度梯度而直接給出床面切應(yīng)力分布。在應(yīng)用已有實驗數(shù)據(jù)及理論公式對數(shù)值模型驗證的基礎(chǔ)上，著重分析了入射波高、植被密度、植被淹沒高度及水流對植被水域波浪作用下床面切應(yīng)力的影響特性。

2 數(shù)值模型

waves2Foam 求解器是由Jacobsen 等[16]基于Open-FOAM 開發(fā)的一個用于模擬波浪的第三方模塊庫，提供了水流、線性波、Stokes 波、孤立波等造波文件，亦可自定義造波文件實現(xiàn)其他類型波浪模擬。

2.1 控制方程

waves2Foam 求解器控制方程基于OpenFOAM 標準，采用氣液兩相流方程。為考慮波浪與植被相互作用，在動量方程和湍流模型中引入植被源項，即拖曳阻力項和慣性力項[17]（稱之為宏觀數(shù)值模型，以下未特別說明的均指該數(shù)值模型）?？刂品匠虨?/p>

傳統(tǒng)湍流模型應(yīng)用于氣液兩相流時，由于氣液界面處較大的密度梯度[19]及湍流動能的過度評估[20]會造成大波陡波傳播過程中的波高非物理衰減。為抑制此種現(xiàn)象，本文采用Larsen 和Fuhrman[20]提出的修正k-ε湍流模型。

表1 經(jīng)驗系數(shù)取值Table 1 Default values for the closure coefficient

waves2Foam 為兩相流求解器，其定義體積分數(shù)α （ 0 ≤α ≤1）來描述每個網(wǎng)格中各相的體積比重。各相的運動滿足的控制方程為

式（8）第三項是用于保持氣液界面清晰的人工壓縮項；ur,i為 壓縮速度；壓縮系數(shù)Cα∈[0，1]，其默認值為1。

2.2 植被拖曳力系數(shù)

植被拖曳力系數(shù)CD是對植被水域水動力特性進行宏觀數(shù)值模擬時的關(guān)鍵參數(shù)，與波浪、水流、植被密切相關(guān)，直接影響數(shù)值模擬的準確性，但其取值尚無普適方法。本文將上述控制方程中的植被源項去除（令CD、CM、Ckp、Cεp4 項經(jīng)驗系數(shù)均取為0），通過直接考慮剛性植被外形對波流的擾動，詳細模擬純波及波流條件下植被水域的水動力特征（稱之為精細數(shù)值模型，在本文中僅用于計算各工況代表CD值），通過提取植被上的波流力及其所在截面的沿植被高度平均流速并使用Morison 公式（式（10））直接計算植被拖曳力系數(shù)。

式中，F(xiàn)D為 植株上的拖曳力；FI為植株上的慣性力；F為植株上的波流力，Dalrymple 等[22]指出阻力與流速之間無相位差，阻力的做功可由阻力和流速的乘積在周期內(nèi)積分獲得，然而慣性力與流速之間存在 π/2的相位差，使得慣性力在周期內(nèi)的做功為0，故上式計算過程中可忽略慣性力[18]，F(xiàn)根據(jù)植被上壓力及附近流場由式（11）計算；u為植株所在截面的平均流速；ρ為水的密度；hv為 植被淹沒高度；bv為植被直徑。

2.3 邊界條件

數(shù)值水槽入口與出口采用松弛區(qū)造波與消波，入口處速度場、壓力場和自由波面由波浪理論公式給出，出口采用自由出流條件；床面采用無滑移邊界條件；兩側(cè)壁采用周期性邊界條件，該邊界條件能夠消除兩側(cè)壁面的影響，可將水槽簡化為較窄的數(shù)值水槽以提高計算效率，且已被成功應(yīng)用于波浪的數(shù)值模擬[23–24]；頂端采用壓力出口邊界。時間步長由庫朗數(shù)自動調(diào)節(jié)，為保證計算穩(wěn)定，最大庫朗數(shù)設(shè)置為0.25。

3 模型驗證

由于缺乏波浪作用下植被水域床面切應(yīng)力實測數(shù)據(jù)，本文將分別驗證植被水域波面演化和無植被時波浪作用下床面切應(yīng)力。采用海岸和近海工程國家重點實驗室開展的植被影響下波浪傳播物理水槽試驗結(jié)果[25]驗證本文模型模擬植被水域波浪傳播的有效性；采用理論公式驗證層流邊界層時的床面切應(yīng)力，采用徐華等[26]的試驗驗證紊流邊界層時的床面切應(yīng)力。

3.1 植被水域波浪傳播驗證

選取表2 中的兩種試驗工況驗證模型模擬植被水域波浪傳播的有效性，其中CD值由精細數(shù)值模型計算，CD值計算過程以表2 中工況1 為例，數(shù)值模擬中設(shè)置數(shù)值水槽長為18 m，高為0.6 m，寬為0.12 m，取植被水域長為1.08 m，其中以圓柱代表剛性植被，整體網(wǎng)格尺寸 Δx=Δy=0 .03 m， Δz=0.01 m，為更好地貼合植被外形，分別加密近植被域及單株植被附近網(wǎng)格。邊界條件同2.3 節(jié)，此外，植被表面設(shè)為無滑移邊界條件。

表2 植被及波浪參數(shù)Table 2 Parameters of vegetation and waves

圖1 為精細模擬中工況1 的波面演化和x=0.6 m（x=0 為植被域起點）處植被上的波浪力及截面平均速度，使用式（10）計算該處周期平均CD值，并利用同樣的方法計算每一排植株處的周期平均CD值后，做平均得出空間–周期平均CD值作為該工況的代表CD值。

圖1 考慮植被外形擾動的模擬結(jié)果Fig. 1 Simulation results considering the disturbance of vegetation shape

將利用精細模型計算的各工況代表CD值代入宏觀模型中模擬計算。波浪經(jīng)過植被水域的波面演化及流速衰減的模擬值與試驗值對比如圖2 和圖3 所示，由圖可知，模型計算結(jié)果和試驗數(shù)據(jù)[25]符合性較好，表明利用精細數(shù)值模型計算CD值準確并且使用宏觀數(shù)值模型能夠準確模擬植被水域波浪及波流水動力變化。

圖2 植被水域波面演化Fig. 2 Free surface evolution along the vegetation zones

圖3 植被水域流速衰減驗證Fig. 3 Verification of velocity attenuation in vegetation zones

3.2 波浪作用下床面切應(yīng)力驗證

本文數(shù)值模型通過求解雷諾平均Navier-Stokes 方程得出波浪作用下全水深速度場及邊界層內(nèi)速度梯度，并由公式（12）計算床面切應(yīng)力。

為驗證模型計算床面切應(yīng)力的準確性，此處設(shè)置3 組驗證工況，其中工況1、工況2 為本試驗工況，工況3 為徐華等[26]的試驗工況，波浪參數(shù)如表3 所示。工況1 和工況2 的邊界層雷諾數(shù)均小于1.26×104，依據(jù)Jonsson[1]提出的判定標準，兩者均屬于層流邊界層，工況3 則屬于光滑紊流邊界層。

表3 驗證工況參數(shù)Table 3 Parameters of verification conditions

邊界層流速剖面及床面切應(yīng)力驗證結(jié)果如圖4至圖6 所示。其中，工況1 采用線性波理論解驗證，工況2 采用五階Stokes 波理論解驗證，工況3 采用試驗測量值[26]驗證。從圖中可以看出，待波浪場穩(wěn)定后，工況1 模擬值與理論值符合性較好；工況2 床面切應(yīng)力在波谷時模擬值較理論值偏大，但差值小于10%；工況3 模擬值與試驗測量值在一個完整周期內(nèi)總體符合良好，表明該數(shù)值模型能夠準確計算波浪作用下床面切應(yīng)力。

圖4 理論值與模擬值對比（工況1）Fig. 4 Comparison between theoretical and simulated values (case 1)

圖5 理論值與模擬值對比（工況2）Fig. 5 Comparison between theoretical and simulated values (case 2)

圖6 床面切應(yīng)力對比（工況3）Fig. 6 Comparison of bed shear stress (case 3)

4 模擬結(jié)果分析

為分析不同入射波高、植被密度、植被淹沒高度及水流流速條件下，植被水域床面切應(yīng)力的分布特征，本文共模擬計算22 種工況，見表4。其中工況1、2、3、12、14、15、18、19、20 為文獻[25]的試驗工況，這9 種工況的波面演化及流速衰減模擬值與試驗值均符合良好。

表4 模擬工況參數(shù)Table 4 Parameters of numerical simulation conditions

純波時，波浪經(jīng)過植被水域后波高衰減的同時床面切應(yīng)力也隨之減小，見圖7（x=0 m 處為植被起點），這是由于植被的阻水作用，近底流速沿植被水域衰減使得邊界層流速梯度減小所致（圖8）。圖7 中植被前后的波面與床面切應(yīng)力均存在 π/4的相位差，說明植被對波高和床面切應(yīng)力的衰減作用是同步的。此外，該相位差與式（13）、式（14）相一致，進一步驗證了該數(shù)值模型的準確性。

圖7 波面及床面切應(yīng)力（工況1）Fig. 7 Free surface and bed shear stress (case 1)

圖8 近底流速剖面（工況1）Fig. 8 The near-bottom velocity profile (case 1)

當植被處于淹沒狀態(tài)時，植被的阻水效應(yīng)僅發(fā)生在冠層以下，使這一區(qū)域流速減小且阻水作用隨著淹沒高度的增大而增大，而冠層以上由于過水斷面的減小使得流速增大且超過無植被時的流速，故在植被頂端形成較大的速度梯度（圖9）；此外，由于邊界層的作用在近底處存在一速度峰值。上述水平速度沿水深的分布規(guī)律與Liu 等[29]的試驗結(jié)果一致。

圖9 淹沒植被沿水深速度剖面Fig. 9 Longitudinal velocity profile of submerged vegetation

同向水流的存在會增加波浪水質(zhì)點速度的非線性，同時增大正向床面切應(yīng)力幅值，減小負向床面切應(yīng)力幅值，使得一個波周期內(nèi)正向切應(yīng)力時間延長，負向切應(yīng)力時間縮短，如圖10 所示。

圖10 不同流速下床面切應(yīng)力變化Fig. 10 Variation of bed shear stress under different current velocity

5 結(jié)論

本文OpenFOAM 中的waves2Foam 求解器建立了含植被水域的波流數(shù)值水槽，模擬分析了波浪、水流、植被要素對植被水域床面切應(yīng)力衰減的影響。模擬結(jié)果表明：當波浪傳播到植被前端時，床面切應(yīng)力會出現(xiàn)增大現(xiàn)象，并隨著入射波高增大、水流流速增大、植被密度增大此種現(xiàn)象愈加明顯；純波時，波高和植被密度的增大均會導(dǎo)致植被水域床面切應(yīng)力衰減幅度的增大，并且床面切應(yīng)力在植被水域前段的衰減幅度較大，后段衰減幅度減??；對于完全淹沒植向床面切應(yīng)力幅值減小，使得一個波周期內(nèi)正向切應(yīng)力時間延長負向切應(yīng)力時間縮短；當水流流速較小時，床面切應(yīng)力的衰減率與無水流時相近，之后隨著流速增大衰減率隨之增大。

圖11 不同工況下植被水域最大床面切應(yīng)力分布Fig. 11 Distribution of maximum bed shear stress in vegetation zones under different conditions

圖12 不同流速下最大床面切應(yīng)力衰減率Fig. 12 Decay rate of maximum bed shear stress at different current velocities