diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files index d45448f27a..b85aac6d2c 100644 --- a/src/fvOptions/Make/files +++ b/src/fvOptions/Make/files @@ -12,6 +12,7 @@ $(generalSources)/semiImplicitSource/semiImplicitSource.C derivedSources=sources/derived $(derivedSources)/actuationDiskSource/actuationDiskSource.C +$(derivedSources)/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C $(derivedSources)/explicitPorositySource/explicitPorositySource.C $(derivedSources)/MRFSource/MRFSource.C $(derivedSources)/pressureGradientExplicitSource/pressureGradientExplicitSource.C @@ -27,7 +28,8 @@ $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModel.C $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModelNew.C $(derivedSources)/rotorDiskSource/trimModel/fixed/fixedTrim.C $(derivedSources)/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C -$(derivedSources)/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C +$(derivedSources)/solidificationMeltingSource/solidificationMeltingSource.C +$(derivedSources)/solidificationMeltingSource/solidificationMeltingSourceIO.C interRegion = sources/interRegion $(interRegion)/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C new file mode 100644 index 0000000000..e17b8c88d5 --- /dev/null +++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C @@ -0,0 +1,358 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013-2014 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "solidificationMeltingSource.H" +#include "fvMatrices.H" +#include "basicThermo.H" +#include "uniformDimensionedFields.H" +#include "zeroGradientFvPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "geometricOneField.H" + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +namespace Foam +{ + template<> + const char* NamedEnum + < + fv::solidificationMeltingSource::thermoMode, + 2 + >::names[] = + { + "thermo", + "lookup" + }; + + namespace fv + { + defineTypeNameAndDebug(solidificationMeltingSource, 0); + + addToRunTimeSelectionTable + ( + option, + solidificationMeltingSource, + dictionary + ); + } +} + +const Foam::NamedEnum + Foam::fv::solidificationMeltingSource::thermoModeTypeNames_; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::fv::solidificationMeltingSource::solveField +( + const word& fieldName +) const +{ + bool result = true; + + switch (mode_) + { + case mdThermo: + { + const basicThermo& thermo = + mesh_.lookupObject("thermophysicalProperties"); + + if (fieldName != thermo.he().name()) + { + result = false; + } + break; + } + case mdLookup: + { + if (fieldName != TName_) + { + result = false; + } + break; + } + default: + { + FatalErrorIn + ( + "bool Foam::fv::solidificationMeltingSource::solveField" + "(" + "const word&" + ") const" + ) + << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_] + << abort(FatalError); + } + } + return result; +} + + +Foam::tmp +Foam::fv::solidificationMeltingSource::Cp() const +{ + switch (mode_) + { + case mdThermo: + { + const basicThermo& thermo = + mesh_.lookupObject("thermophysicalProperties"); + + return thermo.Cp(); + break; + } + case mdLookup: + { + if (CpName_ == "CpRef") + { + scalar CpRef = readScalar(coeffs_.lookup("CpRef")); + + return tmp + ( + new volScalarField + ( + IOobject + ( + name_ + ":Cp", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar + ( + "Cp", + dimEnergy/dimMass/dimTemperature, + CpRef + ), + zeroGradientFvPatchScalarField::typeName + ) + ); + } + else + { + return mesh_.lookupObject(CpName_); + } + + break; + } + default: + { + FatalErrorIn + ( + "Foam::tmp " + "Foam::fv::solidificationMeltingSource::Cp() const" + ) + << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_] + << abort(FatalError); + } + } + + return tmp(NULL); +} + + +Foam::vector Foam::fv::solidificationMeltingSource::g() const +{ + if (mesh_.foundObject("g")) + { + const uniformDimensionedVectorField& value = + mesh_.lookupObject("g"); + return value.value(); + } + else + { + return coeffs_.lookup("g"); + } +} + + +void Foam::fv::solidificationMeltingSource::update(const volScalarField& Cp) +{ + if (curTimeIndex_ == mesh_.time().timeIndex()) + { + return; + } + + if (debug) + { + Info<< type() << ": " << name_ << " - updating phase indicator" << endl; + } + + // update old time alpha1 field + alpha1_.oldTime(); + + const volScalarField& T = mesh_.lookupObject(TName_); + + forAll(cells_, i) + { + label cellI = cells_[i]; + + scalar Tc = T[cellI]; + scalar Cpc = Cp[cellI]; + scalar alpha1New = alpha1_[cellI] + relax_*Cpc*(Tc - Tmelt_)/L_; + + alpha1_[cellI] = max(0, min(alpha1New, 1)); + deltaT_[i] = Tc - Tmelt_; + } + + alpha1_.correctBoundaryConditions(); + + curTimeIndex_ = mesh_.time().timeIndex(); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::fv::solidificationMeltingSource::solidificationMeltingSource +( + const word& sourceName, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh +) +: + option(sourceName, modelType, dict, mesh), + Tmelt_(readScalar(coeffs_.lookup("Tmelt"))), + L_(readScalar(coeffs_.lookup("L"))), + relax_(coeffs_.lookupOrDefault("relax", 0.9)), + mode_(thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"))), + rhoRef_(readScalar(coeffs_.lookup("rhoRef"))), + TName_(coeffs_.lookupOrDefault("TName", "T")), + CpName_(coeffs_.lookupOrDefault("CpName", "Cp")), + UName_(coeffs_.lookupOrDefault("UName", "U")), + phiName_(coeffs_.lookupOrDefault("phiName", "phi")), + Cu_(coeffs_.lookupOrDefault("Cu", 100000)), + q_(coeffs_.lookupOrDefault("q", 0.001)), + beta_(readScalar(coeffs_.lookup("beta"))), + alpha1_ + ( + IOobject + ( + name_ + ":alpha1", + mesh.time().timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("alpha1", dimless, 0.0), + zeroGradientFvPatchScalarField::typeName + ), + curTimeIndex_(-1), + deltaT_(cells_.size(), 0) +{ + fieldNames_.setSize(1, "source"); + applied_.setSize(1, false); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::fv::solidificationMeltingSource::alwaysApply() const +{ + return true; +} + + +void Foam::fv::solidificationMeltingSource::addSup +( + fvMatrix& eqn, + const label fieldI +) +{ + apply(geometricOneField(), eqn); +} + + +void Foam::fv::solidificationMeltingSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + apply(rho, eqn); +} + + +void Foam::fv::solidificationMeltingSource::addSup +( + fvMatrix& eqn, + const label fieldI +) +{ + const volVectorField& U = eqn.psi(); + + if (U.name() != UName_) + { + return; + } + + if (debug) + { + Info<< type() << ": applying source to " << UName_ << endl; + } + + const volScalarField Cp(this->Cp()); + + update(Cp); + + vector g = this->g(); + + scalarField& Sp = eqn.diag(); + vectorField& Su = eqn.source(); + const scalarField& V = mesh_.V(); + + forAll(cells_, i) + { + label cellI = cells_[i]; + + scalar Vc = V[cellI]; + scalar alpha1c = alpha1_[cellI]; + + scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_); + vector Sb = rhoRef_*g*beta_*deltaT_[i]; + + Sp[cellI] += Vc*S; + Su[cellI] += Vc*Sb; + } +} + + +void Foam::fv::solidificationMeltingSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // momentum source uses a Boussinesq approximation - redirect + addSup(eqn, fieldI); +} + + +// ************************************************************************* // diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H new file mode 100644 index 0000000000..a8d6f95143 --- /dev/null +++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H @@ -0,0 +1,280 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013-2014 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::fv::solidificationMeltingSource + +Description + This source is designed to model the effect of solidification and melting + processes, e.g. windhield defrosting. The phase change occurs at the + melting temperature, \c Tmelt. + + The presence of the solid phase in the flow field is incorporated into the + model as a momentum porosity contribution; the energy associated with the + phase change is added as an enthalpy contribution. + + Based on the references: + + 1. V.R. Voller and C. Prakash, A fixed grid numerical modelling + methodology for convection-diffusion mushy phase-change problems, + Int. J. Heat Mass Transfer 30(8):17091719, 1987. + 2. C.R. Swaminathan. and V.R. Voller, A general enthalpy model for + modeling solidification processes, Metallurgical Transactions + 23B:651664, 1992. + + The model generates a field \c \:alpha1 which can be visualised to + to show the melt distribution as a fraction [0-1] + + \heading Source usage + Example usage: + \verbatim + solidificationMeltingSource1 + { + type solidificationMeltingSource; + active on; + selectionMode cellZone; + cellZone iceZone; + + solidificationMeltingSourceCoeffs + { + Tmelt 273; + L 334000; + thermoMode thermo; + beta 50e-6; + rhoRef 800; + } + } + \endverbatim + + Where: + \table + Property | Description | Required | Default value + Tmelt | Melting temperature [K] | yes | + L | Latent heat of fusion [J/kg] | yes | + relax | Relaxation coefficient [0-1] | no | 0.9 + thermoMode | Thermo mode [thermo|lookup] | yes | + rhoRef | Reference (solid) density | yes | + rhoName | Name of density field | no | rho + TName | Name of temperature field | no | T + CpName | Name of specific heat capacity field | no | Cp + UName | Name of velocity field | no | U + phiName | Name of flux field | no | phi + Cu | Model coefficient | no | 100000 + q | Model coefficient | no | 0.001 + beta | Thermal expansion coefficient [1/K] | yes | + g | Accelerartion due to gravity | no | + \endtable + +SourceFiles + solidificationMeltingSource.C + solidificationMeltingSourceIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef solidificationMeltingSource_H +#define solidificationMeltingSource_H + +#include "fvMesh.H" +#include "volFields.H" +#include "fvOption.H" +#include "NamedEnum.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class solidificationMeltingSource Declaration +\*---------------------------------------------------------------------------*/ + +class solidificationMeltingSource +: + public option +{ +public: + + enum thermoMode + { + mdThermo, + mdLookup + }; + + static const NamedEnum thermoModeTypeNames_; + + +private: + + // Private data + + //- Temperature at which melting occurs [K] + scalar Tmelt_; + + //- Latent heat of fusion [J/kg] + scalar L_; + + //- Phase fraction under-relaxation coefficient + scalar relax_; + + //- Thermodynamics mode + thermoMode mode_; + + //- Reference density - typically the solid density + scalar rhoRef_; + + //- Name of temperature field - default = "T" (optional) + word TName_; + + //- Name of specific heat capacity field - default = "Cp" (optional) + word CpName_; + + //- Name of velocity field - default = "U" (optional) + word UName_; + + //- Name of flux field - default = "phi" (optional) + word phiName_; + + //- Mushy region momentum sink coefficient [1/s]; default = 10^5 + scalar Cu_; + + //- Coefficient used in porosity calc - default = 0.001 + scalar q_; + + //- Thermal expansion coefficient [1/K] + scalar beta_; + + //- Phase fraction indicator field + volScalarField alpha1_; + + //- Current time index (used for updating) + label curTimeIndex_; + + //- Temperature change cached for source calculation when alpha1 updated + scalarField deltaT_; + + + // Private Member Functions + + //- Flag to indicate whether to solve for given field + bool solveField(const word& fieldName) const; + + //- Return the specific heat capacity field + tmp Cp() const; + + //- Return the gravity vector + vector g() const; + + //- Update the model + void update(const volScalarField& Cp); + + //- Helper function to apply to the energy equation + template + void apply(const RhoFieldType& rho, fvMatrix& eqn); + + //- Disallow default bitwise copy construct + solidificationMeltingSource(const solidificationMeltingSource&); + + //- Disallow default bitwise assignment + void operator=(const solidificationMeltingSource&); + + +public: + + //- Runtime type information + TypeName("solidificationMeltingSource"); + + + // Constructors + + //- Construct from explicit source name and mesh + solidificationMeltingSource + ( + const word& sourceName, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + + // Member Functions + + //- Flag to bypass the apply flag list checking + virtual bool alwaysApply() const; + + + // Add explicit and implicit contributions + + //- Add explicit contribution to enthalpy equation + virtual void addSup(fvMatrix& eqn, const label fieldI); + + //- Add implicit contribution to momentum equation + virtual void addSup(fvMatrix& eqn, const label fieldI); + + + // Add explicit and implicit contributions to compressible equation + + //- Add explicit contribution to compressible enthalpy equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Add implicit contribution to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + + // I-O + + //- Write the source properties + virtual void writeData(Ostream&) const; + + //- Read source dictionary + virtual bool read(const dictionary& dict); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "solidificationMeltingSourceTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C new file mode 100644 index 0000000000..d6173d1227 --- /dev/null +++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C @@ -0,0 +1,68 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013-2014 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "solidificationMeltingSource.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::fv::solidificationMeltingSource::writeData(Ostream& os) const +{ + os << indent << name_ << endl; + dict_.write(os); +} + + +bool Foam::fv::solidificationMeltingSource::read(const dictionary& dict) +{ + if (option::read(dict)) + { + coeffs_.lookup("Tmelt") >> Tmelt_; + coeffs_.lookup("L") >> L_; + + coeffs_.readIfPresent("relax", relax_); + + mode_ = thermoModeTypeNames_.read(coeffs_.lookup("thermoMode")); + + coeffs_.lookup("rhoRef") >> rhoRef_; + coeffs_.readIfPresent("TName", TName_); + coeffs_.readIfPresent("UName", UName_); + + coeffs_.readIfPresent("Cu", Cu_); + coeffs_.readIfPresent("q", q_); + + coeffs_.readIfPresent("beta", beta_); + + return true; + } + else + { + return false; + } + + return false; +} + + +// ************************************************************************* // diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C new file mode 100644 index 0000000000..2dfddb76c6 --- /dev/null +++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2014 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "fvMatrices.H" +#include "fvcDdt.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::fv::solidificationMeltingSource::apply +( + const RhoFieldType& rho, + fvMatrix& eqn +) +{ + if (!solveField(eqn.psi().name())) + { + return; + } + + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + const volScalarField Cp(this->Cp()); + + update(Cp); + + dimensionedScalar L("L", dimEnergy/dimMass, L_); + + // contributions added to rhs of solver equation + if (eqn.psi().dimensions() == dimTemperature) + { + // isothermal phase change - only include time derivative +// eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); + eqn -= L/Cp*(fvc::ddt(rho, alpha1_)); + } + else + { + // isothermal phase change - only include time derivative +// eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); + eqn -= L*(fvc::ddt(rho, alpha1_)); + } +} + + +// ************************************************************************* //