基于ANSYS CFX的泥石流数值模拟

李战国

(新疆昌吉方汇水电设计有限公司，新疆 昌吉 831100)

泥石流是由山坡和沟渠产生的小石块和大石头构成的水和泥沙固液两相流体。当满足大量松散的固体物质、地形条件和一定强度的降水条件时，则可能会发生滑坡[1]。当山体滑坡发生时，泥石流裹挟的压力会对山谷地形以及生命财产造成一定程度的损害。山体滑坡引起的巨大土石堆积会使肥沃土地变为荒凉戈壁，其对自然的破坏力不言而喻。目前，对泥石流的研究主要包括泥石流流体模型的研究、泥石流的预测预警、泥石流状况评估、泥石流特性的数值模拟、各种泥石流预防结构和系统的研究等[2]。

1 研究区概况

洞峡子沟泥石流位于大川河右岸，南天门路通过沟槽(图1)。

图1 洞峡子沟沟谷范围图

从2011和2012年累计的滑坡总数来看，本次滑坡为中型滑坡。洞峡子沟地区地形通常为高腐蚀性的中型山地地形，其海拔高差在1 100～2 260 m之间，地势陡峭，沟渠形状通常为不规则梯形。洞峡子沟泥石流主要威胁沟阳凯村居民18户68人和移民安置区45户183人住宅共45栋(即将建造)，潜在威胁财产近1 300万元，同时威胁沟渠南天门道路交通安全。

2 模拟方法

2.1 计算流体动力学数值模拟方法

目前在流体计算中使用3种数值方法，即有限差分法、有限元法和有限体积法(也称为控制体积法)。其中，有限体积法具有计算效率高、收敛性好的特点，因此在CFD软件的开发中得到广泛的应用[3]。

有限体积法将分析对象分割成多个彼此不重复的控制体积，然后采用控制方程以积分插值的形式用于每个控制体积，体积界面的物理尺寸近似于物理节点的数量。区域离散化的本质是用有限的离散点代替原来的连续空间。离散区域过程首先将计算出的区域划分为若干个不重叠的细分(即网格)，以确定每个细分节点的位置以及这些点所指示的控制点。这个区域管制流程完成后，4个最基本的要素是：①节点，需要解决的物理数量的几何位置；②应用控制体积、控制方程或规律的最小几何因素；③接口，接口规定了与每个节点对应的控制点的接口位置；④网格线，连接相邻节点形成的曲线簇。在单个进程中，使用节点作为控制点的代表，在该节点上定义和存储每个控制点的物理数量。图2、图3显示了A、B和P是节点编号，φA、φB为节点边界条件值，N、E、S、W是北、东、南、西，n、e、w、s是相应节点单元编号，Δx和Δy是一维、二维有限体积法的P节点的X和Y方向的长度。本文采用ANSYS CFX软件进行相关模拟。

图2 一维问题有限体积法计算网格

图3 二维问题有限体积法计算网格

2.2 泥石流流变模型的选择

研究泥石流流变的特性对于研究其流动规律、了解其性质和数值模拟滑坡具有重要意义。泥石流的流变特性揭示了从泥石流获得的剪切应力与相应的液体流变速度梯度之间的关系。选择流变模型有助于更深入地了解泥石流流变的内在规律。目前，研究人员已经根据实验数据的理论推导和统计分析开发了多种流变模型。最常用的流变模型有牛顿流体模型、空流体模型、假塑性流体模型、膨胀流体模型[4]等。实验表明，典型的泥石流流体特性的数学模型公式如下:

式中：τ为剪切应力；τb为屈服应力；η为黏性系数;du/dy为剪应变。

剪切应力主要是由泥石流中黏性粒子的凝集结构和沉积物粒子之间的摩擦形成的，黏性粒子的凝集结构是宾汉极限剪切应力的主要组成部分，其影响因素与粒度和微粒物质的含量有关，粒子越细，微粒含量越高，产生宾汉极限剪切应力的临界浓度就越低。实验证实，宾汉流体模型更适合一般泥石流流体的特性，因此本次模拟使用宾汉流体模型。

2.3 几何建模

几何建模选择洞峡子沟主渠道大坝建模。结构形式为普通重力块坝。

将选定的研究流域、大坝尺寸和材质与现有研究模拟计算器相结合，根据宾厄姆流体和湍流模型选择单向流固耦合方法。

1) 坝体大，不易变形，变形对整个流场的流体分布影响很小。

2) 若使用双向流固耦合需要设置移动网格，难以确定移动网格计算的准确性，过于依赖手动经验，以及在研究体变形较小的前提下生成移动网格的计算误差可能会更大。

3) 双向流固耦合计算对计算机性能要求很高，计算速度较慢。此模型最好使用单向流固耦合，因为在采用单向-固态耦合的100 s、步长设置为0.1 s的英特尔i7-7800k 32gb ddr4平面下，累积时间超过44 h。

4) 根据经验，单向流固耦合在计算详细应力应变时比双向流固耦合更加精确，前提是变形较小。因此，使用单向流固耦合模拟此模型。流固耦合模拟是两个条件：第一个是大规模泥浆-流池溢流，第二个是中小土流逐渐沉入整个仓库时的情况。流场和坝宽5 m，耦合分析的初始流场数据来自整个流域流场分析的相应部分。使用CFX执行流场计算，然后导入ANSYS执行流固耦合计算。流固耦合模型见图4。如果将流场网格“躯干大小”设置为0.5m并加密FSI面，则“面大小”设置为0.25 m。工作条件2，进口面小，进口面也加密，面大小设置为0.2 m，最终结果节点(nodes)为54 544 607，单位(Element)为511 280(图5)。

图4 流固耦合模型

图5 流体网格FSI细节加密图

需要特别指出的是，FSI面大小必须与坝体网格的大小相匹配。否则，将出现overflow(溢出)错误。最终获得坝网的节点为2 614 341，单位为623 280(图6)。

图6 流体网格细节展示图

“基于SYM面的属性”(Boundary type)是symmetry属性。模拟坝体左右两侧的工作状态与中间的工作状态相匹配，因此选择20 m的坝长进行模拟(图7)。

图7 SYM流场边界图

2.4 流固耦合计算参数设置

见表1。

表1 流固耦合参数设置一览表

3 模拟结果

根据流场数据计算出的流体逐渐渗透到整个水箱中时，自由水平参考图8，图8显示了每个时间速度向量，总液体压力参考图8。

图8 逐渐淤积至满库过流工况泥流流体与拦挡坝相互作用时不同时刻的流速矢量图

t=1 s时，二流流体开始流入二流块水坝，倾斜将使二流流体“水龙头”部分的速度逐渐加快。t=4 s时，泥浆流体“水龙头”部分速度为12.3 m/s。t=8 s，当泥浆流体到达坝底时，泥浆流体“水龙头”速度为9.9 m/s，此时泥浆流开始沉入坝内。t=38.4 s时，上面的泥浆流体沉入水平持续30.4 s到坝顶，当泥浆流体的液体水平足以通过泥浆水坝，第一次流向大流体时，流体速度为t=8 s的“水龙头”速度明显下降，从而逐渐降低液体水平。如果t=54.2 s，则土流水平会持续提高，在整个计算过程中达到最高值，下游流体可能会增加，而大坝前下游速度为4.5 m/s，则此过程可以更准确地说明大坝位置对控制土流的重要性。

将结果流场数据导入耦合计算后，获得内部实体的应力应变和位移等主要参数。当泥浆流再次沉入整个水箱时，通道和水坝的液体最大压力位置在水坝入口底部，压力值为1.891 0×105Pa，水坝前面的最大压力位置是负流冲击区，冲击压力约为1.225 2×105Pa。

大坝的动态响应分析结果显示，最大大坝应力强度位于大坝进水口表面的顶部，最大应力值为4.304 2×105Pa，大坝入口表面的应力强度也较大。块坝的弹性变形强度分布基本上与坝的应力强度分布相匹配，最大应变值为2.224 1×105Pa。

在整个水库过流条件下，块坝是整体拱形层分布。随着坝高的增加，位移值也会增加，最大位移值为4.790 0×10-5m。

4 结 论

泥石流灾害的防治大部分采用经验和规范，很难对所有泥石流流体的管理采用相同的加固标准。因此,要对每条泥石流沟进行单独的数值模拟计算，以便对后防措施提供参考。在此次泥石流的数值模拟中，选择合理的泥石流模拟参数，以便软件模拟泥石流的起始、流动和整个流域流场情况。本文研究成果可为泥石流流场和大坝流固耦合计算、泥石流流场的动态特性和拦截工程的动态响应分析、泥石流流场动力学研究以及泥石流防治工程的理论研究提供参考和借鉴。