剛性植被密度對(duì)岸灘作用的數(shù)值模擬*

肖苡辀，張國(guó)志，馮先導(dǎo)，陳文峰

(1.中交第二航務(wù)工程局有限公司，湖北 武漢 430040；2. 長(zhǎng)大橋梁建設(shè)施工技術(shù)交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430040；3. 交通運(yùn)輸行業(yè)交通基礎(chǔ)設(shè)施智能制造技術(shù)研發(fā)中心，湖北 武漢 430040；4. 中交公路長(zhǎng)大橋建設(shè)國(guó)家工程研究中心有限公司，湖北 武漢 430040)

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

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

(1)

(2)

式中：ρ為流體密度；ui為局部時(shí)均速度分量；u為速度矢量；t為時(shí)間；μ為流體動(dòng)力黏滯系數(shù)；p為壓強(qiáng)；Fi為微元體上的體力。

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

(3)

(4)

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

局部的Shields參數(shù)θi是基于局部床面切應(yīng)力τ計(jì)算的：

(5)

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

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

(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相的體積分?jǐn)?shù)；Φi為無量綱推移質(zhì)輸沙率，假定為0.05。容積床載輸送率qb，i表示為：

(9)

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

(10)

式中：up是在某點(diǎn)j處的主流流速；k為von Karman常數(shù)，為0.418；u為與底部剪切應(yīng)力有關(guān)的剪切速度；cμ為常數(shù)；kj為點(diǎn)j處的紊動(dòng)能；ΔB為粗糙度函數(shù)；zj為點(diǎn)j的垂向高度。

1.2 流場(chǎng)設(shè)置與邊界條件

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

圖1 計(jì)算模型

表1 試驗(yàn)工況設(shè)計(jì)

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

1.3 計(jì)算方法

基于有限差分法采用RNGk-ε紊流模型來求解控制方程，對(duì)流項(xiàng)采用二階迎風(fēng)格式，擴(kuò)散項(xiàng)采用二階中心差分格式。采用TruVOF方法縮減收斂時(shí)間，該方法僅計(jì)算流體單元，不計(jì)算空氣單元，且對(duì)自由液面的描述更為準(zhǔn)確。通過FAVOR法(fractional area volume obstacle representation)對(duì)計(jì)算區(qū)域劃分網(wǎng)格。當(dāng)控制方程的殘差小于1×10-5時(shí)，認(rèn)為計(jì)算結(jié)果收斂。

2 結(jié)果分析

2.1 試驗(yàn)驗(yàn)證

為了更好地研究植被存在時(shí)水動(dòng)力結(jié)構(gòu)和懸沙運(yùn)動(dòng)特性，需要對(duì)數(shù)值模擬的可靠性進(jìn)行驗(yàn)證?；谠恋脑囼?yàn)結(jié)果進(jìn)行數(shù)值模擬，從而對(duì)比分析數(shù)值模擬的可靠性。本研究基于試驗(yàn)設(shè)置和結(jié)果進(jìn)行數(shù)值模擬，結(jié)果見圖2。圖2為垂向流速分布模擬值與試驗(yàn)值對(duì)比，結(jié)果表明模擬值與試驗(yàn)值吻合度較好，該數(shù)值模擬方法能夠有效模擬植被水流工況。

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

2.2 流速分布

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

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

2.3 床面切應(yīng)力

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

(11)

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

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

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

圖4 床面切應(yīng)力分布對(duì)比

2.4 床面形態(tài)變化

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

圖5 各工況沖刷情況

3 結(jié)論

1)與試驗(yàn)結(jié)果對(duì)比表明，過流流速與試驗(yàn)結(jié)果較為吻合，該數(shù)值模擬方法可以有效地進(jìn)行水流流場(chǎng)復(fù)現(xiàn)。

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

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

4)床面形態(tài)變化規(guī)律顯示，在單根植被周圍形成了沖刷孔，在植被下游形成沙脊；當(dāng)植被密度增加到一定程度時(shí)，整體形成沙脊的現(xiàn)象不那么明顯，而是以單根植被下游近區(qū)的局域沙堆出現(xiàn)。