動水壓力作用下堤防邊坡流固耦合滲流特性研究

吳 勇，黃江華，左思賢

(1. 浙江華東工程安全技術有限公司，杭州 310014;2. 鄭州大學 水利與環境學院，鄭州 450001)

堤防邊坡的地基土體多為第四紀松散沉積物，堤身填筑土的組成及密實程度都具有非均質性。水流或波浪沖刷作用下堤岸滲流、變形、穩定性研究對堤防的安全運行和管理具有十分重要的意義。

國內外眾多學者對波浪沖刷作用進行了大量的研究，主要集中于理論分析和數值計算。在理論分析方面，Sollitt 和Cross[1]利用特征函數法求解得到波浪反射系數和透射系數的解析解；Madsen[2]在假定水頭損失與速度為線性的前提下，給出了波浪在堤防土體內的衰減公式；夏艷軍[3]等利用方向譜區分不同方向波浪的特點，提出了一種現場實測防波堤波浪透射系數的方法。在數值計算方面，溫鴻杰[4]基于改進的SPH方法，模擬研究了規則波在出水堤上的破碎和越浪過程；關大瑋[5]基于CFD 軟件Flow-3D建立三維數值波浪水槽，模擬研究了規則波和不規則波作用下不可滲光滑斜坡堤的越浪過程；李東洋等[6]基于Open FOAM建立三維數值波浪水槽生成不規則波，研究了正向入射不規則波與扭王塊體護面斜坡堤的相互作用；張胡等[7]以線性波浪繞射理論為基礎，依據邊界元法，建立計算大尺度結構所受波浪力的三維數學模型。

在一定的滲透壓力作用下，水流或波浪會使堤防土體的土體顆粒產生動態變形，進而影響體積應變，當前關于水流在堤防飽和土體中的流固耦合滲流特性的研究報道較少。將堤防邊坡土體簡化為由固體土骨架、水組成的飽和介質，運用土動力理論，構建了動水壓力作用下堤防邊坡流固耦合滲流動力方程，分析了堤防邊坡中的應力、滲流等特性，對于堤防的防治具有一定的指導意義。

1 動水壓力作用下的流固耦合動態本構模型

堤防土體的流固耦合滲流規律是一個復雜的動力問題，涉及土力學、流體力學、動力學等諸學科，需要綜合考慮堤防土體的運動變形方程、考慮土骨架體積變形的河水滲流方程[8-9]。

1.1 堤防土體的運動變形方程

(1)有效應力原理。

堤防土體屬于多孔飽和介質，由土骨架及顆粒間孔隙共同組成。將堤防土體視為彈性介質，固體骨架的有效應力遵循修正的Terzaghi有效應力原理，即

(1)

(2)物理方程。

根據彈性多孔介質的相關理論可知，堤防土體變形的本構方程(即物理方程)

σij=2Gεij+λεvδij

(2)

(3)幾何方程。

與常規彈性材料相同，堤防土體固體骨架的幾何方程為

(3)

式中：u為堤防土體固體骨架的位移，m。

(4)動力平衡微分方程。

堤防土體固體骨架的位移是空間和時間的參數，即u=u(x,y,z,t)，參照飽和土體的平衡微分方程，給出含堤防土體的平衡微分方程

(4)

式中：w為河水相對于堤防土體固體骨架的位移，m；Fi為應力，Pa；ρ為堤防土體的總密度，kg/m3；ρ=φρs+(1-φ)ρw；ρs為固體土骨架的密度，kg/m3；φ為孔隙率；ρw為流水的密度，kg/m3。

將式(1)和(2)代入式(4)得

(λ+G)uj,ji+Gui,jj+Fi=αpi+ρü+ρw

(5)

根據uj,j=·u、uj,jj=2ui，式(5)可寫成矢量的形式

(λ+G)(·u)+G2u+F=αp+ρü+ρw

(6)

1.2 考慮土骨架體積變形的河水滲流方程

(1)孔隙率動態模型。

因此孔隙率是研究河水運移的重要物理參數之一?？紫堵师针S著堤防土體的運動變形而變化，即φ與其所處的應力狀態有關，對于等溫過程，其計算式如下

(7)

式中：φ0為某一基準溫度下且無荷載作用時堤防土體的初始孔隙率；εv為土骨架的體積應變；△p為河水壓力變化，Pa，△p=p-p0，p0為河水的初始滲透壓力。

(2)滲透率動態模型。

滲透率是反映河水滲流的重要參數指標，常規的河水滲流模型通常將滲透率視為常數，進而忽視了堤防土體骨架變形、孔隙河水壓力變化及孔隙體積變化等對滲透率的影響，忽略整個過程中的熱效應，根據滲流力學的Krozeny-Carman方程，滲透率k的計算式可簡化如下

(8)

式中：k0為某一基準溫度下且無荷載作用時堤防土體的初始滲透率，m2。

(3)河水含量計算模型。

將堤壩土體中的河水全部假定為自由水，即

Q=ρwφ

(9)

(4)Darcy滲流定律。

在壓力梯度△p作用下，河水在孔隙中作線性滲流運動，該過程遵從Darcy滲流定律，關系式為

(10)

式中：v為河水滲流速度矢量，m/s；η為河水的動力黏度，Pa·s。

(5)河水滲流模型。

煤層河水滲流遵從質量守恒定律，對于單位體積的土體，其滲流方程為

(11)

式中：Q為單位體積堤防土體中的河水含量，kg/m3。

(6)堤防土體總的體積應變。

考慮水壓力作用，堤防土體的體積應變εv由兩部分組成

(12)

① 固體骨架的變形引起的體積應變，對于二維平面應變問題，由式(3)可得

(13)

② 河水滲流壓力增量△p引起的堤防土體的體積應變

(14)

將式(3)和(7)～(10)代入式(11)，經過整理，可得孔隙水滲流方程

(15)

其中，S和Qm的表達式如下

(16)

(17)

(18)

圖1 二維計算模型Fig.1 2D computing model

2 數值解答

以某堤防為工程實例，建立二維計算模型，如圖1所示，堤防高10 m，河床厚10 m，河水深7 m，兩側邊坡1:1放坡，兩側各延伸20 m。

初始和邊界條件如下：模型的左右兩側邊界施加水平向約束，底邊界施加固支約束，其他側為自由邊界；孔隙水壓力p0=γwh(h為河水的壓力水頭)。

表1 物理力學參數Tab.1 Physical mechanics parameters

圖2 二維網格劃分模型Fig.2 Two-dimensional meshing model

假定堤防土體為均質的各向同性介質，堤防土體和河水的物理力學參數見表1。

運用COMSOL軟件，對流固耦合動力方程進行數據求解，二維網格劃分如圖2所示。

首先采用模型I，通過數值計算，得到了孔隙水壓力云圖、水平滲流速度云圖、豎向滲流速度云圖，如圖3所示。

3-a Mises應力云圖(MPa)3-b 孔隙水壓力云圖(Pa)

3-c 水平滲流速度云圖(mm/s)3-d 豎向滲流速度云圖(mm/s)圖3 模型I的計算結果圖Fig.3 Simulated results of model I

圖4 水平滲流速度等值線的對比圖(單位：mm/s)Fig.4 Comparison of horizontal seepage velocity contours

從圖3-a可以看出，河水的滲流能夠引起堤防邊坡土體較明顯的應力變化，河水表面與堤防土體的交界面處出現較大的應力集中，因此堤防邊坡在水土界面處添加防沖刷層是非常必要的保護措施；從圖3-b可以看出，河床土體的孔隙水壓力高于堤防邊坡，這與實際情況相符；從圖3-c和3-d可以看出，水平滲流和豎向滲流主要發生在堤防邊坡水土界面處的一定區域內，由于水頭壓力較小，在堤防邊坡遠處的區域內，滲流速度變緩并最終趨于零。

采用模型II，不考慮河水滲流壓力引起土顆粒位移而造成土骨架的體積變形，以水平滲流速度為例，對比了模型I的模型II的計算結果差異，如圖4所示。

圖5 滲流作用下的堤防沉降(單位：m)Fig.5 Embankment settlement under seepage

從圖4可以看出，相較于模型II，模型I(考慮流體壓力引起土骨架壓縮而引起的總體積應變的增大)的滲流影響區域較遠，在河水滲流動壓力的作用下，固體土骨架的壓縮變形增大，土的孔隙率提高，流水的滲流路徑較遠，因此在堤防邊坡防護設計時，應充分考慮水頭、水流速度等引起土體的動位移。

采用模型I，進一步得到了動水壓力作用下堤防邊坡的沉降云圖，如圖5所示。

由圖5所示，堤防邊坡的沉降主要發生在臨水邊坡一側，而滲流作用下的邊坡沉降較小，最大值僅為0.045 mm，可以忽略不計，結合文獻[1-7]的相關研究結果，水流動水壓力引起的沖刷變形是造成邊坡失穩的最主要外因。

3 結論

在外部壓力作用下，土骨架整體會發生壓縮變形，孔隙率降低；孔隙水滲流過程引起的拖曳力(即動水壓力)屬于土體內力，動水壓力差會造成土顆粒自身的壓縮，進而引起整體土骨架孔隙率的增大?？紤]動水壓力引起的土骨架總體積應變、堤防邊坡土體的孔隙率和滲透率受滲流影響的變量，運用飽和土力學和流體力學理論構建了飽和土體的運動方程和Darcy滲流方程，運用COMSOL軟件分析了動水壓力作用下堤防邊坡的流固耦合特性，研究結果表明：(1)模型Ⅰ(充分考慮孔隙水滲流壓力引起的體積變形)比模型Ⅱ(不考慮孔隙水滲流壓力引起的體積變形)更符合實際情況，動水壓力會引起固體土骨架的壓縮變形，土的孔隙率提高，流水的滲流路徑較遠，而這種趨勢隨著高動能水流沖刷力增大而提高；(2)堤防邊坡的沉降主要發生在臨水邊坡一側，滲流作用引起的邊坡沉降較小，最大值僅為0.045 mm，可以忽略不計。