[romant]

Hej,

I have been working with the interPhaseChangeFoam solver for a while now, and have started writing a phase change model based on temperature. Right now the temperature is still given via the dictionary input.

I wanted to add the energy equation for the fluid phase to this solver. The governing equation is


which I tried implementing like this

{

//    const volScalarField rhoAlpha = cp * rho1 * alpha1; // for time derivative
    volScalarField rhoAlpha
    (
        "rhoAlpha",
        cp * rho1 * alpha1
    );


    surfaceScalarField phiAlpha
    (
        "phiAlpha",
        phi * rho1 * cp  * fvc::interpolate(alpha1) //check with rhoPhi from the alphaEqn.H file
//        rhoPhi * cp //* fvc::interpolate(alpha1)
    );

    volScalarField lambdaFAlpha
    (
        "lambdaFAlpha",
        lambdaF * alpha1
    );


    Pair<tmp<volScalarField> > mDotAlphal = twoPhaseProperties->mDotAlphal();
    const volScalarField mDotAlphac = mDotAlphal[0]();

    fvScalarMatrix TEqn
    (
          fvm::ddt(rhoAlpha,T)
        + fvm::div(phiAlpha, T)
        - fvm::laplacian(lambdaFAlpha, T)
       == mDotAlphac * ifg

    );


    TEqn.relax();
    TEqn.solve();

    Info << "Max(T): " << max(T).value() << " K "
    << "Min(T): " << min(T).value() << " K" << endl;

}

Unfortunately, I have had no success in actually using this within the solver, the results are completely bogus. One thing that is obvious is that there might be a problem with the alpha1 attachment to the equation.

My first thought was to maybe be able to extract the mesh that contains water alpha1>0, and impose boundary conditions onto the interface cells. Is there a way to do it, or does somebody else see what could be wrong with the implementation of this equation?