diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C index 4bbaf792af..bef9c2b0ab 100644 --- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C +++ b/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; } } diff --git a/applications/solvers/multiphase/bubbleFoam/kEpsilon.H b/applications/solvers/multiphase/bubbleFoam/kEpsilon.H index 309d4eba80..05cfb08c9d 100644 --- a/applications/solvers/multiphase/bubbleFoam/kEpsilon.H +++ b/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 - - fvm::Sp(C2*beta*epsilon/k, epsilon) + C1*G*epsilon/k + - 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 - - fvm::Sp(beta*epsilon/k, k) + G + - 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; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean index cc138bc068..63db39ff05 100755 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean @@ -6,5 +6,6 @@ wclean libso phaseModel wclean libso interfacialModels wclean libso kineticTheoryModels wclean +wclean MRFtwoPhaseEulerFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake index 29294d166a..faf438d0bd 100755 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake @@ -6,5 +6,6 @@ wmake libso phaseModel wmake libso interfacialModels wmake libso kineticTheoryModels wmake +wmake MRFtwoPhaseEulerFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C new file mode 100644 index 0000000000..357c7e413a --- /dev/null +++ b/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 . + +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; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files new file mode 100644 index 0000000000..45960722ae --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files @@ -0,0 +1,3 @@ +MRFTwoPhaseEulerFoam.C + +EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options new file mode 100644 index 0000000000..b9b19059da --- /dev/null +++ b/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 diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H new file mode 100644 index 0000000000..2aba5b9d3d --- /dev/null +++ b/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(); + } +} diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H new file mode 100644 index 0000000000..0d65353624 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H @@ -0,0 +1,3 @@ + MRFZones mrfZones(mesh); + mrfZones.correctBoundaryVelocity(Ua); + mrfZones.correctBoundaryVelocity(Ub); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H new file mode 100644 index 0000000000..d89136b172 --- /dev/null +++ b/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(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" diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H index 0a590e6bf9..a4d428e176 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H +++ b/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 diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index 95395d97d4..076b14a713 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/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 ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H b/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H deleted file mode 100644 index 8f85e12052..0000000000 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H +++ /dev/null @@ -1,62 +0,0 @@ -if (turbulence) -{ - if (mesh.changing()) - { - y.correct(); - } - - tmp 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; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H b/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H index b614a4dfbe..27519eaf2e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H +++ b/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); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C index 84926910df..33de87c5bc 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C +++ b/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("minCellSizeLimit", 0.0) + ), minLevels_(readLabel(detailsDict().lookup("minLevels"))), maxSizeRatio_(readScalar(detailsDict().lookup("maxSizeRatio"))), volRes_(readLabel(detailsDict().lookup("sampleResolution"))), @@ -899,18 +902,7 @@ List autoDensity::initialPoints() const ); } - // Initialise size of points list. - const scalar volumeBoundBox = Foam::pow3(hierBB.typDim()); - const scalar volumeSmallestCell = Foam::pow3(minCellSizeLimit_); - - const int initialPointEstimate - = min - ( - static_cast(volumeBoundBox/(volumeSmallestCell + SMALL)/10), - 1e6 - ); - - DynamicList initialPoints(initialPointEstimate); + DynamicList initialPoints; Info<< nl << " " << typeName << endl; diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C index aaaa528be2..f1c297837d 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.C +++ b/src/OpenFOAM/meshes/boundBox/boundBox.C @@ -163,6 +163,49 @@ Foam::tmp 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(); diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H index f531283c3c..c59281bdf0 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.H +++ b/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 points() const; + //- Return faces with correct point order + static faceList faces(); + // Manipulate diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C index c2bbf94493..4d4f4dc3ca 100644 --- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C +++ b/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::ParticleErosion } patchIDs_ = uniquePatchIDs.toc(); + + // trigger ther creation of the Q field + preEvolve(); } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C b/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C index b421a92b3f..c96a5c5b3d 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C +++ b/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::update const scalar cps = 0.16432; - const volScalarField& k = *this->kPtr_; - const volScalarField& epsilon = *this->epsilonPtr_; - const volVectorField& gradk = *this->gradkPtr_; + const scalar k = this->kPtr_->internalField()[cellI]; + const scalar epsilon = + this->epsilonPtr_->internalField()[cellI] + ROOTVSMALL; + const vector& gradk = this->gradkPtr_->internalField()[cellI]; const scalar UrelMag = mag(U - Uc - UTurb); - const scalar tTurbLoc = min - ( - k[cellI]/epsilon[cellI], - cps*pow(k[cellI], 1.5)/epsilon[cellI]/(UrelMag + SMALL) - ); + const scalar tTurbLoc = + min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL)); + // Parcel is perturbed by the turbulence if (dt < tTurbLoc) @@ -124,8 +123,8 @@ Foam::vector Foam::GradientDispersionRAS::update { tTurb = 0.0; - scalar sigma = sqrt(2.0*k[cellI]/3.0); - vector dir = -gradk[cellI]/(mag(gradk[cellI]) + SMALL); + scalar sigma = sqrt(2.0*k/3.0); + vector dir = -gradk/(mag(gradk) + SMALL); // Numerical Recipes... Ch. 7. Random Numbers... scalar x1 = 0.0; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C index 73f0cb9553..d3c6c5b739 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C +++ b/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::update const scalar cps = 0.16432; - const volScalarField& k = *this->kPtr_; - const volScalarField& epsilon = *this->epsilonPtr_; + const scalar k = this->kPtr_->internalField()[cellI]; + const scalar epsilon = + this->epsilonPtr_->internalField()[cellI] + ROOTVSMALL; const scalar UrelMag = mag(U - Uc - UTurb); - const scalar tTurbLoc = min - ( - k[cellI]/epsilon[cellI], - cps*pow(k[cellI], 1.5)/epsilon[cellI]/(UrelMag + SMALL) - ); + const scalar tTurbLoc = + min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL)); + // Parcel is perturbed by the turbulence if (dt < tTurbLoc) @@ -92,7 +91,7 @@ Foam::vector Foam::StochasticDispersionRAS::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::one; dir /= mag(dir) + SMALL; diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties index 32a9cb4441..0e47a29149 100644 --- a/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties +++ b/tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties @@ -21,5 +21,6 @@ dragModelb GidaspowSchillerNaumann; dragPhase a; +residualSlip 0; // ************************************************************************* // diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties index 3c4c17acd1..4896ae1bfd 100644 --- a/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties +++ b/tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties @@ -21,5 +21,6 @@ dragModelb GidaspowErgunWenYu; dragPhase a; +residualSlip 0; // ************************************************************************* // diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties index f106b0ba8d..51e5a32fe8 100644 --- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties +++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties @@ -21,5 +21,6 @@ dragModelb SchillerNaumann; dragPhase blended; +residualSlip 0; // ************************************************************************* //