[Daniel_Khazaei]

Dear foamers,

I'm working on a problem dealing with evaporation and condensation of water inside a closed domain. The fixedFluxPressure boundary condition is applied on all the boundaries which necessitate the need for setting a reference pressure before solving the pdEqn.

The evaporation-condensation model adds a mass transfer source in pdEqn as follow:

        fvScalarMatrix pdEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, pd_)
          - (vDotc - vDotv)
        );

The problem is the reference pressure location has an unexpected effect on the simulation which seems totally unrealistic. If I remove the source term from the equation, then everything becomes completely normal. This finally cause the temperature field to become unbounded and drops way below physical temperature range allowed inside the domain. The contours of velocity and temperature shows the problem:





If it helps further investigating the issue, here is the full pEqn function:

void interPhaseChangeFluid::pEqn
(
    const fvMesh& mesh, const fvVectorMatrix& UEqn
)
{
    rAU_ = 1.0/UEqn.A();
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU_));

    volVectorField HbyA("HbyA", U_);
    HbyA = rAU_*UEqn.H();

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU_, rho_, U_, phi_)
    );

    #include "updateRobinFsiInterface.H"

    adjustPhi(phiHbyA, U_, pd_);

    //- calculating surface tension flux
    //- + gravitational flux
    surfaceScalarField phig
    (
        interface_.capillaryFlux(rAUf)
      - ghf_*fvc::snGrad(rho_)*rAUf*mesh.magSf()
    );

    //- adding the resulting flux to the main flux
    phiHbyA += phig;

    Pair<tmp<volScalarField> > vDot = phaseChangeModel_->vDot();
    const volScalarField& vDotc = vDot[0]();
    const volScalarField& vDotv = vDot[1]();

    while (pimple_.correctNonOrthogonal())
    {
        fvScalarMatrix pdEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, pd_)
          - (vDotc - vDotv)
        );

        pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

        pdEqn.solve
        (
            mesh.solutionDict().solver(pd_.select(pimple_.finalInnerIter()))
        );

        if (pimple_.finalNonOrthogonalIter())
        {
            phi_ = phiHbyA + pdEqn.flux();

            U_ = HbyA + rAU_*fvc::reconstruct((phig + pdEqn.flux())/rAUf);
            U_.correctBoundaryConditions();
            gradU_ = fvc::grad(U_);
        }
    }

    //- re-evaluate pressure using additional terms
    p_ = pd_ + rho_*gh_;
    p_.correctBoundaryConditions();
    gradp_ = fvc::grad(p_);

    if (pd_.needReference())
    {
        p_ += dimensionedScalar
        (
           "p",
            p_.dimensions(),
            pRefValue_ - getRefCellValue(p_, pRefCell_)
        );
        pd_ = p_ - rho_*gh_;
    }

    //- Make the fluxes relative to the mesh motion
    fvc::makeRelative(phi_, U_);
}

Regards,
D.Khazaei