diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean new file mode 100755 index 0000000000000000000000000000000000000000..8098c051604f5b67e75d5401178f2be68d974be1 --- /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 0000000000000000000000000000000000000000..dd002ee0698abdaffc004b7a81ee1b57fd996c2e --- /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 0000000000000000000000000000000000000000..8e6a4cb0de4cf91898a21b08e16cbcd9a74b20a6 --- /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 0000000000000000000000000000000000000000..0f47e6979edee7a5294823ff42e883991eb5279e --- /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 0000000000000000000000000000000000000000..a7ee3a7bb4ba946d65f522cf716e0e1d66781e4b --- /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 0000000000000000000000000000000000000000..38cfebde6f3c8412ae037cd96f68907953a30850 --- /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 0000000000000000000000000000000000000000..1fd68ffb790d4bccb99e870e10f9f10577bc0d00 --- /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 <http://www.gnu.org/licenses/>. + +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 0000000000000000000000000000000000000000..09018b603f0089b685ddcca2497877b1eb6db9d1 --- /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 <http://www.gnu.org/licenses/>. + +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 0000000000000000000000000000000000000000..5c7f723a875a0ef6fe8d93b8679f01b3cc5ee7f2 --- /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<compressible::turbulenceModel> 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 0000000000000000000000000000000000000000..93e6eb9e3dd84ac695c3791a6b3fbd45114f73b5 --- /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 0000000000000000000000000000000000000000..eab8cce15dde3c729b68afecb9510af3bb5967e1 --- /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 0000000000000000000000000000000000000000..6872ae0321e74659359ee24d45561f5f01e04fca --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#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<scalar, volMesh>& iF +) +: + zeroGradientFvPatchScalarField(p, iF) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const fvPatchFieldMapper& mapper +) +: + zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper), + thetaProps_(gcpsf.thetaProps_) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const dictionary& dict +) +: + zeroGradientFvPatchScalarField(p, iF), + thetaProps_(dict.lookup("thetaProperties")) +{ + evaluate(); +} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const DimensionedField<scalar, volMesh>& 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 0000000000000000000000000000000000000000..53ce61d24f316391a4941dbeed2415255b92525d --- /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 <http://www.gnu.org/licenses/>. + +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<scalar, volMesh>& + ); + + //- Construct from patch, internal field and dictionary + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const dictionary& + ); + + //- Construct by mapping given alphaContactAngleFvPatchScalarField + // onto a new patch + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp<fvPatchScalarField> clone() const + { + return tmp<fvPatchScalarField> + ( + new alphaContactAngleFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const DimensionedField<scalar, volMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp<fvPatchScalarField> clone + ( + const DimensionedField<scalar, volMesh>& iF + ) const + { + return tmp<fvPatchScalarField> + ( + 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 0000000000000000000000000000000000000000..163955056439cd7615e0acf7e020fb3f026ab643 --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#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<phaseModel>, 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<phaseModel>, phases_, phasei) + { + phasei().correct(); + } + + PtrDictionary<phaseModel>::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<phaseModel>, phases_, phasei) + { + phasei().thermo().rho() += phasei().thermo().psi()*dp; + } +} + + +bool Foam::multiphaseMixtureThermo::incompressible() const +{ + bool ico = true; + + forAllConstIter(PtrDictionary<phaseModel>, phases_, phase) + { + ico &= phase().thermo().incompressible(); + } + + return ico; +} + + +bool Foam::multiphaseMixtureThermo::isochoric() const +{ + bool iso = true; + + forAllConstIter(PtrDictionary<phaseModel>, phases_, phase) + { + iso &= phase().thermo().incompressible(); + } + + return iso; +} + + +Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::he +( + const volScalarField& p, + const volScalarField& T +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> 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::scalarField> Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const labelList& cells +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::scalarField> Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::hc() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> thc(phasei()*phasei().thermo().hc()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + thc() += phasei()*phasei().thermo().hc(); + } + + return thc; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const labelList& cells +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const label patchi +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rho() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> trho(phasei()*phasei().thermo().rho()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trho() += phasei()*phasei().thermo().rho(); + } + + return trho; +} + + +Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cp() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCp() += phasei()*phasei().thermo().Cp(); + } + + return tCp; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::Cv() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCv() += phasei()*phasei().thermo().Cv(); + } + + return tCv; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::gamma() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tgamma() += phasei()*phasei().thermo().gamma(); + } + + return tgamma; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::Cpv() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpv() += phasei()*phasei().thermo().Cpv(); + } + + return tCpv; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::CpByCpv() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpByCpv() += phasei()*phasei().thermo().CpByCpv(); + } + + return tCpByCpv; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::kappa() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappa() += phasei()*phasei().thermo().kappa(); + } + + return tkappa; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappa +( + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::kappaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat); + } + + return tkappaEff; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::alphaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + talphaEff() += phasei()*phasei().thermo().alphaEff(alphat); + } + + return talphaEff; +} + + +Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<scalarField> 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::volScalarField> Foam::multiphaseMixtureThermo::rCv() const +{ + PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin(); + + tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trCv() += phasei()/phasei().thermo().Cv(); + } + + return trCv; +} + + +Foam::tmp<Foam::surfaceScalarField> +Foam::multiphaseMixtureThermo::surfaceTensionForce() const +{ + tmp<surfaceScalarField> 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<phaseModel>, phases_, phase1) + { + const phaseModel& alpha1 = phase1(); + + PtrDictionary<phaseModel>::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<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + solveAlphas(cAlpha); + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_; + } + + rhoPhi_ = rhoPhiSum; + } + else + { + solveAlphas(cAlpha); + } +} + + +Foam::tmp<Foam::surfaceVectorField> 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::surfaceScalarField> 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<alphaContactAngleFvPatchScalarField>(gbf[patchi])) + { + const alphaContactAngleFvPatchScalarField& acap = + refCast<const alphaContactAngleFvPatchScalarField>(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::volScalarField> Foam::multiphaseMixtureThermo::K +( + const phaseModel& alpha1, + const phaseModel& alpha2 +) const +{ + tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2); + + correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField()); + + // Simple expression for curvature + return -fvc::div(tnHatfv & mesh_.Sf()); +} + + +Foam::tmp<Foam::volScalarField> +Foam::multiphaseMixtureThermo::nearInterface() const +{ + tmp<volScalarField> tnearInt + ( + new volScalarField + ( + IOobject + ( + "nearInterface", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("nearInterface", dimless, 0.0) + ) + ); + + forAllConstIter(PtrDictionary<phaseModel>, 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<surfaceScalarField> phiAlphaCorrs(phases_.size()); + int phasei = 0; + + forAllIter(PtrDictionary<phaseModel>, phases_, phase) + { + phaseModel& alpha = phase(); + + phiAlphaCorrs.set + ( + phasei, + new surfaceScalarField + ( + fvc::flux + ( + phi_, + alpha, + alphaScheme + ) + ) + ); + + surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei]; + + forAllIter(PtrDictionary<phaseModel>, 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<phaseModel>, phases_, phase) + { + phaseModel& alpha = phase(); + + surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; + phiAlpha += upwind<scalar>(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<phaseModel>, 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 0000000000000000000000000000000000000000..546048894f494a3cb025f9c569f400df14001485 --- /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 <http://www.gnu.org/licenses/>. + +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<word> + { + public: + + class hash + : + public Hash<interfacePair> + { + 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<word>(alpha1Name, alpha2Name) + {} + + interfacePair(const phaseModel& alpha1, const phaseModel& alpha2) + : + Pair<word>(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<phaseModel> phases_; + + const fvMesh& mesh_; + const volVectorField& U_; + const surfaceScalarField& phi_; + + surfaceScalarField rhoPhi_; + + volScalarField alphas_; + + typedef HashTable<scalar, interfacePair, interfacePair::hash> + 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<surfaceVectorField> nHatfv + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + tmp<surfaceScalarField> nHatf + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + void correctContactAngle + ( + const phaseModel& alpha1, + const phaseModel& alpha2, + surfaceVectorField::GeometricBoundaryField& nHatb + ) const; + + tmp<volScalarField> 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<phaseModel>& phases() const + { + return phases_; + } + + //- Return non-const access to the phases + PtrDictionary<phaseModel>& 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<volScalarField> he + ( + const volScalarField& p, + const volScalarField& T + ) const; + + //- Enthalpy/Internal energy for cell-set [J/kg] + virtual tmp<scalarField> he + ( + const scalarField& p, + const scalarField& T, + const labelList& cells + ) const; + + //- Enthalpy/Internal energy for patch [J/kg] + virtual tmp<scalarField> he + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Chemical enthalpy [J/kg] + virtual tmp<volScalarField> hc() const; + + //- Temperature from enthalpy/internal energy for cell-set + virtual tmp<scalarField> 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<scalarField> 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<volScalarField> rho() const; + + //- Heat capacity at constant pressure [J/kg/K] + virtual tmp<volScalarField> Cp() const; + + //- Heat capacity at constant pressure for patch [J/kg/K] + virtual tmp<scalarField> Cp + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant volume [J/kg/K] + virtual tmp<volScalarField> Cv() const; + + //- Heat capacity at constant volume for patch [J/kg/K] + virtual tmp<scalarField> Cv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- gamma = Cp/Cv [] + virtual tmp<volScalarField> gamma() const; + + //- gamma = Cp/Cv for patch [] + virtual tmp<scalarField> gamma + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant pressure/volume [J/kg/K] + virtual tmp<volScalarField> Cpv() const; + + //- Heat capacity at constant pressure/volume for patch [J/kg/K] + virtual tmp<scalarField> Cpv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity ratio [] + virtual tmp<volScalarField> CpByCpv() const; + + //- Heat capacity ratio for patch [] + virtual tmp<scalarField> 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<volScalarField> kappa() const; + + //- Thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp<scalarField> kappa + ( + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp<volScalarField> kappaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp<scalarField> kappaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp<volScalarField> alphaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp<scalarField> alphaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + + //- Return the phase-averaged reciprocal Cv + tmp<volScalarField> rCv() const; + + tmp<surfaceScalarField> surfaceTensionForce() const; + + //- Indicator of the proximity of the interface + // Field values are 1 near and 0 away for the interface. + tmp<volScalarField> 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 0000000000000000000000000000000000000000..1559f25a483f68ac1d45ac2b64511bb5a1ba2da9 --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#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> Foam::phaseModel::clone() const +{ + notImplemented("phaseModel::clone() const"); + return autoPtr<phaseModel>(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 0000000000000000000000000000000000000000..66d0ac8d63a8d6c008f18a09a75dc98d6c149d99 --- /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 <http://www.gnu.org/licenses/>. + +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<rhoThermo> thermo_; + volScalarField dgdt_; + + +public: + + // Constructors + + //- Construct from components + phaseModel + ( + const word& phaseName, + const volScalarField& p, + const volScalarField& T + ); + + //- Return clone + autoPtr<phaseModel> 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<phaseModel> operator()(Istream& is) const + { + return autoPtr<phaseModel>(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 0000000000000000000000000000000000000000..e42a54aabdac39c54c2db0468736bfda5052dc26 --- /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<fixedFluxPressureFvPatchScalarField> + ( + p_rgh.boundaryField(), + ( + phiHbyA.boundaryField() + - (mesh.Sf().boundaryField() & U.boundaryField()) + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + + PtrList<fvScalarMatrix> p_rghEqnComps(multiphaseProperties.phases().size()); + + label phasei = 0; + forAllConstIter + ( + PtrDictionary<phaseModel>, + 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<fvScalarMatrix> p_rghEqnComp; + + phasei = 0; + forAllConstIter + ( + PtrDictionary<phaseModel>, + multiphaseProperties.phases(), + phase + ) + { + tmp<fvScalarMatrix> 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<phaseModel>, + 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/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T new file mode 100644 index 0000000000000000000000000000000000000000..1dc67cd149c827026b09c4f1738f89f3a7073cc1 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/T @@ -0,0 +1,52 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 293; + +boundaryField +{ + leftWall + { + type fixedValue; + value $internalField; + } + rightWall + { + type fixedValue; + value $internalField; + } + lowerWall + { + type fixedValue; + value uniform 293; + } + atmosphere + { + type inletOutlet; + phi rho*phi; + inletValue $internalField; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U new file mode 100644 index 0000000000000000000000000000000000000000..7ea3a0c32328d93fd01ba3c580c7d73b54473d22 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/U @@ -0,0 +1,51 @@ +/*--------------------------------*- 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 volVectorField; + location "0"; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + leftWall + { + type fixedValue; + value uniform (0 0 0); + } + rightWall + { + type fixedValue; + value uniform (0 0 0); + } + lowerWall + { + type fixedValue; + value uniform (0 0 0); + } + atmosphere + { + type pressureInletOutletVelocity; + value uniform (0 0 0); + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air new file mode 100644 index 0000000000000000000000000000000000000000..be8307fdf27a9e1bedd029b9f857452fd7a09c7e --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.air @@ -0,0 +1,79 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object alpha.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + leftWall + { + type alphaContactAngle; + thetaProperties + ( + ( water air ) 90 0 0 0 + ( oil air ) 90 0 0 0 + ( mercury air ) 90 0 0 0 + ( water oil ) 90 0 0 0 + ( water mercury ) 90 0 0 0 + ( oil mercury ) 90 0 0 0 + ); + value uniform 0; + } + rightWall + { + type alphaContactAngle; + thetaProperties + ( + ( water air ) 90 0 0 0 + ( oil air ) 90 0 0 0 + ( mercury air ) 90 0 0 0 + ( water oil ) 90 0 0 0 + ( water mercury ) 90 0 0 0 + ( oil mercury ) 90 0 0 0 + ); + value uniform 1; + } + lowerWall + { + type alphaContactAngle; + thetaProperties + ( + ( water air ) 90 0 0 0 + ( oil air ) 90 0 0 0 + ( mercury air ) 90 0 0 0 + ( water oil ) 90 0 0 0 + ( water mercury ) 90 0 0 0 + ( oil mercury ) 90 0 0 0 + ); + value uniform 0; + } + atmosphere + { + type inletOutlet; + inletValue uniform 1; + value uniform 1; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury new file mode 100644 index 0000000000000000000000000000000000000000..d224de9509a3b75ab8edca916d9fc9b7cd80b884 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.mercury @@ -0,0 +1,49 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object alpha.mercury; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + leftWall + { + type zeroGradient; + } + rightWall + { + type zeroGradient; + } + lowerWall + { + type zeroGradient; + } + atmosphere + { + type inletOutlet; + inletValue uniform 0; + value uniform 0; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil new file mode 100644 index 0000000000000000000000000000000000000000..bfcff63aedd3359e9ec036b69b1390822c9e548e --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.oil @@ -0,0 +1,49 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object alpha.oil; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + leftWall + { + type zeroGradient; + } + rightWall + { + type zeroGradient; + } + lowerWall + { + type zeroGradient; + } + atmosphere + { + type inletOutlet; + inletValue uniform 0; + value uniform 0; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water new file mode 100644 index 0000000000000000000000000000000000000000..cf96bb9d9f359617557d1c5166fc49bea0b3c8dc --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alpha.water @@ -0,0 +1,49 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object alpha.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + leftWall + { + type zeroGradient; + } + rightWall + { + type zeroGradient; + } + lowerWall + { + type zeroGradient; + } + atmosphere + { + type inletOutlet; + inletValue uniform 0; + value uniform 0; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas new file mode 100644 index 0000000000000000000000000000000000000000..10e3ea7a0724cf102fe0ee5f8d3b7b7c7c3d7ade --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/alphas @@ -0,0 +1,47 @@ +/*--------------------------------*- 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 volScalarField; + location "0"; + object alphas; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + leftWall + { + type zeroGradient; + } + rightWall + { + type zeroGradient; + } + lowerWall + { + type zeroGradient; + } + atmosphere + { + type zeroGradient; + } + defaultFaces + { + type empty; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p new file mode 100644 index 0000000000000000000000000000000000000000..26edbd338c0f88246d8975be0a4ec68723459c10 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p @@ -0,0 +1,53 @@ +/*--------------------------------*- 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 volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + leftWall + { + type calculated; + value $internalField; + } + + rightWall + { + type calculated; + value $internalField; + } + + lowerWall + { + type calculated; + value $internalField; + } + + atmosphere + { + type calculated; + value $internalField; + } + + defaultFaces + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh new file mode 100644 index 0000000000000000000000000000000000000000..16835050d0f6090661d306d9d9e371b29d422ce6 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0.org/p_rgh @@ -0,0 +1,59 @@ +/*--------------------------------*- 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 volScalarField; + object p_rgh; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + leftWall + { + type fixedFluxPressure; + value $internalField; + } + + rightWall + { + type fixedFluxPressure; + value $internalField; + } + + lowerWall + { + type fixedFluxPressure; + value $internalField; + } + + atmosphere + { + type totalPressure; + p0 uniform 1e5; + U U; + phi phi; + rho rho; + psi none; + gamma 1; + value $internalField; + } + + defaultFaces + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean new file mode 100755 index 0000000000000000000000000000000000000000..cd78b26599c3d4f7cd89a2197d3e52eb5f83087b --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean @@ -0,0 +1,11 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +cleanCase + +\rm -rf 0 + +# ----------------------------------------------------------------- end-of-file diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun new file mode 100755 index 0000000000000000000000000000000000000000..887344985c0a184737b5ab1f75f0709c953212b1 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun @@ -0,0 +1,17 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +# Set application name +application=`getApplication` + +\rm -rf 0 +cp -r 0.org 0 + +runApplication blockMesh +runApplication setFields +runApplication $application + +# ----------------------------------------------------------------- end-of-file diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g new file mode 100644 index 0000000000000000000000000000000000000000..e0ac2653b5b370ad62f6770588121d30cac51627 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g @@ -0,0 +1,22 @@ +/*--------------------------------*- 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 uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value ( 0 -9.81 0 ); + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict new file mode 100644 index 0000000000000000000000000000000000000000..11344be6ac848c78c9368b3e09948934cb22b570 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/blockMeshDict @@ -0,0 +1,108 @@ +/*--------------------------------*- 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; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 0.146; + +vertices +( + (0 0 0) + (2 0 0) + (2.16438 0 0) + (4 0 0) + (0 0.32876 0) + (2 0.32876 0) + (2.16438 0.32876 0) + (4 0.32876 0) + (0 4 0) + (2 4 0) + (2.16438 4 0) + (4 4 0) + (0 0 0.1) + (2 0 0.1) + (2.16438 0 0.1) + (4 0 0.1) + (0 0.32876 0.1) + (2 0.32876 0.1) + (2.16438 0.32876 0.1) + (4 0.32876 0.1) + (0 4 0.1) + (2 4 0.1) + (2.16438 4 0.1) + (4 4 0.1) +); + +blocks +( + hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1) + hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1) + hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1) + hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1) + hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1) +); + +edges +( +); + +boundary +( + leftWall + { + type wall; + faces + ( + (0 12 16 4) + (4 16 20 8) + ); + } + rightWall + { + type wall; + faces + ( + (7 19 15 3) + (11 23 19 7) + ); + } + lowerWall + { + type wall; + faces + ( + (0 1 13 12) + (1 5 17 13) + (5 6 18 17) + (2 14 18 6) + (2 3 15 14) + ); + } + atmosphere + { + type patch; + faces + ( + (8 20 21 9) + (9 21 22 10) + (10 22 23 11) + ); + } +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary new file mode 100644 index 0000000000000000000000000000000000000000..1b4dbb60aaeaec2d5cfe40e3d4a35843d35b44a2 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/polyMesh/boundary @@ -0,0 +1,53 @@ +/*--------------------------------*- 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 polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +5 +( + leftWall + { + type wall; + nFaces 50; + startFace 4432; + } + rightWall + { + type wall; + nFaces 50; + startFace 4482; + } + lowerWall + { + type wall; + nFaces 62; + startFace 4532; + } + atmosphere + { + type patch; + nFaces 46; + startFace 4594; + } + defaultFaces + { + type empty; + inGroups 1(empty); + nFaces 4536; + startFace 4640; + } +) + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties new file mode 100644 index 0000000000000000000000000000000000000000..129c3a278ec17aca791224cde4c213984bc7ef22 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties @@ -0,0 +1,33 @@ +/*--------------------------------*- 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; + location "constant"; + object thermophysicalProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +phases (water oil mercury air); + +pMin pMin [ 1 -1 -2 0 0 0 0 ] 10000; + +sigmas +( + (air water) 0.07 + (air oil) 0.07 + (air mercury) 0.07 + (water oil) 0.07 + (water mercury) 0.07 + (oil mercury) 0.07 +); + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air new file mode 100644 index 0000000000000000000000000000000000000000..befc0aeae449cccc50a24e955083b88b914aba67 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.air @@ -0,0 +1,49 @@ +/*--------------------------------*- 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; + location "constant"; + object thermophysicalProperties.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectGas; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + nMoles 1; + molWeight 28.9; + } + thermodynamics + { + Cp 1007; + Hf 0; + } + transport + { + mu 1.84e-05; + Pr 0.7; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury new file mode 100644 index 0000000000000000000000000000000000000000..e90070ef1318deeee4aa18b9912519d548a8eed0 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.mercury @@ -0,0 +1,54 @@ +/*--------------------------------*- 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; + location "constant"; + object thermophysicalProperties.mercury; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectFluid; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + nMoles 1; + molWeight 200.59; + } + equationOfState + { + R 6818; + rho0 13529; + } + thermodynamics + { + Cp 139; + Hf 0; + } + transport + { + mu 1.522e-3; + Pr 0.022; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil new file mode 100644 index 0000000000000000000000000000000000000000..0bcdc33f4cbe6ce39d16ce37169bf10aedde5c9e --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.oil @@ -0,0 +1,54 @@ +/*--------------------------------*- 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; + location "constant"; + object thermophysicalProperties.oil; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectFluid; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + nMoles 1; + molWeight 100.21; + } + equationOfState + { + R 3564; + rho0 684; + } + thermodynamics + { + Cp 2240; + Hf 0; + } + transport + { + mu 3.76e-4; + Pr 6; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water new file mode 100644 index 0000000000000000000000000000000000000000..91e7adc381bccc0fa94958480d65a585a269bee8 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/thermophysicalProperties.water @@ -0,0 +1,54 @@ +/*--------------------------------*- 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; + location "constant"; + object thermophysicalProperties.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectFluid; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + nMoles 1; + molWeight 18.0; + } + equationOfState + { + R 7255; + rho0 1027; + } + thermodynamics + { + Cp 4195; + Hf 0; + } + transport + { + mu 3.645e-4; + Pr 2.289; + } +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties new file mode 100644 index 0000000000000000000000000000000000000000..c2c3b28a1b4e8f4a2cae55f58bd61f9b1a67b488 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/turbulenceProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- 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; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict new file mode 100644 index 0000000000000000000000000000000000000000..c2993ed4911547974bb0345305e85c404919205c --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict @@ -0,0 +1,56 @@ +/*--------------------------------*- 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; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application compressibleMultiphaseInterFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 10; + +deltaT 0.001; + +writeControl adjustableRunTime; + +writeInterval 0.05; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep yes; + +maxCo 0.5; +maxAlphaCo 0.5; + +maxDeltaT 1; + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict new file mode 100644 index 0000000000000000000000000000000000000000..e7e490bf74b1a3f58a34923b4f98569f1b09e483 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict @@ -0,0 +1,45 @@ +/*--------------------------------*- 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; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 4; + +method simple; + +simpleCoeffs +{ + n ( 2 2 1 ); + delta 0.001; +} + +hierarchicalCoeffs +{ + n ( 1 1 1 ); + delta 0.001; + order xyz; +} + +manualCoeffs +{ + dataFile ""; +} + +distributed no; + +roots ( ); + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes new file mode 100644 index 0000000000000000000000000000000000000000..9f3fac1dce597032e84957e7b639b1152fb3cad8 --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes @@ -0,0 +1,64 @@ +/*--------------------------------*- 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; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + div(rho*phi,U) Gauss upwind; + div(phi,alpha) Gauss vanLeer; + div(phirb,alpha) Gauss interfaceCompression; + "div\(phi,.*rho.*\)" Gauss upwind; + div(rho*phi,T) Gauss upwind; + div(rho*phi,K) Gauss upwind; + div(phi,p) Gauss upwind; + div((muEff*dev2(T(grad(U))))) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + pcorr; + p_rgh; + "alpha.*"; +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution new file mode 100644 index 0000000000000000000000000000000000000000..6f88113f81b96e4c2f29a1ba31c9e3ad83287e7a --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution @@ -0,0 +1,120 @@ +/*--------------------------------*- 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; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + pcorr + { + solver PCG; + preconditioner + { + preconditioner GAMG; + tolerance 1e-05; + relTol 0; + smoother GaussSeidel; + nPreSweeps 0; + nPostSweeps 2; + nFinestSweeps 2; + cacheAgglomeration off; + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 2; + } + tolerance 1e-05; + relTol 0; + maxIter 100; + } + + ".*(rho|rhoFinal)" + { + solver diagonal; + } + + p_rgh + { + solver GAMG; + tolerance 1e-07; + relTol 0.05; + smoother GaussSeidel; + nPreSweeps 0; + nPostSweeps 2; + nFinestSweeps 2; + cacheAgglomeration on; + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 1; + } + + p_rghFinal + { + solver PCG; + preconditioner + { + preconditioner GAMG; + tolerance 1e-07; + relTol 0; + nVcycles 2; + smoother GaussSeidel; + nPreSweeps 0; + nPostSweeps 2; + nFinestSweeps 2; + cacheAgglomeration on; + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 1; + } + tolerance 1e-07; + relTol 0; + maxIter 20; + } + + "(U|T|k|B|nuTilda)" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-08; + relTol 0.1; + nSweeps 1; + } + + "(U|T|k|B|nuTilda)Final" + { + $U; + relTol 0; + } +} + +PIMPLE +{ + nCorrectors 2; + nNonOrthogonalCorrectors 0; + nAlphaSubCycles 4; + cAlpha 2; +} + +relaxationFactors +{ + fields + { + } + equations + { + "U.*" 1; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict new file mode 100644 index 0000000000000000000000000000000000000000..0a139a98ecd8a2e9d84aaf2bde4c6bb1306f923e --- /dev/null +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict @@ -0,0 +1,65 @@ +/*--------------------------------*- 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; + location "system"; + object setFieldsDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +defaultFieldValues +( + volScalarFieldValue alpha.air 1 + volScalarFieldValue alpha.water 0 + volScalarFieldValue alpha.oil 0 + volScalarFieldValue alpha.mercury 0 + volVectorFieldValue U ( 0 0 0 ) +); + +regions +( + boxToCell + { + box ( 0 0 -1 ) ( 0.1461 0.292 1 ); + fieldValues + ( + volScalarFieldValue alpha.water 1 + volScalarFieldValue alpha.oil 0 + volScalarFieldValue alpha.mercury 0 + volScalarFieldValue alpha.air 0 + ); + } + boxToCell + { + box ( 0.1461 0 -1 ) ( 0.2922 0.292 1 ); + fieldValues + ( + volScalarFieldValue alpha.water 0 + volScalarFieldValue alpha.oil 1 + volScalarFieldValue alpha.mercury 0 + volScalarFieldValue alpha.air 0 + ); + } + boxToCell + { + box ( 0 0 -1 ) ( 0.1461 0.1 1 ); + fieldValues + ( + volScalarFieldValue alpha.water 0 + volScalarFieldValue alpha.oil 0 + volScalarFieldValue alpha.mercury 1 + volScalarFieldValue alpha.air 0 + ); + } +); + + +// ************************************************************************* //