最新MIKE使用说明.docx
《最新MIKE使用说明.docx》由会员分享,可在线阅读,更多相关《最新MIKE使用说明.docx(29页珍藏版)》请在冰点文库上搜索。
最新MIKE使用说明
MIKE-使用说明
1根本参数与设置
1.1Mike21FlowModule
Mike21流体模型用来模拟二维自由外表流。
适用于水平尺度远大于垂向尺度,垂向流速和垂向加速度可以忽略时,湖泊、河口、海湾、海岸和海洋的水动力、环境现象的模拟。
共分为4个模块:
∙水动力模块〔Hydrodynamic〕;
∙平流扩散模块〔Advection-Dispersion〕;
∙泥沙输运模块〔MudTransport〕;
∙生态过程模块〔ECOLab〕。
其中,水动力模块是根底,为其他三个模块的计算提供动力。
泥沙输运模块可以用来模拟波流共同作用下粉砂、淤泥和粘土的冲刷、输移与沉降。
适用范围:
矩形网格。
1.2MIKE21/3CoupledModelFM
MIKE21/3耦合流体模型是海岸河口地区精确的动态模拟系统,包括以下模块:
∙水动力模块〔HydrodynamicModule〕
∙输运模块〔TransportModule〕
∙生态过程、溢油模块〔ECOLab/OilSpillModule〕
∙淤泥输运模块〔MudTransportModule〕
∙粗砂输运模块〔SandTransportModule〕
∙粒径追踪模块〔ParticleTrackingModule〕
∙波谱模块〔SpectralWaveModule〕
其中,水动力模块与波谱模块是最根本的两种。
适用范围:
三角网格。
1.3求解方法
ADI算法〔AlternativeDirectionImplicitMethod〕:
交替方向隐式方法。
把每一个时间步长分成两步进行,前半步隐式计算x方向流速分量及潮位,显式计算y方向流速分量;后半步隐式计算y方向流速分量及潮位,显式计算x方向流速分量。
1.4Bathymetry
地形中水深采用当地平均海平面以下水深,而非海图水深,故需在海图水深上增加2.9m。
〔timestepinterval〕,使CFLnumber小于1。
在浅水方程、泥沙输运方程计算中,为了使所有计算点的CFL数小于临界值,使用了变时间步长,即对每一步进行了细化计算。
计算对泥沙输运方程的CFL数要求没有浅水方程那么严〔可以大一些〕,因此,在泥沙输运方程计算时可以使用较大的时间步长。
用户可以指定最小时间步长〔0.01s〕和最大时间步长〔10s〕来控制时间步数。
最大时间步长〔maximumtime〕不可超过前面“time〞中定义的时间步长〔timestepinterval〕。
1.5Floodanddry
即干湿判别技术,是为了提高模拟的精度。
假设某网格水深小于dryingdepth那么不参与计算,当水深大于floodingdepth时重新进入计算。
1.6Eddyviscosity〔涡粘性〕
三角网格可选择两种方式:
Ø常数涡流公式〔㎡/s〕
ØSmagorinsky公式〔需指定Smagorinsky系数及最大值、最小值〕。
矩形网格可选择四种方式:
Ø忽略
Ø每个区域给定一个常值
Ø每个区域给定一个dfs2文件
Ø通过Smagorinsky公式进行动态计算〔用户需选择基于流速或基于通量,并指定比例因子〕
1.7Boundaryconditions
水动力与波浪边界条件均为每半小时给一次,意为这半小时内全采用相同的边界条件,这样做的前提是边界条件在半小时之内根本无变化。
边界条件的时间范围必须大于等于模型模拟的时间,边界条件的时间步长不必与模拟时间步长对应,模型会自动根据模拟时间步长插值。
潮流模型边界条件取中潮过程〔2022年9月8日2:
00:
00~9月9日2:
30:
00〕水位、流量。
波浪模型波浪边界条件由东中国海大模型提供,由于波浪边界条件对关注区域影响甚微,故可任意选取;风场为常风天6.1m/s,45°。
进行潮流模型、波浪模型验证的目的是给泥沙输运模型提供边界条件,验证好的潮流模型的水位、流量边界条件为泥沙模型提供潮流边界条件,验证好的波浪模型输出波浪场结果为泥沙输运模型提供波浪场,。
矩形网格的边界条件空间划分必须与边界网格个数一致,而三角网格的边界条件有插值功能,只需给出一个及以上的网格边界条件,沿线其他位置处的边界条件会插值得到。
1.8CFLnumber
浅水方程CFLnumber定义如下:
式中:
h——单元中心点的总水深〔totalwaterdepth〕;
u,v——单元中心点x和y方向的分速度;
Δx,Δy——计算单元x向和y向的特征长度,取边界最小值。
对于三角网格,可采用如下形式:
泥沙输运方程CFLnumber定义如下:
1.9Frequency
Frequency是指间隔多少时间步长输出一次结果。
1.10更换地形
每次更换地形后,系统会提醒是否更新边界条件,选“是〞的话,边界条件会自动恢复为默认值,需从新导入边界条件;建议选“否〞。
1.11输出结果格式
输出结果〔包括水动力、泥沙〕一般为半小时输出一次,因此时间间隔应设为180步。
假设设置过小会因内存缺乏导致运行中断。
1.12Mike21与Mike21/3输出结果项差异
2Flowmodel
2.1Bedresistance〔床面糙率〕
两种形式可供选择:
谢才系数〔Chezynumber〕与曼宁系数〔Manningnumber〕。
床面糙率采用下式计算
Bedresistance=
式中u是速度,C为谢才系数。
假设选择的是曼宁系数〔一般会这么做〕,模型会根据下式将其转换为谢才系数
。
因为多了此步,计算时间会变长。
Mike软件中曼宁系数的表达式为M=1/n,单位是m1/3/s,不同于水力学课本中常用的M=n〔无量纲〕。
在河流数值模拟中,如果某一边界的边界条件不是真实边界条件,计算过程中可能发生不稳定问题,这样建议增大边界区域的糙率,具体做法为沿河道方向将2~4行网格的曼宁系数设置为30。
〔1〕潮差验证时需修改糙率,假设是计算潮差较小,可加大下游曼宁系数〔减小下游糙率〕;〔2〕潮周期验证时要修改水深,因为使用的地形来自测图水深,为可能发生最小值,因此可适当增加水深。
2.2Initialconditions〔初始条件〕
三角网格时
可选择初始潮位或初始水深、流场,dfsu或dfs2文件均可,范围必须覆盖模型区域。
当只有一个时间步,就直接作为初始条件,当有多个时间步时,起始时间必须在第一步与最后一步之间,便于进行插值。
矩形网格时
可将每个区域的初始水位设置成定值或给定dfs2文件。
初始水位必须与边界条件相符合,例如边界条件初始值为0.5m,那么初始水位也应在0.5m附近。
3Wavemodel
3.1Whitecapping〔白浪〕
波高与Cdis成反比,波周期与DELTAdis成正比。
模型考虑白浪时,需指定两个参数:
Cdis和DELTAdis。
Cdis是白浪耗散源函数的比例因子,控制着耗散率的大小;DELTAdis控制着能量耗散在能量谱中的比重。
Cdis与DELTAdis的建议值分别为4.5和0.5。
减小Cdis会使白浪耗散整体减少,波高增大;DELTAdis的取值范围为0~1,DELTAdis小于0.5时,低频波耗散的比重会增加,导致波周期减小;DELTAdis大于0.5时,低频波耗散的比重会减小,导致波周期增大。
3.2Bottomfriction〔底摩阻〕
底摩阻引起的能量耗散会造成波谱频率降低,也就是平均波周期减小。
4Sandmodel
4.1ScientificDoc
Mud是指泥沙粒径小于63微米的细颗粒和粘性泥沙,Clay与Silt均属于Mud。
每层床面可以用下面几项描述:
临界冲刷切应力、冲刷系数〔E〕、冲刷功率、泥沙干密度和冲刷函数。
冲刷系数是控制冲刷速率的比例因子,软泥通常取0.000005~0.00002kg/m2/s,硬泥取0.0001kg/m2/s左右;冲刷功率〔无论软泥还是硬泥〕建议取4~26。
为了描述床层间的交换率,将固结考虑进来。
为了描述波浪引起的床面破坏,考虑了波浪液化作用。
4.2Parameterselection〔参数选择〕
可以选择泥沙粒径组的数量以及床面的层数。
粒径组的数量最多为8种,床层的数量最多为12层,原那么上,床层的数量应反映出床面的抗冲强度变化。
4.3水体参数
4.4初始条件
初始条件包括三项:
Fractionconcentrations〔初始水体含沙量〕、Layerthickness〔初始床面厚度〕、Fractiondistribution〔床面初始粒径组分配〕。
Initialconcentrations是粒径组的含沙量,这里不同的粒径组分别给定含沙量,可是并没有涉及到泥沙粒径的差异,这与利用挟沙力理论推求航道回淤时含沙量公式中没有考虑泥沙粒径相似[4],那么粒径组是以什么划分的?
与干密度有关。
initialbedthickness指每一层床面的厚度;initialfractiondistributioninbed为每一种粒径组在每一层床面所占的百分比。
各组泥沙初始含沙量,及其在每一床层〔一般为1层〕的百分比;
各组泥沙沉速,可考虑絮凝沉降或受阻沉降情况。
4.5Bedshearstress〔床面切应力〕
决定泥沙起动的关键是床面切应力的大小,在模型中,床面切应力采用以下公式进行计算[2]:
1〕当只有水流时,床面平均切应力:
〔5-4〕
2〕当纯波浪作用时,床面平均切应力:
〔5-5〕
3〕波流共同作用时的床面平均切应力:
〔5-6〕
式中:
为水流摩阻系数;
为水流平均速度;
为波浪摩阻系数;
为波浪底部水质点水平运动速度;
为波浪参数。
由以上三式可知,床面切应力只与水流、波浪、水深、床面糙率高度有关,可以从模型中输出结果。
4.6Concentrationprofile〔含沙量剖面〕
在设置临界淤积切应力场之前,必须指定含沙量剖面。
含沙量剖面分布形式有两种:
Rouseprofile和Teeterprofile。
Rouseprofile剖面含沙量不随时间变化,Teeterprofile剖面含沙量会随时间变化〔近底含沙量和深度平均含沙量的关系会根据重力、升力的大小重新计算,深度平均含沙量是如何参与计算的?
〕。
其对泥沙淤积速率的影响表现在cib〔近底含沙量〕与
取值会有所不同。
4.7Settling〔沉降〕
悬泥沉降过程可分为4个阶段:
自由沉降〔速度不变〕、絮凝沉降、受阻沉降、浮泥。
泥沙沉速开展过程线〔考虑受阻沉降〕
各阶段对应的沉速分别为
自由沉降:
絮凝沉降:
受阻沉降:
式中
——泥沙干密度;
——临界絮凝浓度;
——泥沙各层总浓度;
——临界受阻沉降浓度;
——沉速系数;
——功率。
4.8Deposition〔淤积〕
泥沙沉降是指泥沙从水体落到床面,当床面切应力小于临界淤积切应力时发生,第i组泥沙的淤积速率可以用下式表达:
式中,
——为泥沙单位时间、单位面积回淤量〔单位是kg/(m2·s)〕;
ωis——近底泥沙沉速。
软件注释中专门指出,对于细颗粒粘性泥沙〔finegrainedcohesivesediment(<0.004mm)〕,考虑絮凝沉降时,沉速分为三种情况:
含沙量小于临界絮凝浓度、在临界絮凝浓度与临界受阻沉降浓度之间、大于临界受阻沉降浓度,分别采用不同沉速;
cib——近底含沙量〔
,
,к=0.4,Uf为摩阻流速〕;
——淤积概率〔
〕。
可见:
泥沙淤积速率仅与近底泥沙沉速、近底含沙量及床面切应力、临界淤积切应力有关,没考虑泥沙粒径及级配。
注:
近底含沙量的取值取决于深度平均含沙量、淤积概率、近底泥沙沉速及摩阻流速,而近底泥沙沉速的取值完全取决于近底含沙量,两者是相互影响的。
综上,泥沙回淤速率最终取决于:
深度平均含沙量,床面切应力、临界淤积切应力、近底泥沙沉速及摩阻流速。
4.9Erosion〔冲刷〕
床面可选择软泥、局部固结床层〔softmud〕或者硬固结层〔hardmud〕,
软泥、局部固结层的冲刷速率可由下式表示
式中,E——冲刷系数〔erosioncoefficient〕〔表征床面的可冲性大小,单位是kg/(m2·s);
α——冲刷率〔poweroferosion〕;
τce——临界冲刷切应力。
硬固结层的冲刷速率可由下式表示
式中,E0是冲刷系数,Em是冲刷率,PE是冲刷概率。
软泥的抗冲强度沿深度是增加的,即临界冲刷切应力会随着冲刷深度的不同而改变,称为第二种冲刷类型[19];硬泥的抗冲强度沿深度不变,称为第一种冲刷类型。
矩形网格模型中使用的是softmud,进行三角网格计算时,小宝师兄用的是hardmud,之前用的是softmud。
由于连云港海域泥沙干密度在600~640kg/m3[2],属于硬泥〔Hardmud〕,因此取冲刷系数M=0.0001kg/(m2·s)〕;
只有一层河床时,由于每一层的床面参数为定值,所以每一层临界冲刷切应力场也为定值,不存在沿深度变化的问题。
冲刷、淤积速率的单位均是kg/(m2·s),模型中床面地形的改变是以m为单位,因此,上式结果还要除以表层淤积物的干密度。
被冲起的泥沙会根据床面分布,分配到不同的粒径组中。
式中Speedup为加速因子,其他参数意义同上。
总的来说,床面冲淤函数如下
式中:
为水流底部剪切应力;
为临界淤积剪切应力;
为临界冲刷剪切应力;
为淤积概率;M为冲刷系数〔kg/(m2·s)〕;
为泥沙沉速。
这里床面冲刷速率是针对hardmud来说的,softmud的冲刷速率仍为上条公式。
4.10Dispersion
三角网格时,水平扩散以三种不同的方式考虑
Ø无扩散
Ø扩散系数公式〔指定扩散系数㎡/s〕
Ø比例涡粘性公式〔水流涡粘系数乘以比例因子〕
矩形网格时,扩散系数可以设置为如下两种方式
Ø独立于流体〔指定扩散系数㎡/s〕
Ø与计算通量成正比〔每个方向指定一个比例因子〕
注意:
数学模型中的扩散系数取决于网格尺度、时间步长以及问题的物理本质。
4.11Bedroughness〔泥沙床面糙率〕
在计算床面切应力时,需要使用床面糙率。
一般来说,糙率高度kn可定义为2.5倍泥沙粒径,对于泥沙粒径d<0.5mm的较平坦床面,床面形态是引起糙率高度的主要因素,因此有学者建议直接取kn=0.001m[3]。
连云港海域糙率高度取kn=0.001m。
4.12Morphology〔地貌〕
长周期地形模拟时,假设是地形的改变量与水深处在同一量级,可能对水动力造成影响时,需要考虑地貌的变化,特别是浅水区域。
地貌更新方式如下:
式中
——第n步地形;
Zn+1——第n+1步地形;
netsedn——第n步到第n+1步的泥沙净通量。
如果模拟时间非常长,可以考虑使用加速因子提高模拟的精度。
此时:
4.13泥沙模型计算冲淤过程
软件计算航道冲淤时,首先根据床面切应力计算泥沙冲或淤的概率〔
〕,这取决于床面切应力与临界淤积、冲刷切应力的关系,然后利用冲、淤公式计算冲、淤厚度,对于航道而言,一般是只淤不冲。
计算并不区分航道或浅滩。
4.14泥沙模型率定
模型率定的关键为率定临界冲刷、淤积切应力场,其取值范围分别为τce=0.18~0.70N/m2,τcd=0.03~0.15N/m2[1]。
率定步骤:
首先假设临界淤积切应力场所有值均为0.05,根据大西山等测站的实测含沙量以及海床冲淤平衡原那么率定临界冲刷切应力场;然后根据各分段航道实测回淤数据率定航道内临界淤积切应力场。
在进行率定淤积切应力过程中发现,在平衡状态下,临界淤积切应力〔τcd〕与滩槽深度比〔d1/d2〕呈良好的关系,因此,在进行高等级航道回淤预测时,可以直接利用上面的关系率定淤积切应力。
上面是理想情况,实际上,上述关系并不合理。
一方面,建立关系所用数据较少,连云港航道仅考虑w1~外8段,没有考虑更靠外海的情况;另一方面,所用数据中,航道水深根本无变化,连云港航道w1~外5航道水深一致,外5~外8航道水深只有微小变化,徐圩航道各段水深均一致。
因此,与其说是建立临界淤积切应力〔τcd〕与滩槽深度比〔d1/d2〕的关系,不如说是建立临界淤积切应力〔τcd〕与边滩深度的关系。
该公式并不合理,无法推广。
现在,肖天葆[3]提出了一种新的率定方法:
基于床面切应力的临界冲刷、淤积切应力率定方法:
〔1〕对于航道外的临界冲刷切应力及临界淤积切应力,可根据动力塑造海床及海床冲淤平衡的原那么,并采用大西山等测站的实测多年平均含沙量进行率定。
〔2〕对于航道内的临界冲刷切应力,假定其与附近水深相同处的
值一致,而航道内的临界淤积切应力可采用连云港区不同等级航道的实测回淤强度进行率定。
从对模型的测试来说,临界冲刷切应力对含沙量的影响要大很多〔分析原因?
〕。
因此,本文主要根据临界冲刷切应力来确定含沙量,根据航道内临界淤积切应力来确定航道回淤。
在进行计算床面切应力之前,首先需要假设临界冲刷、淤积切应力场,建议将其设为定值,然后根据计算出来的床面切应力场计算临界冲刷、淤积切应力场。
具体方法见“连云港文献阅读情况〞。
注意:
临界冲刷切应力率定时,可以选择7万吨级航道开挖后的地形〔2022年〕,当然也可以选择其他等级航道地形,因为单个航道对大范围海域的影响可以忽略不计,航道外大范围海域的
率定完成以后就不再改变,航道内的
假设与外海水深相同处的
一致,航道内根本不发生冲刷,航道内
根据航道回淤实测资料进行率定。
4.15泥沙模型验证
模型验证时仅验证含沙量与回淤强度,含沙量验证可通过刘家驹公式验证,验证结果较合理;回淤强度验证可通过验证15万吨级回淤强度进行,由于时间限制,此步尚未实施。
4.16冲淤特点
平行于岸线的单元冲淤过程类似,这与潮流、波浪传播方向有关,也是造成沙坝的主要因素。
其中,上面三条曲线是单元t7、t8、t9冲淤过程,中间五条曲线是单元t1、t5、t6、t10、t11,下面三条曲线是单元t2、t3、t4冲淤过程。
岸线走向为西北-东南向。
单元
5
S5
7
S7
2
0.1297
0.1017
-0.1092
0.1017
-0.2118
3
0.3192
0.1095
-1
0.2095
-1
2905
1.0094
0.001
-0.648
0.101
4.16
2907
1.1954
先调临界淤积切应力,让只通过调节
就能调好的单元全部调好,之后再调临界冲刷切应力。
如果单元冲刷大于0.1m,那么
增加0.1,假设果单元淤积大于0.1m,那么
减去0.1,
最小为0.001,不能大于
最大为0.9。
两种情况较难解决:
〔1〕
很小时〔=0.001〕,床面仍旧发生淤积,一种方法是继续调小
,另一种方法是调小
,这里使用后一种方法。
如果si与s(i+1)同号,
增加或减小0.1时,|s〔i+1〕-si|>0.05,并且
>0.2,那么继续调
,否那么直接调
。
如果si与s(i+1)异号,那么根据截距法求临界淤积切应力。
当
=0.001,且s(i+1)<-0.1时,调大
。
对于单元169-175,s5-s7,临界淤积切应力大于临界冲刷切应力,但是在模拟时间段内,始终有冲淤过程。
单元169s5冲淤过程
对于单元3,淤积切应力增加0.1后,床面冲淤量几乎没有变化,仍然冲刷剧烈〔-1m〕,s7m继续将单元3临界淤积切应力增加0.1,效果仍不明显,建议下次调大临界冲刷切应力。
对于单元2,s5与s7m的临界淤积、冲刷切应力均相同,床面切应力过程根本一致〔见图1、2〕,但是冲淤过程却差异很大〔见图3、4〕。
图1单元2S5床面切应力变化过程
图2单元2S7m床面切应力变化过程
图3单元2S5冲淤过程
图4单元2S7m冲淤过程
对于单元2905,s5与s7m的
保持不变,临界淤积切应力增加0.1后,床面发生大量淤积〔冲淤过程见图7、8〕,床面切应力过程变化较大〔见图5、6〕。
图5单元2905S5床面切应力过程
图6单元2905S7m床面切应力过程
图7单元2905S5冲淤过程
图8单元2905S7m冲淤过程
对于单元14879,
为0.001,
为2.5,但是床面淤积达4.8536m,其床面切应力过程〔s11〕如图9,下次直接调
。
图9单元14879s11床面切应力过程
4.17数模预测回淤展望
使用数学模型预测回淤的难点如下:
〔1〕当航道等级变化后,怎么样调整临界淤积切应力场?
临界冲刷切应力场等其它参数用不用改变?
〔2〕假设是水动力与含沙量场也改变的情况下,应如何处理?
〔3〕由于率定参数的多样性,某种情况下的最正确组合可能有多种,怎么调整临界淤积、冲刷切应力场才能使航道冲淤对其敏感性降低,即更好的预测回淤。
理论公式〔刘家驹公式〕计算航道回淤时,依据为浅滩浑水跨越航道时,航道水深大,挟沙能力降低而造成落淤。
囧:
泥沙输运模块输出结果同样是没有勾选bathmetry,矩形网格模式下会输出此项,而三角网格模式下不会输出此项。
原因是什么?
5数据文件
5.1dfs0文件数据特征
dfs0文件即坐标点时间序列值文件。
同一三角网格内部不同坐标点的时间序列值完全相同,而同一矩形网格内部不同坐标点的时间序列值根本相同。
因此,坐标点时间序列值的大小只与其所在的网格有关,而与其所在该网格内部的具体位置无关。
网格值代表点为网格中心点。
5.2Dfs1文件
Dfs1文件为剖面时间序列文件,值随时间与空间变化,默认格式为矩形网格,可应用于三角网格。
时间范围应大于模拟时间;空间范围使用矩形网格时,边界点个数应等于边界网格数,使用三角网格时,边界点个数≥2。
地形信息可不考虑。
5.3矩形网格
(1)文件制作
File——new——file…——Mikezero——Bathymetries——制作绘图区域〔投影坐标系、起始点坐标、长度与高度〕——制作地形
ØWorkarea——backgroundmanagement——导入岸线、地形;
ØWorkarea——bathymetrymanagement——设置网格尺度与个数;
Ø点击“
〞(ImportfromBackground)选中研究区域——再次点击“
〞选中研究区域——bathymetrymanagement——interpolate…插值
Ø保存文件,地形制作完成。
(2)网格嵌套
File——new——file…——Mike21——Mike21toolbox——Hydrodynamics——Boarderadjustment——指定文件、起点坐标。
(3)文件处理
将源文件用Gridseries方式翻开——Tools——Copyfileintodata——导入目标文件,即可对文件进行处理。
5.4Gridseries、dataview翻开方式
用Gridseries方式与dataview方式翻开dfsu文件均可以显示单元格数值,每个单元格对应的值只有一个,这通过dataview下