Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,6 +90,7 @@ int main(int argc, char *argv[])
             if (pimple.turbCorr())
             {
                 #include "kEpsilon.H"
+                nuEffa = sqr(Ct)*nutb + nua;
             }
         }
 

applications/solvers/multiphase/bubbleFoam/kEpsilon.H

@@ -14,16 +14,17 @@ if (turbulence)
     // Dissipation equation
     fvScalarMatrix epsEqn
     (
-        fvm::ddt(beta, epsilon)
+        fvm::ddt(epsilon)
       + fvm::div(phib, epsilon)
+      - fvm::Sp(fvc::div(phib), epsilon)
       - fvm::laplacian
         (
             alphaEps*nuEffb, epsilon,
             "laplacian(DepsilonEff,epsilon)"
         )
       ==
-         C1*beta*G*epsilon/k
+         C1*G*epsilon/k
-       - fvm::Sp(C2*beta*epsilon/k, epsilon)
+       - fvm::Sp(C2*epsilon/k, epsilon)
     );
 
     #include "wallDissipation.H"
@@ -37,16 +38,17 @@ if (turbulence)
     // Turbulent kinetic energy equation
     fvScalarMatrix kEqn
     (
-        fvm::ddt(beta, k)
+        fvm::ddt(k)
       + fvm::div(phib, k)
+      - fvm::Sp(fvc::div(phib), k)
       - fvm::laplacian
         (
             alphak*nuEffb, k,
             "laplacian(DkEff,k)"
         )
       ==
-        beta*G
+        G
-      - fvm::Sp(beta*epsilon/k, k)
+      - fvm::Sp(epsilon/k, k)
     );
     kEqn.relax();
     kEqn.solve();
@@ -59,5 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEffa = sqr(Ct)*nutb + nua;
 nuEffb = nutb + nub;

applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean

@@ -6,5 +6,6 @@ wclean libso phaseModel
 wclean libso interfacialModels
 wclean libso kineticTheoryModels
 wclean
+wclean MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake

@@ -6,5 +6,6 @@ wmake libso phaseModel
 wmake libso interfacialModels
 wmake libso kineticTheoryModels
 wmake
+wmake MRFtwoPhaseEulerFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    twoPhaseEulerFoam
+
+Description
+    Solver for a system of 2 incompressible fluid phases with one phase
+    dispersed, e.g. gas bubbles in a liquid or solid particles in a gas.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "nearWallDist.H"
+#include "wallFvPatch.H"
+#include "Switch.H"
+
+#include "IFstream.H"
+#include "OFstream.H"
+
+#include "dragModel.H"
+#include "phaseModel.H"
+#include "kineticTheoryModel.H"
+
+#include "pimpleControl.H"
+#include "MRFZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "readPPProperties.H"
+    #include "initContinuityErrs.H"
+    #include "createMRFZones.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    pimpleControl pimple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTwoPhaseEulerFoamControls.H"
+        #include "CourantNos.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaEqn.H"
+            #include "liftDragCoeffs.H"
+            #include "UEqns.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+
+                if (correctAlpha && !pimple.finalIter())
+                {
+                    #include "alphaEqn.H"
+                }
+            }
+
+            #include "DDtU.H"
+
+            if (pimple.turbCorr())
+            {
+                #include "kEpsilon.H"
+            }
+        }
+
+        #include "write.H"
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files

@@ -0,0 +1,3 @@
+MRFTwoPhaseEulerFoam.C
+
+EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -I.. \
+    -I../../bubbleFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -IturbulenceModel \
+    -I../kineticTheoryModels/lnInclude \
+    -I../interfacialModels/lnInclude \
+    -I../phaseModel/lnInclude \
+
+EXE_LIBS = \
+    -lEulerianInterfacialModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lincompressibleTransportModels \
+    -lphaseModel \
+    -lkineticTheoryModel

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H

@@ -0,0 +1,99 @@
+fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
+fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
+
+{
+    {
+        volTensorField gradUaT(T(fvc::grad(Ua)));
+
+        if (kineticTheory.on())
+        {
+            kineticTheory.solve(gradUaT);
+            nuEffa = kineticTheory.mua()/rhoa;
+        }
+        else // If not using kinetic theory is using Ct model
+        {
+            nuEffa = sqr(Ct)*nutb + nua;
+        }
+
+        volTensorField Rca
+        (
+            "Rca",
+            (((2.0/3.0)*I)*nuEffa)*tr(gradUaT) - nuEffa*gradUaT
+        );
+
+        if (kineticTheory.on())
+        {
+            Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
+        }
+
+        surfaceScalarField phiRa
+        (
+            -fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
+            /fvc::interpolate(alpha + scalar(0.001))
+        );
+
+        UaEqn =
+        (
+            (scalar(1) + Cvm*rhob*beta/rhoa)*
+            (
+                fvm::ddt(Ua)
+              + fvm::div(phia, Ua, "div(phia,Ua)")
+              - fvm::Sp(fvc::div(phia), Ua)
+            )
+
+          - fvm::laplacian(nuEffa, Ua)
+          + fvc::div(Rca)
+
+          + fvm::div(phiRa, Ua, "div(phia,Ua)")
+          - fvm::Sp(fvc::div(phiRa), Ua)
+          + (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
+         ==
+        //  g                          // Buoyancy term transfered to p-equation
+          - fvm::Sp(beta/rhoa*K, Ua)
+        //+ beta/rhoa*K*Ub             // Explicit drag transfered to p-equation
+          - beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
+        );
+        mrfZones.addCoriolis(UaEqn);
+        UaEqn.relax();
+    }
+
+    {
+        volTensorField gradUbT(T(fvc::grad(Ub)));
+        volTensorField Rcb
+        (
+            "Rcb",
+            (((2.0/3.0)*I)*nuEffb)*tr(gradUbT) - nuEffb*gradUbT
+        );
+
+        surfaceScalarField phiRb
+        (
+            -fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
+            /fvc::interpolate(beta + scalar(0.001))
+        );
+
+        UbEqn =
+        (
+            (scalar(1) + Cvm*rhob*alpha/rhob)*
+            (
+                fvm::ddt(Ub)
+              + fvm::div(phib, Ub, "div(phib,Ub)")
+              - fvm::Sp(fvc::div(phib), Ub)
+            )
+
+          - fvm::laplacian(nuEffb, Ub)
+          + fvc::div(Rcb)
+
+          + fvm::div(phiRb, Ub, "div(phib,Ub)")
+          - fvm::Sp(fvc::div(phiRb), Ub)
+
+          + (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
+         ==
+        //  g                          // Buoyancy term transfered to p-equation
+          - fvm::Sp(alpha/rhob*K, Ub)
+        //+ alpha/rhob*K*Ua            // Explicit drag transfered to p-equation
+          + alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
+        );
+        mrfZones.addCoriolis(UbEqn);
+        UbEqn.relax();
+    }
+}

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H

@@ -0,0 +1,3 @@
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(Ua);
+    mrfZones.correctBoundaryVelocity(Ub);

applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H

@@ -0,0 +1,121 @@
+{
+    surfaceScalarField alphaf(fvc::interpolate(alpha));
+    surfaceScalarField betaf(scalar(1) - alphaf);
+
+    volScalarField rUaA(1.0/UaEqn.A());
+    volScalarField rUbA(1.0/UbEqn.A());
+
+    rUaAf = fvc::interpolate(rUaA);
+    surfaceScalarField rUbAf(fvc::interpolate(rUbA));
+
+    volVectorField HbyAa("HbyAa", Ua);
+    HbyAa = rUaA*UaEqn.H();
+
+    volVectorField HbyAb("HbyAb", Ub);
+    HbyAb = rUbA*UbEqn.H();
+
+    mrfZones.absoluteFlux(phia.oldTime());
+    mrfZones.absoluteFlux(phia);
+
+    mrfZones.absoluteFlux(phib.oldTime());
+    mrfZones.absoluteFlux(phib);
+
+    surfaceScalarField phiDraga
+    (
+        fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
+    );
+
+    if (g0.value() > 0.0)
+    {
+        phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
+    }
+
+    if (kineticTheory.on())
+    {
+        phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
+    }
+
+
+    surfaceScalarField phiDragb
+    (
+        fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
+    );
+
+    // Fix for gravity on outlet boundary.
+    forAll(p.boundaryField(), patchi)
+    {
+        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
+        {
+            phiDraga.boundaryField()[patchi] = 0.0;
+            phiDragb.boundaryField()[patchi] = 0.0;
+        }
+    }
+
+    surfaceScalarField phiHbyAa
+    (
+        "phiHbyAa",
+        (fvc::interpolate(HbyAa) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUaA, Ua, phia)
+      + phiDraga
+    );
+    mrfZones.relativeFlux(phiHbyAa);
+
+    surfaceScalarField phiHbyAb
+    (
+        "phiHbyAb",
+        (fvc::interpolate(HbyAb) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUbA, Ub, phib)
+      + phiDragb
+    );
+    mrfZones.relativeFlux(phiHbyAb);
+
+    surfaceScalarField phiHbyA("phiHbyA", alphaf*phiHbyAa + betaf*phiHbyAb);
+
+    surfaceScalarField Dp
+    (
+        "Dp",
+        alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            surfaceScalarField SfGradp(pEqn.flux()/Dp);
+
+            phia.boundaryField() ==
+                (fvc::interpolate(Ua) & mesh.Sf())().boundaryField();
+            mrfZones.relativeFlux(phia);
+            phia = phiHbyAa - rUaAf*SfGradp/rhoa;
+
+            phib.boundaryField() ==
+                (fvc::interpolate(Ub) & mesh.Sf())().boundaryField();
+            mrfZones.relativeFlux(phib);
+            phib = phiHbyAb - rUbAf*SfGradp/rhob;
+
+            phi = alphaf*phia + betaf*phib;
+
+            p.relax();
+            SfGradp = pEqn.flux()/Dp;
+
+            Ua = HbyAa + fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
+            Ua.correctBoundaryConditions();
+
+            Ub = HbyAb + fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
+            Ub.correctBoundaryConditions();
+
+            U = alpha*Ua + beta*Ub;
+        }
+    }
+}
+
+#include "continuityErrs.H"

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -18,7 +18,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
         volTensorField Rca
         (
             "Rca",
-            ((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT
+            (((2.0/3.0)*I)*nuEffa)*tr(gradUaT) - nuEffa*gradUaT
         );
 
         if (kineticTheory.on())
@@ -62,7 +62,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
         volTensorField Rcb
         (
             "Rcb",
-            ((2.0/3.0)*I)*(k + nuEffb*tr(gradUbT)) - nuEffb*gradUbT
+            (((2.0/3.0)*I)*nuEffb)*tr(gradUbT) - nuEffb*gradUbT
         );
 
         surfaceScalarField phiRb

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -112,7 +112,9 @@
         (
             "phi",
             runTime.timeName(),
-            mesh
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
         ),
         fvc::interpolate(alpha)*phia + fvc::interpolate(beta)*phib
     );
@@ -194,6 +196,17 @@
         }
     }
 
+    dimensionedScalar residualSlip
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "residualSlip",
+            interfacialProperties,
+            0,
+            dimVelocity
+        )
+    );
+
     Info << "dragPhase is " << dragPhase << endl;
     kineticTheoryModel kineticTheory
     (

applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H

@@ -1,62 +0,0 @@
-if (turbulence)
-{
-    if (mesh.changing())
-    {
-        y.correct();
-    }
-
-    tmp<volTensorField> tgradUb = fvc::grad(Ub);
-    volScalarField G(2*nutb*(tgradUb() && dev(symm(tgradUb()))));
-    tgradUb.clear();
-
-    #include "wallFunctions.H"
-
-    // Dissipation equation
-    fvScalarMatrix epsEqn
-    (
-        fvm::ddt(beta, epsilon)
-      + fvm::div(phib, epsilon)
-      - fvm::laplacian
-        (
-            alphaEps*nuEffb, epsilon,
-            "laplacian(DepsilonEff,epsilon)"
-        )
-      ==
-         C1*beta*G*epsilon/k
-       - fvm::Sp(C2*beta*epsilon/k, epsilon)
-    );
-
-    #include "wallDissipation.H"
-
-    epsEqn.relax();
-    epsEqn.solve();
-
-    epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
-
-
-    // Turbulent kinetic energy equation
-    fvScalarMatrix kEqn
-    (
-        fvm::ddt(beta, k)
-      + fvm::div(phib, k)
-      - fvm::laplacian
-        (
-            alphak*nuEffb, k,
-            "laplacian(DkEff,k)"
-        )
-      ==
-        beta*G
-      - fvm::Sp(beta*epsilon/k, k)
-    );
-    kEqn.relax();
-    kEqn.solve();
-
-    k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
-
-    //- Re-calculate turbulence viscosity
-    nutb = Cmu*sqr(k)/epsilon;
-
-    #include "wallViscosity.H"
-}
-
-nuEffb = nutb + nub;

applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H

@@ -1,5 +1,5 @@
     volVectorField Ur(Ua - Ub);
-    volScalarField magUr(mag(Ur));
+    volScalarField magUr(mag(Ur) + residualSlip);
 
     volScalarField Ka(draga->K(magUr));
     volScalarField K(Ka);

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C

@@ -849,7 +849,10 @@ autoDensity::autoDensity
 :
     initialPointsMethod(typeName, initialPointsDict, cvMesh),
     globalTrialPoints_(0),
-    minCellSizeLimit_(readScalar(detailsDict().lookup("minCellSizeLimit"))),
+    minCellSizeLimit_
+    (
+        detailsDict().lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
+    ),
     minLevels_(readLabel(detailsDict().lookup("minLevels"))),
     maxSizeRatio_(readScalar(detailsDict().lookup("maxSizeRatio"))),
     volRes_(readLabel(detailsDict().lookup("sampleResolution"))),
@@ -899,18 +902,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         );
     }
 
-    // Initialise size of points list.
+    DynamicList<Vb::Point> initialPoints;
-    const scalar volumeBoundBox = Foam::pow3(hierBB.typDim());
-    const scalar volumeSmallestCell = Foam::pow3(minCellSizeLimit_);
-
-    const int initialPointEstimate
-        = min
-          (
-              static_cast<int>(volumeBoundBox/(volumeSmallestCell + SMALL)/10),
-              1e6
-          );
-
-    DynamicList<Vb::Point> initialPoints(initialPointEstimate);
 
     Info<< nl << "    " << typeName << endl;
 

src/OpenFOAM/meshes/boundBox/boundBox.C

@@ -163,6 +163,49 @@ Foam::tmp<Foam::pointField> Foam::boundBox::points() const
 }
 
 
+Foam::faceList Foam::boundBox::faces()
+{
+    faceList faces(6);
+
+    forAll(faces, fI)
+    {
+        faces[fI].setSize(4);
+    }
+
+    faces[0][0] = 0;
+    faces[0][1] = 1;
+    faces[0][2] = 2;
+    faces[0][3] = 3;
+
+    faces[1][0] = 2;
+    faces[1][1] = 6;
+    faces[1][2] = 7;
+    faces[1][3] = 3;
+
+    faces[2][0] = 0;
+    faces[2][1] = 4;
+    faces[2][2] = 5;
+    faces[2][3] = 1;
+
+    faces[3][0] = 4;
+    faces[3][1] = 7;
+    faces[3][2] = 6;
+    faces[3][3] = 5;
+
+    faces[4][0] = 3;
+    faces[4][1] = 7;
+    faces[4][2] = 4;
+    faces[4][3] = 0;
+
+    faces[5][0] = 1;
+    faces[5][1] = 5;
+    faces[5][2] = 6;
+    faces[5][3] = 2;
+
+    return faces;
+}
+
+
 void Foam::boundBox::inflate(const scalar s)
 {
     vector ext = vector::one*s*mag();

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -33,6 +33,7 @@ Description
 #define boundBox_H
 
 #include "pointField.H"
+#include "faceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -160,6 +161,9 @@ public:
             //- Return corner points in an order corresponding to a 'hex' cell
             tmp<pointField> points() const;
 
+            //- Return faces with correct point order
+            static faceList faces();
+
 
         // Manipulate
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,6 +101,9 @@ Foam::ParticleErosion<CloudType>::ParticleErosion
     }
 
     patchIDs_ = uniquePatchIDs.toc();
+
+    // trigger ther creation of the Q field
+    preEvolve();
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,17 +103,16 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
 
     const scalar cps = 0.16432;
 
-    const volScalarField& k = *this->kPtr_;
+    const scalar k = this->kPtr_->internalField()[cellI];
-    const volScalarField& epsilon = *this->epsilonPtr_;
+    const scalar epsilon =
-    const volVectorField& gradk = *this->gradkPtr_;
+        this->epsilonPtr_->internalField()[cellI] + ROOTVSMALL;
+    const vector& gradk = this->gradkPtr_->internalField()[cellI];
 
     const scalar UrelMag = mag(U - Uc - UTurb);
 
-    const scalar tTurbLoc = min
+    const scalar tTurbLoc =
-    (
+        min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL));
-        k[cellI]/epsilon[cellI],
+
-        cps*pow(k[cellI], 1.5)/epsilon[cellI]/(UrelMag + SMALL)
-    );
 
     // Parcel is perturbed by the turbulence
     if (dt < tTurbLoc)
@@ -124,8 +123,8 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
         {
             tTurb = 0.0;
 
-            scalar sigma = sqrt(2.0*k[cellI]/3.0);
+            scalar sigma = sqrt(2.0*k/3.0);
-            vector dir = -gradk[cellI]/(mag(gradk[cellI]) + SMALL);
+            vector dir = -gradk/(mag(gradk) + SMALL);
 
             // Numerical Recipes... Ch. 7. Random Numbers...
             scalar x1 = 0.0;

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,16 +72,15 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
 
     const scalar cps = 0.16432;
 
-    const volScalarField& k = *this->kPtr_;
+    const scalar k = this->kPtr_->internalField()[cellI];
-    const volScalarField& epsilon = *this->epsilonPtr_;
+    const scalar epsilon =
+        this->epsilonPtr_->internalField()[cellI] + ROOTVSMALL;
 
     const scalar UrelMag = mag(U - Uc - UTurb);
 
-    const scalar tTurbLoc = min
+    const scalar tTurbLoc =
-    (
+        min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL));
-        k[cellI]/epsilon[cellI],
+
-        cps*pow(k[cellI], 1.5)/epsilon[cellI]/(UrelMag + SMALL)
-    );
 
     // Parcel is perturbed by the turbulence
     if (dt < tTurbLoc)
@@ -92,7 +91,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
         {
             tTurb = 0.0;
 
-            scalar sigma = sqrt(2.0*k[cellI]/3.0);
+            scalar sigma = sqrt(2.0*k/3.0);
             vector dir = 2.0*rnd.sample01<vector>() - vector::one;
             dir /= mag(dir) + SMALL;
 

tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties

@@ -21,5 +21,6 @@ dragModelb      GidaspowSchillerNaumann;
 
 dragPhase       a;
 
+residualSlip    0;
 
 // ************************************************************************* //

tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties

@@ -21,5 +21,6 @@ dragModelb      GidaspowErgunWenYu;
 
 dragPhase       a;
 
+residualSlip    0;
 
 // ************************************************************************* //

tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties

@@ -21,5 +21,6 @@ dragModelb      SchillerNaumann;
 
 dragPhase       blended;
 
+residualSlip    0;
 
 // ************************************************************************* //