浅水数值模拟备受关注^[1-9]. 基于Godunov思想的有限体积法能够自动捕捉激波/间断，在浅水模拟方面得到了长足的发展^[10-13]，从一维到平面二维；从清水定床到泥沙运动；从均匀到非均匀沙；从明渠流到水下异重流^[14-21]. 这类方法的时间步长受到CFL（Courant-Friedrichs-Lewy）稳定性条件的限制. 以往数学模型普遍采用整体时间步长（global-time-step, GTS）理念，取所有单元最小的 $\Delta t_i^{{\rm{CFL}}}$作为整个计算区域的时间步长 $\Delta {t_{\rm{g}}}$（即 $\Delta {t_{\rm{g}}} = {\min _{i = 1,\,2,\,3,\,\cdots,\,N_{\rm c}}}(\Delta t_i^{{\rm{CFL}}})$），大大限制了模型的计算效率^[22-24]. 许多研究人员着眼于局部时间步长（local-time-step, LTS）技术的研究^[22-28]. 按照局部时间步长技术，每个单元采用与CFL条件相适应的足够大的时间步长，从而提高计算效率. Crossley等^[27-28]针对一维圣维南方程，实现了局部时间步长算法. Sanders^[22]提出针对平面二维浅水方程的分级局部时间步长算法. Dazzi等^[24]实现了对于整体计算区域每个单元时间步长选取的动态检查，若该单元当前的时间步长超过CFL条件，则将局部时间步长级数降级. 对浅水模型计算效率的提高还包括并行技术的发展^{[4, 29-31]}. 目前，局部时间步长技术的应用全部针对于清水定床过程. 本文利用三角形网格离散计算区域，采用能够充分反映水-沙-床相互作用的平面二维完整控制方程. 利用有限体积法求解控制方程，用HLLC近似黎曼算子估算界面数值通量，建立平面二维水沙耦合数学模型. 数值更新过程采用局部时间步长技术.

控制方程组包括挟沙水流的质量和动量守恒方程、分粒径组的随流泥沙质量守恒方程、底床所有组份泥沙质量守恒方程（即底床变形方程）、分粒径组份的底床泥沙质量守恒方程（即基于底床活动层概念的非均匀泥沙输移数学模式）. 将控制方程写成向量形式^{[21, 32]}如下：

(1) $\frac{{\partial {{U}}}}{{\partial t}} + \frac{{\partial {{F}}}}{{\partial x}} + \frac{{\partial {{G}}}}{{\partial y}} = {{{S}}_1} + {{{S}}_2}.$

(2a) ${{U}}{\text{ = }}{\left[ {h,hu,hv,h{c_k},{z_{\text{b}}},\delta {f_{{\rm a},k}}} \right]^{\text{T}}}\tag{2a},$

(2b) ${{F}}{\text{ = }}{\left[ {hu,h{u^2} + g{h^2}/2,huv,hu{c_k},0,0} \right]^{\text{T}}}\tag{2b},$

(2c) ${{G}} = {\left[ {hv,huv,h{v^2} + g{h^2}/2,hv{c_k},0,0} \right]^{\text{T}}}\tag{2c},$

(2d) ${{{S}}_1} = \left[ {\begin{array}{*{20}{c}} {({E_{\text{T}}} - {D_{\text{T}}})/(1 - p)} \\ {gh{S_{{\text{b}}x}}} \\ {gh{S_{{\text{b}}y}}} \\ {{E_k} - {D_k}} \\ { - ({E_{\text{T}}} - {D_{\text{T}}})/(1 - p)} \\ {({E_k} - {D_k})/(1 - p) - {f_{{\rm s},k}}\partial \eta /\partial t} \end{array}} \right], \tag{2d}$

(2e) ${{{S}}_2} = \left[ {0} ,\; {{{ - }}gh{S_{{\text{f}}x}}},\; { - gh{S_{{\text{f}}y}}} ,\; {0} ,\; {0},\; {0} \right]^{\rm T}. \tag{2e}$

式中： ${{U}}$ 为守恒变量向量； ${{F}}$ 和 ${{G}}$ 分别为 $x$ 和 $y$ 方向的通量向量； ${{{S}}_1}$、 ${{{S}}_2}$ 为源项向量； $t$ 为时间； $x\text{、} y$ 为空间水平坐标； $h$ 为水深； $u\text{、}v$ 为 $x\text{、}y$ 方向上的深度积分平均流速； $g$ 为重力加速度； ${c_k}$ 为第 $k$ 粒径组沙的深度平均体积含沙量（ $k = 1,2,3,\cdots,{N_{{\rm{sps}}}}$，其中 ${N_{{\rm{sps}}}}$ 为粒径分组数）； ${z_{\rm{b}}}$ 为底床高程； $\delta $ 为活动层厚度； ${f_{{\rm{a}},k}}$ 为活动层泥沙体积百分比； ${E_k}$、 ${D_k}$ 为第 $k$ 粒径组沙在水流与床面之间的上扬和沉降通量， ${E_{\rm{T}}} = \sum {{E_k}} $， ${D_{\rm{T}}} = \sum {{D_k}} $； $p$ 为泥沙孔隙率； ${S_{{\rm{b}}x}} = {\rm{ - }}\partial {z_{\rm b}}/\partial x$ 和 ${S_{{\rm{b}}y}} = {\rm{ - }}\partial {z_{\rm b}}/\partial y$ 为 $x\text{、}y$ 方向的底坡； ${S_{{\rm{f}}x}}$ 和 ${S_{{\rm{f}}y}}$ 为水流在 $x\text{、}y$ 方向的阻力坡度； ${\rho _{\rm{w}}}$ 和 ${\rho _{\rm{s}}}$ 分别为清水和泥沙密度； ${f_{{\rm{s}},k}}$ 为活动层下交界面处的泥沙体积百分比； $\eta = {z_{\rm{b}}} - \delta $ 为存储层与交换层界面的高程.

阻力坡度采用曼宁公式估算：

(3) ${S_{{\rm{f}}x}} = \frac{{{n^2}u\sqrt {{u^2} + {v^2}} }}{{{h^{{4 / 3}}}}}, \;\; {S_{{\rm{f}}y}} = \frac{{{n^2}v\sqrt {{u^2} + {v^2}} }}{{{h^{{4 / 3}}}}}.$

式中： $n$ 为曼宁糙率系数. 分组粒径的泥沙沉降和上扬通量计算采用下式：

(4) ${D_k} = {\alpha _k}{c_k}{\omega _k}(1 - {\alpha _k}{c_k}), \;\; k = 1,2, \cdots ,{N_{{\rm{sps}}}};$

(5) ${E_k} = {\alpha _k}{c_{{\rm e},k}}{\omega _k}(1 - {\alpha _k}{c_{{\rm e},k}}), \;\; k = 1,2, \cdots ,{N_{{\rm{sps}}}}.$

式中: ${\omega _k}$ 为分粒径组沙在静水条件下的沉降速度，采用张瑞瑾沉速公式计算^[33]； ${c_{{\rm{e}},k}}$ 为第 $k$ 粒径组沙的水流挟沙力； ${\alpha _k}$ 为近底泥沙含量与深度平均含沙量的比值. 对于动床溃坝洪水，采用Wu^[33]的全沙经验公式^[33]计算水流挟沙力：

(6) $\begin{split} {c_{{\rm{e}},k}} = & \{ 0.005\,\,3\sqrt {sgd_k^3} {\left[ {{\theta _k}/{\theta _{{\rm c},k}} - 1} \right]^{2.2}} + \\ & 0.000\,\,026\,\,2{\left[ {\left( {{\theta _k}/{\theta _{{\rm c},k}} - 1} \right)U/{\omega _k}} \right]^{1.74}}\} /\left( {hU} \right). \end{split}$

式中： ${\theta _{{\rm{c}},k}} \!=\! {\theta _{\rm{c}}}{({p_{{\rm{e}},k}}/{p_{{\rm{h}},k}})^{ - 0.6}}$，其中 ${p_{{\rm{h}},k}} \!=\! \sum\nolimits_j {[{F_j}{d_k}/({d_k} \!+\! {d_j})]},$ ${p_{{\rm{e}},k}} = \sum\nolimits_j {[{F_j}{d_j}/({d_k} + {d_j})]} $， ${F_j} = [{f_{{\rm{a}},j}}/d_j^{}]/\sum\nolimits_l {({f_{{\rm{a}},l}}/d_l^{})} $； U = $( {{u^2} + {v^2}})^{1/2} $. 对于太平口水道，采用张瑞瑾^[34]水流挟沙力经验公式：

(7) ${c_{{\rm{e}},k}} = K{\left( {\frac{{{U^3}}}{{gh{\omega _k}}}} \right)^m}.$

式中： $K$=0.05~0.15， $m$=0.92，本文取 $K$=0.12.

采用非结构三角网格离散计算区域，将控制方程组（1）在编号为i的三角形单元积分，可得

(8) $\int_{{A_i}} {\frac{{\partial {{U}}}}{{\partial t}}{\rm{d}}A} + \int_{{A_i}} {\nabla \bullet {{E}}{\rm{d}}A} = \int_{{A_i}} {\left( {{{{S}}_{{1}}} + {{{S}}_{{2}}}} \right){\rm{d}}A} .$

式中： ${{E}} = [{{F}},{{G}}]$； ${A_i}$ 为单元面积，下标i为单元编号， $i = 1,2,3,\cdots,N_{\rm{c}}$ 为单元编号, $N_{\rm c}$ 为单元总数. 记 ${{{U}}_i} = \int_{A_{{i}}} {{{U}}{\rm{d}}A} /{A_i},\;$ ${{{S}}_{i1}} = \int_{A_i} {{{{S}}_1}{\rm{d}}A} /{A_{{i}}},\;$ ${{{S}}_{i2}} = \int_{A_{{i}}} {{{{S}}_2}{\rm{d}}A} /{A_i},$对式（8）左边第二项应用高斯定理，将面积分转换为线积分，可得

(9) ${A_i}\frac{{\partial {{{U}}_i}}}{{\partial t}} + \oint\limits_\varGamma {{{{E}}_n}({{U}})} {\rm{d}}\varGamma = {A_i}({{{S}}_{i1}} + {{{S}}_{i2}}).$

式中： ${{{E}}_{{n}}}({{U}}) = {{F}}{{{n}}_x} + {{G}}{{{n}}_y}$，其中 ${{{n}}_x}$ 和 ${{{n}}_y}$ 为积分边界的法向和切向向量. 考虑到三角形由三条边组成，改写式（9）中左边第2项的线积分，可得

(10) $\frac{{\partial {{{U}}_i}}}{{\partial t}} + \frac{1}{{{A_i}}}\sum\limits_{j = 1}^3 {{{{E}}_{nij}}\Delta {L_{ij}}} = {{{S}}_{i1}} + {{{S}}_{i2}};i = 1,2,3,\cdots,N_{\rm{c}}.$

式中： $\Delta {L_{ij}}$ 为单元 $i$ 第 $j$ 条边的边长, 下标 $j = 1,2,3$ 表示3条边的编号； ${{{E}}_{nij}}$ 为通过第 $j$ 条边的对流通量，由HLLC近似黎曼算子估算^[2-3]；床面底坡采用Slope flux method^[2]. 采用显式方法离散式（10）左边第1项，可得

(11) ${{U}}_i^{n + 1} = {{U}}_i^n - \frac{{\Delta {t_{\max }}}}{{{A_i}}}\sum\limits_{j = 1}^3 {{{{E}}_{nij}}\Delta {L_{ij}}} + \Delta {t_{\max }}({{{S}}_{i1}} + {{{S}}_{i2}}).$

式中：上标 $n + 1$ 和 $n$ 分别表示2个相邻时间层， $\Delta {t_{\max }}$ 为最大时间步长. $\Delta {t_{\max }}$ 取值为所有三角形网格局部时间步长的最大值（详见1.3.2节）.

动床水沙耦合数学模型的局部时间步长技术的实现分为以下2步.

1）重构每个单元的局部时间步长 $\Delta {t_i}$. 先根据式（12）计算满足CFL条件的 $\Delta t_i^{{\rm{CFL}}}$：

(12) $ \Delta t_i^{{\rm{CFL}}} = Cr{\min _{j = 1,2,3}}\left( {\frac{{{R_{i,j}}}}{{\sqrt {u_i^2 + v_i^2} + \sqrt {g{h_i}} }}} \right). $

式中： ${R_{i,j}}$ 为单元中心到3条边的距离， $Cr$ 为克朗数. 再计算整体最小时间步长 $\Delta {t_{\rm{g}}} = {\min _{i = 1,2,3,\cdots,N_{\rm{c}}}}$ $(\Delta t_i^{{\rm{CFL}}})$. 根据式（13）计算单元的时间步长级别：

(13) ${m_i} = \operatorname{int}\; \left( {\frac{{\log \;(\Delta t_i^{{\rm{CFL}}}/\Delta {t_{\rm{g}}})}}{{\log\; 2}}} \right).$

式中： $\operatorname{int} ()$ 为取整函数. 按式（14）修改单元级别，计算实际的局部时间步长：

(14) $\left. \begin{gathered} m_i^* = \min \;({m_{\max }},{m_i},{m_{i1}},{m_{i2}},{m_{i3}}), \hfill \\ \Delta {t_i} = {2^{m_i^*}}\Delta {t_{\text{g}}}. \hfill \\ \end{gathered} \right\}$

式中： ${m_{ij}}(j = 1,2,3)$ 为与单元i相邻的3个单元的时间步长级别， ${m_{\max }}$ 为用户设定最大级别. 记重构单元时间步长级别 $m_i^*$ 的最大值为 $m_{\max }^*$. 采用 $L = m_{\max }^* + 1$ 表示局部时间步长算法的级数（ $L = 1$ 代表传统GTS方法）. 式（11）中 $\Delta {t_{\max }} = {2^{m_{\max }^*}}\Delta {t_{\rm{g}}}$.

2）更新单元水沙变量. 每个单元在一个 $\Delta {t_{\max }}$ 内的更新次数为 ${N_{\rm{p}}} = {2^{m_{_{\max }}^{\rm{*}} - m_{_i}^{\rm{*}}}}$. 在一个 $\Delta {t_{\max }}$ 内，需要对地形冲淤进行累加，实现底床冲淤变化的更新. 首先计算2个时间层之间的临时守恒量 ${{U}}_i^*$：

(15) ${{U}}_i^* = {{U}}_i^n - \frac{{\Delta {t_i}}}{{{A_i}}}\sum\limits_{j = 1}^3 {{{{E}}_{nij}}\Delta {L_{ij}}} + \Delta {t_i}{{{S}}_{i1}}.$

根据 ${{U}}_i^*$，计算阻力源项向量 ${{{S}}_{i2}}$，从而得到下一时间层守恒量 ${{U}}_i^{n + 1}$：

(16) ${{U}}_i^{n + 1} = {{U}}_i^* + \Delta {t_i}{{{S}}_{i2}}.$

将式（15）、（16）循环 ${N_{\rm{p}}}$ 次，使所有单元变量更新至同一时间层. 更新变量包括每个单元在 $\Delta {t_i}$ 内地形冲淤的累加，具体如下：

(17) ${{U}}_i^{n + 1} = {{U}}_i^n - \sum\limits^{{N_{\rm{p}}}}_{l=1} {\left[ {\frac{{\Delta {t_i}}}{{{A_i}}}\sum\limits_{j = 1}^3 {{{{E}}_{nij,l}}\Delta {L_{ij}}} {\rm{ - }}\Delta {t_i}({{{S}}_{i1,l}} + {{{S}}_{i2,l}})} \right]} .$

式（17）实现了式（11）的动床局部时间步长计算方法. 如图1所示为局部时间步长技术的实现流程示意图.

为了检验局部时间步长算法的计算精度和计算效率，采用2组算例验证. 第1组算例（2.1节）复演Sanders^[22]针对清水定床模型的局部时间步长算法的结果；第2组算例（2.2节）采用UCL大学动床溃坝水流实验. 为了对比局部时间步长与整体时间步长给计算结果带来的差异，采用两者的均方根误差作为定量指标，定义如下：

(18) ${L_2}(f) = \sqrt {\left( {\sum\limits_{i = 1}^{N_{\rm{c}}} {\frac{1}{{N_{\rm{c}}}}} {{\left( {{f_{i,L}} - {f_{i,1}}} \right)}^2}} \right)} .$

式中： ${L_2}(f)$ 为变量 $f$（如 $h$、 $u$、水位 $z$、底床高程变化量 $\Delta {z_{\rm b}}$ 等变量）的均方根误差，下标 $L$ 表示得到相应计算结果所采用的级数. 采用相对计算总时间来衡量不同 $L$ 下对计算效率的提高：

(19) ${T_{\rm{r}}} = {T_L}/{T_1}.$

式中： ${T_{\rm{r}}}$ 为相对计算时间，数值上等于第 $L$ 级的局部时间步长算例计算的总时间 ${T_L}$ 与GTS算例计算总时间 ${T_1}$ 的比值.

溃坝定床算例采用理想的溃坝过程. 采用1 000 m× 100 m的矩形计算区域，共计3 802个计算单元，平均三角单元网格边长为7.8 m. $x \leqslant 500$ m为坝址上游区域，设置初始水深为1 m. 该组算例包括如下3个工况. 工况1：无摩阻、下游为干床面；工况2：有摩阻（底床曼宁糙率系数 $n = 0.02$）、下游为干床面；工况3：无摩阻、下游为湿床面（下游初始水深为0.1 m）.

计算结果如图2所示. 如图2（a）所示为无摩阻工况. 可以看出，当 $L = 1$ 或2时，计算结果满足要求；当 $L \geqslant 3$ 时，数值结果与解析解在溃坝前锋处产生明显的差异. 这是由于在溃坝前锋流速大，所得到的浅水波速大，前锋界面流速大的单元时间步长小，另一侧单元时间步长大，两者差异明显，导致计算不稳定. 如图2（b）所示为有摩阻工况. 可以看出，当 $L \leqslant 4$ 时，计算结果满足要求. 这是由于溃坝水流前锋由于阻力的存在，浅水波速大大减小，使得大的级数条件下计算保持稳定. 如图2（c）所示为湿床面工况，在 $L \leqslant 4$ 时均能够满足要求.

如表1所示为清水定床算例组在不同 $L$ 下计算时间和计算精度的定量对比. 无摩阻工况表明， $L = 2$ 时的计算结果与GTS的计算结果误差较小；当 $L \geqslant 3$ 时，误差过大. 对于有摩阻工况和湿床面工况， $L$ 均能够在保证精度的情况下达到4，节省计算时间分别为73%和64%. 这表明该模型成功再现了Sanders^[22]的局部时间步长技术.

溃坝动床算例采用比利时UCL大学实验数据^[35-36]. 水槽平面如图3所示，水槽长6 m，宽度在距离左边界4 m的地方展宽，离左侧3 m位置设置水闸；床面铺厚度为0.1 m的泥沙，中值粒径 ${d_{50}} = $ 1.82 mm； ${\rho _{\rm s}}$=2 680 kg/m³； $p$=0.47， $n$=0.024^[34-35]. 采用式（6）计算水流挟沙力， $\alpha = \min \;(10,h/(\omega \Delta {t_i}))$. 初始坝左侧水深 为0.25 m，右侧为干床. 在6个测点（P1~P6，见图3）观测水位；在2个断面（CS1；CS2)测量冲淤. 网格有7 308个三角单元.

如图4所示为P1～P6测点实测和计算水位随时间的变化. 当 $L \leqslant 3$ 时，计算结果与实测吻合很好；当 $L = 4$ 时，P1、P2和P3 3个测点在1 s左右产生了极大值. 这与图2（a）类似，即较高级数（如 $L = 4$）条件下，溃坝前锋两侧浅水波速的剧烈变化会导致计算不稳定. 如图5所示为CS1和CS2断面计算和实测冲淤地形的比较. 结果表明，当 $L \leqslant 3$ 时，计算模拟的结果与实测吻合很好；当 $L = 4$ 时，CS1和CS2的地形冲淤结果与GTS结果产生误差.

如表2所示为不同 $L$ 下计算时间和计算精度的对比. 当 $L \leqslant 3$ 时，计算结果与GTS结果相比误差较小；当 $L = 4$ 时，计算所得的 ${L_2}(\Delta {z_{\rm{b}}})$ 明显增大，这与图3、4一致. $L = 2$ 和 $L = 3$ 能够分别节省47%和68%的计算时间.

将模型应用于长江中游的太平口水道的地形冲淤复演. 如图6所示，太平口水道位于荆州市南面，长约为20 km. 选取太平口水道2012年2月的初始地形，采用2012年2月—2012年10月的来水来沙过程作为上游边界条件（见图7（a）），2012年的下游水位-流量关系作为下游边界条件（见图7（b）），对2012年2月17日—2012年10月5日太平口水道经历一个洪期的地形冲淤演变进行复演，共计231天. 图7中，q_V为体积流量，z为水位. 最大和最小网格边长分别为167和53 m. 采用统一的底床糙率 $n = 0.02$，水流挟沙力采用式（7）.

如图8所示为不同级数下2012年10月5日计算水深分布对比图. 可以看出，在不同级数下，水深分布基本一致. 如图9所示为不同L下2012年2月17日—2012年10月5日计算地形冲淤分布对比图. 对比GTS与实测值计算结果可以看出，该算例复演了太平口水道2012年经历洪峰过程下的地形冲淤演变，实测冲刷和淤积位置与计算保持一致，且从冲淤的量级上来看吻合很好. 对比不同级数条件下的计算结果与GTS的结果可以看出，该模型在不同局部时间步长级数下，地形冲淤的结果图像对比上几乎完全一致. 对于实际河道（如本文太平口水道）冲淤的水沙模拟，存在较大的不确定性. 在没有大规模调整参数的情况下，该模型较好地复演了2012年2月—10月太平口水道的实际冲淤过程.

如表3所示为太平口水道算例组在不同 $L$ 下计算时间和计算精度的对比. 可知，当级数取2~6时， ${L_2}(h)$ 和 ${L_2}(\Delta {z_{\rm{b}}})$ 在同一量级，0.2 m水深误差相较于最大水深24 m约为0.8%，0.2 m地形冲淤误差相较于14 m最大冲淤厚度约为1.4%. 这一误差程度可以接受. $L = 2$ 时可以节省61%的计算时间， $L = 4$ 时可以节省79%的时间， $L = 6$ 时可以节省92%的计算时间，等同于计算效率提高12.5倍.

（1）算例验证表明，现有局部时间步长技术对于无摩阻地形上溃坝水流的模拟要取较小的级数. 这是由于溃坝水流前锋波速大，水流变化剧烈，前锋两侧单元的局部时间步长级数相差较大，导致计算不稳定. 当地形有摩阻时，溃坝前锋遇到阻力，水流变化减弱，可取的局部时间步长级数可以相应增大. 对于UCL大学的溃坝动床算例，该模型在保证精度的前提下，能够节省68%的计算时间. 如何处理水流变化剧烈界面两侧单元的局部时间步长是需要进一步研究的内容.

（2）模型在太平口水道的应用表明，如果水流交替变化不剧烈，较大的局部时间步长级数不会降低计算精度. 本文模拟了太平口水道近7个月的地形冲淤过程，取局部时间步长级数为6，能够节省92%的计算时间（加速12.5倍）.

除局部时间步长技术外，还可以通过并行计算技术（如Open MP、MPI、GPU等）进一步提高计算效率. 本文模型已采用适应于个人计算机的Open MP并行技术.