CFD: 浅水方程算法

浅水方程的假设是高度方向的流动尺度要远远小于水平面的尺度。在这种情况下，可以认为高度方向的速度近似为0。不考虑粘性的浅水方程是标准的双曲偏微分方程系统。其有的时候也会被称之为Saint-Venant方程。由于浅水方程构成了一个双曲系统，因此浅水方程与欧拉方程类似，可以用于研究各类对流格式。

控制方程

浅水方程的控制方程为

(1)\[ \frac{\partial h }{\partial t}+\nabla\cdot(\bfU h) =0 \]

(2)\[ \frac{\partial h\bfU}{\partial t}+\nabla\cdot(h\bfU \bfU) =-\bfF\times\bfU-|\bfg|h\nabla(h+h_0) \]

其中\(h\)表示水深，\(h_0\)表示河底的高度，其为一个用户自定义高度，\(\bfF=(2\boldsymbol{\Omega}\cdot\bfn_g)\bfn_g\)表示科氏力，\(\bfn_g\)表示重力方向的单位矢量。方程(1)表示的是关于水深的连续性方程。方程(2)表示的是关于\(\bfU\)的动量方程。

考虑一维情况，如果忽略科氏力源项，浅水方程可以写为

(3)\[ \frac{\partial h }{\partial t}+\frac{\p hu}{\p x} =0 \]

(4)\[ \frac{\partial h u}{\partial t}+\frac{\p huu}{\p x} =-|\bfg|h\frac{\p( h+h_0)}{\p x} \]

现在观察\(-|\bfg|h\frac{\p( h+h_0)}{\p x}\)，有：

(5)\[\begin{split} -|\bfg|h\frac{\p( h+h_0)}{\p x} =-\frac{1}{2}|\bfg|\frac{\p h^2}{\p h}\frac{\p( h+h_0)}{\p x} =-\frac{1}{2}|\bfg|\frac{\p h^2}{\p h}\frac{\p h}{\p x} - \frac{1}{2}|\bfg|\frac{\p h^2}{\p h}\frac{\p h_0}{\p x} \\\\ =-\frac{1}{2}|\bfg|\frac{\p h^2}{\p x} - |\bfg| h\frac{\p h_0}{\p x} \end{split}\]

将方程(5)代入到(4)有：

(6)\[ \frac{\partial h u}{\partial t}+\frac{\p \left(huu + \frac{1}{2}|\bfg| h^2\right) }{\p x} = -|\bfg| h\frac{\p h_0}{\p x} \]

大量的文献将其写为矢量形式：

(7)\[ \frac{\p \bfU}{\p t}+\frac{\p \mathbf{F}}{\p x}=\mathbf{S} \]

其中

(8)\[\begin{split} \bfU= \begin{bmatrix} h \\ hu \end{bmatrix}, \mathbf{F} = \begin{bmatrix} hu \\ hu^2+\frac{1}{2}|\bfg| h^2 \end{bmatrix}, \mathbf{S} = \begin{bmatrix} 0 \\ -|\bfg| h\frac{\p h_0}{\p x} \end{bmatrix} \end{split}\]

方程(7)可以进一步的用来分析方程的特征值以及特征向量，其被文献证明为一个双曲系统。

离散

回到OpenFOAM中求解的方程(1)与(2)。在求解中动量方程可以使用守恒变量或者原始变量。OpenFOAM求解的为守恒变量，即将\(h\bfU\)看做一个整体。在这里定义：

(9)\[ \mathcal{U}_h=h\bfU \]

因此有：

(10)\[ \frac{\p \mathcal{U}_h}{\partial t}+\nabla\cdot(\bfU \mathcal{U}_h) =-\bfF\times\bfU-|\bfg|h\nabla(h+h_0) \]

如果不考虑科氏力，则可以简化为

(11)\[ \frac{\p \mathcal{U}_h}{\partial t}+\nabla\cdot(\bfU \mathcal{U}_h) = -|\bfg|h\nabla(h+h_0) \]

方程(1)与(11)是shallowWaterFoam求解的浅水方程。(1)可以求解\(h\)，(11)可以求解\(\mathcal{U}_h\)。但由于二者直接互相耦合，如果在一个时间步内顺序求解，存在失耦特性，因此迭代算法必不可少。

现将整个迭代思路进行简述。参考CFD:不可压+瞬态算法，若进行动量预测步骤，则需要针对方程(11)对\(\mathcal{U}_h\)做离散，其中时间项以及对流项均为隐性离散，\(h\)的梯度项需要显性离散，求解后有\(\mathcal{U}_h^*\)。若不进行动量预测步骤，则仅仅使用时间项与对流项进行离散并组建矩阵，\(\mathcal{U}_h\)不更新。两种情况下都可组建\(\bfHbyA_\rP^{*}\)：

(12)\[ \mathbf{HbyA}^{*}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}{\mathcal{U}_h}_\mathrm{N}^{t}} + S_\mathrm{P}^n} \right) \]

注意：

若进行动量预测，方程右侧为\(\mathcal{U}_h^*\)。若不进行动量预测，方程右侧为\(\mathcal{U}_h^t\)。

其中\(\bfHbyA_\rP^*\)并没有考虑\(-|\bfg|h\nabla(h+h_0)\)的贡献。因此在后续组建\(h\)方程的时候，需要考虑体积力（浮力）的附加效应。将\(\bfHbyA_\rP^*\)可以带入到方程(1)中可求解\(h\)，求解\(h\)后，可更新\(\mathcal{U}_h^{**}\)。由于\(h\)与\(\mathcal{U}_h^{**}\)的解耦求解，因此需要迭代2-3次。此即浅水方程求解的思想。

对于连续性方程，将其进行离散且考虑\(h\)的梯度项贡献后有：

(13)\[ V_\rP\int\frac{\partial h }{\partial t}\rd t+ \left(\sum \left(\mathbf{HbyA}_f^{*} -\frac{|\bfg|h_f}{A_{\rP,f}}\nabla(h^t+h_0) \right) \cdot \bfS_f \right) \dt =0 \]

(14)\[ \frac{h^{*}-h^t}{\dt}+ \frac{1}{V_\rP}\sum \left(\mathbf{HbyA}_f^{*} -\frac{\bfg|h_f}{A_{\rP,f}}\nabla_f(h^t+h_0) \right) \cdot \bfS_f =0 \]

(15)\[ \frac{h^{*}-h^t}{\dt}+ \frac{1}{V_\rP}\sum \mathbf{HbyA}_f^{*} \cdot \bfS_f - \frac{1}{V_\rP}\sum \left(\frac{\bfg|h_f}{A_{\rP,f}}\nabla_f(h^t+h_0) \right) \cdot \bfS_f =0 \]

求解方程(15)后，有\(h^*\)。其可用于更新\(\mathcal{U}_h\)：

(16)\[ \mathcal{U}_h^{**}=\mathbf{HbyA}^{*}-\frac{|\bfg|h_f}{A_{\rP,f}}\nabla(h^*+h_0) \]

有了\(\mathcal{U}_h^{**}\)，可以通过方程(12)继续更新下一个迭代步的\(\bfHbyA\)，进而继续更新\(\mathcal{U}_h\)。整个迭代过程迭代2-3次即可进行下一个时间步的推进。

关键代码

fvVectorMatrix hUEqn
(
    fvm::ddt(hU)
  + fvm::div(phiv, hU)
);

上述代码对应无源项的方程(11)。

fvScalarMatrix hEqn
(
    fvm::ddt(h)
  + fvc::div(phiHbyA)
  - fvm::laplacian(ghrAUf, h)
);

上述代码对应无源项的方程(15)。

hU = HbyA - rAU*h*magg*fvc::grad(h + h0);

上述代码对应无源项的方程(16)。