diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H index 2bd8277caa..77a04592e3 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H @@ -5,7 +5,7 @@ surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); - surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, 0.0))); + surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); tmp pPrimeByA; diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index 4e5686f3c3..bb2658a3ef 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -349,7 +349,7 @@ Foam::RASModels::kineticTheoryModel::divDevRhoReff void Foam::RASModels::kineticTheoryModel::correct() { // Local references - volScalarField alpha(max(this->alpha_, 0.0)); + volScalarField alpha(max(this->alpha_, scalar(0))); const volScalarField& rho = phase_.rho(); const surfaceScalarField& alphaPhi = this->alphaPhi_; const volVectorField& U = this->U_; diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C index da1d658f5d..8231cc32eb 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C @@ -88,7 +88,7 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime { volScalarField aByaMax ( - cbrt(min(max(alpha, 1e-3), alphaMinFriction)/alphaMax) + cbrt(min(max(alpha, scalar(1e-3)), alphaMinFriction)/alphaMax) ); return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax)); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C index 2cf5da4021..1f31befdd2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C @@ -52,6 +52,9 @@ makeBaseTurbulenceModel #include "kEpsilon.H" makeRASModel(kEpsilon); +#include "mixtureKEpsilon.H" +makeRASModel(mixtureKEpsilon); + #include "LaheyKEpsilon.H" makeRASModel(LaheyKEpsilon); diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C index b3eb2a9e1f..fe3f48e6be 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C @@ -106,7 +106,7 @@ bool Foam::GeometricField::readIfPresent() ( "GeometricField::readIfPresent()" ) << "read option IOobject::MUST_READ or MUST_READ_IF_MODIFIED" - << "suggests that a read constructor for field " << this->name() + << " suggests that a read constructor for field " << this->name() << " would be more appropriate." << endl; } else if (this->readOpt() == IOobject::READ_IF_PRESENT && this->headerOk()) diff --git a/src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.C b/src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.C index 1ffdd0553a..cadaafc818 100644 --- a/src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.C +++ b/src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.C @@ -200,7 +200,7 @@ NicenoKEqn::phaseTransferCoeff() const return ( - max(alphaInversion_ - alpha, 0.0) + max(alphaInversion_ - alpha, scalar(0)) *rho *min ( diff --git a/src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.C b/src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.C index e3a3b2f58d..6613d81636 100644 --- a/src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.C +++ b/src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.C @@ -137,7 +137,7 @@ continuousGasKEqn::phaseTransferCoeff() const return ( - max(alphaInversion_ - alpha, 0.0) + max(alphaInversion_ - alpha, scalar(0)) *rho *min ( diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.C b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.C index 2194280174..111e5bb8ba 100644 --- a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.C +++ b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.C @@ -83,6 +83,16 @@ LaheyKEpsilon::LaheyKEpsilon ) ), + C3_ + ( + dimensioned::lookupOrAddToDict + ( + "Cp", + this->coeffDict_, + this->C2_.value() + ) + ), + Cmub_ ( dimensioned::lookupOrAddToDict @@ -110,6 +120,7 @@ bool LaheyKEpsilon::read() { alphaInversion_.readIfPresent(this->coeffDict()); Cp_.readIfPresent(this->coeffDict()); + C3_.readIfPresent(this->coeffDict()); Cmub_.readIfPresent(this->coeffDict()); return true; @@ -207,7 +218,7 @@ LaheyKEpsilon::phaseTransferCoeff() const return ( - max(alphaInversion_ - alpha, 0.0) + max(alphaInversion_ - alpha, scalar(0)) *rho *min(gasTurbulence.epsilon()/gasTurbulence.k(), 1.0/U.time().deltaT()) ); @@ -244,7 +255,7 @@ tmp LaheyKEpsilon::epsilonSource() const const volScalarField phaseTransferCoeff(this->phaseTransferCoeff()); return - alpha*rho*this->C2_*this->epsilon_*bubbleG()/this->k_ + alpha*rho*this->C3_*this->epsilon_*bubbleG()/this->k_ + phaseTransferCoeff*gasTurbulence.epsilon() - fvm::Sp(phaseTransferCoeff, this->epsilon_); } diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H index a73255b3e0..7fc8509324 100644 --- a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H +++ b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H @@ -110,6 +110,7 @@ protected: dimensionedScalar alphaInversion_; dimensionedScalar Cp_; + dimensionedScalar C3_; dimensionedScalar Cmub_; diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C index a237a8998c..6fc1832bdf 100644 --- a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C +++ b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C @@ -123,7 +123,7 @@ void continuousGasKEpsilon::correctNut() volScalarField thetal(liquidTurbulence.k()/liquidTurbulence.epsilon()); volScalarField thetag((1.0/(18*liquid.nu()))*sqr(gas.d())); - volScalarField expThetar(exp(min(thetal/thetag, 50.0))); + volScalarField expThetar(exp(min(thetal/thetag, scalar(50)))); volScalarField omega(sqr(expThetar - 1)/(sqr(expThetar) - 1)); nutEff_ = omega*liquidTurbulence.nut(); @@ -163,7 +163,15 @@ continuousGasKEpsilon::nuEff() const { volScalarField blend ( - max(min((this->alpha_ - 0.5)/(alphaInversion_ - 0.5), 1.0), 0.0) + max + ( + min + ( + (this->alpha_ - scalar(0.5))/(alphaInversion_ - 0.5), + scalar(1) + ), + scalar(0) + ) ); return tmp @@ -210,7 +218,7 @@ continuousGasKEpsilon::phaseTransferCoeff() const return ( - max(alphaInversion_ - alpha, 0.0) + max(alphaInversion_ - alpha, scalar(0)) *rho *min ( diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C new file mode 100644 index 0000000000..4cd773db4b --- /dev/null +++ b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C @@ -0,0 +1,698 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "mixtureKEpsilon.H" +#include "bound.H" +#include "twoPhaseSystem.H" +#include "fixedValueFvPatchFields.H" +#include "inletOutletFvPatchFields.H" +#include "fvmSup.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace RASModels +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +mixtureKEpsilon::mixtureKEpsilon +( + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName, + const word& type +) +: + eddyViscosity > + ( + type, + alpha, + rho, + U, + alphaPhi, + phi, + transport, + propertiesName + ), + + liquidTurbulencePtr_(NULL), + + Cmu_ + ( + dimensioned::lookupOrAddToDict + ( + "Cmu", + this->coeffDict_, + 0.09 + ) + ), + C1_ + ( + dimensioned::lookupOrAddToDict + ( + "C1", + this->coeffDict_, + 1.44 + ) + ), + C2_ + ( + dimensioned::lookupOrAddToDict + ( + "C2", + this->coeffDict_, + 1.92 + ) + ), + C3_ + ( + dimensioned::lookupOrAddToDict + ( + "C3", + this->coeffDict_, + C2_.value() + ) + ), + Cp_ + ( + dimensioned::lookupOrAddToDict + ( + "Cp", + this->coeffDict_, + 0.25 + ) + ), + sigmak_ + ( + dimensioned::lookupOrAddToDict + ( + "sigmak", + this->coeffDict_, + 1.0 + ) + ), + sigmaEps_ + ( + dimensioned::lookupOrAddToDict + ( + "sigmaEps", + this->coeffDict_, + 1.3 + ) + ), + + k_ + ( + IOobject + ( + IOobject::groupName("k", U.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ), + epsilon_ + ( + IOobject + ( + IOobject::groupName("epsilon", U.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ) +{ + bound(k_, this->kMin_); + bound(epsilon_, this->epsilonMin_); + + if (type == typeName) + { + this->printCoeffs(type); + } +} + + +template +wordList mixtureKEpsilon::epsilonBoundaryTypes +( + const volScalarField& epsilon +) const +{ + const volScalarField::GeometricBoundaryField& ebf = epsilon.boundaryField(); + + wordList ebt = ebf.types(); + + forAll(ebf, patchi) + { + if (isA(ebf[patchi])) + { + ebt[patchi] = fixedValueFvPatchScalarField::typeName; + } + } + + return ebt; +} + + +template +void mixtureKEpsilon::correctInletOutlet +( + volScalarField& vsf, + const volScalarField& refVsf +) const +{ + volScalarField::GeometricBoundaryField& bf = vsf.boundaryField(); + const volScalarField::GeometricBoundaryField& refBf = + refVsf.boundaryField(); + + forAll(bf, patchi) + { + if + ( + isA(bf[patchi]) + && isA(refBf[patchi]) + ) + { + refCast + (bf[patchi]).refValue() = + refCast + (refBf[patchi]).refValue(); + } + } +} + + +template +void mixtureKEpsilon::initMixtureFields() +{ + if (rhom_.valid()) return; + + // Local references to gas-phase properties + const volScalarField& kg = this->k_; + const volScalarField& epsilong = this->epsilon_; + + // Local references to liquid-phase properties + mixtureKEpsilon& turbc = this->liquidTurbulence(); + const volScalarField& kl = turbc.k_; + const volScalarField& epsilonl = turbc.epsilon_; + + Ct2_.set + ( + new volScalarField + ( + IOobject + ( + "Ct2", + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + Ct2() + ) + ); + + rhom_.set + ( + new volScalarField + ( + IOobject + ( + "rhom", + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + rhom() + ) + ); + + km_.set + ( + new volScalarField + ( + IOobject + ( + "km", + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mix(kl, kg), + kl.boundaryField().types() + ) + ); + correctInletOutlet(km_(), kl); + + epsilonm_.set + ( + new volScalarField + ( + IOobject + ( + "epsilonm", + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mix(epsilonl, epsilong), + epsilonBoundaryTypes(epsilonl) + ) + ); + correctInletOutlet(epsilonm_(), epsilonl); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool mixtureKEpsilon::read() +{ + if (eddyViscosity >::read()) + { + Cmu_.readIfPresent(this->coeffDict()); + C1_.readIfPresent(this->coeffDict()); + C2_.readIfPresent(this->coeffDict()); + C3_.readIfPresent(this->coeffDict()); + Cp_.readIfPresent(this->coeffDict()); + sigmak_.readIfPresent(this->coeffDict()); + sigmaEps_.readIfPresent(this->coeffDict()); + + return true; + } + else + { + return false; + } +} + + +template +void mixtureKEpsilon::correctNut() +{ + this->nut_ = Cmu_*sqr(k_)/epsilon_; + this->nut_.correctBoundaryConditions(); +} + + +template +mixtureKEpsilon& +mixtureKEpsilon::liquidTurbulence() const +{ + if (!liquidTurbulencePtr_) + { + const volVectorField& U = this->U_; + + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + const transportModel& liquid = fluid.otherPhase(gas); + + liquidTurbulencePtr_ = + &const_cast&> + ( + U.db().lookupObject > + ( + IOobject::groupName + ( + turbulenceModel::propertiesName, + liquid.name() + ) + ) + ); + } + + return *liquidTurbulencePtr_; +} + + +template +tmp mixtureKEpsilon::Ct2() const +{ + const mixtureKEpsilon& liquidTurbulence = + this->liquidTurbulence(); + + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + const transportModel& liquid = fluid.otherPhase(gas); + + const volScalarField& alphag = this->alpha_; + + volScalarField magUr(mag(liquidTurbulence.U() - this->U())); + + volScalarField beta + ( + (6*this->Cmu_/(4*sqrt(3.0/2.0))) + *alphag*fluid.drag(gas).K(magUr)/liquid.rho() + *(liquidTurbulence.k_/liquidTurbulence.epsilon_) + ); + volScalarField Ct0((3 + beta)/(1 + beta + 2*gas.rho()/liquid.rho())); + volScalarField fAlphad((180 + (-4.71e3 + 4.26e4*alphag)*alphag)*alphag); + + return sqr(1 + (Ct0 - 1)*exp(-fAlphad)); +} + + +template +tmp mixtureKEpsilon::rholEff() const +{ + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + return fluid.otherPhase(gas).rho(); +} + + +template +tmp mixtureKEpsilon::rhogEff() const +{ + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + return + gas.rho() + + fluid.Cvm()*fluid.otherPhase(gas).rho(); +} + + +template +tmp mixtureKEpsilon::rhom() const +{ + const volScalarField& alphag = this->alpha_; + const volScalarField& alphal = this->liquidTurbulence().alpha_; + + return alphal*rholEff() + alphag*rhogEff(); +} + + +template +tmp mixtureKEpsilon::mix +( + const volScalarField& fc, + const volScalarField& fd +) const +{ + const volScalarField& alphag = this->alpha_; + const volScalarField& alphal = this->liquidTurbulence().alpha_; + + return (alphal*rholEff()*fc + alphag*rhogEff()*fd)/rhom_(); +} + + +template +tmp mixtureKEpsilon::mixU +( + const volScalarField& fc, + const volScalarField& fd +) const +{ + const volScalarField& alphag = this->alpha_; + const volScalarField& alphal = this->liquidTurbulence().alpha_; + + return + (alphal*rholEff()*fc + alphag*rhogEff()*Ct2_()*fd) + /(alphal*rholEff() + alphag*rhogEff()*Ct2_()); +} + + +template +tmp mixtureKEpsilon::mixFlux +( + const surfaceScalarField& fc, + const surfaceScalarField& fd +) const +{ + const volScalarField& alphag = this->alpha_; + const volScalarField& alphal = this->liquidTurbulence().alpha_; + + surfaceScalarField alphalf(fvc::interpolate(alphal)); + surfaceScalarField alphagf(fvc::interpolate(alphag)); + + surfaceScalarField rholEfff(fvc::interpolate(rholEff())); + surfaceScalarField rhogEfff(fvc::interpolate(rhogEff())); + + return + fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_())) + *(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd); +} + + +template +tmp mixtureKEpsilon::bubbleG() const +{ + const mixtureKEpsilon& liquidTurbulence = + this->liquidTurbulence(); + + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + const transportModel& liquid = fluid.otherPhase(gas); + + volScalarField magUr(mag(liquidTurbulence.U() - this->U())); + + // Lahey model + tmp bubbleG + ( + Cp_ + *sqr(liquid)*liquid.rho() + *( + pow3(magUr) + + pow(fluid.drag(gas).K(magUr)*gas.d()/liquid.rho(), 3.0/4.0) + *pow(magUr, 9.0/4.0) + ) + *gas + /gas.d() + ); + + // Simple model + // tmp bubbleG + // ( + // Cp_*sqr(liquid)*gas*fluid.drag(gas).K(magUr)*sqr(magUr) + // ); + + return bubbleG; +} + + +template +tmp mixtureKEpsilon::kSource() const +{ + return fvm::Su(bubbleG(), km_()); +} + + +template +tmp mixtureKEpsilon::epsilonSource() const +{ + return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_()); +} + + +template +void mixtureKEpsilon::correct() +{ + const transportModel& gas = this->transport(); + const twoPhaseSystem& fluid = gas.fluid(); + + // Only solve the mixture turbulence for the gas-phase + if (&gas != &fluid.phase1()) + { + // This is the liquid phase but check the model for the gas-phase + // is consistent + this->liquidTurbulence(); + + return; + } + + if (!this->turbulence_) + { + return; + } + + // Initialise the mixture fields if they have not yet been constructed + initMixtureFields(); + + // Local references to gas-phase properties + const surfaceScalarField& phig = this->phi_; + const volVectorField& Ug = this->U_; + const volScalarField& alphag = this->alpha_; + volScalarField& kg = this->k_; + volScalarField& epsilong = this->epsilon_; + volScalarField& nutg = this->nut_; + + // Local references to liquid-phase properties + mixtureKEpsilon& liquidTurbulence = + this->liquidTurbulence(); + const surfaceScalarField& phil = liquidTurbulence.phi_; + const volVectorField& Ul = liquidTurbulence.U_; + const volScalarField& alphal = liquidTurbulence.alpha_; + volScalarField& kl = liquidTurbulence.k_; + volScalarField& epsilonl = liquidTurbulence.epsilon_; + volScalarField& nutl = liquidTurbulence.nut_; + + // Local references to mixture properties + volScalarField& rhom = rhom_(); + volScalarField& km = km_(); + volScalarField& epsilonm = epsilonm_(); + + eddyViscosity >::correct(); + + // Update the effective mixture density + rhom = this->rhom(); + + // Mixture flux + surfaceScalarField phim("phim", mixFlux(phil, phig)); + + // Mixture velocity divergence + volScalarField divUm + ( + mixU + ( + fvc::div(fvc::absolute(phil, Ul)), + fvc::div(fvc::absolute(phig, Ug)) + ) + ); + + tmp Gc; + { + tmp tgradUl = fvc::grad(Ul); + Gc = tmp + ( + new volScalarField + ( + this->GName(), + nutl*(tgradUl() && dev(twoSymm(tgradUl()))) + ) + ); + tgradUl.clear(); + + // Update k, epsilon and G at the wall + kl.boundaryField().updateCoeffs(); + epsilonl.boundaryField().updateCoeffs(); + } + + tmp Gd; + { + tmp tgradUg = fvc::grad(Ug); + Gd = tmp + ( + new volScalarField + ( + this->GName(), + nutg*(tgradUg() && dev(twoSymm(tgradUg()))) + ) + ); + tgradUg.clear(); + + // Update k, epsilon and G at the wall + kg.boundaryField().updateCoeffs(); + epsilong.boundaryField().updateCoeffs(); + } + + // Mixture turbulence generation + volScalarField Gm(mix(Gc, Gd)); + + // Mixture turbulence viscosity + volScalarField nutm(mixU(nutl, nutg)); + + // Update the mixture k and epsilon boundary conditions + km == mix(kl, kg); + bound(km, this->kMin_); + epsilonm == mix(epsilonl, epsilong); + bound(epsilonm, this->epsilonMin_); + + // Dissipation equation + tmp epsEqn + ( + fvm::ddt(rhom, epsilonm) + + fvm::div(phim, epsilonm) + - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), epsilonm) + - fvm::laplacian(DepsilonEff(rhom*nutm), epsilonm) + == + C1_*rhom*Gm*epsilonm/km + - fvm::SuSp(((2.0/3.0)*C1_)*rhom*divUm, epsilonm) + - fvm::Sp(C2_*rhom*epsilonm/km, epsilonm) + + epsilonSource() + ); + + epsEqn().relax(); + + epsEqn().boundaryManipulate(epsilonm.boundaryField()); + + solve(epsEqn); + bound(epsilonm, this->epsilonMin_); + + + // Turbulent kinetic energy equation + tmp kmEqn + ( + fvm::ddt(rhom, km) + + fvm::div(phim, km) + - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), km) + - fvm::laplacian(DkEff(rhom*nutm), km) + == + rhom*Gm + - fvm::SuSp((2.0/3.0)*rhom*divUm, km) + - fvm::Sp(rhom*epsilonm/km, km) + + kSource() + ); + + kmEqn().relax(); + solve(kmEqn); + bound(km, this->kMin_); + km.correctBoundaryConditions(); + + volScalarField Cc2(rhom/(alphal*rholEff() + alphag*rhogEff()*Ct2_())); + kl = Cc2*km; + kl.correctBoundaryConditions(); + epsilonl = Cc2*epsilonm; + epsilonl.correctBoundaryConditions(); + liquidTurbulence.correctNut(); + + Ct2_() = Ct2(); + kg = Ct2_()*kl; + kg.correctBoundaryConditions(); + epsilong = Ct2_()*epsilonl; + epsilong.correctBoundaryConditions(); + nutg = Ct2_()*(liquidTurbulence.nu()/this->nu())*nutl; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H new file mode 100644 index 0000000000..ccc13b8035 --- /dev/null +++ b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H @@ -0,0 +1,252 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::RASModels::mixtureKEpsilon + +Group + grpRASTurbulence + +Description + Standard k-epsilon turbulence model + + The default model coefficients correspond to the following: + \verbatim + mixtureKEpsilonCoeffs + { + Cmu 0.09; + C1 1.44; + C2 1.92; + C3 -0.33; + sigmak 1.0; + sigmaEps 1.3; + } + \endverbatim + +SourceFiles + mixtureKEpsilon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef mixtureKEpsilon_H +#define mixtureKEpsilon_H + +#include "RASModel.H" +#include "eddyViscosity.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class mixtureKEpsilon Declaration +\*---------------------------------------------------------------------------*/ + +template +class mixtureKEpsilon +: + public eddyViscosity > +{ + // Private data + + mutable mixtureKEpsilon *liquidTurbulencePtr_; + + + // Private Member Functions + + // Disallow default bitwise copy construct and assignment + mixtureKEpsilon(const mixtureKEpsilon&); + mixtureKEpsilon& operator=(const mixtureKEpsilon&); + + //- Return the turbulence model for the other phase + mixtureKEpsilon& liquidTurbulence() const; + + +protected: + + // Protected data + + // Model coefficients + + dimensionedScalar Cmu_; + dimensionedScalar C1_; + dimensionedScalar C2_; + dimensionedScalar C3_; + dimensionedScalar Cp_; + dimensionedScalar sigmak_; + dimensionedScalar sigmaEps_; + + // Fields + + volScalarField k_; + volScalarField epsilon_; + + // Mixture fields + + autoPtr Ct2_; + autoPtr rhom_; + autoPtr km_; + autoPtr epsilonm_; + + + // Protected Member Functions + + wordList epsilonBoundaryTypes(const volScalarField& epsilon) const; + + void correctInletOutlet + ( + volScalarField& vsf, + const volScalarField& refVsf + ) const; + + void initMixtureFields(); + + virtual void correctNut(); + + tmp Ct2() const; + + tmp rholEff() const; + tmp rhogEff() const; + tmp rhom() const; + + tmp mix + ( + const volScalarField& fc, + const volScalarField& fd + ) const; + + tmp mixU + ( + const volScalarField& fc, + const volScalarField& fd + ) const; + + tmp mixFlux + ( + const surfaceScalarField& fc, + const surfaceScalarField& fd + ) const; + + tmp bubbleG() const; + virtual tmp kSource() const; + virtual tmp epsilonSource() const; + + //- Return the effective diffusivity for k + tmp DkEff(const volScalarField& nutm) const + { + return tmp + ( + new volScalarField + ( + "DkEff", + nutm/sigmak_ + ) + ); + } + + //- Return the effective diffusivity for epsilon + tmp DepsilonEff(const volScalarField& nutm) const + { + return tmp + ( + new volScalarField + ( + "DepsilonEff", + nutm/sigmaEps_ + ) + ); + } + +public: + + typedef typename BasicTurbulenceModel::alphaField alphaField; + typedef typename BasicTurbulenceModel::rhoField rhoField; + typedef typename BasicTurbulenceModel::transportModel transportModel; + + + //- Runtime type information + TypeName("mixtureKEpsilon"); + + + // Constructors + + //- Construct from components + mixtureKEpsilon + ( + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName = turbulenceModel::propertiesName, + const word& type = typeName + ); + + + //- Destructor + virtual ~mixtureKEpsilon() + {} + + + // Member Functions + + //- Re-read model coefficients if they have changed + virtual bool read(); + + //- Return the turbulence kinetic energy + virtual tmp k() const + { + return k_; + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp epsilon() const + { + return epsilon_; + } + + //- Solve the turbulence equations and correct the turbulence viscosity + virtual void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "mixtureKEpsilon.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C index cecceba932..15543f1851 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C @@ -205,7 +205,7 @@ void Foam::patchInjectionBase::setPositionAndCell const point pf(tri.randomPoint(rnd)); // position perturbed away from face (into domain) - const scalar a = rnd.position(0.1, 0.5); + const scalar a = rnd.position(scalar(0.1), scalar(0.5)); const vector& pc = mesh.cellCentres()[cellOwner]; const vector d = mag(pf - pc)*patchNormal_[faceI]; diff --git a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C index b4d433d425..7237f30307 100644 --- a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C @@ -197,7 +197,7 @@ updateCoeffs() const vector& myRayId = dom.IRay(rayId).d(); - const scalarField& Ir = dom.Qin(); + const scalarField& Ir = dom.Qin().boundaryField()[patchI]; forAll(Iw, faceI) { diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.air b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.air index 40de4981f5..fe69aba1ea 100644 --- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.air +++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.air @@ -19,10 +19,16 @@ simulationType RAS; RAS { - RASModel continuousGasKEpsilon; + RASModel mixtureKEpsilon; //continuousGasKEpsilon; turbulence on; printCoeffs on; + + // mixtureKEpsilonCoeffs + // { + // Cp 1; + // C3 1; + // } } // ************************************************************************* // diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.water b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.water index 2c0b129e11..b82e81c0cb 100644 --- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.water +++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/constant/turbulenceProperties.water @@ -19,7 +19,7 @@ simulationType RAS; RAS { - RASModel LaheyKEpsilon; + RASModel mixtureKEpsilon; //LaheyKEpsilon; turbulence on; printCoeffs on; diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes index f55c52ee56..d3f22a16e8 100644 --- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes +++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes @@ -40,6 +40,7 @@ divSchemes "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1; "div\(alphaPhi.*,(k|epsilon).*\)" Gauss limitedLinear 1; + "div\(phim,(k|epsilon)m\)" Gauss limitedLinear 1; "div\(\(\(alpha.*nuEff.*\)*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear; }