From 3ba865dadddc9a69986a630ae1a928ec1efe4465 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 08:58:03 +0000 Subject: [PATCH 001/150] Bugfix: correct ramping integral for translation --- .../solidBodyMotion/translation/translation.C | 65 ++++++++++++++----- .../solidBodyMotion/translation/translation.H | 4 +- 2 files changed, 52 insertions(+), 17 deletions(-) diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C index fe5aa5705..037898975 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C @@ -48,24 +48,24 @@ namespace solidBodyMotionFunctions // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -Foam::scalar Foam::solidBodyMotionFunctions::translation::rampFactor() const +Foam::scalar +Foam::solidBodyMotionFunctions::translation::rampFactor() const { - if (rampTime_ > 0) - { - // Calculate ramp-up factor - const scalar t = time_.value(); + const scalar t = time_.value(); - return sin(2*pi/(4*rampTime_)*Foam::min(rampTime_, t)); + if (t < rampTime_) + { + // Ramping region + return sin(pi/(2*rampTime_)*t); } else { + // Past ramping region return 1; } } - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::solidBodyMotionFunctions::translation::translation @@ -77,7 +77,17 @@ Foam::solidBodyMotionFunctions::translation::translation solidBodyMotionFunction(SBMFCoeffs, runTime), velocity_(SBMFCoeffs_.lookup("velocity")), rampTime_(readScalar(SBMFCoeffs_.lookup("rampTime"))) -{} +{ + if (rampTime_ < 0) + { + FatalIOErrorIn + ( + "solidBodyMotionFunctions::translation::translation", + SBMFCoeffs_ + ) << "Negative rampTime not allowed." + << abort(FatalIOError); + } +} // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // @@ -91,12 +101,39 @@ Foam::solidBodyMotionFunctions::translation::~translation() Foam::septernion Foam::solidBodyMotionFunctions::translation::transformation() const { - scalar time = time_.value(); + const scalar t = time_.value(); - septernion TR(rampFactor()*velocity_*time, quaternion::I); + septernion TR; + + if (t < rampTime_) + { + // Ramping region + // Account for ramping using analytical integration from ramped velocity + // distribution in this region + TR = septernion + ( + velocity_*2*rampTime_/pi*(1 - cos(pi/(2*rampTime_)*t)), + quaternion::I + ); + } + else + { + // Past ramping region + TR = septernion + ( + velocity_* + ( + // Displacement during the ramping region + 2*rampTime_/pi + // Displacement during constant velocity after ramping region + + (t - rampTime_) + ), + quaternion::I + ); + } Info<< "solidBodyMotionFunctions::translation::transformation(): " - << "Time = " << time << " velocity = " << rampFactor()*velocity_ + << "Time = " << t << " velocity = " << rampFactor()*velocity_ << " transformation = " << TR << endl; @@ -107,12 +144,10 @@ Foam::solidBodyMotionFunctions::translation::transformation() const Foam::septernion Foam::solidBodyMotionFunctions::translation::velocity() const { - scalar time = time_.value(); - septernion TV(rampFactor()*velocity_, quaternion::zero); Info<< "solidBodyMotionFunctions::translation::transformation(): " - << "Time = " << time << " velocity: " << TV << endl; + << "Time = " << time_.value() << " velocity: " << TV << endl; return TV; } diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.H b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.H index 82109c4c1..e7b4d1678 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.H +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.H @@ -25,7 +25,7 @@ Class Foam::solidBodyMotionFunctions::translation Description - translation motion function + Translation motion function with sine law ramping for velocity. SourceFiles translation.C @@ -70,7 +70,7 @@ class translation void operator=(const translation&); - //- Ramping factor resulting from rampTime_ value + //- Velocity ramping factor resulting from rampTime_ value scalar rampFactor() const; From 20e7bdd79463d587f5e3777c97a188e5d2926842 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:43:24 +0000 Subject: [PATCH 002/150] Tutorial update --- .../shockTube/constant/polyMesh/boundary | 34 +++++++++++++++++++ .../pimpleDyMFoam/wingMotion/Allrun | 4 +-- .../constant/polyMesh/boundary | 22 ++++++------ .../constant/polyMesh/boundary | 22 ++++++------ .../constant/polyMesh/boundary | 24 ++++++++----- 5 files changed, 73 insertions(+), 33 deletions(-) create mode 100644 tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary new file mode 100644 index 000000000..798403436 --- /dev/null +++ b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/polyMesh/boundary @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | For copyright notice see file Copyright | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +2 +( + sides + { + type patch; + nFaces 2; + startFace 999; + } + empty + { + type empty; + nFaces 4000; + startFace 1001; + } +) + +// ************************************************************************* // diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun b/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun index 1ae239a15..22a9194a1 100755 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun @@ -19,8 +19,8 @@ runApplication simpleFoam # Copy the mesh from the steady state case and map the results to a # mesh motion case, then solve transient. cd ../wingMotion2D_pimpleDyMFoam -cp -r ../wingMotion2D_simpleFoam/constant/polyMesh/* constant/polyMesh/ -cp -r 0.org 0 +\cp -r ../wingMotion2D_simpleFoam/constant/polyMesh constant/polyMesh +\cp -r 0.org 0 runApplication mapFields ../wingMotion2D_simpleFoam -sourceTime latestTime -consistent mv 0/pointDisplacement.unmapped 0/pointDisplacement runApplication decomposePar diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/polyMesh/boundary b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/polyMesh/boundary index 62020026c..92ce545c1 100644 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/polyMesh/boundary +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/polyMesh/boundary @@ -3,7 +3,7 @@ | \\ / F ield | foam-extend: Open Source CFD | | \\ / O peration | Version: 3.2 | | \\ / A nd | Web: http://www.foam-extend.org | -| \\/ M anipulation | | +| \\/ M anipulation | For copyright notice see file Copyright | \*---------------------------------------------------------------------------*/ FoamFile { @@ -21,37 +21,37 @@ FoamFile { type patch; nFaces 150; - startFace 76473; + startFace 71816; } inlet { type patch; nFaces 48; - startFace 76623; + startFace 71966; } outlet { type patch; nFaces 48; - startFace 76671; + startFace 72014; } wing { type wall; nFaces 778; - startFace 76719; + startFace 72062; } - back + symBack { type empty; - nFaces 38129; - startFace 77497; + nFaces 35801; + startFace 72840; } - front + symFront { type empty; - nFaces 38129; - startFace 115626; + nFaces 35801; + startFace 108641; } ) diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/constant/polyMesh/boundary b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/constant/polyMesh/boundary index 62020026c..92ce545c1 100644 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/constant/polyMesh/boundary +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/constant/polyMesh/boundary @@ -3,7 +3,7 @@ | \\ / F ield | foam-extend: Open Source CFD | | \\ / O peration | Version: 3.2 | | \\ / A nd | Web: http://www.foam-extend.org | -| \\/ M anipulation | | +| \\/ M anipulation | For copyright notice see file Copyright | \*---------------------------------------------------------------------------*/ FoamFile { @@ -21,37 +21,37 @@ FoamFile { type patch; nFaces 150; - startFace 76473; + startFace 71816; } inlet { type patch; nFaces 48; - startFace 76623; + startFace 71966; } outlet { type patch; nFaces 48; - startFace 76671; + startFace 72014; } wing { type wall; nFaces 778; - startFace 76719; + startFace 72062; } - back + symBack { type empty; - nFaces 38129; - startFace 77497; + nFaces 35801; + startFace 72840; } - front + symFront { type empty; - nFaces 38129; - startFace 115626; + nFaces 35801; + startFace 108641; } ) diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/constant/polyMesh/boundary b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/constant/polyMesh/boundary index c8f40206f..dfb6fb9f4 100644 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/constant/polyMesh/boundary +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/constant/polyMesh/boundary @@ -3,7 +3,7 @@ | \\ / F ield | foam-extend: Open Source CFD | | \\ / O peration | Version: 3.2 | | \\ / A nd | Web: http://www.foam-extend.org | -| \\/ M anipulation | | +| \\/ M anipulation | For copyright notice see file Copyright | \*---------------------------------------------------------------------------*/ FoamFile { @@ -15,37 +15,43 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -5 +6 ( topAndBottom { type patch; nFaces 150; - startFace 7077; + startFace 2087919; } inlet { type patch; nFaces 48; - startFace 7227; + startFace 2088069; } outlet { type patch; nFaces 48; - startFace 7275; + startFace 2088117; } symFront { type symmetryPlane; - nFaces 3600; - startFace 7323; + nFaces 35801; + startFace 2088165; } symBack { type symmetryPlane; - nFaces 3600; - startFace 10923; + nFaces 35792; + startFace 2123966; + } + wing_5degrees.obj_WALL10 + { + type wall; + nFaces 49792; + startFace 2159758; } ) From aa1fa870f4b3ef45ac63376611f47425e7abdf2a Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:44:23 +0000 Subject: [PATCH 003/150] Add div phid upder-relaxation in the pressure equation --- .../solvers/compressible/steadyCompressibleFoam/pEqn.H | 9 +++++++++ .../compressible/steadyCompressibleMRFFoam/pEqn.H | 9 +++++++++ .../compressible/steadyCompressibleSRFFoam/pEqn.H | 9 +++++++++ 3 files changed, 27 insertions(+) diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H index 355dae922..5736d3d6f 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H @@ -21,12 +21,21 @@ p.storePrevIter(); + volScalarField divPhid + ( + "divPhid", + fvc::div(phid) + ); + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psis, p) + fvm::div(phid, p) + // Convective flux relaxation terms + + fvm::SuSp(-divPhid, p) + + divPhid*p + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index b45e1573d..868f74040 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -25,12 +25,21 @@ p.storePrevIter(); + volScalarField divPhid + ( + "divPhid", + fvc::div(phid) + ); + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psis, p) + fvm::div(phid, p) + // Convective flux relaxation terms + + fvm::SuSp(-divPhid, p) + + divPhid*p + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H index f161a3c4d..0718dd8d2 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H @@ -20,12 +20,21 @@ p.storePrevIter(); + volScalarField divPhid + ( + "divPhid", + fvc::div(phid) + ); + for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psis, p) + fvm::div(phid, p) + // Convective flux relaxation terms + + fvm::SuSp(-divPhid, p) + + divPhid*p + fvc::div(phid2) - fvm::laplacian(rho*rUrelA, p) ); From 05fa8214cec17004b1d7fb39e2a12d0892443266 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:44:42 +0000 Subject: [PATCH 004/150] Bugfix: remove zones before writing the mesh --- applications/utilities/mesh/manipulation/splitMesh/splitMesh.C | 3 +++ 1 file changed, 3 insertions(+) diff --git a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C index 7699bfd11..cc66597fa 100644 --- a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C +++ b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C @@ -263,6 +263,9 @@ int main(int argc, char *argv[]) splitter.changeMesh(); + // Remove zones - not needed + mesh.removeZones(); + Info<< "Writing mesh to " << runTime.timeName() << endl; if (!mesh.write()) { From 425569a02d78c1df04b8a82167934164bd9dc1f4 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:45:39 +0000 Subject: [PATCH 005/150] Bugfix: added noMotion --- .../meshMotion/solidBodyMotion/Make/files | 1 + .../solidBodyMotion/noMotion/noMotion.C | 91 +++++++++++++++ .../solidBodyMotion/noMotion/noMotion.H | 107 ++++++++++++++++++ 3 files changed, 199 insertions(+) create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files index 2d8de31b9..1e02ad3a4 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files @@ -1,6 +1,7 @@ solidBodyMotionFunction/solidBodyMotionFunction.C solidBodyMotionFunction/newSolidBodyMotionFunction.C +noMotion/noMotion.C translation/translation.C linearOscillation/linearOscillation.C rotatingOscillation/rotatingOscillation.C diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C new file mode 100644 index 000000000..d1bdc2bf9 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "noMotion.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + defineTypeNameAndDebug(noMotion, 0); + addToRunTimeSelectionTable + ( + solidBodyMotionFunction, + noMotion, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::noMotion::noMotion +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::noMotion::~noMotion() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::noMotion::transformation() const +{ + return septernion(vector::zero, quaternion::I); +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::noMotion::velocity() const +{ + return septernion(vector::zero, quaternion::zero); +} + + +bool Foam::solidBodyMotionFunctions::noMotion::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + return true; +} + + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H new file mode 100644 index 000000000..f2bd0219e --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::solidBodyMotionFunctions::noMotion + +Description + No motion + +SourceFiles + noMotion.C + +\*---------------------------------------------------------------------------*/ + +#ifndef noMotion_H +#define noMotion_H + +#include "solidBodyMotionFunction.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class noMotion Declaration +\*---------------------------------------------------------------------------*/ + +class noMotion +: + public solidBodyMotionFunction +{ + // Private Member Functions + + //- Disallow copy construct + noMotion(const noMotion&); + + //- Disallow default bitwise assignment + void operator=(const noMotion&); + + +public: + + //- Runtime type information + TypeName("noMotion"); + + + // Constructors + + //- Construct from components + noMotion + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + // Destructor + + virtual ~noMotion(); + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 2d4d2879d89f303aa794cbf980599aabd4f1309d Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:45:56 +0000 Subject: [PATCH 006/150] Formatting --- src/foam/containers/Lists/PriorityList/PriorityList.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/foam/containers/Lists/PriorityList/PriorityList.H b/src/foam/containers/Lists/PriorityList/PriorityList.H index 98010378f..6b76f8220 100644 --- a/src/foam/containers/Lists/PriorityList/PriorityList.H +++ b/src/foam/containers/Lists/PriorityList/PriorityList.H @@ -50,7 +50,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class PriorityList Declaration + Class PriorityList Declaration \*---------------------------------------------------------------------------*/ template From 400afe8ec6e618bb4a06272652a905944fff3e0b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:47:03 +0000 Subject: [PATCH 007/150] Formatting --- .../immersedBoundaryFvPatch/immersedBoundaryFvPatch.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C index 556cb9cca..a189ae83c 100644 --- a/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C +++ b/src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C @@ -1035,7 +1035,7 @@ void Foam::immersedBoundaryFvPatch::makeIbPointsAndNormals() const ) << "Can't find nearest triSurface point for cell " << ibc[cellI] << ", " << mesh_.cellCentres()[ibc[cellI]] - << "Hit data = " << pih << nl + << ". Hit data = " << pih << nl << abort(FatalError); } From a55eb64d4cd62cbf9809769afff7191c7dd8cf31 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 10:46:56 +0000 Subject: [PATCH 008/150] Limiters: preparation --- src/foam/Make/files | 1 + .../Fields/VectorNFields/VectorNFields.H | 3 ++ .../fields/Fields/tensorField/tensorField.C | 8 ++- .../fields/Fields/tensorField/tensorField.H | 7 ++- .../fields/Fields/vectorField/vectorField.C | 50 +++++++++++++++++++ .../fields/Fields/vectorField/vectorField.H | 14 ++++++ src/foam/primitives/Tensor/TensorTemplateI.H | 13 +++++ src/foam/primitives/Vector/VectorTemplateI.H | 8 +++ src/foam/primitives/VectorN/TensorN.H | 5 +- src/foam/primitives/VectorN/TensorNI.H | 23 +++++++++ src/foam/primitives/VectorN/VectorN.H | 2 +- src/foam/primitives/VectorN/VectorNI.H | 28 ++++++++--- 12 files changed, 149 insertions(+), 13 deletions(-) create mode 100644 src/foam/fields/Fields/vectorField/vectorField.C diff --git a/src/foam/Make/files b/src/foam/Make/files index 012816086..7cf866868 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -533,6 +533,7 @@ fields/cloud/cloud.C Fields = fields/Fields $(Fields)/labelField/labelField.C $(Fields)/scalarField/scalarField.C +$(Fields)/vectorField/vectorField.C $(Fields)/sphericalTensorField/sphericalTensorField.C $(Fields)/diagTensorField/diagTensorField.C $(Fields)/symmTensorField/symmTensorField.C diff --git a/src/foam/fields/Fields/VectorNFields/VectorNFields.H b/src/foam/fields/Fields/VectorNFields/VectorNFields.H index 839a71c61..e5a1af17d 100644 --- a/src/foam/fields/Fields/VectorNFields/VectorNFields.H +++ b/src/foam/fields/Fields/VectorNFields/VectorNFields.H @@ -52,6 +52,9 @@ UNARY_FUNCTION(CmptType, vectorType, cmptSum) \ BINARY_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ BINARY_TYPE_FUNCTION(vectorType, vectorType, vectorType, cmptMultiply) \ \ +BINARY_FUNCTION(vectorType, vectorType, CmptType, scaleRow) \ +BINARY_TYPE_FUNCTION(vectorType, vectorType, CmptType, scaleRow) \ + \ BINARY_OPERATOR(vectorType, CmptType, vectorType, /, divide) \ BINARY_TYPE_OPERATOR(vectorType, CmptType, vectorType, /, divide) \ \ diff --git a/src/foam/fields/Fields/tensorField/tensorField.C b/src/foam/fields/Fields/tensorField/tensorField.C index abb67e786..7bf8b3344 100644 --- a/src/foam/fields/Fields/tensorField/tensorField.C +++ b/src/foam/fields/Fields/tensorField/tensorField.C @@ -35,7 +35,7 @@ License namespace Foam { -// * * * * * * * * * * * * * * * global functions * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // UNARY_FUNCTION(scalar, tensor, tr) UNARY_FUNCTION(sphericalTensor, tensor, sph) @@ -160,7 +160,11 @@ tmp > transformFieldMask } -// * * * * * * * * * * * * * * * global operators * * * * * * * * * * * * * // +BINARY_FUNCTION(tensor, tensor, vector, scaleRow) +BINARY_TYPE_FUNCTION(tensor, tensor, vector, scaleRow) + + +// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // UNARY_OPERATOR(vector, tensor, *, hdual) UNARY_OPERATOR(tensor, vector, *, hdual) diff --git a/src/foam/fields/Fields/tensorField/tensorField.H b/src/foam/fields/Fields/tensorField/tensorField.H index e85284d02..3831c8b9e 100644 --- a/src/foam/fields/Fields/tensorField/tensorField.H +++ b/src/foam/fields/Fields/tensorField/tensorField.H @@ -51,7 +51,7 @@ namespace Foam typedef Field tensorField; -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // UNARY_FUNCTION(scalar, tensor, tr) UNARY_FUNCTION(sphericalTensor, tensor, sph) @@ -70,8 +70,11 @@ UNARY_FUNCTION(tensor, tensor, hinv) UNARY_FUNCTION(vector, symmTensor, eigenValues) UNARY_FUNCTION(tensor, symmTensor, eigenVectors) +BINARY_FUNCTION(tensor, tensor, vector, scaleRow) +BINARY_TYPE_FUNCTION(tensor, tensor, vector, scaleRow) -// * * * * * * * * * * * * * * * global operators * * * * * * * * * * * * * // + +// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // UNARY_OPERATOR(vector, tensor, *, hdual) UNARY_OPERATOR(tensor, vector, *, hdual) diff --git a/src/foam/fields/Fields/vectorField/vectorField.C b/src/foam/fields/Fields/vectorField/vectorField.C new file mode 100644 index 000000000..9baf788fd --- /dev/null +++ b/src/foam/fields/Fields/vectorField/vectorField.C @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "vectorField.H" + +#define TEMPLATE +#include "FieldFunctionsM.C" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +BINARY_FUNCTION(vector, vector, scalar, scaleRow) +BINARY_TYPE_FUNCTION(vector, vector, scalar, scaleRow) + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "undefFieldFunctionsM.H" + +// ************************************************************************* // diff --git a/src/foam/fields/Fields/vectorField/vectorField.H b/src/foam/fields/Fields/vectorField/vectorField.H index 3f085969e..2042ec01d 100644 --- a/src/foam/fields/Fields/vectorField/vectorField.H +++ b/src/foam/fields/Fields/vectorField/vectorField.H @@ -38,6 +38,9 @@ SourceFiles #include "scalarField.H" #include "vector.H" +#define TEMPLATE +#include "FieldFunctionsM.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -47,12 +50,23 @@ namespace Foam typedef Field vectorField; + +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +BINARY_FUNCTION(vector, vector, scalar, scaleRow) +BINARY_TYPE_FUNCTION(vector, vector, scalar, scaleRow) + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#include "undefFieldFunctionsM.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/foam/primitives/Tensor/TensorTemplateI.H b/src/foam/primitives/Tensor/TensorTemplateI.H index 4e7849f8e..2dcf53582 100644 --- a/src/foam/primitives/Tensor/TensorTemplateI.H +++ b/src/foam/primitives/Tensor/TensorTemplateI.H @@ -567,6 +567,19 @@ inline Cmpt invariantIII(const Tensor& t) } +// Scale row +template +inline Tensor scaleRow(const Tensor& t, const Vector& v) +{ + return Tensor + ( + t.xx()*v.x(), t.xy()*v.x(), t.xz()*v.x(), + t.yx()*v.y(), t.yy()*v.y(), t.yz()*v.y(), + t.zx()*v.z(), t.zy()*v.z(), t.zz()*v.z() + ); +} + + // * * * * * * * * * Mixed Tensor SphericalTensor Operators * * * * * * * * // template diff --git a/src/foam/primitives/Vector/VectorTemplateI.H b/src/foam/primitives/Vector/VectorTemplateI.H index 44c357a0b..fa9772536 100644 --- a/src/foam/primitives/Vector/VectorTemplateI.H +++ b/src/foam/primitives/Vector/VectorTemplateI.H @@ -156,6 +156,14 @@ inline Vector operator^(const Vector& v1, const Vector& v2) } +template +inline Vector scaleRow(const Vector& v, const Cmpt& c) +{ + // Multiply by scalar + return v*c; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/foam/primitives/VectorN/TensorN.H b/src/foam/primitives/VectorN/TensorN.H index dd837cf07..94bc2bca7 100644 --- a/src/foam/primitives/VectorN/TensorN.H +++ b/src/foam/primitives/VectorN/TensorN.H @@ -121,15 +121,16 @@ public: const direction j ); - //- Diagonal + //- Return diagonal inline DiagTensorN diag() const; - //- Transpose + //- Return transpose inline TensorN T() const; //- Negative sum the vertical off-diagonal components inline TensorN negSumDiag() const; + // Member Operators //- Assign to a SphericalTensorN diff --git a/src/foam/primitives/VectorN/TensorNI.H b/src/foam/primitives/VectorN/TensorNI.H index 9e80eec06..360169bdc 100644 --- a/src/foam/primitives/VectorN/TensorNI.H +++ b/src/foam/primitives/VectorN/TensorNI.H @@ -1007,6 +1007,29 @@ inline TensorN negSumDiag(const TensorN& t) return t.negSumDiag(); } + +//- Scale tensor row with a VectorN +template +inline TensorN scaleRow +( + const TensorN& t, + const VectorN& v +) +{ + TensorN result; + + for (label i = 0; i < TensorN::rowLength; i++) + { + for (label j = 0; j < TensorN::rowLength; j++) + { + result(i, j) = t(i, j)*v(j); + } + } + + return result; +} + + //- Return the component sum // template // inline Cmpt sum(const TensorN& t) diff --git a/src/foam/primitives/VectorN/VectorN.H b/src/foam/primitives/VectorN/VectorN.H index dfca611cc..1e4d9e5cc 100644 --- a/src/foam/primitives/VectorN/VectorN.H +++ b/src/foam/primitives/VectorN/VectorN.H @@ -100,7 +100,7 @@ public: //- Return access to ith component inline Cmpt& operator()(const direction i); - //- Componentwise multiply + //- Component-wise multiply inline VectorN cmptMultiply ( const VectorN& diff --git a/src/foam/primitives/VectorN/VectorNI.H b/src/foam/primitives/VectorN/VectorNI.H index a0759b0eb..4595e5370 100644 --- a/src/foam/primitives/VectorN/VectorNI.H +++ b/src/foam/primitives/VectorN/VectorNI.H @@ -73,7 +73,8 @@ inline VectorN::VectorN template inline VectorN::VectorN(const Cmpt& vx) { - VectorSpaceOps::nComponents,0>::eqOpS(*this, vx, eqOp()); + VectorSpaceOps::nComponents,0>::eqOpS + (*this, vx, eqOp()); } @@ -103,10 +104,12 @@ inline Cmpt& VectorN::operator()(const direction i) //- Multiply components of VectorN by VectorN template -inline VectorN VectorN::cmptMultiply(const VectorN& v) +inline VectorN +VectorN::cmptMultiply(const VectorN& v) { VectorN res; - VectorSpaceOps::nComponents,0>::op(res, *this, v, multiplyOp()); + VectorSpaceOps::nComponents,0>::op + (res, *this, v, multiplyOp()); return res; } @@ -175,20 +178,33 @@ inline VectorN operator/(const scalar s, const VectorN& v) { VectorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, v, divideOp3()); + VectorSpaceOps::nComponents,0>::opSV + (res, s, v, divideOp3()); return res; } //- Multiply components of VectorN by VectorN template -inline VectorN cmptMultiply(const VectorN& v1, const VectorN& v2) +inline VectorN +cmptMultiply(const VectorN& v1, const VectorN& v2) { VectorN res; - VectorSpaceOps::nComponents,0>::op(res, v1, v2, multiplyOp()); + VectorSpaceOps::nComponents,0>::op + (res, v1, v2, multiplyOp()); return res; } + +//- Multiply row of VectorN with a scalar +template +inline VectorN +scaleRow(const VectorN& v, const Cmpt& c) +{ + return v*c; +} + + //- Return the component sum // template // inline Cmpt sum(const VectorN& v) From 22c9eccfc8f59dc50bab8233b0718a70ff170423 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 10:47:58 +0000 Subject: [PATCH 009/150] Gradient limiters --- src/dbns/limiter/BarthJespersenLimiter.H | 232 ------------ src/dbns/limiter/MDLimiter.H | 34 +- src/dbns/limiter/VenkatakrishnanLimiter.H | 273 -------------- src/dbns/limiter/firstOrderLimiter.H | 72 +--- src/finiteVolume/Make/files | 4 + .../gradSchemes/gradScheme/gradScheme.H | 144 ++++---- .../barthJespersenGrad/barthJespersenGrad.H | 139 +++++++ .../barthJespersenGrad/barthJespersenGrads.C | 89 +++++ .../limitedGrad/LimitedGrad.C | 347 ++++++++++++++++++ .../limitedGrad/LimitedGrad.H | 150 ++++++++ .../minModGrad/minModGrad.H | 139 +++++++ .../minModGrad/minModGrads.C | 89 +++++ .../venkatakrishnanGrad/venkatakrishnanGrad.H | 139 +++++++ .../venkatakrishnanGrads.C | 89 +++++ .../limitedGradSchemes/wangGrad/wangGrad.H | 139 +++++++ .../limitedGradSchemes/wangGrad/wangGrads.C | 89 +++++ .../gradientLimiters/BarthJespersenLimiter.H | 125 +++++++ .../gradientLimiters/MinModLimiter.H | 126 +++++++ .../gradientLimiters/VenkatakrishnanLimiter.H | 169 +++++++++ .../gradientLimiters/WangLimiter.H | 170 +++++++++ 20 files changed, 2110 insertions(+), 648 deletions(-) delete mode 100644 src/dbns/limiter/BarthJespersenLimiter.H delete mode 100644 src/dbns/limiter/VenkatakrishnanLimiter.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrad.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrads.C create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.C create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrad.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrads.C create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrad.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrads.C create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrad.H create mode 100644 src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrads.C create mode 100644 src/finiteVolume/finiteVolume/gradientLimiters/BarthJespersenLimiter.H create mode 100644 src/finiteVolume/finiteVolume/gradientLimiters/MinModLimiter.H create mode 100644 src/finiteVolume/finiteVolume/gradientLimiters/VenkatakrishnanLimiter.H create mode 100644 src/finiteVolume/finiteVolume/gradientLimiters/WangLimiter.H diff --git a/src/dbns/limiter/BarthJespersenLimiter.H b/src/dbns/limiter/BarthJespersenLimiter.H deleted file mode 100644 index 334d59f1b..000000000 --- a/src/dbns/limiter/BarthJespersenLimiter.H +++ /dev/null @@ -1,232 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / 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 . - -Class - BarthJespersenLimiter - -Description - Barth-Jespersen limiter - -Author - Aleksandar Jemcov, Hrvoje Jasak - -\*---------------------------------------------------------------------------*/ - -#ifndef BarthJespersenLimiter_H -#define BarthJespersenLimiter_H - -#include "vector.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class BarthJespersenLimiter Declaration -\*---------------------------------------------------------------------------*/ - -class BarthJespersenLimiter -{ -public: - - // Constructors - - //- Construct null - BarthJespersenLimiter() - {} - - - // Desctructor - default - - - // Member functions - - //- Return limiter value - inline scalar limiter - ( - const scalar& cellVolume, - const scalar& deltaOneMax, - const scalar& deltaOneMin, - const scalar& deltaTwo - ) - { - if (mag(deltaTwo) < SMALL) - { - return 1; - } - else - { - scalar deltaTwoStab = stabilise(deltaTwo, VSMALL); - - return min - ( - max - ( - max(deltaOneMax/deltaTwoStab, 0), - max(deltaOneMin/deltaTwoStab, 0) - ), - 1 - ); - } - } - - //- Return scalar limiter - inline scalar limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const scalar& deltaOneMax, - const scalar& deltaOneMin, - const vector& gradPhiP, - const vector& gradPhiN, - const vector& d - ) - { - return limiter - ( - cellVolume, - deltaOneMax, - deltaOneMin, - (d & gradPhiP) - ); - } - - //- Return vector limiter - inline vector limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const vector& deltaOneMax, - const vector& deltaOneMin, - const tensor& gradPhiP, - const tensor& gradPhiN, - const vector& d - ) - { - vector deltaTwo = d & gradPhiP; - - return vector - ( - limiter - ( - cellVolume, - deltaOneMax[vector::X], - deltaOneMin[vector::X], - deltaTwo[vector::X] - ), - limiter - ( - cellVolume, - deltaOneMax[vector::Y], - deltaOneMin[vector::Y], - deltaTwo[vector::Y] - ), - limiter - ( - cellVolume, - deltaOneMax[vector::Z], - deltaOneMin[vector::Z], - deltaTwo[vector::Z] - ) - ); - } - - //- Return limiter for a symmTensor - inline symmTensor limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const symmTensor& deltaOneMax, - const symmTensor& deltaOneMin, - const vector& gradPhiPXX, - const vector& gradPhiNXX, - const vector& gradPhiPXY, - const vector& gradPhiNXY, - const vector& gradPhiPXZ, - const vector& gradPhiNXZ, - const vector& gradPhiPYY, - const vector& gradPhiNYY, - const vector& gradPhiPYZ, - const vector& gradPhiNYZ, - const vector& gradPhiPZZ, - const vector& gradPhiNZZ, - const vector& d - ) - { - return symmTensor - ( - limiter - ( - cellVolume, - deltaOneMax[0], - deltaOneMin[0], - d & gradPhiPXX - ), - limiter - ( - cellVolume, - deltaOneMax[1], - deltaOneMin[1], - d & gradPhiPXY - ), - limiter - ( - cellVolume, - deltaOneMax[2], - deltaOneMin[2], - d & gradPhiPXZ - ), - limiter - ( - cellVolume, - deltaOneMax[3], - deltaOneMin[3], - d & gradPhiPYY - ), - limiter - ( - cellVolume, - deltaOneMax[4], - deltaOneMin[4], - d & gradPhiPYZ - ), - limiter - ( - cellVolume, - deltaOneMax[5], - deltaOneMin[5], - d & gradPhiPZZ - ) - ); - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -#endif - -// ************************************************************************* // diff --git a/src/dbns/limiter/MDLimiter.H b/src/dbns/limiter/MDLimiter.H index 122cb9084..5bfb3557f 100644 --- a/src/dbns/limiter/MDLimiter.H +++ b/src/dbns/limiter/MDLimiter.H @@ -191,9 +191,7 @@ public: const GradFieldType& gradPhiIn = gradPhi.internalField(); // Compute limiter values, internal faces - - // Dummy flux - scalar upwindFlux = 1; + Type phiOwnerLimiter, phiNeighbourLimiter; forAll (owner, faceI) { @@ -203,29 +201,27 @@ public: vector deltaRLeft = faceCentre[faceI] - cellCentre[own]; vector deltaRRight = faceCentre[faceI] - cellCentre[nei]; - Type phiOwnerLimiter = limitFunction.limiter + // Owner side + limitFunction.limiter ( + phiOwnerLimiter, cellVolume[own], - upwindFlux, phiMaxIn[own] - phi_[own], phiMinIn[own] - phi_[own], - gradPhiIn[own], - gradPhiIn[nei], - deltaRLeft + (deltaRLeft & gradPhiIn[own]) ); - Type phiNeighbourLimiter = limitFunction.limiter + // Neighbour side + limitFunction.limiter ( + phiNeighbourLimiter, cellVolume[nei], - upwindFlux, phiMaxIn[nei] - phi_[nei], phiMinIn[nei] - phi_[nei], - gradPhiIn[nei], - gradPhiIn[own], - deltaRRight + (deltaRRight & gradPhiIn[nei]) ); - // find minimal limiter value in each cell + // Find minimal limiter value in each cell phiLimiterIn[own] = min(phiLimiterIn[own], phiOwnerLimiter); @@ -255,19 +251,19 @@ public: const GradFieldType gradPhiRight = gradPhi.boundaryField()[patchI].patchNeighbourField(); + Type phiFCLimiter; + forAll (fc, faceI) { const label& curFC = fc[faceI]; - Type phiFCLimiter = limitFunction.limiter + limitFunction.limiter ( + phiFCLimiter, cellVolume[curFC], - upwindFlux, phiMaxIn[curFC] - phi_[curFC], phiMinIn[curFC] - phi_[curFC], - gradPhiLeft[faceI], - gradPhiRight[faceI], - deltaR[faceI] + (deltaR[faceI] & gradPhiLeft[faceI]) ); // Find minimal limiter value in each cell diff --git a/src/dbns/limiter/VenkatakrishnanLimiter.H b/src/dbns/limiter/VenkatakrishnanLimiter.H deleted file mode 100644 index 8f8c6fe1f..000000000 --- a/src/dbns/limiter/VenkatakrishnanLimiter.H +++ /dev/null @@ -1,273 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / 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 . - -Class - VenkatakrishnanLimiter - -Description - Venkatakrishnan differentiable limiter - -Author - Aleksandar Jemcov - Updated by Hrvoje Jasak - -\*---------------------------------------------------------------------------*/ - -#ifndef VenkatakrishnanLimiter_H -#define VenkatakrishnanLimiter_H - -#include "vector.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class VenkatakrishnanLimiter Declaration -\*---------------------------------------------------------------------------*/ - -class VenkatakrishnanLimiter -{ - // Private data - - //- Limiter value - scalar k_; - -public: - - // Constructor - - //- Construct null - VenkatakrishnanLimiter() - : - k_(5) - {} - - //- Return limiter - inline scalar limiter - ( - const scalar& cellVolume, - const scalar& deltaOneMax, - const scalar& deltaOneMin, - const scalar& deltaTwo - ) - { - scalar epsilonSquare = pow3(k_)*cellVolume; - - if (deltaTwo > 0) - { - return max - ( - min - ( - ( - (sqr(deltaOneMax) + epsilonSquare)*deltaTwo - + 2*sqr(deltaTwo)*deltaOneMax - )/ - stabilise - ( - deltaTwo* - ( - sqr(deltaOneMax) - + 2.0*sqr(deltaTwo) - + deltaOneMax*deltaTwo - + epsilonSquare - ), - SMALL - ), - 1 - ), - 0 - ); - } - else if (deltaTwo < 0) - { - return max - ( - min - ( - ( - (sqr(deltaOneMin) + epsilonSquare)*deltaTwo - + 2*sqr(deltaTwo)*deltaOneMin - )/ - stabilise - ( - deltaTwo* - ( - sqr(deltaOneMin) - + 2.0*sqr(deltaTwo) + deltaOneMin*deltaTwo - + epsilonSquare - ), - SMALL - ), - 1 - ), - 0 - ); - } - else - { - return 1; - } - } - - //- Return limiter for a scalar - inline scalar limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const scalar& deltaOneMax, - const scalar& deltaOneMin, - const vector& gradPhiP, - const vector& gradPhiN, - const vector& d - ) - { - return limiter - ( - cellVolume, - deltaOneMax, - deltaOneMin, - d & gradPhiP - ); - } - - //- Return limiter for a vector - inline vector limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const vector& deltaOneMax, - const vector& deltaOneMin, - const tensor& gradPhiP, - const tensor& gradPhiN, - const vector& d - ) - { - vector deltaTwo = d & gradPhiP; - - return vector - ( - limiter - ( - cellVolume, - deltaOneMax[vector::X], - deltaOneMin[vector::X], - deltaTwo[vector::X] - ), - limiter - ( - cellVolume, - deltaOneMax[vector::Y], - deltaOneMin[vector::Y], - deltaTwo[vector::Y] - ), - limiter - ( - cellVolume, - deltaOneMax[vector::Z], - deltaOneMin[vector::Z], - deltaTwo[vector::Z] - ) - ); - } - - //- Return limiter for a symmTensor - inline symmTensor limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const symmTensor& deltaOneMax, - const symmTensor& deltaOneMin, - const vector& gradPhiPXX, - const vector& gradPhiNXX, - const vector& gradPhiPXY, - const vector& gradPhiNXY, - const vector& gradPhiPXZ, - const vector& gradPhiNXZ, - const vector& gradPhiPYY, - const vector& gradPhiNYY, - const vector& gradPhiPYZ, - const vector& gradPhiNYZ, - const vector& gradPhiPZZ, - const vector& gradPhiNZZ, - const vector& d - ) - { - return symmTensor - ( - limiter - ( - cellVolume, - deltaOneMax[0], - deltaOneMin[0], - d & gradPhiPXX - ), - limiter - ( - cellVolume, - deltaOneMax[1], - deltaOneMin[1], - d & gradPhiPXY - ), - limiter - ( - cellVolume, - deltaOneMax[2], - deltaOneMin[2], - d & gradPhiPXZ - ), - limiter - ( - cellVolume, - deltaOneMax[3], - deltaOneMin[3], - d & gradPhiPYY - ), - limiter - ( - cellVolume, - deltaOneMax[4], - deltaOneMin[4], - d & gradPhiPYZ - ), - limiter - ( - cellVolume, - deltaOneMax[5], - deltaOneMin[5], - d & gradPhiPZZ - ) - ); - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -#endif - -// ************************************************************************* // diff --git a/src/dbns/limiter/firstOrderLimiter.H b/src/dbns/limiter/firstOrderLimiter.H index b497c5683..752b2e6b1 100644 --- a/src/dbns/limiter/firstOrderLimiter.H +++ b/src/dbns/limiter/firstOrderLimiter.H @@ -43,7 +43,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class firstOrderLimiter Declaration + Class firstOrderLimiter Declaration \*---------------------------------------------------------------------------*/ class firstOrderLimiter @@ -56,71 +56,37 @@ public: firstOrderLimiter() {} - //- Return limiter - inline scalar limiter + + // Destructor - default + + + // Member functions + + //- Set scalar limiter value + inline void limiter ( + scalar& lim, const scalar& cellVolume, const scalar& deltaOneMax, const scalar& deltaOneMin, const scalar& deltaTwo ) { - return 0; + lim = 0; } - //- Return limiter for a scalar - inline scalar limiter + //- Set Type limiter + template + inline void limiter ( + Type& lim, const scalar& cellVolume, - const scalar& faceFlux, - const scalar& deltaOneMax, - const scalar& deltaOneMin, - const vector& gradPhiP, - const vector& gradPhiN, - const vector& d + const Type& deltaOneMax, + const Type& deltaOneMin, + const Type& extrapolate ) { - return 0; - } - - //- Return limiter for a vector - inline vector limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const vector& deltaOneMax, - const vector& deltaOneMin, - const tensor& gradPhiP, - const tensor& gradPhiN, - const vector& d - ) - { - return vector::zero; - } - - //- Return limiter for a symmTensor - inline symmTensor limiter - ( - const scalar& cellVolume, - const scalar& faceFlux, - const symmTensor& deltaOneMax, - const symmTensor& deltaOneMin, - const vector& gradPhiPXX, - const vector& gradPhiNXX, - const vector& gradPhiPXY, - const vector& gradPhiNXY, - const vector& gradPhiPXZ, - const vector& gradPhiNXZ, - const vector& gradPhiPYY, - const vector& gradPhiNYY, - const vector& gradPhiPYZ, - const vector& gradPhiNYZ, - const vector& gradPhiPZZ, - const vector& gradPhiNZZ, - const vector& d - ) - { - return symmTensor::zero; + lim = pTraits::zero; } }; diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 3fb994320..272471cc8 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -359,6 +359,10 @@ $(limitedGradSchemes)/faceLimitedGrad/faceLimitedGrads.C $(limitedGradSchemes)/cellLimitedGrad/cellLimitedGrads.C $(limitedGradSchemes)/faceMDLimitedGrad/faceMDLimitedGrads.C $(limitedGradSchemes)/cellMDLimitedGrad/cellMDLimitedGrads.C +$(limitedGradSchemes)/minModGrad/minModGrads.C +$(limitedGradSchemes)/barthJespersenGrad/barthJespersenGrads.C +$(limitedGradSchemes)/venkatakrishnanGrad/venkatakrishnanGrads.C +$(limitedGradSchemes)/wangGrad/wangGrads.C snGradSchemes = finiteVolume/snGradSchemes $(snGradSchemes)/snGradScheme/snGradSchemes.C diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H index f7d45f14e..1c66789d7 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H @@ -63,6 +63,25 @@ class gradScheme : public refCount { +public: + + // Public typedefs + + typedef Field FieldType; + typedef GeometricField GeoFieldType; + typedef typename GeoFieldType::GeometricBoundaryField GeoBoundaryFieldType; + + typedef Field::type> GradFieldType; + typedef GeometricField + < + typename outerProduct::type, + fvPatchField, + volMesh + > GeoGradFieldType; + typedef BlockLduSystem::type> + GradMatrixType; + + // Private data const fvMesh& mesh_; @@ -127,75 +146,58 @@ public: return mesh_; } - //- Calculate and return the grad of the given field. - // Used by grad either to recalculate the cached gradient when it is - // out of date with respect to the field or when it is not cached. - virtual tmp - < - GeometricField - ::type, fvPatchField, volMesh> - > calcGrad - ( - const GeometricField&, - const word& name - ) const = 0; - - //- Calculate and return the grad of the given field - // which may have been cached - tmp - < - GeometricField - ::type, fvPatchField, volMesh> - > grad - ( - const GeometricField&, - const word& name - ) const; - - //- Calculate and return the grad of the given field - // with the default name - // which may have been cached - tmp - < - GeometricField - ::type, fvPatchField, volMesh> - > grad - ( - const GeometricField& - ) const; - - //- Calculate and return the grad of the given field - // with the default name - // which may have been cached - tmp - < - GeometricField - ::type, fvPatchField, volMesh> - > grad - ( - const tmp >& - ) const; - - //- Return the BlockLduSystem corresponding to the implicit grad - // discretization. For block coupled systems. - virtual tmp - < - BlockLduSystem::type> - > fvmGrad - ( - const GeometricField& - ) const; - - // Moved from gaussGrad into base class. HJ, 14/Jun/2013 //- Correct the boundary values of the gradient using the patchField // snGrad functions static void correctBoundaryConditions ( - const GeometricField&, - GeometricField - ::type, fvPatchField, volMesh>& + const GeoFieldType&, + GeoGradFieldType& ); + + + // Gradient functions + + //- Calculate and return the grad of the given field + // which may have been cached + tmp grad + ( + const GeoFieldType&, + const word& name + ) const; + + //- Calculate and return the grad of the given field + // with the default name which may have been cached + tmp grad + ( + const GeoFieldType& + ) const; + + //- Calculate and return the grad of the given field + // with the default name which may have been cached + tmp grad + ( + const tmp& + ) const; + + + // Virtual function interface + + //- Calculate and return the grad of the given field. + // Used by grad either to recalculate the cached gradient when + // out of date with respect to the field or when it is not cached. + virtual tmp calcGrad + ( + const GeoFieldType&, + const word& name + ) const = 0; + + //- Return the BlockLduSystem corresponding to the implicit grad + // discretization. For block coupled systems. + virtual tmp fvmGrad + ( + const GeoFieldType& + ) const; }; @@ -211,17 +213,17 @@ public: // Add the patch constructor functions to the hash tables -#define makeFvGradTypeScheme(SS, Type) \ - \ -defineNamedTemplateTypeNameAndDebug(SS, 0); \ - \ -gradScheme::addIstreamConstructorToTable > \ +#define makeFvGradTypeScheme(SS, Type) \ + \ +defineNamedTemplateTypeNameAndDebug(SS, 0); \ + \ +gradScheme::addIstreamConstructorToTable > \ add##SS##Type##IstreamConstructorToTable_; -#define makeFvGradScheme(SS) \ - \ -makeFvGradTypeScheme(SS, scalar) \ +#define makeFvGradScheme(SS) \ + \ +makeFvGradTypeScheme(SS, scalar) \ makeFvGradTypeScheme(SS, vector) diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrad.H new file mode 100644 index 000000000..5b667b1b4 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrad.H @@ -0,0 +1,139 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::fv::barthJespersenGrad + +Description + Barth-Jespersen gradient limiter applied to a runTime selected + base gradient scheme. + +SourceFiles + barthJespersenGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef barthJespersenGrad_H +#define barthJespersenGrad_H + +#include "gradScheme.H" +#include "LimitedGrad.H" +#include "BarthJespersenLimiter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class barthJespersenGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class barthJespersenGrad +: + public fv::gradScheme, + public LimitedGrad +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + barthJespersenGrad(const barthJespersenGrad&); + + //- Disallow default bitwise assignment + void operator=(const barthJespersenGrad&); + + +public: + + //- RunTime type information + TypeName("barthJespersen"); + + + // Constructors + + //- Construct from mesh and schemeData + barthJespersenGrad(const fvMesh& mesh, Istream& schemeData) + : + gradScheme(mesh), + LimitedGrad(mesh, schemeData) + {} + + + // Member Functions + + //- Return the gradient of the given field to the gradScheme::grad + // for optional caching + virtual tmp + < + GeometricField + ::type, fvPatchField, volMesh> + > calcGrad + ( + const GeometricField& vf, + const word& name + ) const + { + return LimitedGrad::gradientField + ( + vf, + name + ); + } + + //- Return the BlockLduSystem corresponding to the implicit cell + // limited grad discretization. For block coupled systems. + virtual tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& vf + ) const + { + return LimitedGrad::gradientMatrix + ( + vf + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrads.C new file mode 100644 index 000000000..c849d86d4 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/barthJespersenGrad/barthJespersenGrads.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "barthJespersenGrad.H" +#include "fvMesh.H" +#include "volMesh.H" +#include "surfaceMesh.H" +#include "volFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFvGradScheme(barthJespersenGrad) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +template<> +tmp +< + BlockLduSystem::type> +> +barthJespersenGrad::fvmGrad +( + const volVectorField& vf +) const +{ + FatalErrorIn + ( + "tmp barthJespersenGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with cell limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.C new file mode 100644 index 000000000..98c889e41 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.C @@ -0,0 +1,347 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "fvMesh.H" +#include "volMesh.H" +#include "surfaceMesh.H" +#include "GeometricField.H" +#include "zeroGradientFvPatchField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +tmp > +LimitedGrad::limiter +( + const GeoFieldType& vf, + const GeoGradFieldType& gradVf +) const +{ + // Get reference to the mesh + const fvMesh& mesh = vf.mesh(); + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + // Calculate min/max of field + + Field maxVf(vf.internalField()); + Field minVf(vf.internalField()); + + const Field& vfIn = vf.internalField(); + + // Internal faces + forAll (owner, faceI) + { + label own = owner[faceI]; + label nei = neighbour[faceI]; + + const Type& vfOwn = vfIn[own]; + const Type& vfNei = vfIn[nei]; + + // Owner side + maxVf[own] = max(maxVf[own], vfNei); + minVf[own] = min(minVf[own], vfNei); + + // Neighbour side + maxVf[nei] = max(maxVf[nei], vfOwn); + minVf[nei] = min(minVf[nei], vfOwn); + } + + const GeoBoundaryFieldType& bf = vf.boundaryField(); + + // Boundary faces + forAll (bf, patchI) + { + const fvPatchField& psf = bf[patchI]; + + const unallocLabelList& pOwner = mesh.boundary()[patchI].faceCells(); + + if (psf.coupled()) + { + // For a coupled boundary, use neighbour field + Field psfNei = psf.patchNeighbourField(); + + forAll (pOwner, pFaceI) + { + label own = pOwner[pFaceI]; + Type vfNei = psfNei[pFaceI]; + + maxVf[own] = max(maxVf[own], vfNei); + minVf[own] = min(minVf[own], vfNei); + } + } + else + { + // For regular boundary, use boundary value + forAll (pOwner, pFaceI) + { + label own = pOwner[pFaceI]; + Type vfNei = psf[pFaceI]; + + maxVf[own] = max(maxVf[own], vfNei); + minVf[own] = min(minVf[own], vfNei); + } + } + } + + // Subtract the cell value to get differences + maxVf -= vf; + minVf -= vf; + + + // Create a limiter + tmp tlimitField + ( + new GeoFieldType + ( + IOobject + ( + "limitField(" + vf.name() + ")", + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensioned("zero", dimless, pTraits::zero), + zeroGradientFvPatchField::typeName + ) + ); + GeoFieldType& limitField = tlimitField(); + + const volVectorField& C = mesh.C(); + const surfaceVectorField& Cf = mesh.Cf(); + const scalarField& vols = mesh.V().field(); + + Field& lfIn = limitField.internalField(); + + const GradFieldType& g = gradVf.internalField(); + + // Apply limiter function + GradientLimiter limitFunction; + + forAll (owner, faceI) + { + label own = owner[faceI]; + label nei = neighbour[faceI]; + + // Owner side + limitFunction.limiter + ( + lfIn[own], + vols[own], + maxVf[own], + minVf[own], + (Cf[faceI] - C[own]) & g[own] + ); + + // Neighbour side + limitFunction.limiter + ( + lfIn[nei], + vols[nei], + maxVf[nei], + minVf[nei], + (Cf[faceI] - C[nei]) & g[nei] + ); + } + + forAll (bf, patchI) + { + const unallocLabelList& pOwner = mesh.boundary()[patchI].faceCells(); + const vectorField& pCf = Cf.boundaryField()[patchI]; + + forAll (pOwner, pFaceI) + { + label own = pOwner[pFaceI]; + + limitFunction.limiter + ( + lfIn[own], + vols[own], + maxVf[own], + minVf[own], + (pCf[pFaceI] - C[own]) & g[own] + ); + } + } + + // Update coupled boundaries for patchNeighbourField + limitField.boundaryField().evaluateCoupled(); + + if (fv::debug) + { + Info<< "gradient limiter for: " << vf.name() + << " max = " << gMax(lfIn) + << " min = " << gMin(lfIn) + << " average: " << gAverage(lfIn) << endl; + } + + return tlimitField; +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +tmp +< + GeometricField + < + typename outerProduct::type, fvPatchField, volMesh + > +> +LimitedGrad::gradientField +( + const GeoFieldType& vf, + const word& name +) const +{ + // Get base gradient + tmp tGrad = basicGradScheme_().grad(vf, name); + GeoGradFieldType& gradVf = tGrad(); + + // Apply the limiter + GeoFieldType limitField = this->limiter(vf, gradVf); + + gradVf.internalField() = + scaleRow(gradVf.internalField(), limitField.internalField()); + + gradVf.correctBoundaryConditions(); + gradScheme::correctBoundaryConditions(vf, gradVf); + + return tGrad; +} + + +template +tmp +< + BlockLduSystem::type> +> +LimitedGrad::gradientMatrix +( + const GeoFieldType& vf +) const +{ + // Calculate base gradient matrix + tmp tbs = basicGradScheme_().fvmGrad(vf); + GradMatrixType& bs = tbs(); + + // Calculate limiter. Using explicit gradient + GeoFieldType limitField = this->limiter + ( + vf, + basicGradScheme_().grad(vf)() + ); + const Field& lfIn = limitField.internalField(); + + typedef typename CoeffField::linearTypeField + linearCoeffType; + + typedef typename outerProduct::type sourceType; + + Field& source = bs.source(); + + // Grab ldu parts of block matrix as linear always + linearCoeffType& d = bs.diag().asLinear(); + + linearCoeffType& u = bs.upper().asLinear(); + + linearCoeffType& l = bs.lower().asLinear(); + + // Limit upper and lower coeffs + + const fvMesh& mesh = vf.mesh(); + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + forAll(u, faceI) + { + u[faceI] = scaleRow(u[faceI], lfIn[owner[faceI]]); + l[faceI] = scaleRow(l[faceI], lfIn[neighbour[faceI]]); + } + + // Limit diag and source coeffs + forAll(d, cellI) + { + d[cellI] = scaleRow(d[cellI], lfIn[cellI]); + + source[cellI] = scaleRow(source[cellI], lfIn[cellI]); + } + + // Limit coupling coeffs + forAll(vf.boundaryField(), patchI) + { + const fvPatchField& pf = vf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + + const labelList& fc = patch.faceCells(); + + if (patch.coupled()) + { + linearCoeffType& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + + linearCoeffType& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const Field lfNei = + limitField.boundaryField()[patchI].patchNeighbourField(); + + forAll(pf, faceI) + { + pcoupleUpper[faceI] = + scaleRow(pcoupleUpper[faceI], lfIn[fc[faceI]]); + + pcoupleLower[faceI] = + scaleRow(pcoupleLower[faceI], lfNei[faceI]); + } + } + } + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.H new file mode 100644 index 000000000..3ea007406 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/limitedGrad/LimitedGrad.H @@ -0,0 +1,150 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::fv::LimitedGrad + +Description + Second-order gradient scheme with a templated limiter. + +SourceFiles + LimitedGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef LimitedGrad_H +#define LimitedGrad_H + +#include "gradScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class LimitedGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class LimitedGrad +{ +public: + + // Public typedefs + typedef typename outerProduct::type OuterProductType; + + typedef Field FieldType; + typedef GeometricField GeoFieldType; + typedef typename GeoFieldType::GeometricBoundaryField GeoBoundaryFieldType; + + typedef Field::type> GradFieldType; + typedef GeometricField + < + typename outerProduct::type, + fvPatchField, + volMesh + > GeoGradFieldType; + typedef BlockLduSystem::type> + GradMatrixType; + + +private: + + // Private Data + + //- Base gradient scheme + tmp > basicGradScheme_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + LimitedGrad(const LimitedGrad&); + + //- Disallow default bitwise assignment + void operator=(const LimitedGrad&); + + + //- Calculate limiter + tmp limiter + ( + const GeoFieldType& vf, + const GeoGradFieldType& gradVf + ) const; + + +public: + + // Constructors + + //- Construct from mesh and schemeData + LimitedGrad(const fvMesh& mesh, Istream& schemeData) + : + basicGradScheme_(fv::gradScheme::New(mesh, schemeData)) + {} + + + // Member Functions + + //- Return the gradient of the given field to the gradScheme::grad + // for optional caching + tmp gradientField + ( + const GeoFieldType& vf, + const word& name + ) const; + + //- Return the BlockLduSystem corresponding to the implicit gradient + // discretization. For block coupled systems. + tmp gradientMatrix + ( + const GeoFieldType& vf + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "LimitedGrad.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrad.H new file mode 100644 index 000000000..0232248fe --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrad.H @@ -0,0 +1,139 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::fv::minModGrad + +Description + MinMod gradient limiter applied to a runTime selected + base gradient scheme. + +SourceFiles + minModGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef minModGrad_H +#define minModGrad_H + +#include "gradScheme.H" +#include "LimitedGrad.H" +#include "MinModLimiter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class minModGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class minModGrad +: + public fv::gradScheme, + public LimitedGrad +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + minModGrad(const minModGrad&); + + //- Disallow default bitwise assignment + void operator=(const minModGrad&); + + +public: + + //- RunTime type information + TypeName("minMod"); + + + // Constructors + + //- Construct from mesh and schemeData + minModGrad(const fvMesh& mesh, Istream& schemeData) + : + gradScheme(mesh), + LimitedGrad(mesh, schemeData) + {} + + + // Member Functions + + //- Return the gradient of the given field to the gradScheme::grad + // for optional caching + virtual tmp + < + GeometricField + ::type, fvPatchField, volMesh> + > calcGrad + ( + const GeometricField& vf, + const word& name + ) const + { + return LimitedGrad::gradientField + ( + vf, + name + ); + } + + //- Return the BlockLduSystem corresponding to the implicit cell + // limited grad discretization. For block coupled systems. + virtual tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& vf + ) const + { + return LimitedGrad::gradientMatrix + ( + vf + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrads.C new file mode 100644 index 000000000..34aa93b73 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/minModGrad/minModGrads.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "minModGrad.H" +#include "fvMesh.H" +#include "volMesh.H" +#include "surfaceMesh.H" +#include "volFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFvGradScheme(minModGrad) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +template<> +tmp +< + BlockLduSystem::type> +> +minModGrad::fvmGrad +( + const volVectorField& vf +) const +{ + FatalErrorIn + ( + "tmp minModGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with cell limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrad.H new file mode 100644 index 000000000..5b57e51ce --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrad.H @@ -0,0 +1,139 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::fv::venkatakrishnanGrad + +Description + Venkatakrishnan gradient limiter applied to a runTime selected + base gradient scheme. + +SourceFiles + venkatakrishnanGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef venkatakrishnanGrad_H +#define venkatakrishnanGrad_H + +#include "gradScheme.H" +#include "LimitedGrad.H" +#include "VenkatakrishnanLimiter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class venkatakrishnanGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class venkatakrishnanGrad +: + public fv::gradScheme, + public LimitedGrad +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + venkatakrishnanGrad(const venkatakrishnanGrad&); + + //- Disallow default bitwise assignment + void operator=(const venkatakrishnanGrad&); + + +public: + + //- RunTime type information + TypeName("venkatakrishnan"); + + + // Constructors + + //- Construct from mesh and schemeData + venkatakrishnanGrad(const fvMesh& mesh, Istream& schemeData) + : + gradScheme(mesh), + LimitedGrad(mesh, schemeData) + {} + + + // Member Functions + + //- Return the gradient of the given field to the gradScheme::grad + // for optional caching + virtual tmp + < + GeometricField + ::type, fvPatchField, volMesh> + > calcGrad + ( + const GeometricField& vf, + const word& name + ) const + { + return LimitedGrad::gradientField + ( + vf, + name + ); + } + + //- Return the BlockLduSystem corresponding to the implicit cell + // limited grad discretization. For block coupled systems. + virtual tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& vf + ) const + { + return LimitedGrad::gradientMatrix + ( + vf + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrads.C new file mode 100644 index 000000000..88863462c --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/venkatakrishnanGrad/venkatakrishnanGrads.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "venkatakrishnanGrad.H" +#include "fvMesh.H" +#include "volMesh.H" +#include "surfaceMesh.H" +#include "volFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFvGradScheme(venkatakrishnanGrad) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +template<> +tmp +< + BlockLduSystem::type> +> +venkatakrishnanGrad::fvmGrad +( + const volVectorField& vf +) const +{ + FatalErrorIn + ( + "tmp venkatakrishnanGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with cell limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrad.H new file mode 100644 index 000000000..3e8b7fa02 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrad.H @@ -0,0 +1,139 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::fv::wangGrad + +Description + Wang gradient limiter applied to a runTime selected + base gradient scheme. + +SourceFiles + wangGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef wangGrad_H +#define wangGrad_H + +#include "gradScheme.H" +#include "LimitedGrad.H" +#include "WangLimiter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class wangGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class wangGrad +: + public fv::gradScheme, + public LimitedGrad +{ + // Private Member Functions + + //- Disallow default bitwise copy construct + wangGrad(const wangGrad&); + + //- Disallow default bitwise assignment + void operator=(const wangGrad&); + + +public: + + //- RunTime type information + TypeName("wang"); + + + // Constructors + + //- Construct from mesh and schemeData + wangGrad(const fvMesh& mesh, Istream& schemeData) + : + gradScheme(mesh), + LimitedGrad(mesh, schemeData) + {} + + + // Member Functions + + //- Return the gradient of the given field to the gradScheme::grad + // for optional caching + virtual tmp + < + GeometricField + ::type, fvPatchField, volMesh> + > calcGrad + ( + const GeometricField& vf, + const word& name + ) const + { + return LimitedGrad::gradientField + ( + vf, + name + ); + } + + //- Return the BlockLduSystem corresponding to the implicit cell + // limited grad discretization. For block coupled systems. + virtual tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& vf + ) const + { + return LimitedGrad::gradientMatrix + ( + vf + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrads.C new file mode 100644 index 000000000..bb49a76b8 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/wangGrad/wangGrads.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "wangGrad.H" +#include "fvMesh.H" +#include "volMesh.H" +#include "surfaceMesh.H" +#include "volFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFvGradScheme(wangGrad) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +template<> +tmp +< + BlockLduSystem::type> +> +wangGrad::fvmGrad +( + const volVectorField& vf +) const +{ + FatalErrorIn + ( + "tmp wangGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with cell limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradientLimiters/BarthJespersenLimiter.H b/src/finiteVolume/finiteVolume/gradientLimiters/BarthJespersenLimiter.H new file mode 100644 index 000000000..c715ba6e9 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradientLimiters/BarthJespersenLimiter.H @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + BarthJespersenLimiter + +Description + Barth-Jespersen limiter + +Author + Aleksandar Jemcov, Hrvoje Jasak + +\*---------------------------------------------------------------------------*/ + +#ifndef BarthJespersenLimiter_H +#define BarthJespersenLimiter_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BarthJespersenLimiter Declaration +\*---------------------------------------------------------------------------*/ + +class BarthJespersenLimiter +{ +public: + + // Constructors + + //- Construct null + BarthJespersenLimiter() + {} + + + // Destructor - default + + + // Member functions + + //- Set scalar limiter value + inline void limiter + ( + scalar& lim, + const scalar& cellVolume, + const scalar& deltaOneMax, + const scalar& deltaOneMin, + const scalar& extrapolate + ) + { + if (mag(extrapolate) < SMALL) + { + lim = 1; + } + else + { + lim = min + ( + max + ( + max(deltaOneMax/extrapolate, 0), + max(deltaOneMin/extrapolate, 0) + ), + 1 + ); + } + } + + //- Set Type limiter + template + inline void limiter + ( + Type& lim, + const scalar& cellVolume, + const Type& deltaOneMax, + const Type& deltaOneMin, + const Type& extrapolate + ) + { + for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) + { + limiter + ( + lim.component(cmpt), + cellVolume, + deltaOneMax.component(cmpt), + deltaOneMin.component(cmpt), + extrapolate.component(cmpt) + ); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradientLimiters/MinModLimiter.H b/src/finiteVolume/finiteVolume/gradientLimiters/MinModLimiter.H new file mode 100644 index 000000000..74def91d5 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradientLimiters/MinModLimiter.H @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + MinModLimiter + +Description + Non-differentiable cell limiter. The limiter limits the gradient such + that the extrapolated face value does not under- or over-shoot the + neighbouring cell values + +Author + Hrvoje Jasak + +\*---------------------------------------------------------------------------*/ + +#ifndef MinModLimiter_H +#define MinModLimiter_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class MinModLimiter Declaration +\*---------------------------------------------------------------------------*/ + +class MinModLimiter +{ +public: + + // Constructors + + //- Construct null + MinModLimiter() + {} + + + // Destructor - default + + + // Member functions + + //- Set scalar limiter value + inline void limiter + ( + scalar& lim, + const scalar& cellVolume, + const scalar& deltaOneMax, + const scalar& deltaOneMin, + const scalar& extrapolate + ) + { + if (mag(extrapolate) < SMALL) + { + lim = 1; + } + else + { + if (extrapolate > deltaOneMax + VSMALL) + { + lim = min(1, deltaOneMax/extrapolate); + } + else if (extrapolate < deltaOneMin - VSMALL) + { + lim = min(1, deltaOneMin/extrapolate); + } + } + } + + //- Set Type limiter + template + inline void limiter + ( + Type& lim, + const scalar& cellVolume, + const Type& deltaOneMax, + const Type& deltaOneMin, + const Type& extrapolate + ) + { + for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) + { + limiter + ( + lim.component(cmpt), + cellVolume, + deltaOneMax.component(cmpt), + deltaOneMin.component(cmpt), + extrapolate.component(cmpt) + ); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradientLimiters/VenkatakrishnanLimiter.H b/src/finiteVolume/finiteVolume/gradientLimiters/VenkatakrishnanLimiter.H new file mode 100644 index 000000000..1b6516753 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradientLimiters/VenkatakrishnanLimiter.H @@ -0,0 +1,169 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + VenkatakrishnanLimiter + +Description + Venkatakrishnan differentiable limiter + +Author + Aleksandar Jemcov + Updated by Hrvoje Jasak + +\*---------------------------------------------------------------------------*/ + +#ifndef VenkatakrishnanLimiter_H +#define VenkatakrishnanLimiter_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class VenkatakrishnanLimiter Declaration +\*---------------------------------------------------------------------------*/ + +class VenkatakrishnanLimiter +{ + // Private data + + //- Limiter value + scalar k_; + + +public: + + // Constructor + + //- Construct null + VenkatakrishnanLimiter() + : + k_(5) + {} + + //- Set scalar limiter value + inline void limiter + ( + scalar& lim, + const scalar& cellVolume, + const scalar& deltaOneMax, + const scalar& deltaOneMin, + const scalar& extrapolate + ) + { + scalar epsilonSquare = pow3(k_)*cellVolume; + + if (extrapolate > SMALL) + { + lim = max + ( + min + ( + ( + (sqr(deltaOneMax) + epsilonSquare)*extrapolate + + 2*sqr(extrapolate)*deltaOneMax + )/ + stabilise + ( + extrapolate* + ( + sqr(deltaOneMax) + + 2.0*sqr(extrapolate) + + deltaOneMax*extrapolate + + epsilonSquare + ), + SMALL + ), + 1 + ), + 0 + ); + } + else if (extrapolate < SMALL) + { + lim = max + ( + min + ( + ( + (sqr(deltaOneMin) + epsilonSquare)*extrapolate + + 2*sqr(extrapolate)*deltaOneMin + )/ + stabilise + ( + extrapolate* + ( + sqr(deltaOneMin) + + 2.0*sqr(extrapolate) + deltaOneMin*extrapolate + + epsilonSquare + ), + SMALL + ), + 1 + ), + 0 + ); + } + else + { + lim = 1; + } + } + + //- Set Type limiter + template + inline void limiter + ( + Type& lim, + const scalar& cellVolume, + const Type& deltaOneMax, + const Type& deltaOneMin, + const Type& extrapolate + ) + { + for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) + { + limiter + ( + lim.component(cmpt), + cellVolume, + deltaOneMax.component(cmpt), + deltaOneMin.component(cmpt), + extrapolate.component(cmpt) + ); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradientLimiters/WangLimiter.H b/src/finiteVolume/finiteVolume/gradientLimiters/WangLimiter.H new file mode 100644 index 000000000..a0b6c81f6 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradientLimiters/WangLimiter.H @@ -0,0 +1,170 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + WangLimiter + +Description + Wang differentiable limiter + +Author + Aleksandar Jemcov + Updated by Hrvoje Jasak + +\*---------------------------------------------------------------------------*/ + +#ifndef WangLimiter_H +#define WangLimiter_H + +#include "vector.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class WangLimiter Declaration +\*---------------------------------------------------------------------------*/ + +class WangLimiter +{ + // Private data + + //- Limiter value + scalar epsilonPrime_; + + +public: + + // Constructor + + //- Construct null + WangLimiter() + : + epsilonPrime_(0.05) + {} + + //- Set scalar limiter value + inline void limiter + ( + scalar& lim, + const scalar& cellVolume, + const scalar& deltaOneMax, + const scalar& deltaOneMin, + const scalar& extrapolate + ) + { + scalar epsilonSquare = + sqr(epsilonPrime_*(deltaOneMax - deltaOneMin)); + + if (extrapolate > SMALL) + { + lim = max + ( + min + ( + ( + (sqr(deltaOneMax) + epsilonSquare)*extrapolate + + 2*sqr(extrapolate)*deltaOneMax + )/ + stabilise + ( + extrapolate* + ( + sqr(deltaOneMax) + + 2.0*sqr(extrapolate) + + deltaOneMax*extrapolate + + epsilonSquare + ), + SMALL + ), + 1 + ), + 0 + ); + } + else if (extrapolate < SMALL) + { + lim = max + ( + min + ( + ( + (sqr(deltaOneMin) + epsilonSquare)*extrapolate + + 2*sqr(extrapolate)*deltaOneMin + )/ + stabilise + ( + extrapolate* + ( + sqr(deltaOneMin) + + 2.0*sqr(extrapolate) + deltaOneMin*extrapolate + + epsilonSquare + ), + SMALL + ), + 1 + ), + 0 + ); + } + else + { + lim = 1; + } + } + + //- Set Type limiter + template + inline void limiter + ( + Type& lim, + const scalar& cellVolume, + const Type& deltaOneMax, + const Type& deltaOneMin, + const Type& extrapolate + ) + { + for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++) + { + limiter + ( + lim.component(cmpt), + cellVolume, + deltaOneMax.component(cmpt), + deltaOneMin.component(cmpt), + extrapolate.component(cmpt) + ); + } + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +#endif + +// ************************************************************************* // From 4c7473ddf447086a3ecc34ed6c29b11e7143ca08 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 10:57:32 +0000 Subject: [PATCH 010/150] Bugfix: wrong directory --- tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun b/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun index 22a9194a1..49182bed1 100755 --- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun +++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/Allrun @@ -19,7 +19,7 @@ runApplication simpleFoam # Copy the mesh from the steady state case and map the results to a # mesh motion case, then solve transient. cd ../wingMotion2D_pimpleDyMFoam -\cp -r ../wingMotion2D_simpleFoam/constant/polyMesh constant/polyMesh +\cp -r ../wingMotion2D_simpleFoam/constant/polyMesh/* constant/polyMesh \cp -r 0.org 0 runApplication mapFields ../wingMotion2D_simpleFoam -sourceTime latestTime -consistent mv 0/pointDisplacement.unmapped 0/pointDisplacement From ec1dfcee2f8b6397f4fc55cf573bdfa5b7ea0e9b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 16:31:51 +0000 Subject: [PATCH 011/150] Bugfix: specified motion and 6-DOF constraints --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 218 +++++++++--- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H | 17 + src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H | 8 + .../meshMotion/solidBodyMotion/Make/files | 7 + .../meshMotion/solidBodyMotion/SDA/SDA.C | 8 +- .../meshMotion/solidBodyMotion/SKA/SKA.C | 8 +- .../constantVelocity/constantVelocity.C | 133 ++++++++ .../constantVelocity/constantVelocity.H | 134 ++++++++ .../solidBodyMotion/graphMotion/graphMotion.C | 281 ++++++++++++++++ .../solidBodyMotion/graphMotion/graphMotion.H | 144 ++++++++ .../graphVelocity/graphVelocity.C | 317 ++++++++++++++++++ .../graphVelocity/graphVelocity.H | 168 ++++++++++ .../harmonicOscillation/harmonicOscillation.C | 216 ++++++++++++ .../harmonicOscillation/harmonicOscillation.H | 151 +++++++++ .../linearOscillation/linearOscillation.C | 2 +- .../rotatingOscillation/rotatingOscillation.C | 15 +- .../rotatingOscillation/rotatingOscillation.H | 43 ++- 17 files changed, 1816 insertions(+), 54 deletions(-) create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.H create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.H create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.H create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.H diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index 62757ffb1..fdfb7bc55 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -74,11 +74,18 @@ Foam::dimensionedVector Foam::sixDOFqODE::A const HamiltonRodriguezRot& rotation ) const { + // Fix the global force for global rotation constraints + dimensionedVector fAbs = force(); + vector& fAbsVal = fAbs.value(); + + // Constrain rotation + fAbsVal = constrainTranslation(fAbsVal); + return ( - (linSpringCoeffs_ & xR) // spring - (linDampingCoeffs_ & uR) // damping - + force() + + fAbs // To absolute + (rotation.invR() & forceRelative()) )/mass_; @@ -95,18 +102,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot dimensionedVector mAbs = moment(); vector& mAbsVal = mAbs.value(); - if (fixedRoll_) - { - mAbsVal.x() = 0; - } - if (fixedPitch_) - { - mAbsVal.y() = 0; - } - if (fixedYaw_) - { - mAbsVal.z() = 0; - } + mAbsVal = constrainRotation(mAbsVal); return inv(momentOfInertia_) @@ -128,6 +124,97 @@ Foam::dimensionedVector Foam::sixDOFqODE::E } +Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const +{ + vector consVec(vector::zero); + + // Constrain the vector in respect to referent or global coordinate system + if (referentMotionConstraints_) + { + consVec = referentRotation_.R() & vec; + + if (fixedRoll_) + { + consVec.x() = 0; + } + if (fixedPitch_) + { + consVec.y() = 0; + } + if (fixedYaw_) + { + consVec.z() = 0; + } + + consVec = referentRotation_.invR() & consVec; + } + else + { + consVec = vec; + + if (fixedRoll_) + { + consVec.x() = 0; + } + if (fixedPitch_) + { + consVec.y() = 0; + } + if (fixedYaw_) + { + consVec.z() = 0; + } + } + + return consVec; +} + + +Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const +{ + vector consVec(vector::zero); + + // Constrain the vector in respect to referent or global coordinate system + if (referentMotionConstraints_) + { + consVec = referentRotation_.R() & vec; + + if (fixedSurge_) + { + consVec.x() = 0; + } + if (fixedSway_) + { + consVec.y() = 0; + } + if (fixedHeave_) + { + consVec.z() = 0; + } + + consVec = referentRotation_.invR() & consVec; + } + else + { + consVec = vec; + + if (fixedSurge_) + { + consVec.x() = 0; + } + if (fixedSway_) + { + consVec.y() = 0; + } + if (fixedHeave_) + { + consVec.z() = 0; + } + } + + return consVec; +} + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -167,7 +254,16 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) fixedHeave_(lookup("fixedHeave")), fixedRoll_(lookup("fixedRoll")), fixedPitch_(lookup("fixedPitch")), - fixedYaw_(lookup("fixedYaw")) + fixedYaw_(lookup("fixedYaw")), + referentMotionConstraints_ + ( + lookupOrDefault + ( + "referentMotionConstraints", + false + ) + ), + referentRotation_(vector::zero, 0) { setCoeffs(); } @@ -224,7 +320,9 @@ Foam::sixDOFqODE::sixDOFqODE fixedHeave_(sd.fixedHeave_), fixedRoll_(sd.fixedRoll_), fixedPitch_(sd.fixedPitch_), - fixedYaw_(sd.fixedYaw_) + fixedYaw_(sd.fixedYaw_), + referentMotionConstraints_(sd.referentMotionConstraints_), + referentRotation_(sd.referentRotation_) {} @@ -299,22 +397,22 @@ void Foam::sixDOFqODE::derivatives dydx[4] = accel.y(); dydx[5] = accel.z(); - // Add translational constraints by setting RHS of given components to zero - if (fixedSurge_) - { - dydx[0] = 0; // Surge velocity - dydx[3] = 0; // Surge acceleration - } - if (fixedSway_) - { - dydx[1] = 0; // Sway velocity - dydx[4] = 0; // Sway acceleration - } - if (fixedHeave_) - { - dydx[2] = 0; // Heave velocity - dydx[5] = 0; // Heave acceleration - } +// // Add translational constraints by setting RHS of given components to zero +// if (fixedSurge_) +// { +// dydx[0] = 0; // Surge velocity +// dydx[3] = 0; // Surge acceleration +// } +// if (fixedSway_) +// { +// dydx[1] = 0; // Sway velocity +// dydx[4] = 0; // Sway acceleration +// } +// if (fixedHeave_) +// { +// dydx[2] = 0; // Heave velocity +// dydx[5] = 0; // Heave acceleration +// } // Set the derivatives for rotation dimensionedVector curOmega @@ -324,6 +422,25 @@ void Foam::sixDOFqODE::derivatives vector(y[6], y[7], y[8]) ); +// dimensionedVector curGlobalOmega = curRotation.invR() & curOmega; +// +// // Add rotational constraints by setting RHS of given components to zero +// if (fixedRoll_) +// { +// curGlobalOmega.value().x() = 0; +// } +// if (fixedPitch_) +// { +// curGlobalOmega.value().y() = 0; +// } +// if (fixedYaw_) +// { +// curGlobalOmega.value().z() = 0; +// } +// +// +// curOmega = curRotation.R() & curGlobalOmega; + const vector omegaDot = OmegaDot(curRotation, curOmega).value(); dydx[6] = omegaDot.x(); @@ -335,19 +452,20 @@ void Foam::sixDOFqODE::derivatives dydx[11] = curRotation.eDot(curOmega.value(), 2); dydx[12] = curRotation.eDot(curOmega.value(), 3); - // Add rotational constraints by setting RHS of given components to zero - if (fixedRoll_) - { - dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) - } - if (fixedPitch_) - { - dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) - } - if (fixedYaw_) - { - dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) - } + +// // Add rotational constraints by setting RHS of given components to zero +// if (fixedRoll_) +// { +// dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) +// } +// if (fixedPitch_) +// { +// dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) +// } +// if (fixedYaw_) +// { +// dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) +// } } @@ -384,6 +502,12 @@ void Foam::sixDOFqODE::update(const scalar delta) Uval.y() = coeffs_[4]; Uval.z() = coeffs_[5]; + // Constrain velocity + Uval = constrainTranslation(Uval); + coeffs_[3] = Uval.x(); + coeffs_[4] = Uval.y(); + coeffs_[5] = Uval.z(); + // Update omega vector& omegaVal = omega_.value(); @@ -391,6 +515,12 @@ void Foam::sixDOFqODE::update(const scalar delta) omegaVal.y() = coeffs_[7]; omegaVal.z() = coeffs_[8]; + // Constrain omega + omegaVal = constrainRotation(omegaVal); + coeffs_[6] = omegaVal.x(); + coeffs_[7] = omegaVal.y(); + coeffs_[8] = omegaVal.z(); + rotation_.updateRotation ( HamiltonRodriguezRot diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H index 6b9eb3d4d..fbe85f95d 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H @@ -151,6 +151,12 @@ class sixDOFqODE //- Fixed yaw (rotation around z) Switch fixedYaw_; + //- Restraints in referent coordinate system + Switch referentMotionConstraints_; + + //- Rotation of referent coordinate system + HamiltonRodriguezRot referentRotation_; + // Private Member Functions @@ -191,6 +197,14 @@ class sixDOFqODE const dimensionedVector& omega ) const; + //- Constrain rotation vector in referent or global coordinate + // system + vector constrainRotation(vector& vec) const; + + //- Constrain translation vector in referent or global coordinate + // system + vector constrainTranslation(vector& vec) const; + public: @@ -286,6 +300,9 @@ public: // coordinate system inline void setOmega(const vector& omega); + //- Set referent coordinate system to apply constraints + inline void setReferentRotation(const HamiltonRodriguezRot& rot); + // Average motion per time-step diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H index 75e878040..2d619a694 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H @@ -166,6 +166,14 @@ void Foam::sixDOFqODE::setOmega(const vector& omega) } +void Foam::sixDOFqODE::setReferentRotation(const HamiltonRodriguezRot& rot) +{ + referentRotation_ = rot; + + referentMotionConstraints_ = true; +} + + const Foam::dimensionedVector& Foam::sixDOFqODE::Uaverage() const { return Uaverage_; diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files index 1e02ad3a4..41f06b6bc 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files @@ -3,9 +3,16 @@ solidBodyMotionFunction/newSolidBodyMotionFunction.C noMotion/noMotion.C translation/translation.C +graphMotion/graphMotion.C + linearOscillation/linearOscillation.C rotatingOscillation/rotatingOscillation.C +harmonicOscillation/harmonicOscillation.C + SDA/SDA.C SKA/SKA.C +constantVelocity/constantVelocity.C +graphVelocity/graphVelocity.C + LIB = $(FOAM_LIBBIN)/libsolidBodyMotion diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/SDA/SDA.C b/src/dynamicMesh/meshMotion/solidBodyMotion/SDA/SDA.C index 49c50fd1d..7a3b35fbf 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/SDA/SDA.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/SDA/SDA.C @@ -116,7 +116,13 @@ Foam::septernion Foam::solidBodyMotionFunctions::SDA::velocity() const scalar t = time_.value(); scalar dt = time_.deltaT().value(); - return (calcTransformation(t + dt) - calcTransformation(t))/dt; + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; } diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/SKA/SKA.C b/src/dynamicMesh/meshMotion/solidBodyMotion/SKA/SKA.C index 3157c4f3a..7fa867869 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/SKA/SKA.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/SKA/SKA.C @@ -129,7 +129,13 @@ Foam::septernion Foam::solidBodyMotionFunctions::SKA::velocity() const scalar t = time_.value(); scalar dt = time_.deltaT().value(); - return (calcTransformation(t + dt) - calcTransformation(t))/dt; + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; } diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C b/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C new file mode 100644 index 000000000..dbb82f7f1 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "constantVelocity.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +defineTypeNameAndDebug(constantVelocity, 0); +addToRunTimeSelectionTable +( + solidBodyMotionFunction, + constantVelocity, + dictionary +); + +} +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::constantVelocity::calcTransformation +( + const scalar t +) const +{ + const vector translation = transVelocity_*(t - startMotionTime_); + const vector rotation = rotVelocity_*(t - startMotionTime_); + + const quaternion R(rotation.x(), rotation.y(), rotation.z()); + const septernion TR + ( + septernion(origin_ + translation)*R*septernion(-origin_) + ); + + return TR; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::constantVelocity::constantVelocity +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime) +{ + read(SBMFCoeffs); +} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::constantVelocity::transformation() const +{ + const scalar t = time_.value(); + + const septernion TR = calcTransformation(t); + + Info<< "solidBodyMotionFunctions::constantVelocity::transformation(): " + << "Time = " << t << " transformation: " << TR << endl; + + return TR; +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::constantVelocity::velocity() const +{ + const scalar t = time_.value(); + const scalar dt = time_.deltaT().value(); + + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; +} + + +bool Foam::solidBodyMotionFunctions::constantVelocity::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + + SBMFCoeffs_.lookup("origin") >> origin_; + SBMFCoeffs_.lookup("translationalVelocity") >> transVelocity_; + SBMFCoeffs_.lookup("rotationalVelocity") >> rotVelocity_; + SBMFCoeffs_.lookup("startMotionTime") >> startMotionTime_; + SBMFCoeffs_.lookup("inDegrees") >> inDegrees_; + + // Convert to radians if necessary + rotVelocity_ *= inDegrees_ ? mathematicalConstant::pi/180.0 : 1; + + return true; +} + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.H b/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.H new file mode 100644 index 000000000..cb18609f3 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.H @@ -0,0 +1,134 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::solidBodyMotionFunction::constantVelocity + +Description + Prescribed constant translational and angular velocity. + +SourceFiles + constantVelocity.C + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef constantVelocity_H +#define constantVelocity_H + +#include "solidBodyMotionFunction.H" +#include "primitiveFields.H" +#include "point.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class constantVelocity Declaration +\*---------------------------------------------------------------------------*/ + +class constantVelocity +: + public solidBodyMotionFunction +{ + // Private data + + //- Centre of gravity + point origin_; + + //- Translational velocity vector + vector transVelocity_; + + //- Rotational velocity vector + vector rotVelocity_; + + //- Start motion time + scalar startMotionTime_; + + //- Is the rotational velocity given in degrees/sec + Switch inDegrees_; + + + // Private Member Functions + + //- Disallow copy construct + constantVelocity(const constantVelocity&); + + //- Disallow default bitwise assignment + void operator=(const constantVelocity&); + + //- Calculate tranformation + septernion calcTransformation(const scalar t) const; + + +public: + + //- Runtime type information + TypeName("constantVelocity"); + + + // Constructors + + //- Construct from components + constantVelocity + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + //- Destructor + virtual ~constantVelocity() + {} + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.C b/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.C new file mode 100644 index 000000000..484f895e1 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.C @@ -0,0 +1,281 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "graphMotion.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" +#include "interpolateXY.H" +#include "IFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +defineTypeNameAndDebug(graphMotion, 0); +addToRunTimeSelectionTable(solidBodyMotionFunction, graphMotion, dictionary); + +} +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::graphMotion::calcTransformation +( + const scalar t +) const +{ + const vector translation + ( + interpolateXY(t, surge_.x(), surge_.y()), + interpolateXY(t, sway_.x(), sway_.y()), + interpolateXY(t, heave_.x(), heave_.y()) + ); + + vector rotation + ( + interpolateXY(t, roll_.x(), roll_.y()), + interpolateXY(t, pitch_.x(), pitch_.y()), + interpolateXY(t, yaw_.x(), yaw_.y()) + ); + + if (inDegrees_) + { + const scalar piBy180 = mathematicalConstant::pi/180.0; + + rotation *= piBy180; + } + + const quaternion R(rotation.x(), rotation.y(), rotation.z()); + const septernion TR + ( + septernion(origin_ + translation)*R*septernion(-origin_) + ); + + return TR; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::graphMotion::graphMotion +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime), + surge_ + ( + "surge", + "t", + "eta1Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("surge")) + )() + ), + sway_ + ( + "sway", + "t", + "eta2Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("sway")) + )() + ), + heave_ + ( + "heave", + "t", + "eta3Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("heave")) + )() + ), + roll_ + ( + "roll", + "t", + "eta4Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("roll")) + )() + ), + pitch_ + ( + "pitch", + "t", + "eta5Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("pitch")) + )() + ), + yaw_ + ( + "yaw", + "t", + "eta6Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("yaw")) + )() + ) +{ + read(SBMFCoeffs); +} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::graphMotion::transformation() const +{ + const scalar t = time_.value(); + const septernion TR = calcTransformation(t); + + Info<< "solidBodyMotionFunctions::graphMotion::transformation(): " + << "Time = " << t << " transformation: " << TR << endl; + + return TR; +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::graphMotion::velocity() const +{ + const scalar t = time_.value(); + const scalar dt = time_.deltaT().value(); + + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; +} + + +bool Foam::solidBodyMotionFunctions::graphMotion::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + + SBMFCoeffs_.lookup("origin") >> origin_; + SBMFCoeffs_.lookup("inDegrees") >> inDegrees_; + + word surge = SBMFCoeffs_.lookup("surge"); + word sway = SBMFCoeffs_.lookup("sway"); + word heave = SBMFCoeffs_.lookup("heave"); + word roll = SBMFCoeffs_.lookup("roll"); + word pitch = SBMFCoeffs_.lookup("pitch"); + word yaw = SBMFCoeffs_.lookup("yaw"); + + surge_ = graph + ( + "surge", + "t", + "eta1Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + surge + )() + ); + sway_ = graph + ( + "sway", + "t", + "eta2Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + sway + )() + ); + heave_ = graph + ( + "heave", + "t", + "eta3Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + heave + )() + ); + roll_ = graph + ( + "roll", + "t", + "eta4Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + roll + )() + ); + pitch_ = graph + ( + "pitch", + "t", + "eta5Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + pitch + )() + ); + yaw_ = graph + ( + "yaw", + "t", + "eta6Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + yaw + )() + ); + + return true; +} + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.H b/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.H new file mode 100644 index 000000000..6bc076f2e --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::solidBodyMotionFunction::graphMotion + +Description + Prescribed translational and rotational motion given graphs for surge, + sway, heave, roll, pitch and yaw. Hence, the motion (parts of trajectory) is + given in graphs, not the velocity. + +SourceFiles + graphMotion.C + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef graphMotion_H +#define graphMotion_H + +#include "solidBodyMotionFunction.H" +#include "graph.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class graphMotion Declaration +\*---------------------------------------------------------------------------*/ + +class graphMotion +: + public solidBodyMotionFunction +{ + // Private data + + //- Centre of gravity + point origin_; + + + // Prescribed translation members + + // - Graph for surge translational motion (x-translation) + graph surge_; + + // - Graph for sway translational motion (y-translation) + graph sway_; + + // - Graph for heave translational motion (z-translation) + graph heave_; + + + // Prescribed rotation members + + // - Graph for roll rotational motion (x-rotation) + graph roll_; + + // - Graph for pitch rotational motion (y-rotation) + graph pitch_; + + // - Graph for yaw rotational motion (z-rotation) + graph yaw_; + + //- Is the rotational velocity given in degrees/sec + Switch inDegrees_; + + + //- Member functions + + //- Calculate transformation septernion + septernion calcTransformation(const scalar t) const; + + +public: + + //- Runtime type information + TypeName("graphMotion"); + + + // Constructors + + //- Construct from components + graphMotion + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + //- Destructor + virtual ~graphMotion() + {} + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.C b/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.C new file mode 100644 index 000000000..1c6b1e1c6 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.C @@ -0,0 +1,317 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "graphVelocity.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" +#include "interpolateXY.H" +#include "IFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +defineTypeNameAndDebug(graphVelocity, 0); +addToRunTimeSelectionTable(solidBodyMotionFunction, graphVelocity, dictionary); + +} +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::vector +Foam::solidBodyMotionFunctions::graphVelocity::translationalVelocity() const +{ + const scalar t = time_.value() - time_.deltaT().value()/2; + + return vector + ( + interpolateXY(t, surge_.x(), surge_.y()), + interpolateXY(t, sway_.x(), sway_.y()), + interpolateXY(t, heave_.x(), heave_.y()) + ); +} + + +Foam::vector +Foam::solidBodyMotionFunctions::graphVelocity::rotationalVelocity() const +{ + const scalar t = time_.value() - time_.deltaT().value()/2; + + scalar rollVelocity = interpolateXY(t, roll_.x(), roll_.y()); + scalar pitchVelocity = interpolateXY(t, pitch_.x(), pitch_.y()); + scalar yawVelocity = interpolateXY(t, yaw_.x(), yaw_.y()); + + if (inDegrees_) + { + const scalar piBy180 = mathematicalConstant::pi/180.0; + + rollVelocity *= piBy180; + pitchVelocity *= piBy180; + yawVelocity *= piBy180; + } + + return vector(rollVelocity, pitchVelocity, yawVelocity); +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::graphVelocity::calcTransformation() const +{ + // Integrate velocity to get position + if(localTimeIndex_ != time_.timeIndex()) + { + // Set old translation and rotation to the previous state + translationOld_ = translation_; + rotationOld_ = rotation_; + + const scalar dt = time_.deltaT().value(); + + translation_ += translationalVelocity()*dt; + rotation_ += rotationalVelocity()*dt; + + localTimeIndex_ = time_.timeIndex(); + } + + const quaternion R(rotation_.x(), rotation_.y(), rotation_.z()); + const septernion TR + ( + septernion(origin_ + translation_)*R*septernion(-origin_) + ); + + return TR; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::graphVelocity::graphVelocity +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime), + translation_(vector::zero), + rotation_(vector::zero), + translationOld_(vector::zero), + rotationOld_(vector::zero), + localTimeIndex_(-1), + surge_ + ( + "surge", + "t", + "eta1Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("surge")) + )() + ), + sway_ + ( + "sway", + "t", + "eta2Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("sway")) + )() + ), + heave_ + ( + "heave", + "t", + "eta3Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("heave")) + )() + ), + roll_ + ( + "roll", + "t", + "eta4Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("roll")) + )() + ), + pitch_ + ( + "pitch", + "t", + "eta5Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("pitch")) + )() + ), + yaw_ + ( + "yaw", + "t", + "eta6Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + word(SBMFCoeffs_.lookup("yaw")) + )() + ) +{ + read(SBMFCoeffs); +} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::graphVelocity::transformation() const +{ + const septernion TR = calcTransformation(); + + Info<< "solidBodyMotionFunctions::graphVelocity::transformation(): " + << "Time = " << time_.value() << " transformation: " << TR << endl; + + return TR; +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::graphVelocity::velocity() const +{ + const scalar dt = time_.deltaT().value(); + + const septernion TR = calcTransformation(); + + const quaternion ROld(rotationOld_.x(), rotationOld_.y(), rotationOld_.z()); + const septernion TROld + ( + septernion(origin_ + translationOld_)*ROld*septernion(-origin_) + ); + + return septernion + ( + (TR.t() - TROld.t())/dt, + (TR.r()/TROld.r())/dt + ); +} + + +bool Foam::solidBodyMotionFunctions::graphVelocity::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + + SBMFCoeffs_.lookup("origin") >> origin_; + SBMFCoeffs_.lookup("inDegrees") >> inDegrees_; + + word surge = SBMFCoeffs_.lookup("surge"); + word sway = SBMFCoeffs_.lookup("sway"); + word heave = SBMFCoeffs_.lookup("heave"); + word roll = SBMFCoeffs_.lookup("roll"); + word pitch = SBMFCoeffs_.lookup("pitch"); + word yaw = SBMFCoeffs_.lookup("yaw"); + + surge_ = graph + ( + "surge", + "t", + "eta1Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + surge + )() + ); + sway_ = graph + ( + "sway", + "t", + "eta2Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + sway + )() + ); + heave_ = graph + ( + "heave", + "t", + "eta3Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + heave + )() + ); + roll_ = graph + ( + "roll", + "t", + "eta4Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + roll + )() + ); + pitch_ = graph + ( + "pitch", + "t", + "eta5Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + pitch + )() + ); + yaw_ = graph + ( + "yaw", + "t", + "eta6Dot", + IFstream + ( + time_.path()/time_.caseConstant()/ + yaw + )() + ); + + return true; +} + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.H b/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.H new file mode 100644 index 000000000..b05e331c6 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.H @@ -0,0 +1,168 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::solidBodyMotionFunction::graphVelocity + +Description + Prescribed translational and rotational velocity given graphs for surge, + sway, heave, roll, pitch and yaw. Hence, the velocity is given in graphs, + not the motion (trajectory). + +SourceFiles + graphVelocity.C + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef graphVelocity_H +#define graphVelocity_H + +#include "solidBodyMotionFunction.H" +#include "graph.H" +#include "primitiveFields.H" +#include "point.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class graphVelocity Declaration +\*---------------------------------------------------------------------------*/ + +class graphVelocity +: + public solidBodyMotionFunction +{ + // Private data + + //- Centre of gravity + point origin_; + + //- Integrated translation + mutable vector translation_; + + //- Integrated rotation + mutable vector rotation_; + + //- Integrated old translation + mutable vector translationOld_; + + //- Integrated old rotation + mutable vector rotationOld_; + + //- Time index, avoid integrating more than once per time step + mutable label localTimeIndex_; + + + // Prescribed translation members + + // - Graph for surge translational velocity (x-translation) + graph surge_; + + // - Graph for sway translational velocity (y-translation) + graph sway_; + + // - Graph for heave translational motion (z-translation) + graph heave_; + + + // Prescribed rotation members + + // - Graph for roll rotational velocity (x-rotation) + graph roll_; + + // - Graph for pitch rotational velocity (y-rotation) + graph pitch_; + + // - Graph for yaw rotational velocity (z-rotation) + graph yaw_; + + //- Is the rotational velocity given in degrees/sec + Switch inDegrees_; + + + //- Member functions + + //- Return translational velocity vector + vector translationalVelocity() const; + + // -Return rotational velocity vector + vector rotationalVelocity() const; + + //- Calculate transformation septernion + septernion calcTransformation() const; + + +public: + + //- Runtime type information + TypeName("graphVelocity"); + + + // Constructors + + //- Construct from components + graphVelocity + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + //- Destructor + virtual ~graphVelocity() + {} + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.C b/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.C new file mode 100644 index 000000000..df90bb122 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.C @@ -0,0 +1,216 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "harmonicOscillation.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +defineTypeNameAndDebug(harmonicOscillation, 0); +addToRunTimeSelectionTable +( + solidBodyMotionFunction, + harmonicOscillation, + dictionary +); + +} +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::harmonicOscillation::calcTransformation +( + const scalar t +) const +{ + const vector translation = + cmptMultiply + ( + transAmplitudes_, + vector + ( + sin(transAngularFreq_.x()*t + transPhaseAngles_.x()), + sin(transAngularFreq_.y()*t + transPhaseAngles_.y()), + sin(transAngularFreq_.z()*t + transPhaseAngles_.z()) + ) + ); + + + const vector rotation = + cmptMultiply + ( + rotAmplitudes_, + vector + ( + sin(rotAngularFreq_.x()*t + rotPhaseAngles_.x()), + sin(rotAngularFreq_.y()*t + rotPhaseAngles_.y()), + sin(rotAngularFreq_.z()*t + rotPhaseAngles_.z()) + ) + ); + + const quaternion R(rotation.x(), rotation.y(), rotation.z()); + const septernion TR + ( + septernion(origin_ + translation)*R*septernion(-origin_) + ); + + return TR; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::harmonicOscillation::harmonicOscillation +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime) +{ + read(SBMFCoeffs); +} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::harmonicOscillation::transformation() const +{ + const scalar t = time_.value(); + + const septernion TR = calcTransformation(t); + + Info<< "solidBodyMotionFunctions::harmonicOscillation::transformation(): " + << "Time = " << t << " transformation: " << TR << endl; + + return TR; +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::harmonicOscillation::velocity() const +{ + const scalar t = time_.value(); + const scalar dt = time_.deltaT().value(); + + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; +} + + +bool Foam::solidBodyMotionFunctions::harmonicOscillation::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + + SBMFCoeffs_.lookup("origin") >> origin_; + SBMFCoeffs_.lookup("translationalAmplitudes") >> transAmplitudes_; + SBMFCoeffs_.lookup("translationalAngularFrequencies") >> transAngularFreq_; + SBMFCoeffs_.lookup("translationalPhaseAngles") >> transPhaseAngles_; + SBMFCoeffs_.lookup("rotationalAmplitudes") >> rotAmplitudes_; + SBMFCoeffs_.lookup("rotationalAngularFrequencies") >> rotAngularFreq_; + SBMFCoeffs_.lookup("rotationalPhaseAngles") >> rotPhaseAngles_; + SBMFCoeffs_.lookup("inDegrees") >> inDegrees_; + + // Sanity check for negative frequencies + if + ( + transAngularFreq_.x() < -SMALL + || transAngularFreq_.y() < -SMALL + || transAngularFreq_.z() < -SMALL + || rotAngularFreq_.x() < -SMALL + || rotAngularFreq_.y() < -SMALL + || rotAngularFreq_.z() < -SMALL + ) + { + FatalErrorIn + ( + "harmonicOscillation::read" + "(\n" + " const fvMesh& mesh\n" + " const dictionary& dict\n" + ")" + ) << "Negative angular frequency detected." + << abort(FatalError); + } + + // Convert to radians if necessary, printing data + if (inDegrees_) + { + const scalar piBy180 = mathematicalConstant::pi/180.0; + + transPhaseAngles_ *= piBy180; + rotAmplitudes_ *= piBy180; + rotPhaseAngles_ *= piBy180; + + Info<< "Data in degrees converted:" << nl + << "translationalPhaseAngles = " << transPhaseAngles_ << nl + << "rotationalAmplitudes = " << rotAmplitudes_ << nl + << "rotationalPhaseAngles = " << rotPhaseAngles_ << endl; + } + + // Count prescribed rotations: only two harmonic rotations can be prescribed + label nRot = 0; + + // Rotation + mag(rotAmplitudes_.x()) > SMALL ? ++nRot : /* null expression */ 0 ; + mag(rotAmplitudes_.y()) > SMALL ? ++nRot : /* null expression */ 0 ; + mag(rotAmplitudes_.z()) > SMALL ? ++nRot : /* null expression */ 0 ; + + // Check if more than two rotations prescribed + if (nRot > 2) + { + FatalErrorIn + ( + "harmonicOscillation::read" + "(\n" + " const fvMesh& mesh\n" + " const dictionary& dict\n" + ")" + ) << "Only two harmonic rotations can be prescribed. Prescribing " + << "more than two yields unphysical results. " + << abort(FatalError); + } + + return true; +} + + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.H b/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.H new file mode 100644 index 000000000..8c1109e7a --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.H @@ -0,0 +1,151 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::solidBodyMotionFunctions::harmonicOscillation + +Description + Prescribed translational and rotational harmonic motion. + The motion is of the following form: x = x_a*sin(omega*t + phi_x). + Hence, velocities go with cosine instead of sine. + +SourceFiles + harmonicOscillation.C + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef harmonicOscillation_H +#define harmonicOscillation_H + +#include "solidBodyMotionFunction.H" +#include "primitiveFields.H" +#include "point.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class harmonicOscillation Declaration +\*---------------------------------------------------------------------------*/ + +class harmonicOscillation +: + public solidBodyMotionFunction +{ + // Private data + + //- Centre of gravity + point origin_; + + + // Prescribed translation members + + //- Translational amplitudes + vector transAmplitudes_; + + //- Translational angular frequencies + vector transAngularFreq_; + + //- Translational phase angles + vector transPhaseAngles_; + + + // Prescribed rotation members + + //- Rotational amplitudes + vector rotAmplitudes_; + + //- Rotational angular frequencies + vector rotAngularFreq_; + + //- Rotational phase angles + vector rotPhaseAngles_; + + //- Is the rotational velocity given in degrees/sec + Switch inDegrees_; + + + // Private Member Functions + + //- Disallow copy construct + harmonicOscillation(const harmonicOscillation&); + + //- Disallow default bitwise assignment + void operator=(const harmonicOscillation&); + + //- Calculate tranformation + septernion calcTransformation(const scalar t) const; + + +public: + + //- Runtime type information + TypeName("harmonicOscillation"); + + + // Constructors + + //- Construct from components + harmonicOscillation + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + //- Destructor + virtual ~harmonicOscillation() + {} + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/linearOscillation/linearOscillation.C b/src/dynamicMesh/meshMotion/solidBodyMotion/linearOscillation/linearOscillation.C index 9762960ec..e4345a8f9 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/linearOscillation/linearOscillation.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/linearOscillation/linearOscillation.C @@ -102,7 +102,7 @@ Foam::solidBodyMotionFunctions::linearOscillation::velocity() const septernion TV ( - (calcPosition(t + dt) - calcPosition(t))/dt, + (calcPosition(t) - calcPosition(t - dt))/dt, quaternion::zero ); diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.C b/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.C index f221dd2a3..3f711cbf7 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.C @@ -57,7 +57,7 @@ Foam::solidBodyMotionFunctions::rotatingOscillation::calcPosition } -Foam::septernion +Foam::septernion Foam::solidBodyMotionFunctions::rotatingOscillation::calcTransformation ( const scalar t @@ -114,15 +114,23 @@ transformation() const return TR; } + Foam::septernion Foam::solidBodyMotionFunctions::rotatingOscillation::velocity() const { scalar t = time_.value(); scalar dt = time_.deltaT().value(); - - return (calcTransformation(t + dt) - calcTransformation(t))/dt; + + const septernion velocity + ( + (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt, + (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt + ); + + return velocity; } + bool Foam::solidBodyMotionFunctions::rotatingOscillation::read ( const dictionary& SBMFCoeffs @@ -137,4 +145,5 @@ bool Foam::solidBodyMotionFunctions::rotatingOscillation::read return true; } + // ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H b/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H index 80ac1c932..16decb950 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H @@ -1,3 +1,37 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Class + Foam::solidBodyMotionFunctions::rotatingOscillation + +Description + Rotating oscillation motion function + +SourceFiles + rotatingOscillation.C + +\*---------------------------------------------------------------------------*/ + #ifndef rotatingOscillation_H #define rotatingOscillation_H @@ -39,11 +73,12 @@ class rotatingOscillation //- Disallow default bitwise assignment void operator=(const rotatingOscillation&); - - //- Calculate position + + + //- Calculate position vector calcPosition(const scalar t) const; - //- Calculate tranformation + //- Calculate tranformation septernion calcTransformation(const scalar t) const; @@ -84,7 +119,7 @@ public: //- Return the solid-body motion transformation septernion virtual septernion transformation() const; - + //- Return the solid-body motion velocity virtual septernion velocity() const; From b3a4b2fad1911167144bc4faf7e1b03fd5fff4fd Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 16:33:03 +0000 Subject: [PATCH 012/150] Formatting --- .../BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H index 1a3efe91a..a48b9c495 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H @@ -87,7 +87,7 @@ class BlockGaussSeidelPrecon //- Calculate inverse diagonal void calcInvDiag(); - // Block Gauss-Seidel sweep, symetric matrix + //- Block Gauss-Seidel sweep, symetric matrix template void BlockSweep ( @@ -97,7 +97,7 @@ class BlockGaussSeidelPrecon const Field& b ) const; - // Block Gauss-Seidel sweep, asymmetric matrix + //- Block Gauss-Seidel sweep, asymmetric matrix template void BlockSweep ( From 76ae588d3b18b5109006fc0424df64f7425a7081 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 12:54:21 +0000 Subject: [PATCH 013/150] Rename small and large in hinv --- src/foam/primitives/Tensor/tensor/tensor.C | 24 +++++++++++----------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/foam/primitives/Tensor/tensor/tensor.C b/src/foam/primitives/Tensor/tensor/tensor.C index 57906b56a..0ce95c46a 100644 --- a/src/foam/primitives/Tensor/tensor/tensor.C +++ b/src/foam/primitives/Tensor/tensor/tensor.C @@ -607,10 +607,10 @@ tensor eigenVectors(const symmTensor& t) // Matrix inversion with singular value decomposition tensor hinv(const tensor& t) { - static const scalar large = 1e10; - static const scalar small = 1e-10; + static const scalar hinvLarge = 1e10; + static const scalar hinvSmall = 1e-10; - if (det(t) > small) + if (det(t) > hinvSmall) { return inv(t); } @@ -627,7 +627,7 @@ tensor hinv(const tensor& t) // Jovani Favero, 18/Nov/2009 // Further bug fix: replace > with == and add SMALL to zeroInv // Dominik Christ, 7/Aug/2012 - if (mag(eig.z()) == large*mag(eig.z())) + if (mag(eig.z()) == hinvLarge*mag(eig.z())) { // Three zero eigen values (z is largest in magnitude). // Return zero inverse @@ -637,13 +637,13 @@ tensor hinv(const tensor& t) // Compare smaller eigen values and if out of range of large // consider them singular - if (mag(eig.z()) > large*mag(eig.x())) + if (mag(eig.z()) > hinvLarge*mag(eig.x())) { // Make a tensor out of symmTensor sqr. HJ, 24/Oct/2009 zeroInv += tensor(sqr(eigVecs.x())); } - if (mag(eig.z()) > large*mag(eig.y())) + if (mag(eig.z()) > hinvLarge*mag(eig.y())) { // Make a tensor out of symmTensor sqr. HJ, 24/Oct/2009 zeroInv += tensor(sqr(eigVecs.y())); @@ -656,10 +656,10 @@ tensor hinv(const tensor& t) symmTensor hinv(const symmTensor& t) { - static const scalar large = 1e10; - static const scalar small = 1e-10; + static const scalar hinvLarge = 1e10; + static const scalar hinvSmall = 1e-10; - if (det(t) > small) + if (det(t) > hinvSmall) { return inv(t); } @@ -676,7 +676,7 @@ symmTensor hinv(const symmTensor& t) // Jovani Favero, 18/Nov/2009 // Further bug fix: replace > with == and add SMALL to zeroInv // Dominik Christ, 7/Aug/2012 - if (mag(eig.z()) == large*mag(eig.z())) + if (mag(eig.z()) == hinvLarge*mag(eig.z())) { // Three zero eigen values (z is largest in magnitude). // Return zero inverse @@ -686,12 +686,12 @@ symmTensor hinv(const symmTensor& t) // Compare smaller eigen values and if out of range of large // consider them singular - if (mag(eig.z()) > large*mag(eig.x())) + if (mag(eig.z()) > hinvLarge*mag(eig.x())) { zeroInv += sqr(eigVecs.x()); } - if (mag(eig.z()) > large*mag(eig.y())) + if (mag(eig.z()) > hinvLarge*mag(eig.y())) { zeroInv += sqr(eigVecs.y()); } From 101e59c0e8cba1d9203df2e123ad7588ae759e8a Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 13:11:19 +0000 Subject: [PATCH 014/150] Bugfix: New way of handling constraints in 6-DOF --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 23 +++++++++-------------- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H | 4 ++-- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index fdfb7bc55..fe292bf08 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -76,10 +76,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::A { // Fix the global force for global rotation constraints dimensionedVector fAbs = force(); - vector& fAbsVal = fAbs.value(); - // Constrain rotation - fAbsVal = constrainTranslation(fAbsVal); + // Constrain translation + constrainTranslation(fAbs.value()); return ( @@ -100,9 +99,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot { // Fix the global moment for global rotation constraints dimensionedVector mAbs = moment(); - vector& mAbsVal = mAbs.value(); - mAbsVal = constrainRotation(mAbsVal); + // Constrain rotation + constrainRotation(mAbs.value()); return inv(momentOfInertia_) @@ -124,7 +123,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::E } -Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const +void Foam::sixDOFqODE::constrainRotation(vector& vec) const { vector consVec(vector::zero); @@ -165,12 +164,10 @@ Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const consVec.z() = 0; } } - - return consVec; } -Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const +void Foam::sixDOFqODE::constrainTranslation(vector& vec) const { vector consVec(vector::zero); @@ -211,8 +208,6 @@ Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const consVec.z() = 0; } } - - return consVec; } @@ -263,7 +258,7 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) false ) ), - referentRotation_(vector::zero, 0) + referentRotation_(vector(1, 0, 0), 0) { setCoeffs(); } @@ -503,7 +498,7 @@ void Foam::sixDOFqODE::update(const scalar delta) Uval.z() = coeffs_[5]; // Constrain velocity - Uval = constrainTranslation(Uval); + constrainTranslation(Uval); coeffs_[3] = Uval.x(); coeffs_[4] = Uval.y(); coeffs_[5] = Uval.z(); @@ -516,7 +511,7 @@ void Foam::sixDOFqODE::update(const scalar delta) omegaVal.z() = coeffs_[8]; // Constrain omega - omegaVal = constrainRotation(omegaVal); + constrainRotation(omegaVal); coeffs_[6] = omegaVal.x(); coeffs_[7] = omegaVal.y(); coeffs_[8] = omegaVal.z(); diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H index fbe85f95d..e61f54f87 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H @@ -199,11 +199,11 @@ class sixDOFqODE //- Constrain rotation vector in referent or global coordinate // system - vector constrainRotation(vector& vec) const; + void constrainRotation(vector& vec) const; //- Constrain translation vector in referent or global coordinate // system - vector constrainTranslation(vector& vec) const; + void constrainTranslation(vector& vec) const; public: From 9705581d89fa6d639f55ea3f6ddb47853e311852 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 13:11:55 +0000 Subject: [PATCH 015/150] Bugfix: force calculation of fluxes after GGI motion --- .../dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C | 5 +++-- src/dynamicMesh/dynamicFvMesh/turboFvMesh/turboFvMesh.C | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/dynamicMesh/dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C b/src/dynamicMesh/dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C index c1b0585e6..1297f6aa2 100644 --- a/src/dynamicMesh/dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C +++ b/src/dynamicMesh/dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C @@ -322,8 +322,9 @@ bool Foam::mixerGgiFvMesh::update() ) ); - // The mesh is not morphing - return false; + // The mesh is not morphing, but flux re-calculation is required + // HJ, 25/Jan/2016 + return true; } diff --git a/src/dynamicMesh/dynamicFvMesh/turboFvMesh/turboFvMesh.C b/src/dynamicMesh/dynamicFvMesh/turboFvMesh/turboFvMesh.C index 1300ea2ce..1a8d26335 100644 --- a/src/dynamicMesh/dynamicFvMesh/turboFvMesh/turboFvMesh.C +++ b/src/dynamicMesh/dynamicFvMesh/turboFvMesh/turboFvMesh.C @@ -235,8 +235,9 @@ bool Foam::turboFvMesh::update() ) ); - // The mesh is not morphing - return false; + // The mesh is not morphing, but flux re-calculation is required + // HJ, 25/Jan/2016 + return true; } From 986bce03e5316fa8be5f7622e3dd85b16205b915 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 13:13:50 +0000 Subject: [PATCH 016/150] Bugfix: corrected velocity transformation --- src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C index d1bdc2bf9..f9a139ee1 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C @@ -74,7 +74,7 @@ Foam::solidBodyMotionFunctions::noMotion::transformation() const Foam::septernion Foam::solidBodyMotionFunctions::noMotion::velocity() const { - return septernion(vector::zero, quaternion::zero); + return septernion(vector::zero, quaternion::I)/time_.deltaT().value(); } From 11122b385ef8fd6de57aebe7ef2e7f0c8a040281 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Feb 2016 12:34:45 +0000 Subject: [PATCH 017/150] Bug fix: memory leak in block AMG --- src/foam/meshes/lduMesh/lduPrimitiveMesh.H | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/foam/meshes/lduMesh/lduPrimitiveMesh.H b/src/foam/meshes/lduMesh/lduPrimitiveMesh.H index 098403d18..f001de9d4 100644 --- a/src/foam/meshes/lduMesh/lduPrimitiveMesh.H +++ b/src/foam/meshes/lduMesh/lduPrimitiveMesh.H @@ -25,7 +25,8 @@ Class Foam::lduPrimitiveMesh Description - Simplest contrete lduMesh which stores the addressing needed bu lduMatrix. + Simplest contrete lduMesh which stores the addressing needed by the + lduMatrix. Addressing is contained locally \*---------------------------------------------------------------------------*/ @@ -115,7 +116,18 @@ public: // Destructor virtual ~lduPrimitiveMesh() - {} + { + // Clear interface pointers: they are not deleted + // by the lduInterfacePtrsList (derived from UPtrList) + // HJ, 5/Feb/2016 + forAll (interfaces_, i) + { + if (interfaces_.set(i)) + { + delete interfaces_(i); + } + } + } // Member Functions From 9710a2246e31799fb6c63aaa73be23f9c0cfc320 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Feb 2016 13:47:27 +0000 Subject: [PATCH 018/150] Bugfix: better signalling --- src/lduSolvers/amg/coarseAmgLevel.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lduSolvers/amg/coarseAmgLevel.C b/src/lduSolvers/amg/coarseAmgLevel.C index c4a9deda5..133edf0e2 100644 --- a/src/lduSolvers/amg/coarseAmgLevel.C +++ b/src/lduSolvers/amg/coarseAmgLevel.C @@ -208,7 +208,7 @@ void Foam::coarseAmgLevel::solve return; } - // Switch of debug in top-level direct solve + // Switch off debug in top-level direct solve label oldDebug = lduMatrix::debug(); if (matrixPtr_->matrix().symmetric()) @@ -261,7 +261,7 @@ void Foam::coarseAmgLevel::solve coarseSolverPerf.print(); } - if (lduMatrix::debug >= 2) + if (lduMatrix::debug >= 3) { coarseSolverPerf.print(); } From 82012c585bdb9cff91cd1b1e41af4c58dc3a6d4e Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Feb 2016 13:47:56 +0000 Subject: [PATCH 019/150] Bugfix: Removed coupled solver settings --- .../transonicBump/system/fvSolution | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution index e499e4b91..65a552736 100644 --- a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution +++ b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution @@ -16,19 +16,6 @@ FoamFile solvers { - // Coupled - Up - { - solver BiCGStab; - preconditioner Cholesky; - - tolerance 1e-09; - relTol 0.0; - - minIter 1; - maxIter 500; - } - // Segregated p { From acb0ba58a43543047680595d0bab65d5638e0cf8 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Feb 2016 15:10:24 +0000 Subject: [PATCH 020/150] Bugfix: simplified laplacian schemes --- .../incompressible/pUCoupledFoam/cavity/system/fvSchemes | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes index 11c47fd1e..c96582ac8 100644 --- a/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes +++ b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes @@ -35,9 +35,7 @@ divSchemes laplacianSchemes { - default none; - laplacian(nuEff,U) Gauss linear corrected; - laplacian(rUAf,p) Gauss linear corrected; + default Gauss linear corrected; } interpolationSchemes From f42adf48abc51b660e07ce8e38cf853307fc871d Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 5 Feb 2016 17:48:36 +0000 Subject: [PATCH 021/150] Bugfix: regular formulation of enthalpy source terms --- .../steadyCompressibleMRFFoam/createFields.H | 2 +- .../compressible/steadyCompressibleMRFFoam/hEqn.H | 7 ++++++- .../compressible/steadyCompressibleMRFFoam/iEqn.H | 15 +++------------ .../compressible/steadyCompressibleMRFFoam/pEqn.H | 3 ++- 4 files changed, 12 insertions(+), 15 deletions(-) diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H index 00073568e..78c32d810 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H @@ -93,4 +93,4 @@ ), h ); - i -= 0.5*magSqr(Urot); + i == h - 0.5*(magSqr(Urot) - magSqr(Urel)); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H index 08392027f..6bf1b5f52 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H @@ -2,7 +2,12 @@ // Solve the enthalpy equation T.storePrevIter(); - surfaceScalarField faceU = phi/fvc::interpolate(rho); + // Calculate face velocity from flux + surfaceScalarField faceU + ( + "faceU", + phi/fvc::interpolate(rho) + ); fvScalarMatrix hEqn ( diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H index ca500d7af..ef6c36c13 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H @@ -26,9 +26,10 @@ fvm::ddt(rho, i) + fvm::div(phi, i) - fvm::laplacian(turbulence->alphaEff(), i) - // u & gradP term (steady-state formulation) - + fvm::SuSp((fvc::div(faceU, p, "div(U,p)") - fvc::div(faceU)*p)/i, i) == + // ddt(p) term removed: steady-state. HJ, 27/Apr/2010 + fvc::div(faceU, p, "div(U,p)") + - p*fvc::div(faceU) // Viscous heating: note sign (devRhoReff has a minus in it) - (turbulence->devRhoReff() && fvc::grad(Urel)) ); @@ -42,16 +43,6 @@ h = i + 0.5*magSqr(Urot); h.correctBoundaryConditions(); - // Bound the enthalpy using TMin and TMax - volScalarField Cp = thermo.Cp(); - - h = Foam::min(h, TMax*Cp); - h = Foam::max(h, TMin*Cp); - h.correctBoundaryConditions(); - - // Re-initialise rothalpy based on limited enthalpy - i = h - 0.5*magSqr(Urot); - thermo.correct(); psis = thermo.psi()/thermo.Cp()*thermo.Cv(); } diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index 868f74040..f66193394 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -4,9 +4,10 @@ surfaceScalarField psisf = fvc::interpolate(psis); surfaceScalarField rhof = fvc::interpolate(rho); - // Needs to be outside of loop since p is changing, but psi and rho are not. + // Needs to be outside of loop since p is changing, but psi and rho are not surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p); + // --- PISO loop for (int corr = 0; corr < nCorr; corr++) { U = rUA*UEqn.H(); From f9698bfc2f47e9b8a1a2b390bc67bad7d95ac8f3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 8 Feb 2016 17:48:59 +0000 Subject: [PATCH 022/150] Formatting --- .../cfdTools/general/porousMedia/porousZone.C | 8 +++++--- .../cfdTools/general/porousMedia/porousZones.C | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C index 0bfcf2e09..8c82e3af2 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C @@ -353,8 +353,8 @@ void Foam::porousZone::addResistance if (correctAUprocBC) { // Correct the boundary conditions of the tensorial diagonal to ensure - // processor boundaries are correctly handled when AU^-1 is interpolated - // for the pressure equation. + // processor boundaries are correctly handled when AU^-1 is + // interpolated for the pressure equation. AU.correctBoundaryConditions(); } } @@ -383,7 +383,9 @@ void Foam::porousZone::writeDict(Ostream& os, bool subDict) const if (dict_.found("porosity")) { - os.writeKeyword("porosity") << porosity() << token::END_STATEMENT << nl; + os.writeKeyword("porosity") + << porosity() << token::END_STATEMENT + << nl; } // powerLaw coefficients diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C index 6caf08430..e9063a6e8 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C @@ -75,7 +75,7 @@ void Foam::porousZones::addResistance ) const { // addResistance for each zone, delaying the correction of the - // precessor BCs of AU + // processor BCs of AU forAll(*this, i) { operator[](i).addResistance(UEqn, AU, false); From cd4951da8842713a6c64119256355ca9f89f6e79 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 16 Feb 2016 10:43:34 +0000 Subject: [PATCH 023/150] New tutorial: MRF mixer with GGI interfaces --- .../MRFSimpleFoam/mixerGgiMRF/0/U | 50 ++++ .../MRFSimpleFoam/mixerGgiMRF/0/p | 48 ++++ .../MRFSimpleFoam/mixerGgiMRF/Allclean | 10 + .../MRFSimpleFoam/mixerGgiMRF/Allrun | 13 + .../MRFSimpleFoam/mixerGgiMRF/boundary | 58 ++++ .../mixerGgiMRF/constant/MRFZones | 31 ++ .../mixerGgiMRF/constant/RASProperties | 24 ++ .../mixerGgiMRF/constant/dynamicMeshDict | 40 +++ .../constant/polyMesh/blockMeshDict | 271 ++++++++++++++++++ .../mixerGgiMRF/constant/transportProperties | 21 ++ .../MRFSimpleFoam/mixerGgiMRF/setBatch.batch | 4 + .../mixerGgiMRF/system/controlDict | 61 ++++ .../mixerGgiMRF/system/decomposeParDict | 35 +++ .../mixerGgiMRF/system/fvSchemes | 60 ++++ .../mixerGgiMRF/system/fvSolution | 54 ++++ 15 files changed, 780 insertions(+) create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/U create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/p create mode 100755 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allclean create mode 100755 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allrun create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/boundary create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/MRFZones create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/RASProperties create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/dynamicMeshDict create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/polyMesh/blockMeshDict create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/transportProperties create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/setBatch.batch create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/controlDict create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/decomposeParDict create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSchemes create mode 100644 tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSolution diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/U b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/U new file mode 100644 index 000000000..f2a729e72 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/U @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + impellerWall + { + type movingWallVelocity; + value uniform (0 0 0); + } + baffleWall + { + type fixedValue; + value uniform (0 0 0); + } + insideSlider + { + type ggi; + value uniform (0 0 0); + } + outsideSlider + { + type ggi; + value uniform (0 0 0); + } + defaultFaces + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/p b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/p new file mode 100644 index 000000000..25ead2d46 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/0/p @@ -0,0 +1,48 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + impellerWall + { + type zeroGradient; + } + baffleWall + { + type zeroGradient; + } + insideSlider + { + type ggi; + value uniform 0; + } + outsideSlider + { + type ggi; + value uniform 0; + } + defaultFaces + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allclean b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allclean new file mode 100755 index 000000000..64dded4a3 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allclean @@ -0,0 +1,10 @@ +#!/bin/sh + +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +cleanCase +\rm -rf VTK +\rm -rf 0/polyMesh +\rm -f constant/polyMesh/boundary +\rm -rf constant/polyMesh/sets diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allrun b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allrun new file mode 100755 index 000000000..5344e7fed --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/Allrun @@ -0,0 +1,13 @@ +#!/bin/sh +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +application="MRFSimpleFoam" + +runApplication blockMesh +\cp -f boundary constant/polyMesh/boundary +runApplication regionCellSets +runApplication setSet -batch setBatch.batch +\rm -f constant/polyMesh/sets/cellRegion* +runApplication setsToZones -noFlipMap +runApplication $application diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/boundary b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/boundary new file mode 100644 index 000000000..4f8c9e7aa --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/boundary @@ -0,0 +1,58 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +5 +( + impellerWall + { + type wall; + nFaces 68; + startFace 1040; + } + baffleWall + { + type wall; + nFaces 84; + startFace 1108; + } + insideSlider + { + type ggi; + nFaces 36; + startFace 1192; + shadowPatch outsideSlider; + zone insideZone; + bridgeOverlap false; + } + outsideSlider + { + type ggi; + nFaces 36; + startFace 1228; + shadowPatch insideSlider; + zone outsideZone; + bridgeOverlap false; + } + defaultFaces + { + type empty; + nFaces 1152; + startFace 1264; + } +) + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/MRFZones b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/MRFZones new file mode 100644 index 000000000..6938817b2 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/MRFZones @@ -0,0 +1,31 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object MRFZones; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +1 +( + rotor + { + // Fixed patches (by default they 'move' with the MRF zone) + nonRotatingPatches (); + + origin origin [0 1 0 0 0 0 0] (0 0 0); + axis axis [0 0 0 0 0 0 0] (0 0 1); + omega omega [0 0 -1 0 0 0 0] 6.2831853; + } +) + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/RASProperties b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/RASProperties new file mode 100644 index 000000000..9ab4ecf1d --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/RASProperties @@ -0,0 +1,24 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object RASProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +RASModel laminar; + +turbulence on; + +printCoeffs on; + + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/dynamicMeshDict b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/dynamicMeshDict new file mode 100644 index 000000000..857f0b55a --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/dynamicMeshDict @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object dynamicMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dynamicFvMeshLib "libtopoChangerFvMesh.so"; +dynamicFvMesh mixerGgiFvMesh; + +mixerGgiFvMeshCoeffs +{ + coordinateSystem + { + type cylindrical; + origin (0 0 0); + axis (0 0 1); + direction (1 0 0); + } + + // Rotational speed in revolutions per minute + rpm 60; + + slider + { + moving ( insideSlider ); + static ( outsideSlider ); + } +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/polyMesh/blockMeshDict b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/polyMesh/blockMeshDict new file mode 100644 index 000000000..a82e9d9c9 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/polyMesh/blockMeshDict @@ -0,0 +1,271 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 0.1; + +vertices +( + (-5 0.01 0) + (-5 -0.01 0) + (-4 0 0) + (-3 0 0) + (-3 0 0) + (-2 0 0) + (-1 0.01 0) + (-1 -0.01 0) + (1 0.01 0) + (1 -0.01 0) + (2 0 0) + (3 0 0) + (3 0 0) + (4 0 0) + (5 0.01 0) + (5 -0.01 0) + (0 -3 0) + (0 -2 0) + (-0.01 -1 0) + (0.01 -1 0) + (-0.01 1 0) + (0.01 1 0) + (0 2 0) + (0 3 0) + (-2.50433 -4.32763 0) + (-2.49567 -4.33263 0) + (-2 -3.4641 0) + (-1.5 -2.59808 0) + (2.49567 -4.33263 0) + (2.50433 -4.32763 0) + (2 -3.4641 0) + (1.5 -2.59808 0) + (-2.50433 4.32763 0) + (-2.49567 4.33263 0) + (-2 3.4641 0) + (-1.5 2.59808 0) + (2.49567 4.33263 0) + (2.50433 4.32763 0) + (2 3.4641 0) + (1.5 2.59808 0) + (-5 0.01 0.5) + (-5 -0.01 0.5) + (-4 0 0.5) + (-3 0 0.5) + (-3 0 0.5) + (-2 0 0.5) + (-1 0.01 0.5) + (-1 -0.01 0.5) + (1 0.01 0.5) + (1 -0.01 0.5) + (2 0 0.5) + (3 0 0.5) + (3 0 0.5) + (4 0 0.5) + (5 0.01 0.5) + (5 -0.01 0.5) + (0 -3 0.5) + (0 -2 0.5) + (-0.01 -1 0.5) + (0.01 -1 0.5) + (-0.01 1 0.5) + (0.01 1 0.5) + (0 2 0.5) + (0 3 0.5) + (-2.50433 -4.32763 0.5) + (-2.49567 -4.33263 0.5) + (-2 -3.4641 0.5) + (-1.5 -2.59808 0.5) + (2.49567 -4.33263 0.5) + (2.50433 -4.32763 0.5) + (2 -3.4641 0.5) + (1.5 -2.59808 0.5) + (-2.50433 4.32763 0.5) + (-2.49567 4.33263 0.5) + (-2 3.4641 0.5) + (-1.5 2.59808 0.5) + (2.49567 4.33263 0.5) + (2.50433 4.32763 0.5) + (2 3.4641 0.5) + (1.5 2.59808 0.5) +); + +blocks +( + hex (5 17 18 7 45 57 58 47) (9 4 1) simpleGrading (1 1 1) + hex (17 10 9 19 57 50 49 59) (9 4 1) simpleGrading (1 1 1) + hex (10 22 21 8 50 62 61 48) (9 4 1) simpleGrading (1 1 1) + hex (22 5 6 20 62 45 46 60) (9 4 1) simpleGrading (1 1 1) + hex (4 16 17 5 44 56 57 45) (9 4 1) simpleGrading (1 1 1) + hex (16 11 10 17 56 51 50 57) (9 4 1) simpleGrading (1 1 1) + hex (11 23 22 10 51 63 62 50) (9 4 1) simpleGrading (1 1 1) + hex (23 4 5 22 63 44 45 62) (9 4 1) simpleGrading (1 1 1) + hex (2 26 27 3 42 66 67 43) (6 4 1) simpleGrading (1 1 1) + hex (26 30 31 27 66 70 71 67) (6 4 1) simpleGrading (1 1 1) + hex (30 13 12 31 70 53 52 71) (6 4 1) simpleGrading (1 1 1) + hex (13 38 39 12 53 78 79 52) (6 4 1) simpleGrading (1 1 1) + hex (38 34 35 39 78 74 75 79) (6 4 1) simpleGrading (1 1 1) + hex (34 2 3 35 74 42 43 75) (6 4 1) simpleGrading (1 1 1) + hex (1 24 26 2 41 64 66 42) (6 4 1) simpleGrading (1 1 1) + hex (25 28 30 26 65 68 70 66) (6 4 1) simpleGrading (1 1 1) + hex (29 15 13 30 69 55 53 70) (6 4 1) simpleGrading (1 1 1) + hex (14 37 38 13 54 77 78 53) (6 4 1) simpleGrading (1 1 1) + hex (36 33 34 38 76 73 74 78) (6 4 1) simpleGrading (1 1 1) + hex (32 0 2 34 72 40 42 74) (6 4 1) simpleGrading (1 1 1) +); + +edges +( + arc 7 18 (-0.707107 -0.707107 0) + arc 19 9 (0.707107 -0.707107 0) + arc 8 21 (0.707107 0.707107 0) + arc 20 6 (-0.707107 0.707107 0) + arc 5 17 (-1.41421 -1.41421 0) + arc 17 10 (1.41421 -1.41421 0) + arc 10 22 (1.41421 1.41421 0) + arc 22 5 (-1.41421 1.41421 0) + arc 4 16 (-2.12132 -2.12132 0) + arc 16 11 (2.12132 -2.12132 0) + arc 11 23 (2.12132 2.12132 0) + arc 23 4 (-2.12132 2.12132 0) + arc 3 27 (-2.59808 -1.5 0) + arc 27 31 (0 -3 0) + arc 31 12 (2.59808 -1.5 0) + arc 12 39 (2.59808 1.5 0) + arc 39 35 (0 3 0) + arc 35 3 (-2.59808 1.5 0) + arc 2 26 (-3.4641 -2 0) + arc 26 30 (0 -4 0) + arc 30 13 (3.4641 -2 0) + arc 13 38 (3.4641 2 0) + arc 38 34 (0 4 0) + arc 34 2 (-3.4641 2 0) + arc 1 24 (-4.33013 -2.5 0) + arc 25 28 (0 -5 0) + arc 29 15 (4.33013 -2.5 0) + arc 14 37 (4.33013 2.5 0) + arc 36 33 (0 5 0) + arc 32 0 (-4.33013 2.5 0) + arc 47 58 (-0.707107 -0.707107 0.5) + arc 59 49 (0.707107 -0.707107 0.5) + arc 48 61 (0.707107 0.707107 0.5) + arc 60 46 (-0.707107 0.707107 0.5) + arc 45 57 (-1.41421 -1.41421 0.5) + arc 57 50 (1.41421 -1.41421 0.5) + arc 50 62 (1.41421 1.41421 0.5) + arc 62 45 (-1.41421 1.41421 0.5) + arc 44 56 (-2.12132 -2.12132 0.5) + arc 56 51 (2.12132 -2.12132 0.5) + arc 51 63 (2.12132 2.12132 0.5) + arc 63 44 (-2.12132 2.12132 0.5) + arc 43 67 (-2.59808 -1.5 0.5) + arc 67 71 (0 -3 0.5) + arc 71 52 (2.59808 -1.5 0.5) + arc 52 79 (2.59808 1.5 0.5) + arc 79 75 (0 3 0.5) + arc 75 43 (-2.59808 1.5 0.5) + arc 42 66 (-3.4641 -2 0.5) + arc 66 70 (0 -4 0.5) + arc 70 53 (3.4641 -2 0.5) + arc 53 78 (3.4641 2 0.5) + arc 78 74 (0 4 0.5) + arc 74 42 (-3.4641 2 0.5) + arc 41 64 (-4.33013 -2.5 0.5) + arc 65 68 (0 -5 0.5) + arc 69 55 (4.33013 -2.5 0.5) + arc 54 77 (4.33013 2.5 0.5) + arc 76 73 (0 5 0.5) + arc 72 40 (-4.33013 2.5 0.5) +); + +boundary +( + impellerWall + { + type wall; + faces + ( + + (7 47 58 18) + (18 58 57 17) + (17 57 59 19) + (19 59 49 9) + (9 49 50 10) + (10 50 48 8) + (8 48 61 21) + (21 61 62 22) + (22 62 60 20) + (20 60 46 6) + (6 46 45 5) + (5 45 47 7) + ); + } + + baffleWall + { + type wall; + faces + ( + (1 24 64 41) + (24 26 66 64) + (26 25 65 66) + (25 28 68 65) + (28 30 70 68) + (30 29 69 70) + (29 15 55 69) + (15 13 53 55) + (13 14 54 53) + (14 37 77 54) + (37 38 78 77) + (38 36 76 78) + (36 33 73 76) + (33 34 74 73) + (34 32 72 74) + (32 0 40 72) + (0 2 42 40) + (2 42 41 1) + ); + } + + insideSlider + { + type patch; + faces + ( + (4 16 56 44) + (16 11 51 56) + (11 23 63 51) + (23 4 44 63) + ); + } + + outsideSlider + { + type patch; + faces + ( + (3 43 67 27) + (27 67 71 31) + (31 71 52 12) + (12 52 79 39) + (39 79 75 35) + (35 75 43 3) + ); + } +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/transportProperties b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/transportProperties new file mode 100644 index 000000000..fe8091e78 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/constant/transportProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +transportModel Newtonian; + +nu nu [0 2 -1 0 0 0 0] 0.01; + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/setBatch.batch b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/setBatch.batch new file mode 100644 index 000000000..540348b33 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/setBatch.batch @@ -0,0 +1,4 @@ +faceSet insideZone new patchToFace insideSlider +faceSet outsideZone new patchToFace outsideSlider +cellSet rotor new cellToCell cellRegion0 +quit diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/controlDict b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/controlDict new file mode 100644 index 000000000..9a299d315 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/controlDict @@ -0,0 +1,61 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application MRFSimpleFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 100; + +deltaT 1; + +writeControl timeStep; + +writeInterval 10; + +cycleWrite 0; + +writeFormat ascii; + +writePrecision 10; + +writeCompression compressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +functions +( + ggiCheck + { + // Type of functionObject + type ggiCheck; + + phi phi; + + // Where to load it from (if not already in solver) + functionObjectLibs ("libcheckFunctionObjects.so"); + } +); + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/decomposeParDict b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/decomposeParDict new file mode 100644 index 000000000..820c4383a --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/decomposeParDict @@ -0,0 +1,35 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 4; + +method simple; + +globalFaceZones ( insideZone outsideZone ); + +simpleCoeffs +{ + n (2 2 1); + delta 0.001; +} + +distributed no; + +roots +( +); + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSchemes b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSchemes new file mode 100644 index 000000000..8ac92b9d2 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSchemes @@ -0,0 +1,60 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; + grad(p) Gauss linear; +} + +divSchemes +{ + default none; + div(phi,U) Gauss upwind; + + div((nuEff*dev(grad(U).T()))) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; + interpolate(HbyA) linear; + interpolate(1|A) linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + pcorr; + p; +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSolution b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSolution new file mode 100644 index 000000000..7eab0a022 --- /dev/null +++ b/tutorials/incompressible/MRFSimpleFoam/mixerGgiMRF/system/fvSolution @@ -0,0 +1,54 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.2 | +| \\ / A nd | Web: http://www.foam-extend.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver CG; + preconditioner Cholesky; + + minIter 0; + maxIter 1000; + tolerance 1e-07; + relTol 0.0; + } + U + { + solver BiCGStab; + preconditioner DILU; + + minIter 0; + maxIter 1000; + tolerance 1e-07; + relTol 0; + } +} + +SIMPLE +{ + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + +relaxationFactors +{ + p 0.3; + U 0.7; +} + +// ************************************************************************* // From d1bbe885820d888bd85cecc4763d3c53a376790d Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 17 Feb 2016 16:33:42 +0000 Subject: [PATCH 024/150] Bugfix: constraint must work on a reference. Inno Gatin --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index fe292bf08..059222465 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -145,23 +145,21 @@ void Foam::sixDOFqODE::constrainRotation(vector& vec) const consVec.z() = 0; } - consVec = referentRotation_.invR() & consVec; + vec = referentRotation_.invR() & consVec; } else { - consVec = vec; - if (fixedRoll_) { - consVec.x() = 0; + vec.x() = 0; } if (fixedPitch_) { - consVec.y() = 0; + vec.y() = 0; } if (fixedYaw_) { - consVec.z() = 0; + vec.z() = 0; } } } @@ -189,23 +187,21 @@ void Foam::sixDOFqODE::constrainTranslation(vector& vec) const consVec.z() = 0; } - consVec = referentRotation_.invR() & consVec; + vec = referentRotation_.invR() & consVec; } else { - consVec = vec; - if (fixedSurge_) { - consVec.x() = 0; + vec.x() = 0; } if (fixedSway_) { - consVec.y() = 0; + vec.y() = 0; } if (fixedHeave_) { - consVec.z() = 0; + vec.z() = 0; } } } From 451aa81a39bb3aceb6bbeb8b3fc22dc42fd43c89 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 8 Mar 2016 12:03:39 +0000 Subject: [PATCH 025/150] Bugfix: rothalpy boundary condition handling. Ilaria De Dominicis --- .../fixedEnthalpyFvPatchScalarField.C | 38 +++++++++++++ .../gradientEnthalpyFvPatchScalarField.C | 56 +++++++++++++++++-- .../mixedEnthalpyFvPatchScalarField.C | 50 +++++++++++++++++ 3 files changed, 139 insertions(+), 5 deletions(-) diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C index 603b94265..7d17ef4d9 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C @@ -115,6 +115,44 @@ void Foam::fixedEnthalpyFvPatchScalarField::updateCoeffs() dimensionedInternalField().name() == db().mangleFileName("i") ) { + // Get access to relative and rotational velocity + const word UrelName("Urel"); + const word UrotName("Urot"); + + if + ( + !this->db().objectRegistry::found(UrelName) + || !this->db().objectRegistry::found(UrotName) + ) + { + // Velocities not available, do not update + InfoIn + ( + "void gradientEnthalpyFvPatchScalarField::" + "updateCoeffs(const vectorField& Up)" + ) << "Velocity fields " << UrelName << " or " + << UrotName << " not found. " + << "Performing enthalpy value update for field " + << this->dimensionedInternalField().name() + << " and patch " << patchi + << endl; + + operator==(thermo.h(Tw, patchi)); + } + else + { + const fvPatchVectorField& Urelp = + lookupPatchField(UrelName); + + const fvPatchVectorField& Urotp = + lookupPatchField(UrotName); + + operator== + ( + thermo.h(Tw, patchi) + - 0.5*(magSqr(Urotp) - magSqr(Urelp)) + ); + } } else { diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C index 2639a4fc6..42eb6be44 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C @@ -121,15 +121,61 @@ void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs() dimensionedInternalField().name() == db().mangleFileName("i") ) { + // Get access to relative and rotational velocity + const word UrelName("Urel"); + const word UrotName("Urot"); + + if + ( + !this->db().objectRegistry::found(UrelName) + || !this->db().objectRegistry::found(UrotName) + ) + { + // Velocities not available, do not update + InfoIn + ( + "void gradientEnthalpyFvPatchScalarField::" + "updateCoeffs(const vectorField& Up)" + ) << "Velocity fields " << UrelName << " or " + << UrotName << " not found. " + << "Performing enthalpy value update for field " + << this->dimensionedInternalField().name() + << " and patch " << patchi + << endl; + + gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad() + + patch().deltaCoeffs()* + ( + thermo.h(Tw, patchi) + - thermo.h(Tw, patch().faceCells()) + ); + } + else + { + const fvPatchVectorField& Urelp = + lookupPatchField(UrelName); + + const fvPatchVectorField& Urotp = + lookupPatchField(UrotName); + + gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad() + + patch().deltaCoeffs()* + ( + thermo.h(Tw, patchi) + - thermo.h(Tw, patch().faceCells()) + ) + - mag(Urotp)*mag(Urotp.snGrad()) + + mag(Urelp)*mag(Urelp.snGrad()); + } } else { gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad() - + patch().deltaCoeffs()* - ( - thermo.hs(Tw, patchi) - - thermo.hs(Tw, patch().faceCells()) - ); + + patch().deltaCoeffs()* + ( + thermo.hs(Tw, patchi) + - thermo.hs(Tw, patch().faceCells()) + ); } fixedGradientFvPatchScalarField::updateCoeffs(); diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C index 6670b3d1d..553f7d1bf 100644 --- a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C +++ b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C @@ -130,6 +130,56 @@ void Foam::mixedEnthalpyFvPatchScalarField::updateCoeffs() dimensionedInternalField().name() == db().mangleFileName("i") ) { + // Get access to relative and rotational velocity + const word UrelName("Urel"); + const word UrotName("Urot"); + + if + ( + !this->db().objectRegistry::found(UrelName) + || !this->db().objectRegistry::found(UrotName) + ) + { + // Velocities not available, do not update + InfoIn + ( + "void gradientEnthalpyFvPatchScalarField::" + "updateCoeffs(const vectorField& Up)" + ) << "Velocity fields " << UrelName << " or " + << UrotName << " not found. " + << "Performing enthalpy value update for field " + << this->dimensionedInternalField().name() + << " and patch " << patchi + << endl; + + refValue() = thermo.h(Tw.refValue(), patchi); + refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad() + + patch().deltaCoeffs()* + ( + thermo.h(Tw, patchi) + - thermo.h(Tw, patch().faceCells()) + ); + } + else + { + const fvPatchVectorField& Urelp = + lookupPatchField(UrelName); + + const fvPatchVectorField& Urotp = + lookupPatchField(UrotName); + + refValue() = thermo.h(Tw.refValue(), patchi) + - 0.5*(magSqr(Urotp) - magSqr(Urelp)); + + refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad() + + patch().deltaCoeffs()* + ( + thermo.h(Tw, patchi) + - thermo.h(Tw, patch().faceCells()) + ) + - mag(Urotp)*mag(Urotp.snGrad()) + + mag(Urelp)*mag(Urelp.snGrad()); + } } else { From fb0b13225ac8d410b00ebcb58bfb52e06d979061 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 8 Mar 2016 14:41:59 +0000 Subject: [PATCH 026/150] Bugfix: reporting y+ in parallel --- .../postProcessing/wall/yPlusLES/yPlusLES.C | 4 ++-- .../postProcessing/wall/yPlusRAS/yPlusRAS.C | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C index 5e9136f25..f1447a5ca 100644 --- a/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C +++ b/applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C @@ -122,8 +122,8 @@ int main(int argc, char *argv[]) Info<< "Patch " << patchi << " named " << currPatch.name() - << " y+ : min: " << min(Yp) << " max: " << max(Yp) - << " average: " << average(Yp) << nl << endl; + << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp) + << " average: " << gAverage(Yp) << nl << endl; } } diff --git a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C index c75871837..33e86a392 100644 --- a/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C +++ b/applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C @@ -91,8 +91,8 @@ void calcIncompressibleYPlus Info<< "Patch " << patchi << " named " << nutPw.patch().name() - << " y+ : min: " << min(Yp) << " max: " << max(Yp) - << " average: " << average(Yp) << nl << endl; + << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp) + << " average: " << gAverage(Yp) << nl << endl; } } @@ -171,8 +171,8 @@ void calcCompressibleYPlus Info<< "Patch " << patchi << " named " << mutPw.patch().name() - << " y+ : min: " << min(Yp) << " max: " << max(Yp) - << " average: " << average(Yp) << nl << endl; + << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp) + << " average: " << gAverage(Yp) << nl << endl; } } @@ -243,8 +243,8 @@ void calcTwoPhaseYPlus Info<< "Patch " << patchi << " named " << nutPw.patch().name() - << " y+ : min: " << min(Yp) << " max: " << max(Yp) - << " average: " << average(Yp) << nl << endl; + << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp) + << " average: " << gAverage(Yp) << nl << endl; } } From 4bb1a7676eb07e15b2ba2e753acd69eb1c0e4d3a Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 8 Mar 2016 20:19:57 +0000 Subject: [PATCH 027/150] Bugfix: translation velocity transformation missing 1/dt --- .../meshMotion/solidBodyMotion/translation/translation.C | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C index 037898975..132453536 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/translation/translation.C @@ -144,7 +144,11 @@ Foam::solidBodyMotionFunctions::translation::transformation() const Foam::septernion Foam::solidBodyMotionFunctions::translation::velocity() const { - septernion TV(rampFactor()*velocity_, quaternion::zero); + septernion TV + ( + rampFactor()*velocity_, + quaternion::I/time_.deltaT().value() + ); Info<< "solidBodyMotionFunctions::translation::transformation(): " << "Time = " << time_.value() << " velocity: " << TV << endl; From 71b947c1db86dfbc3a94853ea2244e7d28ecc9bf Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 14 Mar 2016 11:21:32 +0000 Subject: [PATCH 028/150] Clean-up --- .../transformPoints/transformPoints.C | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C index 43ae87179..3937fb5e2 100644 --- a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C +++ b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C @@ -165,6 +165,14 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + if (args.options().empty()) + { + FatalErrorIn(args.executable()) + << "No options supplied, please use one or more of " + "-translate, -rotate, -scale, or -cylToCart options." + << exit(FatalError); + } + word regionName = polyMesh::defaultRegion; fileName meshDir; @@ -191,14 +199,7 @@ int main(int argc, char *argv[]) ) ); - - if (args.options().empty()) - { - FatalErrorIn(args.executable()) - << "No options supplied, please use one or more of " - "-translate, -rotate, -scale, or -cylToCart options." - << exit(FatalError); - } + // Translation options if (args.optionFound("translate")) { @@ -209,6 +210,8 @@ int main(int argc, char *argv[]) points += transVector; } + // Rotation options + if (args.optionFound("rotate")) { Pair n1n2(args.optionLookup("rotate")()); @@ -298,7 +301,7 @@ int main(int argc, char *argv[]) } } - + // Scale options if (args.optionFound("scale")) { From 291d010e9b9393d72b2db44c850d045f1ba09de8 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 14 Mar 2016 11:22:06 +0000 Subject: [PATCH 029/150] Added rotateAlongVector option --- .../surfaceTransformPoints.C | 27 ++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C index 5bc270880..4966958b4 100644 --- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C +++ b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C @@ -30,6 +30,8 @@ Description - pitch (rotation about y) followed by - yaw (rotation about z) + or -rotateAlongVector (vector and angle in degrees) + The yawPitchRoll does yaw followed by pitch followed by roll. \*---------------------------------------------------------------------------*/ @@ -42,11 +44,11 @@ Description #include "transformField.H" #include "Pair.H" #include "quaternion.H" +#include "RodriguesRotation.H" using namespace Foam; using namespace Foam::mathematicalConstant; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: @@ -59,6 +61,7 @@ int main(int argc, char *argv[]) argList::validArgs.append("output surface file"); argList::validOptions.insert("translate", "vector"); argList::validOptions.insert("rotate", "(vector vector)"); + argList::validOptions.insert("rotateAlongVector", "(vector angleInDegree)"); argList::validOptions.insert("scale", "vector"); argList::validOptions.insert("rollPitchYaw", "(roll pitch yaw)"); argList::validOptions.insert("yawPitchRoll", "(yaw pitch roll)"); @@ -72,7 +75,6 @@ int main(int argc, char *argv[]) Info<< "Writing surf to " << outFileName << " ..." << endl; - if (args.options().empty()) { FatalErrorIn(args.executable()) @@ -85,6 +87,8 @@ int main(int argc, char *argv[]) pointField points(surf1.points()); + // Translation options + if (args.optionFound("translate")) { vector transVector(args.optionLookup("translate")()); @@ -94,6 +98,8 @@ int main(int argc, char *argv[]) points += transVector; } + // Rotation options + if (args.optionFound("rotate")) { Pair n1n2(args.optionLookup("rotate")()); @@ -133,7 +139,6 @@ int main(int argc, char *argv[]) << " pitch " << v.y() << nl << " roll " << v.z() << endl; - // Convert to radians v *= pi/180.0; @@ -148,6 +153,22 @@ int main(int argc, char *argv[]) Info<< "Rotating points by quaternion " << R << endl; points = transform(R, points); } + else if (args.optionFound("rotateAlongVector")) + { + vector rotationAxis; + scalar rotationAngle; + + args.optionLookup("rotateAlongVector")() + >> rotationAxis + >> rotationAngle; + + tensor T = RodriguesRotation(rotationAxis, rotationAngle); + + Info << "Rotating points by " << T << endl; + points = transform(T, points); + } + + // Scale options if (args.optionFound("scale")) { From 73655a655289661b497f487ea7130825a26bfac3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 15 Mar 2016 10:29:30 +0000 Subject: [PATCH 030/150] Added decaying turbulence inlet boundary condition by Kornev --- .../incompressible/LES/Make/files | 3 + .../decayingTurbulenceFvPatchVectorField.C | 584 ++++++++++++++++++ .../decayingTurbulenceFvPatchVectorField.H | 209 +++++++ .../decayingTurbulence/decayingVorton.C | 112 ++++ .../decayingTurbulence/decayingVorton.H | 188 ++++++ 5 files changed, 1096 insertions(+) create mode 100644 src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.C create mode 100644 src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.H create mode 100644 src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.C create mode 100644 src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.H diff --git a/src/turbulenceModels/incompressible/LES/Make/files b/src/turbulenceModels/incompressible/LES/Make/files index 186324fde..ff17b596e 100644 --- a/src/turbulenceModels/incompressible/LES/Make/files +++ b/src/turbulenceModels/incompressible/LES/Make/files @@ -31,6 +31,9 @@ kOmegaSSTSAS/kOmegaSSTSAS.C /* Wall functions */ wallFunctions=derivedFvPatchFields/wallFunctions +derivedFvPatchFields/decayingTurbulence/decayingVorton.C +derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.C + nuSgsWallFunctions=$(wallFunctions)/nuSgsWallFunctions $(nuSgsWallFunctions)/nuSgsWallFunction/nuSgsWallFunctionFvPatchScalarField.C diff --git a/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.C b/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.C new file mode 100644 index 000000000..6cde9de35 --- /dev/null +++ b/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingTurbulenceFvPatchVectorField.C @@ -0,0 +1,584 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "decayingTurbulenceFvPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "volFields.H" +#include "scalarList.H" +#include "vectorList.H" +#include "ListListOps.H" +#include "PstreamReduceOps.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::decayingTurbulenceFvPatchVectorField::updateVortons() +{ + vectorField& patchField = *this; + + vectorField tloc = patch().Cf(); + + List l(Pstream::nProcs()); + List cf(Pstream::nProcs()); + List rf(Pstream::nProcs()); + + l[Pstream::myProcNo()].setSize(LField_.size()); + cf[Pstream::myProcNo()].setSize(this->size()); + rf[Pstream::myProcNo()].setSize(refField_.size()); + + label i = 0; + + forAll (tloc, faceI) + { + l[Pstream::myProcNo()][i] = LField_[faceI]; + cf[Pstream::myProcNo()][i] = tloc[faceI]; + rf[Pstream::myProcNo()][i] = refField_[faceI]; + + ++i; + } + + Pstream::gatherList(l); + Pstream::gatherList(cf); + Pstream::gatherList(rf); + + scalarList L + ( + ListListOps::combine(l, accessOp()) + ); + + vectorList CF + ( + ListListOps::combine(cf, accessOp()) + ); + + vectorList RF + ( + ListListOps::combine(rf, accessOp()) + ); + + List > v(Pstream::nProcs()); + + if (Pstream::master()) + { + forAll (CF, I) + { + if (L[I] > SMALL) + { + scalar x = CF[I].x(); + scalar ls = lSpot(L[I]); + scalar ymin = CF[I].y() - L[I]; + scalar zmin = CF[I].z() - L[I]; + + for (label k = 0; k < nVortexes_; k++) + { + vector v + ( + (direction_ > 0) ? x - ls : x + ls, + ymin + rndGen_.scalar01()*2*L[I], + zmin + rndGen_.scalar01()*2*L[I] + ); + + bool allowed = true; + + for + ( + SLList::const_iterator it = + vortons_.begin(); + it != vortons_.end(); + ++it + ) + { + if (mag(v - it().location()) < vortexOverlap_*ls) + { + allowed = false; + break; + } + } + + if (!allowed) + { + continue; + } + else + { + vortons_.insert + ( + decayingVorton + ( + L[I], + v, + RF[I], + (direction_ > 0) ? x + ls : x - ls + ) + ); + } + } + } + } + + v[Pstream::myProcNo()].setSize(vortons_.size()); + + i = 0; + + for + ( + SLList::iterator it = vortons_.begin(); + it != vortons_.end(); + ++it + ) + { + v[Pstream::myProcNo()][i++] = it(); + } + } + + Pstream::scatterList(v); + + List V + ( + ListListOps::combine > + ( + v, + accessOp >() + ) + ); + + if (!patchField.empty()) + { + vectorField turbulent(refField_.size(), pTraits::zero); + + forAll (patchField, I) + { + vector pt = tloc[I]; + + forAll (V, vI) + { + const decayingVorton& vorton = V[vI]; + + if (mag(vorton.location() - pt) < lSpot(vorton.length())) + { + turbulent[I] += vorton.velocityAt(pt); + } + } + } + + R_ = ((index_ - 1.0)/index_)*R_ + (1.0/index_)*sqr(turbulent); + + tensorField C_(R_.size(), sphericalTensor::I); + + forAll (R_, I) + { + scalar D1 = R_[I].xx(); + + if (D1 > 0) + { + C_[I].xx() = 1.0/sqrt(D1); + } + + scalar D2 = R_[I].xx()*R_[I].yy() - R_[I].xy()*R_[I].xy(); + + if (D1 > 0 && D2 > 0) + { + C_[I].yx() = -R_[I].xy()/sqrt(D1*D2); + C_[I].yy() = sqrt(D1/D2); + } + + scalar D3 = det(R_[I]); + + if (D2 > 0 && D3 > 0) + { + C_[I].zx() = (R_[I].xy()*R_[I].yz() - R_[I].yy()*R_[I].xz())/ + sqrt(D2*D3); + + C_[I].zy() = -(R_[I].xx()*R_[I].yz()-R_[I].xz()*R_[I].xy())/ + sqrt(D2*D3); + + C_[I].zz() = sqrt(D2/D3); + } + } + + index_++; + + turbulent = C_ & turbulent; + turbulent = Lund_ & turbulent; + + fixedValueFvPatchVectorField::operator==(refField_+ turbulent); + } + + if (Pstream::master()) + { + for + ( + SLList::iterator it = vortons_.begin(); + it != vortons_.end(); + ++it + ) + { + it().move(this->db().time().deltaT().value()); + } + + bool modified; + + do + { + modified = false; + + for + ( + SLList::iterator it = vortons_.begin(); + it!=vortons_.end(); + ++it + ) + { + if (direction_ > 0) + { + if (it().location().x() > it().xmax()) + { + vortons_.remove(it); + modified = true; + break; + } + } + else + { + if (it().location().x() < it().xmax()) + { + vortons_.remove(it); + modified = true; + break; + } + } + } + } while (modified); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::decayingTurbulenceFvPatchVectorField:: +decayingTurbulenceFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchVectorField(p, iF), + LField_(p.size()), + refField_(p.size()), + RField_(p.size()), + Lund_(p.size(), pTraits::zero),//HJ??? + R_(p.size()), + direction_(1), + vortexOverlap_(0.3), + nVortexes_(10), + curTimeIndex_(-1), + rndGen_(label(0)), + vortons_(), + index_(2) +{} + + +Foam::decayingTurbulenceFvPatchVectorField:: +decayingTurbulenceFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchVectorField(p, iF), + LField_("LField", dict, p.size()), + refField_("refField", dict, p.size()), + RField_("RField", dict, p.size()), + Lund_(p.size(), pTraits::zero), + R_(RField_), + direction_(readScalar(dict.lookup("direction"))), + vortexOverlap_(dict.lookupOrDefault("vortexOverlap", 0.3)), + nVortexes_(dict.lookupOrDefault