diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean new file mode 100755 index 0000000000..8098c05160 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean libso multiphaseMixtureThermo +wclean + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake new file mode 100755 index 0000000000..dd002ee069 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake libso multiphaseMixtureThermo +wmake + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files new file mode 100644 index 0000000000..8e6a4cb0de --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files @@ -0,0 +1,3 @@ +compressibleMultiphaseInterFoam.C + +EXE = $(FOAM_APPBIN)/compressibleMultiphaseInterFoam diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options new file mode 100644 index 0000000000..0f47e6979e --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options @@ -0,0 +1,17 @@ +EXE_INC = \ + -ImultiphaseMixtureThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lmultiphaseMixtureThermo \ + -lfluidThermophysicalModels \ + -lspecie \ + -linterfaceProperties \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H new file mode 100644 index 0000000000..a7ee3a7bb4 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H @@ -0,0 +1,17 @@ +{ + fvScalarMatrix TEqn + ( + fvm::ddt(rho, T) + + fvm::div(multiphaseProperties.rhoPhi(), T) + - fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T) + + ( + fvc::div(fvc::absolute(phi, U), p) + + fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K) + )*multiphaseProperties.rCv() + ); + + TEqn.relax(); + TEqn.solve(); + + multiphaseProperties.correct(); +} diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H new file mode 100644 index 0000000000..38cfebde6f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H @@ -0,0 +1,27 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(multiphaseProperties.rhoPhi(), U) + + turbulence->divDevRhoReff(U) + ); + + UEqn.relax(); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + multiphaseProperties.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + K = 0.5*magSqr(U); + } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H new file mode 100644 index 0000000000..1fd68ffb79 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Global + CourantNo + +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 + ( + multiphaseProperties.nearInterface()().internalField() + * fvc::surfaceSum(mag(phi))().internalField() + ); + + 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/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C new file mode 100644 index 0000000000..09018b603f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C @@ -0,0 +1,111 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 + compressibleMultiphaseInterFoam + +Description + Solver for n compressible, non-isothermal immiscible fluids using a VOF + (volume of fluid) 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. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "multiphaseMixtureThermo.H" +#include "turbulenceModel.H" +#include "pimpleControl.H" +#include "fixedFluxPressureFvPatchScalarField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + + pimpleControl pimple(mesh); + + #include "readTimeControls.H" + #include "initContinuityErrs.H" + #include "createFields.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()) + { + multiphaseProperties.solve(); + + solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi())); + + #include "UEqn.H" + #include "TEqn.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/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H new file mode 100644 index 0000000000..5c7f723a87 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H @@ -0,0 +1,71 @@ + 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<< "Constructing multiphaseMixtureThermo\n" << endl; + multiphaseMixtureThermo multiphaseProperties(U, phi); + + volScalarField& p = multiphaseProperties.p(); + volScalarField& T = multiphaseProperties.T(); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + multiphaseProperties.rho() + ); + + Info<< max(rho) << min(rho); + + dimensionedScalar pMin(multiphaseProperties.lookup("pMin")); + + + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + // Construct compressible turbulence model + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + multiphaseProperties.rhoPhi(), + multiphaseProperties + ) + ); + + Info<< "Creating field kinetic energy K\n" << endl; + volScalarField K("K", 0.5*magSqr(U)); diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files new file mode 100644 index 0000000000..93e6eb9e3d --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files @@ -0,0 +1,5 @@ +phaseModel/phaseModel.C +alphaContactAngle/alphaContactAngleFvPatchScalarField.C +multiphaseMixtureThermo.C + +LIB = $(FOAM_LIBBIN)/libmultiphaseMixtureThermo diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options new file mode 100644 index 0000000000..eab8cce15d --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options @@ -0,0 +1,8 @@ +EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +LIB_LIBS = \ + -lfluidThermophysicalModels \ + -lspecie \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C new file mode 100644 index 0000000000..6872ae0321 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "alphaContactAngleFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps +( + Istream& is +) +: + theta0_(readScalar(is)), + uTheta_(readScalar(is)), + thetaA_(readScalar(is)), + thetaR_(readScalar(is)) +{} + + +Istream& operator>> +( + Istream& is, + alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp +) +{ + is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_; + return is; +} + + +Ostream& operator<< +( + Ostream& os, + const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp +) +{ + os << tp.theta0_ << token::SPACE + << tp.uTheta_ << token::SPACE + << tp.thetaA_ << token::SPACE + << tp.thetaR_; + + return os; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + zeroGradientFvPatchScalarField(p, iF) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper), + thetaProps_(gcpsf.thetaProps_) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + zeroGradientFvPatchScalarField(p, iF), + thetaProps_(dict.lookup("thetaProperties")) +{ + evaluate(); +} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const DimensionedField& iF +) +: + zeroGradientFvPatchScalarField(gcpsf, iF), + thetaProps_(gcpsf.thetaProps_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void alphaContactAngleFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + os.writeKeyword("thetaProperties") + << thetaProps_ << token::END_STATEMENT << nl; + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + alphaContactAngleFvPatchScalarField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H new file mode 100644 index 0000000000..53ce61d24f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H @@ -0,0 +1,215 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::alphaContactAngleFvPatchScalarField + +Description + Contact-angle boundary condition for multi-phase interface-capturing + simulations. Used in conjuction with multiphaseMixture. + +SourceFiles + alphaContactAngleFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef alphaContactAngleFvPatchScalarField_H +#define alphaContactAngleFvPatchScalarField_H + +#include "zeroGradientFvPatchFields.H" +#include "multiphaseMixtureThermo.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class alphaContactAngleFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class alphaContactAngleFvPatchScalarField +: + public zeroGradientFvPatchScalarField +{ +public: + + class interfaceThetaProps + { + //- Equilibrium contact angle + scalar theta0_; + + //- Dynamic contact angle velocity scale + scalar uTheta_; + + //- Limiting advancing contact angle + scalar thetaA_; + + //- Limiting receeding contact angle + scalar thetaR_; + + + public: + + // Constructors + interfaceThetaProps() + {} + + interfaceThetaProps(Istream&); + + + // Member functions + + //- Return the equilibrium contact angle theta0 + scalar theta0(bool matched=true) const + { + if (matched) return theta0_; + else return 180.0 - theta0_; + } + + //- Return the dynamic contact angle velocity scale + scalar uTheta() const + { + return uTheta_; + } + + //- Return the limiting advancing contact angle + scalar thetaA(bool matched=true) const + { + if (matched) return thetaA_; + else return 180.0 - thetaA_; + } + + //- Return the limiting receeding contact angle + scalar thetaR(bool matched=true) const + { + if (matched) return thetaR_; + else return 180.0 - thetaR_; + } + + + // IO functions + + friend Istream& operator>>(Istream&, interfaceThetaProps&); + friend Ostream& operator<<(Ostream&, const interfaceThetaProps&); + }; + + typedef HashTable + < + interfaceThetaProps, + multiphaseMixtureThermo::interfacePair, + multiphaseMixtureThermo::interfacePair::hash + > thetaPropsTable; + + +private: + + // Private data + + thetaPropsTable thetaProps_; + + +public: + + //- Runtime type information + TypeName("alphaContactAngle"); + + + // Constructors + + //- Construct from patch and internal field + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given alphaContactAngleFvPatchScalarField + // onto a new patch + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new alphaContactAngleFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new alphaContactAngleFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Return the contact angle properties + const thetaPropsTable& thetaProps() const + { + return thetaProps_; + } + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C new file mode 100644 index 0000000000..1639550564 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -0,0 +1,1129 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "multiphaseMixtureThermo.H" +#include "alphaContactAngleFvPatchScalarField.H" +#include "Time.H" +#include "subCycle.H" +#include "MULES.H" +#include "fvcDiv.H" +#include "fvcGrad.H" +#include "fvcSnGrad.H" +#include "fvcFlux.H" +#include "fvcMeshPhi.H" +#include "surfaceInterpolate.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(multiphaseMixtureThermo, 0); +} + + +const Foam::scalar Foam::multiphaseMixtureThermo::convertToRad = + Foam::constant::mathematical::pi/180.0; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::multiphaseMixtureThermo::calcAlphas() +{ + scalar level = 0.0; + alphas_ == 0.0; + + forAllIter(PtrDictionary, phases_, phase) + { + alphas_ += level*phase(); + level += 1.0; + } + + alphas_.correctBoundaryConditions(); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::multiphaseMixtureThermo::multiphaseMixtureThermo +( + const volVectorField& U, + const surfaceScalarField& phi +) +: + psiThermo(U.mesh(), word::null), + phases_(lookup("phases"), phaseModel::iNew(p_, T_)), + + mesh_(U.mesh()), + U_(U), + phi_(phi), + + rhoPhi_ + ( + IOobject + ( + "rho*phi", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("rho*phi", dimMass/dimTime, 0.0) + ), + + alphas_ + ( + IOobject + ( + "alphas", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh_, + dimensionedScalar("alphas", dimless, 0.0), + zeroGradientFvPatchScalarField::typeName + ), + + sigmas_(lookup("sigmas")), + dimSigma_(1, 0, -2, 0, 0), + deltaN_ + ( + "deltaN", + 1e-8/pow(average(mesh_.V()), 1.0/3.0) + ) +{ + calcAlphas(); + alphas_.write(); + correct(); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::multiphaseMixtureThermo::correct() +{ + forAllIter(PtrDictionary, phases_, phasei) + { + phasei().correct(); + } + + PtrDictionary::iterator phasei = phases_.begin(); + + psi_ = phasei()*phasei().thermo().psi(); + mu_ = phasei()*phasei().thermo().mu(); + alpha_ = phasei()*phasei().thermo().alpha(); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + psi_ += phasei()*phasei().thermo().psi(); + mu_ += phasei()*phasei().thermo().mu(); + alpha_ += phasei()*phasei().thermo().alpha(); + } +} + + +void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp) +{ + forAllIter(PtrDictionary, phases_, phasei) + { + phasei().thermo().rho() += phasei().thermo().psi()*dp; + } +} + + +bool Foam::multiphaseMixtureThermo::incompressible() const +{ + bool ico = true; + + forAllConstIter(PtrDictionary, phases_, phase) + { + ico &= phase().thermo().incompressible(); + } + + return ico; +} + + +bool Foam::multiphaseMixtureThermo::isochoric() const +{ + bool iso = true; + + forAllConstIter(PtrDictionary, phases_, phase) + { + iso &= phase().thermo().incompressible(); + } + + return iso; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const volScalarField& p, + const volScalarField& T +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the(phasei()*phasei().thermo().he(p, T)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += phasei()*phasei().thermo().he(p, T); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const labelList& cells +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the + ( + scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the + ( + phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += + phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::hc() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp thc(phasei()*phasei().thermo().hc()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + thc() += phasei()*phasei().thermo().hc(); + } + + return thc; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const labelList& cells +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const label patchi +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::rho() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp trho(phasei()*phasei().thermo().rho()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trho() += phasei()*phasei().thermo().rho(); + } + + return trho; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cp() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCp(phasei()*phasei().thermo().Cp()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCp() += phasei()*phasei().thermo().Cp(); + } + + return tCp; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cp +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCp + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCp() += + phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi); + } + + return tCp; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCv(phasei()*phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCv() += phasei()*phasei().thermo().Cv(); + } + + return tCv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCv + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCv() += + phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi); + } + + return tCv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::gamma() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tgamma(phasei()*phasei().thermo().gamma()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tgamma() += phasei()*phasei().thermo().gamma(); + } + + return tgamma; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::gamma +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tgamma + ( + phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tgamma() += + phasei().boundaryField()[patchi] + *phasei().thermo().gamma(p, T, patchi); + } + + return tgamma; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cpv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpv(phasei()*phasei().thermo().Cpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpv() += phasei()*phasei().thermo().Cpv(); + } + + return tCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpv + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpv() += + phasei().boundaryField()[patchi] + *phasei().thermo().Cpv(p, T, patchi); + } + + return tCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpByCpv(phasei()*phasei().thermo().CpByCpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpByCpv() += phasei()*phasei().thermo().CpByCpv(); + } + + return tCpByCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpByCpv + ( + phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpByCpv() += + phasei().boundaryField()[patchi] + *phasei().thermo().CpByCpv(p, T, patchi); + } + + return tCpByCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappa() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappa(phasei()*phasei().thermo().kappa()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappa() += phasei()*phasei().thermo().kappa(); + } + + return tkappa; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappa +( + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappa + ( + phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappa() += + phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi); + } + + return tkappa; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappaEff(phasei()*phasei().thermo().kappaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat); + } + + return tkappaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappaEff + ( + phasei().boundaryField()[patchi] + *phasei().thermo().kappaEff(alphat, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappaEff() += + phasei().boundaryField()[patchi] + *phasei().thermo().kappaEff(alphat, patchi); + } + + return tkappaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::alphaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp talphaEff(phasei()*phasei().thermo().alphaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + talphaEff() += phasei()*phasei().thermo().alphaEff(alphat); + } + + return talphaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::alphaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp talphaEff + ( + phasei().boundaryField()[patchi] + *phasei().thermo().alphaEff(alphat, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + talphaEff() += + phasei().boundaryField()[patchi] + *phasei().thermo().alphaEff(alphat, patchi); + } + + return talphaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::rCv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp trCv(phasei()/phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trCv() += phasei()/phasei().thermo().Cv(); + } + + return trCv; +} + + +Foam::tmp +Foam::multiphaseMixtureThermo::surfaceTensionForce() const +{ + tmp tstf + ( + new surfaceScalarField + ( + IOobject + ( + "surfaceTensionForce", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar + ( + "surfaceTensionForce", + dimensionSet(1, -2, -2, 0, 0), + 0.0 + ) + ) + ); + + surfaceScalarField& stf = tstf(); + + forAllConstIter(PtrDictionary, phases_, phase1) + { + const phaseModel& alpha1 = phase1(); + + PtrDictionary::const_iterator phase2 = phase1; + ++phase2; + + for (; phase2 != phases_.end(); ++phase2) + { + const phaseModel& alpha2 = phase2(); + + sigmaTable::const_iterator sigma = + sigmas_.find(interfacePair(alpha1, alpha2)); + + if (sigma == sigmas_.end()) + { + FatalErrorIn + ( + "multiphaseMixtureThermo::surfaceTensionForce() const" + ) << "Cannot find interface " << interfacePair(alpha1, alpha2) + << " in list of sigma values" + << exit(FatalError); + } + + stf += dimensionedScalar("sigma", dimSigma_, sigma()) + *fvc::interpolate(K(alpha1, alpha2))* + ( + fvc::interpolate(alpha2)*fvc::snGrad(alpha1) + - fvc::interpolate(alpha1)*fvc::snGrad(alpha2) + ); + } + } + + return tstf; +} + + +void Foam::multiphaseMixtureThermo::solve() +{ + const Time& runTime = mesh_.time(); + + const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE"); + + label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); + + scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha"))); + + + volScalarField& alpha = phases_.first(); + + if (nAlphaSubCycles > 1) + { + surfaceScalarField rhoPhiSum(0.0*rhoPhi_); + dimensionedScalar totalDeltaT = runTime.deltaT(); + + for + ( + subCycle alphaSubCycle(alpha, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + solveAlphas(cAlpha); + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_; + } + + rhoPhi_ = rhoPhiSum; + } + else + { + solveAlphas(cAlpha); + } +} + + +Foam::tmp Foam::multiphaseMixtureThermo::nHatfv +( + const volScalarField& alpha1, + const volScalarField& alpha2 +) const +{ + /* + // Cell gradient of alpha + volVectorField gradAlpha = + alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2); + + // Interpolated face-gradient of alpha + surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha); + */ + + surfaceVectorField gradAlphaf + ( + fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1)) + - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2)) + ); + + // Face unit interface normal + return gradAlphaf/(mag(gradAlphaf) + deltaN_); +} + + +Foam::tmp Foam::multiphaseMixtureThermo::nHatf +( + const volScalarField& alpha1, + const volScalarField& alpha2 +) const +{ + // Face unit interface normal flux + return nHatfv(alpha1, alpha2) & mesh_.Sf(); +} + + +// 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::multiphaseMixtureThermo::correctContactAngle +( + const phaseModel& alpha1, + const phaseModel& alpha2, + surfaceVectorField::GeometricBoundaryField& nHatb +) const +{ + const volScalarField::GeometricBoundaryField& gbf + = alpha1.boundaryField(); + + const fvBoundaryMesh& boundary = mesh_.boundary(); + + forAll(boundary, patchi) + { + if (isA(gbf[patchi])) + { + const alphaContactAngleFvPatchScalarField& acap = + refCast(gbf[patchi]); + + vectorField& nHatPatch = nHatb[patchi]; + + vectorField AfHatPatch + ( + mesh_.Sf().boundaryField()[patchi] + /mesh_.magSf().boundaryField()[patchi] + ); + + alphaContactAngleFvPatchScalarField::thetaPropsTable:: + const_iterator tp = + acap.thetaProps().find(interfacePair(alpha1, alpha2)); + + if (tp == acap.thetaProps().end()) + { + FatalErrorIn + ( + "multiphaseMixtureThermo::correctContactAngle" + "(const phaseModel& alpha1, const phaseModel& alpha2, " + "fvPatchVectorFieldField& nHatb) const" + ) << "Cannot find interface " << interfacePair(alpha1, alpha2) + << "\n in table of theta properties for patch " + << acap.patch().name() + << exit(FatalError); + } + + bool matched = (tp.key().first() == alpha1.name()); + + scalar theta0 = convertToRad*tp().theta0(matched); + scalarField theta(boundary[patchi].size(), theta0); + + scalar uTheta = tp().uTheta(); + + // Calculate the dynamic contact angle if required + if (uTheta > SMALL) + { + scalar thetaA = convertToRad*tp().thetaA(matched); + scalar thetaR = convertToRad*tp().thetaR(matched); + + // Calculated the component of the velocity parallel to the wall + vectorField Uwall + ( + U_.boundaryField()[patchi].patchInternalField() + - U_.boundaryField()[patchi] + ); + Uwall -= (AfHatPatch & Uwall)*AfHatPatch; + + // Find the direction of the interface parallel to the wall + vectorField nWall + ( + nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch + ); + + // Normalise nWall + nWall /= (mag(nWall) + SMALL); + + // Calculate Uwall resolved normal to the interface parallel to + // the interface + scalarField uwall(nWall & Uwall); + + theta += (thetaA - thetaR)*tanh(uwall/uTheta); + } + + + // Reset nHatPatch to correspond to the contact angle + + scalarField a12(nHatPatch & AfHatPatch); + + scalarField b1(cos(theta)); + + scalarField b2(nHatPatch.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); + + nHatPatch = a*AfHatPatch + b*nHatPatch; + + nHatPatch /= (mag(nHatPatch) + deltaN_.value()); + } + } +} + + +Foam::tmp Foam::multiphaseMixtureThermo::K +( + const phaseModel& alpha1, + const phaseModel& alpha2 +) const +{ + tmp tnHatfv = nHatfv(alpha1, alpha2); + + correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField()); + + // Simple expression for curvature + return -fvc::div(tnHatfv & mesh_.Sf()); +} + + +Foam::tmp +Foam::multiphaseMixtureThermo::nearInterface() const +{ + tmp tnearInt + ( + new volScalarField + ( + IOobject + ( + "nearInterface", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("nearInterface", dimless, 0.0) + ) + ); + + forAllConstIter(PtrDictionary, phases_, phase) + { + tnearInt() = max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase())); + } + + return tnearInt; +} + + +void Foam::multiphaseMixtureThermo::solveAlphas +( + const scalar cAlpha +) +{ + static label nSolves=-1; + nSolves++; + + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phic(mag(phi_/mesh_.magSf())); + phic = min(cAlpha*phic, max(phic)); + + PtrList phiAlphaCorrs(phases_.size()); + int phasei = 0; + + forAllIter(PtrDictionary, phases_, phase) + { + phaseModel& alpha = phase(); + + phiAlphaCorrs.set + ( + phasei, + new surfaceScalarField + ( + fvc::flux + ( + phi_, + alpha, + alphaScheme + ) + ) + ); + + surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei]; + + forAllIter(PtrDictionary, phases_, phase2) + { + phaseModel& alpha2 = phase2(); + + if (&alpha2 == &alpha) continue; + + surfaceScalarField phir(phic*nHatf(alpha, alpha2)); + + phiAlphaCorr += fvc::flux + ( + -fvc::flux(-phir, alpha2, alpharScheme), + alpha, + alpharScheme + ); + } + + MULES::limit + ( + geometricOneField(), + alpha, + phi_, + phiAlphaCorr, + zeroField(), + zeroField(), + 1, + 0, + 3, + true + ); + + phasei++; + } + + MULES::limitSum(phiAlphaCorrs); + + rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0); + + volScalarField sumAlpha + ( + IOobject + ( + "sumAlpha", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("sumAlpha", dimless, 0) + ); + + + volScalarField divU(fvc::div(fvc::absolute(phi_, U_))); + + + phasei = 0; + + forAllIter(PtrDictionary, phases_, phase) + { + phaseModel& alpha = phase(); + + surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; + phiAlpha += upwind(mesh_, phi_).flux(alpha); + + volScalarField::DimensionedInternalField Sp + ( + IOobject + ( + "Sp", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0) + ); + + volScalarField::DimensionedInternalField Su + ( + IOobject + ( + "Su", + mesh_.time().timeName(), + mesh_ + ), + // Divergence term is handled explicitly to be + // consistent with the explicit transport solution + divU*min(alpha, scalar(1)) + ); + + { + const scalarField& dgdt = alpha.dgdt(); + + forAll(dgdt, celli) + { + if (dgdt[celli] < 0.0 && alpha[celli] > 0.0) + { + Sp[celli] += dgdt[celli]*alpha[celli]; + Su[celli] -= dgdt[celli]*alpha[celli]; + } + else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0) + { + Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]); + } + } + } + + forAllConstIter(PtrDictionary, phases_, phase2) + { + const phaseModel& alpha2 = phase2(); + + if (&alpha2 == &alpha) continue; + + const scalarField& dgdt2 = alpha2.dgdt(); + + forAll(dgdt2, celli) + { + if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0) + { + Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]); + Su[celli] += dgdt2[celli]*alpha[celli]; + } + else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0) + { + Sp[celli] += dgdt2[celli]*alpha2[celli]; + } + } + } + + MULES::explicitSolve + ( + geometricOneField(), + alpha, + phiAlpha, + Sp, + Su + ); + + rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha; + + Info<< alpha.name() << " volume fraction, min, max = " + << alpha.weightedAverage(mesh_.V()).value() + << ' ' << min(alpha).value() + << ' ' << max(alpha).value() + << endl; + + sumAlpha += alpha; + + phasei++; + } + + Info<< "Phase-sum volume fraction, min, max = " + << sumAlpha.weightedAverage(mesh_.V()).value() + << ' ' << min(sumAlpha).value() + << ' ' << max(sumAlpha).value() + << endl; + + calcAlphas(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H new file mode 100644 index 0000000000..546048894f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H @@ -0,0 +1,434 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::multiphaseMixtureThermo + +Description + +SourceFiles + multiphaseMixtureThermo.C + +\*---------------------------------------------------------------------------*/ + +#ifndef multiphaseMixtureThermo_H +#define multiphaseMixtureThermo_H + +#include "phaseModel.H" +#include "PtrDictionary.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "rhoThermo.H" +#include "psiThermo.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class multiphaseMixtureThermo Declaration +\*---------------------------------------------------------------------------*/ + +class multiphaseMixtureThermo +: + public psiThermo +{ +public: + + class interfacePair + : + public Pair + { + public: + + class hash + : + public Hash + { + public: + + hash() + {} + + label operator()(const interfacePair& key) const + { + return word::hash()(key.first()) + word::hash()(key.second()); + } + }; + + + // Constructors + + interfacePair() + {} + + interfacePair(const word& alpha1Name, const word& alpha2Name) + : + Pair(alpha1Name, alpha2Name) + {} + + interfacePair(const phaseModel& alpha1, const phaseModel& alpha2) + : + Pair(alpha1.name(), alpha2.name()) + {} + + + // Friend Operators + + friend bool operator== + ( + const interfacePair& a, + const interfacePair& b + ) + { + return + ( + ((a.first() == b.first()) && (a.second() == b.second())) + || ((a.first() == b.second()) && (a.second() == b.first())) + ); + } + + friend bool operator!= + ( + const interfacePair& a, + const interfacePair& b + ) + { + return (!(a == b)); + } + }; + + +private: + + // Private data + + //- Dictionary of phases + PtrDictionary phases_; + + const fvMesh& mesh_; + const volVectorField& U_; + const surfaceScalarField& phi_; + + surfaceScalarField rhoPhi_; + + volScalarField alphas_; + + typedef HashTable + sigmaTable; + + sigmaTable sigmas_; + dimensionSet dimSigma_; + + //- Stabilisation for normalisation of the interface normal + const dimensionedScalar deltaN_; + + //- Conversion factor for degrees into radians + static const scalar convertToRad; + + + // Private member functions + + void calcAlphas(); + + void solveAlphas(const scalar cAlpha); + + tmp nHatfv + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + tmp nHatf + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + void correctContactAngle + ( + const phaseModel& alpha1, + const phaseModel& alpha2, + surfaceVectorField::GeometricBoundaryField& nHatb + ) const; + + tmp K + ( + const phaseModel& alpha1, + const phaseModel& alpha2 + ) const; + + +public: + + //- Runtime type information + TypeName("multiphaseMixtureThermo"); + + + // Constructors + + //- Construct from components + multiphaseMixtureThermo + ( + const volVectorField& U, + const surfaceScalarField& phi + ); + + + //- Destructor + virtual ~multiphaseMixtureThermo() + {} + + + // Member Functions + + //- Return the phases + const PtrDictionary& phases() const + { + return phases_; + } + + //- Return non-const access to the phases + PtrDictionary& phases() + { + return phases_; + } + + //- Return the velocity + const volVectorField& U() const + { + return U_; + } + + //- Return the volumetric flux + const surfaceScalarField& phi() const + { + return phi_; + } + + const surfaceScalarField& rhoPhi() const + { + return rhoPhi_; + } + + //- Update properties + virtual void correct(); + + //- Update densities for given pressure change + void correctRho(const volScalarField& dp); + + //- Return true if the equation of state is incompressible + // i.e. rho != f(p) + virtual bool incompressible() const; + + //- Return true if the equation of state is isochoric + // i.e. rho = const + virtual bool isochoric() const; + + + // Access to thermodynamic state variables + + //- Enthalpy/Internal energy [J/kg] + // Non-const access allowed for transport equations + virtual volScalarField& he() + { + notImplemented("multiphaseMixtureThermo::he()"); + return phases_[0]->thermo().he(); + } + + //- Enthalpy/Internal energy [J/kg] + virtual const volScalarField& he() const + { + notImplemented("multiphaseMixtureThermo::he() const"); + return phases_[0]->thermo().he(); + } + + //- Enthalpy/Internal energy + // for given pressure and temperature [J/kg] + virtual tmp he + ( + const volScalarField& p, + const volScalarField& T + ) const; + + //- Enthalpy/Internal energy for cell-set [J/kg] + virtual tmp he + ( + const scalarField& p, + const scalarField& T, + const labelList& cells + ) const; + + //- Enthalpy/Internal energy for patch [J/kg] + virtual tmp he + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Chemical enthalpy [J/kg] + virtual tmp hc() const; + + //- Temperature from enthalpy/internal energy for cell-set + virtual tmp THE + ( + const scalarField& h, + const scalarField& p, + const scalarField& T0, // starting temperature + const labelList& cells + ) const; + + //- Temperature from enthalpy/internal energy for patch + virtual tmp THE + ( + const scalarField& h, + const scalarField& p, + const scalarField& T0, // starting temperature + const label patchi + ) const; + + + // Fields derived from thermodynamic state variables + + //- Density [kg/m^3] + virtual tmp rho() const; + + //- Heat capacity at constant pressure [J/kg/K] + virtual tmp Cp() const; + + //- Heat capacity at constant pressure for patch [J/kg/K] + virtual tmp Cp + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant volume [J/kg/K] + virtual tmp Cv() const; + + //- Heat capacity at constant volume for patch [J/kg/K] + virtual tmp Cv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- gamma = Cp/Cv [] + virtual tmp gamma() const; + + //- gamma = Cp/Cv for patch [] + virtual tmp gamma + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant pressure/volume [J/kg/K] + virtual tmp Cpv() const; + + //- Heat capacity at constant pressure/volume for patch [J/kg/K] + virtual tmp Cpv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity ratio [] + virtual tmp CpByCpv() const; + + //- Heat capacity ratio for patch [] + virtual tmp CpByCpv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + + // Fields derived from transport state variables + + //- Thermal diffusivity for temperature of mixture [J/m/s/K] + virtual tmp kappa() const; + + //- Thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp kappa + ( + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp kappaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp kappaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp alphaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp alphaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + + //- Return the phase-averaged reciprocal Cv + tmp rCv() const; + + tmp surfaceTensionForce() const; + + //- Indicator of the proximity of the interface + // Field values are 1 near and 0 away for the interface. + tmp nearInterface() const; + + //- Solve for the mixture phase-fractions + void solve(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C new file mode 100644 index 0000000000..1559f25a48 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "phaseModel.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseModel::phaseModel +( + const word& phaseName, + const volScalarField& p, + const volScalarField& T +) +: + volScalarField + ( + IOobject + ( + IOobject::groupName("alpha", phaseName), + p.mesh().time().timeName(), + p.mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + p.mesh() + ), + name_(phaseName), + p_(p), + T_(T), + thermo_(NULL), + dgdt_ + ( + IOobject + ( + IOobject::groupName("dgdt", phaseName), + p.mesh().time().timeName(), + p.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + p.mesh(), + dimensionedScalar("0", dimless/dimTime, 0) + ) +{ + { + volScalarField Tp(IOobject::groupName("T", phaseName), T); + Tp.write(); + } + + thermo_ = rhoThermo::New(p.mesh(), phaseName); + thermo_->validate(phaseName, "e"); + + correct(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::autoPtr Foam::phaseModel::clone() const +{ + notImplemented("phaseModel::clone() const"); + return autoPtr(NULL); +} + + +void Foam::phaseModel::correct() +{ + thermo_->he() = thermo_->he(p_, T_); + thermo_->correct(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H new file mode 100644 index 0000000000..66d0ac8d63 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H @@ -0,0 +1,156 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::phaseModel + +Description + Single incompressible phase derived from the phase-fraction. + Used as part of the multiPhaseMixture for interface-capturing multi-phase + simulations. + +SourceFiles + phaseModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef phaseModel_H +#define phaseModel_H + +#include "rhoThermo.H" +#include "volFields.H" +#include "dictionaryEntry.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class phaseModel Declaration +\*---------------------------------------------------------------------------*/ + +class phaseModel +: + public volScalarField +{ + // Private data + + word name_; + const volScalarField& p_; + const volScalarField& T_; + autoPtr thermo_; + volScalarField dgdt_; + + +public: + + // Constructors + + //- Construct from components + phaseModel + ( + const word& phaseName, + const volScalarField& p, + const volScalarField& T + ); + + //- Return clone + autoPtr clone() const; + + //- Return a pointer to a new phaseModel created on freestore + // from Istream + class iNew + { + const volScalarField& p_; + const volScalarField& T_; + + public: + + iNew + ( + const volScalarField& p, + const volScalarField& T + ) + : + p_(p), + T_(T) + {} + + autoPtr operator()(Istream& is) const + { + return autoPtr(new phaseModel(is, p_, T_)); + } + }; + + + // Member Functions + + const word& name() const + { + return name_; + } + + const word& keyword() const + { + return name(); + } + + //- Return const-access to phase rhoThermo + const rhoThermo& thermo() const + { + return thermo_(); + } + + //- Return access to phase rhoThermo + rhoThermo& thermo() + { + return thermo_(); + } + + //- Return const-access to phase divergence + const volScalarField& dgdt() const + { + return dgdt_; + } + + //- Return access to phase divergence + volScalarField& dgdt() + { + return dgdt_; + } + + void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H new file mode 100644 index 0000000000..e42a54aabd --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H @@ -0,0 +1,142 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + + surfaceScalarField phig + ( + ( + multiphaseProperties.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + // Update the fixedFluxPressure BCs to ensure flux consistency + setSnGrad + ( + p_rgh.boundaryField(), + ( + phiHbyA.boundaryField() + - (mesh.Sf().boundaryField() & U.boundaryField()) + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + + PtrList p_rghEqnComps(multiphaseProperties.phases().size()); + + label phasei = 0; + forAllConstIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + const rhoThermo& thermo = phase().thermo(); + const volScalarField& rho = thermo.rho()(); + + p_rghEqnComps.set + ( + phasei, + ( + fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh)) + + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho) + ).ptr() + ); + + phasei++; + } + + // Cache p_rgh prior to solve for density update + volScalarField p_rgh_0(p_rgh); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqnIncomp + ( + fvc::div(phiHbyA) + - fvm::laplacian(rAUf, p_rgh) + ); + + tmp p_rghEqnComp; + + phasei = 0; + forAllConstIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + tmp hmm + ( + (max(phase(), scalar(0))/phase().thermo().rho()) + *p_rghEqnComps[phasei] + ); + + if (phasei == 0) + { + p_rghEqnComp = hmm; + } + else + { + p_rghEqnComp() += hmm; + } + + phasei++; + } + + solve + ( + p_rghEqnComp + + p_rghEqnIncomp, + mesh.solver(p_rgh.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + // p = max(p_rgh + multiphaseProperties.rho()*gh, pMin); + // p_rgh = p - multiphaseProperties.rho()*gh; + + phasei = 0; + forAllIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + phase().dgdt() = + pos(phase()) + *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho(); + } + + phi = phiHbyA + p_rghEqnIncomp.flux(); + + U = HbyA + + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf); + U.correctBoundaryConditions(); + } + } + + p = max(p_rgh + multiphaseProperties.rho()*gh, pMin); + + // Update densities from change in p_rgh + multiphaseProperties.correctRho(p_rgh - p_rgh_0); + rho = multiphaseProperties.rho(); + + K = 0.5*magSqr(U); + + Info<< "max(U) " << max(mag(U)).value() << endl; + Info<< "min(p_rgh) " << min(p_rgh).value() << endl; +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index f347de82a0..995d07d7e4 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -126,11 +126,25 @@ void Foam::multiphaseSystem::solveAlphas() } // Ensure that the flux at inflow BCs is preserved - phiAlphaCorr.boundaryField() = min - ( - phase1.phi().boundaryField()*alpha1.boundaryField(), - phiAlphaCorr.boundaryField() - ); + forAll(phiAlphaCorr.boundaryField(), patchi) + { + fvsPatchScalarField& phiAlphaCorrp = + phiAlphaCorr.boundaryField()[patchi]; + + if (!phiAlphaCorrp.coupled()) + { + const scalarField& phi1p = phase1.phi().boundaryField()[patchi]; + const scalarField& alpha1p = alpha1.boundaryField()[patchi]; + + forAll(phiAlphaCorrp, facei) + { + if (phi1p[facei] < 0) + { + phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei]; + } + } + } + } MULES::limit ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H index 3127450f00..f53973f0cc 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H @@ -23,7 +23,7 @@ )/rho1 //***HGW- fvm::laplacian(alpha1*turbulence1->alphaEff(), he1) - - fvm::laplacian(alpha1*turbulence1->nuEff(), he1) + - fvm::laplacian(alpha1*phase1.turbulence().nuEff(), he1) == heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1 + heatTransferCoeff*he1/Cpv1/rho1 @@ -46,7 +46,7 @@ )/rho2 //***HGW- fvm::laplacian(alpha2*turbulence2->alphaEff(), he2) - - fvm::laplacian(alpha2*turbulence2->nuEff(), he2) + - fvm::laplacian(alpha2*phase2.turbulence().nuEff(), he2) == heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2 + heatTransferCoeff*he2/Cpv2/rho2 diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H index dea48be335..fad73cfeb8 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H @@ -26,7 +26,7 @@ volScalarField dragCoeff(fluid.dragCoeff()); - fvm::Sp(fvc::div(phi1), U1) ) - + turbulence1->divDevReff(U1) + + phase1.turbulence().divDevReff(U1) == - fvm::Sp(dragCoeff/rho1, U1) - alpha1*alpha2/rho1*(liftForce - fluid.Cvm()*rho2*DDtU2) @@ -50,7 +50,7 @@ volScalarField dragCoeff(fluid.dragCoeff()); + fvm::div(phi2, U2) - fvm::Sp(fvc::div(phi2), U2) ) - + turbulence2->divDevReff(U2) + + phase2.turbulence().divDevReff(U2) == - fvm::Sp(dragCoeff/rho2, U2) + alpha1*alpha2/rho2*(liftForce + fluid.Cvm()*rho2*DDtU1) diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H deleted file mode 100644 index be420a36e1..0000000000 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H +++ /dev/null @@ -1,165 +0,0 @@ -{ - word alphaScheme("div(phi," + alpha1.name() + ')'); - word alpharScheme("div(phir," + alpha1.name() + ')'); - - alpha1.correctBoundaryConditions(); - - surfaceScalarField phic("phic", phi); - surfaceScalarField phir("phir", phi1 - phi2); - - surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); - - tmp pPrimeByA; - - if (implicitPhasePressure) - { - pPrimeByA = - fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime()) - + fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime()); - - surfaceScalarField phiP - ( - pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh.magSf() - ); - - phic += alpha1f*phiP; - phir += phiP; - } - - for (int acorr=0; acorr 0.0 && alpha1[celli] > 0.0) - { - Sp[celli] -= dgdt[celli]*alpha1[celli]; - Su[celli] += dgdt[celli]*alpha1[celli]; - } - else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0) - { - Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]); - } - } - - dimensionedScalar totalDeltaT = runTime.deltaT(); - if (nAlphaSubCycles > 1) - { - alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0); - } - - for - ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); - !(++alphaSubCycle).end(); - ) - { - surfaceScalarField alphaPhic1 - ( - fvc::flux - ( - phic, - alpha1, - alphaScheme - ) - + fvc::flux - ( - -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), - alpha1, - alpharScheme - ) - ); - - // Ensure that the flux at inflow BCs is preserved - forAll(alphaPhic1.boundaryField(), patchi) - { - fvsPatchScalarField& alphaPhic1p = - alphaPhic1.boundaryField()[patchi]; - - if (!alphaPhic1p.coupled()) - { - const scalarField& phi1p = phi1.boundaryField()[patchi]; - const scalarField& alpha1p = alpha1.boundaryField()[patchi]; - - forAll(alphaPhic1p, facei) - { - if (phi1p[facei] < 0) - { - alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei]; - } - } - } - } - - MULES::explicitSolve - ( - geometricOneField(), - alpha1, - phi, - alphaPhic1, - Sp, - Su, - 1, - 0 - ); - - if (nAlphaSubCycles > 1) - { - alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1; - } - else - { - alphaPhi1 = alphaPhic1; - } - } - - if (implicitPhasePressure) - { - fvScalarMatrix alpha1Eqn - ( - fvm::ddt(alpha1) - fvc::ddt(alpha1) - - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded") - ); - - alpha1Eqn.relax(); - alpha1Eqn.solve(); - - alphaPhi1 += alpha1Eqn.flux(); - } - - alphaPhi2 = phi - alphaPhi1; - alpha2 = scalar(1) - alpha1; - - Info<< "Dispersed phase volume fraction = " - << alpha1.weightedAverage(mesh.V()).value() - << " Min(alpha1) = " << min(alpha1).value() - << " Max(alpha1) = " << max(alpha1).value() - << endl; - } -} - -rho = fluid.rho(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index f279445cf6..a382b769b5 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H @@ -10,19 +10,13 @@ volVectorField& U1 = phase1.U(); surfaceScalarField& phi1 = phase1.phi(); - surfaceScalarField alphaPhi1 - ( - IOobject::groupName("alphaPhi", phase1.name()), - fvc::interpolate(alpha1)*phi1 - ); + surfaceScalarField& alphaPhi1 = phase1.phiAlpha(); volVectorField& U2 = phase2.U(); surfaceScalarField& phi2 = phase2.phi(); - surfaceScalarField alphaPhi2 - ( - IOobject::groupName("alphaPhi", phase2.name()), - fvc::interpolate(alpha2)*phi2 - ); + surfaceScalarField& alphaPhi2 = phase2.phiAlpha(); + + surfaceScalarField& phi = fluid.phi(); dimensionedScalar pMin ( @@ -55,19 +49,6 @@ fluid.U() ); - surfaceScalarField phi - ( - IOobject - ( - "phi", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - fluid.phi() - ); - volScalarField rho ( IOobject @@ -97,10 +78,6 @@ - fvc::div(phi2)*U2 ); - - Info<< "Calculating field g.h\n" << endl; - volScalarField gh("gh", g & mesh.C()); - volScalarField rAU1 ( IOobject @@ -133,13 +110,6 @@ scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue); - - volScalarField dgdt - ( - pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001)) - ); - - Info<< "Creating field dpdt\n" << endl; volScalarField dpdt ( @@ -157,29 +127,3 @@ Info<< "Creating field kinetic energy K\n" << endl; volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1)); volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2)); - - autoPtr > - turbulence1 - ( - PhaseIncompressibleTurbulenceModel::New - ( - alpha1, - U1, - alphaPhi1, - phi1, - phase1 - ) - ); - - autoPtr > - turbulence2 - ( - PhaseIncompressibleTurbulenceModel::New - ( - alpha2, - U2, - alphaPhi2, - phi2, - phase2 - ) - ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index 9e99c5cc61..ea3f6bdb33 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -31,17 +31,19 @@ surfaceScalarField phiP1 ( "phiP1", - fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime()) + fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime()) *fvc::snGrad(alpha1)*mesh.magSf() ); + phiP1.boundaryField() == 0; // Phase-2 pressure flux (e.g. due to particle-particle pressure) surfaceScalarField phiP2 ( "phiP2", - fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime()) + fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime()) *fvc::snGrad(alpha2)*mesh.magSf() ); + phiP2.boundaryField() == 0; surfaceScalarField phiHbyA1 ( @@ -126,9 +128,11 @@ ); pEqnComp1 = - fvc::ddt(rho1) - + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1) - + correction + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1) + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1/rho1)*correction ( psi1*fvm::ddt(p) + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p) @@ -137,9 +141,11 @@ pEqnComp1().relax(); pEqnComp2 = - fvc::ddt(rho2) - + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2) - + correction + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2) + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + (alpha2/rho2)*correction ( psi2*fvm::ddt(p) + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p) @@ -150,12 +156,18 @@ else { pEqnComp1 = - fvc::ddt(rho1) + psi1*correction(fvm::ddt(p)) - + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1); + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1) + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1*psi1/rho1)*correction(fvm::ddt(p)); pEqnComp2 = - fvc::ddt(rho2) + psi2*correction(fvm::ddt(p)) - + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2); + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2) + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + (alpha2*psi2/rho2)*correction(fvm::ddt(p)); } // Cache p prior to solve for density update @@ -171,11 +183,7 @@ solve ( - ( - (alpha1/rho1)*pEqnComp1() - + (alpha2/rho2)*pEqnComp2() - ) - + pEqnIncomp, + pEqnComp1() + pEqnComp2() + pEqnIncomp, mesh.solver(p.select(pimple.finalInnerIter())) ); @@ -199,10 +207,10 @@ phi = alpha1f*phi1 + alpha2f*phi2; - dgdt = + fluid.dgdt() = ( - pos(alpha2)*(pEqnComp2 & p)/rho2 - - pos(alpha1)*(pEqnComp1 & p)/rho1 + pos(alpha2)*(pEqnComp2 & p)/max(alpha2, scalar(1e-3)) + - pos(alpha1)*(pEqnComp1 & p)/max(alpha1, scalar(1e-3)) ); p.relax(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H deleted file mode 100644 index 29353a8fa1..0000000000 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H +++ /dev/null @@ -1,7 +0,0 @@ - #include "readTimeControls.H" - #include "alphaControls.H" - - Switch implicitPhasePressure - ( - alphaControls.lookupOrDefault("implicitPhasePressure", false) - ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index 1e7cbafea7..69682a9a3d 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C @@ -31,13 +31,10 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" -#include "subCycle.H" -#include "rhoThermo.H" #include "twoPhaseSystem.H" +#include "PhaseIncompressibleTurbulenceModel.H" #include "dragModel.H" #include "heatTransferModel.H" -#include "PhaseIncompressibleTurbulenceModel.H" #include "pimpleControl.H" #include "IOMRFZoneList.H" #include "fixedFluxPressureFvPatchScalarField.H" @@ -66,7 +63,7 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readTwoPhaseEulerFoamControls.H" + #include "readTimeControls.H" #include "CourantNos.H" #include "setDeltaT.H" @@ -76,7 +73,10 @@ int main(int argc, char *argv[]) // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { - #include "alphaEqn.H" + fluid.solve(); + rho = fluid.rho(); + fluid.correct(); + #include "EEqns.H" #include "UEqns.H" @@ -90,8 +90,7 @@ int main(int argc, char *argv[]) if (pimple.turbCorr()) { - turbulence1->correct(); - turbulence2->correct(); + fluid.correctTurbulence(); } } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files index fb49c3cef7..0de9e8939c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files @@ -4,6 +4,12 @@ diameterModels/diameterModel/newDiameterModel.C diameterModels/constantDiameter/constantDiameter.C diameterModels/isothermalDiameter/isothermalDiameter.C +diameterModels/IATE/IATE.C +diameterModels/IATE/IATEsources/IATEsource/IATEsource.C +diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C +diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C +diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C + twoPhaseSystem.C LIB = $(FOAM_LIBBIN)/libcompressibleTwoPhaseSystem diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options index ab3c396f57..eec60d23e7 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options @@ -3,7 +3,10 @@ EXE_INC = \ -I../interfacialModels/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/transportModels/incompressible/lnInclude + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude LIB_LIBS = \ -lincompressibleTransportModels \ diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C new file mode 100644 index 0000000000..ef4bc12d0c --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C @@ -0,0 +1,188 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "IATE.H" +#include "IATEsource.H" +#include "twoPhaseSystem.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvmSup.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcAverage.H" +#include "mathematicalConstants.H" +#include "fundamentalConstants.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ + defineTypeNameAndDebug(IATE, 0); + + addToRunTimeSelectionTable + ( + diameterModel, + IATE, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::diameterModels::IATE::IATE +( + const dictionary& diameterProperties, + const phaseModel& phase +) +: + diameterModel(diameterProperties, phase), + kappai_ + ( + IOobject + ( + IOobject::groupName("kappai", phase.name()), + phase_.U().time().timeName(), + phase_.U().mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + phase_.U().mesh() + ), + dMax_("dMax", dimLength, diameterProperties_.lookup("dMax")), + dMin_("dMin", dimLength, diameterProperties_.lookup("dMin")), + d_ + ( + IOobject + ( + IOobject::groupName("d", phase.name()), + phase_.U().time().timeName(), + phase_.U().mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + dsm() + ), + sources_ + ( + diameterProperties_.lookup("sources"), + IATEsource::iNew(*this) + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::diameterModels::IATE::~IATE() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::diameterModels::IATE::dsm() const +{ + return max(6/max(kappai_, 6/dMax_), dMin_); +} + +// Placeholder for the nucleation/condensation model +// Foam::tmp Foam::diameterModels::IATE::Rph() const +// { +// const volScalarField& T = phase_.thermo().T(); +// const volScalarField& p = phase_.thermo().p(); +// +// scalar A, B, C, sigma, vm, Rph; +// +// volScalarField ps(1e5*pow(10, A - B/(T + C))); +// volScalarField Dbc +// ( +// 4*sigma*vm/(constant::physicoChemical::k*T*log(p/ps)) +// ); +// +// return constant::mathematical::pi*sqr(Dbc)*Rph; +// } + +void Foam::diameterModels::IATE::correct() +{ + // Initialise the accumulated source term to the dilatation effect + volScalarField R + ( + ( + (2.0/3.0) + /max + ( + fvc::average(phase_ + phase_.oldTime()), + phase_.fluid().residualPhaseFraction() + ) + ) + *(fvc::ddt(phase_) + fvc::div(phase_.phiAlpha())) + ); + + // Accumulate the run-time selectable sources + forAll(sources_, j) + { + R -= sources_[j].R(); + } + + // Construct the interfacial curvature equation + fvScalarMatrix kappaiEqn + ( + fvm::ddt(kappai_) + fvm::div(phase_.phi(), kappai_) + - fvm::Sp(fvc::div(phase_.phi()), kappai_) + == + - fvm::SuSp(R, kappai_) + //+ Rph() // Omit the nucleation/condensation term + ); + + kappaiEqn.relax(); + kappaiEqn.solve(); + + // Update the Sauter-mean diameter + d_ = dsm(); +} + + +bool Foam::diameterModels::IATE::read(const dictionary& phaseProperties) +{ + diameterModel::read(phaseProperties); + + diameterProperties_.lookup("dMax") >> dMax_; + diameterProperties_.lookup("dMin") >> dMin_; + + // Re-create all the sources updating number, type and coefficients + PtrList + ( + diameterProperties_.lookup("sources"), + IATEsource::iNew(*this) + ).transfer(sources_); + + return true; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H new file mode 100644 index 0000000000..f599fb9959 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATE + +Description + IATE (Interfacial Area Transport Equation) bubble diameter model. + + Solves for the interfacial curvature per unit volume of the phase rather + than interfacial area per unit volume to avoid stability issues relating to + the consistency requirements between the phase fraction and interfacial area + per unit volume. In every other respect this model is as presented in the + paper: + + \verbatim + "Development of Interfacial Area Transport Equation" + M. Ishii, + S. Kim, + J Kelly, + Nuclear Engineering and Technology, Vol.37 No.6 December 2005 + \endverbatim + +SourceFiles + IATE.C + +\*---------------------------------------------------------------------------*/ + +#ifndef IATE_H +#define IATE_H + +#include "diameterModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ + +// Forward declaration of classes +class IATEsource; + +/*---------------------------------------------------------------------------*\ + Class IATE Declaration +\*---------------------------------------------------------------------------*/ + +class IATE +: + public diameterModel +{ + // Private data + + //- Interfacial curvature (alpha*interfacial area) + volScalarField kappai_; + + //- Maximum diameter used for stabilisation in the limit kappai->0 + dimensionedScalar dMax_; + + //- Minimum diameter used for stabilisation in the limit kappai->inf + dimensionedScalar dMin_; + + //- The Sauter-mean diameter of the phase + volScalarField d_; + + //- IATE sources + PtrList sources_; + + + // Private member functions + + tmp dsm() const; + + +public: + + friend class IATEsource; + + //- Runtime type information + TypeName("IATE"); + + + // Constructors + + //- Construct from components + IATE + ( + const dictionary& diameterProperties, + const phaseModel& phase + ); + + + //- Destructor + virtual ~IATE(); + + + // Member Functions + + //- Return the interfacial curvature + const volScalarField& kappai() const + { + return kappai_; + } + + //- Return the interfacial area + tmp a() const + { + return phase_*kappai_; + } + + //- Return the Sauter-mean diameter + virtual tmp d() const + { + return d_; + } + + //- Correct the diameter field + virtual void correct(); + + //- Read phaseProperties dictionary + virtual bool read(const dictionary& phaseProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C new file mode 100644 index 0000000000..4d83431331 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "IATEsource.H" +#include "twoPhaseSystem.H" +#include "fvMatrix.H" +#include "PhaseIncompressibleTurbulenceModel.H" +#include "uniformDimensionedFields.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ + defineTypeNameAndDebug(IATEsource, 0); + defineRunTimeSelectionTable(IATEsource, dictionary); +} +} + + +// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * // + +Foam::autoPtr +Foam::diameterModels::IATEsource::New +( + const word& type, + const IATE& iate, + const dictionary& dict +) +{ + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(type); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "IATEsource::New" + "(const word& type, const IATE&, const dictionary&)" + ) << "Unknown IATE source type " + << type << nl << nl + << "Valid IATE source types : " << endl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(iate, dict)); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::diameterModels::IATEsource::Ur() const +{ + const uniformDimensionedVectorField& g = + phase().U().db().lookupObject("g"); + + return + sqrt(2.0) + *pow025 + ( + fluid().sigma()*mag(g) + *(otherPhase().rho() - phase().rho()) + /sqr(otherPhase().rho()) + ) + *pow(max(1 - phase(), scalar(0)), 1.75); +} + +Foam::tmp Foam::diameterModels::IATEsource::Ut() const +{ + return sqrt(2*otherPhase().turbulence().k()); +} + +Foam::tmp Foam::diameterModels::IATEsource::Re() const +{ + return max(Ur()*phase().d()/otherPhase().nu(), scalar(1.0e-3)); +} + +Foam::tmp Foam::diameterModels::IATEsource::CD() const +{ + const volScalarField Eo(this->Eo()); + const volScalarField Re(this->Re()); + + return + max + ( + min + ( + (16/Re)*(1 + 0.15*pow(Re, 0.687)), + 48/Re + ), + 8*Eo/(3*(Eo + 4)) + ); +} + +Foam::tmp Foam::diameterModels::IATEsource::Mo() const +{ + const uniformDimensionedVectorField& g = + phase().U().db().lookupObject("g"); + + return + mag(g)*pow4(otherPhase().nu())*sqr(otherPhase().rho()) + *(otherPhase().rho() - phase().rho()) + /pow3(fluid().sigma()); +} + +Foam::tmp Foam::diameterModels::IATEsource::Eo() const +{ + const uniformDimensionedVectorField& g = + phase().U().db().lookupObject("g"); + + return + mag(g)*sqr(phase().d()) + *(otherPhase().rho() - phase().rho()) + /fluid().sigma(); +} + +Foam::tmp Foam::diameterModels::IATEsource::We() const +{ + return otherPhase().rho()*sqr(Ur())*phase().d()/fluid().sigma(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H new file mode 100644 index 0000000000..f21f4f3876 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H @@ -0,0 +1,192 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATEsource + +Description + IATE (Interfacial Area Transport Equation) bubble diameter model + run-time selectable sources. + +SourceFiles + IATEsource.C + +\*---------------------------------------------------------------------------*/ + +#ifndef IATEsource_H +#define IATEsource_H + +#include "IATE.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ + +/*---------------------------------------------------------------------------*\ + Class IATEsource Declaration +\*---------------------------------------------------------------------------*/ + +class IATEsource +{ + +protected: + + // Protected data + + //- Reference to the IATE this source applies to + const IATE& iate_; + + +public: + + //- Runtime type information + TypeName("IATEsource"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + IATEsource, + dictionary, + ( + const IATE& iate, + const dictionary& dict + ), + (iate, dict) + ); + + + //- Class used for the read-construction of + // PtrLists of IATE sources + class iNew + { + const IATE& iate_; + + public: + + iNew(const IATE& iate) + : + iate_(iate) + {} + + autoPtr operator()(Istream& is) const + { + word type(is); + dictionary dict(is); + return IATEsource::New(type, iate_, dict); + } + }; + + + // Constructors + + IATEsource(const IATE& iate) + : + iate_(iate) + {} + + autoPtr clone() const + { + notImplemented("autoPtr clone() const"); + return autoPtr(NULL); + } + + + // Selectors + + static autoPtr New + ( + const word& type, + const IATE& iate, + const dictionary& dict + ); + + + //- Destructor + virtual ~IATEsource() + {} + + + // Member Functions + + const phaseModel& phase() const + { + return iate_.phase(); + } + + const twoPhaseSystem& fluid() const + { + return iate_.phase().fluid(); + } + + const phaseModel& otherPhase() const + { + return phase().otherPhase(); + } + + scalar phi() const + { + return 1.0/(36*constant::mathematical::pi); + } + + //- Return the bubble relative velocity + tmp Ur() const; + + //- Return the bubble turbulent velocity + tmp Ut() const; + + //- Return the bubble Reynolds number + tmp Re() const; + + //- Return the bubble drag coefficient + tmp CD() const; + + //- Return the bubble Morton number + tmp Mo() const; + + //- Return the bubble Eotvos number + tmp Eo() const; + + //- Return the bubble Webber number + tmp We() const; + + virtual tmp R() const = 0; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C new file mode 100644 index 0000000000..d442373b0c --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C @@ -0,0 +1,66 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "dummy.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + defineTypeNameAndDebug(dummy, 0); + addToRunTimeSelectionTable(IATEsource, dummy, word); +} +} +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::diameterModels::IATEsources::dummy::R() const +{ + return tmp + ( + new volScalarField + ( + IOobject + ( + "R", + iate_.phase().U().time().timeName(), + iate_.phase().mesh() + ), + iate_.phase().U().mesh(), + dimensionedScalar("R", dimless/dimTime, 0) + ) + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H new file mode 100644 index 0000000000..9fcdbf5325 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H @@ -0,0 +1,96 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATEsources::dummy + +Description + +SourceFiles + dummy.C + +\*---------------------------------------------------------------------------*/ + +#ifndef dummy_H +#define dummy_H + +#include "IATEsource.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + +/*---------------------------------------------------------------------------*\ + Class dummy Declaration +\*---------------------------------------------------------------------------*/ + +class dummy +: + public IATEsource +{ + +public: + + //- Runtime type information + TypeName("dummy"); + + // Constructors + + dummy + ( + const word& name, + const IATE& iate, + const dictionary& dict + ) + : + IATEsource(iate) + {} + + + //- Destructor + virtual ~dummy() + {} + + + // Member Functions + + virtual tmp R() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace IATEsources +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C new file mode 100644 index 0000000000..3ddfd70713 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C @@ -0,0 +1,109 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "randomCoalescence.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + defineTypeNameAndDebug(randomCoalescence, 0); + addToRunTimeSelectionTable(IATEsource, randomCoalescence, dictionary); +} +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::diameterModels::IATEsources::randomCoalescence:: +randomCoalescence +( + const IATE& iate, + const dictionary& dict +) +: + IATEsource(iate), + Crc_("Crc", dimless, dict.lookup("Crc")), + C_("C", dimless, dict.lookup("C")), + alphaMax_("alphaMax", dimless, dict.lookup("alphaMax")) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::diameterModels::IATEsources::randomCoalescence::R() const +{ + tmp tR + ( + new volScalarField + ( + IOobject + ( + "R", + iate_.phase().U().time().timeName(), + iate_.phase().mesh() + ), + iate_.phase().U().mesh(), + dimensionedScalar("R", dimless/dimTime, 0) + ) + ); + + volScalarField R = tR(); + + scalar Crc = Crc_.value(); + scalar C = C_.value(); + scalar alphaMax = alphaMax_.value(); + volScalarField Ut(this->Ut()); + const volScalarField& alpha = phase(); + const volScalarField& kappai = iate_.kappai(); + scalar cbrtAlphaMax = cbrt(alphaMax); + + forAll(R, celli) + { + if (alpha[celli] < alphaMax - SMALL) + { + scalar cbrtAlphaMaxMAlpha = cbrtAlphaMax - cbrt(alpha[celli]); + + R[celli] = + 12*phi()*kappai[celli]*alpha[celli] + *Crc + *Ut[celli] + *(1 - exp(-C*cbrt(alpha[celli]*alphaMax)/cbrtAlphaMaxMAlpha)) + /(cbrtAlphaMax*cbrtAlphaMaxMAlpha); + } + } + + return tR; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H new file mode 100644 index 0000000000..d8c03b47c3 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATEsources::randomCoalescence + +Description + Random coalescence IATE source as defined in paper: + + \verbatim + "Development of Interfacial Area Transport Equation" + M. Ishii, + S. Kim, + J Kelly, + Nuclear Engineering and Technology, Vol.37 No.6 December 2005 + \endverbatim + + +SourceFiles + randomCoalescence.C + +\*---------------------------------------------------------------------------*/ + +#ifndef randomCoalescence_H +#define randomCoalescence_H + +#include "IATEsource.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + +/*---------------------------------------------------------------------------*\ + Class randomCoalescence Declaration +\*---------------------------------------------------------------------------*/ + +class randomCoalescence +: + public IATEsource +{ + // Private data + + dimensionedScalar Crc_; + dimensionedScalar C_; + dimensionedScalar alphaMax_; + + +public: + + //- Runtime type information + TypeName("randomCoalescence"); + + // Constructors + + randomCoalescence + ( + const IATE& iate, + const dictionary& dict + ); + + + //- Destructor + virtual ~randomCoalescence() + {} + + + // Member Functions + + virtual tmp R() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace IATEsources +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C new file mode 100644 index 0000000000..39bd0daf89 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "turbulentBreakUp.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + defineTypeNameAndDebug(turbulentBreakUp, 0); + addToRunTimeSelectionTable(IATEsource, turbulentBreakUp, dictionary); +} +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::diameterModels::IATEsources::turbulentBreakUp:: +turbulentBreakUp +( + const IATE& iate, + const dictionary& dict +) +: + IATEsource(iate), + Cti_("Cti", dimless, dict.lookup("Cti")), + WeCr_("WeCr", dimless, dict.lookup("WeCr")) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::diameterModels::IATEsources::turbulentBreakUp::R() const +{ + tmp tR + ( + new volScalarField + ( + IOobject + ( + "R", + iate_.phase().U().time().timeName(), + iate_.phase().mesh() + ), + iate_.phase().U().mesh(), + dimensionedScalar("R", dimless/dimTime, 0) + ) + ); + + volScalarField R = tR(); + + scalar Cti = Cti_.value(); + scalar WeCr = WeCr_.value(); + volScalarField Ut(this->Ut()); + volScalarField We(this->We()); + const volScalarField& d(iate_.d()()); + + forAll(R, celli) + { + if (We[celli] > WeCr) + { + R[celli] = + (1.0/3.0) + *Cti/d[celli] + *Ut[celli] + *sqrt(1 - WeCr/We[celli]) + *exp(-WeCr/We[celli]); + } + } + + return tR; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H new file mode 100644 index 0000000000..c900b5d592 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATEsources::turbulentBreakUp + +Description + Turbulence-induced break-up IATE source as defined in paper: + + \verbatim + "Development of Interfacial Area Transport Equation" + M. Ishii, + S. Kim, + J Kelly, + Nuclear Engineering and Technology, Vol.37 No.6 December 2005 + \endverbatim + +SourceFiles + turbulentBreakUp.C + +\*---------------------------------------------------------------------------*/ + +#ifndef turbulentBreakUp_H +#define turbulentBreakUp_H + +#include "IATEsource.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + +/*---------------------------------------------------------------------------*\ + Class turbulentBreakUp Declaration +\*---------------------------------------------------------------------------*/ + +class turbulentBreakUp +: + public IATEsource +{ + // Private data + + dimensionedScalar Cti_; + dimensionedScalar WeCr_; + + +public: + + //- Runtime type information + TypeName("turbulentBreakUp"); + + // Constructors + + turbulentBreakUp + ( + const IATE& iate, + const dictionary& dict + ); + + + //- Destructor + virtual ~turbulentBreakUp() + {} + + + // Member Functions + + virtual tmp R() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace IATEsources +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C new file mode 100644 index 0000000000..66c324c761 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C @@ -0,0 +1,72 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "wakeEntrainmentCoalescence.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + defineTypeNameAndDebug(wakeEntrainmentCoalescence, 0); + addToRunTimeSelectionTable + ( + IATEsource, + wakeEntrainmentCoalescence, + dictionary + ); +} +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence:: +wakeEntrainmentCoalescence +( + const IATE& iate, + const dictionary& dict +) +: + IATEsource(iate), + Cwe_("Cwe", dimless, dict.lookup("Cwe")) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp +Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence::R() const +{ + return (-12)*phi()*Cwe_*cbrt(CD())*iate_.a()*Ur(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H new file mode 100644 index 0000000000..7c2df39912 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence + +Description + Bubble coalescence due to wake entrainment IATE source as defined in paper: + + \verbatim + "Development of Interfacial Area Transport Equation" + M. Ishii, + S. Kim, + J Kelly, + Nuclear Engineering and Technology, Vol.37 No.6 December 2005 + \endverbatim + +SourceFiles + wakeEntrainmentCoalescence.C + +\*---------------------------------------------------------------------------*/ + +#ifndef wakeEntrainmentCoalescence_H +#define wakeEntrainmentCoalescence_H + +#include "IATEsource.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace diameterModels +{ +namespace IATEsources +{ + +/*---------------------------------------------------------------------------*\ + Class wakeEntrainmentCoalescence Declaration +\*---------------------------------------------------------------------------*/ + +class wakeEntrainmentCoalescence +: + public IATEsource +{ + // Private data + + dimensionedScalar Cwe_; + + +public: + + //- Runtime type information + TypeName("wakeEntrainmentCoalescence"); + + // Constructors + + wakeEntrainmentCoalescence + ( + const IATE& iate, + const dictionary& dict + ); + + + //- Destructor + virtual ~wakeEntrainmentCoalescence() + {} + + + // Member Functions + + virtual tmp R() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace IATEsources +} // End namespace diameterModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C index 315f160893..233b9bf49e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C @@ -48,12 +48,12 @@ namespace diameterModels Foam::diameterModels::constant::constant ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ) : - diameterModel(dict, phase), - d_("d", dimLength, dict.lookup("d")) + diameterModel(diameterProperties, phase), + d_("d", dimLength, diameterProperties_.lookup("d")) {} @@ -84,4 +84,14 @@ Foam::tmp Foam::diameterModels::constant::d() const } +bool Foam::diameterModels::constant::read(const dictionary& phaseProperties) +{ + diameterModel::read(phaseProperties); + + diameterProperties_.lookup("d") >> d_; + + return true; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H index 071820c319..6f616c6f73 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H @@ -22,7 +22,7 @@ License along with OpenFOAM. If not, see . Class - Foam::constant + Foam::diameterModels::constant Description Constant dispersed-phase particle diameter model. @@ -69,7 +69,7 @@ public: //- Construct from components constant ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ); @@ -80,7 +80,11 @@ public: // Member Functions - tmp d() const; + //- Return the diameter as a field + virtual tmp d() const; + + //- Read diameterProperties dictionary + virtual bool read(const dictionary& diameterProperties); }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C index 2618815d24..55225147ac 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C @@ -38,11 +38,11 @@ namespace Foam Foam::diameterModel::diameterModel ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ) : - dict_(dict), + diameterProperties_(diameterProperties), phase_(phase) {} @@ -53,4 +53,18 @@ Foam::diameterModel::~diameterModel() {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::diameterModel::correct() +{} + + +bool Foam::diameterModel::read(const dictionary& phaseProperties) +{ + diameterProperties_ = phaseProperties.subDict(type() + "Coeffs"); + + return true; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H index c14f65e364..b40a6cb8ca 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H @@ -36,26 +36,27 @@ SourceFiles #ifndef diameterModel_H #define diameterModel_H -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - #include "dictionary.H" #include "phaseModel.H" #include "runTimeSelectionTables.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + namespace Foam { /*---------------------------------------------------------------------------*\ - Class diameterModel Declaration + Class diameterModel Declaration \*---------------------------------------------------------------------------*/ class diameterModel { + protected: // Protected data - const dictionary& dict_; + dictionary diameterProperties_; const phaseModel& phase_; @@ -73,10 +74,10 @@ public: diameterModel, dictionary, ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ), - (dict, phase) + (diameterProperties, phase) ); @@ -84,7 +85,7 @@ public: diameterModel ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ); @@ -97,15 +98,33 @@ public: static autoPtr New ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ); // Member Functions + //- Return the phase diameter properties dictionary + const dictionary& diameterProperties() const + { + return diameterProperties_; + } + + //- Return the phase + const phaseModel& phase() const + { + return phase_; + } + //- Return the phase mean diameter field virtual tmp d() const = 0; + + //- Correct the diameter field + virtual void correct(); + + //- Read phaseProperties dictionary + virtual bool read(const dictionary& phaseProperties) = 0; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C index 4b9e571503..aab9b469c7 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C @@ -48,13 +48,13 @@ namespace diameterModels Foam::diameterModels::isothermal::isothermal ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ) : - diameterModel(dict, phase), - d0_("d0", dimLength, dict.lookup("d0")), - p0_("p0", dimPressure, dict.lookup("p0")) + diameterModel(diameterProperties, phase), + d0_("d0", dimLength, diameterProperties_.lookup("d0")), + p0_("p0", dimPressure, diameterProperties_.lookup("p0")) {} @@ -77,4 +77,15 @@ Foam::tmp Foam::diameterModels::isothermal::d() const } +bool Foam::diameterModels::isothermal::read(const dictionary& phaseProperties) +{ + diameterModel::read(phaseProperties); + + diameterProperties_.lookup("d0") >> d0_; + diameterProperties_.lookup("p0") >> p0_; + + return true; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H index 0d155ad16e..6e997a210a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H @@ -22,7 +22,7 @@ License along with OpenFOAM. If not, see . Class - Foam::isothermal + Foam::diameterModels::isothermal Description Isothermal dispersed-phase particle diameter model. @@ -72,7 +72,7 @@ public: //- Construct from components isothermal ( - const dictionary& dict, + const dictionary& diameterProperties, const phaseModel& phase ); @@ -83,7 +83,11 @@ public: // Member Functions - tmp d() const; + //- Return the diameter field + virtual tmp d() const; + + //- Read phaseProperties dictionary + virtual bool read(const dictionary& phaseProperties); }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C index 697371c044..da43a49881 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C @@ -26,6 +26,8 @@ License #include "phaseModel.H" #include "twoPhaseSystem.H" #include "diameterModel.H" +#include "fvMatrix.H" +#include "PhaseIncompressibleTurbulenceModel.H" #include "dragModel.H" #include "heatTransferModel.H" #include "fixedValueFvPatchFields.H" @@ -73,6 +75,17 @@ Foam::phaseModel::phaseModel IOobject::AUTO_WRITE ), fluid.mesh() + ), + phiAlpha_ + ( + IOobject + ( + IOobject::groupName("alphaPhi", name_), + fluid.mesh().time().timeName(), + fluid.mesh() + ), + fluid.mesh(), + dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0) ) { thermo_->validate("phaseModel " + name_, "h", "e"); @@ -153,6 +166,16 @@ Foam::phaseModel::phaseModel phaseDict_, *this ); + + turbulence_ = + PhaseIncompressibleTurbulenceModel::New + ( + *this, + U_, + phiAlpha_, + phi(), + *this + ); } @@ -164,10 +187,39 @@ Foam::phaseModel::~phaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +const Foam::phaseModel& Foam::phaseModel::otherPhase() const +{ + return fluid_.otherPhase(*this); +} + + Foam::tmp Foam::phaseModel::d() const { return dPtr_().d(); } +Foam::PhaseIncompressibleTurbulenceModel& +Foam::phaseModel::turbulence() +{ + return turbulence_(); +} + +const Foam::PhaseIncompressibleTurbulenceModel& +Foam::phaseModel::turbulence() const +{ + return turbulence_(); +} + +void Foam::phaseModel::correct() +{ + return dPtr_->correct(); +} + +bool Foam::phaseModel::read(const dictionary& phaseProperties) +{ + phaseDict_ = phaseProperties.subDict(name_); + return dPtr_->read(phaseDict_); +} + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H index ed2a611546..7ca4255207 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H @@ -48,6 +48,9 @@ namespace Foam class twoPhaseSystem; class diameterModel; +template +class PhaseIncompressibleTurbulenceModel; + /*---------------------------------------------------------------------------*\ Class phaseModel Declaration @@ -74,12 +77,18 @@ class phaseModel //- Velocity volVectorField U_; - //- Fluxes + //- Volumetric flux of the phase + surfaceScalarField phiAlpha_; + + //- Volumetric flux of the phase autoPtr phiPtr_; //- Diameter model autoPtr dPtr_; + //- turbulence model + autoPtr > turbulence_; + public: @@ -99,19 +108,47 @@ public: // Member Functions + //- Return the name of this phase + const word& name() const + { + return name_; + } + //- Return the twoPhaseSystem to which this phase belongs const twoPhaseSystem& fluid() const { return fluid_; } - const word& name() const - { - return name_; - } + //- Return the other phase in this two-phase system + const phaseModel& otherPhase() const; + //- Return the Sauter-mean diameter tmp d() const; + //- Return the turbulence model + const PhaseIncompressibleTurbulenceModel& + turbulence() const; + + //- Return non-const access to the turbulence model + // for correction + PhaseIncompressibleTurbulenceModel& + turbulence(); + + //- Return the thermophysical model + const rhoThermo& thermo() const + { + return thermo_(); + } + + //- Return non-const access to the thermophysical model + // for correction + rhoThermo& thermo() + { + return thermo_(); + } + + //- Return the laminar viscosity tmp nu() const { return thermo_->nu(); @@ -123,57 +160,71 @@ public: return thermo_->nu(patchi); } + //- Return the thermal conductivity tmp kappa() const { return thermo_->kappa(); } + //- Return the specific heat capacity tmp Cp() const { return thermo_->Cp(); } + //- Return the density const volScalarField& rho() const { return thermo_->rho(); } - const rhoThermo& thermo() const - { - return thermo_(); - } - - rhoThermo& thermo() - { - return thermo_(); - } - + //- Return the velocity const volVectorField& U() const { return U_; } + //- Return non-const access to the velocity + // Used in the momentum equation volVectorField& U() { return U_; } + //- Return the volumetric flux const surfaceScalarField& phi() const { return phiPtr_(); } + //- Return non-const access to the volumetric flux surfaceScalarField& phi() { return phiPtr_(); } - //- Dummy correct - void correct() - {} + //- Return the volumetric flux of the phase + const surfaceScalarField& phiAlpha() const + { + return phiAlpha_; + } - //- Dummy read - bool read() + //- Return non-const access to the volumetric flux of the phase + surfaceScalarField& phiAlpha() + { + return phiAlpha_; + } + + //- Correct the phase properties + // other than the thermodynamics and turbulence + // which have special treatment + void correct(); + + //- Read phaseProperties dictionary + virtual bool read(const dictionary& phaseProperties); + + //- Dummy Read for transportModel + virtual bool read() { return true; } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 83902a7b7d..8ba3382b86 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -24,9 +24,19 @@ License \*---------------------------------------------------------------------------*/ #include "twoPhaseSystem.H" +#include "fvMatrix.H" +#include "PhaseIncompressibleTurbulenceModel.H" #include "surfaceInterpolate.H" -#include "fixedValueFvsPatchFields.H" +#include "MULES.H" +#include "subCycle.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcSnGrad.H" +#include "fvcFlux.H" #include "fvcCurl.H" +#include "fvmDdt.H" +#include "fvmLaplacian.H" +#include "fixedValueFvsPatchFields.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -63,6 +73,37 @@ Foam::twoPhaseSystem::twoPhaseSystem wordList(lookup("phases"))[1] ), + phi_ + ( + IOobject + ( + "phi", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->calcPhi() + ), + + dgdt_ + ( + IOobject + ( + "dgdt", + mesh.time().timeName(), + mesh + ), + pos(phase2_)*fvc::div(phi_)/max(phase2_, scalar(0.0001)) + ), + + sigma_ + ( + "sigma", + dimensionSet(1, 0, -2, 0, 0), + lookup("sigma") + ), + Cvm_ ( "Cvm", @@ -170,7 +211,7 @@ Foam::tmp Foam::twoPhaseSystem::U() const } -Foam::tmp Foam::twoPhaseSystem::phi() const +Foam::tmp Foam::twoPhaseSystem::calcPhi() const { return fvc::interpolate(phase1_)*phase1_.phi() @@ -366,18 +407,242 @@ Foam::tmp Foam::twoPhaseSystem::heatTransferCoeff() const } +void Foam::twoPhaseSystem::solve() +{ + const Time& runTime = mesh_.time(); + + volScalarField& alpha1 = phase1_; + volScalarField& alpha2 = phase2_; + + const surfaceScalarField& phi1 = phase1_.phi(); + const surfaceScalarField& phi2 = phase2_.phi(); + + const dictionary& alphaControls = mesh_.solverDict + ( + alpha1.name() + ); + + label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); + label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr"))); + Switch implicitPhasePressure + ( + alphaControls.lookupOrDefault("implicitPhasePressure", false) + ); + + word alphaScheme("div(phi," + alpha1.name() + ')'); + word alpharScheme("div(phir," + alpha1.name() + ')'); + + alpha1.correctBoundaryConditions(); + + + surfaceScalarField phic("phic", phi_); + surfaceScalarField phir("phir", phi1 - phi2); + + surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); + + tmp pPrimeByA; + + if (implicitPhasePressure) + { + const volScalarField& rAU1 = mesh_.lookupObject + ( + IOobject::groupName("rAU", phase1_.name()) + ); + const volScalarField& rAU2 = mesh_.lookupObject + ( + IOobject::groupName("rAU", phase2_.name()) + ); + + pPrimeByA = + fvc::interpolate((1.0/phase1_.rho()) + *rAU1*phase1_.turbulence().pPrime()) + + fvc::interpolate((1.0/phase2_.rho()) + *rAU2*phase2_.turbulence().pPrime()); + + surfaceScalarField phiP + ( + pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf() + ); + + phic += alpha1f*phiP; + phir += phiP; + } + + for (int acorr=0; acorr 0.0 && alpha1[celli] > 0.0) + { + Sp[celli] -= dgdt_[celli]*alpha1[celli]; + Su[celli] += dgdt_[celli]*alpha1[celli]; + } + else if (dgdt_[celli] < 0.0 && alpha1[celli] < 1.0) + { + Sp[celli] += dgdt_[celli]*(1.0 - alpha1[celli]); + } + } + + dimensionedScalar totalDeltaT = runTime.deltaT(); + if (nAlphaSubCycles > 1) + { + phase1_.phiAlpha() = + dimensionedScalar("0", phase1_.phiAlpha().dimensions(), 0); + } + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + surfaceScalarField alphaPhic1 + ( + fvc::flux + ( + phic, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + alpha1, + alpharScheme + ) + ); + + // Ensure that the flux at inflow BCs is preserved + forAll(alphaPhic1.boundaryField(), patchi) + { + fvsPatchScalarField& alphaPhic1p = + alphaPhic1.boundaryField()[patchi]; + + if (!alphaPhic1p.coupled()) + { + const scalarField& phi1p = phi1.boundaryField()[patchi]; + const scalarField& alpha1p = alpha1.boundaryField()[patchi]; + + forAll(alphaPhic1p, facei) + { + if (phi1p[facei] < 0) + { + alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei]; + } + } + } + } + + MULES::explicitSolve + ( + geometricOneField(), + alpha1, + phi_, + alphaPhic1, + Sp, + Su, + 1, + 0 + ); + + if (nAlphaSubCycles > 1) + { + phase1_.phiAlpha() += (runTime.deltaT()/totalDeltaT)*alphaPhic1; + } + else + { + phase1_.phiAlpha() = alphaPhic1; + } + } + + if (implicitPhasePressure) + { + fvScalarMatrix alpha1Eqn + ( + fvm::ddt(alpha1) - fvc::ddt(alpha1) + - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded") + ); + + alpha1Eqn.relax(); + alpha1Eqn.solve(); + + phase1_.phiAlpha() += alpha1Eqn.flux(); + } + + phase2_.phiAlpha() = phi_ - phase1_.phiAlpha(); + alpha2 = scalar(1) - alpha1; + + Info<< alpha1.name() << " volume fraction = " + << alpha1.weightedAverage(mesh_.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + } +} + + +void Foam::twoPhaseSystem::correct() +{ + phase1_.correct(); + phase2_.correct(); +} + + +void Foam::twoPhaseSystem::correctTurbulence() +{ + phase1_.turbulence().correct(); + phase2_.turbulence().correct(); +} + + bool Foam::twoPhaseSystem::read() { if (regIOobject::read()) { bool readOK = true; - readOK &= phase1_.read(); - readOK &= phase2_.read(); + readOK &= phase1_.read(*this); + readOK &= phase2_.read(*this); + lookup("sigma") >> sigma_; lookup("Cvm") >> Cvm_; lookup("Cl") >> Cl_; + // drag1_->read(*this); + // drag2_->read(*this); + + // heatTransfer1_->read(*this); + // heatTransfer2_->read(*this); + + lookup("dispersedPhase") >> dispersedPhase_; + lookup("residualPhaseFraction") >> residualPhaseFraction_; + lookup("residualSlip") >> residualSlip_; + return readOK; } else diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H index 3a7896ea02..0b501152c1 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H @@ -74,6 +74,12 @@ class twoPhaseSystem phaseModel phase1_; phaseModel phase2_; + surfaceScalarField phi_; + + volScalarField dgdt_; + + dimensionedScalar sigma_; + dimensionedScalar Cvm_; dimensionedScalar Cl_; @@ -88,6 +94,10 @@ class twoPhaseSystem dimensionedScalar residualSlip_; + //- Return the mixture flux + tmp calcPhi() const; + + public: // Constructors @@ -140,6 +150,28 @@ public: return phase2_; } + //- Return the mixture flux + const surfaceScalarField& phi() const + { + return phi_; + } + + //- Return the mixture flux + surfaceScalarField& phi() + { + return phi_; + } + + const volScalarField& dgdt() const + { + return dgdt_; + } + + volScalarField& dgdt() + { + return dgdt_; + } + const dragModel& drag1() const { return drag1_(); @@ -193,8 +225,11 @@ public: //- Return the mixture velocity tmp U() const; - //- Return the mixture flux - tmp phi() const; + //- Return the surface tension coefficient + dimensionedScalar sigma() const + { + return sigma_; + } //- Return the virtual-mass coefficient dimensionedScalar Cvm() const @@ -208,9 +243,14 @@ public: return Cl_; } - //- Dummy correct - void correct() - {} + //- Solve for the two-phase-fractions + void solve(); + + //- Correct two-phase properties other than turbulence + void correct(); + + //- Correct two-phase turbulence + void correctTurbulence(); //- Read base phaseProperties dictionary bool read(); diff --git a/applications/test/HashingSpeed/Test-HashingSpeed.C b/applications/test/HashingSpeed/Test-HashingSpeed.C index 1b42b5d34d..c3cbda702d 100644 --- a/applications/test/HashingSpeed/Test-HashingSpeed.C +++ b/applications/test/HashingSpeed/Test-HashingSpeed.C @@ -1016,7 +1016,6 @@ int i; // found somewhere in the 2nd addition uint32_t stroustrup (const char * s, int len) { uint32_t h; - int i; for (register int i=0; i < len; ++s, ++i) { diff --git a/applications/test/List/Test-List.C b/applications/test/List/Test-List.C index ca303d53b7..719a5ebc5c 100644 --- a/applications/test/List/Test-List.C +++ b/applications/test/List/Test-List.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -137,7 +137,7 @@ int main(int argc, char *argv[]) Info<< " std::list constructed from Foam::List: "; std::list::iterator it; - for (it=stlList.begin(); it != stlList.end(); it++) + for (it=stlList.begin(); it != stlList.end(); ++it) { Info<< *it << " "; } diff --git a/applications/test/string/Test-string.C b/applications/test/string/Test-string.C index 8494e8a104..f39be8c206 100644 --- a/applications/test/string/Test-string.C +++ b/applications/test/string/Test-string.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -79,7 +79,7 @@ int main(int argc, char *argv[]) while (iter != iter2) { Info<< *iter; - iter++; + ++iter; } Info<< "<\n"; diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index bf8eba6e0d..9db10ea85e 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -474,10 +474,12 @@ meshQualityControls // Advanced // Flags for optional output -// 0 : only write final meshes -// 1 : write intermediate meshes -// 2 : write volScalarField with cellLevel for postprocessing -// 4 : write current intersections as .obj files +// 0 : only write final meshes +// 1 : write intermediate meshes +// 2 : write volScalarField with cellLevel for postprocessing +// 4 : write current mesh intersections as .obj files +// 8 : write information about explicit feature edge refinement +// 16 : write information about layers debug 0; // Merge tolerance. Is fraction of overall bounding box of initial mesh. diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C index 7e395be879..b0c4d37c03 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C @@ -1,3 +1,28 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + #include "checkTopology.H" #include "polyMesh.H" #include "Time.H" @@ -9,6 +34,8 @@ #include "emptyPolyPatch.H" #include "processorPolyPatch.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + Foam::label Foam::checkTopology ( const polyMesh& mesh, @@ -287,6 +314,37 @@ Foam::label Foam::checkTopology rs ); ctr.write(); + + + // write cellSet for each region + PtrList cellRegions(rs.nRegions()); + for (label i = 0; i < rs.nRegions(); i++) + { + cellRegions.set + ( + i, + new cellSet + ( + mesh, + "region" + Foam::name(i), + mesh.nCells()/100 + ) + ); + } + + forAll(rs, i) + { + cellRegions[rs[i]].insert(i); + } + + for (label i = 0; i < rs.nRegions(); i++) + { + Info<< " <()) + << " cells to cellSet " << cellRegions[i].name() << endl; + + cellRegions[i].write(); + } } } @@ -437,3 +495,6 @@ Foam::label Foam::checkTopology return noFailedChecks; } + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun deleted file mode 100755 index 11f1cb983e..0000000000 --- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/sh -cd ${0%/*} || exit 1 # run from this directory - -# Source tutorial run functions -. $WM_PROJECT_DIR/bin/tools/RunFunctions - -# Get application name -application=`getApplication` - -runApplication blockMesh -# Make random cell numbering to get faceZone orientation random -runApplication renumberMesh -dict system/renumberMeshDict - -# Construct faceZone -runApplication topoSet - -runApplication decomposePar -force -cellDist -runParallel orientFaceZone 2 f0 '(10 0 0)' -runApplication reconstructParMesh -latestTime -mergeTol 1e-6 - - -# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt deleted file mode 100644 index df42877302..0000000000 --- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt +++ /dev/null @@ -1 +0,0 @@ -2013-09-16 To test orientFaceZone. Make random orientation. Orient afterwards. diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict deleted file mode 100644 index 60454a0a87..0000000000 --- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict +++ /dev/null @@ -1,135 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: dev | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - note "mesh decomposition control dictionary"; - object decomposeParDict; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -numberOfSubdomains 2; - -//- Keep owner and neighbour on same processor for faces in zones: -// preserveFaceZones (heater solid1 solid3); - -//- Keep owner and neighbour on same processor for faces in patches: -// (makes sense only for cyclic patches) -//preservePatches (cyclic_half0 cyclic_half1); - -//- Keep all of faceSet on a single processor. This puts all cells -// connected with a point, edge or face on the same processor. -// (just having face connected cells might not guarantee a balanced -// decomposition) -// The processor can be -1 (the decompositionMethod chooses the processor -// for a good load balance) or explicitly provided (upsets balance). -//singleProcessorFaceSets ((f0 -1)); - - -//- Use the volScalarField named here as a weight for each cell in the -// decomposition. For example, use a particle population field to decompose -// for a balanced number of particles in a lagrangian simulation. -// weightField dsmcRhoNMean; - -//method scotch; -method hierarchical; -// method simple; -// method metis; -// method manual; -// method multiLevel; -// method structured; // does 2D decomposition of structured mesh - -multiLevelCoeffs -{ - // Decomposition methods to apply in turn. This is like hierarchical but - // fully general - every method can be used at every level. - - level0 - { - numberOfSubdomains 64; - //method simple; - //simpleCoeffs - //{ - // n (2 1 1); - // delta 0.001; - //} - method scotch; - } - level1 - { - numberOfSubdomains 4; - method scotch; - } -} - -// Desired output - -simpleCoeffs -{ - n (2 1 1); - delta 0.001; -} - -hierarchicalCoeffs -{ - n (2 1 1); - delta 0.001; - order xyz; -} - -metisCoeffs -{ - /* - processorWeights - ( - 1 - 1 - 1 - 1 - ); - */ -} - -scotchCoeffs -{ - //processorWeights - //( - // 1 - // 1 - // 1 - // 1 - //); - //writeGraph true; - //strategy "b"; -} - -manualCoeffs -{ - dataFile "decompositionData"; -} - -structuredCoeffs -{ - // Patches to do 2D decomposition on. Structured mesh only; cells have - // to be in 'columns' on top of patches. - patches (bottomPatch); -} - -//// Is the case distributed? Note: command-line argument -roots takes -//// precedence -//distributed yes; -//// Per slave (so nProcs-1 entries) the directory above the case. -//roots -//( -// "/tmp" -// "/tmp" -//); - -// ************************************************************************* // diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict deleted file mode 100644 index addb04236e..0000000000 --- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict +++ /dev/null @@ -1,110 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: dev | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - note "mesh renumbering dictionary"; - object renumberMeshDict; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -// Write maps from renumbered back to original mesh -writeMaps true; - -// Optional entry: sort cells on coupled boundaries to last for use with -// e.g. nonBlockingGaussSeidel. -sortCoupledFaceCells false; - -// Optional entry: renumber on a block-by-block basis. It uses a -// blockCoeffs dictionary to construct a decompositionMethod to do -// a block subdivision) and then applies the renumberMethod to each -// block in turn. This can be used in large cases to keep the blocks -// fitting in cache with all the the cache misses bunched at the end. -// This number is the approximate size of the blocks - this gets converted -// to a number of blocks that is the input to the decomposition method. -//blockSize 1000; - -// Optional entry: sort points into internal and boundary points -//orderPoints false; - - - -//method CuthillMcKee; -//method Sloan; -//method manual; -method random; -//method structured; -//method spring; -//method zoltan; // only if compiled with zoltan support - -//CuthillMcKeeCoeffs -//{ -// // Reverse CuthillMcKee (RCM) or plain -// reverse true; -//} - -manualCoeffs -{ - // In system directory: new-to-original (i.e. order) labelIOList - dataFile "cellMap"; -} - - -// For extruded (i.e. structured in one direction) meshes -structuredCoeffs -{ - // Patches that mesh was extruded from. These determine the starting - // layer of cells - patches (movingWall); - // Method to renumber the starting layer of cells - method random; - - // Renumber in columns (depthFirst) or in layers - depthFirst true; - - // Optional: reverse ordering - //reverse false; -} - - -springCoeffs -{ - // Maximum jump of cell indices. Is fraction of number of cells - maxCo 0.01; - - // Limit the amount of movement; the fraction maxCo gets decreased - // with every iteration - freezeFraction 0.999; - - // Maximum number of iterations - maxIter 1000; -} - - -blockCoeffs -{ - method scotch; - //method hierarchical; - //hierarchicalCoeffs - //{ - // n (1 2 1); - // delta 0.001; - // order xyz; - //} -} - - -zoltanCoeffs -{ - ORDER_METHOD LOCAL_HSFC; -} - - -// ************************************************************************* // diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict b/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict deleted file mode 100644 index fca42b45f8..0000000000 --- a/applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict +++ /dev/null @@ -1,64 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: 2.2.1 | -| \\ / A nd | Web: www.OpenFOAM.org | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - object topoSetDict; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -actions -( - // Get 'T' of faces - { - name f0; - type faceSet; - action new; - source boxToFace; - sourceInfo - { - boxes - ( - (0.0499 -100 -100)(0.0501 100 100) - (0.03 0.0499 -100)(0.05 0.0501 100) - ); - } - } - -// // Pick up some cells -// { -// name c0; -// type cellSet; -// action new; -// source boxToCell; -// sourceInfo -// { -// boxes -// ( -// (-100 -100 -100)(0.05 0.05 100) -// (0.05 0.05 -100)(100 100 100) -// ); -// } -// } - - // Convert faceSet to faceZone - { - name f0; - type faceZoneSet; - action new; - source setToFaceZone; - sourceInfo - { - faceSet f0; - } - } -); - -// ************************************************************************* // diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict index 9dee4d1df9..5f0e5c0b39 100644 --- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict +++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict @@ -368,6 +368,8 @@ FoamFile // { // faceSet f0; // name of faceSet // cellSet c0; // name of cellSet of slave side +// flip false; // optional: flip the faceZone (so now the cellSet +// // is the master side) // } // // // Select based on surface. Orientation from normals on surface diff --git a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C index 3e9f8cd159..e0dda32b2a 100644 --- a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C +++ b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C @@ -46,8 +46,8 @@ int main(int argc, char *argv[]) ); argList::noBanner(); - argList::noParallel(); argList::addBoolOption("times", "list available times"); + argList::addBoolOption("latestTime", "list last time"); argList::addBoolOption ( "keywords", @@ -75,6 +75,15 @@ int main(int argc, char *argv[]) Info<< times[i].name() << endl; } } + else if (args.optionFound("latestTime")) + { + instantList times + ( + Foam::Time::findTimes(args.rootPath()/args.caseName()) + ); + + Info<< times.last().name() << endl; + } if (args.optionFound("dict")) { diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C index 21ba1cdfa1..b569601734 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C +++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C @@ -99,6 +99,40 @@ Usage // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +const labelIOList& procAddressing +( + const PtrList& procMeshList, + const label procI, + const word& name, + PtrList& procAddressingList +) +{ + const fvMesh& procMesh = procMeshList[procI]; + + if (!procAddressingList.set(procI)) + { + procAddressingList.set + ( + procI, + new labelIOList + ( + IOobject + ( + name, + procMesh.facesInstance(), + procMesh.meshSubDir, + procMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ) + ); + } + return procAddressingList[procI]; +} + + + int main(int argc, char *argv[]) { argList::addNote @@ -805,74 +839,29 @@ int main(int argc, char *argv[]) } const fvMesh& procMesh = procMeshList[procI]; + const labelIOList& faceProcAddressing = procAddressing + ( + procMeshList, + procI, + "faceProcAddressing", + faceProcAddressingList + ); - if (!faceProcAddressingList.set(procI)) - { - faceProcAddressingList.set - ( - procI, - new labelIOList - ( - IOobject - ( - "faceProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ) - ); - } - const labelIOList& faceProcAddressing = - faceProcAddressingList[procI]; + const labelIOList& cellProcAddressing = procAddressing + ( + procMeshList, + procI, + "cellProcAddressing", + cellProcAddressingList + ); - - if (!cellProcAddressingList.set(procI)) - { - cellProcAddressingList.set - ( - procI, - new labelIOList - ( - IOobject - ( - "cellProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ) - ); - } - const labelIOList& cellProcAddressing = - cellProcAddressingList[procI]; - - - if (!boundaryProcAddressingList.set(procI)) - { - boundaryProcAddressingList.set - ( - procI, - new labelIOList - ( - IOobject - ( - "boundaryProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ) - ); - } - const labelIOList& boundaryProcAddressing = - boundaryProcAddressingList[procI]; + const labelIOList& boundaryProcAddressing = procAddressing + ( + procMeshList, + procI, + "boundaryProcAddressing", + boundaryProcAddressingList + ); // FV fields @@ -959,27 +948,13 @@ int main(int argc, char *argv[]) || pointTensorFields.size() ) { - if (!pointProcAddressingList.set(procI)) - { - pointProcAddressingList.set - ( - procI, - new labelIOList - ( - IOobject - ( - "pointProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ) - ); - } - const labelIOList& pointProcAddressing = - pointProcAddressingList[procI]; + const labelIOList& pointProcAddressing = procAddressing + ( + procMeshList, + procI, + "pointProcAddressing", + pointProcAddressingList + ); const pointMesh& procPMesh = pointMesh::New(procMesh); diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C index 4eda411482..923d9ae1cf 100644 --- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C +++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C @@ -38,6 +38,7 @@ License #include "cellSet.H" #include "faceSet.H" #include "pointSet.H" +#include "uniformDimensionedFields.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -195,6 +196,60 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets) } + autoPtr cellLevelPtr; + { + IOobject io + ( + "cellLevel", + facesInstance(), + polyMesh::meshSubDir, + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE + ); + if (io.headerOk()) + { + Info<< "Reading hexRef8 data : " << io.name() << endl; + cellLevelPtr.reset(new labelIOList(io)); + } + } + autoPtr pointLevelPtr; + { + IOobject io + ( + "pointLevel", + facesInstance(), + polyMesh::meshSubDir, + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE + ); + if (io.headerOk()) + { + Info<< "Reading hexRef8 data : " << io.name() << endl; + pointLevelPtr.reset(new labelIOList(io)); + } + } + autoPtr level0EdgePtr; + { + IOobject io + ( + "level0Edge", + facesInstance(), + polyMesh::meshSubDir, + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE + ); + if (io.headerOk()) + { + Info<< "Reading hexRef8 data : " << io.name() << endl; + level0EdgePtr.reset(new uniformDimensionedScalarField(io)); + } + } + + + label maxProcCells = 0; label totProcFaces = 0; label maxProcPatches = 0; @@ -767,8 +822,28 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets) procMesh.write(); + // Write points if pointsInstance differing from facesInstance + if (facesInstancePointsPtr_.valid()) + { + pointIOField pointsInstancePoints + ( + IOobject + ( + "points", + pointsInstance(), + polyMesh::meshSubDir, + procMesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + xferMove(procPoints) + ); + pointsInstancePoints.write(); + } + // Decompose any sets if (decomposeSets) { forAll(cellSets, i) @@ -813,25 +888,67 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets) } - // Write points if pointsInstance differing from facesInstance - if (facesInstancePointsPtr_.valid()) + // hexRef8 data + if (cellLevelPtr.valid()) { - pointIOField pointsInstancePoints + labelIOList ( IOobject ( - "points", - pointsInstance(), + cellLevelPtr().name(), + facesInstance(), polyMesh::meshSubDir, procMesh, IOobject::NO_READ, - IOobject::NO_WRITE, - false + IOobject::AUTO_WRITE ), - xferMove(procPoints) - ); - pointsInstancePoints.write(); + UIndirectList