diff --git a/applications/solvers/multiphase/interIsoFoam/Make/files b/applications/solvers/multiphase/interIsoFoam/Make/files new file mode 100644 index 0000000000..10f61a6885 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/Make/files @@ -0,0 +1,3 @@ +interIsoFoam.C + +EXE = $(FOAM_APPBIN)/interIsoFoam diff --git a/applications/solvers/multiphase/interIsoFoam/Make/options b/applications/solvers/multiphase/interIsoFoam/Make/options new file mode 100644 index 0000000000..cfe4a92072 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/Make/options @@ -0,0 +1,20 @@ +EXE_INC = \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -limmiscibleIncompressibleTwoPhaseMixture \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling diff --git a/applications/solvers/multiphase/interIsoFoam/UEqn.H b/applications/solvers/multiphase/interIsoFoam/UEqn.H new file mode 100644 index 0000000000..77d1dcd83e --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/UEqn.H @@ -0,0 +1,33 @@ + MRF.correctBoundaryVelocity(U); + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(rho, U) + == + fvOptions(rho, U) + ); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + fvOptions.correct(U); + } diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H new file mode 100644 index 0000000000..db77d94af4 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H @@ -0,0 +1,3 @@ +const dictionary& alphaControls = mesh.solverDict(alpha1.name()); + +label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); diff --git a/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H new file mode 100644 index 0000000000..a084155c8a --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2016 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 + alphaCourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar maxAlphaCo +( + readScalar(runTime.controlDict().lookup("maxAlphaCo")) +); + +scalar alphaCoNum = 0.0; +scalar meanAlphaCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + mixture.nearInterface()().primitiveField() + *fvc::surfaceSum(mag(phi))().primitiveField() + ); + + alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanAlphaCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Interface Courant Number mean: " << meanAlphaCoNum + << " max: " << alphaCoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H new file mode 100644 index 0000000000..cefabfda73 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H @@ -0,0 +1,35 @@ +// If there are more than one outer corrector, we use a mixture of old and +// new U and phi for propagating alpha1 in all but the first outer iteration +if (!pimple.firstIter()) +{ + // We are recalculating alpha1 from the its old time value + alpha1 = alpha1.oldTime(); + // Temporarily storing new U and phi values in prevIter storage + U.storePrevIter(); + phi.storePrevIter(); + + // Overwriting new U and phi values with mixture of old and new values + phi = 0.5*(phi + phi.oldTime()); + U = 0.5*(U + U.oldTime()); +} + +// Update alpha1 +advector.advect(); + +// Update rhoPhi +rhoPhi = advector.getRhoPhi(rho1, rho2); + +alpha2 = 1.0 - alpha1; + +if (!pimple.firstIter()) +{ + // Restoring new U and phi values temporarily saved in prevIter() above + U = U.prevIter(); + phi = phi.prevIter(); +} + +Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(" << alpha1.name() << ") = " << min(alpha1).value() + << " Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1 + << endl; diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H new file mode 100644 index 0000000000..8f0af80e0d --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H @@ -0,0 +1,35 @@ +if (nAlphaSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum + ( + IOobject + ( + "rhoPhiSum", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("0", rhoPhi.dimensions(), 0) + ); + + tmp trSubDeltaT; + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ + #include "alphaEqn.H" +} + +rho == alpha1*rho1 + alpha2*rho2; diff --git a/applications/solvers/multiphase/interIsoFoam/correctPhi.H b/applications/solvers/multiphase/interIsoFoam/correctPhi.H new file mode 100644 index 0000000000..9afcd58a66 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/correctPhi.H @@ -0,0 +1,11 @@ +CorrectPhi +( + U, + phi, + p_rgh, + dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1), + geometricZeroField(), + pimple +); + +#include "continuityErrs.H" diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H new file mode 100644 index 0000000000..8c3a05e3fa --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/createFields.H @@ -0,0 +1,123 @@ +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 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; +immiscibleIncompressibleTwoPhaseMixture mixture(U, phi); + +volScalarField& alpha1(mixture.alpha1()); +volScalarField& alpha2(mixture.alpha2()); + +const dimensionedScalar& rho1 = mixture.rho1(); +const dimensionedScalar& rho2 = mixture.rho2(); + + +// Need to store rho for ddt(rho, U) +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + alpha1*rho1 + alpha2*rho2 +); +rho.oldTime(); + + +// Mass flux +surfaceScalarField rhoPhi +( + IOobject + ( + "rhoPhi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::interpolate(rho)*phi +); + + +// Construct incompressible turbulence model +autoPtr turbulence +( + incompressible::turbulenceModel::New(U, phi, mixture) +); + + +#include "readGravitationalAcceleration.H" +#include "readhRef.H" +#include "gh.H" + + +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, + pimple.dict(), + pRefCell, + pRefValue +); + +if (p_rgh.needReference()) +{ + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; +} + +mesh.setFluxRequired(p_rgh.name()); +mesh.setFluxRequired(alpha1.name()); + +#include "createMRF.H" +#include "createIsoAdvection.H" diff --git a/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H new file mode 100644 index 0000000000..889af7228f --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H @@ -0,0 +1,10 @@ +// Defining isoAdvection +isoAdvection advector(alpha1, phi, U); + +bool frozenFlow = pimple.dict().lookupOrDefault("frozenFlow", false); +if (frozenFlow) +{ + Info<< "Employing frozen-flow assumption: " + << "pressure-velocity system will not be updated" + << endl; +} diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C new file mode 100644 index 0000000000..1514c5c0cb --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C @@ -0,0 +1,143 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016 DHI +------------------------------------------------------------------------------- +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 + interIsoFoam + +Group + grpMultiphaseSolvers + +Description + Solver derived from interFoam for 2 incompressible, isothermal immiscible + fluids using the iso-advector phase-fraction based interface capturing + approach. + + The momentum and other fluid properties are of the "mixture" and a single + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + For a two-fluid approach see twoPhaseEulerFoam. + + Reference: + \verbatim + Roenby, J., Bredmose, H. and Jasak, H. (2016). + A computational method for sharp interface advection + Royal Society Open Science, 3 + doi 10.1098/rsos.160405 + \endverbatim + + isoAdvector code supplied by Johan Roenby, DHI (2016) + +\*---------------------------------------------------------------------------*/ + +#include "isoAdvection.H" +#include "fvCFD.H" +#include "subCycle.H" +#include "immiscibleIncompressibleTwoPhaseMixture.H" +#include "turbulentTransportModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "CorrectPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "correctPhi.H" + + turbulence->validate(); + + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "alphaControls.H" + #include "alphaEqnSubCycle.H" + + mixture.correct(); + + if (frozenFlow) + { + continue; + } + + #include "UEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + 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/interIsoFoam/pEqn.H b/applications/solvers/multiphase/interIsoFoam/pEqn.H new file mode 100644 index 0000000000..17c3ae9748 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/pEqn.H @@ -0,0 +1,67 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::flux(HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + MRF.makeRelative(phiHbyA); + adjustPhi(phiHbyA, U, p_rgh); + + surfaceScalarField phig + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() +// - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf + ); + + phiHbyA += phig; + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF); + + 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/interIsoFoam/setDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H new file mode 100644 index 0000000000..9cc860c032 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 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 + setDeltaT + +Description + Reset the timestep to maintain a constant maximum courant Number. + Reduction of time-step is immediate, but increase is damped to avoid + unstable oscillations. + +\*---------------------------------------------------------------------------*/ + +if (adjustTimeStep) +{ + scalar maxDeltaTFact = + min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL)); + + scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); + + runTime.setDeltaT + ( + min + ( + deltaTFact*runTime.deltaTValue(), + maxDeltaT + ) + ); + + Info<< "deltaT = " << runTime.deltaTValue() << endl; +} + +// ************************************************************************* // diff --git a/applications/test/parallelOverset/heatTransfer/0.orig/T b/applications/test/parallelOverset/heatTransfer/0.orig/T index a60edc4bb6..c775c6f0c5 100644 --- a/applications/test/parallelOverset/heatTransfer/0.orig/T +++ b/applications/test/parallelOverset/heatTransfer/0.orig/T @@ -20,8 +20,7 @@ internalField uniform 273; boundaryField { - //- Set patchGroups for constraint patches - #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes" + #includeEtc "caseDicts/setConstraintTypes" "(walls|hole)" { diff --git a/applications/test/parallelOverset/heatTransfer/0.orig/zoneID b/applications/test/parallelOverset/heatTransfer/0.orig/zoneID index 1a49886c28..3e3f18697d 100644 --- a/applications/test/parallelOverset/heatTransfer/0.orig/zoneID +++ b/applications/test/parallelOverset/heatTransfer/0.orig/zoneID @@ -21,8 +21,7 @@ internalField uniform 0; boundaryField { - //- Set patchGroups for constraint patches - #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes" + #includeEtc "caseDicts/setConstraintTypes" overset { diff --git a/applications/utilities/mesh/generation/snappyHexMesh/meshQualityDict b/applications/utilities/mesh/generation/snappyHexMesh/meshQualityDict index 9cf0a3e506..3d43ee3a28 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/meshQualityDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/meshQualityDict @@ -15,7 +15,7 @@ FoamFile // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Include defaults parameters from master dictionary -#include "$WM_PROJECT_DIR/etc/caseDicts/meshQualityDict" +#includeEtc "caseDicts/meshQualityDict" // ************************************************************************* // diff --git a/applications/utilities/preProcessing/setAlphaField/Make/files b/applications/utilities/preProcessing/setAlphaField/Make/files new file mode 100755 index 0000000000..23145059e2 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/Make/files @@ -0,0 +1,3 @@ +setAlphaField.C + +EXE = $(FOAM_APPBIN)/setAlphaField diff --git a/applications/utilities/preProcessing/setAlphaField/Make/options b/applications/utilities/preProcessing/setAlphaField/Make/options new file mode 100755 index 0000000000..b4ad3d8bdf --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C new file mode 100644 index 0000000000..6fe8001e31 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C @@ -0,0 +1,191 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016-2017 DHI +------------------------------------------------------------------------------- +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 + setAlphaField + +Description + Uses isoCutCell to create a volume fraction field from either a cylinder, + a sphere or a plane. + + Original code supplied by Johan Roenby, DHI (2016) + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "isoCutFace.H" +#include "isoCutCell.H" +#include "Enum.H" +#include "mathematicalConstants.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class shapeSelector +{ + public: + + enum class shapeType + { + PLANE, + SPHERE, + CYLINDER, + SIN + }; + + static const Foam::Enum shapeTypeNames; +}; + +const Foam::Enum shapeSelector::shapeTypeNames +{ + { shapeSelector::shapeType::PLANE, "plane" }, + { shapeSelector::shapeType::SPHERE, "sphere" }, + { shapeSelector::shapeType::CYLINDER, "cylinder" }, + { shapeSelector::shapeType::SIN, "sin" }, +}; + + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + Info<< "Reading setAlphaFieldDict\n" << endl; + + IOdictionary dict + ( + IOobject + ( + "setAlphaFieldDict", + runTime.system(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + const shapeSelector::shapeType surfType + ( + shapeSelector::shapeTypeNames.read(dict.lookup("type")) + ); + const vector centre(dict.lookup("centre")); + const word fieldName(dict.lookup("field")); + + Info<< "Reading field " << fieldName << "\n" << endl; + volScalarField alpha1 + ( + IOobject + ( + fieldName, + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + scalar f0 = 0.0; + scalarField f(mesh.points().size()); + + Info<< "Processing type '" << shapeSelector::shapeTypeNames[surfType] + << "'" << endl; + + switch (surfType) + { + case shapeSelector::shapeType::PLANE: + { + const vector direction(dict.lookup("direction")); + + f = -(mesh.points() - centre) & (direction/mag(direction)); + f0 = 0.0; + break; + } + case shapeSelector::shapeType::SPHERE: + { + const scalar radius(readScalar(dict.lookup("radius"))); + + f = -mag(mesh.points() - centre); + f0 = -radius; + break; + } + case shapeSelector::shapeType::CYLINDER: + { + const scalar radius(readScalar(dict.lookup("radius"))); + const vector direction(dict.lookup("direction")); + + f = -sqrt + ( + sqr(mag(mesh.points() - centre)) + - sqr(mag((mesh.points() - centre) & direction)) + ); + f0 = -radius; + break; + } + case shapeSelector::shapeType::SIN: + { + const scalar period(readScalar(dict.lookup("period"))); + const scalar amplitude(readScalar(dict.lookup("amplitude"))); + const vector up(dict.lookup("up")); + const vector direction(dict.lookup("direction")); + + const scalarField xx + ( + (mesh.points() - centre) & direction/mag(direction) + ); + const scalarField zz((mesh.points() - centre) & up/mag(up)); + + f = amplitude*Foam::sin(2*mathematical::pi*xx/period) - zz; + f0 = 0; + break; + } + } + + + // Define function on mesh points and isovalue + + // Calculating alpha1 volScalarField from f = f0 isosurface + isoCutCell icc(mesh, f); + icc.volumeOfFluid(alpha1, f0); + + // Writing volScalarField alpha1 + ISstream::defaultPrecision(18); + alpha1.write(); + + Info<< nl << "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(" << alpha1.name() << ") = " << min(alpha1).value() + << " Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1 + << nl << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict b/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict new file mode 100644 index 0000000000..6bbb616b54 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict @@ -0,0 +1,24 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +field alpha.water; +type cylinder; +radius 0.25; +direction (0 1 0); +centre (0.5 0 0.5); + +// ************************************************************************* // diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C index 5f19867486..1d8e807045 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -234,25 +234,16 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate { const label celli = patch.faceCells()[facei]; - const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei]; - const scalar w = cornerWeights[facei]; - if (yPlus > yPlusLam_) - { - epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]); + epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]); - G0[celli] += - w - *(nutw[facei] + nuw[facei]) - *magGradUw[facei] - *Cmu25*sqrt(k[celli]) - /(kappa_*y[facei]); - } - else - { - epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]); - } + G0[celli] += + w + *(nutw[facei] + nuw[facei]) + *magGradUw[facei] + *Cmu25*sqrt(k[celli]) + /(kappa_*y[facei]); } } diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 5bed4060c5..59145dc988 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -238,8 +238,6 @@ void omegaWallFunctionFvPatchScalarField::calculate { const label celli = patch.faceCells()[facei]; - const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei]; - const scalar w = cornerWeights[facei]; const scalar omegaVis = 6*nuw[facei]/(beta1_*sqr(y[facei])); @@ -257,27 +255,17 @@ void omegaWallFunctionFvPatchScalarField::calculate omega0[celli] += w*sqrt(sqr(omegaVis) + sqr(omegaLog)); } - if (yPlus > yPlusLam_) + if (!blended_) { - if (!blended_) - { - omega0[celli] += w*omegaLog; - } + omega0[celli] += w*omegaLog; + } - G0[celli] += - w - *(nutw[facei] + nuw[facei]) - *magGradUw[facei] - *Cmu25*sqrt(k[celli]) - /(kappa_*y[facei]); - } - else - { - if (!blended_) - { - omega0[celli] += w*omegaVis; - } - } + G0[celli] += + w + *(nutw[facei] + nuw[facei]) + *magGradUw[facei] + *Cmu25*sqrt(k[celli]) + /(kappa_*y[facei]); } } diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 1087c3e378..ca9bf31ff6 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -238,6 +238,9 @@ fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C fvMatrices/solvers/MULES/MULES.C fvMatrices/solvers/MULES/CMULES.C +fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C +fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C +fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C interpolation = interpolation/interpolation diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C new file mode 100644 index 0000000000..33f63c2b6c --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C @@ -0,0 +1,1165 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016-2017 DHI +------------------------------------------------------------------------------- + +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 . + +\*---------------------------------------------------------------------------*/ + +#include "isoAdvection.H" +#include "volFields.H" +#include "interpolationCellPoint.H" +#include "volPointInterpolation.H" +#include "fvcSurfaceIntegrate.H" +#include "fvcGrad.H" +#include "upwind.H" +#include "cellSet.H" +#include "meshTools.H" +#include "OBJstream.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(isoAdvection, 0); +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::isoAdvection::isoAdvection +( + volScalarField& alpha1, + const surfaceScalarField& phi, + const volVectorField& U +) +: + // General data + mesh_(alpha1.mesh()), + dict_(mesh_.solverDict(alpha1.name())), + alpha1_(alpha1), + alpha1In_(alpha1.ref()), + phi_(phi), + U_(U), + dVf_ + ( + IOobject + ( + "dVf_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimVol, 0) + ), + advectionTime_(0), + + // Interpolation data + ap_(mesh_.nPoints()), + + // Tolerances and solution controls + nAlphaBounds_(dict_.lookupOrDefault