CN102359862A - 粉沙质和淤泥质海岸泥沙运动数值模拟方法 - Google Patents

Abstract

本发明公开了一种粉沙质和淤泥质海岸泥沙运动数值模拟方法。主要步骤包括：步骤一、建立基本方程，所述基本方程包括水流连续方程、水流运动方程、悬移质输移扩散方程、床面冲淤变化方程；步骤二，计算悬沙沉速、泥沙冲淤函数、波流共同作用下的床面剪切应力；步骤三，泥沙数学模型的建立和求解；步骤四，模型验证，包括水流验证、含沙量验证；步骤五，根据泥沙数学模型判定悬沙沉速场。本发明在考虑淤泥质和粉沙质泥沙特性的基础上，构建泥沙数学模型，可以同时模拟两种不同性质泥沙运动；不仅如此，利用所建立的泥沙数学模型还可以计算悬沙的沉速场，据此分析研究区域的特征和影响因素。

Description

粉沙质和淤泥质海岸泥沙运动数值模拟方法

技术领域

本发明涉及一种粉沙质和淤泥质海岸泥沙运动数值模拟方法，属于泥沙运动数值模拟的技术领域。

背景技术

近年来，随着我国经济实力迅速提升，港口航道的建设和发展面临着巨大的需求和机遇。自然条件好的岸线已经得到充分开发，港口航道工程的建设重点逐渐转向淤泥质海岸和粉沙质海岸。由于坡缓水浅，动力复杂，淤泥质海岸和粉沙质海岸的泥沙回淤问题较为突出，给港口航道工程的建设和维护带来了一定的困难，“浅水深用”也成为工程界主要关心的问题。

目前，国内外对于淤泥质浅滩泥沙运动的研究较多，建港技术也相对成熟，如连云港港口等。至于粉沙质海岸，在波浪、潮流等动力作用下，泥沙易悬易沉，运动十分活跃，航道经常会出现大风骤淤的现象，是工程界关注的热点。而对于同时包括淤泥质和粉沙质泥沙的海域，由于水沙环境更为复杂，国内外相关成果很少，有待于进一步研究。

例如，位于江苏连云港的埒子口～灌河口海域，海床底质较为复杂，其中，埒子口以西海域的泥沙颗粒较细，为典型的淤泥质海岸；灌河口附近海域的底床泥沙粒径略粗，为粉沙质海岸。因此，对于连云港大范围海域来说，由于淤泥质海岸和粉沙质海岸并存，水沙运动十分复杂，是海岸泥沙运动研究的难点。

发明内容

本发明所要解决的技术问题是在考虑淤泥质和粉沙质泥沙特性的基础上，构建泥沙数学模型，模拟计算粉沙质和淤泥质同时存在海域的泥沙运动，并利用所建立的泥沙数学模型计算悬沙的沉速场。

本发明为解决上述技术问题采用以下技术方案：

一种粉沙质和淤泥质海岸泥沙运动数值模拟方法，其特征在于，包括如下步骤：

步骤一：建立基本方程，所述基本方程包括水流连续方程、水流运动方程、悬移质输移扩散方程、床面冲淤变化方程；具体如下：

①水流连续方程：

式中：x、y为直角坐标系坐标；t为时间；h为平均水深；ζ为相对于平均海平面的潮位；U_x、U_y分别为x、y方向上的垂线平均速度；

②水流运动方程：

$\frac{\∂U_x}{\∂t}+U_x\frac{{\∂U}_x}{\∂x}+U_y\frac{\∂U_x}{\∂y}={\mathrm{fU}}_y-g\frac{\∂\mathrm{\ζ}}{\∂x}-\frac{{\mathrm{\τ}}_x}{\mathrm{\ρ}(h+\mathrm{\ζ})}+N_x(\frac{{\∂}^2{U}_x}{\∂x^2}+\frac{{\∂}^2{U}_x}{\∂y^2})$

$\frac{\∂U_y}{\∂t}+U_x\frac{{\∂U}_y}{\∂x}+U_y\frac{\∂U_y}{\∂y}={-\mathrm{fU}}_x-g\frac{\∂\mathrm{\ζ}}{\∂y}-\frac{{\mathrm{\τ}}_y}{\mathrm{\ρ}(h+\mathrm{\ζ})}+N_y(\frac{{\∂}^2{U}_y}{\∂x^2}+\frac{{\∂}^2{U}_y}{\∂y^2})---\left(2\right)$

式中：ρ为水体密度；g为重力加速度；N_x、N_y分别为x、y方向的水平紊动粘性系数；f为科氏参数；τ_x、τ_y分别为波流共同作用下床面剪切应力在x、y方向的分量；

③悬移质输移扩散方程：

式中：S为垂线平均含沙量；D_x、D_y分别为x、y方向的泥沙扩散系数；F_s为泥沙冲淤函数；

④床面冲淤变化方程：

${\mathrm{\γ}}_d\frac{{\∂\mathrm{\η}}_b}{\∂t}-F_s=0---\left(4\right)$

式中：γ_d为床沙干密度，η_b为海底床面的竖向位移，即冲淤变化量；F_s为泥沙冲淤函数。

步骤二，计算悬沙沉速、泥沙冲淤函数、波流共同作用下的床面剪切应力；具体如下：

(1)悬沙沉速：

针对中值粒径小于63μm和含沙量小于10kg/m³的泥沙，关于沉速与含沙量之间的关系式：

ω＝kS^α                            (5)

式中k、a为经验系数，ω为泥沙沉降速度，S为含沙量；

(2)泥沙冲淤函数：

泥沙冲淤函数F_s与底部剪切应力及泥沙特征的关系由下式确定：

式中：τ为水流底部剪切应力；τ_d为不淤临界剪切应力；τ_e为起动临界剪切应力；α为淤积概率；M为冲刷系数；ω为泥沙沉降速度；

(3)波流共同作用下的床面剪切应力：

(3)-a，当只考虑水流作用时，床面平均剪切应力：

${\mathrm{\τ}}_c=\frac{1}{2}{\mathrm{\ρf}}_c{U}_c^2$

式中：U_c为水流平均速度；f_c为水流摩阻系数，可用下式计算

$f_c=2{\left\{2.5\right[\ln \left(\frac{30h}{k}\right)-1\left]\right\}}^{-2}---\left(7\right)$

其中k为床面阻力系数，h为平均水深；

(3)-b，当只考虑波浪作用时，床面平均剪切应力：

${\mathrm{\τ}}_w=\frac{1}{2}{\mathrm{\ρf}}_w{U}_b^2---\left(8\right)$

式中：U_b为波浪底部水质点水平运动速度；f_w为波浪摩阻系数，其中

$\begin{array}{cc}{f}_w=0.00251\exp \left[{5.21(A/k_s)}^{-0.19}\right]& A/k_s>1.57\\ =0.3& A/k_s\≤1.57\end{array}---\left(9\right)$

其中k_s为Nikuradse糙率，取30z₀；z₀为波浪摩阻高度；A为波浪近底质点运动的振幅，按照线性波理论近似表示为：

$A=\frac{H}{\mathrm{\π}}\frac{1}{\mathrm{sh}(2\mathrm{\πh}/L)}---\left(10\right)$

式中，H代表波浪高度，L代表波浪的长度；

(3)-c，当波流共同作用，采用波浪摩阻系数概念

$\frac{{\mathrm{\τ}}_{\mathrm{mean}}}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}(1+b{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^p{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^p)---\left(11\right)$

$\frac{{\mathrm{\τ}}_{\max }}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=1+a{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^m{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^n---\left(12\right)$

式中：τ_c为考虑纯水流作用下的切应力，τ_w为考虑纯波浪作用的切应力，b、p、q、a、m、n为综合表达式。

步骤三，泥沙数学模型的建立和求解：以MIKE21为工具，结合步骤一的基本方程、步骤二的计算结果建立泥沙数学模型，然后采用有限体积法求解泥沙数学模型，计算网格为三角网格，最大网格尺度1500m，最小20m；

步骤四，泥沙数学模型验证：

4-1，采用实测中潮资料对泥沙数学模型进行水流验证；

4-2，采用平均含沙量作为检验的标准，相对应的动力条件为典型中潮资料加上年均有效波高波浪场，对泥沙数学模型进行含沙量验证；

步骤五，根据泥沙数学模型判定悬沙沉速场：

5-1，根据悬沙沉速来判断泥沙性质：当沉速大于1mm/s的泥沙属于粉沙质泥沙，沉速在0.6mm/s以下的泥沙判定为淤泥质泥沙，而既包括粉沙质泥沙又包括淤泥质泥沙的粉沙淤泥质泥沙，其沉速范围为0.6mm/s～1mm/s；

5-2，根据模型图中潮周期内不同沉速线的外包络线判断悬沙沉速场与底质分布。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)建立了波流作用下的泥沙数学模型，可以用来同时模拟包含粉沙质和淤泥质在内的泥沙运动。

(2)本发明建立的二维潮流数学模型和二维泥沙数学模型范围覆盖了整个工程海域，且模型经过实测资料验证，计算结果合理可信，利用所建立的泥沙数学模型计算悬沙的沉速场，并由此分析不同性质泥沙的分布特征。

附图说明

图1是2005年9月现场实测海床底质中值粒径分布图；

图2是模型范围示意图；

图3是2005年9月水文测点布置图；

图4是2005年9月中潮潮位验证图；

图5是2005年9月中潮潮流部分测点验证图；

图6是年均含沙量计算结果(kg/m3)；

图7是常风天沉速计算结果(m/s)；

图8是本发明的模拟方法流程图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

实施例：针对连云港埒子口～灌河口海域浅滩泥沙运动情况进行分析说明。

首先如图8所示，本发明的模拟方法流程如下：

(1)建立数学模型：

建立水流连续方程、水流运动方程、悬移质输移扩散方程以及床面充淤变化方程等，如式(1)、(2)、(3)、(4)所示。

(2)关键问题的处理：

沉速的计算

悬沙沉速不仅是影响含沙量计算的重要参数，也是反映水体中泥沙性质的依据。根据以往研究，含沙量是影响沉速的最关键因素。针对中值粒径小于63μm和含沙量小于10kg/m³的泥沙，Burt(1986)曾提出关于沉速与含沙量之间的关系式：

ω＝kS^α                (5)

式中k、α为经验系数，ω_s为沉速，S为含沙量

Shi(2010)曾根据长江口悬沙颗粒分析，认为悬沙中大部分颗粒(95％左右)小于63μm，而大于63μm的泥沙基本上以推移质或者近底高含沙水体的形态运移，并根据Rouse含沙量公式率定计算了长江口絮凝体的沉速为ω＝2.37S^0.84，取得了较好的结果，计算所得沉速为0.4～4.1mm/s，而这一沉速范围几乎包括了淤泥和部分粉沙的沉速范围。

虽然连云港大范围海域既有淤泥质海岸，又有粉沙质海岸，但该海域悬沙颗粒分析表明，水体中粒径小于0.063mm的泥沙占90％以上，与长江口流域泥沙情况十分类似，因此可以应用Burt公式计算连云港海域大范围悬沙沉速，k、α参数的取值需要在计算中通过率定给出。

底部冲淤函数

近底通量反映了底床泥沙与上层水体的物质交换程度，也是泥沙运动数值模拟中的底部边界条件，包括床面冲刷(erosion)和淤积(deposition)两个方面。在近底通量的描述方式上，国际范围内主要采用以床面切应力为判据的模式。其中，底部冲淤函数Fs与底部剪切应力及泥沙特征有关，由式(6)确定。

式中：τ为水流底部剪切应力；τ_d为不淤临界剪切应力；τ_e为起动临界剪切应力；α为淤积概率；M为冲刷系数；ω为泥沙沉降速度。

波流共同作用下的床面剪切应力

床面剪应力即可由潮流产生，也可由波浪产生，可以分别加考虑。

(a)当只考虑水流作用时，床面平均剪切应力用式(7)计算：

${\mathrm{\τ}}_c=\frac{1}{2}{\mathrm{\ρf}}_c{U}_c^2$

式中：U_c为水流平均速度；f_c为水流摩阻系数，可用下式计算

$f_c=2{\left\{2.5\right[\ln \left(\frac{30h}{k}\right)-1\left]\right\}}^{-2}---\left(7\right)$

其中k为床面阻力系数，h为水深。

当只考虑潮流对泥沙起动的影响时，床面阻力对泥沙起动的影响最为主要，不仅影响着水流的计算，也与床面剪切应力的计算有关。床面阻力包括沙粒阻力和沙波阻力，沙粒阻力系床面沙粒的摩阻而引起，主要受粒径大小的影响，而沙波阻力属于形状阻力；由于粉沙质海岸床面形态较为复杂，床面阻力明显不同于淤泥质海岸，具体参数需通过率定得出。

(b)当考虑纯波浪作用时，床面平均剪切应力以式(8)计算。

${\mathrm{\τ}}_w=\frac{1}{2}{\mathrm{\ρf}}_w{U}_b^2---\left(8\right)$

式中：U_b为波浪底部水质点水平运动速度；f_w为波浪摩阻系数，采用Swart(1974)提出的表达式：

$\begin{array}{cc}{f}_w=0.00251\exp \left[{5.21(A/k_s)}^{-0.19}\right]& A/k_s>1.57\\ =0.3& A/k_s\≤1.57\end{array}---\left(9\right)$

其中k_s为Nikuradse糙率，一般取为30z₀；z₀为波浪摩阻高度；A为波浪近底质点运动的振幅，可以按照线性波理论近似表示为：

$A=\frac{H}{\mathrm{\π}}\frac{1}{\mathrm{sh}(2\mathrm{\πh}/L)}---\left(10\right)$

(c)当波流共同作用，以式(11)和式(12)计算。

$\frac{{\mathrm{\τ}}_{\mathrm{mean}}}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}(1+b{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^p{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^p)---\left(11\right)$

$\frac{{\mathrm{\τ}}_{\max }}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=1+a{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^m{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^n---\left(12\right)$

式中：τ_c为考虑纯流作用下的切应力，τ_w为考虑纯波浪作用的切应力，b、p、q、a、m、n为综合表达式。

(3)模型的建立和求解：

数学模型北起日照(35°22′30″N，119°33′E)，东至(35°22′30″N，120°17′E)，南至废黄河口附近(34°17′00″N，120°17′E)，东西宽约99.7km，南北长约119.3km，模型范围内水域面积约8648km²，详见图2。模型还包括部分灌河河段，上游至响水大桥附近。

采用有限体积法(Finite Volume Method)求解数学模型。计算网格为三角网格，最大网格尺度1500m，最小20m。

数学模型中的外海边界条件由东中国海潮波模型提供；灌河边界位于响水附近，由灌河长河段模型提供。

(4)模型的验证：

(a)水流验证

采用2005年9月实测中潮资料对于水流数学模型进行验证，流速测点和潮位站位置详见图3。各潮位站的潮位过程验证结果见图4，潮位基面为当地平均海平面，各潮流测点的流速过程验证结果见图5。图中“o”表示为实测值，实线代表计算结果；代表流向的角度以北向为度量起点，顺时针方向为正。验证结果表明：连云港大范围海域潮流计算结果与实测资料吻合较好，能够反映连云港大范围海域潮流运动情况，如图5所示。

(b)含沙量验证

对于淤泥质海岸来说，泥沙运动的机理可以概括为“风浪掀沙、潮流输沙”；而对粉沙质海岸而言，泥沙运动的机理是“波流共同掀沙，潮流输沙”。但风浪因素往往具有很强的随机性，这就造成了同一地点在不同时间段内的含沙量值可能会相差较大。因此，采用多年平均含沙量作为检验的标准，相对应的动力条件为典型中潮加上年均有效波高波浪场。在连云港海域，有关大范围泥沙运动趋势已有较为清晰的认识，而且大西山、临洪河口、羊山岛和灌河口都有较多的实测资料，可以用来对于泥沙数学模型进行验证。

图6为常风天情况下含沙量场的计算结果，表1则为主要站点的含沙量比较。图6的结果表明，计算所得的含沙量等值线与等深线走势相近，近岸水体含沙量较高，向外海逐渐降低，外海含沙量大致在0.1kg/m³左右，这种横向分布形式与金镠(1990)和刘家驹(1990)等对连云港大范围含沙量场的研究结果是一致的；灌河口海域含沙量较高，并呈向西北方向扩散的趋势，这也与以前的研究成果相符。灌河口海域是高含沙区域，口内含沙量较高，口外高含沙水体沿着西水道向埒子口海域逐渐扩散。从表1可以看出，几个主要验证站点的计算含沙量与实测值吻合较好，特别是灌河口内陈家港的含沙量超过4.0kg/m³，而连云港海域洋山岛附近的含沙量仅0.2kg/m³左右，两者相差很大，能够在一个数学模型中同时进行模拟实属不易。

表1主要站点年均含沙量验证情况(kg/m³)

(5)悬沙沉速场：

利用所建立的泥沙数学模型，不仅可以计算悬沙的含沙量场，还可以计算悬沙的沉速场，并由此分析不同性质泥沙的分布特征。

悬沙特性与沉速

粉沙质海岸的特点：(1)从岸坡角度分析，粉沙质海岸岸坡变幅很大，一般岸坡i＝1/3000～1/50(李孟国，曹祖德，2009)，既不同于淤泥质海岸(i＝1/1000～1/2000)，也与沙质海岸(i＞1/500)有所区别；(2)从泥沙运移型态来看，粉沙质泥沙的基本特性是易悬、易沉、易板结，既有悬移质，也有推移质，还有混移质(近底高含沙水体)，但是泥沙运动仍以悬移质为主(刘家驹，2004)。

悬沙沉速是能够反映水体中泥沙性质的重要参数：孔令双等(2004)认为对于淤泥质海岸，不管泥沙粒径多大，最终都以絮凝速度沉降，絮凝沉速范围一般0.45～0.55mm/s；杨华(2008)曾分析泥沙粒径与沉速的关系，如表2所示，结果表明在泥沙颗粒大于0.03mm时，絮凝现象已不再明显，中值粒径大于0.03mm的泥沙沉速在1mm/s以上。而对于粉沙质海岸的泥沙，泥沙中值粒径一般都大于0.03mm，由于絮凝作用较弱，其沉速也比较大，例如黄骅港的泥沙沉降速度为1.20mm/s，如东港的为2.42mm/s，一般都要大于1.0mm/s。

既然淤泥质泥沙和粉沙质泥沙的沉速是不同的，因此，也可以通过悬沙沉速来判断泥沙性质：当沉速大于1mm/s的泥沙属于粉沙质泥沙，沉速在0.6mm/s以下的泥沙可以判定为淤泥质泥沙，而既包括粉沙质泥沙又包括淤泥质泥沙的粉沙淤泥质泥沙，其沉速范围一般为0.6mm/s～1mm/s。

表2泥沙粒径与沉速的关系

沉速场与底质分布

图7为1个潮周期内0.6mm/s和1.0mm/s等沉速线的外包络线，由图可见：计算的悬沙沉速场与海床底质分布图(见图1)有一定的相似之处，特别是在埒子口～灌河口海域，沉速场与底质分布走向几乎完全一致，说明两者存在一定的内在联系：埒子口以西海域，悬沙沉速小于0.6mm/s，多以絮凝沉速沉降，属于淤泥质海岸；埒子口以东海域，悬沙沉速大于1mm/s，多以粉沙质泥沙沉速沉降，属于粉沙质海岸；埒子口海域悬沙沉速在0.6mm/s～1mm/s之间，与该海域处于粉沙质海岸向淤泥质海岸的过渡特征相一致。

灌河为潮汐河流，一般情况下涨落潮水流主要从西水道进出灌河，也就是说，涨潮时，外海水流主要由口门以西海域沿西水道进入灌河；落潮时，灌河内的水流同样由西水道返回西侧海域。由于灌河的含沙量较高，且多为粉沙质，因此，在落潮过程中，灌河口内高含沙水体随落潮流向口门以西水域扩散，沿程落淤，形成了埒子口～灌河口海域为淤泥质向粉沙质的过渡带。应该指出的是，落淤在水深较小浅滩上的粉沙质泥沙，易沉也易悬，一遇大风浪天气，再次被掀起悬浮，部分泥沙随着涨潮流又返回灌河，如此循环往复。

本发明利用所建立的泥沙数学模型，不仅可以计算悬沙的含沙量场，还可以计算悬沙的沉速场，并由此分析不同性质泥沙的分布特征，根据计算的悬沙沉速场可以看出：埒子口以西海域，悬沙沉速小于0.6mm/s，多以絮凝沉速沉降，属于淤泥质海岸；埒子口以东海域，悬沙沉速大于1mm/s，多以粉沙质泥沙沉速沉降，属于粉沙质海岸；至于埒子口海域，悬沙沉速在0.6mm/s～1mm/s之间，与该海域处于粉沙质海岸向淤泥质海岸的过渡特征相一致。

综上所述，对本实施例总结如下：

(1)连云港大范围海域底质比较复杂，既有粉沙质泥沙，又有淤泥质泥沙，其中，埒子口以西的底质沉积物以d₅₀＜0.031mm的细颗粒泥沙为主，而埒子口以东的底质沉积物以d₅₀介于0.031～0.075mm之间的粉沙质泥沙为主，粉沙质泥沙包络范围大致位于埒子口海域附近。

(2)建立了波流作用下的泥沙数学模型，可以用来模拟包含粉沙质和淤泥质在内的泥沙运动，特别是灌河口内陈家港的含沙量超过4.0kg/m³，而连云港海域洋山岛附近的含沙量仅0.2kg/m³左右，尽管两者相差很大，但所建模型均能给出较为满意的结果。

(3)粉沙质和淤泥质泥沙运动都以悬移质运动为主，连云港海域悬沙粒径小于0.063mm的细颗粒泥沙占90％以上，符合Burt沉速公式的应用条件，可以将其应用于连云港海域大范围悬沙沉速以及含沙量的模拟计算中。

(4)当只考虑潮流对泥沙起动的影响时，床面阻力对泥沙起动的影响最为主要，不仅影响着水流的计算，也与床面剪切应力的计算有关。床面阻力包括沙粒阻力和沙波阻力，沙粒阻力系床面沙粒的摩阻而引起，主要受粒径大小的影响，而沙波阻力属于形状阻力；由于粉沙质海岸床面形态较为复杂，床面阻力明显不同于淤泥质海岸，具体参数需通过率定得出。

本发明未进一步展开详细说明的内容均属于本技术领域的公知常识，以上实施例是本发明技术方案的的一个实例，本发明的公开范围并不局限于此，任何对于本发明的技术方案的等同替换，都落在本发明的保护范围之内。

Claims (3)

1.一种粉沙质和淤泥质海岸泥沙运动数值模拟方法，其特征在于，包括如下步骤：

步骤一：建立基本方程，所述基本方程包括水流连续方程、水流运动方程、悬移质输移扩散方程、床面冲淤变化方程；

步骤二，计算悬沙沉速、泥沙冲淤函数、波流共同作用下的床面剪切应力；

步骤三，泥沙数学模型的建立和求解：以MIKE21为工具，结合步骤一的基本方程、步骤二的计算结果建立泥沙数学模型，然后采用有限体积法求解泥沙数学模型，计算网格为三角网格，最大网格尺度1500m，最小20m；

步骤四，泥沙数学模型验证：

4-1，采用实测中潮资料对泥沙数学模型进行水流验证；

4-2，采用平均含沙量作为检验的标准，相对应的动力条件为典型中潮资料加上年均有效波高波浪场，对泥沙数学模型进行含沙量验证；

步骤五，根据泥沙数学模型判定悬沙沉速场：

5-1，根据悬沙沉速来判断泥沙性质：当沉速大于1mm/s的泥沙属于粉沙质泥沙，沉速在0.6mm/s以下的泥沙判定为淤泥质泥沙，而既包括粉沙质泥沙又包括淤泥质泥沙的粉沙淤泥质泥沙，其沉速范围为0.6mm/s～1mm/s；

5-2，根据模型图中潮周期内不同沉速线的外包络线判断悬沙沉速场与底质分布。

2.根据权利要求1所述的粉沙质和淤泥质海岸泥沙运动数值模拟方法，其特征在于，步骤一所述基本方程具体如下：

①水流连续方程：

式中：x、y为直角坐标系坐标；t为时间；h为平均水深；ζ为相对于平均海平面的潮位；U_x、U_y分别为x、y方向上的垂线平均速度；

②水流运动方程：

$\frac{\∂U_x}{\∂t}+U_x\frac{{\∂U}_x}{\∂x}+U_y\frac{\∂U_x}{\∂y}={\mathrm{fU}}_y-g\frac{\∂\mathrm{\ζ}}{\∂x}-\frac{{\mathrm{\τ}}_x}{\mathrm{\ρ}(h+\mathrm{\ζ})}+N_x(\frac{{\∂}^2{U}_x}{\∂x^2}+\frac{{\∂}^2{U}_x}{\∂y^2})$

$\frac{\∂U_y}{\∂t}+U_x\frac{{\∂U}_y}{\∂x}+U_y\frac{\∂U_y}{\∂y}={-\mathrm{fU}}_x-g\frac{\∂\mathrm{\ζ}}{\∂y}-\frac{{\mathrm{\τ}}_y}{\mathrm{\ρ}(h+\mathrm{\ζ})}+N_y(\frac{{\∂}^2{U}_y}{\∂x^2}+\frac{{\∂}^2{U}_y}{\∂y^2})---\left(2\right)$

式中：ρ为水体密度；g为重力加速度；N_x、N_y分别为x、y方向的水平紊动粘性系数；f为科氏参数；τ_x、τ_y分别为波流共同作用下床面剪切应力在x、y方向的分量；

③悬移质输移扩散方程：

式中：S为垂线平均含沙量；D_x、D_y分别为x、y方向的泥沙扩散系数；F_s为泥沙冲淤函数；

④床面冲淤变化方程：

${\mathrm{\γ}}_d\frac{{\∂\mathrm{\η}}_b}{\∂t}-F_s=0---\left(4\right)$

式中：γ_d为床沙干密度，η_b为海底床面的竖向位移，即冲淤变化量；F_s为泥沙冲淤函数。

3.根据权利要求1所述的粉沙质和淤泥质海岸泥沙运动数值模拟方法，其特征在于，步骤二的计算具体如下：

(1)悬沙沉速：

针对中值粒径小于63μm和含沙量小于10kg/m³的泥沙，关于沉速与含沙量之间的关系式：

ω＝kS^α                                (5)

式中k、a为经验系数，ω为泥沙沉降速度，S为含沙量；

(2)泥沙冲淤函数：

泥沙冲淤函数F_s与底部剪切应力及泥沙特征的关系由下式确定：

式中：τ为水流底部剪切应力；τ_d为不淤临界剪切应力；τ_e为起动临界剪切应力；α为淤积概率；M为冲刷系数；ω为泥沙沉降速度；

(3)波流共同作用下的床面剪切应力：

(3)-a，当只考虑水流作用时，床面平均剪切应力：

${\mathrm{\τ}}_c=\frac{1}{2}{\mathrm{\ρf}}_c{U}_c^2$

式中：U_c为水流平均速度；f_c为水流摩阻系数，可用下式计算

$f_c=2{\left\{2.5\right[\ln \left(\frac{30h}{k}\right)-1\left]\right\}}^{-2}---\left(7\right)$

其中k为床面阻力系数，h为平均水深；

(3)-b，当只考虑波浪作用时，床面平均剪切应力：

${\mathrm{\τ}}_w=\frac{1}{2}{\mathrm{\ρf}}_w{U}_b^2---\left(8\right)$

式中：U_b为波浪底部水质点水平运动速度；f_w为波浪摩阻系数，其中

$\begin{array}{cc}{f}_w=0.00251\exp \left[{5.21(A/k_s)}^{-0.19}\right]& A/k_s>1.57\\ =0.3& A/k_s\≤1.57\end{array}---\left(9\right)$

其中k_s为Nikuradse糙率，取30z₀；z₀为波浪摩阻高度；A为波浪近底质点运动的振幅，按照线性波理论近似表示为：

$A=\frac{H}{\mathrm{\π}}\frac{1}{\mathrm{sh}(2\mathrm{\πh}/L)}---\left(10\right)$

式中，H代表波浪高度，L代表波浪的长度；

(3)-c，当波流共同作用，采用波浪摩阻系数概念

$\frac{{\mathrm{\τ}}_{\mathrm{mean}}}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}(1+b{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^p{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^p)---\left(11\right)$

$\frac{{\mathrm{\τ}}_{\max }}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}=1+a{\left(\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w}\right)}^m{(1-\frac{{\mathrm{\τ}}_c}{{\mathrm{\τ}}_c+{\mathrm{\τ}}_w})}^n---\left(12\right)$

式中：τ_c为考虑纯水流作用下的切应力，τ_w为考虑纯波浪作用的切应力，b、p、q、a、m、n为综合表达式。

Images