基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型.pdf

泥石流断面内流速垂向分布是研究其流量、冲击力和沟床侵蚀过程的关键。然而，受限于测量装置布设条件，泥石流现场实测及水槽试验中常用的分层流速仪等设备仅能采集断面内少量样本点的流速数据，导致基于实测结果拟合回归的线性分布模型难以准确描述泥石流速度分布规律。对此，本文依托大比例尺泥石流水槽试验开展研究，利用所构建的基于 本构的光滑粒子流体动力学（）数值模型反演泥石流三维动力过程，通过分层统计算法对大量粒子速度数据进行分析处理，获得了断面内速度垂向分布规律，并据此提出了基于对数函数的泥石流流速垂向非线性分布模型。为验证所提出模型的准确性，利用其他多组水槽试验数据进行了对比分析，结果表明，本文提出的对数分布模型比传统线性分布模型能够更准确地拟合速度剖面，并在模型参数敏感性方面具有更强鲁棒性。关键词泥石流；流速；垂向分布；对数模型；数值模拟中图分类号：文献标识码：收稿日期：；修回日期：基金项目：国家重点研发计划专项（资助号：），国家自然科学基金面上项目（资助号：），湖南省自然科学基金优秀青年项目（资助号：）（），（）（）第一作者简介：韩征（），男，博士，教授，博士生导师，主要从事地质灾害防治与数值模拟方面的研究工作 ：通讯作者简介：李艳鸽（），女，博士，副教授，博士生导师，主要从事灾害风险评价与遥感解译方面的研究工作 ：（，）（，）（，）（，）（，），（），；引言泥石流是由山区松散堆积物和水二相混合的非牛顿流体，与水流的特性相似，其流速具有显著的垂向分布特征（，；，；，）。不少研究（，；，；王东坡等，）认为，表征流深和流速相关性的流速垂向分布是理解和预测泥石流侵蚀行为和冲击力大小的关键。对于泥石流的侵蚀行为，理论公式表明计算泥石流侵蚀量需要获取流体的基底流速（，），而基底流速难以直接通过试验测得，只能根据流速的垂向分布以及整体的平均流速进行计算；对于冲击力的计算，由于泥石流的冲击力与流速的平方呈正比关系，流速的垂向分布造成了冲击力深度方向上的不均匀（马宗源等，）。目前设计拦砂坝和桥墩等结构物时尚未考虑冲击力的分布，使得结构在泥石流作用下受力不均，导致近年来防治工程结构失效致灾的案例时有发生。因此，明确泥石流流速的垂向分布对于泥石流动力过程分析和防治工程设计具有重要意义。目前，泥石流流速垂向分布的研究主要采用野外原型观测与模型试验两种方法。（）和 （）通过现场观测数据提出了泥石流流体表层速度的横向不对称分布，即沟道中泓线速度大于两侧的速度。康志成等（）和 （）也观测到泥石流的流速在流体自由面比底部高得多的垂向分布特点。但由于泥石流的发生往往没有预兆并且持续过程较短，采用现场观测中较难采集泥石流流速，因此许多学者通过实验室水槽实验对泥石流的流速分布进行研究（韦方强等，）。典型研究有：（，）首次在水槽试验中测量出泥石流的速度剖 工程地质学报 面，并且基于 模型推导出流速分布呈现出幂律 关 系；（）、（）和 （）的试验证实了泥石流自由面流速高于底部流速的分布特点，并提出了不同的线性分布模型来拟合测量结果。但总体而言，泥石流流速垂向分布方面的相关试验研究较少，这一方面是由于泥石流动力过程复杂，作为一种二相混合流体，固相的颗粒粒径繁杂，试验中往往使用理想级配的砂石进行模拟，难以反映复杂的泥石流流变特性；另一方面，在试验中由于仪器与布设条件限制，很难详细准确地采集不同流深位置的流体速度，例如，适用于洪水等透明介质的高速摄影测量技术（）很难观测到泥石流这种泥浆与砂石的不透明流体的内部流态，而常用的分层流速仪由于仪器体积较大，一个断面内仅能布设少量样本点采集数据，因此难以准确还原泥石流流速垂向分布规律。鉴于上述难点，在已有的流速垂向分布模型中，目前采纳较多的是 （）提出的线性分布模型，该模型在一些研究（，）中得到了应用。然而该模型将断面上的流速剖面简化为流深方向上的线性分布，因此为了得到更好的速度拟合结果，通常需要在该模型中引入一个经验性分布系数，的建议取值范围为 （，）。由于不同深度的流速值与该系数线性相关，因此敏感性较强。在前期研究（，）中，我们发现 时能够获得拟合最好的流速分布情况。实质上，引入经验性分布系数 的泥石流流速垂向线性分布模型是在断面平均流速的基础上沿流深方向进行的流速调整，虽然在流体中部拟合结果较好，但在流体自由面和底部则会产生较大的误差（，）。所以总体看来，泥石流流体在流深方向上的速度分布特征更近似于非线性分布，这也是近期一些关于二相流、碎屑流流速分布研究（，）所观察到的一致现象。对泥石流流速垂向分布进行非线性回归分析需要流深方向上大量的流速样本数据。而基于大量离散粒子的三维光滑粒子流体动力学方法（）能够从三维角度对泥石流动力过程进行模拟（胡凯衡等，；乔成等，；，；，），也为获取大量的流速样本数据提供了有效手段。传统的泥石流数值研究大多在欧拉网格内，利用浅水波假设对纳维 斯托克斯方程进行深度平均的二维简化，并进行有限差分求解，尽管可以获得复杂三维地形上泥石流流速及流深的时变分布，但无法解析流体的速度剖面（，；，；，）。而三维的非连续大变形分析方法（）、颗粒流分析方法（）等方法受限于本构模型，只能将泥石流作为碎屑流进行模拟，无法反映其真实的流态特征（胡明鉴等，；王学良等，；郑博宁等，；刘伦杰等，；潘青等，）。相比于传统网格方法，三维 方法因其粒子特性，在解析流体断面速度场的方面具有很大优势，由于泥石流流体质量被视为一组离散的颗粒，其行为可以通过直接求解纳维 斯托克斯方程来描述，不需要深度积分和浅水假设。典型研究如：（）最先将其用于泥石流、滑坡问题中；缪吉伦等（）从纵向及深度方向对泥石流运动纳维 斯托克斯方程进行离散，建立了二维 堆积形态数值模型；（）也运用 方法从三维尺度上成功实现了泥石流运动的初步模拟。我们也从泥石流流变模型角度，构建了基于 （）本构的三维 数值模型，实现了泥石流动力过程（韩征等，）及沟床松散堆积物侵蚀过程的有效模拟（，）。上述研究表明，三维 方法能够有效模拟泥石流动力过程及流态，为解析泥石流垂向速度分布规律提供了支撑。本文在前期研究所构建的三维 数值模型（韩征等，）的基础上，对 年美国地调局泥石流水槽试验（，）进行了模拟，反演泥石流的动力过程，并且利用粒子分层统计算法提取了断面流深方向的流速分布特征，通过大量流速样本数据的回归分析，提出了基于对数函数的流速垂向非线性分布模型，并通过其他多组试验数据将其与线性分布模型进行了对比分析。结果表明本文提出的分布模型能够更好地拟合流速剖面，且具有鲁棒性。模型和方法 流速垂向线性分布模型如前所述，不少研究已经建立了泥石流垂向分布的线性分布模型。典型研究有：（）通过大型水槽试验对泥石流的三维流速分布进行了研究，通过颜色标记块和高速摄影技术测 （）韩征等：基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型量了泥石流的表面流速，之后建立了关于流动的形状和速度场函数的数学模型。为了从深度积分速度确定三维速度场 ，需要对速度剖面进行描述。确定水平速度分量 和 需要通过：（，）（）（）式中：（）表示一个假定的速度分布，速度的垂直分量通过在相对于 的移动坐标系中计算质量守恒方程来获得：（，）()（）速度分布图 （）虽然很难在试验中直接观察得到，但是可以从以下几个方面加以约束。在整体静止的系统中，深度平均速度 （，珋）满足：珋 （）因此，速度分布 （）需要满足：（）（）从物理规律上，也期望流速在流动表面最大，并随着流动深度而减小，将 约束为非递减函数。因此考虑单参数的流速剖面如下：()()（）该公式的详细推导可参考 （）的研究，式中：参数（）用于控制大部分流体内的剪切量，当 时，（，），表示流体上下等速流动；当 时，（，），表示整个流体中速度随深度呈线性变化。为了进一步利用速度分布计算泥石流的侵蚀率，将基底流速 与更容易观察到的深度方向平均流速珋 联系起来很有必要，根据式（）可以得到如下关系：珋 （）珋（）当取 时，即为基底流速 。虽然这一线性分布模型在流体中部拟合结果较好，但是由于泥石流过程的复杂性和偶然性，该模型在流体自由面和底部会产生较大的误差。数值模拟方法在 方法中，泥石流等流体被视为连续体和不可压缩的流体，并用一组离散化的粒子表征，其行为可以通过求解纳维 斯托克斯方程来描述，从而提供了一种获取 个维度上的流速场分量的数值解。为了模拟泥石流的动态特性，本文采用了基于 本构的 模型（，；韩征等，）（即三维 模型）。流变模型表达式如下：（）（）（）（）（）式中：为剪切应力张量；表示等效黏性系数；为局部应变率张量；表示 黏性系数；和 分别为控制不同剪切率下应力增长的常数和幂律指数；为 屈服准则的屈服应力；表示土体的黏聚力；表示内摩擦角；表示正应力；表示剪切应变率，定义为：槡（）在拉格朗日形式下，三维 框架整合 流变模型，由动量守恒方程组成的纳维 斯托克斯方程有效地描述了泥石流的动力过程。()()（）式中：表示核函数；、分别表示粒子速度和重力。我们前期研究（，；韩征等，）已表明，模型在泥石流动态过程分析中具有较好的适用性，能够更好的模拟真实泥石流的流速分布规律。粒子分层统计算法利用 方法对泥石流动力过程进行模拟后，可以获得不同时间步下大量粒子的空间位置和速度矢量，为拟合分析泥石流不同流深的速度分布情况提供了大量数据样本。其中：粒子的空间位置坐标是位于整体坐标系中的 ，个方向的位置，速度矢量也可以解析为整体坐标系下 ，个方向的速度分量。由于泥石流流动方向大多不与整体坐标系的 工程地质学报 轴平行，因此直接导出的粒子数据不可以直接用来绘制剖面图。为了方便研究泥石流的速度垂向分布规律，往往需要对坐标系进行转换操作，将所得粒子数据的坐标系转为：轴平行于泥石流的流动方向并以泥石流的前进方向为正，垂直流动方向向上为 轴方向。我们采用以下公式进行坐标转换操作。（）（）（）（）此处，针对一般的倾斜水槽试验，有如下计算公式。（）（）（）（）（）（）（）（）（）（）（）（）类似的，采用以下公式进行速度的方向变换，由于速度是一种矢量，不需要对其进行平移变换，只需要对其方向进行变换。（）（）（）（）之后需要筛选位于沟道某一剖面范围内的粒子数据，即选定某一个断面的顺泥石流方向的 坐标以及选取的统计粒子的宽度范围 ，通过坐标条件 判断粒子是否位于所选取的断面范围内，以此确定位于作用域范围内的粒子编号及其数据并保存。对于断面范围内的粒子数据，由于粒子并不是真实存在于一个断面上，同一高度处可能存在多个粒子速度，因此无法直接绘制出流速剖面图，需要先按照 的流深方向上按深度位置进行分层统计，计算每一层的流深和平均流速，这需要首先获取作用域范围内粒子的最大深度和最小深度。（，）（，）（）为了数据的清晰和方便进行分析，将最高点的坐标值和最低点的坐标值按照单层高度 进行取整，取整后的最高点和最低点坐标记为 和 。（）则第 层的平均高度为：（，）（）式中：表示粒子分层的单层高度；表示粒子分层的数目。（）为了计算第 层的平均速度，需要统计第 层粒子的总数 ，。，()（）计算第 层粒子沿流动方向上速度分量的总和 ，：，（）计算第 层粒子的平均速度：，（）通过处理后得到不同流深层的粒子 方向的平均速度及其对应的高度绘制速度剖面图，其中以平均流速作为横坐标，以高度作为纵坐标。该分层统计算法已经通过 实现。粒子分层统计方法的示意图和程序流程图如图 、图 所示。大型泥石流水槽试验数值模拟 数值模拟为分析泥石流速度的垂向分布规律，我们采用基于 本构模型的 方法对美国地质调查局大尺度水槽试验案例开展了数值模拟实验。（）开展了一系列大尺度泥石流水槽试验，用于探究泥石流的复杂动力过程。其通过改变沉积物的体积含水量，共计完成了 组试验，包括两个不设置沉积物层的对照试验，并且在观测点 处安装了孔隙水压力传感器、流深传感器和侵蚀传感器，在 到 处还设置了多组侵蚀传感器。上述 组试验中，各组的初始泥石流浆体体积都是，泥石流的浆体是由 砾石、沙 （）韩征等：基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型图 粒子分层统计算法示意图 图 流速剖面分析计算流程图 子和 泥浆细颗粒混合而成。由于该试验对真实尺度下的泥石流进行了多次试验，获得了大量可靠性高并且可以对照的试验数据，为数值模拟研究提供了有力的数据支持。本文采用 方法对该试验进行了数值模拟，数值模型如图 所示，在 数值模拟中，水槽被视为由固定颗粒组成的边界，泥石流体也被离散成大量可以运动的粒子。采取初始粒子距离为 的长度建立粒子模型，离散结果产生了 个边界粒子和 个流体粒子。数值模拟采用图 泥石流水槽试验数值模型 高性能服务器进行计算，服务器配置有 颗 核心 、块 （包含 显存 ，个 核心，理论单精度浮点数运算能力：）、内存。使用服务器计算一组 时间步的水槽试验需耗时 ，本文一共完成了 组数值模拟。为了验证速度模拟的准确性，我们将模拟泥石流前端位置坐标的时间序列与 （）的实验数据进行比较。如图 中所示，数值模拟中的前缘位置与实验数据中的位置非常接近，这意味着泥石流体的总体速度在数值模拟中与实验结果吻合良好。图 还展示了不同时间（，）沉积物铺床段不同位置的基底流速和流深云图。工程地质学报 图 流体前端位置随时间的变化 图 不同时间点的流深和基底流速云图 ；流速分布特征流速分布的结果如图 所示，图中分别显示了采用 的数值模型分析得到的 断面处的粒子数分布和 个典型时间步的流速剖面。由于这些时间点的流体粒子数量较多，分层统图 流速剖面分析结果 不同时间步的断面内粒子数目；时间步 ；时间步 ；时间步 ；时间步 计后的流速误差更小。结果表明，在时间步 时间步，流深达到峰值 ，流体的速度剖面呈现出流速由流体自由面沿流深方向非线性减小的分布形式，基底流速较小，约为 ，流体自由面流速较大，约为 ，此时间步的流速的最大值出现在 的高度处，并不位于流体表面。随着时间的推移，流速剖面逐渐由非线性分布向线性分布演变，并且最大流速相比于之前呈现逐渐减小。从时间步 到时间步 的两个时间步可以看出流速剖面的弯曲半径逐渐变大，逐渐接近于直线分布，且表层流速和基底流速的差异逐渐减小。流速剖面呈现的规律表明流体的切向流速基本上符合从基底到表层随着高度增大而增大的规律。但在时间步 和时间步 的速度剖面中观察到了表面流速减小的现象，对于这种表面流速降低的问题有待进一步研究。如图 所示是时间步 时间步对应的流速剖面中 到 的 个不同深度处的流速分布。该时间步中，所选断面共计分布粒子数为 个。为了说明每一层的粒子速度的分布情况，对 深度层的粒子速度进行了单独统计分析，绘制出速度的频率分布直方图。由于单层粒子数量较多，可以近似将粒子数看作频率，并且采用正态分布 （）韩征等：基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型图 单层粒子流速箱线图 曲线进行拟合。可以看出每一层的速度均值都分布在概率最大的位置，采用速度均值代表单层粒子速度较为合理。图中速度均值的连接线形成的曲线就是此断面位置的速度剖面。线性分布模型 演化规律采用 （）提出的垂向流速的线性分布模型 进行描述时，由式（）可知 存在以下关系：珋（）式中：表示泥石流基底平行于水槽方向的流速，在本次研究中，我们采用最底下一层的流速平均值作为基底流速 ；珋 表示泥石流垂向断面的平均流速，采用如下公式进行计算。珋 ，（）由于在泥石流峰值流量集中在前 个时间步以内，后续经过监测点的都是流量较小的流体，为了避免粒子数量过少带来的误差影响，只取前 个时间步的速度剖面进行分析。由监测点 处的速度分布数据计算得到的 随时间变化的关系如图 所示。图中可以看出 的最佳拟合值是变化的，在泥石流洪峰刚到达监测点时，迅速增大到 ，在泥石流洪峰流过之后，随时间逐渐减小，并且最终稳定在 之间。在前 个时间的范围内，的变化范围为 ，可以看出，在这样大的变化范围内想找到一个恒定的的值来拟合整个流动过程十分困难，因此线性分布模型描图 随时间变化的关系 述泥石流垂向流速分布存在显著局限。而从 图 可以看出流速剖面呈现非常明显的非线性分布。因此，构建能够更加精确拟合流速垂向分布规律的非线性模型具有重要意义。流速垂向分布非线性模型 流速分布拟合方程为了明确泥石流断面流速垂向分布的规律，对 处的流速剖面进行了非线性回归分析，（）提出流速剖面曲线可以用自然数为底的指数函数描述，但本文为了便于描述不同深度的流速，需要采用形如 （）的函数模型，因此选用对数函数较为合理。而在绘制流速剖面时为了直观表现速度的分布规律，需要交换横纵坐标后进行绘图，仍然选择将流速作为横坐标，流深作为纵坐标。我们针对时间步为 、的 个流速剖面进行了回归分析，结果如图 所示。可以看出，基于对数形式的函数方程很好地拟合了 个时间点的速度剖面。为了简化分析，个剖面均采用同一对数函数形式：()（）式中：和表示流速和深度进行归一化后的无量纲值，其取值范围均为 ，为了方便叙述，以下分别将其记作 、；表示此断面上的最大流速；表示此断面上的流深；是引入的拟合控制参数。该模型假定流速最大值出现在流体自由面，即函数始终满足当 时，。参数 的取值和回归指数 如图 所示，结果表明随着时间步从 工程地质学报 图 流速剖面的拟合回归分析 时间步 ；时间步 ；时间步 的增大，函数关系式的拟合参数 的取值从 向 逐渐减小，并且拟合曲线弯曲弧度逐渐增大。在时间步 ，时，函数图像在 的范围接近于线性分布，这与实际分布规律相吻合。图像还可以看出，在 位于 的区间，拟合曲线与数值模拟数据吻合度很高，在流体靠近自由面和底部的区域拟合程度相对较差。整体上看，采用对数函数进行分析得到的拟合曲线匹配了不同时间步的速度剖面，通过参数 的不同表现出流速分布从非线性逐渐向线性变化的趋势。拟合方程的参数分析为了进一步探明参数 对流速剖面的影响，本文对不同取值得到的剖面进行研究，采用 、以及 的取值作了分析，结果如图 所示。可以看出随着参数 的增大，曲线的曲率半径逐渐增大，曲线弯曲的程度逐渐减小，相同的自变量对应的因变量的值减小，流速剖面的上半部分逐渐接近于线性。讨论回归分析得到的非线性分布函数很好地拟合了模拟结果，为了验证本文所提出的对数函数分布模型（如式（）所示）的有效性，并给出参数值 的建议取值范围，我们利用 （）及 （）前期研究所公开的 组水槽试验中数据进行拟合，并将得到的拟合结果与已有的线性模型进行比较。在 （）的试验中，水槽宽 ，长 ，水槽底部人工覆盖一层由直径 的沙粒组成的薄砂层，水槽的倾斜角图 拟合函数参数 分析 度可以调节，以便产生不同的速度。采用控制水量的方式控制最大流量和水流深度，水槽旁架设用于记录不同深度流速的高速摄像系统。为了更加准确评估拟合函数的性能，我们将拟合结果与试验测量值进行比较，并采用误差平方和（）来评价拟合结果。的计算过程如式（）所示，越小表示拟合结果与试验值差别越小，拟合程度越好。（）（）式中：表示实测的速度值；表示拟合得到的速度值；表示试验数据点的个数。我们针对每组试验数据进行同样的迭代计算，将 从 的区间离散，分别计算每个离散值 对应的 值并绘制成曲线。如图 所示是以第 组数据为例迭代计算的结果，图中最小的 值 对应的 就是此组数据拟合最优的参数 值。组数据拟合的结果汇总于表 ，并且选择了数据 和数据 的拟合结果绘制成图，如图 所示。由表中数据可以看出，组数据中有 组数据 （）韩征等：基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型图 随参数 变化的关系 （数据 、）都是非线性模型拟合的 最小，说明了非线性分布模型的有效性。另外，为了对比非线性分布模型与线性分布模型（如式（）所示）的鲁棒性差异，我们采用同样的迭代方法分析了每组数据拟合最佳的 值，结果如表格最后一列所示，然后分别对 组试验数据得到的 和非线图 非线性分布和线性分布对比图 数据 对比分析结果；数据 对比分析结果性分布参数 绘制了箱线图，如图 所示。由图像可以看出参数 的取值更为稳定，其箱体的变化范围大小为 ，远小于参数 的变化范围 ，后者约为前者的 倍。因此，基于对数形式的非线性分布对比线性分布具有更好的拟合精度。同时，从图 所示的箱线图中可以看出，组试验数据拟合得到的参数 的中值为 ，最大值为 ，最小值为 。考虑 组数据的概率分布，我们根据研究（，）所用的方法，选取上四分位数 与下四分位数 的近似值作为非线性分布的拟合参数 的建议取值范围，即所有数据中，参数 有 的拟合结果分布于 ，的区间内。而对比图 中线性分布参数 的上、下四分位数范围 ，表 （）及 （）试验数据的拟合结果对比 （）（）数据编号误差平方和（）误差平方和（）误差平方和（）误差平方和（）线性最优的 图 参数 和 箱线图 非线性分布的拟合参数 的建议范围是更加靠近于中值的，这也表明基于对数函数的非线性流速分布 工程地质学报 模型具有更强的鲁棒性。结论本文采用前期研究所构建的基于 本构模型的光滑粒子动力学（）模型，对美国地调局大比例尺泥石流水槽试验进行模拟，反演其三维动力过程，然后利用提出的粒子分层统计算法分析了断面内流速垂向分布规律，通过大量数据的回归分析提出了基于对数函数的流速垂向非线性分布模型，依托其他多组水槽试验实测数据，与基于 的线性分布模型进行了比较研究。本文取得主要结论如下：（）利用 对泥石流水槽模型试验进行了模拟，通过粒子分层统计得到了流深方向上的流速剖面，断面内的流速分布呈现出流速由泥石流流体自由面沿流深方向非线性减小的分布特征，基底处的流速最小，最大流速一般出现在流体自由面，部分时间最大流速出现在流体自由面以下一定距离位置。并且随着时间推移，流速剖面出现由非线性分布向线性分布变化的趋势。（）建立了基于对数函数的非线性分布模型，该模型能够较好地拟合流体深度中部区域的流速剖面，在流体靠近自由面和基底的深度拟合程度相对较差。非线性分布模型引入了参数 ，对参数 的敏感性分析表明，随着 增大，流速分布的曲率逐渐增大，弯曲程度减小，逐渐向线性分布模型变化；通过对试验数据的拟合分析得到参数 的建议取值范围为 。（）通过多组水槽试验实测数据对非线性分布的准确性和鲁棒性进行了分析，并且与线性分布模型进行了对比，结果表明基于对数函数的非线性分布模型能够更准确地拟合断面内的速度剖面，且模型参数敏感性方面具有更好的鲁棒性，参数 的变化范围仅为线性分布参数 的 。参考文献 ，（）：，：，：，（）：，（）：，：，：，：，：，（）：，：，：，（）：，：，：，（）：，：，（）：，（）：，（）：，（）：，（）韩征等：基于流体粒子分层统计算法的泥石流流速垂向非线性分布模型 ，（）：，（）：，（）：，：，：，（）：，（）：，（）：，（）：，?（?）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?（）?工程地质学报