# Problem with Diffusion Equation Implementation

 Register Blogs Members List Search Today's Posts Mark Forums Read

February 22, 2019, 11:45
Problem with Diffusion Equation Implementation
#1
Member

Raunak Bardia
Join Date: Jan 2015
Posts: 32
Rep Power: 11
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

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

mdot1 = mdot0; // mdot0 stores the initial condition
(
fvm::laplacian(mdotDiffusionConstant, mdot1)
);
The mdot fields were created as follows.

createFields.H
Code:
```// Calculating Source Term at Interfacial Cells
volScalarField mdot0
(
IOobject
(
"mdot0",
runTime.timeName(),
mesh,
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::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.*"
{
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.
Attached Images
 InitialCondition.jpg (25.9 KB, 15 views) Diffused.jpg (32.4 KB, 14 views)

 February 25, 2019, 03:04 #2 Senior Member   Robert Join Date: May 2015 Location: Bremen, GER Posts: 292 Rep Power: 12 Did you allready take a look at the scalarTransportFoam solver? It solves an advection-diffusion equation. Maybe you can take parts of that eqn. and implement it into your solver. __________________ If you liked my answer to your question, please consider leaving a "Like" in return

 Tags diffusion equation, interface advection, openfoam 1706, solver