基于LBM 水力裂缝中支撑剂颗粒团簇升阻力系数研究

张 涛，吴春燕，孙 堃，卢 聪，李 聪

(1.西南石油大学石油与天然气工程学院，成都 610000；2.空气动力学国家重点实验室，绵阳 621000)

石油工程中，油气藏储层水力压裂时压裂液携支撑剂在缝内输送是一个狭窄裂缝中的固液两相流过程[1-3]，此过程中支撑剂颗粒会聚集成团形成不同构型的颗粒团簇[4]，从而影响颗粒阻力的大小[5]。同时，裂缝壁面的存在会影响颗粒团簇受力，从而对颗粒的运移铺置产生影响[6]。因此，明确颗粒团簇在窄缝内运移过程中的受力规律，对水力压裂过程中支撑剂输送问题的研究有重要意义。

格子玻尔兹曼方法(LBM, Lattice Boltzmann Method)在介观尺度将流体离散为微团进行模拟计算。相对于物理实验和基于分子动力学的微观尺度模拟以及基于N-S 方程宏观尺度模拟，LBM 没有连续性假设且计算量小，对复杂流固边界附近的流场刻画较为精细[7]，容易获得颗粒在流场中的受力信息。多位学者采用该方法研究流固多相流动过程。BEETSTRA 等[5]基于LBM 推导不同几何形状颗粒群的曳力系数公式，其研究表明颗粒群中任意单颗粒的曳力系数低于孤立的单颗粒曳力系数。GONG 等[8]采用LBM 研究颗粒簇的阻力系数与颗粒雷诺数在沉降过程中的关系，随着颗粒雷诺数的增加，团簇中颗粒间相互作用减少。HILL等[9]基于LBM 研究得到曳力随孔隙率及雷诺数变化的关系式。YIN 等[10]基于LBM 对多分散相系下颗粒曳力系数进行了研究，并提出了曳力系数公式。SARKAR 等[11]模拟了液体中的颗粒曳力，扩大了固液多分散系相修正因子的应用范围。SOMMERFELD 等[12]应用LBM 确定不规则形状颗粒流动阻力系数的平均标准差随雷诺数增加而减小，表明颗粒形状的不规则性对较高的雷诺数有更强的影响。王春雨等[13]模拟了非均匀多孔颗粒的升力随达西数变化规律。杨佼[14]采用D2Q9模型，研究了双颗粒、四颗粒在不同颗粒间距、不同颗粒排布方式和不同雷诺数下的绕流模式与受力特性。

从上述文献看出，国内外学者基于LBM 对单颗粒及颗粒群的阻力问题进行了大量的研究，但关于壁面对颗粒团簇受力影响的研究较少。由于LBM 能实现对流固边界的精细刻画及对单个颗粒受力的计算，因此，本文基于LBM 研究颗粒团簇在窄缝壁面影响下的升力、阻力变化规律，为支撑剂输送研究提供理论依据。

1 LBM 数值模型

1.1 LBM 模型

格子玻尔兹曼模型包含三个要素[15]：流体粒子的离散速度集合、格子演化方程和平衡态分布函数。其中，离散速度集合eα选用三维19 个速度方向的D3Q19 速度模型，如下：

格子演化方程描述了具有离散速度的流体粒子分布函数fi(r,t) 在一个固定格子上不同t时刻的运动过程：

式中：r/m 为格子上节点坐标矢量；ci/(m/s)为速度； Δt/s 为时间步长； Δi为碰撞算子。

碰撞模型采用LADD 等[16]提出的双松弛模型，对一个立方体的网格来说：

式中：ωi为权重系数，无量纲；cs/(m/s)为理想气体方程中的声速；1 是单位矩阵。对D3Q19 模型，ω0=1/3，ω1-6=1/18，ω7-18=1/36，c2s=c2/3。

式中： ρ/(kg/m3)为流体密度；j/(kg/(m2·s))为动量密度；u/(m/s)为速度矢量。

由质量守恒和动量守恒得：

平衡态分布函数常用形式如下[18]：

1.2 边界条件

固体壁面边界采用半步长反弹格式[19]，固体边界格点不在流体网格点上，而是位于流体网格中间，即 (xf+xb)/2处，其表达式为：

LADD[20]对于固液两相流的移动边界，提出基于半步长反弹格式加入反映颗粒运动速度的项。假定流体质点位置r恰好位于颗粒表面的外部，并选定离散速度方向cb使r+cbΔt位于颗粒内部格点处，从而对移动边界进行处理。

流体在流体点和固体点连线的中点发生碰撞，将这些中点依次连接，曲折的边界线随着网格的细化逐渐接近颗粒的实际表面边界线，如图1 所示。

图1 不同位置球形颗粒边界Fig.1 Spherical particle boundary at different positions

1.3 流体-颗粒作用力

采用动量交换法计算颗粒受力[20]，颗粒表面的局部速度可以表示为：

局部速度由颗粒线速度U、角速度 Ω和颗粒质心r确定，其中rb=r+1/2ciΔt是颗粒边界点的位置。

移动边界处理导致在流体点和固体点之间产生局部动量交换，但颗粒和流体总的动量保持守恒。边界上流体对固体颗粒的作用力可以通过动量交换计算出来：

流体作用在颗粒上总的合力与合力矩通过对边界格点上的f(rb)与rb×f(rb)求和累积得到。

单个球形颗粒的绕流阻力与升力表示为：

式中：CD为曳力系数，无因次；CL为升力系数，无因次； ρ/(kg/m3)为流体的密度；U∞/(m/s)为来流的速度；A/m2为固体在来流方向的投影表面积；D/m 为颗粒直径；B/m 为壁面缝宽。

1.4 模型求解

采用LADD[20]教授开发的悬浮颗粒流开源软件Susp3D 进行问题求解。LBM 模拟需要将物理单位转换到格子单位[21-22]进行求解计算，且颗粒间相对位置固定以保证颗粒团簇形状参数不变，使用周期边界条件以有效地模拟颗粒团簇在无限静止流体中平行于壁面的等速运动。其中：Re=0.1、1 时，计算域为100×27×100；Re=10、100 时，计算域为100×62×100。模拟循环1500 个，共计计算150 000个格子时间步，从而获取颗粒曳力、升力。

计算采用Intel E5-2670 处理器，工作主频2.3 GHZ，24 线程并行。若在120 h 的计算时间后，颗粒受力未达到稳定状态，则通过调整颗粒速度、流体粘度、颗粒粒径的参数重新计算至收敛。单颗粒壁面影响下(颗粒中心到壁面距离与颗粒粒径比值H/d=1)曳力收敛曲线如图2 所示。

图2 曳力收敛曲线图Fig.2 Drag convergence curve

2 颗粒团簇物理模型

通过高清相机观测支撑剂运移过程中颗粒微观结构[23]如图3 所示，存在多个颗粒紧密连接聚集成团[24-25]的现象，形成不同构型的支撑剂颗粒团簇。

图3 颗粒团簇微观结构Fig.3 Microstructure of particle clusters

根据颗粒团簇运移迎流面的形状选取三种具有代表性的构型：星型、三棱柱型、长方体型，如图4 所示。颗粒团簇在无限静止流体中沿x轴正方向，平行于壁面做等速运动。

图4 颗粒团簇构型Fig.4 Configuration of particle clusters

3 模拟结果及分析

3.1 模型验证

不同雷诺数Re下星型颗粒团簇在无壁面流场中运移的曳力系数Cd与BEETSTRA等[5]的模拟结果、TRAN-CONG 等[26]的实验结果以及CLIFT 等[27]和STOKES[28]的曳力系数公式相吻合如图5 所示，从而验证了模型的准确性。对比发现，大雷诺数时模拟结果与文献存在一定的偏差，这是由于后颗粒受到沿来流方向前颗粒尾流的影响[29]以及颗粒间距的变化会干扰颗粒受力[30]。

图5 无壁面时团簇曳力系数对比Fig.5 Comparison of drag coefficients of particle cluster without wall

模拟与ZENG 等[31]相同的工况下(颗粒中心到壁面距离与颗粒粒径比值H/d=1)，有壁面影响时，单颗粒曳力系数Cd与Re关系如图6 所示，通过与ZENG 等[31]的数值模拟结果对比可以发现，最大误差小于1.8%，从而验证了在壁面影响下该模型的有效性。

图6 壁面影响下单颗粒曳力系数对比Fig.6 Comparison of drag coefficient of single particle with wall

3.2 壁面影响下颗粒团簇曳力与升力模拟

壁面会影响颗粒团簇受力，从而对颗粒沉降和运移产生干扰作用，在狭缝中这种现象更为明显。为研究此问题，在垂直于y方向增加一组平行壁面，如图4 所示，并保持缝间距为5 倍颗粒粒径。模拟壁面影响下，不同Re(颗粒团簇位于壁面中心)、不同颗粒团簇与壁面相对距离H/d时颗粒团簇的运移(H/d是最外侧颗粒中心到壁面距离与颗粒粒径比值，星型及三棱柱型取0.7、0.9、1.1、1.3、1.4，长方体型取0.7、0.9、1.1、1.5、1.9)，从而获取颗粒的曳力系数、升力系数。

3.2.1 星型团簇

星型团簇曳力系数CD与雷诺数Re的关系如图7 所示。与单颗粒曳力系数相同，CD随着Re增大，CD逐渐减小。

图7 不同雷诺数时星型团簇曳力系数Fig.7 Drag coefficient of star clusters at different Re

不同雷诺数Re时星型颗粒团簇各颗粒曳力系数与团簇平均曳力系数的比值CD/CD,s，如图8 所示。中心颗粒和前颗粒处于流场中心位置，与单颗粒曳力相似，其曳力主要受来流Re影响，CD随Re增大而减小，CD/CD,s随之减小。雷诺数较小时，后颗粒随Re的增大受壁面边界影响剧烈，CD/CD,s随之增大；雷诺数增大至100 时，壁面边界效应对后颗粒影响减弱，曳力主要受流场Re影响，CD/CD,s随Re增大而减小。上下左右(近壁面)四颗粒曳力受壁面影响，CD/CD,s随Re增大而增大。就各颗粒间的曳力大小而言，壁面影响下近壁面颗粒曳力占团簇总曳力比重较大，中心颗粒占比小，说明壁面不仅会加剧对外侧颗粒曳力受力的影响，且外侧颗粒对内部颗粒的受力有屏障作用。星型团簇各颗粒升力系数CL与雷诺数Re的关系，如图9 所示。随着Re的增大，CL总体呈下降趋势；壁面影响下，边界效应对上下左右(近壁面)颗粒影响作用明显，CL数量级较大；由于垂直于z轴方向未设置壁面，上下颗粒无壁面(z)方向CL总小于壁面(y)方向CL；其它颗粒位于流场中心，受到的升力小及对应数量级较小，CL基本保持不变。

图8 不同雷诺数时星型团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.8 Ratio of each particle drag coefficient of star cluster to average drag coefficient of cluster at different Re

图9 不同雷诺数时星型团簇各颗粒升力系数Fig.9 Lift coefficient of each particle of star cluster at different Re

星型颗粒团簇不同H/d时曳力系数CD，如图10所示。随着H/d的增大，壁面影响逐渐减小，颗粒团簇CD逐渐减小；在H/d=0.7～H/d=0.9 阶段，由于壁面边界层效应的影响较为明显，颗粒团簇在远离边界层的过程中CD下降较快；而在H/d=0.9～H/d=1.4 阶段，颗粒团簇从近壁面流场远离至中心流场处，在这一区域距离边界层较远，受壁面影响较小，CD下降程度较小。不同H/d时星型团簇各颗粒曳力系数与团簇平均曳力系数的比值CD/CD,s，如图11 所示。上下及右颗粒随H/d增大逐渐靠近壁面，边界效应增强，CD/CD,s随之增大；左(近壁面)颗粒随H/d增大逐渐远离壁面，边界效应CD/CD,s随之减小；H/d=1.4时，左右两颗粒处于流场对称位置，CD/CD,s相等；中心颗粒及前后颗粒，曳力主要受来流的流场影响，CD/CD,s随H/d增大保持平稳且占比较小，说明外侧颗粒对内部颗粒受力有屏蔽作用。

图10 不同H/d 时星型团簇曳力系数Fig.10 Drag coefficient of star clusters at different H/d

图11 不同H/d 时星型团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.11 Ratio of each particle drag coefficient of star cluster to average drag coefficient of cluster at different H/d

不同H/d时星型颗粒团簇各颗粒壁面方向升力系数CL，如图12 所示。团簇内除前后颗粒外所有颗粒在周围颗粒间屏蔽作用下，CL都随着H/d的变大而小幅度减小；前后侧颗粒由于没有其它颗粒的屏蔽作用，受壁面边界效应影响明显，CL大幅度减小。

图12 不同H/d 时星型团簇各颗粒壁面方向升力系数Fig.12 Lift coefficient of each particle wall direction of star cluster at different H/d

3.2.2 三棱柱型

三棱柱型颗粒团簇曳力系数CD与雷诺数Re关系，如图13 所示，壁面影响下CD随Re增大而减小。

图13 不同雷诺数时三棱柱团簇曳力系数Fig.13 Drag coefficient of triangular prism cluster under different Re

不同雷诺数Re时三棱柱颗粒团簇各颗粒曳力系数与团簇平均曳力系数比值CD/CD,s，如图14 所示。中心及前颗粒处于流场中心位置，壁面影响较弱，曳力主要与Re有关，CD/CD,s随Re增大而减小；后侧近壁面颗粒在边界效应影响下，CD/CD,s随Re增大而增大；后排中心颗粒，在小雷诺数时CD受壁面影响先增大，Re较大时CD主要受流场Re影响骤减，CD/CD,s随之先增大再减小；就各颗粒曳力大小而言，由于壁面影响下外侧颗粒对内部颗粒受力有屏障作用，后侧近壁面颗粒曳力比重较其他颗粒大。

图14 不同雷诺数时的三棱柱团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.14 Ratio of each particle drag coefficient of triangular prism cluster to average drag coefficient of cluster at different Re

三棱柱型颗粒团簇升力系数CL与雷诺数Re的关系，如图15 所示，前颗粒及后排中间颗粒位于流场中心，受到升力及相应CL量数级较小。其它颗粒位于团簇外侧，壁面(y)方向及无壁面(z)方向升力系数CL均随Re的增大而减小；且壁面方向CL总小于无壁面方向CL。

图15 不同雷诺数时三棱柱团簇升力系数Fig.15 Lift coefficient of triangular prism cluster under different Re

三棱柱型颗粒团簇不同H/d下曳力系数CD与雷诺数Re的关系，如图16 所示。壁面影响下三棱柱型颗粒团簇CD随着H/d的增大而减小；在H/d=0.7～H/d=0.9 阶段，颗粒团簇逐渐远离壁面，由于边界层效应的影响较为明显，CD下降较快；在H/d=0.9～H/d=1.4 阶段，颗粒团簇远离壁面至中心流场处，距边界层较远，受壁面影响较小，CD下降程度较小。

图16 不同H/d 时三棱柱团簇曳力系数Fig.16 Drag coefficient of triangular prism cluster at different H/d

三棱柱型颗粒团簇不同H/d时，各颗粒曳力系数与团簇平均曳力系数比值CD/CD,s，如图17 所示。随H/d的增大，后排近壁面颗粒逐渐远离壁面，边界效应减弱，CD/CD,s显著减小；后排远壁面颗粒逐渐靠近壁面，边界效应加强，CD/CD,s显著增大；H/d=1.4 时，颗粒团簇位于流场中间位置，后排近壁面与远离壁面颗粒处于流场对称位置，CD/CD,s相同；在外侧颗粒的屏蔽作用下，中间及后排中间颗粒CD/CD,s保持平稳；壁面边界效应会加剧颗粒受力，使后排近壁面颗粒曳力占团簇总曳力比重较大。三棱柱型颗粒团簇不同H/d下各颗粒壁面方向升力系数CL结果，如图18 所示。分析发现，不同H/d时，近壁面颗粒随H/d的增大逐渐远离边界层壁面方向CL减小，远壁面颗粒随着H/d的增大逐渐靠近壁面边界，壁面方向CL增大；壁面边界效应加剧颗粒受力，使近壁面颗粒CL总大于远壁面颗粒CL。

图17 不同H/d 时三棱柱团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.17 Ratio of each particle drag coefficient of triangular prism cluster to average drag coefficient of cluster at different H/d

图18 不同H/d 时三棱柱团簇各颗粒壁面方向升力系数Fig.18 Lift coefficient of each particle wall direction of triangular prism cluster at different H/d

3.2.3 长方体型

长方体型颗粒团簇曳力系数CD与雷诺数Re的关系，如图19 所示，与星型及三棱柱型颗粒团簇类似，壁面影响下长方体型颗粒团簇CD随Re的增大。

图19 不同雷诺数时长方体型团簇曳力系数Fig.19 Drag coefficient of cuboid clusters under different Re

不同雷诺数Re时长方体型颗粒团簇各颗粒曳力系数与团簇平均曳力系数比值CD/CD,s，如图20所示。可以发现，前排颗粒曳力受来流Re影响，CD/CD,s随着Re的增大而减小；中间颗粒CD/CD,s基本保持不变；后排颗粒在Re较小时曳力受壁面边界效应影响，CD/CD,s先增大，Re较大时对曳力的影响增强，CD/CD,s减小。就各颗粒间的曳力大小而言，中间颗粒受力在外侧颗粒屏蔽作用下，曳力系数占比较小。

图20 不同雷诺数时长方体团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.20 Ratio of each particle drag coefficient of cuboid cluster to average drag coefficient of cluster at different Re

不同雷诺数Re时长方体型颗粒团簇各颗粒升力系数CL，如图21 所示。通过分析发现，壁面影响下各颗粒壁面(y)方向和无壁面(z)方向CL均随着Re的增大而减小；且壁面影响下，颗粒的无壁面方向CL总体上小于壁面方向的结果。

图21 不同雷诺数时长方体型团簇各颗粒升力系数Fig.21 Lift coefficient of each particle of cuboid cluster at different Re

长方体型颗粒团簇不同H/d下曳力系数CD的结果，如图22 所示。随着H/d的增大，壁面影响减小，颗粒团簇CD逐渐减小；颗粒团簇在远离边界层的H/d=0.7～H/d=0.9 阶段，由于壁面边界层效应的影响较为明显，CD下降较快；在H/d=0.9～H/d=1.9阶段，颗粒团簇从靠近壁面流场处远离至中心流场，距离壁面边界层较远，受壁面影响较小，CD下降程度较小。

图22 不同H/d 时长方体团簇曳力系数Fig.22 Drag coefficient of cuboid clusters at different H/d

不同H/d时长方体型颗粒团簇各颗粒曳力系数与团簇平均曳力系数比值CD/CD,s，如图23 所示。随着H/d的增大，近壁面前中后三个颗粒逐渐远离壁面，壁面边界效应减弱，CD/CD,s逐渐减小；远离壁面前中后三颗粒逐渐靠近壁面，壁面边界效益增强，CD/CD,s随之增大；就各颗粒曳力大小而言，边界效应会加剧近壁面颗粒受力，使远壁面颗粒CD/CD,s较小；内侧颗粒受力在团簇外侧颗粒的屏蔽作用下，中间颗粒曳力系数占比较小。

图23 不同H/d 时长方体型团簇各颗粒曳力系数与团簇平均曳力系数比值Fig.23 Ratio of each particle drag coefficient of cuboid cluster to average drag coefficient of cluster at different H/d

长方体型颗粒团簇不同H/d下各颗粒升力系数CL结果，如图24 所示。团簇内各颗粒壁面方向CL均随H/d的变大而减小；近壁面颗粒的CL总大与远壁面颗粒；中间两颗粒壁面方向升力在外侧颗粒屏蔽作用下，随着H/d的变大基本保持不变。

图24 不同H/d 时长方体型团簇各颗粒壁面方向升力系数Fig.24 Lift coefficient of each particle wall direction of cuboid cluster at different H/d

4 结论

本文构建了星型、三棱柱型、长方体型颗粒团簇模型，并基于LBM 研究了窄缝中颗粒团簇在不同雷诺数Re、不同颗粒团簇与壁面相对距离H/d下的升阻力系数变化规律，得到以下结论：

(1)星型、三棱柱型、长方体型颗粒团簇的曳力系数CD随着Re、H/d的增大而减小；且颗粒团簇在H/d=0.7～0.9 阶段，距边界层较近，受壁面影响明显，CD下降较快；在颗粒团簇继续逐渐远离至中心流场阶段，距边界层较远，受壁面影响程度小，CD下降较慢。

(2)星型、三棱柱型、长方体型颗粒团簇，颗粒曳力系数占总团簇平均曳力系数比值CD/CD,s，随Re增大而减小及随壁面边界效应减弱而增大的程度不同，前排颗粒曳力系数比值减小、中间颗粒先增大后减小、近壁面颗粒逐渐增大。随着H/d的增大，近壁面颗粒远离壁面使边界效应减弱，比值随之减小；远壁面颗粒反之；中间颗粒位于流场中心，比值保持稳定。

(3)星型、三棱柱型、长方体型颗粒团簇，各颗粒升力系数CL随着Re、H/d的增大总体呈减小趋势；且壁面影响下随Re的增大，颗粒壁面方向CL均小于无壁面方向的结果。

(4)星型、三棱柱型、长方体型颗粒团簇，近壁面颗粒曳力系数与总团簇平均曳力系数比值CD/CD,s及升力系数CL，总大于远壁面颗粒的结果；且外侧颗粒受力总大于中间颗粒。壁面影响下，边界效应会加剧近壁面颗粒受力，且颗粒团簇外侧颗粒对内部颗粒受力有屏蔽作用。