数值计算方法

采用边界拟合能力强和易于局部网格加密的三角形网格剖分计算域，利用基于水位-体积关系的斜底三角单元模型，有效解决了小尺度线状地形模拟难题；以能够有效捕获激波的 Godunov 型有限体积法为框架，运用在时间上和空间上均具有二阶精度的 MUSCL-Hancock 预测-校正格式离散洪水控制方程，采用 HLLC 近似 Riemann 算子计算对流数值通量，采用直接近似方法计算扩散数值通量，并结合斜率限制器以保证模型的高分辨率特性，避免在间断或大梯度解附近产生非物理虚假振荡；基于单元中心型底坡项近似方法，在不使用任何额外动量通量校正项的前提下模型能保持通量梯度与底坡项之间的平衡，即模型具有和谐性质；采用半隐式格式处理摩阻项，该半隐式格式既能保证不改变流速分量的方向，也能避免小水深引起的非物理大流速问题，有利于计算稳定性；实现了固壁、水位、流量、自由出流等边界条件；基于 CFL 稳定条件实现了数值模型的自适应时间步长技术。

• 数值计算方法
  • 计算网格
  • 斜底三角单元
  • 有限体积离散
  • 数值通量计算
  • 变量重构和斜率限制函数
  • 底坡项近似
  • 摩阻项处理
  • 柯氏力项处理
  • 时间二阶积分
    • （1）预测步
    • （2）校正步
  • 边界条件
    • （1）急流开边界条件
    • （2）缓流开边界条件
  • 稳定条件

计算网格

鉴于三角形单元具有复杂边界拟合能力强、便于网格生成和局部加密等特点，本模型采用非结构三角形单元作为计算网格，如图3.2-1所示。其中，为待计算单元，其顶点1-顶点2-顶点3排序服从逆时针方向；与顶点相对的边为，其外法向单位向量为；单元的邻接单元中，与顶点相对的单元为。

对于每一个节点而言，其网格拓扑信息包括：节点序号，节点坐标，节点周围的单元；对于每一条边而言，其网格拓扑信息包括：边序号，边的左、右单元，边的始、末节点；对于每一个单元而言，其网格拓扑信息包括：单元序号，单元的三个顶点，单元的三条边，单元的三个邻接单元。在模型开始计算之前，需要根据网格文件构造包含上述拓扑信息的计算网格系统。

斜底三角单元

本模型采用斜底三角单元模型，即将底高程定义于单元顶点处。这种方式在地形表达精度方面为二阶精度。

模型计算时，水深为守恒变量，表征单元的水量，即 $V_i=h_i\Omega_i$，其中 $V_i$ 为水量。图2中，斜线阴影面为水面，三角形 $\Delta123$ 为单元底面，（ａ）中3个顶点的水深均大于零，（ｂ）和（ｃ）存在水深为零的顶点。假设底高程满足 $b_i^{n,1} \leqslant b_i^{n,2} \leqslant b_i^{n,3}$，基于此，单元水位 $z_i$ 与水量存在如下关系： $$\tag{1} V_i/\Omega_i = h_i = \begin{cases} (z_i-b_i^{n,1})^3/[3(b_i^{n,2}-b_i^{n,1})(b_i^{n,3}-b_i^{n,1})] &\text{if} \ b_i^{n,1} < z_i \leqslant b_i^{n,2} \\ (z_i^2+z_ib_i^{n,3}-3z_ib_i^{n,1}-b_i^{n,2}b_i^{n,3}+b_i^{n,1}b_i^{n,2}+b_i^{n,1}b_i^{n,1})/[3(b_i^{n,3}-b_i^{n,1})] &\text{if} \ b_i^{n,2} < z_i \leqslant b_i^{n,3} \\ (z_i-b_i) &\text{if} \ z_i > b_i^{n,3} \end{cases}$$ 式中 $b_i$ 为单元形心处的底高程，$b_i=(b_i^{n,1}+b_i^{n,2}+b_i^{n,3})/3\geqslant b_i^{n,1}$。模型初始化时，若以水位为初始条件，则通过式（1）可计算出单元水深值；模型计算过程中，对单元水深值进行更新后，通过式（1）可计算出相应的单元水位值。与平底单元相比，基于水位～体积关系建立的斜底单元模型，使单元水位的值域由 $[b_i, +\infty)$ 扩展为 $[b_i^{n,1}, +\infty)$，避免了干湿界面计算时重构干单元底高程，提高了模型的动边界处理能力，有利于设计具有和谐性的计算格式。

有限体积离散

在任意控制体上对式（1）所示的控制方程进行积分得： $$\tag{2} \frac{\partial}{\partial t}\int_\Omega{\bf{U}d\Omega}+ \int_\Omega(\frac{\partial{\bf{E}^{adv}}}{\partial x}+\frac{\partial{\bf{G}^{adv}}}{\partial y})d\Omega= \int_\Omega(\frac{\partial{\bf{E}^{diff}}}{\partial x}+\frac{\partial{\bf{G}^{diff}}}{\partial y})d\Omega + \int_\Omega{\bf{S}d\Omega}$$

根据格林公式，将上式转化为沿其边界的线积分，可得： $$\tag{3} \frac{\partial}{\partial t}\int_\Omega{\bf{U}d\Omega}+ \oint_{\partial \Omega}{\bf{F}^{adv} \cdot \bf{n}}d\ell = \oint_{\partial \Omega}{\bf{F}^{diff} \cdot \bf{n}}d\ell+ \int_\Omega{\bf{S}d\Omega}$$ 式中，$\partial \Omega$为控制体$\Omega$的边界；$n$为边界的外法向单位向量；$d\Omega$、$d\ell$分别为面积微元和弧微元；$F^{adv}= [E^{adv},G^{adv}]$;$F^{diff}= [E^{diff},G^{diff}]$。

在平面二维三角形网格中，线积分可分别由下式计算： $$\tag{4} \oint_{\partial \Omega}{\bf{F}^{adv} \cdot \bf{n}}d\ell = \sum_1^3 \bf{F}^{adv} \cdot \bf{n}_k L_k$$ $$\tag{5} \oint_{\partial \Omega}{\bf{F}^{diff} \cdot \bf{n}}d\ell = \sum_1^3 \bf{F}^{diff} \cdot \bf{n}_k L_k$$ 假设$U_i$为$U$在单元$C_i$内的平均值，即： $$\tag{6} \bf{U}_i = \frac{1}{\Omega_i}\int_{C_i}\bf{U}d\Omega$$

则有： $$\tag{7} \Omega_i \frac{d \bf{U}_i}{dt}= -\sum_{k=1}^{k=3} \bf{F}_{i,k}^{adv} \cdot \bf{n}_{i,k} L_{i,k} +\sum_{k=1}^{k=3} \bf{F}_{i,k}^{diff} \cdot \bf{n}_{i,k} L_{i,k} +\bf{S}_i$$ 式中，$\Omega_i$为单元$C_i$的面积；$F_{i,k}^{adv}$ 、$F_{i,k}^{diff}$、$n_{i,k}$、$L_{i,k}$分别代表单元第条边的对流数值通量、扩散数值通量、外法向单位向量和长度；$S_i$为源项近似。

数值通量计算

由于积分平均，物理变量在每个单元内部为常数，在整个计算域内形成阶梯状分布，因此在单元界面处物理量存在间断，即界面左、右两侧的物理量不相等，故而在界面处构成了一个局部 Riemann 问题。通过 Riemann 问题的求解可得到界面处的对流数值通量，即：

$$\tag{8} \bf{F}^{adv} \cdot \bf{n} = [\bf{E}^{adv},\bf{G}^{adv}]\cdot [n_x,n_y] = \bf{E}^{adv}n_x +\bf{G}^{adv}n_y =\begin{bmatrix} h(un_x+vny) \\ hu(un_x+vn_y)+\frac 12 g(h^2-b^2)n_x \\ hv(un_x+vn_y)+\frac12g(h^2-b^2)n_y \end{bmatrix}$$ 令 $$\tag{9} u_\perp = un_x+vn_y \\ u_{//} = -un_y+vn_x$$ 式中，$u_\perp$、$u_{//}$分别为与界面垂直和平行的流速分量；定义： $$\tag{10} \hat{\bf{U}} = \bf{T} (\bf{n}) \cdot \bf{U} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & n_x & n_y \\ 0 & -n_y & n_x \end{bmatrix} \cdot\bf{U} =\begin{bmatrix} h \\ hu_\perp \\ hu_{//} \end{bmatrix}$$ 代入(12)可得： $$\tag{11} \bf{F}^{adv} \cdot \bf{n} =\begin{bmatrix} hu_\perp \\ huu_\perp +\frac 12 g(h^2-b^2)n_x \\ hvu_\perp +\frac12g(h^2-b^2)n_y \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ 0 & n_x & -n_y \\ 0 & n_y & n_x \end{bmatrix} \bf{E}( \hat{\bf{U}} )$$ 可知，二维浅水方程的对流数值通量计算可转化为界面处一维 Riemann 问题求解。

如图2所示，在平面上，二维浅水方程 Riemann 解中存在三个波，分为三类：激波（shock wave）、稀疏波(rarefaction wave)和接触波（contact wave）。三个波将初始的两个常状态区分割为四个常状态区，其中左、右两个区域的状态分别为初始的左、右状态值，中间两个区域的状态则由 Riemann 解的结构确定。

$$\tag{12} \bf{F}^{adv}(\bf{U}_L,\bf{U}_R) \cdot \bf{n}= \begin{cases} \bf{F}_L^{adv} &\text{if} \ s_1 \geqslant 0 \\ \bf{F}_{*L}^{adv} &\text{if} \ s_1 <0 \leqslant s_2 \\ \bf{F}_{*R}^{adv} &\text{if} \ s_2 <0 < s_3 \\ \bf{F}_R^{adv} &\text{if} \ s_3 \leqslant 0 \end{cases}$$ 式中，$\bf{F}_L^{adv} = \bf{F}^{adv}(\bf{U}_L)\cdot \bf{n}$、$\bf{F}_R^{adv} = \bf{F}^{adv}(\bf{U}_R)\cdot \bf{n}$ ，均有式(3.2-15)计算，其中$\bf{U}_L$、$\bf{U}_R$分别为局部 Riemann 问题所在边界左侧与右侧的守恒变量；$\bf{F}_{*L}^{adv}$、$\bf{F}_{*R}^{adv}$分别为局部 Riemann 解中间区域接触波左右侧的数值通量；$s_1$、$s_2$、$s_3$分别为左波、接触波和右波的波速。 $$\tag{13} \bf{F}_{*L}^{adv}=\begin{bmatrix} (\bf{E}_{HLL}^{adv})^1 \\ (\bf{E}_{HLL}^{adv})^2n_x-u_{//,L}(\bf{E}_{HLL}^{adv})^1n_y\\ (\bf{E}_{HLL}^{adv})^2n_y+u_{//,L}(\bf{E}_{HLL}^{adv})^1n_x \end{bmatrix}$$ $$\tag{14} \bf{F}_{*R}^{adv}=\begin{bmatrix} (\bf{E}_{HLL}^{adv})^1 \\ (\bf{E}_{HLL}^{adv})^2n_x-u_{//,R}(\bf{E}_{HLL}^{adv})^1n_y\\ (\bf{E}_{HLL}^{adv})^2n_y+u_{//,R}(\bf{E}_{HLL}^{adv})^1n_x \end{bmatrix}$$

式中$(E_{HLL}^{adv})^1$、$(E_{HLL}^{adv})^2$分别为运用 HLL 格式计算得到的法向数值通量的第一、二个分量： $$\tag{15} \bf{E}_{HLL}^{adv}=\frac{s_3 \bf{E}^{adv}(\hat{\bf{U}}_L) -s_1\bf{E}^{adv}(\hat{\bf{U}}_R) +s_1s_3((\hat{\bf{U}}_R)-(\hat{\bf{U}}_L))}{s_3-s_1}$$ 本模型采用 Einfeldt 波速近似计算波速： $$\tag{16} s_1 = \begin{cases} min(u_{\perp,L}-\sqrt{gh_L},u_{\perp ,*}-\sqrt{gh_{*}}) &\text{if} \ h_L>0 \\ u_{\perp,R}-2\sqrt{gh_R} &\text{ if} \ h_L=0 \end{cases}$$

$$\tag{17} s_3 = \begin{cases} min(u_{\perp,L}+\sqrt{gh_R},u_{\perp ,*}+\sqrt{gh_{*}}) &\text{if} \ h_R>0 \\ u_{\perp,L}-2\sqrt{gh_L} &\text{ if} \ h_R=0 \end{cases}$$ 其中$h_*$、$u_{\perp,*}$为 Roe 平均: $$\tag{18} h_* =\frac 12(h_L+h_R)$$ $$\tag{19} u_{\perp,*}=\frac{\sqrt{h_L}u_{\perp,L}+\sqrt{h_R}u_{\perp,R}}{\sqrt{h_L}+\sqrt{h_R}}$$

由于激波的波速小于激波后面区域的特征波速，此时 Einfeldt 波速即为激波速的 Roe 近似，因此使用 Einfeldt 波速可以获得在激波附近更加准确的数值解。接触波波速： $$\tag{20} s_2=\frac{s_1h_R(u_{\perp,R}-s_3)-s_3h_L(u_{\perp,L}-s_1)}{h_R(u_{\perp,R}-s_3)-h_L(u_{\perp,L}-s_1)}$$

扩散数值通量采用下式计算： $$\tag{21} \bf{F}_{i,k}^{diff} \cdot \bf{n}_{i,k} = \begin{bmatrix} 0 \\ h_{i,k} \nu_t [2 \frac{\partial u}{\partial x}|_{i,k} n_{i,k}^{x}+(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})|_{i,k} n_{i,k}^{y}] \\ h_{i,k} \nu_t [(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})|_{i,k} n_{i,k}^{x}+2\frac{\partial v}{\partial y}|_{i,k} n_{i,k}^{y} ] \end{bmatrix}$$ 式中，$h_{i,k}$为单元$C_i$第$k$条边所在界面处的水深；$\frac{\partial u}{\partial x}|_{i,k}$、$\frac{\partial u}{\partial y}|_{i,k}$、$\frac{\partial v}{\partial x}|_{i,k}$、$\frac{\partial v}{\partial y}|_{i,k}$为流速分量在单元$C_i$第$k$条边所在界面处的斜率。

变量重构和斜率限制函数

如果物理量在单元内近似为常数分布，并且以各单元形心处的值作为界面处局部 Riemann 问题的初始间断条件，则相应的计算格式在空间上仅具有一阶精度，存在较大的数值耗散。对水流模拟而言，一阶精度的计算格式能基本满足工程应用的要求。然而，对泥沙运动或者污染物扩散模拟而言，一阶格式的过度数值耗散往往导致计算结果失真。为了提高格式的空间精度，在构造界面处局部 Riemann 问题时，需要采用比分片常数近似函数更高精度的重构方法对界面左、右两侧的变量进行重构，并基于重构变量求解界面处的数值通量。本模型针对不同的单元干湿状态，采用不同的变量重构方法。

若单元处于局部淹没状态，则采用采用该单元形心处的值为重构值，否则若单元处于全淹状态，则采用线性重构并结合梯度限制器计算水位、流速的重构值： $$\tag{22} p(x,y)=p_i+\overline{\nabla p_i}\cdot \bf{r}$$ 式中，$p$为重构变量，如水位、流速；$\bf{r}$为界面中点相对于单元形心处的位置矢量；$\overline{\nabla p_i}$为限制梯度： $$\tag{23} \overline{\nabla p_i}=\varphi_i \nabla p_i$$ 其中$\varphi_i$为限制函数： $$\tag{24} \varphi_i = \min_{k=1,2,3}(\varphi_k) \ \ \varphi_k= \begin{cases} min(1,\frac{p_i^{max}-p_i}{p_{i,k}^{Rec}-p_i} )&\text{if} \ p_{i,k}^{Rec}-p_i>0 \\ min(1,\frac{p_i^{min}-p_i}{p_{i,k}^{Rec}-p_i} )&\text{if} \ p_{i,k}^{Rec}-p_i<0 \\ 1 &\text{if} \ p_{i,k}^{Rec}-p_i=0 \end{cases}$$ 其中$p_{i,k}^{Rec}$为单元顶点处未受限制的重构值，$p_i^{max}=\max p_{i,k}$,$p_i^{min}= \min p_{i,k} ,\ k=1,2,3$。

底坡项近似

在每一个计算步，本模型均对所有单元进行连续方程和动量方程的更新，即同时更新单元的水深和流速。因此，需要对所有单元进行底坡项处理。 若单元为全淹没状态，本模型采用单元中心型近似方法处理底坡项： $$\tag{25} S_{i,0x}=-\int_{C_1}g(h+b)\frac{\partial b}{\partial x}d\Omega=-g(h+b)\Omega\frac{\partial b}{\partial x}|_i$$

否则，本模型采用 DFB（Divergence Form for Bed slope source term）技术处理底坡项： $$\tag{26} S_{i,0x}= \sum_{k=1}^3\frac12g(h_{i,k}^{Rec})^2L_{i,k}n_{i,k}^x- \sum_{k=1}^3\frac12g(b_{i,k})^2L_{i,k}n_{i,k}^x$$

式中，$b_{i,k}$为单元$C_i$第$k$条边的中点位置底高程。 采用上述底坡项近似方法，保证了模型的和谐性。

摩阻项处理

一般地，摩阻项通过算子分裂法进行处理： $$\tag{27} \frac{d\bf{U}}{dt}=S_f(\bf{U}) \implies \frac {d}{dt} \begin{bmatrix} h \\ hu \\ hv \end{bmatrix} = \begin{bmatrix} 0 \\ {-gn^2hu\sqrt{u^2+v^2}}/{h^{4/3}} \\ {-gn^2hv\sqrt{u^2+v^2}}/{h^{4/3}} \end{bmatrix}$$

单元水深保持不变 式（31）可简化为： $$\tag{28} \frac {d}{dt} \begin{bmatrix} u \\ v \end{bmatrix} =-gn^2 h^{-{4}/{3}} \begin{bmatrix} u\sqrt{u^2+v^2} \\ v\sqrt{u^2+v^2} \end{bmatrix}$$ 复杂地形的陡峭坡面，使局部区域的水深较小、流速较大，使得式（3.2-32）对应的常微分方程系统的 Lipschitz 常数很大，因此，摩阻项可能引起刚性问题。此时，若采用一般的显式数值方法处理摩阻项，将显著影响数值计算的稳定性，或将极大减小全局时间步长，从而严重降低计算效率。为解决上述问题，需要采用隐式或半隐式格式处理摩阻项。然而，由于水深变量位于摩阻项的分母，一般的隐式或半隐式计算格式仍面临一些问题，如产生错误的大流速、改变流速分量的方向等。

综合考虑摩阻项处理的稳定性和计算效率，本模型采用如下半隐式格式处理摩阻项:

$$\tag{29} \frac{u^{n+1}-u^{n}}{\Delta t}=(-gn^2\sqrt{u^2+v^2} h^{-{4}/{3}})|^nu^{n+1}$$

柯氏力项处理

与魔阻项类似，柯氏力项通过算子分裂法进行处理： $$\tag{30} \frac{\partial\bf{U}}{\partial t}=S_k(\bf{U}) \implies \frac {\partial}{\partial t} \begin{bmatrix} h \\ hu \\ hv \end{bmatrix} = \begin{bmatrix} 0\\ fhv\\ -fhu \end{bmatrix}$$

由于在柯氏力项处理过程中， 式（34）可简化为： $$\tag{31} \frac {\partial}{\partial t} \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} fv \\ -fu \end{bmatrix}$$ 采用如下隐式差分格式：

$$\tag{32} \frac{u^{n+1}-u^{n}}{\Delta t}=f^nv^{n+1}\\ \frac{v^{n+1}-v^{n}}{\Delta t}=f^nu^{n+1}$$

时间二阶积分

在二维浅水方程数值求解过程中，时间上离散可采用 Runge-Kutta、Hancock 预测-校正等多步格式，以提高模型的时间精度。从计算效率看，在一个时间步长内，两步 Runge-Kutta 法分别对所有单元界面计算两次黎曼问题，而 Hancock 预测-校正法只需要计算一次黎曼问题，因此，Hancock 预测-校正法的计算效率相对较高；从格式稳定性看，Runge-Kutta 法可满足 TVD 性质，而 MUSCL-Hancock 法在合理选择时间步长的前提下满足 L1 稳定。考虑模型的计算效率，本模型采用 Hancock 预测-校正格式实现二维浅水方程数值求解的时间二阶积分。

（1）预测步

在预测步，计算域的水流状态由$t$时刻更新至$t+\frac 12 \Delta t$时刻，其中$\Delta t$为计算时间步长。$\Delta t$可以设定为常数，也可以根据稳定条件进行自适应调整。预测步不考虑紊动扩散项。

若单元为全干或局部淹没状态，则以当前时刻的水位、流速作为预测步的结果；若单元为全淹没状态，则以水位、流速作为预测变量： $$\tag{33} \eta _i^{t+{\Delta t} /{2}} =\eta _i^t \frac{\Delta t}{2}( h\overline{\partial _x}u +h\overline{\partial _y}v +u\overline{\partial _x}h +v\overline{\partial _y}h)|_i^t$$ $$\tag{34} u_i^{t+{\Delta t} /{2}} =\frac {1}{1+gn^2h^{-{4}/{3}}\sqrt{u^2+v^2}{\Delta t}/{2}}[ u -\frac {\Delta t}{2}(u\overline{\partial _x}u+g\overline{\partial_x}\eta+v\overline{\partial_y}u)] |_i^t$$ $$\tag{35} v_i^{t+{\Delta t} /{2}} =\frac {1}{1+gn^2h^{-{4}/{3}}\sqrt{u^2+v^2}{\Delta t}/{2}}[ v -\frac {\Delta t}{2}(u\overline{\partial _x}v+g\overline{\partial_y}\eta+v\overline{\partial_y}v)] |_i^t$$ 由于静水条件下流速为零、水位为常数，因此，由式（3.2-33）~（3.2-35）可知，预测步的水流静止状态得以维持。

（2）校正步

在校正步，计算域的水流状态由$t$时刻更新至$t+ \Delta t$时刻。基于守恒形式浅水方程，由下式对单元值进行更新：

$$\tag{36} U_i^{t+\Delta t} =U_i^t+\frac{\Delta t}{\Omega_i}[\sum_{k=1}^{3}(-F_{i,k}^{adv}+F_{i,k}^{diff}) \cdot n_{i,k}\ell_{i,k}+S_i]^{t+{\Delta t}/{2}}$$

边界条件

一般情况下，模型的边界条件实现方式有两种：镜像单元法和直接计算数值通量法。其中，前者在基于结构网格的数学模型中应用较广，而后者被广泛运用于基于非结构网格的数值模型。本模型采用直接计算数值通量的方式实现边界条件。 边界处的对流数值通量为： $$\tag{37} \bf{F}^{adv}\cdot \bf{n} = \begin{bmatrix} h_*U_{\perp,*} \\ h_*u_*u_{\perp,*}+\frac 12g(h_*^2-b^2)n_x\\ h_*v_*u_{\perp,*}+\frac 12g(h_*^2-b^2)n_y \end{bmatrix}$$ 式中，$h_*$、$u_*$、$u_{\perp,*}$分别代表边界边中点处的水深、$x$方向流速分量、界面外法向方向流速分量；$b$为边界边中点处的底高程值；$n_x$、$n_y$为分别为边界边的外法向单位向量在$x$和$y$方向分量。

由式（37）可知，计算边界边对流数值通量的关键在于估算边界边中点处的水深、流速等物理量。

（1）急流开边界条件

由浅水方程特征理论可知，水流为急流状态时信息沿下游传播，因此急流出口边界的扰动不会对计算域内的水流状态产生影响。故急流边界处的水深、流速等水力要素值采用边界边所在内单元的水力要素值： $$\tag{38} h_*=h_L$$ $$\tag{39} u_*=u_L, \ v_*=v_L$$ 在数值模拟中，常采用自由出流边界条件，其实现方式与急流开边界条件相同。

（2）缓流开边界条件

假设边界边的左侧位于计算域之内、右侧位于计算域之外，基于浅水方程特征理论，由一维浅水方程的 Riemann 不变量可得： $$\tag{40}u_{\perp,*}+2\sqrt{gh_*}=u_{\perp,L}+2\sqrt{gh_L}$$

式中，$h_L$、$u_{\perp,L}$分别为边界边所在内单元形心处的水深和界面外法向方向流速分量。 式（40）中存在两个未知变量，因此需要结合边界条件建立关于$h_*$与$u_{\perp,*}$的等式，进而通过求解方程组计算和。

对于水位边界，可根据水位求出边界上的水深$h_*$，从而可以求出流速$u_{\perp,*}$；对于流量边界，可求出边界边的单宽流量$q_{\perp,*}$，从而求出$h_*$与$u_{\perp,*}$的一个关系式：$q_{\perp,*}=h_*u_{\perp,*}$,结合式（3.2-44）可以求出$h_*$与$u_{\perp,*}$；对于固壁边界，即边界上的法向流速为零，且假设边界上的水深等于边界边所在单元形心处的水深值，即$u_{\perp,*}=0.0$， $h_*=h_L$。

稳定条件

为保持显示格式的稳定，通常采用 CFL 稳定性条件限制时间步长$\Delta t$： $$\tag{41} \Delta t=C_r\cdot \min_{i,k}[(\frac{\Omega}{(|u_\perp|+\sqrt{gh})L_k})_i] ,i=1,2,\cdots,N;k=1,2,3$$

式中，$\Delta t$为 Courant 数，$0<C_r<1$,$u_\perp$和$h$为界面的 Roe 平均；$N$为计算网格单元的个数。