壩岸根石沖刷坍塌數(shù)值模擬

，，，

(南昌大學(xué) 建筑工程學(xué)院，南昌 330031)

1 研究背景

壩垛是為抵御水流直接沖刷堤身而修建的控導(dǎo)工程，是黃河防洪工程體系中的重要組成部分[1]。根石走失的方式主要是沖刷坍塌和沖揭走失2種[2]。沖刷坍塌即壩垛根基周圍河床由于長(zhǎng)期受到含沙水流的淘刷形成沖刷坑，隨著沖刷坑的不斷加大，壩岸根石自身穩(wěn)定遭到破壞，根石發(fā)生坍塌[3]。自1954年以來，黃河河務(wù)部門對(duì)根石走失開展探索試驗(yàn)與研究，分析了走失原因及規(guī)律，列舉了根石防護(hù)整治措施，為實(shí)際工程中根石防護(hù)積累了經(jīng)驗(yàn)[4]。田德本[5]從水流條件、根石斷面形態(tài)等方面對(duì)丁壩走失原因進(jìn)行了分析和探討，提出壩岸根石受旋渦條件影響最大。崔占峰等[6]引入VOF和空度法，模擬了丁壩附近水流，揭示了丁壩附近水流流動(dòng)的規(guī)律。張紅武等[7]通過局部動(dòng)床模型試驗(yàn)對(duì)不同的護(hù)根措施進(jìn)行模擬研究，提出了網(wǎng)罩護(hù)根措施。張志軍[8]從丁壩沖刷坑的形成及根石走失原因入手，分析了不同工況下各計(jì)算公式的適應(yīng)情況，得出根石走失規(guī)律。緱元有[9]通過對(duì)河道整治工程中根石受力情況，建立了根石走失的物理模型，得出根石的臨界起動(dòng)尺寸隨著含沙量和水流流速的增大而增大。孫桂環(huán)等[10]對(duì)壩垛沖刷坑的形成機(jī)理進(jìn)行了分析，對(duì)黃河下游河道整治工程的沖刷坑進(jìn)行了公式分析，并對(duì)實(shí)際河道段沖刷坑的計(jì)算參數(shù)進(jìn)行了選取驗(yàn)證。光治河等[11]根據(jù)山東黃河河道整治工程現(xiàn)狀提出，水流作用、根石粒徑和工程布局等都直接影響著壩岸中根石的走失。

由于黃河水流為高含沙水流，并且目前的根石探測(cè)技術(shù)局限性較大，對(duì)根石受沖揭影響的變形量及走失范圍多數(shù)仍停留在經(jīng)驗(yàn)性判斷[12-14]，根石走失尚無系統(tǒng)的理論分析方法，缺少可用于計(jì)算與分析水流特性與根石走失的數(shù)值模型。本文基于Flow-3D軟件建立根石沖刷坍塌走失特性研究的數(shù)值模型，計(jì)算水流對(duì)黃河壩垛工程的河床沖刷，對(duì)沖刷完成后壩垛周圍不同區(qū)域的泥沙高程面進(jìn)行分析，對(duì)壩垛周圍四個(gè)位置的邊坡進(jìn)行穩(wěn)定性分析，得到實(shí)際工程中容易發(fā)生沖刷坍塌的部位，為河道整治工程設(shè)計(jì)、防汛搶險(xiǎn)預(yù)案的制定提供理論依據(jù)與技術(shù)支撐。

2 數(shù)學(xué)模型與基本理論

2.1 水流控制方程

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

(3)

(4)

式中：u，v，w為x，y，z3個(gè)方向上的流速分量；Gx，Gy，Gz為x，y，z3個(gè)方向上的重力加速度；fx，fy，fz為x，y，z3個(gè)方向上的黏滯力加速度；Ax，Ay，Az為x，y，z3個(gè)方向上流體面積分?jǐn)?shù)；VF為流體的體積分?jǐn)?shù)；ρ為流體密度；p為作用在流體微元上的壓強(qiáng)。

2.2 泥沙控制方程

懸移質(zhì)方程：

(5)

基于梅耶-彼得的推移質(zhì)方程：

(6)

2.3 RNG k-ε模型

RNGk-ε模型控制方程和標(biāo)準(zhǔn)的k-ε模型相似，該模型在模擬計(jì)算低強(qiáng)度紊流和具有強(qiáng)剪切力區(qū)域時(shí)能力較好，因此適用于泥沙沖刷模型?？刂品匠贪ㄎ蓜?dòng)能k方程和紊動(dòng)能耗散率ε方程。

紊動(dòng)能k方程為

(7)

紊動(dòng)能耗散率ε方程為

(8)

2.4 FAVOR技術(shù)

Flow-3D采用結(jié)構(gòu)化矩形網(wǎng)格的有限差分法[15-16]，在此基礎(chǔ)上應(yīng)用FAVOR技術(shù)，通過計(jì)算該區(qū)域內(nèi)物體的面積和體積分?jǐn)?shù)定義復(fù)雜幾何實(shí)體模型。將該技術(shù)應(yīng)用于泥沙模型，并通過FAVOR技術(shù)在每個(gè)計(jì)算時(shí)間步里計(jì)算出網(wǎng)格內(nèi)泥沙的面積和體積分?jǐn)?shù)，對(duì)沖刷沉積過程中泥沙曲面的變化描述是準(zhǔn)確的，解決了在湍流剪切應(yīng)力作用下，泥沙面變化復(fù)雜難以準(zhǔn)確模擬的難題。

3 泥沙模型的數(shù)學(xué)驗(yàn)證

3.1 槽內(nèi)潰壩水流沖沙試驗(yàn)

泥沙模型采用Fraccarollo等[17]進(jìn)行的槽內(nèi)潰壩水流沖沙的物理模型試驗(yàn)進(jìn)行驗(yàn)證，模型布置見圖1(a)。試驗(yàn)在水槽內(nèi)進(jìn)行，水槽長(zhǎng)為2.5 m，寬為0.1 m，試驗(yàn)開始前在水槽內(nèi)平鋪滿0.06 m厚的細(xì)沙，泥沙面上以上設(shè)置高0.1 m、長(zhǎng)1.0 m的水體并拿閘門攔擋，試驗(yàn)?zāi)M閘門快速提起后形成的潰壩水流對(duì)泥沙面的沖刷情況。按1∶1比例建立數(shù)學(xué)模型，如圖1(b)。采用0.005 m×0.005 m×0.005 m的網(wǎng)格對(duì)計(jì)算區(qū)域劃分見圖1(c)，計(jì)算總網(wǎng)格數(shù)為40萬個(gè)。

圖1模型布置及網(wǎng)格劃分
Fig.1Modelandgriddivision

邊界條件設(shè)定：側(cè)面與底面均設(shè)為固壁邊界，計(jì)算區(qū)域上方設(shè)一個(gè)標(biāo)準(zhǔn)大氣壓。泥沙設(shè)定：泥沙參數(shù)設(shè)定參考原文獻(xiàn)中的泥沙屬性，取泥沙中值粒徑為3.5 mm，密度為1 540 kg/m3，泥沙休止角取30°。重力加速度為9.81 m/s2。

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

閘門的瞬間提起引起水體形成潰壩水流在槽內(nèi)傳播，造成閘門下游泥沙顆粒的起動(dòng)，在閘門下游處形成沖刷坑。本文由此分析深泓線處的水面線和泥沙面變化。由圖2中物理模型試驗(yàn)與數(shù)值模擬計(jì)算結(jié)果中水面線變化中可知，物理模型試驗(yàn)實(shí)測(cè)的水面線水舌處水位起伏較大，而數(shù)值模擬計(jì)算的水面線較為平滑，但兩者在數(shù)值上基本接近。結(jié)合0.25，0.50，0.75，1.00 s時(shí)的水面線可以看出物理模型試驗(yàn)結(jié)果中水面線均略高于數(shù)值模擬的計(jì)算值，考慮到物理試驗(yàn)中影響因素較多，卷氣現(xiàn)象和泥沙對(duì)水流影響較大，認(rèn)為存在該誤差合理。結(jié)合圖2中泥沙線可以看出，數(shù)值模擬結(jié)果中泥沙面計(jì)算值和物理模型試驗(yàn)實(shí)測(cè)值接近，說明軟件中的泥沙模型計(jì)算水流作用下泥沙床面變化準(zhǔn)確度高。

圖2不同時(shí)刻的泥沙線及水面線
Fig.2Profilesofsedimentandwaterlevelatdifferentinstances

4 根石沖刷坍塌數(shù)值模擬研究分析

4.1 根石沖刷坍塌模型建立

利用泥沙模型計(jì)算過壩洪水對(duì)河床的沖刷，對(duì)沖刷完成后壩垛周圍不同區(qū)域的泥沙高程面進(jìn)行分析。建立壩垛模型如圖3(a)和圖3(b)，將所建的實(shí)體模型按照1∶1比例導(dǎo)入軟件中，計(jì)算區(qū)域長(zhǎng)190 m、寬35 m、高2 m，根石護(hù)坡埋入泥沙內(nèi)部1 m，護(hù)坡的坡度為1∶1，在出口處設(shè)尾墩控制流量。采用結(jié)構(gòu)化矩形網(wǎng)格劃分如圖3(c)，網(wǎng)格總數(shù)約為110萬個(gè)。試驗(yàn)?zāi)M洪水對(duì)壩垛周圍泥沙床面的淘刷，A-A，B-B，C-C，D-D斷面分別為離入水口處60，70，80，95 m的監(jiān)測(cè)斷面，如圖3(d)。

邊界條件設(shè)定：入口邊界采用水位上升邊界，設(shè)定入口邊界為在300 s以內(nèi)水深由0.5 m上升到1.84 m隨后保持在1.84 m；出口邊界設(shè)定為自由出流邊界；壁面及底面均為固壁邊界；計(jì)算區(qū)域上方為一個(gè)標(biāo)準(zhǔn)大氣壓。計(jì)算時(shí)間為900 s。泥沙設(shè)置：泥沙的中值粒徑為0.2 mm，密度為2 630 kg/m3，推移質(zhì)系數(shù)取8，攜沙系數(shù)取0.018，泥沙休止角為32°。

圖3壩垛模型計(jì)算網(wǎng)格及泥沙斷面監(jiān)測(cè)布置
Fig.3Computationmodelofdambuttressandarrangementofsedimentmonitoringsections

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

4.2.1 流速場(chǎng)分析

圖4為不同時(shí)刻水流過壩垛后的水流矢量圖。

圖4不同時(shí)刻水流矢量圖
Fig.4Flowvectorsatdifferentinstances

水流過壩后由于河床的束窄作用壩前流速變大，并且在壩襠間形成漩渦。由圖4中可以看出隨著上游來水水位的上漲，壩垛將主流挑向壩對(duì)岸，流速較大的區(qū)域?yàn)閴吻昂蛪螌?duì)岸處，此處水流的沖刷能力較大，而壩襠間的渦流形態(tài)變化復(fù)雜，泥沙在此區(qū)域可能沉積也可能被帶動(dòng)。

4.2.2 流速分析

圖5為不同時(shí)刻水流平均流速的等值線圖。

圖5不同時(shí)刻水流平均流速等值線
Fig.5Contoursofflowvelocityatdifferentinstances

由圖5可知：t=100 s時(shí)水流過壩，由于壩垛的束窄作用，水流流速在壩前頭區(qū)域附近最大，迎水面與壩襠間處水流流速小；t=200 s時(shí)入口邊界水深加高，水流在壩前受擠壓作用并被挑向壩垛對(duì)岸后一分為二，匯入壩襠間的水流在該處形成漩渦，靠近第2個(gè)壩段迎水面根部處水流流速較大和漩渦強(qiáng)度大；t=300 s時(shí)入流邊界已經(jīng)達(dá)到穩(wěn)定，水位保持在洪水位1.84 m，水流已漫過整個(gè)壩垛段的根石護(hù)坡，水流在壩前頭處流速達(dá)到最大，主流被壩垛挑向?qū)Π叮髟诘?個(gè)壩垛迎水面根部處流速較大；t=600 s和t=900 s時(shí)，由于上游來水已經(jīng)穩(wěn)定，因此水流流態(tài)變化不大，整體區(qū)域的水流流速變化小，壩襠間存在較大的水流漩渦。由系列等值線圖可以看出整個(gè)過程中迎水面處水流流速小且變化不大，壩前頭處水流流速最大，上跨角次之，下跨角處水流流速較小，壩襠間主要是漩渦水流的運(yùn)動(dòng)。

4.2.3 泥沙分析

圖6是壩垛河床的泥沙輸運(yùn)圖，結(jié)合圖5可知：t=100 s時(shí)，水流在壩前頭處流速較大，水流帶動(dòng)該處泥沙形成懸移質(zhì)運(yùn)輸，而其他區(qū)域基本無泥沙起動(dòng)；t=200 s時(shí)，壩前頭處被帶動(dòng)的泥沙在壩襠間運(yùn)動(dòng)，結(jié)合圖4可知此時(shí)壩襠間有水流漩渦的存在；t=300 s時(shí)，水流中懸沙濃度較大的區(qū)域?yàn)閴吻邦^主流區(qū)，結(jié)合圖5可知此時(shí)壩前頭處水流流速達(dá)到最大，過壩水流帶動(dòng)壩前頭大部分泥沙形成高含沙水流，被攜帶的泥沙主要隨壩垛對(duì)岸主流運(yùn)動(dòng)，少部分泥沙進(jìn)入壩襠間隨水流漩渦運(yùn)動(dòng)；t=600 s和t=900 s時(shí)水流已經(jīng)穩(wěn)定，水中的含沙量變化不大，由于沖刷坑的形成加大了壩前頭處水深，此時(shí)水流的挾沙能力保持穩(wěn)定減小的趨勢(shì)。由圖6可知：在壩垛沖刷的整個(gè)過程中壩前頭處水流挾沙能力最大，當(dāng)水位達(dá)到設(shè)計(jì)洪水位時(shí)水流的挾沙能力最強(qiáng)，隨著壩前沖坑的形成與上游來水水位的穩(wěn)定，水流的挾沙能力逐漸減小并繼續(xù)淘刷周圍河床，壩前頭處泥沙大部分被主流帶往下游，少部分泥沙進(jìn)入壩襠間隨漩渦水流運(yùn)動(dòng)。

圖6壩垛周圍河床泥沙輸運(yùn)圖
Fig.6Sedimenttransportalongtheriverbed

4.2.4 泥沙床面變化分析

圖7是泥沙床面高程等值線，圖8是各監(jiān)測(cè)斷面的泥沙高程線。

圖7 泥沙沖刷高程線變化Fig.7 Variation of sediment line elevation

圖8各監(jiān)測(cè)斷面泥沙高程線
Fig.8Elevationsofsedimentatmonitoringsections

由圖7和圖8可知：t=100 s時(shí)僅壩前頭處有微小沖刷，泥沙床面基本保持平整；t=200 s時(shí)壩前頭處泥沙沖刷范圍變大，結(jié)合圖4可知壩襠間的漩渦水流對(duì)連壩段處造成輕微的淘刷；t=300 s時(shí)泥沙床面有輕微變化，壩前頭出現(xiàn)沖刷坑，結(jié)合圖8中C-C監(jiān)測(cè)斷面可知壩前沖坑深達(dá)到-0.25 m，B-B監(jiān)測(cè)斷面顯示上跨角處也有沖坑的存在，從圖7可知上跨角和壩垛對(duì)岸都出現(xiàn)了相應(yīng)的沖坑，在下跨角前部出現(xiàn)了泥沙淤積現(xiàn)象，D-D斷面出現(xiàn)了泥沙的淤積現(xiàn)象，且泥沙淤積高度最大，約為0.07 m；t=600 s時(shí)由B-B監(jiān)測(cè)斷面可知上跨角處泥沙沖坑深度加深，達(dá)-0.25 m，相比于壩前頭沖刷坑的縱向發(fā)展速度更快，此時(shí)水流在上跨角處的泥沙淘刷能力強(qiáng)于壩前頭處，D-D斷面顯示由于水流作用下跨角處的泥沙沖坑逐漸加深，但沖坑的發(fā)展速度緩慢；t=900 s時(shí)壩前頭處的沖坑深度達(dá)到約-0.55 m，B-B斷面顯示水流在600 s以后并沒有對(duì)沖坑形成縱向的淘刷，而是橫向的淘刷，此時(shí)沖坑的深度為-0.3m，D-D斷面顯示600 s以后由于沖刷范圍的變大，水流開始淘刷下跨角，下跨角處的沖坑快速發(fā)展，深度約為-0.3 m。在整個(gè)計(jì)算中A-A斷面泥沙床面變化程度小，說明該處河床受水流作用影響小。

5 結(jié) 論

本文基于Flow-3D軟件建立了根石沖刷坍塌走失的數(shù)值模型，對(duì)實(shí)際壩垛工程進(jìn)行計(jì)算得到結(jié)論如下：

(1)水流經(jīng)過壩垛時(shí)流速沿壩身逐漸增大，在壩前頭處達(dá)到最大，上跨角次之，迎水面和下跨角處水流流速較小，壩襠間區(qū)域存在較大的水流漩渦。

(2)壩前頭處水流挾沙能力最大，該處沖坑出現(xiàn)的時(shí)間最早，最終沖坑深最大，該處沖坑前期主要以縱向淘刷為主，隨后逐漸向上跨角與下跨角區(qū)域橫向發(fā)展，此兩處水流沖刷程度相當(dāng)，水流淘刷上跨角的時(shí)間早于下跨角，此幾處被挾帶的泥沙在壩襠間由于水流漩渦作用部分沉積形成堆丘。

(3)對(duì)各斷面的泥沙沖坑進(jìn)行了對(duì)比分析，沖坑深度由大到小依次為壩前頭、上跨角、下跨角和迎水面，推斷出實(shí)際工程中根石護(hù)坡在水流沖刷作用下最易出現(xiàn)沖刷坍塌的部位為壩前頭處，其次為上跨角與下跨角，上跨角出險(xiǎn)時(shí)間早于下跨角，迎水面處泥沙沖刷程度很小，該處發(fā)生沖刷坍塌險(xiǎn)情的幾率較小。