陳 明，婁 廈，2，劉曙光，2，RADNAEVA Larisa Dorzhievna，NIKITINA Elena，CHALOV Sergey Romanovich

（1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.同濟(jì)大學(xué)長(zhǎng)江水環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；3.俄羅斯科學(xué)院西伯利亞分院貝加爾湖自然管理研究所，烏蘭 烏德 670047，俄羅斯；4.莫斯科羅蒙諾索夫國(guó)立大學(xué)地理學(xué)院，莫斯科 119991，俄羅斯）

海岸帶、濱海濕地等區(qū)域普遍存在的水生植物，對(duì)水動(dòng)力和水環(huán)境具有重要影響。一方面，植物可以減小水流流速，削弱波浪對(duì)岸灘的沖擊；另一方面，植物也會(huì)改變水體的紊動(dòng)特征，影響泥沙輸運(yùn)和污染物遷移過(guò)程。因此，研究含植物水流水動(dòng)力特征及泥沙輸運(yùn)過(guò)程對(duì)岸灘保護(hù)和水環(huán)境修復(fù)具有重要意義。

目前，大多數(shù)含植物水流數(shù)學(xué)模型將植物群設(shè)置為多體結(jié)構(gòu)或多孔介質(zhì)進(jìn)行模擬。多體結(jié)構(gòu)模型是將每根植物視為獨(dú)立的個(gè)體，進(jìn)而可將整體植物群等效為柱體群結(jié)構(gòu)，利用網(wǎng)格讀取每根植物的邊界。如Chang等［1］利用分離渦模擬（DES），構(gòu)建了圓柱群陣列，對(duì)每根植物的幾何形狀進(jìn)行了精細(xì)模擬，分析了淹沒(méi)剛性植物對(duì)水流流速和紊動(dòng)的影響，在植物尾跡區(qū)觀察到U形渦旋。多孔介質(zhì)模型是用一個(gè)帶有孔隙率的整體代替植物群，通過(guò)在計(jì)算方程中添加概化固體阻力源項(xiàng)來(lái)模擬。如Zhan等［2］采用多孔介質(zhì)模擬方法，通過(guò)設(shè)置統(tǒng)一的孔隙率將植物群概化為等間距圓柱群，并引入一個(gè)非恒定的慣性阻力系數(shù)來(lái)模擬不均勻的植物阻力分布。與多體結(jié)構(gòu)模型相比，多孔介質(zhì)模型雖然可以節(jié)省大量計(jì)算時(shí)間，但可能存在植物群內(nèi)部水體結(jié)構(gòu)模擬失真的問(wèn)題［3］。如Yu等［4］利用Fluent模型對(duì)比了兩種方法下剛性挺水植物與水流相互作用，發(fā)現(xiàn)多體結(jié)構(gòu)模型可以相對(duì)較好地模擬出植物影響下的平均流速場(chǎng)和紊流場(chǎng)，多孔介質(zhì)模型只能模擬植物群整體產(chǎn)生的大尺寸紊動(dòng)而不能很好地模擬單根植物引起的的小尺寸紊動(dòng)。因此，為了更加精確地模擬植物對(duì)水動(dòng)力特征的影響，獲得植物群內(nèi)部的水流結(jié)構(gòu)和泥沙運(yùn)動(dòng)規(guī)律，本研究選用多體結(jié)構(gòu)方法模擬剛性沉水植物。

近年來(lái)越來(lái)越多的研究發(fā)現(xiàn)，植物的存在使得水體中懸浮泥沙濃度顯著增加，含植物床面中植物引起的紊動(dòng)對(duì)泥沙起動(dòng)及懸浮起主導(dǎo)作用［5-6］。因此，僅考慮平均流速的床面切應(yīng)力并不適用于含植物水流中泥沙起動(dòng)與輸運(yùn)的模擬。Tinoco和Coco［5，7］用紊動(dòng)強(qiáng)度代替平均流速，采用物理模型試驗(yàn)數(shù)據(jù)，對(duì)傳統(tǒng)臨界希爾茲（Shields）數(shù)進(jìn)行修正，提出了適用于波浪、水流條件下的含植物床面泥沙起動(dòng)判斷方法，但該方法是基于特定的物理模型試驗(yàn)結(jié)果進(jìn)行的推導(dǎo)，其應(yīng)用范圍仍需進(jìn)一步探討。Lou等［8］進(jìn)一步修正Shields數(shù)并應(yīng)用于NHWAVE模型，模擬了植物影響下的泥沙懸浮過(guò)程，并獲得了較為精確的結(jié)果。然而，上述改進(jìn)的模型需要對(duì)較多參數(shù)進(jìn)行校正，其中一些參數(shù)的物理意義難以得到合理解釋。因此，本文在前期物理模型試驗(yàn)的基礎(chǔ)上，基于當(dāng)前被廣泛認(rèn)可的利用流速和植物特征計(jì)算植物尾流紊動(dòng)的方法［9-10］，從理論上采用紊動(dòng)能對(duì)傳統(tǒng)Shields數(shù)進(jìn)行修正，從而獲得可應(yīng)用于含植物床面的泥沙起動(dòng)與懸浮的預(yù)測(cè)方法。該方法具有較強(qiáng)的物理意義，且并未引入新的參數(shù)。應(yīng)用該方法，改進(jìn)了Flow-3D模型中的泥沙模塊，開(kāi)展了波浪通過(guò)剛性沉水植物區(qū)時(shí)泥沙懸浮的模擬。

1 數(shù)學(xué)模型建立

1.1 基本方程

1.1.1 水流方程

Flow-3D模型中笛卡爾坐標(biāo)系下的控制方程如下：

連續(xù)性方程：

動(dòng)量方程：

式（1）—（2）中：VF為流體的體積分?jǐn)?shù)；ρ為流體密度；t為時(shí)間；xi為笛卡爾坐標(biāo)（x，y，z）；ui為速度分量；Ai為i向流體的面積分?jǐn)?shù)；p為壓力；Gi為i向體加速度；fi為i向粘滯力加速度。

采用Yakhot和Orszag［11］提出的RNGk-ε模型作為紊流模型。紊動(dòng)能kT和紊動(dòng)耗散率εT滿足如下方程：

式（3）—（4）中：PT為速度梯度引起的紊動(dòng)能產(chǎn)生項(xiàng)；GT為浮力引起的紊動(dòng)能產(chǎn)生項(xiàng)；Dk和Dε為擴(kuò)散項(xiàng)；C1ε、C2ε和C3ε為模型參數(shù)。

1.1.2 泥沙運(yùn)動(dòng)方程

Flow-3D模型中，以臨界切應(yīng)力作為判斷泥沙起動(dòng)的條件，采用Soulsby-Whitehouse公式［12］計(jì)算量綱一化臨界Shields數(shù)θcr：

式中：d*為量綱一化泥沙粒徑；ds為泥沙粒徑；ρs為泥沙密度；||g||為重力加速度的大??；μ為流體動(dòng)力粘度。

泥沙挾帶上升速度ulift用于計(jì)算床沙轉(zhuǎn)化為懸浮狀態(tài)的量，采用Mastbergen和Berg［13］的理論公式：

式中：α為泥沙挾帶系數(shù)；ns為床面法向向量；θ為Shields數(shù)；τ為床面切應(yīng)力。

懸移質(zhì)泥沙擴(kuò)散方程為

式中：Cs為懸浮泥沙的質(zhì)量濃度為水沙混合物的速度，利用ulift和水流速度相乘的通量計(jì)算得到；v為流體運(yùn)動(dòng)粘度；D為擴(kuò)散系數(shù)。

1.2 改進(jìn)的泥沙模塊

由式（7）和（8）可知，F(xiàn)low-3D原始泥沙模塊是基于流速計(jì)算的床面切應(yīng)力來(lái)判斷和計(jì)算泥沙起動(dòng)及懸浮的。本研究考慮植物引起的紊動(dòng)作用對(duì)泥沙輸運(yùn)的影響，對(duì)該模型進(jìn)行了改進(jìn)。

對(duì)于裸床面（無(wú)植物），當(dāng)單向流通過(guò)時(shí)，由床面引起的紊動(dòng)kTb與床面剪切τc之間存在線性關(guān)系［14］：

同時(shí)，τc也可用平均流速U和床面摩擦系數(shù)Cf直接計(jì)算［12］：

由式（10）和（11）可得：

Tang等［15］基于水槽試驗(yàn)，初步證明了式（12）同樣適用于波浪條件，得出：

式中：Cw為比例因子，取值為Cw=0.02±0.01；Uw為波浪水平振蕩速度幅值。

將式（13）代入Flow-3D原始泥沙模塊Shields數(shù)的計(jì)算公式（8），將其改寫(xiě)為基于紊動(dòng)能的Shields數(shù)θ0：

式中：fw為波浪條件下底部摩擦系數(shù)。

Tanino和Nepf［9］根據(jù)水槽試驗(yàn)結(jié)果提出了單向流條件下植物尾流紊動(dòng)kTv與水流平均流速U的關(guān)系公式：

式中：δ為比例因子；CD為植物拖曳力系數(shù)；?為冠層體積分?jǐn)?shù)，對(duì)圓柱群陣列，?=nπd2，其中n為植物密度，d為植物直徑。當(dāng)n較小使得植物間距S＞d/0.56時(shí)，l=d，δ=1.1；當(dāng)n較大使得S≤d/0.56時(shí)，l=S，δ=0.88。Zhang等［10］利用物理模型試驗(yàn)證明了該公式可應(yīng)用于波浪條件，但需滿足Aw/S＞1，其中Aw為波浪軌道偏移量，此時(shí)δ將減小為0.76。結(jié)合式（13）和（15），波浪條件下，含植物尾流紊動(dòng)的床面總紊動(dòng)能為

將式（16）代入式（14），可得到含植物床面的Shields數(shù)θ?為

將式（17）代入式（7），可得到波浪條件下含植物床面的泥沙挾帶上升速度ulift，?的計(jì)算公式（18）。?=0時(shí)，θ?=θ0。用ulift，?替換Flow-3D原始泥沙模塊中泥沙挾帶上升速度ulift，即可得到改進(jìn)后的泥沙模塊，以用于計(jì)算含植物床面上泥沙起動(dòng)及懸浮過(guò)程。

1.3 模型設(shè)定

本文建立基于Flow-3D 的含剛性沉水植物水槽數(shù)學(xué)模型，利用作者前期物理模型試驗(yàn)數(shù)據(jù)［16-17］對(duì)數(shù)值模擬結(jié)果進(jìn)行驗(yàn)證。計(jì)算區(qū)域布置如圖1 所示，長(zhǎng)20 m（x向），寬0.2 m（y向），高1 m（z向）。水槽入口采用波浪邊界（Wave），出口采用出流邊界（Outflow），底面采用無(wú)滑移壁面邊界（No-slip Wall），其他三面均采用對(duì)稱邊界（Symmetry）。

模型中采用直徑為8 mm 的圓柱體構(gòu)建長(zhǎng)6 m，寬0.2 m 的植物區(qū)（圖1），共設(shè)置3組植物布置方案（圖2），其中工況a為裸床對(duì)比工況（無(wú)植物）；工況b中植物密度n為200 株·m-2，高度hv為0.4 m；工況c中n為400株·m-2，hv為0.4 m；為模擬真實(shí)生態(tài)環(huán)境中植物高度的空間變化特征，在工況d中布置4種不同高度的植物（0.2、0.4、0.6、0.8 m），其密度和工況c相同。在模擬植物區(qū)底部鋪設(shè)0.1 m厚的沙床，泥沙粒徑為0.16 mm，密度為1.45 g·cm-3。

圖1 計(jì)算區(qū)域布置示意圖（單位：m；尺寸不按比例）Fig.1 Schematic of computational domain（unit:m,not to scale）

圖2 模擬植物分布示意圖（單位：cm）Fig.2 Schematic of mimic vegetation canopy（unit:cm）

模型中設(shè)置4 組波浪工況W1～W4，波高分別為0.06、0.08、0.10、0.12 m，波周期均為1.8 s。采用五階斯托克斯波進(jìn)行模擬，水槽末端設(shè)置海綿層用于消波。為節(jié)省計(jì)算時(shí)間、確保計(jì)算精度，本模型采用雙重網(wǎng)格嵌套的方式，設(shè)置粗網(wǎng)格覆蓋整個(gè)水槽區(qū)域，耦合細(xì)網(wǎng)格對(duì)植物區(qū)進(jìn)行加密。粗網(wǎng)格的水平間距均為0.04 m，垂向平均間距為0.025 m，在水面和底部對(duì)網(wǎng)格進(jìn)行局部加密。細(xì)網(wǎng)格垂向間距和粗網(wǎng)格相同，水平間距選擇三種不同尺寸進(jìn)行網(wǎng)格敏感性分析。結(jié)果表明，與間距為0.001 m 網(wǎng)格的模擬結(jié)果相比，0.002 5 m 間距網(wǎng)格相差小于10%，而0.004 m間距網(wǎng)格相差大于20%。兼顧計(jì)算精度和效率，選擇細(xì)網(wǎng)格間距為0.002 5 m。

采用上述尺寸的雙重網(wǎng)格嵌套模式，對(duì)模型的造波能力進(jìn)行驗(yàn)證。無(wú)植物時(shí)W4 條件下x=3.5 m 和x=19.98 m 處（入口邊界定為x=0）的計(jì)算值與理論值對(duì)比如圖3所示。在x=3.5 m處，可以生成穩(wěn)定的波浪，平均波高約0.12 m，計(jì)算值與五階斯托克斯波理論值及實(shí)測(cè)值均吻合較好。在x=19.98 m 處，波面升高接近于0，說(shuō)明波浪在海綿層內(nèi)充分衰減，不會(huì)在出口邊界產(chǎn)生反射波影響模擬結(jié)果。因此，本模型能夠準(zhǔn)確模擬波浪傳播，可用于研究植物-波浪-泥沙之間的相互作用。

圖3 不同位置處波面升高計(jì)算值與五階斯托克斯波理論值及實(shí)測(cè)值對(duì)比Fig.3 Comparison of numerical, fifth-order Stokes theoretical and measured results of wave trains at different locations

1.4 模型精度評(píng)價(jià)

本文采用均方根誤差（RMSE）和一致性指數(shù)（IA）定量評(píng)估數(shù)值模型的模擬效果。

式中：Pi為模擬值；Oi為實(shí)測(cè)值；n為樣本個(gè)數(shù)為實(shí)測(cè)值的平均值。RMSE值越小，模擬值與實(shí)測(cè)值越接近，模擬效果越好。IA值在0-1 之間，越接近1，模擬效果越好。

2 計(jì)算結(jié)果分析

2.1 波浪動(dòng)力特性分析

2.1.1 流速分析

圖4為斷面x=7.5 m處波浪水平振蕩速度幅值Uw（=0.5（Umax-Umin），Umax和Umin分別為波浪最大流速和最小流速）在各工況下的對(duì)比結(jié)果，其中hv為工況b和c中植物高度，其值為0.4 m，水平虛線表示植物冠層頂部，后同。無(wú)植物（工況a）的4組波浪條件下，模擬值與實(shí)測(cè)值較為接近，實(shí)測(cè)與模擬的Uw均自水面向底部逐漸減小，說(shuō)明該模型可以較好地模擬波浪流速（RMSE＜0.005 m·s-1，IA＞0.91）。三組植物工況下（工況b～工況d），模擬值與實(shí)測(cè)值分布趨勢(shì)一致，但在個(gè)別條件下出現(xiàn)了一定程度的誤差（RMSE＜0.014 m·s-1，IA＞0.75），如工況c-W1和c-W4中，冠層下方模擬值與實(shí)測(cè)值存在一定偏差，這可能是由于物理模型試驗(yàn)中測(cè)點(diǎn)位置偏差引起的。

圖4 x=7.5 m處波浪水平振蕩速度幅值Uw垂向分布Fig.4 Vertical distributions of horizontal wave orbital velocity amplitude Uw at x=7.5 m

與無(wú)植物工況相比，不同波浪條件下Uw在工況b和c的植物冠層內(nèi)的垂向平均模擬值分別衰減了3%和2%（W1），6%和8%（W2），10%和15%（W3），9%和17%（W4）。植物的存在對(duì)波浪具有明顯的阻礙作用，且植物引起的波浪速度衰減幅度與植物密度和波高均呈正相關(guān)關(guān)系。當(dāng)波高較小時(shí)，水動(dòng)力相對(duì)較弱，植物密度對(duì)速度衰減的影響較弱。反之，隨著波高的增大，水動(dòng)力增強(qiáng)，增大植物密度會(huì)引起更顯著的速度衰減幅度。

垂向均勻密度植物影響下，如圖4b、4c 所示，植物的阻礙作用使得下層水體流速顯著小于上層。在冠層頂部（z/hv=1），存在顯著的流速梯度，出現(xiàn)Uw峰值，且該值隨著植物密度和波高的增大而增大。為進(jìn)一步分析該現(xiàn)象，圖5繪制了工況a-W3和c-W3條件下x=7.5 m 處的波浪最大流速Umax和最小流速Umin的垂向分布情況。在工況c中，垂向上Umax在冠層頂部附近取得最大值，該值甚至大于無(wú)植物工況，說(shuō)明植物的存在對(duì)該處的正向流速有一定的增幅。在冠頂處一個(gè)完整的波浪運(yùn)動(dòng)周期內(nèi)，水分子的運(yùn)動(dòng)軌跡是在冠頂上方與波浪傳播方向相同，為正向流速；而冠頂下方與波浪傳播方向相反，為負(fù)向流速。冠頂處的植物拖曳力和水體慣性力會(huì)引起較大的雷諾剪切應(yīng)力［18］，增大了冠頂上方的正向流速，表現(xiàn)出冠層頂部的Umax極值；而冠頂下方的負(fù)向流速被植物拖曳力減小，使得冠頂附近的水分子運(yùn)動(dòng)出現(xiàn)周期不對(duì)稱的現(xiàn)象，表現(xiàn)出正向的平均流速。該現(xiàn)象可等效為斯托克斯漂移，與粗糙底床上波浪邊界層的流速類似，這也使得很多數(shù)學(xué)模型中將植物冠層概化為用拖曳力和慣性力模擬的粗糙多孔介質(zhì)。

圖5 工況a-W3和c-W3條件下x = 7.5 m處波浪最大流速Umax和最小流速Umin垂向分布Fig.5 Vertical distributions of maximum and minimum wave velocity (Umax and Umin) at x = 7.5 m in test a-W3 and c-W3

在工況d 垂向不均勻密度植物的影響下，植物密度由水體上層向底部逐層增加，導(dǎo)致水體流速隨之逐層減小，且在每個(gè)冠層頂部附近均有較大的流速梯度（圖4d）。同時(shí)，在工況d中，各冠頂上方均有更高植物的存在，此處正向雷諾剪切應(yīng)力對(duì)流速的增強(qiáng)作用會(huì)被上方植物的拖曳力作用抵消，因此在圖4d中，各冠頂附近沒(méi)有出現(xiàn)流速極值。

2.1.2 紊動(dòng)能分析

圖6 為斷面x=7.5 m 處紊動(dòng)能其中u′，v′，w′分別為x，y，z向的紊動(dòng)流速）的垂向分布。無(wú)植物（工況a）模擬中，RMSE較?。ǎ?.35 cm2·s-2），但由于紊動(dòng)強(qiáng)度較弱，模擬值與實(shí)測(cè)值的差距相對(duì)較大（IA＞0.62）。與之相比，植物工況b～工況d模擬中，對(duì)波高較小的W1和W2，模擬精度較高（如工況c-W2中，RMSE=0.87 cm2·s-2，IA=0.96）；而對(duì)波高較大的W3和W4，模擬精度較低（如工況d-W4中，RMSE=3.06 cm2·s-2，IA=0.42），且模擬值小于實(shí)測(cè)值。這可能是由于波高較大、水動(dòng)力作用較強(qiáng)時(shí)，物理模型試驗(yàn)中的植物會(huì)隨波浪的傳播發(fā)生輕微的晃動(dòng)；而在數(shù)值模擬中的植物被設(shè)定為固定不動(dòng)的理想圓柱體，未考慮其發(fā)生的晃動(dòng)。

圖6a 中，kT由上層水體至下層逐漸減小。無(wú)植物影響時(shí)，波浪傳播過(guò)程中，水面附近產(chǎn)生相對(duì)較大的紊動(dòng)，且紊動(dòng)由水面至底面逐漸減弱。與無(wú)植物條件相比，植物的存在使得水體的紊動(dòng)明顯增強(qiáng)，表現(xiàn)為工況b～工況d的kT顯著大于工況a。

如圖6b、6c所示，在統(tǒng)一高度植物影響下，尾流紊動(dòng)在植物的中部偏上位置處取得最大值，然后向冠層頂部和底部減小。在冠層下方，所有波浪條件下，工況c的kT模擬值均比工況b大30%～70%，說(shuō)明尾流紊動(dòng)強(qiáng)度與植物密度正相關(guān)。而在植物垂向密度發(fā)生變化的工況d中，一方面，由水面至底部，植物密度逐漸增加，尾流紊動(dòng)也隨之增強(qiáng)；另一方面，尾流紊動(dòng)的最大值出現(xiàn)在不同高度植物的中上部。兩種效果的疊加使得圖6d中kT模擬值的垂向梯度較小。

圖6 x=7.5 m處紊動(dòng)能kT垂向分布Fig.6 Vertical distributions of turbulent kinetic energy kT at x=7.5 m

選取工況c-W3條件下兩根植物之間的三個(gè)測(cè)點(diǎn)，分析一個(gè)波周期內(nèi)kT和流速u(mài)的相位變化。如圖7d所示，植物B位于植物A正前方5 cm，測(cè)點(diǎn)1位于植物A正前方7 mm，測(cè)點(diǎn)2位于兩根植物中點(diǎn)，測(cè)點(diǎn)3位于植物B正后方7 mm（植物A和B半徑均為4 mm，規(guī)定波浪傳播方向?yàn)檎?。?duì)比三個(gè)測(cè)點(diǎn)的縱向流速u(mài)（θ），測(cè)點(diǎn)2的流速幅值顯著大于測(cè)點(diǎn)1和3。當(dāng)水流流向植物時(shí)，在植物前緣由于其阻流作用水流流速迅速減小并向兩側(cè)分散開(kāi)來(lái)，而后在植物后方匯聚，形成尾流區(qū)，在尾流區(qū)至下一植物前緣區(qū)域流速逐漸增大至穩(wěn)定。由圖7d可知，當(dāng)波浪流速為正向時(shí)，測(cè)點(diǎn)1位于植物A的尾流區(qū)，測(cè)點(diǎn)3位于植物B的前緣；當(dāng)波浪流速為負(fù)向時(shí)，測(cè)點(diǎn)1位于植物A的前緣，測(cè)點(diǎn)3位于植物B的尾流區(qū)。因此，測(cè)點(diǎn)1和3始終處于流速減小區(qū)域，其流速幅值小于測(cè)點(diǎn)2。

如圖7a、7c所示，在測(cè)點(diǎn)1和3處，波周期內(nèi)紊動(dòng)能kT（θ）出現(xiàn)兩個(gè)峰值，且最大峰值顯著大于第二峰值。最大峰值與縱向流速u(mài)（θ）的峰值/谷值之間存在相位差，而第二峰值基本與u（θ）的峰值/谷值同時(shí)出現(xiàn)。在植物前緣（波浪流速為正時(shí)測(cè)點(diǎn)3，為負(fù)時(shí)測(cè)點(diǎn)1），植物的阻流作用使得該區(qū)域出現(xiàn)較大的流速梯度，進(jìn)而產(chǎn)生較大的紊動(dòng)，表現(xiàn)為測(cè)點(diǎn)1的u（θ）谷值和測(cè)點(diǎn)3 的u（θ）峰值對(duì)應(yīng)的kT（θ）第二峰值。而尾流區(qū)（波浪流速為正時(shí)測(cè)點(diǎn)1，為負(fù)時(shí)測(cè)點(diǎn)3）的紊動(dòng)需要經(jīng)歷在植物周?chē)纬赏暾曳€(wěn)定的渦旋后再脫落至水體中的過(guò)程，因此尾流區(qū)的kT（θ）峰值滯后于流速極值，表現(xiàn)為測(cè)點(diǎn)1 中kT（θ）最大峰值滯后于u（θ）峰值約30°，測(cè)點(diǎn)3中kT（θ）最大峰值滯后于u（θ）谷值約50°。

圖7 工況c-W3 條件下，一個(gè)波周期內(nèi)波浪流速u(mài)、紊動(dòng)能kT在不同測(cè)點(diǎn)處的相位圖和三個(gè)測(cè)點(diǎn)位置示意圖（y= 0.125 m，z=0.3 m，尺寸不按比例）Fig.7 Wave velocity u and turbulent kinetic energy kT versus wave phase at different measurement locations in test c-W3,and schematic of three measurement locations（y=0.125 m,z=0.3 m,not to scale）

2.2 泥沙懸浮分析

2.2.1 Flow-3D原始泥沙模塊模擬結(jié)果

采用Flow-3D原始泥沙模塊（見(jiàn)1.1.2節(jié)），模擬W3波浪條件下泥沙懸浮過(guò)程，結(jié)果如圖8所示。工況a 中，模擬的垂向泥沙濃度分布與實(shí)測(cè)值吻合較好（RMSE=0.019 kg·m-3，IA=0.99），說(shuō)明該模型可以較好地模擬出裸床（無(wú)植物）條件下的泥沙起動(dòng)及懸浮過(guò)程。然而，植物影響下的模擬結(jié)果卻與實(shí)測(cè)值存在極大差別。工況b～工況d的模擬值遠(yuǎn)小于實(shí)測(cè)值，三組的評(píng)估指標(biāo)RMSE＞0.21 kg·m-3，IA＜0.61。含植物水流數(shù)值模擬中，植物的存在減小了流速（圖4），由式（7）和（8）可知，床面切應(yīng)力、Shields數(shù)和泥沙挾帶上升速度均減小，泥沙起動(dòng)量降低，表現(xiàn)出圖8中含植物工況泥沙濃度模擬值小于無(wú)植物工況。而在物理模型試驗(yàn)中，植物的存在顯著增強(qiáng)了水體紊動(dòng)（圖6），進(jìn)而引起大量的泥沙懸浮，使得含植物工況泥沙濃度實(shí)測(cè)值大于無(wú)植物工況。由此可見(jiàn)，以切應(yīng)力進(jìn)行底部泥沙起動(dòng)及懸浮判斷的Flow-3D原始泥沙模塊無(wú)法精確模擬植物影響下的泥沙運(yùn)動(dòng)。

圖8 W3條件下x=7.5 m處泥沙濃度C垂向分布Fig.8 Vertical distributions of sediment concentration C at x=7.5 m in test W3

2.2.2 改進(jìn)的Flow-3D泥沙模塊模擬結(jié)果

為精確模擬含植物床面泥沙懸浮過(guò)程，利用基于水體紊動(dòng)能改進(jìn)的泥沙模塊（見(jiàn)1.2節(jié)），對(duì)各工況下的泥沙運(yùn)動(dòng)過(guò)程進(jìn)行計(jì)算。在計(jì)算過(guò)程中，基于對(duì)本研究水槽試驗(yàn)數(shù)據(jù)的分析，對(duì)式（17）中參數(shù)δ在床面附近取值為0.63較為合理［16］，該值與Zhang等［10］的推薦值（0.76）相近；對(duì)式（17）中參數(shù)Cw取值為0.015，該值在Tang等［15］的推薦范圍內(nèi)（0.02±0.01）。

圖9 為泥沙改進(jìn)模塊的計(jì)算結(jié)果，各工況下模擬結(jié)果均較好（RMSE＜0.1 kg·m-3，IA＞0.89），說(shuō)明該模型能夠合理預(yù)測(cè)波浪條件下含植物床面泥沙起動(dòng)及懸浮過(guò)程。雖然流速Uw和紊動(dòng)能kT的模擬結(jié)果存在部分偏差（圖4、圖6），但泥沙模擬值與實(shí)測(cè)值吻合較好，IA值相對(duì)較高。這主要有兩方面原因：一方面，物理模型試驗(yàn)中，用三維聲學(xué)多普勒測(cè)速儀（ADV）測(cè)量流速參數(shù)時(shí)（采樣頻率64 Hz）對(duì)環(huán)境噪聲的敏感性較高，且易受植物晃動(dòng)影響；而光學(xué)后向散射濁度傳感器（OBS）測(cè)量泥沙濃度時(shí)（采樣頻率8 Hz）不易受環(huán)境噪聲影響，精度較高。另一方面，泥沙輸運(yùn)的模擬效果很大程度取決于近床面泥沙的起動(dòng)及懸浮，而物模試驗(yàn)中植物晃動(dòng)主要影響上層水體的紊動(dòng)，對(duì)近床面的紊動(dòng)影響很微弱。因此，數(shù)值模擬中，在式（17）中選擇合理的參數(shù)δ、CD和Cw使得改進(jìn)后的泥沙模型充分考慮到床面剪切紊動(dòng)和植物尾流紊動(dòng)后，可得到較好的泥沙輸運(yùn)模擬效果。

圖9 x=7.5 m處泥沙濃度C垂向分布Fig.9 Vertical distributions of sediment concentration C at x=7.5 m

圖9 中，模擬含植物床面（工況b～工況d）的泥沙濃度顯著大于無(wú)植物床面（工況a），且波高越大越明顯，與實(shí)測(cè)數(shù)據(jù)一致，證明剛性植物引起的尾流紊動(dòng)對(duì)泥沙起動(dòng)及懸浮的促進(jìn)作用十分顯著。在近底層（z/hv＜0.5），工況c的泥沙濃度顯著大于工況b，主要是因?yàn)楦呙芏戎参镆鸬奈擦魑蓜?dòng)更強(qiáng)，使得更多的泥沙起動(dòng)。而在水體上部（z/hv＞0.5），工況b和c的泥沙濃度垂向分布結(jié)構(gòu)相近，主要是因?yàn)閮煞N工況下紊動(dòng)能kT垂向分布結(jié)構(gòu)相似（圖6b、6c）。在工況d 中，底層水體的泥沙濃度略大于工況c，而上層水體則相反。雖然工況c和d植物密度相同，但兩者的kT垂向分布結(jié)構(gòu)差異較大（圖6c、6d），導(dǎo)致了不同的泥沙濃度垂向分布規(guī)律。

3 結(jié)語(yǔ)

本文基于Flow-3D 構(gòu)建了含剛性沉水植物的三維水槽模型，改進(jìn)了泥沙計(jì)算模塊，模擬了不同波高的波浪通過(guò)植物區(qū)時(shí)水動(dòng)力特征及泥沙輸運(yùn)過(guò)程，利用均方根誤差（RMSE）和一致性指數(shù)（IA）證明模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)較為吻合。主要結(jié)論如下：

（1）植物的存在使得波浪流速減小，其衰減幅度與植物密度正相關(guān)。在垂向密度均勻的植物群中，冠層頂部正向流速被增大；而在垂向密度不均勻的植物群中，由于不同高度植物的相互影響，該現(xiàn)象并不明顯。

（2）植物的存在顯著增強(qiáng)了水體紊動(dòng)，其增幅與植物密度正相關(guān)。植物尾流紊動(dòng)在植物中上部取得最大值。對(duì)單根植物而言，在其植物前緣和后方尾流區(qū)均屬于流速減小區(qū)。當(dāng)波浪通過(guò)植物區(qū)時(shí)，在每根植物附近一個(gè)周期內(nèi)水體紊動(dòng)能存在兩個(gè)峰值，最大峰值與尾流區(qū)的渦旋有關(guān)，且與流速的峰值/谷值之間存在相位差；第二峰值與植物前緣的阻流作用有關(guān)，且與流速的峰值/谷值同步。

（3）在Flow-3D原始泥沙模塊中，僅利用床面切應(yīng)力計(jì)算泥沙挾帶上升速度。將其應(yīng)用于含植物床面時(shí)，沒(méi)有考慮到植物尾流紊動(dòng)對(duì)泥沙懸浮的促進(jìn)作用，因此模擬出的泥沙懸浮量顯著偏小。本研究從紊動(dòng)能的角度改進(jìn)了泥沙模塊，在模型中即考慮植物尾流紊動(dòng)，也考慮床面剪切紊動(dòng)。改進(jìn)后對(duì)泥沙懸浮的模擬評(píng)價(jià)指標(biāo)由RMSE＞0.21 kg·m-3，IA＜0.61 提升為RMSE＜0.1 kg·m-3，IA＞0.89，顯著提高了含植物床面泥沙懸浮的模擬精度。