From a6ef8b902736dd99fcc845795c73f369706e28d5 Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Tue, 20 Jun 2017 14:36:15 +0100 Subject: [PATCH 1/5] INT: Integration of isoAdvector and supporting material Community contribution from Johan Roenby, DHI IsoAdvector is a geometric Volume-of-Fluid method for advection of a sharp interface between two incompressible fluids. It works on both structured and unstructured meshes with no requirements on cell shapes. IsoAdvector is as an alternative choice for the interface compression treatment with the MULES limiter implemented in the interFoam family of solvers. The isoAdvector concept and code was developed at DHI and was funded by a Sapere Aude postdoc grant to Johan Roenby from The Danish Council for Independent Research | Technology and Production Sciences (Grant-ID: DFF - 1337-00118B - FTP). Co-funding is also provided by the GTS grant to DHI from the Danish Agency for Science, Technology and Innovation. The ideas behind and performance of the isoAdvector scheme is documented in: Roenby J, Bredmose H, Jasak H. 2016 A computational method for sharp interface advection. R. Soc. open sci. 3: 160405. [http://dx.doi.org/10.1098/rsos.160405](http://dx.doi.org/10.1098/rsos.160405) Videos showing isoAdvector's performance with a number of standard test cases can be found in this youtube channel: https://www.youtube.com/channel/UCt6Idpv4C8TTgz1iUX0prAA Project contributors: * Johan Roenby (Inventor and main developer) * Hrvoje Jasak (Consistent treatment of boundary faces including processor boundaries, parallelisation, code clean up * Henrik Bredmose (Assisted in the conceptual development) * Vuko Vukcevic (Code review, profiling, porting to foam-extend, bug fixing, testing) * Tomislav Maric (Source file rearrangement) * Andy Heather (Integration into OpenFOAM for v1706 release) See the integration repository below to see the full set of changes implemented for release into OpenFOAM v1706 https://develop.openfoam.com/Community/Integration-isoAdvector --- .../multiphase/interIsoFoam/Make/files | 3 + .../multiphase/interIsoFoam/Make/options | 20 + .../solvers/multiphase/interIsoFoam/UEqn.H | 33 + .../multiphase/interIsoFoam/alphaControls.H | 4 + .../multiphase/interIsoFoam/alphaCourantNo.H | 57 + .../multiphase/interIsoFoam/alphaEqn.H | 35 + .../interIsoFoam/alphaEqnSubCycle.H | 35 + .../multiphase/interIsoFoam/correctPhi.H | 11 + .../multiphase/interIsoFoam/createFields.H | 123 ++ .../interIsoFoam/createIsoAdvection.H | 10 + .../multiphase/interIsoFoam/interIsoFoam.C | 143 ++ .../solvers/multiphase/interIsoFoam/pEqn.H | 67 + .../multiphase/interIsoFoam/setDeltaT.H | 53 + .../preProcessing/setAlphaField/Make/files | 3 + .../preProcessing/setAlphaField/Make/options | 9 + .../setAlphaField/setAlphaField.C | 191 +++ .../setAlphaField/setAlphaFieldDict | 24 + src/finiteVolume/Make/files | 3 + .../isoAdvection/isoAdvection/isoAdvection.C | 1165 +++++++++++++++++ .../isoAdvection/isoAdvection/isoAdvection.H | 388 ++++++ .../isoAdvection/isoAdvectionTemplates.C | 111 ++ .../isoAdvection/isoCutCell/isoCutCell.C | 709 ++++++++++ .../isoAdvection/isoCutCell/isoCutCell.H | 200 +++ .../isoAdvection/isoCutFace/isoCutFace.C | 629 +++++++++ .../isoAdvection/isoCutFace/isoCutFace.H | 200 +++ src/functionObjects/field/Make/files | 2 + src/functionObjects/field/setFlow/setFlow.C | 448 +++++++ src/functionObjects/field/setFlow/setFlow.H | 199 +++ .../multiphase/interIsoFoam/damBreak/0.orig/U | 51 + .../interIsoFoam/damBreak/0.orig/alpha.water | 51 + .../interIsoFoam/damBreak/0.orig/p_rgh | 59 + .../multiphase/interIsoFoam/damBreak/Allclean | 9 + .../multiphase/interIsoFoam/damBreak/Allrun | 12 + .../interIsoFoam/damBreak/Allrun-parallel | 13 + .../damBreak/constant/dynamicMeshDict | 21 + .../interIsoFoam/damBreak/constant/g | 22 + .../damBreak/constant/transportProperties | 38 + .../damBreak/constant/turbulenceProperties | 21 + .../damBreak/system/blockMeshDict | 108 ++ .../interIsoFoam/damBreak/system/controlDict | 56 + .../damBreak/system/decomposeParDict | 45 + .../interIsoFoam/damBreak/system/fvSchemes | 60 + .../interIsoFoam/damBreak/system/fvSolution | 76 ++ .../damBreak/system/setFieldsDict | 36 + .../interIsoFoam/discInConstantFlow/0.orig/U | 39 + .../discInConstantFlow/0.orig/alpha.water | 37 + .../discInConstantFlow/0.orig/p_rgh | 37 + .../interIsoFoam/discInConstantFlow/Allclean | 8 + .../interIsoFoam/discInConstantFlow/Allrun | 11 + .../discInConstantFlow/constant/g | 22 + .../constant/transportProperties | 37 + .../constant/turbulenceProperties | 21 + .../discInConstantFlow/system/blockMeshDict | 85 ++ .../discInConstantFlow/system/controlDict | 60 + .../discInConstantFlow/system/fvSchemes | 60 + .../discInConstantFlow/system/fvSolution | 44 + .../system/setAlphaFieldDict | 24 + .../discInReversedVortexFlow/0.orig/U | 39 + .../0.orig/alpha.water | 38 + .../discInReversedVortexFlow/0.orig/p_rgh | 37 + .../discInReversedVortexFlow/Allclean | 7 + .../discInReversedVortexFlow/Allrun | 17 + .../discInReversedVortexFlow/constant/g | 22 + .../constant/transportProperties | 86 ++ .../constant/turbulenceProperties | 21 + .../system/blockMeshDict | 85 ++ .../system/controlDict | 81 ++ .../system/decomposeParDict | 61 + .../discInReversedVortexFlow/system/fvSchemes | 51 + .../system/fvSolution | 43 + .../system/refineMeshDict | 62 + .../system/setAlphaFieldDict | 24 + .../system/topoSetDict | 37 + .../system/topoSetDict2 | 35 + .../notchedDiscInSolidBodyRotation/0.orig/U | 39 + .../0.orig/alpha.water | 38 + .../0.orig/p_rgh | 37 + .../notchedDiscInSolidBodyRotation/Allclean | 9 + .../notchedDiscInSolidBodyRotation/Allrun | 19 + .../notchedDiscInSolidBodyRotation/constant/g | 22 + .../constant/transportProperties | 37 + .../constant/turbulenceProperties | 21 + .../system/blockMeshDict | 85 ++ .../system/controlDict | 76 ++ .../system/fvSchemes | 52 + .../system/fvSolution | 43 + .../system/setAlphaFieldDict | 24 + .../system/setFieldsDict | 31 + .../sphereInReversedVortexFlow/0.orig/U | 31 + .../0.orig/alpha.water | 30 + .../sphereInReversedVortexFlow/0.orig/p_rgh | 30 + .../sphereInReversedVortexFlow/Allclean | 7 + .../sphereInReversedVortexFlow/Allrun | 13 + .../sphereInReversedVortexFlow/constant/g | 22 + .../constant/transportProperties | 37 + .../constant/turbulenceProperties | 21 + .../system/blockMeshDict | 71 + .../system/controlDict | 81 ++ .../system/decomposeParDict | 31 + .../system/fvSchemes | 51 + .../system/fvSolution | 43 + .../system/setAlphaFieldDict | 23 + .../interIsoFoam/standingWave/0.orig/U | 51 + .../standingWave/0.orig/alpha.water | 49 + .../interIsoFoam/standingWave/0.orig/p_rgh | 53 + .../interIsoFoam/standingWave/Allclean | 9 + .../interIsoFoam/standingWave/Allrun | 17 + .../standingWave/constant/dynamicMeshDict | 21 + .../interIsoFoam/standingWave/constant/g | 22 + .../standingWave/constant/transportProperties | 37 + .../constant/turbulenceProperties | 21 + .../standingWave/system/blockMeshDict | 106 ++ .../standingWave/system/controlDict | 56 + .../standingWave/system/decomposeParDict | 45 + .../standingWave/system/fvSchemes | 59 + .../standingWave/system/fvSolution | 75 ++ .../standingWave/system/refineMeshDict1 | 62 + .../standingWave/system/refineMeshDict2 | 62 + .../standingWave/system/setAlphaFieldDict | 26 + .../standingWave/system/topoSetDict1 | 32 + .../standingWave/system/topoSetDict2 | 32 + 121 files changed, 8643 insertions(+) create mode 100644 applications/solvers/multiphase/interIsoFoam/Make/files create mode 100644 applications/solvers/multiphase/interIsoFoam/Make/options create mode 100644 applications/solvers/multiphase/interIsoFoam/UEqn.H create mode 100644 applications/solvers/multiphase/interIsoFoam/alphaControls.H create mode 100644 applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H create mode 100644 applications/solvers/multiphase/interIsoFoam/alphaEqn.H create mode 100644 applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H create mode 100644 applications/solvers/multiphase/interIsoFoam/correctPhi.H create mode 100644 applications/solvers/multiphase/interIsoFoam/createFields.H create mode 100644 applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H create mode 100644 applications/solvers/multiphase/interIsoFoam/interIsoFoam.C create mode 100644 applications/solvers/multiphase/interIsoFoam/pEqn.H create mode 100644 applications/solvers/multiphase/interIsoFoam/setDeltaT.H create mode 100755 applications/utilities/preProcessing/setAlphaField/Make/files create mode 100755 applications/utilities/preProcessing/setAlphaField/Make/options create mode 100644 applications/utilities/preProcessing/setAlphaField/setAlphaField.C create mode 100644 applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvectionTemplates.C create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.H create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C create mode 100644 src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H create mode 100644 src/functionObjects/field/setFlow/setFlow.C create mode 100644 src/functionObjects/field/setFlow/setFlow.H create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/damBreak/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/damBreak/Allrun create mode 100755 tutorials/multiphase/interIsoFoam/damBreak/Allrun-parallel create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/constant/dynamicMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/decomposeParDict create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/damBreak/system/setFieldsDict create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/discInConstantFlow/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/discInConstantFlow/Allrun create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/discInConstantFlow/system/setAlphaFieldDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allrun create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/decomposeParDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/refineMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/setAlphaFieldDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict create mode 100644 tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict2 create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allrun create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setAlphaFieldDict create mode 100644 tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setFieldsDict create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allrun create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/decomposeParDict create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/setAlphaFieldDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/0.orig/U create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/0.orig/alpha.water create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/0.orig/p_rgh create mode 100755 tutorials/multiphase/interIsoFoam/standingWave/Allclean create mode 100755 tutorials/multiphase/interIsoFoam/standingWave/Allrun create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/constant/dynamicMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/constant/g create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/constant/transportProperties create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/constant/turbulenceProperties create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/blockMeshDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/controlDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/decomposeParDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/fvSchemes create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/fvSolution create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict1 create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict2 create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/setAlphaFieldDict create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict1 create mode 100644 tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict2 diff --git a/applications/solvers/multiphase/interIsoFoam/Make/files b/applications/solvers/multiphase/interIsoFoam/Make/files new file mode 100644 index 0000000000..10f61a6885 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/Make/files @@ -0,0 +1,3 @@ +interIsoFoam.C + +EXE = $(FOAM_APPBIN)/interIsoFoam diff --git a/applications/solvers/multiphase/interIsoFoam/Make/options b/applications/solvers/multiphase/interIsoFoam/Make/options new file mode 100644 index 0000000000..cfe4a92072 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/Make/options @@ -0,0 +1,20 @@ +EXE_INC = \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -limmiscibleIncompressibleTwoPhaseMixture \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling diff --git a/applications/solvers/multiphase/interIsoFoam/UEqn.H b/applications/solvers/multiphase/interIsoFoam/UEqn.H new file mode 100644 index 0000000000..77d1dcd83e --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/UEqn.H @@ -0,0 +1,33 @@ + MRF.correctBoundaryVelocity(U); + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(rho, U) + == + fvOptions(rho, U) + ); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + fvOptions.correct(U); + } diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H new file mode 100644 index 0000000000..92a8fdf951 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H @@ -0,0 +1,4 @@ +const dictionary& alphaControls = mesh.solverDict(alpha1.name()); + +label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); + diff --git a/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H new file mode 100644 index 0000000000..a084155c8a --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + alphaCourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar maxAlphaCo +( + readScalar(runTime.controlDict().lookup("maxAlphaCo")) +); + +scalar alphaCoNum = 0.0; +scalar meanAlphaCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + mixture.nearInterface()().primitiveField() + *fvc::surfaceSum(mag(phi))().primitiveField() + ); + + alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanAlphaCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Interface Courant Number mean: " << meanAlphaCoNum + << " max: " << alphaCoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H new file mode 100644 index 0000000000..cefabfda73 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H @@ -0,0 +1,35 @@ +// If there are more than one outer corrector, we use a mixture of old and +// new U and phi for propagating alpha1 in all but the first outer iteration +if (!pimple.firstIter()) +{ + // We are recalculating alpha1 from the its old time value + alpha1 = alpha1.oldTime(); + // Temporarily storing new U and phi values in prevIter storage + U.storePrevIter(); + phi.storePrevIter(); + + // Overwriting new U and phi values with mixture of old and new values + phi = 0.5*(phi + phi.oldTime()); + U = 0.5*(U + U.oldTime()); +} + +// Update alpha1 +advector.advect(); + +// Update rhoPhi +rhoPhi = advector.getRhoPhi(rho1, rho2); + +alpha2 = 1.0 - alpha1; + +if (!pimple.firstIter()) +{ + // Restoring new U and phi values temporarily saved in prevIter() above + U = U.prevIter(); + phi = phi.prevIter(); +} + +Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(" << alpha1.name() << ") = " << min(alpha1).value() + << " Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1 + << endl; diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H new file mode 100644 index 0000000000..8f0af80e0d --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H @@ -0,0 +1,35 @@ +if (nAlphaSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum + ( + IOobject + ( + "rhoPhiSum", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("0", rhoPhi.dimensions(), 0) + ); + + tmp trSubDeltaT; + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ + #include "alphaEqn.H" +} + +rho == alpha1*rho1 + alpha2*rho2; diff --git a/applications/solvers/multiphase/interIsoFoam/correctPhi.H b/applications/solvers/multiphase/interIsoFoam/correctPhi.H new file mode 100644 index 0000000000..9afcd58a66 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/correctPhi.H @@ -0,0 +1,11 @@ +CorrectPhi +( + U, + phi, + p_rgh, + dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1), + geometricZeroField(), + pimple +); + +#include "continuityErrs.H" diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H new file mode 100644 index 0000000000..8c3a05e3fa --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/createFields.H @@ -0,0 +1,123 @@ +Info<< "Reading field p_rgh\n" << endl; +volScalarField p_rgh +( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "createPhi.H" + + +Info<< "Reading transportProperties\n" << endl; +immiscibleIncompressibleTwoPhaseMixture mixture(U, phi); + +volScalarField& alpha1(mixture.alpha1()); +volScalarField& alpha2(mixture.alpha2()); + +const dimensionedScalar& rho1 = mixture.rho1(); +const dimensionedScalar& rho2 = mixture.rho2(); + + +// Need to store rho for ddt(rho, U) +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + alpha1*rho1 + alpha2*rho2 +); +rho.oldTime(); + + +// Mass flux +surfaceScalarField rhoPhi +( + IOobject + ( + "rhoPhi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::interpolate(rho)*phi +); + + +// Construct incompressible turbulence model +autoPtr turbulence +( + incompressible::turbulenceModel::New(U, phi, mixture) +); + + +#include "readGravitationalAcceleration.H" +#include "readhRef.H" +#include "gh.H" + + +volScalarField p +( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rho*gh +); + +label pRefCell = 0; +scalar pRefValue = 0.0; +setRefCell +( + p, + p_rgh, + pimple.dict(), + pRefCell, + pRefValue +); + +if (p_rgh.needReference()) +{ + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; +} + +mesh.setFluxRequired(p_rgh.name()); +mesh.setFluxRequired(alpha1.name()); + +#include "createMRF.H" +#include "createIsoAdvection.H" diff --git a/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H new file mode 100644 index 0000000000..889af7228f --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H @@ -0,0 +1,10 @@ +// Defining isoAdvection +isoAdvection advector(alpha1, phi, U); + +bool frozenFlow = pimple.dict().lookupOrDefault("frozenFlow", false); +if (frozenFlow) +{ + Info<< "Employing frozen-flow assumption: " + << "pressure-velocity system will not be updated" + << endl; +} diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C new file mode 100644 index 0000000000..1514c5c0cb --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C @@ -0,0 +1,143 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016 DHI +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + interIsoFoam + +Group + grpMultiphaseSolvers + +Description + Solver derived from interFoam for 2 incompressible, isothermal immiscible + fluids using the iso-advector phase-fraction based interface capturing + approach. + + The momentum and other fluid properties are of the "mixture" and a single + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + For a two-fluid approach see twoPhaseEulerFoam. + + Reference: + \verbatim + Roenby, J., Bredmose, H. and Jasak, H. (2016). + A computational method for sharp interface advection + Royal Society Open Science, 3 + doi 10.1098/rsos.160405 + \endverbatim + + isoAdvector code supplied by Johan Roenby, DHI (2016) + +\*---------------------------------------------------------------------------*/ + +#include "isoAdvection.H" +#include "fvCFD.H" +#include "subCycle.H" +#include "immiscibleIncompressibleTwoPhaseMixture.H" +#include "turbulentTransportModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "CorrectPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "correctPhi.H" + + turbulence->validate(); + + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "alphaControls.H" + #include "alphaEqnSubCycle.H" + + mixture.correct(); + + if (frozenFlow) + { + continue; + } + + #include "UEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interIsoFoam/pEqn.H b/applications/solvers/multiphase/interIsoFoam/pEqn.H new file mode 100644 index 0000000000..17c3ae9748 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/pEqn.H @@ -0,0 +1,67 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::flux(HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + MRF.makeRelative(phiHbyA); + adjustPhi(phiHbyA, U, p_rgh); + + surfaceScalarField phig + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() +// - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf + ); + + phiHbyA += phig; + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqn + ( + fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) + ); + + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA - p_rghEqn.flux(); + + p_rgh.relax(); + + U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); + fvOptions.correct(U); + } + } + + #include "continuityErrs.H" + + p == p_rgh + rho*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rho*gh; + } +} diff --git a/applications/solvers/multiphase/interIsoFoam/setDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H new file mode 100644 index 0000000000..9cc860c032 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + setDeltaT + +Description + Reset the timestep to maintain a constant maximum courant Number. + Reduction of time-step is immediate, but increase is damped to avoid + unstable oscillations. + +\*---------------------------------------------------------------------------*/ + +if (adjustTimeStep) +{ + scalar maxDeltaTFact = + min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL)); + + scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); + + runTime.setDeltaT + ( + min + ( + deltaTFact*runTime.deltaTValue(), + maxDeltaT + ) + ); + + Info<< "deltaT = " << runTime.deltaTValue() << endl; +} + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/setAlphaField/Make/files b/applications/utilities/preProcessing/setAlphaField/Make/files new file mode 100755 index 0000000000..23145059e2 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/Make/files @@ -0,0 +1,3 @@ +setAlphaField.C + +EXE = $(FOAM_APPBIN)/setAlphaField diff --git a/applications/utilities/preProcessing/setAlphaField/Make/options b/applications/utilities/preProcessing/setAlphaField/Make/options new file mode 100755 index 0000000000..b4ad3d8bdf --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C new file mode 100644 index 0000000000..6fe8001e31 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C @@ -0,0 +1,191 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016-2017 DHI +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + setAlphaField + +Description + Uses isoCutCell to create a volume fraction field from either a cylinder, + a sphere or a plane. + + Original code supplied by Johan Roenby, DHI (2016) + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "isoCutFace.H" +#include "isoCutCell.H" +#include "Enum.H" +#include "mathematicalConstants.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class shapeSelector +{ + public: + + enum class shapeType + { + PLANE, + SPHERE, + CYLINDER, + SIN + }; + + static const Foam::Enum shapeTypeNames; +}; + +const Foam::Enum shapeSelector::shapeTypeNames +{ + { shapeSelector::shapeType::PLANE, "plane" }, + { shapeSelector::shapeType::SPHERE, "sphere" }, + { shapeSelector::shapeType::CYLINDER, "cylinder" }, + { shapeSelector::shapeType::SIN, "sin" }, +}; + + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + Info<< "Reading setAlphaFieldDict\n" << endl; + + IOdictionary dict + ( + IOobject + ( + "setAlphaFieldDict", + runTime.system(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + const shapeSelector::shapeType surfType + ( + shapeSelector::shapeTypeNames.read(dict.lookup("type")) + ); + const vector centre(dict.lookup("centre")); + const word fieldName(dict.lookup("field")); + + Info<< "Reading field " << fieldName << "\n" << endl; + volScalarField alpha1 + ( + IOobject + ( + fieldName, + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + scalar f0 = 0.0; + scalarField f(mesh.points().size()); + + Info<< "Processing type '" << shapeSelector::shapeTypeNames[surfType] + << "'" << endl; + + switch (surfType) + { + case shapeSelector::shapeType::PLANE: + { + const vector direction(dict.lookup("direction")); + + f = -(mesh.points() - centre) & (direction/mag(direction)); + f0 = 0.0; + break; + } + case shapeSelector::shapeType::SPHERE: + { + const scalar radius(readScalar(dict.lookup("radius"))); + + f = -mag(mesh.points() - centre); + f0 = -radius; + break; + } + case shapeSelector::shapeType::CYLINDER: + { + const scalar radius(readScalar(dict.lookup("radius"))); + const vector direction(dict.lookup("direction")); + + f = -sqrt + ( + sqr(mag(mesh.points() - centre)) + - sqr(mag((mesh.points() - centre) & direction)) + ); + f0 = -radius; + break; + } + case shapeSelector::shapeType::SIN: + { + const scalar period(readScalar(dict.lookup("period"))); + const scalar amplitude(readScalar(dict.lookup("amplitude"))); + const vector up(dict.lookup("up")); + const vector direction(dict.lookup("direction")); + + const scalarField xx + ( + (mesh.points() - centre) & direction/mag(direction) + ); + const scalarField zz((mesh.points() - centre) & up/mag(up)); + + f = amplitude*Foam::sin(2*mathematical::pi*xx/period) - zz; + f0 = 0; + break; + } + } + + + // Define function on mesh points and isovalue + + // Calculating alpha1 volScalarField from f = f0 isosurface + isoCutCell icc(mesh, f); + icc.volumeOfFluid(alpha1, f0); + + // Writing volScalarField alpha1 + ISstream::defaultPrecision(18); + alpha1.write(); + + Info<< nl << "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(" << alpha1.name() << ") = " << min(alpha1).value() + << " Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1 + << nl << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict b/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict new file mode 100644 index 0000000000..6bbb616b54 --- /dev/null +++ b/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict @@ -0,0 +1,24 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +field alpha.water; +type cylinder; +radius 0.25; +direction (0 1 0); +centre (0.5 0 0.5); + +// ************************************************************************* // diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 1087c3e378..ca9bf31ff6 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -238,6 +238,9 @@ fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C fvMatrices/solvers/MULES/MULES.C fvMatrices/solvers/MULES/CMULES.C +fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C +fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C +fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C interpolation = interpolation/interpolation diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C new file mode 100644 index 0000000000..33f63c2b6c --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C @@ -0,0 +1,1165 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- + isoAdvector | Copyright (C) 2016-2017 DHI +------------------------------------------------------------------------------- + +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "isoAdvection.H" +#include "volFields.H" +#include "interpolationCellPoint.H" +#include "volPointInterpolation.H" +#include "fvcSurfaceIntegrate.H" +#include "fvcGrad.H" +#include "upwind.H" +#include "cellSet.H" +#include "meshTools.H" +#include "OBJstream.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(isoAdvection, 0); +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::isoAdvection::isoAdvection +( + volScalarField& alpha1, + const surfaceScalarField& phi, + const volVectorField& U +) +: + // General data + mesh_(alpha1.mesh()), + dict_(mesh_.solverDict(alpha1.name())), + alpha1_(alpha1), + alpha1In_(alpha1.ref()), + phi_(phi), + U_(U), + dVf_ + ( + IOobject + ( + "dVf_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimVol, 0) + ), + advectionTime_(0), + + // Interpolation data + ap_(mesh_.nPoints()), + + // Tolerances and solution controls + nAlphaBounds_(dict_.lookupOrDefault