一种基于生物絮凝的水沙动力学模拟方法与流程

本发明涉及水利工程技术领域，特别涉及生物絮凝的水沙动力学领域。

背景技术：
我国河道含沙量高、输沙量大，泥沙输移对河床演变、水利工程建设等均产生重要影响，因此泥沙研究在我国一直受到高度重视，并已形成一套相对成熟的研究体系。但传统的泥沙研究，以干净泥沙为研究对象，干净泥沙主要由岩石风化产生，其表面特性、颗粒特性和群体特性均以岩石颗粒为基础，忽略了微生物等的影响。近年来随着社会经济的快速发展，生产和生活用水的大量排放，点源和面源污染导致河道水体污染日趋严重、营养化程度增高，河流中的干净泥沙颗粒作为水环境中的重要载体，吸附了大量污染物，为微生物在其表面的吸附成膜提供了更为有利的条件，从而进一步加剧了水环境中泥沙颗粒、生物膜生长及污染物吸附的相互影响问题。这一问题传统的泥沙研究无法回答。目前已有的水沙动力学模型，只能模拟干净泥沙颗粒的絮凝及输移，以及由此所导致的河床冲淤演变及流量、水位、糙率等水沙动力特性变化。微生物是自然界中个体最小、数目最多、分布范围最广的生命体，在生长环境中大多是以附着状态而非游离状态存在，其通过代谢活动分泌称之为生物膜的胞外聚合物质(EPS)粘附于固体基底表面。在河道、湖泊水体中，有大量的微生物。观测发现，微生态系统对河流、湖泊、潮间带地区等水沙环境具有重要作用，微生物通过新陈代谢活动分泌生物膜等作用于泥沙颗粒，并形成新的泥沙生物絮凝体，又称之为“生物絮凝泥沙”或“生物絮凝体”。与传统模型研究中的普通泥沙颗粒相比，生物絮凝泥沙的理化性质均发生很大变化，生物絮凝泥沙的输移规律也发生改变。微生物吸附及生物膜的形成会强烈地改变泥沙颗粒本身的物理化学特征和运动输移规律等，并且导致泥沙颗粒的新的环境效应问题。但是对于生物絮凝影响下的水沙特性，由于其作用机理不同，影响因子由于侧重于具有生物活性的微生态系统而更显复杂化，单纯的物理化学理论已无法准确全面的描述其实际情况。因此，为了更加准确全面的描述泥沙生物絮凝体的输移状况，亟需建立一种耦合泥沙颗粒和微生态系统之间相关关系的模型系统。

技术实现要素：
本发明主要针对实际工程需要，建立了一种基于泥沙生物絮凝的水沙动力学模型系统，该系统涵盖水位流量预报、河道演进、泥沙生物絮凝变化等模拟功能，具有功能齐全，方便实用的特点。据此，本发明提供了一种基于生物絮凝的水沙动力学模拟方法，所述方法包括如下步骤：A)给定初始河道断面资料，利用上游边界条件、下游边界条件及河道糙率计算出各个断面的水力要素，进而得到整个河道的水力要素分布；B)利用所述水力要素计算出水体紊动剪切速率，然后利用所述水体紊动剪切速率和测定的水质参数计算出悬移泥沙生物絮凝体的特征数据；C)利用上游来沙过程线、所述各个断面的水力要素以及所述悬移泥沙生物絮凝体的特征计算相应断面泥沙生物絮凝体的冲刷淤积及运动输移，给出冲淤之后的河道断面情况；D)在相应的时间节点上，重复步骤A)-C)直到完成整个时段的计算，模拟河道冲淤的全过程，给出河道的演变规律。其中，步骤A)对应于水动力学模块，步骤B)对应于泥沙生物絮凝模块，步骤C)对应于泥沙动力学模块。由于传统的水沙动力学模型系统忽略生物絮凝对水沙动力学的影响，因而不能真实的反应河道中实际情况。然而，发明人通过研究发现，生物絮凝对水沙动力学具有一定的影响，其不可作为一个忽略因素来用于实际的水沙动力学评价当中。因此，将生物絮凝的因素考虑到水沙动力学系统当中，无疑使河道相关测报的数据更加科学、全面，从而更精准的指导实际现场。在一个具体实施方式中，所述水力要素包括流量、流速、水位和水深。在一个具体实施方式中，所述水质参数可以包括水体营养浓度、生物膜最大指定生长速率、营养物质的半饱和浓度、悬移泥沙生物絮凝体携带生物膜的能力系数、生物絮凝计算系数和破坏计算系数。在一个具体实施方式中，所述悬移泥沙生物絮凝体的特征可以包括悬移泥沙生物絮凝体的粒径；优选所述悬移泥沙生物絮凝体的特征数据还包括悬移泥沙生物絮凝体的密度；特别优选所述悬移泥沙生物絮凝体的特征数据还包括悬移泥沙生物絮凝体的组成和/或其结构。其中，悬移泥沙生物絮凝体的组成主要包括悬移泥沙和附着于悬移泥沙上的生物膜两部分。附着于悬移泥沙上的生物膜可以简称为悬移生物膜；悬移泥沙生物絮凝体的结构主要指悬移泥沙生物絮凝体的孔隙率等。在一个优选的实施方式中，所述悬移泥沙生物絮凝体的特征包括悬移泥沙生物絮凝体的粒径和密度。为了更精准地将生物絮凝应用于泥沙动力学模拟方法中，并且为了使具体的操作更具可行性，在本发明的一个实施例中，将所述悬移泥沙生物絮凝体分为悬移泥沙和悬移生物膜两部分，并引入分形维数理论。由此可根据絮体的分形几何特征，结合悬移泥沙生物絮凝体的体积来推导求出其粒径。在一个具体实施方式中，根据方程求出所述悬移泥沙生物絮凝体的粒径L，其中，V为所述悬移泥沙生物絮凝体的总体积，d为所述悬移泥沙生物絮凝体的分形维数，LP为所述悬移泥沙生物絮凝体中的所述悬移泥沙的泥沙颗粒的粒径，是所述悬移泥沙生物絮凝体中的所述悬移泥沙的泥沙颗粒的体积；以及任选地，ρf＝ζρB+(1-ζ)ρS，其中，ρf为悬移泥沙生物絮凝体的密度，ρB为所述悬移泥沙生物絮凝体中的所述悬移生物膜的密度，ρS为所述悬移泥沙生物絮凝体中的所述悬移泥沙的密度。在一个具体实施方式中，根据方程V＝VS+VB＝(1-ζ)V+ζV或求出所述悬移泥沙生物絮凝体的总体积V；其中，VS为所述悬移泥沙生物絮凝体中的所述悬移泥沙的总体积，VB为所述悬移泥沙生物絮凝体中的所述悬移生物膜的总体积，ζ表示所述悬移泥沙生物絮凝体中的所述悬移生物膜的体积含量，且ζ＝VB/V，(1-ζ)表示所述悬移泥沙生物絮凝体中的所述悬移泥沙的体积含量，t为时间。在此求出的V可用于方程在一个具体实施方式中，用于方程中的通过方程求得，其中，Cs是悬移泥沙生物絮凝体中的所述悬移泥沙的浓度，且CS＝(1-Ω)C，其中，Ω是总悬移泥沙生物絮凝体中的所述悬移生物膜的质量含量，C是悬移泥沙生物絮凝体的浓度；(1-ζ)是悬移泥沙生物絮凝体中的所述悬移泥沙的体积含量；G为水体紊动剪切速率；ka为悬移泥沙生物絮凝体的絮凝速率，kb为悬移泥沙生物絮凝体的破坏速率，且ka和kb计算方程分别为和其中，k'a为无量纲悬移泥沙生物絮凝体的絮凝计算系数，kb'为无量纲悬移泥沙生物絮凝体的破坏计算系数，μ为水的动力粘度系数，Fy为悬移泥沙生物絮凝体的强度。在此求出的可用于方程在一个具体实施方式中，用于方程中的通过方程求得，其中，CB是所述悬移泥沙生物絮凝体中的所述悬移生物膜的浓度，且CB＝ΩC；η为所述悬移泥沙生物絮凝体中的所述悬移生物膜的生长速率，且ηmax为所述悬移泥沙生物絮凝体中的所述悬移生物膜最大指定生长速率，N为水体的营养浓度，Km为营养物质的半饱和浓度；K为悬移泥沙生物絮凝体携带生物膜的能力，K＝βVp＝β(L3-V)，β为悬移泥沙生物絮凝体携带生物膜的能力系数，Vp为悬移泥沙生物絮凝体中的孔隙体积。在本发明中，最后得到的是河道冲淤的全过程，这其中包括河道的河床高程、河道的断面形状、河道断面的平均挟沙量、河道断面的平均含沙量以及泥沙的沉降速度等。在一个优选的实施方式中，所述河道冲淤全过程包括河道的河床高程和河道的断面形状。在本发明中，术语“泥沙颗粒”是指生长生物膜之前分散的干净的泥沙颗粒；“悬移泥沙”是指泥沙生物絮凝体中除去(或者说是不包括)生物膜部分的部分。其中，一般来说，泥沙生物絮凝体包括两部分，即泥沙和生物膜两部分。例如，悬移泥沙是指悬移在水体中的泥沙生物絮凝体除去生物膜部分后的部分。附图说明图1模拟系统结构基本框架。图2Preissmann差分格式网格。图3Zeebrugge处测得的悬浮絮凝体(SPM)的浓度、水体紊动剪切速率G和平均絮体大小L。图4BFLOC模型的率定验证结果。图5实验水槽示意图。图6干净泥沙和泥沙生物絮凝体对比。其中，(a)为干净泥沙；(b)为泥沙生物絮凝体。图7基于生物絮凝的水沙动力学模型率定结果。其中，(a)水位率定结果；(b)含沙量率定结果。图8基于生物絮凝的水沙动力学模型验证结果。其中，(a)水位验证结果；(b)含沙量验证结果。具体实施方式以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。实施例1基于生物絮凝的水沙动力学模型系统总思路基于生物絮凝的水沙动力学模型系统总思路如图1所示。其主要包括如下几部分：(1)给定初始河道断面资料，利用上游边界条件、下游边界条件以及河道糙率，通过水动力学模块计算河道沿程水力要素的大小及变化情况，得到水面线(水位)分布、水深分布、流速分布、流量分布。其中，河道断面资料包括河道断面各测点的位置和高程。上游边界条件可以包括上游流量过程线、上游水位过程线，在本实施例中仅使用上游流量过程线。下游边界条件包括下游水位过程线、下游流量过程线，在本实施例中仅使用下游流量过程线。因此，在本实施例中利用河道断面资料，上游流量过程线、下游水位过程线及河道糙率计算出水力要素分布数据。水力要素分布包括水位、水深、流速、流量等。(2)根据河道水动力要素计算水体紊动剪切速率，同时给定测得的水质参数，通过泥沙生物絮凝模块(简称生物絮凝模块)计算泥沙颗粒生物絮凝后的絮凝体的密度、絮凝体的粒径等理化特性的变化情况。其中，需要测定的水质参数包括水体营养浓度、生物膜最大指定生长速率、营养物质的半饱和浓度、悬移泥沙生物絮凝体携带生物膜的能力系数、生物絮凝和破坏计算系数。(3)根据上游来沙过程线，结合计算得到的水力要素及悬移泥沙生物絮凝体的理化特性，通过泥沙动力学模块计算各断面泥沙生物絮凝体的冲刷淤积及运动输移，给出冲淤之后的河道断面情况。(4)在相应的时间节点上，重复步骤(1)-(3)至完成整个时段的计算，模拟河道冲淤的全过程，给出河道的演变规律。实施例2模型方程及原理1.水流方程明渠流动常被考虑为一维流动，所以一维流动连续方程和运动方程是明渠水流问题的控制方程。使用一维非恒定流数学模型计算的目的是为了正确估计这些水力现象发生时的水力要素(如水位、流量、流速等)的大小及变化情况。一维非恒定流的控制方程为圣维南方程组：式中，Q为水流流量，x为断面间河段长度，t为时间，A为过水断面面积，q为支流入流单宽流量，g为重力加速度，Z为水位，C为谢才阻力系数n为糙率，R为水力半径。2.泥沙方程对于悬移质泥沙的输移，采用非饱和非均匀输沙方法计算(其中由于河道一般以对流运动为主，这里忽略扩散项的影响)。这种方法的主要特点是，每一断面的含沙量不一定刚好等于其水流挟沙力，即不一定处于饱和输沙状态。泥沙连续方程：式中S为断面平均含沙量，x为断面间河段长度，S*为断面平均挟沙力，Z0为河床高程，B为河宽，α为泥沙恢复饱和系数，ω为干净泥沙和/或泥沙生物絮凝体的沉速。3河床变形方程按河床冲淤守恒关系式中，ρ'S是冲刷或淤积泥沙的干密度。实施例3主要参数及相应的生物膜影响修正1悬移质泥沙的冲淤计算悬移质泥沙的挟沙力，引入张瑞瑾挟沙力公式进行计算：式中，S*为断面平均挟沙力；U为断面平均速度；R为水力半径，对于宽浅河道，式(5)中的R可用平均水深H代替；K(＝0.245)及m(＝0.92)分别为挟沙力系数与指数。公式(3)～(5)中的ω对单纯的泥沙颗粒或泥沙生物絮凝体，或它们的组合都适用，因此在此处未使用下标标注。本模型系统着重考虑生物絮凝对泥沙输移的影响，故将生物絮凝视为导致泥沙絮凝变化的主要因素。关于泥沙生物絮凝体沉速的计算，采用方红卫等推导的生长有生物膜的泥沙生物絮凝体的沉降公式(方红卫等.泥沙颗粒生长生物膜后沉降的实验研究Ⅱ：沉降速度计算[J].水利学报,43(4):336-391.)：式中，ω0为泥沙颗粒在无限静水中的沉速；d0为泥沙颗粒的粒径；ωf为泥沙生物絮凝体的沉速；df为泥沙生物絮凝体的粒径；ρW为水体密度；ρS为泥沙颗粒密度；ρf为泥沙生物絮凝体的密度；CD0为光滑球体阻力系数；CDf为泥沙生物絮凝体阻力系数。2生物絮凝粒径及密度计算生长有生物膜的泥沙生物絮凝体的沉速公式中，涉及到泥沙生物絮凝体粒径和密度因生物絮凝而发生的变化，引入专门的生物絮凝模型，将悬移泥沙生物絮凝体分成悬移泥沙和悬移生物膜两部分，结合分形维数的理论，模拟粘性泥沙絮凝模型，分别计算两部分的絮凝变化。具体计算方法如下：总的悬移泥沙生物絮凝体包括悬移泥沙部分和悬移生物膜部分，各部分体积存在如下关系：V＝VS+VB＝(1-ζ)V+ζV(7)其中，V为所述悬移泥沙生物絮凝体的总体积，VS为所述悬移泥沙生物絮凝体中悬移泥沙部分的总体积，VB为所述悬移泥沙生物絮凝体中悬移生物膜部分的总体积，ζ为所述悬移泥沙生物絮凝体中悬移生物膜的体积含量，且ζ＝VB/V，(1-ζ)表示述悬移泥沙生物絮凝体中的悬移泥沙的体积含量。泥沙生物絮凝体的密度ρf可表示为ρf＝ζρB+(1-ζ)ρS，其中ρB和ρS分别为泥沙生物絮凝体中的生物膜和泥沙的密度。假设泥沙生物絮凝体具有分形几何特征，则其中，d为悬移泥沙生物絮凝体的分形维数，LP为所述悬移泥沙生物絮凝体中的所述悬移泥沙的泥沙颗粒的粒径，是所述悬移泥沙生物絮凝体中的所述悬移泥沙的泥沙颗粒的体积。总的泥沙生物絮凝体的变化包括了泥沙部分的变化和生物膜部分的变化，即由于泥沙颗粒和生物膜组分具有粘结性，和项均有聚集和破裂的过程，即泥沙生物絮凝体体积由于悬移泥沙颗粒和生物膜的聚集而增加，由于它们的分离而减小。此外，泥沙生物絮凝体体积也会因微生物自身的生长和死亡而增加或减小，这点包含于项中。和定义如下：项的变化速率可通过差分链的转化通过泥沙生物絮凝体的粒径L来表达。其中，差分式(7)可得差分式(8)可得为由于泥沙颗粒的聚集和破裂而导致的泥沙生物絮凝体粒径的变化速率，采用粘性泥沙絮凝模型(Winterwerp,J.C..Asimplemodelforturbulenceinducedflocculationofcohesivesediment[J].JournalofHydraulicResearch,1998,36(3):309-326.)计算。将以上各项代入差分链中可得：式中，CS＝(1-Ω)C是悬移泥沙生物絮凝体中的泥沙部分的浓度，Ω是总的悬移泥沙生物絮凝体中悬移生物膜的质量含量，C是悬移泥沙生物絮凝体的浓度；(1-ζ)是悬移泥沙生物絮凝体中悬移泥沙的体积含量；G为水体紊动剪切速率；ka为悬移泥沙生物絮凝体的絮凝速率，kb为悬移泥沙生物絮凝体的破坏速率，ka和kb可按下式计算：式中，k'a和kb'分别为无量纲的悬移泥沙生物絮凝体的絮凝和破坏计算系数；μ为水的动力粘度系数；Fy为悬移泥沙生物絮凝体强度；ρf为悬移泥沙生物絮凝体的密度。生物膜的变化速率可写成与推导方程(10)所用的差分链相似的形式，即其中，项用于描述生物膜体积由于微生物细胞自身的生长和死亡而产生的变化，采用Logistic方程来描述(Murray,J.D..MathematicalBiologyI.AnIntroduction[M].Thirded.German:Springer,2002)。则有：式中，CB是悬移泥沙生物絮凝体中的悬移生物膜部分的浓度，CB＝ΩC；η为生物膜生长速率，ηmax为最大指定生长速率，N为水体的营养浓度，Km为营养物质的半饱和浓度；K为悬移泥沙生物絮凝体携带生物膜的能力，K＝βVp＝β(L3-V)，β为悬移泥沙生物絮凝体携带生物膜的能力系数，Vp为悬移泥沙生物絮凝体中的孔隙体积。实施例4模型方程的数值离散1圣维南方程的数值离散对于上述非线性偏微分方程组，只有在极少数理想化的情形下，才能得到解析解。在实际情况下，一般都必须采取数值方法，用有限个离散的网格点来逼近连续区域中无限个点，用这些结点上离散的近似解来逼近准确解，常用的数值方法有差分、特征线、有限元、有限体积、边界元等，但对于一维非恒定流问题，差分方法仍然应用得最普遍，其中差分格式存在许多改进的变形的方法，这里采用四点偏心Preissmann格式。四点偏心隐式差分格式的特点是围绕矩形网格中的一点M来取偏导数并进行差商逼近。其网格的距离步长Δs可以是不等距的，而时间步长一般是等间距的，点M距离已知时层为θΔt，距未知时层为(1-θ)Δt，网格形式如图2所示，M处于距离步长中间。设每一矩形网格内函数f呈直线变化，则M点的函数值f可以由四个顶点的函数值表示如下：式中，f为连续方程和运动方程中的相关变量。如果写出fn+1＝fn+Δf，上面的表达式可以写成：微分方程的形式是：这两个方程的差分格式按下列过程分别建立。利用(14)～(16)式，可得到连续方程(17)的差分形式：运动方程(18)的差分形式为：利用下面的关系将式(19)和(20)线性化，在线性化的过程中，略去增量的乘积项，例如ΔAΔQ，……等，线性化后的式(19)可写成，A1jΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1＝E1j(21)这里，对于梯形断面，即边坡系数。线性化后的式(20)可以写成，A2jΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1＝E2j(22)这里，假设如下的两个线性关系式，ΔQj＝FjΔZj+Gj(23)ΔZj＝HjΔQj+1+IjΔZj+1+Jj(24)将式(23)代入(21)，A1j(FjΔZj+Gj)+B1jΔZj+C1jΔQj+1+D1jΔZj+1＝E1j(A1jFj+B1j)ΔZj＝-C1jΔQj+1-D1jΔZj+1+(E1j-A1jGj)(25)比较式(25)和(24)，求得，将式(23)，(24)代入式(22)整理后得，比较式(29)和(23)，设α＝A2jFj+B2j，则得，利用循环计算式(26)-(31)，在追的过程中可求得系数Hj，Ij，Jj，Fj，Gj，而后在赶的过程中求出和2.泥沙方程的数值离散对于非恒定输沙的泥沙连续方程(3)，采用特征线方法求解该非恒定一阶偏微分方程可得：式中，下标0表示河段进口断面变量，Δx为断面间河段长度。实施例5模型的率定和验证1生物絮凝模型(BFLOC)的率定和验证对于计算生物絮凝泥沙粒径变化的BFLOC模型，直接采用文献中的实测数据进行率定和验证。实验数据于2003年9月8日取自比利时Zeebrugge处的北海区，采用光散射浊度传感器(OBS)量测悬浮泥沙生物絮凝体的浓度，Acustic多普勒海流剖面仪(ADCP)量测流速以计算水体紊动剪切流速，激光原位散射和透射仪(LISST)量测泥沙生物絮凝体的尺寸，测得数据如图3所示。计算时，测得的泥沙生物絮凝体浓度SPM和水体紊动剪切速率G同时作为输入参数以确定模型系统，泥沙生物絮凝体大小L用于率定和验证。计算时采用的相关参数值如表1所示：表1BFLOC模型参数值参数单位数值k'a[-]0.189kb'[s-1]×10-411.41ηmax[M]×10-66.586km[-]×10-61.159β[-]0.226Lp[μm]2Ω(0)[-]0.04N[-]×10-620ρW[kgm-3]1015ρS[kgm-3]2650ρB[kgm-3]1025Fy[N]×10-113d[-]2VS(0)[mm3]2×10-8VB(0)[mm3]1×10-8率定验证结果如图4所示，可知BFLOC模型可以较好的模拟泥沙生物絮凝体的粒径变化情况。2.基于生物絮凝的水沙动力学模型的率定和验证对于基于生物絮凝的水沙动力学模型，采用如下实验数据对其进行率定和验证。实验水槽结构如图5所示，实验水槽主要分成窄段、渐变段和宽段三部分，可分别侧重于泥沙的冲刷和淤积。其中，窄段长L1＝6m，宽b1＝33cm；过渡段长L2＝2m；宽段L3＝6m，宽b3＝50cm。实验所用泥沙用按照一定的比例配置的营养液浸泡十天，其表面及孔隙之间生长有生物膜，形成了泥沙生物絮凝体，其与干净泥沙的对比如图6所示，可知泥沙生物絮凝体与传统泥沙运动理论研究中的泥沙絮凝具有较大的不同，其运动输移规律亦会有所改变。实验分别按照31.1m3/h和21.8m3/h两种流量控制水槽水体进行循环，依次取不同位置处的8个断面，待稳定后分别测其相应的水位、含沙量等相应数据，实验数据如表2所示：表2泥沙生物絮凝体输移实验数据采用实验1的相关数据对模型相关物理参数进行率定，调整相应参数进行计算，提取不同断面的水位和含沙量进行计算值和实测值对比，率定结果如图7所示，二者拟合结果较好，由此初步确定模型参数。采用实验2的相关数据对确定参数后的模型进行验证，仍提取不同断面的水位和含沙量进行计算值和实测值对比，验证结果如图8所示，由验证结果可见，水位拟合结果较好，含沙量略有偏差，主要是由于泥沙生物絮凝粒径的变化对泥沙沉降及水流挟沙力影响较为直接，从而导致含沙量变化较为明显，即含沙量对模型参数的敏感性更高一些。