刘家峡洮河口排沙洞水下岩塞爆通后多相流运动数值模拟

徐小武，薛立梅，李俊富，吴希华

（1.中水东北勘测设计研究有限责任公司，吉林 长春 130061；2.水利部松辽水利委员会，吉林 长春 130021）

1 概 述

刘家峡水电站位于甘肃省永靖县的黄河干流上，主坝为混凝土重力坝，最大坝高147 m，总库容57.01亿m3，总装机容量1390 MW，为一等大（1）型工程。

洮河水少沙多，在距刘家峡大坝上游1.5 km的右岸汇入黄河干流。洮河口段死库容于1987年淤满，在河口附近黄河干流形成沙坎，且淤积面逐年抬高，现有排沙设施已不能解决洮河泥沙淤积并向坝前推移的问题，给电站的安全运行和渡汛造成了严重危害。

为解决泥沙淤积对大坝安全运行带来的严重危害和对大坝安全渡汛构成的威胁及电站机组磨损等问题，在刘家峡水电站左岸增建排沙洞工程，排沙洞的进水口采用水下岩塞爆破方案。

岩塞爆破口位于洮河出口，黄河左岸；排沙洞全长1486 m，岩塞进口底板高程为1664.53 m；岩塞内口为圆形，内径10 m，外口近似椭圆，外口尺寸约20.3 m×27.84 m，岩塞最小厚度12.3 m，岩塞进口轴线与水平面夹角45°，塞体方量约2474 m3。

2 计算目的和内容

排沙洞进口段采用水下岩塞爆破施工。在排沙洞岩塞口爆破区范围内，淤积层厚度变化较大，11～58 m；另外，淤积层内有522 m3的水下人工堆渣，渣块块径一般5～10 cm，最大块径100 cm。厚淤泥层的存在对岩塞口爆破会带来不利影响；另外，在排沙洞运行过程中，若渣块进入洞内，也会影响排沙洞的正常运行，成为工程的隐患。为了进一步掌握排沙洞水下岩塞爆破后的水沙运动规律，有必要开展排沙洞水下岩塞爆通后多相流运动数值模拟计算，为排沙洞水下岩塞爆破方案设计提供参考。

为了解水下岩塞爆破后厚覆盖层对排沙洞内水沙运动的影响程度，本文基于当前多相流计算技术的发展水平，假定在排沙洞水下岩塞爆破后，爆破岩塞所产生的岩渣瞬间全部堆满预挖的集渣坑，忽略岩塞岩渣和集渣坑的影响。采用Euler多相流模型，以排沙洞岩塞口区及排沙洞为研究对象，建立三维排沙洞多相流非恒定多相流计算模型；计算正常蓄水位、淤泥层不同颗粒组成等条件下，排沙洞水下岩塞爆破后排沙洞中多相流运动过程；对比分析不同时间过程排沙洞内各相体积率及其流动速度、压力分布等。

3 计算模型原理和方法

3.1 多相流模型

排沙洞岩塞口区覆盖层较深且夹杂人工堆渣，在排沙洞岩塞爆通后,岩塞口区的淤泥将夹带人工堆渣流入排沙洞。夹带堆渣的淤泥在排沙洞内的流动是一个多相流问题，需要建立排沙洞多相流数学模型来研究分析排沙洞内淤泥的流动规律。目前，处理多相流多采用欧拉模型的计算方法。

欧拉模型是一种非常复杂的多相流模型，它包含有n个动量方程和连续方程，来求解每一相；压力项和各界面交换系数是耦合在一起的。耦合的方式则依赖于所含相的情况，颗粒流(流—固)的处理与非颗粒流(流—流)是不同的。对于颗粒流，可应用分子运动理论来求得流动特性。不同相之间的动量交换也依赖于混合物的类别。

3.2 流体流动的基本控制方程

流体流动受物理守恒定律的支配，要满足质量守恒方程、动量守恒方程。对于定常流动，密度为常数的情况下，质量守恒方程可表述为：

对于不可压缩常黏度的流体，在不考虑流体黏性的情况下，欧拉方程可表述为：

3.3 湍流模型

湍流是一种高度复杂的三维非稳态、带旋转的不规则流动。湍流中流动的各个物理量参数，如流速、压力、温度等都随时间和空间发生随机变化，特点是物理量的脉动，非稳态的N-S方程对湍流运动仍是适用的。湍流模型[1]分为标准k-ε湍流模型、RNGk-ε模型、Realizablek-ε模型3种，这些模型均是针对湍流发展非常充分的湍流流动来创建的。

3.4 计算模型的求解参数

计算模型的参数选取对于数学模型的计算结果至关重要。正确设置求解的控制方程、压力速度耦合方法、松弛因子、离散格式等参数，有利于提高模型计算的仿真精度。

1）求解的控制方程

在求解参数设置中，可选择的方程包括流动方程、湍流方程、能量方程、体积分数方程等。在求解过程中，有时为了得到收敛的解，等一些简单的方程收敛后，再开启复杂的方程一起计算。

2）压力速度耦合方法选择

在基于压力求解器中，一般有SIMPLE，SIMPLEC，PISO以及Coupled等4种压力速度耦合的方法。定常状态计算一般使用SIMPLE或SIMPLEC方法，但对于许多问题如果使用SIMPLEC可能会得到更好的结果。

3）亚松弛因子

亚松弛因子是基于压力求解器所使用的加速收敛参数，用于控制每个迭代步内所计算的场变量的更新。除耦合方程之外的所有方程，包括耦合隐式求解器中的非耦合方程（如湍流方程），均有与之相关的亚松弛因子。一般，压力、动量、k和ε的亚松弛因子分别为0.2，0.5，0.5和0.5。

4 计算模型建立

4.1 计算范围及计算网格

按照设计要求，模型计算水位选用正常蓄水位1735 m；库区水平范围以岩塞进口中心点为中心、按75 m爆破影响半径范围选取；计算假定在排沙洞水下岩塞爆破后，爆破岩塞所产生的岩渣瞬间全部堆满预挖的集渣坑，忽略岩塞岩渣和集渣坑的影响。

选取的计算库区尺寸为99 m×150 m×78.82 m(长×宽×高)；岩塞为圆台型，上口直径为22.3 m、下口直径为10 m、高为12.3 m，圆台轴线倾角为45°，岩塞进口中心点高程为1673.18 m；岩塞的下面接直径为10 m、高为3 m的圆柱；圆柱后面接一直径为10 m、中心角为51°的圆环；圆环下游接一坡度为6°的逆坡段圆柱，圆柱直径10 m，长度116.86 m；最后与逆坡段圆柱连接的是一段长度为1286.417 m、直径为10 m的顺坡圆柱，坡角为1.488°，圆柱出口高程为1631.6 m。参照黄河刘家峡洮河口排沙洞工程岩塞地质勘测工程地质勘查报告及剖面图（H-H，I-I），模型水库计算域分为上、中、下三层：上层为空气层，高度为5 m（高程从1735.0～1740.0 m）；中层为库水层，高度为41.5 m（高程从1693.5～1735.0 m）；下层为淤泥层，淤泥层高度为 32.32 m（高程 1661.18～1693.5 m）、含水量为 33.9%、孔隙比为1.114。

4.2 建立求解模型

1）模型求解器设置：①开启标准-湍流模型；②选择k-epsilon作为 Model；③在 k-epsilon Model下，保留其默认选择Standard不变；④保持Near-Wall Treatment下的默认设置Standard Wall Functions不变。

2）设置多相流模型，开启欧拉模型选项。

3）设置用于计算相间动量传递的拖曳规律，即将各相间相互影响的拖拽作用选择为schiller-naumann。

4）设置重力加速度。

4.3 离散格式

对流项可选择不同的离散格式：

1）一般，当使用基于压力求解器时，所有方程中的对流项均用一阶精度格式进行离散。

2）当使用基于密度求解器时，流动方程使用二阶精度格式，其他方程使用一阶精度格式进行离散。

3）此外，当使用分离式求解器时，可为压力选择插值方式。

4.4 流体材料

模型计算材料分为空气、水和沙三类，相应的计算参数见表1。

4.5 计算相

为了研究水下岩塞爆破后淤泥颗粒组成及堆渣含量对排沙洞中水沙运动规律的影响，该多相流模型的计算相大体分为空气、水和沙粒三类。针对工程特点并结合多相流模型的要求，模型中将空气作为主相，其他类均设为副相。由于淤泥层孔隙比为1.114，所以淤泥层中水的体积率为0.527。

表1 计算参数表

4.6 计算方案

按淤泥层中不同沙粒的组成，模型计算分4个方案，见表2。

表2 计算方案表

4.7 计算边界

边界条件的正确设置是多相流计算的关键。岩塞爆通后，岩塞口区淤泥将流向排沙洞，同时岩塞口区外的淤泥和水又会流向岩塞口区，使得岩塞口区的淤泥和水得到补充。

由于库区边界处各相的流速或流量均未知，为了便于模型计算，同时又能真实模拟水沙运动现象，特将模型水库计算域上游侧设置为压力进口边界，从水库由上到下依次为空气、库水和淤泥3种介质，按照上述水库计算域的3个分层，进口边界相应分为空气压力进口、库水压力进口和淤泥压力进口三部分；排沙洞出口设置为压力出口边界，压力按自由出流边界控制。

4.8 堆渣不冲流速

对岩块来说，最常用的不冲流速公式为伊兹巴士公式：

式中：V为不冲平均流速，m/s；D为粒径，m；g为重力加速度，m/s2；γs，γ 分别为岩块和水的容重；C 为反映岩块稳定状况的无量纲系数，取0.86时块石是处于滑动临界失稳状态，取1.2时块石是处于倾覆失稳临界状态；K一般为5～7 m1/2/s。

0.6m 堆渣，γs为 2.73 g/cm3，当C取 0.86～1.2 时，由（3）式计算的不冲流速为3.88～5.42 m/s。说明当流速大于3.88 m/s时，0.6 m堆渣处于滑动状态；当流速大于5.42 m/s时，0.6 m堆渣处于滚动状态。20 mm的沙粒，相应的不冲流速为0.71～0.99 m/s。

5 计算成果分析

本文采用Euler多相流模型，建立三维排沙洞水沙非恒定多相流计算模型，计算刘家峡排沙洞正常蓄水位、淤泥层不同颗粒组成等条件下，水下岩塞爆破后排沙洞中多相流运动过程，计算成果包括不同时间过程排沙洞内各相体积率及其流动速度、压力分布等。

1）对于不考虑岩塞口区淤泥层的情况，参照有压管道流速计算公式计算的岩塞爆破后排沙洞内水的运动速度为25.08 m/s。

2）方案1（3相流，不含堆渣)各相流速过程线见图1（a）。爆破后160 s，水相和淤泥相均流出排沙洞出口，淤泥（1 mm沙粒）相在逆坡段中心的流速为12.7 m/s，是排沙洞内的最小流速。

图1 排沙洞各相流速过程线图

3）方案2（3相流）各相流速过程线见图1（b）。爆破后175 s，水相和堆渣相均流出排沙洞出口，但在逆坡段中心处，堆渣（0.6 m堆渣）相的流速为5.13 m/s，是排沙洞内的最小流速，小于0.6 m堆渣5.42 m/s的不冲流速，存在0.6 m堆渣停止流动的可能。

4）方案3（4相流，不含堆渣）各相流速过程线见图1（c）。爆破后173 s，水相流至排沙洞全长97.66%的地方，淤泥相流至排沙洞全长97.41%的地方，堆渣（20 mm沙粒）相在逆坡段中心的最小流速为7.28 m/s远远大于20 mm沙粒0.99 m/s的不冲流速。

5）方案4（4相流，淤泥含0.6 m堆渣）各相流速过程线见图1（d）。岩塞爆破后198 s，在逆坡段中心，水相、堆渣(0.6 m堆渣)相及淤泥(1 mm沙粒)相均达到最小流速，依次为 5.38，5.08，5.16 m/s，均小于 0.6 m 堆渣 5.42 m/s的不冲流速，存在0.6 m堆渣在逆坡段停止运动的可能。

[1]计算流体动力学分析[M].北京：清华大学出版社，2004.