CN105260571A - 一种滩浅海水动力模拟方法 - Google Patents

Abstract

本发明公开了一种滩浅海水动力模拟方法，包括以下步骤：（1）、获取目标海域的水深数据；（2）、网格划分步骤，将目标海域分别采用主网格和子网格两套网格进行划分；（3）、获取曼宁系数的分布；（4）、建立边界条件和初始条件；（5）、基于改进型浅水方程，在所述主网格上执行改进型浅水方程的离散和计算，在所述子网格上计算孔隙率Θ、等效水深Y和等效底摩阻系数Cd，最终计算输出目标海域的流速、水位结果。本滩浅海水动力模拟方法，同时采用两套网格，即较粗的主网格和较细的子网格。主网格计算效率高，而子网格的信息精度高，本发明结合了这两点优势，提高了计算速度的同时增加了计算精度。

Description

一种滩浅海水动力模拟方法

技术领域

本发明涉及一种涉及海岸动力模拟技术，尤其涉及一种滩浅海水动力模拟方法。

背景技术

滩浅海是指水深较浅的近岸滩涂区域。我国一些滩浅海区域不仅是重要的滩涂湿地自然保护区，同时也蕴藏着丰富的石油、天然气等资源，如黄河三角洲区域。在海岸带的开发过程中，需要对工程所在海域进行数值模拟，以掌握该海域的水动力规律，预测工程建设对水动力要素、生态环境、地貌形态等的影响，为海岸的合理开发、科学保护提供指导依据。

在滩浅海区域，潮沟多级发育，各级潮沟相互连接形成复杂的潮沟网系。潮沟网系是滩浅海区域滩涂、湿地部分与外海进行水体交换、泥沙输运的主要通道；潮沟-滩涂交错的地形特点，使得水深在很小的距离上就发生较大的变化。这些特点决定了滩浅海水动力的复杂性，也构成了对其进行数值模拟的主要困难：模型中需要以很高的精度分辨出精细的潮沟-滩涂地形。对于传统的海岸水动力模型，精度提高必然导致网格数增加、时间步长减小，计算量大大增加。根据作者实际经验，对于10平方千米的海域，采用传统浅水方程模型NearCoM，以2m精度模拟10天过程，采用800核并行计算时需要耗时3天。如此庞大的计算量无论是在时间上还是经济上都难以承受，因而限制了其在大范围、长时间、高精度的滩浅海水动力模拟中的应用。

在此背景下，如何发明一种适用于滩浅海动力环境的多尺度、高效率模拟技术，是本发明主要解决的技术问题。

发明内容

本发明为了解决现有模拟技术在应对滩浅海特殊的动力环境时，存在着计算量极其庞大、干-湿算法不足和缺乏精确的底摩阻处理方法等问题，提出了一种滩浅海水动力模拟方法，可以解决上述问题。

为了解决上述技术问题，本发明采用以下技术方案予以实现：

一种滩浅海水动力模拟方法，包括以下步骤：

(1)、获取目标海域的水深数据；

(2)、网格划分步骤，将目标海域分别采用主网格和子网格两套网格进行划分，所述主网格的宽度大于所述子网格的宽度，每一个主网格内包括若干子网格；

(3)、获取曼宁系数的分布，根据目标海域的海床泥沙特性、潮沟潮滩的分布以及植被的空间分布，建立曼宁系数的分布；

(4)、建立边界条件和初始条件；

(5)、基于改进型浅水方程，在所述主网格上执行改进型浅水方程的离散和计算，在所述子网格上计算孔隙率Θ、等效水深Y和等效底摩阻系数C_d，最终计算输出目标海域的流速、水位结果，所述改进型浅水方程为：

$\mathrm{\Θ}\frac{\∂\mathrm{\η}}{\∂t}+\frac{\∂P}{\∂x}+\frac{\∂Q}{\∂y}=0$

$\frac{\∂P}{\∂t}+\frac{\∂}{\∂x}\left({\∈}_{xx}\frac{P^2}{Y}\right)+\frac{\∂}{\∂y}\left({\∈}_{xy}\frac{PQ}{Y}\right)+gY\frac{\∂\mathrm{\η}}{\∂x}-\frac{\∂R_{xx}}{\∂x}-\frac{\∂R_{xy}}{\∂y}-\frac{{\mathrm{\τ}}_{sx}}{\mathrm{\ρ}}+\frac{{\mathrm{\τ}}_{bx}}{\mathrm{\ρ}}=0$

$\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left({\∈}_{xy}\frac{PQ}{Y}\right)+\frac{\∂}{\∂y}\left({\∈}_{yy}\frac{Q^2}{Y}\right)+gY\frac{\∂\mathrm{\η}}{\∂y}-\frac{\∂R_{xy}}{\∂x}-\frac{\∂R_{yy}}{\∂y}-\frac{{\mathrm{\τ}}_{sy}}{\mathrm{\ρ}}+\frac{{\mathrm{\τ}}_{by}}{\mathrm{\ρ}}=0$

式中η为主网格水位，P和Q分别为所述主网格x和y方向上的流量，Θ为孔隙率，所述孔隙率为所述主网格内湿的子网格个数占其所在主网格的总子网格个数的比例，Y为等效水深，也即所述主网格上所有湿的子网格水深的平均，C_d为等效底摩阻，也即所述主网格的底摩阻，R为水体黏性导致的剪切应力张量。

与现有技术相比，本发明的优点和积极效果是：本发明的滩浅海水动力模拟方法，不同于传统水动力模型中地形数据和计算网格采用同一套网格的做法，本发明同时采用两套网格，即较粗的主网格和较细的子网格。每个主网格中包含若干子网格，浅水方程的离散和计算在主网格上执行，而高精度的地形和底摩阻信息在子网格上给出，计算过程中，高精度子网格上的信息以孔隙率、等效水深、等效底摩阻的形式反馈到主网格上。主网格计算效率高，而子网格的信息精度高，本发明结合了这两点优势，提高了计算速度的同时增加了计算精度。

结合附图阅读本发明实施方式的详细描述后，本发明的其他特点和优点将变得更加清楚。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明所提出的滩浅海水动力模拟方法一种实施例中所划分的主网格和子网格的示意图；

图2是本发明所提出的滩浅海水动力模拟方法一种实施例中主网格上的变量MASK定义示意图；

图3是本发明所提出的滩浅海水动力模拟方法一种实施例中滩浅海区域图；

图4是图3中计算结果与实测数据的对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一，为了对本技术方案更清楚、完整的说明，首先介绍一下本技术方案中涉及的一些概念，如图1所示，为本实施例中所划分的主网格和子网格的示意图，其中，主网格的大小为Δx×Δy，每个主网格内的子网格大小为δx×δy，A2平面为平均海平面位置，A1平面为某一时刻的水面所在位置，η为水位，hs为定义在子网格水深，也即子网格的静止水深，Hs为该子网格的总水深。

本实施例提出了一种滩浅海水动力模拟方法，包括以下步骤：

S1、获取目标海域的水深数据；

滩浅海可以分为几个特征显著的部分，即陆地、潮间带、潮沟和外海。潮间带区域在大潮低潮期间大部分露出水面，因此可以采用机载激光雷达(LiDAR)进行扫描测量，获取露出区域和陆地区域的高程。大部分潮沟被水所覆盖，潮沟可以利用GPS定位-高程仪测量进行水深走航调查。对于外海区域，可以采用该地区的海图数据，将上述地形数据统一至该海域平均海平面的基准。

S2、网格划分步骤，将目标海域分别采用主网格和子网格两套网格进行划分，所述主网格的宽度大于所述子网格的宽度，每一个主网格内包括若干子网格；

根据步骤S1中所获取的水深数据覆盖区域，确定计算范围。根据该海区地形的特征(如主潮沟的宽度、进潮口的大小等)，结合所能使用的计算资源，选取合适的主网格大小。一般来说，关键地形特征应该有2-3个主网格来分辨，如主潮沟的进潮口处宽100米，则主网格大小最好能达到30-50米，根据计算范围和主网格的尺寸，则可以建立其覆盖目标海域的计算主网格。

每个主网格内包含若干个子网格，子网格的尺寸取决于步骤S1中所获取目标海域的水深数据的分辨率，目前的LiDAR测量精度一般可达1米左右，故可以将子网格尺寸定为1米。这个精度足以分辨出大部分滩浅海的盐沼湿地、滩涂的潮沟。将水深数据采用加权反距离方法插值到子网格上，获得精细子网格上的水深数据。

S3、获取曼宁系数的分布，根据目标海域的海床泥沙特性、潮沟潮滩的分布以及植被的空间分布，建立曼宁系数的分布；

根据目标海域的海床泥沙特性(淤泥质或砂质)、潮沟潮滩的分布以及植被的空间分布，建立曼宁系数的分布。曼宁系数在模型中将被用来计算底摩阻系数，属于可率定参数，计算时需要采用实测数据对该参数及其空间分布采用“trialanderror”方法进行率定。

S4、建立边界条件和初始条件；边界条件反应了模型范围以外的海洋动力对模型内部的驱动作用。

S5、基于改进型浅水方程，将步骤S1-S4中的数据输入至改进型浅水方程，在所述主网格上执行改进型浅水方程的离散和计算，在所述子网格上计算孔隙率Θ、等效水深Y和等效底摩阻系数C_d，最终计算输出目标海域的流速、水位结果，所述改进型浅水方程为：

$\mathrm{\Θ}\frac{\∂\mathrm{\η}}{\∂t}+\frac{\∂P}{\∂x}+\frac{\∂Q}{\∂y}=0$

$\frac{\∂P}{\∂t}+\frac{\∂}{\∂x}\left({\∈}_{xx}\frac{P^2}{Y}\right)+\frac{\∂}{\∂y}\left({\∈}_{xy}\frac{PQ}{Y}\right)+gY\frac{\∂\mathrm{\η}}{\∂x}-\frac{\∂R_{xx}}{\∂x}-\frac{\∂R_{xy}}{\∂y}-\frac{{\mathrm{\τ}}_{sx}}{\mathrm{\ρ}}+\frac{{\mathrm{\τ}}_{bx}}{\mathrm{\ρ}}=0$

$\frac{\∂Q}{\∂t}+\frac{\∂}{\∂x}\left({\∈}_{xy}\frac{PQ}{Y}\right)+\frac{\∂}{\∂y}\left({\∈}_{yy}\frac{Q^2}{Y}\right)+gY\frac{\∂\mathrm{\η}}{\∂y}-\frac{\∂R_{xy}}{\∂x}-\frac{\∂R_{yy}}{\∂y}-\frac{{\mathrm{\τ}}_{sy}}{\mathrm{\ρ}}+\frac{{\mathrm{\τ}}_{by}}{\mathrm{\ρ}}=0$

式中η为主网格水位，P和Q分别为所述主网格x和y方向上的流量，Θ为所述主网格内湿的子网格个数占其所在主网格的总子网格个数的比例，也即孔隙率，Y为等效水深，也即所述主网格上所有湿的子网格水深的平均，C_d为等效底摩阻，也即所述主网格的底摩阻，R为水体黏性导致的剪切应力张量。

本实施例的滩浅海水动力模拟方法，不同于传统水动力模型中地形数据和计算网格采用同一套网格的做法，本发明同时采用两套网格，即较粗的主网格和较细的子网格。每个主网格中包含若干子网格，浅水方程的离散和计算在主网格上执行，而高精度的地形和底摩阻信息在子网格上给出，计算过程中，高精度子网格上的信息以孔隙率、等效水深、等效底摩阻的形式反馈到主网格上。主网格计算效率高，而子网格的信息精度高，本发明结合了这两点优势，提高了计算速度的同时增加了计算精度。本实施例的滩浅海多尺度水动力模拟方法，假设主网格内水位一致、水流方向一致，且摩阻坡度在子网格上的变化可以忽略，考虑水深、底摩阻、流速在子网格尺度上空间变化的等效底摩阻系数。

本实施例中采用改进型浅水方程，与传统的浅水方程相比，改进方程中出现了系数Θ，该变量在物理意义上可将其看成孔隙率，以考虑主网格为部分干部分湿的状态，可定义为湿的子网格个数占所在主网格的总子网格点数的比例，

所述孔隙率Θ的计算方法如下：

其中，A为主网格的面积，hs为子网格水深，z为子网格实际水面位置，Hs为该子网格的总水深，对于空间任意一点x，定义其状态函数ψ(x)，如果该点为湿时其值为1，反之则为0，表达式：

所述等效水深Y的计算方法如下：

所述等效底摩阻系数Cd的计算方法如下：

$C_d=\frac{Y}{H_f}$

$H_f={\[\frac{\underset{A}{\∫\∫}{H}_s\sqrt{H_s/C_{ds}}dA}{AY}\]}^2.$

所述步骤S5中，离散后的改进型浅水方程为：

$P^{n+1}=-A^n\frac{\∂\mathrm{\η}}{\∂x}+B^n$

$Q^{n+1}=-A^n\frac{\∂\mathrm{\η}}{\∂x}+C^n$

其中：

$A^n={\left(\frac{{\mathrm{gY}}^2\mathrm{\Δ}t}{Y+C_d\left|u\right|\mathrm{\Δ}t}\right)}^n$

$B^n=\frac{Y\mathrm{\Δ}t}{Y+C_d\left|u\right|\mathrm{\Δ}t}{\[-\frac{\∂}{\∂x}\left(\frac{P^2}{Y}\right)-\frac{\∂}{\∂y}\left(\frac{PQ}{Y}\right)+\frac{\∂R_{xx}}{\∂x}+\frac{\∂R_{xy}}{\∂y}+{\mathrm{\τ}}_{sx}+\frac{P}{\mathrm{\Δ}t}\]}^n$

$C^n=\frac{Y\mathrm{\Δ}t}{Y+C_d\left|u\right|\mathrm{\Δ}t}{\[-\frac{\∂}{\∂x}\left(\frac{PQ}{Y}\right)-\frac{\∂}{\∂y}\left(\frac{Q^2}{Y}\right)+\frac{\∂R_{xy}}{\∂x}+\frac{\∂R_{yy}}{{\∂}_y}+{\mathrm{\τ}}_{sy}+\frac{Q}{\mathrm{\Δ}t}\]}^n$

将该式代入到连续性方程，得到一组混合微分-差分方程：

$\frac{\∂\mathrm{\η}}{\∂t}=\frac{1}{\mathrm{\Θ}}\[\frac{\∂}{\∂x}\left(A\frac{\∂\mathrm{\η}}{\∂x}\right)+\frac{\∂}{\∂y}\left(A\frac{\∂\mathrm{\η}}{\∂y}\right)-\frac{\∂B}{\∂x}-\frac{\∂C}{\∂y}\].$

混合微分-差分方程优选采用奇偶跳点法分两步进行求解，可以求解出主网格水位η，进而可以通过改进型浅水方程求得P和Q，也即获得模拟滩浅海水动力的主要参数，实现对滩浅海水动力模拟。

步骤S4中，边界条件反应了模型范围以外的海洋动力对模型内部的驱动作用。对于不同的资料收集程度，本步骤给出不同的边界条件获取方法。对于平面二维水动力模型，边界条件仅需给出水位一个变量即可满足数值模型的适定性。当计算范围较大、边界位于外海处时，很难保证边界上每个网格节点都有实测数据，因此可以采用所在海域已有的其他大区域海洋模型进行嵌套计算，即大模型在小模型的边界处输出水位时间序列，并采用空间上线性插值、时间上二次插值的方法获得边界上各个节点的水位边界条件。当计算范围较小、边界附近设有潮位观测站时，可以根据实测的水位时间序列，采用空间上线性插值、时间上二次插值的方法，获得模型边界每个网格节点上的水位边界条件。对于初始水位场和流速场，由于初始条件对水动力的影响在几个潮周期后即可消去，故可设为零，即冷启动。

所述步骤S5中，还包括定义主网格上的变量MASK的步骤，其位置位于主网格的中心，MASK值与所述主网格的等效水深有关，如图2所示，其判断方法如下：

MASK_i，j＝1ifY_i，j＞0

MASK_i，j＝0ifY_i，j≤0

其中，MASK＝0代表该主网格为干，MASK＝1代表该主网格为湿，i,j分别为主网格在x和y方向上的网格编号；

参考MASK的值，给出干湿条件，当同时满足如下两条件时，关闭x向的流量，即令P等于0：

a、该x向的流量所通过的网格面两侧形成干湿界限，也即MASK_i,j＝1且MASK_i+1,j＝0；

b、该x向的流量所通过的网格面两侧的主网格至少一侧的等效水深小于临界值；

当同时满足如下两条件时，关闭y向的流量，即令Q等于0：

a、该y向的流量所通过的网格面两侧形成干湿界限，也即MASK_i,j＝1且MASK_i,j+1＝0；

b、该y向的流量所通过的网格面两侧的主网格至少一侧的等效水深小于临界值。

孔隙率Θ、等效水深Y和等效底摩阻系数Cd的计算均需要在子网格尺度上进行积分得到，且每一个时间步都需要更新，因此仍然需要耗费较多的计算时间。在这几个变量的计算中，只有主网格水位η是未知的。因此对于预定的一系列离散水位值，我们可以预先计算出相应的孔隙率Θ、等效水深Y和等效底摩阻系数Cd的值；将这些值存储在内存中，在模型运行的每个时间步上，采用插值或者多项式拟合的方法可以获得连续的变量值，从而避免了每步都进行繁琐的积分计算。总结出来即为“预先计算——储存查找表——插值或多项式拟合”，下面详细给出预先计算并存储孔隙率Θ、等效水深Y、以及等效底摩阻系数Cd的步骤的具体方法。包括：

首先，设定水位值η’的变化步长Δη，得到一组离散水位值，η_min<η’<η_max，定义离散水位值：

η_k＝η_min+(k-1)Δη；k＝1,…,K，其中Δη＝(η_max-η_min)/(K-1)；

其次，预先计算出上述离散水位值η_k下的孔隙率Θ_k、等效水深Y_k和等效底摩阻系数C_d，k，并存储在内存中；当采用多项式拟合时，则存储多项式的系数。

最后，在每一个时间步上，采用插值或者多项式拟合的方法获取连续的变量值Θ(η)，Y(η)和Cd(η)。

在采用插值或者多项式拟合的方法获取连续的变量值步骤中，以孔隙率的多项式拟合方法为例，其计算公式为：

$\mathrm{\&Theta;}\left(\mathrm{\&eta;}\right)=\left\{\begin{array}{cc}0.0& {\mathrm{\&eta;}}_{\min }\&le;\mathrm{\&eta;}<z_{low}\\ \underset{i=0}{\overset{n}{\&Sigma;}}{a}_i{\mathrm{\&eta;}}^i& z_{low}\&le;\mathrm{\&eta;}\&le;z_{high}\\ 1.0& z_{high}<\mathrm{\&eta;}\&le;{\mathrm{\&eta;}}_{\max }\end{array}\right.$

其中z_low是该主网格中底高程的最小值，z_high为该主网格中底高程的最大值。

为了验证本模拟方法的精确性，还包括模型率定和验证的步骤。模型率定的主要参数为曼宁系数，采用方法一般为试算法，将步骤S5中计算的结果与实际测量数据的对比，根如果模拟的流速偏大、水位偏高，则可以调高曼宁系数；反之则可以调低曼宁系数。一般来说，当模拟水位与实测水位的相关系数达到0.9，流速的相关系数达到0.7，即可认为模拟结果可以较好的反映真实的水动力环境。

模拟结果处理和分析。模型计算可以获得整个区域的水位场、流速场、底部剪切应力等水动力要素。采用Matlab对计算结果进行处理、绘图，分析潮汐水位的涨落规律，以及涨急、落急、涨憩、落憩流场的分布，根据底部剪切应力的分布，可以预测容易发生侵蚀或泥沙淤积的地区。这些结果可以为滩浅海的工程建设、环境生态预测提供有意义的依据。

为了验证本实施例的滩浅海水动力模拟方法，在实际滩浅海区域的模拟能力，特选取美国特拉华州的一处滩浅海区域，包括滩涂、盐沼和近岸浅水地区，其地形如图3所示。为了展示本发明提出的基于子网格概念的模拟技术相对于传统模拟技术在计算效率方面的提高，我们设计了5组模拟，分别称之为A，B，B1，B2，C组，见表1。

表1

A组直接采用2m×2m高精度的地形数据进行计算，因此不需要开启子网格技术，我们称A组为全网格模型组；B组的主网格精度设置为8m×8m，开启子网格技术，每个主网格内含有4×4个子网格。C组中原始的2m×2m精度地形被平滑到8m×8m的计算网格上，但不开启子网格技术。C组在这里主要作为对比组，用以说明采用8m的精度不能准确模拟该地区的水动力特点。在B组中，关闭预存储技术，每个时间上变量孔隙率Θ、等效水深Y和等效底摩阻系数Cd均需要在子网格上进行积分计算；B1和B2组中则开启预存储计算，分别采用插值方法和多项式拟合方法。

模型边界条件采用实测水位时间序列进行驱动，对于滩涂盐沼区域的底摩阻曼宁系数取为0.25，对于潮沟和近岸区域取为0.02。模型运行10天，图3中7个站位的实测数据用来验证本模拟方法的有效性。计算结果与实测数据的对比在图4中给出，其中虚线为实测数据，点画线为A组全网格模型结果，点线为B组开启子网格技术的结果，实线为C组粗网格直接计算的结果。B1和B2组也开启了子网格技术，不同的是同时采用了预存储技术，其结果与B组几乎完全一致，所以在图中不予展示。

可以看出，由于采用了本发明提出的计算方法，B组采用8m×8m的计算精度，其结果与A组2m×2m计算精度的结果几乎完全一样，且与实测数据吻合非常好；作为对比，C组虽然也是8m×8m的计算精度，但没有采用子网格技术，其模拟结果差了很多，尤其是在离外海更远的细窄潮沟的站位处。

表1中最后一栏给出了各组模拟的计算耗时，可以看出同样都采用20核并行计算的前提下，A组模拟耗时763.3小时，而B组只需要20.7小时，计算效率提高了36倍；当在B组基础上采用预存储技术时(B1、B2组)，时间进一步缩短到7.0小时，计算效率提高了109倍；而B、B1、B2的模拟结果却是跟A组几乎完全一样的。综上可知，本发明给出的模拟技术在维持了计算准确性的同时，计算效率可以大幅提升。本技术兼顾高精度与高效性，是大范围、高精度、长时间的复杂滩浅海地区模拟的一个新选择。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的普通技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (7)

1.一种滩浅海水动力模拟方法，其特征在于，包括以下步骤：

(1)、获取目标海域的水深数据；

(2)、网格划分步骤，将目标海域分别采用主网格和子网格两套网格进行划分，所述主网格的宽度大于所述子网格的宽度，每一个主网格内包括若干子网格；

(3)、获取曼宁系数的分布，根据目标海域的海床泥沙特性、潮沟潮滩的分布以及植被的空间分布，建立曼宁系数的分布；

(4)、建立边界条件和初始条件；

(5)、基于改进型浅水方程，在所述主网格上执行改进型浅水方程的离散和计算，在所述子网格上计算孔隙率Θ、等效水深Y和等效底摩阻系数C_d，最终计算输出目标海域的流速、水位结果，所述改进型浅水方程为：

$\mathrm{\&Theta;}\frac{\&part;\mathrm{\&eta;}}{\&part;t}+\frac{\&part;P}{\&part;x}+\frac{\&part;Q}{\&part;y}=0$

$\frac{\&part;P}{\&part;t}+\frac{\&part;}{\&part;x}\left({\&Element;}_{xx}\frac{P^2}{Y}\right)+\frac{\&part;}{\&part;y}\left({\&Element;}_{xy}\frac{PQ}{Y}\right)+gY\frac{\&part;\mathrm{\&eta;}}{\&part;x}-\frac{\&part;R_{xx}}{\&part;x}-\frac{\&part;R_{xy}}{\&part;y}-\frac{{\mathrm{\&tau;}}_{sx}}{\mathrm{\&rho;}}+\frac{{\mathrm{\&tau;}}_{bx}}{\mathrm{\&rho;}}=0$

$\frac{\&part;Q}{\&part;t}+\frac{\&part;}{\&part;x}\left({\&Element;}_{xy}\frac{PQ}{Y}\right)+\frac{\&part;}{\&part;y}\left({\&Element;}_{yy}\frac{Q^2}{Y}\right)+gY\frac{\&part;\mathrm{\&eta;}}{\&part;y}-\frac{\&part;R_{yx}}{\&part;x}-\frac{\&part;R_{yy}}{\&part;y}-\frac{{\mathrm{\&tau;}}_{sy}}{\mathrm{\&rho;}}+\frac{{\mathrm{\&tau;}}_{by}}{\mathrm{\&rho;}}=0$

式中η为主网格水位，P和Q分别为所述主网格x和y方向上的流量，Θ为孔隙率，所述空隙率为所述主网格内湿的子网格个数占其所在主网格的总子网格个数的比例，Y为等效水深，也即所述主网格上所有湿的子网格水深的平均，C_d为等效底摩阻，也即所述主网格的底摩阻，R为水体黏性导致的剪切应力张量。

2.根据权利要求1所述的滩浅海水动力模拟方法，其特征在于，

所述步骤(5)中，离散后的改进型浅水方程为：

$P^{n+1}=-A^n\frac{\&part;\mathrm{\&eta;}}{\&part;x}+B^n$

$Q^{n+1}=-A^n\frac{\&part;\mathrm{\&eta;}}{\&part;x}+C^n$

其中：

$A^n={\left(\frac{{\mathrm{gY}}^2\mathrm{\&Delta;}t}{Y+C_d\left|u\right|\mathrm{\&Delta;}t}\right)}^n$

$B^n=\frac{Y\mathrm{\&Delta;}t}{Y+C_d\left|u\right|\mathrm{\&Delta;}t}{\&lsqb;-\frac{\&part;}{\&part;x}\left(\frac{P^2}{Y}\right)-\frac{\&part;}{\&part;y}\left(\frac{PQ}{Y}\right)+\frac{\&part;R_{xx}}{\&part;x}+\frac{\&part;R_{xy}}{\&part;y}+{\mathrm{\&tau;}}_{sx}+\frac{P}{\mathrm{\&Delta;}t}\&rsqb;}^n$

$C^n=\frac{Y\mathrm{\&Delta;}t}{Y+C_d\left|u\right|\mathrm{\&Delta;}t}{\&lsqb;-\frac{\&part;}{\&part;x}\left(\frac{PQ}{Y}\right)-\frac{\&part;}{\&part;y}\left(\frac{Q^2}{Y}\right)+\frac{\&part;R_{xy}}{\&part;x}+\frac{\&part;R_{yy}}{\&part;y}+{\mathrm{\&tau;}}_{sy}+\frac{Q}{\mathrm{\&Delta;}t}\&rsqb;}^n$

$\frac{\&part;\mathrm{\&eta;}}{\&part;t}=\frac{1}{\mathrm{\&Theta;}}\&lsqb;\frac{\&part;}{\&part;x}\left(A\frac{\&part;\mathrm{\&eta;}}{\&part;x}\right)+\frac{\&part;}{\&part;y}\left(A\frac{\&part;\mathrm{\&eta;}}{\&part;y}\right)-\frac{\&part;B}{\&part;x}-\frac{\&part;C}{\&part;y}\&rsqb;.$

3.根据权利要求1所述的滩浅海水动力模拟方法，其特征在于，所述孔隙率Θ的计算方法如下：

其中，A为主网格的面积，hs为子网格水深，z为子网格实际水面位置，Hs为该子网格的总水深，对于空间任意一点x，定义其状态函数ψ(x)，如果该点为湿时其值为1，反之则为0，表达式：

4.根据权利要求3所述的滩浅海水动力模拟方法，其特征在于，所述等效水深Y的计算方法如下：

5.根据权利要求4所述的滩浅海水动力模拟方法，其特征在于，所述等效底摩阻系数Cd的计算方法如下：

$C_d=\frac{Y}{H_f}$

$H_f={\&lsqb;\frac{\underset{A}{\&Integral;\&Integral;}{H}_s\sqrt{H_s/C_{ds}}dA}{AY}\&rsqb;}^2.$

6.根据权利要求5所述的滩浅海水动力模拟方法，其特征在于，所述步骤(5)中，还包括定义主网格上的变量MASK的步骤，其位置位于主网格的中心，MASK值与所述主网格的等效水深有关，其判断方法如下：

MASK_i，j＝1ifY_i，j＞0

MASK_i，j＝0ifY_i，j≤0

其中，MASK＝0代表该主网格为干，MASK＝1代表该主网格为湿，i,j分别为主网格在x和y方向上的网格编号；

参考MASK的值，给出干湿条件，当同时满足如下两条件时，关闭x向的流量，即令P等于0：

a、该x向的流量所通过的网格面两侧形成干湿界限，也即MASK_i,j＝1且MASK_i+1,j＝0；

b、该x向的流量所通过的网格面两侧的主网格至少一侧的等效水深小于临界值；

当同时满足如下两条件时，关闭y向的流量，即令Q等于0：

a、该y向的流量所通过的网格面两侧形成干湿界限，也即MASK_i,j＝1且MASK_i,j+1＝0；

b、该y向的流量所通过的网格面两侧的主网格至少一侧的等效水深小于临界值。

7.根据权利要求5所述的滩浅海水动力模拟方法，其特征在于，所述孔隙率Θ、等效水深Y、以及等效底摩阻系数Cd的计算方法中，还包括预先计算并存储孔隙率Θ、等效水深Y、以及等效底摩阻系数Cd的步骤，具体包括：

首先，设定水位值η’的变化步长Δη，得到一组离散水位值，η_min<η’<η_max，定义离散水位值：

η_k＝η_min+(k-1)Δη；k＝1,…,K，其中Δη＝(η_max-η_min)/(K-1)；

其次，预先计算出上述离散水位值η_k下的孔隙率Θ_k、等效水深Y_k和等效底摩阻系数C_d，k，并存储在内存中；

最后，在每一个时间步上，采用插值或者多项式拟合的方法获取连续的变量值Θ(η)，Y(η)和Cd(η)。