葛昭佩，唐軍*，趙楚嫣

( 1. 大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023)

1 引言

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

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

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

2 數(shù)值模型

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

2.1 控制方程

waves2Foam 求解器控制方程基于OpenFOAM 標(biāo)準(zhǔn)，采用氣液兩相流方程。為考慮波浪與植被相互作用，在動(dòng)量方程和湍流模型中引入植被源項(xiàng)，即拖曳阻力項(xiàng)和慣性力項(xiàng)[17]（稱之為宏觀數(shù)值模型，以下未特別說明的均指該數(shù)值模型）。控制方程為

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

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

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

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

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

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

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

2.3 邊界條件

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

3 模型驗(yàn)證

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

3.1 植被水域波浪傳播驗(yàn)證

選取表2 中的兩種試驗(yàn)工況驗(yàn)證模型模擬植被水域波浪傳播的有效性，其中CD值由精細(xì)數(shù)值模型計(jì)算，CD值計(jì)算過程以表2 中工況1 為例，數(shù)值模擬中設(shè)置數(shù)值水槽長(zhǎng)為18 m，高為0.6 m，寬為0.12 m，取植被水域長(zhǎng)為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 為精細(xì)模擬中工況1 的波面演化和x=0.6 m（x=0 為植被域起點(diǎn)）處植被上的波浪力及截面平均速度，使用式（10）計(jì)算該處周期平均CD值，并利用同樣的方法計(jì)算每一排植株處的周期平均CD值后，做平均得出空間–周期平均CD值作為該工況的代表CD值。

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

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

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

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

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

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

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

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

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

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

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

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

4 模擬結(jié)果分析

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

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

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

圖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)

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

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

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

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

5 結(jié)論

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

圖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