From a1558b88f5a28b292cd5f56cc1dd919899f56e6c Mon Sep 17 00:00:00 2001 From: Andrew Heather <> Date: Thu, 6 Dec 2018 22:28:43 +0000 Subject: [PATCH] ENH: Added new utility to create a box of isotropic turbulence --- .../preProcessing/createBoxTurb/Make/files | 3 + .../preProcessing/createBoxTurb/Make/options | 12 ++ .../createBoxTurb/createBlockMesh.H | 104 ++++++++++ .../createBoxTurb/createBoxTurb.C | 184 ++++++++++++++++++ .../createBoxTurb/createFields.H | 33 ++++ 5 files changed, 336 insertions(+) create mode 100644 applications/utilities/preProcessing/createBoxTurb/Make/files create mode 100644 applications/utilities/preProcessing/createBoxTurb/Make/options create mode 100644 applications/utilities/preProcessing/createBoxTurb/createBlockMesh.H create mode 100644 applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C create mode 100644 applications/utilities/preProcessing/createBoxTurb/createFields.H diff --git a/applications/utilities/preProcessing/createBoxTurb/Make/files b/applications/utilities/preProcessing/createBoxTurb/Make/files new file mode 100644 index 0000000000..718509c756 --- /dev/null +++ b/applications/utilities/preProcessing/createBoxTurb/Make/files @@ -0,0 +1,3 @@ +createBoxTurb.C + +EXE = $(FOAM_APPBIN)/createBoxTurb diff --git a/applications/utilities/preProcessing/createBoxTurb/Make/options b/applications/utilities/preProcessing/createBoxTurb/Make/options new file mode 100644 index 0000000000..8e62b99f63 --- /dev/null +++ b/applications/utilities/preProcessing/createBoxTurb/Make/options @@ -0,0 +1,12 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/mesh/blockMesh/lnInclude \ + -I$(LIB_SRC)/fileFormats/lnInclude + + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lblockMesh \ + -lfileFormats diff --git a/applications/utilities/preProcessing/createBoxTurb/createBlockMesh.H b/applications/utilities/preProcessing/createBoxTurb/createBlockMesh.H new file mode 100644 index 0000000000..a4c482a954 --- /dev/null +++ b/applications/utilities/preProcessing/createBoxTurb/createBlockMesh.H @@ -0,0 +1,104 @@ +const cellModel& hex = cellModel::ref(cellModel::HEX); + +cellShapeList cellShapes; +faceListList boundary; +pointField points; +{ + Info<< "Creating block" << endl; + + block b + ( + cellShape(hex, identity(8), false), + pointField + ( + { + point(0, 0, 0), + point(L.x(), 0, 0), + point(L.x(), L.y(), 0), + point(0, L.y(), 0), + point(0, 0, L.z()), + point(L.x(), 0, L.z()), + point(L.x(), L.y(), L.z()), + point(0, L.y(), L.z()) + } + ), + blockEdgeList(), + blockFaceList(), + N, + List(12) + ); + + Info<< "Creating cells" << endl; + + List> bCells(b.cells()); + cellShapes.setSize(bCells.size()); + forAll(cellShapes, celli) + { + cellShapes[celli] = + cellShape(hex, labelList(bCells[celli]), false); + } + + Info<< "Creating boundary faces" << endl; + + boundary.setSize(b.boundaryPatches().size()); + forAll(boundary, patchi) + { + faceList faces(b.boundaryPatches()[patchi].size()); + forAll(faces, facei) + { + faces[facei] = face(b.boundaryPatches()[patchi][facei]); + } + boundary[patchi].transfer(faces); + } + + points.transfer(const_cast(b.points())); +} + +Info<< "Creating patch dictionaries" << endl; +wordList patchNames(boundary.size()); +forAll(patchNames, patchi) +{ + patchNames[patchi] = "patch" + Foam::name(patchi); +} + +PtrList boundaryDicts(boundary.size()); +forAll(boundaryDicts, patchi) +{ + boundaryDicts.set(patchi, new dictionary()); + dictionary& patchDict = boundaryDicts[patchi]; + word nbrPatchName; + if (patchi % 2 == 0) + { + nbrPatchName = "patch" + Foam::name(patchi + 1); + } + else + { + nbrPatchName = "patch" + Foam::name(patchi - 1); + } + + patchDict.add("type", cyclicPolyPatch::typeName); + patchDict.add("neighbourPatch", nbrPatchName); +} + +Info<< "Creating polyMesh" << endl; +polyMesh mesh +( + IOobject + ( + polyMesh::defaultRegion, + runTime.constant(), + runTime, + IOobject::NO_READ + ), + std::move(points), + cellShapes, + boundary, + patchNames, + boundaryDicts, + "defaultFaces", + cyclicPolyPatch::typeName, + false +); + +Info<< "Writing polyMesh" << endl; +mesh.write(); diff --git a/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C b/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C new file mode 100644 index 0000000000..9da5d76795 --- /dev/null +++ b/applications/utilities/preProcessing/createBoxTurb/createBoxTurb.C @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + createBoxTurb + +Description + Creates a box of isotropic turbulence based on a user-specified + energy spectrum. + + Based on the reference + \verbatim + Saad, T., Cline, D., Stoll, R., Sutherland, J.C. + "Scalable Tools for Generating Synthetic Isotropic Turbulence with + Arbitrary Spectra" + AIAA Journal, Vol. 55, No. 1 (2017), pp. 327-331. + \endverbatim + + The \c -createBlockMesh option creates a block mesh and exits, which + can then be decomposed and the utility run in parallel. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "block.H" +#include "mathematicalConstants.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::vector randomUnitVector(Random& rndGen) +{ + // Sample point on a sphere + scalar t = rndGen.globalPosition(-1, 1); + scalar phim = rndGen.globalSample01()*mathematical::twoPi; + scalar thetam = Foam::acos(t); + + return vector + ( + Foam::sin(thetam*Foam::cos(phim)), + Foam::sin(thetam*Foam::sin(phim)), + Foam::cos(thetam) + ); +} + + +int main(int argc, char *argv[]) +{ + argList::addBoolOption + ( + "createBlockMesh", + "create the block mesh and exit" + ); + + #include "setRootCase.H" + + #include "createTime.H" + #include "createFields.H" + + if (args.found("createBlockMesh")) + { + // Create a box block mesh with cyclic patches + #include "createBlockMesh.H" + return 0; + } + + #include "createMesh.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + // Minimum wave number + scalar kappa0 = mathematical::twoPi/cmptMin(L); + + // Maximum wave number + scalar kappaMax = mathematical::pi/cmptMin(delta); + + Info<< "Wave number min/max = " << kappa0 << ", " << kappaMax << endl; + + Info<< "Generating velocity field" << endl; + + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector(dimVelocity, Zero) + ); + + vectorField& Uc = U.primitiveFieldRef(); + const scalar deltaKappa = (kappaMax - kappa0)/scalar(nModes - 1); + const vectorField& C(mesh.C()); + for (label modei = 1; modei <= nModes; ++modei) + { + // Equidistant wave mode + scalar kappaM = kappa0 + deltaKappa*(modei-1); + + Info<< "Processing mode:" << modei << " kappaM:" << kappaM << endl; + + // Energy + scalar E = Ek->value(kappaM); + + // Wave amplitude + scalar qm = Foam::sqrt(E*deltaKappa); + + // Wave number unit vector + const vector kappaHatm(randomUnitVector(rndGen)); + + vector kappaTildem(0.5*kappaM*cmptMultiply(kappaHatm, delta)); + for (direction i = 0; i < 3; ++i) + { + kappaTildem[i] = 2/delta[i]*Foam::sin(kappaTildem[i]); + } + + // Intermediate unit vector zeta + const vector zetaHatm(randomUnitVector(rndGen)); + + // Unit vector sigma + vector sigmaHatm(zetaHatm^kappaTildem); + sigmaHatm /= mag(kappaTildem); + + // Phase angle + scalar psim = 0.5*rndGen.position(-mathematical::pi, mathematical::pi); + + // Add the velocity contribution per mode + Uc += 2*qm*cos(kappaM*(kappaHatm & C) + psim)*sigmaHatm; + } + + U.write(); + + { + Info<< "Generating kinetic energy field" << endl; + volScalarField k("k", 0.5*magSqr(U)); + k.write(); + Info<< "min/max/average k = " + << gMin(k) << ", " << gMax(k) << ", " << gAverage(k) + << endl; + } + + { + Info<< "Generating div(U) field" << endl; + volScalarField divU(fvc::div(U)); + divU.write(); + Info<< "min/max/average div(U) = " + << gMin(divU) << ", " << gMax(divU) << ", " << gAverage(divU) + << endl; + } + + Info<< nl; + runTime.printExecutionTime(Info); + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/createBoxTurb/createFields.H b/applications/utilities/preProcessing/createBoxTurb/createFields.H new file mode 100644 index 0000000000..2b0fb5e643 --- /dev/null +++ b/applications/utilities/preProcessing/createBoxTurb/createFields.H @@ -0,0 +1,33 @@ +IOdictionary dict +( + IOobject + ( + "createBoxTurbDict", + runTime.constant(), + runTime, + IOobject::MUST_READ + ) +); + +// Extents in x, y, z directions +const vector L(dict.get("L")); + +// Number of cells in x, y, z directions +const Vector