基于SPH的泥石流重力式拦砂坝削峰效应研究

彭章强 周伟 许佳路 杨立

摘要：拦砂坝削峰效应是衡量拦砂坝防治效益的一个重要指标。为研究泥石流重力式拦砂坝的削峰效果，利用SPH建立三维模型，对比分析数值模拟与物理试验两种不同方法的模拟结果，论证SPH模型的可用性;并以四川省绵竹市清平乡文家沟为例，估算不同泥石流事件下拦砂坝对过坝泥石流的削峰比例。研究结果表明：泥石流削峰效应受总库比和泥石流重度的影响，削峰率与总库比呈对数关系，与泥石流重度呈线性关系;在数值模拟结果与物理试验结果的对比中，两种方法的相对误差均值为12%，采用SPH模型研究泥石流削峰效应可行;对文家沟泥石流事件进行削峰分析时，获得了近年来文家沟泥石流的削峰率均值为53%，较好地量化了拦砂坝对泥石流的削峰效应。

关 键 词：拦砂坝; 削峰效应; SPH模型; 文家沟; 四川省

中图法分类号： P642   文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2022.01.007

0 引 言

2008年汶川地震诱发了大量的崩塌滑坡，在短历时强降雨触发下，堆积于坡面或沟道的松散物源极易转化形成泥石流。为降低泥石流造成的人员伤亡和经济损失，对人类活动区内的泥石流流域实施了工程治理。重力式拦砂坝（简称重力坝）是泥石流防治工程中较为常用的结构，占西南地区拦挡坝的88.89%[1]，其主要作用是拦截泥石流中的固体物质，削减泥石流峰值流量，大幅度降低堆积扇上的泥石流规模。削峰效应特指过坝后泥石流峰值流量的削弱程度，是评价重力坝防治效益的一个重要指标，对拦砂坝设计和泥石流防治效益评价具有重要意义，可为后续拦砂坝的设计提供参考。

目前，拦砂坝的研究集中于破坏模式[2]、防治功能[3-4]和坝体作用力[5]等方面，对重力坝削峰效应的研究不多。陈宁生等[1]对云南省大盈江浑水沟和东川市小水沟拦砂坝的削峰情况进行了分析，得出浑水沟与小沟拦砂坝的削峰比例分别为94.75%和66.7%。熊木齐等[6]指出在百年一遇条件下，甘肃省马槽沟一号坝的削峰比例为53.3%。还有部分研究者采用物理试验方法，研究了拦砂坝的削峰效应。贾世涛等[7]通过试验确定了拦砂坝开孔率为2.2%时，过坝前后流量比为开孔率为6.6%时的两倍。孙昊等[8]将试验数据进行拟合，分析了格栅坝在不同条件下峰值流量削减率与相对开度的关系。胡伟[9]利用物理试验，得到泥石流在拦砂坝后的回淤坡度和削峰效应的拟合曲线。

物理模拟试验可用于研究拦砂坝削峰效应，但物理试验周期长，投入的人力物力多，且最终结果易被试验误差影响。因此，本文采用数值模拟的方法，解决物理试验所遇到的问题。SPH是一种三维有限元模型，能完整模拟拦砂坝与泥石流相互作用时的三维应力状态，反映拦砂坝的削峰过程[10]。本文采用SPH模型建模，将模拟结果与削峰效应物理试验结果进行对比，判断方法可行性。最后利用该方法分析了四川省绵竹市清平乡文家沟在不同工况下的拦砂坝削峰效应。

1 SPH模型基本原理

光滑粒子流体动力学（SPH）模型是一种无网格的拉格朗日算法[11]，能描述泥石流的流速场及冲击力分布，可解决大变形及复杂边界等问题[12]。SPH模型以核近似和粒子近似为原则，将连续介质离散为粒子，采用离散元的方式分析粒子受力与运动状态。核近似是指将连续介质流体动力学的守恒定律从偏微分形式转化为适用于粒子模拟的形式。粒子近似指将连续介质的积分形式转换为支持域内粒子相关变量的求和[11]。

2.2 试验数值模型

数值模拟试验的三维模型如图1所示，尺寸和位置参考削峰物理试验模型。模型由水箱、水槽、拦砂坝（无泄水孔）和收集箱4个部分组成。流体从水箱流出，经过水槽加速后，越过拦砂坝，最终流入收集箱。

数值模拟计算过程与物理模拟试验相同。首先，根据水槽坡度与坝高确定拦砂坝的库容，根据总库比计算出流体总量，将得到的流体总量与拦砂坝尺寸一起输入三维数值模型进行模拟。在计算过坝前后峰值流量时，忽略拦砂坝宽度及流体回淤范围。以坝前流体稳定时，入库处的最大速度与此时的坝前流体截面面积计算坝前峰值流量。以流体越坝后，流体稳定时落地处的流体截面面积与该位置的最大流速计算坝后峰值流量。过坝前后峰值流量差值与坝前峰值流量之比即为拦砂坝的削峰率。

2.3 试验结果

数值模拟试验结果如表2所示。模拟结果表明：泥石流削峰效应受泥石流重度、水槽坡度和总库比的影响。试验共进行16组，对比前8组试验结果（泥石流重度为1 690 kg/m3），在水槽坡度一定的情况下，总库比越大，削峰率越低。分析原因，当一次泥石流总量远大于拦砂坝库容，拦砂坝对泥石流的拦挡作用不明显，削峰效果较差，削峰率低。在总库比一定的条件下，水槽坡度对泥石流削峰率的影响较小。分析原因，削峰率为过坝前后泥石流峰值流量之比，而过坝前后泥石流峰值流量均考虑水槽坡度，导致最终结果忽略了坡度对泥石流削峰率的影响。综合对比1～16組试验结果，在控制总库比和水槽坡度不变的条件下，分析泥石流重度对泥石流削峰效应的影响。发现重度越小，削峰率越大。泥石流重度减小时，其到达拦砂坝处的动能也将减小，容易在坝前淤积，拦砂坝对泥石流的削峰作用更加明显，削峰率大。

为进一步分析总库比与削峰率、泥石流重度与削峰率之间的关系，对模拟数据进行拟合，如图2所示。结果显示：总库比与削峰率呈对数关系，削峰率随总库比的增加而减小。当总库比足够小时（泥石流未越过拦砂坝），削峰率为100%;当总库比无限大时，削峰率为0。泥石流重度与削峰率呈线性关系，随泥石流重度的增加，削峰率减小。但当总库比过大时（总库比为100），泥石流重度对削峰率影响较小，此时拦砂坝对泥石流的削峰效应主要由总库比控制。

2.4 SPH模拟结果与物理模拟结果对比

数值模拟试验初始条件与物理试验（参照组）条件基本相同。对比两种方法的模拟结果（见图3），发现数值模拟结果略小于物理试验。分析原因，数值模拟忽略了水槽粗糙度的影响，导致泥石流削峰效应减小，削峰率降低。对两种方法进行相对误差分析，最大相对误差为20%（第2组试验），最小相对误差为1%（第12组试验），平均相对误差为12%，误差均在合理范围内。由此可认为采用SPH模型研究泥石流削峰效应是可行的。

3 文家沟拦砂坝削峰效应分析

在判断SPH模型研究泥石流削峰效应可行后，采用SPH模拟文家沟修建拦砂坝后泥石流的运动堆积，分析空库条件下拦砂坝的削峰效应。

3.1 研究区概况

文家沟地处四川省绵竹市清平乡，流域面积7.73 km2，相对高差1 540 m，主沟长度4.5 km，主沟纵比降425.8%。2008年汶川地震后，文家沟发生了一起高位顺层滑坡。2 750万m3的岩体从山顶滑下，巨大的滑体不断侵蚀和堆积，最终在文家沟形成约3 000万m3的松散堆积物。2010年8月12日晚至13日凌晨，清平乡遭受局部强降雨，累计雨量达227 mm。在强降雨的触发下，大量的松散堆积物起动[14]，形成泥石流。13日00：30左右，泥石流冲毁拦砂坝，形成溃坝型泥石流[15]。泥石流总量达310万m3，最大淤积厚度超过15 m[16]。

在文家沟“8·13”泥石流发生后，文家沟仍存在大量的松散堆积物，在暴雨条件下极易起动形成泥石流。为降低文家沟泥石流造成的危害，在该沟重新修建了泥石流防治工程，设计频率为20 a一遇。文家沟拦砂坝位置分布如图4所示，设计参数如表3所列。

3.2 文家沟拦砂坝削峰效应

在重新修建泥石流治理工程后，文家沟没有发生可以冲到沟口的特大规模泥石流，故本文以2010年文家沟已发生的4次泥石流为例，分析在空库条件下文家沟新建拦砂坝对泥石流的削峰效应。SPH模拟参数取值如表4所列，其中泥石流重度取自余斌等[16]计算结果。

模拟结果显示，拦砂坝对泥石流有明显的淤停和拦阻效应。图5为2010年7月31日文家沟泥石流建坝后的三维流速分布图，图中拦砂坝对泥石流有明显的淤停和阻挡作用，泥石流流速在过坝前后变化明显。本文以泥石流过坝后落地处的最大速度与截面面积计算泥石流坝后峰值流量。将坝后峰值流量与坝前峰值流量进行对比，计算每个拦砂坝对过坝泥石流的削峰比例。对4次泥石流事件[17]的模拟结果进行统计，如表5所列。在文家沟4次泥石流事件中，重力式拦砂坝对过坝泥石流的削峰率均值为53%，其中在同一泥石流事件中，3号坝库容最小，削峰率最低;在同库容条件下，2010年9月18日一次泥石流总量最小，削峰率最高。结合文家沟4次泥石流事件的体积和4个拦砂坝库容，发现在空库条件下，总库比越大，削峰率越低。随总库比的变化，削峰率在0～100%之间变化（当总库比无限大时，削峰率为0，泥石流未漫坝时，削峰率为100%）。

4 结 论

拦砂坝削峰效应分析是泥石流防治工程效益评价的重要内容之一。本文将数值模拟与物理试验进行对比，分析了采用数值模拟评价拦砂坝削峰效应的可行性，利用SPH模型对文家沟4次泥石流事件进行建模，计算了在空库条件下文家沟拦砂坝的削峰率，得出以下結论：

（1） 泥石流削峰率受总库比和泥石流重度影响，总库比越大，削峰率越低;泥石流重度越大，削峰率越低。

（2） 数值模拟结果与物理试验结果的相对误差在合理范围内，表明采用SPH模型对过坝泥石流进行削峰分析可行。

（3） 对文家沟4次泥石流事件进行模拟，得到空库条件下文家沟拦砂坝对泥石流峰值流量的削峰率，其范围为27%～100%，文家沟拦砂坝的削峰效应得以准确量化。

除总库比之外，影响拦砂坝削峰效应的因素还包括地形、拦砂坝位置和初始库容量等条件。本文仅分析了在空库条件下文家沟拦砂坝对过坝泥石流的削峰效应，对于半库及满库等其他工况，拦砂坝的削峰效应还需进一步探索。

参考文献：

[1] 陈宁生，周海波，卢阳，等.西南山区泥石流防治工程效益浅析[J].成都理工大学学报（自然科学版），2013，40（1）：50-58.

[2] 黄海，石胜伟，杨顺，等.2017年“8·8”九寨沟地震对景区泥石流治理工程影响机制研究[J].岩石力学与工程学报，2020，39（9）：1773-1786.

[3] 袁东，柳金峰，游勇，等.格子坝拦挡黏性泥石流闭塞过程及其闭塞判识试验研究[J].岩石力学与工程学报，2019，38（增1）：2881-2888.

[4] 曾庆利，岳中琦，杨志法，等.谷坊在泥石流防治中的作用：以云南蒋家沟2条支沟的对比为例[J].岩石力学与工程学报，2005（17）：3137-3145.

[5] 张莉，游勇，柳金峰，等.泥石流拦砂坝坝基土颗粒级配对扬压力影响试验研究[J].岩石力学与工程学报，2018，37（1）：190-198.

[6] 熊木齐，郭富赟，崔志杰，等.甘肃省武都区马槽沟泥石流特征及其治理工程效应[J].兰州大学学报（自然科学版），2015，51（6）：831-836.

[7] 贾世涛，崔鹏，陈晓清，等.拦沙坝调节泥石流拦挡与输移性能的试验研究[J].岩石力学与工程学报，2011，30（11）：2338-2345.

[8] 孙昊，游勇，柳金峰.泥石流梁式格栅坝调控性能试验研究[J].铁道学报，2019，41（3）：163-168.

[9] 胡伟.泥石流在拦沙坝后的回淤坡度和削峰效应试验研究[D].成都：成都理工大学，2017.

[10] CRESPO A J C，DOMíNGUEZ J M，ROGERS B D，et al.DualSPHysics：Open-source parallel CFD solver based on Smoothed Particle Hydrodynamics（SPH）[J].Computer Physics Communications，2015，187：204-216.

[11] 许波，谢谟文，胡嫚.基于GIS空间数据的滑坡SPH粒子模型研究[J].岩土力学，2016，37（9）：2696-2705.

[12] 韩征，粟滨，李艳鸽，等.基于HBP本构模型的泥石流动力过程SPH数值模拟[J].岩土力学，2019，40（增1）：477-485，510.

[13] LIANG H，HE S M，LIU W.Dynamic simulation of rockslide-debris flow based on an elastic-plastic framework using the SPH method[J].Bulletin of Engineering Geology & the Environment，2019，79：451-465.

[14] NI H Y，ZHENG W M，TIE Y B，et al.Formation and characteristics of post-earthquake debris flow：a case study from Wenjia gully in Mianzhu，Sichuan，SW China[J].Natural Hazards，2012，61（2）：317-335.

[15] YU B，MA Y，WU Y.Case study of a giant debris flow in theWenjia Gully，Sichuan Province，China[J].Natural Hazards，2013，65（1）：835-849.

[16] 余斌，马煜，吴雨夫.汶川地震后四川省绵竹市清平乡文家沟泥石流灾害调查研究[J].工程地质学报，2010，18（6）：827-836.

[17] LIU F Z，XU Q，DONG X J，et al.Design and performance of a novel multi-function debris flow mitigation system in Wenjia Gully，Sichuan[J].Landslides，2017，14：2089-2104.

（编辑：刘 媛）