37、堤防决口建模：原理、方法与数值模拟

堤防决口建模：原理、方法与数值模拟

1. 引言

堤防通常沿河流、溪流和渠道修建，用于防洪保护，或围绕蓄水设施建造，一般采用可侵蚀材料。这些土堤结构可能因漫顶、管涌或边坡失稳而失效，其中漫顶是最常见的失效原因。堤防决口可能导致灾难性的财产损失和人员伤亡，如2005年卡特里娜飓风后新奥尔良市的破坏、2012年桑迪飓风对美国东海岸的影响，以及2015年华金飓风对南卡罗来纳州哥伦比亚市的破坏。

准确预测决口的发展、决口出流以及由此产生的流场，对于应急准备和减轻洪水灾害至关重要。堤防决口过程复杂，涉及水流、泥沙输移和相应的地貌变化之间的相互作用。

堤防决口可分为意外决口和有计划的工程决口。对于准确的洪水模拟，需要知道决口的形状和大小。工程决口的这些参数是预先已知的，但对于意外决口，决口的形状、大小和随时间的发展必须与水流计算同时确定，因为它们相互依赖。目前大多数洪水模拟模型需要用户指定决口大小作为输入，而适当的动态建模应将决口发展视为一个时间演化过程，取决于决口内外的水流、堤防材料和堤防断面。

2. 决口流量估算

侧向出流通过堤防决口可以通过数值求解一维空间变化流方程或二维浅水方程来估算。与一维模型相比，二维模型需要更多的地形信息、更多的离散化工作，并且需要将决口作为内部边界包含在内。而一维模型通常更易于应用。

以下是一些不同学者开发的模型：
- Han等人（1998） ：开发了一个结合一维和二维的水动力模型，一维模型求解主河道的动态波方程，二维模型求解漫滩的扩散波方程，并应用于汉江下游的实际堤防决口案例。
- Roger等人（2009） ：开发了两个堤防决口水流模型，分别使用总变差递减的Runge - Kutta间断Galerkin有限元方法和有限体积法求解浅水方程。
- Van Emelen等人（2012） ：数值求解二维浅水方程，其计算结果与新奥尔良第17街运河决口的1:50比例模型测量结果吻合良好。

Cheong（1991）利用能量和动量方程以及实验结果估算棱柱形梯形渠道的侧向出流，建议不要使用下游水深来估算侧堰的流量系数。Singh等人（1994）的实验表明，矩形侧堰的流量系数是上游弗劳德数和堰顶高度与上游水深之比的函数。

Elafy等人（2018）将堤防决口的侧向出流视为矩形和梯形主河道中宽顶侧堰的侧向出流，求解空间变化流方程，并引入了修正系数。

一维空间变化流方程（有侧向出流）可写为：
[
\frac{dh}{dx} = \frac{S_0 - S_f - \frac{Qq_b}{gA^2}}{1 - \frac{Q^2}{gA^2D}}
]
其中，(\frac{dh}{dx}) 是沿决口宽度的水流深度梯度，(S_0) 是渠道底部坡度，(S_f) 是能量梯度线坡度，(Q) 是主河道流量，(q_b) 是决口单位宽度的侧向出流强度，(A) 和 (D) 分别是决口处的过流面积和水力平均深度。对于水平棱柱形渠道上的宽顶侧堰，(q_b) 可根据Hager（1987）的数据计算。

为了修正一维模型计算的堤防决口流量，Elafy等人（2018）引入了修正系数 (C_f)：
[
C_f = \frac{Q_{bs} - Q_{b2}}{Q_i}
]
其中，(Q_{bs}) 是使用空间变化流模型计算的决口流量，(Q_{b2}) 是使用二维模型计算的决口流量。通过比较一维和二维模型的结果并进行多元回归分析，得到不同主河道形状下 (C_f) 的表达式：
- 矩形主河道 ：
[
C_f = 0.48S_r^{6.2}b_r^{0.72}F_{app}^{-0.15}F_{out}^{-0.13}
]
- 梯形主河道 ：
[
C_f = 0.3S_r^{0.8}b_r^{0.72}F_{app}^{-0.15}F_{out}^{-0.13}
]
其中，(S_r = \frac{h_d - h_b}{h_u})，(h_u) 和 (h_d) 分别是主河道决口上下游的水深，(h_b) 是决口顶部高于主河道底部的高度，(b_r = \frac{b_b}{b_c}) 是决口宽度与河道底宽之比。(S_r) 的范围是0.71 - 0.90，(b_r) 的范围是0.12 - 0.86，(F_{app})（决口上游主河道的弗劳德数）的范围是0.46 - 0.9，(F_{out})（模型出口的弗劳德数）的范围是0.01 - 0.62。

从这些方程可以看出，修正系数受淹没比的影响比相对决口宽度更大，并且与接近弗劳德数和出口弗劳德数呈反比关系。

3. 漫顶导致的堤防决口

对于非粘性土堤的漫顶破坏，不同学者进行了相关研究：
- Faeh（2007） ：开发了一个二维数值模型，使用有限体积法求解浅水方程和泥沙质量守恒方程，应用于基于侵蚀的大坝决口和易北河堤防决口案例，指出决口边坡角度是影响侧向侵蚀的最敏感参数。
- Kakinuma和Shimizu（2014） ：进行了四次大规模的河道堤防决口实验，根据入流率、堤防材料和堤防形状的变化，将堤防破坏过程分为四个阶段：前两个阶段是垂直方向的决口形成，接着是主河道上下游方向的决口拓宽，最后是决口拓宽减速阶段。他们基于实验开发了一个二维深度平均流模型，成功预测了决口拓宽阶段的决口体积。

Elafy等人（2018）开发了一个使用修正泥沙质量守恒方程的二维深度平均模型，在经典泥沙质量守恒方程中加入了一个新的源项，以考虑决口边坡的侧向侵蚀。该模型通过与实验室测量的稳定和非稳定流通过堤防决口的结果进行比较来验证，然后应用于非粘性土堤的漫顶破坏，以预测决口的演变，包括加深和拓宽阶段以及决口水文过程线。

4. 数值模型

4.1 控制方程

• 水动力方程 ：深度平均浅水方程可以写成向量形式：
  [
  U_t + E_x + F_y + S = 0
  ]
  其中，
  [
  U =
  \begin{pmatrix}
  h \
  q_x \
  q_y
  \end{pmatrix};
  E =
  \begin{pmatrix}
  q_x \
  \frac{q_x^2}{h} + \frac{1}{2}gh^2 \
  \frac{q_xq_y}{h}
  \end{pmatrix}
  ]
  [
  F =
  \begin{pmatrix}
  q_y \
  \frac{q_xq_y}{h} \
  \frac{q_y^2}{h} + \frac{1}{2}gh^2
  \end{pmatrix};
  S =
  \begin{pmatrix}
  0 \
  -gh(S_{ox} - S_{fx}) \
  -gh(S_{oy} - S_{fy})
  \end{pmatrix}
  ]
  这里，(U) 是守恒变量向量，(E) 和 (F) 是通量向量函数，(S) 是源项向量，(t) 是时间，(x) 和 (y) 是水平笛卡尔坐标，(h) 是水流深度，(q_x) 和 (q_y) 分别是 (x) 和 (y) 方向的单位宽度流量，(g) 是重力加速度，(S_{ox}) 和 (S_{oy}) 分别是 (x) 和 (y) 方向的床面坡度，(S_{fx}) 和 (S_{fy}) 分别是 (x) 和 (y) 方向的摩擦坡度。

基于曼宁方程，摩擦坡度 (S_{fx}) 和 (S_{fy}) 可写为：
[
S_{f(x;y)} = \frac{n^2q_{(x;y)}}{\sqrt{q_x^2 + q_y^2}} \left(\frac{h}{10}\right)^{\frac{3}{2}}
]
其中，(n) 是曼宁糙率系数。

• 泥沙方程 ：包含新源项以考虑侧向侵蚀的二维泥沙质量守恒方程可写为：
  [
  (1 - \lambda_p)\frac{\partial z}{\partial t} = -\frac{\partial q_{bx}}{\partial x} - \frac{\partial q_{by}}{\partial y} - \frac{\partial q_{sf}}{\partial x}
  ]
  其中，床面孔隙率 (\lambda_p = 0.43)，(z) 是床面高程，(q_{sf}) 是由边坡失稳引起的泥沙入流率，(q_{bx}) 和 (q_{by}) 分别是 (x) 和 (y) 方向的单位宽度体积床沙输移率。

由于通过堤防决口的水流包括在 (x) 和 (y) 方向的斜坡床面上的流动，基于平板假设获得的临界剪应力需要修正。根据Van Rijn（1993），修正系数可写为：
[
k_1 = \cos\beta_1 \left[1 - \left(\frac{\tan^2\beta_1}{\tan^2\varphi}\right)\right]^{0.5}
]
[
k_2 =
\begin{cases}
\frac{\sin(\varphi - \beta_2)}{\sin\varphi} & \text{对于向下坡的流动} \
\frac{\sin(\varphi + \beta_2)}{\sin\varphi} & \text{对于向上坡的流动}
\end{cases}
]
[
\tau_c = \frac{k_1}{k_2}\tau_{c,f}
]
其中，(k_1) 和 (k_2) 分别是 (x) 和 (y) 方向的修正系数，(\beta_1) 和 (\beta_2) 分别是 (x) 和 (y) 方向的床面坡度角，(\varphi) 是床面材料的休止角，(\overline{\tau} c) 是无量纲修正临界剪应力，(\tau {c,f}) 是平板的无量纲临界剪应力。

剪应力分量可计算为：
[
\tau_{(x;y)} = \rho ghS_{f(x;y)}
]
总剪应力为：
[
\tau = \sqrt{\tau_x^2 + \tau_y^2}
]
无量纲总剪应力表示为：
[
\overline{\tau} = \frac{\tau}{\rho gRD}
]
其中，(R) 是泥沙的淹没比重，(D) 是中值粒径。

单位宽度体积床沙输移率的无量纲分量可表示为：
[
\overline{q} {bx} = \overline{q}_b\cos\theta
]
[
\overline{q} {by} = \overline{q}_b\sin\theta
]
其中，(\theta) 是泥沙输移角度，考虑了主流剪切和 (x) 和 (y) 方向床面坡度的影响，(\overline{q}_b) 是单位宽度的无量纲总体积床沙输移率，根据Meyer - Peter和Müller（1948）（MPM）方程可表示为：
[
\overline{q}_b = \alpha(\overline{\tau} - \overline{\tau}_c)^\beta
]
其中，(\alpha) 和 (\beta) 分别是MPM公式的系数和指数。

泥沙输移角度 (\theta) 可根据Van Bendegom公式（1947，经Talmon等人（1995）修正）计算：
[
\theta = \tan^{-1}
\begin{pmatrix}
\frac{\sin\delta - \frac{1}{1.7}\sqrt{\overline{\tau}}\frac{\partial z}{\partial y}}{\cos\delta - \frac{1}{1.7}\sqrt{\overline{\tau}}\frac{\partial z}{\partial x}}
\end{pmatrix}
]
其中，(\delta) 是床面剪应力的方向角，表示为：
[
\delta = \tan^{-1}\left(\frac{V}{U}\right)
]
其中，(U) 和 (V) 分别是深度平均流速在流向（(x)）和横向（(y)）方向的分量。

(x) 和 (y) 方向的单位宽度体积床沙输移率可计算为：
[
q_{b(x;y)} = \sqrt{RgD^3}\overline{q}_{b(x;y)}
]

4.2 由于坍塌破坏引起的侧向泥沙负荷

由于水流通过决口造成的垂直侵蚀，决口边坡变得陡峭。如果坡度超过土壤休止角，不稳定的土壤段可能会向下坍塌。

Spinewine等人（2002）区分了淹没区（水面以下）和露出区（水面以上）的土壤休止角，指出露出区的休止角 (\varphi_e) 通常大于淹没区的休止角 (\varphi_s)，这是由于水面以上部分饱和土壤的表观凝聚力。他们通过一个小水箱实验测量了这些角度，使用的土壤是粒径均匀为1.8毫米的粗砂，估算的 (\varphi_s) 和 (\varphi_e) 分别为35°和87°。Wu等人（2009）开发了一个二维深度平均模型，包括一个雪崩算法来模拟非粘性土坝因漫顶而决口，当 (\varphi_s) 和 (\varphi_e) 分别设置为33°和79°时，数值模拟结果与实验测量结果吻合良好。

由这种边坡失稳引起的泥沙入流率 (q_{sf}) 可计算为：
[
q_{sf} = \frac{\Delta A\zeta}{\Delta t}
]
其中，(\Delta A) 是向下坍塌的不稳定面积，(\zeta) 是松弛系数。决口右侧淹没区的不稳定面积可估算为：
[
\Delta A = 0.5\Delta x(z_{(i,j)} - z_{(i - 1,j)}) - 0.5(\Delta x)^2\tan\varphi_s
]
其中，(\Delta x) 是 (x) 方向的网格间距，(i) 表示 (x) 方向的节点编号，(j) 表示 (y) 方向的节点编号。松弛系数可表示为：
[
\zeta = (1 - e^{(-\Delta t/T_b)})
]
其中，(\Delta t) 是计算时间步长，(T_b) 是坍塌材料向下坍塌的时间。

4.3 数值解法

• 有限差分格式 ：对于稳定流通过堤防决口，水动力方程（式16 - 5）使用MacCormack显式有限差分格式求解。该格式在时间和空间上都是二阶精确的，允许求解渐变流和急变流，无需对陡峭波前进行特殊处理。求解过程包括一个两步预测 - 校正序列，预测步骤中的空间导数用向前有限差分代替，校正步骤中用向后有限差分代替，这个序列可以每隔一个时间步交替进行。

对于非粘性土堤的漫顶破坏模拟，泥沙质量守恒方程（式16 - 7）与水动力方程（式16 - 5）以耦合方式求解。在预测和校正步骤开始时，检查土堤和漫滩区域的边坡稳定性，以计算 (q_{sf})。预测的床面高程 (z^ ) 用于校正步骤中检查边坡稳定性和评估 (\frac{\partial z}{\partial x}) 项，预测的水流变量 (h^ )、(q_x^ )、(q_y^ ) 用于评估泥沙输移能力 (q_{b(x,y)})，这些值在校正步骤中使用。因此，每个时间步每个变量的最终值考虑了所有其他变量的变化。

数值格式用于计算内部节点的水流变量，而在边界处，水流变量根据每个边界的水流条件进行计算。

• 初始和边界条件 ：对于非稳定流，初始条件由等于入口流量的恒定流量和主河道中相应的恒定水深指定。对于稳定流，初始条件是任意的，因为只关注最终的稳定状态解。

计算域的边界可分为开放边界和固体边界。开放边界包括主河道的上下游端和漫滩出口，固体边界是主河道的墙壁。由于上游边界的水流是亚临界流，所有时间步的 (q_x) 和 (q_y) 都被指定，(h) 从内部节点外推。主河道下游边界的水流也是亚临界流，因此 (q_x) 和 (q_y) 从内部节点外推，对于稳定流，(h) 指定为恒定值，对于非稳定流，(h) 等于传感器2测量的水位过程线。漫滩出口的边界条件是自由落体，因此由水流类型决定：对于亚临界流，(h) ......（此处原文未完整）

以下是一个简单的mermaid流程图，展示了整个堤防决口建模的大致流程：

graph LR
    A[开始] --> B[确定初始条件和边界条件]
    B --> C[选择数值解法（有限差分格式）]
    C --> D[求解水动力方程]
    D --> E[求解泥沙方程]
    E --> F[检查边坡稳定性，计算q_sf]
    F --> G[更新水流变量和泥沙输移能力]
    G --> H{是否达到稳定状态或结束时间}
    H -- 否 --> D
    H -- 是 --> I[输出结果（决口演变、流量等）]
    I --> J[结束]

通过以上的分析和模型，我们可以更深入地理解堤防决口的过程，并进行准确的模拟和预测，为防洪减灾提供有力的支持。

堤防决口建模：原理、方法与数值模拟

5. 模型验证与应用案例

为了确保所建立的堤防决口模型的准确性和可靠性，需要进行模型验证。通常会将模型计算结果与实验室测量数据或实际案例中的观测数据进行对比。

例如，Elafy等人（2018）开发的模型就通过与实验室测量的稳定和非稳定流通过堤防决口的结果进行比较来验证。在验证过程中，会关注多个关键指标，如决口流量、决口宽度变化、淹没区域等。以下是一个简单的验证指标对比表格：
| 验证指标 | 模型计算结果 | 实验室测量结果 | 误差（%） |
| — | — | — | — |
| 决口流量 | [具体数值1] | [具体数值2] | [误差值1] |
| 决口宽度 | [具体数值3] | [具体数值4] | [误差值2] |
| 淹没面积 | [具体数值5] | [具体数值6] | [误差值3] |

在实际应用方面，模型可以用于预测不同情况下的堤防决口后果，为防洪减灾决策提供依据。比如，在面临洪水威胁时，可以模拟不同水位、流量条件下堤防决口的可能性和影响范围，从而提前制定疏散计划、加固堤防等措施。

以下是一个简单的应用流程mermaid流程图：

graph LR
    A[获取洪水预报信息（水位、流量等）] --> B[输入模型参数（堤防材料、断面等）]
    B --> C[运行堤防决口模型]
    C --> D[输出决口可能性、影响范围等结果]
    D --> E{是否需要采取措施}
    E -- 是 --> F[制定防洪减灾措施（疏散、加固等）]
    E -- 否 --> G[继续监测]
    F --> G

6. 影响因素分析

堤防决口过程受到多种因素的影响，了解这些因素对于准确建模和预测至关重要。以下是一些主要影响因素的分析：
- 水流条件 ：包括流量、水位、流速等。流量和水位的增加会加大堤防承受的压力，增加决口的风险。流速则会影响泥沙输移和侵蚀作用，进而影响决口的发展。例如，高流速可能导致更严重的侵蚀，加速决口的形成和扩大。
- 堤防材料 ：不同的堤防材料具有不同的抗侵蚀能力和稳定性。非粘性土堤相对更容易受到侵蚀和破坏，而粘性土堤则具有较好的抗侵蚀性能。此外，材料的颗粒大小、孔隙率等也会影响堤防的性能。
- 堤防断面形状 ：合理的堤防断面形状可以提高堤防的稳定性。例如，较宽的堤顶和缓坡的堤坡可以分散水流压力，减少侵蚀的可能性。相反，过于陡峭的堤坡可能会导致边坡失稳，增加决口的风险。
- 地质条件 ：堤防所在的地质条件，如地基的承载能力、地下水位等，也会对堤防的稳定性产生影响。如果地基承载能力不足，可能会导致堤防下沉、开裂，从而引发决口。

以下是一个影响因素的重要性排序表格：
| 影响因素 | 重要性排序 | 说明 |
| — | — | — |
| 水流条件 | 1 | 直接影响堤防承受的压力和侵蚀作用 |
| 堤防材料 | 2 | 决定堤防的抗侵蚀能力和稳定性 |
| 堤防断面形状 | 3 | 影响水流压力分布和边坡稳定性 |
| 地质条件 | 4 | 对堤防的整体稳定性有潜在影响 |

7. 模型的局限性与改进方向

尽管目前的堤防决口模型已经取得了一定的成果，但仍然存在一些局限性：
- 数据获取困难 ：准确的模型需要大量的输入数据，如地形信息、堤防材料参数、水流条件等。然而，在实际应用中，这些数据的获取可能存在困难，导致模型的准确性受到影响。
- 复杂过程模拟不足 ：堤防决口过程涉及水流、泥沙输移、边坡失稳等多个复杂过程，目前的模型可能无法完全准确地模拟这些过程之间的相互作用。
- 不确定性因素考虑不足 ：实际情况中存在许多不确定性因素，如气候变化、人类活动等，这些因素可能会对堤防决口产生影响，但目前的模型往往难以充分考虑这些因素。

针对这些局限性，可以从以下几个方面进行改进：
- 加强数据采集和处理 ：利用先进的测量技术和数据采集系统，获取更准确、详细的输入数据。同时，开发数据处理方法，提高数据的质量和可用性。
- 改进模型算法 ：不断改进模型的算法，提高对复杂过程的模拟能力。例如，引入更先进的数值方法和物理模型，更好地模拟水流、泥沙输移和边坡失稳等过程。
- 考虑不确定性因素 ：在模型中引入不确定性分析方法，考虑气候变化、人类活动等不确定性因素的影响，提高模型的可靠性和适应性。

以下是一个改进方向的实施步骤列表：
1. 建立数据采集网络，定期收集地形、水流、堤防材料等数据。
2. 对采集到的数据进行清洗、预处理，提高数据质量。
3. 研究和应用新的数值方法和物理模型，改进模型算法。
4. 引入不确定性分析方法，如蒙特卡罗模拟等，评估不确定性因素的影响。
5. 定期对模型进行验证和更新，确保模型的准确性和可靠性。

通过不断改进模型，我们可以更准确地预测堤防决口的发生和发展，为防洪减灾工作提供更有力的支持。

综上所述，堤防决口建模是一个复杂但重要的研究领域。通过深入研究其原理、方法，建立准确可靠的模型，并不断改进和完善，我们可以更好地应对洪水灾害，保护人民生命财产安全和社会经济的稳定发展。