水流作用下填埋場(chǎng)開(kāi)裂黏土襯墊數(shù)值仿真計(jì)算

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

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

?

水流作用下填埋場(chǎng)開(kāi)裂黏土襯墊數(shù)值仿真計(jì)算

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

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

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

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

1引言

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

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

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

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

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

2.1　流體流速方程

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

(1)

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

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

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

(2)

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

對(duì)不可壓縮流體來(lái)說(shuō)，ρf可以從等式兩端約去，因此控制式2可化為：

(3)

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

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

以時(shí)間為變量的土體基質(zhì)中的水流流動(dòng)方程可以通過(guò)達(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ì)裂縫中水流流動(dòng)方程

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

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

(7)

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

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

(8)

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

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

p=p0.

(9)

p=p0-t·10Pa/s.

(10)

(11)

2.4　平均流速方程

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

uα=u/θs.

(12)

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

2.5　相似系數(shù)

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

(13)

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

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

表1　計(jì)算參數(shù)

3水滲流仿真模型的建立

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

圖1　一條裂縫的幾何模型

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

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

圖2　三條裂縫的幾何模型

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

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

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

圖3　裂縫中水流流動(dòng)等壓面圖

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

圖4　裂縫中水流流動(dòng)的流速場(chǎng)的分布

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

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

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

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

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

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

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

5結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

[1]浦辛剛,陳新民.填埋場(chǎng)襯墊系統(tǒng)的評(píng)價(jià)[J]. 南京工業(yè)大學(xué)學(xué)報(bào),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]劉建國(guó),聶永豐.填埋場(chǎng)黏土襯里破壞機(jī)理研究[J].城市環(huán)境與城市生態(tài), 2000, 13(6): 51-53.

[6]錢(qián)學(xué)德,郭志平.填埋場(chǎng)黏土襯墊的設(shè)計(jì)與施工[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)過(guò)程中膨脹土的脹縮變形特征[J].巖土工程學(xué)報(bào)，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]徐國(guó)建,沈揚(yáng),劉漢龍. 孔隙率、級(jí)配參數(shù)對(duì)粉土雙軸壓縮性狀影響的顆粒流分析[J]. 巖土力學(xué), 2013, 34(11):3321-3328.

[13]趙強(qiáng), 李皋, 肖貴林,等. 單裂隙滲流有限元數(shù)值仿真研究[J]. 水資源與水工程學(xué)報(bào), 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

中圖分類(lèi)號(hào)：TP 391.9；TV139.14

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