MIKE-使用说明.docx

1 基本参数与设置1.1 Mike 21 Flow ModuleMike 21流体模型用来模拟二维自由表面流。适用于水平尺度远大于垂向尺度，垂向流速和垂向加速度可以忽略时，湖泊、河口、海湾、海岸和海洋的水动力、环境现象的模拟。共分为4个模块： 水动力模块（Hydrodynamic）； 平流扩散模块（Advection-Dispersion）； 泥沙输运模块（Mud Transport）； 生态过程模块（ECO Lab）。其中，水动力模块是基础，为其他三个模块的计算提供动力。泥沙输运模块可以用来模拟波流共同作用下粉砂、淤泥和粘土的冲刷、输移与沉降。适用范围：矩形网格。1.2 MIKE 21/3 Coupled Model FMMIKE 21/3耦合流体模型是海岸河口地区精确的动态模拟系统，包括以下模块： 水动力模块（Hydrodynamic Module） 输运模块（Transport Module） 生态过程、溢油模块（ECO Lab/Oil Spill Module） 淤泥输运模块（Mud Transport Module） 粗砂输运模块（Sand Transport Module） 粒径追踪模块（Particle Tracking Module） 波谱模块（Spectral Wave Module）其中，水动力模块与波谱模块是最基本的两种。适用范围：三角网格。1.3 求解方法ADI算法（Alternative Direction Implicit Method）：交替方向隐式方法。把每一个时间步长分成两步进行，前半步隐式计算x方向流速分量及潮位，显式计算y方向流速分量；后半步隐式计算y方向流速分量及潮位，显式计算x方向流速分量。1.4 Bathymetry地形中水深采用当地平均海平面以下水深，而非海图水深，故需在海图水深上增加2.9m。1.5 Time指模拟的时间段，包括起始时间（simulation start time）、总步数（no. of time steps）及时间步长（time step interval）。时间步长一方面决定了结果文件的最大输出频率，另一方面实现了不同模块的同步耦合。具体计算时的时间步长则是在solution technique中定义的。水动力、平流扩散和波谱模型的计算步数是动态的，只要满足稳定性要求即可；对于泥沙输运和水质模型，进程可能多个步长更新一回。1.6 Solution technique对于时间积分（time integration）与空间离散（space discretization）两项，一般选作低阶快速模式（lower order，fast algorithm），这种模式没有高阶模式（higher order）计算精确；对于平面二维计算，浅水方程、泥沙输运方程和平流扩散方程的时间积分采用的是显式格式，三维计算采用的是半隐式格式（水平方向显式，竖直方向隐式）。由于使用显式计算时的稳定性要求，必须指定时间间隔（time step interval），使CFL number小于1。在浅水方程、泥沙输运方程计算中，为了使所有计算点的CFL数小于临界值，使用了变时间步长，即对每一步进行了细化计算。计算对泥沙输运方程的CFL数要求没有浅水方程那么严（可以大一些），因此，在泥沙输运方程计算时可以使用较大的时间步长。用户可以指定最小时间步长（0.01s）和最大时间步长（10s）来控制时间步数。最大时间步长（maximum time）不可超过前面“time”中定义的时间步长（time step interval）。1.7 Flood and dry 即干湿判别技术，是为了提高模拟的精度。若某网格水深小于drying depth 则不参与计算，当水深大于flooding depth 时重新进入计算。1.8 Eddy viscosity（涡粘性）三角网格可选择两种方式： 常数涡流公式（/s） Smagorinsky公式（需指定Smagorinsky系数及最大值、最小值）。矩形网格可选择四种方式： 忽略 每个区域给定一个常值 每个区域给定一个dfs2文件 通过Smagorinsky公式进行动态计算（用户需选择基于流速或基于通量，并指定比例因子）1.9 Boundary conditions水动力与波浪边界条件均为每半小时给一次，意为这半小时内全采用相同的边界条件，这样做的前提是边界条件在半小时之内基本无变化。边界条件的时间范围必须大于等于模型模拟的时间，边界条件的时间步长不必与模拟时间步长对应，模型会自动根据模拟时间步长插值。潮流模型边界条件取中潮过程（2005年9月8 日2:00:009月9日 2:30:00）水位、流量。波浪模型波浪边界条件由东中国海大模型提供，由于波浪边界条件对关注区域影响甚微，故可任意选取；风场为常风天6.1m/s，45。进行潮流模型、波浪模型验证的目的是给泥沙输运模型提供边界条件，验证好的潮流模型的水位、流量边界条件为泥沙模型提供潮流边界条件，验证好的波浪模型输出波浪场结果为泥沙输运模型提供波浪场，。矩形网格的边界条件空间划分必须与边界网格个数一致，而三角网格的边界条件有插值功能，只需给出一个及以上的网格边界条件，沿线其他位置处的边界条件会插值得到。1.10 CFL number浅水方程CFL number定义如下：式中： h单元中心点的总水深（total water depth）； u,v单元中心点x和y方向的分速度；x，y计算单元x向和y向的特征长度，取边界最小值。对于三角网格，可采用如下形式：泥沙输运方程CFL number定义如下：1.11 FrequencyFrequency 是指间隔多少时间步长输出一次结果。1.12 更换地形每次更换地形后，系统会提醒是否更新边界条件，选“是”的话，边界条件会自动恢复为默认值，需从新导入边界条件；建议选“否”。1.13 输出结果格式输出结果（包括水动力、泥沙）一般为半小时输出一次，因此时间间隔应设为180步。若设置过小会因内存不足导致运行中断。1.14 Mike21与Mike21/3输出结果项差别2 Flow model2.1 Bed resistance（床面糙率）两种形式可供选择：谢才系数（Chezy number）与曼宁系数（Manning number）。床面糙率采用下式计算Bed resistance=式中u是速度，C为谢才系数。若选择的是曼宁系数（一般会这么做），模型会根据下式将其转换为谢才系数。因为多了此步，计算时间会变长。Mike软件中曼宁系数的表达式为M=1/n，单位是m1/3/s，不同于水力学课本中常用的M=n（无量纲）。在河流数值模拟中，如果某一边界的边界条件不是真实边界条件，计算过程中可能发生不稳定问题，这样建议增大边界区域的糙率，具体做法为沿河道方向将24行网格的曼宁系数设置为30。（1）潮差验证时需修改糙率，若是计算潮差较小，可加大下游曼宁系数（减小下游糙率）；（2）潮周期验证时要修改水深，因为使用的地形来自测图水深，为可能发生最小值，因此可适当增加水深。2.2 Initial conditions（初始条件）三角网格时可选择初始潮位或初始水深、流场，dfsu或dfs2文件均可，范围必须覆盖模型区域。当只有一个时间步，就直接作为初始条件，当有多个时间步时，起始时间必须在第一步与最后一步之间，便于进行插值。矩形网格时可将每个区域的初始水位设置成定值或给定dfs2文件。初始水位必须与边界条件相符合，例如边界条件初始值为0.5m，则初始水位也应在0.5m附近。3 Wave model3.1 White capping（白浪）波高与Cdis成反比，波周期与DELTAdis成正比。模型考虑白浪时，需指定两个参数：Cdis 和 DELTAdis。Cdis是白浪耗散源函数的比例因子，控制着耗散率的大小；DELTAdis控制着能量耗散在能量谱中的比重。Cdis与DELTAdis的建议值分别为4.5和0.5。减小Cdis 会使白浪耗散整体减少，波高增大；DELTAdis的取值范围为01，DELTAdis小于0.5时，低频波耗散的比重会增加，导致波周期减小；DELTAdis大于0.5时，低频波耗散的比重会减小，导致波周期增大。3.2 Bottom friction（底摩阻）底摩阻引起的能量耗散会造成波谱频率降低，也就是平均波周期减小。4 Sand model4.1 Scientific DocMud是指泥沙粒径小于63微米的细颗粒和粘性泥沙，Clay与Silt均属于Mud。每层床面可以用下面几项描述：临界冲刷切应力、冲刷系数（E）、冲刷功率、泥沙干密度和冲刷函数。冲刷系数是控制冲刷速率的比例因子，软泥通常取0.0000050.00002 kg/m2/s，硬泥取0.0001 kg/m2/s左右；冲刷功率（无论软泥还是硬泥）建议取426。 为了描述床层间的交换率，将固结考虑进来。为了描述波浪引起的床面破坏，考虑了波浪液化作用。4.2 Parameter selection（参数选择）可以选择泥沙粒径组的数量以及床面的层数。粒径组的数量最多为8种，床层的数量最多为12层，原则上，床层的数量应反映出床面的抗冲强度变化。4.3 水体参数4.4 初始条件初始条件包括三项：Fraction concentrations（初始水体含沙量）、Layer thickness（初始床面厚度）、Fraction distribution（床面初始粒径组分配）。Initial concentrations是粒径组的含沙量，这里不同的粒径组分别给定含沙量，可是并没有涉及到泥沙粒径的差别，这与利用挟沙力理论推求航道回淤时含沙量公式中没有考虑泥沙粒径相似4，那么粒径组是以什么划分的？与干密度有关。 initial bed thickness指每一层床面的厚度；initial fraction distribution in bed为每一种粒径组在每一层床面所占的百分比。各组泥沙初始含沙量，及其在每一床层（一般为1层）的百分比；各组泥沙沉速，可考虑絮凝沉降或受阻沉降情况。4.5 Bed shear stress（床面切应力）决定泥沙起动的关键是床面切应力的大小，在模型中，床面切应力采用以下公式进行计算2：1）当只有水流时，床面平均切应力：（5-4）2）当纯波浪作用时，床面平均切应力：（5-5）3）波流共同作用时的床面平均切应力：（5-6）式中：为水流摩阻系数；为水流平均速度；为波浪摩阻系数；为波浪底部水质点水平运动速度；为波浪参数。由以上三式可知，床面切应力只与水流、波浪、水深、床面糙率高度有关，可以从模型中输出结果。4.6 Concentration profile（含沙量剖面）在设置临界淤积切应力场之前，必须指定含沙量剖面。含沙量剖面分布形式有两种：Rouse profile和Teeter profile。Rouse profile剖面含沙量不随时间变化，Teeter profile剖面含沙量会随时间变化（近底含沙量和深度平均含沙量的关系会根据重力、升力的大小重新计算，深度平均含沙量是如何参与计算的？）。其对泥沙淤积速率的影响表现在cib（近底含沙量）与取值会有所不同。4.7 Settling（沉降）悬泥沉降过程可分为4个阶段：自由沉降（速度不变）、絮凝沉降、受阻沉降、浮泥。泥沙沉速发展过程线（考虑受阻沉降）各阶段对应的沉速分别为自由沉降：絮凝沉降：受阻沉降：式中 泥沙干密度； 临界絮凝浓度；泥沙各层总浓度；临界受阻沉降浓度；沉速系数；功率。4.8 Deposition（淤积）泥沙沉降是指泥沙从水体落到床面，当床面切应力小于临界淤积切应力时发生，第i组泥沙的淤积速率可以用下式表达：式中，为泥沙单位时间、单位面积回淤量（单位是kg/(m2s)）；is近底泥沙沉速。软件注释中专门指出，对于细颗粒粘性泥沙（fine grained cohesive sediment (0.004 mm)），考虑絮凝沉降时，沉速分为三种情况：含沙量小于临界絮凝浓度、在临界絮凝浓度与临界受阻沉降浓度之间、大于临界受阻沉降浓度，分别采用不同沉速；cib近底含沙量（，=0.4，Uf为摩阻流速）；淤积概率（）。可见：泥沙淤积速率仅与近底泥沙沉速、近底含沙量及床面切应力、临界淤积切应力有关，没考虑泥沙粒径及级配。注：近底含沙量的取值取决于深度平均含沙量、淤积概率、近底泥沙沉速及摩阻流速，而近底泥沙沉速的取值完全取决于近底含沙量，两者是相互影响的。综上，泥沙回淤速率最终取决于：深度平均含沙量，床面切应力、临界淤积切应力、近底泥沙沉速及摩阻流速。4.9 Erosion（冲刷）床面可选择软泥、部分固结床层（soft mud）或者硬固结层（hard mud），软泥、部分固结层的冲刷速率可由下式表示式中，E冲刷系数（erosion coefficient）（表征床面的可冲性大小，单位是kg/(m2s)； 冲刷率（power of erosion）；ce临界冲刷切应力。硬固结层的冲刷速率可由下式表示式中，E0是冲刷系数，Em是冲刷率，PE是冲刷概率。软泥的抗冲强度沿深度是增加的，即临界冲刷切应力会随着冲刷深度的不同而改变，称为第二种冲刷类型19；硬泥的抗冲强度沿深度不变，称为第一种冲刷类型。矩形网格模型中使用的是soft mud，进行三角网格计算时，小宝师兄用的是hard mud，之前用的是soft mud。由于连云港海域泥沙干密度在600640kg/m32，属于硬泥（Hard mud），因此取冲刷系数M=0.0001 kg/(m2s)）；只有一层河床时，由于每一层的床面参数为定值，所以每一层临界冲刷切应力场也为定值，不存在沿深度变化的问题。冲刷、淤积速率的单位均是kg/(m2s)，模型中床面地形的改变是以m为单位，因此，上式结果还要除以表层淤积物的干密度。被冲起的泥沙会根据床面分布，分配到不同的粒径组中。式中Speedup为加速因子，其他参数意义同上。总的来说，床面冲淤函数如下式中：为水流底部剪切应力；为临界淤积剪切应力；为临界冲刷剪切应力；为淤积概率；M为冲刷系数（kg/(m2s)）；为泥沙沉速。这里床面冲刷速率是针对hard mud来说的，soft mud的冲刷速率仍为上条公式。4.10 Dispersion 三角网格时，水平扩散以三种不同的方式考虑 无扩散 扩散系数公式（指定扩散系数/s） 比例涡粘性公式（水流涡粘系数乘以比例因子）矩形网格时，扩散系数可以设置为如下两种方式 独立于流体（指定扩散系数/s） 与计算通量成正比（每个方向指定一个比例因子）注意：数学模型中的扩散系数取决于网格尺度、时间步长以及问题的物理本质。4.11 Bed roughness（泥沙床面糙率）在计算床面切应力时，需要使用床面糙率。一般来说，糙率高度kn可定义为2.5倍泥沙粒径，对于泥沙粒径d0.05，并且ce-cd0.2，则继续调cd，否则直接调ce。如果si与s(i+1)异号，则根据截距法求临界淤积切应力。当ce-cd=0.001，且s(i+1)-0.1时，调大ce。对于单元169-175，s5-s7，临界淤积切应力大于临界冲刷切应力，但是在模拟时间段内，始终有冲淤过程。单元169 s5冲淤过程对于单元3，淤积切应力增加0.1后，床面冲淤量几乎没有变化，仍然冲刷剧烈（-1m），s7 m继续将单元3临界淤积切应力增加0.1，效果仍不明显，建议下次调大临界冲刷切应力。对于单元2，s5与s7 m的临界淤积、冲刷切应力均相同，床面切应力过程基本一致（见图1、2），但是冲淤过程却差别很大（见图3、4）。图1 单元2 S5 床面切应力变化过程图2 单元2 S7 m 床面切应力变化过程图3 单元2 S5 冲淤过程图4 单元2 S7 m 冲淤过程对于单元2905，s5与s7 m的 ce保持不变，临界淤积切应力增加0.1后，床面发生大量淤积（冲淤过程见图7、8），床面切应力过程变化较大（见图5、6）。图5 单元2905 S5 床面切应力过程图6 单元2905 S7 m 床面切应力过程图7 单元2905 S5 冲淤过程图8 单元2905 S7 m 冲淤过程对于单元14879，cd为0.001，ce为2.5，但是床面淤积达4.8536m，其床面切应力过程（s11）如图9，下次直接调ce。图9 单元14879 s11 床面切应力过程4.17 数模预测回淤展望使用数学模型预测回淤的难点如下：（1）当航道等级变化后，怎么样调整临界淤积切应力场？临界冲刷切应力场等其它参数用不用改变？（2）若是水动力与含沙量场也改变的情况下，应如何处理？（3）由于率定参数的多样性，某种情况下的最佳组合可能有多种，怎么调整临界淤积、冲刷切应力场才能使航道冲淤对其敏感性降低，即更好的预测回淤。理论公式（刘家驹公式）计算航道回淤时，依据为浅滩浑水跨越航道时，航道水深大，挟沙能力降低而造成落淤。囧：泥沙输运模块输出结果同样是没有勾选bathmetry，矩形网格模式下会输出此项，而三角网格模式下不会输出此项。原因是什么？5 数据文件5.1 dfs0文件数据特征dfs0文件即坐标点时间序列值文件。同一三角网格内部不同坐标点的时间序列值完全相同，而同一矩形网格内部不同坐标点的时间序列值基本相同。因此，坐标点时间序列值的大小只与其所在的网格有关，而与其所在该网格内部的具体位置无关。网格值代表点为网格中心点。5.2 Dfs1文件Dfs1文件为剖面时间序列文件，值随时间与空间变化，默认格式为矩形网格，可应用于三角网格。时间范围应大于模拟时间；空间范围使用矩形网格时，边界点个数应等于边界网格数，使用三角网格时，边界点个数2。地形信息可不考虑。5.3 矩形网格（1） 文件制作FilenewfileMike zeroBathymetries制作绘图区域（投影坐标系、起始点坐标、长度与高度）制作地形 Work areabackground management导入岸线、地形； Work areabathymetry management设置网格尺度与个数； 点击“”( Import from Background)选中研究区域再次点击“”选中研究区域bathymetry managementinterpolate插值 保存文件，地形制作完成。（2） 网格嵌套FilenewfileMike 21Mike 21 toolboxHydrodynamicsBoarder adjustment指定文件、起点坐标。 （3） 文件处理将源文件用Grid series 方式打开ToolsCopy file into data导入目标文件，即可对文件进行处理。5.4 Grid series、data view打开方式用Grid series 方式与data view方式打开dfsu文件均可以显示单元格数值，每个单元格对应的值只有一个，这通过data view下box contour显示方式可以发现，但我们一般选择shaded contour显示方式，在这种方式下，大范围场值会自动插值，从而每个单元的值有无数个。此外，Grid series下，除了显示插值后的单元值外，还会显示地形值。5.5 dfsu文件制作可用于制作糙率文件，临界冲刷、淤积切应力场文件等。FilenewfileMike zerodata manager。首先自文件导入网格，可选择mesh文件或者dfs2、dfs3文件。mesh文件数据值是网格结点值，而.dfs2、.dfs3文件数据值是网格中心点值，均可用来制作原始网格。制作糙率文件时，需要根据运行结果中的still water depth，利用公式计算，然后赋值到空白网格内；制作临界切应力场文件时，需要用到切应力过程最大值或最小值，乘以系数，进行赋值。5.6 结果提取矩形网格结果提取方式：FilenewfileMike zeroMike zero toolboxStatisticstxstatsetup name（可填名称或选默认，点下一步）statistics on matrixselect grid file（选取源文件，点下一步）设置提取结果的时间范围（起始步和终止步）设置提取结果的空间范围设置提取结果的类型设置输出文件的位置及名称。三角网格结果提取方式：FilenewfileMike zeroData extraction FM在input与outputs里设置参数。输出结果文件的时间项包括起始步数、最终步数以及输出频率（frequency）。其中，输出频率是指每隔多少步输出一次。例如：时间步长为30s，输出频率为60，输出结果为每隔半小时保存一次。之所以设定输出频率，是为了节省计算机存储空间。每一步输出结果均是基于这半小时内的值，一般取这半小时之内的平均值。例如：连云港地区波周期为3s，半小时内共发生600次波峰，输出的有效波高就是取自这600个波高中前1/3大波高的平均值。有效波高周期同潮周期相似，均为12h。5.7 不同网格数量dfsu文件插值转换 Toolsexport导出网格值、水深值； 将其制成xyz文件（删掉水深值，保留网格值）； 将该文件导入新的网格文件中（在sms里），转化为mesh文件； 在模型中运行，输出still water depth，即为制作好的新dfs2文件。6 试验方案6.1 小宝、刘杰率定大范围冲淤平衡模型参数比较小宝刘杰水动力边界40570床面泥沙密度(kg/m)620600Dispersion（scaled eddy viscosity formulation）1.10.56.2 问：模型计算要求波浪场网格必须与地形网格一致吗？答：不必要。波浪场网格可以不一样，但初始场（水动力、含沙量）、临界切应力场的网格必须一样；见huyan-5calculation-15wan-“15wan mud lj 3(changing mesh).mfm”。6.3 问：航道内临界淤积切应力是否影响床面切应力？答：不影响，见huyan-8result-15wan-“15w-mud-45d-6 lj 4 航段4中心点stress.dfs0 VS 15w-mud-45d-6 lj航段4中心点stress.dfs0”。6.4 问：影响床面切应力的主要因素是什么？答：水深。首先，边滩上切应力明显大于航道内切应力，见“15w-mud-45d-6 lj 4 航段7航槽及边滩stress.dfs0”；其次，随着航道等级提高（航深加大），床面切应力减小，尤其是最大值的减小，最小值有可能基本不变或增大，见“航道内床面切应力比较.xlsx外6、9段15wan、25wan比较”；最后，随着边滩水深的增加，航道内床面切应力增加，见“航道内床面切应力比较.xlsx15万沿程stress”。6.5 问：初始水动力场、波浪场对航道回淤或含沙量场的影响大吗？初始含沙量场呢？答：初始水动力场、波浪场对航道回淤或含沙量场的影响并不大，见huyan-8result-15wan-“15w-mud-45d-6 lj 4-1、4-3”，初始含沙量场若不小心设成负数，则会大幅度影响含沙量场及航道回淤。6.6 矩形、三角网格计算航道回淤时，对临界淤积切应力的敏感性一样吗？答：不一样。矩形网格基本呈线性关系，三角网格敏感性远大于矩形网格，见huyan-8result-25wan- Sensitivity test“Sensitivity.xlsx” VS nisha-7结果- 25wan tcd敏感性测试。6.7 三角网格、矩形网格航道回淤与床面切应力的关系对于25万吨级外4段、外8段，矩形网格床面切应力过程极小值小于三角网格，但是根据回淤量率定临界淤积切应力却大于三角网格，见“航道内床面切应力比较.xlsx25wan外4段stress三角、矩形比较”。鉴于此，推断当临界淤积切应力大于床面切应力较多时，敏感性会低。6.8 航道回淤敏感性是否与临界淤积切应力相对床面切应力大小有关？答：有关。对25万吨级航道各段临界淤积切应力乘以1.3，求其对临界淤积切应力的敏感性，见huyan-8result-25wan- Sensitivity test“S