CN106096212A - 水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法 - Google Patents

Abstract

本发明提供一种水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于，包括以下步骤：步骤1. 漂浮物迁移过程模拟：基于漂浮物集合F中各个漂浮物在水体运动过程中的受力情况，将每个漂浮物都分别概化为具有一个质心的实体，求取每个漂浮物实体的六个属性的变化的过程；步骤2. 漂浮物堆积过程模拟：将具有相同离散网格编号Q _ij 的漂浮物合并为一个新的漂浮物，重新求取该新的漂浮物的质量、运动速度和质心坐标，并将新的漂浮物纳入漂浮物集合中；步骤3. 清漂过程模拟：根据漂浮物清理计划，将漂浮物进行剔除。本方法能够快速、准确地模拟漂浮物在水体中的迁移规律，从而能够提高海水体漂浮物的综合管理效率。

Description

水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法

技术领域

本发明属于水力学、环境管理和环境模拟领域，具体涉及一种水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法。

技术背景

漂浮物是依赖于浮力漂浮在水体表面的固体废弃物，包括原木、树枝、稻草、桔杆、塑料制品和人畜尸体等。漂浮物的存在，对水体水质，水面景观，供水，水产，航运等构成不利影响，不仅破坏生态环境，还会威胁饮水安全。尤其是在具有发电功能的水库区域，漂浮物顺流而下，易于聚集在河道凹岸、拦河坝前，不仅影响发电水头，减小水电枢纽的发电效益，对枢纽运行安全构成威胁。漂浮物的存在对水体正常功能的发挥造成了极大的影响，开展漂浮物的监测、清理和排除，是减少漂浮物危害的重要手段。而开展漂浮物的治理，需要对水体中漂浮物的运动规律有全面的把握。漂浮物在水体表面的迁移模拟技术，能提供漂浮物运动的实时信息，有助于漂浮物治理方案的构建。基于漂浮物模拟技术提供的实时信息，在合适的位置和时间开展清漂工作，进行漂浮物监控，有助于快速有效的清理流域漂浮物、降低排漂清漂的费用，节省人力物力。

近年来，随着漂浮物对水质、水生态、水工建筑的不利影响逐渐受到重视，针对漂浮物的研究不断发展。目前，主要在海洋区域针对漂浮物的相关研究较多，国外学者通过拉格朗日和空间建模技术，构建了漂浮物预测模型，模型能够模拟海洋中漂浮物的运移路径和轨迹，但对漂浮物的堆积以及质量缺乏量化和考虑，且缺少对河流、水库等水体漂浮物的迁移运动过程的模拟研究。我国学者在对漂浮物的研究中，主要采用采样和实验分析方法，包括研究漂浮物的丰度、粒径、空间分布及其组成成分，研究漂浮物对水工建筑的影响以及漂浮物的治理措施等。关于漂浮物的定量研究较少，如通过分析漂浮物的水流特性，提出漂浮物模拟的相似条件，应用正态模型或小变率的变态模型来模拟漂浮物的堆积形态；通过假定漂浮物与水流质点运动的轨迹一致，基于平面二维水流数值模型用水流速度模拟漂浮物运动，分析漂浮物的堆积问题。目前已有的研究在一定程度上为漂浮物数值模拟的进一步开展提供了技术借鉴，但现有研究仍然不能快速、准确的模拟漂浮物从产生到堆积或清理的过程。

漂浮物在水体中受水和风的双重动力作用，单纯用水的流速代替漂浮物的速度进行模拟，考虑不够全面，且水流数学模型采用离散化的网格形式，模拟的是网格的速度，格点的速度与漂浮物质心的速度存在一定的差别；而在漂浮物的堆积问题中，清漂措施会极大的影响漂浮物的总量和运移情况。因此，在现有研究的基础上，通过合理概化漂浮物受水流、风速的影响，综合考虑漂浮物的入汇、堆积以及清漂过程，模拟漂浮物的运移过程。为了更加有效的进行漂浮物治理，建立漂浮物模型，快速、准确地模拟漂浮物在水体中的迁移规律，对大尺度流域的产漂量进行定量研究，对漂浮物的去除，提高海水体漂浮物的综合管理有重要意义。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种能够快速、准确地模拟漂浮物在水体中的迁移规律的体表面漂浮物迁移、堆积和清漂过程的数值模拟方法。

本发明为了实现上述目的，采用了以下方案：

本发明提供一种水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于，包括以下步骤：

步骤1.漂浮物迁移过程模拟：

基于漂浮物集合F中各个漂浮物在水体运动过程中的受力情况，将每个漂浮物都分别概化为具有一个质心的实体，求取每个漂浮物实体的六个属性的变化的过程，这六个属性分别为质量Qf，正东方向流速uf，正北方向流速vf，质心在正东和正北方向的大地坐标xf和yf，以及质心所在的离散网格编号Q_ij；

步骤2.漂浮物堆积过程模拟：将具有相同离散网格编号Q_ij的漂浮物合并为一个新的漂浮物，重新求取该新的漂浮物的质量、运动速度和质心坐标，并将新的漂浮物纳入漂浮物集合中；

步骤3.清漂过程模拟：

根据水体的漂浮物清理计划，将处于清理计划所在的清漂时间段和清理区域的漂浮物进行剔除：求取当前时刻下漂浮物的六个属性，对比漂浮物所在的网格是否处于清漂区域，若是，则进一步对比当前时刻是否处于清漂时间段内，若是，则将漂浮物清理出漂浮物集合F。

进一步地，本发明提供的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法中，步骤1.漂浮物迁移过程模拟，还包括以下子步骤：

步骤1-1.漂浮物入汇条件给定：漂浮物入汇条件是漂浮物计算的前提，给定在时刻t所汇入的漂浮物的质量Qf₀，排放位置xf₀和yf₀，流速uf₀＝0和vf₀＝0，以及排放处的网格编号Q_ij0，设漂浮物实体为float，漂浮物组成的漂浮物集合为F＝{float|float为漂浮物实体}，则每排入一个漂浮物，F中的漂浮物的数量增加1个；

步骤1-2.水动力模型构建：基于二维规则网格构建技术，构建用于漂浮物模拟的二维规则网格体系G[ij]及其水动力模型，该水动力模型基于二维圣维南方程，如下：

连续方程：

$\frac{\∂h}{\∂t}+\frac{\∂\left(uh\right)}{\∂x}+\frac{\∂\left(vh\right)}{\∂y}=0,$

动量方程：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=fv-g\frac{\∂z}{\∂x}-\frac{gu\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_x{\▿}^{2}u+\frac{{\mathrm{\τ}}_x}{\mathrm{\ρ}h},$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-fu-g\frac{\∂z}{\∂y}-\frac{gv\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_y{\▿}^{2}v+\frac{{\mathrm{\τ}}_y}{\mathrm{\ρ}h},$

式中：h为水深，u为正东方向的流速，v为正北方向的流速，x为正东方向的投影坐标，y为正北方向的投影坐标，t为时间，g为重力加速度，ρ为水体密度，c为谢才系数，f为柯氏力常数， 为纬度，Ω为地转角速度，约为2π/(24×3600)弧度每秒，ξ_x、ξ_y为分别为x、y方向上的涡动粘滞系数；τ_x、τ_y为分别为x、y方向上的风切应力，其表达形式为：

${\mathrm{\τ}}_x=C_a{\mathrm{\ρ}}_a{W}_x{(W_x^2+W_y^2)}^{1/2},$

${\mathrm{\τ}}_y=C_a{\mathrm{\ρ}}_a{W}_y{(W_x^2+W_y^2)}^{1/2},$

式中：C_a为风阻力系数；ρ_a为空气密度；W_x，W_y为分别为x、y方向上的风速；

通过二维水动力模型，模拟任意时刻每个网格上的正东方向和正北方向的表面流速，u_ij和v_ij；

步骤1-3.气象插值：构建二维规则网格内及其周边的气象站集合W，W中包括有每个气象站的位置(xw,yw)，该气象站t时刻监测到的正东方向的风速wu和正北方向的风速wv，通过反距离权重法，将集合W中所有气象站的风速插值到每个二维规则网格，得到t时刻每个网格上的正东方向和正北方向的风速为wu_ij，wv_ij；反距离权重插值的方程如下：

${\mathrm{wu}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wu\left(s_i\right),$

${\mathrm{wv}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wv\left(s_i\right),$

${\mathrm{\λ}}_i=\frac{d_{i0}}{\underset{i=1}{\overset{N}{\Σ}}{d}_{i0}},$

式中，λ_i为第s_i个气象站所占权重，wu(s_i)，wv(s_i)为第s_i个气象站实测的正东和正北方向的风速值，d_i0为网格ij与其第s_i个气象站之间的距离，N为气象站个数；

步骤1-4.漂浮物运移速度计算：漂浮物在水体中受水和风的作用力而发生变化，将漂浮物所受的力概化成流速和风速对漂浮物运动速度的影响，漂浮物质心所在位置的运动速度uf和vf通过下式求取：

uf＝c_w×wub+c_a×ub，

vf＝c_w×wvb+c_a×vb，

式中，c_w为漂浮物运动的水阻系数，c_a为漂浮物运动的风阻系数，wub，wvb，ub，vb分别为漂浮物质心所在位置正东和正北方向的风速和流速；在所述步骤1-2和所述步骤1-3中，求取到的速度为网格的流速(u_ij和v_ij)与风速(wu_ij和wv_ij)；wub，wvb，ub，vb通过漂浮物所在网格1及其相邻网格2～9的流速和风速，利用反距离权重法插值获得：

$\begin{array}{c}wub=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{\mathrm{wu}}_{ij}\\ wvb=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{\mathrm{wv}}_{ij}\\ ub=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{u}_{ij}\\ vb=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{v}_{ij}\\ {\mathrm{\λ}}_{ij}=\frac{d_{ij}}{\underset{ij=1}{\overset{9}{\Σ}}{d}_{ij}}\\ ij=\left[\begin{array}{ccc}8& 5& 9\\ 2& 1& 3\\ 6& 4& 7\end{array}\right]\end{array},$

其中，i为与漂浮物相关的网格编号，1为漂浮物质心所在的网格，2,3,4,5分别为1的西、东、南、北方向的网格，6,7为4号网格的西、东方向的网格、8,9为5号网格的西、东方向的网格，d_ij为漂浮物的质心(xf，yf)到网格质心ij(x，y)的距离；

步骤1-5.漂浮物位置计算：漂浮物在时刻t所处的位置为xf，yf，则其位置满足：

xf_t＝xf_(t-1)+uf×Δt

yf_t＝yf_(t-1)+vf×Δt，

其中，xf_t，yf_t为漂浮物当前时刻t的位置，xf_t-1，yf_t-1为上一时刻的位置，Δt为当前时刻与上一时刻的时间差，等于一个计算时间步长；若t-1时刻为漂浮物排入的时刻，xf_t-1＝xf₀，yf_t-1＝yf₀；

步骤1-6，漂浮物所在网格计算：

在漂浮物模拟计算中，一个计算时间步长内，漂浮物运移的速度可能超过一个网格的范围，采用模式分离法，通过动步长模式，计算漂浮物所在的网格；设定当前漂浮物所在的网格为ij，网格的x方向的长度为dx，y方向的长度为dy，模型模拟步长为dt，则用于计算漂浮物运移的步长为dtf＝max(dt/(dx/uf)，dt/(dy/vf))；

在一个漂浮物模拟步长后，漂浮物的移动距离为dxf和dyf，其当前的位置为xf_t＝xf_t-1+dxf，yf_t＝yf_t-1+dyf，设运动前漂浮物所在的网格号为ij(i,j)，在ij网格的东西南北方向外扩至少1层网格进行搜索；

若漂浮物质心所在的位置Q与网格ij_N四个顶点(A,B,C,D)所组成的向量及网格四条边的向量，满足关系式：

AB×AQ＞0

BC×BQ＞0

CD×CQ＞0

DA×DQ＞0，

则Q在网格ij_N(ABCD)的范围内，且有漂浮物所在的网格编号Q_ij＝ij_N。

进一步地，本发明提供的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法中，步骤2.漂浮物堆积过程模拟，包括以下子步骤：

步骤2-1.判断漂浮物集合中是否存在具有相同编号的漂浮物：

在完成一个步长dt的计算后，检索漂浮物集合F中所有漂浮物所在的网格编号，在集合F中，从第1个漂浮物向第N个漂浮物依次扫描，这里N为F中，漂浮物的总数，以当前编号为n∈[1，N]的漂浮物的所在网格编号Qij(n)为判断依据，从编号为第N的漂浮物向第n+1个漂浮物，编号nk依次减1进行扫描，判断Qij(n)与Qij(nk)的大小；

若不存在网格编号相同的漂浮物，则继续下一个步长的计算；

若存在具有相同编号的漂浮物，则转入步骤2-2；

步骤2-2.合并相同网格内的漂浮物：

在漂浮物体系F中，当Qij(n)＝Qij(nk)时，合并漂浮物n与nk为新的漂浮物，并将该新的漂浮物的编号设为n，它的坐标与速度满足：

$\begin{array}{c}xf\left(n\right)=\frac{xf\left(n\right)Qf\left(n\right)+xf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ yf\left(n\right)=\frac{yf\left(n\right)Qf\left(n\right)+yf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ uf\left(n\right)=\frac{uf\left(n\right)Qf\left(n\right)+uf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ vf\left(n\right)=\frac{vf\left(n\right)Qf\left(n\right)+vf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\end{array},$

新的漂浮物的质量为两个合并的漂浮物质量之和：Qf(n)＝Qf(n)+Qf(nk)；

步骤2-3.对漂浮物集合F进行重新排列：

在将漂浮物n和漂浮物nk合并后，重新排列漂浮物体系，漂浮物1～nk-1保持序号不变，漂浮物nk+1～N序号依次减1，得到合并后总数为N-1的漂浮物集合F；

重复执行步骤2-1～步骤2-3，直到具有相同ij号的漂浮物都被合并并且新的漂浮物都被重新排列完成为止。

进一步地，本发明提供的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法中，步骤3.清漂过程模拟，还包括以下子步骤：

步骤3-1.定义清漂方案：按照水体中清漂的规划方案，定义清漂的时间段为TC[t1,t2]，清漂位置集为FC＝{Fij|Fij为清漂的网格编号}；

步骤3-2.清理满足条件的漂浮物：

在模拟过程中，若满足条件：当前所模拟的时间t满足t1≤t≤t2，且当前漂浮物所在的网格Qij(n)∈FC，则将编号为n的漂浮物清理出漂浮物集合F；

并将漂浮物1～n保持不变，漂浮物n+1～N序号依次减1，得到清出漂浮物编号为n后的新的漂浮物集合F；

重复执行步骤3-1至步骤3-2，直到所有满足条件的漂浮物都被清理完为止。

进一步地，本发明提供的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法还可以具有这样的特征：在步骤1-6中，是在网格的东西南北方向外扩3层网格进行搜索。

发明的作用与效果

本发明与现有技术相比，具有以下优点和效果：

(1)将漂浮物概化为具有6个属性的实体，漂浮物移动、堆积和清理的过程，即是对漂浮物6个属性的定量求解的过程。本技术相比于现有漂浮物的相关技术，具有明显的创新性，技术可行，计算速度快。

(2)基于水力学领域成熟的平面二维水动力模拟技术，以离散网格的形式获得水体表面的流速。分析漂浮物在水体所受的水流作用力和风作用力，通过流速和风速的概化模式，基于漂浮物的质心与网格体系的距离关系，利用反距离权重法，获得漂浮物的运移速度。通过该方法的处理，避免了对漂浮物的形态，在水体中所受复杂的力学结构的解析，具有更好的可行性。

(3)利用模式分裂法，通过构建模型模拟步长和漂浮物运移计算步长的方式，进行漂浮物移动的模拟，确保漂浮物在移动的过程中，每次计算步长内其运移的距离不会超过一个网格的边长，从而确保对漂浮物所在网格的搜索更准确和可靠。

(4)基于构建的漂浮物迁移模拟技术，能对任何水体中浮在水面的漂浮物进行模拟和预测，确定漂浮物在任意时刻所处的位置以及水体中任意位置所含有的漂浮物质量，分析漂浮物的运移轨迹，为漂浮物的监控、清理、处置提供数据支撑。

综上所述，本发明的体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，能够快速、准确地模拟漂浮物在水体中的迁移规律，对大尺度流域的产漂量进行定量研究，对漂浮物的去除，从而能够切实提高海水体漂浮物的综合管理效率。

附图说明

图1为本发明实施例的模拟范围网格体系的示意图；

图2为本发明实施例的漂浮物入汇位置和过程的示意图；

图3为本发明实施例中通过气象站插值获得网格上的风速的示意图；

图4为本发明实施例中漂浮物质心与周边9个网格的反距离权重插值的示意图；

图5为本发明实施例中漂浮物运移后搜索其所在网格ij号的示意图；

图6为本发明实施例中合并具有相同网格号的漂浮物的流程图；

图7为本发明实施例中模拟2天后漂浮物运移状态的效果图。

具体实施方式

以下结合附图对本发明涉及的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法的具体实施方案进行详细地说明。

实施例

本实施例的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法包括以下步骤：

步骤1.漂浮物迁移过程模拟

1-1.漂浮物水动力模型构建：利用规则网格绘制软件，如delft3d等软件绘制要模拟区域的二维规则网格体系G[ij]，如图1所示。基于二维圣维南方程，构建水动力模型，水动力模型基于二维圣维南方程，如下：

连续方程：

$\frac{\∂h}{\∂t}+\frac{\∂\left(uh\right)}{\∂x}+\frac{\∂\left(vh\right)}{\∂y}=0---\left(1\right)$

动量方程：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=fv-g\frac{\∂z}{\∂x}-\frac{gu\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_x{\▿}^{2}u+\frac{{\mathrm{\τ}}_x}{\mathrm{\ρ}h}---\left(2\right)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-fu-g\frac{\∂z}{\∂y}-\frac{gv\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_y{\▿}^{2}v+\frac{{\mathrm{\τ}}_y}{\mathrm{\ρ}h}---\left(3\right)$

式中：式中：h为水深；u为正东方向的流速；v为正北方向的流速，x为正东方向的投影坐标，y为正北方向的投影坐标；g为重力加速度；ρ为水体密度；c为谢才系数；f为柯氏力常数， 为纬度，Ω为地转角速度，约为2π/(24×3600)弧度每秒；ξ_x、ξ_y为分别为x、y方向上的涡动粘滞系数；τ_x、τ_y为分别为x、y方向上的风切应力，其表达形式为：

${\mathrm{\τ}}_x=C_a{\mathrm{\ρ}}_a{W}_x{(W_x^2+W_y^2)}^{1/2}---\left(4\right)$

${\mathrm{\τ}}_y=C_a{\mathrm{\ρ}}_a{W}_y{(W_x^2+W_y^2)}^{1/2}---\left(5\right)$

式中：C_a为风阻力系数；ρ_a为空气密度；W_x，W_y为分别为x、y方向上的风速。通过二维水动力模型，模拟任意时刻每个网格上的正东方向和正北方向的表面流速，u_ij和v_ij。在本实施例中，设定模拟时间为30天，模拟步长为360秒，模拟的第1天的流速如附图1所示。

1-2.给定入汇条件：在漂浮物的计算中，漂浮物入汇条件是漂浮物计算的前提，在本实施例中给定从时刻0到第30天的漂浮物的入汇量。如图2中，为30天内在入汇点S1，S2处的漂浮物的入汇位置及其入汇过程线。

1-3.气象插值:构建二维规则网格内及其周边的气象站集合W，在本实施例中W包括有2个气象站，分别W1和W2。其中，W1的位置为(x_w1，y_w1)，W2的位置为(x_w2，y_w2)。每个气象站t时刻监测到的正东方向的风速(wu)和正北方向的风速(wv)。通过反距离权重法，将集合W中所有气象站的风速插值到每个二维规则网格，得到t时刻每个网格上的正东方向和正北方向的风速为wu_ij，wv_ij。反距离权重插值的方程如下：

$\begin{array}{c}{\mathrm{wu}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wu\left(s_i\right)\\ {\mathrm{wv}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wv\left(s_i\right)\end{array}---\left(6\right)$

${\mathrm{\λ}}_i=\frac{d_{i0}}{\underset{i=1}{\overset{N}{\Σ}}{d}_{i0}}=\frac{\sqrt{{(x_{ij}-x_{w_i})}^2+{(y_{ij}-y_{w_i})}^2}}{\underset{i=1}{\overset{2}{\Σ}}\sqrt{{(x_{ij}-x_{w_i})}^2+{(y_{ij}-y_{w_i})}^2}}---\left(7\right)$

式中，λ_i为第s_i个气象站所占权重，x_ij，y_ij为网格ij的x，y值；wu(s_i)，wv(s_i)为第s_i个气象站实测的正东和正北方向的风速值；d_i0为网格ij与其第s_i个气象站之间的距离。插值获得的第1天的网格上的风速如图3所示。

1-4.漂浮物运移速度计算：漂浮物在水体中受水和风的作用力而发生变化。将漂浮物所受的力概化成流速和风速对漂浮物运动速度的影响。漂浮物质心所在位置的运动速度uf，vf通过下式求取：

$\begin{array}{c}uf=c_w\×wub+c_a\×ub\\ vf=c_w\×wvb+c_a\×vb\end{array}---\left(8\right)$

其中，c_w为漂浮物运动的水阻系数，c_a为漂浮物运动的风阻系数；在本例中，水阻系数为0.8，风阻系数为0.2。wub，wvb，ub，vb分别为漂浮物质心所在位置正东和正北方向的风速和流速。利用漂浮物所在的位置与周边9个网格ij1～ij9的距离，通过反距离权重法，获得漂浮物质心到9个网格的权重(如图4所示)，将9个网格的流速(u_ij和v_ij)与风速(wu_ij和wv_ij)插值到漂浮物所在的位置，通过公式(8)求取到漂浮物所在位置的速度uf，wf。

1-5.漂浮物位置计算：漂浮物在时刻t所处的位置为xf，yf，则其位置满足：

$\begin{array}{c}{\mathrm{xf}}_t={\mathrm{xf}}_{(t-1)}+uf\×\mathrm{\Δ}t\\ {\mathrm{yf}}_t={\mathrm{yf}}_{(t-1)}+vf\×\mathrm{\Δ}t\end{array}---\left(9\right)$

其中，xf_t，yf_t为漂浮物当前时刻t的位置，xf_t-1，yf_t-1为上一时刻的位置，Δt为当前时刻与上一时刻的时间差，等于一个计算时间步长。若t-1时刻为漂浮物排入的时刻，则有xf_t-1＝xf₀，yf_t-1＝yf₀。

1-6.漂浮物所在网格计算：在漂浮物模拟计算中，一个时间步长60秒内，漂浮物运移的速度可能超过一个网格的范围，采用模式分离法，通过动步长模式，计算漂浮物所在的网格。设定当前漂浮物所在的网格为ij，网格的x方向的长度为dx，y方向的长度为dy，模型模拟步长为dt，则用于计算漂浮物运移的步长为dtf＝max(dt/(dx/uf)，dt/(dy/vf))。在一个漂浮物模拟步长后，漂浮物的移动距离为dxf和dyf，其当前的位置为xf_t＝xf_t-1+dxf，yf_t＝yf_t-1+dyf。设运动前漂浮物所在的网格号为ij(i,j)，在ij网格的东西南北方向外扩3层网格，即[i-3，i+3]与[j-3，j+3]的范围内搜索。若漂浮物质心所在的位置Q与网格ij_N四个顶点(A，B，C，D)所组成的向量及网格四条边的向量，满足关系式：

$\begin{array}{c}AB\×AQ>0\\ BC\×BQ>0\\ CD\×CQ>0\\ DA\×DQ>0\end{array}---\left(10\right)$

则Q在网格ij_N(ABCD)的范围内，且有漂浮物所在的网格编号Q_ij＝ij_N。漂浮物由ij处运移到Q处过程中，网格搜索的范围和确定漂浮物所在网格的ij如图5所示。

步骤2.漂浮物堆积过程模拟

2-1.漂浮物集合中，判断所有漂浮物所在网格：在完成一个步长dt的计算后，检索漂浮物集合F中，所有漂浮物所在的网格编号。在集合F中，从第1个漂浮物向第N个漂浮物依次扫描(N为F中，漂浮物的总数)，以当前编号为n∈[1，N]的漂浮物的所在网格编号Qij(n)为判断依据，从编号为第N的漂浮物向第n+1个漂浮物，编号nk依次减1进行扫描，判断Qij(n)与Qij(nk)的大小。若不存在若不存在网格编号相同的漂浮物，则继续下一个步长的计算。若存在具有相同编号的漂浮物，则转入步骤2-2。

2-2.合并相同网格内的漂浮物：当漂浮物体系F中，Qij(n)＝Qij(nk)时，合并漂浮物n与nk为新的漂浮物，其编号为n，其坐标与速度满足：

$\begin{array}{c}xf\left(n\right)=\frac{xf\left(n\right)Qf\left(n\right)+xf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ yf\left(n\right)=\frac{yf\left(n\right)Qf\left(n\right)+yf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ uf\left(n\right)=\frac{uf\left(n\right)Qf\left(n\right)+uf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ vf\left(n\right)=\frac{vf\left(n\right)Qf\left(n\right)+vf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\end{array}---\left(11\right)$

漂浮物质量为两个合并的漂浮物质量之和，Qf(n)＝Qf(n)+Qf(nk)

2-3.漂浮物体系F重新排列：在将漂浮物n和漂浮物nk合并后，重新排列漂浮物体系，漂浮物1～nk-1保持序号不变，漂浮物nk+1～N序号依次减1，得到合并后总数为N-1的漂浮物集合F。重复执行步骤2-1～步骤2-3，直到具有所有相同ij号的漂浮物合并完成为止。漂浮物堆积的计算流程如图6所示。

步骤3.漂浮物清理过程模拟

3-1.定义清漂方案：按照水体中清漂的规划方案，定义清漂的时间段为TC[t1，t2]和清漂位置集FC＝{Fij|Fij为清漂的网格编号}。

3-2.清理满足条件的漂浮物：在模拟过程中，若当前所模拟的时间t满足t1≤t≤t2，且当前漂浮物所在的网格Qij(n)∈FC，则将编号为n的漂浮物清理出F。将漂浮物1～n保持不变，漂浮物n+1～N序号依次减1，则得到清出漂浮物编号为n后的漂浮物体系F。重复执行，直到所有满足条件的漂浮物被清理完为止。经过计算，模拟的2天后计算范围内漂浮物的状态如图7所示，从图中可以看出，漂浮物的移动轨迹与水流流场的轨迹相似。

实施例的作用与效果：

(1)将漂浮物概化为具有6个属性的实体，漂浮物移动、堆积和清理的过程，即是对漂浮物6个属性的定量求解的过程，技术可行，计算速度快。

(2)基于水力学领域成熟的平面二维水动力模拟技术，以离散网格的形式获得水体表面的流速。分析漂浮物在水体所受的水流作用力和风作用力，通过流速和风速的概化模式，基于漂浮物的质心与网格体系的距离关系，利用反距离权重法，获得漂浮物的运移速度。通过该方法的处理，避免了对漂浮物的形态，在水体中所受复杂的力学结构的解析，具有更好的可行性。

(3)利用模式分裂法，通过构建模型模拟步长和漂浮物运移计算步长的方式，进行漂浮物移动的模拟，确保漂浮物在移动的过程中，每次计算步长内其运移的距离不会超过一个网格的边长，从而确保对漂浮物所在网格的搜索更准确和可靠。

(4)基于构建的漂浮物迁移模拟技术，能对任何水体中浮在水面的漂浮物进行模拟和预测，确定漂浮物在任意时刻所处的位置以及水体中任意位置所含有的漂浮物质量，分析漂浮物的运移轨迹，为漂浮物的监控、清理、处置提供数据支撑。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

Claims (5)

1.一种水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于，包括以下步骤：

步骤1.漂浮物迁移过程模拟：

基于漂浮物集合F中各个漂浮物在水体运动过程中的受力情况，将每个漂浮物都分别概化为具有一个质心的实体，求取每个漂浮物实体的六个属性的变化的过程，这六个属性分别为质量Qf，正东方向流速uf，正北方向流速vf，质心在正东和正北方向的大地坐标xf和yf，以及质心所在的离散网格编号Q_ij；

步骤2.漂浮物堆积过程模拟：

将具有相同离散网格编号Q_ij的漂浮物合并为一个新的漂浮物，重新求取该新的漂浮物的质量、运动速度和质心坐标，并将新的漂浮物纳入漂浮物集合中；

步骤3.清漂过程模拟：

根据水体的漂浮物清理计划，将处于清理计划所在的清漂时间段和清理区域的漂浮物进行剔除：求取当前时刻下漂浮物的六个属性，对比漂浮物所在的网格是否处于清漂区域，若是，则进一步对比当前时刻是否处于清漂时间段内，若是，则将漂浮物清理出漂浮物集合F。

2.根据权利要求1所述的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于：

步骤1.漂浮物迁移过程模拟，包括以下子步骤：

步骤1-1.漂浮物入汇条件给定：给定在时刻t所汇入的漂浮物的质量Qf₀，排放位置xf₀和yf₀，流速uf₀＝0和vf₀＝0，以及排放处的网格编号Q_ij0，设漂浮物实体为float，漂浮物组成的漂浮物集合为F＝{float|float为漂浮物实体}，则每排入一个漂浮物，F中的漂浮物的数量增加1个；

步骤1-2.水动力模型构建：基于二维规则网格构建技术，构建用于漂浮物模拟的二维规则网格体系G[ij]及其水动力模型，该水动力模型基于二维圣维南方程，如下：

连续方程：

$\frac{\∂h}{\∂t}+\frac{\∂\left(uh\right)}{\∂x}+\frac{\∂\left(vh\right)}{\∂y}=0,$

动量方程：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=fv-g\frac{\∂z}{\∂x}-\frac{gu\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_x{\▿}^{2}u+\frac{{\mathrm{\τ}}_x}{\mathrm{\ρ}h},$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-fu-g\frac{\∂z}{\∂y}-\frac{gv\sqrt{u^2+v^2}}{\left(c^{2}h\right)}+{\mathrm{\ξ}}_y{\▿}^{2}v+\frac{{\mathrm{\τ}}_y}{\mathrm{\ρ}h},$

式中：h为水深，u为正东方向的流速，v为正北方向的流速，x为正东方向的投影坐标，y为正北方向的投影坐标，g为重力加速度，t为时间，ρ为水体密度，c为谢才系数，f为柯氏力常数， 为纬度，Ω为地转角速度，约为2π/(24×3600)弧度每秒，ξ_x、ξ_y为分别为x、y方向上的涡动粘滞系数；τ_x、τ_y为分别为x、y方向上的风切应力，其表达形式为：

${\mathrm{\τ}}_x=C_a{\mathrm{\ρ}}_a{W}_x{(W_x^2+W_y^2)}^{1/2},$

${\mathrm{\τ}}_y=C_a{\mathrm{\ρ}}_a{W}_y{(W_x^2+W_y^2)}^{1/2},$

式中：C_a为风阻力系数；ρ_a为空气密度；W_x，W_y为分别为x、y方向上的风速；

通过二维水动力模型，模拟任意时刻每个网格上的正东方向和正北方向的表面流速，u_ij和v_ij；

步骤1-3.气象插值：构建二维规则网格内及其周边的气象站集合W，W中包括有每个气象站的位置(xw,yw)，该气象站t时刻监测到的正东方向的风速wu和正北方向的风速wv，通过反距离权重法，将集合W中所有气象站的风速插值到每个二维规则网格，得到t时刻每个网格上的正东方向和正北方向的风速为wu_ij，wv_ij，反距离权重插值的方程如下：

${\mathrm{wu}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wu\left(s_i\right),$

${\mathrm{wv}}_{ij}=\underset{i=0}{\overset{N}{\Σ}}{\mathrm{\λ}}_{i}wv\left(s_i\right),$

${\mathrm{\λ}}_i=\frac{d_{i0}}{\underset{i=1}{\overset{N}{\Σ}}{d}_{i0}},$

式中，λ_i为第s_i个气象站所占权重，wu(s_i)，wv(s_i)为第s_i个气象站实测的正东和正北方向的风速值，d_i0为网格ij与其第s_i个气象站之间的距离，N为气象站个数；

步骤1-4.漂浮物运移速度计算：漂浮物在水体中受水和风的作用力而发生变化，将漂浮物所受的力概化成流速和风速对漂浮物运动速度的影响，漂浮物质心所在位置的运动速度uf和vf通过下式求取：

uf＝c_w×wub+c_a×ub，

vf＝c_w×wvb+c_a×vb，

式中，c_w为漂浮物运动的水阻系数，c_a为漂浮物运动的风阻系数，wub，wvb，ub，vb分别为漂浮物质心所在位置正东和正北方向的风速和流速；在所述步骤1-2和所述步骤1-3中，求取到的速度为网格的流速(u_ij和v_ij)与风速(wu_ij和wv_ij)；wub，wvb，ub，vb通过漂浮物所在网格1及其相邻网格2～9的流速和风速，利用反距离权重法插值获得：

$wub=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{\mathrm{wu}}_{ij}$

$wvb=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{\mathrm{wv}}_{ij}$

$ub=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{u}_{ij}$

$\begin{array}{c}vb=\underset{ij=1}{\overset{9}{\Σ}}{\mathrm{\λ}}_{ij}{v}_{ij}\\ {\mathrm{\λ}}_{ij}=\frac{d_{ij}}{\underset{ij=1}{\overset{9}{\Σ}}{d}_{ij}}\end{array},$

$ij=\left[\begin{array}{ccc}8& 5& 9\\ 2& 1& 3\\ 6& 4& 7\end{array}\right]$

其中，i为与漂浮物相关的网格编号，1为漂浮物质心所在的网格，2,3,4,5分别为1的西、东、南、北方向的网格，6,7为4号网格的西、东方向的网格、8,9为5号网格的西、东方向的网格，d_ij为漂浮物的质心(xf，yf)到网格质心ij(x，y)的距离；

步骤1-5.漂浮物位置计算：漂浮物在时刻t所处的位置为xf，yf，则其位置满足：

xf_t＝xf_(t-1)+uf×Δt

yf_t＝yf_(t-1)+vf×Δt，

其中，xf_t，yf_t为漂浮物当前时刻t的位置，xf_t-1，yf_t-1为上一时刻的位置，Δt为当前时刻与上一时刻的时间差，等于一个计算时间步长，若t-1时刻为漂浮物排入的时刻，xf_t-1＝xf₀，yf_t-1＝yf₀；

步骤1-6，漂浮物所在网格计算：

在漂浮物模拟计算中，一个计算时间步长内，漂浮物运移的速度可能超过一个网格的范围，采用模式分离法，通过动步长模式，计算漂浮物所在的网格；设定当前漂浮物所在的网格为ij，网格的x方向的长度为dx，y方向的长度为dy，模型模拟步长为dt，则用于计算漂浮物运移的步长为dtf＝max(dt/(dx/uf)，dt/(dy/vf))；

在一个漂浮物模拟步长后，漂浮物的移动距离为dxf和dyf，其当前的位置为xf_t＝xf_t-1+dxf，yf_t＝yf_t-1+dyf，设运动前漂浮物所在的网格号为ij(i,j)，在ij网格的东西南北方向外扩至少1层网格进行搜索；

若漂浮物质心所在的位置Q与网格ij_N四个顶点(A,B,C,D)所组成的向量及网格四条边的向量，满足关系式：

AB×AQ＞0

BC×BQ＞0

CD×CQ＞0

DA×DQ＞0，

则Q在网格ij_N(ABCD)的范围内，且有漂浮物所在的网格编号Q_ij＝ij_N。

3.根据权利要求1所述的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于：

步骤2.漂浮物堆积过程模拟，包括以下子步骤：

步骤2-1.判断漂浮物集合中是否存在具有相同编号的漂浮物：

在完成一个步长dt的计算后，检索漂浮物集合F中所有漂浮物所在的网格编号，在集合F中，从第1个漂浮物向第N个漂浮物依次扫描，这里N为F中，漂浮物的总数，以当前编号为n∈[1，N]的漂浮物的所在网格编号Qij(n)为判断依据，从编号为第N的漂浮物向第n+1个漂浮物，编号nk依次减1进行扫描，判断Qij(n)与Qij(nk)的大小；

若不存在网格编号相同的漂浮物，则继续下一个步长的计算；

若存在具有相同编号的漂浮物，则转入步骤2-2；

步骤2-2.合并相同网格内的漂浮物：

在漂浮物体系F中，当Qij(n)＝Qij(nk)时，合并漂浮物n与nk为新的漂浮物，并将该新的漂浮物的编号设为n，它的坐标与速度满足：

$\begin{array}{c}xf\left(n\right)=\frac{xf\left(n\right)Qf\left(n\right)+xf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ yf\left(n\right)=\frac{yf\left(n\right)Qf\left(n\right)+yf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ uf\left(n\right)=\frac{uf\left(n\right)Qf\left(n\right)+uf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\\ vf\left(n\right)=\frac{vf\left(n\right)Qf\left(n\right)+vf\left(nk\right)Qf\left(nk\right)}{Qf\left(n\right)+Qf\left(nk\right)}\end{array},$

新的漂浮物的质量为两个合并的漂浮物质量之和：Qf(n)＝Qf(n)+Qf(nk)；

步骤2-3.对漂浮物集合F进行重新排列：

在将漂浮物n和漂浮物nk合并后，重新排列漂浮物体系，漂浮物1～nk-1保持序号不变，漂浮物nk+1～N序号依次减1，得到合并后总数为N-1的漂浮物集合F；

重复执行步骤2-1～步骤2-3，直到具有相同ij号的漂浮物都被合并并且新的漂浮物都被重新排列完成为止。

4.根据权利要求1所述的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于：

步骤3.清漂过程模拟，包括以下子步骤：

步骤3-1.定义清漂方案：按照水体中清漂的规划方案，定义清漂的时间段为TC[t1,t2]，清漂位置集为FC＝{Fij|Fij为清漂的网格编号}；

步骤3-2.清理满足条件的漂浮物：

在模拟过程中，若满足条件：当前所模拟的时间t满足t1≤t≤t2，且当前漂浮物所在的网格Qij(n)∈FC，则将编号为n的漂浮物清理出漂浮物集合F；

并将漂浮物1～n保持不变，漂浮物n+1～N序号依次减1，得到清出漂浮物编号为n后的新的漂浮物集合F；

重复执行步骤3-1至步骤3-2，直到所有满足条件的漂浮物都被清理完为止。

5.根据权利要求2所述的水体表面漂浮物迁移、堆积和清漂过程的数值模拟方法，其特征在于：

其中，在步骤1-6中，是在网格的东西南北方向外扩3层网格进行搜索。