2019, Vol. 51引用本文 |
2. 中国科学院 水利部 成都山地灾害与环境研究所，四川 成都 610041;
3. 中国科学院大学，北京 100049;
4. 香港科技大学 土木与环境工程学系，香港
2. Inst. of Mountain Hazards and Environment, CAS, Chengdu 610044, China;
3. Univ. of Chinese Academy of Sciences, Beijing 100049, China;
4. Dept. of Civil and Environmental Eng., The Hong Kong Univ. of Sci. and Technol., Hong Kong, China
松散土体（有学者称为宽级配弱固结土）是发生滑坡、泥石流的典型土体类型之一[1]，尤其在地震影响区内该类土体分布更加广泛。2008年汶川地震后形成了大量的地震滑坡堆积体[2]，在降雨诱发条件下，坡面发生失稳进而转化为泥石流。例如：2010年，位于文家沟的滑坡松散堆积体在降雨作用下，形成大规模泥石流，造成了巨大灾害[3]。
已有研究表明，对于松散土坡体，其内部细颗粒在渗流作用下的迁移是诱发该类土体滑坡的主要原因之一[1]。Wang等[4]利用小型水槽进行不同细颗粒含量的降雨诱发滑坡试验，发现滑坡过程中细颗粒的存在将产生高孔隙水压力，细颗粒含量越高孔隙水压力越易升高，且在坡体滑动过程中仍保持此压力。国内学者通过室内与野外坡体降雨试验发现：细颗粒的存在不仅会引起内部孔隙水压力增高，同时其在土体内部的迁移会造成孔隙通道堵塞或贯通，将导致孔隙水压力随时间发生变化[5–6]；细颗粒总体沿渗流方向运移至坡体中部、坡脚及坡底处聚集[7–8]，由此而形成的高孔隙水压力面进一步发展为滑坡主滑移面[9]。孔隙通道堵塞和贯通与孔隙直径直接相关，当孔喉直径小于颗粒直径3倍或4倍时，拱桥效应将易导致孔隙被颗粒堵塞[10]。土体内部孔隙堵塞将造成局部孔隙水压力增高，土骨架有效应力降低，达到临界破坏面；若造成大面积孔隙堵塞将形成滞水面，为坡体内形成土体溜滑提供破坏面。当孔隙直径高于细颗粒直径5倍以上时，细颗粒群的运动造成孔隙堵塞的概率较低[10]，在坡体内渗流场作用下沿孔隙通道运移并最终流失，造成土体结构松散而易于坍塌。周小军等[11]通过人工降雨水槽试验定量分析降雨强度、斜坡坡度和降雨时间对松散土体内部细颗粒运移过程的影响程度，发现降雨历时与强度是影响细颗粒迁移积聚的主要因素。松散土内部所形成大孔隙为细颗粒迁移提供了充足的空间条件，降雨为细颗粒运移提供了动力条件，孔隙的堵塞和再贯通改变坡体内渗流场并导致不同土体破坏形式。以上研究对于认识松散土滑坡发生机制具有重要价值，但目前研究鲜见阐释此类土体内部细颗粒的分布特征、运移与孔隙堵塞过程，制约了进一步认识细颗粒迁移的机理和定量描述其迁移过程。
近年来，X射线断层扫描技术在医学、生物、材料及土壤等领域得到了广泛应用，此项技术能够在保证样品完整性的基础上快速获得样品3维微观结构特征[12]。利用上海同步辐射光源X射线成像及生物医学应用光束线站，对历经不同渗流时间的试样进行CT扫描以获得微观尺度下土体内部细颗粒分布、孔隙变化和侵蚀特征；结合基于离散元和Darcy流耦合模型的数值模拟，分析粗细颗粒相互作用机制和整体侵蚀过程；通过分析细颗粒动能随时间变化特征提出细颗粒堵塞孔隙和运移的机制，为深入理解降雨过程中松散土坡体破坏机制提供有力参考。
1 研究方法 1.1 试验设计研究表明：松散土体内部发生运移的细颗粒直径小于0.05 mm的含量达到85%[5]；CT技术所能分辨的颗粒直径与样品尺寸呈反比；根据现场测试，直径为10 mm左右砂土样品在9 μm成像分辨率平台上成像效果较好。利用筛分法和水分法难以精确分离直径小于0.05 mm左右的砂土颗粒并确定最小颗粒直径，故试验采用直径为0.075～0.100 mm的砂土颗粒作为细颗粒，样品直径为8.5 mm。依据微观土力学和微观流体力学研究尺度[13–14]，本研究中细颗粒运动过程中流体与颗粒及颗粒与颗粒相互作用的微观尺度范围为9～100 μm。采用不同粗细颗粒粒径之比表征孔隙特征对细颗粒运移过程的影响，粗颗粒粒径选取2种粒径区间（表1）：0.6～0.8 mm、0.8～0.9 mm。保证水力梯度一致，水流方向为自上而下。
| 表1 渗流试验设计方案 Tab. 1 Design scheme of seepage experiment |
 |
由于细小的样品直径不利于进行流体加压和铺设底部滤网，故设计制作可拆解铁质容器，通过在其中部增加橡胶弹性膜以增加壁厚，扩大滤网与渗流仪底部可接触面积，利用塑料软管将焊接于仪器顶部、直径为5 mm的导水管与供水槽连接。为防止流出颗粒在底部聚集而堵塞出流通道，将滤网底部设计为圆锥漏斗状，利用颗粒收集器皿获得试验过程中流出的细颗粒（图1（a））。
 |
| 图1 渗流试验仪器及CT扫描设施示意图 Fig. 1 Schematic diagram of seepage experimental instrument and CT scan equipments |
同时，3维CT技术的较长扫描周期导致无法将渗流试验和扫描过程同步进行，故采用拟静态研究方法，即当试验进行到0、4、8、12、16、20 h时终止（表1），将样品冻结运至扫描实验室，融化后再进行X射线扫描，试验设施如图1（b）所示。虽然每一个时间段内采用不同试样，颗粒分布存在差异，但能保证每个样品达到基本均匀混合，故初始状态下同一组试验的颗粒空间分布和孔隙特征具有一定相似性。
1.2 数值计算方法1971年，Cundall[15]为解决由破碎岩体造成的边坡稳定性问题，首次提出离散单元法（DEM）的概念，即一种能够模拟物体大变形、断裂等非连续介质力学问题的数值模拟方法。离散元颗粒流是一种目前发展较为完善的离散元方法，在颗粒流方法中，将单元体假设为一个3维球体或2维圆盘，颗粒与颗粒、颗粒与墙体之间遵循力与位移定律和牛顿第二定律[16]。
数值模拟是计算连续时间尺度上物质受力状态、速度及相互作用等动力学特征的有效手段，离散单元法与流体耦合计算是岩土工程领域近几十年发展起来的一种研究方法。相比于纳维–斯托克方程传统数值解法[17]、格子玻尔兹曼方法[18]及光滑粒子流[19]等流体力学计算方法，Darcy流模型具有计算简单、运算速度快的优势，结合颗粒流计算模型，能够有效地解决流体对颗粒的作用问题。因此，本研究将采用Darcy流与离散元进行耦合的方法开展计算。
根据Darcy流质量守恒方程，在饱和介质中的流体方程为[20]：
| $ \frac{\partial }{{\partial x}}\left({k_x}\frac{{\partial H}}{{\partial x}}\right) + \frac{\partial }{{\partial y}}\left({k_y}\frac{{\partial H}}{{\partial y}}\right) + \frac{\partial }{{\partial\textit{z}}}\left({k_{\textit{z}}}\frac{{\partial H}}{{\partial\textit{z}}}\right) +F = {C_{\rm{w}}}n{\gamma _{\rm{w}}}\frac{{\partial H}}{{\partial t}} $ | (1) |
式中：
水头边界条件：
| $ H = {H_{\rm{b}}}(t) $ | (2) |
流动边界条件：
| $ \left({k_x}\frac{{\partial H}}{{\partial x}}\right){l_x} + \left({k_y}\frac{{\partial H}}{{\partial y}}\right){l_y} + \left({k_{\textit{z}}}\frac{{\partial H}}{{\textit{z}}}\right){l_{\textit{z}}} + q(t) = 0 $ | (3) |
式中，
采用有限差分法求解方程（1），水头
| $\begin{aligned}[b] & \frac{\partial }{{\partial x}}\left({k_x}\frac{{\partial H}}{{\partial x}}\right) = {k_{x + \frac{{\Delta x}}{2}}}\frac{{H(x + \Delta x,y,{\textit{z}},t) - H(x,y,{\textit{z}},t)}}{{\Delta {x^2}}}- \\ & \quad\quad{k_{x - \frac{{\Delta x}}{2}}} \frac{{H(x,y,{\textit{z}},t) - H(x - \Delta x,y,{\textit{z}},t)}}{{\Delta {x^2}}} \frac{{\partial H}}{{\partial t}} =\\ & \quad\quad\frac{{H(x,y,{\textit{z}},t + \Delta t) - H(x,y,{\textit{z}},t)}}{{\Delta t}}\\[-15pt] \end{aligned}$ | (4) |
根据Darcy定律，通过渗透系数和水力梯度可计算每一个节点上的流速，采用形函数插值计算任意点的流速，根据此流速进一步获得流体拖拽力并施加到细颗粒上[21]。拖拽力大小与颗粒体积及流体速度有关，采用式（5）计算拖拽力
| $ {{{f}}_{\rm d}} = \frac{n}{{1 - n}}\beta {V_{\rm{p}}}{{{u}}_{\rm r}} $ | (5) |
式中：
| $ \beta=\left\{\begin{aligned} & {150 \mu \dfrac{(1-n)^{2}}{d_{\rm{p}}^{2} n^{2}}+1.75 \dfrac{(1-n) \rho_{0}\left|{{u}}_{\rm r}\right|}{d_{\rm{p}} n^{2}}}, { n<0.8}; \\ & {\dfrac{3}{4} C_{\rm{d}} \dfrac{|{{u}}_{\rm r}| \rho_{0}(1-n) n^{-2.7}}{d_{\rm{p}}}}, { n \ge 0.8}\end{aligned}\right. $ | (6) |
式中：
| $ C_{\mathrm{d}}=\left\{\begin{aligned} & {\dfrac{24\left(1+0.15 R e^{0.67}\right)}{R e}} , { Re<1\;000} ;\\ &{0.43}, { Re \ge 1\;000}\end{aligned}\right. $ | (7) |
因离散元计算中的时间步长与颗粒质量成正比，若取细颗粒直径为0.075～0.100 mm，将因颗粒质量过小导致计算效率低，故采用细颗粒粒径为0.5 mm，并将粗颗粒粒径按一定比例放大。模拟中所采用的球体颗粒相比于真实土体中不规则颗粒更有利于细颗粒在孔隙通道中运移，故粗细颗粒粒径之比需相对于真实土颗粒较小。综合考虑以上因素及不同粗细颗粒粒径之比的模拟计算情况，最终选取粗细颗粒粒径分别为3.0和0.5 mm、孔隙率为0.4、边长为20 mm的立方体模型，流体及离散元中材料的具体参数如表2所示。
| 表2 数值模拟参数值 Tab. 2 Parameters values of numerical simulation |
 |
2 结果与分析 2.1 CT扫描结果
将试样CT扫描结果进行重构获取颗粒分布的横截面图像，取未渗流和渗流时间20 h后的试样底部、中部和顶部图像进行分析，如图2～3所示。
 |
| 图2 粗颗粒粒径为0.8～0.9 mm时不同渗流时间CT结果对比 Fig. 2 CT results of coarse grain diameter 0.8～0.9 mm with different seepage time |
 |
| 图3 粗颗粒粒径为0.6～0.8 mm时不同渗流时间CT结果对比 Fig. 3 CT results of coarse grain diameter 0.6～0.8 mm with different seepage time |
在渗流20 h、粗粒直径为0.8～0.9 mm的图像中，顶部区域细颗粒被完全侵蚀（图2（a）和（d）），中部仍聚集大量细颗粒。试样底部虽然遭受一定程度的侵蚀作用，但其上部细颗粒的补充作用使其侵蚀程度相对较低（图2（c）和（f））。当粗颗粒直径为0.6～0.8 mm时，顶部区域内细颗粒虽被侵蚀但程度相对较低（图3（a）和（d）），但其底部位置侵蚀强烈，仅有局部区域存在少量细颗粒（图3（c）和（f））。
两种不同粗颗粒粒径样品的内部侵蚀特征存在明显差异。第1种存在顶部和底部优先侵蚀且上部颗粒对底部流失颗粒进行补充；第2种以从底部逐渐向上侵蚀为主，上部颗粒补充速度相对底部流失颗粒速度较小，主要因为相对较小的粗颗粒形成孔隙通道尺寸较小，同时通道内存在的大量细颗粒使颗粒运移空间进一步降低，从而被侵蚀的底部颗粒为上部细颗粒运移提供足够运移空间，在渗流力作用下上部颗粒随之发生运移。
2.2 细颗粒运移数值模拟结果沿
 |
| 图4 细颗粒迁移速度和侵蚀过程剖面示意图 Fig. 4 Schematic diagram of fine particle migration velocity and erosion process |
 |
| 图5 第15 s时沿渗流方向不同位置处横截面图 Fig. 5 Cross-sectional view at different locations in the seepage direction at 15 s |
在垂直于
 |
| 图6 不同位置处细颗粒数量和平均动能分布 Fig. 6 Distribution of fine particle number and average kinetic energy at different locations |
根据式（8）计算颗粒的平均动能
| $ {E_{\rm{a}}} = \frac{1}{{{n_{\rm{p}}}}}\sum\limits_{i = 1}^{{n_{\rm{p}}}} {\frac{1}{2}} {m_i}{v_i}^2 $ | (8) |
式中，
同时，截取此长方体区域内的区域①和②（图4）进一步分析细颗粒数量与平均动能随时间的变化特征（图（7）），将每一时刻的颗粒数量同区域内计算过程中颗粒最多的数量值相除，进行归一化处理，比值用
 |
| 图7 区域①和②内细颗粒数量和平均动能随时间变化统计分析 Fig. 7 Statistical analysis of the number of fine particles and the average kinetic energy over time in zone ① and zone ② |
细颗粒迁移的CT扫描和数值计算结果是人工降雨试验中松散土坡破坏微观机制的一种展现。顶部颗粒的快速流失与降雨过程中坡体顶部粗化现象相一致，底部颗粒的补充与流失平衡变化解释了坡脚处颗粒的聚集与流失问题。当坡体底部有大量的细颗粒快速流失时，虽有上部颗粒进行补充，但侵蚀程度通常大于补充速度，甚至出现较大粒径的颗粒流出，造成坡脚土体松散进而发生坍塌[25]。数值计算过程中孔隙堵塞存在长期与临时性两种现象，临时性堵塞导致局部孔隙水压力升高一定程度后，孔隙堵塞结构被打破而出现颗粒速度瞬间增高，在坡体内滞水面位置形成的高孔隙水压力即是冲坏表层土体的主要动力之一。
3 结论及展望本文通过CT扫描技术获得土体内部细颗粒运移特征，并结合数值模拟计算研究微观尺度动态流固耦合，阐明了试验中无法直接测量的颗粒位移、流速及颗粒间相互作用随时间的变化特征。进而结合这两方面研究，从微观尺度分析了坡体内部颗粒运移过程及其导致土体破坏的机理，得到如下主要结论：
1）Darcy流与离散元耦合模型对于计算土体内部侵蚀具有很好的效果，能与CT扫描结果和坡体人工降雨试验在侵蚀特征上得到互相验证。同时，数值模拟方法可较好地描述坡体内部细颗粒迁移及其导致孔隙堵塞与破坏的机制。
2）分析CT扫描和数值模拟计算结果，发现优先侵蚀发生在入流口和出流口，其驱动作用分别受流体压力和运移空间的影响，这与人工降雨试验中坡体顶部粗化和坡脚逐渐变得松散的现象一致，从而得到试验结果的印证。当孔隙尺寸较大时，流失的细颗粒能够获得其他颗粒补给；当孔隙尺寸较小时，流失的颗粒不易得到其他颗粒补给，表现出沿渗流方向空隙逐渐扩大的特征。
3）细颗粒在孔隙通道内会形成长期和临时堵塞的现象。通过颗粒的平均动能描述颗粒运移潜势，用以在数值模拟中描述颗粒位移特征并判断堵塞性质。长期堵塞区域内颗粒的平均动能随时间逐渐减小，并切断对其下方区域流失颗粒的补给；在临时堵塞区域，存在颗粒平均动能增高后颗粒数量继而增多的滞后效应。
本文针对尚未引起足够关注的松散土坡体内部细颗粒运移的微观过程进行研究，得到关于土体内部侵蚀的重要特征，发现CT技术与数值模拟方法对研究松散土破坏的微观机制具有重要作用。虽然Darcy流与离散元耦合模型能够在一定程度上反映土体内部侵蚀过程，但仍需进一步改进方法以体现真实情况下的颗粒不规则特征和微观扰流的影响。尽管CT扫描结果足以分辨出细颗粒整体侵蚀特征，但无法获取渗流过程中单个颗粒随时间的运动过程。同时，受水和气体影响，成像效果并不十分完美，需要进一步改进试验方法和图像处理技术。
[致谢]感谢上海同步辐射光源X射线成像及生物医学应用光束线站为研究提供试验平台，以及工作人员为试验顺利开展提供的指导和帮助。感谢长安大学博士生郭剑、中科院成都山地灾害与环境研究所博士生李尧对试验开展提供的帮助。
| [1] |
郭朝旭, 崔鹏. 宽级配弱固结土体内细颗粒迁移规律研究评述[J]. 山地学报, 2017, 35(2): 179-186. |
| [2] |
黄润秋. 汶川地震地质灾害后效应分析[J]. 工程地质学报, 2011, 19(2): 145-151. DOI:10.3969/j.issn.1004-9665.2011.02.001 |
| [3] |
刘传正. 汶川地震区文家沟泥石流成因模式分析[J]. 地质论评, 2012, 58(4): 709-716. DOI:10.3969/j.issn.0371-5736.2012.04.012 |
| [4] |
Wang G, Sassa K. Pore-pressure generation and movement of rainfall-induced landslides: Effects of grain size and fine-particle content[J]. Engineering Geology, 2003, 69(1/2): 109-125. |
| [5] |
Cui Peng,Guo Chaoxu,Zhou Jiawen,et al. The mechanisms behind shallow failures in slopes comprised of landslide deposits[J]. Engineering Geology, 2014, 180: 34-44. DOI:10.1016/j.enggeo.2014.04.009 |
| [6] |
陈晓清,崔鹏,冯自立,等. 滑坡转化泥石流起动的人工降雨试验研究[J]. 岩石力学与工程学报, 2006, 25(1): 106-116. DOI:10.3321/j.issn:1000-6915.2006.01.018 |
| [7] |
矫滨田,鲁晓兵,王淑云,等. 土体降雨滑坡中细颗粒运移及效应[J]. 地下空间与工程学报, 2005, 1(S1): 36-38. |
| [8] |
王志兵,李凯,汪稔,等. 细粒含量对泥石流斜坡失稳模式与规模的影响[J]. 水利水电科技进展, 2016(2): 35-41. DOI:10.3880/j.issn.1006-7647.2016.02.007 |
| [9] |
Cui Yifei,Zhou Xiaojun,Guo Chaoxu. Experimental study on the moving characteristics of fine grains in wide grading unconsolidated soil under heavy rainfall[J]. Journal of Mountain Science, 2017, 14(3): 417-431. DOI:10.1007/s11629-016-4303-x |
| [10] |
Valdes J R,Santamarina J C. Clogging: Bridge formation and vibration-based destabilization[J]. Canadian Geotechnical Journal, 2008, 45(2): 177-184. DOI:10.1139/T07-088 |
| [11] |
Zhou Xiaojun,Cui Peng,Jia Shitao,et al. Flume test study on the movement of fine grains based on orthogonal design [J].Journal of Sichuan University (Engineering Science edition),2012,44(Supp 1):83–88 周小军,崔鹏,贾世涛,等.基于正交设计的土体细颗粒迁移积聚水槽实验研究[J].四川大学学报(工程科学版),2012,44(增刊1): 83–88. |
| [12] |
周虎, 李文昭, 张中彬, 等. 利用X射线CT研究多尺度土壤结构[J]. 土壤学报, 2013, 50(6): 1226-1230. |
| [13] |
尹振宇. 土体微观力学解析模型: 进展及发展[J]. 岩土工程学报, 2013, 35(6): 993-1009. |
| [14] |
Zhang Kai.Research on the fluid flow and sample mixing in the microfluidic devices[D].Hangzhou:Zhejiang University, 2007. 张凯.微器件中流体的流动与混合研究[D].浙江:浙江大学, 2007. |
| [15] |
Cundall P A.A computer model for simulating progressive, large-scale movements in block rock systems[C]//Proceedings of International Symposium Fracture.Nancy,1971.
|
| [16] |
Potyondy D O,Cundall P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329-1364. DOI:10.1016/j.ijrmms.2004.09.011 |
| [17] |
Zou Yuhua,Chen Qun,Chen Xiaoqing,et al. Discrete numerical modeling of particle transport in granular filters[J]. Computers and Geotechnics, 2013, 47: 48-56. DOI:10.1016/j.compgeo.2012.06.002 |
| [18] |
Galindo-Torres S A,Scheuermann A,Mühlhaus H B,et al. A micro-mechanical approach for the study of contact erosion[J]. Acta Geotechnica, 2013, 10(3): 357-368. |
| [19] |
Wang Zhichao.Research on droplet impact to discrete particle based on a SPH-DEM coupling method [D].Tianjin: Tianjin University,2015. 王志超.基于SPH-DEM耦合方法的液滴冲击散粒体运动机理研究[D].天津: 天津大学, 2015. |
| [20] |
Tang Y,Chan D H,Zhu D Z. A coupled discrete element model for the simulation of soil and water flow through an orifice[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2017, 41(14): 1477-1493. DOI:10.1002/nag.v41.14 |
| [21] |
Cui Y F,Nouri A,Chan D,et al. A new approach to DEM simulation of sand production[J]. Journal of Petroleum Science and Engineering, 2016, 147: 56-67. DOI:10.1016/j.petrol.2016.05.007 |
| [22] |
Tsuji Y,Kawaguchi T,Tanaka T. Discrete particle simulation of two-dimensional fluidized bed[J]. Powder Technology, 1993, 77(1): 79-87. DOI:10.1016/0032-5910(93)85010-7 |
| [23] |
Wen C Y,Yu Y H. Mechanics of fluidization[J]. Chemical Engineering Progress Symposium Series, 1966, 62: 100-111. |
| [24] |
Ergun S. Fluid flow through packed columns[J]. Chemical Engineering Progress, 1952, 48(2): 6. |
| [25] |
Maeda K,Wood D M,Kondo A.Micro and macro modeling of internal erosion and scouring with fine particle dynamics[C]//Proceedings of 6th International Conference on Scour and Erosion.Paris:ISSMGE,2012:321–328.
|