CFD Online Discussion Forums - Problem with Diffusion Equation Implementation

Hi all,

I am using OpenFOAM v1706 and trying to make some additions to the interFlow solver.

I have initialized a field that is concentrated at the interfacial cells of my two-phase system. It is shown by the attached image (InitialCondition.jpg). I use a simple diffusion equation to diffuse this field on the rest of the grid. The result is shown in the second image (Diffused.jpg). Clearly, the circular symmetry is lost and it has to be due to some error in my implementation of the diffusion equation. But the volume integration of the initial source is conserved to the tolerance level that was specified.

Here is my source code used for the diffusion equation

spreadSource.H

Code:

const dimensionedScalar mdotDiffusionConstant("diff",dimArea,dict_.lookupOrDefault<scalar>("mdotDiffusionConstant",1000));



mdot1 = mdot0; // mdot0 stores the initial condition

fvScalarMatrix mdotSpreadEqn

(

    fvm::laplacian(mdotDiffusionConstant, mdot1)

);

mdotSpreadEqn.solve();

The mdot fields were created as follows.

createFields.H

Code:

// Calculating Source Term at Interfacial Cells

volScalarField mdot0

(

    IOobject

    (   

        "mdot0",

        runTime.timeName(),

        mesh,

        IOobject::NO_READ,

        IOobject::AUTO_WRITE

    ),  

    mesh,

    dimensionedScalar("zero",dimensionSet(1,-3,-1,0,0,0,0),0.0)

);



// Distributing Source Term in the domain

volScalarField mdot1

(

    IOobject

    (   

        "mdot1",

        runTime.timeName(),

        mesh,

        IOobject::NO_READ,

        IOobject::AUTO_WRITE

    ),  

    mesh,

    dimensionedScalar("zero",dimensionSet(1,-3,-1,0,0,0,0),0.0),

    mdotBoundaries

);

Here is the fvSolution file I am using:

fvSolution

Code:

solvers

{

    "alpha.water.*"

    {

        interfaceMethod "isoAdvector";

        isoFaceTol      1e-8;

        surfCellTol     1e-8;

        snapAlphaTol    1e-8;

        nAlphaBounds    3;

        clip            true;

        writeIsoFaces   false;

        writeInterfaceDetails   true;

        phaseChange true;

        prescribedmdot  true;

        massfluxvalue   -100.0;



        nAlphaCorr      2;

        nAlphaSubCycles 1;

        cAlpha          1;  



        MULESCorr       yes;

        nLimiterIter    3;  



        solver          smoothSolver;

        smoother        symGaussSeidel;

        tolerance       1e-8;

        relTol          0;  

        mdotDiffusionConstant   6.25e-08;

    }   



    pcorr

    {   

        solver          PCG;

        preconditioner  DIC;

        tolerance       1e-10;

        relTol          0;  

    }   



    p_rgh

    {   

        solver          PCG;

        preconditioner  DIC;

        tolerance       1e-07;

        relTol          0.05;

    }   



    p_rghFinal

    {   

        $p_rgh;

        tolerance       1e-07;

        relTol          0;  

    }   

    mdot1

    {   

        solver          PCG;

        preconditioner  DIC;

        tolerance       1e-8;

        relTol          0.1;

    }

    mdot1Final

    {

        $mdot1;

    }

    U

    {

        solver          PBiCG;

        preconditioner  DILU;

        tolerance       1e-06;

        relTol          0;

    }

}



PIMPLE

{

    momentumPredictor no;

    nCorrectors     -1;

    nNonOrthogonalCorrectors -1;

    nAlphaCorr      1;

    nAlphaSubCycles 1;

    cAlpha          1;

    pRefCell        0;

    pRefValue       0;

}



SIMPLE

{

    nNonOrthogonalCorrectors 0;

}

I believe that implementation of a Diffusion equation should be straightforward in OpenFOAM. Any help will be much appreciated.

Thanks.