Code of the Week #21: MULES

OpenFOAM
1

#201908

MULES (Multidimensional Universal Limiter with Explicit Solution)，

\frac{\partial \alpha}{\partial t} + \nabla\cdot(\alpha \boldsymbol{U})=0 \tag{1}

利用 Gauss 定理，把上式写成积分形式，有

\int_{\Omega}\frac{\partial \alpha}{\partial t} \mathrm{d}V + \int_{\partial\Omega}\alpha\boldsymbol{U}\cdot n\mathrm{d} S =0 \tag{2}

其在网格上的离散形式为

\frac{\alpha_i^{n+1}-\alpha_i^n}{\Delta t}=\frac{-1}{\Omega_i}\sum_{f\in\partial\Omega_i}(F_u+\lambda_M F_c)^n \tag{3}

上式的左边为时间的离散，右边为对流项的离散，写成该网格每个面上的通量之和。其中，

F_u = \phi_f \alpha_{f,\text{upwind}} \tag{4}

为低阶格式的离散通量， F_c 为

F_c = \phi_f \alpha_f +\phi_{rf} \alpha_{rf} (1-\alpha)_{rf} - F_u \tag{5}

F_c 的第一项是高阶格式的通量，可以看到， \lambda_M 是一个控制因子，在非界面的网格上， \lambda_M=0 ，在界面所在的网格上， \lambda_M=1 ，

当 \lambda_M=0 时，通量为 F_u ，当 \lambda_M=1 时，为 \phi_f \alpha_f +\phi_{rf} \alpha_{rf} (1-\alpha)_{rf} ，第二项 \phi_{rf} \alpha_{rf} (1-\alpha)_{rf} 是界面压缩通量 (interfacial compression flux)，这样处理的目的就是降低界面的数值耗散，并减小计算 cost，把高阶格式限定在界面上。

而这一段的程序实现可见 src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,                 // alpha
    const surfaceScalarField& phi,       // U_f * S_f
    surfaceScalarField& phiPsi,          // Fc 前两项
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    Info<< "MULES: Solving for " << psi.name() << endl;
    const fvMesh& mesh = psi.mesh();

    psi.correctBoundaryConditions();

    surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi); // F_u

    surfaceScalarField& phiCorr = phiPsi;            // 高阶的部分
    phiCorr -= phiBD;                                // 把低阶的部分拿掉

    scalarField allLambda(mesh.nFaces(), 1.0);       // 预设 lamda_M

    slicedSurfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimless,
        allLambda,
        false   // Use slices for the couples
    );

    limiter
    (
        allLambda,
        rho,
        psi,
        phiBD,
        phiCorr,
        Sp.field(),
        Su.field(),
        psiMax,
        psiMin,
        3
    );

    phiPsi = phiBD + lambda*phiCorr;         // Fu + lamda_M * Fc

    scalarField& psiIf = psi;
    const scalarField& psi0 = psi.oldTime();
    const scalar deltaT = mesh.time().deltaT().value();

    psiIf = 0.0;
    fvc::surfaceIntegrate(psiIf, phiPsi);

    if (mesh.moving())
    {
        psiIf =
        (
            mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
          + Su.field()
          - psiIf
        )/(rho/deltaT - Sp.field());
    }
    else
    {
        psiIf =
        (
            rho.oldTime()*psi0/deltaT
          + Su.field()
          - psiIf
        )/(rho/deltaT - Sp.field());
    }

    psi.correctBoundaryConditions();
}

limiter() 函数就在这个文件里，实现对 \alpha\in[\alpha_{\min},\alpha_{\max}] 、 \lambda_M 的限制。 F_c 的前两项在 alphaEqn.H 中定义，

        surfaceScalarField phiAlpha =
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
                alpha1,
                alpharScheme
            );

phir 为

    surfaceScalarField phic = mag(phi/mesh.magSf());
    phic = min(interface.cAlpha()*phic, max(phic));
    surfaceScalarField phir = phic*interface.nHatf();

alphaScheme 定义了离散格式

    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");
4

image 前辈，您好，能谈一下，这些参数的具体用法嘛，我在处理比较复杂的案例时在确定网格和时间步以及边界条件没有问题的情况下，怎么调整这些参数呢，