initial commit for mass transfer

applications/solvers/cfdemSolverMultiphaseScalar/CEqn.H

@@ -0,0 +1,26 @@
+        // get mixture properties
+        D  = mixture.D();
+	Cs  = mixture.Cs();
+        //keff=particleCloud.thermCondM().thermCond(); 
+
+        // get scalar source from DEM        
+        //particleCloud.energyContributions(Qsource);
+        //particleCloud.energyCoefficients(QCoeff);
+
+        fvScalarMatrix CEqn
+        (
+	 fvm::ddt(voidfraction,C)
+		 + fvm::div(phi,C)
+	 - fvc::div(D*voidfraction*fvc::grad(C))
+         ==
+	 Csource + fvm::Sp(CCoeff,C)
+        );
+
+
+        CEqn.relax();
+        fvOptions.constrain(CEqn);
+        CEqn.solve();
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");

applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C

@@ -111,6 +111,7 @@ int main(int argc, char *argv[])
             rhoEps = rho * voidfraction;
 
 	    #include "EEqn.H"
+	    #include "CEqn.H"
             
              // --- Pressure-velocity PIMPLE corrector loop
             while (pimple.loop())

applications/solvers/cfdemSolverMultiphaseScalar/createFields.H

@@ -190,6 +190,75 @@ volScalarField keff
     dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
 );
 
+Info<< "Reading/creating concentration fields\n" << endl;
+volScalarField C
+(
+    IOobject
+    (
+        "C",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField Csource
+(
+    IOobject
+    (
+        "Csource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
+);
+
+volScalarField CCoeff
+(
+    IOobject
+    (
+        "Csource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(0,0,-1,0,0,0,0), 0.0)
+);
+
+
+volScalarField D
+(
+    IOobject
+    (
+        "D",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.D()
+);
+
+volScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.Cs()
+);
+
 
 //========================
 

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C

@@ -416,6 +416,45 @@ Foam::multiphaseMixture::kf() const
     return tkf;
 }
 
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::D() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // 1/D
+    tmp<volScalarField> tDInv = iter()/iter().D();
+    volScalarField& DInv = tDInv.ref();
+
+    // D
+    tmp<volScalarField> tD = iter()*iter().D();
+    volScalarField& D = tD.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+      DInv += iter()/iter().D();
+    }
+
+    D = 1/DInv;
+    return tD;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::Cs() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // Cs
+    tmp<volScalarField> tCs = iter()*iter().Cs();
+    volScalarField& Cs = tCs.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        Cs += iter()*iter().Cs();
+    }
+
+    return tCs;
+}
+
 void Foam::multiphaseMixture::solve()
 {
     correct();

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H

@@ -259,6 +259,12 @@ public:
         //- Return the thermal conductivity
         tmp<volScalarField> kf() const;
 
+        //- Return the diffusion coefficient
+        tmp<volScalarField> D() const;
+
+        //- Return the solubility
+        tmp<volScalarField> Cs() const;
+
         tmp<surfaceScalarField> surfaceTensionForce() const
         {
             return surfaceTensionForce_;

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C

@@ -61,7 +61,9 @@ Foam::phase::phase
     ),
     rho_("rho", dimDensity, phaseDict_),
     Cp_("Cp", (dimSpecificHeatCapacity), phaseDict_),
-    kf_("kf", (dimPower/dimLength/dimTemperature), phaseDict_)
+    kf_("kf", (dimPower/dimLength/dimTemperature), phaseDict_),
+    D_("D", dimViscosity, phaseDict_),
+    Cs_("Cs", dimDensity, phaseDict_)
 {}
 
 
@@ -86,6 +88,8 @@ bool Foam::phase::read(const dictionary& phaseDict)
 
     phaseDict_.lookup("Cp") >> Cp_;
     phaseDict_.lookup("kf") >> kf_;
+    phaseDict_.lookup("D") >> D_;
+    phaseDict_.lookup("Cs") >> Cs_;
 
     if (nuModel_->read(phaseDict_))
     {

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H

@@ -62,6 +62,8 @@ class phase
         dimensionedScalar rho_;
         dimensionedScalar Cp_;
         dimensionedScalar kf_;
+        dimensionedScalar D_;
+        dimensionedScalar Cs_;
 
 public:
 
@@ -154,6 +156,18 @@ public:
             return kf_;
         }
 
+         //- Return const-access to phase1 diffusion coefficient
+        const dimensionedScalar& D() const
+        {
+            return D_;
+        }
+
+        //- Return const-access to phase1 solubility
+        const dimensionedScalar& Cs() const
+        {
+            return Cs_;
+        } 
+
         //- Correct the phase properties
         void correct();