update cfdemSolverPimple for periodic box simulations

applications/solvers/cfdemSolverPimple/UEqn.H

@@ -4,7 +4,7 @@ tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
   + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + particleCloud.divVoidfractionTau(U, voidfraction)  // in case of doing "periodic box" simulations viscous term can be commented and the solver can be compiled 							       // again. However, the effect of this term on the results is negligible (about 1-2%). 
   - fOther/rho
   ==
     fvOptions(U)

applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C

@@ -58,6 +58,12 @@ int main(int argc, char *argv[])
     cfdemCloud particleCloud(mesh);
     #include "checkModelType.H"
 
+    // switch for periodic box simulations
+    Switch periodicBoxSwitch
+    (
+        pimple.dict().lookupOrDefault<Switch>("periodicBoxSwitch", false)
+    );    
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
@@ -117,6 +123,11 @@ int main(int argc, char *argv[])
             Info << "skipping flow solution." << endl;
         }
 
+        if (periodicBoxSwitch) 
+        {         
+           #include "periodicBoxProperties.H"
+        }
+         
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverPimple/createFields.H

@@ -139,6 +139,30 @@
         linearInterpolate(U) & mesh.Sf()
     );
 
+//Periodic box
+    volScalarField unity
+    (
+     	IOobject
+     	(
+	   "unity",
+	   runTime.timeName(),
+	   mesh,
+	   IOobject::NO_READ,
+	   IOobject::NO_WRITE
+     	),
+     	mesh,
+     	dimensionedScalar("unity", dimless, 1.0)
+    );
+
+// gravity vector for periodic boundary conditions
+    dimensionedVector gN
+    (
+       "gN",
+       dimensionSet(0,1,-2,0,0,0,0),
+       vector(0,0,-9.81)
+    );
+
+
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;

applications/solvers/cfdemSolverPimple/periodicBoxProperties.H

@@ -0,0 +1,65 @@
+// gravity direction should always be -z!
+dimensionedScalar volume = fvc::domainIntegrate(unity);
+
+Info<< "particle_ENSTROPHY: "
+    << (
+            fvc::domainIntegrate( 0.5*magSqr(fvc::curl(Us)))
+            /volume
+       ).value()
+    << endl;
+
+Info<< "air_ENSTROPHY: "
+    << (
+           fvc::domainIntegrate( 0.5*magSqr(fvc::curl(U)))
+           /volume
+       ).value()
+    << endl;
+
+Info<< "slip_velocity: "
+    << - ((
+		fvc::domainIntegrate(voidfraction*(U&gN)).value()
+	       /fvc::domainIntegrate(voidfraction*mag(gN)).value()
+	     )
+	   - (
+		fvc::domainIntegrate((1.0-voidfraction)*(Us&gN)).value()
+	       /fvc::domainIntegrate((1.0-voidfraction)*mag(gN)).value()
+         ))
+   << endl;
+   
+dimensionedVector alpha1Us = fvc::domainIntegrate((1.0-voidfraction)*(Us))/fvc::domainIntegrate((1.0-voidfraction));
+dimensionedVector alpha2U = fvc::domainIntegrate(voidfraction*(U))/fvc::domainIntegrate(voidfraction);
+dimensionedScalar alpha1M  = fvc::domainIntegrate((1.0-voidfraction))/volume;
+dimensionedScalar alpha2M  = scalar(1.0) - alpha1M;
+
+Info<< "TKE gas: "
+    << 0.5
+      *(
+          fvc::domainIntegrate(voidfraction*(U&U)).value()
+         /fvc::domainIntegrate(voidfraction).value()
+       )
+     - 0.5
+      *(
+          alpha2U.value()
+         &alpha2U.value()
+       )
+    << endl;
+
+Info<< "TKE solid: "
+    << 0.5
+      *(
+          fvc::domainIntegrate((1.0-voidfraction)*(Us&Us)).value()
+         /fvc::domainIntegrate(1.0-voidfraction).value()
+       )
+     - 0.5
+      *(
+          alpha1Us.value()
+         &alpha1Us.value()
+       )
+    << endl;
+
+Info<< "PhiP2: "
+    << fvc::domainIntegrate((1.0-voidfraction)*(1.0-voidfraction)).value()
+       /fvc::domainIntegrate(unity).value()
+       - alpha1M.value()*alpha1M.value()
+    << endl;
+