Problem adding a source term in interfoam's alphaEqn.H
 

Hi foamers,
Related to this thread concerning "numerical ventilation" with interFoam, i'n now trying to add a source term to the alpha equation in interFoam.

As I don't perfectly understand how to handle implicit/explicit source term and how it is treated in MULES::correct and MULES::explicitsolve, i try to impelment my source term totally explicitely.

My source term would have this form (with an arbitrary distance from the wall, the distance from the cell center to the wall, an arbitrary number of iterations and the water phase fraction of the considered cell) :


if

if //mostly air in the cell / above interface

//i'm removing water from the cell

else //mosty water in the cell /under interface

//i'm adding water to the cell

else

I try to implement this in interFoam solver :

1. calculating the wall distance in createField.H
   Code:

   ---

   Info<< "Calculating wall distance field" <<endl;
    volScalarField y = wallDist(mesh).y();

   ---

2. creating the source terms in alphaEqn.H
   Code:

   ---

   //Creating Source Terms
    volScalarField::DimensionedInternalField Su
    (
        IOobject
        (
            "Su",
            runTime.timeName(),
            mesh
        ),
    mesh,
        dimensionedScalar("Su", dimensionSet(0,0,-1,0,0,0,0), scalar(0.0))
    );

   
    volScalarField::DimensionedInternalField Source
    (
        IOobject
        (
            "Source",
            runTime.timeName(),
            mesh
        ),
    mesh,
        dimensionedScalar("Source", dimensionSet(0,0,-1,0,0,0,0), scalar(0.0))
    );

   ---

3. evaluating the source term in alphaEqn.H ; as i want to getan only explicitly handled source term, i let Su to 0.
   Code:

   ---

   //evaluating Source Term
    forAll( mesh.C(), celli)
    {
        if(y[celli] <= 0.02)
        {
            if(alpha1[celli] > 0.5)
            {
                Source[celli] = ( 1.0 - alpha1[celli])/10.0;
            }
            else
            {
                Source[celli] = - alpha1[celli]/10.0;
            }
        }
        else
        {
            Source[celli] = 0.0;
        }
    }

   ---

4. adding the source term the equation
   Code:

   ---

   fvScalarMatrix alpha1Eqn
        (
            #ifdef LTSSOLVE
            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
            #else
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            #endif
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)
          ==
            Source
        );
   
        alpha1Eqn.solve();

   ---

5. adding the source terms in MULES::correct and MULES::explicitSolve (I don't understand how it is working, so i just copied the function calls of interPhaseChangeFoam ...). there are three call of those MULES fonction in the code, i just quote the first one of them :
   Code:

   ---

   volScalarField alpha10(alpha1);
            #ifdef LTSSOLVE
            MULES::LTScorrect
            (
                geometricOneField(),
                alpha1,
                tphiAlpha(),
                tphiAlphaCorr0(),
                Source,
                (-Su*alpha10)(),
                1,
                0
            );
            #else
            MULES::correct
            (
                geometricOneField(),
                alpha1,
                tphiAlpha(),
                tphiAlphaCorr0(),
                Source,
                (-Su*alpha10)(),
                1,
                0
            );
            #endif

   ---

This is compiling correctly ; but it seems that there is no difference between this solver and the classic LTSInterFoam one. I changed a bit my code to be able to see my Source term in paraview and it is calculated.


I don't know where I am wrong - in fact if I do I would not ask fo some help.


Any Idea would be very welcome.

I read quite a lot of threads on this forum related to the source term, but didn't get any more clue concerning my implementation.

solved
 

Solution here.