spatially-varied time-constant diffusivity in interMixingFoam
 

I'm modeling diffusion/perfusion of two miscible chemical solutions through human tissue with interMixingFoam in OpenFOAM 7. The physical components in the domain (eg. cells, membranes, intercellular matrix) all are fixed in position during the simulation, but each component has its own diffusion coefficients in the real world. So, instead of using the constant diffusion coefficient D23, it would be more realistic to set up a volScalarField for DiffCoef (just picking a name) in the time 0 folder like is done for alpha2, use topoSetDict to define cellZones for the different physical components, and set their unique values of DiffCoef in setFieldsDict. But I suspect this would require modification of the solver and I'm better at Fortran77 and Scilab than C++ so I could use some hints.


I see that in alphaEqn.H for interMixingFoam, the constant D23 is used to construct the volScalarFields Dc23 and Dc32. So I assume these lines would need to be altered to handle D23 as a volScalarField (ie. my DiffCoef). It's not clear to me how D23 got to this namespace, but volScalarField DiffCoef would need to find its way here. Is there another solver that has a volScalarField diffusivity that I could emulate? I'm fairly good at monkey see, monkey do.


Thanks for your help, -Ron.

update to spatially-varied time-constant diffusivity in interMixingFoam
 

This turned out to be easier than I imagined. In my test case, I set the diffusion coefficient to zero in a region and diffusion was prevented there as expected. Here is what I did, for future reference.

To create the new version of interMixingFoam, I used:

mkdir -p $WM_PROJECT_USER_DIR/applications/solvers
cp -r interMixingFoam $WM_PROJECT_USER_DIR/applications/solvers/interMixingDFieldFoam
cd $WM_PROJECT_USER_DIR/applications/solvers/interMixingDFieldFoam
mv interMixingFoam.C interMixingDFieldFoam.C
rm interMixingFoam.dep

Then I edited Make/files to update the last lines:

interMixingDFieldFoam.C

EXE = $(FOAM_USER_APPBIN)/interMixingDFieldFoam

I removed the old binaries:

rm -rf Make/linux64GccDPint32Opt

I tested the compilation process with:
wmake

Then confirmed the new solver binary was in:
ls $FOAM_USER_APPBIN

I modified createFields.H to replace the scalar D23 with the volScalarField Diff. I doesn't need a value specified after the second "mesh", because it will be reading the field of values set in setFieldsDict. The dimensions are set in transportProperties. It is set to NO-WRITE because the initialized values are never changed.

//dimensionedScalar D23("D23", dimViscosity, mixture);
volScalarField Diff
(
IOobject
(
"Diff",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);

In alphaEqn.H, I replaced D23 with Diff. No other changes needed.

// Create the diffusion coefficients for alpha2<->alpha3
// volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2));
// volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3));
volScalarField Dc23(Diff*max(alpha3, scalar(0))*pos0(alpha2));
volScalarField Dc32(Diff*max(alpha2, scalar(0))*pos0(alpha3));

I used wclean and then wmake to compile the new solver, then focused on trying it out.

In the time 0.orig folder in my case folder (run/NTB30), I added a file named Diff with all the boundary conditions set to zeroGradient. (I don't know if these are even needed):

FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object Diff;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}

In the transportProperties file I replaced D23 with Diff and did not specify a value since they are initialized in setFieldsDict:

// Diffusivity between miscible phases,
//D23 80.0e-5 ;
Diff Diff [0 2 -1 0 0 0 0];

In setFieldsDict, I added Diff to the defaultFieldValues and any region (which are defined in topoSetDict) in which the value is different:

volScalarFieldValue Diff 102.0e-5

regions
(
zoneToCell
{
name "meshLayer";
fieldValues
(
volScalarFieldValue alpha.two 0
volScalarFieldValue alpha.three 1
volScalarFieldValue Diff 0.008
);
}