郝由之，假冬冬,2*，張幸農(nóng)，陳長英，吳 磊，楊 俊

(1.南京水利科學(xué)研究院 港口航道泥沙工程交通行業(yè)重點實驗室，江蘇 南京 210029；2.長江保護(hù)與綠色發(fā)展研究院，江蘇 南京 210098)

天然河流大多是由主槽和漫灘組成的復(fù)式河道，枯水期主槽承擔(dān)泄流功能，洪水期水流溢出主槽，由主槽和漫灘共同泄流[1]。而在自然河流近岸臨水區(qū)，即主槽岸坡，經(jīng)常會觀察到喬木和灌木等植被。傳統(tǒng)觀點認(rèn)為，岸坡植被不僅可以減緩水流速度，減弱河岸沖刷侵蝕，避免崩岸頻發(fā)，從而穩(wěn)定河勢，還能起到凈化水體、美化環(huán)境的作用，是一種既經(jīng)濟(jì)又有效的生態(tài)治河措施[2]。然而在洪水期，岸坡植被在一定程度上會改變?yōu)┎蹆?nèi)部原有水流結(jié)構(gòu)，同時植被和水流間的相互作用會使河道中部分水流的能量轉(zhuǎn)化為植被附近產(chǎn)生的紊動能，從而降低河道行洪能力。因此，了解復(fù)式河道岸坡植被對水流的擾動影響，進(jìn)而分析岸坡植被的利害關(guān)系，對于河道防洪、護(hù)岸工程設(shè)計及水生態(tài)保護(hù)與修復(fù)均具有重要的研究意義[3]。

關(guān)于植被水流研究，在水槽試驗方面，Nepf[4]利用示蹤劑顯示了圓形非淹沒植被陣列內(nèi)的紊流形態(tài)變化。López[5]、Nezu[6]和閆靜[7]等以木棒或條狀薄片模擬淹沒剛性植被，研究了水流的紊動特性。趙芳[8]和楊克君[9]等研究了不同植被類型條件下河道的水流特性。隨著數(shù)值計算理論和計算機(jī)技術(shù)的發(fā)展，解析方法和數(shù)值模擬技術(shù)得到成功的應(yīng)用，且與試驗相比，易得到整個計算域的水流特性。Huai等[10]提出了預(yù)測平均流速垂向分布的分析模型，很好地模擬了含雙層剛性植被河道水流的垂向分布。Wu[11]和Zhang[12]等在動量方程和k-ε模型中考慮水流對植被的阻力影響，建立了植被水流2維模型。但是，含植被的河道水流具有復(fù)雜的3維各向異性，因此，3維模型的建立非常必要。吳夢瑤等[13]將多孔介質(zhì)理論應(yīng)用于自由表面水流3維模型，較好地反映了有限植被群的尾流區(qū)流動特征。劉誠[14]將Wu[11]等將模型擴(kuò)展到非正交曲線坐標(biāo)系下建立了3維模型。以上3維模擬方法可實現(xiàn)植被對水流的宏觀效應(yīng)模擬，但對植被群內(nèi)紊流特性的捕捉仍存在一定困難?；蔽男臶15]和Nadaoka[16]等采用大渦模擬技術(shù)建立了植被水流模型，較好地模擬了植被尾流和植被群內(nèi)細(xì)部水流結(jié)構(gòu)和大渦的發(fā)展過程，但當(dāng)計算域較大，而細(xì)部渦流復(fù)雜區(qū)域占比相對較小時，采用大渦模擬運算量較大。

復(fù)式河道受自身斷面形態(tài)特點影響，水流特性復(fù)雜，再加上河道內(nèi)植被與水體間相互作用，水體流動更加復(fù)雜。然而，當(dāng)前研究大多以河床或漫灘植被為研究對象，如：米云彤[17]和任俊韜[18]等考慮了河床植被對水流特性的影響；蔣北寒[19]和袁素勤[20]等探究了漫灘植被對流速分布的影響。而實際河流中，植被在復(fù)式河道岸坡廣泛分布，其對河道行洪的影響尚不明確。以往研究也多以單個植被或植被斑塊為研究對象，而對順?biāo)鞣较蛑脖谎爻谈采w的研究較少，且現(xiàn)有研究結(jié)果主要集中于植被對河道流速分布和紊流特性的影響，對二次流結(jié)構(gòu)、河道流量分配及河床剪應(yīng)力分布影響的研究尚不多見，而同時考慮這些影響能夠更加全面認(rèn)識植被對水流的擾動。因此，岸坡植被對復(fù)式河道水動力特性的影響還需進(jìn)一步研究。

本文將植被作為固壁邊界，并采用雷諾應(yīng)力模型封閉紊流控制方程，建立3維植被河道水動力數(shù)學(xué)模型，相比于2維模型，可以獲取河道橫斷面的水流特性；同時，與傳統(tǒng)研究在動量方程中添加植被阻力項相比，本文方法可有效捕捉植被區(qū)域內(nèi)部的紊流結(jié)構(gòu)，同時又避免了使用大渦模擬時的高運算量。為了驗證模型的準(zhǔn)確性，對復(fù)式河道的水流特性進(jìn)行了模擬，并將結(jié)果與Tominaga和Nezu[21]的試驗數(shù)據(jù)進(jìn)行比較。隨后，利用該模型模擬了復(fù)式河道岸坡種植不同淹沒狀態(tài)的植被和不同河岸坡比條件下的河道水流結(jié)構(gòu)，以揭示岸坡植被對復(fù)式河道水流擾動的影響。

1 方 法

1.1 數(shù)學(xué)模型構(gòu)建

1.1.1 控制方程

利用CFD計算流體力學(xué)軟件建立基于有限體積法的3維植被河道水動力數(shù)學(xué)模型，其中紊流控制方程采用連續(xù)性方程和雷諾時均納維-斯托克斯方程[8]，并通過雷諾應(yīng)力模型對紊流模型進(jìn)行封閉。

對于不可壓縮均勻流體，連續(xù)性方程與雷諾時均納維-斯托克斯方程形式如下：

式（1）～（2）中：i,j= 1, 2, 3為空間指數(shù)，x1、x2、x3分別代表x、y、z方向，即橫向、縱向、垂向；ρ為水的密度，取998.2 kg/m3；μ為水的動力黏滯系數(shù)，取1.002×10-3Pa·s；p為時均壓力；Ui（Uj）和ui（uj）分別為在xi（xj）方向上的時均流速和瞬時流速；為雷諾應(yīng)力項；gi為重力加速項。

紊流封閉采用雷諾應(yīng)力模型，該模型通過求解雷諾應(yīng)力輸運方程對方程組進(jìn)行封閉[22]。雷諾應(yīng)力輸運方程具體形式如下：

其中：

對流項為：

式中：xk為垂直于xi和xj的方向；uk為xk方向的瞬時流速。

紊動擴(kuò)散項可通過Daly等[23]給出的廣義梯度擴(kuò)散模型進(jìn)行計算：

分子黏性擴(kuò)散項為：

剪應(yīng)力產(chǎn)生項為：

壓力應(yīng)變項為：

式中，p′為單位時間內(nèi)的壓力。

黏性耗散項為：

雷諾應(yīng)力輸運方程中，i、j取1、2、3時，得到3個正應(yīng)力方程，將3個方程相加得到紊動能的準(zhǔn)確方程；而方程式中除了剪應(yīng)項外，均含有2階或3階相關(guān)矩，因此，需給定邊界條件使雷諾應(yīng)力輸運方程封閉，才能對紊動能、紊動耗散率等參數(shù)進(jìn)行求解。

1.1.2 邊界條件

模型計算域及邊界條件如圖1所示，邊界條件類型主要包含速度進(jìn)口邊界、自由流出口邊界、對稱邊界和壁面邊界。復(fù)式河道中，水面線幾乎與底坡平行，自由水面高程的起伏變化很小，在水深較大的情況下可以用對稱邊界條件近似代替自由水面，因此，自由水面采用剛蓋假定法做近似處理。模型進(jìn)口處設(shè)定為水流速度進(jìn)口邊界條件，并給定平均流速Uave。出口處設(shè)為自由出流邊界條件，即出口處沿水流流向的速度梯度和紊動強(qiáng)度梯度為零。由于河道形態(tài)沿主槽中心對稱分布，為了減少計算工作量，計算區(qū)域選取整個斷面的一半，將主槽中心縱剖面設(shè)為對稱邊界。河床底部、灘地邊界及剛性植被均設(shè)為壁面邊界，黏性流動中壁面默認(rèn)為靜止無滑移邊界條件，即Ui= 0。由于Tominaga等[21]研究發(fā)現(xiàn)相對于斷面形態(tài)與水力條件，壁面糙率對水流結(jié)構(gòu)的影響十分微弱[24]，因此，該模型中未考慮壁面糙率的影響。

圖1 模型計算域及邊界條件示意圖Fig.1 Schematic diagram of model calculation domain and boundary conditions

1.1.3 植被固壁邊界處理方法及網(wǎng)格劃分

網(wǎng)格劃分方法和網(wǎng)格質(zhì)量對模擬結(jié)果有很大的影響，因此合理的網(wǎng)格劃分可以提高模型的計算精度和計算效率。針對本文模型的特點，要保證在靠近河床和植被邊壁的流速大梯度變化區(qū)域網(wǎng)格有較細(xì)的間距，而在主槽和漫灘中心的低速梯度區(qū)域網(wǎng)格有更寬的間距。此外，主槽和漫灘區(qū)水平面上劃分為均勻的笛卡爾網(wǎng)格塊，而岸坡斜面上劃分為向植被圓柱體拉伸的“O”型網(wǎng)格塊。此劃分方法能夠在保證捕捉植被區(qū)細(xì)部水流結(jié)構(gòu)的同時減小網(wǎng)格數(shù)量，節(jié)省計算資源。模擬時采用了不同分辨率的網(wǎng)格驗證了模型的收斂性和模擬精度，進(jìn)而選取適當(dāng)?shù)挠嬎憔W(wǎng)格。以岸坡種植非淹沒剛性植被為例，其3維幾何模型不同剖面的網(wǎng)格劃分如圖2所示。

圖2 幾何模型在不同剖面的網(wǎng)格劃分Fig.2 Grid division of the geometric model in different sections

圖2中：近床面和岸坡植被區(qū)網(wǎng)格細(xì)密，主槽和漫灘中心網(wǎng)格較粗；同時，在x-z剖面和植被區(qū)水平剖面圖中，x1∶x2、z1∶z2及y1∶y2應(yīng)小于網(wǎng)格尺度比1∶5。

1.1.4 求解方法

采用基于有限體積法的計算流體力學(xué)軟件建立含植被河道3維紊流數(shù)值模型，方程離散時，壓力速度耦合方式選擇SIMPLE算法，壓力的插值格式采用體積力加權(quán)格式，動量等其余量的離散格式采用2階迎風(fēng)插值格式。當(dāng)各計算參量殘差小于0.000 01且監(jiān)測某點的參量不隨迭代次數(shù)變化時，可認(rèn)為計算收斂，迭代終止。

1.2 模型驗證

采用Tominaga和Nezu[21]矩形復(fù)式河道紊流試驗中工況 S-2的結(jié)果對所建數(shù)學(xué)模型進(jìn)行驗證。該試驗中，水槽長12.5 m，寬0.4 m，高0.4 m；水槽底板為鐵板，邊壁為有機(jī)玻璃。河道主槽水深hm為0.08 m，灘地水深hf為0.04 m，主槽寬度b為0.2 m，總寬度B為0.4 m，水力參數(shù)見表1。

表1 驗證工況相關(guān)水力參數(shù)Tab.1 Related hydraulic parameters of verification conditions

選擇河道中心橫斷面作為代表斷面，繪制了二次流、主流速及河床剪切應(yīng)力分布，如圖3～5所示。圖3～5中，橫、縱坐標(biāo)軸均使用主槽水深hm進(jìn)行無量綱處理，主流速使用最大流速Umax進(jìn)行無量綱處理，河床剪切應(yīng)力 τ采用平均剪切應(yīng)力 τ進(jìn)行無量綱處理。

圖3 二次流分布驗證Fig.3 Secondary flow distribution verification

由圖3可知，試驗和模擬結(jié)果在復(fù)式河道灘槽交接處均存在一對方向相反的二次流漩渦，二次流漩渦的影響范圍從灘槽交界處開始延伸至自由表面，水平方向的影響范圍也一致且漩渦的位置、大小和強(qiáng)度十分接近，故擬合效果較好。

由圖4可知，試驗與模擬結(jié)果規(guī)律相近，只在自由水面處略有誤差，這是由于數(shù)值計算時將自由水面作為對稱面邊界處理，無法較真實地模擬自由表面渦，導(dǎo)致近水面流速誤差相對較大。但數(shù)值模擬與試驗結(jié)果均顯示，復(fù)式河道灘槽交界處的主流速等值線向主槽一側(cè)的自由水面凸起，灘槽交界區(qū)主流速有明顯減小趨勢，說明數(shù)值模擬結(jié)果仍能較好地反映斷面主流速分布的主要特點。

圖4 主流速分布驗證Fig.4 Main flow velocity distribution verification

由圖5可知，復(fù)式河道河床剪切應(yīng)力分布的數(shù)值模擬與試驗結(jié)果基本一致，均在灘槽交界處剪應(yīng)力迅速減小并達(dá)到最小值，并在交界處向主河道區(qū)和漫灘區(qū)兩個方向過渡時剪切應(yīng)力值迅速增加至最大值，然后再分別減小，直至邊界處達(dá)到最小值。

圖5 河床剪切應(yīng)力分布驗證Fig.5 River bed shear stress distribution verification

1.3 工況設(shè)置

模型設(shè)置為復(fù)式梯形河道，如圖2所示，模擬類似洪水期漫灘水流河道，設(shè)置主槽寬度為1.5 m，漫灘寬度為2.0 m，岸坡高度為0.4 m，水深0.8 m，流向長度30.0 m，河道坡降為0.1%。岸坡順?biāo)鞣较蜓爻滩荚O(shè)剛性植被，植被以平行方式排布，間距為0.14 m，植被直徑為20 mm。模擬工況設(shè)置了不同的植被淹沒方式及不同的河岸坡比（河岸橫向?qū)挾扰c高度之比）進(jìn)行對比，所有模擬工況中的進(jìn)口平均流速Uave均取0.5 m/s，具體工況設(shè)置見表2。

表2 數(shù)值模擬工況Tab.2 Numerical simulation conditions

2 結(jié)果與討論

2.1 流速分布與流量分配

圖6為岸坡種植剛性非淹沒植被時（z=0.6 m）處的河道水平面流速分布。由圖6可知，左側(cè)進(jìn)口處流場為非充分發(fā)展階段，該區(qū)域流速分布與單個圓柱繞流相似[25]，也說明了模擬方法可以很好地捕捉植被區(qū)細(xì)部水流特點，具有一定的精度。由于本文考慮植被在岸坡沿程覆蓋的情況，因此重點關(guān)注水流充分發(fā)展階段的水流結(jié)構(gòu)，故下文河道斷面水流特性數(shù)據(jù)也均取自充分發(fā)展階段的斷面進(jìn)行分析。

圖6 水平面流速分布Fig.6 Velocity distribution on the horizontal surface

圖7為8種工況下的復(fù)式河道橫斷面流速分布。

圖7 斷面主流速分布Fig.7 Distribution of main flow velocity in section

對比圖7（a）～7（d）發(fā)現(xiàn)：岸坡植被使得岸坡區(qū)流速顯著減小，且流速減小區(qū)域與植被分布相關(guān)；其中，剛性非淹沒植被所占岸坡范圍最大，因此對灘槽交界區(qū)流速的減小作用最強(qiáng)。

圖7（e）～7（h）分別為1∶1.5和1∶2.0坡比條件下，岸坡有、無植被時的斷面主流速分布，與圖7（a）和7（b）對比可發(fā)現(xiàn)：無植被情況下，主流速在灘槽交界處凸起，這與Xiao等[26]對不同坡比下復(fù)式河道水流研究得出的結(jié)果相似，這是由于此處邊界凸起，引起了二次流的產(chǎn)生；但當(dāng)岸坡植被存在時，灘槽交互區(qū)流速仍顯著減小。

岸坡植被改變了河道斷面流速分布，進(jìn)而影響流量在主槽、河岸和漫灘3個區(qū)域的分配，圖8為工況1～4的斷面流量分配占比。由圖8可知，在來流量相同條件下，岸坡植被的存在大大削弱了岸坡區(qū)過流量，但主槽和漫灘區(qū)過流量增加，這在一定程度上增加了河道的過流負(fù)荷。其中，河岸種植非淹沒剛性植被影響最顯著，與岸坡無植被相比，岸坡區(qū)流量分配由9.21%縮減到2.07%，減小了80%，使得主槽區(qū)流量分配增加了6%，漫灘區(qū)流量分配增加了13%，這些結(jié)果可為確定河道行洪負(fù)荷提供依據(jù)。

圖8 斷面流量分配Fig.8 Cross-section flow distribution

2.2 二次流分布

圖9給出8種工況下復(fù)式河道的二次流分布，橫、縱坐標(biāo)軸均使用主槽水深hm進(jìn)行無量綱處理。

圖9 斷面二次流分布Fig.9 Secondary flow distribution in section

圖9（a）～（d）顯示：同一坡比條件下，岸坡無植被時，灘槽交界區(qū)存在一對方向相反的二次流渦旋。當(dāng)岸坡種植非淹沒剛性植被時，二次流渦旋數(shù)量增加，灘槽交界區(qū)有兩對方向相反的二次流渦旋；當(dāng)岸坡交替種植淹沒與非淹沒剛性植被時，灘槽交界區(qū)有3個相鄰且方向相反的二次流渦旋；當(dāng)岸坡種植淹沒剛性植被時，灘槽交界區(qū)有一對方向相反的二次流渦旋，但相比無植被時，二次流渦旋水平影響范圍增加。此外，根據(jù)二次流流速矢量圖例可發(fā)現(xiàn)，岸坡有植被時的二次流約為無植被情況的3倍，說明岸坡植被可以增加二次流渦旋的數(shù)量、強(qiáng)度及影響范圍，其中非淹沒剛性植被影響最顯著。

由圖9（a）、（e）及（g）可見，當(dāng)河岸坡比逐漸減小時，渦旋水平方向的影響范圍進(jìn)一步減小，二次流的形狀也逐漸模糊，這是由于坡角變緩后，主槽與灘地之間存在一個流速過渡區(qū)，灘槽間橫向動量交換減弱，二次流強(qiáng)度也隨之減弱。與圖9（b）、（f）及（h）對比同樣可發(fā)現(xiàn)，岸坡植被能夠促使二次流渦旋的產(chǎn)生，并擴(kuò)大二次流的影響范圍，增加灘槽間橫向動量交換，但隨著坡比的減小，植被的影響會隨之減弱。

2.3 河床剪應(yīng)力分布

圖10給出了8種工況下，對平均剪應(yīng)力 τ進(jìn)行無量綱化后的河床剪應(yīng)力分布。圖10中，河床剪切應(yīng)力τ可由CFD計算后直接得出，再將 τ沿床面平均即可求得平均剪應(yīng)力 τ。

圖10 河床剪應(yīng)力分布Fig.10 River bed shear stress distribution

由圖10可見：各工況下河床剪應(yīng)力分布趨勢基本一致，均在主槽和灘地中心接近均勻分布，主槽和漫灘向岸坡過渡時河床剪應(yīng)力逐漸增大，進(jìn)入河岸邊坡區(qū)又發(fā)生突變減??；但不同工況下，主槽和灘地中心河床剪應(yīng)力值及灘槽交界區(qū)河床剪應(yīng)力的變化幅度有所不同。

對比工況1～4可以發(fā)現(xiàn)：同一坡比條件下，岸坡無植被時，主槽和漫灘河床剪應(yīng)力最小，無量綱化剪應(yīng)力值約為1；岸坡種植植被后，主槽和漫灘河床剪應(yīng)力均有所增大，可增大至1.1～1.3。而對于岸坡區(qū)無植被時，河床剪應(yīng)力最大，無量綱化剪應(yīng)力值約為0.7～1.0，該結(jié)果與Xiao等[26]對坡比1∶1復(fù)式河道河床剪應(yīng)力的模擬結(jié)果基本吻合；有植被時，河床剪應(yīng)力大幅減小，甚至趨近于零。在主槽和漫灘向岸坡過渡區(qū)，岸坡種植非淹沒剛性植被時，河床剪應(yīng)力增加幅度最大，交替種植淹沒與非淹沒剛性植被時次之，種植淹沒剛性植被時最小。這說明非淹沒剛性植被能夠引起強(qiáng)烈的灘槽動量交換，從而使得該處河床剪應(yīng)力急劇增加。

對比工況1～2和5～8可以發(fā)現(xiàn)，在岸坡有、無植被情況下，隨著坡比的減小，主槽、漫灘及灘槽交界區(qū)的河床剪應(yīng)力均逐漸增大，說明隨著坡比減小，二次流引起的動量交換強(qiáng)度變小，對邊界剪切應(yīng)力分布的影響減弱，河床剪應(yīng)力分布向均勻化趨近。

2.4 雷諾應(yīng)力分布

由圖11（a）可知，岸坡無植被時，雷諾應(yīng)力在岸坡區(qū)呈凸起分布。與圖11（b）～（d）對比可以發(fā)現(xiàn)：岸坡種植植被后，灘槽交界區(qū)雷諾應(yīng)力分布與無植被時有很大差異，雷諾應(yīng)力分布不再呈凸起狀態(tài)，而是與岸坡植被分布相關(guān)；岸坡主槽側(cè)存在一個緊鄰植被的雷諾應(yīng)力值增大區(qū)，且為正值；岸坡漫灘側(cè)存在一個緊鄰植被的雷諾應(yīng)力值減小區(qū)，且為負(fù)值。這說明岸坡植被的存在顯著影響了雷諾應(yīng)力分布的變化，使得臨近岸坡的主槽和漫灘與岸坡之間的橫向動量交換更加劇烈。

圖11 斷面雷諾應(yīng)力分布Fig.11 Reynolds stress distribution of section

2.5 紊動能分布

圖12為1∶1坡比條件下4種工況的復(fù)式河道橫斷面紊動能分布，其中，紊動能k由平均摩阻流速U*進(jìn)行無量綱化處理。

圖12 斷面紊動能分布Fig.12 Distribution of section turbulent kinetic energy

由圖12可知：無植被時，灘槽交界區(qū)紊動能呈凸起分布；而岸坡種植植被后，同雷諾應(yīng)力分布相似，該區(qū)域紊動能不再呈凸起狀態(tài)，而是與岸坡植被分布相關(guān)。岸坡種植非淹沒剛性植被時，紊動能在岸坡植被區(qū)兩側(cè)顯著增加，無量綱化的紊動能值可達(dá)到20；岸坡交替種植淹沒與非淹沒剛性植被時，紊動能也在岸坡植被區(qū)兩側(cè)顯著增加，但幅度比前者略??；岸坡種植淹沒剛性植被時，紊動能在岸坡植被區(qū)兩側(cè)及植被頂部顯著增加，增加幅度也小于前者。這說明岸坡植被的存在能夠增加緊鄰植被周圍區(qū)域的紊動交換，而植被區(qū)域內(nèi)部紊動能會顯著減小。米云彤等[17]對植被群周圍水流特性進(jìn)行測定，其植被群處斷面紊動能結(jié)果與本文結(jié)果趨勢相同。

3 結(jié) 論

本文采用雷諾應(yīng)力模型封閉紊流模型，建立了3維植被河道水動力數(shù)學(xué)模型，在模型驗證相似的基礎(chǔ)上，利用所建模型研究了岸坡剛性植被對水流的影響，主要研究結(jié)論為：

1）岸坡植被增加灘槽間動量交換，使得岸坡區(qū)的流速、流量及河床剪應(yīng)力減小，而二次流渦旋數(shù)量、強(qiáng)度和范圍均有所增加；此外，隨著河岸坡比減小，無植被河道灘槽間動量交換逐漸衰減，但在坡面種植非淹沒剛性植被后，灘槽間動量交換的衰減程度有所減小。植被區(qū)域內(nèi)紊動有所減弱，而臨近植被的周邊區(qū)域紊動交換劇烈，且岸坡種植非淹沒剛性植被時影響最顯著。

2）根據(jù)斷面流量分配結(jié)果，可確定岸坡植被對行洪負(fù)荷的影響，而河床剪應(yīng)力結(jié)果又可與河岸土體起動建立聯(lián)系，進(jìn)一步為坡腳沖刷及河岸穩(wěn)定提供判據(jù)，因此該模擬方法和結(jié)果有助于評估和權(quán)衡植被的利弊影響，進(jìn)而制定適宜的河岸帶生態(tài)護(hù)岸植被配置方案，使得植被發(fā)揮到最大的生態(tài)效益和工程安全效益。