流体体积法（VOF）求解两相流问题的步骤

北京理工大学 | 李明健

控制方程：

对于两相流动，采用 one-fluid 模型，即用统一的物理场描述两种流体，用体积分数的方式区分。例如，密度场可表示为如下形式：

\begin{matrix} (1) & ρ = α ρ_{1} + (1 - α) ρ_{2} \end{matrix}

$\alpha$ $\rho_1$ $\rho_2$ 分别为两种流体的密度。对于粘度等其他物理量，用同样的方式处理。

连续性方程

不可压流体问题的连续性方程通常写为速度无散形式如下：

\begin{matrix} (2) & \nabla \cdot u = 0 \end{matrix}

动量方程

Navier-Stokes方程采用守恒形式如下：

\begin{matrix} (3) & ρ \frac{\partial u}{\partial t} + ρ \nabla \cdot (u u) = - \nabla p + \nabla \cdot μ [\nabla u + (\nabla u)^{T}] + f \end{matrix}

$\boldsymbol u$ $p$ $\rho$ $\mu$ $\boldsymbol f$ 分别为速度、密度、压力、粘度、体力。

求解方法：

采用 Chorin 投影法，我们把N-S方程分裂成以下形式：

\begin{matrix} (4) & \begin{matrix} \frac{u^{*} - u^{n}}{Δ t} = - C^{n} + \frac{1}{ρ^{n}} (D^{n} + f^{n}) \\ \frac{u^{n + 1} - u^{*}}{Δ t} = - \frac{1}{ρ^{n}} \nabla p \end{matrix} \end{matrix}

$\boldsymbol C$ $\boldsymbol D$ 分别为对流项和扩散项，两式相加即为N-S方程。

对第二式取散度得到压力泊松方程：

\begin{matrix} (5) & \nabla \cdot (\frac{1}{ρ^{n}} \nabla p) = \frac{1}{Δ t} \nabla \cdot u^{*} \end{matrix}

采用交错网格，在格心求解压力，面心求解速度。

对流项

对流项具体展开形式为：

\begin{matrix} (6) & \begin{matrix} C = \nabla \cdot (u u) = [\begin{array}{lll} \frac{\partial (u u)}{\partial x} + \frac{\partial (u v)}{\partial y} + \frac{\partial (u w)}{\partial z} \\ \frac{\partial (u v)}{\partial x} + \frac{\partial (v v)}{\partial y} + \frac{\partial (v w)}{\partial z} \\ \frac{\partial (u w)}{\partial x} + \frac{\partial (v w)}{\partial y} + \frac{\partial (w w)}{\partial z} \end{array}] \end{matrix} \end{matrix}

在有限体积法中，根据奥-高公式，将体积分变成面积分，对流项表示如下：

\begin{matrix} (7) & C = \frac{1}{V} \int_{V} \nabla \cdot (u u) d V = \frac{1}{V} \oint_{A} u (u \cdot n) d A \end{matrix}

$x$ 方向面心对流项分量为例：

\begin{matrix} (8) & \begin{matrix} (C_{x})_{i + 1 / 2, j, k} = \frac{1}{V} {[(u u)_{i + 1, j, k} - (u u)_{i, j, k}] A_{x} \\ + [(u v)_{i + 1 / 2, j + 1 / 2, k} - (u v)_{i + 1 / 2, j - 1 / 2, k}] A_{y} \\ + [(u w)_{i + 1 / 2, j, k + 1 / 2} - (u w)_{i + 1 / 2, j, k - 1 / 2}] A_{z}} \end{matrix} \end{matrix}

vof1

$A_x, A_y, A_z$ $x$ $u$ 为例，插值形式如下：

\begin{matrix} (9) & \begin{matrix} u_{i, j, k} = {\begin{cases} \frac{3}{8} u_{i + 1 / 2, j, k} + \frac{3}{4} u_{i - 1 / 2, j, k} - \frac{1}{8} u_{i - 3 / 2, j, k}, \frac{1}{2} (u_{i + 1 / 2, j, k} + u_{i - 1 / 2, j, k}) > 0 \\ \frac{3}{8} u_{i - 1 / 2, j, k} + \frac{3}{4} u_{i + 1 / 2, j, k} - \frac{1}{8} u_{i + 3 / 2, j, k}, \frac{1}{2} (u_{i + 1 / 2, j, k} + u_{i - 1 / 2, j, k}) < 0 \end{cases} \end{matrix} \end{matrix}

vof1

$x$ $u$ 为例，插值形式如下：

\begin{matrix} (10) & \begin{matrix} u_{i + 1 / 2, j + 1 / 2, k} = {\begin{cases} \frac{3}{8} u_{i + 1 / 2, j + 1, k} + \frac{3}{4} u_{i + 1 / 2, j, k} - \frac{1}{8} u_{i + 1 / 2, j - 1, k}, \frac{1}{2} (v_{i, j + 1 / 2, k} + v_{i + 1, j + 1 / 2, k}) > 0 \\ \frac{3}{8} u_{i + 1 / 2, j - 1, k} + \frac{3}{4} u_{i + 1 / 2, j, k} - \frac{1}{8} u_{i + 1 / 2, j + 1, k}, \frac{1}{2} (v_{i, j + 1 / 2, k} + v_{i + 1, j + 1 / 2, k}) < 0 \end{cases} \end{matrix} \end{matrix}

vof3

类似可以得到所有速度分量，并得到对流项所有分量。

扩散项

粘性扩散项具体展开形式为：

\begin{matrix} (11) & \begin{matrix} D = \nabla \cdot τ = [\begin{array}{lll} \frac{\partial τ^{x x}}{\partial x} + \frac{\partial τ^{x y}}{\partial y} + \frac{\partial τ^{x z}}{\partial z} \\ \frac{\partial τ^{y x}}{\partial x} + \frac{\partial τ^{y y}}{\partial y} + \frac{\partial τ^{y z}}{\partial z} \\ \frac{\partial τ^{z x}}{\partial x} + \frac{\partial τ^{z y}}{\partial y} + \frac{\partial τ^{z z}}{\partial z} \end{array}] \end{matrix} \end{matrix}

其中粘性应力具体展开形式为：

\begin{matrix} (12) & \begin{matrix} τ = [\begin{matrix} 2 \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} & \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \\ \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} & 2 \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \\ \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} & \frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} & 2 \frac{\partial w}{\partial z} \end{matrix}] \end{matrix} \end{matrix}

扩散项写成面积分的形式为：

\begin{matrix} (13) & D = \frac{1}{V} \int_{V} \nabla \cdot τ d V = \frac{1}{V} \oint_{A} τ \cdot n d A, τ = μ [\nabla u + (\nabla u)^{T}] \end{matrix}

$x$ 方向面心分量为例：

\begin{matrix} (14) & (D_{x})_{i + 1 / 2, j, k} = \frac{τ_{i + 1, j, k}^{x x} - τ_{i, j, k}^{x x}}{Δ x} + \frac{τ_{i + 1 / 2, j + 1 / 2, k}^{x y} - τ_{i + 1 / 2, j - 1 / 2, k}^{x y}}{Δ y} + \frac{τ_{i + 1 / 2, j, k + 1 / 2}^{x z} - τ_{i + 1 / 2, j, k - 1 / 2}^{x z}}{Δ z} \end{matrix}

vof1

粘性应力分量采用二阶中心差分得到，以其中两个分量为例：

\begin{matrix} (15) & \begin{matrix} τ_{i, j, k}^{x x} = 2 μ_{i, j, k} \frac{u_{i + 1 / 2, j, k} - u_{i - 1 / 2, j, k}}{Δ x} \\ τ_{i + 1 / 2, j + 1 / 2, k}^{x y} = μ_{i + 1 / 2, j + 1 / 2, k} (\frac{u_{i + 1 / 2, j + 1, k} - u_{i + 1 / 2, j, k}}{Δ x} + \frac{v_{i + 1, j + 1 / 2, k} - v_{i, j + 1 / 2, k}}{Δ y}) \end{matrix} \end{matrix}

vof1

其余分量同理，最终得到扩散项所有分量。

速度预测

Chorin 投影法第一步是预测一个中间速度：

\begin{matrix} (16) & u^{*} = u^{n} - Δ t C^{n} + Δ t \frac{1}{ρ^{n}} (D^{n} + f^{n}) \end{matrix}

预测速度在面上求解：

\begin{matrix} (17) & u_{i + 1 / 2, j, k}^{*} = u_{i + 1 / 2, j, k}^{n} - Δ t {C_{i + 1 / 2, j, k}^{n} + \frac{1}{ρ_{i + 1 / 2, j, k}^{n}} (D_{i + 1 / 2, j, k}^{n} + f_{i + 1 / 2, j, k}^{n})} \end{matrix}

压力泊松方程

Chorin 投影法第二步是求解压力泊松方程：

\begin{matrix} (18) & \nabla \cdot (\frac{1}{ρ^{n}} \nabla p) = \frac{1}{Δ t} \nabla \cdot u^{*} \end{matrix}

在格心求解格式为：

\begin{matrix} (19) & \begin{matrix} \frac{1}{Δ x^{2}} (\frac{p_{i + 1, j, k} - p_{i, j, k}}{ρ_{i + 1 / 2, j, k}} - \frac{p_{i, j, k} - p_{i - 1, j, k}}{ρ_{i - 1 / 2, j, k}}) \\ + \frac{1}{Δ y^{2}} (\frac{p_{i, j + 1, k} - p_{i, j, k}}{ρ_{i, j + 1 / 2, k}} - \frac{p_{i, j, k} - p_{i, j - 1, k}}{ρ_{i, j - 1 / 2, k}}) \\ + \frac{1}{Δ z^{2}} (\frac{p_{i, j, k + 1} - p_{i, j, k}}{ρ_{i, j, k + 1 / 2}} - \frac{p_{i, j, k} - p_{i, j, k - 1}}{ρ_{i, j, k - 1 / 2}}) \\ = \frac{1}{Δ t} (\frac{u_{i + 1 / 2, j, k}^{*} - u_{i - 1 / 2, j, k}^{*}}{Δ x} + \frac{v_{i, j + 1 / 2, k}^{*} - v_{i, j - 1 / 2, k}^{*}}{Δ y} + \frac{w_{i, j, k + 1 / 2}^{*} - w_{i, j, k - 1 / 2}^{*}}{Δ z}) \end{matrix} \end{matrix}

vof1

压力泊松方程可以采用SOR（successive over-relaxation）迭代法求解，或共轭梯度等其他方法。

速度修正

在面心用压力梯度修正中间速度，得到最终的速度场：

\begin{matrix} (20) & u^{n + 1} = u^{*} - \frac{Δ t}{ρ^{n}} \nabla p \end{matrix}

在面心求解格式为：

\begin{matrix} (21) & u_{i + 1 / 2, j, k}^{n + 1} = u_{i + 1 / 2, j, k}^{*} - \frac{Δ t}{ρ_{i + 1 / 2, j, k}} \frac{p_{i + 1, j, k} - p_{i, j, k}}{Δ x} \end{matrix}

vof1

密度和粘度

为了处理大密度比和粘性比，首先在面心通过PLIC中的切割等方法得到体积分数，随后，在面上用体积分数更新密度和粘度：

\begin{matrix} (22) & \begin{matrix} ρ_{i + 1 / 2, j, k} = α_{i + 1 / 2, j, k} ρ_{1} + (1 - α_{i + 1 / 2, j, k}) ρ_{2} \\ μ_{i + 1 / 2, j, k} = α_{i + 1 / 2, j, k} μ_{1} + (1 - α_{i + 1 / 2, j, k}) μ_{2} \end{matrix} \end{matrix}