Merge branch 'feature/cfdemSolverRhoPimple' of https://github.com/ParticulateFlow/CFDEMcoupling into feature/cfdemSolverRhoPimple

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -137,17 +137,15 @@ int main(int argc, char *argv[])
 	        #include "molConc.H"
                 #include "pEqn.H"
             }
-            m=gSum(rhoeps*1.0*rhoeps.mesh().V());
+           
-            if(counter==0) m0=m;
-	    counter++;
-            Info << "\ncurrent gas mass = " << m << "\n" << endl;
-	    Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl;
             if (pimple.turbCorr())
             {
                 turbulence->correct();
             }
         }
 
+        #include "monitorMass.H"
+        
         particleCloud.clockM().start(31,"postFlow");
         particleCloud.postFlow();
         particleCloud.clockM().stop("postFlow");

applications/solvers/cfdemSolverRhoPimpleChem/monitorMass.H

@@ -0,0 +1,7 @@
+{
+    m=gSum(rhoeps*1.0*rhoeps.mesh().V());
+    if(counter==0) m0=m;
+    counter++;
+    Info << "\ncurrent gas mass = " << m << "\n" << endl;
+    Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl;
+}

applications/solvers/cfdemSolverRhoPimpleChem/pEqn.H

@@ -3,12 +3,6 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-rhoeps = rho * voidfraction;
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution - done in 2 parts. Part 1:
-thermo.rho() -= psi*p;
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
 if (modelType=="A")
@@ -44,22 +38,20 @@ else
 
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
-
+    
     volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
         (
-         //   fvm::ddt(psi*voidfraction, p)
+          fvm::ddt(voidfraction, psi, p)
-          fvc::ddt(rhoeps) + psi*correction(fvm::ddt(voidfraction,p))
           + fvc::div(phi)
           - fvm::laplacian(rhorAUf, p)
           ==
-     //      particleCloud.chemistryM(0).Sm()
             fvm::Sp(SmbyP, p)
           + fvOptions(psi, p, rho.name())
-
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -76,21 +68,17 @@ else
 
 // Explicitly relax pressure for momentum corrector
 p.relax();
-// Second part of thermodynamic density update
-thermo.rho() += psi*p;
 
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
-
-rhoeps = rho * voidfraction;
-
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
+rhoeps = rho * voidfraction;
+
 if (modelType=="A")
 {
     U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);

applications/solvers/cfdemSolverRhoPimpleChem/pEqn_alternative.H

@@ -0,0 +1,109 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+rhoeps = rho * voidfraction;
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution - done in 2 parts. Part 1:
+thermo.rho() -= psi*p;
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
+if (modelType=="A")
+{
+    rhorAUf *= fvc::interpolate(voidfraction);
+}
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+
+if (pimple.nCorrPISO() <= 1)
+{
+    tUEqn.clear();
+}
+
+if (pimple.transonic())
+{
+//  transonic version not implemented yet
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rhoeps*HbyA)
+    //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+        )
+    );
+    
+    // flux without pressure gradient contribution
+    phi = phiHbyA + phiUs;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rhoeps, U, phi, rhorAUf);
+
+    volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+         //   fvm::ddt(psi*voidfraction, p)
+          fvc::ddt(rhoeps) + psi*correction(fvm::ddt(voidfraction,p))
+          + fvc::div(phi)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+     //      particleCloud.chemistryM(0).Sm()
+            fvm::Sp(SmbyP, p)
+          + fvOptions(psi, p, rho.name())
+
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrsPU.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+// Second part of thermodynamic density update
+thermo.rho() += psi*p;
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+rhoeps = rho * voidfraction;
+
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+if (modelType=="A")
+{
+    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+}
+else
+{
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+}
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(voidfraction,p);
+}

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C

@@ -334,7 +334,7 @@ void species::execute()
         }
         massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * timestep);
 	massSourceTot_ += massSourceCurr_;
-	Info << "total conversion of mass:\tcurrent source = " << massSourceCurr_ << ", total source = " << massSourceTot_ << "\n" << endl;
+	Info << "total conversion of mass:\n\tcurrent source = " << massSourceCurr_ << "\n\ttotal source = " << massSourceTot_ << "\n" << endl;
         Info << "get data done" << endl;
 }
 

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/0/CO

@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 0 0 0 0 0];
 
-internalField   uniform 0.3;
+internalField   uniform 0.0;
 
 boundaryField
 {

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/0/O2

@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 0 0 0 0 0];
 
-internalField   uniform 0.7;
+internalField   uniform 1.0;
 
 boundaryField
 {

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/0/T

@@ -23,27 +23,32 @@ boundaryField
 {
     top
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1293.15;
     }
 
     bottom
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1293.15;
     }
 
     side-walls
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1293.15;
     }
 
     inlet
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1293.15;
     }
 
     outlet
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 1293.15;
     }
 }
 

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/constant/couplingProperties

@@ -31,7 +31,7 @@ FoamFile
 
 modelType "A"; // A or B
 
-couplingInterval 100;//1000;
+couplingInterval 10;//1000;
 
 voidFractionModel divided;//centre;//
 

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/system/controlDict

@@ -25,11 +25,11 @@ stopAt          endTime;
 
 endTime         5.0;
 
-deltaT          0.0001;
+deltaT          0.00001;
 
 writeControl    timeStep;
 
-writeInterval   50;
+writeInterval   1000;
 
 purgeWrite      0;
 
@@ -104,6 +104,7 @@ functions
         (
             rhoeps
             rho
+            molarConc
         );
     }
 

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/system/fvOptions

@@ -22,7 +22,7 @@ source1
     {
         active          yes;
         selectionMode   all;
-        Tmin     1293.1;
+        Tmin     1293.15;
-        Tmax     1293.2;
+        Tmax     2000;
     }
 }

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/system/fvSolution

@@ -58,14 +58,14 @@ solvers
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
-        tolerance       1e-05;
+        tolerance       1e-06;
         relTol          0.1;
     }
 
     "(U|h|R|k|epsilon)Final"
     {
         $U;
-        tolerance       1e-05;
+        tolerance       1e-07;
         relTol          0;
     }
 
@@ -77,7 +77,7 @@ solvers
     "(Yi|CO2|O2)Final"
     {
         $Yi;
-        tolerance       1e-06;
+        tolerance       1e-08;
         relTol          0;
     }
 }
@@ -98,13 +98,13 @@ relaxationFactors
 {
     fields
     {
-        ".*"            1;
+          p               0.3;
     }
     equations
     {
         ".*"            1;
     }
 }
-
+*/
 
 // ************************************************************************* //

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/CFD/system/probesDict

@@ -21,7 +21,7 @@ fields
     rho
     p
     T
-    molC
+    molarConc
     O2
     O
     CO2

tutorials/cfdemSolverRhoPimpleChem/Shrinking_particle_model_multispecies/DEM/in.liggghts_run

@@ -45,21 +45,21 @@ fix     zwalls2 all wall/gran model hertz tangential history primitive type 1 zp
 
 ###############################################
 # cfd coupling
-fix     cfd all couple/cfd couple_every 100 mpi
+fix     cfd all couple/cfd couple_every 10 mpi
 fix     cfd2 all couple/cfd/force
 
 # this should invoke chemistry
 fix     cfd3 all couple/cfd/chemistry n_species 5 species_names O2 CO2 N2 CO O
 
 # this should shrink the particle
-#fix     cfd4 all chem/shrink speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 k 2.5e1 rmin 0.005 nevery 100 screen yes
+fix     cfd4 all chem/shrink speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 k 0.1 rmin 0.005 nevery 10 screen yes
-#fix     cfd5 all chem/shrink speciesA CO2 molMassA 44.01 speciesC CO molMassC 28.01 nuC 2 molMassB 12.01 k 2.5e1 rmin 0.005 nevery 100
+fix     cfd5 all chem/shrink speciesA CO2 molMassA 44.01 speciesC CO molMassC 28.01 nuC 2 molMassB 12.01 k 0.02 rmin 0.005 nevery 10 screen yes
 
 # values for k0 and T0 according to Shen et al., Fuel (2011) with a coke density of 800 kg / m^3
 # C + O2 -> CO2		k0 = 4e3 m/s		T0 = 10855 K 
 # C + CO2 -> 2CO	k0 = 6e6 m/s		T0 = 29018 K 
-fix     cfd4 all chem/shrink/Arrhenius speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 rmin 0.005 nevery 100 screen yes k 4.0e3 T 10855
+#fix     cfd4 all chem/shrink/Arrhenius speciesA O2 molMassA 31.99 speciesC CO2 molMassC 44.01 molMassB 12.01 rmin 0.005 nevery 10 screen yes k 4.0e3 T 10855
-fix	cfd5 all chem/shrink/Arrhenius speciesA CO2 molMassA 44.01 speciesC CO molMassC 28.01 nuC 2 molMassB 12.01 rmin 0.005 nevery 100 screen yes k 6e6 T 29018
+#fix	cfd5 all chem/shrink/Arrhenius speciesA CO2 molMassA 44.01 speciesC CO molMassC 28.01 nuC 2 molMassB 12.01 rmin 0.005 nevery 10 screen yes k 6e6 T 29018
 
 
 # apply nve integration to all particles that are inserted as single particles