• 返回主页
  • twoLiquidMixingFoam解析

    Reviewed by:
    求解器类型:俩种不可压缩可溶液体的混合流动。

    首先进入createFields.H,在createFields.H中,定义:

    Info<< "Reading transportProperties\n" << endl;
    incompressibleTwoPhaseMixture mixture(U, phi);
    

    其声明一个incompressibleTwoPhaseMixture类型,并称为mixtureUphi传入构造函数。

    dimensionedScalar Dab("Dab", dimViscosity, mixture);
    
    // Read the reciprocal of the turbulent Schmidt number
    dimensionedScalar alphatab("alphatab", dimless, mixture);
    
    // Need to store rho for ddt(rho, U)
    volScalarField rho("rho", alpha1*rho1 + alpha2*rho2);
    rho.oldTime();
    

    然后,定义一个带量纲的标量(dimensionedScalarDab以及alphatab,其值为mixture,但是mixture为一个incompressibleTwoPhaseMixture类型。实际上incompressibleTwoPhaseMixture类型中定义了一个字典,因此dimensionedScalar则会从算例的字典中进行读取。

    dimensionedScalar Dab("Dab", dimViscosity, mixture);
    
    // Need to store rho for ddt(rho, U)
    volScalarField rho("rho", alpha1*rho1 + alpha2*rho2);
    rho.oldTime();
    

    接下来我们定义密度rho,其中rho.oldTime()函数执行后会存储rho的旧值。

    #include "readGravitationalAcceleration.H"
    #include "readhRef.H"
    #include "gh.H"
    

    我们把头文件展开有下述代码:

        Info<< "\nReading g" << endl;
        uniformDimensionedVectorField g
        (
            IOobject
            (
                "g",
                runTime.constant(),
                mesh,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        );//创建重力场
    
        Info<< "\nReading hRef" << endl;
        uniformDimensionedScalarField hRef
        (
            IOobject
            (
                "hRef",
                runTime.constant(),
                mesh,
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE
            ),
            dimensionedScalar("hRef", dimLength, 0)
        );//创立参考高度
        Info<< "Calculating field g.h\n" << endl;
        dimensionedScalar ghRef
        (
            mag(g.value()) > SMALL
          ? g & (cmptMag(g.value())/mag(g.value()))*hRef
          : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
        );//如果g大于SMALL(1e-15),则将ghRef初始化,否则为0。假设g为(0,0,-9.81),g.value()为(0,0,9.81),cmptMag(g.value())为(0,0,9.81)。
    mag(g.value())为9.81。
    volScalarField gh("gh", (g & mesh.C()) - ghRef); surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);

    至此,createFields.H分析完毕。下面我们进入alphaCourantNo.H,

    scalar maxAlphaCo
    (
        readScalar(runTime.controlDict().lookup("maxAlphaCo"))
    );//从字典中读取manAlphaCo
    
    scalar alphaCoNum = 0.0;
    scalar meanAlphaCoNum = 0.0;
    
    if (mesh.nInternalFaces())
    {
        scalarField sumPhi
        (
            pos(alpha1 - 0.01)*pos(0.99 - alpha1)
           *fvc::surfaceSum(mag(phi))().internalField()
        );//公式(1)
    
        alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();//公式(2)
    
        meanAlphaCoNum =
            0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();//公式(3)
    }
    
    Info<< "Interface Courant Number mean: " << meanAlphaCoNum
        << " max: " << alphaCoNum << endl;
    

    首先定义sumPhi: \begin{equation} \mathrm{sumPhi}=\left\{\begin{matrix} \sum |\phi_f| &0.01<\alpha<0.99 \\ 0, &\mathrm{else} \end{matrix}\right. \end{equation} 然后定义最大相库朗数alphaCoNum: \begin{equation} \mathrm{alphaCoNum}=0.5*\mathrm{max}\left( \frac{\sum |\phi_f|}{V} \Delta t \right) \end{equation} 然后定义平均相库朗数meanAlphaCoNum: \begin{equation} \mathrm{meanAlphaCoNum}=0.5* \frac{\sum (\sum |\phi_f|)}{\sum V} \Delta t \end{equation} 然后通过#include "setDeltaT.H"设定库朗数。在设定库朗数之后,执行:

    Info<< "Time = " << runTime.timeName() << nl << endl;
    
    mixture.correct();
    

    其中的mixture.correct()函数请参考incompressibleTwoPhaseMixture.H,此函数通过速度场来返回粘度。如果为牛顿流体则返回流体的物理粘度,如果为非牛顿流体则对粘度进行计算。下面,进入相方程求解部分:

    #include "alphaEqnSubCycle.H"//求解alpha方程的时候,分割时间步,比如:每个时间步压力求解一次,但是alpha求解两次,相应地时间步减半
    #include "alphaDiffusionEqn.H"
    

    展开alphaEqn.H有:

    {
        word alphaScheme("div(phi,alpha)");//定义alphaScheme,可以调用通量限制器离散格式
    
        surfaceScalarField alphaPhi
        (
            phi.name() + alpha1.name(),
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
        );//公式(4)
    
        MULES::explicitSolve(alpha1, phi, alphaPhi, 1, 0);//通过MULES求解相方程,MULES求解理论是OpenFOAM基金会研发的一种有效的解决相场有界性的求解方法。
    对应公式(6)
    rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;//公式(5) Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; }

    其中的alphaPhi定义为: \begin{equation} \mathrm{alphaPhi}=\alpha_1 u \end{equation} 其中的rhoPhi定义为(整理后即对应代码中的定义): \begin{equation} \mathrm{rhoPhi}=u (\alpha_1 \rho_1 + \alpha_2 \rho_2) \end{equation} MULES求解的相方程对应于: \begin{equation} \frac{\partial \alpha_1}{\partial t}+\nabla \cdot (\alpha_1 u) = 0 \end{equation} 然后进入alphaDiffusionEqn.H,求解相场的扩散方程:

    {
        fvScalarMatrix alpha1Eqn
        (
            fvm::ddt(alpha1)
          - fvc::ddt(alpha1)
          - fvm::laplacian
            (
                volScalarField("Dab", Dab + alphatab*turbulence->nut()),
                alpha1
            )
        );//公式(7)
    
        alpha1Eqn.solve();
    
        alpha2 = 1.0 - alpha1;//alpha_2通过alpha_1减得
        //rhoPhi += alpha1Eqn.flux()*(rho1 - rho2);
        rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
    }
    
    rho = alpha1*rho1 + alpha2*rho2;
    

    其中相场的扩散方程定义为: \begin{equation} \frac{\partial \alpha_1}{\partial t} - \nabla \cdot \left( \left( \nu_m + \nu_t \right)u \right)= 0 \end{equation} 其中$\nu_m$为分子扩散系数,其对应为Dab,$\nu_t$为湍流粘度,alphatab为一个用户自定义的常数。

    在求得相场之后,就可以进行NS方程求解了。首先我们组建速度矩阵:

        fvVectorMatrix UEqn
        (
            fvm::ddt(rho, U)
          + fvm::div(rhoPhi, U)
          + turbulence->divDevRhoReff(rho, U)
        );//公式(8)左边
    
        UEqn.relax();
    
        if (pimple.momentumPredictor())
        {
            solve
            (
                UEqn
             ==
                fvc::reconstruct
                (
                    (
                      - ghf*fvc::snGrad(rho)
                      - fvc::snGrad(p_rgh)
                    ) * mesh.magSf()
                )
            );////公式(8)右边
        }
    

    其对应速度方程: \begin{equation} \frac{\partial \rho u}{\partial t}+\nabla \cdot (\rho uu) - \nabla \cdot(\nu_t \nabla \rho u)=-g \cdot h \nabla p - \nabla p_{rgh} \end{equation} 其中方程的右边项请参考interFoam解析。组建速度矩阵之后则进入压力方程求解,求解思路请参考icoFoam解析,不再赘述。
    如何引用:
    李东岳. twoLiquidMixingFoam解析.[OL]. http://dyfluid.com/twoLiquidMixingFoam.html, 2016.3.10

    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国