返回主页

twoLiquidMixingFoam解析

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

首先进入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}\label{eight} \frac{\partial \rho u}{\partial t}+\nabla \cdot (\rho uu) - \nabla \cdot(\nu_t \nabla \rho u)=-g \cdot h \nabla \rho - \nabla p_{rgh} \end{equation} 其中方程的右边项请参考interFoam解析。组建速度矩阵之后则进入压力方程求解,求解思路请参考icoFoam解析,不再赘述。

更新历史

2018.03.30:依据Dingcy的建议更正方程\eqref{eight}

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