杨 晨 ，陈 鑫 ，薛 亮 ，孙 健

（1. 华北水利水电大学 水利学院，河南 郑州 450046；2. 中国农业大学 水利与土木工程学院，北京 100083；3. 中国石油大学（北京） 石油工程学院，北京 102249；4. 清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084）

相对于欧美国家来看，我国河流含沙量普遍偏高，特别是北方地区，如黄河流域。另外，我国有广大的河口海岸地区如环渤海湾、江苏中部沿长江口至浙江北部都属于淤泥质海岸，悬移质含沙量比较高。由于生产活动的需要，我国许多水利与海洋工程都不可避免地需要直接应用高含沙水，例如使用高含沙水吹填岛礁或者用于农业灌溉，整治深水航道的疏浚维护等，即使某些应用会带来如增加灌溉成本、损害水力机械等很多危害，因此研究高含沙条件下的水沙运动具有重要意义。随着计算机软硬件技术的发展，工程中比较普遍使用ANSYS 软件模拟高含沙流动，以其中的Fluent 通用两相流模型的应用为例，涉及的模型控制方程和紊流模型众多，边界条件复杂，参数不确定性极大，本文对这些问题进行了探讨与研究。

1 控制方程

根据Elghobashi S E[1]的总结，水沙两相的耦合模式按照级别从低到高可以分为单向耦合、双向耦合和四向耦合三种类型。单向耦合模式下，泥沙颗粒对流体运动的影响被完全忽略，只需要描述颗粒在紊流中的对流扩散运动。当泥沙体积浓度达到10-6～10-3时，需要用双向耦合模式，此时考虑了颗粒相对流体相紊动的影响。特别是颗粒自身对周边水流绕流涡的影响，以及这种影响对颗粒运动的作用。当泥沙相浓度更高的时候，相间的相互作用需要考虑；泥沙颗粒之间的相互碰撞即碰撞效应应力变得非常重要，并且这种碰撞还强烈影响到自身以及流体的紊动特别是绕流涡，以及大量泥沙悬浮后对更多泥沙悬浮的再抑制的影响，这些因素都考虑的就是四向耦合模式。然而现有的Fluent 软件多相流模型却并非针对高含沙流动构建，故而需要在运用中针对此问题进行大量的二次开发，如层移输沙问题、岸滩底层泥沙运动问题和高含沙泥浆抽取及输送问题等。Fluent 软件包的通用两相流模型包括VOF 模型、混合模型和欧拉模型。VOF 模型所考虑的问题，两相之间并没有充分耦合，而是彼此之间有很明确的连续整段分界面分层或者自由表面流。混合模型跟VOF 方法一样仍然使用单一流体的方法，求解混合物各种物理量的方程，不过考虑相间的滑移，由于不考虑相间作用力和颗粒粒间应力，不适合于高含沙的情况。此处讨论的高含沙数值模拟是体积浓度达到10-3时的情形，适合于欧拉模型，是一种四向耦合模式。

紊流时均的相连续方程和动量方程，软件中略去了其中紊动的时间导数项、三阶紊动项、压力梯度紊动项、应力相关紊动项，得到如下控制方程组：

式中：下标 k 指代某一项，k=f，s；而 f 代表流体相；s 代表泥沙相；ρ 为相密度；t 代表时间；x 代表空间坐标轴；下标i,j=1，2 分别代表水平和垂直方向上的坐标，遵循爱因斯坦求和约定；u 为流速；p 为压力；τ 为应力张量；g 为单位质量力；F 为相间作用力，它由相界面上应力的积分所得。

运用Fluent 连续介质模型，需要强调的是，模型的建立基于双流体模型连续介质理论，只适用于高含沙浓度的情况，Hsu T 等[2]认为是泥沙体积浓度大于10-4。低浓度条件下的连续性假设不是客观成立的，在实际问题中强行将模型运用于低浓度问题不一定能获得收敛解，计算精度极有可能低于混合模型甚至是普通的扩散模型。此外，这个模型得到的只是两相的浓度和运动场，并没有方程表达单个粒子的运动，所以不能直接从此模型得到泥沙颗粒运动轨迹。

在一般浓度两相流问题中，相界面不连续且相间作用力太过复杂，很难对这些基本方程进行封闭求解。对模型各项的物理意义理解不一样而导致模型封闭不同，从而使得模型的最终表达式存在差异。Fluent 的紊动关联项中，两相体积分数与两相速度相关项可以根据梯度输移假定由式（3）给出，雷诺应力项可以根据Boussinesq 假设由式（4）给出。

式中：κ=νst/σ 为泥沙紊动扩散系数，σ 为泥沙Schmidt 数；νft和 νst分别为流体相和泥沙相紊动粘性系数。

2 基本应力项

Fluent 通用多项模型的基础是分子运动论，大量使用动理学理论描述泥沙相运动。Fluent 软件完整的固体动量方程式（2）包含了固体应力项τ，只考虑平移和碰撞从颗粒的动量交换中产生的剪切和体积应力，源于 Chapman S 等[3]、Syamlal M 等[4]、Gidaspow D 等[5]学者的理论。基础是分子运动论，多用动理学模型建模，主要运用于气固两相流和颗粒流。与颗粒温度有关，随颗粒温度增加而增大。颗粒温度物理意义与紊动能相同，与颗粒随机运动的动能成比例。应力粘性包括碰撞粘性（Collisional viscosity）、动力粘度（Kinetic viscosity）、体积粘度（Bulk viscosity）和摩擦粘度（Friction viscosity）。

这些源于分子运动论和动理学理论的应力项运用于气固两相流和颗粒流获得了较好的效果，但是在高含沙水流中使用仍然有待完善，因为即使紊动微弱甚至无紊动的高含沙流动中颗粒的碰撞和摩擦也是存在的。在河流海洋极限含沙区域，如底床区域，流动微弱，其适用性值得进一步商榷。此外体积粘度是颗粒压缩和扩张的抵抗力；摩擦粘度与可压缩机制下颗粒流动的固体压力正相关，而固体压力随颗粒温度或者是紊动能增加而增大。水沙两相流动中的假设是泥沙颗粒不可压缩，因此在高含沙的计算中不适合使用这两项。

Fluent 软件并不包含颗粒粒间应力 (Intergranular stress)，也叫离散应力(Dispersive stress)，这个是高含沙模拟的关键。为了获得更贴合实际的高精度模拟结果，可考虑在二次开发中将此项编译于源项中。Bagnold R A[6-7]最早对它进行了相关研究，其后还有 Ahilan R V 等[8]和 Fredsoe 等[9]，其等价粘性的一个基本形式如式（5）和式（6）：

式中：λ 是泥沙相的线性浓度；αsm是泥沙相最大体积浓度，一般取值是 0.6～0.65。Liu H 等[10]、Li M等[11]、Bakhtyar R 等[12]在高浓度水沙两相模型中充分使用并发展了粒间应力，Fluent 模拟中可考虑编译此项。

3 紊流模型

单相流紊流模型仅有速度紊动相关项，形式相对简单；而两相流动量方程由于包含体积分数、压力、速度的紊动相关项，使用的紊流模型远远要复杂。基于Fluent 可以使用或者二次开发的两相紊流模型跟单相紊流一样有很多种。在模拟高含沙水流时所有的常规紊流模型以及较为先进的大涡模拟和直接数值模拟都可以选择使用。由于紊流理论及先进的紊流模型本身也一直在发展中，综合考虑计算成本、稳定性、成熟性以及精度，其中最常用的是基于双方程紊流模型的模化。无论使用何种紊流模型，在高含沙状态下均需要合理考虑水沙的特性，特别是大量泥沙运动对水流的影响、颗粒碰撞的影响和颗粒绕流涡的影响。以经典常用的k-ε 模型为基础，Fluent 提供了3 种方法模拟多相流中的紊流。不过其既有两相紊流模型主要源于气固两相流动，有关固体相的紊动理论主要依靠分子运动论，在水沙两相流动中的适用性值得进一步完善。

最完整和容易理解的是每相紊流模型（Turbulence Model for Each Phase）。它为水相和泥沙相都求解一套k 和ε 输运方程。同时考虑体积分数、压力和流速的紊动后，两相流动方程按雷诺展开平均后推导所得到的紊动能方程和紊动能耗散率项数非常大[13]，分别具有三十多项和六十多项。当紊流传递在相间起重要作用时，这个紊流模型是合适的选择。然而此模型需要模化参数太多，计算成本巨大并且不易收敛，还有待进一步发展以提高结果的可靠性，周立行[14]进行了大量研究。考虑到两相最大速度差一般不会超过最大流速2%[11]，在高含沙水流动模拟中的应用需要慎重考虑。

Fluent 默认的是混合紊流模型（Mixture Turbulence Model）。它源于混合动量方程的时均化，仅求解混合物共同的k-ε 方程，一般认为它适合于相分离、分层或接近分层的多相流，或者两相之间的密度比接近1 的情况。不去考虑泥沙相跟水相紊动区别的情形下，使用混合属性和混合速度捕获紊流的重要特征是足够的。比较理想的运用情况是泥沙相的浓度不会太高，在一般的水力机械问题里面应当是合适的；而在河流海洋具有低床的情形，由于床面附近含沙浓度达到了极限状态，模拟效果会差很多。

最普通的是分散紊流模型（Dispersed Turbulence Model）。对水相使用k-ε 方程求解，在此基础上采用颗粒追随理论获得泥沙相紊动特征。这种情形下，对泥沙相随机运动的起支配作用的是水相紊流的影响，颗粒间的碰撞可忽略。此时泥沙相的紊动量根据水相的平均特征、颗粒弛豫时间和粒子相互作用时间的旋涡给出。当明显地水相是主要连续相而泥沙相是分散稀释的第二相时，这个模型是比较适用的。总体上看水力机械中很多高含沙水流都是属于这种情形，再考虑到此模型相对成熟和完善，计算成本不大，因此采用此模型用于水力机械高含沙水流计算是比较合适的。

关于紊流模型，还需要注意到高含沙浓度对水流紊动的抑制作用。Elghobashi S E 等[13]最初紊流模型包含体积紊动相关量，但并没有具体应用于高含沙问题，其参数也没有含沙影响。大量两相模型不直接计算高浓度处的紊动能量和耗散率，而像单相流动利用壁函数在底部一定高度给定高含沙条件下的紊动能和紊动能耗散率[12，15-16]。

4 边界条件

水相的边界条件中，最需要注意的是壁面边界条件。传统的壁函数表达式仅仅适合于清水状况，与泥沙（粒径、浓度）没有关系。当大量泥沙存在时，粗糙高度、粘性底层和过渡层的厚度都会发生较大的变化；流速分布跟清水条件下的分布也有很大不同，如图 1 所示[17]。其中：uf与式（1）一样代表水相流速；u*代表摩阻流速；y' 代表相对水深，y'=0 为水底，y'=1 为水面。而壁面边界条件的选择又极大影响紊流模型及动量方程的求解精度，然后影响泥沙扩散系数和浓度分布计算的精度，最终影响的是泥沙输送率的估计。关于这方面的论著已有不少[18-20]，可以考虑将最先进的成果用于Fluent 高含沙模拟中，改善数值结果。

泥沙相的边界条件，最大的问题是现有的泥沙运动理论体系都源于恒定流动，在考虑时间变化时并非总是适用的，极有可能给计算结果带来巨大误差。以周期性的流动为例，见图2 二阶Stokes 波流动过程中层移输沙问题泥沙参考浓度的变化过程：点代表实验值，点画线和虚线代表两个经典恒定流动公式，实线代表非恒定流动理论推导的理论解[21]，t 即式（1）的时间；T 代表二阶 Stokes 波流动的周期。传统恒定流动的边界条件不受相位差影响，在流动转向时候t/T～0.42 瞬时流速很小，根据经典的Engelund F 等[22]和 Zyserman J A 等[23]公式（即图 2中EF76 和ZF94）获得的参考浓度很小，而实际情况下已悬浮泥沙含量来不及沉降会保持一定的数值；不考虑大量已悬浮泥沙对再悬浮的抑制，两个经典公式的解随着流动速度增加一直增长到很大的水平t/T～0.21，与实验背离。扬沙率或者挟沙能力等泥沙参数的变化也跟参考浓度类似。由于现实中大部分流动都是时间变化的非恒定流动，因此可以考虑在运用Fluent 软件时重新编译泥沙边界条件。

5 结论

图1 清水和饱和含沙水流速分布

图2 非恒定流动中泥沙浓度边界条件

高含沙两相流动的精确数值模拟在我国水利工程中具有重要的意义。本文针对较为普遍使用的Fluent 商业软件通用两相模型，在控制方程、基本应力项、紊流模型和边界条件这几个要素的比选上展开了讨论。在此基础上本文提出了在软件模拟仿真中可能出现问题的原因及改进的方向，并期望Fluent 软件在高含沙流动模拟这方面能得到更好的发展和应用。