diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/files b/applications/solvers/incompressible/porousSimpleFoam/Make/files new file mode 100644 index 0000000000..e72f7ed78c --- /dev/null +++ b/applications/solvers/incompressible/porousSimpleFoam/Make/files @@ -0,0 +1,3 @@ +porousSimpleFoam.C + +EXE = $(FOAM_APPBIN)/porousSimpleFoam diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/options b/applications/solvers/incompressible/porousSimpleFoam/Make/options new file mode 100644 index 0000000000..f490133ba6 --- /dev/null +++ b/applications/solvers/incompressible/porousSimpleFoam/Make/options @@ -0,0 +1,14 @@ +EXE_INC = \ + -I../simpleFoam \ + -I$(LIB_SRC)/turbulenceModels \ + -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lincompressibleRASModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H new file mode 100644 index 0000000000..c1925c7708 --- /dev/null +++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H @@ -0,0 +1,47 @@ + // Construct the Momentum equation + + tmp UEqn + ( + fvm::div(phi, U) + - fvm::Sp(fvc::div(phi), U) + + turbulence->divDevReff(U) + ); + + UEqn().relax(); + + // Include the porous media resistance and solve the momentum equation + // either implicit in the tensorial resistance or transport using by + // including the spherical part of the resistance in the momentum diagonal + + tmp trAU; + tmp trTU; + + if (pressureImplicitPorosity) + { + tmp tTU = tensor(I)*UEqn().A(); + pZones.addResistance(UEqn(), tTU()); + trTU = inv(tTU()); + trTU().rename("rAU"); + + volVectorField gradp = fvc::grad(p); + + for (int UCorr=0; UCorr turbulence + ( + incompressible::RASModel::New(U, phi, laminarTransport) + ); + + + porousZones pZones(mesh); + Switch pressureImplicitPorosity(false); + + int nUCorr = 0; + if (pZones.size()) + { + // nUCorrectors for pressureImplicitPorosity + if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors")) + { + nUCorr = readInt + ( + mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors") + ); + } + + if (nUCorr > 0) + { + pressureImplicitPorosity = true; + } + } diff --git a/applications/solvers/incompressible/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H new file mode 100644 index 0000000000..e6797cd860 --- /dev/null +++ b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H @@ -0,0 +1,59 @@ +if (pressureImplicitPorosity) +{ + U = trTU()&UEqn().H(); +} +else +{ + U = trAU()*UEqn().H(); +} + +UEqn.clear(); +phi = fvc::interpolate(U) & mesh.Sf(); +adjustPhi(phi, U, p); + +for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) +{ + tmp tpEqn; + + if (pressureImplicitPorosity) + { + tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi)); + } + else + { + tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi)); + } + + tpEqn().setReference(pRefCell, pRefValue); + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = tpEqn().solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + tpEqn().solve(); + } + + if (nonOrth == nNonOrthCorr) + { + phi -= tpEqn().flux(); + } +} + +#include "continuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +if (pressureImplicitPorosity) +{ + U -= trTU()&fvc::grad(p); +} +else +{ + U -= trAU()*fvc::grad(p); +} + +U.correctBoundaryConditions(); diff --git a/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C new file mode 100644 index 0000000000..b78f1688fe --- /dev/null +++ b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C @@ -0,0 +1,85 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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 + +Application + porousSimpleFoam + +Description + Steady-state solver for incompressible, turbulent flow with + implicit or explicit porosity treatment + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" +#include "porousZones.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createFields.H" + #include "initContinuityErrs.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + #include "readSIMPLEControls.H" + #include "initConvergenceCheck.H" + + p.storePrevIter(); + + // Pressure-velocity SIMPLE corrector + { + #include "UEqn.H" + #include "pEqn.H" + } + + turbulence->correct(); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + #include "convergenceCheck.H" + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/Make/files b/applications/solvers/multiphase/interMixingFoam/Make/files new file mode 100644 index 0000000000..488cd77e56 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/Make/files @@ -0,0 +1,6 @@ +incompressibleThreePhaseMixture/threePhaseMixture.C +threePhaseInterfaceProperties/threePhaseInterfaceProperties.C +interMixingFoam.C + +EXE = $(FOAM_APPBIN)/interMixingFoam + diff --git a/applications/solvers/multiphase/interMixingFoam/Make/options b/applications/solvers/multiphase/interMixingFoam/Make/options new file mode 100644 index 0000000000..d8e4da2313 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/Make/options @@ -0,0 +1,17 @@ +INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam + +EXE_INC = \ + -I$(INTERFOAM) \ + -IincompressibleThreePhaseMixture \ + -IthreePhaseInterfaceProperties \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/transportModels + +EXE_LIBS = \ + -linterfaceProperties \ + -lincompressibleTransportModels \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H new file mode 100644 index 0000000000..f9bad0a705 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H @@ -0,0 +1,164 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phir + ( + IOobject + ( + "phir", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf() + ); + + for (int gCorr=0; gCorr(mesh, phi).flux(alpha1); + + // Calculate the flux correction for alpha1 + phiAlpha1 -= phiAlpha1BD; + + // Calculate the limiter for alpha1 + MULES::limiter + ( + allLambda, + oneField(), + alpha1, + phiAlpha1BD, + phiAlpha1, + zeroField(), + zeroField(), + 1, + 0, + 3 + ); + + // Create the complete flux for alpha2 + surfaceScalarField phiAlpha2 = + fvc::flux + ( + phi, + alpha2, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(phir, alpha1, alpharScheme), + alpha2, + alpharScheme + ); + + // Create the bounded (upwind) flux for alpha2 + surfaceScalarField phiAlpha2BD = + upwind(mesh, phi).flux(alpha2); + + // Calculate the flux correction for alpha2 + phiAlpha2 -= phiAlpha2BD; + + // Further limit the limiter for alpha2 + MULES::limiter + ( + allLambda, + oneField(), + alpha2, + phiAlpha2BD, + phiAlpha2, + zeroField(), + zeroField(), + 1, + 0, + 3 + ); + + // Construct the limited fluxes + phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1; + phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2; + + // Solve for alpha1 + solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1)); + + // Create the diffusion coefficients for alpha2<->alpha3 + volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2); + volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3); + + // Add the diffusive flux for alpha3->alpha2 + phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1); + + // Solve for alpha2 + fvScalarMatrix alpha2Eqn + ( + fvm::ddt(alpha2) + + fvc::div(phiAlpha2) + - fvm::laplacian(Dc23 + Dc32, alpha2) + ); + alpha2Eqn.solve(); + + // Construct the complete mass flux + rhoPhi = + phiAlpha1*(rho1 - rho3) + + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3) + + phi*rho3; + + alpha3 = 1.0 - alpha1 - alpha2; + } + + Info<< "Air phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + + Info<< "Liquid phase volume fraction = " + << alpha2.weightedAverage(mesh.V()).value() + << " Min(alpha2) = " << min(alpha2).value() + << " Max(alpha2) = " << max(alpha2).value() + << endl; +} diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H new file mode 100644 index 0000000000..765087a183 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H @@ -0,0 +1,43 @@ +label nAlphaCorr +( + readLabel(piso.lookup("nAlphaCorr")) +); + +label nAlphaSubCycles +( + readLabel(piso.lookup("nAlphaSubCycles")) +); + +if (nAlphaSubCycles > 1) +{ + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + dimensionedScalar totalDeltaT = runTime.deltaT(); + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { +# include "alphaEqns.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ +# include "alphaEqns.H" +} + +interface.correct(); + +{ + volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3; + + //solve(fvm::ddt(rho) + fvc::div(rhoPhi)); + //Info<< "density error = " + // << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl; + + rho == rhoNew; +} diff --git a/applications/solvers/multiphase/interMixingFoam/createFields.H b/applications/solvers/multiphase/interMixingFoam/createFields.H new file mode 100644 index 0000000000..9615e6c186 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/createFields.H @@ -0,0 +1,132 @@ + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 + ( + IOobject + ( + "alpha1", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Reading field alpha2\n" << endl; + volScalarField alpha2 + ( + IOobject + ( + "alpha2", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Reading field alpha3\n" << endl; + volScalarField alpha3 + ( + IOobject + ( + "alpha3", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + alpha3 == 1.0 - alpha1 - alpha2; + + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +# include "createPhi.H" + + threePhaseMixture threePhaseProperties(U, phi); + + const dimensionedScalar& rho1 = threePhaseProperties.rho1(); + const dimensionedScalar& rho2 = threePhaseProperties.rho2(); + const dimensionedScalar& rho3 = threePhaseProperties.rho3(); + + dimensionedScalar D23(threePhaseProperties.lookup("D23")); + + // Need to store rho for ddt(rho, U) + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + alpha1*rho1 + alpha2*rho2 + alpha3*rho3, + alpha1.boundaryField().types() + ); + rho.oldTime(); + + + // Mass flux + // Initialisation does not matter because rhoPhi is reset after the + // alpha solution before it is used in the U equation. + surfaceScalarField rhoPhi + ( + IOobject + ( + "rho*phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho1*phi + ); + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); + + + // Construct interface from alpha distribution + threePhaseInterfaceProperties interface(threePhaseProperties); + + + // Construct incompressible turbulence model + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, threePhaseProperties) + ); diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C new file mode 100644 index 0000000000..f27e17923f --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C @@ -0,0 +1,204 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseMixture + +\*---------------------------------------------------------------------------*/ + +#include "threePhaseMixture.H" +#include "addToRunTimeSelectionTable.H" +#include "surfaceFields.H" +#include "fvc.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +//- Calculate and return the laminar viscosity +void Foam::threePhaseMixture::calcNu() +{ + // Average kinematic viscosity calculated from dynamic viscosity + nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::threePhaseMixture::threePhaseMixture +( + const volVectorField& U, + const surfaceScalarField& phi +) +: + transportModel(U, phi), + + phase1Name_("phase1"), + phase2Name_("phase2"), + phase3Name_("phase3"), + + nuModel1_ + ( + viscosityModel::New + ( + "nu1", + subDict(phase1Name_), + U, + phi + ) + ), + nuModel2_ + ( + viscosityModel::New + ( + "nu2", + subDict(phase2Name_), + U, + phi + ) + ), + nuModel3_ + ( + viscosityModel::New + ( + "nu3", + subDict(phase2Name_), + U, + phi + ) + ), + + rho1_(nuModel1_->viscosityProperties().lookup("rho")), + rho2_(nuModel2_->viscosityProperties().lookup("rho")), + rho3_(nuModel3_->viscosityProperties().lookup("rho")), + + U_(U), + phi_(phi), + + alpha1_(U_.db().lookupObject ("alpha1")), + alpha2_(U_.db().lookupObject ("alpha2")), + alpha3_(U_.db().lookupObject ("alpha3")), + + nu_ + ( + IOobject + ( + "nu", + U_.time().timeName(), + U_.db() + ), + U_.mesh(), + dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0), + calculatedFvPatchScalarField::typeName + ) +{ + calcNu(); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::threePhaseMixture::mu() const +{ + return tmp + ( + new volScalarField + ( + "mu", + alpha1_*rho1_*nuModel1_->nu() + + alpha2_*rho2_*nuModel2_->nu() + + alpha3_*rho3_*nuModel3_->nu() + ) + ); +} + + +Foam::tmp Foam::threePhaseMixture::muf() const +{ + surfaceScalarField alpha1f = fvc::interpolate(alpha1_); + surfaceScalarField alpha2f = fvc::interpolate(alpha2_); + surfaceScalarField alpha3f = fvc::interpolate(alpha3_); + + return tmp + ( + new surfaceScalarField + ( + "mu", + alpha1f*rho1_*fvc::interpolate(nuModel1_->nu()) + + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu()) + + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu()) + ) + ); +} + + +Foam::tmp Foam::threePhaseMixture::nuf() const +{ + surfaceScalarField alpha1f = fvc::interpolate(alpha1_); + surfaceScalarField alpha2f = fvc::interpolate(alpha2_); + surfaceScalarField alpha3f = fvc::interpolate(alpha3_); + + return tmp + ( + new surfaceScalarField + ( + "nu", + ( + alpha1f*rho1_*fvc::interpolate(nuModel1_->nu()) + + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu()) + + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu()) + )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_) + ) + ); +} + + +bool Foam::threePhaseMixture::read() +{ + if (transportModel::read()) + { + if + ( + nuModel1_().read(*this) + && nuModel2_().read(*this) + && nuModel3_().read(*this) + ) + { + nuModel1_->viscosityProperties().lookup("rho") >> rho1_; + nuModel2_->viscosityProperties().lookup("rho") >> rho2_; + nuModel3_->viscosityProperties().lookup("rho") >> rho3_; + + return true; + } + else + { + return false; + } + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H new file mode 100644 index 0000000000..2712b5d0cb --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H @@ -0,0 +1,197 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseMixture + +Description + +SourceFiles + threePhaseMixture.C + +\*---------------------------------------------------------------------------*/ + +#ifndef threePhaseMixture_H +#define threePhaseMixture_H + +#include "incompressible/transportModel/transportModel.H" +#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H" +#include "dimensionedScalar.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class threePhaseMixture Declaration +\*---------------------------------------------------------------------------*/ + +class threePhaseMixture +: + public transportModel +{ + // Private data + + word phase1Name_; + word phase2Name_; + word phase3Name_; + + autoPtr nuModel1_; + autoPtr nuModel2_; + autoPtr nuModel3_; + + dimensionedScalar rho1_; + dimensionedScalar rho2_; + dimensionedScalar rho3_; + + const volVectorField& U_; + const surfaceScalarField& phi_; + + const volScalarField& alpha1_; + const volScalarField& alpha2_; + const volScalarField& alpha3_; + + volScalarField nu_; + + + // Private Member Functions + + //- Calculate and return the laminar viscosity + void calcNu(); + + +public: + + // Constructors + + //- Construct from components + threePhaseMixture + ( + const volVectorField& U, + const surfaceScalarField& phi + ); + + + // Destructor + + ~threePhaseMixture() + {} + + + // Member Functions + + //- Return const-access to phase1 viscosityModel + const viscosityModel& nuModel1() const + { + return nuModel1_(); + } + + //- Return const-access to phase2 viscosityModel + const viscosityModel& nuModel2() const + { + return nuModel2_(); + } + + //- Return const-access to phase3 viscosityModel + const viscosityModel& nuModel3() const + { + return nuModel3_(); + } + + //- Return const-access to phase1 density + const dimensionedScalar& rho1() const + { + return rho1_; + } + + //- Return const-access to phase2 density + const dimensionedScalar& rho2() const + { + return rho2_; + }; + + //- Return const-access to phase3 density + const dimensionedScalar& rho3() const + { + return rho3_; + }; + + const volScalarField& alpha1() const + { + return alpha1_; + } + + const volScalarField& alpha2() const + { + return alpha2_; + } + + const volScalarField& alpha3() const + { + return alpha3_; + } + + //- Return the velocity + const volVectorField& U() const + { + return U_; + } + + //- Return the dynamic laminar viscosity + tmp mu() const; + + //- Return the face-interpolated dynamic laminar viscosity + tmp muf() const; + + //- Return the kinematic laminar viscosity + tmp nu() const + { + return nu_; + } + + //- Return the face-interpolated dynamic laminar viscosity + tmp nuf() const; + + //- Correct the laminar viscosity + void correct() + { + calcNu(); + } + + //- Read base transportProperties dictionary + bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C new file mode 100644 index 0000000000..a18df8d12a --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Application + interMixingFoam + +Description + Solver for 3 incompressible fluids, two of which are miscible, + using a VOF method to capture the interface. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "MULES.H" +#include "subCycle.H" +#include "threePhaseMixture.H" +#include "threePhaseInterfaceProperties.H" +#include "turbulenceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + #include "readPISOControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + #include "correctPhi.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readPISOControls.H" + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + threePhaseProperties.correct(); + + #include "alphaEqnsSubCycle.H" + + #define twoPhaseProperties threePhaseProperties + #include "UEqn.H" + + // --- PISO loop + for (int corr=0; corrcorrect(); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "\n end \n"; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C new file mode 100644 index 0000000000..1ab10bd743 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C @@ -0,0 +1,209 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Application + threePhaseInterfaceProperties + +Description + Properties to aid interFoam : + 1. Correct the alpha boundary condition for dynamic contact angle. + 2. Calculate interface curvature. + +\*---------------------------------------------------------------------------*/ + +#include "threePhaseInterfaceProperties.H" +#include "alphaContactAngleFvPatchScalarField.H" +#include "mathematicalConstants.H" +#include "surfaceInterpolate.H" +#include "fvcDiv.H" +#include "fvcGrad.H" + +// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * // + +const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad = + Foam::constant::mathematical::pi/180.0; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Correction for the boundary condition on the unit normal nHat on +// walls to produce the correct contact angle. + +// The dynamic contact angle is calculated from the component of the +// velocity on the direction of the interface, parallel to the wall. + +void Foam::threePhaseInterfaceProperties::correctContactAngle +( + surfaceVectorField::GeometricBoundaryField& nHatb +) const +{ + const volScalarField::GeometricBoundaryField& alpha1 = + mixture_.alpha1().boundaryField(); + const volScalarField::GeometricBoundaryField& alpha2 = + mixture_.alpha2().boundaryField(); + const volScalarField::GeometricBoundaryField& alpha3 = + mixture_.alpha3().boundaryField(); + const volVectorField::GeometricBoundaryField& U = + mixture_.U().boundaryField(); + + const fvMesh& mesh = mixture_.U().mesh(); + const fvBoundaryMesh& boundary = mesh.boundary(); + + forAll(boundary, patchi) + { + if (isA(alpha1[patchi])) + { + const alphaContactAngleFvPatchScalarField& a2cap = + refCast + (alpha2[patchi]); + + const alphaContactAngleFvPatchScalarField& a3cap = + refCast + (alpha3[patchi]); + + scalarField twoPhaseAlpha2 = max(a2cap, scalar(0)); + scalarField twoPhaseAlpha3 = max(a3cap, scalar(0)); + + scalarField sumTwoPhaseAlpha = + twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL; + + twoPhaseAlpha2 /= sumTwoPhaseAlpha; + twoPhaseAlpha3 /= sumTwoPhaseAlpha; + + fvsPatchVectorField& nHatp = nHatb[patchi]; + + scalarField theta = + convertToRad + *( + twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp)) + + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp)) + ); + + vectorField nf = boundary[patchi].nf(); + + // Reset nHatPatch to correspond to the contact angle + + scalarField a12 = nHatp & nf; + + scalarField b1 = cos(theta); + + scalarField b2(nHatp.size()); + + forAll(b2, facei) + { + b2[facei] = cos(acos(a12[facei]) - theta[facei]); + } + + scalarField det = 1.0 - a12*a12; + + scalarField a = (b1 - a12*b2)/det; + scalarField b = (b2 - a12*b1)/det; + + nHatp = a*nf + b*nHatp; + + nHatp /= (mag(nHatp) + deltaN_.value()); + } + } +} + + +void Foam::threePhaseInterfaceProperties::calculateK() +{ + const volScalarField& alpha1 = mixture_.alpha1(); + + const fvMesh& mesh = alpha1.mesh(); + const surfaceVectorField& Sf = mesh.Sf(); + + // Cell gradient of alpha + volVectorField gradAlpha = fvc::grad(alpha1); + + // Interpolated face-gradient of alpha + surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha); + + // Face unit interface normal + surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_); + correctContactAngle(nHatfv.boundaryField()); + + // Face unit interface normal flux + nHatf_ = nHatfv & Sf; + + // Simple expression for curvature + K_ = -fvc::div(nHatf_); + + // Complex expression for curvature. + // Correction is formally zero but numerically non-zero. + //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_); + //nHat.boundaryField() = nHatfv.boundaryField(); + //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties +( + const threePhaseMixture& mixture +) +: + mixture_(mixture), + cAlpha_ + ( + readScalar + ( + mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha") + ) + ), + sigma12_(mixture.lookup("sigma12")), + sigma13_(mixture.lookup("sigma13")), + + deltaN_ + ( + "deltaN", + 1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0) + ), + + nHatf_ + ( + ( + fvc::interpolate(fvc::grad(mixture.alpha1())) + /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_) + ) & mixture.alpha1().mesh().Sf() + ), + + K_ + ( + IOobject + ( + "K", + mixture.alpha1().time().timeName(), + mixture.alpha1().mesh() + ), + -fvc::div(nHatf_) + ) +{ + calculateK(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H new file mode 100644 index 0000000000..3a402feeed --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseInterfaceProperties + +Description + Properties to aid interFoam : + 1. Correct the alpha boundary condition for dynamic contact angle. + 2. Calculate interface curvature. + +SourceFiles + threePhaseInterfaceProperties.C + +\*---------------------------------------------------------------------------*/ + +#ifndef threePhaseInterfaceProperties_H +#define threePhaseInterfaceProperties_H + +#include "threePhaseMixture.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class threePhaseInterfaceProperties Declaration +\*---------------------------------------------------------------------------*/ + +class threePhaseInterfaceProperties +{ + // Private data + + const threePhaseMixture& mixture_; + + //- Compression coefficient + scalar cAlpha_; + + //- Surface tension 1-2 + dimensionedScalar sigma12_; + + //- Surface tension 1-3 + dimensionedScalar sigma13_; + + //- Stabilisation for normalisation of the interface normal + const dimensionedScalar deltaN_; + + surfaceScalarField nHatf_; + volScalarField K_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct and assignment + threePhaseInterfaceProperties(const threePhaseInterfaceProperties&); + void operator=(const threePhaseInterfaceProperties&); + + //- Correction for the boundary condition on the unit normal nHat on + // walls to produce the correct contact dynamic angle + // calculated from the component of U parallel to the wall + void correctContactAngle + ( + surfaceVectorField::GeometricBoundaryField& nHat + ) const; + + //- Re-calculate the interface curvature + void calculateK(); + + +public: + + //- Conversion factor for degrees into radians + static const scalar convertToRad; + + + // Constructors + + //- Construct from volume fraction field alpha and IOdictionary + threePhaseInterfaceProperties(const threePhaseMixture& mixture); + + + // Member Functions + + scalar cAlpha() const + { + return cAlpha_; + } + + const dimensionedScalar& deltaN() const + { + return deltaN_; + } + + const surfaceScalarField& nHatf() const + { + return nHatf_; + } + + const volScalarField& K() const + { + return K_; + } + + tmp sigma() const + { + volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0)); + volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0)); + + return + (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_) + /(limitedAlpha2 + limitedAlpha3 + SMALL); + } + + tmp sigmaK() const + { + return sigma()*K_; + } + + void correct() + { + calculateK(); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/test/CompactListList/CompactListListTest.C b/applications/test/CompactListList/CompactListListTest.C index 136b5f4583..cd3270a169 100644 --- a/applications/test/CompactListList/CompactListListTest.C +++ b/applications/test/CompactListList/CompactListListTest.C @@ -32,6 +32,9 @@ Description #include "CompactListList.H" #include "IOstreams.H" +#include "OStringStream.H" +#include "IStringStream.H" +#include "faceList.H" using namespace Foam; @@ -40,7 +43,29 @@ using namespace Foam; int main(int argc, char *argv[]) { - CompactListList