diff --git a/applications/solvers/basic/potentialFoam/potentialFoam.C b/applications/solvers/basic/potentialFoam/potentialFoam.C index e2e8702280..cdada0c353 100644 --- a/applications/solvers/basic/potentialFoam/potentialFoam.C +++ b/applications/solvers/basic/potentialFoam/potentialFoam.C @@ -132,6 +132,7 @@ int main(int argc, char *argv[]) "Execute functionObjects" ); + #include "addRegionOption.H" #include "addCheckCaseOptions.H" #include "setRootCaseLists.H" #include "createTime.H" diff --git a/applications/solvers/basic/potentialFoam/potentialFoam.C.orig b/applications/solvers/basic/potentialFoam/potentialFoam.C.orig new file mode 100644 index 0000000000..ec1c198493 --- /dev/null +++ b/applications/solvers/basic/potentialFoam/potentialFoam.C.orig @@ -0,0 +1,254 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- + | Copyright (C) 2011-2016 OpenFOAM Foundation +------------------------------------------------------------------------------- +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 + potentialFoam + +Group + grpBasicSolvers + +Description + Potential flow solver which solves for the velocity potential, to + calculate the flux-field, from which the velocity field is obtained by + reconstructing the flux. + + \heading Solver details + The potential flow solution is typically employed to generate initial fields + for full Navier-Stokes codes. The flow is evolved using the equation: + + \f[ + \laplacian \Phi = \div(\vec{U}) + \f] + + Where: + \vartable + \Phi | Velocity potential [m2/s] + \vec{U} | Velocity [m/s] + \endvartable + + The corresponding pressure field could be calculated from the divergence + of the Euler equation: + + \f[ + \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0 + \f] + + but this generates excessive pressure variation in regions of large + velocity gradient normal to the flow direction. A better option is to + calculate the pressure field corresponding to velocity variation along the + stream-lines: + + \f[ + \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0 + \f] + where the flow direction tensor \f$\vec{F}\f$ is obtained from + \f[ + \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}} + \f] + + \heading Required fields + \plaintable + U | Velocity [m/s] + \endplaintable + + \heading Optional fields + \plaintable + p | Kinematic pressure [m2/s2] + Phi | Velocity potential [m2/s] + | Generated from p (if present) or U if not present + \endplaintable + + \heading Options + \plaintable + -writep | write the Euler pressure + -writePhi | Write the final velocity potential + -initialiseUBCs | Update the velocity boundaries before solving for Phi + \endplaintable + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "pisoControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::addNote + ( + "Potential flow solver which solves for the velocity potential" + ); + + argList::addOption + ( + "pName", + "pName", + "Name of the pressure field" + ); + + argList::addBoolOption + ( + "initialiseUBCs", + "Initialise U boundary conditions" + ); + + argList::addBoolOption + ( + "writePhi", + "Write the final velocity potential field" + ); + + argList::addBoolOption + ( + "writep", + "Calculate and write the Euler pressure field" + ); + + argList::addBoolOption + ( + "withFunctionObjects", + "Execute functionObjects" + ); + + #include "addCheckCaseOptions.H" + #include "setRootCaseLists.H" + #include "createTime.H" + #include "createMesh.H" + + pisoControl potentialFlow(mesh, "potentialFlow"); + + #include "createFields.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< nl << "Calculating potential flow" << endl; + + // Since solver contains no time loop it would never execute + // function objects so do it ourselves + runTime.functionObjects().start(); + + MRF.makeRelative(phi); + adjustPhi(phi, U, p); + + // Non-orthogonal velocity potential corrector loop + while (potentialFlow.correctNonOrthogonal()) + { + fvScalarMatrix PhiEqn + ( + fvm::laplacian(dimensionedScalar("1", dimless, 1), Phi) + == + fvc::div(phi) + ); + + PhiEqn.setReference(PhiRefCell, PhiRefValue); + PhiEqn.solve(); + + if (potentialFlow.finalNonOrthogonalIter()) + { + phi -= PhiEqn.flux(); + } + } + + MRF.makeAbsolute(phi); + + Info<< "Continuity error = " + << mag(fvc::div(phi))().weightedAverage(mesh.V()).value() + << endl; + + U = fvc::reconstruct(phi); + U.correctBoundaryConditions(); + + Info<< "Interpolated velocity error = " + << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value() + << endl; + + // Write U and phi + U.write(); + phi.write(); + + // Optionally write Phi + if (args.found("writePhi")) + { + Phi.write(); + } + + // Calculate the pressure field from the Euler equation + if (args.found("writep")) + { + Info<< nl << "Calculating approximate pressure field" << endl; + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell + ( + p, + potentialFlow.dict(), + pRefCell, + pRefValue + ); + + // Calculate the flow-direction filter tensor + volScalarField magSqrU(magSqr(U)); + volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU))); + + // Calculate the divergence of the flow-direction filtered div(U*U) + // Filtering with the flow-direction generates a more reasonable + // pressure distribution in regions of high velocity gradient in the + // direction of the flow + volScalarField divDivUU + ( + fvc::div + ( + F & fvc::div(phi, U), + "div(div(phi,U))" + ) + ); + + // Solve a Poisson equation for the approximate pressure + while (potentialFlow.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(p) + divDivUU + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + } + + p.write(); + } + + runTime.functionObjects().end(); + + runTime.printExecutionTime(Info); + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/include/createMesh.H b/src/OpenFOAM/include/createMesh.H index 06d07591b2..da6fbbccd9 100644 --- a/src/OpenFOAM/include/createMesh.H +++ b/src/OpenFOAM/include/createMesh.H @@ -1,4 +1,5 @@ Foam::autoPtr meshPtr(nullptr); +Foam::word regionName = Foam::fvMesh::defaultRegion; if (args.found("dry-run") || args.found("dry-run-write")) { @@ -6,10 +7,16 @@ if (args.found("dry-run") || args.found("dry-run-write")) << "Operating in 'dry-run' mode: case will run for 1 time step. " << "All checks assumed OK on a clean exit" << Foam::endl; + // Allow region in combination with dry-run + args.readIfPresent("region", regionName); + Foam::FieldBase::allowConstructFromLargerSize = true; // Create a simplified 1D mesh and attempt to re-create boundary conditions - meshPtr.reset(new Foam::simplifiedMeshes::columnFvMesh(runTime)); + meshPtr.reset + ( + new Foam::simplifiedMeshes::columnFvMesh(runTime, regionName) + ); // Stop after 1 iteration of the simplified mesh @@ -30,9 +37,18 @@ if (args.found("dry-run") || args.found("dry-run-write")) } else { - Foam::Info - << "Create mesh for time = " - << runTime.timeName() << Foam::nl << Foam::endl; + if (args.readIfPresent("region", regionName)) + { + Foam::Info + << "Create mesh " << regionName << " for time = " + << runTime.timeName() << Foam::nl << Foam::endl; + } + else + { + Foam::Info + << "Create mesh for time = " + << runTime.timeName() << Foam::nl << Foam::endl; + } meshPtr.reset ( @@ -40,7 +56,7 @@ else ( Foam::IOobject ( - Foam::fvMesh::defaultRegion, + regionName, runTime.timeName(), runTime, Foam::IOobject::MUST_READ diff --git a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C index 31a0ea8ab4..df47e836be 100644 --- a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C +++ b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C @@ -67,14 +67,14 @@ Foam::simplifiedMeshes::simplifiedDynamicFvMeshBase::New if (cstrIter.found()) { Info<< "Selecting simplified mesh model " << modelType << endl; - return autoPtr(cstrIter()(io.time())); + return autoPtr(cstrIter()(io.time(), io.name())); } } Info<< "Selecting simplified mesh model " << staticFvMesh::typeName << endl; return autoPtr ( - new SimplifiedDynamicFvMesh(io.time()) + new SimplifiedDynamicFvMesh(io.time(), io.name()) ); } diff --git a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.H b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.H index f8ac3c722d..4925fd1e12 100644 --- a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.H +++ b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.H @@ -67,9 +67,10 @@ public: dynamicFvMesh, time, ( - const Time& runTime + const Time& runTime, + const word& regionName ), - (runTime) + (runTime, regionName) ); @@ -101,7 +102,7 @@ public: ClassNameNoDebug(DynamicMeshType::typeName_.c_str()); //- Constructor - SimplifiedDynamicFvMesh(const Time& runTime); + SimplifiedDynamicFvMesh(const Time& runTime, const word& regionName); //- Update the mesh for both mesh motion and topology change virtual bool update() diff --git a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMeshTemplates.C b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMeshTemplates.C index 37209b20bd..ac9429597f 100644 --- a/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMeshTemplates.C +++ b/src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMeshTemplates.C @@ -29,16 +29,17 @@ template Foam::simplifiedMeshes::SimplifiedDynamicFvMesh:: SimplifiedDynamicFvMesh ( - const Time& runTime + const Time& runTime, + const word& regionName ) : simplifiedDynamicFvMeshBase(), - columnFvMeshInfo(runTime), + columnFvMeshInfo(runTime, regionName), DynamicMeshType ( IOobject ( - fvMesh::defaultRegion, + regionName, runTime.constant(), runTime, IOobject::NO_READ, // Do not read any existing mesh diff --git a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C index 36f7283e0e..bcca22f366 100644 --- a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C +++ b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C @@ -62,7 +62,7 @@ bool Foam::simplifiedMeshes::columnFvMeshInfo::setPatchEntries ( "boundary", localInstance_, - polyMesh::meshSubDir, + regionPrefix_ + polyMesh::meshSubDir, runTime, IOobject::MUST_READ, IOobject::NO_WRITE, @@ -97,7 +97,12 @@ bool Foam::simplifiedMeshes::columnFvMeshInfo::setPatchEntries else { // No boundary file - try reading from a field - IOobjectList objects(runTime, runTime.timeName()); + IOobjectList objects + ( + runTime, + runTime.timeName(), + (regionName_ == fvMesh::defaultRegion ? "" : regionName_) + ); if (objects.empty()) { @@ -202,7 +207,7 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::initialise(const Time& runTime) ( "points", localInstance_, - polyMesh::meshSubDir, + regionPrefix_ + polyMesh::meshSubDir, runTime, IOobject::MUST_READ, IOobject::NO_WRITE, @@ -371,12 +376,7 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::addLocalPatches if (debug) { - Pout<< "patches:" << nl << endl; - forAll(patches, patchi) - { - Pout<< "patch: " << patches[patchi]->name() << nl - << *patches[patchi] << endl; - } + Pout<< "patches:" << nl << mesh.boundaryMesh() << endl; } } @@ -400,13 +400,24 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::initialiseZones(fvMesh& mesh) // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo(const Time& runTime) +Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo +( + const Time& runTime, + const word& regionName +) : + regionName_(regionName), + regionPrefix_ + ( + regionName_ == fvMesh::defaultRegion + ? "" + : regionName_ + '/' + ), localInstance_ ( runTime.findInstance ( - polyMesh::meshSubDir, + regionPrefix_ + polyMesh::meshSubDir, "boundary", IOobject::READ_IF_PRESENT ) @@ -428,14 +439,18 @@ Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo(const Time& runTime) } -Foam::simplifiedMeshes::columnFvMesh::columnFvMesh(const Time& runTime) +Foam::simplifiedMeshes::columnFvMesh::columnFvMesh +( + const Time& runTime, + const word& regionName +) : - columnFvMeshInfo(runTime), + columnFvMeshInfo(runTime, regionName), simplifiedFvMesh ( IOobject ( - fvMesh::defaultRegion, + regionName, runTime.constant(), runTime, IOobject::NO_READ, // Do not read any existing mesh diff --git a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H index 01b65bdd9f..6347b5c8ed 100644 --- a/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H +++ b/src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2018-2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -73,6 +73,12 @@ protected: // Protected data + //- Region of existing mesh + const word regionName_; + + //- Additional prefix for region. Empty if default region + const fileName regionPrefix_; + //- Location of existing mesh (if present) const word localInstance_; @@ -118,7 +124,11 @@ public: ClassName("columnFvMeshInfo"); // Constructor - columnFvMeshInfo(const Time& runTime); + columnFvMeshInfo + ( + const Time& runTime, + const word& regionName + ); }; @@ -134,7 +144,11 @@ public: TypeName("columnFvMesh"); //- Constructor - columnFvMesh(const Time& runTime); + columnFvMesh + ( + const Time& runTime, + const word& regionName = fvMesh::defaultRegion + ); };