diff --git a/README.md b/README.md index aec529e09f..fed8a6e758 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ OpenCFD Ltd grants use of the OpenFOAM trademark by Third Parties on a licence b Please [contact OpenCFD](http://www.openfoam.com/contact) if you have any questions on the use of the OpenFOAM trademark. -Violations of the Trademark are continuously monitored, and will be duly prosecuted. +Violations of the Trademark are continuously monitored, and will be duly prosecuted. # Useful Links diff --git a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/files b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/files new file mode 100644 index 0000000000..0c6257a79d --- /dev/null +++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/files @@ -0,0 +1,3 @@ +laplacianDyMFoam.C + +EXE = $(FOAM_APPBIN)/overLaplacianDyMFoam diff --git a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options new file mode 100644 index 0000000000..149f9694b6 --- /dev/null +++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options @@ -0,0 +1,8 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -loverset diff --git a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H new file mode 100644 index 0000000000..d14935caab --- /dev/null +++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H @@ -0,0 +1,53 @@ + Info<< "Reading field T\n" << endl; + + volScalarField T + ( + IOobject + ( + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // Add overset specific interpolations + { + dictionary oversetDict; + oversetDict.add("T", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); + } + + + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + + Info<< "Reading diffusivity DT\n" << endl; + + dimensionedScalar DT + ( + transportProperties.lookup("DT") + ); diff --git a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/laplacianDyMFoam.C b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/laplacianDyMFoam.C new file mode 100644 index 0000000000..7bcd6b6cb6 --- /dev/null +++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/laplacianDyMFoam.C @@ -0,0 +1,110 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + laplacianFoam + +Group + grpBasicSolvers + +Description + Laplace equation solver for a scalar quantity. + + \heading Solver details + The solver is applicable to, e.g. for thermal diffusion in a solid. The + equation is given by: + + \f[ + \ddt{T} = \div \left( D_T \grad T \right) + \f] + + Where: + \vartable + T | Scalar field which is solved for, e.g. temperature + D_T | Diffusion coefficient + \endvartable + + \heading Required fields + \plaintable + T | Scalar field which is solved for, e.g. temperature + \endplaintable + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "fvOptions.H" +#include "simpleControl.H" +#include "dynamicFvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createNamedDynamicFvMesh.H" + + simpleControl simple(mesh); + + #include "createFields.H" + #include "createFvOptions.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nCalculating temperature distribution\n" << endl; + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + mesh.update(); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix TEqn + ( + fvm::ddt(T) - fvm::laplacian(DT, T) + == + fvOptions(T) + ); + + fvOptions.constrain(TEqn); + TEqn.solve(); + fvOptions.correct(T); + } + + #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/basic/laplacianFoam/overLaplacianDyMFoam/write.H b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/write.H new file mode 100644 index 0000000000..47aa182c0a --- /dev/null +++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/write.H @@ -0,0 +1,46 @@ + if (runTime.outputTime()) + { + volVectorField gradT(fvc::grad(T)); + + volScalarField gradTx + ( + IOobject + ( + "gradTx", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::X) + ); + + volScalarField gradTy + ( + IOobject + ( + "gradTy", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::Y) + ); + + volScalarField gradTz + ( + IOobject + ( + "gradTz", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::Z) + ); + + + runTime.write(); + } diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/files b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/files new file mode 100644 index 0000000000..3d96a548d0 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/files @@ -0,0 +1,3 @@ +rhoPimpleDyMFoam.C + +EXE = $(FOAM_APPBIN)/overRhoPimpleDyMFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options new file mode 100644 index 0000000000..ca5b141675 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options @@ -0,0 +1,25 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +EXE_LIBS = \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -loverset \ + -lfvOptions \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H new file mode 100644 index 0000000000..37072312ff --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H @@ -0,0 +1,11 @@ +CorrectPhi +( + U, + phi, + p, + rho, + psi, + dimensionedScalar("rAUf", dimTime, 1), + divrhoU, + pimple +); diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H new file mode 100644 index 0000000000..aed0e76956 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H @@ -0,0 +1,11 @@ +#include "createTimeControls.H" + +bool correctPhi +( + pimple.dict().lookupOrDefault("correctPhi", true) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H new file mode 100644 index 0000000000..27f568bce7 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H @@ -0,0 +1,117 @@ +Info<< "Reading thermophysical properties\n" << endl; + +autoPtr pThermo +( + psiThermo::New(mesh) +); +psiThermo& thermo = pThermo(); +thermo.validate(args.executable(), "h", "e"); + +volScalarField& p = thermo.p(); +const volScalarField& psi = thermo.psi(); + +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "compressibleCreatePhi.H" + +dimensionedScalar rhoMax +( + dimensionedScalar::lookupOrDefault + ( + "rhoMax", + pimple.dict(), + dimDensity, + GREAT + ) +); + +dimensionedScalar rhoMin +( + dimensionedScalar::lookupOrDefault + ( + "rhoMin", + pimple.dict(), + dimDensity, + 0 + ) +); + +Info<< "Creating turbulence model\n" << endl; +autoPtr turbulence +( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) +); + +mesh.setFluxRequired(p.name()); + +Info<< "Creating field dpdt\n" << endl; +volScalarField dpdt +( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("dpdt", p.dimensions()/dimTime, 0) +); + +Info<< "Creating field kinetic energy K\n" << endl; +volScalarField K("K", 0.5*magSqr(U)); + + +//- Overset specific + +// Add solver-specific interpolations +{ + dictionary oversetDict; + oversetDict.add("U", true); + oversetDict.add("p", true); + oversetDict.add("HbyA", true); + oversetDict.add("grad(p)", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); +} + +// Mask field for zeroing out contributions on hole cells +#include "createCellMask.H" diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H new file mode 100644 index 0000000000..5f1fc7451c --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H @@ -0,0 +1,123 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU)); +volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p)); +//mesh.interpolate(HbyA); + +if (pimple.nCorrPISO() <= 1) +{ + tUEqn.clear(); +} + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + fvc::flux(HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho) + ) + ); + + fvc::makeRelative(phid, psi, U); + MRF.makeRelative(fvc::interpolate(psi), phid); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi == pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::flux(rho*HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, rhoUf) + ); + + fvc::makeRelative(phiHbyA, rho, U); + MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + + while (pimple.correctNonOrthogonal()) + { + // Pressure corrector + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); +Info<< "rho max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + +volVectorField gradP(fvc::grad(p)); +//mesh.interpolate(gradP); +U = HbyA - rAU*cellMask*gradP; +U.correctBoundaryConditions(); +fvOptions.correct(U); +K = 0.5*magSqr(U); + +{ + rhoUf = fvc::interpolate(rho*U); + surfaceVectorField n(mesh.Sf()/mesh.magSf()); + rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf)); +} + +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); + + if (mesh.moving()) + { + dpdt -= fvc::div(fvc::meshPhi(rho, U), p); + } +} diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H new file mode 100644 index 0000000000..ed2db49fb4 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H @@ -0,0 +1,6 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", true); + +checkMeshCourantNo = + pimple.dict().lookupOrDefault("checkMeshCourantNo", false); diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/rhoPimpleDyMFoam.C new file mode 100644 index 0000000000..0bbcefa14b --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/rhoPimpleDyMFoam.C @@ -0,0 +1,173 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + rhoPimpleFoam + +Group + grpCompressibleSolvers grpMovingMeshSolvers + +Description + Transient solver for laminar or turbulent flow of compressible fluids + for HVAC and similar applications. + + Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and + pseudo-transient simulations. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "psiThermo.H" +#include "turbulentFluidThermoModel.H" +#include "bound.H" +#include "pimpleControl.H" +#include "CorrectPhi.H" +#include "fvOptions.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" +#include "cellCellStencilObject.H" +#include "localMin.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + + pimpleControl pimple(mesh); + + #include "createRDeltaT.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createMRF.H" + #include "createFvOptions.H" + #include "createRhoUf.H" + #include "createControls.H" + + turbulence->validate(); + + if (!LTS) + { + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + } + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readControls.H" + + { + // Store divrhoU from the previous mesh so that it can be mapped + // and used in correctPhi to ensure the corrected phi has the + // same divergence + volScalarField divrhoU + ( + "divrhoU", + fvc::div(fvc::absolute(phi, rho, U)) + ); + + if (LTS) + { + #include "setRDeltaT.H" + } + else + { + #include "compressibleCourantNo.H" + #include "setDeltaT.H" + } + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // Store momentum to set rhoUf for introduced faces. + volVectorField rhoU("rhoU", rho*U); + + // Do any mesh changes + mesh.update(); + + if (mesh.changing()) + { + #include "setCellMask.H" + } + + if (mesh.changing() && correctPhi) + { + // Calculate absolute flux from the mapped surface velocity + phi = mesh.Sf() & rhoUf; + + #include "correctPhi.H" + + // Make the fluxes relative to the mesh-motion + fvc::makeRelative(phi, rho, U); + } + } + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + + #include "rhoEqn.H" + Info<< "rhoEqn max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "EEqn.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/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H new file mode 100644 index 0000000000..6353e7908f --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + surfaceScalarField phiMask(localMin(mesh).interpolate(cellMask)); + + scalarField sumPhi(fvc::surfaceSum(mag(phiMask*phi))().internalField()); + + CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/files b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/files new file mode 100644 index 0000000000..7c557aba24 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/files @@ -0,0 +1,3 @@ +pimpleDyMFoam.C + +EXE = $(FOAM_APPBIN)/overPimpleDyMFoam diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/options b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/options new file mode 100644 index 0000000000..c538a99b57 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/Make/options @@ -0,0 +1,24 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +EXE_LIBS = \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lfvOptions \ + -lsampling \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -loverset diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H new file mode 100644 index 0000000000..aef0116f7b --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H @@ -0,0 +1,24 @@ +// Solve the Momentum equation + +MRF.correctBoundaryVelocity(U); + +tmp tUEqn +( + fvm::ddt(U) + fvm::div(phi, U) + + MRF.DDt(U) + + turbulence->divDevReff(U) + == + fvOptions(U) +); +fvVectorMatrix& UEqn = tUEqn.ref(); + +UEqn.relax(); + +fvOptions.constrain(UEqn); + +if (pimple.momentumPredictor()) +{ + solve(UEqn == -cellMask*fvc::grad(p)); + + fvOptions.correct(U); +} diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/continuityErrs.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/continuityErrs.H new file mode 100644 index 0000000000..c787055623 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/continuityErrs.H @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + continuityErrs + +Description + Calculates and prints the continuity errors. + +\*---------------------------------------------------------------------------*/ + +{ + volScalarField contErr(interpolatedCells*cellMask*fvc::div(phi)); + + scalar sumLocalContErr = runTime.deltaTValue()* + mag(contErr)().weightedAverage(mesh.V()).value(); + + scalar globalContErr = runTime.deltaTValue()* + contErr.weightedAverage(mesh.V()).value(); + cumulativeContErr += globalContErr; + + Info<< "time step continuity errors : sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr + << endl; +} + +// ************************************************************************* // diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/correctPhi.H new file mode 100644 index 0000000000..1951316f46 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/correctPhi.H @@ -0,0 +1,113 @@ +if (mesh.changing()) +{ + volVectorField::Boundary& bfld = U.boundaryFieldRef(); + forAll(bfld, patchi) + { + if (bfld[patchi].fixesValue()) + { + bfld[patchi].initEvaluate(); + } + } + + surfaceScalarField::Boundary& phiBfld = phi.boundaryFieldRef(); + forAll(bfld, patchi) + { + if (bfld[patchi].fixesValue()) + { + bfld[patchi].evaluate(); + + phiBfld[patchi] = bfld[patchi] & mesh.Sf().boundaryField()[patchi]; + } + } +} + // Initialize BCs list for pcorr to zero-gradient + wordList pcorrTypes + ( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + // Set BCs of pcorr to fixed-value for patches at which p is fixed + forAll(p.boundaryField(), patchi) + { + if (p.boundaryField()[patchi].fixesValue()) + { + pcorrTypes[patchi] = fixedValueFvPatchScalarField::typeName; + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0.0), + pcorrTypes + ); + +{ + dimensionedScalar rAUf("rAUf", dimTime, 1.0); + + const cellCellStencilObject& overlap = Stencil::New(mesh); + const labelList& cellTypes = overlap.cellTypes(); + const labelIOList& zoneIDs = overlap.zoneID(); + + while (pimple.correctNonOrthogonal()) + { + label nZones = gMax(zoneIDs)+1; + + //label refCellI2 = -1; + labelList refCells(nZones, -1); + labelList refZones(nZones, -1); + + forAll(zoneIDs, cellI) + { + label zoneId = zoneIDs[cellI]; + if + ( + refCells[zoneId] == -1 + && cellTypes[cellI] == cellCellStencil::CALCULATED + && refZones[zoneId] == -1 + ) + { + refCells[zoneId] = cellI; + refZones[zoneId] = zoneId; + } + } + + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(rAUf, pcorr) == fvc::div(phi) + ); + + //pcorrEqn.setReference(refCellI2, 0.0, true); + scalarList values(nZones, 0.0); + pcorrEqn.setReferences(refCells, values, true); + + const dictionary& d = mesh.solver + ( + pcorr.select + ( + pimple.finalInnerIter() + ) + ); + mesh.fvMesh::solve(pcorrEqn, d); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + } + + if (runTime.outputTime()) + { + volScalarField("contPhiPcorr", fvc::div(phi)).write(); + pcorr.write(); + } +} diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H new file mode 100644 index 0000000000..53ed2009c7 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H @@ -0,0 +1,26 @@ +#include "createTimeControls.H" + +bool correctPhi +( + pimple.dict().lookupOrDefault("correctPhi", false) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); + +bool massFluxInterpolation +( + pimple.dict().lookupOrDefault("massFluxInterpolation", false) +); + +bool adjustFringe +( + pimple.dict().lookupOrDefault("oversetAdjustPhi", false) +); + +bool ddtCorr +( + pimple.dict().lookupOrDefault("ddtCorr", true) +); diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H new file mode 100644 index 0000000000..ba991bfd48 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H @@ -0,0 +1,70 @@ +Info<< "Reading field p\n" << endl; +volScalarField p +( + IOobject + ( + "p", + 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" + + +label pRefCell = 0; +scalar pRefValue = 0.0; +setRefCell(p, pimple.dict(), pRefCell, pRefValue); +mesh.setFluxRequired(p.name()); + + +//- Overset specific + +// Add solver-specific interpolations +{ + dictionary oversetDict; + oversetDict.add("U", true); + oversetDict.add("p", true); + oversetDict.add("HbyA", true); + oversetDict.add("grad(p)", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); +} + +// Mask field for zeroing out contributions on hole cells +#include "createCellMask.H" + +// Create bool field with interpolated cells +#include "createInterpolatedCells.H" + +singlePhaseTransportModel laminarTransport(U, phi); + +autoPtr turbulence +( + incompressible::turbulenceModel::New(U, phi, laminarTransport) +); diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H new file mode 100644 index 0000000000..e13a2cc0a7 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H @@ -0,0 +1,269 @@ +// Interpolation used +interpolationCellPoint UInterpolator(HbyA); + +// Determine faces on outside of interpolated cells +PackedBoolList isOwnerInterpolatedFace(mesh.nInternalFaces()); +PackedBoolList isNeiInterpolatedFace(mesh.nInternalFaces()); + +// Determine donor cells +labelListList donorCell(mesh.nInternalFaces()); + +scalarListList weightCellCells(mesh.nInternalFaces()); + +// Interpolated HbyA faces +vectorField UIntFaces(mesh.nInternalFaces(), vector::zero); + +// Determine receptor neighbourd cells +labelList receptorNeigCell(mesh.nInternalFaces(), -1); + +{ + const cellCellStencilObject& overlap = Stencil::New(mesh); + const labelList& cellTypes = overlap.cellTypes(); + const labelIOList& zoneID = overlap.zoneID(); + + label nZones = gMax(zoneID)+1; + PtrList meshParts(nZones); + labelList nCellsPerZone(nZones, 0); + + forAll(nCellsPerZone, zoneI) + { + meshParts.set(zoneI, new fvMeshSubset(mesh)); + meshParts[zoneI].setLargeCellSubset(zoneID, zoneI); + } + + for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + { + label ownType = cellTypes[mesh.faceOwner()[faceI]]; + label neiType = cellTypes[mesh.faceNeighbour()[faceI]]; + if + ( + ownType == cellCellStencil::INTERPOLATED + && neiType == cellCellStencil::CALCULATED + ) + { + isOwnerInterpolatedFace[faceI] = true; + + const vector& fc = mesh.faceCentres()[faceI]; + + for (label zoneI = 0; zoneI < nZones; zoneI++) + { + if (zoneI != zoneID[mesh.faceOwner()[faceI]]) + { + const fvMesh& partMesh = meshParts[zoneI].subMesh(); + const labelList& cellMap = meshParts[zoneI].cellMap(); + label cellI = partMesh.findCell(fc); + + if (cellI != -1) + { + // Determine weights + labelList stencil(partMesh.cellCells()[cellI]); + + stencil.append(cellI); + + label st = stencil.size(); + + donorCell[faceI].setSize(st); + + weightCellCells[faceI].setSize(st); + + scalarField weights(st); + + forAll(stencil, i) + { + scalar d = mag + ( + partMesh.cellCentres()[stencil[i]] + - fc + ); + weights[i] = 1.0/d; + donorCell[faceI][i] = cellMap[stencil[i]]; + } + weights /= sum(weights); + + weightCellCells[faceI] = weights; + + forAll(stencil, i) + { + UIntFaces[faceI] += + weightCellCells[faceI][i] + *UInterpolator.interpolate + ( + fc, + donorCell[faceI][i] + ); + } + + break; + } + } + } + + receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI]; + } + else if + ( + ownType == cellCellStencil::CALCULATED + && neiType == cellCellStencil::INTERPOLATED + ) + { + isNeiInterpolatedFace[faceI] = true; + + const vector& fc = mesh.faceCentres()[faceI]; + for (label zoneI = 0; zoneI < nZones; zoneI++) + { + if (zoneI != zoneID[mesh.faceNeighbour()[faceI]]) + { + const fvMesh& partMesh = meshParts[zoneI].subMesh(); + const labelList& cellMap = meshParts[zoneI].cellMap(); + label cellI = partMesh.findCell(fc); + + if (cellI != -1) + { + // Determine weights + labelList stencil(partMesh.cellCells()[cellI]); + + stencil.append(cellI); + + label st = stencil.size(); + + donorCell[faceI].setSize(st); + + weightCellCells[faceI].setSize(st); + + scalarField weights(st); + + forAll(stencil, i) + { + scalar d = mag + ( + partMesh.cellCentres()[stencil[i]] + - fc + ); + weights[i] = 1.0/d; + donorCell[faceI][i] = cellMap[stencil[i]]; + } + weights /= sum(weights); + + weightCellCells[faceI] = weights; + + forAll(stencil, i) + { + UIntFaces[faceI] += + weightCellCells[faceI][i] + *UInterpolator.interpolate + ( + fc, + donorCell[faceI][i] + ); + } + + break; + } + } + } + + receptorNeigCell[faceI] = mesh.faceOwner()[faceI]; + } + } +} + +// contravariant U +vectorField U1Contrav(mesh.nInternalFaces(), vector::zero); + +surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf()); + +forAll(isNeiInterpolatedFace, faceI) +{ + label cellId = -1; + if (isNeiInterpolatedFace[faceI]) + { + cellId = mesh.faceNeighbour()[faceI]; + } + else if (isOwnerInterpolatedFace[faceI]) + { + cellId = mesh.faceOwner()[faceI]; + } + + if (cellId != -1) + { + const vector& n = faceNormals[faceI]; + vector n1 = vector::zero; + // 2-D cases + if (mesh.nSolutionD() == 2) + { + for (direction cmpt=0; cmpt mag(n.y()) && mag(n.x()) > mag(n.z())) + { + n1 = vector(n.y(), -n.x(), 0); + } + else if (mag(n.y()) > mag(n.z())) + { + n1 = vector(0, n.z(), -n.y()); + } + else + { + n1 = vector(-n.z(), 0, n.x()); + } + } + + n1 /= mag(n1); + + vector n2 = n ^ n1; + n2 /= mag(n2); + + tensor rot = + tensor + ( + n.x() ,n.y(), n.z(), + n1.x() ,n1.y(), n1.z(), + n2.x() ,n2.y(), n2.z() + ); + +// tensor rot = +// tensor +// ( +// n & x ,n & y, n & z, +// n1 & x ,n1 & y, n1 & z, +// n2 & x ,n2 & y, n2 & z +// ); + + U1Contrav[faceI].x() = + 2*transform(rot, UIntFaces[faceI]).x() + - transform(rot, HbyA[receptorNeigCell[faceI]]).x(); + + U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y(); + + U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z(); + + HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]); + } +} diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H new file mode 100644 index 0000000000..028dfe4306 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H @@ -0,0 +1,109 @@ +volScalarField rAU(1.0/UEqn.A()); + +// Option 1: interpolate rAU, do not block out rAU on blocked cells +//mesh.interpolate(rAU, false); +//surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + +// Option 2: do not interpolate rAU but block out rAU +//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU)); + + +// Option 3: do not interpolate rAU but zero out rAUf on faces on holes +// But what about: +// +// H +// H I C C C C +// H +// + +surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); +surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU)); + +volVectorField HbyA("HbyA", U); +HbyA = constrainHbyA(rAU*UEqn.H(), U, p); + +//mesh.interpolate(HbyA); +if (massFluxInterpolation) +{ + #include "interpolatedFaces.H" +} + +if (pimple.nCorrPISO() <= 1) +{ + tUEqn.clear(); +} + +surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); + +if (ddtCorr) +{ + phiHbyA += rAUf*fvc::ddtCorr(U, Uf); +} + +MRF.makeRelative(phiHbyA); + +if (p.needReference()) +{ + fvc::makeRelative(phiHbyA, U); + adjustPhi(phiHbyA, U, p); + fvc::makeAbsolute(phiHbyA, U); +} + +if (adjustFringe) +{ + fvc::makeRelative(phiHbyA, U); + oversetAdjustPhi(phiHbyA, U); + fvc::makeAbsolute(phiHbyA, U); +} + +if (runTime.outputTime()) +{ + volScalarField + ( + "div(phiHbyA)", + fvc::div(phiHbyA) + //interpolatedCells*cellMask*fvc::div(phiHbyA) + ).write(); +} + +while (pimple.correctNonOrthogonal()) +{ + fvScalarMatrix pEqn + ( + fvm::laplacian(rAUf, p) == fvc::div(phiHbyA) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA - pEqn.flux(); + } +} + +#include "continuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); +volVectorField gradP(fvc::grad(p)); +//mesh.interpolate(gradP); + +// Option 1: leave velocity intact on blocked out cells +//U = HbyA - rAU*gradP; + +// Option 2: zero out velocity on blocked out cells +U = (HbyA - rAU*cellMask*gradP); +U.correctBoundaryConditions(); + +fvOptions.correct(U); + +{ + Uf = fvc::interpolate(U); + surfaceVectorField n(mesh.Sf()/mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); +} + +// Make the fluxes relative to the mesh motion +fvc::makeRelative(phi, U); diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pimpleDyMFoam.C new file mode 100644 index 0000000000..f3dbd6a989 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pimpleDyMFoam.C @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + pimpleDyMFoam.C + +Group + grpIncompressibleSolvers grpMovingMeshSolvers + +Description + Transient solver for incompressible, flow of Newtonian fluids + on a moving mesh using the PIMPLE (merged PISO-SIMPLE) algorithm. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" + +#include "cellCellStencilObject.H" +#include "zeroGradientFvPatchFields.H" +#include "localMin.H" +#include "interpolationCellPoint.H" +#include "transform.H" +#include "fvMeshSubset.H" +#include "oversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::addNote + ( + "Experimental version of pimpleDyMFoam with support for overset meshes" + ); + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "initContinuityErrs.H" + + pimpleControl pimple(mesh); + + #include "createFields.H" + #include "createUf.H" + #include "createMRF.H" + #include "createFvOptions.H" + #include "createControls.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readControls.H" + #include "CourantNo.H" + + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool changed = mesh.update(); + + if (changed) + { + #include "setCellMask.H" + #include "setInterpolatedCells.H" + } + + // Calculate absolute flux from the mapped surface velocity + phi = mesh.Sf() & Uf; + + if (runTime.outputTime()) + { + volScalarField + ( + "contPhi", + interpolatedCells*cellMask*fvc::div(phi) + ).write(); + } + + if (mesh.changing() && correctPhi) + { + #include "correctPhi.H" + } + + // Make the flux relative to the mesh motion + fvc::makeRelative(phi, U); + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + laminarTransport.correct(); + 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/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H new file mode 100644 index 0000000000..0141c85225 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H @@ -0,0 +1,10 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", false); + +checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false); + +massFluxInterpolation = + pimple.dict().lookupOrDefault("massFluxInterpolation", false); + +ddtCorr = pimple.dict().lookupOrDefault("ddtCorr", true); diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/files b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/files new file mode 100644 index 0000000000..3294949edd --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/files @@ -0,0 +1,3 @@ +simpleFoam.C + +EXE = $(FOAM_APPBIN)/overSimpleFoam diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/options b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/options new file mode 100644 index 0000000000..639b94132a --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/Make/options @@ -0,0 +1,25 @@ +EXE_INC = \ + -I. \ + -I$(FOAM_SOLVERS)/incompressible/pimpleFoam/overPimpleDyMFoam \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +EXE_LIBS = \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lfvOptions \ + -lsampling \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -loverset diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H new file mode 100644 index 0000000000..5e9861ccc3 --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H @@ -0,0 +1,21 @@ + // Momentum predictor + + MRF.correctBoundaryVelocity(U); + + tmp tUEqn + ( + fvm::div(phi, U) + + MRF.DDt(U) + + turbulence->divDevReff(U) + == + fvOptions(U) + ); + fvVectorMatrix& UEqn = tUEqn.ref(); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + solve(UEqn == -cellMask*fvc::grad(p)); + + fvOptions.correct(U); diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/continuityErrs.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/continuityErrs.H new file mode 100644 index 0000000000..c787055623 --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/continuityErrs.H @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + continuityErrs + +Description + Calculates and prints the continuity errors. + +\*---------------------------------------------------------------------------*/ + +{ + volScalarField contErr(interpolatedCells*cellMask*fvc::div(phi)); + + scalar sumLocalContErr = runTime.deltaTValue()* + mag(contErr)().weightedAverage(mesh.V()).value(); + + scalar globalContErr = runTime.deltaTValue()* + contErr.weightedAverage(mesh.V()).value(); + cumulativeContErr += globalContErr; + + Info<< "time step continuity errors : sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr + << endl; +} + +// ************************************************************************* // diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H new file mode 100644 index 0000000000..b3d40776f4 --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createFields.H @@ -0,0 +1,47 @@ +Info<< "Reading field p\n" << endl; +volScalarField p +( + IOobject + ( + "p", + 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" + + +label pRefCell = 0; +scalar pRefValue = 0.0; +setRefCell(p, simple.dict(), pRefCell, pRefValue); +mesh.setFluxRequired(p.name()); + + +singlePhaseTransportModel laminarTransport(U, phi); + +autoPtr turbulence +( + incompressible::turbulenceModel::New(U, phi, laminarTransport) +); + +#include "createMRF.H" + +#include "createOversetFields.H" diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H new file mode 100644 index 0000000000..7eda905f50 --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H @@ -0,0 +1,34 @@ +//- Overset specific + +// Add solver-specific interpolations +{ + dictionary oversetDict; + oversetDict.add("U", true); + oversetDict.add("p", true); + oversetDict.add("HbyA", true); + oversetDict.add("grad(p)", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); +} + +// Mask field for zeroing out contributions on hole cells +#include "createCellMask.H" + +#include "createInterpolatedCells.H" + +bool adjustFringe +( + simple.dict().lookupOrDefault("oversetAdjustPhi", false) +); +bool massFluxInterpolation +( + simple.dict().lookupOrDefault("massFluxInterpolation", false) +); diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H new file mode 100644 index 0000000000..db58f90828 --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/createUpdatedDynamicFvMesh.H @@ -0,0 +1,22 @@ + Info<< "Create dynamic mesh for time = " + << runTime.timeName() << nl << endl; + + autoPtr meshPtr + ( + dynamicFvMesh::New + ( + IOobject + ( + dynamicFvMesh::defaultRegion, + runTime.timeName(), + runTime, + IOobject::MUST_READ + ) + ) + ); + + dynamicFvMesh& mesh = meshPtr(); + + // Calculate initial mesh-to-mesh mapping. Note that this should be + // done under the hood, e.g. as a MeshObject + mesh.update(); diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H new file mode 100644 index 0000000000..c64f92e33f --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H @@ -0,0 +1,57 @@ +{ + surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); + + volScalarField rAU(1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = constrainHbyA(rAU*UEqn.H(), U, p); + + //mesh.interpolate(HbyA); + if (massFluxInterpolation) + { + #include "interpolatedFaces.H" + } + + tUEqn.clear(); + + surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); + + adjustPhi(phiHbyA, U, p); + + if (adjustFringe) + { + oversetAdjustPhi(phiHbyA, U); + } + + // Non-orthogonal pressure corrector loop + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rAUf, p) == fvc::div(phiHbyA) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA - pEqn.flux(); + } + } + + #include "continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + // Momentum corrector + volVectorField gradP(fvc::grad(p)); + //mesh.interpolate(gradP); + + U = HbyA - rAU*cellMask*gradP; + U.correctBoundaryConditions(); + fvOptions.correct(U); +} diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/simpleFoam.C new file mode 100644 index 0000000000..4a7cf1315b --- /dev/null +++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/simpleFoam.C @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + simpleFoam + +Group + grpIncompressibleSolvers + +Description + Steady-state solver for incompressible flows with turbulence modelling. + + \heading Solver details + The solver uses the SIMPLE algorithm to solve the continuity equation: + + \f[ + \div \vec{U} = 0 + \f] + + and momentum equation: + + \f[ + \div \left( \vec{U} \vec{U} \right) - \div \gvec{R} + = - \grad p + \vec{S}_U + \f] + + Where: + \vartable + \vec{U} | Velocity + p | Pressure + \vec{R} | Stress tensor + \vec{S}_U | Momentum source + \endvartable + + \heading Required fields + \plaintable + U | Velocity [m/s] + p | Kinematic pressure, p/rho [m2/s2] + \ | As required by user selection + \endplaintable + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "simpleControl.H" +#include "fvOptions.H" + +#include "dynamicFvMesh.H" +#include "cellCellStencilObject.H" +#include "localMin.H" +#include "interpolationCellPoint.H" +#include "fvMeshSubset.H" +#include "transform.H" +#include "oversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #define CREATE_MESH createUpdatedDynamicFvMesh.H + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createUpdatedDynamicFvMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity SIMPLE corrector + { + #include "UEqn.H" + #include "pEqn.H" + } + + laminarTransport.correct(); + 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/interFoam/overInterDyMFoam/CourantNo.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.H new file mode 100644 index 0000000000..4182f4e209 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.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 + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + surfaceScalarField phiMask(localMin(mesh).interpolate(cellMask)); + + scalarField sumPhi + ( + fvc::surfaceSum(mag(phiMask*phi))().internalField() + ); + + CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/files b/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/files new file mode 100644 index 0000000000..e229fbaed0 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/files @@ -0,0 +1,3 @@ +interDyMFoam.C + +EXE = $(FOAM_APPBIN)/overInterDyMFoam diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/options b/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/options new file mode 100644 index 0000000000..3b8ac67af5 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/Make/options @@ -0,0 +1,31 @@ +EXE_INC = \ + -I. \ + -I.. \ + -I../../VoF \ + -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)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(FOAM_SOLVERS)/incompressible/pimpleFoam/overPimpleDyMFoam \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -limmiscibleIncompressibleTwoPhaseMixture \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lfiniteVolume \ + -ldynamicMesh \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -loverset \ + -lfvOptions \ + -lsampling \ + -lwaveModels diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H new file mode 100644 index 0000000000..600c94829c --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/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 + == + cellMask*fvc::reconstruct + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + fvOptions.correct(U); + } diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/continuityErrs.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/continuityErrs.H new file mode 100644 index 0000000000..c787055623 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/continuityErrs.H @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + continuityErrs + +Description + Calculates and prints the continuity errors. + +\*---------------------------------------------------------------------------*/ + +{ + volScalarField contErr(interpolatedCells*cellMask*fvc::div(phi)); + + scalar sumLocalContErr = runTime.deltaTValue()* + mag(contErr)().weightedAverage(mesh.V()).value(); + + scalar globalContErr = runTime.deltaTValue()* + contErr.weightedAverage(mesh.V()).value(); + cumulativeContErr += globalContErr; + + Info<< "time step continuity errors : sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr + << endl; +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H new file mode 100644 index 0000000000..81948b54da --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/correctPhi.H @@ -0,0 +1,125 @@ +{ + if (mesh.changing()) + { + volVectorField::Boundary& bfld = U.boundaryFieldRef(); + forAll(bfld, patchi) + { + if (bfld[patchi].fixesValue()) + { + bfld[patchi].initEvaluate(); + } + } + + surfaceScalarField::Boundary& phiBfld = phi.boundaryFieldRef(); + forAll(bfld, patchi) + { + if (bfld[patchi].fixesValue()) + { + bfld[patchi].evaluate(); + + phiBfld[patchi] = + bfld[patchi] + & mesh.Sf().boundaryField()[patchi]; + } + } + } + + wordList pcorrTypes + ( + p_rgh.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + for (label i=0; i("correctPhi", true) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); + +bool moveMeshOuterCorrectors +( + pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false) +); + + +bool massFluxInterpolation +( + pimple.dict().lookupOrDefault("massFluxInterpolation", false) +); + +bool ddtCorr +( + pimple.dict().lookupOrDefault("ddtCorr", true) +); diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H new file mode 100644 index 0000000000..4a83ab5850 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H @@ -0,0 +1,171 @@ +#include "createRDeltaT.H" + +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" + + + + //- Overset specific + + // Add solver-specific interpolations + { + dictionary oversetDict; + oversetDict.add("U", true); + oversetDict.add("p", true); + oversetDict.add("HbyA", true); + oversetDict.add("p_rgh", true); + oversetDict.add("alpha1", true); + oversetDict.add("minGradP", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); + } + + // Mask field for zeroing out contributions on hole cells + #include "createCellMask.H" + + // Create bool field with interpolated cells + #include "createInterpolatedCells.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()); + +// MULES compressed flux is registered in case scalarTransport FO needs it. +surfaceScalarField alphaPhiUn +( + IOobject + ( + "alphaPhiUn", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", phi.dimensions(), 0.0) +); + +#include "createMRF.H" diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C new file mode 100644 index 0000000000..7bb7f8af6b --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C @@ -0,0 +1,209 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 + overInterDyMFoam + +Group + grpMultiphaseSolvers grpMovingMeshSolvers + +Description + Solver for 2 incompressible, isothermal immiscible fluids using a VOF + (volume of fluid) phase-fraction based interface capturing approach, + with optional mesh motion and mesh topology changes including adaptive + re-meshing. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "CMULES.H" +#include "EulerDdtScheme.H" +#include "localEulerDdtScheme.H" +#include "CrankNicolsonDdtScheme.H" +#include "subCycle.H" +#include "immiscibleIncompressibleTwoPhaseMixture.H" +#include "turbulentTransportModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "CorrectPhi.H" +#include "fvcSmooth.H" +#include "cellCellStencilObject.H" +#include "localMin.H" +#include "interpolationCellPoint.H" +#include "transform.H" +#include "fvMeshSubset.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "initContinuityErrs.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "createDyMControls.H" + #include "createFields.H" + #include "createAlphaFluxes.H" + #include "createFvOptions.H" + + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0) + ); + + #include "correctPhi.H" + #include "createUf.H" + + turbulence->validate(); + + if (!LTS) + { + #include "CourantNo.H" + #include "setInitialDeltaT.H" + } + + #include "setCellMask.H" + #include "setInterpolatedCells.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readControls.H" + + if (LTS) + { + #include "setRDeltaT.H" + } + else + { + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" + } + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + if (pimple.firstIter() || moveMeshOuterCorrectors) + { + scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); + + mesh.update(); + + if (mesh.changing()) + { + Info<< "Execution time for mesh.update() = " + << runTime.elapsedCpuTime() - timeBeforeMeshUpdate + << " s" << endl; + + // Do not apply previous time-step mesh compression flux + // if the mesh topology changed + if (mesh.topoChanging()) + { + talphaPhiCorr0.clear(); + } + + gh = (g & mesh.C()) - ghRef; + ghf = (g & mesh.Cf()) - ghRef; + + // Update cellMask field for blocking out hole cells + #include "setCellMask.H" + #include "setInterpolatedCells.H" + } + + if ((mesh.changing() && correctPhi) || mesh.topoChanging()) + { + // Calculate absolute flux from the mapped surface velocity + // Note: temporary fix until mapped Uf is assessed + Uf = fvc::interpolate(U); + + // Calculate absolute flux from the mapped surface velocity + phi = mesh.Sf() & Uf; + + #include "correctPhi.H" + + // Make the flux relative to the mesh motion + fvc::makeRelative(phi, U); + + mixture.correct(); + } + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + } + + #include "alphaControls.H" + #include "alphaEqnSubCycle.H" + + mixture.correct(); + + #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/interFoam/overInterDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H new file mode 100644 index 0000000000..0a5b0296b8 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H @@ -0,0 +1,96 @@ +{ + rAU = 1.0/UEqn.A(); + surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); + + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + //HbyA = rAU*UEqn.H(); + HbyA = constrainHbyA(rAU*UEqn.H(), U, p_rgh); + + if (massFluxInterpolation) + { + #include "interpolatedFaces.H" + } + + surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); + + if (ddtCorr) + { + phiHbyA += fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf); + } + MRF.makeRelative(phiHbyA); + + if (p_rgh.needReference()) + { + fvc::makeRelative(phiHbyA, U); + adjustPhi(phiHbyA, U, p_rgh); + fvc::makeAbsolute(phiHbyA, U); + } + + surfaceScalarField phig + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*faceMask*rAUf*mesh.magSf() + ); + + 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(faceMask*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(); + + // Reconstruct body forces (-grad(p) and gh etc) + volVectorField minGradP + ( + "minGradP", + fvc::reconstruct((phig - p_rghEqn.flux())/rAUf) + ); + U = HbyA + rAU*cellMask*minGradP; + U.correctBoundaryConditions(); + fvOptions.correct(U); + } + } + + #include "continuityErrs.H" + + { + Uf = fvc::interpolate(U); + surfaceVectorField n(mesh.Sf()/mesh.magSf()); + Uf += n*(phi/mesh.magSf() - (n & Uf)); + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + 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/interFoam/overInterDyMFoam/readControls.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H new file mode 100644 index 0000000000..0df2a74f30 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H @@ -0,0 +1,15 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", true); + +checkMeshCourantNo = + pimple.dict().lookupOrDefault("checkMeshCourantNo", false); + +moveMeshOuterCorrectors = + pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false); + +massFluxInterpolation = + pimple.dict().lookupOrDefault("massFluxInterpolation", false); + +ddtCorr = + pimple.dict().lookupOrDefault("ddtCorr", true); \ No newline at end of file diff --git a/applications/test/pTraits/Test-pTraits.C b/applications/test/pTraits/Test-pTraits.C index d300ddc87a..68d6a01527 100644 --- a/applications/test/pTraits/Test-pTraits.C +++ b/applications/test/pTraits/Test-pTraits.C @@ -92,14 +92,14 @@ int main() std::cout<< "max = " << pTraits