重力影响下黏性细颗粒泥沙絮凝的三维模拟

柴朝晖1，李会云1,2，王茜1,2，刘同宦1,2

(1.长江科学院河流研究所,湖北　武汉　430010；2.水利部江湖治理与防洪重点实验室，湖北　武汉　430010)

【摘　要】　考虑重力作用是静水下黏性细颗粒泥沙絮凝与胶体絮凝主要不同点，因此研究重力对黏性细颗粒泥沙絮凝的影响是十分必要的。基于分形聚集生长理论，对黏性细颗粒泥沙的絮凝过程进行了三维数值模拟，分析了不同边界条件和含沙量下重力对絮团分形维数及粒径分布的影响，模拟结果表明：周期性和封闭性边界条件下分形维数的变化规律大致相同，即初始阶段分形维数的变化与布朗运动下的变化相似，而后分形维数先增加后减小，最终达到稳定值，且封闭性边界条件下的稳定值大于周期性边界下的稳定值；对于絮团粒径分布而言，考虑重力后，泥沙絮团粒径分布范围变广，异质性变强；并且随着泥沙浓度的增加，粒径较大的絮团比重变大。

【关键词】　黏性细颗粒泥沙；重力；三维模拟；边界条件；含沙量

基金项目：国家自然科学基金(51609012,51339001)；长江科学院中央级公益性科研院所基本科研业务费(CKSF2016010/HL)。

作者简介：柴朝晖(1986—　)，男，河北魏县人，高级工程师，博士，主要从事黏性泥沙运动特性研究。

E-mail：a3515522@126.com

1　引言

黏性细颗粒泥沙由于其表面电化学性质的作用，一定条件下，在水中会絮凝形成絮团或絮网。对于水流条件较弱的水库、湖泊等区域，重力作用引起的差速沉降是黏性细颗粒泥沙絮凝的主要作用力。泥沙絮凝改变了泥沙的沉速，同时使泥沙悬浊液由牛顿流体向宾汉流体转变，造成淤积物启动时具有成团的特性，导致泥沙大量淤积，影响水库和湖泊的蓄洪能力[1,2]；同时，细颗粒泥沙絮凝形成的网状结构改变了水库、湖泊底泥的水分构成，增加了内部封闭水的含量，给底泥的疏浚和处置带来了很大的困难[3]，并且，泥沙絮团物理吸附作用较强，是污染物的主要载体，影响污染物的迁移[4]，因此，研究重力作用对黏性细颗粒泥沙絮凝的影响具有重要意义。

沉降试验和图像分析法是絮凝研究中常用的方法[5,6]，但由于观测条件的限制，这两种方法无法观测絮团形成过程及其空间结构的变化。随着分形理论在各领域的广泛应用和计算机性能的提高，学者们开始通过计算机模拟黏性细颗粒泥沙的絮凝，如Kim等[7]基于扩散限制聚集模型(DLA)和扩散受限絮团聚集模型(DLCA)建立了细颗粒的差速絮凝模型，研究了絮团的形成和碰撞效率；杨铁笙等[8]基于DLCA模拟了黏性细颗粒泥沙静水絮凝沉降过程；张金凤等[9]基于Lattice Boltzmann方法模拟了不同粒径颗粒之间的碰撞和黏结。

上述模拟工作虽然考虑了由重力引起的差速絮凝，但未对由重力作用引起的絮团分形维数、粒径的变化等作深入的研究，而这些参数正有助于深入了解细颗粒泥沙的絮凝过程及絮团的性质，因此，本文通过三维数值模拟研究了不同浓度和边界条件下，重力作用对细颗粒泥沙静水絮凝形成的絮团分形维数的影响，同时利用多重分形分析了絮团粒径分布情况。

2　数值模拟方法

本文基于絮团分形聚集模型，模拟了黏性细颗粒泥沙的碰撞黏结过程。为了简化模拟过程，文中做以下假设：泥沙颗粒之间碰撞后会永久性结合；初始泥沙颗粒粒径是均一的[7-9]；模拟过程如下：泥沙颗粒(絮团)在布朗及重力作用下运动，经过单位时间步长运动后，判断颗粒(絮团)之间是否发生碰撞黏结，然后形成的新絮团作为一个整体进入下一次运动，重复这一过程直至模拟结束。这一过程通过在MATLAB平台上编程实现。

2.1　基本参数和边界条件

由于每次判断颗粒之间是否发生碰撞时，均需计算任意两个颗粒之间的距离，初始泥沙颗粒较多后，会形成一个巨大的矩阵，程序运行成本较高，同时考虑到泥沙絮团具有自相似性，能用絮团的某一部分反映整个絮团的性质，因此，本文选择在一个较小的区域[40μm×40μm×100μm(长×宽×高)]内进行模拟试验。根据《全国主要河流中水文特征统计》中河流含沙量的统计结果，初始含沙量分别取1.7kg/m3、4.3kg/m3和6.9kg/m3。考虑到初始泥沙颗粒的布朗运动，初始泥沙颗粒粒径取1.0μm，时间步长取1s。初始泥沙颗粒的分布位置由单独的子程序生成，每次运行主程序时，首先调用此数据文件，这样可确保每次模拟试验中初始泥沙分布位置是相同的，提高结果的可比性。

水平X和Y方向的均采用周期性边界条件，所谓周期性边界是指当泥沙颗粒运动出边界时，会从其相反界面重新进入。对于竖直Z方向的边界条件，本文采用周期性边界模拟远离水面和底面内的细颗粒泥沙絮凝过程，采用封闭性边界模拟近底面区域内细颗粒泥沙的絮凝，所谓封闭性边界是指泥沙絮团运动到底部时，会在底部发生沉积，由于泥沙浓度较小，并不考虑泥沙底部絮团的压缩变形。

2.2　颗粒运动方式及其计算

模拟中泥沙颗粒(絮团)的运动方式包括布朗运动和重力沉降，对于布朗运动，根据爱因斯坦布朗位移公式可知[10]：泥沙颗粒(絮团)在t时间内沿某一方向的平均位移:

(1)

式中：R为摩尔气体常数；NA为阿伏伽德罗常数，η为介质动力黏度；T为绝对温度；r为泥沙颗粒(絮团)的半径，r由絮团长、中和短轴的平均值确定。

由式(1)可得泥沙颗粒(絮团)在t时间内沿x、y、z方向的位移:

(2)

式中：φ,θ为球面坐标参数，模拟时由MATLAB随机函数随机产生。

对于泥沙颗粒(絮团)的重力沉降速度，由于絮团结构复杂、形状不规则、内部孔隙多，不宜采用传统的泥沙公式计算。在假设絮团孔隙中充满水和考虑絮团分形特性的前提下，笔者采用以下公式计算絮团沉速[11]：

(3)

式中：g为重力加速度；ρa为颗粒的质量密度；ρl为悬液的质量密度；d0为初始泥沙粒径；df为泥沙絮团粒径；DF为泥沙絮团分形维数。

3　研究参数

3.1　絮团分形维数

分形维数是分析絮体性质的重要参数，文中采用改进的回转半径法[12]计算絮团的分形维数DF，其基本原理为

N(l)∝lDF

(4)

式中：l为距离絮团重心的距离；N(l)为距离中心为l的空间内初始颗粒的数目。

通过计算N(l)与l的双对数坐标曲线的斜率可得到絮团分形维数DF。

3.2　粒径分布多重分形维

采用多重分形分析泥沙絮团粒径的分布情况，即用尺度为λ的“盒子”对整个絮团粒径范围进行划分，所需的盒子总数为N，由每个盒子的概率测度μi(λ)(该盒子内絮团的体积百分比)、尺度λ和q(文中取-5～5)，可得到广义分形维D(q)[13]为:

(5)

当q=0时，D(q)为粒径分布的容量维数D0，也就是集合维数，用于反映絮团粒径分布的最基本信息；当q=1时，D(q)为信息维D1,其主要提供粒径分布的异质性(均匀性)信息；当q=2时，D(q)为关联维数D2，其主要反映粒径分布的聚集程度。当D(q)随q的变化是一条直线时，则该分布不具有多重分形的性质，且为均匀分布[14]。

4　模拟结果与讨论

4.1　模型验证

为了验证此模拟方法的正确性，选用800颗初始泥沙颗粒，在布朗运动和周期性边界条件下进行模拟试验。根据Smoluchowski絮凝动力学理论可知：初始颗粒由于只能与其他颗粒或絮团碰撞黏结，其数量应一直减小；k级絮团(k为絮团中包含初始颗粒的数目)的数目则遵循先增加后减小的规律。图1是模拟过程中1级、2级、3级、5级和8级絮团随时间的变化情况，由图可知：初始颗粒一直在减少；不同级别的絮团出现的时间也不同，且其数目基本遵循先增加后减小的规律，符合Smoluchowski絮凝动力学理论。为了进一步验证模型，笔者计算了不同时刻的絮团分形维数，如图2，由于初始阶段絮团中颗粒的数量较小，且颗粒之间组合的随机性较强，絮团的分形维数波动较大，但当模拟时间超过200s后，絮团分形维数则稳定在1.85附近，此值与文献[15]～[17]中通过试验或模拟得到的分形维数相近，以上结果表明：此模拟方法可用于模拟细颗粒的絮凝过程。

4.2　分形维数的变化

使用含沙量为6.9kg/m3时的模拟结果分析重力作用对絮团分形维数的影响，分形维数为含有初始颗粒数大于10的絮团分形维数的平均值。图3和图4分别为周期性和封闭性边界条件下泥沙絮团分形维数随时间的变化情况，由图中可以看出：考虑重力作用后，无论在哪种边界条件下，其随时间的变化均与仅考虑布朗运动时有较大的差别，但两种边界条件下的变化规律大致相同，即初始阶段(0～140s)，分形维数的变化与布朗作用下的变化相似，而后分形维数先增加后减小，最终达到稳定值(周期性边界下为2.04，封闭性边界下为2.12)，且均大于布朗运动下的稳定值1.85。主要原因是：初始阶段絮凝形成的絮团粒径较小，重力沉降作用不明显，絮团与初始颗粒之间的差速絮凝作用下，分形维数的变化与仅考虑布朗运动时相差不大；随着时间的推移，絮团粒径逐渐变大，大絮团与颗粒或小絮团之间的相对沉速也随之变大，其在沉降过程中卷扫周围的小颗粒(絮团)，且这些小颗粒(絮团)易进入絮团的内部[18]，使絮团变密实，导致絮团分形维数增大；但当模拟区域内颗粒浓度降到一定程度后，大絮团在沉降过程中只能卷扫到其边缘附近的颗粒，絮团结构又变疏松，分形维数随之下降，这一过程中絮团结构形态的变化见图5。不同边界条件下，重力对分形维数的影响也有所不同，见图6。封闭性边界条件下分形维数开始增加的时间早于周期性边界条件下，且增幅较大，而开始减小的时间晚于周期性边界条件下，并且最终稳定时封闭性边界条件下的分形维数(2.12)大于周期性边界下的2.05。这主要是由于：沉积性边界条件下，絮团会在底部发生沉积，使模拟区域下部泥沙颗粒浓度较大，造成模拟区域中部的絮团沉降时能卷捕更多的小颗粒(絮团)，导致其分形维数开始增加的时刻较早，并且增幅较大；而后，模拟区域中上部泥沙颗粒浓度降低，絮团在沉降过程中很难卷扫到小颗粒(絮团)，分形维数在一段时间内保持不变，但当接近模拟区域底部时，又与底部沉积的絮团发生结合，使絮团分形维数降低；而周期性边界条件下则不存在絮团在模拟区域底部沉积的情况。

4.3　絮团粒径分布的变化

泥沙絮团粒径分布情况对于研究细颗粒泥沙沉降过程中清混面的形成机理有重要的意义[19,20]，因此笔者利用多重分形分析了不同边界条件和含沙量下泥沙絮团粒径分布的变化。

通过式(5)可得到絮团粒径分布的广义分形维D(q)，结果见表1(表中只给出了D0、D1和D2)。由表中的数据可知：无论哪种含沙量和边界条件，考虑重力作用后，D0、D1和D2均变大，也就意味着细颗粒泥沙絮凝形成的絮团与仅考虑布朗作用下相比，粒径分布范围更广、异质性较大，其主要原因是：考虑重力作用后，絮团与颗粒(絮团)之间发生差速絮凝，使部分粒径较大的絮团会卷扫周围的颗粒(絮团)，形成粒径更大的絮团，使絮团粒径分布较广，见图7。考虑重力作用后，含沙量对絮团粒径分布也有一定的影响。当含沙量从1.7kg/m3增至6.9kg/m3时，D0、D1和D2均增大，即随含沙量的增加，絮团粒径分布变宽，异质性变大，同时，D1/D0的值随含沙量的增加逐渐减小，如周期性边界条件下，含沙量为1.7kg/m3时D1/D0为0.969、含沙量为4.3kg/m3时D1/D0为0.922、含沙量为6.9kg/m3时D1/D0为0.814。由于D1/D0的值与粗颗粒含量成反比[21]，由此可知，随泥沙浓度的增加，大絮团含量增大。从表1中的数据还可发现：周期性边界条件下的D0、D1和D2均大于封闭性边界条件下的(6.9kg/m3下的数据D1和D2出现异常)，其原因是：周期性边界条件下，絮团周而复始地在模拟区域内运动，能卷扫更多的颗粒或絮团，然而由于絮凝的影响因素较多，具体机理还有待进一步研究。

5　结论

本文通过三维数值模拟研究了重力作用对絮团分形维数及粒径分布影响。最终得到以下初步认识：①周期性和封闭性边界条件下的絮团分形维数的变化规律大致相同，即初始阶段与布朗运动下的变化相似，而后分形维数先增加后减小，最终达到稳定值，且封闭性边界条件下的稳定值(2.12)大于周期性边界下2.04；②考虑重力作用后，泥沙絮团粒径分布范围变广，异质性变大；随泥沙浓度的增加，较大的絮团数目变多，并且周期性边界条件下的粒径分布更广。

参考文献：

[1]　王家生，陈立，刘林，等.阳离子浓度对泥沙沉速影响实验研究[J].水科学进展，2005,16(2):169-173.

[2]　洪国军，杨铁笙.黏性细颗粒泥沙絮凝及沉降三维模拟[J].水利学报，2006,37(2):172-177.

[3]　FU J J,CAI W M.Application of a well-designed cationic Polyelectrolyte for activated sludge dewatering [J].Journal of Chemical Engineering of Japan,2007,40(12):1113-1120.

[4]　陈明洪，方红卫，陈志和.泥沙颗粒表面磷分布吸附的实验研究[J].泥沙研究，2009,4:51-57.

[5]　LEE D G,BONNER J S,GARTON L,et al.Modeling coagulation kinetics incorporating fractal theories:comparison with observed data [J].Water Research,2002,36(4):1056-1066.

[6]　陈洪松,邵明安.NaCl对细颗粒泥沙静水絮凝沉降动力学模式的影响[J].水利学报，2002,8:63-67.

[7]　KIM A S,STOLZENBACH K D.Aggregate formation and collision efficiency in differential settling [J].Journal of Colloid and Interface Science,2004,271(1):110-119.

[8]　杨铁笙，李富根，梁朝皇.黏性细颗粒泥沙静水絮凝沉降生长的计算机模拟[J].泥沙研究，2005,4:14-19.

[9]　张金凤，张庆河.黏性泥沙不等速沉降絮凝的格子Boltzmann模拟[J].水利学报，2009,40(4):385-390.

[10]　胡纪华，杨兆禧，郑忠.胶体与界面化学[M].广州：华南理工大学出版社，1997:1-5.

[11]　柴朝晖，杨国录，陈萌，等.黏性细颗粒泥沙静水絮凝-沉降模拟[J].四川大学学报(工程科学版)，2012,44(1):48-53.

[12]　金鹏康，王晓昌，郭坤.絮凝体的DLA分形模拟及其分形维数的计算方法 [J].环境化学,2007,26(1):5-9.

[13]　MONTERO E.Rényi dimensions analysis of soil particle-size distributions [J].Ecological Modelling,2005,182(3-4):305-315.

[14]　DATHE A,TARQUIS A M,Perrier E.Multifractal analysis of the pore and solid phases in binary two-dimensional images of natural porous structures [J].Geoderma,2006,134(3-4):318-326.

[15]　LI X Y,PASSOW U,LOGAN B E.Fractal dimensions of small particle in Eastern Pacific Coastal Waters [J].Deep-Sea Research,1998,1(45):115-131.

[16]　LIN M Y,LINFDSY H M,WEITZ D A,et al.Universality in colloid aggregation [J].Nature,1989,339:360-362.

[17]　GARDNER K H,THEIS T L,YOUNG T C.Colloid aggregation:numerical solution and measurements [J].Colloids and surfaces A:Phy-sico-Chemical and Engineering Aspects,1998,141(2):137-252.

[18]　GONZLEZ A E.Colloidal Aggtrgation with sedimentation:computer simulations [J].Physical Review Letters,2001,86(7):1243-1246.

[19]　ALLAIN C,CLOITRE M,PARISSE F.Settling by Cluster Deposition in Aggregating Colloidal Suspensions [J].Journal of Colloid and Interface Science,1996,178(2):411-416.

[20]　WINTERWERP J C.On the flocculation and settling velocity of estuarine mud[J].Continental Shelf Research,2002,22(9):1339-1360.

[21]　POSADAS A N D,GIMÉNEZ D,BITTELLI M,et al.Multifractal Characterization of soil Particle-size Distributions [J].Soil Science Society of America Journal,2001,65(5):1361-136.

3-D Simulation of Flocculation of Cohesive Fine Sediment Induced by Gravitational Sedimentation

CHAI Zhaohui1,LI Huiyun1,2,WANG Xi1,2,LIU Tonghuan1,2

(1.River Department,Changjiang River Scientific Research Institute,Wuhan Hubei Province　430010;2.Key Laboratory of rivers and lakes harnessing and flood control,Wuhan Hubei Province　430010)

Abstract:The major difference between the flocculation cohesive fine sediment and colloid flocculation in still water is the gravitational sedimentation.Hence,flocculation process of cohesive fine sediment induced by gravitational sedimentation was simulated in 3-D space based on fractal cluster growing theory.The results indicate that changes of fractal dimension were similar to these of Brownian motion at early stage.With time went on,the fractal dimension increased and then decreased.Finally,it reached a stable value,which is 2.12 in closed boundary and 2.05 in periodic boundary.Moreover,when the gravitational sedimentation was considered,the range of floc size became larger and its heterogeneity was bigger and the ratio of large-sized flocs increased with the increase of sediment concentration.

Key words:Cohesive Fine Sediment;Gravity;3-D Simulation;Boundary Condition;Sediment Concentration