mike21c泥沙模块科学手册中文版

3 河流地貌学（河床变形模型）

河流地貌学模型是水动力学模型和泥沙输移模型的结合。根据河床地形的变化不断更新水流流场。

河流地貌学模型传统上分为耦合模型和非耦合模型两种。对于耦合模型，水流控制方程和泥沙输移方程合并到一个方程组中，同时求解。对于非耦合模型，水动力学流控制方程的求解在泥沙输移方程求解的前一个时间步里完成，然后计算出新的床面高程，并对下一时间步的水流方程进行求解。本模型采用的就是后一种方法（非耦合）。本模型包括的其他子模型，诸如河岸侵蚀，岸线更新，动床阻力，床面形态，床沙级配等，都将在本章中加以描述。

3.1 泥沙连续方程

根据计算的床沙质（推移质和悬移质）输移（状况），可以从以下方程求得床面高程的变化：

式中：Sx 为X方向输沙量 Sy 为Y方向输沙量 n 床沙孔隙率 t 时间

（x，y） 笛卡儿坐标系

ΔSe 由于河岸侵蚀引起的侧向泥沙补给。

总输沙量等于推移质和悬移质之和，对于曲线网格（s，n），方程3.1将略有不同，数值处理中在每一个网格单元的入流和出流边界上使用不同的Δs和Δn，见图3.2。

本模型采用空间中心差分和时间上的向前差分格式。时间步长受Courant准则的限制，即Courant数必须小于1。波数可以用式3.1的一维版本加以估算：

通过假定床

波的波速就可以表示为：

2

和泥沙输移是仅是希尔兹切应力θ的函数，那么沙

如果谢才系数假定为常数，那么θ与h2成反比，即& Hansen输沙公式，有：

，对Engelund

式中： S 输沙率 h 水深

cbw 粗估的沙波速

基于图3.2的符号系统，曲线网格的泥沙连续方程可以表达为如下差分格式：

式中：

Ss s方向的输沙率 Sn n方向的输沙率 n 床沙孔隙率 t 时间

（s，n） 曲线坐标系 Δs s方向的空间步长 Δn n方向的空间步长 （j，k） 网格坐标

为了系统的封闭，在上边界上需要提供边界条件。可选择两种（边界条件），指定床面高程变化率

或简单的指定进入系统的输沙率。

理论上讲，只需指定上边界泥沙输移条件。然而，由于模型允许模拟过程中水流方向的变化，所以在所有的边界上都要指定输沙率。每一指定的边界条件只有在入流情况下才是有效的。

3

泥沙输移边界

边界条件可以在每一边界点上指定，也可以在整个边界上指定。如果指定沿河宽的泥沙总量（可以是常数，也可以是时间序列），模型单位为m3/s。模型将自动按局部输沙能力的比例对总输沙量进行分配如下：

式中：

qs,i 边界线上i点的实际输沙率（m/s）(即等价于Ss或Sn) Qs 边界上总输沙率

qes,i边界线上i点的输沙能力（m2/s） n 边界线上的网格点数

Δxi 沿边界线上的网格点间距（即等于Δs或Δn）

输沙能力（也叫平衡挟沙能力，因为泥沙输移在局部水流条件下属于平衡输沙）通过选用的的泥沙输移公式计算出来。计算的输沙率可能会大于（或小于）挟沙能力，如果遇到这种情况，将在边界线上发生冲刷（或淤积）。

在指定泥沙输移边界条件时，对悬移质和推移质是不加区分的，边界条件直接应用到方程（3.1）和（3.2）。然而，当求解描述悬移质运动的对流扩散方程时，悬沙模型需要单独的含沙量边界条件。

含沙量边界条件通过推移质和全沙输沙量按如下方式计算自动得到：

2

式中：

q s,sus i点处实际悬移质输沙率 qs,bed i点处推移质输沙率 C 悬移质质量含沙量(g/m3) q i点流量（单宽流量）（m2/s）

4

ps 泥沙密度（kg/m3）

实际推移质与按局部水力条件计算得到的推移质输沙能力相同。只有悬沙由于时空的滞后可以超饱和或欠饱和输沙。这就意味作如果指定的边界条件超过全沙挟沙能力，多余的泥沙将作为额外的悬沙进入模型。

床面高程变化

如果边界条件是以dz/dt的形式指定，则方程（3.2）将不会求解，因为指定的边界条件可以直接计算边界处的新的床面高程。

也是在这种情况下，要求提供单独的质量含沙量，以求解反映悬移质运动的对流扩散方程。边界条件是按如下方式推导的：

式中：

C 悬沙质量含沙量(g/m3) Ce 悬沙平衡含沙量(g/m3)

dz/dt 指定的床面高程变化率(m/day) 86400 一天的秒数 q 点i的流量（单宽流量） Δy 边界线上网格点的间距（m）

上式是基于这样的假定：推移质输沙率总是局部平衡的，且泥沙保持连续性。给定的边界河床变化率将导致边界处的悬移质超饱和输移。作为泥沙输移边界条件，指定的值将直接带入泥沙连续方程，上式仅在对悬沙子模型指定边界条件时使用。

Flood and Dry Points（淹没和干涸点）

在数学模型中那些水深较浅，时常遭遇淹没和干涸的点需要特殊的处理。 在水动力学模型中，淹没水深和干涸水深被应用以处理水深变化所导致的洪泛区的不同蔓延这样的水力条件：

干涸水深：如果局部水深比它小，则此点作为陆地处理。 淹没水深：如果局部水深比它大，则此点作为水域处理。

为了确保水动力学模型的数值稳定性，有必要为淹没水深和干涸水深指定一

5

个确定的数值（通常在0.2m～0.6m）。而且，必须使淹没水深和干涸水深的差值不大，以避免振荡交替所带来的不稳定性。

浅水中的地貌变化特点与深水中情况略有不同。比如排水渠道的形态就比主槽中小，因此要求模型有更高的分辨率。然而这些对总体形态影响较小而没有包括。

采用的方法是干单元不产生泥沙输移，除非是湿单元，否则床面高程不改变。仅有的限制为开边界线上的网格点必须保持要么是干单元要么是湿单元。然而，如果采用准恒定流求解方法来求解水动力学模型，模型的也允许开边界单元在模拟过程中有干湿变化。

Availability of Sediment 泥沙补给量

本模型可以应用到泥沙可用量减少的情况。例如当冲刷到基岩时。 对于沙质河床，没有必要使用这一工具。如要使用，需指定河床上泥沙的初始量。

3.2 河岸侵蚀

泥沙输移模型的连续方程可以包括河岸侵蚀：

Eb 河岸侵蚀率(m/s), z 局部床面高程 S 近岸输沙率 h 局部水深

αβγ 模型中指定的率定参数

当河岸侵蚀时释放的额外泥沙加入泥沙平衡方程（3.1），其贡献可从以下关系得到：

式中：

hb 水面以上河岸高度 Δn 靠岸线单元宽度

6

式（3.10）第一项包括了近岸坡度有关的影响，第二项间接包含了水流对河岸施加的切应力。最后一项，代表固定侵蚀率，与水流条件无关。

第一项假定横断面形态不变。因此，从几何上考虑河岸侵蚀与河床侵蚀成正比。系数α与河岸处河床横向坡度相应。天然河道中沿河岸线河床的横向坡度的典型量级为5～20。这就意味作系数α应该在0～20的范围内指定。

第2项是基于如下考虑：侵蚀量被水流带向下游确定了输沙率，或者由比例为β的侵蚀量所为。这就意味作β应在0～1的范围内选择。注意河岸侵蚀模块假定河岸物质与河床上的泥沙一样。

观测到河岸侵蚀率被用来率定系数和α和β。如果没有发现t与S的明显的相关关系，可以把河岸侵蚀率γ指定为常数。

每一时间步都要计算河岸侵蚀，河岸冲刷量将在泥沙连续性方程中包括，见式（3.1）。如果某一部位河岸累积侵蚀量超过了预先定义的临界宽度（与沿河岸的计算网格宽度有关），曲线网格将用新的河岸来根新。这就意味作模型中的如下步骤：

由新的河岸线产生新的正交网格。

从新计算网格参数，比如空间步长?s，?n，和网格线曲率半径Rs和Rn。 将模型地形从老网格转换到新网格。

初始化河岸累积侵蚀记录（累积侵蚀记录在网格更新时被用来作检验）。

7

如果河岸侵蚀只影响侵蚀河岸附近的地貌发展而不是整个河道，则平面更新将排除在外。

3.3 线性稳定性分析（本小节翻译中问题较多）

为了理解地貌模型输出的物理背景，对水流和泥沙输移方程进行稳定性分析是有用的。下面将作简要介绍。更纤细的信息参考Olesen(1987)的文献。方法如下：

与二维流有关的长度量级（length scale）的推导。 与河床地貌学有关的长度量级（length scale）的推导。 为了得到长度量级，假定了恒定流条件。

最后，两种方法结合起来，从线性化的方程中推导出河流的整个地貌特性。 纵向的动量方程可以写成：

方程（3.12）已经线性化，所有因变量都已经分界成平均值与脉动值之和：

将式（3.13）～（3.16）带入式（3.12）得到：

忽略超过一阶的高阶项，并假定v0=0（按水流沿流线运动考虑），式（3.17）变为：

式（3.18）第5项可以如下线性化：

8

式（3.19）带入式（3.18）并作时间平均，零阶项之和等于零而何以去掉。这样得到如下脉动项的线性微分方程：

假定压力可以忽略，这就意味推导出方程仅在水流惯性和河道阻力控制的摩擦控制的水流（friction dominated flow）中才是合理的。由式（3.20）得：

式中长度比尺λw可以由下式计算：

参数λw是长度比尺的度量为流场适应河床高程的变化，参见示意图4.3。

长度比尺在不同的方面都很重要。如果一个心滩在发育中，水流的惯性趋向于冲平河床，而摩擦力会增强沙洲的进一步生长。在数值项中，浅水点的流速既能由水深减小而单宽流量相同而增加，又能由于沙洲的分流而减小。

不难预料，自适应长度与局部水深成比例。由谢才系数代表的河床阻力也会

9

影响到长度比尺。假定沙洲上的床面阻力较高，则长度比尺小，水流易于偏离。相反（若阻力小），则自适应长度较大（或水流偏离的强度会减弱）。

床面的形成不仅受控于水力条件，而且受控于泥沙输移和泥沙连续方程：

方程（3.23）是建立在曲线坐标系中流线s方向上的输沙率公式。横向的输沙率Sn 由螺旋流强度和横向床面坡度来确定，见第2节式（2.67）：

近似的，流线上的泥沙输移假定依赖于下式：

式中比例因子a1和指数b都假定为常数。式（3.25）对曲线坐标s求偏导数，得到：

如果式（3.24）中希尔兹数θ可以忽略，结合式（3.23），（3.24）和（3.26）并加以线性化，可得下式：

流体连续方程有：

用线性化格式公式（3.28）可写成：

结合式（3.27）和（3.29），可推导以下共轭微分方程：

除了微分方程式（3.30）外，还有边界条件。对于过上游边界水深为正弦变化的情况，解可以写成如下形式：

10

式中：i 虚数 w 河宽

n 横向坐标（单位向量） s 纵向坐标（单位向量） 带入式（3.30），得：

与水动力学得耦合通过定义以下参数来实现：

而床面扰动的长度变换的长度比尺定义为：

式（3.32）的解为：

一个求解的图例见图3.5。

自适应长度λs预示作在通过泥沙输移夷平之前，下游一段距离的任何扰动都将被感受到。横向坡度因子G，结合螺旋流强度，控制床面切应力的偏向（进

11

而控制泥沙输移的方向），看来对自适应长度比尺相当重要。理论上讲，如果床面切应力不受河床横向坡度的影响，长度比尺将变成无穷大。不然，任何干扰被保留。如果G值较大，横向坡度的重力效应将对床面切应力产生强烈的偏向，自适应长度将变小，干扰将迅速消失。

从式（3.34）可以看到，宽深比河单独的水深对长度比尺都很重要。事实上，求解的边界条件确定了宽深比的重要性。

如果将式（3.33）带入作为对干扰在在流场中对河流地貌产生相映的式（3.32），那么对于坐标s，我们有：

作为河床地貌变化对流场的相映的式（3.21），也可以对坐标s写成：

将

项从式（3.32）中分离出来，

项从式（3.36）中分离出来，带入式（3.37），

可以得到：

式（3.38）的复数解为：

式中：

j 解的个数，对于二阶微分方程有两个根存在。 复波数可以分解成代表波长的实部，和代表波阻的虚部：

把式（3.39）带入式（3.37）得到确定波数的多项式：

复行列式为：

12

式（3.41）有复数解：

对于

，波数的实部为：

而波数的虚部为：

求解草图见图3.6。

对于

大于1时，解变得不稳定，因为阻尼变成负数，线性化方程时采用

的假定不再正确。

3.4 动床阻力

有两种方式可以更新谢才系数和满宁系数，要么为水深的简单函数，要么为计算的沙波尺度的函数。

由水深更新阻力系数可表达如下：

式中：A 阻力系数

13

b 阻力指数

可变阻力的应用将影响冲淤模式。采用了阻力更新，水流将更加偏离浅水区，从而造成床面切应力加大，泥沙输移增强。通常，这将导致河床冲刷的超条效应，沙洲的地貌将变得更加圆滑，正如1992年Talmon所述。

3.5 补给受限的泥沙输移

补给受限的泥沙输移意味作，局部的挟沙力由两个因数的乘积确定，一个是局部床沙层厚度与平衡床沙厚度的相对比值，另一个是由选定的公式计算出的输沙能力。补给限制的泥沙输移是在有一层床沙时激活，而不是缺省的0层。输沙公式修改为：

式中对输沙能力

进行了缩放得到输沙率S。缩放是通过函数f，以及床沙

时，输沙

厚度?相对于指定的可以充足补给泥沙的床沙厚度?eq的比值。当率小于挟沙力。函数f变为

。当床沙厚度大时输沙率不受此影响。

平衡厚度作为水深比例来计算的，缺省为h/12，因为典型的沙波高度为水深的1/6。可以指定平衡厚度相对于水深的最大和最小值。如果最大和最小值相同，也可以指定厚度为常数。注意，很薄的床沙会很快移动，将要求较小的地貌（河床变形）时间步长。

补给受限的泥沙输移模型可以用来模拟薄层泥沙，比如下水道、粗化河床表面、局部抗冲区。补给受限模型将跟踪局部床沙厚度，能够再生易于冲刷的薄层床沙。补给受限的泥沙输移模型有两个主要的应用：

局部抗冲区域

例如河流中被植被覆盖的小岛，就可以用补给受限的泥沙输移模型。当指定初始床沙厚度的分布时，岛上可以指定为0，而其他地方用一个较大数值，比如100m。这就意味作岛是不可冲的，但多余的泥沙可以在岛上沉积（以后也可以被冲走）。

对地貌系统的补给比输沙能力低（不饱和输沙）

这是补给受限的泥沙输移模型最重要的应用。这里泥沙将在粗化河床（或混凝土管道）上形成厚度可变的床沙层。床沙厚度由进入的泥沙补给确定。这对

14

于粗化河道很重要，也可以用来模拟沙洲。

15

本文来源：https://www.bwwdw.com/article/j5vp.html