山洪输移漂木演进特征数值试验研究

Numerical Investigation of Large Wood Evolution Characteristics Transported by Flash Floods

摘要

洪水期间河道附近的木材在水流作用下被输运到下游，并产生堆积，造成水位上升，进而加剧山洪灾害。现阶段漂木运动模型主要复现漂木的二维运动状态，但对漂木的拥挤运输、三维堆积行为的模拟存在缺陷，因此，此类模型不适于应用在宽深比较小、蜿蜒崎岖的山区河流。计算流体力学与离散元（CFD‒DEM）耦合的三维流场与漂木运动的数值模型（3D‒3D）能模拟漂木三维输移与堆积行为，将其与水槽试验对比，可支持后续山洪输移漂木的演进特征研究。基于“8·20”三江镇发生的山洪事件，采用3D‒3D模型模拟漂木输移与堆积，结果显示：在不同漂木补给量下，分布在距河段上游100～200 m、400～500 m的漂木分别占总补给量的14%～31%、51%～71%。漂木分布的空间位置并未受补给量的变化而发生显著变化。漂木补给量超过200根时，漂木长度对输移距离的影响并不显著。漂木长度超过3.6 m时，淹没率随漂木补给量的增加而降低。相同漂木长度下，不同漂木补给量对淹没率有影响，但相关关系并不显著。漂木易在河道缩窄处与较大弯道处外侧堆积，导致水位上升，前者堆积行为受河道缩窄处的河岸阻碍影响，后者受弯道处离心力或二次流的影响。上游的河道缩窄处漂木长度从短到长（分别为2.4、3.0、3.6 m）形成堆积体的数量为47、50和52根，形成堆积体的漂木越长越易拦截后续漂木，从而扩大堆积范围。

Abstract

Objective During flood events, wood scattered near river channels is subjected to flowing water and transported downstream. During this process, structures and locations such as bridges, the outer banks of river bends, and channel constrictions are impacted by large wood, which can accumulate and form logjams. The formation of these logjams alters local flow conditions, resulting in increased upstream water levels and enhanced sediment deposition. These conditions directly affect the stability and safety of bridges and dams. Currently, most numerical models primarily consider the coupling of two-dimensional flow fields with two-dimensional large wood motion models (2D‒2D) or three-dimensional flow fields with two-dimensional large wood motion models (3D‒2D). However, two-dimensional large wood motion models cannot adequately represent the vertical structure of logjams, effectively limiting the simulation of large wood movement and accumulation to shallow-water or low-Froude-number conditions. This modeling approach primarily applies to rivers with large width-to-depth ratios in plain regions. However, such rivers generally experience fewer large wood blockage problems. In contrast, mountainous rivers characterized by small width-to-depth ratios and pronounced channel bends present a greater risk of large wood blockage and therefore require primary consideration. Therefore, most existing numerical models exhibit significant limitations that can lead to the underestimation of the impacts associated with large wood blockage hazards. A two-way coupled numerical model integrating three-dimensional flow fields and three-dimensional large wood motion (3D‒3D) is required to enhance the capability of numerical models to accurately reproduce large wood movement and accumulation processes. Methods The three-dimensional flow field and three-dimensional large wood motion two-way coupling numerical model (3D‒3D) presented in this study was based on the finite volume method and employed a second-order nonlinear k-ε model for turbulence calculations. Each piece of large wood was constructed from multiple rigidly connected discrete particles, providing improved collision calculation performance and greater flexibility. The two-way coupling algorithm reflected the influence of large wood on water flow, enabling the analysis of backwater rise caused by logjams. The drag force acting on the large wood was represented by the linear superposition of the forces exerted by the discrete particles constituting the large wood on the water flow. The control equation for large wood considered the combined effects of drag force, water flow acceleration, gravity, virtual mass force, buoyancy, particle collisions, and bed interaction, closely approximating the actual conditions governing large wood dynamics. Results and Discussions The 3D‒3D motion and accumulation model was utilized to simulate the significantly large wood transport and accumulation phenomena in the Xihe River of Sanjiang Town during the Wenchuan "8·20" heavy rainfall-induced mountain torrent and debris flow event. Numerical simulation results showed that the transport distance of large wood under different recruitment scenarios exhibited similar distribution patterns within the river channel. Specifically, large wood distributed within the 100～200 m and 400～500 m sections accounted for 14%～31% and 51%～71% of the total recruitment, respectively. However, the results obtained under larger recruitment conditions were statistically more significant. The findings confirmed that large wood length was only one of several influencing factors, indicating that future studies should investigate additional variables that can affect large wood transport. When large wood recruitment exceeded 200 pieces, the influence of large wood length on transport distance was not significant. The submergence state of large wood varied over time, and a statistical description was obtained using the time-averaged submergence state. Although changes in large wood length affected the number and proportion of submerged wood, the correlation was not significant, indicating that submergence behavior was more strongly influenced by external flow conditions. When the length of the large wood exceeded 3.6 m, the submergence rate decreased as recruitment increased. In the study reach of the Xihe River in Sanjiang Town, no bridge piers were present. Therefore, the study did not consider the direct effects of bridge piers on large wood transport and accumulation. Large wood accumulation tended to occur at channel constrictions and along the outer banks of large river bends, with a stronger tendency toward the latter under the examined conditions. Logjams at channel constrictions were likely caused by large wood being blocked by upstream obstacles under flowing water conditions, eventually reaching a relatively stable state. At larger bends, logjams were more likely influenced by centrifugal forces or secondary flows, causing large wood to deposit along the riverbank or become compressed, increasing friction with the riverbank. Upstream of the logjams, the water depth was found to be 2～3 cm greater than under conditions without logjams. The backwater rise phenomenon can be related to the three-dimensional accumulation pattern and the number of large wood pieces within the logjam, requiring further investigation. At the upstream channel constriction, the numbers of large wood pieces forming logjams from the shortest to the longest length categories were 47, 50, and 52, respectively. Longer pieces of large wood were more likely to intercept subsequent pieces, expanding the accumulation area. Conclusions The 3D‒3D model demonstrated strong capability in reproducing the movement and accumulation of large wood. Simulations of large wood distribution, submergence, and accumulation in the Xihe River segment of Sanjiang Town revealed the evolutionary characteristics of large wood transport during mountain torrent events. The transport of large wood was not entirely confined to the vicinity of the water-flow interface but exhibited distinct three-dimensional movement characteristics over time. Large wood primarily accumulated near channel constrictions and major bends, where logjams were formed. These logjams exhibited a pronounced three-dimensional structure and caused a measurable increase in backwater levels.

Graphical abstract

本刊网刊

漂木对河流生态系统有着重要的保护与恢复作用^[1‒2]。地震、滑坡等因素会导致大量木材散落在河道附近，当山洪发生时，这部分木材受到水流作用向下游输移，可能直接冲击并破坏桥梁等建筑物^[3]。漂木在建筑物处形成堆积体，造成水位上升^[4‒5]或加剧河床冲刷与河岸侵蚀^[6]。此外，漂木堆积体的溃决使流量瞬间放大，影响下游居民生命财产安全^[7‒8]。

De‒Cicco等^[9]利用室内试验研究了桥墩形状对漂木拦截能力的影响，旨在减轻桥梁处漂木堆积引发的相关灾害。Schalko等^[5]开展了一系列水槽试验，认为来流的弗劳德数、漂木堆积的致密性和细颗粒的含量主导了回水上升行为。Ruiz‒Villanueva等^[10]利用录像研究了含漂木水流的主要特征，并拓宽了漂木运输行为的定义范围。Wyss等^[11]在瑞士的格拉特河进行了原位试验，认为漂木堆积过程可以用冲击、旋转和分离3个阶段描述。Kramer等^[12]基于现场调查结果，提出漂木尺寸主导河道漂木运输机制的河流分类方法，进一步识别了漂木运输机制。数值模拟方法相较于上述方法具有更高的灵活性。Mazzorana等^[13]利用ArcGIS软件，基于事件数据集提出了漂木输运的二维建模方法。Panici^[14]提出了仅考虑拖曳力对漂木运动影响的控制方程，用于预测下一个时间点漂木的运动位置，但忽略了其他作用力的影响，还需进一步改进。Persi等^[15]提出了单向耦合的水动力对流‒扩散模型，用以模拟漂木拥挤运输，但是漂木的拥挤运输涉及小范围内大量颗粒对水流的作用，单向耦合这一方法不够完善。Ruiz‒Villanueva等^[10,16]耦合了二维水动力学模型与漂木运动模型，再现了漂木的运动过程，但该耦合模型基于浅水方程，在漂木拥挤运输与三维堆积问题上存在方法的局限性。浅水条件下，部分垂直力对漂木的影响被忽略，漂木运动主要受垂直方向受力引起的二维平面效应控制，因此，漂木主要表现为二维平面运动行为^[17‒18]。上述数值模型均不能反映漂木的三维输移与堆积特性，仅适用于弗劳德数较小或平原地区宽深比大的河流，对宽深比较小、河道蜿蜒崎岖的山区河流存在不足。

本文采用计算流体力学与离散元（CFD‒DEM）双向耦合的数值模拟方法，构建了适用于弗劳德数较大或宽深比较小的山区弯曲河流的漂木输移与堆积模型。利用室内试验数据验证模型准确性后，分析不同漂木长度和补给量下三江镇西河漂木分布、淹没、堆积情况，揭示了山区河流中漂木输移与堆积特性，为山区河流漂木风险评估与管理提供参考。

1 数值模型与验证

1.1 数值模型

1.1.1 水流控制方程

采用三维水流计算软件Nays CUBE为流场计算提供支持，求解方法基于有限体积法（FVM），这一方法已在Fluent、OpenFOAM等计算流体力学软件中得到广泛的验证^[19‒21]。湍流模型采用二阶非线性k‒ε模型（non-linear k‒ε model^[22‒23]，其中：k为湍动能，m²/s²；ε为湍流耗散率，m³/s²）。

连续性方程如下：

\frac{\partial U_i}{\partial x_i}=\mathrm{0}（1）

式中：U_i 为水流在i方向的流速，m/s；x_i 为水流在i方向上的空间坐标，m。

动量方程如下：

\begin{array}{l}\frac{\partial U_i}{\partial t}+\frac{\partial U_i{U}_j}{\partial x_j}=-\frac{\mathrm{1}}{\rho }·\frac{\partial p}{\partial x_i}+{\delta }_{i\mathrm{3}}g-\frac{\mathrm{1}}{\rho }{F}_{\mathrm{d}i}+\\ \frac{\partial }{\partial x_j}\left(\nu \frac{\partial U_i}{\partial x_j}-\overline{u_i\mathrm{\text{'}}{u}_j\mathrm{\text{'}}}\right)\end{array}（2）

式中：t为时间，s；g为重力加速度，m/s²；ρ为流体密度，kg/m³；p为压强，Pa；i、j为方向分量；-\overline{u_i\mathrm{\text{'}}{u}_j\mathrm{\text{'}}}为雷诺应力，m²/s²；ν为运动黏滞系数，m²/s；F_di 为漂木在i方向的拖曳力，N；δ_ij 为克罗内克尔符号（Kronecker symbol），当i和j为同方向时，{\delta }_{ij}=1，反之则为0，其中，j=3时，{\delta }_{ij}为{\delta }_{i\mathrm{3}}，即当i=\mathrm{3}时，{\delta }_{i\mathrm{3}}=1，i\ne \mathrm{3}时，{\delta }_{i\mathrm{3}}=0。

湍动能k计算方程如下：

\frac{\partial k}{\partial t}+\frac{\partial k{U}_j}{\partial x_j}=-\overline{u_i\mathrm{\text{'}}{u}_j\mathrm{\text{'}}}\frac{\partial U_i}{\partial x_j}+\frac{\partial }{\partial x_j}\left[\left(\frac{{\nu }_{\mathrm{t}}}{{\sigma }_k}+\nu \right)\frac{\partial k}{\partial x_j}\right]-\epsilon（3）

式中：σ_k 为模型常数，通常σ_k =1.0；ν_t为湍流运动黏滞系数，m²/s。

湍动能耗散率ε计算方程如下：

\begin{array}{l}\frac{\partial \epsilon }{\partial t}+\frac{\partial \epsilon U_j}{\partial x_j}=-C_{\mathrm{1}}\frac{\epsilon }{k}\overline{u_i\mathrm{\text{'}}{u}_j\mathrm{\text{'}}}\frac{\partial U_i}{\partial x_j}+\\ \frac{\partial }{\partial x_j}\left(\left(\frac{{\nu }_{\mathrm{t}}}{{\sigma }_{\epsilon }}+\nu \right)\frac{\partial \epsilon }{\partial x_j}\right)-C_{\mathrm{2}}\frac{{\epsilon }^{\mathrm{2}}}{k}\end{array}（4）

式中，σ_ε 、C₁与C₂均为模型常数，通常σ_ε =1.3，C₁=1.44，C₂=1.92。

湍流模型的非线性本构方程如下：

\begin{array}{l}-\overline{u_i{u}_j}={\nu }_{\mathrm{t}}{S}_{ij}-\frac{\mathrm{2}}{\mathrm{3}}k{\delta }_{ij}-\frac{k}{\epsilon }{\nu }_{\mathrm{t}}({\alpha }_{\mathrm{1}}(S_{il}{K}_{lj}+S_{jl}{K}_{li})+\\ {\alpha }_{\mathrm{2}}(S_{il}{S}_{lj}-\frac{\mathrm{1}}{\mathrm{3}}{S}_{em}{S}_{me}{\delta }_{ij})+{\alpha }_{\mathrm{3}}(K_{il}{K}_{lj}-\frac{\mathrm{1}}{\mathrm{3}}{K}_{em}{K}_{me}{\delta }_{ij}))\end{array}（5）

{\nu }_{\mathrm{t}}=c_{\mathrm{\mu }}\frac{k^{\mathrm{2}}}{\epsilon }（6）

\begin{array}{c}{c}_{\mathrm{\mu }}=(c_{\mu \mathrm{0}}(\mathrm{1}+c_{nS}{S}^{\mathrm{2}}+c_{nK}{K}^{\mathrm{2}}))/(\mathrm{1}+c_{dS}{S}^{\mathrm{2}}+c_{dK}{K}^{\mathrm{2}}+\\ c_{\mathrm{d}SK}SK+c_{\mathrm{d}S\mathrm{1}}{S}^{\mathrm{4}}+c_{\mathrm{d}K\mathrm{1}}{K}^{\mathrm{4}}+c_{\mathrm{d}SK\mathrm{1}}{S}^{\mathrm{2}}{K}^{\mathrm{2}})\end{array}（7）

\left\{\begin{array}{l}{\alpha }_{\mathrm{1}}=-\mathrm{0.132}\mathrm{5}f(M),\\ {\alpha }_{\mathrm{2}}=\mathrm{0.067}\mathrm{5}f(M),\\ {\alpha }_{\mathrm{2}}=-\mathrm{0.067}\mathrm{5}f(M)\end{array}\right.（8）

f(M)=\frac{\mathrm{1}}{\mathrm{1}+\mathrm{0.01}{M}^{\mathrm{2}}}（9）

M=\mathrm{m}\mathrm{a}\mathrm{x}(S,K)（10）

S=\frac{k}{\epsilon }\sqrt[]{\frac{\mathrm{1}}{\mathrm{2}}{S}_{ij}{S}_{ij}},K=\frac{k}{\epsilon }\sqrt[]{\frac{\mathrm{1}}{\mathrm{2}}{K}_{ij}{K}_{ij}}（11）

S_{ij}=\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i},K_{ij}=\frac{\partial U_i}{\partial x_j}-\frac{\partial U_j}{\partial x_i}（12）

式（5）～（12）中，α₁、α₂与α₃均为非线性项系数，f(M)与M为中间变量，S_ij 、K_ij 分别为应变张量与旋转张量，e、m、l为表示方向的中间变量，S与K分别为应变参数与旋转参数，c_μ为非线性模型常数^[23]，c_μ0、c_nS 、c_nK 、c_dS 、c_dK 、c_dSK 、c_dS1、c_dK1与c_dSK1均为计算c_μ所需模型常数。

1.1.2 漂木控制方程

漂木在水流中受到多种力的共同作用，颗粒的运动方程如式（13）所示：

\begin{array}{l}\sigma A_{\mathrm{3}}{d}^{\mathrm{3}}\frac{\mathrm{d}{U}_{\mathrm{p}i}}{\mathrm{d}t}=\frac{\mathrm{1}}{\mathrm{2}}{C}_{\mathrm{d}}\rho {\lambda }_{\mathrm{A}\_\mathrm{s}\mathrm{u}\mathrm{b}}{A}_{\mathrm{2}}{d}^{\mathrm{2}}\left|U_i-U_{\mathrm{p}i}\right|(U_i-U_{\mathrm{p}i})+\\ \rho {\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{A}_{\mathrm{3}}{d}^{\mathrm{3}}\frac{\mathrm{d}{U}_i}{\mathrm{d}t}+\sigma A_{\mathrm{3}}{d}^{\mathrm{3}}{\delta }_{i\mathrm{3}}g+\\ C_{\mathrm{m}}\rho {\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{A}_{\mathrm{3}}{d}^{\mathrm{3}}\left(\frac{\mathrm{d}{U}_i}{\mathrm{d}t}-\frac{\mathrm{d}{U}_{\mathrm{p}i}}{\mathrm{d}t}\right)-\\ \rho {\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{A}_{\mathrm{3}}{d}^{\mathrm{3}}{\delta }_{i\mathrm{3}}g+F_{\mathrm{p}i}+\varphi F_{\mathrm{b}i}\end{array}（13）

式中：\sigma为离散颗粒密度，kg/m³；U_pi 为颗粒在i方向的速度，m/s；C_d为拖曳力系数^[24]；C_m为虚质量系数；A₂为颗粒的二维形状系数，A₂=π/4；A₃为颗粒的三维形状系数，A₃=π/6；λ_{A_sub}为颗粒的投影面积在流向的淹没比；λ_{V_sub}为颗粒体积的淹没比；d为耦合形成漂木的颗粒直径，m；γ为颗粒的比重；F_pi 为不同颗粒间的碰撞力，N；F_bi 为床层对漂木的作用力，N；ϕ为评估床层作用的系数，接触床面时ϕ=1，反之ϕ=0。

床层作用力一般可以分为漂木与床层的摩擦力、床层梯度的重力效应。床层小梯度假设下，床层作用力可以简化为：

F_{\mathrm{b}i}=\frac{\mu A_{\mathrm{3}}{d}^{\mathrm{3}}g(\rho {\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}-\sigma )({\delta }_{i\mathrm{1}}+{\delta }_{i\mathrm{2}})U_{\mathrm{p}i}}{{‖U_{\mathrm{p}i}‖}_{\mathrm{2}}}+F_{\mathrm{g}i}（14）

式中：μ为床层摩擦系数^[25]；{\delta }_{i\mathrm{1}}与{\delta }_{i\mathrm{2}}均为克罗内克尔符号，下标1和2为两正交方向；{‖U_{\mathrm{p}i}‖}_{\mathrm{2}}为颗粒速度的2阶范数；F_{\mathrm{g}i}为重力分量，N。

式（13）经过上述简单替换可以得到离散颗粒加速度与影响因素的关系，表达式为：

\begin{array}{l}\frac{\mathrm{d}{U}_{\mathrm{p}i}}{\mathrm{d}t}=\frac{\mathrm{1}}{\mathrm{2}}{C}_{\mathrm{d}}{\lambda }_{\mathrm{A}\_\mathrm{s}\mathrm{u}\mathrm{b}}\frac{A_{\mathrm{2}}}{A_{\mathrm{3}}}·\frac{\mathrm{1}}{d}·\frac{\mathrm{1}}{\gamma +{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{C}_{\mathrm{m}}}\left|U_i-U_{\mathrm{p}i}\right|(U_i-U_{\mathrm{p}i})+\\ \frac{{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}(\mathrm{1}+C_{\mathrm{m}})}{\gamma +{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{C}_{\mathrm{m}}}·\frac{\mathrm{d}{U}_i}{\mathrm{d}t}+\frac{\gamma -{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}}{\gamma +{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{C}_{\mathrm{m}}}{\delta }_{i\mathrm{3}}g+\\ \frac{\mathrm{1}}{A_{\mathrm{3}}{d}^{\mathrm{3}}\rho (\gamma +{\lambda }_{\mathrm{V}\_\mathrm{s}\mathrm{u}\mathrm{b}}{C}_{\mathrm{m}})}(F_{\mathrm{p}i}+\varphi F_{\mathrm{b}i})\end{array}（15）

刚性连接多个离散颗粒构成单根漂木，该方法碰撞计算效果好、灵活性强。双向耦合方法额外考虑了颗粒对水流的反向作用，在CFD‒DEM耦合时具有更高的理论模拟精度。漂木对水流的拖曳力采用构成漂木的多个离散颗粒对水流阻力作用的线性叠加表示，表达式为：

F_{\mathrm{d}i}=\frac{\mathrm{1}}{\mathrm{2}}\rho C_{\mathrm{d}}{\lambda }_{\mathrm{A}}{\displaystyle \sum _{q=\mathrm{1}}^{N_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}}}{A_q}_{,\mathrm{a}\mathrm{v}\mathrm{g}}\left|U_i-U_{\mathrm{p}i}\right|(U_i-U_{\mathrm{p}i})（16）

式中：λ_A为颗粒的投影面积顺流向的淹没比；A_q,avg为第q个颗粒顺流向的投影面积，m²；N_cell为网格内粒子数量。

1.2 模型验证

Kimura等^[26]利用室内水槽试验验证了二维流场与漂木运动耦合模型（2D‒2D）和三维流场与二维漂木运动耦合模型（3D‒2D）的模拟精度。本文在相同工况下验证三维流场与漂木运动耦合模型（3D‒3D）的模拟精度，并将其与3D‒2D模型的模拟效果进行对比。工况Ⅰ、Ⅱ分别用于验证漂木运动模型与堆积模型的准确率，工况详情如表1所示。

水槽形态为正弦函数型，水槽宽为0.2 m。来流中漂木补给数量共110根，待水流流态稳定后，在水槽前端位置按水平方向内随机角度自由落入水槽中，漂木密度为600 kg/m³。工况Ⅱ系列中4个模型桥墩宽度均为1 cm，除弯道出口宽度为4 cm，其余出口均为3 cm。漂木的运动采用漂木理查德松数（DRI，记为I_d）描述^[26]，表达式为：

I_{\mathrm{d}}=\frac{g(\rho -{\rho }_{\mathrm{t}})h}{\rho u^{\mathrm{2}}},{\rho }_{\mathrm{w}}\ge {\rho }_{\mathrm{t}}（17）

式中：ρ_w为水的密度，kg/m³；ρ_t为漂木密度，kg/m³；u为流速，m/s；h为水深，m。

1.2.1 运动模型验证

工况Ⅰ用于检验模型结果是否遵循试验结果中漂木的基本运动状态。在相同时间点下对比3D‒3D、3D‒2D与水槽试验结果，如图1所示。

由图1可知，与试验结果相比，3D‒2D与3D‒3D模型均能模拟出漂木群的整体运动趋势，对再现漂木在弯道中因惯性而向水道外侧移动的行为具有较强的模拟能力。但两类模型在细节模拟上具有明显的差别。3D‒3D模型与试验结果得到的漂木运动情况表现出较高的一致性，较好地模拟了部分漂木的拥挤运输行为。

3D‒2D模型由于缺少对漂木垂向方向的运动模拟，不能体现漂木的淹没情况，漂木聚团的现象并不属于拥挤运输行为的体现。定义漂木通量Φ为不同时刻通过水槽某一特定断面（图1中红线）的漂木数量与总补给量之比，以此量化漂木的运动情况。试验与数值模拟条件下漂木通量的变化如图2所示。由图2可知，3D‒3D数值模型得出的漂木通量较3D‒2D数值模型更加接近水槽试验结果。

1.2.2 堆积模型验证

堆积模型的验证采用工况Ⅱ系列条件，模型桥墩的设置与Kimura等^[26]的水槽试验一致，主要从桥墩对受困漂木沉没率的影响论证堆积模型的合理性。沉没率R_s能体现漂木堆积后的三维结构特性，表达式为：

R_{\mathrm{s}}=\frac{N_{\mathrm{d}}}{N_{\mathrm{t}}}（18）

式中，N_t与N_d分别表示被困漂木数量与被困漂木中沉入水下（未接触自由水面）的漂木数量。

沉没率与理查德松数的关系如图3所示。由图3可知，本文3D‒3D模型数值模拟结果与Kimura等^[26]的水槽试验结果表现出良好的一致性。当I_d较大时，漂木的浮力起主导作用，受困漂木只有少部分处于沉没状态，堆积结构呈现二维形态；当I_d较小时，漂木的惯性力起主导作用，受困漂木倾向沉没状态，堆积结构呈现三维形态^[6]，这一复杂堆积结构使R_s变化范围更大。

漂木堆积将导致上游水位上涨^[27]。图4为基准线处水深变化情况。图4中，定义上游基准线为漂木堆积体前缘位置（图4中红线），1、2、3、4和5分别为障碍物拦截后形成的5个出水断面的序号。图4（b）中，无漂木时采用相同位置作为对比。漂木在运动过程中受到模型桥墩的阻挡形成了堆积体。出口2、5上游基准线处水深如图5所示，当漂木在上游桥墩前发生堆积后，水深分别加深1.7、1.5 mm。在距离基准线约0.12、0.15 m处，出口2、5恢复至正常水深。不论是否发生漂木堆积，最大水面坡降均集中于桥墩间隔范围内。

2 漂木演进特征模拟

2.1 研究区域概况

研究区域为四川省汶川县三江镇西河某一河段，如图6所示^[28]。研究河段范围内共有3座桥梁，计算域网格划分如图7所示。其中：C桥建造于2020年之后，本次研究可以忽略该桥的影响；A桥跨度较小，主河道中没有桥墩的存在；B桥为悬索桥，桥面较窄、跨度较大，主河道中无桥墩存在，几乎不会发生漂木堆积问题。研究区3座桥梁均可排除桥墩的影响导致漂木在桥梁上游的堆积问题，因而可以简化河道模型。研究区段的上游控制断面位于A桥上游约170 m处。

2019年8月19日晚至20日晨，四川省汶川县突发强降雨，诱发特大山洪泥石流（简称“8·20”山洪）。水流携带漂木在多处形成了堆积体，堆积体溃决后产生了流量放大效应，造成下游多处严重灾害^[29‒30]。在此次事件中，三江镇受灾严重，共有6人遇难、9人失联。

2.2 数值模型条件

图8为汶川“8·20”三江镇西河洪水流量过程^[31]。由图8可知，在“8·20”山洪发生时段内，三江镇西河最大洪峰流量为524.14 m³/s。将实际地形高程与无人机拍摄画面对比估计得到河段下游水位高程为1 016 m。由于发生漂木输运耗时较短，故可视为河流流量与水位高程不发生变化。本次数值模拟采用的CFD‒DEM参数设置如表2所示。

上下游边界条件分别由来流流量与水位高程控制，上游来流流量为500 m³/s，下游水位高程为1 016 m。探究漂木长度L与漂木补给量N变化条件下西河部分河段漂木输移与堆积情况。漂木补给量N分别为100根（A）、200根（B）、300根（C）（A、B、C为相对应的工况编号），漂木长度分别为2.4 m（1）、3.0 m（2）、3.6 m（3）（1、2、3为对应的工况编号），例如工况B2表示上游有200根3.0 m长的漂木补给。通过分析现场影像资料，获得漂木尺寸数据，通过无人机航测获取了三江镇河道及沿河村镇的数字地表模型（DSM），分辨率为5 m。由于模拟基于水位与流量不发生变化这一假设，计算时长的选择不宜过长。当计算时长大于360 s后，漂木的运动与堆积变化并不显著，不影响最终结果，综合考虑计算量与模拟效果后，总计算时长确定为360 s。

2.3 数值试验分析

对宽深比较大的山区河流或平原河流，由于漂木的输移与堆积通常不涉及三维状态，所以采用二维模型能降低计算量且不影响实际模拟效果。但在宽深比不大的山区河流中，漂木受河道地形与流动复杂性的影响，将会表现为三维输运与堆积状态，现有二维模型描述此类问题时存在缺陷。3D‒3D模型已被证实能捕捉到漂木拥挤运输和三维堆积现象，这在以往漂木数值模拟中是很少考虑的内容。

2.3.1 漂木分布

漂木补给量较少时，其运动过程相对简单；但当补给量较大时，直接研究其运动过程就十分复杂。在上游有大量漂木补给时，漂木的运动过程通常不是关注的核心内容，其分布和堆积状况更直观地反映了漂木造成的影响。

采用漂木分布率P量化分析漂木在河段中的分布情况，代表一定范围河段内漂木数量N_p与当前研究区中漂木总数N_a的比值，表达式为：

P=\frac{N_{\mathrm{p}}}{N_{\mathrm{a}}}\times \mathrm{100}\mathrm{\%}（19）

不同漂木补给量下，漂木沿程分布情况如图9所示。由图9可知：不同工况下漂木的分布情况相似，在距离上游控制断面100～200 m、400～500 m范围内漂木分布率分别处于14%～31%与51%～71%，这两段河段漂木分布率总计达80%以上。同一种漂木补给量下，漂木长度变化带来的影响仅在补给量较低（100根）时较明显，而在较大补给量时这一现象并不显著，且补给量的变化并没有改变研究区中漂木的空间分布模式。上述结果与Picco等^[32]对皮亚韦河的漂木的输移研究结论一致，证实漂木的运输是一个复杂的机制，漂木长度对输移距离的影响较弱。采用更大的漂木补给量得出的结果更具有统计意义，后续主要对漂木补给量较大工况（工况C）的堆积情况进行分析。

2.3.2 漂木淹没

漂木的淹没状态体现了漂木在水中的垂向动态变化过程。不同时间点漂木的淹没状态受自身运动状态、水流状态等多种条件的控制，单纯研究某一时间点的漂木淹没情况有较大的偶然性。本文考虑计算60、120、180、240、300与360 s时漂木的淹没情况，在一定程度上减弱偶然性带来的影响，得到更具说服力的统计性结果。一般情况下，漂木密度均小于水的密度，所以漂木难以达到完全淹没状态。仅当其形成三维堆积结构时，受到结构控制并处于堆积体下方的漂木才可能处于完全淹没状态。由于三维结构的复杂性，以漂木的淹没情况量化堆积结构的三维特性。不同漂木补给量、长度下漂木淹没情况分别如图10、11所示。

由图10可知：漂木补给量相同情况下，较长的漂木更难以被水淹没；当漂木补给量增大时，复杂的堆积结构导致漂木淹没不符合上述规律。工况C2中淹没漂木数量的平均值均大于工况C1与C3，但是工况C1与C2的淹没漂木根数平均值仅相差4.4根，也就是每100根相差约1.4根淹没漂木，这是在数值模拟中可以接受的误差大小，所以可以认为在同种漂木补给条件下，漂木长度越长，漂木淹没数量N_s越少。为考虑相同漂木长度下不同漂木补给量对漂木淹没行为的控制结果，本文采用淹没率S_d来进行分析，其定义式为：

S_{\mathrm{d}}=\frac{N_{\mathrm{s}}}{N_{\mathrm{a}}}\times \mathrm{100}\mathrm{\%}（20）

由图11可知：相同漂木长度下，不同漂木补给量对淹没率有影响，但相关性并不显著。漂木较短（2.4 m）时，中等漂木补给量（200根）条件下出现最大淹没率均值；漂木长度适中（3.0 m）时，中等漂木补给量条件下出现最小淹没率均值；当漂木足够长（3.6 m）时，淹没率随着漂木补给量的增大而减小。漂木长度的增加导致其受到外界的作用更复杂，而复杂的作用得到的结果往往更贴近于现实情况，所以当漂木长度超过一定阈值时，淹没率随着漂木补给量的增加而降低的结论更加具有现实意义。

2.3.3 漂木堆积

西河漂木堆积情况如图12所示。以工况C3为例，漂木在河道缩窄处（CS37）与弯道处（CS100）发生了两处明显的堆积，其他C类工况均发生类似堆积现象。河道缩窄处漂木堆积现象可能是河道左岸缩窄导致外侧漂木撞击缩窄处前缘，后在水流的推力作用下在此处形成了一个相对稳定的堆积体。弯道处漂木堆积现象可能是漂木受到离心力或二次流的作用导致倾向于弯道外侧运动，运动过程中受到岸滩的摩擦、部分凸起石块的阻挡进而形成了堆积体。

漂木在河道缩窄处、弯道处的堆积情况如图13所示。

由图13可知，漂木长度明显影响堆积体的形成。较长的漂木倾向于形成更大的堆积区域，从而更容易拦截并形成含有更多漂木的堆积体。这一现象表明漂木的尺寸直接影响后续形成堆积体的漂木数量，以最靠近上游的河道缩窄处形成的漂木堆积体为例，漂木长度从短到长形成漂木堆积体的数量为47、50和52根。漂木在河道缩窄处（CS37）形成堆积体后导致周围一定范围内水深增加。无漂木最大水深为5.15 m；有漂木最大水深范围为5.17～5.18 m。由于堆积体并未在桥墩或涵洞等过水断面较窄处形成，所以水深增加较小。上游回水深度的增加可能与形成堆积体的漂木数量与堆积模式有关。

3 结论

三维流场与漂木运动耦合模型（3D‒3D）被证实具有复现河道中漂木三维输移与堆积的能力。基于汶川“8·20”特大山洪泥石流灾害事件，模拟得到了山区河流中漂木可能存在的分布、淹没与堆积情况，分析了山洪输移漂木演进特征，主要结论如下：

1）研究漂木在河流中的分布是针对较多漂木补给时的可行方法。西河在不同工况下漂木分布行为相似，均处于距上游断面100～200 m、400～500 m区段，漂木分布率分别为14%～31%与51%～71%。当漂木补给量超过200根时，漂木长度对输移距离的影响并不显著。

2）漂木淹没情况一定程度上反映了漂木的三维运动特性。工况C1、C2和C3中，淹没漂木数量每百根相差约1.4根淹没漂木，在数值模拟误差接受范围内，故可以认为，在相同漂木补给条件下，漂木长度越长，淹没漂木数量越少。当漂木长度超过3.6 m时，更多的漂木补给导致更小的漂木淹没率。

3）数值模型中，易发生漂木堆积的位置主要在河道缩窄处与较大弯道处外侧，前者被认为可能是由于障碍物的直接阻挡作用，而后者被认为可能是弯道离心力或二次流导致漂木在河岸的沉积。形成堆积体的漂木越长，越易拦截后续漂木，进而形成更大的堆积体。堆积体上游有回水深度增加的现象，由于堆积位置不处于桥梁或涵洞等过水断面较窄处，所以有漂木时的水深较无漂木时仅加深2～3 cm。

参考文献

[1] Gurnell A, England J, Burgess‒Gamble L.Trees and wood:Working with natural river processes[J].Water and Environment Journal,2019,33(3):342‒352. doi:10.1111/wej.12426

[2] Wohl E, Bledsoe B P, Fausch K D,et al.Management of large wood in streams:An overview and proposed framework for hazard evaluation[J].JAWRA Journal of the American Water Resources Association,2016,52(2):315‒335. doi:10.1111/1752-1688.12388

[3] Schmocker L, Hager W H.Probability of drift blockage at bridge decks[J].Journal of Hydraulic Engineering,2011,137(4):470‒479. doi:10.1061/(asce)hy.1943-7900.0000319

[4] Ruiz‒Villanueva V, Díez‒Herrero A, Bodoque J M,et al.Large wood in rivers and its influence on flood hazard[J].Cuadernos de Investigación Geográfica,2014,40(1):229‒246. doi:10.18172/cig.2523

[5] Schalko I, Schmocker L, Weitbrecht V,et al.Backwater rise due to large wood accumulations[J].Journal of Hydraulic Engineering,2018,144(9):04018056. doi:10.1061/(asce)hy.1943-7900.0001501

[6] Okamoto T, Takebayashi H, Sanjou M,et al.Log jam formation at bridges and the effect on floodplain flow:A flume experiment[J].Journal of Flood Risk Management,2020,13(S1):e12562. doi:10.1111/jfr3.12562

[7] Chen Jiangang, Wang Daozheng, Zhao Wanyu,et al.Laboratory study on the characteristics of large wood and debris flow processes at slit-check dams[J].Landslides,2020,17(7):1703‒1711. doi:10.1007/s10346-020-01409-3

[8] Wang Daozheng, Wang Xingang, Chen Xiaoqing,et al.Analysis of factors influencing the large wood transport and block-outburst in debris flow based on physical model experiment[J].Geomorphology,2022,398:108054. doi:10.1016/j.geomorph.2021.108054

[9] De‒Cicco P N, Paris E, Solari L,et al.Bridge pier shape influence on wood accumulation:Outcomes from flume experiments and numerical modelling[J].Journal of Flood Risk Management,2020,13(2):e12599. doi:10.1111/jfr3.12599

[10] Ruiz‒Villanueva V, Mazzorana B, Bladé E,et al.Characterization of wood-laden flows in rivers[J].Earth Surface Processes and Landforms,2019,44(9):1694‒1709. doi:10.1002/esp.4603

[11] Wyss A, Schalko I, Weitbrecht V.Field study on wood accumulation at a bridge pier[J].Water,2021,13(18):2475. doi:10.3390/w13182475

[12] Kramer N, Wohl E.Rules of the road:A qualitative and quantitative synthesis of large wood transport through drainage networks[J].Geomorphology,2017,279:74‒97. doi:10.1016/j.geomorph.2016.08.026

[13] Mazzorana B, Hübl J, Zischg A,et al.Modelling woody material transport and deposition in alpine rivers[J].Natural Hazards,2011,56(2):425‒449. doi:10.1007/s11069-009-9492-y

[14] Panici D.An experimental and numerical approach to modeling large wood displacement in rivers[J].Water Resources Research,2021,57(7):e2021WR029860. doi:10.1029/2021wr029860

[15] Persi E, Petaccia G, Sibilla S,et al.A one-way coupled hydrodynamic advection-diffusion model to simulate congested large wood transport[J].Hydrology,2021,8(1):21. doi:10.3390/hydrology8010021

[16] Ruiz‒Villanueva V, Bladé E, Sánchez‒Juny M,et al.Two-dimensional numerical modeling of wood transport[J].Journal of Hydroinformatics,2014,16(5):1077‒1096. doi:10.2166/hydro.2014.026

[17] Kang T, Kimura I.Computational modeling for large wood dynamics with root wad and anisotropic bed friction in shallow flows[J].Advances in Water Resources,2018,121:419‒431. doi:10.1016/j.advwatres.2018.09.006

[18] Kang T, Kimura I, Shimizu Y.Numerical simulation of large wood deposition patterns and responses of bed morphology in a braided river using large wood dynamics model[J].Earth Surface Processes and Landforms,2020,45(4):962‒977. doi:10.1002/esp.4789

[19] Ma Jianxin, Shi Xi, Xie Xiaogang,et al.Study on influence of FLUENT rigid vegetation diameter on water flow characteristics of open channel[J].Water Resources and Power,2023,41(9):116‒119.

[20] 马健鑫,石喜,谢晓刚,等.基于FLUENT刚性植被直径对明渠水流特性的影响[J].水电能源科学,2023,41(9):116‒119.

[21] Liu Hanyue, Ren Chunping, Huang Xuzhong.Numerical study of two-phase flow based on OpenFOAM for dynamic characteristics of deflecting flow of a certain spillway[J].Water Resources and Power,2024,42(1):115‒119.

[22] 刘晗月,任春平,黄旭中.基于OpenFOAM对某溢洪道挑流水动力特性的两相流数值研究[J].水电能源科学,2024,42(1):115‒119.

[23] Shaheed R, Mohammadian A, Kheirkhah Gildeh H.A comparison of standard k-ε and realizable k-ε turbulence models in curved and confluent channels[J].Environmental Fluid Mechanics,2019,19(2):543‒568. doi:10.1007/s10652-018-9637-1

[24] Kimura I, Hosoda T.A non-linear k-ε model with realizability for prediction of flows around bluff bodies[J].International Journal for Numerical Methods in Fluids,2003,42(8):813‒837. doi:10.1002/fld.540

[25] Ali S, Hosoda T, Kimura I.A non-linear k-.EPSILON.Model to predict the spatial change of turbulent structures in large scale vortices[J].Journal of Applied Mechanics,2007,10:723‒732. doi:10.2208/journalam.10.723

[26] Schiller L, Naumann A.A drag coefficient correlation[J].Zeits-chrift des Vereins Deutscher Ingenieure,1935,77:318‒320.

[27] Kimura I, Kang T, Kato K.3D‒3D computations on submerged-driftwood motions in water flows with large wood density around driftwood capture facility[J].Water,2021,13(10):1406. doi:10.3390/w13101406

[28] Kimura I, Kitazono K.Effects of the driftwood Richardson number and applicability of a 3D‒2D model to heavy wood jamming around obstacles[J].Environmental Fluid Mechanics,2020,20(3):503‒525. doi:10.1007/s10652-019-09709-6

[29] Follett E, Schalko I, Nepf H.Momentum and energy predict the backwater rise generated by a large wood jam[J].Geophysical Research Letters,2020,47(17):e2020GL089346. doi:10.1029/2020gl089346

[30] Chen Jiangang, Fei Gaogao, Wang Xi'an,et al.Advances on disaster effects of drift wood in flash flood debris flows[J].Advances in Science and Technology of Water Resources,2022,42(3):104‒111. doi:10.3880/j.issn.1006-7647.2022.03.018

[31] 陈剑刚,费高高,王喜安,等.漂木对山洪泥石流运动致灾影响研究进展[J].水利水电科技进展,2022,42(3):104‒111. doi:10.3880/j.issn.1006-7647.2022.03.018

[32] Fei Gaogao, Wang Xiekang.A review of large wood dynamics relevant to hazard characteristics for built structures[J].Geomorphology,2024,453:109152. doi:10.1016/j.geomorph.2024.109152

[33] 大众网.汶川三江镇一农贸市场被冲垮,村委会、派出所成孤岛[EB/OL].(2019‒08‒22)[2024‒07‒01].

[34] Yang Qingyuan, Wang Xiekang, Sun Yi,et al.Numerical investigation on a flash flood disaster in streams with confluence and bifurcation[J].Water,2022,14(10):1646. doi:10.3390/w14101646

[35] Picco L, Tonon A, Rainato R,et al.Bank erosion and large wood recruitment along a gravel bed river[J].Journal of Agricultural Engineering,2016,47(2):72. doi:10.4081/jae.2016.488