compressibleInterFoam解析


如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次

1. 引言

compressibleInterFoam是一个可压缩版本的VOF模型求解器,可用于对气液界面现象进行模拟。可以将它理解为interFoam求解器的可压缩版。compressibleInterFoam可用于模拟可压缩界面形变问题(如气泡溃灭),以及涉及传热的界面问题。在将不可压缩求解器拓展到可压缩求解器的情况下,主要需要附加考虑密度的变化和能量方程。有关不可压缩VOF模型的详细推导请参考interFoam解析。本文在interFoam解析的基础上,讨论界面类问题可压缩性的植入。首先,不可压缩流体VOF模型的控制方程为

\begin{equation}\label{vof1} \frac{\partial \alpha}{\partial t}+\bfU\cdot\nabla\alpha=0 \end{equation}
\begin{equation}\label{vof2} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =-\nabla p_\mathrm{rgh}-\bfg\cdot \bfh \nabla \rho+\sigma \kappa \nabla \alpha \end{equation}
\begin{equation}\label{vof3} \nabla\cdot\bfU=0 \end{equation}
下一节将讨论可压缩VOF模型的控制方程。

2. 控制方程与算法
2.1. 可压缩相方程

可压缩VOF模型中的动量方程与不可压缩VOF模型中的动量方程相同。连续性方程\eqref{vof3}则需要考虑密度的变化进行适配:

\begin{equation}\label{drhodt} \frac{\p \rho}{\p t}+\nabla\cdot\left(\rho\bfU\right)=0 \end{equation}
同时,可压缩VOF模型的相方程为:
\begin{equation}\label{alpha1} \frac{\p \alpha_1 \rho_1}{\partial t}+\nabla \cdot (\rho_1 \alpha_1 \mathbf{U})=0 \end{equation}
将方程\eqref{drhodt}减去\eqref{alpha1}有
\begin{equation}\label{alpha2} \frac{\partial \alpha_2 \rho_2}{\partial t}+\nabla \cdot (\rho_2 \alpha_2 \mathbf{U})=0 \end{equation}
对方程\eqref{alpha1}展开有:
\begin{equation} \alpha_1 \frac{\partial \rho_1}{\partial t}+\rho_1 \frac{\partial \alpha_1}{\partial t} +\rho_1 \nabla \cdot (\alpha_1 \mathbf{U}) + \alpha_1 \mathbf{U} \nabla \rho_1=0 \end{equation}
即:
\begin{equation}\label{alpha1dt} \frac{\partial \alpha_1}{\partial t}+\nabla \cdot (\alpha_1 \mathbf{U})=-\frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \end{equation}
同理有:
\begin{equation}\label{alpha2dt} \frac{\partial \alpha_2}{\partial t}+\nabla \cdot (\alpha_2 \mathbf{U})=-\frac{\alpha_2}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} \end{equation}
将方程\eqref{alpha1dt}、\eqref{alpha2dt}相加有:
\begin{equation}\label{divU} \nabla \cdot \mathbf{U}=-\frac{\alpha_2}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t}-\frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \end{equation}
同时对方程\eqref{alpha1dt}进一步展开有:
\begin{equation}\label{divU2} \frac{\partial \alpha_1}{\partial t}+\mathbf{U} \cdot \nabla \alpha_1+\alpha_1\nabla \cdot \mathbf{U}=-\frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \end{equation}
将方程\eqref{divU}代入到方程\eqref{divU2}有:
\begin{equation}\label{alpha1d} \frac{\partial \alpha_1}{\partial t}+\mathbf{U} \cdot \nabla \alpha_1=\alpha_1 \alpha_2 \left( \frac{1}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} - \frac{1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \right) \end{equation}
在这里,用dgdt表示可压缩项:
\begin{equation} \mathrm{dgdt}= \left( \frac{1}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} - \frac{1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \right) \end{equation}
方程\eqref{alpha1d}可以化为:
\begin{equation}\label{alpha1dgdtpre} \frac{\partial \alpha_1}{\partial t}+\mathbf{U} \cdot \nabla \alpha_1=\alpha_1 \alpha_2 \mathrm{dgdt} \end{equation}
将方程\eqref{alpha1dgdtpre}左右加上$\alpha_1 \nabla \cdot \mathbf{U}$有守恒形式:
\begin{equation}\label{alpha1dgdt0} \frac{\partial \alpha_1}{\partial t}+ \nabla \cdot (\alpha_1 \mathbf{U})=\alpha_1 \alpha_2 \mathrm{dgdt}+\alpha_1 \nabla \cdot \mathbf{U} \end{equation}
参考interFoam解析中的界面压缩技术,其可以写为
\begin{equation}\label{alpha1dgdt} \frac{\partial \alpha_1}{\partial t}+ \nabla \cdot (\alpha_1 \mathbf{U})+\nabla\cdot\left(\alpha_1(1-\alpha_1)\bfU_c\right)=\alpha_1\alpha_2 \mathrm{dgdt}+\alpha_1 \nabla \cdot \mathbf{U} \end{equation}
方程\eqref{alpha1dgdt}为compressibleInterFoam中需要求解的相方程的最终形式。

方程\eqref{alpha1dgdt}在求解的过程中,为了保证相分数的有界性。需要对对流项以及非线性源项$\alpha_1\alpha_2\mathrm{dgdt}$进行特殊的数值处理。首先看源项部分,为了更好的理解数值内容,考虑已经对除了非线性源项进行欧拉显性离散后的方程:

\begin{equation}\label{test} \frac{\alpha_1^{t+\Delta t}-\alpha_1^t}{\Delta t}+\sum\phi_f^{t}\alpha_{1,f}^{t}=\alpha_1\alpha_2\mathrm{dgdt}+\alpha_1^{t} \nabla \cdot \mathbf{U}^{t} \end{equation}
其中$\alpha_1^{t} \nabla \cdot \mathbf{U}^{t}$应该进行显性离散,因为其来源于方程左边的对流项。同时,若$\mathrm{dgdt}<0$,$\alpha_1\alpha_2\mathrm{dgdt}$中可对$\alpha_1$隐性离散,$\alpha_2$显性离散为$1-\alpha_1^t$:
\begin{equation}\label{test3} \frac{\alpha_1^{t+\Delta t}-\alpha_1^t}{\Delta t}+\sum\phi_f^t\alpha_{1,f}^t=\alpha_1^{t+\Delta t}\mathrm{Sp}+\mathrm{Su} \end{equation}
其中
\begin{equation} \mathrm{Sp}=(1-\alpha_1^t)\mathrm{dgdt} \end{equation}
\begin{equation} \mathrm{Su}=\alpha_1^{t} \nabla \cdot \mathbf{U}^t \end{equation}
反之,若$\mathrm{dgdt}>0$,$\alpha_1\alpha_2\mathrm{dgdt}$中可对$\alpha_1$显性离散,$\alpha_2$隐性离散为$1-\alpha_1^{t+\Delta t}$:
\begin{equation}\label{test4} \frac{\alpha_1^{t+\Delta t}-\alpha_1^t}{\Delta t}+\sum\phi_f^t\alpha_{1,f}^t = -\alpha_1^{t+\Delta t}\mathrm{Sp}+\mathrm{Su} \end{equation}
其中
\begin{equation} \mathrm{Sp}=\alpha_1^t\mathrm{dgdt} \end{equation}
\begin{equation} \mathrm{Su}=\alpha_1^{t}\mathrm{dgdt}+ \alpha_1^{t}\nabla \cdot \mathbf{U}^t \end{equation}
方程\eqref{test4}和\eqref{test3}中等号右侧第一项一直小于0,在移项之后可保证$\alpha^{t+\Delta t}$表达式中右侧的分母($1/\Delta t -\mathrm{Sp}$)一直大于0,保证解符合物理、$\alpha_1$和$\alpha_2$的耦合以及有界性。

2.2. 压力方程

方程\eqref{divU}可以写为

\begin{equation}\label{pressure0} \nabla \cdot \mathbf{U}+\frac{\alpha_2}{\rho_2}\left(\frac{\p\rho_1}{\p t}+\bfU\cdot\nabla\rho_1\right)+\frac{\alpha_1}{\rho_1}\left(\frac{\p\rho_2}{\p t}+\bfU\cdot\nabla\rho_2\right)=0 \end{equation}
依据链条法则:
\begin{equation}\label{chain} \frac{\p\rho}{\p t}=\frac{\p\rho}{\p p}\frac{\p p}{\p t}=\psi\frac{\p p}{\p t} \end{equation}
代入有
\begin{equation}\label{pressure1} \nabla \cdot \mathbf{U}+\frac{\alpha_2\psi_2}{\rho_2}\left(\frac{\p p}{\p t}+\bfU\cdot\nabla p\right)+\frac{\alpha_1\psi_1}{\rho_1}\left(\frac{\p p}{\p t}+\bfU\cdot\nabla p\right)= \\\\ \nabla \cdot \mathbf{U}+\left(\frac{\alpha_2\psi_2}{\rho_2}+\frac{\alpha_1\psi_1}{\rho_1}\right)\left(\frac{\p p}{\p t}+\bfU\cdot\nabla p\right)=0 \end{equation}
其中$\nabla \cdot \mathbf{U}$的处理参考不可压缩VOF模型中的处理方法,可用来组建压力泊松方程。$\frac{\p p}{\p t}$可用来隐性离散。在近音速情况下,$\bfU\cdot\nabla p$隐性离散。

2.3. 能量方程

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网