CN107798176B - 一种高低浓度自适应的泥沙运动数值模拟方法 - Google Patents

Abstract

本发明公开了一种高低浓度自适应的泥沙运动数值模拟方法：(1)建立高低浓度自适应的泥沙运动数学模型；(2)计算边界条件、悬沙沉速、床面切应力、稳定性函数及适应于高低浓度泥沙条件的修正；(3)求解高低浓度自适应的泥沙运动数学模型；(4)采用实测资料对泥沙运动数学模型进行水流验证，含沙量验证，航道淤积分布验证；(5)根据泥沙运动数学模型结果，分析各因素对河口航道泥沙淤积的影响。本发明为一种三维泥沙输移模拟技术，合理考虑高、低浓度泥沙条件下水沙相互作用机制和泥沙悬浮机制的泥沙运动数值模拟技术，对河口航道泥沙淤积等工程问题进行数值模拟研究，解析泥沙输移的作用过程和机制，为工程实践提供科技支持。

Description

一种高低浓度自适应的泥沙运动数值模拟方法

技术领域

本发明涉及泥沙运动数值模拟，特别涉及一种在同一海域水体中能有效地模拟同时存在低浓度泥沙和高浓度泥沙现象的泥沙运动数值模拟方法。

背景技术

海岸河口地区的动力条件复杂多变，人类活动频繁，海岸泥沙运动活跃。海岸河口泥沙运动与海岸的演变、航道与港池的淤积、污染物输移以及海岸工程建筑物的安全等问题有密切联系，其在海岸河口动力过程研究中占有重要地位，一直是港口海岸及近海工程学科中的研究热点。

目前，在国内外泥沙运动数值模拟研究中，已经形成开源软件、商业软件、自主软件等多种形式的众多的三维水沙输移数学模型，并且在海岸河口泥沙研究和工程泥沙实践中起到了重要的支撑作用。这些模型采用的数值方法、网格形式、理论模式各有不同，形成了各具特色、百花齐放的局面。但总体上看，对于海岸河口泥沙输移数值模拟问题，现有的泥沙运动数值模拟技术大多只适用于低浓度泥沙条件，而天然状态下，在同一海域水体中低浓度泥沙和高浓度泥沙现象是同时存在的。因此现有的方法很难合理反映高浓度泥沙条件下泥沙悬浮机制与输运规律，具体来说，还存在以下几方面的不足。

(1)现有的泥沙运动数值模拟技术大多没有考虑泥沙对床面切应力的影响，个别考虑该影响的方法在近底泥沙浓度垂向梯度大时存在计算稳定性等问题。

(2)现有的泥沙运动数值模拟技术没有合理考虑高浓度泥沙对水体紊动结构的影响。

发明内容

本发明的目的是克服现有技术中的不足，提供一种高低浓度自适应的泥沙运动数值模拟方法，本发明方法为一种三维泥沙输移模拟技术，合理考虑高、低浓度泥沙条件下水沙相互作用机制和泥沙悬浮机制的泥沙运动数值模拟技术，对河口航道泥沙淤积等工程问题进行数值模拟研究，解析泥沙输移的作用过程和机制，为工程实践提供科技支持。

本发明所采用的技术方案是：一种高低浓度自适应的泥沙运动数值模拟方法，包括如下步骤：

步骤一：建立高低浓度自适应的泥沙运动数学模型，所述模型的基本方程包括三维水流运动方程、三维泥沙输移扩散方程与床面冲淤变化方程；

步骤二：计算边界条件、悬沙沉速、床面切应力、稳定性函数及适应于高低浓度泥沙条件的修正；

步骤三：求解高低浓度自适应的泥沙运动数学模型；

步骤四：采用实测资料对泥沙运动数学模型进行水流验证，含沙量验证，航道淤积分布验证；

步骤五：根据泥沙运动数学模型结果，分析上游流量、外海潮差、泥沙沉降速度、航道深度、外海平均潮位差与高低浓度泥沙自适应技术对河口航道泥沙淤积的影响。

进一步的，步骤一中，所述模型的基本方程具体如下：

①三维水流运动方程：

式中，x,y,z为笛卡尔坐标系下的三维坐标，分别为东西方向、南北方向以及垂向的坐标；u,v分别为水平方向的东分量速度和北分量速度；w为垂向速度；g为重力加速度；t为时间；ρ₀为水体参考密度；ρ为水体密度；p为静水压力；h为水深；f为科氏参数；K_m为垂向涡粘系数；S_xx为x向波浪作用在x方向上的波浪辐射应力张量，S_xy为x向波浪作用在y方向上的波浪辐射应力张量，S_yx为y向波浪作用在x方向上的波浪辐射应力张量，S_yy为y向波浪作用在y方向上的波浪辐射应力张量；F_u为x方向的动量扩散系数，F_v为y方向的动量扩散系数；

②三维泥沙输移扩散方程：

式中，c为悬沙浓度；u_x，u_y和u_z分别为x，y和z方向的流速分量；ω_s为悬沙沉降速度；ε_x为x方向泥沙紊动扩散系数，ε_y为y方向泥沙紊动扩散系数；ε_z为垂向泥沙紊动扩散系数；

③床面冲淤变化方程：

式中，F为底床局部泥沙净通量；γ_s为泥沙干容重；z_b为计算时间T内单位长度上的冲淤强度。

进一步的，步骤二中，计算边界条件具体如下：

边界条件：

(1)表面边界条件

(2)床面边界条件

式中，z为笛卡尔坐标系下的垂向坐标；ε_z为垂向泥沙紊动扩散系数；c为悬沙浓度；E为冲刷通量；D为沉积通量；

冲刷通量E表示为：

式中，E₀为泥沙的床面冲刷强度，τ_b是床面切应力，τ_e是临界冲刷应力；

沉积通量D表示为：

式中，τ_d为临界淤积应力，ω_s为悬沙沉降速度。

进一步的，步骤二中，计算悬沙沉速具体如下：

悬沙沉速：

(1)粘性沙沉速公式如下：

式中，ω_s为悬沙沉降速度；ω_max为最大絮凝沉速；S为盐度；S_max为最佳絮凝盐度；c为悬沙浓度；

ε和υ分别为紊动能量耗散率和水体运动粘滞性系数；a、b、n、m、B₁和B₂均为经验系数，根据泥沙特性确定；温度对泥沙沉速的影响，主要通过不同温度条件下运动粘滞系数差异导致的单颗粒泥沙静水沉降速度差异来体现；

(2)非粘性沙沉速公式如下：

式中，ω_s0为在静水中的泥沙颗粒沉降速度；c_v为悬沙体积浓度；n₁为非粘性沙沉降速度制约系数；

ρ_s和ρ分别为泥沙和水的密度；

g为重力加速度；υ为水体运动粘滞性系数；d为泥沙颗粒粒径。

进一步的，步骤二中，计算床面切应力具体如下：

床面切应力：

式中，

为摩阻流速，ρ为水的密度，τ_b是床面切应力；u_b为近床面流速；C_d为底摩阻系数；

式中，κ为卡门常数，取0.41；z₀取底层网格单元高度的一半；z_r为底部粗糙率，取0.001m～0.01m。

进一步的，步骤二中，计算稳定性函数具体如下：

在近底泥沙浓度大于10kg/m³的情况下，为了保证方法的稳定性，将垂向各层的沉积通量分解为对流通量conv和耗散通量diss，具体如下：

式中，k为垂向单元层数；conv_k为第k层垂向单元的流通量；diss_k为第k层垂向单元的耗散通量；ω_s,k为第k层垂向单元的悬沙沉降速度；ω_s,k-1为第k-1层垂向单元的悬沙沉降速度；c_k为第k层垂向单元的悬沙浓度；c_k-1为第k-1层垂向单元的悬沙浓度；c_k+1为第k+1层垂向单元的悬沙浓度；c_k-2为第k-2层垂向单元的悬沙浓度；令a₁＝c_k+1-c_k，b₁＝c_k-1-c_k-2，则：

其中，

其中，指数q的取值为：

对于表层和底层沉积通量的计算，采用如下边界条件：

其中，kb为垂向层面数；c_-1为表层之外第-1层面的悬沙浓度；c₀为表层之外第0层面的悬沙浓度；c₁为表层的悬沙浓度；c_k+2为第k+2层面的悬沙浓度；ω_s,-1为表层之外第-1层面的悬沙沉降速度；ω_s,0为表层之外第0层面的悬沙沉降速度；ω_s,1为表层的悬沙沉降速度；ω_s,k+1为第k+1层面的悬沙沉降速度；ω_s,k+2为第k+2层面的悬沙沉降速度。

进一步的，步骤二中，计算适应于高低浓度泥沙条件的修正具体如下：

(1)泥沙对床面切应力的影响：

含沙水体摩阻流速：

式中，

为摩阻流速，τ_b是床面切应力；ρ_s和ρ分别为泥沙和水的密度，C_m为高浓度层内泥沙平均浓度，α₁为无量纲系数。

结合床面切应力τ_b和含沙水体摩阻流速u′_*，即反映泥沙对床面切应力的影响。

(2)泥沙对水体紊动结构的影响：

含沙水体垂向泥沙紊动扩散系数表示为掺混长度l_c和掺混速度w_mc的表达式：

式中，ε_z为垂向泥沙紊动扩散系数；ω_s为悬沙沉降速度；κ_s为泥沙掺混长度系数，根据实验确定；v′为水体垂向脉动速度，u′_*为含沙水体摩阻流速，h为水深；w_mc为掺混速度；z_z为坐标系下z坐标轴上数值，随位置变化而变化。

进一步的，步骤三具体为：以FVCOM为工具，结合步骤一的基本方程、步骤二的计算方程，采用有限体积法对所建立的高低浓度自适应的泥沙运动数学模型行求解，计算网格为三角形网格，网格尺度根据研究的区域而定。

本发明的有益效果是：

(1)本发明建立了高低浓度自适应的泥沙运动数学模型，可以用来模拟同时存在低浓度泥沙和高浓度泥沙现象的三维泥沙运动，能合理反映高、低浓度泥沙条件下泥沙悬浮机制与输运规律。

(2)本发明建立的全动力过程数学模型和三维泥沙数学模型范围覆盖了整个工程海域，且模型经过实测资料验证，计算结果合理可信，并利用该发明对河口航道泥沙淤积的影响因素进行了分析。

附图说明

图1为大模型和小模型计算范围示意图；

图2为2012年8月水文测点布置图；

图3为2012年8月全潮潮位验证图；

图3-a为鸡骨礁潮位站全潮潮位验证图；

图3-b为北槽中潮位站全潮潮位验证图；

图3-c为横沙潮位站全潮潮位验证图；

图3-d为长兴潮位站全潮潮位验证图；

图3-e为石洞口潮位站全潮潮位验证图；

图3-f为徐六泾潮位站全潮潮位验证图；

图4为2012年8月全潮潮流部分测点验证图；

图4-a为CS10S潮流站全潮流速验证图；

图4-b为CS10S潮流站全潮流向验证图；

图4-c为CS7S潮流站全潮流速验证图；

图4-d为CS7S潮流站全潮流向验证图；

图4-e为CSWS潮流站全潮流速验证图；

图4-f为CSWS潮流站全潮流向验证图；

图4-g为CS0S潮流站全潮流速验证图；

图4-h为CS0S潮流站全潮流向验证图；

图5为2012年8月全潮含沙量部分测点验证图；

图5-a为CS10S全潮表层处含沙量验证图；

图5-b为CS10S全潮0.2倍水深处含沙量验证图；

图5-c为CS10S全潮0.4倍水深处含沙量验证图；

图5-d为CS10S全潮0.6倍水深处含沙量验证图；

图5-e为CS10S全潮0.8倍水深处含沙量验证图；

图5-f为CS10S全潮底层水深处含沙量验证图；

图5-g为CS3N全潮表层水深处含沙量验证图；

图5-h为CS3N全潮0.2倍水深处含沙量验证图；

图5-i为CS3N全潮0.4倍水深处含沙量验证图；

图5-j为CS3N全潮0.6倍水深处含沙量验证图；

图5-k为CS3N全潮0.8倍水深处含沙量验证图；

图5-l为CS3N全潮底层水深处含沙量验证图；

图5-m为CSWS全潮表层水深处含沙量验证图；

图5-n为CSWS全潮0.2倍水深处含沙量验证图；

图5-o为CSWS全潮0.4倍水深处含沙量验证图；

图5-p为CSWS全潮0.6倍水深处含沙量验证图；

图5-q为CSWS全潮0.8倍水深处含沙量验证图；

图5-r为CSWS全潮底层水深处含沙量验证图；

图5-s为CS0S全潮表层水深处含沙量验证图；

图5-t为CS0S全潮0.2倍水深处含沙量验证图；

图5-u为CS0S全潮0.4倍水深处含沙量验证图；

图5-v为CS0S全潮0.6倍水深处含沙量验证图；

图5-w为CS0S全潮0.8倍水深处含沙量验证图；

图5-x为CS0S全潮底层水深处含沙量验证图；

图6为2012年全年航道常态淤积量验证图；

图7为本发明的整体思路流程图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步的描述。

本发明整体思路如图7所示，为一种三维泥沙输移模拟技术，合理考虑高、低浓度泥沙条件下水沙相互作用机制和泥沙悬浮机制的泥沙运动数值模拟技术，对河口航道泥沙淤积等工程问题进行数值模拟研究，解析泥沙输移的作用过程和机制，为工程实践提供科技支持。

实施例：

针对长江河口泥沙运动情况进行分析说明。

一种高低浓度自适应的泥沙运动数值模拟方法，包括以下步骤：

(一)建立数学模型：

建立高低浓度自适应的泥沙运动数学模型，所述模型的基本方程包括三维水流运动方程、三维泥沙输移扩散方程与床面冲淤变化方程。

①三维水流运动方程：

式中，x,y,z为笛卡尔坐标系下的三维坐标，分别为东西方向、南北方向以及垂向的坐标；u,v分别为水平方向的东分量速度和北分量速度；w为垂向速度；g为重力加速度；t为时间；ρ₀为水体参考密度；ρ为水体密度；p为静水压力；h为水深；f为科氏参数；K_m为垂向涡粘系数；S_xx为x向波浪作用在x方向上的波浪辐射应力张量，S_xy为x向波浪作用在y方向上的波浪辐射应力张量，S_yx为y向波浪作用在x方向上的波浪辐射应力张量，S_yy为y向波浪作用在y方向上的波浪辐射应力张量；F_u为x方向的动量扩散系数，F_v为y方向的动量扩散系数。

②三维泥沙输移扩散方程：

式中，c为悬沙浓度；u_x，u_y和u_z分别为x，y和z方向的流速分量；ω_s为悬沙沉降速度；ε_x为x方向泥沙紊动扩散系数，ε_y为y方向泥沙紊动扩散系数；ε_z为垂向泥沙紊动扩散系数；

③床面冲淤变化方程：

式中，F为底床局部泥沙净通量；γ_s为泥沙干容重；z_b为计算时间T内单位长度上的冲淤强度；

(二)关键问题的处理：

计算边界条件、悬沙沉速、床面切应力、稳定性函数及适应于高低浓度泥沙条件的修正。

①边界条件

(1)表面边界条件

(2)床面边界条件

式中，z为笛卡尔坐标系下的垂向坐标；ε_z为垂向泥沙紊动扩散系数；c为悬沙浓度；E为冲刷通量；D为沉积通量；

冲刷通量E表示为：

式中，E₀为泥沙的床面冲刷强度，τ_b是床面切应力，τ_e是临界冲刷应力；

沉积通量表示为：

式中，τ_d为临界淤积应力。

②悬沙沉速：

泥沙沉速与泥沙性质有关，泥沙分为粘性泥沙和非粘性泥沙，本发明针对不同的泥沙，给出了对应的沉速计算公式。

(1)粘性细颗粒泥沙沉降的形成既有化学原因(属于胶体化学性质，主要为电化学性质)，又有物理原因(布朗运动、不等速沉降、水流紊动)。国内外学者对沉速进行了大量研究，提出的沉速计算公式也会因考虑影响因素以及研究方法的不同而不同，计算沉速公式通常涉及絮凝沉降段和制约沉降段。影响悬浮泥沙絮凝沉降速度的因素很多，主要有：泥沙浓度、水流紊动、盐度、含离子浓度、电解质浓度、阳离子化合价、颗粒粒径、絮团强度、分形结构、泥沙组成、絮团形成时间、PH值、有机质等。例如针对长江口粘性细颗粒泥沙沉降规律，主要考虑含沙浓度、盐度、水流紊动、温度因素。

本发明综合Deflt3D、Hwang(1989)以及Van Leussen(1994)相关研究，构建粘性细颗粒泥沙沉降速度公式如式(11)所示：

式中，ω_s为悬沙沉降速度；ω_max为最大絮凝沉速；S为盐度；S_max为最佳絮凝盐度；c为悬沙浓度；

ε和υ分别为紊动能量耗散率和水体运动粘滞性系数；a、b、n、m、B₁和B₂均为经验系数，根据泥沙特性确定；温度对泥沙沉速的影响，主要通过不同温度条件下运动粘滞系数差异导致的单颗粒泥沙静水沉降速度差异来体现；

(2)非粘性泥沙沉降主要与含沙量有关，针对非粘性泥沙，Richardson和Zaki(1954)曾提出沉降速度与含沙量之间的关系式：

式中，ω_s0为在静水中的泥沙颗粒沉降速度；c_v为悬沙体积浓度；n₁为非粘性沙沉降速度制约系数；

对于上式中制约系数的确定。Richardson和Zaki(1954)认为其与颗粒雷诺数有关，随着颗粒雷诺数的减小而增大。Cheng(1997)提出了一个与雷诺数和悬沙体积浓度有关的公式：

式中，

ρ_s和ρ分别为泥沙和水的密度；

g为重力加速度；υ为水体运动粘滞性系数；d为泥沙颗粒粒径。

③床面切应力：

(1)当水体中不含泥沙时，床面切应力用式(14)计算：

式中，

为摩阻流速，ρ为水的密度，τ_b是床面切应力；u_b为近床面流速；C_d为底摩阻系数；

式中，κ为卡门常数，取0.41；z₀取底层网格单元高度的一半；z_r为底部粗糙率，取0.001m～0.01m。

(2)当水体中含有泥沙时，床面切应力计算公式形式与式(14)保持一致，但其中的摩阻流速发生改变，含沙水体摩阻流速按式(22)计算。

④稳定性函数

近底通量放映了床面泥沙与上层水体的物质交换程度，即三维泥沙运动数学模型中的床面边界条件，包括床面冲刷和沉积两方面。以国际通用的切应力模式为例，近底通量与床面切应力及泥沙特征有关，由式(9)和(10)确定。

在本实施例中，在近底泥沙浓度大于10kg/m³的情况下，为了保证方法的稳定性，将垂向各层的沉积通量分解为对流通量conv和耗散通量diss，具体如下：

式中，k为垂向单元层数；conv_k为第k层垂向单元的流通量；diss_k为第k层垂向单元的耗散通量；ω_s,k为第k层垂向单元的悬沙沉降速度；ω_s,k-1为第k-1层垂向单元的悬沙沉降速度；c_k为第k层垂向单元的悬沙浓度；c_k-1为第k-1层垂向单元的悬沙浓度；c_k+1为第k+1层垂向单元的悬沙浓度；c_k-2为第k-2层垂向单元的悬沙浓度；令a₁＝c_k+1-c_k，b₁＝c_k-1-c_k-2，则：

其中，

其中，指数q的取值为：

对于表层和底层沉积通量的计算，采用如下边界条件：

其中，kb为垂向层面数；c_-1为表层之外第-1层面的悬沙浓度；c₀为表层之外第0层面的悬沙浓度；c₁为表层的悬沙浓度；c_k+2为第k+2层面的悬沙浓度；ω_s,-1为表层之外第-1层面的悬沙沉降速度；ω_s,0为表层之外第0层面的悬沙沉降速度；ω_s,1为表层的悬沙沉降速度；ω_s,k+1为第k+1层面的悬沙沉降速度；ω_s,k+2为第k+2层面的悬沙沉降速度。

⑤适应于高低浓度泥沙条件的修正

(1)泥沙对床面切应力的影响：

含沙水体摩阻流速：

式中，

为摩阻流速，τ_b是床面切应力；ρ_s和ρ分别为泥沙和水的密度，C_m为高浓度层内泥沙平均浓度，α₁为无量纲系数，本发明中取α₁＝70。

结合床面切应力τ_b和含沙水体摩阻流速u′_*，即反映泥沙对床面切应力的影响。

(2)泥沙对水体紊动结构的影响：

本实施例中，为了使泥沙垂向结构同时适应于高低浓度泥沙条件，采用改进的垂向泥沙紊动扩散系数，与传统的垂向泥沙紊动扩散系数相比，改进的系数中引入了掺混长度和掺混速度，以同时满足高、低浓度泥沙条件。

改进的垂向泥沙紊动扩散系数计算方法如式(23)所示：

式中，ε_z为垂向泥沙紊动扩散系数；ω_s为悬沙沉降速度；l_c为掺混长度；w_mc为掺混速度。

本实施例中分别采用Nezu和Nakagawa(1993)的明渠水流紊动强度分布公式和卡门紊流相似假说来表示垂向掺混速度w_mc和掺混长度l_c：

式中，κ_s为泥沙掺混长度系数，需根据实验来确定；v′为水体垂向脉动速度，u′_*为含沙水体摩阻流速，h为水深，z_z为坐标系下z坐标轴上数值，随位置变化而变化。

(三)模型的建立与求解：

以FVCOM为工具，结合步骤一的基本方程、步骤二的计算方程，采用有限体积法对所建立的高低浓度自适应的泥沙运动数学模型行求解，计算网格为三角形网格，网格尺度根据研究的区域而定。

本实施例采用大、小模型嵌套的计算模式，大模型主要用来调试外海潮汐边界和向小模型提供长江口上游流量边界。大模型和小模型采用相同的外海边界范围，东至东经125°，北至江苏盐城(北纬34°)，南至浙江温州(北纬28°)，模型外海海域南北宽约667km，从长江口外至东边界约280km。大模型长江口上游边界至安徽省大通水文观测站(距离徐六泾约500km)，小模型上游边界至江苏省江阴(距离徐六径约95km，距北槽下口约137km)，详见图1。

采用有限体积法求解数学模型。计算网格水平方向为非结构化的三角形网格，垂向采用σ坐标，最大网格尺度10000m，位于外海开边界处，最小网格尺度20m。

模型外海边界由Mike全球潮汐预报系统提供；大模型上游边界为流量边界，直接输入大通水文站实测流量数据；小模型上游边界同样为流量边界，由大模型计算结果提供流量数据。

陆域边界位置根据提供的工程布置图、卫星遥感影像、海图、航道图等资料确定。

(四)模型的验证：

(1)水流验证

选用长江口在2012年8月12日～8月19日期间的实测全潮数据对水流数学模型进行验证，潮位测站和潮流测站位置如图2所示。各潮位测站潮位过程验证结果见图3，各潮流测站的流速流向过程验证结果见图4。图中“_◇”表示实测值，实线表示计算值。验证结果表明：模型计算值与实测值吻合良好，说明率定工作良好，三维水动力数学模型合理反映了长江口深水航道及其附近海域的潮流运动情况。

(2)含沙量验证

含沙量验证时间和测点布置与水动力一致。本实施例比较了2012年8月12日-8月19日长江口深水航道北槽中各测点的垂向分层含沙量历时曲线计算和实测结果，如图5所示。从图中可知，长江口深水航道沿程各含沙量测点的含沙量在数值上有很大差异，有的时段含沙量低于0.5kg/m³，而有的时段含沙量高于20kg/m³，本实施例利用高低浓度自适应的泥沙运动数值模拟方法对上述现象进行了很好的模拟，验证结果表明：计算值与实测值吻合良好，说明了含沙量验证工作良好，三维泥沙运动数学模型合理反映了长江口深水航道及其附近海域泥沙场的变化。

(3)航道淤积分布验证

长江口南港-北槽深水航道2012年全年常态淤积量计算值与实测值比较见图6。可见，模型计算得到的淤积分布基本反映了长江口北槽深水航道淤积分布趋势，从全航道淤积量看，计算结果与实测结果吻合较好，全年航道实测常态回淤总量为10390万方，计算结果为9610万方，与实测结果相差-7.5％。

(五)航道淤积影响因素分析

本实施例中采用本发明计算分析了上游流量、外海潮差、泥沙沉降速度、航道深度(地形变化)、外海平均潮位差与高低浓度泥沙自适应技术六方面对航道回淤的影响。计算条件及航道回淤量变化见表1。从表中可知：高低浓度泥沙自适应技术、泥沙沉降速度、航道深度(地形变化)、外海潮差对航道淤积量起到重要影响，其中前两项因素导致的淤积量变化可达30％以上，可见，合理考虑两项因素，会很大程度提高数值模拟方法的准确性。上述两项因素也是本发明的关键点和创新点。

表1不同因素对航道回淤量的影响

尽管上面结合附图对本发明的优选实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，并不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护范围之内。

Claims (6)

1.一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，包括以下步骤：

步骤一：建立高低浓度自适应的泥沙运动数学模型，所述模型的基本方程包括三维水流运动方程、三维泥沙输移扩散方程与床面冲淤变化方程；

步骤二：计算边界条件、悬沙沉速、床面切应力、稳定性函数及适应于高低浓度泥沙条件的修正；

其中，计算稳定性函数具体如下：

在近底泥沙浓度大于10kg/m³的情况下，为了保证方法的稳定性，将垂向各层的沉积通量分解为对流通量conv和耗散通量diss，具体如下：

式中，k为垂向单元层数；conv_k为第k层垂向单元的流通量；diss_k为第k层垂向单元的耗散通量；ω_s,k为第k层垂向单元的悬沙沉降速度；ω_s,k-1为第k-1层垂向单元的悬沙沉降速度；c_k为第k层垂向单元的悬沙浓度；c_k-1为第k-1层垂向单元的悬沙浓度；c_k+1为第k+1层垂向单元的悬沙浓度；c_k-2为第k-2层垂向单元的悬沙浓度；令a₁＝c_k+1-c_k，b₁＝c_k-1-c_k-2，则：

其中，

其中，指数q的取值为：

对于表层和底层沉积通量的计算，采用如下边界条件：

其中，kb为垂向层面数；c_-1为表层之外第-1层面的悬沙浓度；c₀为表层之外第0层面的悬沙浓度；c₁为表层的悬沙浓度；c_k+2为第k+2层面的悬沙浓度；ω_s,-1为表层之外第-1层面的悬沙沉降速度；ω_s,0为表层之外第0层面的悬沙沉降速度；ω_s,1为表层的悬沙沉降速度；ω_s,k+1为第k+1层面的悬沙沉降速度；ω_s,k+2为第k+2层面的悬沙沉降速度；

计算适应于高低浓度泥沙条件的修正具体如下：

(1)泥沙对床面切应力的影响：

含沙水体摩阻流速：

式中，

为摩阻流速，τ_b是床面切应力；ρ_s和ρ分别为泥沙和水的密度，C_m为高浓度层内泥沙平均浓度，α₁为无量纲系数；

结合床面切应力τ_b和含沙水体摩阻流速u′_*，即反映泥沙对床面切应力的影响；

(2)泥沙对水体紊动结构的影响：

含沙水体垂向泥沙紊动扩散系数表示为掺混长度l_c和掺混速度w_mc的表达式：

式中，ε_z为垂向泥沙紊动扩散系数；ω_s为悬沙沉降速度；κ_s为泥沙掺混长度系数，根据实验确定；v′为水体垂向脉动速度，u′_*为含沙水体摩阻流速，h为水深；w_mc为掺混速度；z_z为坐标系下z坐标轴上数值，随位置变化而变化；

步骤三：求解高低浓度自适应的泥沙运动数学模型；

步骤四：采用实测资料对泥沙运动数学模型进行水流验证，含沙量验证，航道淤积分布验证；

步骤五：根据泥沙运动数学模型结果，分析上游流量、外海潮差、泥沙沉降速度、航道深度、外海平均潮位差与高低浓度泥沙自适应技术对河口航道泥沙淤积的影响。

2.根据权利要求1所述的一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，步骤一中，所述模型的基本方程具体如下：

①三维水流运动方程：

式中，x,y,z为笛卡尔坐标系下的三维坐标，分别为东西方向、南北方向以及垂向的坐标；u,v分别为水平方向的东分量速度和北分量速度；w为垂向速度；g为重力加速度；t为时间；ρ₀为水体参考密度；ρ为水体密度；p为静水压力；h为水深；f为科氏参数；K_m为垂向涡粘系数；S_xx为x向波浪作用在x方向上的波浪辐射应力张量，S_xy为x向波浪作用在y方向上的波浪辐射应力张量，S_yx为y向波浪作用在x方向上的波浪辐射应力张量，S_yy为y向波浪作用在y方向上的波浪辐射应力张量；F_u为x方向的动量扩散系数，F_v为y方向的动量扩散系数；

②三维泥沙输移扩散方程：

式中，c为悬沙浓度；u_x，u_y和u_z分别为x，y和z方向的流速分量；ω_s为悬沙沉降速度；ε_x为x方向泥沙紊动扩散系数，ε_y为y方向泥沙紊动扩散系数；ε_z为垂向泥沙紊动扩散系数；

③床面冲淤变化方程：

式中，F为底床局部泥沙净通量；γ_s为泥沙干容重；z_b为计算时间T内单位长度上的冲淤强度。

3.根据权利要求1所述的一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，步骤二中，计算边界条件具体如下：

边界条件：

(1)表面边界条件

(2)床面边界条件

式中，z为笛卡尔坐标系下的垂向坐标；ε_z为垂向泥沙紊动扩散系数；c为悬沙浓度；E为冲刷通量；D为沉积通量；

冲刷通量E表示为：

式中，E₀为泥沙的床面冲刷强度，τ_b是床面切应力，τ_e是临界冲刷应力；

沉积通量D表示为：

式中，τ_d为临界淤积应力，ω_s为悬沙沉降速度。

4.根据权利要求1所述的一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，步骤二中，计算悬沙沉速具体如下：

悬沙沉速：

(1)粘性沙沉速公式如下：

式中，ω_s为悬沙沉降速度；ω_max为最大絮凝沉速；S为盐度；S_max为最佳絮凝盐度；c为悬沙浓度；

ε和υ分别为紊动能量耗散率和水体运动粘滞性系数；a、b、n、m、B₁和B₂均为经验系数，根据泥沙特性确定；温度对泥沙沉速的影响，主要通过不同温度条件下运动粘滞系数差异导致的单颗粒泥沙静水沉降速度差异来体现；

(2)非粘性沙沉速公式如下：

式中，ω_s0为在静水中的泥沙颗粒沉降速度；c_v为悬沙体积浓度；n₁为非粘性沙沉降速度制约系数；

ρ_s和ρ分别为泥沙和水的密度；

g为重力加速度；υ为水体运动粘滞性系数；d为泥沙颗粒粒径。

5.根据权利要求1所述的一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，步骤二中，计算床面切应力具体如下：

床面切应力：

式中，

为摩阻流速，ρ为水的密度，τ_b是床面切应力；u_b为近床面流速；C_d为底摩阻系数；

式中，κ为卡门常数，取0.41；z₀取底层网格单元高度的一半；z_r为底部粗糙率，取0.001m～0.01m。

6.根据权利要求1所述的一种高低浓度自适应的泥沙运动数值模拟方法，其特征在于，步骤三具体为：以FVCOM为工具，结合步骤一的基本方程、步骤二的计算方程，采用有限体积法对所建立的高低浓度自适应的泥沙运动数学模型行求解，计算网格为三角形网格，网格尺度根据研究的区域而定。

Images