高阶有限差分法求解可压缩流体的典型步骤

北京理工大学 | 李明健

有限差分法（Finite difference method，FDM）是计算流体力学中常用的数值方法，下面介绍其求解可压缩流体问题的典型过程（感谢我的学生于荣岳对本文公式整理做出的贡献）。

考虑无粘、可压缩的理想流体，给出以下欧拉方程：

\begin{matrix} (1) & \frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} + \frac{\partial H}{\partial z} = 0 \end{matrix}

其中：

\begin{matrix} (2) & \begin{matrix} U = [\begin{matrix} ρ \\ ρ u \\ ρ v \\ ρ w \\ E \end{matrix}], F = [\begin{matrix} ρ u \\ ρ u^{2} + p \\ ρ u v \\ ρ u w \\ u (E + p) \end{matrix}], G = [\begin{matrix} ρ v \\ ρ u v \\ ρ v^{2} + p \\ ρ v w \\ v (E + p) \end{matrix}], H = [\begin{matrix} ρ w \\ ρ w u \\ ρ w v \\ ρ w^{2} + p \\ w (E + p) \end{matrix}] \end{matrix} \end{matrix}

$\rho$ $u$ $v$ $\ w$ $x$ $y$ $z$ $p$ $E$ 为单位质量的总能量，写为内能和动能之和：

\begin{matrix} (3) & E = ρ e + \frac{1}{2} ρ (u^{2} + v^{2} + w^{2}) \end{matrix}

以上方程需要状态方程进行封闭，对于气、水混合介质，可采用以下状态方程：

\begin{matrix} (4) & p = (γ - 1) ρ e - γ B \end{matrix}

$B=0$ $\gamma = 1.4$ $\gamma = 7.5$ $B=3.309\times 10^8$ Pa 时退化为水的 Tait 状态方程。

流体声速满足以下形式：

\begin{matrix} (5) & c = \sqrt{γ (p + B) / ρ} \end{matrix}

$j$ 节点上的半离散格式如下：

\begin{matrix} (6) & \frac{\partial U_{j}}{\partial t} + \frac{F_{j + \frac{1}{2}} - F_{j - \frac{1}{2}}}{Δ x} + \frac{G_{j + \frac{1}{2}} - G_{j - \frac{1}{2}}}{Δ y} + \frac{H_{j + \frac{1}{2}} - H_{j - \frac{1}{2}}}{Δ z} = 0 \end{matrix}

采用通量分裂方法，进一步得到以下形式：

\begin{matrix} (7) & \begin{matrix} \frac{\partial U_{j}}{\partial t} + \frac{({\hat{F}}_{j + \frac{1}{2}}^{+} + {\hat{F}}_{j + \frac{1}{2}}^{-}) - ({\hat{F}}_{j - \frac{1}{2}}^{+} + {\hat{F}}_{j - \frac{1}{2}}^{-})}{Δ x} + \frac{({\hat{G}}_{j + \frac{1}{2}}^{+} + {\hat{G}}_{j + \frac{1}{2}}^{-}) - ({\hat{G}}_{j - \frac{1}{2}}^{+} + {\hat{G}}_{j - \frac{1}{2}}^{-})}{Δ y} \\ + \frac{({\hat{h}}_{j + \frac{1}{2}}^{+} + {\hat{H}}_{j + \frac{1}{2}}^{-}) - ({\hat{H}}_{j - \frac{1}{2}}^{+} + {\hat{H}}_{j - \frac{1}{2}}^{-})}{Δ z} = 0 \end{matrix} \end{matrix}

${\hat{\boldsymbol{F}}} ^{+}$ ${\hat{\boldsymbol{F}}} ^{-}$ $\boldsymbol F$ 分裂得到的正、负通量，其他同理。

Lax-Friedrichs 通量分裂算法 $x$ $\boldsymbol{F}$ 分裂为：

\begin{matrix} (8) & \begin{matrix} {\hat{F}}^{+} = (F + λ^{max} U) / 2 \\ {\hat{F}}^{-} = (F - λ^{max} U) / 2 \end{matrix} \end{matrix}

$\lambda^{\text{max}}$ ${\partial \boldsymbol F}/{\partial \boldsymbol U}$ 的最大特征值，通过：

\begin{matrix} (9) & λ^{max} = max (| u |, | u - c |, | u + c |) \end{matrix}

$y$ $z$ $v$ $w$ $\boldsymbol{G}$ $\boldsymbol{H}$ 进行分裂。

5 阶 WENO（Weighted Essentially Non-Oscillatory）格式 $x$ ${\hat{\boldsymbol{F}}}_{j + \frac{1}{2}}^{+}$ $j-2$ $j-1$ $j$ $j+1$ $j+2$ 五个模板点的值：

\begin{matrix} (10) & {\hat{F}}_{j + \frac{1}{2}}^{+} = \sum_{k = 0}^{2} ω_{k}^{+} q_{k}^{+} (F_{j - 2}^{+}, F_{j - 1}^{+}, F_{j}^{+}, F_{j + 1}^{+}, F_{j + 2}^{+}) \end{matrix}

${\hat{\boldsymbol{F}}}_{j + \frac{1}{2}}^{-}$ $j-1$ $j$ $j+1$ $j+2$ $j+3$ 五个模板点的值：

\begin{matrix} (11) & {\hat{F}}_{j + \frac{1}{2}}^{-} = \sum_{k = 0}^{2} ω_{k}^{-} q_{k}^{-} (F_{j - 1}^{-}, F_{j}^{-}, F_{j + 1}^{-}, F_{j + 2}^{-}, F_{j + 3}^{-}) \end{matrix}

${{\boldsymbol{F}}}_{j + \frac{1}{2}}$ 作为总通量。同理，得到所有通量后，即得到式（6）完整的半离散格式。式（10）和（11）中的权系数具体取值参考舒其望的论文。

$\boldsymbol U$ $\boldsymbol{F}$ $\boldsymbol{G}$ $\boldsymbol{H}$ $\boldsymbol U$ 表示，因此可将半离散的方程（7）简写为：

\begin{matrix} (12) & \frac{\partial U_{j}}{\partial t} = L (U) \end{matrix}

$\mathcal L$ 为空间离散算子。

采用三阶 Runge-Kutta 法对时间进行离散，即可得到全离散格式：

\begin{matrix} (13) & \begin{matrix} U_{j}^{(1)} = U_{j}^{n} + Δ t L (U^{n}) \\ U_{j}^{(2)} = \frac{3}{4} U_{j}^{n} + \frac{1}{4} U_{j}^{(1)} + \frac{1}{4} Δ t L (U^{(1)}) \\ U_{j}^{(3)} = \frac{1}{3} U_{j}^{n} + \frac{2}{3} U_{j}^{(2)} + \frac{2}{3} Δ t L (U^{(2)}) \end{matrix} \end{matrix}

每个时间步分为三个子步，每个子步内的求解流程如下：

$\rho$ $p$ $c$ ；
$\rho$ $p$ $e$ $\rho$ $p$ $u$ $v$ $w$ $e$ $E$ ；
$\rho$ $p$ $u$ $v$ $w$ $E$ $\boldsymbol U$ $\boldsymbol{F}$ $\boldsymbol{G}$ $\boldsymbol{H}$ ；
$u$ $v$ $w$ $c$ $\lambda^\text{max}$ $\lambda^\text{max}$ $\boldsymbol U$ $\boldsymbol{F}$ $\boldsymbol{G}$ $\boldsymbol{H}$ $\hat{\boldsymbol{F}}^{+}, \hat{\boldsymbol{F}}^{-}, \hat{\boldsymbol{G}}^{+}, \hat{\boldsymbol{G}}^{-}, \hat{\boldsymbol{H}}^{+}, \hat{\boldsymbol{H}}^{-}$ ；
$\hat{\boldsymbol{F}}^{+}, \hat{\boldsymbol{F}}^{-}, \hat{\boldsymbol{G}}^{+}, \hat{\boldsymbol{G}}^{-}, \hat{\boldsymbol{H}}^{+}, \hat{\boldsymbol{H}}^{-}$ $\mathcal{L}$ ；
$\boldsymbol U$ $\mathcal{L}$ $\boldsymbol U$ ；
$\boldsymbol U$ $\rho$ $u$ $v$ $w$ $e$ $E$ $p$ 。

重复以上步骤，即可完成求解。