43、一维非恒定变形床模型：原理、应用与验证

一维非恒定变形床模型：原理、应用与验证

1. 一维非恒定泥沙输移模型分类

一维非恒定泥沙输移模型主要分为两类：
- 非耦合模型 ：在给定时间步长内，水流方程和泥沙连续性方程不耦合。
- 准恒定流模型 ：将能量方程与泥沙连续性方程联立求解。

Lyn（1987）使用摄动技术识别控制方程中的多个时间尺度，建议在边界条件快速变化的情况下，将完全非恒定水流方程与泥沙连续性方程完全耦合，他提出使用Preissmann线性隐式格式联立求解控制方程。Park和Jain（1986）在其非恒定、非耦合模型中使用该方法分析泥沙过载引起的淤积，但为避免不稳定，当河床高程变化的空间梯度过大时，需迭代求解。

2. 一维非恒定耦合变形床模型

这里介绍一个一维非恒定耦合变形床模型（Bhallamudi和Chaudhry，1991），该模型使用MacCormack显式格式（MacCormack，1969）联立求解完整的圣维南水流方程和泥沙连续性方程，并将计算结果与实验数据对比验证模型。

2.1 控制方程

描述宽矩形可变形河床非恒定流的方程如下：
- 水流连续性方程 ：
[
\frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=0 \quad(18 - 36)
]
- 水流动量方程 ：
[
\frac{\partial q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{q^{2}}{h}+\frac{1}{2} g h^{2}\right)+g h \frac{\partial z}{\partial x}+g h S_{f}=0 \quad(18 - 37)
]
- 泥沙连续性方程 ：
[
\frac{\partial}{\partial t}\left((1 - p) z+\frac{q_{s} h}{q}\right)+\frac{\partial q_{s}}{\partial x}=0 \quad(18 - 38)
]
其中，$q$ 为单位宽度的水流量；$h$ 为水流深度；$z$ 为河床底部高程；$q_{s}$ 为单位泥沙流量；$p$ 为河床层的孔隙率。

泥沙流量可通过水流速度的经验幂函数估算：
[
q_{s}=a\left(\frac{q}{h}\right)^{b} \quad(18 - 39)
]
其中，$a$ 和 $b$ 是经验常数，其值取决于泥沙特性。此泥沙流量关系主要为简化模型，由于控制方程采用显式数值格式求解，可轻松包含更复杂的关系。

实验研究表明，淤积河道非均匀条件下的曼宁系数 $n$ 小于均匀流时的值。在非均匀条件下，若使用局部摩擦坡度而非底部坡度，可使用均匀流的阻力方程。但恒定 - 均匀泥沙输移公式在非均匀条件下的适用性有待进一步研究。

2.2 数值格式

方程（18 - 36）至（18 - 38）是一组非线性双曲型方程，仅在理想情况下有解析解，这里使用MacCormack格式进行数值求解。该格式易于实现，无需特殊处理即可捕捉激波，且便于纳入粗糙度和泥沙流量的一般经验方程。

使用前向有限差分计算预测步的空间偏导数，使用后向有限差分计算校正步的空间偏导数，且该顺序可每两个时间步反转一次。

• 预测步 ：
  [
  h_{i}^{ }=h_{i}^{k}-\frac{\Delta t}{\Delta x}\left(q_{i + 1}^{k}-q_{i}^{k}\right)
  ]
  [
  q_{i}^{ }=q_{i}^{k}-\frac{\Delta t}{\Delta x}\left(\frac{\left(q_{i + 1}^{k}\right)^{2}}{h_{i + 1}^{k}}-\frac{\left(q_{i}^{k}\right)^{2}}{h_{i}^{k}}+\frac{g}{2}\left(\left(h_{i + 1}^{k}\right)^{2}-\left(h_{i}^{k}\right)^{2}\right)\right)-g h_{i}^{k} \frac{\Delta t}{\Delta x}\left(z_{i + 1}^{k}-z_{i}^{k}\right)-g h_{i}^{k} \Delta t \frac{\left(q_{i}^{k} n\right)^{2}}{\left(h_{i}^{k}\right)^{3.33}}
  ]
  [
  z_{i}^{ }=z_{i}^{k}+\frac{1}{1 - p}\left(\left(\frac{q_{s} h}{q}\right) {i}^{k}-\left(\frac{q {s} h}{q}\right)_{i}^{ }\right)-\frac{\Delta t}{(1 - p) \Delta x}\left((q_{s}) {i + 1}^{k}-(q {s}) {i}^{k}\right)
  ]
  [
  (q {s}) {i}^{ }=a\left(\frac{q_{i}^{ }}{h {i}^{*}}\right)^{b} \quad(18 - 40)
  ]
• 校正步 ：
  [
  h_{i}^{ }=h_{i}^{ }-\frac{\Delta t}{\Delta x}\left(q_{i}^{ }-q_{i - 1}^{*}\right)
  ]
  [
  q_{i}^{ }=q_{i}^{ }-\frac{\Delta t}{\Delta x}\left(\frac{\left(q_{i}^{ }\right)^{2}}{h_{i}^{ }}-\frac{\left(q_{i - 1}^{ }\right)^{2}}{h_{i - 1}^{ }}+\frac{g}{2}\left(\left(h_{i}^{ }\right)^{2}-\left(h_{i - 1}^{ }\right)^{2}\right)\right)-g h_{i}^{ } \frac{\Delta t}{\Delta x}\left(z_{i}^{ }-z_{i - 1}^{ }\right)-g h_{i}^{ } \Delta t \frac{\left(q_{i}^{ } n\right)^{2}}{\left(h_{i}^{ }\right)^{3.33}}
  ]
  [
  z_{i}^{ }=z_{i}^{ }+\frac{1}{1 - p}\left(\left(\frac{q_{s} h}{q}\right) {i}^{ }-\left(\frac{q_{s} h}{q}\right)_{i}^{ }\right)-\frac{\Delta t}{(1 - p) \Delta x}\left((q_{s})_{i}^{ }-(q {s})_{i - 1}^{ }\right)
  ]
  [
  (q_{s})_{i}^{ }=a\left(\frac{q_{i}^{ }}{h_{i}^{ *}}\right)^{b} \quad(18 - 41)
  ]

时间层 $k + 1$ 的未知量值由下式给出：
[
h_{i}^{k + 1}=\frac{1}{2}\left(h_{i}^{k}+h_{i}^{ }\right)
]
[
q_{i}^{k + 1}=\frac{1}{2}\left(q_{i}^{k}+q_{i}^{ }\right)
]
[
z_{i}^{k + 1}=\frac{1}{2}\left(z_{i}^{k}+z_{i}^{**}\right) \quad(18 - 42)
]

边界节点的因变量值通过边界条件确定。对于亚临界流，根据特征理论，上游边界需指定两个边界条件，下游边界需指定一个边界条件。未通过边界条件指定的因变量值可由特征方程确定，也可通过内节点的已知值插值得到。

2.3 稳定性与人工粘性

• 稳定性 ：MacCormack格式需满足Courant - Friedrichs - Lewy（CFL）条件以保证稳定性。由于水波传播速度远高于河床瞬变速度，该条件为：
  [
  C_{n}=\left(\frac{q}{h}+\sqrt{g h}\right) \frac{\Delta t}{\Delta x} \leq 1 \quad(18 - 43)
  ]
  其中，$C_{n}$ 为Courant数。为保证格式稳定，每个网格点都需满足该条件。
• 人工粘性 ：可通过引入人工粘性来抑制陡峭波前附近的数值振荡，可使用Jameson方法。在研究的几个案例中，仅在节点迁移案例中需要进行平滑处理。

2.4 计算步骤

假设已知时间层 $k$ 时，将河道划分为 $N$ 个河段的所有网格点（$i = 1, \cdots, N + 1$）上的 $h_{i}^{k}$、$q_{i}^{k}$ 和 $z_{i}^{k}$ 值，要确定未知时间层 $k + 1$ 的值，步骤如下：
1. 使用方程（18 - 40）计算内节点（$i = 2, \cdots, N$）的 $h_{i}^{ }$、$q_{i}^{ }$ 和 $z_{i}^{ }$ 值，使用适当的边界条件计算边界节点（$i = 1$ 和 $i = N + 1$）的值。
2. 使用方程（18 - 41）确定内节点（$i = 2, \cdots, N$）的 $h_{i}^{ }$、$q_{i}^{ }$ 和 $z_{i}^{ *}$ 值，从边界条件确定边界节点的值。
3. 使用方程（18 - 42）确定时间间隔 $\Delta t$ 结束时的 $h$、$q$ 和 $z$ 值，即 $h_{i}^{k + 1}$、$q_{i}^{k + 1}$ 和 $z_{i}^{k + 1}$。
4. 如有必要，使用Jameson人工粘性方法修改步骤3中确定的值，以抑制高频振荡。
5. 将下一个时间间隔的 $h_{i}^{k}$、$q_{i}^{k}$ 和 $z_{i}^{k}$ 值设置为 $h_{i}^{k + 1}$、$q_{i}^{k + 1}$ 和 $z_{i}^{k + 1}$；根据方程（18 - 43）确定下一个时间步的时间间隔 $\Delta t$；重复该过程，直到完成指定时间的模拟。

该过程中水流方程和泥沙连续性方程是耦合的，因为采用了两层预测 - 校正方法。严格来说，预测步没有耦合，但预测的河床高程值用于确定方程（18 - 37）中流量的正确值，预测的 $h$ 和 $q$ 值用于确定 $q_{s}$ 并计算方程（18 - 38）中的空间导数项。因此，每个时间步结束时计算的因变量值考虑了所有其他变量的变化。

3. 模型应用与验证

为说明该模型的应用，将计算结果与两个案例的实验结果进行比较。

3.1 泥沙过载引起的淤积

实验室渠道淤积过程的测试结果（Soni等，1980）与计算结果进行对比。渠道宽0.2 m，长30 m，河床砂和注入材料的平均直径为0.32 mm。根据均匀流测量，计算得到：$a = 1.45×10^{-3}$；$b = 5.0$；曼宁系数 $n = 0.022$；泥沙床层孔隙率 $p = 0.4$。初始水流量 $q_{0}=0.020 m^{2}/s$；均匀流深度 $h_{0}=0.05 m$；初始河床坡度 $S_{0}=0.00356$；平衡泥沙流量 $q_{s0}=4$。

在数学模型中，将均匀单位流量、均匀流深度和根据 $S_{0}$ 计算的初始河床高程指定为每个节点的初始条件。通过在上游端增加泥沙流量速率 $\Delta q_{s}$ 启动瞬态状态。指定下游端一个边界条件和上游端两个边界条件。上游边界代表恒定流量，通过指定 $q(0, t)=q_{0}$（$t\geq0$）实现。第二个边界 $q_{s}(0, t)=q_{s0}+\Delta q_{s}$ 的处理不那么直接，需转化为方程以计算上游端的河床高程，通过假设节点1上游的虚拟节点并指定该节点的泥沙流量等于 $q_{s0}+\Delta q_{s}$，使用泥沙连续性方程并对空间微分项应用后向有限差分，得到：
[
\left((1 - p) z+\frac{q_{s} h}{q}\right) {1}^{k + 1}=\left((1 - p) z+\frac{q {s} h}{q}\right) {1}^{k}+\frac{\Delta t}{\Delta x}\left(\left(q {s0}+\Delta q_{s}\right)-\left(q_{s}\right) {1}^{k}\right) \quad(18 - 44)
]
由于时间层 $k$ 右侧各项已知，可计算未知时间层 $k + 1$ 方程（18 - 44）左侧的值，进而得到 $z$。边界处的水流深度由特征方程和已知流量确定。下游边界指定为恒定深度 $h(N\Delta x, t)=h {0}$（$t\geq0$），下游端的流量和河床高程通过内节点值外推确定。

数学模型应用于50 m长的渠道，分为50个河段（$\Delta x = 1.0 m$），选择计算时间步 $\Delta t$ 以满足稳定性的Courant条件（$C_{n}=0.9$）。图18 - 6比较了 $t = 40$ 分钟时计算的瞬态河床和水面轮廓与测量值，数学模型能较好地模拟渠道河床的淤积和瞬态水面轮廓。

graph TD;
    A[开始] --> B[设置初始条件];
    B --> C[计算预测步];
    C --> D[计算校正步];
    D --> E[更新变量值];
    E --> F[检查稳定性条件];
    F -- 满足 --> G[是否完成指定时间];
    G -- 否 --> C;
    G -- 是 --> H[结束];
    F -- 不满足 --> I[调整时间步长];
    I --> C;

3.2 节点迁移

节点定义为河道纵向底部轮廓的突然变化（图18 - 7）。在非侵蚀性河床的河道中，节点保持不变；在侵蚀性河床的河道中，节点向上游迁移并逐渐消失。Brush和Wolman（1960）解释节点迁移是由于坡度变化处的侵蚀潜力（$hS_{f}$）达到最大。节点上游为亚临界流，下游为超临界流，底部坡度突变处的边界剪切力（$\tau_{0}=\gamma hS_{f}$）最大，下游随水流深度减小而减小，上游虽水流深度增加但能量梯度线变平也减小。由于泥沙输移与剪切应力成正比，节点处较高的剪切应力导致该位置带走的物质比上下游河段多，因此节点向上游迁移，侵蚀物质在下游沉积，最终过陡河段变平。

Brush和Wolman（1960）在15.8 m长、1.2 m宽的水槽中获得的实验数据与计算结果进行对比。实验前，在中值粒径为0.67 mm的非粘性砂中塑造了一个宽0.21 m、带圆角的梯形渠道，在距上游端10.8 m处设置约0.0305 m的落差模拟过陡河段或节点，该点上下游渠道坡度约为0.00125。然后开启水流（$q_{0}=0.0028 m^{2}/s$，$h_{0}=0.0305 m$），并记录不同时间的河床高程。

为模拟该实验，将渠道分为52个河段（$\Delta x = 0.3048 m$），初始和边界条件与前一案例相同，Courant数为0.85，使用Jameson方法添加人工粘性。泥沙流量由方程（18 - 39）确定，$b = 4.2$，系数 $a$ 假设为能量梯度线 $S_{e}$ 的函数。每个新时间步开始时，使用后向有限差分和已知时间层的 $q$、$h$ 和 $z$ 值计算每个网格点的能量梯度线坡度：
[
S_{e}=\frac{1}{\Delta x}\left(z_{i - 1}+h_{i - 1}+\frac{q_{i - 1}^{2}}{2 g h_{i - 1}^{2}}\right)-\frac{1}{\Delta x}\left(z_{i}+h_{i}+\frac{q_{i}^{2}}{2 g h_{i}^{2}}\right) \quad(18 - 45)
]
然后由 $a = S_{e}^{1.71}$ 确定 $a$ 的值。方程（18 - 39）中的指数根据泥沙输移测量值（Brush和Wolman，1960）估算。除渠道坡度外，两个实验条件相同，泥沙输移速率的差异可通过幂律与能量梯度线坡度相关。

图18 - 8比较了 $t = 2.67$ 小时时的计算结果和测量结果，尽管泥沙输移方程存在固有不确定性，但计算的河床轮廓与测量的河床轮廓匹配良好。实验中节点下游的河床高程高于模型预测值，但上游侵蚀量与下游沉积量几乎相等。实验中观察到节点下游河道展宽，实验中较大的沉积量可能是由于河岸侵蚀提供了额外的泥沙。

综上所述，该一维非恒定耦合变形床模型在模拟渠道淤积和节点迁移等现象方面具有较好的效果，能够为实际工程中的河道演变分析提供有力的工具。通过与实验数据的对比验证，模型的可靠性得到了一定程度的证明，但在泥沙输移方程等方面仍存在改进空间。

一维非恒定变形床模型：原理、应用与验证

4. 模型优势与局限性分析

4.1 模型优势

• 耦合性优势 ：该模型采用两层预测 - 校正方法实现了水流方程和泥沙连续性方程的耦合。在每个时间步结束时，计算的因变量值能考虑到所有其他变量的变化，相较于非耦合模型，能更准确地模拟河道中水流与泥沙相互作用的复杂过程。例如在泥沙过载引起的淤积模拟中，能更合理地反映出河床高程变化对水流流量以及泥沙输移的影响。
• 数值格式优势 ：使用MacCormack显式格式进行数值求解，此格式易于实现，无需特殊处理即可捕捉激波，并且便于纳入粗糙度和泥沙流量的一般经验方程。这使得模型在处理复杂的水流和泥沙问题时具有较高的灵活性，能够适应不同的实际情况。
• 验证效果较好 ：通过与两个案例（泥沙过载引起的淤积和节点迁移）的实验结果进行对比，计算结果与实验数据匹配良好，说明该模型在模拟河道淤积和节点迁移等现象方面具有较好的效果，能够为实际工程中的河道演变分析提供有力的工具。

4.2 模型局限性

• 泥沙输移方程的不确定性 ：虽然模型中使用的泥沙流量经验幂函数关系简单，但存在一定的局限性。实际的泥沙输移过程非常复杂，受到多种因素的影响，如泥沙颗粒的形状、水流的紊动特性等。因此，即使通过实验数据估算了经验常数，泥沙输移方程仍存在固有不确定性，可能导致模拟结果与实际情况存在一定偏差。例如在节点迁移案例中，实验中节点下游的河床高程高于模型预测值，可能与泥沙输移方程的不完善有关。
• 边界条件处理的复杂性 ：对于亚临界流，需要根据特征理论在上游边界指定两个边界条件，下游边界指定一个边界条件。并且部分边界条件的处理不那么直接，需要进行转化和额外的计算，如上游边界的第二个边界 $q_{s}(0, t)=q_{s0}+\Delta q_{s}$ 需要转化为方程以计算上游端的河床高程。这增加了模型应用的复杂性，并且边界条件的准确性对模拟结果的影响较大。
• 未考虑一些实际因素 ：模型在模拟过程中没有考虑到一些实际因素，如河道的侧向侵蚀、河岸的坍塌等。在实际的河道演变过程中，这些因素可能会对河床的变化产生重要影响。例如在节点迁移案例中，实验中观察到节点下游河道展宽，而模型未考虑这一因素，可能导致实验中较大的沉积量无法被准确模拟。

5. 模型的改进方向与展望

5.1 改进泥沙输移方程

• 引入更复杂的物理模型 ：可以考虑引入基于物理过程的泥沙输移模型，如考虑泥沙颗粒的受力分析、水流的紊动扩散等因素，以提高泥沙输移方程的准确性。例如，采用两相流模型来描述水流和泥沙的相互作用，能够更详细地模拟泥沙的运动过程。
• 结合机器学习方法 ：利用机器学习算法对大量的实验数据和实际观测数据进行分析，建立更准确的泥沙输移预测模型。通过训练模型，使其能够自动学习泥沙输移与各种影响因素之间的复杂关系，从而提高模拟结果的精度。

5.2 优化边界条件处理

• 开发通用的边界条件处理方法 ：针对不同类型的边界条件，开发一套通用的处理方法，减少边界条件处理的复杂性。例如，建立边界条件数据库，根据不同的河道情况和水流条件，自动选择合适的边界条件处理方式。
• 提高边界条件的准确性 ：通过更多的现场观测和实验数据，对边界条件进行更准确的确定。例如，在确定上游泥沙流量边界条件时，可以结合实际的泥沙来源和输移情况进行精确估算。

5.3 考虑更多实际因素

• 纳入侧向侵蚀和河岸坍塌模型 ：在模型中加入侧向侵蚀和河岸坍塌的模拟模块，考虑河道的侧向变形对水流和泥沙输移的影响。可以采用经验公式或物理模型来描述侧向侵蚀和河岸坍塌的过程，提高模型对实际河道演变的模拟能力。
• 考虑气候变化和人类活动的影响 ：气候变化和人类活动（如水利工程建设、土地利用变化等）会对河道的水流和泥沙条件产生重要影响。在模型中考虑这些因素，能够更全面地模拟河道的长期演变过程，为水资源管理和河道治理提供更准确的决策依据。

6. 实际工程应用建议

6.1 前期数据收集与分析

• 地形数据 ：收集详细的河道地形数据，包括河床高程、河道宽度、坡度等信息，为模型的初始化提供准确的基础数据。可以采用地形测量、遥感等技术获取高精度的地形数据。
• 水流和泥沙数据 ：测量河道的水流流量、流速、泥沙含量等数据，确定模型中的参数，如曼宁系数、泥沙输移经验常数等。同时，收集不同季节和工况下的水流和泥沙数据，以反映实际情况的变化。
• 边界条件数据 ：确定合适的边界条件，如上游的水流流量和泥沙流量边界条件、下游的水位或流量边界条件。通过现场观测和历史数据统计，准确确定边界条件的取值。

6.2 模型参数校准与验证

• 参数校准 ：在模型应用前，利用收集到的实际数据对模型中的参数进行校准。通过调整参数值，使模型的计算结果与实际观测数据尽可能匹配。可以采用试错法、优化算法等方法进行参数校准。
• 模型验证 ：使用独立的实验数据或实际观测数据对校准后的模型进行验证，评估模型的准确性和可靠性。如果模型的计算结果与验证数据存在较大偏差，需要进一步分析原因并调整模型。

6.3 模拟结果的应用与决策支持

• 河道演变预测 ：利用校准和验证后的模型对河道的未来演变进行预测，分析不同工况下河道的淤积和侵蚀情况，为河道治理和水资源管理提供科学依据。例如，预测在不同的水流流量和泥沙输入条件下，河道的河床高程和水面线变化。
• 工程方案评估 ：对不同的水利工程方案（如修建水库、河道整治等）进行模拟评估，比较不同方案对河道水流和泥沙输移的影响，选择最优的工程方案。例如，评估水库建设对下游河道淤积和生态环境的影响。

7. 总结

本研究介绍了一维非恒定耦合变形床模型，详细阐述了其原理、数值格式、计算步骤以及应用案例。该模型在模拟渠道淤积和节点迁移等现象方面具有较好的效果，通过与实验数据的对比验证，证明了其可靠性。然而，模型也存在一些局限性，如泥沙输移方程的不确定性、边界条件处理的复杂性以及未考虑一些实际因素等。为了提高模型的准确性和适用性，提出了改进泥沙输移方程、优化边界条件处理和考虑更多实际因素等改进方向。在实际工程应用中，建议进行充分的前期数据收集与分析、模型参数校准与验证，并将模拟结果应用于河道演变预测和工程方案评估等方面。

以下是一个总结模型应用流程的表格：
|步骤|具体内容|
| ---- | ---- |
|前期准备|收集地形、水流、泥沙和边界条件数据|
|模型设置|确定模型参数，设置初始条件和边界条件|
|模型计算|按照计算步骤进行迭代计算，更新变量值|
|结果分析|对比计算结果与实验数据，评估模型准确性|
|应用决策|利用模拟结果进行河道演变预测和工程方案评估|

graph LR;
    A[前期准备] --> B[模型设置];
    B --> C[模型计算];
    C --> D[结果分析];
    D -- 不满意 --> B;
    D -- 满意 --> E[应用决策];

通过不断地改进和完善模型，以及合理地应用模型结果，该一维非恒定耦合变形床模型将在河道演变分析和水利工程设计中发挥更大的作用。