[Phigo90]

Hallo everyone,


I am working on the InterCondensatingEvaporatingFoam, which uses the "alphaEqnSubCycle.H" and "alphaEqn.H" from the InterphaseChangeFoam.


The original solver is written like follows:

        fvScalarMatrix alpha1Eqn
        (
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)
          - fvm::Sp(divU, alpha1)
         ==
            fvm::Sp(vDotvmcAlphal, alpha1)
            + vDotcAlphal
        );

In mathematical notation, this can be written as:


d(alpha)/dt+v*div(alpha)=(S_vap-S_con)*alpha+S_con


My frist question is about the RHS.
The expression can be rewritten as:
(S_vap-S_con)*alpha+S_con=S_vap*alpha+S_con(1-alpha).


However, I want to change the linear dependency on the right HS like:
S_vap*sqrt(alpha)+S_con(1-sqrt(alpha)).


Does anyone of you know, how it is possible to change the source Code to use a non linear source term? The source term will decrease the values of alpha, so I think this is the reason why it is treated in an implicit way.



The second question is about MULEs.


Original Code:

    volScalarField alpha10("alpha10", alpha1);

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        tmp<surfaceScalarField> talphaPhiCorr
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

        if (MULESCorr)
        {
            talphaPhiCorr.ref() -= talphaPhi();
            volScalarField alpha100("alpha100", alpha10);
            alpha10 = alpha1;

            MULES::correct
            (
                geometricOneField(),
                alpha1,
                talphaPhi(),
                talphaPhiCorr.ref(),
                vDotvmcAlphal,
                (
                    divU*(alpha10 - alpha100)
                  - vDotvmcAlphal*alpha10
                )(),
                1,
                0
            );

MULES is called after solving the equation mentioned above to guarantee boundedness. I think one problem is, that I do not completly understand what MULES is doing, but in my opinion it doesn't work. If I get for example a value for alpha=-1*e^-5 (which is lower than the physical limit of alpha, which is in the region of 0 to 1), the value after running MULES does not change.