水流交汇区底泥影响下微塑料输运模拟方法

1.本发明属于水流与颗粒两相耦合运动模拟技术领域，具体涉及水流交汇区底泥影响下微塑料输运模拟方法。


背景技术：

2.微塑料是尺寸小于5mm的塑料碎片、纤维或颗粒，对于生态环境具有严重的威胁。微塑料颗粒体积小，比表面积大，疏水性强，颗粒本身具有吸附性，能够吸附有毒物质且难以降解，会伴随着物理行为在环境中迁移输运，最终危害整个生态系统。目前对于微塑料迁移输运问题的研究，大部分集中在海洋甚至外海中，主要探究海洋表层及深层微塑料的运移，近年来随着人们对微塑料的认识逐渐加深，学者们对于微塑料的研究方向从海洋转移至河口与内陆河流，相较于海洋环境，内陆河流尤其是河道中微塑料的迁移运输情况与我们的生活联系更为紧密，密布交错的河流水网是微塑料不断迁移污染的主要载体，因此对于内陆河流中微塑料输运过程的研究极为必要。
3.天然内陆河流中普遍存在水流交汇现象，在河道交汇处，水流紊动混掺剧烈，对微塑料的运动有着较强的干扰、驱动作用。交汇区流速、湍动能及结构的变化对微塑料输运、扩散有较大影响，因此本方法主要是针对交汇区微塑料输运规律的研究。其次，受水土流失及自然地理的影响，天然河道底部普遍有泥沙颗粒层的存在，水流中泥沙的运动特性显著影响微塑料颗粒的聚沉、再悬浮和渗透等行为，而微塑料颗粒与悬浮泥沙发生聚集会导致微塑料聚集体密度增加，从而影响微塑料在水中垂直分布和长期迁移。泥沙自身的运动特性(粒径分布、组成、启动等)、湍流的扰动时间与扰动强度等是影响河流中微塑料迁移的关键影响因素。
4.目前国内外学者对于交汇区、微塑料的输运规律的模拟分别作了大量研究，传统的交汇区水流的模拟方法基本以rngk-ε紊流模型为基础，相比于标准k-ε模型，rngk-ε紊流模型修正了湍动粘度，考虑了平均流动中的旋转及旋转波动情况，同时还在ε方程中增加一项，反映了主流时均应变率，因此rngk-ε模型可以更好的处理高应变率及流线弯曲程度较大的流动。采用vof方法模拟两相流，得到交汇区的水流特性。对于含有泥沙或污染物颗粒的水流，目前主要采用dpm、dem方法进行研究，dpm模型可以模拟水流中污染物的分布情况、沉降规律等。但dpm模型其存在着以下的问题：首先，dpm模型相比于dem模型只能设置颗粒的密度、粒径与形状因子，无法具体准确地设置微塑料及泥沙颗粒的物性、形状、颗粒的相互作用系数、剪切/弹性模量等；其次dpm模型只适用于颗粒相体积分数小于10％的情况，同时不考虑颗粒体积以及颗粒和颗粒之间相互作用力的情况，不适用与内陆污染程度大、微塑料赋存含量高的河流。


技术实现要素：

5.本发明的目的是提供水流交汇区底泥影响下微塑料输运模拟方法，通过准确设置微塑料及泥沙颗粒的性质，耦合模拟微塑料颗粒在含有底泥的交汇区水流中的运动，得到
水流交汇区底泥影响下微塑料颗粒的输运规律与特性，为后续处理河流中微塑料污染的问题提供相关的参考依据。
6.本发明所采用的技术方案是，水流交汇区底泥影响下微塑料输运模拟方法，具体按照以下步骤实施：
7.步骤1、根据内陆河水流交汇区泥沙颗粒分布情况，查询泥沙颗粒的物理性质以及泥沙颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域泥沙颗粒的形状，构建泥沙颗粒模型；
8.步骤2、根据内陆河流水流交汇区微塑料颗粒分布情况，查询微塑料颗粒的物理性质以及微塑料颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域微塑料颗粒的形状，构建微塑料颗粒模型；
9.步骤3、根据干流、直流尺寸数据建立内陆河水流交汇区模型，对内陆河水流交汇区模型进行网格划分，根据设定的最小网格尺寸生成相应的网格数量；
10.步骤4、设置泥沙颗粒、微塑料颗粒相相关计算参数，建立静态泥沙颗粒工厂，设定动态微塑料颗粒入射参数；
11.步骤5、根据干支流的来流情况计算流体域边界条件；
12.步骤6、基于内陆河水流交汇区模型与泥沙颗粒模型、微塑料颗粒模型列出流体域与颗粒相的控制方程组，将流体域与颗粒相进行耦合计算，在流体域计算中，将流体域所有六面体单元的控制方程组联立，采用有限体积法建立控制方程组离散代数方程组，将流体域边界条件代入，并将全流体域压强值、流速值初始化，流体域使用三维瞬态方法计算；在颗粒相的计算中，每个颗粒都包含6个自由度，两种运动方式：平动和旋转，在计算时，运用牛顿第二定律计算平动加速度和旋转加速度，在一个时间步上对他们进行数值积分，以此更新每个粒子的速度和位置，借助计算机进行数值计算求解，每隔一定时间保存迭代的计算数据，通常流体域计算数据保存时间为颗粒相计算数据保存时间的100倍；
13.步骤7、将保存的计算结果输入ensight软件进行后处理，得到流场、微塑料颗粒及泥沙颗粒运动参数。
14.泥沙颗粒的物理性质、微塑料颗粒的物理性质包括密度、泊松比、弹性模量、剪切模量、恢复系数、静摩擦系数和动摩擦系数。
15.步骤4中相关计算参数包括泥沙颗粒的物理性质、材料、分布类别，微塑料的物理性质、材料、分布类别。
16.建立静态泥沙颗粒工厂包括查阅相关研究文献确定底泥覆盖层厚度、在底泥覆盖区添加静态泥沙颗粒；动态微塑料颗粒入射参数包括微塑料颗粒入射类型、微塑料颗粒入射位置及微塑料颗粒入射颗粒数目。
17.流体域的边界条件包括干流速度进口边界条件、支流速度进口边界条件、出口边界条件、以及河流边界条件，河流边界条件定义为粗糙度无滑移边界，粗糙度根据实际河流情况决定，采用标准壁面函数法进行计算，出口边界条件为压力出口，干流速度进口边界条件、支流进口边界条件均包括进口流速，其中：
18.进口流速计算式为：
19.20.式(1)中，u
进口
为进口流速，m/s；q
进口
为进口流量，m3/s；b为进口断面河宽，m；h
进口
为进口断面平均水深，m。
21.流体域使用三维瞬态方法计算的方程如下：
22.质量守恒方程：
[0023][0024]
式(2)中，v表示控制体，控制体为流场单个非结构网格的体积；a表示控制面；ρ为流体密度；v为流体流速；n为控制面数；等式左边第一项表示控制体v内部质量的增量，第二项表示通过控制表面流入控制体的净通量；
[0025]
动量守恒方程：
[0026][0027]
式(3)中：f
bx
，f
by
，f
bz
分别是单位质量流体上的质量力在三个方向上的分量；p
xy
，p
yz
，p
xz
等是该流体内应力张量的分量；
[0028]
为保证赋予的重力加速度方向与实际情况保持一致，重力加速度方向依据网格本身坐标系位置确定，湍流模型使用rngk-ε：
[0029]
k方程：
[0030][0031]
ε方程：
[0032][0033]
式(4)与式(5)中，gk表示由于平均速度变化梯度引起的湍动能产生；gb是由于浮力影响引起的湍动能产生；ym为可压缩湍流脉动膨胀对总耗散率的影响；μ为湍流粘性系数。
[0034]
步骤6中微塑料颗粒的控制方程如下：
[0035]
旋转运动：
[0036][0037]
式(6)中，i是转动惯性矩，单位为kg
·
m2；ω是角速度，单位为rad/s；m是作用在质量上的总扭矩，单位为n
·
m；t是时间，单位为s；
[0038]
平移运动：
[0039]
[0040]
式(7)中，v是颗粒的平移速度，单位为m/s；m是颗粒的重量，单位为kg；fg是作用在颗粒上的联合重力，单位为n；fc和f
nc
是颗粒与周围颗粒或壁面之间的接触和非接触力，单位为n；
[0041]
数值积分：
[0042][0043]
式(8)中，v(t)是速度，单位为m/s；x(t)是位置，单位为m；a(t)是颗粒在给定步长的加速度，单位为m/s2；
△
t是时间步长，单位为s；同理更新颗粒的角速度和方向。
[0044]
步骤6中流体域与颗粒相耦合计算的具体过程为，流体域计算完一个时间步长后,将计算的最新计算结果传递给泥沙颗粒、微塑料颗粒相，计算相互间作用力，将相互间作用力引入微塑料颗粒的控制方程计算微塑料颗粒速度与位置，最后将计算体积分数连同相互作用力传递回内陆河水流交汇区模型计算，从流体域开始计算到最终传递回流体域，这是一个循环过程。
[0045]
步骤6中采用有限体积法建立控制方程组离散代数方程组具体过程包括将每个六面体网格中心节点的控制方程在该六面体网格进行积分，积分得到六面体网格的代数方程组，该方程组中每个六面体网格边界上的参数值由二阶迎风方程得到。
[0046]
步骤7中微塑料颗粒及泥沙颗粒运动参数包括速度、压强、湍动能、颗粒位置、颗粒留存时间。
[0047]
本发明的有益效果是：
[0048]
本发明中水流交汇区底泥影响下微塑料输运模拟方法，对泥沙颗粒和微塑料颗粒的密度、泊松比、弹性模量、剪切模量、恢复系数、静摩擦系数、动摩擦系数、相互作用进行设定，使建立的颗粒与实际颗粒更接近，并将颗粒相与流场作用耦合，考虑颗粒相体积分数，模拟了底泥影响下不同流速、不同性质的三种微塑料颗粒在交汇区附近的流动，有利于探究底泥影响下微塑料颗粒的输运规律。
附图说明
[0049]
图1是本发明水流交汇区底泥影响下微塑料输运模拟方法流程图；
[0050]
图2是fluent求解流程图；
[0051]
图3是edem软件结构框架图；
[0052]
图4是fluent耦合edem流程图；
[0053]
图5是建立的不同泥沙颗粒形状图；
[0054]
图6是建立的不同微塑料颗粒形状图；
[0055]
图7是颗粒相及流场网格图；
[0056]
图8是流场说明图；
[0057]
图9是单体泥沙颗粒图；
[0058]
图10是单体微塑料颗粒图；
[0059]
图11是本工况下各时刻微塑料颗粒分布图；
[0060]
图12是本工况下微塑料颗粒悬浮或沉降淤积示意图；
[0061]
图13是本工况下微塑料颗粒在渠道中驻留时间示意图。
具体实施方式
[0062]
下面结合附图及具体实施方式对本发明进行详细说明。
[0063]
如图1所示，本发明水流交汇区底泥影响下微塑料输运模拟方法，具体按照以下步骤实施：
[0064]
步骤1、根据内陆河水流交汇区泥沙颗粒分布情况，查询泥沙颗粒的物理性质以及泥沙颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域泥沙颗粒的形状，构建泥沙颗粒模型；
[0065]
步骤2、根据内陆河流水流交汇区微塑料颗粒分布情况，查询微塑料颗粒的物理性质以及微塑料颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域微塑料颗粒的形状，构建微塑料颗粒模型；
[0066]
泥沙颗粒的物理性质、微塑料颗粒的物理性质包括密度、泊松比、弹性模量、剪切模量、恢复系数、静摩擦系数和动摩擦系数。
[0067]
步骤3、根据干流、直流尺寸数据建立内陆河水流交汇区模型，对内陆河水流交汇区模型进行网格划分，根据设定的最小网格尺寸生成相应的网格数量；
[0068]
步骤4、设置泥沙颗粒、微塑料颗粒相相关计算参数，建立静态泥沙颗粒工厂，设定动态微塑料颗粒入射参数；其中，相关计算参数包括泥沙颗粒的物理性质、材料、分布类别，微塑料的物理性质、材料、分布类别；具体包括：
[0069]
步骤4.1、根据步骤1得到的泥沙颗粒与微塑料颗粒的物理性质，定义颗粒相中不同形状泥沙颗粒、不同形状且不同类别的微塑料颗粒对应的物理性质；添加微塑料颗粒、泥沙颗粒、河流边界三者间的相互作用的恢复系数、静摩擦系数、动摩擦系数；
[0070]
步骤4.2、对步骤3得到的六面体网格的各特征面进行命名。
[0071]
步骤4.3、设置颗粒工厂，包括静态颗粒工厂及动态颗粒工厂。
[0072]
静态颗粒为底泥覆盖区域的泥沙颗粒，查阅相关文献确定底泥覆盖层厚度，在底泥覆盖区添加步骤4.1所设置的泥静态沙颗粒，动态颗粒为微塑料颗粒，选择步骤4.2的干流与支流入口为微塑料颗粒入射面，入射的微塑料颗粒为步骤4.1中所设置的微塑料颗粒，选择微塑料颗粒注入的模式，确定每秒入射的颗粒数量，微塑料颗粒入射方向垂直虚体面，给定微塑料颗粒一定的入射速度，入射位置于整个虚体面随机分布；因涉及不同材料及尺寸的颗粒，故考虑-z方向加速度为9.81m/s2。
[0073]
特征面中除bottom、side为实体，使用环境材料，其余均为虚体，以保证颗粒自由流动。
[0074]
步骤5、根据干支流的来流情况计算流体域边界条件；其中，流体域的边界条件包括干流速度进口边界条件、支流速度进口边界条件、出口边界条件、以及河流边界条件，河流边界条件定义为粗糙度无滑移边界，粗糙度根据实际河流情况决定，采用标准壁面函数法进行计算，出口边界条件为压力出口，干流速度进口边界条件、支流进口边界条件均包括进口流速，其中：
[0075]
进口流速计算式为：
[0076][0077]
式(1)中，u
进口
为进口流速，m/s；q
进口
为进口流量，m3/s；b为进口断面河宽，m；h
进口
为进口断面平均水深，m。
[0078]
步骤6、基于内陆河水流交汇区模型与泥沙颗粒模型、微塑料颗粒模型列出流体域与颗粒相的控制方程组，控制方程组具体过程为：
[0079]
流体域使用三维瞬态方法计算，基本控制方程如下：
[0080]
质量守恒方程：
[0081][0082]
式(2)中，v表示控制体，控制体为流场单个非结构网格的体积；a表示控制面；ρ为流体密度；v为流体流速；n为控制面数；等式左边第一项表示控制体v内部质量的增量，第二项表示通过控制表面流入控制体的净通量；
[0083]
动量守恒方程：
[0084][0085]
式(3)中：f
bx
，f
by
，f
bz
分别是单位质量流体上的质量力在三个方向上的分量；p
xy
，p
yz
，p
xz
等是该流体内应力张量的分量；
[0086]
为保证赋予的重力加速度方向与实际情况保持一致，重力加速度方向依据网格本身坐标系位置确定，湍流模型使用rngk-ε：
[0087]
k方程：
[0088][0089]
ε方程：
[0090][0091]
式(4)与式(5)中，gk表示由于平均速度变化梯度引起的湍动能产生；gb是由于浮力影响引起的湍动能产生；ym为可压缩湍流脉动膨胀对总耗散率的影响；μ为湍流粘性系数。
[0092]
微塑料颗粒的控制方程如下：
[0093]
旋转运动：
[0094][0095]
式(6)中，i是转动惯性矩，单位为kg
·
m2；ω是角速度，单位为rad/s；m是作用在质
量上的总扭矩，单位为n
·
m；t是时间，单位为s；
[0096]
平移运动：
[0097][0098]
式(7)中，v是颗粒的平移速度，单位为m/s；m是颗粒的重量，单位为kg；fg是作用在颗粒上的联合重力，单位为n；fc和f
nc
是颗粒与周围颗粒或壁面之间的接触和非接触力，单位为n；
[0099]
数值积分：
[0100][0101]
式(8)中，v(t)是速度，单位为m/s；x(t)是位置，单位为m；a(t)是颗粒在给定步长的加速度，单位为m/s2；
△
t是时间步长，单位为s；同理更新颗粒的角速度和方向。
[0102]
通过读取耦合接口将流体域与颗粒相进行耦合计算，在流体域计算中，将流体域所有六面体单元的控制方程组联立，采用有限体积法建立控制方程组离散代数方程组，将流体域边界条件代入，并将流体域进行数据初始化，在颗粒相的计算中，每个颗粒都包含6个自由度，因此有两种运动方式：平动和旋转，在计算时，运用牛顿第二定律计算平动加速度和旋转加速度，在一个时间步上对他们进行数值积分，以此更新每个粒子的速度和位置。借助计算机进行数值计算求解，每隔一定时间保存迭代的计算数据，通常流体域计算数据保存时间为颗粒相计算数据保存时间的100倍。
[0103]
利用计算机进行数值求解中使用的算法为simple算法，迭代收敛条件为残差小于10-6
。
[0104]
步骤7、将保存的计算结果输入ensight软件进行后处理，得到流场、微塑料颗粒及泥沙颗粒运动参数。微塑料颗粒及泥沙颗粒运动参数包括速度、压强、湍动能、颗粒位置、颗粒留存时间。
[0105]
实施例
[0106]
如图2、图3、图4所示，水流交汇区底泥影响下微塑料输运模拟方法，具体按照以下步骤实施：
[0107]
查阅相关文献，查找内陆河流中泥沙颗粒的赋存特征，重点记录推移质泥沙颗粒的材料、尺寸、形状以及其物性(颗粒物性包含密度、泊松比、弹性模量、剪切模量、恢复系数、静摩擦系数、动摩擦系数)，如表1所示，泥沙颗粒形状图见图5，主要以球状、块状为代表，同时查找内陆河流中微塑料颗粒的赋存特征，对颗粒的材料、尺寸、形状等分类并记录微塑料颗粒的物性(同泥沙颗粒)，如表2所示，微塑料颗粒形状图见图6，主要以发泡状/球状、薄膜状、纤维状、片状为代表，本模拟方法选择三种不同密度的微塑料来进行输移规律的研究，分别是pe(聚乙烯)、ps(聚苯乙烯)和pa(聚碳酸酯)；
[0108]
表1
[0109][0110]
表2
[0111][0112][0113]
采用有限体积法建立颗粒相与流场计算使用的六面体网格，颗粒相网格图与流场的六面体网格数量均在100000左右，颗粒相及流场网格图见图7，流场说明图见图8；
[0114]
建立泥沙颗粒形状，利用图9所示的单体颗粒，以0.1mm-0.5mm为半径的单体颗粒大量拼接组成图5所示的泥沙颗粒形状，基本包括球状和块状，其中以球状居多，组合形状的尺寸分别为0.5mm、1mm和1.5mm，在底泥覆盖层注入泥沙颗粒时，使这三种粒径的泥沙颗粒各占三分之一；
[0115]
建立微塑料颗粒形状，利用图10所示的单体颗粒，以0.1mm-0.5mm为半径的单体颗粒大量拼接能够组成图6所示的微塑料颗粒非规则形状，在组合过程中单体数目需小于500，以提高计算效率，组合后的颗粒尺寸均保持在2mm，为保证不同密度的微塑料颗粒在尺寸相同的情况下进行对比；
[0116]
对泥沙颗粒定义其名字、泊松比、密度、弹性模量、其次对不同种类的微塑料命名、定义各自的泊松比、密度、剪切/弹性模量、导入由单体颗粒所组合成的泥沙及微塑料颗粒；添加微塑料颗粒、泥沙颗粒、河流边界三者间的相互作用的恢复系数、静摩擦系数、动摩擦系数；
[0117]
设置颗粒相六面体网格，命名干流入口为inlet1与inlet2，支流入口为inlet3和inlet4，流场出口outlet，流场上顶面top，底面bottom，侧面side。除bottom及side为实体，使用环境材料外，其余均为虚体。
[0118]
设置颗粒工厂，包括静态和动态颗粒两种颗粒工厂，静态颗粒为底泥覆盖区的泥沙颗粒，通过查阅文献，总结相关经验，最终设定泥沙层厚度为0.03m，分别建立干流和支流的泥沙分布区域，选择静态颗粒区域均为虚体，以保证泥沙可以正常流动，选择填充满静态颗粒分布区域，颗粒的初始速度为0；动态颗粒为微塑料颗粒，由于inlet2与inlet4是与底泥区域所重合的部分，故该部分不需要设置动态的微塑料颗粒入射，因此将inlet1与
inlet3添加动态颗粒工厂，根据不同的入射需求选择颗粒的注入模式，设定每秒从虚体面射入的颗粒数目为50，颗粒入射的方向垂直于虚体面，给定颗粒以初速度，防止颗粒在入射口淤积，颗粒入射位置于整个虚体面随机分布。因为涉及不同的泥沙颗粒与微塑料颗粒，考虑泥沙颗粒与微塑料颗粒在-z方向加速度为9.81m/s。
[0119]
通过查阅相关文献并且前期进行不同流量比的工况的数值模拟计算，最终通过公式计算出流体域边界条件：干流进口流速为0.3m/s，支流进口流速为0.2m/s，河流边界粗糙度为0.003m。
[0120]
基于内陆水流交汇区模型与颗粒模型列出流体域与颗粒相的控制方程组，控制方程组具体过程为：
[0121]
流体域使用三维瞬态方法计算，基本控制方程如下：
[0122]
质量守恒方程：
[0123][0124]
式(1)中，v表示控制体，控制体为流场单个非结构网格的体积；a表示控制面；ρ为流体密度；v为流体流速；n为控制面数；等式左边第一项表示控制体v内部质量的增量，第二项表示通过控制表面流入控制体的净通量；
[0125]
动量守恒方程：
[0126][0127]
式(2)中：f
bx
，f
by
，f
bz
分别是单位质量流体上的质量力在三个方向上的分量；p
xy
，p
yz
，p
xz
等是该流体内应力张量的分量；
[0128]
为保证赋予的重力加速度方向与实际情况保持一致，重力加速度方向依据网格本身坐标系位置确定，湍流模型使用rngk-ε：
[0129]
k方程：
[0130][0131]
ε方程：
[0132][0133]
式(3)与式(4)中，gk表示由于平均速度变化梯度引起的湍动能产生；gb是由于浮力影响引起的湍动能产生；ym为可压缩湍流脉动膨胀对总耗散率的影响；μ为湍流粘性系数。
[0134]
微塑料颗粒的控制方程如下：
[0135]
旋转运动：
[0136][0137]
式(5)中，i是转动惯性矩，单位为kg
·
m2；ω是角速度，单位为rad/s；m是作用在质量上的总扭矩，单位为n
·
m；t是时间，单位为s；
[0138]
平移运动：
[0139][0140]
式(6)中，v是颗粒的平移速度，单位为m/s；m是颗粒的重量，单位为kg；fg是作用在颗粒上的联合重力，单位为n；fc和f
nc
是颗粒与周围颗粒或壁面之间的接触和非接触力，单位为n；
[0141]
数值积分：
[0142][0143]
式(7)中，v(t)是速度，单位为m/s；x(t)是位置，单位为m；a(t)是颗粒在给定步长的加速度，单位为m/s2；
△
t是时间步长，单位为s；同理更新颗粒的角速度和方向。
[0144]
通过读取耦合接口将流体域与颗粒相进行耦合计算，在流体域计算中，将流体域所有六面体单元的控制方程组联立，采用有限体积法建立控制方程组离散代数方程组，离散格式为二阶迎风格式，将流体域边界条件代入，并将流体域进行数据初始化，在颗粒相的计算中，每个颗粒都包含6个自由度，因此有两种运动方式：平动和旋转，在计算时，运用牛顿第二定律计算平动加速度和旋转加速度，在一个时间步上对他们进行数值积分，以此更新每个粒子的速度和位置。借助计算机进行数值计算求解，求解方法为simple算法，迭代收敛条件设置为残差小于10-6
，流体域计算设置每0.02s保存一次，颗粒相计算设置每0.0002s保存一次计算数据。
[0145]
利用ensight软件进行后处理，依据计算结果分析交汇区流场速度、压强、湍动能特性，同时研究在流场与泥沙的作用下，不同形状及粒径的微塑料颗粒的速度、压强、湍动能、受力情况、运动轨迹及驻留时间。
[0146]
同等流场的计算参数下，计算工况如表3所示：
[0147]
表3
[0148][0149]
下面以上述球形微塑料颗粒的运动工况为例，说明该方法的发明创新点，具体如下：
[0150]
(1)该方法既考虑了天然河道的交汇河段对微塑料输运的影响，又考虑了河床中底泥覆盖层这一影响因素，在以往的研究及公开专利中，研究者大多仅考虑天然河道的交
汇情况，而不考虑河道中泥沙颗粒对于微塑料颗粒输运的影响；
[0151]
(2)本发明方法涉及的泥沙颗粒及微塑料颗粒类别不一、形状大小各不相同，较为贴切的模拟了天然河道中微塑料及泥沙颗粒的赋存情况；
[0152]
(3)本发明方法采用dem模型，建立了动态颗粒及静态颗粒工厂，添加了微塑料与泥沙颗粒、微塑料与微塑料颗粒、泥沙与泥沙颗粒之间的相互作用，而dpm模型则难以实现这一要求，所以本发明相比于以往的发明研究更贴近天然河道的真实情况。
[0153]
图11为该工况下各时刻三种微塑料颗粒分布示意图，对比三种颗粒的分布情况可以发现，pa颗粒几乎完全沉降至泥沙层，ps颗粒沉降较多，但有少部分颗粒有悬浮-再沉降现象，pe颗粒沉降较少，几乎全部漂浮至液体表面，在下游侧pe颗粒发生淤积。
[0154]
图12为该工况下微塑料颗粒悬浮或沉降淤积示意图，pa颗粒在随干支流水流运动过程中，在还未到达交汇区时，就已经沉降完全。ps颗粒在随干支流水流运动过程中，沉降现象发生在交汇区及交汇区下游侧，结合图11可以发现，有少部分ps颗粒在沉降后产生悬浮-再沉降现象。由图12可知pe颗粒在干流进口不远处便开始上浮，且在交汇之前pe颗粒基本全部上浮至水面，结合图11发现，受交汇区水流的影响，pe颗粒集中分布在交汇区左岸，且交汇后pe颗粒沉降极少，下游水面淤积的颗粒基本为pe颗粒。
[0155]
图13为该工况下微塑料颗粒在渠道中驻留时间示意图，通过对比驻留时间发现，仅有pe颗粒持续向下游测流动，且随着时间的增加，pe颗粒全部聚集在下游的两近岸侧，沉降-悬浮-再沉降的ps颗粒小部分在泥沙表层做推移运动，大部分被泥沙所裹挟不再运动，pa颗粒沉降后全部被泥沙裹挟，不再继续向下游侧运动。
[0156]
通过上述方式，在研究交汇区微塑料颗粒输运规律时，底泥覆盖层对颗粒的运移过程有着至关重要的影响，大部分学者在研究交汇区微塑料输运这一问题时，只考虑了汇流比及物性对其运动的影响，未考虑多形状、多种类、河道底泥覆盖层等关键影响因素，本方法通过耦合颗粒与流体运动，考虑了流场对颗粒的影响，微塑料颗粒间、微塑料颗粒与泥沙颗粒间的作用影响，微塑料颗粒自身特性的影响(形状、粒径等)，使模拟的结果更加符合实际情况，更具有说服力。

技术特征：
1.水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，具体按照以下步骤实施：步骤1、根据内陆河水流交汇区泥沙颗粒分布情况，查询泥沙颗粒的物理性质以及泥沙颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域泥沙颗粒的形状，构建泥沙颗粒模型；步骤2、根据内陆河流水流交汇区微塑料颗粒分布情况，查询微塑料颗粒的物理性质以及微塑料颗粒的材料、尺寸、形状、分布类别，选择半径为0.1mm-0.5mm的单体颗粒创建研究区域微塑料颗粒的形状，构建微塑料颗粒模型；步骤3、根据干流、直流尺寸数据建立内陆河水流交汇区模型，对内陆河水流交汇区模型进行网格划分，根据设定的最小网格尺寸生成相应的网格数量；步骤4、设置泥沙颗粒、微塑料颗粒相相关计算参数，建立静态泥沙颗粒工厂，设定动态微塑料颗粒入射参数；步骤5、根据干支流的来流情况计算流体域边界条件；步骤6、基于内陆河水流交汇区模型与泥沙颗粒模型、微塑料颗粒模型列出流体域与颗粒相的控制方程组，将流体域与颗粒相进行耦合计算，在流体域计算中，将流体域所有六面体单元的控制方程组联立，采用有限体积法建立控制方程组离散代数方程组，将流体域边界条件代入，并将全流体域压强值、流速值初始化，流体域使用三维瞬态方法计算；在颗粒相的计算中，每个颗粒都包含6个自由度，两种运动方式：平动和旋转，在计算时，运用牛顿第二定律计算平动加速度和旋转加速度，在一个时间步上对他们进行数值积分，以此更新每个粒子的速度和位置，借助计算机进行数值计算求解，每隔一定时间保存迭代的计算数据，通常流体域计算数据保存时间为颗粒相计算数据保存时间的100倍；步骤7、将保存的计算结果输入ensight软件进行后处理，得到流场、微塑料颗粒及泥沙颗粒运动参数。2.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，所述泥沙颗粒的物理性质、微塑料颗粒的物理性质包括密度、泊松比、弹性模量、剪切模量、恢复系数、静摩擦系数和动摩擦系数。3.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，所述步骤4中相关计算参数包括泥沙颗粒的物理性质、材料、分布类别，微塑料的物理性质、材料、分布类别。4.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，所述建立静态泥沙颗粒工厂包括查阅相关研究文献确定底泥覆盖层厚度、在底泥覆盖区添加静态泥沙颗粒；所述动态微塑料颗粒入射参数包括微塑料颗粒入射类型、微塑料颗粒入射位置及微塑料颗粒入射颗粒数目。5.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，所述流体域的边界条件包括干流速度进口边界条件、支流速度进口边界条件、出口边界条件、以及河流边界条件，河流边界条件定义为粗糙度无滑移边界，粗糙度根据实际河流情况决定，采用标准壁面函数法进行计算，出口边界条件为压力出口，干流速度进口边界条件、支流进口边界条件均包括进口流速，其中：进口流速计算式为：
式(1)中，u
进口
为进口流速，m/s；q
进口
为进口流量，m3/s；b为进口断面河宽，m；h
进口
为进口断面平均水深，m。6.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，所述流体域使用三维瞬态方法计算的方程如下：质量守恒方程：式(2)中，v表示控制体，控制体为流场单个非结构网格的体积；a表示控制面；ρ为流体密度；v为流体流速；n为控制面数；等式左边第一项表示控制体v内部质量的增量，第二项表示通过控制表面流入控制体的净通量；动量守恒方程：式(3)中：f
bx
，f
by
，f
bz
分别是单位质量流体上的质量力在三个方向上的分量；p
xy
，p
yz
，p
xz
等是该流体内应力张量的分量；为保证赋予的重力加速度方向与实际情况保持一致，重力加速度方向依据网格本身坐标系位置确定，湍流模型使用rngk-ε：k方程：ε方程：式(4)与式(5)中，g
k
表示由于平均速度变化梯度引起的湍动能产生；g
b
是由于浮力影响引起的湍动能产生；y
m
为可压缩湍流脉动膨胀对总耗散率的影响；μ为湍流粘性系数。7.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，步骤6中所述微塑料颗粒的控制方程如下：旋转运动：式(6)中，i是转动惯性矩，单位为kg
·
m2；ω是角速度，单位为rad/s；m是作用在质量上
的总扭矩，单位为n
·
m；t是时间，单位为s；平移运动：式(7)中，v是颗粒的平移速度，单位为m/s；m是颗粒的重量，单位为kg；f
g
是作用在颗粒上的联合重力，单位为n；f
c
和f
nc
是颗粒与周围颗粒或壁面之间的接触和非接触力，单位为n；数值积分：式(8)中，v(t)是速度，单位为m/s；x(t)是位置，单位为m；a(t)是颗粒在给定步长的加速度，单位为m/s2；
△
t是时间步长，单位为s；同理更新颗粒的角速度和方向。8.根据权利要求7所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，步骤6中流体域与颗粒相耦合计算的具体过程为，流体域计算完一个时间步长后，将计算的最新计算结果传递给泥沙颗粒、微塑料颗粒相，计算相互间作用力，将相互间作用力引入微塑料颗粒的控制方程计算微塑料颗粒速度与位置，最后将计算体积分数连同相互作用力传递回内陆河水流交汇区模型计算，从流体域开始计算到最终传递回流体域，这是一个循环过程。9.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，步骤6中采用有限体积法建立控制方程组离散代数方程组具体过程包括将每个六面体网格中心节点的控制方程在该六面体网格进行积分，积分得到六面体网格的代数方程组，该方程组中每个六面体网格边界上的参数值由二阶迎风方程得到。10.根据权利要求1所述水流交汇区底泥影响下微塑料输运模拟方法，其特征在于，步骤7中微塑料颗粒及泥沙颗粒运动参数包括速度、压强、湍动能、颗粒位置、颗粒留存时间。

技术总结
本发明公开了水流交汇区底泥影响下微塑料输运模拟方法，在研究交汇区微塑料颗粒输运规律时，底泥覆盖层对颗粒的运移过程有着至关重要的影响，大部分学者在研究交汇区微塑料输运这一问题时，只考虑了汇流比及物性对其运动的影响，未考虑多形状、多种类、河道底泥覆盖层等关键影响因素，本方法通过耦合颗粒与流体运动，考虑了流场对颗粒的影响，微塑料颗粒间、微塑料颗粒与泥沙颗粒间的作用影响，微塑料颗粒自身特性的影响(形状、粒径等)，使模拟的结果更加符合实际情况，更具有说服力。更具有说服力。更具有说服力。


技术研发人员：谌霞 李德鸿 高玮峥 蔡焕杰 李昇 袁旭年 卿登科 曹力玮
受保护的技术使用者：西北农林科技大学
技术研发日：2023.04.11
技术公布日：2023/7/20