高山峡谷区山洪沟口泥沙堆积特征数值试验

聂一品 兰玲 王协康

聂一品, 兰玲, 王协康. 高山峡谷区山洪沟口泥沙堆积特征数值试验 [J]. 工程科学与技术, 2024, 56(1): 237-244. doi: 10.15961/j.jsuese.202200524
引用本文: 聂一品, 兰玲, 王协康. 高山峡谷区山洪沟口泥沙堆积特征数值试验 [J]. 工程科学与技术, 2024, 56(1): 237-244. doi: 10.15961/j.jsuese.202200524
NIE Yipin, LAN Ling, WANG Xiekang. Numerical Experiment on Sediment Accumulation Characteristics at the Mouth of Flash Flood Gully in Alpine Canyon Area [J]. Advanced Engineering Sciences, 2024, 56(1): 237-244. doi: 10.15961/j.jsuese.202200524
Citation: NIE Yipin, LAN Ling, WANG Xiekang. Numerical Experiment on Sediment Accumulation Characteristics at the Mouth of Flash Flood Gully in Alpine Canyon Area [J]. Advanced Engineering Sciences, 2024, 56(1): 237-244. doi: 10.15961/j.jsuese.202200524

高山峡谷区山洪沟口泥沙堆积特征数值试验

基金项目: 国家自然科学基金重点项目(52239006;51639007)
详细信息
    • 收稿日期:  2022-05-26
    • 网络出版时间:  2023-11-24 12:00:00
  • 作者简介:

    聂一品(1998—),男,博士生. 研究方向:水力学及河流动力学. E-mail:nieyipin@163.com

    通信作者:

    王协康, 教授,E-mail:wangxiekang@scu.edu.cn

  • 中图分类号: TV142

Numerical Experiment on Sediment Accumulation Characteristics at the Mouth of Flash Flood Gully in Alpine Canyon Area

  • 摘要: 山洪泥石流是高山峡谷区一种常见的自然现象,其运动过程通常伴随剧烈的沟道侵蚀,并在沟口形成显著的泥沙堆积体。由于高山峡谷地区山高谷深,难以通过实地调查及模型实验等传统的研究手段追踪沟床泥沙受山洪泥石流冲刷的运动过程,无法揭示山洪泥石流挟带泥沙的堆积形态及堆积体中泥沙粒径分布特征,易产生对灾害影响范围估计不足导致的区域人员伤亡和财产损失。为探究山洪沟口泥沙堆积形态和粒度特征,以概化的高山峡谷区山洪沟为研究对象,采用耦合的计算流体力学(CFD)与离散单元法(DEM)数学模型,模拟了受到不同流变特性山洪泥石流冲刷的泥沙颗粒在沟口的堆积过程,重点分析了山洪泥石流体积浓度对堆积区中不同粒径泥沙颗粒堆积发展过程的影响。结果表明:泥沙堆积体的发展速度随着山洪泥石流体积浓度增加呈现先增加后减小的趋势,过高的体积浓度促使堆积体形态发生改变。泥沙平均堆积距离随时间变化可以分为高速增加、初次减速、增速恢复和稳定发展4个阶段。随着粒径的增大,泥沙与堆积体中心的平均距离也增大但堆积分散程度减小。泥沙颗粒的分散程度随时间变化过程与泥沙粒径和体积浓度密切相关,可以使用具有3个参数的幂函数对这个过程进行预测。泥沙颗粒粒径的增大加快其在堆积体中的分散速度而山洪泥石流体积浓度的增大使泥沙颗粒的分散速度先增加后减小。研究可为进一步理解山洪泥石流堆积致灾机理提供科学依据。

     

    Abstract: In alpine canyon, flash floods and debris flows are common. Usually, they are accompanied by strong gully erosion, and consequently create significant sediment accumulation in the gully mouth, which easily causes casualties and property losses in the region. It is difficult to track the movement process of gully bed sediment washed by flash floods and debris flows through field surveys, model experiments, and other traditional methods in alpine canyon areas due to high mountains and deep valleys. Therefore, it is impossible to provide information regarding the accumulation forms of sediment carried by mountain floods and debris flows as well as the distribution properties of sediment particle size in the accumulation body. To investigate the sediment accumulation morphology and particle size characteristics of mountain flood and debris flow, the generalized mountain gully area was taken as the research object. The coupled numerical methods of computational fluid dynamics (CFD) and discrete element method (DEM) were applied to simulation of sediment deposition process of mountain flood debris flow with different rheological characteristics. The effect of the volume concentration of mountain flood and debris flow on the accumulation process of sediment particles of different sizes in the accumulation area is analyzed in detail. As the volume concentration of mountain flood and debris flow increases, sediment accumulation rate first increases and then decreases. Furthermore, excessive volume concentration causes sediment to accumulate in a different form. The variations in sediment accumulation distance with time can be categorized into four stages: rapid increase, initial deceleration, recovery of growth and stable development. The variation of particle dispersion with movement time is of a power function, and the development rate is closely correlated with sediment particle size and volume concentration. There is a strong correlation between the degree of dispersion of sediment particles with time and the size and volume concentration of sediment particles. A power function with three parameters can be used to predict this process. Increased sediment particle size leads to an increase in the dispersion of sediment particles in the accumulation body. However, the increase in mountain floods and debris flows volume concentration causes the dispersion of sediment particles to increase at first and then decrease. The findings in this study can serve as a scientific basis for further understanding of mountain floods and debris flows.

     

  • 中国有近一半的地区被划分为山洪灾害防治区[1]。山洪泥石流作为一种典型的山区地质现象,不仅有着广泛的分布[2],且频繁成灾。在以西藏东部为典型代表的高山峡谷地区,常存在高差大、海拔高且狭窄的山洪泥石流沟道[3]。在地质活动和地区气候环境的影响下,高原冰川发育,局地性暴雨强盛,松散堆积物厚度大[4],山洪泥石流活动频繁。受国民经济发展和可居住条件影响,高山峡谷沟口附近仍有大量的人类活动。山洪泥石流冲出沟道发生堆积会严重威胁地区人民生命财产和重大工程安全[5]。因此,众多学者利用现场调查、物理实验以及数值模拟的手段对山洪泥石流堆积体的物质组成、发展过程和致灾机理进行研究。Cenderelli等[6]通过对历史山洪泥石流事件的调查,发现堆积体的前缘和表面的泥沙颗粒的粒径通常是堆积体中最大的。针对2019年汶川“8·20”山洪泥石流事件,Li等[7]通过野外调查并结合影像分析,指出不当的人类活动加剧了山洪泥石流损失。Zheng等[8]基于系列泥沙堆积实验,表明了山洪泥石流释放流量的增加会显著地增大堆积体厚度和长度比值。王协康等[9]以2010年“8·13”文家沟山洪泥石流试验表明,受主流挤压作用,堆积体易在沟口形成多元堆积结构。Chen等[10]基于有限元特征的分裂算法建立了山洪泥石流堆积与干流交汇的模型,模拟了文家沟山洪泥石流的堆积过程。Han等[11]引入HBP流体模型,采用SPH方法模拟了2010年日本Yohutagawa山洪泥石流运动,获得了与现场调查一致的堆积体延伸范围。

    综上,山洪泥石流堆积数值模拟大都是将多相流动视为单相流动[12],进而应用库伦塑性理论或黏塑性理论对流体流动进行描述[13],但忽略了不同大小颗粒流动过程的相互作用及其运动属性变化[14],难以揭示颗粒堆积过程致灾机理。离散单元法(DEM)作为一种能准确描述颗粒系统中颗粒之间的复杂的相互作用以及颗粒与周围环境之间的相互作用的数值模拟方法[15],能直接得到颗粒在运动过程中的动态信息[16]。计算流体力学方法(CFD)因其具有计算效率高、适用性好等优点,在地球表面的各种流体的计算中发挥着重要的作用[17]。为此,耦合CFD–DEM描述流体–粒子系统运动,已在研究颗粒系统的演进过程中成为了现实[18],例如滑坡坝破坏[19]以及山区河道泥沙输移[20]等。Li等[21]采用耦合CFD–DEM方法,研究了山洪泥石流对柔性防护网的影响,指出防护效率与山洪泥石流浆体弗劳德数密切相关。Fang等[22]结合CFD–DEM方法,评估了山洪泥石流对刚性防护结构的冲击特征,发现山洪泥石流固体体积分数的变化决定了防护结构上静荷载的分布。然而,耦合CFD–DEM方法对山洪泥石流的研究大都没有考虑其浆体流变特性的变化对颗粒运动的影响。为此,基于耦合CFD–DEM方法开展数值模拟,拟探究山洪泥石流体积浓度对不同粒径泥沙颗粒堆积过程中平均堆积距离和分散程度的差异,揭示了颗粒堆积分散程度随时间的变化规律,为理解山洪泥石流堆积致灾机理提供科学依据。

    采用CFDEM[23]耦合引擎耦合CFD开源软件OpenFOAM和DEM开源软件LIGGGHTS对高山峡谷区山洪沟口泥沙堆积形态及粒度特征进行计算。计算过程中使用的控制方程如下:

    1)颗粒运动控制方程。颗粒运动主要考虑平动和转动。运动过程颗粒可能与邻近的颗粒和边壁发生碰撞或者与周围的流体产生相互作用进而实现动量和能量交换。因此,单个颗粒的运动可以通过牛顿第二定律进行描述,其控制方程可以写成:

    $$ {m_i}\frac{{{\text{d}}{{\boldsymbol{v}}_i}}}{{{\text{d}}t}} = \sum\limits_j {{\boldsymbol{F}}_{i,j}^{\text{c}}} + {\boldsymbol{F}}_i^{\text{f}} + {\boldsymbol{F}}_i^{\text{g}} $$ (1)
    $$ {I_i}\frac{{{\text{d}}{{\boldsymbol{\omega }}_i}}}{{{\text{d}}t}} = \sum\limits_j {{{\boldsymbol{M}}_{i,j}}} $$ (2)

    式(1)~(2)中,${m_i}$为颗粒$i$的质量,$ {{\boldsymbol{v}}_i} $为颗粒$i$的平动速度,$ {I_i} $为颗粒的转动惯量,$ {{\boldsymbol{\omega }}_i} $为颗粒$i$的角速度,$ {\boldsymbol{F}}_{i,j}^{\text{c}} $为其它颗粒或边壁作用于颗粒$i$上的接触力,$ {\boldsymbol{F}}_i^{\text{f}} $为颗粒$i$通过颗粒–流体相互作用受到的作用力,$ {\boldsymbol{F}}_i^{\text{g}} $为颗粒$i$的重力,$ {{\boldsymbol{M}}_{i,j}} $为其它颗粒或边壁作用于颗粒$i$上的力矩。

    2)流体运动控制方程。黏性不可压缩流体运动通过连续性方程和动量方程进行描述,其形式如下:

    $$ \frac{{\partial {\alpha _{\text{f}}}}}{{\partial t}} + \nabla \cdot \left( {{\alpha _{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}} \right) = 0 $$ (3)
    $$ \begin{split} & \frac{\partial }{{\partial t}}\left( {{\alpha _{\text{f}}}{\rho _{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}} \right) + \nabla \cdot \left( {{\alpha _{\text{f}}}{\rho _{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}} \right) =\\&\qquad - {\alpha _{\text{f}}}\nabla p + \nabla\cdot \left( {{\alpha _{\text{f}}}{\tau _{\text{f}}}} \right) + {\alpha _{\text{f}}}{\rho _{\text{f}}}g - {{\boldsymbol{f}}_{{\text{fluid}}}} \end{split}$$ (4)

    式中,$ {\alpha _{\text{f}}} $为流体的体积分数,$ {{\boldsymbol{u}}_{\text{f}}} $为流体的运动速度,$ {\rho _{\text{f}}} $为流体的密度,$ p $为流体的压强,$ {\tau _{\text{f}}} $为流体的应力张量,$g$为重力加速度,$ {{\boldsymbol{f}}_{{\text{fluid}}}} $为流体平均下的颗粒通过颗粒–流体相互作用受到的阻力。

    3)颗粒–流体间相互作用力控制方程。颗粒–流体相互作用力是由颗粒和流体之间相互运动而产生的黏性阻力,其大小与流体的流动状态和颗粒雷诺数有关。表达式为:

    $$ {{\boldsymbol{f}}_{{\text{fluid}}}}{\text{ = }}\frac{1}{2}{C_{\mathrm{d}}}{\rho _{\mathrm{f}}}{\text{π}} {r^2}\left| {{{\boldsymbol{u}}_{\mathrm{f}}} - {{\boldsymbol{v}}_i}} \right|\left( {{{\boldsymbol{u}}_{\mathrm{f}}} - {{\boldsymbol{v}}_i}} \right) \cdot {n^{1 - \chi }} $$ (5)

    式中:$ {C_{\mathrm{d}}} $为阻力系数;$\chi $为经验修正因子,$\chi {\text{ = }}3.7 - 0.65\exp ( {{{{{ - }}{{\left( {1.5 - \lg R{e_{\text{p}}}} \right)}^2}} / 2}} )$,其中,$ R{e_{\text{p}}} $为颗粒雷诺数。

    1)模型设置

    采用的概化高山峡谷区山洪泥石流沟道数值模型如图1所示。通过控制水箱中流体的高度,使不同体积浓度山洪泥石流在挡板下的出口处具有相同的压强。将粒径为1至2 mm的泥沙颗粒视为泥石流浆体的组成部分,使用单相流体模拟泥石流浆体的运动并结合流体体积分数(VOF)方法捕捉泥石流与空气的界面。挡板处于x = 0 的平面上。流通区底部铺有厚度为0.1 m的卵砾石层,其粒径为6至24 mm。流通区和加速区比降为10%,长度分别为3 m和1 m。堆积区采用1%的比降,大小为6 m×6 m。模拟的总时间为20 s。

    图  1  典型山洪泥石流沟道数值试验模型
    Fig.  1  Numerical test model of a typical landslide gully
    下载: 全尺寸图片

    2)计算参数设置

    已有研究表明,山洪泥石流浆体的流变特性可以使用牛顿流体、宾汉流体、幂律流体以及二项式模式进行表述[24],而宾汉流体因其形式简单且具有较高的精度已经在山洪泥石流3维模拟中有着广泛的应用[11]。此外,已有大量研究对宾汉流体模型的流变参数进行了研究[12],故在本数值试验中使用宾汉流体模型对山洪泥石流的流变行为进行表达。试验中采用4种典型的山洪泥石流体积浓度(${C_{\mathrm{v}}}$)分别为15%、30%、45%和60%。在屈服应力的计算中,当体积浓度较低时采用文献[24]中的Kang和Zhang公式($0.15 < {C_{\mathrm{v}}} < 0.5$);而当体积浓度较高时采用文献[25]中的Major和Pierson公式($0.44 < {C_{\mathrm{v}}} < 0.66$)。刚度系数通过计算不同体积浓度山洪泥石流的相对黏度,再将其乘以常温下水的动力黏度得到。相对黏度采用文献[26]中的Thomas经验公式计算。为了解决宾汉流体模型存在的低剪切情况下刚度系数数值发散问题[27],采用双黏度正则化模型,在低剪切状态下引入等效黏性系数解决该问题[28]。经过计算,4种典型山洪泥石流浆体的模拟参数见表1

    表  1  山洪泥石流浆体模拟参数
    Table  1  Simulation parameters of debris flow
    体积
    浓度/%
    密度/
    (kg·m–3)
    屈服应
    力/Pa
    刚度系数/
    (Pa·s)
    等效黏度
    系数/(Pa·s)
    15 1248 0.566 1.63×10–3 1.634
    30 1495 1.83 3.05×10–3 3.052
    45 1743 5.91 8.95×10–3 8.950
    60 1990 44.4 6.39×10–2 63.89

    耦合模型中颗粒计算参数见表2。其中,颗粒之间的计算参数与颗粒和壁面之间的计算参数保持一致。

    表  2  耦合模型中颗粒计算参数
    Table  2  Particle calculation parameters
    颗粒参数 数值
    密度/(kg·m–3) 2 650
    泊松比 0.3
    杨氏模量/(N·m–1) 5×106
    碰撞恢复系数 0.3
    滑动摩擦系数 0.4
    滚动摩擦系数 0.1

    使用泥浆溃坝实验来检验耦合模型对两相流体的交界面捕捉能力以验证模型的适用性。图2(a)绘制了Komatina[29]实验中第0.6 s时的泥浆表面和模拟结果。由图2(a)可见,模型对两相流体交界面的捕捉是较为准确的。随后计算颗粒在宾汉流体中沉降时阻力系数–雷诺数的关系来验证模型的准确性。经过模拟的阻力系数–雷诺数关系如图2(b)所示,其中理论结果选用Cheng[30]提出的公式。由图2(b)可见,模型能准确计算颗粒运动过程中颗粒–流体相互作用力。因此,耦合CFD–DEM模型可以准确地计算宾汉流体中颗粒的运动情况。

    图  2  CFD–DEM耦合模型验证
    Fig.  2  CFD–DEM coupling model validation results
    下载: 全尺寸图片

    计算开始后,床面泥沙颗粒在山洪泥石流的冲刷下发生运动。当泥沙颗粒运动到堆积区后,由于山洪泥石流速度降低,泥沙颗粒开始发生沉积。泥沙颗粒在堆积区中不同时刻的堆积的形态如图3所示,图3中,体积分数表明该时刻堆积区上山洪泥石流的运动痕迹。从图3(a)(c)可以看出:体积浓度较低时(${C_{\mathrm{v}}} \leqslant 45\% $),山洪泥石流冲刷的颗粒在初始堆积阶段(t=5.0 s),由于堆积平面足够大且没有阻挡物,泥沙呈现扇形堆积;随着山洪泥石流将更多的颗粒携带到堆积区中(t=7.5 s),堆积体堆积范围不断扩大,但依然呈现扇形分布。当${C_{\mathrm{v}}} \leqslant $45%时,山洪泥石流对流通区床面的冲刷随着体积浓度的增加而增加[31],因此,在同一时刻下泥沙颗粒的堆积范围随着${C_{\mathrm{v}}}$的增加不断增加。在t=10 s时,由于$ {C_{\mathrm{v}}}= $30%和45%的模型中的床面已经被侵蚀完毕,因此在堆积区入口处,泥沙不断受到后续山洪泥石流的推动而产生空隙。在图3(d)中,沟床颗粒在${C_{\mathrm{v}}} = $60%山洪泥石流冲刷下形成的堆积体形态在堆积区与流通区交界处有一段很短的梯形堆积区域,随后扩展成半圆状。此时,由于山洪泥石流体积浓度很高,在运动过程中需要更多地克服黏性阻力,因此,堆积体发展的速度有所减缓。

    图  3  泥沙颗粒堆积过程模拟结果
    Fig.  3  Simulation results of sediment particle accumulation process
    下载: 全尺寸图片

    为了探究颗粒的粒径对其堆积行为的影响,分别计算了4种体积浓度下不同粒径颗粒与堆积体几何中心之间的距离($\Delta r$)随时间变化情况,计算结果如图4所示。从图4中可以看出,各粒径泥沙颗粒$\Delta r$随时间变化可分为4个阶段,分别为高速增加阶段、初次减速阶段、增速恢复阶段和稳定发展阶段。

    图  4  不同粒径颗粒距堆积体几何中心的平均距离
    Fig.  4  Average distance between particles of different sizes and geometric center of accumulation
    下载: 全尺寸图片

    第1阶段为高速增加阶段,该阶段中山洪泥石流从水箱中流出后在床面产生沿程冲刷,在流通区中产生了“龙头”,大量泥沙颗粒在聚集在山洪泥石流的前缘。当山洪泥石流到达堆积区时,“龙头”中的颗粒以很高的速度进入,导致$\Delta r$的高速增加。这一阶段中,山洪泥石流的冲刷使粒径较大泥沙颗粒运动速度显著快于较小的泥沙颗粒,因而当其进入堆积区后,大粒径泥沙的$\Delta r$增速大于粒径较小的泥沙。第2个阶段为初次减速阶段。在山洪泥石流的持续冲刷下,流通区床面颗粒不断受到侵蚀。但紧随“龙头”的后续浆体中泥沙含量相对较低,导致了$\Delta r$的增速减缓。

    当流通区床面泥沙被山洪泥石流完全侵蚀后,处于堆积区中的泥沙颗粒在后续浆体的推动下产生整体性运动,$\Delta r$随时间的变化进入到第3阶段即增速恢复阶段。这一阶段中,更大粒径泥沙颗粒的$\Delta r$增速依然大于较小的泥沙。此时,山洪泥石流浆体受到外围颗粒的阻挡,后续的浆体尚未从堆积体中流出。但随着模拟的进行,在某一方向上颗粒的堆积宽度减小,导致浆体从堆积体中流出,对堆积体的整体推动作用减小,形成了稳定发展的第4个阶段。${C_{\mathrm{v}}}= $15%时,稳定发展阶段的堆积体中颗粒的$\Delta r$最终不变,但随着${C_{\mathrm{v}}}$的增加,$\Delta r$不再保持不变,其增速显著增大。在该阶段中,$\Delta r$的增速随着粒径的减小逐渐增加,表明了在沟床的颗粒完全进入到堆积平面中后,由于山洪泥石流的持续作用,各个粒径的颗粒逐渐地混合在一起。Cv=60%的山洪泥石流的流动速度较低,到模拟结束时,沟床泥沙颗粒没有被完全侵蚀,因此,模拟中$\Delta r$随时间的变化仅有高速增加阶段和初次减速阶段。当${C_{\mathrm{v}}} \leqslant $45%,随着${C_{\mathrm{v}}}$的增加,相同粒径的泥沙颗粒在同一时刻下的$\Delta r$逐渐增加,堆积体对堆积区的侵占更加严重。

    模拟计算了大量(共计30577个)泥沙颗粒的运动数据,因此从统计的角度分析堆积区内颗粒运动的内在规律是可行的。将每一组模拟中颗粒的$\Delta r$按照绘制箱线图的方法进行统计。如图5所示,dQ1dQ3分别代表经过统计得到的$\Delta r$的下和上四分位数。dIQR为四分位距,等于dQ1dQ3的差值,反映了统计意义上$\Delta r$的集中趋势,可以表征颗粒在堆积体中分布的分散程度。上边界使用dQ3+1.5×dIQR表示,代表了非正态分布的$\Delta r$异常值的判别标准,用于表示堆积体外边界范围。经过统计分析,模拟过程中任何时刻异常$\Delta r$的数量都不多于8.5%,因此使用统计的方法得到的dQ3dIQRdQ3+1.5×dIQR来分别判断山洪泥石流堆积体中粒度特征和堆积范围是可以接受的。

    图  5  箱线图组成部分
    Fig.  5  Components of box plot
    下载: 全尺寸图片

    图6为堆积体中泥沙颗粒分散程度(dIQR)与泥沙颗粒的上四分位数(dQ3)的关系随粒径变化趋势,其中,dIQR/dQ3是将模拟过程中统计得到的dIQRdQ3的数值相除后进行平均。

    图  6  颗粒堆积体中dIQRdQ3之间的关系
    Fig.  6  Relationship between dIQR and dQ3
    下载: 全尺寸图片

    图6可知,当山洪泥石流体积浓度较低(${C_{\mathrm{v}}} \leqslant 45{\text{%}} $)时,dIQR/dQ3的值会随泥沙颗粒粒径的增加而不断的减小。这一现象说明,随着粒径的增大,泥沙颗粒堆积更加集中。此外,不同的dIQR/dQ3代表了各粒径颗粒在堆积平面上所处位置的不同。由图3(a)(c)可知:在${C_{\mathrm{v}}} \leqslant 45{\text{%}} $的山洪泥石流冲刷形成的泥沙堆积体中,较大的泥沙颗粒显著地位于堆积体的外侧且各粒径颗粒在堆积体中的分选良好;${C_{\mathrm{v}}} = 60{\text{%}} $时,堆积体中不同粒径颗粒的dIQR/dQ3不随粒径发生显著的变化,大致保持在0.45附近。这个现象表明,随着山洪泥石流体积浓度增加,各种粒径的颗粒较为均匀的混杂在堆积体当中。结合图3(d)可发现,即使当t=10 s时,堆积体的中心仍然可以明显看见大粒径泥沙的存在,泥沙颗粒之间的分选性较差。除粒径为6 mm的颗粒外,当${C_{\mathrm{v}}} \leqslant 45{\text{%}} $时,dIQR/dQ3${C_{\mathrm{v}}}$的增加逐渐减小。这个现象表明随着山洪泥石流体积浓度的增加,泥沙颗粒的分选性逐渐变差。综合以上分析,堆积体中dIQR/dQ3的比值可以较为准确地代表各组分颗粒在堆积体内的相对位置,同时不同颗粒dIQR/dQ3数值的变化也表明了堆积体中颗粒之间具有良好的分选性。

    图7${C_{\mathrm{v}}} = 15\% $山洪泥石流冲刷下,粒径为6 mm的颗粒在堆积过程中dIQR随时间的变化过程。由图7可知,随着时间的增加,dIQR的增速逐渐减小,且有接近一个最大值的趋势。由于被冲刷的泥沙颗粒需要一定的时间才能进入堆积区之中,故假设dIQR随时间的增长满足形如$ a \times {t^b} + c $的幂函数关系。根据模拟数据点的分布情况,$a$$b$都小于0,而$c$大于0。图7的预测带代表了dIQR在95%置信度上可能的范围,且几乎所有的模拟数据点都在95%预测带以内,这代表了幂函数关系对dIQR随时间的变化过程较为精准。结合图6的分析结果,可以在任意时刻通过dIQR的拟合结果,推断出dQ3的数值,从而估算山洪泥石流堆积体影响范围时间的变化过程。

    图  7  泥沙颗粒分散程度(dIQR)随时间变化过程
    Fig.  7  Variation in the dispersion of sediment particles (dIQR) over time
    下载: 全尺寸图片

    所有模拟泥沙颗粒分散程度(dIQR)随时间变化过程的拟合参数与体积浓度(${C_{\mathrm{v}}}$)和粒径($d$)之间的关系如图8所示。

    图  8  拟合指标ifit与碰撞指标icollision的关系
    Fig.  8  Relationship between the ifit and icollision
    下载: 全尺寸图片

    图8中箭头的方向指示了不同模拟中同一粒径颗粒的拟合指标(ifit,即$a \times {b / c}$)变化趋势。${C_{\mathrm{v}}}$的增加代表了山洪泥石流浆体内颗粒碰撞的几率的增加[32],而离散的颗粒与浆体中颗粒碰撞几率会随着颗粒粒径$d$增大而增大,故碰撞指标(icollision,即${C_{\mathrm{v}}} \times d$)能代表离散的颗粒通过颗粒碰撞获得能量的几率大小。而$a \times {b / c}$的增加则代表了颗粒的dIQR极值的减小和其增速的增大,也即代表着颗粒达到极限分散程度所需时间也较小。

    粒径越大的泥沙颗粒可以通过碰撞能够获得更多的能量,其进入到堆积区后能在更短时间内克服较多摩擦阻力,实现更快的分散。随着${C_{\mathrm{v}}}$的逐渐减小,粒径较小的颗粒与浆体发生能量交换的几率大大减小,因此ifit$d$的增益也越大。而当颗粒完成一定程度的分散之后,堆积体的存在导致了浆体的运动方向发生偏转,处于堆积体前缘的大颗粒难以被浆体推动,导致其分散状态难以继续改变,达到了极限分散状态。对相同粒径的泥沙颗粒来说,随着${C_{\mathrm{v}}}$的增加,浆体在运动过程中因其黏性消耗的能量不断增加,抵消了一部分传递给离散颗粒的能量,抑制了颗粒的分散。因此,如图8所示,同一${C_{\mathrm{v}}}$下,随着$d$的增大,ifit也随之增加,且${C_{\mathrm{v}}}$越小,$a \times {b \mathord{\left/ {\vphantom {b c}} \right. } c}$$d$的增益也越大。而在同一粒径下,$a \times {b \mathord{\left/ {\vphantom {b c}} \right. } c}$会随着${C_{\mathrm{v}}}$的增加而呈现出先增加后减小的趋势。

    基于高山峡谷区山洪泥石流沟道的特点,借助耦合CFD–DEM方法对沟口泥沙颗粒的运动进行了数值计算,重点研究了堆积体整体的形态及不同粒径颗粒的分布特征随时间的变化,主要结论如下:

    1)随着山洪泥石流体积浓度的增加,泥沙堆积体的发展速度呈现先增加后减小的趋势,当体积浓度较高时堆积体的形态也随之发生改变。

    2)不同粒径泥沙颗粒与堆积体几何中心之间的平均堆积距离($\Delta r$)随着模拟时间的变化呈现出4个不同的阶段。粒径的增大会导致堆积距离的增速增大,但随着运动的发展,各粒径颗粒堆积距离之间的差距会逐渐减小。

    3)通过统计得到颗粒堆积距离的四分位距(dIQR)和上边界(dQ3+1.5×dIQR)可以表示堆积体中颗粒的分散程度以及堆积范围。不同粒径颗粒在堆积体中dIQR/dQ3的变化则反映了堆积体中分选性的良好与否。dIQR随时间变化的趋势的符合幂函数关系,通过拟合结果与dIQR/dQ3的数值变化规律可以得到任意时刻的泥石流堆积体的堆积范围。

    4)相同体积浓度下,随着颗粒粒径的增加,颗粒的极限分散程度减小但能更快地分散。在颗粒粒径相同时,随着体积浓度的增加,颗粒的极限分散程度先减小后增加而分散速度先增加后减小。

  • 图  1   典型山洪泥石流沟道数值试验模型

    Fig.  1   Numerical test model of a typical landslide gully

    下载: 全尺寸图片

    图  2   CFD–DEM耦合模型验证

    Fig.  2   CFD–DEM coupling model validation results

    下载: 全尺寸图片

    图  3   泥沙颗粒堆积过程模拟结果

    Fig.  3   Simulation results of sediment particle accumulation process

    下载: 全尺寸图片

    图  4   不同粒径颗粒距堆积体几何中心的平均距离

    Fig.  4   Average distance between particles of different sizes and geometric center of accumulation

    下载: 全尺寸图片

    图  5   箱线图组成部分

    Fig.  5   Components of box plot

    下载: 全尺寸图片

    图  6   颗粒堆积体中dIQRdQ3之间的关系

    Fig.  6   Relationship between dIQR and dQ3

    下载: 全尺寸图片

    图  7   泥沙颗粒分散程度(dIQR)随时间变化过程

    Fig.  7   Variation in the dispersion of sediment particles (dIQR) over time

    下载: 全尺寸图片

    图  8   拟合指标ifit与碰撞指标icollision的关系

    Fig.  8   Relationship between the ifit and icollision

    下载: 全尺寸图片

    表  1   山洪泥石流浆体模拟参数

    Table  1   Simulation parameters of debris flow

    体积
    浓度/%
    密度/
    (kg·m–3)
    屈服应
    力/Pa
    刚度系数/
    (Pa·s)
    等效黏度
    系数/(Pa·s)
    15 1248 0.566 1.63×10–3 1.634
    30 1495 1.83 3.05×10–3 3.052
    45 1743 5.91 8.95×10–3 8.950
    60 1990 44.4 6.39×10–2 63.89

    表  2   耦合模型中颗粒计算参数

    Table  2   Particle calculation parameters

    颗粒参数 数值
    密度/(kg·m–3) 2 650
    泊松比 0.3
    杨氏模量/(N·m–1) 5×106
    碰撞恢复系数 0.3
    滑动摩擦系数 0.4
    滚动摩擦系数 0.1
  • [1] 王协康,刘兴年,周家文.泥沙补给突变下的山洪灾害研究构想和成果展望[J].工程科学与技术,2019,51(4):1–10. doi: 10.15961/j.jsuese.201900261

    Wang Xiekang,Liu Xingnian,Zhou Jiawen.Research framework and anticipated results of flash flood disasters under the mutation of sediment supply[J].Advanced Engineering Sciences,2019,51(4):1–10 doi: 10.15961/j.jsuese.201900261
    [2] 费祥俊,舒安平.泥石流运动机理与灾害防治[M].北京:清华大学出版社,2004.
    [3] 陈龙,余斌,黄海,等.西藏天摩沟泥石流物源演变特征[J].地质通报,2021,40(12):2089–2097. doi: 10.16030/j.cnki.issn.1000-3665.2019.05.19

    Chen Long,Yu Bin,Huang Hai,et al.Characteristics of provenance evolution of debris flow in Tianmogou,Tibet[J].Geological Bulletin of China,2021,40(12):2089–2097 doi: 10.16030/j.cnki.issn.1000-3665.2019.05.19
    [4] 吴积善,程尊兰,耿学勇.西藏东南部泥石流堵塞坝的形成机理[J].山地学报,2005,23(4):399–405. doi: 10.3969/j.issn.1008-2786.2005.04.003

    Wu Jishan,Cheng Zunlan,Geng Xueyong.Formation of dam from debris flow in the southeast Tibet[J].Journal of Mountain Science,2005,23(4):399–405 doi: 10.3969/j.issn.1008-2786.2005.04.003
    [5] 崔鹏,邹强.山洪泥石流风险评估与风险管理理论与方法[J].地理科学进展,2016,35(2):137–147. doi: 10.18306/dlkxjz.2016.02.001

    Cui Peng,Zou Qiang.Theory and method of risk assessment and risk management of debris flows and flash floods[J].Progress in Geography,2016,35(2):137–147 doi: 10.18306/dlkxjz.2016.02.001
    [6] Cenderelli D A,Kite J S.Geomorphic effects of large debris flows on channel morphology at North Fork Mountain,eastern West Virginia,USA[J].Earth Surface Processes and Landforms,1998,23(1):1–19. doi: 10.1002/(SICI)1096-9837(199801)23:1<1::AID-ESP814>3.0.CO;2-3
    [7] Li Yu,Liu Xingnian,Gan Binrui,et al.Formation-evolutionary mechanism analysis and impacts of human activities on the 20 August 2019 clustered debris flows event in Wenchuan County,southwestern China[J].Frontiers in Earth Science,2021,9:616113. doi: 10.3389/feart.2021.616113
    [8] Zheng Hongchao,Shi Zhenming,Hanley K J,et al.Deposition characteristics of debris flows in a lateral flume considering upstream entrainment[J].Geomorphology,2021,394:107960. doi: 10.1016/j.geomorph.2021.107960
    [9] 王协康,叶龙,王玉林,等.陡坡沟道强输沙汇口堆积形态演变试验研究[J].四川大学学报(工程科学版),2014,46(5):7–13. doi: 10.15961/j.jsuese.2014.05.032

    Wang Xiekang,Ye Long,Wang Yulin,et al.Experimental study on the evolution of the accumulation body in the confluence due to steep gully of the high intensity sediment transport[J].Journal of Sichuan University(Engineering Science,2014,46(5):7–13 doi: 10.15961/j.jsuese.2014.05.032
    [10] Chen Ridong,Liu Xingnian,Cao Shuyou,et al.Numerical simulation of deposit in confluence zone of debris flow and mainstream[J].Science China Technological Sciences,2011,54(10):2618–2628. doi: 10.1007/s11431-011-4510-1
    [11] Han Zheng,Su Bin,Li Yange,et al.Numerical simulation of debris-flow behavior based on the SPH method incorporating the Herschel-Bulkley-Papanastasiou rheology model[J].Engineering Geology,2019,255:26–36. doi: 10.1016/j.enggeo.2019.04.013
    [12] 崔鹏,唐金波,林鹏智.泥石流运动阻力特性及其研究进展[J].四川大学学报(工程科学版),2016,48(3):1–11. doi: 10.15961/j.jsuese.2016.03.001

    Cui Peng,Tang Jinbo,Lin Pengzhi.Research progress of resistance character of debris-flow[J].Journal of Sichuan University(Engineering Science Edition),2016,48(3):1–11 doi: 10.15961/j.jsuese.2016.03.001
    [13] Ancey C.Plasticity and geophysical flows:A review[J].Journal of Non-Newtonian Fluid Mechanics,2007,142(1/2/3):4–35. doi: 10.1016/j.jnnfm.2006.05.005
    [14] Chen H,Lee C F.A dynamic model for rainfall-induced landslides on natural slopes[J].Geomorphology,2003,51(4):269–288. doi: 10.1016/S0169-555X(02)00224-6
    [15] Zhu H P,Zhou Z Y,Yang R Y,et al.Discrete particle simulation of particulate systems:A review of major applications and findings[J].Chemical Engineering Science,2008,63(23):5728–5770. doi: 10.1016/j.ces.2008.08.006
    [16] Lei Ming,Yang Po,Wang Yikui,et al.Numerical analyses of the influence of baffles on the dynamics of debris flow in a gully[J].Arabian Journal of Geosciences,2020,13(19):1052. doi: 10.1007/s12517-020-06016-z
    [17] von Boetticher A,Turowski J M,McArdell B W,et al.DebrisInterMixing-2.3:Afinite volume solver for three-dimensional debris–flow simulations with two calibration parameters–Part 1:Model description[J].Geoscientific Model Development,2016,9(9):2909–2923. doi: 10.5194/gmd-9-2909-2016
    [18] Zhao T,Utili S,Crosta G B.Rockslide and impulse wave modelling in the vajont reservoir by DEM–CFD analyses[J].Rock Mechanics and Rock Engineering,2016,49(6):2437–2456. doi: 10.1007/s00603-015-0731-0
    [19] 雷明,余海逖,许泽星,等.山区变比降河道卵砾石溯源淤积的数值模拟[J].工程科学与技术,2019,51(1):45–51. doi: 10.15961/j.jsuese.201700981

    Lei Ming,Yu Haiti,Xu Zexing,et al.Numerical simulation of retrogressive pebble deposition in the changing slope zone of a mountainous river[J].Advanced Engineering Sciences,2019,51(1):45–51 doi: 10.15961/j.jsuese.201700981
    [20] Nie Yipin,Lan Ling,Yan Xufeng,et al.Effect of gradation variation on particle transport process in a generalized flash flood gully via CFD–DEM method[J].Acta Geophysica,2023,71(1):391–404. doi: 10.1007/s11600-022-00862-z
    [21] Li Xingyue,Zhao Jidong,Kwan J S H.Assessing debris flow impact on flexible ring net barrier:A coupled CFD–DEM study[J].Computers and Geotechnics,2020,128:103850. doi: 10.1016/j.compgeo.2020.103850
    [22] Fang Jun,Wang Lizhong,Hong Yi,et al.Influence of solid–fluid interaction on impact dynamics against rigid barrier:CFD–DEM modelling[J].Géotechnique,2022,72(5):391–406. doi: 10.1680/jgeot.19.p.160
    [23] Goniva C,Kloss C,Deen N G,et al.Influence of rolling friction on single spout fluidized bed simulation[J].Particuology,2012,10(5):582–591.
    [24] 王裕宜,詹钱登,严璧玉.泥石流体的流变特性与运移特征[M].长沙:湖南科学技术出版社,2014.
    [25] Major J J,Pierson T C.Debris flow rheology:Experimental analysis of fine-grained slurries[J].Water Resources Research,1992,28(3):841–857. doi: 10.1029/91wr02834
    [26] Thomas D G.Transport characteristics of suspension:Ⅷ.A note on the viscosity of Newtonian suspensions of uniform spherical particles[J].Journal of Colloid Science,1965,20(3):267–277. doi: 10.1016/0095-8522(65)90016-4
    [27] Papanastasiou T C.Flows of materials with yield[J].Journal of Rheology,1987,31(5):385–404. doi: 10.1122/1.549926
    [28] Hojeij A,Jossic L,Séchet P,et al.Experimental study and numerical modeling of mixing by air injection in yield stress fluids using the OpenFOAM software[J].AIChE Journal,2022,68(2):17442–1-17442-19. doi: 10.1002/aic.17442
    [29] Komatina D,Jovanovic M.Experimental study of steady and unsteady free surface flows with water-clay mixtures[J].Journal of Hydraulic Research,1997,35(5):579–590. doi: 10.1080/00221689709498395
    [30] Cheng Niansheng.Comparison of formulas for drag coefficient and settling velocity of spherical particles[J].Powder Technology,2009,189(3):395–398. doi: 10.1016/j.powtec.2008.07.006
    [31] 王兆印,张新玉.水流冲刷沉积物生成泥石流的条件及运动规律的试验研究[J].地理学报,1989,44(3):291–301. doi: 10.3321/j.issn:0375-5444.1989.03.006

    Wang Zhaoyin,Zhang Xinyu.Experimental study of formation and laws ofmotion of debris flow[J].Acta Geographica Sinica,1989,44(3):291–301 doi: 10.3321/j.issn:0375-5444.1989.03.006
    [32] Zhong Deyu,Wang Guangqian,Zhang Lei.A Bed-load function based on kinetic theory[J].International Journal of Sediment Research,2012,27(4):460–472. doi: 10.1016/S1001-6279(13)60005-0
图(8)  /  表(2)

本文结构

    /

    返回文章
    返回