一种高海况下宽级配沙坝潟湖海域动力地貌数值模拟方法

1.本发明涉及海域动力地貌数值模拟领域，具体涉及一种高海况下宽级配沙坝潟湖海域动力地貌数值模拟方法。


背景技术：

2.在海岸滩涂动力地貌的模拟上，往往将不断变化的潮汐简化为有代表性的潮汐分潮，其研究大多注重较长时间内对地貌演变有净效果的泥沙余输运，很难考虑到长期演变过程中风、浪等驱动力的随机性，对于极端气候或动力条件引起的短期内的冲淤演变则难以把握。在一般的动力条件下海床易形成沙波，海岸泥沙输沙量并不明显。但在强风浪作用尤其是风暴潮来临的情况下，波浪作用会使泥沙大量悬扬，可能会导致短期内较大的冲淤变化。在沿海港口航道建设中，大量的实测资料和研究实践表明，一次大风或风暴潮往往造成港口航道大量淤积，带来巨大损失。
3.高海况(风暴潮或者大风天)下海岸动力地貌模拟准确的重要前提就是能够正确求解波流共同作用下的挟沙能力。目前对于波流挟沙能力求解方法主要包括经验分析法、因次分析法、能量平衡法和基于紊动猝发理论等。经验分析法通过对水流挟沙能力的各要素进行分析，寻找主要影响因素，采用实测资料进行参数拟合，从而建立适合特定海岸区域的经验公式。不同地区经验系数不同，因此经验公式具有简便和针对性强的特点，但受资料和海岸区域限制，适用性不广泛。因次分析方法是对描述各物理过程内在联系进行求解从而建立关系式，该类关系式具有结构合理以及物理意义清晰的无量纲关系，最具代表性的为刘家驹、罗肇森公式和胡春宏公式等。能量平衡法具有一定的理论基础，得出的挟沙能力公式较为可靠，因此广泛运用于波流泥沙运动研究。国内外学者对水流及波浪的挟沙能力进行了大量的研究。窦国仁等、王尚毅等、乐陪九、曹祖德等分别得到了波流共同作用下的挟沙能力。其中代表性的为窦国仁依据能量迭加原理，将潮流和波浪用于悬浮泥沙的能量相加，从而导出潮流和波浪共同作用下的挟沙能力公式。从机理上看，能量平衡方法在一定程度上也回避了潮流和波浪作用下泥沙运动的微观机制，且存在挟沙能力与周期成反比的不合理现象。
4.目前已有应用的波流挟沙能力公式大多数是半经验半理论公式，没有充分考虑泥沙运动的微观机制，还存在较大的弊端。特别是在高海况条件下，河床底部存在高含沙层运动，沙质泥沙或者淤泥质泥沙的层移运动下含沙量剖面均有一定的计算方法，但对于宽级配泥沙海域的波浪挟沙能力计算方法还不明确。


技术实现要素：

5.本发明的目的在于克服上述现有技术存在的不足，提供一种高海况下宽级配沙坝潟湖海域动力地貌模拟方法，该方法能够较为准确的模拟风暴潮或者极端天气下宽级配海域的泥沙骤淤，以期为保障风暴潮来袭时港口、航道、取排水口等水工建筑物的正常运行提供理论依据。
6.为实现上述技术目的，本发明采用如下方案：
7.一种高海况下宽级配沙坝潟湖海域动力地貌数值模拟方法，包括：
8.s1利用层移条件下悬沙层含沙量剖面的计算公式表征波浪作用下悬沙层含沙量分布；
9.s2建立高海况下波流挟沙能力的计算公式：
10.波流挟沙能力采用水流挟沙能力和波浪挟沙能力之和表征，其中：水流挟沙能力通过实测数据和水力因子拟合获取；波浪挟沙能力根据床面形态划分：
11.层移床面下，将所述含沙量剖面的计算公式根据水深积分，用以表征波浪挟沙能力；
12.沙波床面下采用沙波床面波浪挟沙能力公式计算波浪挟沙能力；
13.s3建立沙坝潟湖型海岸大范围潮波模型及局部动力地貌数学模型，并引入波浪模块建立耦合模型，利用s2获取的计算公式表征波流挟沙能力参数，模拟高海况下宽级配沙坝潟湖海域动力地貌。
14.作为一种优选的实施方式，所述s2中，基于波浪移动数判断层移临界条件，确定床面形态；
15.当ψ≥226.21β
c-0.89
red
*-0.274
时，利用根据水深积分的含沙量剖面表征波浪挟沙能力；ψ＜226.21β
c-0.89
red
*-0.274
时，采用沙波床面波浪挟沙能力公式计算波浪挟沙能力；
16.其中，ψ为波浪移动数，βc为密实度系数，red
*
为泥沙雷诺数。
17.作为一种优选的实施方式，对所述含沙量剖面的计算公式在层移层顶部高度至水面的水深范围内积分。
18.作为一种优选的实施方式，层移床面下波浪挟沙能力计算公式如下：
[0019][0020]
式中，
[0021][0022][0023][0024]
h为水深；α
fi
为考虑层化效应的床面平坦时的平均rouse系数，i＝1，2，3；和由每层的平均含沙量分别通过层化效应和阻滞沉速公式计算得到；zi为高度，c
zi
为zi高度处含沙量；za为层移层顶部高度；ca为层移层顶部含沙量；δw为波浪边界层厚度；rs＝0.8；af＝1.17。
[0025]
作为一种优选的实施方式，所述s3中，所述大范围潮波模型以基于浅水假设及boussinesq假设的不可压缩流体navier-stokes方程作为水动力控制方程。
[0026]
作为一种优选的实施方式，大范围潮波模型的开边界设置于断面平均水深＞2m处，且开边界与研究区域距离＞10km。
[0027]
作为一种优选的实施方式，所述s3中，局部动力地貌数学模型模拟时，将泥沙组分按照中值粒径范围a～b划分为三组粘性沙和三组非粘性沙，其中三组粘性沙的中值粒径分别为a，60，三组非粘性沙的中值粒径分别为100，b；单位均为μm。
[0028]
作为一种优选的实施方式，所述s3中，局部动力地貌数学模型模拟时，各泥沙组分的厚度基于如下方式确定：
[0029]
将海床实测点的中值粒径在海域上根据模型网格进行插值，得到海域内海床上任意一点的中值粒径；
[0030]
根据各点中值粒径和泥沙组分中值粒径数值范围，确定代表泥沙组分并获取泥沙组分厚度：
[0031][0032]
d2＝δ-d1[0033]
其中，d1为组分1厚度，c1为组分1中值粒径，d2为组分2厚度，c2为组分2中值粒径，c为海床上任意一点的中值粒径，c1＜c＜c2，δ为床面可动层厚度。
[0034]
作为一种优选的实施方式，局部动力地貌数学模型计算时，基于沙坝潟湖型海岸近n年间最大一次风暴潮发生后持续时间内最大冲刷深度确定床面可动层厚度，n＞40。
[0035]
作为一种优选的实施方式，所述s3中，局部动力地貌数学模型包含一个离岸边界和两个与海岸垂直的接岸边界；所述离岸边界采用水位边界，所述接岸边界采用neumann边界。
[0036]
本发明具有如下有益效果：
[0037]
1、建立高海况下宽级配海域波流挟沙能力公式，可以同时适用于沙波及层移质运动，并且同时适用于粉沙-沙的宽级配海域，采用本发明提出的挟沙能力公式可以更准确地模拟高海况下宽级配海域的动力地貌过程。
[0038]
2、本发明对高海况下沙坝潟湖海域初始海床地形进行精细化设置，根据不同泥沙粒径大小给出相适应的不同级配的海床初始厚度，使得初始海床泥沙更加符合实际情况，提升了高海况下宽级配海域的数值模拟计算精度。
[0039]
3、给出了高海况宽级配沙坝潟湖海域动力地貌过程预测的验证方法，并提出沙坝外侧的风暴剖面规律，从而根据验证结果可以直观的判断模型准确度。
附图说明
[0040]
图1为se向偏东6度10年一遇波浪有效波高分布图(中潮位：0.06m)。
[0041]
图2为10年一遇波浪下落急时刻含沙量场。
[0042]
图3为10年一遇波浪下老龙沟海区冲淤变化。
[0043]
图4为大风天前后沙坝外侧床面剖面高程对比所取剖面示意图。
[0044]
图5为大风天前后沙坝外侧床面剖面高程对比。
具体实施方式
[0045]
下面结合说明书附图和具体的实施例，对本发明作详细描述。
[0046]
实施例1
[0047]
(1)建立高海况下宽级配泥沙悬沙层含沙量剖面计算公式。
[0048]
波浪作用下悬沙层含沙量分布可采用层移条件下悬沙层含沙量剖面的计算公式，适用于细粉沙到粗沙较为广泛的粒径：
[0049][0050]
式中，αf为考虑层化效应的床面平坦时的rouse系数，其中u
′
*
为有效平均波浪剪切速度，u
′
*
＝0.5u
*w
，u
*w
为波浪最大剪切速度，φd为阻尼系数，κ为卡门系数，ws为阻滞泥沙沉速；δw为波浪边界层厚度，δw＝κu
*
/ω，ω为波浪频率，u
*
为波浪摩阻流速；rs＝0.8；af＝1.17；z1＝max(za，0.5δw)；z2＝max(za，2.5δw)；ca为层移层顶部含沙量，c
z1
为z1高度处含沙量，c
z2
为z2高度处含沙量，h
′
＝az2/b，a和b为拟合系数。
[0051]
鉴于底部高含沙层影响，公式中考虑了阻滞沉速效应和制紊效应。具体的，当含沙量较大时，沉速会减小，即所谓的阻滞沉速效应，对于细沙和粉沙制约沉速机理有所不同，这是因为d＜100μm时颗粒周围的流动为层流，而d＞100μm的沙颗粒周围的流动为紊流。细沙和粉沙处于stokes区，几何形状的影响不大。粉沙和细沙的阻滞泥沙沉速可采用：
[0052][0053]
式中：w
s，0
为泥沙初始沉速；cv为体积含沙量；φ
s，struct
＝0.5为结构密度，即固体含量所能达到的床面结构密度；φ
s，max
＝0.65为最大密度，即固体含量所能达到的最大体积含量；φ为泥沙密度；m＝1-2为系数，代表颗粒尾流的影响。
[0054]
而在浑水中，泥沙浓度梯度将导致密度变化。若含沙量梯度较大，其制紊效应可能改变流速剖面和输沙率，对细颗粒泥沙更为显著。因此，含沙量层化效应在高含沙量层模拟中是一个重要参数。
[0055]
为了考虑层化效应对沙波漩涡层涡黏系数的影响，基于下式计算阻尼系数：
[0056]
φd＝φ
fs
[1+(cv/c
gel，s
)
0.8-2(cv/c
gel，s
)
0.4
]；
[0057]
其中φ
fs
＝d
50
/(1.5d
sand
)；当d
50
≥1.5d
sand
时，φ
fs
＝1；c
gel，s
＝0.67为最大体积含沙量d
sand
＝0.12mm。
[0058]
(2)建立高海况下宽级配海域波流挟沙能力公式
[0059]
二维悬移质不平衡输沙基本方程为：
[0060]
[0061]
式中，h为水深，c为悬沙浓度，ci为不同泥沙组分的悬沙浓度；u
x
和uy分别是x和y方向上垂线流速分量；ε
s，x
和ε
s，y
是泥沙成分的涡动扩散系数，fs为源汇项。
[0062]
源汇项fs是泥沙模拟的关键问题，利用挟沙能力的概念，将悬浮和沉降通量表示为，
[0063]
e＝ωss
*
，d＝ωss
ꢀꢀꢀ
(4)
[0064]
考虑恢复饱和系数αs，源汇项可写成：
[0065]fs
＝asωs(s
*-s)
ꢀꢀꢀ
(5)
[0066]
式中e为悬浮通量，d为沉降通量，s
*
为挟沙能力，s为垂线平均含沙量，ωs为泥沙沉速。
[0067]
底床变形采用公式(6)计算：
[0068][0069]
式中γs为泥沙容重。
[0070]
底床冲淤采用含沙量与挟沙能力对比的判别条件。即当s＞s
*
，含沙量大于挟沙能力，底床淤积；当s≤s
*
，且u≥u0或τ≥τ0，含沙量小于挟沙能力，且流速(剪切力)大于起动流速(剪切力)，底床冲刷。
[0071]
底床更新通过求解底床变形方程得到底床冲淤厚度，反馈给地形，再进行水流、波浪及泥沙模拟。在一次风暴潮计算过程中，短时间内可能造成较大的冲淤变化，而计算时长有限，因此，从底床影响和计算效率两方面考虑，采用实时更新，逐步计算。
[0072]
从能量观点来看，波流挟沙能力表现为水流和波浪挟沙能力之和。鉴于水流挟沙能力研究相对成熟，在河口海岸一般为水流的二次方关系，因此采用广泛应用的水流挟沙能力关系与上述所得的波浪挟沙能力叠加，从而求解波流挟沙能力。
[0073]s*cw
＝s
*c
+s
*w
ꢀꢀꢀ
(7)
[0074]
式中，s
*cw
为波流挟沙能力，s
*w
为波浪挟沙能力，波浪因子考虑水流的影响，s
*c
为水流挟沙能力，通过实测资料和水力因子拟合，采用式(8)计算：
[0075][0076]
式中k0为系数，s0为背景含沙量。
[0077]s*w
为波浪挟沙能力，本技术对其进行改进，在高海况下海域底部床面发生层移质运动，在层移条件下，将含沙量剖面根据水深进行积分，可以得到垂线平均含沙量s
*
，即代表波浪挟沙能力，h为水深。
[0078]
对含沙量剖面计算表达式(式(1))进行积分可得：
[0079][0080]
式中，cz为z高度处含沙量。
[0081]
因此层移条件下s
*
可表示为：
[0082]
[0083]
式中，
[0084][0085][0086]
α
fi
为考虑层化效应的床面平坦时的平均rouse系数，i＝1，2，3；和由每层的平均含沙量分别通过层化效应和阻滞沉速公式计算而得。
[0087]
一次风暴潮过程包括发展、强波和衰退的过程，在波高较小时底床可能出现沙波床面，沙波床面波浪挟沙能力采用richardson j f(1997)发表的公式计算，如下：
[0088][0089]
式中，
[0090][0091][0092][0093]
α
ri
是考虑层化效应的沙波床面下平均rouse系数，i＝0，1，2。
[0094]dh
＝4(h-2η)-3(z
3r-2η)，每层的βi由下式计算，
[0095][0096]
式中，δz＝z
i+1-zi。
[0097]
加入密实度系数后的层移临界条件下移动数计算公式为(参见cn114611425a)：
[0098]
ψ＝226.21β
c-0.89
red
*-0.274
(13)
[0099]
其中，ψ为波浪移动数，βc为密实度系数，red
*
为泥沙雷诺数，βc＝(ρ0/ρ
0*
)
2.5
其中，ρ0和ρ
0*
分别是泥沙干容重和稳定干容重，ρ
0*
＝0.68ρs(d
50
/d0)n，n＝0.08+0.014(d
50
/d
25
)为系数；d0为0.001m；d
25
为体积百分比为25％的较细的泥沙所对应的粒径大小。
[0100]
根据层移临界条件计算公式(13)可判断床面处于何种形态。
[0101]
则由式(10)和式(11)可得到沙波及层移条件下的垂线平均含沙量计算公式：
[0102]
[0103]
(3)建立模拟高海况下宽级配沙坝潟湖海域动力地貌的耦合模型；
[0104]
①
建立渤海湾大范围海域潮波模型，开边界设置在大连至烟台，该区域水深较大且远离研究区域。在开边界上，水位由8个天文潮汐分潮组成(m2，s2，n2，k2，k1，o1，p1，q1)，这些分潮的振幅和相位是从东中国海潮波预报模型中获取，该模型范围包含渤海、黄海和东海且经过较好验证。
[0105]
根据以下公式计算开放边界处的水位：
[0106][0107]
其中，η(t)表示关于时间t的水位；a0表示一个潮周期的平均水位；n表示相关分潮的个数；表示每个分潮的振幅；fi表示节点振幅系数；ωi表示频率；表示相位。
[0108]
对于研究区域内的浅水区(如潮滩)，该模型考虑了干湿过程，其主要方法是在每个时间步计算之前对水位点和速度点(由于使用正交网格)进行一系列检查。对于水位点，使用周围四个深度点的最大值(d
max
)作为水深。在速度点，床面高程是单元顶点(d
mean
)中指定值的算术平均值。对于速度点处的水位，采用所谓的逆风方法。通过比较水位点和速度点的总水深与阈值深度
±
(本文中
±
＝0.01m)来确定计算网格是湿还是干。该模型在零初始条件(冷启动)下启动，运行3个月(2006年2月1日至4月1日)，涵盖渤海湾实测数据过程。前7天被视为螺旋上升期，在分析中被忽略。
[0109]
②
采用neumann开边界建立沙坝潟湖型潮汐汊道海域中长时间尺度动力地貌数学模型；
[0110]
模型中的粘性沙计算通过用源和汇项来解决平流扩散方程：
[0111][0112]
式中，h是水深(m)；ci是第i个泥沙成分的悬沙浓度(kg/m3)；u
x
和uy分别是x和y方向上垂线流速分量(m/s)；ε
s，x
和ε
s，y
是泥沙成分的涡动扩散系数(m2/s)；ei和di分别是第i个泥沙成分的悬浮和沉降通量(kg/m2/s)。
[0113]
非粘性沙的底部参考高度：
[0114][0115]
式中，fac是自定义的比例因子，ks是与水流相关的粗糙高度，δr是波浪引起的沙纹高度(0.025m)，h是水深。
[0116]
波流共同作用下的推移质输沙公式采用：
[0117][0118]
式中，sb是推移质输沙率(km/m/s)；d
50
是泥沙中值粒径；s是泥沙容重；ws是泥沙沉速；u
eff
是基于有效波高的垂线平均流速和近底计算层速度联合计算出的速度值；u
cr
是基于谢尔兹曲线的非粘性泥沙的临界垂线平均速度；g是重力加速度。
[0119]
模型中地貌形态演变在水流计算中进行更新，泥沙和动力地貌的计算均与水流同时计算。该模型可以通过使用代表性的水动力设置和加速因子来模拟从数秒到数千年时间尺度上的动力地貌演化过程。
[0120]
动力地貌模型有三个开边界，一个是离岸边界，另外两个是与海岸垂直的接岸边
界。在应用动力地貌模型研究近海工程及其演变时，边界条件在设置和验证过程中通常会遇到困难。主要问题是开边界如何确定合适的边界条件。这是由于作用在模型域上的过程的组合使得水位和流速在沿岸方向变化。如果边界条件与这种分布不相匹配，则会产生边界扰动，模型计算时位于边界处的结果将不再准确。此外，由于动力地貌模拟的时间尺度一般为几十年到上千年，若采用常用的水位时间过程边界，则数据量过大，且容易出现边界上流线紊乱的问题，本模型中a-a’和b-b’两个接岸边界采用neumann边界，只需提供频率、相位和振幅三个参数，即可代表任意时刻潮位变化数据，不仅可以解决在潮波在传递过程中对开边界的扰动，而且解决了中长时间尺度计算时，开边界数据量过大的复杂难题。
[0121]
根据下列公式给出接岸边界的水位过程：
[0122][0123]
其中：为振幅(m)；表示频率(rad/h)；表示波浪数(rad/m)；d
ab
代表a和b之间的距离，则a和b之间相位的差别可表示为lj是波长，u是流速，对于浅水流速可以根据求得。
[0124]
对上式方程两边求导可得接岸边界设置为纽曼型边界的动力方程：
[0125][0126]
开边界上的泥沙通量通设置为给定零含沙量梯度。糙率根据沉积物粒径分布和水深计算，其值在0.014～0.018之间。
[0127]
③
将波浪计算模块与已建立的动力地貌模型进行耦合，在潮波和泥沙模块中加入波浪的影响
[0128]
波浪的影响很难在动力地貌模型中加入，因为它们的出现往往不能被准确预测。本发明采用波作用量守恒的波浪相位平均模型作为控制方程，它主要描述波浪能量谱或波作用量的变化，是一种统计的方法。基于离散波作用平衡方程，通过耦合其它计算模块，在水流和泥沙模块中加入波浪的影响。
[0129]
波作用量与波能密度之间关系为：
[0130][0131]
式中，σ，α分别为角频率和波向。
[0132]
控制方程，即波作用量守恒方程在笛卡尔坐标系下表达式为：
[0133][0134]
式中，第一项代表波作用量随时间的变化率；第二项和第三项分别代表了波作用量在坐标空间x、y上的变化；第四项为水深和流速的变化引起的波作用量在频率上的变化；第五项表示由于水深和流场变化引起的折射，即波作用量在空间方向的传播；右侧s表示波作用量的源汇项，包括风输入项、波与波之间的非线性相互作用、白帽耗散等。c
x
、cy、c
σ
、c
α
分别表示在x、y、σ、α空间的波浪传播速度，可以采用线性波理论进行计算。
[0135]
风能的输入，即风能向波浪的转化主要采用共振机制和反馈机制来描述，可以表
示为线性增长和指数增长之和：
[0136]
sw＝g+pe(σ，α)
ꢀꢀꢀ
(23)
[0137]
式中，g、p与波频、波向、风速、风向有关。g为线性增长项。p为指数增长项。
[0138]
波能耗散机理主要包括三个方面：白帽耗散s
ds，w
、底摩擦耗散s
ds，b
以及由于水深变浅引起的波浪破碎s
ds，br
。深水中波浪表面的破碎现象称之为白帽耗散，主要由波陡决定。目前的白帽耗散公式主要基于脉动平均模型进行求解：
[0139][0140]
式中，γ为波陡系数；为平均角频率；为平均波数；k为各组分波数。
[0141]
低摩擦模型则包括三种模型：jonswap经验模型、collins拖拽理论模型、madsen等的涡黏模型。三折均可表示为：
[0142][0143]
式中，c
bottom
为与波浪质点底部轨道速度有关的底摩擦系数。
[0144]
由于水深变浅导致的波浪破碎的过程较为复杂，采用单位时间内谱分量的耗散表达式：
[0145][0146]
式中，e
tot
为总波能；d
tot
为单位水平面积上由波浪破碎导致的波能平均耗散率。
[0147]
不同频率间的波浪由于非线性相互作用发生能量交换，深水中四阶波-波非线性相互作用控制着波谱的变化，可以使波能从谱峰向低频(使谱峰频率减小)和高频(产生白帽耗散)转移；在浅水中，三阶波-波非线性相互作用可以使波能从低频向高频转移。
[0148]
四阶波-波非线性相互作用主要采用离散作用近似(dia)方法计算；三阶波-波非线性相互作用主要基于双线性波的离散三阶近似(dta)模型进行计算。
[0149]
由于计算时边界水流条件从水动力模块中获得，其结果从水动力网格插入波浪网格，所以通常波浪模型的网格必须大于水动力模型网格。
[0150]
重现波浪采用多代表波方法来确定代表波的波要素。将观察到的波浪数据进行分组，分为等距离的8个方向*8个波高。通过加权平均计算出每个等级中相应的有效波高和方向，然后在模型中将时间分成恒定长度的若干段，并根据波玫瑰图中的出现频率用代表波要素与其对应。波段长度设置为一天(两个潮汐周期)，并且将这些分好的段重新随即排序以减少由潮位变化和波浪条件的系统排序引起的误差。最后，根据重建的波要素时间序列与动力地貌模型相结合，每计算一个小时，计算文件将互相交换信息，作为下一步计算的初始条件。
[0151]
(4)高海况下沙坝潟湖海域初始海床地形精细化设置
[0152]
若某沙坝潟湖型海域床沙中值粒径为a～b微米(一般a＜60表示存在粘性沙，b＞100表示存在非粘性沙)。模型可以设置泥沙组分为，粘性沙三组m1，m2，m3，中值粒径分别为a，60微米；非粘性沙三组s1，s2，s3，中值粒径分为100，b微米。
[0153]
沙坝潟湖海域包括潮滩、潟湖、沙坝及潮汐深槽等地貌结构，其中潮滩与潟湖内以粉沙和淤泥为主，易在大风浪条件下发生剧烈运动，导致较大的冲淤变化，沙岛周边和深槽
存在沙质泥沙，较难起动。为了准确模拟沙坝潟湖海域动力地貌过程，必需对海床初始泥层厚度进行精细化设置。
[0154]
首先需要确定泥层总厚度δ，δ采用计算海域内近50年来最大一次风暴潮(大风天)发生后持续时间内最大冲刷深度。
[0155]
其次需要确定泥层中每种粒径的泥沙厚度，将海床实测点的中值粒径在海域上根据模型网格进行差值得到海域内海床上任意一点的中值粒径，假设海床某处实测中值粒径为c微米，若c介于a和微米之间，则对应选择m1，m2两种泥沙组分进行代表，m1的厚度设置为m2的厚度设置为若c介于和60微米之间，则对应选择m2，m3两种泥沙组分进行代表，m2的厚度设置为ms的厚度设置为若c介于60～100微米之间，则对应选择m3，s1两种泥沙组分进行代表，m3的厚度设置为s1的厚度设置为若c介于微米之间，则对应选择s1，s2两种泥沙组分进行代表，s1的厚度设置为s2的厚度设置为若c介于和b微米之间，则对应选择s2，s3两种泥沙组分进行代表，s2的厚度设置为s3的厚度设置为
[0156]
(5)高海况下沙坝潟湖海域动力地貌演变预测
[0157]
为了采用建立的模型预测高海况下沙坝潟湖海域动力地貌演变，需给出高海况条件。
[0158]
老龙沟潮滩与潟湖内以粉沙和沙质粉沙为主，易在风浪条件下发生剧烈运动，导致较大的冲淤变化，沙岛周边和深槽存在沙质泥沙。老龙沟海域经常发生7级以上大风，例如2008年12月发生8～9级大风，最大风速达16m/s，平均风速约10m/s。为了研究大风浪天气引起的滩槽冲淤演变，采用常浪向为se向偏东6度10年一遇波浪和大、中潮组合进行了规划方案实施后老龙沟冲淤演变计算，考虑到风速的实际过程，最大瞬时风速一般作用时间有限，波浪作用时间为48h，计算地形采用2010年实测地形。
[0159]
在中潮位为0.06m(理论基面)时，10年一遇有效波高分布见图2。由图2可见，10年一遇有效波高在外海最大可达4m左右，至老龙沟口门处衰减为3m左右，进入口门后波高迅速减小至1m以内，沙岛北侧浅滩范围波高小于0.5m。同时得到落急时刻含沙量场(2，可见10年一遇波浪情形下，潮汐汊道系统口门最大含沙量可达2kg/m3，进入口门后含沙量在0.5～1kg/m3。根据对不同风浪条件下老龙沟口门处含沙量资料可知，在无风浪时含沙量一般小于0.1kg/m3，多年平均波浪下含沙量在0.3～0.5kg/m3，常风天时最大可达1～1.5kg/m3，可见风浪对水体含沙量影响十分显著。
[0160]
分别计算中潮位和高潮位下10年一遇波浪作用下海域动力地貌演变结果，结果(图3)。可见，由东至西沿着月坨、东坑坨以及曹妃甸陆域围垦区形成了一条较为明显的冲刷带，也是该海域冲淤变化最为剧烈的区域。最大冲刷发生在月坨～东坑坨沙岛链一线，最大冲深达1m以上。东坑坨沙岛外侧以及口门西侧浅滩冲刷深度在1m左右，拦门沙冲刷幅度较小，基本在0.5m以内。而深槽以及口门-5m等深线外侧海域内出现淤积，包括东支以及西支深槽，西槽和东槽以及外海的淤积厚度大多在0.3m以内。
[0161]
为验证模型计算的可靠度，取沙坝外侧0m～10m等深线的断面，根据计算前后床面高程画出对比图。若剖面从沙坝外侧向海方向先冲刷后淤积，则与风暴剖面实验结果相符，结果(图4、5)表明预测的高海况下沙坝潟湖海域动力地貌演变结果可信。

技术特征：
1.一种高海况下宽级配沙坝潟湖海域动力地貌数值模拟方法，其特征在于，包括：s1利用层移条件下悬沙层含沙量剖面的计算公式表征波浪作用下悬沙层含沙量分布；s2建立高海况下波流挟沙能力的计算公式：波流挟沙能力采用水流挟沙能力和波浪挟沙能力之和表征，其中：水流挟沙能力通过实测数据和水力因子拟合获取；波浪挟沙能力根据床面形态划分：层移床面下，将所述含沙量剖面的计算公式根据水深积分，用以表征波浪挟沙能力；沙波床面下采用沙波床面波浪挟沙能力公式计算波浪挟沙能力；s3建立沙坝潟湖型海岸大范围潮波模型及局部动力地貌数学模型，并引入波浪模块建立耦合模型，利用s2获取的计算公式表征波流挟沙能力参数，模拟高海况下宽级配沙坝潟湖海域动力地貌。2.根据权利要求1所述的方法，其特征在于，所述s2中，基于波浪移动数判断层移临界条件，确定床面形态；当ψ≥226.21β
c-0.89
red
*-0.274
时，利用根据水深积分的含沙量剖面表征波浪挟沙能力；ψ＜226.21β
c-0.89
red
*-0.274
时，采用沙波床面波浪挟沙能力公式计算波浪挟沙能力；其中，ψ为波浪移动数，β
c
为密实度系数，red
*
为泥沙雷诺数。3.根据权利要求1所述的方法，其特征在于，对所述含沙量剖面的计算公式在层移层顶部高度至水面的水深范围内积分。4.根据权利要求1或3所述的方法，其特征在于，层移床面下波浪挟沙能力计算公式如下：式中，式中，式中，h为水深；α
fi
为考虑层化效应的床面平坦时的平均rouse系数，i＝1,2,3；和由每层的平均含沙量分别通过层化效应和阻滞沉速公式计算得到；z
i
为高度，c
zi
为z
i
高度处含沙量；z
a
为层移层顶部高度；c
a
为层移层顶部含沙量；δ
w
为波浪边界层厚度；r
s
＝0.8；a
f
＝1.17。5.根据权利要求1所述的方法，其特征在于，所述s3中，所述大范围潮波模型以基于浅水假设及boussinesq假设的不可压缩流体navier-stokes方程作为水动力控制方程。6.根据权利要求1或5所述的方法，其特征在于，大范围潮波模型的开边界设置于断面平均水深＞2m处，且开边界与研究区域距离＞10km。
7.根据权利要求1所述的方法，其特征在于，所述s3中，局部动力地貌数学模型模拟时，将泥沙组分按照中值粒径范围a～b划分为三组粘性沙和三组非粘性沙，其中三组粘性沙的中值粒径分别为a，60，三组非粘性沙的中值粒径分别为100，b；单位均为μm。8.根据权利要求1所述的方法，其特征在于，所述s3中，局部动力地貌数学模型模拟时，各泥沙组分的厚度基于如下方式确定：将海床实测点的中值粒径在海域上根据模型网格进行插值，得到海域内海床上任意一点的中值粒径；根据各点中值粒径和泥沙组分中值粒径数值范围，确定代表泥沙组分并获取泥沙组分厚度：d2＝δ-d1其中，d1为组分1厚度，c1为组分1中值粒径，d2为组分2厚度，c2为组分2中值粒径，c为海床上任意一点的中值粒径，c1＜c＜c2，δ为床面可动层厚度。9.根据权利要求8所述的方法，其特征在于，局部动力地貌数学模型计算时，基于沙坝潟湖型海岸近n年间最大一次风暴潮发生后持续时间内最大冲刷深度确定床面可动层厚度，n＞40。10.根据权利要求1所述的方法，其特征在于，所述s3中，局部动力地貌数学模型包含一个离岸边界和两个与海岸垂直的接岸边界；所述离岸边界采用水位边界，所述接岸边界采用neumann边界。

技术总结
本发明涉及一种高海况下宽级配沙坝潟湖海域动力地貌数值模拟方法，包括利用层移条件下悬沙层含沙量剖面的计算公式表征波浪作用下悬沙层含沙量分布；采用水流挟沙能力和波浪挟沙能力之和表征波流挟沙能力，其中波浪挟沙能力根据床面形态划分：层移床面下，将所述含沙量剖面的计算公式根据水深积分，用以表征波浪挟沙能力；建立沙坝潟湖型海岸大范围潮波模型及局部动力地貌数学模型，并引入波浪模块建立耦合模型，模拟高海况下宽级配沙坝潟湖海域动力地貌。本发明的方法可以同时适用于沙波及层移质运动，并且同时适用于粉沙-沙的宽级配海域，采用本发明提出的挟沙能力公式可以更准确地模拟高海况下宽级配海域的动力地貌过程。确地模拟高海况下宽级配海域的动力地貌过程。确地模拟高海况下宽级配海域的动力地貌过程。


技术研发人员：朱昊 王茂枚 刘嘉琦 张玉龙 王宏伟 李寿千 徐毅 姓海涛 陈颖 寇依诺 高业何敏 王振兵
受保护的技术使用者：集美大学
技术研发日：2023.07.11
技术公布日：2023/9/9