2. 河北大学岩土工程研究所, 河北 保定 071000
2. Institute of Geotechnical Engineering, Hebei University, Baoding 071000, Hebei, China
0 引言
管涌、潜蚀、河床冲刷、泥石流等诸多岩土工程问题中往往涉及颗粒-流体耦合，在对其中机理进行的研究中，近些年发展起来的颗粒尺度流固耦合方法正发挥着越来越重要的作用[1]。利用离散元(discrete element method, DEM)[2]建模岩土体，结合常用于流体流动建模的计算流体动力学(computational fluid dynamics, CFD)方法，Tsuji等[3]提出了一种颗粒尺度的流固耦合方法——CFD-DEM耦合方法，并成功将其用于模拟流体流动引起的球床流化问题。随后，Sharmy等[4-5]、Zhao等[6]、蒋明镜等[7]、洪勇[8]等以及周健等[9]论证了CFD-DEM耦合方法在岩土工程流固耦合问题模拟分析中的有效性。
CFD-DEM耦合模拟的准确性在很大程度上取决于能否合理对颗粒-流体相互作用建模。需要指出的是，按照流场求解的精细程度，目前的CFD-DEM耦合方法又细分为非解流(un-resolved)、全解流(fully-resolved)和半解流(semi-resolved)CFD-DEM耦合方法[10]，不同耦合方法中处理颗粒-流体相互作用的方式不同[10]。本文采用的是非解流CFD-DEM耦合方法，即通过引入颗粒拖曳力模型[11]来处理颗粒-流体相互作用，因此，拖曳力模型的精度很大程度上决定了CFD-DEM耦合模拟的准确性。
目前常用于验证拖曳力模型精度的基准算例包括水中单颗粒沉降、土柱中向上渗流和固结等。其中，水中单颗粒沉降涉及到水中单颗粒在自重和流体拖曳力作用下的沉降过程，该物理过程能够直接且简洁地反映单颗粒与流体的相互作用，因此常作为验证拖曳力模型的基准算例[7-12]。然而，当前关于拖曳力模型的验证，多是针对单一拖曳力模型，缺少不同拖曳力模型之间的对比，以3种典型拖曳力模型——Ergun、Wen和Yu模型[11]，Di Felice模型[13]及Hill和Koch[14]模型为例：蒋明镜等[7]，Guo等[15]对Ergun、Wen和Yu模型进行了验证，吴野[12]，Zhao等[16]，Hu等[17]对Di Felice模型进行了验证，Vaibhav等[18]，Wdfa等[19]对Hill和Koch模型进行了验证。已知在单颗粒沉降基准算例中，颗粒粒径越大，颗粒最终沉降速度越大，意味着颗粒雷诺数(Rep)越大，而且一般情况下，拖曳力模型的精度会随着颗粒雷诺数的变化而改变；因此，在基于单颗粒沉降基准算例验证拖曳力模型精度时，有必要考虑颗粒雷诺数的影响。
本文将基于水中单颗粒沉降的CFD-DEM耦合模拟，研究不同颗粒雷诺数时3种典型拖曳力模型的模拟精度，以期探讨颗粒雷诺数对其模拟精度的影响规律。
1 CFD-DEM耦合理论CFD-DEM耦合方法是一种从颗粒尺度求解颗粒-流体系统的耦合方法。在该方法中，通过离散元(DEM)建模固体颗粒，采用牛顿运动定律描述颗粒运动[1]，颗粒p的运动方程为
(1) 
(2) 式中：mp，up，ρp，Vp，ωp和Jp分别为颗粒p的质量、平移速度、密度、体积、角速度和转动惯量，粗细颗粒均适用，单位分别为kg, m/s, kg/m3, m3, m/s, kg·m2；t为时间，s；g为重力加速度，m/s2；
为了考虑颗粒-流体相互作用，在颗粒运动方程式(1)(2)中引入了流体对颗粒的作用力，主要包括拖曳力、压力梯度力、黏性力和浮力；其中，拖曳力是相间作用力的重要组成部分，模拟时由拖曳力模型计算获得，关于拖曳力模型的介绍将在第2节中给出。
流体建模时基于计算流体动力学(CFD)方法，采用局部平均化(locally averaged)Navier-Stokes(N-S)方程描述流体流动：
(3) 
(4) 式中：ε为网格孔隙率；▽表示对函数在各个正交方向上求导数以后再分别乘以各个方向上的单位向量; uf为流体速度，m/s；ρf为流体密度, kg/m3；p和τ分别为局部平均压力和剪应力, Pa；Ffp为颗粒作用于单位体积流体的相间作用力, N；fFD为虚体力的局部平均值, N。
在流体控制方程中，通过引入相间作用力项来考虑颗粒对流体的作用。
2 拖曳力模型介绍CFD-DEM耦合方法中采用拖曳力模型描述颗粒-流体相互作用，拖曳力模型的精度对于CFD-DEM耦合模拟的准确性来说至关重要。本文考虑了目前常用的3种拖曳力模型。
1) Ergun、Wen和Yu模型
基于Ergun公式[11]以及Wen和Yu公式[20]，结合计算密集颗粒堆积体和稀疏颗粒堆积体的两种公式，得到了Ergun、Wen和Yu模型[11]：
(5) 
(6) 式中：fp为流体对颗粒产生的拖曳力, N；dp为颗粒直径, m；μf为流体动力黏度, Pa·s；ufi为颗粒球心处的流体速度, m/s，通常采用颗粒球心所在流体网格中的流速；Cp为拖曳力系数；Rep为颗粒雷诺数，
(7) 2) Di Felice模型
Di Felice[13]通过对文献中颗粒堆积体内渗流试验数据进行整理和拟合，得到了Di Felice模型公式，此模型认为颗粒堆积体内流体作用于颗粒的拖曳力可以根据单颗粒在流体内平稳沉降时受到的拖曳力fp0乘以一个关于堆积体孔隙率ε的函数f(ε)来获得：
(8) 式中：
(9) 其中，Cp0为此模型的拖曳力系数，
(10) f(ε)为幂函数形式，
(11) 其中，
(12) 1) 和2)中的模型是都是基于渗流试验数据拟合所得。
3) Hill和Koch模型
另外一种常用的拖曳力模型是基于颗粒堆积体的完全求解方法获得的Hill和Koch模型[14]：
(13) 其中，
(14) 
(15) 本文采用CFD-DEM耦合方法建立了水中单颗粒沉降数值模型。其中，采用开源CFD软件OpenFOAM[21]建模流体(水)运动，通过开源离散元软件LIGGGHTS[22]模拟颗粒运动；同时采用CFDEMcoupling[23]处理CFD和DEM之间的数据交换，实现颗粒-流体耦合。
数值模型如图 1所示。CFD计算域为高0.750 m、边长0.225 m的长方体；颗粒初始位置(图 1黄色球体)位于CFD计算域中轴线上，且距离底面0.600 m。在CFD计算中，CFD计算域的网格划分情况如图 1a，c所示，其中，x和y向分别划分9个网格，z向划分30个网格，网格尺寸为25 mm。CFD模拟设置参数为：流体密度1 000 kg/m3，动力黏度10-3 Pa·s，时间步长0.001 s。DEM模拟设置参数为：颗粒密度2 650 kg/m3，杨氏模量2×1010 Pa，泊松比0.3，摩擦系数0.5，时间步长0.000 01 s。DEM程序每运行100个时间步后与CFD程序进行一次信息交换，以实现CFD-DEM耦合模拟。模拟时，颗粒从图示位置由静止开始下落，在重力、浮力以及流体拖曳力作用下发生自由沉降。
|
| a. 网格俯视图；b. 3D视图；c. 网格正视图。黄色球体代表颗粒。 图 1 数值模型示意图 Fig. 1 Numerical model |
|
|
采用Ergun、Wen和Yu拖曳力模型(式(5))，以粒径为1.0 mm颗粒的沉降为例，模拟得到了颗粒在水中的沉降过程。图 2为颗粒在水中沉降的速度-时间曲线，即沉降曲线。
|
| 颗粒粒径为1.0 mm，模拟中采用了Ergun、Wen和Yu拖曳力模型。 图 2 颗粒速度和拖曳力随时间变化曲线 Fig. 2 Variations of particle velocity and drag force with time |
|
|
由图 2可以看出：初始沉降时，颗粒速度随时间增加逐渐增大；随着时间继续增加，颗粒速度的增加率逐渐减小，最后趋于恒定值。这是由于随着颗粒的沉降，颗粒速度的增大先引起颗粒所受的拖曳力增加，导致颗粒拖曳力越来越接近颗粒浮重力(即重力减去浮力)；之后，颗粒拖曳力等于颗粒浮重力，颗粒受力平衡，沉降速度不再变化，此时的沉降速度称为颗粒最终沉降速度ut。
为了验证不同拖曳力模型的精度，将模拟得到的结果与试验数据对比。因Brown等[24]根据单颗粒沉降试验结果提出的预测单颗粒拖曳力经验公式在预测单颗粒与流体相互作用方面相对准确；因此，本文选取Brown和Lawler经验公式[24]作为标准，检验3种典型拖曳力模型在模拟单颗粒-流体相互作用时的精度。Brown和Lawler经验公式以颗粒拖曳力系数的形式给出：
(16) 考虑颗粒所受拖曳力
下文将根据模拟得到的颗粒沉降曲线以及最终沉降速度，验证和对比不同拖曳力模型的精度，并分析颗粒粒径对拖曳力模型精度的影响。
4.1 网格与颗粒尺寸比的敏感性在CFD-DEM耦合模拟中，网格与颗粒尺寸比(以下简称尺寸比)对模拟结果会产生重要影响，因此，在考察拖曳力模型精度之前，首先应分析尺寸比对模拟结果的影响。模拟中颗粒粒径为1.0 mm，考虑了5种不同的尺寸比(1.0，2.0，5.0，10.0，12.5)，并针对每一种尺寸比分别采用3种拖曳力模型进行模拟，得到了颗粒最终沉降速度ut随尺寸比的变化曲线(图 3)。
|
| 图 3 颗粒最终沉降速度随网格颗粒尺寸比变化曲线 Fig. 3 Variations of particle final sedimentation velocity with the ratio of mesh size to particle diameter |
|
|
由图 3可以看出：当尺寸比相对较小(< 2)时，ut与尺寸比具明显的相关性，即ut随着尺寸比的增加而增大，说明尺寸比确实会对模拟结果产生影响，但这在物理学上是不合理的；但当尺寸比≥5时，随着尺寸比的继续增加，ut不再发生明显变化，基于各个拖曳力模型得到的模拟结果趋于收敛。参考上述敏感性分析结果，下文模拟中采用的网格尺寸为25 mm，保证模拟中尺寸比均大于5.0。
4.2 3种经典拖曳力模型精度针对粒径为0.2 mm的颗粒分别采用3种拖曳力模型进行模拟，得到颗粒的沉降曲线，如图 4所示。由图 4可以看出，当采用不同拖曳力模型时，颗粒均先做加速度逐渐减小的加速运动，之后速度趋于恒定，即达到最终沉降速度ut。同时，将模拟结果与式(16)预测结果对比发现，采用3种拖曳力模型的模拟精度存在差异。其中：采用Hill和Koch模型模拟得到的ut明显大于经验公式的预测结果，采用Di Felice模型模拟得到的ut明显小于经验公式的预测结果，采用Ergun、Wen和Yu模型模拟得到的ut与Brown和Lawler公式的预测结果较为接近。可以看出，Ergun、Wen和Yu模型的精度相对较高。
|
| 图 4 颗粒速度随时间变化曲线 Fig. 4 Variations of particle velocity with time |
|
|
下文针对3种拖曳力模型，分析颗粒雷诺数对拖曳力模型精度的影响规律。需要说明的是，本节采用最终沉降速度ut(即将式(7)中的流体速度uf替换为ut)计算颗粒雷诺数Rep, 随着颗粒粒径的增大，颗粒最终沉降速度增加，颗粒雷诺数增大; 因此，可以通过变化颗粒粒径得到不同的颗粒雷诺数，进而用于考虑颗粒雷诺数对拖曳力模型精度的影响。本节考虑了13组不同的颗粒粒径(表 1)，模拟得到了每组颗粒粒径对应的最终沉降速度和颗粒雷诺数，如图 5所示。由图 5可见，ut和Rep均随着颗粒粒径的增大而增大。
| 粒径/mm | 网格尺寸/mm | 尺寸比 |
| 0.2 | 25 | 125.0 |
| 0.3 | 25 | 83.3 |
| 0.4 | 25 | 62.5 |
| 0.5 | 25 | 50.0 |
| 0.6 | 25 | 41.6 |
| 0.7 | 25 | 35.7 |
| 0.8 | 25 | 31.3 |
| 0.9 | 25 | 27.7 |
| 1.0 | 25 | 25.0 |
| 1.3 | 25 | 19.2 |
| 1.5 | 25 | 16.6 |
| 1.7 | 25 | 14.7 |
| 2.0 | 25 | 12.5 |
|
| 模拟中采用Ergun、Wen和Yu拖曳力模型。 图 5 最终沉降速度与颗粒雷诺数随粒径变化曲线 Fig. 5 Variations of particle final sedimentation velocity and particle Reynolds number with particle size |
|
|
为了分析Rep对拖曳力模型精度的影响，选取3个Rep值(5、26和150)，模拟得到了Rep在不同值时采用不同拖曳力模型的颗粒沉降曲线，并将模拟得到的ut与式(16)预测结果进行对比(图 6)。由图 6可知：当Rep = 5时，采用Ergun、Wen和Yu拖曳力模型模拟得到的ut与经验公式结果较为接近，模拟精度最高，Hill和Koch模型模拟得到的ut偏大，而Di Felice模型得到的模拟结果偏小；当Rep增大到26时，Hill和Koch模型模拟得到的ut与经验公式结果最为接近，Ergun、Wen和Yu模型和Di Felice模型得到的模拟结果均偏小；当Rep增大到150时，采用3种拖曳力模型模拟得到的ut均偏小，其中Ergun、Wen和Yu拖曳力模型的模拟精度相对较高。
|
| 图 6 不同颗粒雷诺数下各拖曳力模型的颗粒沉降曲线对比 Fig. 6 Sedimentation curves of particle with particle Reynolds number using different drag force models |
|
|
图 7中给出了采用3种拖曳力模型模拟得到的ut随Rep变化的曲线。由图 7可以看出: 当Rep≤14时，采用Ergun、Wen和Yu模型模拟得到的ut与经验公式预测结果最为接近，3种拖曳力模型的的精度由高到低为Ergun、Wen和Yu模型> Di Felice模型> Hill和Koch模型；当14 < Rep≤40时，Hill和Koch拖曳力模型的精度最高，Di Felice拖曳力模型的精度最低; 当40 < Rep≤72时，采用Ergun、Wen和Yu拖曳力模型模拟得到的ut与经验公式预测结果最为接近，采用Hill和Koch拖曳力模型所得结果次之，Di Felice拖曳力模型的精度最低；当Rep>72时，采用Ergun、Wen和Yu拖曳力模型模拟得到的ut与预测结果最为接近，采用Di Felice拖曳力模型所得结果次之，Hill和Koch拖曳力模型的精度最低。
|
| 图 7 各拖曳力模型下颗粒最终沉降速度随颗粒雷诺数变化曲线 Fig. 7 Variations of particle final sedimentation velocity with particle Reynolds number using different drag force models |
|
|
由拖曳力模型误差随Rep变化曲线(图 8)可以看出：随着Rep增大，采用Ergun、Wen和Yu模型及Di Felice模型模拟结果的误差也持续增大；而采用Hill和Koch模型模拟结果的误差在Rep小于26时先迅速减小，至Rep=26时达到最小值，接下来随着Rep增大，误差快速增大并超过其他2种模型。也就是说，Ergun、Wen和Yu模型以及Di Felice模型的精度均随着Rep的增大而降低，而Hill和Koch模型的精度随着Rep的增大出现先升高后减低的趋势。此外，当Rep > 150时，采用3种拖曳力模型模拟得到的最终沉降速度的误差均已经超过了15%。
|
| 图 8 最终沉降速度误差随颗粒雷诺数变化曲线 Fig. 8 Variations of particle final sedimentation velocity errors with particle Reynolds number |
|
|
1) Ergun、Wen和Yu模型以及Di Felice模型的精度均随着Rep的增大而降低，而Hill和Koch模型的精度随着Rep的增大出现先升高后减低的趋势，在Rep = 26时精度最高。当Rep > 150时，采用3种拖曳力模型模拟得到的最终沉降速度的误差均已经超过了15%。
2) 一般情况下，Rep≤14和Rep＞72时，3种拖曳力模型精度从大到小的顺序为Ergun、Wen和Yu模型> Di Felice模型> Hill和Koch模型，而当14 < Rep≤40时，Hill和Koch模型的精度最高，Di Felice模型的精度最低。
本研究主要基于相对简单的颗粒-流体耦合算例，即单颗粒沉降算例，初步探讨3种常用拖曳力模型的模拟精度及颗粒雷诺数对其模拟精度的影响规律。由于这3种拖曳力模型均是基于试验或全解流模拟试验数据拟合得出的经验公式，其精度受到试验条件、模拟精度、拟合误差等多种复杂因素的影响，其不同精度产生的原因尚有待进一步研究。
| [1] |
季顺迎. 计算颗粒力学及工程应用[M]. 北京: 科学出版社, 2018. Ji Shunying. Computational Granular Mechanics and Its Engineering Applications[M]. Beijing: Science Press, 2018. |
| [2] |
Cundall P A, Strack O D L. A Discrete Numerical Model for Granular Assemblies[J]. Géotechnique, 1979, 29(1): 47-65. DOI:10.1680/geot.1979.29.1.47 |
| [3] |
Tsuji Y, Kawaguchi T, Tanaka T. Discrete Particle Simulation of Two-Dimensional Fluidized Bed[J]. Powder Technology, 1993, 77(1): 79-87. DOI:10.1016/0032-5910(93)85010-7 |
| [4] |
El Shamy U, Zeghal M. Coupled Continum Discrete Model for Saturated Granular Soils[J]. Journal of Engineering Mechanics, 2005, 131(4): 413-426. DOI:10.1061/(ASCE)0733-9399(2005)131:4(413) |
| [5] |
El Shamy U, Zeghal M. A Micro-Mechanical Investigation of the Dynamic Response and Liquefaction of Saturated Granular Soils[J]. Soil Dynamics and Earthquake Engineering, 2007, 27: 712-729. DOI:10.1016/j.soildyn.2006.12.010 |
| [6] |
Zhao J, Shan T. Coupled CFD-DEM Simulation of Fluid-Particle Interaction in Geomechanics[J]. Powder Technology, 2013, 239: 248-258. DOI:10.1016/j.powtec.2013.02.003 |
| [7] |
蒋明镜, 张望城. 一种考虑流体状态方程的土体CFD-DEM耦合数值方法[J]. 岩土工程学报, 2013, 36(5): 793-801. Jiang Mingjing, Zhang Wangcheng. Coupled CFD-DEM Method for Soils Incorporating Equation of State for Liquid[J]. Journal of Geotechnical Engineering, 2013, 36(5): 793-801. |
| [8] |
洪勇, 李子睿, 唐少帅, 等. 平均粒径对砂土剪切特性的影响及细观机理[J]. 吉林大学学报(地球科学版), 2020, 50(6): 1814-1822. Hong Yong, Li Zirui, Tang Shaoshuai, et al. Effect of Average Particle Size on Shear Properties of Sand and Its Mesomechanical Analysis[J]. Journal of Jilin University (Earth Science Edition), 2020, 50(6): 1814-1822. |
| [9] |
周健, 周凯敏, 姚志雄, 等. 砂土管涌-滤层防治的离散元数值模拟[J]. 水利学报, 2010, 41(1): 17-24. Zhou Jian, Zhou Kaimin, Yao Zhixiong, et al. Numerical Simulation of Piping-Filter Prevention in Sandy Soil by Discrete Element Method[J]. Chinese Journal of Water Resources, 2010, 41(1): 17-24. |
| [10] |
程旷. 基于离散元的渗流致断级配土颗粒运移数值分析方法研究[D]. 大连: 大连理工大学, 2019. Cheng Kuang. Study on DEM-Based Numerical Methods for Seepage-Induced Fine Particle Migration in Gap-Graded Soils[D]. Dalian: Dalian University of Technology, 2019. |
| [11] |
Ergun S. Fluid Flow Through Packed Columns[J]. Chemical Engineering Progress, 1952, 48: 89-94. |
| [12] |
吴野. 天然岩土颗粒材料液体中拖曳力系数试验研究[D]. 大连: 大连理工大学, 2017. Wu Ye. Experimental Study on the Drag Force Coefficient of Natural Geotechnical Particle Materials in the Liquid[D]. Dalian: Dalian University of Technology, 2017. |
| [13] |
Di Felice R. The Voidage Function for Fluid-Particle Interaction Systems[J]. International Journal of Multiphase Flow, 1994, 20(1): 153-159. DOI:10.1016/0301-9322(94)90011-6 |
| [14] |
Hill R J, Koch D L, Ladd A J C. Moderate-Reynolds-Number Flows in Ordered and Random Arrays of Spheres[J]. Journal of Fluid Mechanics, 2001, 448: 243-278. DOI:10.1017/S0022112001005936 |
| [15] |
Guo Y, Yu X B. Comparison of the Implementation of Three Common Types of Coupled CFD-DEM Model for Simulating Soil Surface Erosion[J]. International Journal of Multiphase Flow, 2017, 91: 89-100. DOI:10.1016/j.ijmultiphaseflow.2017.01.006 |
| [16] |
Zhao J, Shan T. Coupled CFD-DEM Simulation of Fluid-Particle Interaction in Geomechanics[J]. Powder Technology, 2013, 239: 248-258. DOI:10.1016/j.powtec.2013.02.003 |
| [17] |
Hu Z, Zhang Y, Yang Z. Suffusion-Induced Deformation and Microstructural Change of Granular Soils: A Coupled CFD-DEM Study[J]. Acta Geotechnica, 2019, 14(1): 795-814. |
| [18] |
Vaibhav A, Yogesh S, hah M T, et al. Effect of Drag Models on CFD-DEM Predictions of Bubbling Fluidized Beds with Geldart D Particles[J]. Advanced Powder Technology, 2018, 29: S0921883118303212. |
| [19] |
Wdfa B, Jm A. CFD-DEM Solution Verification: Fixed-Bed Studies[J]. Powder Technology, 2018, 339: 760-764. DOI:10.1016/j.powtec.2018.08.044 |
| [20] |
Wen C Y, Yu Y H. A Generalized Method for Predicting the Minimum Fluidization Velocity[J]. AIChE Journal, 1966, 12(3): 610-612. DOI:10.1002/aic.690120343 |
| [21] |
Cheng K, Wang Y, Yang Q, et al. Determination of Microscopic Parameters of Quartz Sand Through Tri-Axial Test Using the Discrete Element Method[J]. Computers and Geotechnics, 2017, 92: 22-40. DOI:10.1016/j.compgeo.2017.07.017 |
| [22] |
Cheng K, Wang Y, Yang Q. A Semi-Resolved CFD-DEM Model for Seepage-Induced Fine Particle Migration in Gap-Graded Soils[J]. Computers and Geotechnics, 2018, 100: 30-51. DOI:10.1016/j.compgeo.2018.04.004 |
| [23] |
Kloss C, Goniva C, Hager A, et al. Models, Algorithms and Validation for Opensource DEM and CFD-DEM[J]. Progress in Computational Fluid Dynamics: An International Journal, 2012, 12(2/3): 140-152. DOI:10.1504/PCFD.2012.047457 |
| [24] |
Brown P P, Lawler D F. Sphere Drag and Settling Velocity Revisited[J]. Journal of Environmental Engineering, 2003, 129(3): 222-231. DOI:10.1061/(ASCE)0733-9372(2003)129:3(222) |