数值解法

二维对流扩散方程的求解，依托于二维水动力模型的求解出的流场，采用有限体积法进行求解。

有限体积离散

在空间上，在任意有限控制体上进行积分得：
$$\tag{1} \frac{\partial}{\partial t} \int_{\Omega}hcd\Omega+\int_{\Omega}(\frac{\partial huc}{\partial x}+\frac{\partial hvc}{\partial x})d\Omega =\int_\Omega(\frac{\partial}{\partial x}(hE_x\frac{\partial c}{\partial x})+\frac{\partial }{\partial y}(hE_y\frac{\partial c}{\partial y}))d\Omega +\int_\Omega sd\Omega$$

根据格林公式将面积分转为线积分：
$$\tag{2} \frac{\partial}{\partial t} \int_{\Omega}hcd\Omega+\int_{\partial \Omega}[huc,hvc]\cdot {\bf{n}} d\ell =\int_{\partial \Omega}([hE_x\frac{\partial c}{\partial x},hE_y\frac{\partial c}{\partial y}]\cdot {\bf{n}} d\ell +\int_\Omega sd\Omega$$

其中$\partial \Omega$为控制体的边界，${\bf{n}}$为边界的外法向量，$d\Omega$、$d\ell$分别为面积微元和弧微元。针对采用三角形网格系统的二维水动力模型，每个控制体为三角形，线积分有：
$$\tag{3} \int_{\partial \Omega}[huc,hvc]\cdot {\bf{n}} d\ell=\sum_{k=1}^3[huc,hvc]\cdot[n_{k,x},n_{k,y}]L_k$$ $$\tag{4} \int_{\partial \Omega}([hE_x\frac{\partial c}{\partial x},hE_y\frac{\partial c}{\partial y}]\cdot {\bf{n}} d\ell =\sum_{k=1}^3[hE_x\frac{\partial c}{\partial x},hE_y\frac{\partial c}{\partial y}]\cdot[n_{k,x},n_{k,y}]L_k$$
其中$L_k$为第条边的长度，$[n_{k,x},n_{k,y}]$为第条边的外法向量。

令$Q_c$为在单元$C_i$内的平均值，即：
$$\tag{5} Q_c=\frac 1{\Omega_i}\int_{C_i}hcd\Omega_i$$ 则可以得到半离散式常微方程：

$$\tag{6} \Omega_i\frac{dQ_c}{dt}=-\sum_{k=1}^3[huc,hvc]\cdot[n_{k,x},n_{k,y}]L_k+\sum_{k=1}^3[hE_x\frac{\partial c}{\partial x},hE_y\frac{\partial c}{\partial y}]\cdot[n_{k,x},n_{k,y}]L_k+S_c$$

其中$\Omega_i$为单元$C_i$的面积,$S_c$为源项。

在实际计算时，单元每条边的对质量流通量已经由二维潜水方程采用 HLLC 黎曼求解器求出，水质的通量采用一阶迎风格式：
$$\tag{7} \text{Flux}C=\text{flux}_0\begin{cases} c_L& \text{if} \ \text{flux}_0>0.0\\ c_R& \text{else} \end{cases}$$ 计算，其中$c_L$、$c_R$分别为边左右水质浓度的重构值。重构的方法同流速的重构。

扩散通量采用：
$$\tag{8} F_{diff}=\frac 12(h_L+h_R)(E_x\partial_x c n_x+E_y\partial_y c n_y)$$

其中分$h_L$、$h_R$别为边左右水深的重构值，$\partial_x c$、$\partial_y c$分别为边左右水质浓度的原始梯度。

时间二阶积分

在时间上，采用 Hancock 预估校正法，在采用合理的时间步长的条件下可以稳定运行，虽然不具有采用四阶龙格库塔法具有 TVD 性质的优点，但只需要计算一次黎曼问题，效率能够大大提高。

1.预测步

在预测步，计算区域内的水质状态由$t$时刻更新到$t+\frac{\Delta t}{2}$时刻，其中$\Delta t$为时间步长。 若单元为全淹没状态，浓度的预测量为：
$$\tag{9} c_{i}^{t+\frac{\Delta t}{2}}=c_i^{\Delta t}+\frac{\Delta t}{2}(u\partial_xc+v\partial_yc)|_i^t$$ 否则采用当前时刻的浓度作为预测量。

2.校正步

在校正步，水质状态由$t$时刻更新到$t+\Delta t$时刻，以$t+\frac{\Delta t}{2}$时刻水流、水质的状态计算式(5.1-7)右端的通量，时间上采用一阶差分向前差分，即：
$$\tag{10} (hc)_i^{t+\Delta t} =(hc)_i^{t} +\frac{\Delta t}{\Omega_i} \{ -\sum_{k=1}^3[huc,hvc]\cdot[n_{k,x},n_{k,y}]L_k+\sum_{k=1}^3[hE_x\frac{\partial c}{\partial x},hE_y\frac{\partial c}{\partial y}]\cdot[n_{k,x},n_{k,y}]L_k+S_c \}_i^{t+\frac{\Delta t}{2}}$$