release on 2012-10-23_10-36-40

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -93,7 +93,7 @@ int main(int argc, char *argv[])
 
             fvVectorMatrix UEqn
             (
-                fvm::ddt(U)
+                fvm::ddt(voidfraction,U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
             );
@@ -123,7 +123,7 @@ int main(int argc, char *argv[])
 
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
+                        fvm::laplacian(rUA, p) == fvc::div(phi) + fvc::ddt(voidfraction)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -157,7 +157,8 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         Info << "particleCloud.calcVelocityCorrection() " << endl;
-        particleCloud.calcVelocityCorrection(p,U,phiIB);
+        volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
+        particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
         runTime.write();
 

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
 
-        particleCloud.clockM().start(10,"Flow");
+        particleCloud.clockM().start(16,"Flow");
         // Pressure-velocity PISO corrector
         {
             // Momentum predictor
@@ -94,16 +94,15 @@ int main(int argc, char *argv[])
               - fvm::Sp(Ksl/rho,U)
             );
 
+            if (modelType=="B")
+                UEqn == - fvc::grad(p) + Ksl/rho*Us;
+            else
+                UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
+
             UEqn.relax();
 
             if (momentumPredictor)
-            {
-                //solve UEqn
-                if (modelType=="B")
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
-            }
+                solve(UEqn);
 
             // --- PISO loop
 
@@ -114,15 +113,15 @@ int main(int argc, char *argv[])
             {
                 volScalarField rUA = 1.0/UEqn.A();
                 surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
 
                 U = rUA*UEqn.H();
 
-                phi = fvc::interpolate(U*voidfraction) & mesh.Sf();
-                                      //+ fvc::ddtPhiCorr(rUA, U, phi)
+                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                      + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
                 surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
                 surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
 
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
                 if (modelType=="A")
                     rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
 
@@ -156,7 +155,7 @@ int main(int argc, char *argv[])
 
                 } // end non-orthogonal corrector loop
 
-                #include "continuityErrs.H"
+                #include "continuityErrorPhiPU.H"
 
                 if (modelType=="B")
                     U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;

applications/solvers/cfdemSolverPisoMS/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverPisoMS.C
+
+EXE=$(FOAM_USER_APPBIN)/cfdemSolverPisoMS

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(CFDEM_SRC_DIR)/lnInclude \
+    -I ../cfdemSolverPiso \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(CFDEM_SRC_DIR)/cfdTools \
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN)\
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -l$(CFDEM_LIB_NAME)

applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPisoMS
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, 
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+
+#include "cfdemCloudMS.H"
+#include "implicitCouple.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+
+    #include "initContinuityErrs.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudMS particleCloud(mesh);
+
+    #include "checkModelType.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+
+        // do particle stuff
+        Info << "- evolve()" << endl;
+        particleCloud.evolve(voidfraction,Us,U);
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl.internalField() = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        // debug info
+        //Info << "totaldragforceEuler calculus" << endl;
+        //vector totaldragforceEuler(0,0,0);
+        //forAll(Ksl,cellI)
+        //{
+        //    totaldragforceEuler += Ksl[cellI]*(Us[cellI]-U[cellI])/rho[cellI] * Ksl.mesh().V()[cellI];
+        //}
+        //Pout <<"totaldragforceEuler = "<< mag(totaldragforceEuler) << endl;
+        //Pout << "dv/dt =" << sum(fvc::ddt(voidfraction)) << endl;
+        //-----------
+
+        // Pressure-velocity PISO corrector
+        {
+            // Momentum predictor
+            fvVectorMatrix UEqn
+            (
+                fvm::ddt(voidfraction,U)
+              + fvm::div(phi, U)
+//              + turbulence->divDevReff(U)
+              + particleCloud.divVoidfractionTau(U, voidfraction)
+             == 
+              - fvm::Sp(Ksl/rho,U) 
+            );
+
+            UEqn.relax();
+
+            if (momentumPredictor)
+            {
+                //solve UEqn
+                if (modelType=="B")
+                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+                else
+                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+            }
+
+            // --- PISO loop
+
+            //for (int corr=0; corr<nCorr; corr++)
+            int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
+            for (int corr=0; corr<nCorrSoph; corr++)
+            {
+                volScalarField rUA = 1.0/UEqn.A();
+                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+
+                U = rUA*UEqn.H();
+
+                phi = fvc::interpolate(U*voidfraction) & mesh.Sf();
+                                      //+ fvc::ddtPhiCorr(rUA, U, phi)
+                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
+
+                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                if (modelType=="A")
+                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
+
+                // Non-orthogonal pressure corrector loop
+                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                {
+                    // Pressure corrector
+                    fvScalarMatrix pEqn
+                    (
+                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + fvc::ddt(voidfraction)
+                    );
+                    pEqn.setReference(pRefCell, pRefValue);
+
+                    if
+                    (
+                        corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                    {
+                        pEqn.solve(mesh.solver("pFinal"));
+                    }
+                    else
+                    {
+                        pEqn.solve();
+                    }
+
+                    if (nonOrth == nNonOrthCorr)
+                    {
+                        phiGes -= pEqn.flux();
+                    }
+
+                } // end non-orthogonal corrector loop
+
+                #include "continuityErrorPhiPU.H"
+
+                if (modelType=="B")
+                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                else
+                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                
+                U.correctBoundaryConditions();
+
+            } // end piso loop
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -126,14 +126,14 @@ int main(int argc, char *argv[])
             {
                 volScalarField rUA = 1.0/UEqn.A();
                 surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
 
                 U = rUA*UEqn.H();
 
-                phi = fvc::interpolate(U*voidfraction) & mesh.Sf();
-                                      //+ fvc::ddtPhiCorr(rUA, U, phi)
+                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                      + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
                 surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
                 surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
                 if (modelType=="A")
                     rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
 
@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
 
                 } // end non-orthogonal corrector loop
 
-                #include "continuityErrs.H"
+                #include "continuityErrorPhiPU.H"
 
                 if (modelType=="B")
                     U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;