围压条件下原状花岗岩残积土细观渗流数值模拟

蔡沛辰， 阙云， 李显

(福州大学土木工程学院， 福建 福州 350108)

0 引言

花岗岩残积土主要分布于我国东南沿海地区， 属于部分连续的多孔介质材料[1]. 饱和土体中的水分在孔隙中发生不均衡流动的现象， 称为渗流. 渗流现象普遍存在于实际工程项目中， 并常伴有失稳、 渗透变形等不利影响. 传统的渗流研究多集中于宏观尺度， 用试验手段测得渗透系数， 来表征土体的渗流特性， 但这种研究无法直观准确地描述孔隙结构特征、 流体渗流路径以及孔隙结构与渗流特性之间的内在关系[2]. 因此， 在细观尺度上对土体进行渗流特性研究显得尤为必要.

目前， 细观渗流的研究方法主要有物理试验和数值模拟两种. 在物理试验方面， 如张明鸣等[3]采用渗透系数试验， 对淤泥土的细观渗流本质进行了研究； 梁健伟[4]通过对人工土采用微电场效应试验， 研究了微电场变化对土体细观渗透特性的影响. 虽然物理试验研究方法可以对多孔介质渗流的一些细观机理定性研究， 但缺少科学的理论对其进行精细的定量表述. 近年来， 数值模拟在细观渗流方面的研究主要采用微孔流理论、 孔隙网络模型和格子Boltzmann等方法， 如： ① 基于微孔流理论， 赵延林等[5]、 李晶晶等[6]和冯上鑫等[7]对岩石和土石混合体从细观角度进行了渗流场数值模拟， 分析了细观渗流的特性； ② 基于孔隙网络模型方法， Qren等[8]通过构建大量的砂岩孔隙网络模型， 进行了渗流模拟， 并通过模拟结果预测了相对渗透率； ③ 基于格子Boltzmann方法， Shen等[9]、 崔冠哲等[10]建立了水压力作用下土体和重构土体细观渗流场的二维模型， 分析了整体和局部渗流速度的变化规律. 上述相关研究已取得了丰硕成果， 但一方面多侧重于岩石细观渗流的研究， 且研究对象多为随机建立的孔隙结构模型， 与实际孔隙结构有一定差异， 对实践不能起到很好的指导意义. 另一方面， 只在岩石细观渗流研究中有考虑围压的影响， 而事实上， 任何土体四周不可避免地都存在着围压， 围压将如何影响土体细观渗流， 还有待探究.

本研究以原状花岗岩残积土为研究对象， 对其进行工业CT扫描， 建立真实反映土体孔隙结构的几何模型. 采用COMSOL Multiphyscis高级数值模拟仿真软件进行数值模拟， 直观展现各孔隙区域的渗流速度分布情况， 并分析考虑和不考虑围压时的细观渗流特性， 以期探究和认识不同条件下花岗岩残积土的细观渗流规律.

1 CT扫描及细观模型构建

1.1 CT扫描及图像处理

图1 现场取样Fig.1 Field sampling

试验原状土选自福州市某地山坡， 现场取样如图1所示.

刘勇等[11]研究了在扫描土壤孔隙结构特征上医学CT和工业CT的区别， 结果表明， 在土壤横断面图像扫描质量方面， 工业CT要明显高于医学CT. 故本研究对处理后的原状土采用英华公司工业CT扫描仪进行XZ向扫描， 得到能真实体现土壤孔隙分布的CT扫描图像. 通过Image J软件中的Threshold功能， 对图像进行二值化处理， 形成只包含黑白两色的图像， 再利用中值滤波对二值化图像进行降噪处理， 除去图像中独立的噪点， 最终得到的一系列15 cm×15 cm的二维扫描切片， 如图2所示. 选取中间位置的切片， 截取其中孔隙连通性较好的区域作为模拟对象， 局部放大如图3所示， 大小为6.0 mm×3.3 mm.

图2 二维扫描切片Fig.2 2D scanning slice

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

1.2 细观几何模型

本研究所采用的COMSOL Multiphyscis是一款基于偏微分方程的高级数值模拟仿真软件， 具有多场耦合、 2D和3D模型直接导入、 模型网格变形、 粒子追踪及动态可视化结果处理等优点. 为转化为COMSOL中可导入的二维模型文件格式， 对选取的二值化图像进行一系列处理[12]， 过程参见图4. 构建的细观几何模型， 如图5所示， 灰色为基质区域， 蓝色为孔隙区域.

图4 几何模型构建路线Fig.4 Geometric model and generation route

图5 细观几何模型Fig.5 Mesoscopic geometric model

2 细观渗流仿真模拟

2.1 控制方程

2.1.1 细观渗流控制方程

鉴于CT扫描精度较高， 故不考虑流体在土体基质中的流动， 可认为流体只在孔隙和喉道构成的通道中进行流动[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 应力场控制方程

研究需要得到变形后孔隙结构， 故采用几何非线性步骤， 应力为第二类Piola-Kirchhoff应力， 控制方程如下

Δ·(FaS)T+FV=0

(2)

Ip+Δr=Fa

(3)

其中：FV为相对于未变形体积给出的体积力分量；Fa为变形梯度；S为第二类Piola-Kirchhoff应力张量；Ip为单位张量；r为位移场.

2.2 材料属性

采用水作为渗流流体， 为模拟围压作用下的土体变形， 分别用两种线弹性材料区分孔隙和基质区域， 具体材料属性如表1所示. 其中， 基质区域杨氏模量数值通过查阅相关文献[15-16]获得.

表1 材料属性

2.3 边界及初始条件

细观渗流模拟中， 依据构建的几何模型， 进行边界条件设定： 上下边界为进出口， 左右边界为不透水层， 壁为无滑移条件， 同时考虑重力影响因素. 应力场模拟中， 模型上边界和左右边界施加围压， 下边界设为固定约束.

3 结果与讨论

3.1 数值模型验证

为验证文中数值模型的正确性， 采用文献[17]的方法， 验证几何模型中渗流是否符合Darcy定律. 流体通过水压力进行驱动， 从上边界进入孔隙区域， 到下边界流出， 并在孔隙区域内保持恒定的水压差. 整个孔隙区域的速度大小变化如图6所示(以100 Pa水压差为例). 从图6中可以发现： 在渗流路径上存在一些高速区域， 计算可知整个区域的平均速度大小为2.98 mm·s-1， 平均体流速大小为19.1 mm3·s-1.

通过改变模型边界的水压差， 计算相应的体积流量， 对其进行参数化研究. 图7表示水压差和流量的变化情况， 从图7中可以观察到： 流量随施加的水压差呈线性变化， 符合Darcy定律， 也验证了该数值模型的正确性.

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

图7 体积流量随所施加水压差的变化图Fig.7 Variation of volume flow rate with applied water pressure

3.2 模型网格变形及应力场分布

根据土力学原理[18]可知， 土体作用竖向均布荷载时， 其侧压力p0与竖向力pz之间的关系可表述如下

p0=K0pz+K0γz

(4)

K0=1-sinφ

(5)

图8 围压作用下XZ模型应力云图Fig.8 Stress cloud diagram of XZ model under confining pressure

式中：K0为静止土压力系数；γ为土体重度， 此处γ=ρg=17.938 kN·m-3；z为竖向深度， mm；φ(°)为内摩擦角， 查阅相关文献[19]， 取其平均值26.7°.

在模型上边界， 施加均布竖向力1 kPa时， 为简化计算， 且考虑到土样深度仅3.3 mm， 故忽略深度对侧压力的影响， 近似取550 Pa， 应力分布如图8所示. 从图8中可以看出： 孔隙区域内部应力在1 kPa左右， 但在孔隙介质处, 应力集中在3～5 kPa， 局部更高达206 kPa. 此外， 由受力边界可看出： 土体在受到围压作用后， 存在孔隙的部分变形程度明显大于基质部分， 计算知孔隙率为31.71%， 降低了0.51%.

3.3 土体细观渗流分析

3.3.1 计算渗透率

渗透率是指在一定压差下， 岩石或土体本身允许流体通过的能力[20]. 依据达西定律， 计算渗透率K表述为

(6)

其中：K为计算渗透率， mD；u为平均流速， m·s-1；μ为流体动力粘度， Pa·s；L为渗流路径长度， m， 此处为3.3 mm， 确定方法同文献[13]； Δp为水压差， Pa.

不考虑围压时， 细观渗流结果分析如下. 在出口边界(Z=0)处画截线， 沿线渗流速度如图9所示. 再将速度积分后除以整个区域宽度得到平均速度u， 最后根据公式(6)得到计算渗透率.

对出入口边界设置不同的水压力， 模拟不同水压差下的渗流速度场分布情况， 按上述步骤得到计算渗透率， 水压差与计算渗透率的关系如图10所示， 从图10中可以发现： 计算渗透率随水压差的增大逐渐减小， 且近似成线性关系， 拟合表达式为y=0.002 48-1.830 73×10-8x， 相关系数为0.994 55.

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

图10 不同水压差下计算渗透率Fig.10 Calculated permeability under different water pressure differences

图11 围压作用下的计算渗透率Fig.11 Calculated permeability under confining pressure

考虑围压时， 细观渗流结果分析如下. 在边界依次分级施加围压， 并用前文方法求得计算渗透率. 围压与计算渗透率之间的关系如图11所示. 从图11中可以发现： 6组直线呈先急剧到平缓下降， 再急剧上升， 最后趋于平缓下降的变化规律. 竖向力0～100 Pa时， 计算渗透率降低最为显著， 达52%左右， 分析原因是： 孔道受围压作用而变窄， 从而通过孔隙截面的流量减少. 竖向力300～400 Pa时， 计算渗透率反而明显增大， 分析原因是： 模型中的部分孤立孔隙受围压的挤压作用， 使其与原本连通孔道相贯通， 进而致使通过孔隙截面的流量增加.

3.3.2 速度场下的粒子追踪

在细观渗流速度场的基础上， 运用COMSOL后处理模块中的粒子追踪方法， 粒子起点控制参数表达式为X： range(0, 0.005, 6)，Z： 3.3， 单位： mm， 绘制考虑和不考虑围压作用下， 随时间变化速度场下的粒子追踪图像， 如图12、 13所示(1 kPa水压差). 不考虑围压条件下， 由图12(a)可知， 流体在渗流路径节点O处， 有孔道a和b可进入， 但流体优先经孔道b进入到孔隙B中， 相类似的还有图13(a)中节点C和C′， 流体优先渗流进入到孔隙D和D′中； 从图13(a)可以发现， 渗流路径2-2′近乎笔直， 相比其他路径孔道弯曲程度更低， 渗流速度相应较大. 考虑围压条件下， 如图12(b)和13(b)也可发现上述相类似现象. 综合上述现象表明： 考虑和不考虑围压时， 流体渗流速度都受孔道弯曲程度控制， 且流体优先选择连通性较好的孔道进行渗流.

分别对比图12(a)、 (b)和图13(a)、 (b)可知： 同一时刻， 不考虑围压时的实际渗流路径长度明显大于考虑围压时的实际渗流长度， 即前者渗流速度大于后者. 经计算， 考虑围压时渗流长度减少了16%左右， 速度大小降低了13%左右.

图12 0.01 s时刻速度场下的粒子追踪图像Fig.12 Particle tracking image in the velocity field at 0.01 s moment

图13 0.05 s时刻速度场下的粒子追踪图像Fig.13 Particle tracking image in the velocity field at 0.05 s moment

3.3.3 应力敏感性指数

土体对外界的敏感程度用敏感性指数进行评价. 本研究以花岗岩残积土渗透率为例， 用应力敏感性指数来评价其受围压作用时的敏感程度[21]， 定义为

(7)

图14 应力敏感性指数变化图Fig.14 Stress sensitivity index change chart

式中： SIp为应力敏感性指数；K0为围压为0时的计算渗透率；Kp为围压为p时的计算渗透率.

图14总结了不同围压作用下应力敏感性指数的变化情况. 从图中可以看出： 6组直线变化差别不大， 且呈“Z”字型变化. 竖向力300～400 Pa时， 应力敏感性指数降低最为显著， 降低幅度为82.5%. 上述现象表明： 水压差对应力敏感性指数影响甚微， 而围压对它影响较大. 分析原因是： 土体具有可塑性， 对土体进行外部加载， 使其内部形成了新的孔隙结构. 再结合图11分析可知： 应力敏感性指数越大， 土体孔道变窄的程度越大. 反之， 应力敏感性指数越小， 土体中孤立孔隙与连通孔道相贯通的程度越大.

4 结语

1) 两种条件下， 计算渗透率都随水压差的增加， 逐渐降低， 且基本成线性关系. 考虑围压时， 计算渗透率随竖向力呈先急剧后平缓下降， 再急剧上升， 最后趋于平缓下降的变化规律. 急剧下降段的原因为孔道受围压作用而变窄， 从而通过孔隙截面的流量减少； 而急剧上升段的原因为模型中的部分孤立孔隙受围压的挤压作用， 使得其与原本连通孔道相贯通， 进而致使通过孔隙截面的流量增加.

2) 同一时刻， 考虑围压(1 kPa竖向力)时的渗流长度较未考虑时的减少了16%左右， 速度大小降低了13%左右. 两种条件下， 流体渗流速度都受孔道弯曲程度控制， 且流体优先选择连通性较好的孔道进行渗流.

3) 应力敏感性指数主要受围压控制， 水压差对其影响甚微. 不同围压作用下应力敏感性指数呈“Z”字型变化. 其中， 竖向力300～400 Pa时， 应力敏感性指数降低最为显著， 降低幅度为82.5%. 应力敏感性指数增加， 土体孔道变窄的程度越大； 而应力敏感性指数减小， 土体中孤立孔隙与连通孔道相贯通的程度越大.