一種模擬黏土心墻壩漫頂潰決的水沙耦合數(shù)學模型

王 廷，王亞林

(三峽大學 水利與環(huán)境學院，湖北 宜昌 443002)

1 研究背景

黏土心墻壩具有土石壩就地取材，對地形、地質適應能力強等特點，與均質壩和斜墻壩相比，其壩坡陡，相對工程量小，基礎防滲線短，抗震性能好，因此長期以來黏土心墻壩都是壩工界設計人員優(yōu)先選擇的一種壩型。在世界范圍內黏土心墻壩數(shù)量眾多，建設規(guī)模也不斷突破，目前在建和已建的超高黏土心墻壩有雙江口黏土心墻壩(314.0 m)，努列克黏土心墻壩(300.0 m)以及兩河口黏土心墻壩(295.0 m)等。這些工程在帶來巨大的經(jīng)濟效益與社會效益的同時也存在著潰壩風險。據(jù)統(tǒng)計1954—2014年我國的土質心墻壩潰壩數(shù)量達到183座，造成了大量的人員傷亡和財產(chǎn)損失，其中漫頂導致的潰壩超過50%[1]。因此，為了提高黏土心墻壩的安全控制水平，減輕黏土心墻壩漫頂潰決損失，需要對黏土心墻壩的漫頂潰決機理、潰口發(fā)展過程、潰口洪水流量過程、洪水演進過程以及評估潰壩風險[2]進行深入研究。

國內外學者通過黏土心墻壩漫頂潰決試驗發(fā)現(xiàn)[3-4]：隨著下游壩殼逐漸被沖蝕，心墻在上游動水壓力和土壓力的作用下，發(fā)生滑移或傾倒破壞(如圖1、圖2所示，其中Fa為壩殼作用在心墻上的土壓力，F(xiàn)w為庫水作用在心墻上的水壓力，F(xiàn)sb為沿破壞面底部作用的摩擦力，F(xiàn)cb為沿破壞面底部的黏結力)，導致潰口流量迅速增加，沖蝕加劇，直至大壩完全潰決。

圖1 水流淘刷下游壩殼示意圖Fig.1 Schematic diagram of scouring and breaking of downstream dam shell

圖2 黏土心墻滑動失穩(wěn)壩示意圖Fig.2 Schematic diagram of clay core sliding

根據(jù)這一潰決機理，任強等[5](2010)建立的數(shù)學模型能夠較好地描述心墻失穩(wěn)所引起的潰口間歇性增大現(xiàn)象。陳生水等[3](2011)提出了漫壩水流對下游壩殼料沖蝕的經(jīng)驗公式。鐘啟明等[6](2017)模擬了心墻傾倒破壞和剪切破壞。這些成果對于潰口發(fā)展的模擬研究均頗有意義，然而上述模型中水流計算采用的一維平均流速很難體現(xiàn)潰口水流復雜的水力特性，進而無法精確模擬下游壩殼的沖蝕過程。而下游壩殼沖蝕過程對于心墻失穩(wěn)破壞具有重要影響。

為了實現(xiàn)對潰口復雜沖蝕過程的模擬，本文建立了黏土心墻壩漫頂潰決的動床水沙耦合數(shù)學模型。其中水流運動方程采用固液兩相湍流方程，并通過PLIC-VOF法(分段線性界面重構-流體體積法)捕捉洪水自由液面的位置，精確預測洪水運動過程。對于壩殼非黏性沙輸移的計算則采用了不飽和全沙輸移方程。黏土心墻的破壞通過剪切力進行判斷。

2 數(shù)學模型

2.1 水-沙耦合的固液兩相湍流方程

水動力學模型采用了固液兩相湍流方程，可以較好地反映泥沙運移與床面變形對于水流運動的影響。

連續(xù)方程：

(1)

動量方程：

(2)

懸沙對流擴散方程：

流體湍流動能k方程：

(4)

其中：

μmt=ρmCμk2/ε。

式中：t為時間；ρm為混合物密度；um,i為混合物在xi方向上的流速；um,j為混合物在xj方向上的流速；αp為泥沙體積分數(shù)；δj3為重力方向Kronecker delta符號；ω為泥沙沉降速度；μm為混合流體黏性系數(shù)；p為動水壓強；νmt為湍渦黏度；σp為臨界正應力；Cμ和σk均為系數(shù)，取Cμ=0.09，σk=1.0；ε為動能耗散率。

流體湍流動能耗散率ε方程：

式中：C1、C2、σε均為系數(shù)，C1=1.44，C2=1.92，σε=1.3。

考慮含沙量變化的混合流體黏性系數(shù)為

(6)

式中：αs為泥沙體積分數(shù)；αcr為臨界底沙體積分數(shù)；αco為給定的常數(shù)；μf為清水有效黏性系數(shù)。

2.2 自由面處理

自由面重構采用PLIC-VOF法，F(xiàn)為流體體積函數(shù)，滿足輸運方程，即

(7)

PLIC-VOF法采用斜線段來模擬單個網(wǎng)格內的交界面。該方法主要涉及兩方面的問題[7]：

(1)已知網(wǎng)格體積分數(shù)F值條件下，計算自由表面的法向矢量n，從而精確重構自由表面。

(2)相鄰網(wǎng)格內的對流體積計算。

2.3 泥沙侵蝕

泥沙剪應力τb為

(8)

式中：τbo,i,j為編號網(wǎng)格剪應力值；τbo,i±1,j為i方向相鄰網(wǎng)格剪應力值；τbo,i,j±1為j方向相鄰網(wǎng)格剪應力值；|xn|、|zn|分別為底沙表面的單位法向向量在x、z方向的投影。

單個網(wǎng)格的剪應力為

(9)

其中：

(10)

μtot=μm+μmt。

(11)

泥沙減少量Cremoved為

Cremoved=αcrρsLtot。

(12)

式中：ρs為泥沙密度；Ltot為泥沙通量。Ltot可通過下式計算，即

(13)

式中：Δt為時間步長；τc為臨界剪應力。τc可表述為

τc=θcgi(ρs-ρf)d50。

(14)

式中：gi為重力加速度；d50為中值粒徑；ρf為流體密度；θc為校正斜率。

θc可以表述為

式中：η為休止角；ξ為河床坡度；θc,0為斜率初始值。

2.4 心墻滑移

流體計算過程中，黏土心墻作為固壁邊界，采用基于離散力的浸入邊界法來模擬。垂直邊界的速度滿足對流方程[8]為

(16)

式中u為垂直邊界速度。

完成每一步的流體計算，需進行心墻結構受力計算。心墻的破壞面位于壩殼垂向下切深度的水平面，計算破壞面上每個節(jié)點的法相應力σ和剪切應力τ，通過沿破壞面積分計算出整個破壞面上總的滑動力和抗滑力，心墻發(fā)生滑移破壞[9]時，作用在心墻上的力滿足式(17)，即

(17)

式中：c為黏土心墻的內聚力；φ為黏土心墻的內摩擦角；τxy為滑移面剪應力。

數(shù)值計算采用結構網(wǎng)格和控制體積法對計算區(qū)域進行離散，詳細計算過程可以參看文獻[10]。

3 數(shù)值模擬

數(shù)值計算的模型尺寸，材料參數(shù)與開展的室內驗證試驗一致，壩高0.6 m，頂寬0.1 m，上下游壩坡坡比1∶2，心墻頂寬0.06 m，心墻坡比1∶0.2。壩殼堆石體含水量2.7%，特征粒徑d50=6.5 mm，干重度20.7 kN/m3，無黏聚力，內摩擦角47.8°；心墻含礫黏土含水量6.1%，d50=0.9 mm，干重度19.5 kN/m3，黏聚力36 kPa，內摩擦角24.9°。初始水頭0.05 m，試驗水流流量0.01 m3/s。

數(shù)值模型由160 000個單元網(wǎng)格組成，其中單元尺寸為2 cm×2 cm×4 cm。前、后邊界及底面邊界均為壁邊界，頂面為對稱邊界。左、右邊界分別為流量邊界和出流邊界。

3.1 數(shù)值模擬潰口發(fā)展過程

整個潰壩過程可分為以下2個階段。

(1)逐漸沖蝕階段。在這一階段中，一方面水流結構，與下游壩殼的床面形態(tài)息息相關，而下游壩殼的床面形態(tài)又隨著水流結構的變化而不斷改變。另一方面隨著下游壩殼被沖蝕，心墻下游側臨空面逐漸加大，當心墻發(fā)生滑移破壞，將導致下泄流量突增。如圖3所示，由于局部擾動引起近壁層流層波動，下游壩坡形成沙紋，沙紋尺度較小，此時壩面起伏幅度不大(圖3(a)、圖3(b))。隨著水流沖蝕，壩面形成與水面波同相位的沙浪形態(tài)，這時的壩面沙浪與水面波浪之間有強烈的相互干涉作用(圖3(c)—圖3(g))。此后，水流在迎水面產(chǎn)生局部水躍，背水面則產(chǎn)生水跌，致使水面波動超過壩面的起伏，形成“陡坎”(圖3(h)—圖3(l))。“陡坎”逐漸發(fā)展，使得壩頂下緣位置心墻臨空面逐漸加大，當土壓力和水壓力的共同作用超過心墻的極限抗剪強度，心墻發(fā)生第一次滑移破壞。

圖3 潰壩過程逐漸沖蝕階段Fig.3 Gradual erosion stage

(2)劇烈潰決階段。如圖4所示，心墻發(fā)生第一次滑移破壞之后，泄流流量突增，漩渦水流向內劇烈淘刷，引起第二次心墻滑移破壞，流量進一步增大，在達到峰值之后又快速減小，直至潰決結束。

圖4 潰壩過程劇烈沖蝕階段Fig.4 Violent outburst stage

3.2 數(shù)值模擬與室內試驗潰口流量結果對比

為了驗證模擬的精度，開展了室內漫頂潰壩試驗。對于正態(tài)模型，非黏性土在滿足砂礫雷諾數(shù)Re*>70的前提下，符合式(18)、式(19)，即可滿足模型的水流流態(tài)相似[11]。

λd=λh=λL，

(18)

λρs=λρ=1 。

(19)

式中：λL為幾何比尺；λd為砂粒粒徑比尺；λρs為砂粒密度比尺；λh為水深比尺；λρ為水密度比尺。

試驗系統(tǒng)工作原理如圖5所示，試驗系統(tǒng)平面布置示意圖如圖6所示。試驗中使用高清攝像機對潰口發(fā)展過程進行記錄。粒子圖像測速系統(tǒng)采集平均流速，并結合斷面形態(tài)可以得到潰壩流量。全自動水位激光觀測儀，進行水位變化的監(jiān)測。初始水頭為0.05 m，試驗水流流量為0.01 m3/s。

圖5 試驗系統(tǒng)工作原理Fig.5 Working principle diagram of test system

圖6 試驗系統(tǒng)平面布置示意圖Fig.6 Plane layout diagram of test system

根據(jù)圖7可知，計算與試驗中的流量過程線較為接近。第一次心墻滑移導致流量迅速增大的時間點計算值比試驗提前42 s，模擬的心墻滑移面較試驗高0.03 m。第二次心墻滑移導致流量迅速增大的時間點計算值比試驗滯后27 s，模擬的心墻滑移面較試驗值低0.05 m。計算值與試驗值的潰口峰值流量分別為0.142、0.149 m3/s，數(shù)值模擬中峰值流量出現(xiàn)的時間點較試驗值推后58 s。

圖7 潰口流量過程對比曲線Fig.7 Computed and test discharge curves at the breach

4 結 論

(1)黏土心墻壩在漫頂潰決過程中，潰壩水流沖刷壩殼，壩殼發(fā)生變形又反過來影響水流波動。下游壩殼的變化過程為沙紋、沙浪、陡坎3種形態(tài)，直到黏土心墻發(fā)生失穩(wěn)滑移破壞，導致潰壩水流突增，大壩進入劇烈潰決階段。

(2)基于黏土心墻壩的漫頂潰決機理，建立的黏土心墻壩潰壩水沙耦合模型，能夠在一定程度上模擬整個黏土心墻壩漫頂潰決過程，經(jīng)室內水槽試驗檢驗，預測的潰口流量過程與實測值也較為吻合，為研究黏土心墻壩漫頂潰口發(fā)展過程提供一種可選工具。