不同粒径组合块石群水下漂移数值模拟

李 小 超， 张 盖， 李 卫 国， 张 伟， 李 明 益， 常 留 红

( 1.长沙理工大学 水利工程学院， 湖南 长沙 410114；2.水沙科学与水灾害防治湖南省重点实验室， 湖南 长沙 410114；3.洞庭湖水环境治理与生态修复湖南省重点实验室， 湖南 长沙 410114；4.湖南省环境保护河湖污染控制工程技术中心， 湖南 长沙 410114；5.湖南省水运建设投资集团有限公司， 湖南 长沙 410011；6.中交天津航道局有限公司， 天津 300461 )

0 引 言

水上抛石施工作业时块石在水流作用下会发生漂移，块石落点位置不易确定，抛石不到位将造成石料浪费及作业效率降低．块石在动水中的漂距与抛投方式、块石自身特性、水流条件等因素有关，既有一定的规律性，又呈现出一定的随机性，目前很多学者针对单颗粒块石漂距提出了不同的经验公式[1-5]，但这些公式大多适用于较小水深和较小流速条件，随着水深和流速的增加，经验公式得到的块石漂距与实测值偏差逐渐增大[5]．此外，实际抛石施工中多采用抓斗、网兜等方式进行机械抛投，更大型的抛投方式还有底开驳船、侧抛船等船舶直接抛投，即实际工程中抛石为石料群抛，而非单颗粒抛石．对于群抛石料，下落过程中块石之间的碰撞增加了确定落点位置的难度．目前对于混合石料群抛漂移特性的研究还非常少，且大部分采用实验室水槽进行物理模型试验研究[6-7]．由于影响因素众多，难以用物理模型试验对各种影响因素条件下的块石漂移特性进行研究，且物理模型试验的时间周期和成本均较高，需要建立较为实用的数值模型用于群抛石料的漂移计算．

抛石水下漂移运动过程可以看作一种颗粒流现象，其中块石颗粒为固体相，水流为液体相．在颗粒流数值模拟方面，Sun等[8]将程序OpenFOAM和LAMMPS进行耦合，开发了用于求解固液两相流的SediFoam开源求解器，流场和颗粒运动的计算分别在OpenFOAM和LAMMPS中进行，其研究重点关注水流作用下的泥沙输移特性．张强强[9]将Fluent和EDEM进行耦合用于块石水中沉降分析，计算信息通过UDF和API接口在两个软件之间传递．刘卡等[10]采用MFIX-DEM软件对水下抛石运动进行了数值模拟．以上关于块石漂移的数值模拟研究主要针对单颗粒块石进行，且缺乏足够的试验数据验证．

航道整治工程中很多抛石工程面临工期紧、任务重、级配种类多、验收标准高等问题，迫切需要对混合石料群抛漂移特性开展系统研究，以指导精准抛投施工，提升水下抛石质量控制和成本控制．本文采用计算流体力学(CFD)与离散元方法(DEM)耦合的方法对块石群水下漂移进行计算，考虑到OpenFOAM和LAMMPS的开源性，颗粒粒径组合、亚格子尺度应力模型、块石碰撞模型等更易于实现，计算在SediFoam求解器中进行，重点对大水深、大流速条件下的石料群抛漂移特性及级配的影响进行探讨．

1 数值方法

1.1 水流运动方程

不可压缩流体的连续性方程和动量方程为

∇·(εsus+εfuf)=0

(1)

ρfεfg+Ffp)

(2)

式中：εs为固体体积分数；εf=1-εs，为流体体积分数；us为固体颗粒速度矢量；uf为流体速度矢量；p为压力；τ为流体黏性应力张量；ρf为流体密度；g为重力加速度矢量；Ffp为平均化的水流-颗粒相互作用力．

采用大涡模拟对流体方程进行求解，考虑块石漂移具有较高的雷诺数，亚格子尺度应力模型采用单方程亚格子涡黏模型[11]，其表达式为

(3)

1.2 颗粒运动方程

颗粒的平移和转动方程为

(4)

(5)

式中：m为颗粒质量；u为颗粒运动速度；fcol为相邻颗粒之间的接触力，块石之间的接触具有较强的非线性，采用考虑黏弹性影响的非线性弹簧-阻尼模型[12]计算；flub是由于流体从颗粒之间的间隙中挤出而产生的润滑力，采用Ball-Melrose模型[13]计算，该力与颗粒间的相对速度成正比，与颗粒间的距离成反比；ffp是水流-颗粒相互作用力，包括浮力fp、拖曳力fd、升力fl、附加质量力fa、Basset历史应力fB，即ffp=fp+fd+fl+fa+fB，拖曳力、升力、附加质量力和Basset历史应力的计算方法可分别参见文献[14]、[15]、[16]、[17]；I为颗粒转动惯性矩；ω为颗粒角速度；Mcol为接触力扭矩；Mlub为润滑力扭矩；Mfp为水流-颗粒相互作用力扭矩．

1.3 数值求解

计算域如图1所示，为一长方体区域，尺寸为18 m×10 m×15 m(长×宽×高)，原点设置于0号点处，水流沿X方向流动．采用有限体积法对计算域进行离散，网格划分采用结构化六面体网格进行，根据不同网格尺寸下的颗粒沉速和漂距试算结果，网格尺寸与颗粒直径的比值确定为2～3．河流的尺度非常大，若对整个区域进行计算，计算量将非常大．块石漂距最远只有十几米，与整个河流相比，只是非常微小的一个河段．对于此问题，本文数值模拟在X方向上设置为周期性边界条件，利用长度小很多(10～20 m)的河段代替实际河流，以减小模拟尺度，可在不失精确度的前提下，大大降低计算量．流速的垂向分布是影响块石漂移的主要因素，计算采用工程上常用的指数分布：

图1 数值水槽Fig.1 Numerical tank

u=um(Y/H)n

(6)

式中：um为水流表面流速，Y为距离床面的垂直高度，H为水深，n为指数．底部为无滑移边界，顶部为滑移边界．在Z方向上亦设置为周期性边界．

(7)

式中：r为颗粒半径，ρs为颗粒密度，G为剪切模量，vRay为Rayleigh波波速．上式确定的时间步长小于接触模型中弹簧-质量振子的半个周期．河道治理工程所用石料一般为石灰岩和花岗岩，块石质地较为均匀，根据文献[19]取块石密度为2 650 kg/m3，泊松比为0.3，弹性模量为3.755×104MPa，剪切模量为1.444 2×104MPa，Rayleigh 波波速为2 000 m/s．块石质量范围为1～300 kg，根据最大块石质量对应的等容粒径确定最大计算时间步长为2.467×10-6s．

2 结果与分析

2.1 数值模型验证

块石形状千变万化，很难对复杂的块石形状进行模拟，计算时将块石模拟为球体，采用原型块石进行，粒径为块石质量对应的等容粒径．水的密度为1 000 kg/m3，运动黏度系数为1×10-6m2/s．本文重点关注大水深、大流速下的块石群漂移特性，计算水深为15 m、表面流速为1.5 m/s，指数n为1/12．图2将计算流速与水槽模型试验测试数据进行了比较．图3为质量25 kg的单颗粒块石漂移位移随时间的变化，可以看出，块石在初始阶段有一个短暂的加速过程，之后由于流体与块石之间的相互作用，块石漂移速度基本保持稳定．漂距为块石初始下落点与块石顺水流下沉到底床着落点之间的X方向水平距离．图4给出了块石漂距数值模拟与模型试验结果[5]，两者吻合较好，数值模拟得到的漂距变化规律与模型试验结果整体上一致，均表现为随着块石质量(或块石粒径)的增大而逐渐减小．

图2 垂向流速分布Fig.2 Vertical flow velocity distribution

图3 质量25 kg的单颗粒块石漂移位移随时间的变化Fig.3 Drift displacement against time for single particle rock with mass of 25 kg

图4 单颗粒块石漂距Fig.4 Drift distance for single particle rock

文献[7]中作者在实验室水槽内开展了混合石料群抛试验，采用一次抛投的块石总体积来表征铲斗、网兜和抓斗等不同机械抛投方式．本文对抓斗群体抛石漂移进行计算，一次抛投块石总体积为5 m3．初始时刻块石群为一球体颗粒矩形空间阵列，最上层的颗粒紧贴水面，在X、Y、Z3个方向上颗粒等中心间距排列，间距大小根据试算确定，且不小于组合中颗粒最大直径．实际工程中块石质量范围比较大，但根据比尺(质量比尺为27 000)换算为模型后，其范围变得非常小．此外，在从碎石堆中选取块石时，实验人员会习惯性地参考之前已称取好的块石，使得选取到的块石质量大都比较接近，因此计算采用等粒径颗粒进行模拟，块石个数与试验时块石个数保持一致，块石群参数见表1，表中块石规格为实际工程中使用的抛投石料规格，排列方式为X、Y、Z3个方向上的块石排列个数．图5为块石群在底床上的落点分布形态，图中m*为块石平均质量．块石落至底床后没有出现堆积现象，且分布形状较为不规则．下落过程中水流紊动以及块石之间的碰撞使块石运动方向或速度大小发生改变，使得集中在一起的块石逐渐分散开．可以看出，在一次抛投块石总体积不变的情况下，平均质量越大，块石数量越少，块石落至底床后越分散．图6给出了相应的模型试验块石群落点分布照片，对比图5和图6可见，数值模拟和模型试验得到的块石群落点分布形态较为接近．

表1 等粒径块石群参数Tab.1 Parameters for rock group with equal particle diameter

(a) m*=33.129 kg

(a) m*=33.129 kg

漂距和块石群在底床上的落点分布面积是施工过程中比较关心的两个参数．图7为块石群平均漂距计算值与试验值的比较，平均漂距为块石群中所有块石漂距的平均值．数值模拟与模型试验结果较为一致，均随平均质量的增大而减小．块石群落点分布的形状较为不规则，密集区域与空白区域交错分布，分布面积存在一定的不确定性，本文采用如图8所示不规则图形的面积对块石群落点分布面积进行估算．选取块石群落点最外层的A、B、C、D4个点，分布面积为(ZC-ZA)×(XB-XD)/2，ZA、ZC、XB、XD为各点的坐标值，相当于图中虚线围成的三角形ABC和三角形ACD的面积之和．计算时主要考虑较为集中的范围，对于个别离分布集中区域较远的块石，计算不予考虑，例如图5(c)中右上角最远处的块石，计算时没有考虑．块石群落点分布面积计算结果见表2，表中模型试验值为多次抛投的平均值，偏差为相对偏差，即模拟计算值和试验值之差的绝对值与试验值的比值．在抛投总体积以及块石数量一定的情况下，落点分布面积随着平均质量的增大而减小．表中，除1～200 kg的块石计算值与试验值存在一定偏差外，其他计算值均与试验值吻合较好，偏差的产生与模型试验块石形状不规则等因素有关．

图7 块石群平均漂距Fig.7 Average drift distance for rock group

图8 块石群落点分布面积计算示意图Fig.8 Calculation diagram for rock group falling point distribution area

表2 块石群落点分布面积计算值与试验值的比较Tab.2 Comparison of rock group falling point distribution area between simulations and experiments

2.2 级配的影响

为了考察级配的影响，对4种不同规格的块石分别采用4种不同粒径组合进行模拟计算．对于每种块石规格，进行不同粒径组合时选取了5种块石，分别表示为块石1、块石2、块石3、块石4、块石5，其质量和粒径见表3，块石1～5质量是以平均质量(块石3的质量)为中心，相应地增加或减少相同的质量得到，以保证平均质量不变．对于不同块石规格，块石1、块石2、块石3、块石4、块石5所代表块石并不相同．计算时，对于同一规格块石，块石群的总体积、总质量、块石数量、块石平均质量及排列方式均保持不变．计算采用的粒径组合方式及各组合中采用的块石数量见表4，组合1、组合2、组合3、组合4分别为2、3、5、1种块石的组合，例如块石规格1～50 kg的组合1为100个块石1和100个块石5组成的块石群，组合4为上文讨论的等粒径组合．从组合1到组合4，其不均匀系数逐渐变小，见表5，即块石群级配组成逐渐变得更均匀；从连续性上看，从组合1到组合3，级配组成逐渐变得更连续，组合4由于缺少其他粒径的块石，其连续性最差．

表3 块石质量和粒径Tab.3 Mass and particle diameters of rocks

表4 不同粒径组合中块石数量Tab.4 Number of rocks for different particle diameter combinations

表5 不同块石组合的不均匀系数Tab.5 Non-uniformity coefficients for different rock combinations

计算均在水深15 m、表面流速1.5 m/s、指数为1/12的条件下进行．图9～12为4种不同规格的块石采用4种不同粒径组合计算得到的块石群在底床上的落点分布形态．从图中可以看出，不管是何种组合方式，块石群落点分布形态都呈现为不规则图形，随着组合方式的变化，分布形态的变化没有呈现出任何规律性．

表6为不同粒径组合方式下不同块石规格的块石群落至底床上的扩散带宽，水流向扩散带宽和横向扩散带宽分别为块石群落点在水流方向和垂直于水流方向上最大坐标和最小坐标之间的差值．对于相同规格的块石群，随着级配的变化，水流向扩散带宽和横向扩散带宽均没有出现单调的变化规律．

(a) 组合1

(a) 组合1

(a) 组合1

(a) 组合1

表6 扩散带宽Tab.6 Diffusion breadth

块石群平均漂距及在底床上的落点分布面积计算结果分别见图13、14．与水流向扩散带宽的变化规律类似，对于相同规格的块石群，从组合1变化到组合3，随着块石群级配组成的均匀性和连续性变好，平均漂距和分布面积逐渐增大．但从组合3变化到组合4，由于连续性变差，平均漂距和分布面积有所减小．对于不同的粒径组合方式，平均漂距随平均质量的变化趋势是相似的，均表现为随着平均质量(块石规格)的增大而减小．

图13 不同粒径组合的块石群平均漂距Fig.13 Average drift distance of rock groups with different particle diameter combinations

图14 不同粒径组合方式下块石群落点分布面积Fig.14 Falling point distribution area for rock group with different particle diameter combinations

3 结 语

本文基于CFD-DEM耦合方法，建立了用于计算抛石水下漂移的数值模型．模型采用大涡模拟方法对流场进行计算，对块石之间的接触力、块石与水流之间相互作用时的各种力进行了详细模拟，采用物理模型试验结果对数值模型的正确性进行了验证．对不同规格、不同粒径组合的块石群水下漂移进行了计算，分析了块石群在底床上的落点分布形态和扩散范围、漂距、分布面积等抛石施工过程中比较关心的块石漂移特性．结果表明，级配对块石群的漂移与扩散有较大的影响，块石群级配组成的均匀性和连续性越好，平均漂距和分布面积越大．抛石施工过程中抓取得到的块石粒径和质量分布具有很强的随机性，针对具体工程计算时，可参考所有块石(从石料场运来的块石)的级配构造块石组合．