一种基于随机方法的粘性泥沙絮凝群体平衡数值模拟方法

1.本发明涉及泥沙絮凝的数值模拟技术领域，尤其涉及一种基于随机方法的粘性泥沙絮凝群体平衡数值模拟方法。


背景技术：

2.自然界的水陆域环境中的粘性泥沙由于其粒径较小，受水体中复杂的物理、化学及生物过程的影响，不断发生聚并与破碎过程，产生絮凝现象，主要以生物-矿物絮团的形式存在于水体中，进而改变粘性泥沙的粒径分布、沉降速度。对粘性泥沙絮凝机理的认识制约着航道疏浚和污染物输移等工程和环境问题的研究。
3.如何准确地预测絮团的粒径分布，是揭示其运动特性的关键难题，当前粘性泥沙絮凝模型往往集中在其矿物成分和物理性质上，对生物絮凝的影响和粒径分布的模拟缺乏认识，这种处理方法一定程度上制约了大尺度泥沙输运模型的预测精度与应用。
4.因此，为了更加充分的考虑粘性泥沙的絮凝特性，亟需建立一种耦合生物颗粒影响，包含絮团聚并、破碎过程的数值模型。


技术实现要素：

5.发明目的：本发明旨在提供一种提高模拟结果精确度的基于随机方法的粘性泥沙絮凝群体平衡数值模拟方法。
6.技术方案：本发明的基于随机方法的粘性泥沙絮凝群体平衡数值模拟方法，所述数值模拟方法包括以下步骤：
7.(1)生成拟随机数序列；
8.(2)根据水体环境条件、粘性泥沙性质和生物质颗粒性质计算颗粒聚并概率与破碎概率；
9.(3)根据过大絮团数进行聚并事件与破碎事件的选择；
10.(4a)若步骤(3)中选择了聚并事件，随机选絮团颗粒i和j，对其聚并概率采用接受-拒绝法检验，若检验通过则两颗粒发生聚并，若未通过则重新随机选择，直至聚并事件成功发生；
11.(4b)若步骤(3)中选择了破碎事件，随机选絮团颗粒i，对其破碎概率采用接受-拒绝法检验，若检验通过则该颗粒发生破碎，若未通过则重新随机选择，直至破碎事件成功发生；
12.(5)在步骤(4a)或(4b)后计算更新絮团颗粒性质，包括粒径、分形维数、聚并概率、破碎概率；
13.(6)计算推进模拟时间，模拟将重复步骤(3)～(6)直到完成计算时段；
14.(7)模拟在模拟时间满足计算时段后，计算包括颗粒粒径分布、特征粒径、生物组分的絮团性质。
15.进一步地，所述步骤(2)～(5)基于双变量群体平衡方程，双变量群体平衡方程为：
[0016][0017]
式中：v和vb分别为絮团体积和生物组分体积，n(v,vb,t)为拥有体积(v,vb)的絮团在t时刻的数量浓度；β(v,vb；v
′
,v
′b)和α(v,vb；v
′
,v
′b)分别表示体积分别为(v,vb)和(v
′
,v
′b)的两絮团的碰撞频率和聚并效率函数；a(v,vb)为体积为(v,vb)的絮团破碎概率函数；b(v,vb|v
′
,v
′b)为体积为(v
′
,v
′b)的絮团破碎成体积为(v,vb)絮团的分布函数。
[0018]
进一步地，所述步骤(2)中，水体环境条件包括紊动剪切率、水体密度、水体动力粘度、温度；粘性泥沙性质和生物质颗粒性质包括初始粒径分布、颗粒密度。
[0019]
进一步地，所述步骤(2)中，聚并概率为碰撞频率与聚并效率的乘积，碰撞频率考虑布朗运动、不等速沉降、湍流剪切作用以及生物粘性作用。
[0020]
更进一步地，碰撞频率计算公式为：
[0021][0022][0023][0024][0025][0026]
式中：β
ij
为絮团碰撞频率，β1和β2分别为粒径相关碰撞频率和生物组分相关碰撞频率，di和dj分别为絮团i和j的粒径，k为玻尔兹曼常数，t为绝对温度，μ为水体动力粘度，g为紊动剪切率，g为重力加速度，ωi和ωj分别为颗粒i和j的沉速，和分别为絮团i和j的生物体积组分，η
mb
为矿物-生物组分对碰撞频率的影响参数。
[0027]
进一步地，所述步骤(3)中，设定过大絮团粒径与柯尔莫哥洛夫微尺度成正比，考虑生物粘性对絮团最大粒径的增长影响；破碎事件发生概率为：
[0028][0029]
式中：nb为过大絮团的数目；当p
frag
大于所生成拟随机数，选择破碎事件发生，反之则选择聚并事件。
[0030]
进一步地，所述步骤(4a)中，若选中的两个絮团的聚并概率与所有可能的絮团对
中最大聚并概率的比值大于所生成拟随机数，则选中的两个絮团成功聚并，并在系统中随机复制一个絮团填补因聚并产生的颗粒空缺，在模拟过程中保持颗粒数目不变，反之则重新随机寻找絮团对直到事件成功发生；最大聚并概率的技术公式为：
[0031][0032]
式中：na和nr分别为聚并成功数与失败数。
[0033]
进一步地，所述步骤(4b)中，若选中的絮团的破碎概率与所有可能的絮团中最大破碎概率的比值大于所生成拟随机数，则选中的絮团成功破碎，产生的额外新絮团将随机替换系统中的等量絮团，在模拟过程中保持颗粒数目不变，反之则重新随机寻找絮团直到事件成功发生。
[0034]
进一步地，根据不同事件的发生率求解时间步长并推进，时间步长的计算公式为：
[0035][0036]
式中：δt为时间步长，《a
ij
》为所有颗粒对的平均聚并概率，c0为颗粒初始浓度，和分别为事件发生前后的系统平均质量，为初始系统平均质量，《b
ij
》为所有颗粒的平均破碎概率，κ为成功发生事件数。
[0037]
进一步地，所述步骤(7)中，若需要计算平衡态下絮团粒径特征，则需在达到拟平衡态后对结果进行平均。
[0038]
有益效果：与现有技术相比，本发明具有如下显著优点：显著提高模拟结果的精确性，立足于粘性泥沙絮凝机理，引入生物组分体积，将传统的絮凝群体平衡模型扩展为双变量群体平衡模型，并采用随机方法求解，实现对粘性泥沙-生物絮团性质的精确模拟，为深入研究矿物-生物絮凝机理，以及为粘性泥沙运动的精确模拟提供一种重要的技术手段，弥补了此前忽略生物组分的絮凝模型中存在的一些不足。
附图说明
[0039]
图1为本发明的流程图；
[0040]
图2为本发明的单变量模型验证结果；
[0041]
图3是本发明的双变量模型验证结果。
具体实施方式
[0042]
下面结合附图对本发明的技术方案作进一步说明。
[0043]
本发明在考虑粘性泥沙絮凝机理的基础上，构建矿物-生物粘性泥沙絮凝数学模型，模拟泥沙/生物颗粒间、絮团与流体环境间的固-固两相、固-液两相相互作用的颗粒絮凝过程，计算粘性泥沙絮团粒径分布、分形维数、生物组分等絮团性质随时间的演变；本发明的矿物-生物粘性泥沙絮凝数学模型，即双变量群体平衡方程如下：
[0044][0045]
式中：v和vb分别为絮团体积和生物组分体积，n(v,vb,t)为拥有体积(v,vb)的絮团在t时刻的数量浓度；β(v,vb；v
′
,v
′b)和α(v,vb；v
′
,v
′b)分别表示体积分别为(v,vb)和(v
′
,v
′b)的两絮团的碰撞频率和聚并效率函数；a(v,vb)为体积为(v,vb)的絮团破碎概率函数；b(v,vb|v
′
,v
′b)为体积为(v
′
,v
′b)的絮团破碎成体积为(v,vb)絮团的分布函数；方程右边的第一项第一项为由于聚并导致体积为(v,vb)的絮团生成的速率，第二项)的絮团生成的速率，第二项为由于聚并导致体积为(v,vb)的絮团减少的速率，第三项-n(v,vb,t)
·
a(v,vb)为由于破碎导致体积为(v,vb)的絮团减少的速率，第四项为由于破碎导致体积为(v,vb)的絮团生成的速率。
[0046]
具体的，如图1所示，本发明的基于随机方法的粘性泥沙絮凝群体平衡数值模拟方法，包括以下步骤：
[0047]
(1)生成拟随机数(低差异数)序列；
[0048]
(2)根据水体环境条件，包括紊动剪切率、水体密度ρw、水体动力粘度、温度；泥沙颗粒密度ρs、生物质颗粒性质，包括初始粒径分布，颗粒密度，计算颗粒聚并概率与破碎概率；
[0049]
其中，聚并概率为碰撞频率β
ij
与聚并效率的乘积，所述碰撞频率考虑布朗运动、不等速沉降、湍流剪切作用以及生物粘性作用，聚并效率一般取为常数1。
[0050]
具体的，碰撞频率函数β
ij
计算公式为：
[0051][0052][0053][0054][0055][0056]
式中：β
ij
为絮团碰撞频率，β1和β2分别为粒径相关碰撞频率和生物组分相关碰撞频率，di和dj分别为絮团i和j的粒径，k为玻尔兹曼常数，t为绝对温度，μ为水体动力粘度，g为紊动剪切率，g为重力加速度，ωi和ωj分别为颗粒i和j的沉速，和分别为絮团i和j的
生物体积组分，η
mb
为矿物-生物组分对碰撞频率的影响参数。
[0057]
其中，絮团的破碎概率函数采用其强度相关的计算公式，破碎概率函数ai的计算公式为：
[0058][0059]
式中：ai为破碎概率函数，c1为校正系数，μ为水体动力粘度，g为紊动剪切率，fy为絮团强度，nf为絮团分形维数，di为絮团i的粒径，d
p
为初始颗粒粒径。絮团发生破碎后，采用二元破碎分布，即絮团破碎成两个子絮团，质量比由拟随机数决定。
[0060][0061]
式中：df为絮团粒径，δ为拟合常数。
[0062]
(3)根据过大絮团数进行聚并事件与破碎事件的选择；
[0063]
过大絮团粒径与柯尔莫哥洛夫微尺度成正比，考虑生物粘性对絮团最大粒径的增长影响；破碎事件发生概率p
frag
的计算公式为：
[0064][0065]
式中：nb为过大絮团的数目。当p
frag
大于所生成拟随机数，选择破碎事件发生，反之则选择聚并事件。
[0066]
(4)采用接受-拒绝法检验絮团碰撞频率、聚并效率、破碎概率，模拟絮团聚并和破碎过程；
[0067]
(4a)若步骤(3)中选择了聚并事件，随机选絮团颗粒i和j，对其聚并概率采用接受-拒绝法检验，若检验通过则两颗粒发生聚并，若未通过则重新随机选择，直至聚并事件成功发生；
[0068]
具体的，聚并事件中，若选中的两个絮团的聚并概率与所有可能的絮团对中最大聚并概率的比值大于所生成拟随机数，则选中的两个絮团成功聚并，并在系统中随机复制一个絮团填补因聚并产生的颗粒空缺，在模拟过程中保持颗粒数目不变，反之则重新随机寻找絮团对直到事件成功发生。在计算最大概率时，可采用近似系数与平均聚并概率的乘积来近似最大聚并概率，以提高模型运算效率；
[0069]
最大聚并概率的计算公式为
[0070][0071]
式中：na和nr分别为聚并成功数与失败数。
[0072]
(4b)若步骤(3)中选择了破碎事件，随机选絮团颗粒i，对其破碎概率采用接受-拒绝法检验，若检验通过则该颗粒发生破碎，若未通过则重新随机选择，直至破碎事件成功发生；
[0073]
具体的，破碎事件中，若选中的絮团的破碎概率与所有可能的絮团中最大破碎概率的比值大于所生成拟随机数，则选中的絮团成功破碎，产生的额外新絮团将随机替换系统中的等量絮团，在模拟过程中保持颗粒数目不变，反之则重新随机寻找絮团直到事件成
功发生。
[0074]
(5)在步骤(4a)或(4b)后计算更新絮团颗粒性质，包括粒径、分形维数、聚并概率、破碎概率；
[0075]
(6)计算推进模拟时间，模拟将重复步骤(3)～(6)直到完成计算时段；
[0076]
具体的，每一个事件的模拟时间步长为隐式，即根据不同事件的发生率求解时间步长并推进，根据事件驱动蒙特卡罗方法计算时间，对聚并和破碎事件中时间步长δt为：
[0077][0078]
式中：δt为时间步长，《a
ij
》为所有颗粒对的平均聚并概率，c0为颗粒初始浓度，和分别为事件发生前后的系统平均质量，为初始系统平均质量，《b
ij
》为所有颗粒的平均破碎概率，κ为成功发生事件数。
[0079]
(7)模拟在模拟时间满足计算时段后，计算颗粒粒径分布、特征粒径、生物组分等絮团性质。
[0080]
其中，若需要计算平衡态下絮团粒径特征，则需在达到拟平衡态后对结果进行平均以减少随机方法带来的误差。
[0081]
模型验证：
[0082]
(1)单变量模型验证(不考虑生物颗粒体积组分)
[0083]
采用文献中的数据进行验证，试验中泥沙矿物单颗粒粒径为0.1～5μm，泥沙浓度为0.5g/l，水体紊动剪切率为g＝5,10,20,40s-1
。
[0084]
模型中初始颗粒粒径与试验中采用的泥沙矿物颗粒粒径分布相同，泥沙浓度、水体紊动剪切率均与试验条件一致，采用泥沙絮团粒径分布与中值粒径进行验证，所采用的相关参数值如表1所示。
[0085]
表1单变量模拟实例验证相关参数
[0086][0087]
验证的粒径分布和中值粒径结果如图2所示，可知本模型可较好的模拟单变量(纯
泥沙)情况下絮凝过程中絮团粒径分布、中值粒径的变化情况。
[0088]
(2)双变量模型验证(考虑生物颗粒体积组分)
[0089]
采用微藻-高岭土试验的数据进行验证，试验中采用的微藻为梅尼小环藻，高岭土矿物颗粒粒径为4μm，微藻细胞粒径为13μm，泥沙浓度为0.1g/l，微藻浓度为3
×
109cells/m3，水体紊动剪切率为g＝40,60,80s-1
。
[0090]
模型中采用的环境条件、颗粒性质与试验条件一致，过大絮团粒径与柯尔莫哥洛夫微尺度成正比，比例系数取3.1，采用泥沙絮团粒径分布与中值粒径进行验证，所采用的相关参数值如表2所示。
[0091]
表2双变量模拟实例验证相关参数
[0092][0093]
双变量模型验证的粒径分布和中值粒径结果如图3所示，可知本模型亦可较好地模拟双变量(泥沙 微藻)情况下絮凝过程中絮团粒径分布、中值粒径的变化情况。