• 返回主页
  • interFoam解析

    2017.07.19:添加有关相方程的连接,小修格式

    VOF模型推导

    在OpenFOAM中,VOF模型被植入到了interFoam中,用于处理表面跟踪(界面捕获)问题。例如:溃坝、极少量的气泡界面捕获、液体晃荡等问题。本文从最基本的液膜压力降开始分析,一步一步推导VOF模型。

    VOF模型中定义了一个示踪函数,$\alpha$,其在流体1处为1,流体2处为0,0到1之间的值则表示自由表面。可见,CFD后处理中VOF的界面并不是直接算出来的,而是通过计算每个网格内的相场值而后处理出来的。因此,网格越密,界面越薄。我们从下图的小界面微元开始分析:我们定义表面张力,$\sigma$,是一种为了保持界面平衡而作用在每单位长度上的力,它是一种和界面相切的界面张力,表面张力的大小主要和流体的特性有关。在曲面中,均衡的流体表面张力和界面压力降,$\Delta p$,均衡。这个界面压力降主要取决于表面张力和曲面弧度。


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

    如图(1)所示的一个曲面微元。其中$p_1$为曲面微元上方液体向下施加的力,$p_2$为曲面微元下方液体向上施加的力,$s_1$和$s_2$分别表示曲面微元的俩个边长。那么我们有曲面微元上的压力为: \begin{equation} F_p=(p_1 - p_2) \mathrm{d}s_1 \mathrm{d}s_2 \end{equation} 基于表面张力的定义,其作用在曲面微元的四个边上,且合力向上和曲面微元上的压力均衡。如果表面张力系数定义为:$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} 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} VOF中的Young-laplacian方程:我们用$n$来表示曲面微元的单位法向。我们有: \begin{equation} \frac{\partial n}{\partial x} \mathrm{d}x= \mathrm{mag}\left(n_x|_{x+\mathrm{d}x} \right) - \mathrm{mag} \left(n_x|_x\right) =\mathrm{sind}\alpha \approx \mathrm{d}\alpha \end{equation} 也即: \begin{equation} \frac{\partial n}{\partial x} = \frac{\mathrm{d}\alpha}{\mathrm{d} x} = \frac{1}{R_1} \end{equation} \begin{equation}\label{six} \nabla \cdot n = \frac{\partial n}{\partial x} + \frac{\partial n}{\partial y} + \frac{\partial n}{\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 n \end{equation} 在VOF中,我们用$\alpha$表示相的体积分数。当$\alpha=1$,表示此控制体的第一相的体积分数为1,当$\alpha=0$,表示此控制体的第一相的体积分数为0。当$0<\alpha<1$的时候,表示处于界面状态。


    图2. 流体界面示意图

    如图(2)所示,对于其中的自由表面,某一点的法向量为$\nabla\alpha$,其为一个连续函数,其在流体1、2区为0,在自由表面存在值。上图中,自由表面的单位法向量可表示为: \begin{equation} n=\frac{\nabla \alpha}{|\nabla \alpha|} \end{equation} 进一步的,我们定义界面处的曲率为: \begin{equation} \kappa = - \nabla \cdot n \end{equation} CSF模型: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} 如图(3)所示:其中$c$为界面处的位置函数,其可表示为$c=\alpha_\mathrm{t} - \alpha_\mathrm{a}$,$C$为表面张力系数。方程\eqref{ten}将界面处间断的压力函数表示为一个连续函数。


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

    关于$c$的定义可以参考图(4):


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

    考虑整个区域内,方程\eqref{ten}可以表示为: \begin{equation} \nabla p=C \kappa \nabla \alpha \end{equation} 其中$\nabla p$仅在界面处有值。

    CFD中的VOF

    下面我们来看VOF模型,我们有连续性方程: \begin{equation} \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 g \end{equation} 方程\eqref{mom}中的密度$\rho$可以表示为: \begin{equation} \rho = \alpha_1 \rho_1 + (1-\alpha_1)\rho_2 \end{equation} 下面我们来推导相方程,由连续性方程我们有: \begin{equation} \nabla \cdot \mathbf{U}= - \frac{1}{\rho} \frac{\mathrm{D}\rho}{\mathrm{D}t} \end{equation} 假设不可压缩流体有: \begin{equation}\label{cont} \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{mom}和\eqref{s},\eqref{alpha}即为不可压缩VOF模型中的连续性方程、动量方程以及相方程。需要注意的是,VOF中的相分数可以看作为一种被动的标量(Passive Scalar),其需要严格有界与$0-1$之间。在求解方程\eqref{alpha}的时候,需要更进一步的处理,请参考:interPhaseChangeFoam解析

    代码分析

    下面我们进入interFoam求解器分析,首先进入interFoam.C文件:

    #include "fvCFD.H"
    #include "MULES.H"//使用MULES求解
    #include "subCycle.H"//调用子时间步,参考twoLiquidMixingFoam
    #include "interfaceProperties.H"//调用其中定义的界面张力相关内容
    #include "incompressibleTwoPhaseMixture.H"//alpha方程在
    #include "turbulenceModel.H"//湍流
    #include "pimpleControl.H"//pimple算法
    #include "fvIOoptionList.H"//源项
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"//算例目录
        #include "createTime.H"//创建时间对象
        #include "createMesh.H"//创建网格对象
    
        pimpleControl pimple(mesh);//pimple控制
    
        #include "initContinuityErrs.H"//连续性误差
        #include "createFields.H"
    

    至此,我们进入createField.H文件(略去部分声明场的内容):

        Info<< "Reading transportProperties\n" << endl;
        incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);//创建incompressibleTwoPhaseMixture类型,并称为twoPhaseProperties
    
        volScalarField& alpha1(twoPhaseProperties.alpha1());//声明alpha1引用
        volScalarField& alpha2(twoPhaseProperties.alpha2());//声明alpha2引用
    
        const dimensionedScalar& rho1 = twoPhaseProperties.rho1();//rho1引用,const类型,表明不会更改rho1的值,有关何时使用const,
    何时不使用const,请参考《C++ Primer Plus》
    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();//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解析中的内容
    

    回到interFoam.C代码:

        #include "readTimeControls.H"//读取时间控制
        #include "correctPhi.H"//如果从动网格算例结果初始场,通过pcorr进行通量修正,对于interFoam,可忽略
        #include "CourantNo.H"//库朗数计算
        #include "setInitialDeltaT.H"//设置时间步
    
    
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
        Info<< "\nStarting time loop\n" << endl;
    
        while (runTime.run())
        {
            #include "readTimeControls.H"
            #include "CourantNo.H"
            #include "alphaCourantNo.H"
            #include "setDeltaT.H"
    
            runTime++;
    
            Info<< "Time = " << runTime.timeName() << nl << endl;
    
            twoPhaseProperties.correct();//进入twoPhaseProperties.correct()函数
    

    下面我们看twoPhaseProperties.correct()是怎样执行的。首先,correct()函数来自twoPhasePropertiestwoPhaseProperties是一个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();
    

    同样的,我们分析interface.correct()函数。进入interfaceProperties.H文件,我们有:

        void correct()
        {
            calculateK();
        }
    

    进入interfaceProperties.C文件,我们有:

    void Foam::interfaceProperties::calculateK()
    {
        const fvMesh& mesh = alpha1_.mesh();
        const surfaceVectorField& Sf = mesh.Sf();//面法向矢量
    
        // Cell gradient of alpha,计算alpha的梯度
        const volVectorField gradAlpha(fvc::grad(alpha1_));//$\nabla \alpha_1$
    
        // Interpolated face-gradient of alpha,对alpha的梯度再面上插值
        surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));//$\nabla \alpha_{1,f}$
    
        // 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转化为面标量n
    
        // Simple expression for curvature
        K_ = -fvc::div(nHatf_);//公式(9),fvc::div()的处理相同于icoFoam中的fvc::div(phiHbyA)
    
    }
    

    回到VOF的动量方程,把公式(13),(14)结合有: \begin{equation} \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 g \end{equation} 在interFoam中通常我们求解大名鼎鼎的p_rgh,这样可以使得边界条件的定义更加简单(在最新版本的twoPhaseEulerfoam中,Henry表示在非正交网格是求解p_rgh更加稳定)。我们定义: \begin{equation} \mathrm{p_{rgh}}=p-\rho g h \end{equation} 也即: \begin{equation} \nabla \mathrm{p_{rgh}}=\nabla p-g h \nabla \rho - \rho g \end{equation} 把公式(24)回代到公式(22)中,我们有动量方程: \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 \mathrm{p_{rgh}} \end{equation} 下面看我们的速度方程:

        fvVectorMatrix UEqn
        (
            fvm::ddt(rho, U)//公式(25)第1项
          + fvm::div(rhoPhi, U)//公式(25)第2项
          + turbulence->divDevRhoReff(rho, U)//公式(25)第3项
        );
    
        UEqn.relax();
    
        if (pimple.momentumPredictor())
        {
            solve
            (
                UEqn
             ==
                fvc::reconstruct
                (
                    (
                        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)//公式(25)第4项
                      - ghf*fvc::snGrad(rho)//公式(25)第5项
                      - fvc::snGrad(p_rgh)//公式(25)第6项
                    ) * mesh.magSf()
                )
            );
        }
    

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

    p == p_rgh + rho*gh;
    

    来更新压力p
    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国