董祎挈，陸海軍，戴　睿，何　麗

(1.武漢輕工大學(xué) 多孔介質(zhì)力學(xué)研究所，湖北 武漢 430023；2.江漢大學(xué)文理學(xué)院 機(jī)電與建筑工程學(xué)部，湖北 武漢 430056)

?

水流作用下填埋場開裂黏土襯墊數(shù)值仿真計算

董祎挈1，陸海軍1，戴睿1，何麗2

(1.武漢輕工大學(xué) 多孔介質(zhì)力學(xué)研究所，湖北 武漢 430023；2.江漢大學(xué)文理學(xué)院 機(jī)電與建筑工程學(xué)部，湖北 武漢 430056)

摘要：針對垃圾填埋場壓實(shí)黏土封場系統(tǒng)在干旱與半干旱地區(qū)出現(xiàn)的干燥開裂的問題，基于流體動力學(xué)流動方程，建立了多孔介質(zhì)不連續(xù)裂縫中水流流動數(shù)學(xué)模型，采用COMSOL模擬開裂黏土水滲流規(guī)律。數(shù)值模擬計算結(jié)果表明，流體的壓強(qiáng)所出現(xiàn)的最大值和最小值分別出現(xiàn)在入口和出口的邊緣處，流體在裂縫中流速較大的區(qū)域主要分布在水流的入口和出口處的邊緣，越靠近裂縫的兩端，裂縫中的流體的運(yùn)動狀態(tài)與基質(zhì)中的流動狀態(tài)差異性越大。

關(guān)鍵詞：數(shù)值仿真；壓實(shí)黏土；干燥開裂

1引言

壓實(shí)黏土封場系統(tǒng)即填埋場頂部隔斷層，由1m厚壓實(shí)黏土層、1mm厚的HDPE防滲膜和400 g/m2的無紡布復(fù)合構(gòu)成，用來阻擋地表降水滲入廢物層，以減少滲瀝液的滲出量。在自然降雨入滲條件下封場系統(tǒng)的開裂是威脅填埋場運(yùn)行安全的重要因素之一。封場系統(tǒng)內(nèi)的水體流動易致使壓實(shí)黏土層滲透系數(shù)、水導(dǎo)熱系數(shù)等物性參數(shù)變化，進(jìn)而改變封場系統(tǒng)的水分遷移和土體變形等特性，最終影響封場系統(tǒng)的開裂[1-3]。

水流流動是造成壓實(shí)黏土開裂的原因之一，水流在壓實(shí)黏土裂縫和孔隙中流動過程中所引起黏土顆粒表面摩擦系數(shù)、顆粒所受壓力等物性參數(shù)的變化對黏土的抗壓、抗剪強(qiáng)度等工程特性造成一定的影響，進(jìn)而影響壓實(shí)黏土中的應(yīng)力場與位移的變化，最終導(dǎo)致壓實(shí)黏土封場系統(tǒng)的開裂[4-6]。由于流體在多孔介質(zhì)中流動會受到孔隙壁摩擦力的影響，流體的流速會變得非常小，故一般會采用達(dá)西滲流模型來模擬流體通過多孔介質(zhì)間隙的流動[7]。這種模型也被廣泛應(yīng)用于模擬水流在蓄水層和河岸里的遷移過程、油井中石油的遷移過程以及地底巖漿向火山口的遷移過程等。

國內(nèi)外科學(xué)工作者[8-9]大多是通過干濕循環(huán)試驗對填埋場黏土開裂情況進(jìn)行宏觀觀測，或采用離散元模擬了土壤干燥收縮及開裂特性[10]，而水流對開裂黏土裂縫的數(shù)學(xué)模型研究相對較少。

為了準(zhǔn)確預(yù)測水流作用下的填埋場壓實(shí)黏土封場系統(tǒng)的物性參數(shù)的變化過程，通過COMSOL開展了流體在三條裂縫中流動的仿真研究，建立多孔介質(zhì)不連續(xù)裂縫中水滲流數(shù)學(xué)模型，模擬研究壓實(shí)黏土水流作用下的物性參數(shù)變化對填埋場壓實(shí)黏土封場系統(tǒng)開裂，為預(yù)測水流作用下壓實(shí)黏土封場系統(tǒng)物性參數(shù)的變化過程提供了依據(jù)。

2開裂黏土水滲流數(shù)學(xué)模型

2.1　流體流速方程

達(dá)西定律的應(yīng)用條件為流體在多孔介質(zhì)中流動的驅(qū)動力為水力勢能梯度，通過考慮壓力和位置勢能的不同，而從流線起點(diǎn)到終點(diǎn)將水力勢能場形象化。根據(jù)達(dá)西定律，通過單位孔隙表面的凈通量[11]為：

(1)

式中：u為達(dá)西流速(m/s)；κ為多孔介質(zhì)的滲透系數(shù)(m2)；η為流體的動力黏度(Pa·s)；p為流體的壓力(Pa)；ρf為體密度(kg/m3)；g為重力加速度(m/s2)；▽D為g方向上的單位矢量。其中滲透系數(shù)κ表示單位水力梯度下的單位流量，通常表征流體通過孔隙骨架的難易程度。

方程中的水力勢能，來自壓力p和重度ρfgD。定義g=9.82 m/s2，D=0。D的選擇對模擬的結(jié)果和其中包括的物理參數(shù)有重要影響。例如，如果D是垂直z軸的方向，同時流體流向是完全與xy平面平行的，隨著D梯度的變化，流體的驅(qū)動力就只有單獨(dú)的壓力梯度。

將達(dá)西定律代入連續(xù)性方程后得到以下廣義控制方程：

(2)

式中：θs為流體的體積部分；Qs為流體的流量。

對不可壓縮流體來說，ρf可以從等式兩端約去，因此控制式2可化為：

(3)

其中S為相關(guān)系數(shù)，可用來表示包括在模型中設(shè)定固體變形方程或是從其他對溫度和濃度的分析結(jié)果。COMSOL在其表達(dá)式區(qū)域內(nèi)將S定義為利用流體或固體壓縮系數(shù)求得的指定的系數(shù)。應(yīng)用型達(dá)西定律模型采用了式(3)，并同時要求流體是不可壓縮的(即流體密度ρf是一個定值)。

2.2　土體基質(zhì)中水滲流方程

以時間為變量的土體基質(zhì)中的水流流動方程可以通過達(dá)西定律求得：

(4)

式中：p為流體壓強(qiáng)(Pa)；θs為孔隙率；χf、χs分別為流體和固體的壓縮系數(shù)(m·s2/kg)；κm為土體滲透系數(shù)(m2)；η為黏滯系數(shù)(kg/m·s)。

基質(zhì)模型中已經(jīng)設(shè)定好的流速變量為uesdl(m/s)，給出了達(dá)西流速表達(dá)式：

(5)

基質(zhì)塊的所有面上的流量平衡方程為：

n·uesdl=0.

(6)

式中：n為垂直邊界向外的單位方向向量。

2.3　基質(zhì)裂縫中水流流動方程

模型中的裂縫是一系列的內(nèi)部邊界。傳統(tǒng)的邊界模型定義流體的流動是垂直于邊界的而不是沿著邊界的。而在這個模型中，可通過從COMSOL軟件中特定的切向選項來定義流體是可以沿著內(nèi)部邊界和裂縫中流動的。

因此，為了使定義的這種模型可以有效地實(shí)行，必須建立裂縫邊界方程來求解和土體基質(zhì)方程所求得的相同的變量p(壓力)。在此模型中，裂縫中流體流速方程是基質(zhì)流動流速方程(即達(dá)西定律)的變形?？梢酝ㄟ^改變方程中的相關(guān)系數(shù)來改變流體在裂縫中流動的阻力系數(shù)以此適應(yīng)達(dá)西滲流方程。下式給出了裂縫和土體基質(zhì)之間的連續(xù)性方程：

(7)

式中：Sf為裂縫存儲系數(shù)(m·s2/kg)；κf為裂縫的滲透系數(shù)(m2)；dfrac為裂縫的厚度(m)；

由于流動方程中出現(xiàn)了裂縫的厚度，故單位長度上的流量為：

(8)

式中：▽T為裂縫切平面上的梯度算子；裂縫中的線流速為ulin=uesdl/dfrac(m/s)

裂縫的邊界條件及假設(shè)：(1)出口和入口的壓力是已知的；(2)入口的壓力是不變的常數(shù)；(3)出口的壓力是隨著時間而線性變化的；(4)其他邊界沒有流體流入或流出。方程表達(dá)式如下：

p=p0.

(9)

p=p0-t·10Pa/s.

(10)

(11)

2.4　平均流速方程

由于流體一般只占多孔介質(zhì)總體積的10%到50%，因此流體在孔隙通道內(nèi)的流速會超過達(dá)西流速u(大約等于它的2—10倍)，為了明白起見，設(shè)定流體在給定的孔隙空間內(nèi)的流速為平均線流速uα(也叫滲流流速)為：

uα=u/θs.

(12)

式中：θs為流體的體積部分，對于飽和介質(zhì)而言，θs=0；

2.5　相似系數(shù)

應(yīng)用型達(dá)西滲流模型有可供選擇的相似系數(shù)用來進(jìn)一步簡化分析和減少迭代次數(shù)或是進(jìn)行參數(shù)化的運(yùn)算。這種相似系數(shù)的分析模式能夠采用雙子域計算系統(tǒng)，對包括流體在裂縫中的相對快速流動，多相問題和單獨(dú)密度等進(jìn)行耦合計算。由于相似系數(shù)的加入，控制方程可采用以下形式：

(13)

式中：δs為S的相似系數(shù)；δK為滲透系數(shù)κ的相似系數(shù)；δQ為流量Q的相似系數(shù)；

由式1—式5構(gòu)成水流作用下的多孔介質(zhì)不連續(xù)裂縫中水流流動數(shù)學(xué)模型。計算參數(shù)[12-13]見表1。該模型可開展水流作用下壓實(shí)黏土封場系統(tǒng)物性參數(shù)的變化過程預(yù)測與評價研究，對垃圾填埋場封場系統(tǒng)開裂失效的研究具有一定的理論意義。

表1　計算參數(shù)

3水滲流仿真模型的建立

此模型采用了一種精確高效的方法來模擬開裂基質(zhì)環(huán)境中的裂縫和孔隙中的流動。如圖1所示，通過在相鄰基質(zhì)塊邊界處構(gòu)造裂縫來建立模型。當(dāng)設(shè)定流體在邊界上狹窄的裂縫中流動時，由達(dá)西定律來控制基質(zhì)塊中的流速把裂縫作為基質(zhì)塊內(nèi)部的邊界這種模型是非常有效的，這樣減少了創(chuàng)建一個幾何形狀長而窄即具有較高的長寬比的裂縫的必要。否則，需要建立一個由大量微小單元構(gòu)成的高密度的幾何模型。

圖1　一條裂縫的幾何模型

圖1所示為一個長寬高各為1 m的開裂的多孔介質(zhì)基質(zhì)塊。裂縫厚度為0.000 1 m，遠(yuǎn)小于基質(zhì)塊的厚度。水流在裂縫中的滲透性能遠(yuǎn)大于基質(zhì)中滲透性能。除了裂縫的出入口所在的面上以外，基質(zhì)塊其他的平面都是不可滲透的。模型中的流體從基質(zhì)塊的裂縫右上方邊緣流向左下方的邊緣。在數(shù)值計算過程中，最初狀態(tài)的流體在基質(zhì)中是不能流動的，僅當(dāng)裂縫的入口那一側(cè)出現(xiàn)壓力梯度時，流體才會在其作用下流動，出口處的壓力隨時間逐漸減小，入口處的壓力仍保持恒定。計算過程約1 000 s。

為了模擬較復(fù)雜工況下的壓實(shí)黏土裂縫中水流流動模型，人工添加了三條裂縫，對其中水流流動進(jìn)行數(shù)值分析及模擬如圖2所示。

圖2　三條裂縫的幾何模型

4計算結(jié)果與分析

水流在壓實(shí)黏土裂縫和孔隙中流動過程中所引起黏土物性參數(shù)的變化對黏土的抗壓，抗剪強(qiáng)度等工程特性造成一定的影響。此處借助多孔介質(zhì)不連續(xù)裂縫水流流動模型對水流在裂縫中流動時的壓強(qiáng)和流速變化進(jìn)行數(shù)值計算。

圖3表示了在最終輸出時間1 000 s時的流體壓力等值面圖?？梢钥闯觯黧w的壓力在裂縫和基質(zhì)之間是連續(xù)的，等壓面成彎曲狀表示流體在裂縫中的流動狀態(tài)與其在土體基質(zhì)中的流動不同。這是由于裂縫和基質(zhì)(土顆粒)具有不同的滲透系數(shù)所造成的。根據(jù)達(dá)西滲流定律可知，當(dāng)水力梯度相同時，流體在介質(zhì)中的滲流流速與介質(zhì)的滲透系數(shù)成正比[14]。同時，在此過程中流體的壓強(qiáng)所出現(xiàn)的最大值和最小值分別出現(xiàn)在入口和出口的邊緣處，大小分別為99 kPa和91.5 kPa。

圖3　裂縫中水流流動等壓面圖

由圖4可以看出，流體在裂縫中流速較大的區(qū)域主要分布在水流的入口和出口處的邊緣，這是因為流體在裂縫中流動過程中受到流體與基質(zhì)邊界的摩擦力作用，同時裂縫的寬度對流體流動起到一定的約束作用，進(jìn)而影響流體在裂縫中的流速場分布。

圖4　裂縫中水流流動的流速場的分布

由圖5分析可得，流體在基質(zhì)中和裂縫中的流速之比的分布規(guī)律在裂縫寬度方向上表現(xiàn)為各向同性。同時，由右側(cè)的數(shù)據(jù)可以看出，流體在基質(zhì)中和裂縫中的流速之比最大為0.010 1，并且都出現(xiàn)在流體的出入口兩側(cè)，表明水流在這兩處的流速與裂縫中的流速相比都比較大。這是由于流體流動方向的兩側(cè)都存在壓力梯度所導(dǎo)致的。根據(jù)達(dá)西滲流定律，同種介質(zhì)中，滲透流速與水力梯度成正比，故在兩側(cè)的都存在的水力梯度下，使得出口和入口處的滲透流速遠(yuǎn)大于裂縫中的滲透流速。同時，由于滲透系數(shù)的不同，導(dǎo)致沿著水流方向上的基質(zhì)滲透流速的變化量遠(yuǎn)大于裂縫中中流體的流速。所以在上圖中會出現(xiàn)基質(zhì)塊的內(nèi)部區(qū)域內(nèi)，兩者流速之比會出現(xiàn)最小值。

圖5　流體在基質(zhì)中和裂縫中的流速之比(平均線流速)

圖6所表現(xiàn)的是以上所有物理量在1 000 s后的運(yùn)算結(jié)果。圖中的箭頭表示的是流體在基質(zhì)和裂縫中的流動方向。形象清晰地表達(dá)了最終時間狀態(tài)下的各物理量在空間上的分布情況。由于在建立多孔介質(zhì)流動模型的時候，假設(shè)了流體在其出入口兩側(cè)的其他基質(zhì)塊表面沒有流入或流出。圖6中所出現(xiàn)的箭頭表示流體只是從入口流向出口的方向。

圖6　終態(tài)流速-壓強(qiáng)分布圖

為了準(zhǔn)確地動態(tài)分析水流在裂縫中流動狀態(tài)，采用時間長度為1 000 s的數(shù)值計算過程，其計算結(jié)果如圖7所示：

通過圖7可以看出，通過COMSOL形象的利用多孔介質(zhì)不連續(xù)裂縫水流流動模型模擬壓實(shí)黏土裂縫中水流流動條件下的流速壓強(qiáng)分布。在不同的階段，基質(zhì)中的流體流速和等壓面的壓強(qiáng)隨著時間的變化呈現(xiàn)出各自不同的規(guī)律。流體的平均線流速場與裂縫中水流流速場在開始很短一段時間內(nèi)就已經(jīng)達(dá)到穩(wěn)定狀態(tài)，其變化規(guī)律和之前分析的1 000 s時的規(guī)律一樣。等壓面壓強(qiáng)隨著時間的推移，由入口向出口依次遞減。且等壓面的曲率由基質(zhì)塊軸線向裂縫兩端遞增。表明越靠近裂縫的兩端，裂縫中的流體的運(yùn)動狀態(tài)與基質(zhì)中的流動狀態(tài)差異性越大。水流的流動狀態(tài)的變化主要體現(xiàn)在裂縫形狀的變化上。

圖7　不同時間時流速-壓強(qiáng)分布圖

5結(jié)論

通過建立多孔介質(zhì)不連續(xù)裂縫水流流動模型和設(shè)定其邊界條件來對壓實(shí)黏土裂縫中水流流動進(jìn)行數(shù)值分析，以評價水滲流條件下對壓實(shí)黏土開裂的影響，對垃圾填埋場封場系統(tǒng)的安全運(yùn)行提供理論依據(jù)，由仿真計算得出以下結(jié)論。

(1)過程中流體的壓強(qiáng)所出現(xiàn)的最大值和最小值分別出現(xiàn)在入口和出口的邊緣處，大小分別為99 kPa和91.5 kPa；

(2)流體在裂縫中流速較大的區(qū)域主要分布在水流的入口和出口處的邊緣；由于流體與基質(zhì)邊界的摩擦力作用，以及受到裂縫寬度的影響，對流體流動起到一定的約束作用，進(jìn)而影響流體在裂縫中的流速場分布.

(3)流體在基質(zhì)中和裂縫中的流速之比的分布規(guī)律在裂縫寬度方向上表現(xiàn)為各向同性，流體在基質(zhì)中和裂縫中的流速的比值最大為0.010 1，由于滲透系數(shù)的不同，導(dǎo)致沿著水流方向上的基質(zhì)滲透流速的變化量遠(yuǎn)大于裂縫中中流體的流速。

(4)等壓面壓強(qiáng)隨著時間的推移，由入口向出口依次遞減，且等壓面的曲率由基質(zhì)塊軸線向裂縫兩端遞增，表明越靠近裂縫的兩端，裂縫中的流體的運(yùn)動狀態(tài)與基質(zhì)中的流動狀態(tài)差異性越大。

參考文獻(xiàn)：

[1]浦辛剛,陳新民.填埋場襯墊系統(tǒng)的評價[J]. 南京工業(yè)大學(xué)學(xué)報,2004,26 (2):33-38.

[2]Baer J U, Kent T F, Anderson S H. Image analysis and fractal geometry to characterize soil desiccation cracks[J].Geoderma,2009.

[3]Tang C, Shi B, Liu C,et al. Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[J]. Engineering Geology, 2008, 101: 204-217.

[4]袁俊平.非飽和膨脹土的裂隙概化模型與邊坡穩(wěn)定研究[D].南京：河海大學(xué),2003.

[5]劉建國,聶永豐.填埋場黏土襯里破壞機(jī)理研究[J].城市環(huán)境與城市生態(tài), 2000, 13(6): 51-53.

[6]錢學(xué)德,郭志平.填埋場黏土襯墊的設(shè)計與施工[J]. 水利水電科技發(fā)展, 1997, 17(4): 55-59.

[7]Mathias S A, Todman L C,Step-Drawdown Tests and the forchheimer Equation[J]. Water Resources Resrarch,2010,46(7):1-9.

[8]唐朝生，施斌.干濕循環(huán)過程中膨脹土的脹縮變形特征[J].巖土工程學(xué)報，2011,33(9):1376-1384.

[9]Pzbek A.Investigation of the effects of wetting-drying and freezing-thawing cycles on some physical and mechanical properties of selected ignimbrites[J].Bulletin of Engineering Geology and the Environment,2014,73(2):595-609.

[10]Peron H, Delenne J Y, Laloui L,et al. Discrete element modelling of drying shrinkage and cracking of soils[J].Computers and Geotechnics,2009, 36:61-69.

[11]王剛,安林. COMSOL Multiphysics 工程實(shí)踐與理論仿真[M].北京:電子工業(yè)出版社, 2012.

[12]徐國建,沈揚(yáng),劉漢龍. 孔隙率、級配參數(shù)對粉土雙軸壓縮性狀影響的顆粒流分析[J]. 巖土力學(xué), 2013, 34(11):3321-3328.

[13]趙強(qiáng), 李皋, 肖貴林,等. 單裂隙滲流有限元數(shù)值仿真研究[J]. 水資源與水工程學(xué)報, 2014, 25(2):195-199.

[14]翟云芳. 滲流力學(xué)[M], 北京:石油工業(yè)出版社, 2011.

Mathematical simulation of cracking clay in landfill liner under water seepage

DONGYi-qie1,LUHai-jun1，DAIRui1,HELi2

(1.Institute of Poromechanics, Wuhan Polytechnic University, Wuhan 430023，China;2.Department of

Mechatronic Engineering College,College of Arts and Science of Jianghan University,Wuhan 430056,China)

Abstract：Aiming at the problem that the landfill compacted clay cover system is vulnerable to appear desiccation cracking in drought or semi-arid area, based on the Hydrodynamic flow equation, simulate to analyze the mathematical model of water flow in discontinuous cracks in the porous medium, simulate the water seepage rule in crack clay by COMSOL. Based on the results of mathematical simulation, the maximum and minimum values of fluid pressure appear in the inlet and outlet edges, the large area of flow velocity in fluid appear in the inlet and outlet edge in the flow, and closer to the end of cracks, greater the difference in state of motion in fracture fluid and current state of the matrix.

Key words：mathematical simulation; compacted clay; drying and cracking

中圖分類號：TP 391.9；TV139.14

文獻(xiàn)標(biāo)識碼：A