diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H index ca29fcd0..c4806145 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H @@ -3,9 +3,11 @@ { fvScalarMatrix hEqn ( - thermo.isotropic() - ? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)") - : -fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)") + ( + thermo.isotropic() + ? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)") + : -fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)") + ) == fvOptions(rho, h) ); diff --git a/applications/solvers/multiphase/driftFluxFoam/Make/files b/applications/solvers/multiphase/driftFluxFoam/Make/files new file mode 100644 index 00000000..adbf0798 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/Make/files @@ -0,0 +1,3 @@ +driftFluxFoam.C + +EXE = $(FOAM_APPBIN)/driftFluxFoam diff --git a/applications/solvers/multiphase/driftFluxFoam/Make/options b/applications/solvers/multiphase/driftFluxFoam/Make/options new file mode 100644 index 00000000..ce089874 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/Make/options @@ -0,0 +1,11 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lfvOptions diff --git a/applications/solvers/multiphase/driftFluxFoam/UEqn.H b/applications/solvers/multiphase/driftFluxFoam/UEqn.H new file mode 100644 index 00000000..b6d87564 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/UEqn.H @@ -0,0 +1,38 @@ + // Solve the Momentum equation + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(rhoPhi, U) + + fvc::div + ( + (alpha/(scalar(1.001) - alpha))*((rhoc*rhod)/rho)*Vdj*Vdj, + "div(phiVdj,Vdj)" + ) + - fvm::laplacian(muEff, U) + - fvc::div(muEff*dev2(T(fvc::grad(U)))) + == + fvOptions(rho, U) + ); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ) + ); + + fvOptions.correct(U); + } diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaControls.H b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H new file mode 100644 index 00000000..5a3d10bd --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/alphaControls.H @@ -0,0 +1,4 @@ +const dictionary& alphaControls = mesh.solverDict(alpha.name()); + +label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr"))); +label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H new file mode 100644 index 00000000..89c139df --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H @@ -0,0 +1,31 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + for (int aCorr=0; aCorr 1) + { + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField phiAlphaSum + ( + IOobject + ( + "phiAlphaSum", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("0", phi.dimensions(), 0) + ); + + for + ( + subCycle alphaSubCycle(alpha, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqn.H" + phiAlphaSum += (runTime.deltaT()/totalDeltaT)*phiAlpha; + } + + phiAlpha = phiAlphaSum; + } + else + { + #include "alphaEqn.H" + } + + // Apply the diffusion term separately to allow implicit solution + // and boundedness of the explicit advection + { + fvScalarMatrix alphaEqn + ( + fvm::ddt(alpha) - fvc::ddt(alpha) + - fvm::laplacian(mut/rho, alpha) + ); + + alphaEqn.solve(); + + phiAlpha += alphaEqn.flux(); + } + + Info<< "Phase-1 volume fraction = " + << alpha.weightedAverage(mesh.Vsc()).value() + << " Min(alpha) = " << min(alpha).value() + << " Max(alpha) = " << max(alpha).value() + << endl; + + rhoPhi = phiAlpha*(rhod - rhoc) + phi*rhoc; + rho == alpha*rhod + (scalar(1) - alpha)*rhoc; +} diff --git a/applications/solvers/multiphase/driftFluxFoam/calcVdj.H b/applications/solvers/multiphase/driftFluxFoam/calcVdj.H new file mode 100644 index 00000000..7a921ef7 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/calcVdj.H @@ -0,0 +1,20 @@ +if (VdjModel == "general") +{ + Vdj = V0* + ( + exp(-a*max(alpha - alphaMin, scalar(0))) + - exp(-a1*max(alpha - alphaMin, scalar(0))) + ); +} +else if (VdjModel == "simple") +{ + Vdj = V0*pow(10.0, -a*max(alpha, scalar(0))); +} +else +{ + FatalErrorIn(args.executable()) + << "Unknown VdjModel : " << VdjModel + << abort(FatalError); +} + +Vdj.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H b/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H new file mode 100644 index 00000000..0bf086a4 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/correctViscosity.H @@ -0,0 +1,39 @@ +{ + mul = muc + + plasticViscosity + ( + plasticViscosityCoeff, + plasticViscosityExponent, + alpha + ); + + if (BinghamPlastic) + { + volScalarField tauy = yieldStress + ( + yieldStressCoeff, + yieldStressExponent, + yieldStressOffset, + alpha + ); + + mul = + tauy/ + ( + mag(fvc::grad(U)) + + 1.0e-4* + ( + tauy + + dimensionedScalar + ( + "deltaTauy", + tauy.dimensions(), + 1.0e-15 + ) + )/mul + ) + + mul; + } + + mul = min(mul, muMax); +} diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H new file mode 100644 index 00000000..950235cf --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H @@ -0,0 +1,384 @@ + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field alpha\n" << endl; + volScalarField alpha + ( + IOobject + ( + "alpha", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "createPhi.H" + + + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + + dimensionedScalar rhoc(transportProperties.lookup("rhoc")); + + dimensionedScalar rhod(transportProperties.lookup("rhod")); + + dimensionedScalar muc(transportProperties.lookup("muc")); + dimensionedScalar muMax(transportProperties.lookup("muMax")); + + dimensionedScalar plasticViscosityCoeff + ( + transportProperties.lookup("plasticViscosityCoeff") + ); + + dimensionedScalar plasticViscosityExponent + ( + transportProperties.lookup("plasticViscosityExponent") + ); + + dimensionedScalar yieldStressCoeff + ( + transportProperties.lookup("yieldStressCoeff") + ); + + dimensionedScalar yieldStressExponent + ( + transportProperties.lookup("yieldStressExponent") + ); + + dimensionedScalar yieldStressOffset + ( + transportProperties.lookup("yieldStressOffset") + ); + + Switch BinghamPlastic(transportProperties.lookup("BinghamPlastic")); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + alpha*rhod + (scalar(1) - alpha)*rhoc + ); + rho.oldTime(); + + // Mass flux + surfaceScalarField rhoPhi + ( + IOobject + ( + "rhoPhi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::interpolate(rho)*phi + ); + + Info<< "Calculating field mul\n" << endl; + volScalarField mul + ( + IOobject + ( + "mul", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + muc + + plasticViscosity + ( + plasticViscosityCoeff, + plasticViscosityExponent, + alpha + ) + ); + + + Info<< "Initialising field Vdj\n" << endl; + volVectorField Vdj + ( + IOobject + ( + "Vdj", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("0.0", U.dimensions(), vector::zero), + U.boundaryField().types() + ); + + + Info<< "Selecting Drift-Flux model " << endl; + + const word VdjModel(transportProperties.lookup("VdjModel")); + + Info<< tab << VdjModel << " selected\n" << endl; + + const dictionary& VdjModelCoeffs + ( + transportProperties.subDict(VdjModel + "Coeffs") + ); + + dimensionedVector V0(VdjModelCoeffs.lookup("V0")); + + dimensionedScalar a(VdjModelCoeffs.lookup("a")); + + dimensionedScalar a1(VdjModelCoeffs.lookup("a1")); + + dimensionedScalar alphaMin(VdjModelCoeffs.lookup("alphaMin")); + + + IOdictionary RASProperties + ( + IOobject + ( + "RASProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + + Switch turbulence(RASProperties.lookup("turbulence")); + + dictionary kEpsilonDict(RASProperties.subDictPtr("kEpsilonCoeffs")); + + dimensionedScalar Cmu + ( + dimensionedScalar::lookupOrAddToDict + ( + "Cmu", + kEpsilonDict, + 0.09 + ) + ); + + dimensionedScalar C1 + ( + dimensionedScalar::lookupOrAddToDict + ( + "C1", + kEpsilonDict, + 1.44 + ) + ); + + dimensionedScalar C2 + ( + dimensionedScalar::lookupOrAddToDict + ( + "C2", + kEpsilonDict, + 1.92 + ) + ); + + dimensionedScalar C3 + ( + dimensionedScalar::lookupOrAddToDict + ( + "C3", + kEpsilonDict, + 0.85 + ) + ); + + dimensionedScalar sigmak + ( + dimensionedScalar::lookupOrAddToDict + ( + "sigmak", + kEpsilonDict, + 1.0 + ) + ); + + dimensionedScalar sigmaEps + ( + dimensionedScalar::lookupOrAddToDict + ( + "sigmaEps", + kEpsilonDict, + 1.3 + ) + ); + + dictionary wallFunctionDict(RASProperties.subDictPtr("wallFunctionCoeffs")); + + dimensionedScalar kappa + ( + dimensionedScalar::lookupOrAddToDict + ( + "kappa", + wallFunctionDict, + 0.41 + ) + ); + + dimensionedScalar E + ( + dimensionedScalar::lookupOrAddToDict + ( + "E", + wallFunctionDict, + 9.8 + ) + ); + + if (RASProperties.lookupOrDefault("printCoeffs", false)) + { + Info<< "kEpsilonCoeffs" << kEpsilonDict << nl + << "wallFunctionCoeffs" << wallFunctionDict << endl; + } + + + nearWallDist y(mesh); + + Info<< "Reading field k\n" << endl; + volScalarField k + ( + IOobject + ( + "k", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field epsilon\n" << endl; + volScalarField epsilon + ( + IOobject + ( + "epsilon", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Calculating field mut\n" << endl; + volScalarField mut + ( + IOobject + ( + "mut", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + Cmu*rho*sqr(k)/epsilon + ); + + + Info<< "Calculating field muEff\n" << endl; + volScalarField muEff + ( + IOobject + ( + "muEff", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mut + mul + ); + + + Info<< "Calculating field (g.h)f\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rho*gh + ); + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell + ( + p, + p_rgh, + mesh.solutionDict().subDict("PIMPLE"), + pRefCell, + pRefValue + ); + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; + } diff --git a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C new file mode 100644 index 00000000..d0aafc6e --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C @@ -0,0 +1,114 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2014 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 + driftFluxFoam + +Description + Solver for 2 incompressible fluids using the mixture approach with the + drift-flux approximation for relative motion of the phases. + + Used for simulating the settling of the dispersed phase and other similar + separation problems. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "CMULES.H" +#include "subCycle.H" +#include "nearWallDist.H" +#include "wallFvPatch.H" +#include "bound.H" +#include "Switch.H" +#include "plasticViscosity.H" +#include "yieldStress.H" +#include "pimpleControl.H" +#include "fvIOoptionList.H" +#include "fixedFluxPressureFvPatchScalarField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + + pimpleControl pimple(mesh); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "alphaControls.H" + + #include "calcVdj.H" + #include "alphaEqnSubCycle.H" + #include "correctViscosity.H" + + #include "UEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + #include "kEpsilon.H" + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H new file mode 100644 index 00000000..5aaccab8 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H @@ -0,0 +1,78 @@ +if (turbulence) +{ + if (mesh.changing()) + { + y.correct(); + } + + dimensionedScalar k0("k0", k.dimensions(), 0); + dimensionedScalar kMin("kMin", k.dimensions(), SMALL); + dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0); + dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL); + + volScalarField divU(fvc::div(rhoPhi/fvc::interpolate(rho))); + + tmp tgradU = fvc::grad(U); + volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU())))); + tgradU.clear(); + + volScalarField Gcoef + ( + Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin) + ); + + #include "wallFunctions.H" + + // Dissipation equation + fvScalarMatrix epsEqn + ( + fvm::ddt(rho, epsilon) + + fvm::div(rhoPhi, epsilon) + - fvm::laplacian + ( + mut/sigmaEps + muc, epsilon, + "laplacian(DepsilonEff,epsilon)" + ) + == + C1*G*epsilon/(k + kMin) + - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon) + - fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon) + ); + + #include "wallDissipation.H" + + epsEqn.relax(); + epsEqn.solve(); + + bound(epsilon, epsilon0); + + + // Turbulent kinetic energy equation + fvScalarMatrix kEqn + ( + fvm::ddt(rho, k) + + fvm::div(rhoPhi, k) + - fvm::laplacian + ( + mut/sigmak + muc, k, + "laplacian(DkEff,k)" + ) + == + G + - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k) + - fvm::Sp(rho*epsilon/(k + kMin), k) + ); + + kEqn.relax(); + kEqn.solve(); + + bound(k, k0); + + + //- Re-calculate viscosity + mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin); + + #include "wallViscosity.H" +} + +muEff = mut + mul; diff --git a/applications/solvers/multiphase/driftFluxFoam/pEqn.H b/applications/solvers/multiphase/driftFluxFoam/pEqn.H new file mode 100644 index 00000000..d32683bc --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H @@ -0,0 +1,73 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + adjustPhi(phiHbyA, U, p_rgh); + fvOptions.makeRelative(phiHbyA); + + surfaceScalarField phig + ( + ( + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + // Update the fixedFluxPressure BCs to ensure flux consistency + setSnGrad + ( + p_rgh.boundaryField(), + ( + phiHbyA.boundaryField() + - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField()) + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqn + ( + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) + ); + + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA - p_rghEqn.flux(); + + p_rgh.relax(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); + fvOptions.correct(U); + } + } + + #include "continuityErrs.H" + + p == p_rgh + rho*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; + } +} diff --git a/applications/solvers/multiphase/driftFluxFoam/plasticViscosity.H b/applications/solvers/multiphase/driftFluxFoam/plasticViscosity.H new file mode 100644 index 00000000..4be4e6ea --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/plasticViscosity.H @@ -0,0 +1,21 @@ +volScalarField plasticViscosity +( + const dimensionedScalar& plasticViscosityCoeff, + const dimensionedScalar& plasticViscosityExponent, + const volScalarField& Alpha +) +{ + tmp tfld + ( + plasticViscosityCoeff* + ( + pow + ( + 10.0, + plasticViscosityExponent*Alpha + SMALL + ) - scalar(1) + ) + ); + + return tfld(); +} diff --git a/applications/solvers/multiphase/driftFluxFoam/wallDissipation.H b/applications/solvers/multiphase/driftFluxFoam/wallDissipation.H new file mode 100644 index 00000000..c8b88597 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/wallDissipation.H @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2014 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 . + +Global + wallDissipation + +Description + Set wall dissipation in the epsilon matrix + +\*---------------------------------------------------------------------------*/ + +{ + const fvPatchList& patches = mesh.boundary(); + + forAll(patches, patchi) + { + const fvPatch& p = patches[patchi]; + + if (isA(p)) + { + epsEqn.setValues + ( + p.faceCells(), + epsilon.boundaryField()[patchi].patchInternalField() + ); + } + } +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H new file mode 100644 index 00000000..149d787c --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H @@ -0,0 +1,85 @@ +{ + labelList cellBoundaryFaceCount(epsilon.size(), 0); + + const scalar Cmu25 = ::pow(Cmu.value(), 0.25); + const scalar Cmu75 = ::pow(Cmu.value(), 0.75); + const scalar kappa_ = kappa.value(); + const scalar muc_ = muc.value(); + + const fvPatchList& patches = mesh.boundary(); + + //- Initialise the near-wall P field to zero + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isA(curPatch)) + { + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + epsilon[faceCelli] = 0.0; + G[faceCelli] = 0.0; + } + } + } + + //- Accumulate the wall face contributions to epsilon and G + // Increment cellBoundaryFaceCount for each face for averaging + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isA(curPatch)) + { + const scalarField& mutw = mut.boundaryField()[patchi]; + + scalarField magFaceGradU + ( + mag(U.boundaryField()[patchi].snGrad()) + ); + + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + // For corner cells (with two boundary or more faces), + // epsilon and G in the near-wall cell are calculated + // as an average + + cellBoundaryFaceCount[faceCelli]++; + + epsilon[faceCelli] += + Cmu75*::pow(k[faceCelli], 1.5) + /(kappa_*y[patchi][facei]); + + G[faceCelli] += + (mutw[facei] + muc_) + *magFaceGradU[facei] + *Cmu25*::sqrt(k[faceCelli]) + /(kappa_*y[patchi][facei]); + } + } + } + + + // perform the averaging + + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isA(curPatch)) + { + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + G[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + cellBoundaryFaceCount[faceCelli] = 1; + } + } + } +} diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H new file mode 100644 index 00000000..d63f7e79 --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H @@ -0,0 +1,38 @@ +{ + const scalar Cmu25 = ::pow(Cmu.value(), 0.25); + const scalar kappa_ = kappa.value(); + const scalar E_ = E.value(); + const scalar muc_ = muc.value(); + const scalar nuc_ = muc_/rhoc.value(); + + const fvPatchList& patches = mesh.boundary(); + + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isA(curPatch)) + { + scalarField& mutw = mut.boundaryField()[patchi]; + + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + scalar yPlus = + Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) + /nuc_; + + if (yPlus > 11.6) + { + mutw[facei] = + muc_*(yPlus*kappa_/::log(E_*yPlus) - 1); + } + else + { + mutw[facei] = 0.0; + } + } + } + } +} diff --git a/applications/solvers/multiphase/driftFluxFoam/yieldStress.H b/applications/solvers/multiphase/driftFluxFoam/yieldStress.H new file mode 100644 index 00000000..eda55d2e --- /dev/null +++ b/applications/solvers/multiphase/driftFluxFoam/yieldStress.H @@ -0,0 +1,27 @@ +volScalarField yieldStress +( + const dimensionedScalar& yieldStressCoeff, + const dimensionedScalar& yieldStressExponent, + const dimensionedScalar& yieldStressOffset, + const volScalarField& alpha +) +{ + tmp tfld + ( + yieldStressCoeff* + ( + pow + ( + 10.0, + yieldStressExponent*(max(alpha, scalar(0)) + yieldStressOffset) + ) + - pow + ( + 10.0, + yieldStressExponent*yieldStressOffset + ) + ) + ); + + return tfld(); +} diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C index 37c543a8..f83cdd55 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/noBlending/noBlending.C @@ -88,7 +88,7 @@ Foam::tmp Foam::blendingMethods::noBlending::f1 ( "f", dimless, - phase1.name() == continuousPhase_ + phase2.name() != continuousPhase_ ) ) ); @@ -101,7 +101,28 @@ Foam::tmp Foam::blendingMethods::noBlending::f2 const phaseModel& phase2 ) const { - return f1(phase1, phase2); + const fvMesh& mesh(phase1.mesh()); + + return + tmp + ( + new volScalarField + ( + IOobject + ( + "f", + mesh.time().timeName(), + mesh + ), + mesh, + dimensionedScalar + ( + "f", + dimless, + phase1.name() == continuousPhase_ + ) + ) + ); } diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C index 267b8dc8..bae7d395 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C +++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C @@ -195,41 +195,44 @@ void Foam::ensightMesh::correct() { forAll(mesh_.boundary(), patchi) { - const polyPatch& p = mesh_.boundaryMesh()[patchi]; - - labelList& tris = boundaryFaceSets_[patchi].tris; - labelList& quads = boundaryFaceSets_[patchi].quads; - labelList& polys = boundaryFaceSets_[patchi].polys; - - tris.setSize(p.size()); - quads.setSize(p.size()); - polys.setSize(p.size()); - - label nTris = 0; - label nQuads = 0; - label nPolys = 0; - - forAll(p, faceI) + if (mesh_.boundary()[patchi].size()) { - const face& f = p[faceI]; + const polyPatch& p = mesh_.boundaryMesh()[patchi]; - if (f.size() == 3) + labelList& tris = boundaryFaceSets_[patchi].tris; + labelList& quads = boundaryFaceSets_[patchi].quads; + labelList& polys = boundaryFaceSets_[patchi].polys; + + tris.setSize(p.size()); + quads.setSize(p.size()); + polys.setSize(p.size()); + + label nTris = 0; + label nQuads = 0; + label nPolys = 0; + + forAll(p, faceI) { - tris[nTris++] = faceI; - } - else if (f.size() == 4) - { - quads[nQuads++] = faceI; - } - else - { - polys[nPolys++] = faceI; + const face& f = p[faceI]; + + if (f.size() == 3) + { + tris[nTris++] = faceI; + } + else if (f.size() == 4) + { + quads[nQuads++] = faceI; + } + else + { + polys[nPolys++] = faceI; + } } + + tris.setSize(nTris); + quads.setSize(nQuads); + polys.setSize(nPolys); } - - tris.setSize(nTris); - quads.setSize(nQuads); - polys.setSize(nPolys); } } @@ -240,9 +243,12 @@ void Foam::ensightMesh::correct() if (patchNames_.empty() || patchNames_.found(patchName)) { - nfp.nTris = boundaryFaceSets_[patchi].tris.size(); - nfp.nQuads = boundaryFaceSets_[patchi].quads.size(); - nfp.nPolys = boundaryFaceSets_[patchi].polys.size(); + if (mesh_.boundary()[patchi].size()) + { + nfp.nTris = boundaryFaceSets_[patchi].tris.size(); + nfp.nQuads = boundaryFaceSets_[patchi].quads.size(); + nfp.nPolys = boundaryFaceSets_[patchi].polys.size(); + } } reduce(nfp.nTris, sumOp