diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean b/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean new file mode 100755 index 0000000000..99c52e59f9 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean libso phaseChangeTwoPhaseMixtures +wclean +wclean interPhaseChangeDyMFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake b/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake new file mode 100755 index 0000000000..414ed904a0 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake libso phaseChangeTwoPhaseMixtures +wmake +wmake interPhaseChangeDyMFoam + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/files b/applications/solvers/multiphase/interPhaseChangeFoam/Make/files index 5f65a3a4c2..4358afef8b 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/files +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Make/files @@ -1,8 +1,3 @@ interPhaseChangeFoam.C -phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C -phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C -phaseChangeTwoPhaseMixtures/Kunz/Kunz.C -phaseChangeTwoPhaseMixtures/Merkle/Merkle.C -phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C EXE = $(FOAM_APPBIN)/interPhaseChangeFoam diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options index d4ead14721..4721c88961 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options @@ -4,13 +4,14 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ - -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \ + -IphaseChangeTwoPhaseMixtures/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/fvOptions/lnInclude\ -I$(LIB_SRC)/sampling/lnInclude EXE_LIBS = \ + -lphaseChangeTwoPhaseMixtures \ -ltwoPhaseMixture \ -linterfaceProperties \ -ltwoPhaseProperties \ diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/files b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/files new file mode 100644 index 0000000000..2d4521425f --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/files @@ -0,0 +1,3 @@ +interPhaseChangeDyMFoam.C + +EXE = $(FOAM_APPBIN)/interPhaseChangeDyMFoam diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options new file mode 100644 index 0000000000..3bf8a5dfb6 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options @@ -0,0 +1,31 @@ +EXE_INC = \ + -I.. \ + -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/incompressible/turbulenceModel \ + -I../phaseChangeTwoPhaseMixtures/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude\ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lphaseChangeTwoPhaseMixtures \ + -ltwoPhaseMixture \ + -linterfaceProperties \ + -ltwoPhaseProperties \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume \ + -ldynamicMesh \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -lmeshTools \ + -lfvOptions \ + -lsampling diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C new file mode 100644 index 0000000000..da11697426 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2013 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 . + +Application + interPhaseChangeDyMFoam + +Description + Solver for 2 incompressible, isothermal immiscible fluids with phase-change + (e.g. cavitation). Uses a VOF (volume of fluid) phase-fraction based + interface capturing approach, with optional mesh motion and mesh topology + changes including adaptive re-meshing. + + The momentum and other fluid properties are of the "mixture" and a + single momentum equation is solved. + + The set of phase-change models provided are designed to simulate cavitation + but other mechanisms of phase-change are supported within this solver + framework. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "MULES.H" +#include "subCycle.H" +#include "interfaceProperties.H" +#include "phaseChangeTwoPhaseMixture.H" +#include "turbulenceModel.H" +#include "pimpleControl.H" +#include "fvIOoptionList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "readGravitationalAcceleration.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "readTimeControls.H" + + pimpleControl pimple(mesh); + + surfaceScalarField phiAbs("phiAbs", phi); + fvc::makeAbsolute(phiAbs, U); + + #include "../interFoam/interDyMFoam/correctPhi.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "../interFoam/interDyMFoam/readControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); + + { + // Ensure old-time U exists for mapping + U.oldTime(); + + // Calculate the relative velocity used to map the relative flux phi + volVectorField Urel("Urel", U); + + if (mesh.moving()) + { + Urel -= fvc::reconstruct(fvc::meshPhi(U)); + } + + // Do any mesh changes + mesh.update(); + } + + if (mesh.changing()) + { + Info<< "Execution time for mesh.update() = " + << runTime.elapsedCpuTime() - timeBeforeMeshUpdate + << " s" << endl; + + gh = g & mesh.C(); + ghf = g & mesh.Cf(); + } + + if (mesh.changing() && correctPhi) + { + #include "../interFoam/interDyMFoam/correctPhi.H" + } + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + + twoPhaseProperties->correct(); + + #include "alphaEqnSubCycle.H" + interface.correct(); + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #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/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H new file mode 100644 index 0000000000..afb6478b9d --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H @@ -0,0 +1,77 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phiAbs) + ); + + if (p_rgh.needReference()) + { + fvc::makeRelative(phiHbyA, U); + adjustPhi(phiHbyA, U, p_rgh); + fvc::makeAbsolute(phiHbyA, U); + } + + phiAbs = phiHbyA; + + surfaceScalarField phig + ( + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + Pair > vDotP = twoPhaseProperties->vDotP(); + const volScalarField& vDotcP = vDotP[0](); + const volScalarField& vDotvP = vDotP[1](); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqn + ( + fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh) + - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh) + ); + + p_rghEqn.setReference(pRefCell, pRefValue); + + p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + p_rghEqn.flux(); + + U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf); + U.correctBoundaryConditions(); + fvOptions.correct(U); + } + } + + phiAbs = phi; + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + 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/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files new file mode 100644 index 0000000000..ac0e4f4787 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files @@ -0,0 +1,7 @@ +phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C +phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C +Kunz/Kunz.C +Merkle/Merkle.C +SchnerrSauer/SchnerrSauer.C + +LIB = $(FOAM_LIBBIN)/libphaseChangeTwoPhaseMixtures diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options new file mode 100644 index 0000000000..2e61e4ba77 --- /dev/null +++ b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options @@ -0,0 +1,11 @@ +EXE_INC = \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +LIB_LIBS = \ + -ltwoPhaseMixture \ + -ltwoPhaseProperties \ + -lincompressibleTransportModels \ + -lfiniteVolume diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H index 9be5cd0420..e8ed67a9d8 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/CoBlended/CoBlended.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,11 +28,31 @@ Description Two-scheme Courant number based blending differencing scheme. Similar to localBlended but uses a blending factor computed from the - face-based Courant number and the alpha factor supplied: + face-based Courant number and the lower and upper Courant number limits + supplied: + \f[ + weight = 1 - max(min((Co - Co1)/(Co2 - Co1), 1), 0) + \f] + where + \vartable + Co1 | Courant number below which scheme1 is used + Co2 | Courant number above which scheme2 is used + \endvartable - weight = 1 - Co/alpha + The weight applies to the first scheme and 1-weight to the second scheme. - The weight applies to the first scheme and 1-factor to the second scheme. + Example of the CoBlended scheme specification using LUST for Courant numbers + less than 1 and linearUpwind for Courant numbers greater than 10: + \verbatim + divSchemes + { + . + . + div(phi,U) Gauss CoBlended 1 LUST grad(U) 10 linearUpwind grad(U); + . + . + } + \endverbatim SourceFiles CoBlended.C @@ -60,13 +80,15 @@ class CoBlended { // Private data - const scalar alpha_; - - // Private Member Functions + //- Courant number below which scheme1 is used + const scalar Co1_; //- Scheme 1 tmp > tScheme1_; + //- Courant number above which scheme2 is used + const scalar Co2_; + //- Scheme 2 tmp > tScheme2_; @@ -74,6 +96,8 @@ class CoBlended const surfaceScalarField& faceFlux_; + // Private Member Functions + //- Disallow default bitwise copy construct CoBlended(const CoBlended&); @@ -99,11 +123,12 @@ public: ) : surfaceInterpolationScheme(mesh), - alpha_(readScalar(is)), + Co1_(readScalar(is)), tScheme1_ ( surfaceInterpolationScheme::New(mesh, is) ), + Co2_(readScalar(is)), tScheme2_ ( surfaceInterpolationScheme::New(mesh, is) @@ -116,11 +141,11 @@ public: ) ) { - if (alpha_ <= 0 ) + if (Co1_ < 0 || Co2_ < 0 || Co1_ >= Co2_) { FatalIOErrorIn("CoBlended(const fvMesh&, Istream&)", is) - << "coefficient = " << alpha_ - << " should be > 0" + << "coefficients = " << Co1_ << " and " << Co2_ + << " should be > 0 and Co2 > Co1" << exit(FatalIOError); } } @@ -135,22 +160,23 @@ public: ) : surfaceInterpolationScheme(mesh), - alpha_(readScalar(is)), + Co1_(readScalar(is)), tScheme1_ ( surfaceInterpolationScheme::New(mesh, faceFlux, is) ), + Co2_(readScalar(is)), tScheme2_ ( surfaceInterpolationScheme::New(mesh, faceFlux, is) ), faceFlux_(faceFlux) { - if (alpha_ <= 0) + if (Co1_ < 0 || Co2_ < 0 || Co1_ >= Co2_) { FatalIOErrorIn("CoBlended(const fvMesh&, Istream&)", is) - << "coefficient = " << alpha_ - << " should be > 0" + << "coefficients = " << Co1_ << " and " << Co2_ + << " should be > 0 and Co2 > Co1" << exit(FatalIOError); } } @@ -165,11 +191,18 @@ public: return ( + scalar(1) - max ( - scalar(1) - - mesh.time().deltaT()*mesh.deltaCoeffs()*mag(faceFlux_) - /(mesh.magSf()*alpha_), + min + ( + ( + mesh.time().deltaT()*mesh.deltaCoeffs() + *mag(faceFlux_)/mesh.magSf() + - Co1_ + )/(Co2_ - Co1_), + scalar(1) + ), scalar(0) ) ); @@ -185,6 +218,8 @@ public: { surfaceScalarField bf(blendingFactor()); + Info<< "weights " << max(bf) << " " << min(bf) << endl; + return bf*tScheme1_().weights(vf) + (scalar(1.0) - bf)*tScheme2_().weights(vf);