From 4a24c35c679cccd007e8120084b4dcf84075fe07 Mon Sep 17 00:00:00 2001
From: ekinaci <efekinaci@yahoo.com>
Date: Fri, 3 Feb 2017 12:20:08 +0100
Subject: [PATCH] modifed solver

---
 .../solvers/cfdemSolverRhoPimpleChem/EEqn.H   | 58 +++++++++----------
 .../cfdemSolverRhoPimpleChem.C                | 15 +++--
 2 files changed, 35 insertions(+), 38 deletions(-)

diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
index 606d9fe2..6dbf1d50 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H
@@ -1,45 +1,46 @@
 // contributions to internal energy equation can be found in
 // Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
 {
-    // dim he = J / kg
-    volScalarField& he = thermo.he();
-    particleCloud.energyContributions(Qsource);
-    particleCloud.energyCoefficients(QCoeff);
+// dim he = J / kg
+volScalarField& he = thermo.he();
+particleCloud.energyContributions(Qsource);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();   
-    thCond=particleCloud.thermCondM().thermCond(); 
+thDiff=particleCloud.thermCondM().thermDiff();
+// BC for thDiff are calculated from its mathematical expression
+// to enforce user-defined values, change thDiff to MUST_READ and do
+// thDiff.correctBoundaryConditions();
 
-    addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
-      + (
-            he.name() == "e"
-          ? fvc::div
+fvScalarMatrix EEqn
+(
+    fvm::ddt(rhoeps, he) + fvm::div(phi, he)
+  + fvc::ddt(rhoeps, K) + fvc::div(phi, K)
+    + (
+        he.name() == "e"
+        ? fvc::div
             (
                 fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
                 p,
-                "div(phiv,p)"
-            )
+              "div(phiv,p)"
+           )
           : -dpdt
-        );
-
-    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
-
-
-    fvScalarMatrix EEqn
-    (
-        fvm::ddt(rhoeps, he) + fvm::div(phi, he)
-      + addSource
+       )
      // net heat transfer from particles to fluid
       - Qsource
-      - fvm::Sp(QCoeff/Cpv, he)
-     // thermal conduction of the fluid with effective conductivity
-     // - fvm::laplacian(rhoeps*thDiff,he)
-      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+    // thermal conduction of the fluid with effective conductivity
+      - fvm::laplacian(rhoeps*thDiff,he)
+
       // + particle-fluid energy transfer due to work
-      // + fluid energy dissipation due to shearing
+     // + fluid energy dissipation due to shearing
      ==
         fvOptions(rho, he)
     );
-    
+
+
+    // work in progress:
+    // if (particleCloud.implicitEnergyModel)
+    // volScalarField& tpart (lookup ...)
+    // particleCloud.energyCoefficients(Qcoeff);
+    // EEqn += ...
 
     EEqn.relax();
 
@@ -53,7 +54,4 @@
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
-    particleCloud.clockM().start(31,"postFlow");
-    particleCloud.postFlow();
-    particleCloud.clockM().stop("postFlow");
 }
diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C b/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
index 602d11ae..42e6d31f 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
@@ -53,10 +53,13 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
 
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createFieldRefs.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     
@@ -97,19 +100,15 @@ int main(int argc, char *argv[])
 
        //Force Checks
        vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
+       vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
        reduce(fImpTotal, sumOp<vector>());
        Info << "TotalForceExp: " << fTotal << endl;
        Info << "TotalForceImp: " << fImpTotal << endl;
 
-       Info << "number of particles = " << particleCloud.numberOfParticles() << nl << endl;
+       #include "solverDebugInfo.H"
+       particleCloud.clockM().stop("Coupling");
 
-        #include "solverDebugInfo.H"
-        particleCloud.clockM().stop("Coupling");
-
-
-
-        particleCloud.clockM().start(26,"Flow");
+       particleCloud.clockM().start(26,"Flow");
 	
         if (pimple.nCorrPIMPLE() <= 1)
         {