shallowWaterFoam解析

不考虑粘性的浅水方程是标准的双曲偏微分方程系统。其有的时候也会被称之为Saint-Venant方程。浅水方程的假设是高度方向的流动尺度要远远小于水平面的尺度。在这种情况下，可以认为高度方向的速度近似为0。

控制方程与算法

(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) \]

其中\(\bfF=(2\boldsymbol{\Omega}\cdot\bfn_g)\bfn_g\)，\(\bfn_g\)表示重力方向的单位矢量。方程(1)表示的是关于水深\(h\)的连续性方程。方程(2)表示的是关于\(\bfU\)的动量方程。在求解中动量方程可以使用守恒变量或者原始变量。OpenFOAM求解的为守恒变量，即将\(h\bfU\)看做一个整体。在这里定义：

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

因此有：

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

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

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

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

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

(6)\[ \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\)的梯度项贡献后有：

(7)\[ 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 \]

(8)\[ \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 \]

(9)\[ \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 \]

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

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

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

关键代码

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

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

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

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

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

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