diff --git a/applications/solvers/overset/icoDyMOversetFoam/Make/files b/applications/solvers/overset/icoDyMOversetFoam/Make/files new file mode 100644 index 000000000..3543e4382 --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/Make/files @@ -0,0 +1,3 @@ +icoDyMOversetFoam.C + +EXE = $(FOAM_APPBIN)/icoDyMOversetFoam diff --git a/applications/solvers/overset/icoDyMOversetFoam/Make/options b/applications/solvers/overset/icoDyMOversetFoam/Make/options new file mode 100644 index 000000000..5668d862d --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/Make/options @@ -0,0 +1,14 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -ldynamicFvMesh \ + -ldynamicMesh \ + -loversetMesh \ + -llduSolvers diff --git a/applications/solvers/overset/icoDyMOversetFoam/UEqn.H b/applications/solvers/overset/icoDyMOversetFoam/UEqn.H new file mode 100644 index 000000000..e1e3824ce --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/UEqn.H @@ -0,0 +1,12 @@ +// Momentum equation +fvVectorMatrix UEqn +( + fvm::ddt(U) + + fvm::div(phi, U) + - fvm::laplacian(nu, U) +); + +if (piso.momentumPredictor()) +{ + solve(UEqn == -fvc::grad(p)); +} diff --git a/applications/solvers/overset/icoDyMOversetFoam/correctPhi.H b/applications/solvers/overset/icoDyMOversetFoam/correctPhi.H new file mode 100644 index 000000000..53eaae5a7 --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/correctPhi.H @@ -0,0 +1,57 @@ +{ + wordList pcorrTypes(p.boundaryField().types()); + + for (label i = 0; i < p.boundaryField().size(); i++) + { + if (p.boundaryField()[i].fixesValue()) + { + pcorrTypes[i] = fixedValueFvPatchScalarField::typeName; + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0), + pcorrTypes + ); + + // Initialise flux with interpolated velocity + phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf()); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment + + while (piso.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(rAU, pcorr) == fvc::div(phi) + ); + + // Adjust non-orthogonal fluxes if necessary + om.correctNonOrthoFluxes(pcorrEqn, U); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (piso.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + + // Fluxes are corrected to absolute velocity and further corrected + // later. HJ, 6/Feb/2009 + } + +# include "oversetContinuityErrs.H" +} diff --git a/applications/solvers/overset/icoDyMOversetFoam/createControls.H b/applications/solvers/overset/icoDyMOversetFoam/createControls.H new file mode 100644 index 000000000..63d39deef --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/createControls.H @@ -0,0 +1,11 @@ +#include "createTimeControls.H" + +bool correctPhi +( + piso.dict().lookupOrDefault("correctPhi", false) +); + +bool checkMeshCourantNo +( + piso.dict().lookupOrDefault("checkMeshCourantNo", false) +); diff --git a/applications/solvers/overset/icoDyMOversetFoam/createFields.H b/applications/solvers/overset/icoDyMOversetFoam/createFields.H new file mode 100644 index 000000000..77193f386 --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/createFields.H @@ -0,0 +1,79 @@ + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + dimensionedScalar nu + ( + transportProperties.lookup("nu") + ); + + 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" + + // Set reference cell taking into account whether implicit or explicit + // overset aglorithm is used for field p. + const oversetMesh& om = oversetMesh::New(mesh); + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, piso.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(p.name()); + if (correctPhi) + { + mesh.schemesDict().setFluxRequired("pcorr"); + } + + Info<< "Reading field rAU if present\n" << endl; + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dt", dimTime, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + +# include "writeOversetMasks.H" diff --git a/applications/solvers/overset/icoDyMOversetFoam/icoDyMOversetFoam.C b/applications/solvers/overset/icoDyMOversetFoam/icoDyMOversetFoam.C new file mode 100644 index 000000000..7b0ed13b1 --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/icoDyMOversetFoam.C @@ -0,0 +1,167 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + icoDyMOversetFoam + +Description + Transient solver for incompressible, laminar flow of Newtonian fluids + with dynamic mesh and overset mesh support. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" +#include "pisoControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" +# include "createTime.H" +# include "createDynamicFvMesh.H" + + pisoControl piso(mesh); + +# include "initContinuityErrs.H" +# include "initTotalVolume.H" +# include "createControls.H" +# include "createFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readControls.H" +# include "checkTotalVolume.H" + + // Make the fluxes absolute + fvc::makeAbsolute(phi, U); + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + reduce(meshChanged, orOp()); + +# include "createOversetMasks.H" + + if (correctPhi && (mesh.moving() || meshChanged)) + { +# include "correctPhi.H" + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "oversetCourantNo.H" +# include "setDeltaT.H" + + if (mesh.moving() && checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + +# include "UEqn.H" + + // --- PISO loop + + while (piso.correct()) + { + p.boundaryField().updateCoeffs(); + + rAU = 1.0/UEqn.A(); + oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update + + U = rAU*UEqn.H(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + + phi = fvc::interpolate(U) & mesh.Sf(); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + while (piso.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rAU, p) == fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pEqn, U); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve + ( + mesh.solutionDict().solver(p.select(piso.finalInnerIter())) + ); + + if (piso.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + } + +# include "oversetContinuityErrs.H" + + U -= rAU*fvc::grad(p); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the + // interpolation twice (once in correctBoundaryConditions and once + // in oversetInterpolate). Reorganize. VV, 4/Oct/2016. + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/overset/icoDyMOversetFoam/readControls.H b/applications/solvers/overset/icoDyMOversetFoam/readControls.H new file mode 100644 index 000000000..7cc9b3071 --- /dev/null +++ b/applications/solvers/overset/icoDyMOversetFoam/readControls.H @@ -0,0 +1,5 @@ +#include "readTimeControls.H" + +correctPhi = piso.dict().lookupOrDefault("correctPhi", false); + +checkMeshCourantNo = piso.dict().lookupOrDefault("checkMeshCourantNo", false); diff --git a/applications/solvers/overset/icoOversetFoam/Make/files b/applications/solvers/overset/icoOversetFoam/Make/files new file mode 100644 index 000000000..5523b30a3 --- /dev/null +++ b/applications/solvers/overset/icoOversetFoam/Make/files @@ -0,0 +1,3 @@ +icoOversetFoam.C + +EXE = $(FOAM_APPBIN)/icoOversetFoam diff --git a/applications/solvers/overset/icoOversetFoam/Make/options b/applications/solvers/overset/icoOversetFoam/Make/options new file mode 100644 index 000000000..2bb879819 --- /dev/null +++ b/applications/solvers/overset/icoOversetFoam/Make/options @@ -0,0 +1,10 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -loversetMesh \ + -llduSolvers diff --git a/applications/solvers/overset/icoOversetFoam/createFields.H b/applications/solvers/overset/icoOversetFoam/createFields.H new file mode 100644 index 000000000..cc0702878 --- /dev/null +++ b/applications/solvers/overset/icoOversetFoam/createFields.H @@ -0,0 +1,74 @@ + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + dimensionedScalar nu + ( + transportProperties.lookup("nu") + ); + + 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" + + // Set reference cell taking into account whether implicit or explicit + // overset aglorithm is used for field p. + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, piso.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(p.name()); + + Info<< "Reading field rAU if present\n" << endl; + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dt", dimTime, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + +# include "writeOversetMasks.H" diff --git a/applications/solvers/overset/icoOversetFoam/icoOversetFoam.C b/applications/solvers/overset/icoOversetFoam/icoOversetFoam.C new file mode 100644 index 000000000..009831081 --- /dev/null +++ b/applications/solvers/overset/icoOversetFoam/icoOversetFoam.C @@ -0,0 +1,150 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + icoOversetFoam + +Description + Transient solver for incompressible, laminar flow of Newtonian fluids + with overset mesh support. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "pisoControl.H" + +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + pisoControl piso(mesh); + +# include "createOversetMasks.H" +# include "createFields.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + +# include "oversetCourantNo.H" + + // Momentum equation + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + - fvm::laplacian(nu, U) + ); + + if (piso.momentumPredictor()) + { + solve(UEqn == -fvc::grad(p)); + } + + // --- PISO loop + + while (piso.correct()) + { + p.boundaryField().updateCoeffs(); + + rAU = 1.0/UEqn.A(); + oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update + + U = rAU*UEqn.H(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + + phi = fvc::interpolate(U) & mesh.Sf(); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + while (piso.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rAU, p) == fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pEqn, U); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve + ( + mesh.solutionDict().solver(p.select(piso.finalInnerIter())) + ); + + if (piso.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + } + +# include "oversetContinuityErrs.H" + + U -= rAU*fvc::grad(p); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the + // interpolation twice (once in correctBoundaryConditions and once + // in oversetInterpolate). Reorganize. VV, 4/Oct/2016. + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/overset/interOversetFoam/Make/files b/applications/solvers/overset/interOversetFoam/Make/files new file mode 100644 index 000000000..b757bc5c0 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/Make/files @@ -0,0 +1,3 @@ +interOversetFoam.C + +EXE = $(FOAM_APPBIN)/interOversetFoam diff --git a/applications/solvers/overset/interOversetFoam/Make/options b/applications/solvers/overset/interOversetFoam/Make/options new file mode 100644 index 000000000..42ca88a42 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/Make/options @@ -0,0 +1,21 @@ +EXE_INC = \ + -I$(LIB_SRC)/turbulenceModels \ + -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude + +EXE_LIBS = \ + -linterfaceProperties \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume \ + -loversetMesh \ + -llduSolvers diff --git a/applications/solvers/overset/interOversetFoam/UEqn.H b/applications/solvers/overset/interOversetFoam/UEqn.H new file mode 100644 index 000000000..49404121b --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/UEqn.H @@ -0,0 +1,33 @@ + surfaceScalarField muEff + ( + "muEff", + twoPhaseProperties.muf() + + fvc::interpolate(rho*turbulence->nut()) + ); + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(rhoPhi, U) + - fvm::laplacian(muEff, U) + - (fvc::grad(U) & fvc::grad(muEff)) + ); + + UEqn.relax(); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + - fvc::snGrad(pd) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/overset/interOversetFoam/alphaEqn.H b/applications/solvers/overset/interOversetFoam/alphaEqn.H new file mode 100644 index 000000000..ff09cd698 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/alphaEqn.H @@ -0,0 +1,60 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + // New formulation of phir: Co condition + // Note faceOversetMask. HJ, 16/Jun/2015 + surfaceScalarField phir + ( + "phir", + faceOversetMask*interface.cAlpha()*interface.nHatf()* + min + ( + scalar(0.5)/ // This is ref compression Co number for cAlpha = 1 + ( + mesh.surfaceInterpolation::deltaCoeffs()* + mesh.time().deltaT() + ), + mag(phi/mesh.magSf()) + + dimensionedScalar("small", dimVelocity, 1e-3) + ) + ); + + fvScalarMatrix alpha1Eqn + ( + fvm::ddt(alpha1) + + fvm::div(phi, alpha1, alphaScheme) + + fvm::div + ( + -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + alpha1, + alpharScheme + ) + ); + + alpha1Eqn.relax(); + alpha1Eqn.solve(); + + // Update mass flux + rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2; + + Info<< "Liquid phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(cellOversetMask*alpha1).value() + << " Max(alpha1) = " << max(cellOversetMask*alpha1).value() + << endl; + + // Limit alpha to handle overset cells + alpha1.max(0); + alpha1.min(1); + + // Update of interface and density + interface.correct(); + twoPhaseProperties.correct(); + + // Calculate density using limited alpha1 + rho = twoPhaseProperties.rho(); + + // Parallel update in snGrad. HJ, 8/Dec/2012 + rho.correctBoundaryConditions(); +} diff --git a/applications/solvers/overset/interOversetFoam/correctPhi.H b/applications/solvers/overset/interOversetFoam/correctPhi.H new file mode 100644 index 000000000..d117b375c --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/correctPhi.H @@ -0,0 +1,60 @@ +{ + // Set flux required + mesh.schemesDict().setFluxRequired("pcorr"); + + wordList pcorrTypes(p.boundaryField().types()); + + for (label i = 0; i < p.boundaryField().size(); i++) + { + if (p.boundaryField()[i].fixesValue()) + { + pcorrTypes[i] = fixedValueFvPatchScalarField::typeName; + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0), + pcorrTypes + ); + + // Initialise flux with interpolated velocity + phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf()); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(rAU, pcorr) == fvc::div(phi) + ); + + // Adjust non-orthogonal fluxes if necessary + om.correctNonOrthoFluxes(pcorrEqn, U); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + + // Fluxes are corrected to absolute velocity and further corrected + // later. HJ, 6/Feb/2009 + } + +# include "oversetContinuityErrs.H" +} diff --git a/applications/solvers/overset/interOversetFoam/createFields.H b/applications/solvers/overset/interOversetFoam/createFields.H new file mode 100644 index 000000000..350c0c72e --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/createFields.H @@ -0,0 +1,137 @@ + Info<< "Reading field pd\n" << endl; + volScalarField pd + ( + IOobject + ( + "pd", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 + ( + IOobject + ( + "alpha1", + 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; + twoPhaseMixture twoPhaseProperties(U, phi, "alpha1"); + + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); + const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); + + + // Need to store rho for ddt(rho, U) + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + alpha1*rho1 + (scalar(1) - alpha1)*rho2, + alpha1.boundaryField().types() + ); + rho.oldTime(); + + + // Mass flux + // Initialisation does not matter because rhoPhi is reset after the + // alpha1 solution before it is used in the U equation. + surfaceScalarField rhoPhi + ( + IOobject + ( + "rho*phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho1*phi + ); + + + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + pd + rho*gh + ); + + // Set reference cell taking into account whether implicit or explicit + // overset aglorithm is used for field p. + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, pimple.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(pd.name()); + mesh.schemesDict().setFluxRequired(alpha1.name()); + + // Construct interface from alpha1 distribution + interfaceProperties interface(alpha1, U, twoPhaseProperties); + + // Construct incompressible turbulence model + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) + ); + + Info<< "Reading field rAU if present\n" << endl; + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + runTime.deltaT()/rho1, + zeroGradientFvPatchScalarField::typeName + ); + +# include "writeOversetMasks.H" diff --git a/applications/solvers/overset/interOversetFoam/interOversetFoam.C b/applications/solvers/overset/interOversetFoam/interOversetFoam.C new file mode 100644 index 000000000..a232a74a4 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/interOversetFoam.C @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + interOversetFoam + +Description + Solver for 2 incompressible, isothermal immiscible fluids using a VOF + (volume of fluid) phase-fraction based interface capturing approach, + with overset mesh support + + 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. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "interfaceProperties.H" +#include "twoPhaseMixture.H" +#include "turbulenceModel.H" +#include "pimpleControl.H" + +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + pimpleControl pimple(mesh); + +# include "readGravitationalAcceleration.H" +# include "initContinuityErrs.H" +# include "createOversetMasks.H" +# include "createFields.H" +# include "createTimeControls.H" +# include "correctPhi.H" +# include "CourantNo.H" +# include "setInitialDeltaT.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readTimeControls.H" +# include "oversetCourantNo.H" +# include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // Pressure-velocity corrector + while (pimple.loop()) + { + twoPhaseProperties.correct(); + +# include "alphaEqn.H" + +# include "UEqn.H" + + // --- PISO loop + while (pimple.correct()) + { +# include "pdEqn.H" + } + +# include "oversetContinuityErrs.H" + +# include "limitU.H" + + p = pd + cellOversetMask*rho*gh; + + 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/overset/interOversetFoam/limitU.H b/applications/solvers/overset/interOversetFoam/limitU.H new file mode 100644 index 000000000..7ac6f7b44 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/limitU.H @@ -0,0 +1,29 @@ +{ + scalar limitMagU = readScalar(pimple.dict().lookup("limitMagU")); + + volScalarField magU(mag(U)); + + scalar maxMagU = max(magU).value(); + + Info<< "mag(U): max: " << maxMagU + << " avg: " << magU.weightedAverage(mesh.V()).value(); + + if (maxMagU > limitMagU) + { + U.internalField() *= + neg(magU.internalField() - limitMagU) + + pos(magU.internalField() - limitMagU)* + limitMagU/(magU.internalField() + SMALL); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the + // interpolation twice (once in correctBoundaryConditions and once in + // oversetInterpolate). Reorganize. VV, 4/Oct/2016. + + Info << " ...limiting" << endl; + } + else + { + Info << endl; + } +} diff --git a/applications/solvers/overset/interOversetFoam/pdEqn.H b/applications/solvers/overset/interOversetFoam/pdEqn.H new file mode 100644 index 000000000..df1f35aa2 --- /dev/null +++ b/applications/solvers/overset/interOversetFoam/pdEqn.H @@ -0,0 +1,65 @@ +{ + pd.boundaryField().updateCoeffs(); + + rAU = 1.0/UEqn.A(); + oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update + surfaceScalarField rAUf = fvc::interpolate(rAU); + + U = rAU*UEqn.H(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + + surfaceScalarField phiU + ( + "phiU", + fvc::interpolate(U) & mesh.Sf() + ); + + phi = phiU + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf(); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pdEqn + ( + fvm::laplacian(rAUf, pd) == fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pdEqn, U); + + pdEqn.setReference(pRefCell, pRefValue); + + pdEqn.solve + ( + mesh.solutionDict().solver(pd.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pdEqn.flux(); + } + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + } + + // Explicitly relax pressure for momentum corrector + if (!pimple.finalIter()) + { + p.relax(); + } + + U += rAU*fvc::reconstruct((phi - phiU)/rAUf); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the interpolation + // twice (once in correctBoundaryConditions and once in oversetInterpolate) + // Reorganize. VV, 4/Oct/2016. +} diff --git a/applications/solvers/overset/laplacianOversetFoam/Make/files b/applications/solvers/overset/laplacianOversetFoam/Make/files new file mode 100644 index 000000000..1c684c0fc --- /dev/null +++ b/applications/solvers/overset/laplacianOversetFoam/Make/files @@ -0,0 +1,3 @@ +laplacianOversetFoam.C + +EXE = $(FOAM_APPBIN)/laplacianOversetFoam diff --git a/applications/solvers/overset/laplacianOversetFoam/Make/options b/applications/solvers/overset/laplacianOversetFoam/Make/options new file mode 100644 index 000000000..60140c729 --- /dev/null +++ b/applications/solvers/overset/laplacianOversetFoam/Make/options @@ -0,0 +1,8 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -loversetMesh \ + -llduSolvers diff --git a/applications/solvers/overset/laplacianOversetFoam/createFields.H b/applications/solvers/overset/laplacianOversetFoam/createFields.H new file mode 100644 index 000000000..3359f2896 --- /dev/null +++ b/applications/solvers/overset/laplacianOversetFoam/createFields.H @@ -0,0 +1,37 @@ + Info<< "Reading field T\n" << endl; + + volScalarField T + ( + IOobject + ( + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + + Info<< "Reading diffusivity DT\n" << endl; + + dimensionedScalar DT + ( + transportProperties.lookup("DT") + ); diff --git a/applications/solvers/overset/laplacianOversetFoam/laplacianOversetFoam.C b/applications/solvers/overset/laplacianOversetFoam/laplacianOversetFoam.C new file mode 100644 index 000000000..e776e213c --- /dev/null +++ b/applications/solvers/overset/laplacianOversetFoam/laplacianOversetFoam.C @@ -0,0 +1,99 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + laplaciantOversetFoam + +Description + Solves a simple Laplace equation, e.g. for thermal diffusion in a solid + with support for overset meshes. + + Experimental: no overset mesh changes required + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "simpleControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" + +# include "createTime.H" +# include "createMesh.H" + + simpleControl simple(mesh); + +# include "createFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nCalculating temperature distribution\n" << endl; + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + T.correctBoundaryConditions(); + Pout<< "T: " << T.internalField() << endl; + return 0; + + fvScalarMatrix TEqn + ( + fvm::laplacian(DT, T) + ); + Pout<< "TEqn: " << TEqn << endl; +// TEqn.solve(); + + volScalarField residual + ( + "residual", + T + ); + + residual.internalField() = TEqn.residual(); +// residual.boundaryField() == 0; + residual.write(); + + Info<< "residual " << gSumMag(residual.internalField()) << endl; + +# 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/overset/laplacianOversetFoam/write.H b/applications/solvers/overset/laplacianOversetFoam/write.H new file mode 100644 index 000000000..a3cec5419 --- /dev/null +++ b/applications/solvers/overset/laplacianOversetFoam/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/overset/pimpleDyMOversetFoam/Make/files b/applications/solvers/overset/pimpleDyMOversetFoam/Make/files new file mode 100644 index 000000000..5980944cd --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/Make/files @@ -0,0 +1,3 @@ +pimpleDyMOversetFoam.C + +EXE = $(FOAM_APPBIN)/pimpleDyMOversetFoam diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/Make/options b/applications/solvers/overset/pimpleDyMOversetFoam/Make/options new file mode 100644 index 000000000..9810903fc --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/Make/options @@ -0,0 +1,25 @@ +EXE_INC = \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetDynamicFvMesh/lnInclude + +EXE_LIBS = \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume \ + -llduSolvers \ + -L$(MESQUITE_LIB_DIR) -lmesquite \ + -loversetMesh \ + -loversetDynamicFvMesh diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/UEqn.H b/applications/solvers/overset/pimpleDyMOversetFoam/UEqn.H new file mode 100644 index 000000000..d87bb7db1 --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/UEqn.H @@ -0,0 +1,19 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + + turbulence->divDevReff() + ); + + UEqn.relax + ( + mesh.solutionDict().equationRelaxationFactor + ( + U.select(pimple.finalIter()) + ) + ); + + if (pimple.momentumPredictor()) + { + solve(UEqn == -fvc::grad(p)); + } diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/correctPhi.H b/applications/solvers/overset/pimpleDyMOversetFoam/correctPhi.H new file mode 100644 index 000000000..d289344da --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/correctPhi.H @@ -0,0 +1,57 @@ +{ + wordList pcorrTypes(p.boundaryField().types()); + + for (label i = 0; i < p.boundaryField().size(); i++) + { + if (p.boundaryField()[i].fixesValue()) + { + pcorrTypes[i] = fixedValueFvPatchScalarField::typeName; + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0), + pcorrTypes + ); + + // Initialise flux with interpolated velocity + phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf()); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(rAU, pcorr) == fvc::div(phi) + ); + + // Adjust non-orthogonal fluxes if necessary + om.correctNonOrthoFluxes(pcorrEqn, U); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + + // Fluxes are corrected to absolute velocity and further corrected + // later. HJ, 6/Feb/2009 + } + +# include "oversetContinuityErrs.H" +} diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/createControls.H b/applications/solvers/overset/pimpleDyMOversetFoam/createControls.H new file mode 100644 index 000000000..ca5e25906 --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/createControls.H @@ -0,0 +1,11 @@ +#include "createTimeControls.H" + +bool correctPhi +( + pimple.dict().lookupOrDefault("correctPhi", false) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/createFields.H b/applications/solvers/overset/pimpleDyMOversetFoam/createFields.H new file mode 100644 index 000000000..79a42639d --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/createFields.H @@ -0,0 +1,68 @@ + 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" + + // Set reference cell taking into account whether implicit or explicit + // overset aglorithm is used for field p. + const oversetMesh& om = oversetMesh::New(mesh); + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, pimple.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(p.name()); + if (correctPhi) + { + mesh.schemesDict().setFluxRequired("pcorr"); + } + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + + Info<< "Reading field rAU if present\n" << endl; + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + runTime.deltaT(), + zeroGradientFvPatchScalarField::typeName + ); + +# include "writeOversetMasks.H" diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/pEqn.H b/applications/solvers/overset/pimpleDyMOversetFoam/pEqn.H new file mode 100644 index 000000000..8e5c887f2 --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/pEqn.H @@ -0,0 +1,76 @@ +{ + p.boundaryField().updateCoeffs(); + + rAU = 1.0/UEqn.A(); + oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update + surfaceScalarField rAUf = fvc::interpolate(rAU); + + U = rAU*UEqn.H(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + + surfaceScalarField phiU + ( + "phiU", + fvc::interpolate(U) & mesh.Sf() + ); + + phi = phiU; + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rAUf, p) == fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pEqn, U); + + pEqn.setReference(pRefCell, pRefValue); + + if + ( + pimple.corrPISO() == pimple.nCorrPISO() + && pimple.finalNonOrthogonalIter() + ) + { + pEqn.solve(mesh.solutionDict().solver("pFinal")); + } + else + { + pEqn.solve(); + } + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + } + +# include "oversetContinuityErrs.H" + + // Explicitly relax pressure for momentum corrector + if (!pimple.finalIter()) + { + p.relax(); + } + +# include "movingMeshContinuityErrs.H" + + U += rAU*fvc::reconstruct((phi - phiU)/rAUf); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the interpolation + // twice (once in correctBoundaryConditions and once in oversetInterpolate) + // Reorganize. VV, 4/Oct/2016. + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); +} diff --git a/applications/solvers/overset/pimpleDyMOversetFoam/pimpleDyMOversetFoam.C b/applications/solvers/overset/pimpleDyMOversetFoam/pimpleDyMOversetFoam.C new file mode 100644 index 000000000..626881d21 --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/pimpleDyMOversetFoam.C @@ -0,0 +1,129 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + pimpleDyMOversetFoam.C + +Description + Transient solver for incompressible, flow of Newtonian fluids + with with dynamic mesh and overset 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 "singlePhaseTransportModel.H" +#include "turbulenceModel.H" +#include "dynamicFvMesh.H" +#include "pimpleControl.H" + +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" + +# include "createTime.H" +# include "createDynamicFvMesh.H" + + pimpleControl pimple(mesh); + +# include "initContinuityErrs.H" +# include "createControls.H" +# include "createFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readControls.H" + + // Make the fluxes absolute + fvc::makeAbsolute(phi, U); + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + reduce(meshChanged, orOp()); + +# include "createOversetMasks.H" +# include "volContinuity.H" + + if (checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + + // Mesh motion update + if (correctPhi && meshChanged) + { + // Fluxes will be corrected to absolute velocity + // HJ, 6/Feb/2009 +# include "correctPhi.H" + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "oversetCourantNo.H" +# include "setDeltaT.H" + + // --- PIMPLE loop + while (pimple.loop()) + { +# include "UEqn.H" + + // --- PISO loop + while (pimple.correct()) + { +# include "pEqn.H" + } + + 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/overset/pimpleDyMOversetFoam/readControls.H b/applications/solvers/overset/pimpleDyMOversetFoam/readControls.H new file mode 100644 index 000000000..9f982e260 --- /dev/null +++ b/applications/solvers/overset/pimpleDyMOversetFoam/readControls.H @@ -0,0 +1,5 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", false); + +checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false); diff --git a/applications/solvers/overset/potentialDyMOversetFoam/Make/files b/applications/solvers/overset/potentialDyMOversetFoam/Make/files new file mode 100644 index 000000000..55b579515 --- /dev/null +++ b/applications/solvers/overset/potentialDyMOversetFoam/Make/files @@ -0,0 +1,3 @@ +potentialDyMOversetFoam.C + +EXE = $(FOAM_APPBIN)/potentialDyMOversetFoam diff --git a/applications/solvers/overset/potentialDyMOversetFoam/Make/options b/applications/solvers/overset/potentialDyMOversetFoam/Make/options new file mode 100644 index 000000000..47a5fd0fb --- /dev/null +++ b/applications/solvers/overset/potentialDyMOversetFoam/Make/options @@ -0,0 +1,18 @@ +EXE_INC = \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetDynamicFvMesh/lnInclude + +EXE_LIBS = \ + -ldynamicFvMesh \ + -ldynamicMesh \ + -lengine \ + -lmeshTools \ + -lfiniteVolume \ + -loversetMesh \ + -loversetDynamicFvMesh \ + -llduSolvers \ + -L$(MESQUITE_LIB_DIR) -lmesquite diff --git a/applications/solvers/overset/potentialDyMOversetFoam/createFields.H b/applications/solvers/overset/potentialDyMOversetFoam/createFields.H new file mode 100644 index 000000000..ff90e4893 --- /dev/null +++ b/applications/solvers/overset/potentialDyMOversetFoam/createFields.H @@ -0,0 +1,66 @@ + 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 + ); + + surfaceScalarField phi + ( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fvc::interpolate(U) & mesh.Sf() + ); + + // Write overset div phi for visualisation purposes + volScalarField oversetDivPhi + ( + IOobject + ( + "oversetDivPhi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fvc::div(phi) + ); + + // Set reference cell taking into account whether implicit or explicit + // overset algorithm is used for field p. + const oversetMesh& om = oversetMesh::New(mesh); + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, piso.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(p.name()); + + // Write overset masks +# include "writeOversetMasks.H" diff --git a/applications/solvers/overset/potentialDyMOversetFoam/potentialDyMOversetFoam.C b/applications/solvers/overset/potentialDyMOversetFoam/potentialDyMOversetFoam.C new file mode 100644 index 000000000..b85c0ae58 --- /dev/null +++ b/applications/solvers/overset/potentialDyMOversetFoam/potentialDyMOversetFoam.C @@ -0,0 +1,190 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + potentialDyMOversetFoam + +Description + Transient solver for potential flow with dynamic overset mesh. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "pisoControl.H" + +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::validOptions.insert("reconstructU", ""); + +# include "setRootCase.H" +# include "createTime.H" +# include "createDynamicFvMesh.H" + + pisoControl piso(mesh); + +# include "createFields.H" +# include "initTotalVolume.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { +# include "checkTotalVolume.H" + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + reduce(meshChanged, orOp()); + +# include "createOversetMasks.H" + + // Update moving wall velocity boundary condition and calculate the flux + U.correctBoundaryConditions(); + + // Forced overset update: make sure the overset interpolation is + // performed regardless of whether the coupledFringe is specified. + oversetFvPatchVectorField::oversetInterpolate(U); + + phi == (linearInterpolate(U) & mesh.Sf()); + + // Resetting pressure field + p.internalField() = 0; + +# include "volContinuity.H" +# include "meshCourantNo.H" + + // Solve potential flow equations + + // Adjust fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + while (piso.correctNonOrthogonal()) + { + p.storePrevIter(); + + Info<< "Initial flux contour continuity error = " + << mag(sum(phi.boundaryField())) + << endl; + + fvScalarMatrix pEqn + ( + fvm::laplacian + ( + dimensionedScalar + ( + "1", + dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0), + 1 + ), + p + ) + == + fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pEqn, U); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + if (piso.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); +# include "oversetContinuityErrs.H" + } + else + { + p.relax(); + } + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + } + + // Update div phi field for visualisation purposes + oversetDivPhi = cellOversetMask*fvc::div(phi); + + if (args.optionFound("reconstructU")) + { + U = fvc::reconstruct(phi); + U.correctBoundaryConditions(); + } + + Info<< "Interpolated U error = " + << ( + sqrt + ( + sum + ( + sqr + ( + faceOversetMask* + ( + (fvc::interpolate(U) & mesh.Sf()) + - phi + ) + ) + ) + )/sum(mesh.magSf()) + ).value() + << endl; + + // Calculate velocity magnitude + { + volScalarField magU = mag(U); + + Info<< "mag(U): max: " << gMax(magU.internalField()) + << " min: " << gMin(magU.internalField()) << endl; + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/overset/potentialOversetFoam/Make/files b/applications/solvers/overset/potentialOversetFoam/Make/files new file mode 100644 index 000000000..1e33100e3 --- /dev/null +++ b/applications/solvers/overset/potentialOversetFoam/Make/files @@ -0,0 +1,3 @@ +potentialOversetFoam.C + +EXE = $(FOAM_APPBIN)/potentialOversetFoam diff --git a/applications/solvers/overset/potentialOversetFoam/Make/options b/applications/solvers/overset/potentialOversetFoam/Make/options new file mode 100644 index 000000000..2bb879819 --- /dev/null +++ b/applications/solvers/overset/potentialOversetFoam/Make/options @@ -0,0 +1,10 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/oversetMesh/oversetMesh/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -loversetMesh \ + -llduSolvers diff --git a/applications/solvers/overset/potentialOversetFoam/createFields.H b/applications/solvers/overset/potentialOversetFoam/createFields.H new file mode 100644 index 000000000..82e16d476 --- /dev/null +++ b/applications/solvers/overset/potentialOversetFoam/createFields.H @@ -0,0 +1,56 @@ + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh + ); + + // Do not reset p + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + if (args.optionFound("resetU")) + { + U.internalField() = vector::zero; + } + + surfaceScalarField phi + ( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fvc::interpolate(U) & mesh.Sf() + ); + + + // Set reference cell taking into account whether implicit or explicit + // overset algorithm is used for field p. + label pRefCell = 0; + scalar pRefValue = 0.0; + om.setRefCell(p, simple.dict(), pRefCell, pRefValue); + + mesh.schemesDict().setFluxRequired(p.name()); diff --git a/applications/solvers/overset/potentialOversetFoam/potentialOversetFoam.C b/applications/solvers/overset/potentialOversetFoam/potentialOversetFoam.C new file mode 100644 index 000000000..6df89804f --- /dev/null +++ b/applications/solvers/overset/potentialOversetFoam/potentialOversetFoam.C @@ -0,0 +1,262 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + potentialOversetFoam + +Description + Potential flow solver with overset mesh support. + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "simpleControl.H" + +#include "oversetMesh.H" +#include "oversetFvPatchFields.H" +#include "oversetAdjustPhi.H" +#include "globalOversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::validOptions.insert("resetU", ""); + argList::validOptions.insert("writep", ""); + +# include "setRootCase.H" + +# include "createTime.H" +# include "createMesh.H" + + simpleControl simple(mesh); + +# include "createOversetMasks.H" +# include "createFields.H" + +# include "writeOversetMasks.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< nl << "Calculating potential flow" << endl; + + // Do correctors over the complete set + while (simple.correctNonOrthogonal()) + { + phi = (linearInterpolate(U) & mesh.Sf()); + + // Adjust overset fluxes + oversetAdjustPhi(phi, U); // Fringe flux adjustment + globalOversetAdjustPhi(phi, U, p); // Global flux adjustment + + Info<< "Initial flux contour continuity error = " + << mag(sum(phi.boundaryField())) + << endl; + + p.storePrevIter(); + + fvScalarMatrix pEqn + ( + fvm::laplacian + ( + dimensionedScalar + ( + "1", + dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0), + 1 + ), + p + ) + == + fvc::div(phi) + ); + + // Adjust non-orthogonal fringe fluxes if necessary + om.correctNonOrthoFluxes(pEqn, U); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + // Correct the flux + phi -= pEqn.flux(); + + volScalarField divFlux + ( + "divFlux", + cellOversetMask*fvc::div(phi) + ); + divFlux.write(); + + // Perform overset interpolation (after flux reconstruction) + oversetFvPatchScalarField::oversetInterpolate(p); + + if (!simple.finalNonOrthogonalIter()) + { + p.relax(); + } + + Info<< "p min " << gMin(p.internalField()) + << " max " << gMax(p.internalField()) + << " masked min " + << gMin(cellOversetMask.internalField()*p.internalField()) + << " max " + << gMax(cellOversetMask.internalField()*p.internalField()) + << endl; + + Info<< "continuity error = " + << mag + ( + cellOversetMask*fvc::div(phi) + )().weightedAverage(mesh.V()).value() + << endl; + + Info<< "Contour continuity error = " + << mag(sum(phi.boundaryField())) + << endl; + + U = fvc::reconstruct(phi); + U.correctBoundaryConditions(); + oversetFvPatchVectorField::oversetInterpolate(U); // Overset update + // Note: if implicit fringe is switched on, we are doing the + // interpolation twice (once in correctBoundaryConditions and once + // in oversetInterpolate). Reorganize. VV, 4/Oct/2016. + + Info<< "Interpolated U error = " + << ( + sqrt + ( + sum + ( + sqr + ( + faceOversetMask* + ( + (fvc::interpolate(U) & mesh.Sf()) + - phi + ) + ) + ) + )/sum(mesh.magSf()) + ).value() + << endl; + } + + // Calculate velocity magnitude + { + volScalarField magU = cellOversetMask*mag(U); + + Info << "Masked mag(U): max: " << gMax(magU.internalField()) + << " min: " << gMin(magU.internalField()) << endl; + } + + // Force the write + U.write(); + p.write(); + phi.write(); + cellOversetMask.write(); + + if (args.optionFound("writep")) + { + // Find reference patch + label refPatch = -1; + scalar maxMagU = 0; + + // Go through all velocity patches and find the one that fixes + // velocity to the largest value + + forAll (U.boundaryField(), patchI) + { + const fvPatchVectorField& Upatch = U.boundaryField()[patchI]; + + if (Upatch.fixesValue()) + { + // Calculate mean velocity + scalar u = sum(mag(Upatch)); + label patchSize = Upatch.size(); + + reduce(u, sumOp()); + reduce(patchSize, sumOp