• 返回主页
  • interPhaseChangeFoam解析

    于海东(清华大学航空航天学院)
    李东岳(亥姆霍茨计算流体力学研究所)
    求解器类型interPhaseChangeFoam是一个基于interFoam,附加相变(如空化)的不可压缩两相流求解器。其中界面处理方式是利用VOF方法进行捕获。在interFoam求解器的基础之上,interPhaseChangeFoam附加了由于相变引起质量传递的子模型。

    在推导interPhaseChangeFoam的相方程中,由于附加质量传递,首先考虑每一相各自的连续性方程(质量守恒方程): \begin{equation} \frac{\partial \alpha_\rl\rho_\rl}{\partial t}+\nabla \cdot (\alpha_\rl\rho_\rl \mathbf{U}_\rl)=\dot{m} \label{alphal} \end{equation} \begin{equation} \frac{\partial \alpha_\rv\rho_\rv}{\partial t}+\nabla \cdot (\alpha_\rv\rho_\rv \mathbf{U}_\rv)=-\dot{m} \label{alphav} \end{equation} 其中下标$\rl$表示液相,$\rv$表示气相,$\alpha_\rl$表示液相的体积分数,$\alpha_\rv$表示气相的体积分数,$\rho_\rl$表示液相密度,$\rho_\rv$表示气相密度,$\bf{U}_\rl$表示液相的速度矢量,$\bf{U}_\rv$表示气相的速度矢量,$\dot{m}$表示每单位体积相的质量交换。考虑不可压缩流体,提出方程中的密度,方程($\ref{alphal}$)及($\ref{alphav}$)可以化为: \begin{equation} \frac{\partial \alpha_\rl}{\partial t}+\nabla \cdot (\alpha_\rl\mathbf{U}_\rl)=\frac{\dot{m}}{\rho_\rl} \label{alphal_incomp} \end{equation} \begin{equation} \frac{\partial \alpha_\rv}{\partial t}+\nabla \cdot (\alpha_\rv\mathbf{U}_\rv)=-\frac{\dot{m}}{\rho_\rv} \label{alphav_incomp} \end{equation} 在这里需要注意的是,VOF方法通常求解一个相方程(如液相的方程),气相则由$\alpha_\rv=1-\alpha_\rl$得出。因此方程($\ref{alphal_incomp}$)即为interPhaseChangeFoam需要求解的相方程,在后续我们将具体展开讨论。

    将方程($\ref{alphal_incomp}$)和($\ref{alphav_incomp}$)加和有: \begin{equation} \nabla\cdot\bf{U}=\frac{\dot{m}}{\rho_l}-\frac{\dot{m}}{\rho_v} \label{U} \end{equation} 其中$\bf{U}=\alpha_l\bfU_l+\alpha_v\bfU_v$。方程($\ref{U}$)即为连续性方程,其主要用于组件压力泊松方程后对速度场附加连续性限制。

    interPhaseChangeFoam求解的动量方程和interFoam相同即(此处推导略): \begin{equation} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau = C \kappa \nabla \alpha -g h \nabla \rho -\nabla p_\mathrm{rgh} \label{mom} \end{equation}

    相方程

    在VOF模型中,相分数作为一个特殊的被动标量,具有明显的数值特性。首先是相分数$\alpha$在数值上应该保证在$0-1$之间,任何越界的相分数都不符合物理。其次,如果考虑网格非常细密的情况下,相分数应该只存在两个数值:0或者1。通常情况下这种非常致密的网格难以达到,因此会存在$\alpha$在$0-1$之间的情况。最后,在界面重构框架下的VOF方法应该使得相界面保证足够的尖锐。在界面跟踪框架下,相界面可以尖锐的通过动网格或者示踪颗粒来获得。比如在通过网格的界面跟踪框架下,界面外的相分数只存在两个数值。目前免费CFD软件OpenFOAM中植入的是界面重构框架下的VOF模型。由于界面重构框架并不会重构尖锐的相界面,因此在界面重构框架下又分为不同的数值方法。这些数值方法在保证有界的同时,还可以使用高分辨率格式(如HRIC,CICSAM)。但是并没有采用几何重构(有关界面重构请参考:链接)来进行相界面的尖锐,因此在求解相方程的过程中,需要一种数值的方法来保证相界面的尖锐。为了降低相界面附近的相分数的数值耗散,Weller提出通过一种人工的对流项来对相界面附近的相分数进行挤压来抗衡这种数值耗散带来的相界面模糊性,并且这一项在数值上需要保证非相界面处为零。

    如果使用$\mathbf{U}_\mathrm{r}$表示相对速度且$\bfU_\rr=\bfU_\rl-\bfU_\rv$,Weller在双流体模型中将方程\eqref{alphal_incomp}化为(此处不考虑源项): \begin{equation} \frac{\partial \alpha_\rl}{\partial t}+\nabla \cdot (\alpha_\rl\mathbf{U})+\nabla\cdot(\alpha_\rl\alpha_\rv\mathbf{U}_\mathrm{r})=0 \label{wellerTFM} \end{equation} 由于在VOF模型中并不存在多个相速度,因此Weller将方程\eqref{wellerTFM}中的相对速度进行模化: \begin{equation} \bfU_\rr=c_\alpha|\bfU| \end{equation} 其中$c_\alpha$为相界面压缩因子,其可以进一步的进行模化来保证数值稳定。这样,无源项的方程\eqref{alphal_incomp}可以表示为: \begin{equation} \frac{\partial \alpha_l}{\partial t}+\nabla \cdot (\alpha_l\mathbf{U})+\nabla\cdot(\alpha_l\alpha_v\mathbf{U}_\rr)=0 \label{wellerVOF} \end{equation} 方程\eqref{wellerVOF}中模化后的$\mathbf{U}_\rr$用来保证界面的压缩性。方程\eqref{wellerVOF}为OpenFOAM中VOF方法下(如interFoam)求解的相方程。 在interPhaseChangeFoam中,附加源项的方程\eqref{wellerVOF}需要进一步重构,观察其中的$\nabla \cdot (\alpha_l\mathbf{U})$,其可以重构为: \begin{equation} \nabla \cdot (\alpha_l\mathbf{U})=\alpha_\rl\nabla\cdot\bfU+\bfU\nabla\alpha_\rl=\alpha_\rl\nabla\cdot\bfU \end{equation} 将其带入到附加源项的方程\eqref{wellerVOF}有: \begin{equation} \frac{\partial \alpha_l}{\partial t}+\alpha_\rl\nabla\cdot\bfU+\bfU\nabla\alpha_\rl+\nabla\cdot(\alpha_l\alpha_v\mathbf{U}_\rr)=\frac{\dot{m}}{\rho_l} \label{fummy} \end{equation} 将方程\eqref{U}代入到方程\eqref{fummy}有: \begin{equation} \frac{\partial \alpha_\rl}{\partial t}+\bfU\nabla\alpha_\rl+\nabla\cdot(\alpha_\rl\alpha_\rv\mathbf{U}_\rr)=-\alpha_\rl\left(\frac{\dot{m}}{\rho_l}-\frac{\dot{m}}{\rho_v}\right)+\frac{\dot{m}}{\rho_l} \end{equation} \begin{equation} \frac{\partial \alpha_\rl}{\partial t}+\nabla\cdot(\alpha_\rl\bfU)+\nabla\cdot(\alpha_l\alpha_v\mathbf{U}_\rr)=\alpha_\rl\nabla\cdot\bfU-\alpha_\rl\left(\frac{\dot{m}}{\rho_l}-\frac{\dot{m}}{\rho_v}\right)+\frac{\dot{m}}{\rho_l} \label{alphaLast} \end{equation} 方程\eqref{alphaLast}即为最终需要求解的相方程。

    动量方程

    方程\eqref{mom}即为求解的动量方程。在代码中,其左边的部分用来组建UEqn,右边的源项部分移至压力方程来处理,最终代码中求解的动量方程为: \begin{equation} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau = 0 \label{momPre} \end{equation} 动量方程的求解也即动量预测步骤。

    压力方程

    在组建方程\eqref{momPre}后,结合方程\eqref{U}可以用来组建压力方程。由于方程\eqref{momPre}并没有考虑方程\eqref{mom}中的右边部分,因此在组建压力方程的时候,需要更新通量。我们令$\bfHbyA^n$为求解方程\eqref{mom}后的预测速度,在考虑动量方程源项之后的预测速度为: \begin{equation} \bfU^{n+1} = \bfHbyA^n + C \kappa \nabla \alpha -g h \nabla \rho \end{equation} 其中的$\bfU^{n+1}$,在考虑压力梯度的贡献之后,应该符合方程\eqref{mom}。

    空化模型

    不同于interFoaminterPhaseChangeFoam需要附加质量传递模型,本文选取Sauer空化模型进行分析。 \begin{equation} \dot{m}=\frac{3\rho_l\rho_v}{\rho}\alpha(1-\alpha)\frac{1}{R_b}\sqrt{\frac{2|p_v-p|}{3\rho_l}} \end{equation} 其中$p_v$对应温度下的饱和蒸汽压力,$R_b$为平均泡半径: \begin{equation} R_b=\left(\frac{1-\alpha}{\alpha}\frac{3}{4\pi n}\right)^{1/3} \end{equation} $n$为液相中的空泡数目。

    对数值求解过程有个整体的认识,介绍流程如下:首先,对相方程进行求解,分为显式和隐式,更新$\alpha$及其通量和$\rho_l$;其次,由动量方程给出速度预测值;最后,求解泊松方程进行压力值更新与速度修正,并更新通量。较之interFoam求解器而言,该求解器最大的不同在于空化模型的加入和处理。

    空化模型需要通过四个参数,在计算设置时,需要用户在transportProperties内指定:

    参数 n dNuc Cc Cv
    物理意义 来流成核数 成核直径 冷凝率系数 蒸发率系数
    一般取值 1.6e13 2e-6 1.0 1.0


    上述参数均定义在SchnerrSauer.H中,在此基础上,延伸出重要的二级中间参数有三个,它们定义在SchnerrSauer.C的成员函数里:

    参数 alphaNuc() rRb pCoeff
    物理意义 蒸汽泡体积分数$\gamma$ 气泡半径倒数 公共系数
    一般取值 $\gamma=\frac{n\pi d^3/6}{1+n\pi d^3/6}$ $\sqrt[3]{\frac{4n\pi}{3}\alpha\frac{1}{1+\gamma-\alpha}}$ $\frac{3\rho_l\rho_v}{\rho}\mathrm{rRb}\sqrt{\frac{2}{3\rho_l}}\frac{1}{\sqrt{|p_v-p|+0.01p_v}}$


    最后给出质量传递源项,该项分别在求解相方程和泊松方程中加入。具体分别是:SchnerrSauer.C中的mDotAlphal()和mDotP()的质量传递源项;phaseChangeTwoPhaseMixture.C中的vDotAlphal()和vDotP()体积传递源项。这四项返回值均为Pair>,即为一个数组,代表了净空化传质率和净冷凝传质率,它们具体的形式在相方程和泊松方程的求解中会继续分析,下面给出具体形式:
    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国