圍壓條件下原狀花崗巖殘積土細(xì)觀滲流數(shù)值模擬

蔡沛辰， 闕云， 李顯

(福州大學(xué)土木工程學(xué)院， 福建 福州 350108)

0 引言

花崗巖殘積土主要分布于我國東南沿海地區(qū)， 屬于部分連續(xù)的多孔介質(zhì)材料[1]. 飽和土體中的水分在孔隙中發(fā)生不均衡流動的現(xiàn)象， 稱為滲流. 滲流現(xiàn)象普遍存在于實(shí)際工程項(xiàng)目中， 并常伴有失穩(wěn)、 滲透變形等不利影響. 傳統(tǒng)的滲流研究多集中于宏觀尺度， 用試驗(yàn)手段測得滲透系數(shù)， 來表征土體的滲流特性， 但這種研究無法直觀準(zhǔn)確地描述孔隙結(jié)構(gòu)特征、 流體滲流路徑以及孔隙結(jié)構(gòu)與滲流特性之間的內(nèi)在關(guān)系[2]. 因此， 在細(xì)觀尺度上對土體進(jìn)行滲流特性研究顯得尤為必要.

目前， 細(xì)觀滲流的研究方法主要有物理試驗(yàn)和數(shù)值模擬兩種. 在物理試驗(yàn)方面， 如張明鳴等[3]采用滲透系數(shù)試驗(yàn)， 對淤泥土的細(xì)觀滲流本質(zhì)進(jìn)行了研究； 梁健偉[4]通過對人工土采用微電場效應(yīng)試驗(yàn)， 研究了微電場變化對土體細(xì)觀滲透特性的影響. 雖然物理試驗(yàn)研究方法可以對多孔介質(zhì)滲流的一些細(xì)觀機(jī)理定性研究， 但缺少科學(xué)的理論對其進(jìn)行精細(xì)的定量表述. 近年來， 數(shù)值模擬在細(xì)觀滲流方面的研究主要采用微孔流理論、 孔隙網(wǎng)絡(luò)模型和格子Boltzmann等方法， 如： ① 基于微孔流理論， 趙延林等[5]、 李晶晶等[6]和馮上鑫等[7]對巖石和土石混合體從細(xì)觀角度進(jìn)行了滲流場數(shù)值模擬， 分析了細(xì)觀滲流的特性； ② 基于孔隙網(wǎng)絡(luò)模型方法， Qren等[8]通過構(gòu)建大量的砂巖孔隙網(wǎng)絡(luò)模型， 進(jìn)行了滲流模擬， 并通過模擬結(jié)果預(yù)測了相對滲透率； ③ 基于格子Boltzmann方法， Shen等[9]、 崔冠哲等[10]建立了水壓力作用下土體和重構(gòu)土體細(xì)觀滲流場的二維模型， 分析了整體和局部滲流速度的變化規(guī)律. 上述相關(guān)研究已取得了豐碩成果， 但一方面多側(cè)重于巖石細(xì)觀滲流的研究， 且研究對象多為隨機(jī)建立的孔隙結(jié)構(gòu)模型， 與實(shí)際孔隙結(jié)構(gòu)有一定差異， 對實(shí)踐不能起到很好的指導(dǎo)意義. 另一方面， 只在巖石細(xì)觀滲流研究中有考慮圍壓的影響， 而事實(shí)上， 任何土體四周不可避免地都存在著圍壓， 圍壓將如何影響土體細(xì)觀滲流， 還有待探究.

本研究以原狀花崗巖殘積土為研究對象， 對其進(jìn)行工業(yè)CT掃描， 建立真實(shí)反映土體孔隙結(jié)構(gòu)的幾何模型. 采用COMSOL Multiphyscis高級數(shù)值模擬仿真軟件進(jìn)行數(shù)值模擬， 直觀展現(xiàn)各孔隙區(qū)域的滲流速度分布情況， 并分析考慮和不考慮圍壓時(shí)的細(xì)觀滲流特性， 以期探究和認(rèn)識不同條件下花崗巖殘積土的細(xì)觀滲流規(guī)律.

1 CT掃描及細(xì)觀模型構(gòu)建

1.1 CT掃描及圖像處理

圖1 現(xiàn)場取樣Fig.1 Field sampling

試驗(yàn)原狀土選自福州市某地山坡， 現(xiàn)場取樣如圖1所示.

劉勇等[11]研究了在掃描土壤孔隙結(jié)構(gòu)特征上醫(yī)學(xué)CT和工業(yè)CT的區(qū)別， 結(jié)果表明， 在土壤橫斷面圖像掃描質(zhì)量方面， 工業(yè)CT要明顯高于醫(yī)學(xué)CT. 故本研究對處理后的原狀土采用英華公司工業(yè)CT掃描儀進(jìn)行XZ向掃描， 得到能真實(shí)體現(xiàn)土壤孔隙分布的CT掃描圖像. 通過Image J軟件中的Threshold功能， 對圖像進(jìn)行二值化處理， 形成只包含黑白兩色的圖像， 再利用中值濾波對二值化圖像進(jìn)行降噪處理， 除去圖像中獨(dú)立的噪點(diǎn)， 最終得到的一系列15 cm×15 cm的二維掃描切片， 如圖2所示. 選取中間位置的切片， 截取其中孔隙連通性較好的區(qū)域作為模擬對象， 局部放大如圖3所示， 大小為6.0 mm×3.3 mm.

圖2 二維掃描切片F(xiàn)ig.2 2D scanning slice

圖3 文中所采用的二值化圖像Fig.3 Binary image used in this paper

1.2 細(xì)觀幾何模型

本研究所采用的COMSOL Multiphyscis是一款基于偏微分方程的高級數(shù)值模擬仿真軟件， 具有多場耦合、 2D和3D模型直接導(dǎo)入、 模型網(wǎng)格變形、 粒子追蹤及動態(tài)可視化結(jié)果處理等優(yōu)點(diǎn). 為轉(zhuǎn)化為COMSOL中可導(dǎo)入的二維模型文件格式， 對選取的二值化圖像進(jìn)行一系列處理[12]， 過程參見圖4. 構(gòu)建的細(xì)觀幾何模型， 如圖5所示， 灰色為基質(zhì)區(qū)域， 藍(lán)色為孔隙區(qū)域.

圖4 幾何模型構(gòu)建路線Fig.4 Geometric model and generation route

圖5 細(xì)觀幾何模型Fig.5 Mesoscopic geometric model

2 細(xì)觀滲流仿真模擬

2.1 控制方程

2.1.1 細(xì)觀滲流控制方程

鑒于CT掃描精度較高， 故不考慮流體在土體基質(zhì)中的流動， 可認(rèn)為流體只在孔隙和喉道構(gòu)成的通道中進(jìn)行流動[13]. 用不可壓縮N-S方程描述流體流動[14]， 控制方程如下

(1)

其中：ρ為流體密度， kg·m-3；p為壓力， Pa；I為單位矩陣；μ為流體動力粘度， Pa·s；F為體積力， N·m-3；u為流速， m·s-1；Δ為梯度算子；g為重力加速度， m·s-2.

2.1.2 應(yīng)力場控制方程

研究需要得到變形后孔隙結(jié)構(gòu)， 故采用幾何非線性步驟， 應(yīng)力為第二類Piola-Kirchhoff應(yīng)力， 控制方程如下

Δ·(FaS)T+FV=0

(2)

Ip+Δr=Fa

(3)

其中：FV為相對于未變形體積給出的體積力分量；Fa為變形梯度；S為第二類Piola-Kirchhoff應(yīng)力張量；Ip為單位張量；r為位移場.

2.2 材料屬性

采用水作為滲流流體， 為模擬圍壓作用下的土體變形， 分別用兩種線彈性材料區(qū)分孔隙和基質(zhì)區(qū)域， 具體材料屬性如表1所示. 其中， 基質(zhì)區(qū)域楊氏模量數(shù)值通過查閱相關(guān)文獻(xiàn)[15-16]獲得.

表1 材料屬性

2.3 邊界及初始條件

細(xì)觀滲流模擬中， 依據(jù)構(gòu)建的幾何模型， 進(jìn)行邊界條件設(shè)定： 上下邊界為進(jìn)出口， 左右邊界為不透水層， 壁為無滑移條件， 同時(shí)考慮重力影響因素. 應(yīng)力場模擬中， 模型上邊界和左右邊界施加圍壓， 下邊界設(shè)為固定約束.

3 結(jié)果與討論

3.1 數(shù)值模型驗(yàn)證

為驗(yàn)證文中數(shù)值模型的正確性， 采用文獻(xiàn)[17]的方法， 驗(yàn)證幾何模型中滲流是否符合Darcy定律. 流體通過水壓力進(jìn)行驅(qū)動， 從上邊界進(jìn)入孔隙區(qū)域， 到下邊界流出， 并在孔隙區(qū)域內(nèi)保持恒定的水壓差. 整個(gè)孔隙區(qū)域的速度大小變化如圖6所示(以100 Pa水壓差為例). 從圖6中可以發(fā)現(xiàn)： 在滲流路徑上存在一些高速區(qū)域， 計(jì)算可知整個(gè)區(qū)域的平均速度大小為2.98 mm·s-1， 平均體流速大小為19.1 mm3·s-1.

通過改變模型邊界的水壓差， 計(jì)算相應(yīng)的體積流量， 對其進(jìn)行參數(shù)化研究. 圖7表示水壓差和流量的變化情況， 從圖7中可以觀察到： 流量隨施加的水壓差呈線性變化， 符合Darcy定律， 也驗(yàn)證了該數(shù)值模型的正確性.

圖6 XZ向模型速度云圖Fig.6 XZ model velocity cloud

圖7 體積流量隨所施加水壓差的變化圖Fig.7 Variation of volume flow rate with applied water pressure

3.2 模型網(wǎng)格變形及應(yīng)力場分布

根據(jù)土力學(xué)原理[18]可知， 土體作用豎向均布荷載時(shí)， 其側(cè)壓力p0與豎向力pz之間的關(guān)系可表述如下

p0=K0pz+K0γz

(4)

K0=1-sinφ

(5)

圖8 圍壓作用下XZ模型應(yīng)力云圖Fig.8 Stress cloud diagram of XZ model under confining pressure

式中：K0為靜止土壓力系數(shù)；γ為土體重度， 此處γ=ρg=17.938 kN·m-3；z為豎向深度， mm；φ(°)為內(nèi)摩擦角， 查閱相關(guān)文獻(xiàn)[19]， 取其平均值26.7°.

在模型上邊界， 施加均布豎向力1 kPa時(shí)， 為簡化計(jì)算， 且考慮到土樣深度僅3.3 mm， 故忽略深度對側(cè)壓力的影響， 近似取550 Pa， 應(yīng)力分布如圖8所示. 從圖8中可以看出： 孔隙區(qū)域內(nèi)部應(yīng)力在1 kPa左右， 但在孔隙介質(zhì)處, 應(yīng)力集中在3～5 kPa， 局部更高達(dá)206 kPa. 此外， 由受力邊界可看出： 土體在受到圍壓作用后， 存在孔隙的部分變形程度明顯大于基質(zhì)部分， 計(jì)算知孔隙率為31.71%， 降低了0.51%.

3.3 土體細(xì)觀滲流分析

3.3.1 計(jì)算滲透率

滲透率是指在一定壓差下， 巖石或土體本身允許流體通過的能力[20]. 依據(jù)達(dá)西定律， 計(jì)算滲透率K表述為

(6)

其中：K為計(jì)算滲透率， mD；u為平均流速， m·s-1；μ為流體動力粘度， Pa·s；L為滲流路徑長度， m， 此處為3.3 mm， 確定方法同文獻(xiàn)[13]； Δp為水壓差， Pa.

不考慮圍壓時(shí)， 細(xì)觀滲流結(jié)果分析如下. 在出口邊界(Z=0)處畫截線， 沿線滲流速度如圖9所示. 再將速度積分后除以整個(gè)區(qū)域?qū)挾鹊玫狡骄俣萿， 最后根據(jù)公式(6)得到計(jì)算滲透率.

對出入口邊界設(shè)置不同的水壓力， 模擬不同水壓差下的滲流速度場分布情況， 按上述步驟得到計(jì)算滲透率， 水壓差與計(jì)算滲透率的關(guān)系如圖10所示， 從圖10中可以發(fā)現(xiàn)： 計(jì)算滲透率隨水壓差的增大逐漸減小， 且近似成線性關(guān)系， 擬合表達(dá)式為y=0.002 48-1.830 73×10-8x， 相關(guān)系數(shù)為0.994 55.

圖9 XZ向模型出口邊界速度圖Fig.9 XZ direction model exit boundary velocity map

圖10 不同水壓差下計(jì)算滲透率Fig.10 Calculated permeability under different water pressure differences

圖11 圍壓作用下的計(jì)算滲透率Fig.11 Calculated permeability under confining pressure

考慮圍壓時(shí)， 細(xì)觀滲流結(jié)果分析如下. 在邊界依次分級施加圍壓， 并用前文方法求得計(jì)算滲透率. 圍壓與計(jì)算滲透率之間的關(guān)系如圖11所示. 從圖11中可以發(fā)現(xiàn)： 6組直線呈先急劇到平緩下降， 再急劇上升， 最后趨于平緩下降的變化規(guī)律. 豎向力0～100 Pa時(shí)， 計(jì)算滲透率降低最為顯著， 達(dá)52%左右， 分析原因是： 孔道受圍壓作用而變窄， 從而通過孔隙截面的流量減少. 豎向力300～400 Pa時(shí)， 計(jì)算滲透率反而明顯增大， 分析原因是： 模型中的部分孤立孔隙受圍壓的擠壓作用， 使其與原本連通孔道相貫通， 進(jìn)而致使通過孔隙截面的流量增加.

3.3.2 速度場下的粒子追蹤

在細(xì)觀滲流速度場的基礎(chǔ)上， 運(yùn)用COMSOL后處理模塊中的粒子追蹤方法， 粒子起點(diǎn)控制參數(shù)表達(dá)式為X： range(0, 0.005, 6)，Z： 3.3， 單位： mm， 繪制考慮和不考慮圍壓作用下， 隨時(shí)間變化速度場下的粒子追蹤圖像， 如圖12、 13所示(1 kPa水壓差). 不考慮圍壓條件下， 由圖12(a)可知， 流體在滲流路徑節(jié)點(diǎn)O處， 有孔道a和b可進(jìn)入， 但流體優(yōu)先經(jīng)孔道b進(jìn)入到孔隙B中， 相類似的還有圖13(a)中節(jié)點(diǎn)C和C′， 流體優(yōu)先滲流進(jìn)入到孔隙D和D′中； 從圖13(a)可以發(fā)現(xiàn)， 滲流路徑2-2′近乎筆直， 相比其他路徑孔道彎曲程度更低， 滲流速度相應(yīng)較大. 考慮圍壓條件下， 如圖12(b)和13(b)也可發(fā)現(xiàn)上述相類似現(xiàn)象. 綜合上述現(xiàn)象表明： 考慮和不考慮圍壓時(shí)， 流體滲流速度都受孔道彎曲程度控制， 且流體優(yōu)先選擇連通性較好的孔道進(jìn)行滲流.

分別對比圖12(a)、 (b)和圖13(a)、 (b)可知： 同一時(shí)刻， 不考慮圍壓時(shí)的實(shí)際滲流路徑長度明顯大于考慮圍壓時(shí)的實(shí)際滲流長度， 即前者滲流速度大于后者. 經(jīng)計(jì)算， 考慮圍壓時(shí)滲流長度減少了16%左右， 速度大小降低了13%左右.

圖12 0.01 s時(shí)刻速度場下的粒子追蹤圖像Fig.12 Particle tracking image in the velocity field at 0.01 s moment

圖13 0.05 s時(shí)刻速度場下的粒子追蹤圖像Fig.13 Particle tracking image in the velocity field at 0.05 s moment

3.3.3 應(yīng)力敏感性指數(shù)

土體對外界的敏感程度用敏感性指數(shù)進(jìn)行評價(jià). 本研究以花崗巖殘積土滲透率為例， 用應(yīng)力敏感性指數(shù)來評價(jià)其受圍壓作用時(shí)的敏感程度[21]， 定義為

(7)

圖14 應(yīng)力敏感性指數(shù)變化圖Fig.14 Stress sensitivity index change chart

式中： SIp為應(yīng)力敏感性指數(shù)；K0為圍壓為0時(shí)的計(jì)算滲透率；Kp為圍壓為p時(shí)的計(jì)算滲透率.

圖14總結(jié)了不同圍壓作用下應(yīng)力敏感性指數(shù)的變化情況. 從圖中可以看出： 6組直線變化差別不大， 且呈“Z”字型變化. 豎向力300～400 Pa時(shí)， 應(yīng)力敏感性指數(shù)降低最為顯著， 降低幅度為82.5%. 上述現(xiàn)象表明： 水壓差對應(yīng)力敏感性指數(shù)影響甚微， 而圍壓對它影響較大. 分析原因是： 土體具有可塑性， 對土體進(jìn)行外部加載， 使其內(nèi)部形成了新的孔隙結(jié)構(gòu). 再結(jié)合圖11分析可知： 應(yīng)力敏感性指數(shù)越大， 土體孔道變窄的程度越大. 反之， 應(yīng)力敏感性指數(shù)越小， 土體中孤立孔隙與連通孔道相貫通的程度越大.

4 結(jié)語

1) 兩種條件下， 計(jì)算滲透率都隨水壓差的增加， 逐漸降低， 且基本成線性關(guān)系. 考慮圍壓時(shí)， 計(jì)算滲透率隨豎向力呈先急劇后平緩下降， 再急劇上升， 最后趨于平緩下降的變化規(guī)律. 急劇下降段的原因?yàn)榭椎朗車鷫鹤饔枚冋?從而通過孔隙截面的流量減少； 而急劇上升段的原因?yàn)槟Ｐ椭械牟糠止铝⒖紫妒車鷫旱臄D壓作用， 使得其與原本連通孔道相貫通， 進(jìn)而致使通過孔隙截面的流量增加.

2) 同一時(shí)刻， 考慮圍壓(1 kPa豎向力)時(shí)的滲流長度較未考慮時(shí)的減少了16%左右， 速度大小降低了13%左右. 兩種條件下， 流體滲流速度都受孔道彎曲程度控制， 且流體優(yōu)先選擇連通性較好的孔道進(jìn)行滲流.

3) 應(yīng)力敏感性指數(shù)主要受圍壓控制， 水壓差對其影響甚微. 不同圍壓作用下應(yīng)力敏感性指數(shù)呈“Z”字型變化. 其中， 豎向力300～400 Pa時(shí)， 應(yīng)力敏感性指數(shù)降低最為顯著， 降低幅度為82.5%. 應(yīng)力敏感性指數(shù)增加， 土體孔道變窄的程度越大； 而應(yīng)力敏感性指數(shù)減小， 土體中孤立孔隙與連通孔道相貫通的程度越大.