剛性植被密度對岸灘作用的數值模擬*

肖苡辀，張國志，馮先導，陳文峰

(1.中交第二航務工程局有限公司，湖北 武漢 430040；2. 長大橋梁建設施工技術交通行業重點實驗室，湖北 武漢 430040；3. 交通運輸行業交通基礎設施智能制造技術研發中心，湖北 武漢 430040；4. 中交公路長大橋建設國家工程研究中心有限公司，湖北 武漢 430040)

岸灘修復工程通常會利用植被來增加岸坡的穩定性，為水生動物提供遮蔭和避難所。水生植被常存在于河岸、海岸邊，可顯著改變局部流場、泥沙輸移、河岸穩定性和水生棲息地[1-5]。了解植被對水流、泥沙輸移的影響是保護港口海岸岸線的關鍵。

近年來，學者對植被層內外的水流懸移質泥沙輸移規律及床面推移質泥沙起動特性進行了系統研究。張英豪等[6]通過試驗研究了植物對水流和泥沙運動的影響、植物的消波與捕沙作用，認為合理地布置植被可以減少泥沙輸運；Vastila等[7]通過試驗發現清水中的泥沙沉降與蘆葦植被水流中的泥沙沉降基本一致；薛萬云等[8]對剛性植被環境中床面形態變化進行了試驗研究，結果表明淹沒植被條件下動床床面泥沙運動可分為靜止、部分起動及全部起動3個階段。

綜上，植被對岸灘的穩定性分析主要采取試驗的研究方法，缺乏數值模擬植被密度對促淤保灘作用方面的研究。本文選取圓柱代替剛性挺水植被，采用數值模擬方法探索植被水流中床面變化形態規律，分析植株周圍水流流速、床面切應力與植被密度的關系，對促進植被-水流-泥沙研究的發展具有重要價值。

1 數值模擬方法

1.1 控制方程

假定水流為不可壓縮黏性流體，則水流計算可運用Navier-Stokes方程描述[9]：

(1)

(2)

式中：ρ為流體密度；ui為局部時均速度分量；u為速度矢量；t為時間；μ為流體動力黏滯系數；p為壓強；Fi為微元體上的體力。

本研究采用的泥沙沖刷模型基于文獻[10]。Soulsby-Whitehouse方程[11]用于計算臨界Shields數θcr，i：

(3)

(4)

式中：ρi為第i相泥沙密度；ρf為流體密度；di為顆粒直徑；μf為流體動力黏滯系數；g為重力加速度。

局部的Shields參數θi是基于局部床面切應力τ計算的：

(5)

其中τ考慮床面粗糙度，采用壁面定律進行計算。

利用夾帶系數(0.005)來衡量沖刷速率，并與試驗數據進行擬合。沉降速度控制著Soulsby沉降方程。公式(6)～(8)用來計算河床單寬的體積輸沙率：

(6)

Φi=βMPM，i(θi-θcr，i)1.5cb，i

(7)

(8)

式中：βVR，i、βMPM，i和βNie，i一般為0.053、8.0和12.0；cb，i為第i相的體積分數；Φi為無量綱推移質輸沙率，假定為0.05。容積床載輸送率qb，i表示為：

(9)

根據粗糙度修正后的壁面平均速度定律為：

(10)

式中：up是在某點j處的主流流速；k為von Karman常數，為0.418；u為與底部剪切應力有關的剪切速度；cμ為常數；kj為點j處的紊動能；ΔB為粗糙度函數；zj為點j的垂向高度。

1.2 流場設置與邊界條件

本研究基于袁代亮[12]的試驗設置建立模型，通過CFD軟件計算流場，旨在驗證模型的準確性。該驗證模型長3.0 m、寬0.5 m，水深0.24 m，流速0.227 8 m/s，主要監測其垂向流速分布。在模型得到驗證后進行不同植被密度下的水流-泥沙計算(圖1)，具體工況參數見表1。計算區域原點設于水槽進口底板中心點，整體計算區域長3.0 m、寬0.5 m、高0.5 m。剛性植被采用圓柱代替，直徑為0.1 m，第1排植被(植被圓心)位于距離水槽進口1 m，植被縱向間距為0.1 m，植被橫向間距依據植被密度不同而不同(從工況1～4，橫向間距分別為0.125、0.100、0.083、0.071 m)。泥沙中值粒徑為0.006 m，密度為2 650 kg/m3，厚度為0.23 m，鋪設在水槽床面上。

圖1 計算模型

表1 試驗工況設計

模型采用流量進口邊界作為進口邊界條件，壓力出口邊界作為出口邊界條件。壁面邊界用于水槽底板和邊壁。渠道頂部空氣入口設定為對稱邊界，即默認無流體穿過該邊界。

1.3 計算方法

基于有限差分法采用RNGk-ε紊流模型來求解控制方程，對流項采用二階迎風格式，擴散項采用二階中心差分格式。采用TruVOF方法縮減收斂時間，該方法僅計算流體單元，不計算空氣單元，且對自由液面的描述更為準確。通過FAVOR法(fractional area volume obstacle representation)對計算區域劃分網格。當控制方程的殘差小于1×10-5時，認為計算結果收斂。

2 結果分析

2.1 試驗驗證

為了更好地研究植被存在時水動力結構和懸沙運動特性，需要對數值模擬的可靠性進行驗證。基于袁代亮的試驗結果進行數值模擬，從而對比分析數值模擬的可靠性。本研究基于試驗設置和結果進行數值模擬，結果見圖2。圖2為垂向流速分布模擬值與試驗值對比，結果表明模擬值與試驗值吻合度較好，該數值模擬方法能夠有效模擬植被水流工況。

注：z為垂向高度，H為該工況水深。圖2 垂向流速分布模擬值與試驗值對比

2.2 流速分布

基于速度場的平均水力特性研究植被對水流結構的影響。圖3為各工況水平截面z/D=15的速度分布。由圖3可知，植被的存在阻礙了水流的流動，而植被間的流速則增大。值得注意的是，從工況1～3，隨著植被密度的增加，流速逐漸增大；但比較工況3和4，植被密度的增加反而降低了流速。這個現象是在先前研究中出現過的[13]，其原因是工況1～3中，卡門渦街在動量運輸中占有主導作用；而在工況3和4中，卡門渦街已經被密度較高的植被耗散了，以至于渦街并未達到下一排植被。

圖3 各工況水平截面zD=15的速度分布

2.3 床面切應力

床面切應力影響著泥沙輸移，因此在工程中具有重要意義。床面切應力τ計算公式如下：

(11)

式中：Ud是水深平均流速；Cf是表面摩擦系數。

圖4呈現了各工況橫斷面和工況4沿程橫斷面床面切應力分布情況，其中圖4a)為工況4的監測斷面分布，由虛線表示。由圖4可以看到，植被水流中的床面切應力分布呈現波浪形狀，這說明植被的存在降低了平均流速，消耗了水流動量，增加了阻力，從而降低床面切應力。由圖4b)可知，從工況1～4，床面切應力隨著植被密度的增加而出現先增大再減小的趨勢，其中工況3的床面切應力最大。由此可知，一味地增大植被密度并不一定能夠保護好灘地，而是要到一定密度之后才會降低床面切應力從而保護床面。

由圖4c)可知，工況4斷面的床面切應力呈現沿程逐漸增加的規律。在植被群上游斷面(斷面x/D=93)，床面切應力分布較為均勻；當水流經過第一排植被后，床面切應力開始呈現波浪形狀分布；隨后在植被群內，床面切應力沿程先增大后減小；直至水流經過所有植被后，床面切應力逐漸減小。比較圖4c)的橫斷面x/D=93和x/D=175床面切應力沿程變化規律顯示，植被的存在減小了床面大部分的切應力，具有促淤保灘的效果。值得一提的是，在本研究中，橫斷面靠邊壁處的床面切應力經過植被群之后是增加的，這是由于計算區域的邊界受到限制而導致的。

圖4 床面切應力分布對比

2.4 床面形態變化

前文已對植被對水流結構影響進行了定量分析，現對床面泥沙的影響進行直觀展現。圖5展示了各工況床面形態變化，較為直接地呈現了植被水流對床面泥沙的作用與影響。由圖5可見，在單根植被周圍形成了沖刷孔，在單根植被下游形成沙脊。這是由于水流在遇到植被時形成圓柱繞流，其生成的馬蹄渦向下淘刷形成沖刷孔；單根植被下游近區存在滯水區，滯水區流速較低，泥沙容易在此沉積，所以形成了沙脊。同時，在圖5中發現，從工況1至工況3，隨著植被密度的增加，植被區域形成沙脊逐漸明顯；而從工況3到工況4時，植被密度的增加使整體形成沙脊的現象不那么明顯，而是以單根植被下游近區的局域沙堆出現。這說明了渦街被較高密度的植被耗散了，以至于渦街并未達到下一排植被，只能在單根植被下游近區局部形成沙堆，從而減少了植被群下游床面泥沙的沖刷。

圖5 各工況沖刷情況

3 結論

1)與試驗結果對比表明，過流流速與試驗結果較為吻合，該數值模擬方法可以有效地進行水流流場復現。

2)隨著植被密度的增加，流速呈現先增大后減小的規律，其中原因是由于卡門渦街在較高密度的植被中容易耗散，以至于渦街并未達到下一排植被。

3)床面切應力沿程變化規律表明，植被的存在減小了床面大部分的切應力，具有促淤保灘的效果；隨著植被密度的增加，床面切應力呈現先增大后減小的規律。

4)床面形態變化規律顯示，在單根植被周圍形成了沖刷孔，在植被下游形成沙脊；當植被密度增加到一定程度時，整體形成沙脊的現象不那么明顯，而是以單根植被下游近區的局域沙堆出現。