return 0;
wmake

interFoam解析

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


VOF(volume of fraction)模型为一种多相流模型。其直接捕获相见的界面,可以看作为多相流的直接模拟。在OpenFOAM中,VOF模型被植入到了interFoam中,用于求解多相流的界面问题,如溃坝、极少量的气泡界面、液体晃荡等。本文从最基本的液膜压力降开始分析,一步一步推导VOF模型。

VOF模型中定义了一个体分数变量$\alpha$,其表示流体的相分数。考虑某一个网格单元的气液两相系统,如果此网格单元内充满了流体,则$\alpha=1$;如果此网格单元内充满了气体,则则$\alpha=0$。如果$\alpha$的值介于$0$和$1$之间,则此网格单元内为气液混合。可见,CFD后处理中VOF的界面并不是直接算出来的,而是通过计算每个网格内的相场值而后处理出来的。网格越密,VOF后处理出来的界面越薄。考虑下图的小界面微元,定义表面张力为$\sigma$,其是一种为了保持界面平衡而作用在每单位长度上的力,是一种和界面相切的张力。表面张力的大小主要和流体的物理特性有关。在曲面中,均衡的流体表面张力和界面压力降$\Delta p$均衡。这个界面压力降主要取决于表面张力和曲面弧度。


图1. 两相界面处的曲面微元受力示意图

如图(1)所示的一个曲面微元。其中$p_1$为曲面微元上方液体向下施加的压强,$p_2$为曲面微元下方液体向上施加的压强,$s_1$和$s_2$分别表示曲面微元的俩个边长。曲面微元上的总压力的值为:

\begin{equation}\label{Fp} F_p=(p_1 - p_2) \mathrm{d}s_1 \mathrm{d}s_2 \end{equation}

基于表面张力的定义,其作用在曲面微元的四个边上,且合力向上和曲面微元上的总压力$F_p$均衡。如果表面张力系数定义为$C$。那么作用在$\mathrm{d}s_1$上的表面张力为$C\mathrm{d}s_1$,其在$y$方向上的分量为:$C\mathrm{d}s_1 \mathrm{sin} \frac{\mathrm{d}\alpha}{2} \approx C\mathrm{d}s_1 \frac{\mathrm{d}\alpha}{2} $。以此类推,我们对四个边加和之后,有作用在曲面微元上的表面张力之和的大小为:

\begin{equation}\label{Fsigma} F_\sigma=C\left(\frac{1}{R_1}+\frac{1}{R_2} \right) \mathrm{d}s_1 \mathrm{d}s_2 \end{equation}

其中$R_1$,$R_2$为$\mathrm{d}s_1$,$\mathrm{d}s_2$的主曲率半径。均衡状态下有$F_p = F_\sigma$,即:

\begin{equation}\label{three} p_1 - p_2=C\left(\frac{1}{R_1}+\frac{1}{R_2} \right) \end{equation}

Young-laplacian方程

我们用$\bfn$来表示曲面微元的单位法向,有:

\begin{equation}\label{pn} \frac{\p \bfn}{\p x} \mathrm{d}x= \mathrm{mag}\left(\bfn_x|_{x+\mathrm{d}x} \right) - \mathrm{mag}\left(\bfn_x|_x\right) =\mathrm{sind}\alpha \approx \mathrm{d}\alpha \end{equation}
也即:
\begin{equation} \frac{\partial \bfn}{\partial x} = \frac{\mathrm{d}\alpha}{\mathrm{d} x} = \frac{1}{R_1} \end{equation}
\begin{equation}\label{six} \nabla \cdot \bfn = \frac{\partial \bfn}{\partial x} + \frac{\partial \bfn}{\partial y} + \frac{\partial \bfn}{\partial z} =\left( \frac{1}{R_1}+ \frac{1}{R_2} \right) \end{equation}

将方程\eqref{six}代入到方程\eqref{three}有Young-laplacian方程:

\begin{equation} p_1 - p_2 = C \nabla \cdot \bfn \end{equation}

图2. 流体界面示意图

如图(2)所示,对于CFD中的自由表面,某一点的法向量为$\nabla\alpha$,其为一个连续函数,其在流体1和2区内部为0,在自由表面存在非零值。图(2)自由表面的单位法向量可表示为:

\begin{equation} \bfn=\frac{\nabla \alpha}{|\nabla \alpha|} \end{equation}

进一步的,我们定义界面处的曲率为:

\begin{equation} \kappa = - \nabla \cdot n \end{equation}

CSF模型

如图(3)所示,CSF模型将界面处的压力表示为压力的连续函数:

\begin{equation}\label{ten} p_t=\left\{\begin{matrix} p_b & c=1\\ p_a+C c \nabla \cdot n & 0< c< 0\\ p_a & c=0 \end{matrix}\right. \end{equation}

其中$C$为表面张力系数,$c$为界面处的位置函数,其可表示为$c=\alpha_\mathrm{t} - \alpha_\mathrm{a}$。方程\eqref{ten}将界面处间断的压力函数表示为一个连续函数,这为CSF模型的模型基础。


图3. CSF模型下连续的压力函数

具体的,$c$的定义可以参考图(4):


图4. 两相体系下CSF模型界面处$c$的定义

考虑整个区域内,方程\eqref{ten}可以表示为:

\begin{equation} \nabla p=C \kappa \nabla \alpha \end{equation}

其中$\nabla p$仅在界面处有值。

OpenFOAM中的VOF模型

下面我们来看VOF模型,首先我们有连续性方程以及动量方程:

\begin{equation}\label{cont} \frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \mathbf{U})=0 \end{equation}
\begin{equation}\label{mom} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =- \nabla \mathit{p} + \textit{S} \end{equation}

其中$\tau$为剪切应力,$S$为源项,其进一步可以表示为:

\begin{equation}\label{s} S=F_\sigma+\rho g=C \kappa \nabla \alpha+\rho \bfg \end{equation}

方程\eqref{mom}中的密度$\rho$可以表示为:

\begin{equation} \rho = \alpha_1 \rho_1 + (1-\alpha_1)\rho_2 \end{equation}

把公式\eqref{mom},\eqref{s}结合有:

\begin{equation}\label{momemtum} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot (\rho \mathbf{U} \mathbf{U} ) - \nabla \cdot \tau =- \nabla \mathit{p} + C \kappa \nabla \alpha+\rho \bfg \end{equation}

下面我们来推导相方程,由连续性方程\eqref{cont}我们有:

\begin{equation} \nabla \cdot \mathbf{U}= - \frac{1}{\rho} \frac{\mathrm{D}\rho}{\mathrm{D}t} \end{equation}

假设不可压缩流体有:

\begin{equation} \nabla \cdot \mathbf{U}=0 \end{equation}

也即:

\begin{equation} - \frac{1}{\rho} \frac{\mathrm{D}\rho}{\mathrm{D}t} = - \frac{1}{\rho} \frac{\mathrm{D} \left(\alpha\left(\rho_1 - \rho_2 \right)+\rho_2 \right)}{\mathrm{D}t}=-\frac{\rho_1-\rho_2}{\rho} \frac{\mathrm{D}\alpha}{\mathrm{D}t}=0 \end{equation}

进一步的可表示为:

\begin{equation}\label{alpha} \frac{\mathrm{D}\alpha}{\mathrm{D}t}=\frac{\partial \alpha}{\partial t}+\nabla \cdot (\alpha \mathbf{U})=0 \end{equation}

方程\eqref{cont},\eqref{momemtum},\eqref{alpha}即为不可压缩VOF模型中的连续性方程、动量方程以及相方程。

在interFoam中,方程\eqref{momemtum}需要做进一步的变化,以使得边界条件的定义更加简单(在双流体模型中也可以进行类似的数学操作)。首先定义

\begin{equation}\label{prgh} p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh \end{equation}

其中$\bfh$表示网格单元体心的位置矢量。对方程\eqref{prgh}进行梯度操作:

\begin{equation}\label{nablaprgh} \nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg \end{equation}

把公式\eqref{nablaprgh}回代到公式\eqref{momemtum}中有最终的动量方程:

\begin{equation}\label{momentum2} \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} \end{equation}

需要注意的是,VOF中的相分数可以看作为一种被动的标量(Passive Scalar),其需要严格有界与$0-1$之间。在求解方程\eqref{alpha}的时候,需要更进一步的处理,请参考interPhaseChangeFoam解析

代码分析

首先进入createFields.H文件:

    Info<< "Reading transportProperties\n" << endl;
    incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
    //创建incompressibleTwoPhaseMixture类型,并称为twoPhaseProperties

    volScalarField& alpha1(twoPhaseProperties.alpha1());
    //声明alpha1为一个volScalarField引用
    volScalarField& alpha2(twoPhaseProperties.alpha2()); 
    //声明alpha2为一个volScalarField引用

    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
    //声明rho1为一个dimensionedScalar的常量引用
    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();


    // Need to store rho for ddt(rho, U)
    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT
        ),
        alpha1*rho1 + alpha2*rho2, //公式(15)
        alpha1.boundaryField().types()
        //rho的边界条件和alpha1的边界条件相同
    );
    rho.oldTime(); //存储rho的旧值

    ...

    // Construct interface from alpha1 distribution
    interfaceProperties interface(alpha1, U, twoPhaseProperties);
    //创建interfaceProperties类型,并命名为interface

    // Construct incompressible turbulence model
    autoPtr<incompressible::turbulenceModel> turbulence
    (
        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
    );
    //创建湍流

    #include "readGravitationalAcceleration.H"
    //创建重力场,下面省略有关重力场的内容,请参考twoLiquidMixingFoam解析中的内容

然后我们看主要的函数twoPhaseProperties.correct()是怎样执行的。首先,correct()函数来自twoPhaseProperties,twoPhaseProperties是一个incompressibleTwoPhaseMixture类型,在incompressibleTwoPhaseMixture.H中,其定义:

    //- Correct the laminar viscosity
    virtual void correct()
    {
        calcNu();
    }

再看calcNu()函数:

//- Calculate and return the laminar viscosity
void Foam::incompressibleTwoPhaseMixture::calcNu()
{
    nuModel1_->correct();
    //对第一相计算粘度,interFoam支持非牛顿流体以及牛顿流体
    nuModel2_->correct();
    //对第二相计算粘度

    const volScalarField limitedAlpha1
    (
        "limitedAlpha1",
        min(max(alpha1_, scalar(0)), scalar(1))
    );
    //对alpha进行有界限定在0-1之间

    // Average kinematic viscosity calculated from dynamic viscosity
    nu_ = mu()/(limitedAlpha1*rho1_ + (scalar(1) - limitedAlpha1)*rho2_);
    //公式(20)
}

其通过对俩相计算运动粘度来更新平均运动粘度:

\begin{equation} \nu=\frac{\mu}{\alpha_1 \rho_1 + \alpha_2 \rho_2} \end{equation}

其中的$\mu$为俩相的平均动力粘度,其这样计算:

\begin{equation} \mu=\alpha_1 \rho_1 \nu_1 + \alpha_2 \rho_2 \nu_2 \end{equation}

回到interFoam.C文件有:

    #include "alphaEqnSubCycle.H"//参考twoLiquidMixingFoam解析
    interface.correct();

同样的,分析

    void correct()
    {
        calculateK();
    }

进入interfaceProperties.C文件有:

void Foam::interfaceProperties::calculateK()
{
    const fvMesh& mesh = alpha1_.mesh();
    const surfaceVectorField& Sf = mesh.Sf();
    //面法向矢量

    // Cell gradient of alpha
    const volVectorField gradAlpha(fvc::grad(alpha1_));
    //∇α1,计算alpha的梯度

    // Interpolated face-gradient of alpha
    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
    //∇α1,f 对alpha的梯度在面上插值

    // Face unit interface normal
    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));//公式(8)
    correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());

    // Face unit interface normal flux
    nHatf_ = nHatfv & Sf;
    //将面矢量n乘以网格面矢量

    // Simple expression for curvature
    K_ = -fvc::div(nHatf_);
    //公式(9),fvc::div()的处理相同于icoFoam中的fvc::div(phiHbyA),可参考OpenFOAM对标量求散度?
}

下面看我们的速度方程:

fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) //公式(23)第1项
      + fvm::div(rhoPhi, U) //公式(23)第2项
      + turbulence->divDevRhoReff(rho, U) //公式(23)第3项
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                  //公式(23)第4项
                  - ghf*fvc::snGrad(rho)
                  //公式(23)第5项
                  - fvc::snGrad(p_rgh)
                  //公式(23)第6项
                ) * mesh.magSf()
            )
        );
    }

组建速度方程之后,进入压力方程,求解思想请参考其他求解器解析。需要提及的是,在求解压力p_rgh之后,需要通过代码:

    p == p_rgh + rho*gh;

来更新压力p。

2018.04.20依据齐雪宇的建议修正方程\eqref{prgh}、\eqref{nablaprgh}中的$\bfg \bfh$为$\bfg\cdot\bfh$ | 2018.04.08增加fvc::div()的讨论 | 2018.04.06重新排版,小修内容

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