吳 勇，黃江華，左思賢

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

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

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

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

1 動水壓力作用下的流固耦合動態(tài)本構(gòu)模型

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

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

(1)有效應(yīng)力原理。

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

(1)

(2)物理方程。

根據(jù)彈性多孔介質(zhì)的相關(guān)理論可知，堤防土體變形的本構(gòu)方程(即物理方程)

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

(2)

(3)幾何方程。

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

(3)

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

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

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

(4)

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

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

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

(5)

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

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

(6)

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

(1)孔隙率動態(tài)模型。

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

(7)

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

(2)滲透率動態(tài)模型。

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

(8)

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

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

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

Q=ρwφ

(9)

(4)Darcy滲流定律。

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

(10)

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

(5)河水滲流模型。

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

(11)

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

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

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

(12)

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

(13)

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

(14)

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

(15)

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

(16)

(17)

(18)

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

2 數(shù)值解答

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

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

表1 物理力學參數(shù)Tab.1 Physical mechanics parameters

圖2 二維網(wǎng)格劃分模型Fig.2 Two-dimensional meshing model

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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