diff --git a/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.C b/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.C index d15a8cf332..3c97e6eb03 100644 --- a/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.C +++ b/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -60,7 +60,7 @@ Giesekus::Giesekus type ), - alphaG_("alphaG", dimless, this->coeffDict_) + alphaGs_(this->readModeCoefficients("alphaG", dimless)) { if (type == typeName) { @@ -76,7 +76,7 @@ bool Giesekus::read() { if (Maxwell::read()) { - alphaG_.read(this->coeffDict()); + alphaGs_ = this->readModeCoefficients("alphaG", dimless); return true; } @@ -89,12 +89,16 @@ bool Giesekus::read() template tmp -Giesekus::sigmaSource() const +Giesekus::sigmaSource +( + const label modei, + volSymmTensorField& sigma +) const { return fvm::Su ( this->alpha_*this->rho_ - *alphaG_*innerSqr(this->sigma_)/this->nuM_, this->sigma_ + *alphaGs_[modei]*innerSqr(sigma)/this->nuM_, sigma ); } diff --git a/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.H b/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.H index a56364d756..5fd1d2807f 100644 --- a/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.H +++ b/src/TurbulenceModels/turbulenceModels/laminar/Giesekus/Giesekus.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,8 @@ Class Foam::laminarModels::Giesekus Description - Giesekus model for viscoelasticity. + Giesekus model for viscoelasticity using the upper-convected time + derivative of the stress tensor with support for multiple modes. Reference: \verbatim @@ -68,14 +69,18 @@ protected: // Protected data - // Model coefficients + // Mode coefficients - dimensionedScalar alphaG_; + PtrList alphaGs_; // Protected Member Functions - virtual tmp sigmaSource() const; + virtual tmp sigmaSource + ( + const label modei, + volSymmTensorField& sigma + ) const; public: diff --git a/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C b/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C index e667f43a9f..905c440c27 100644 --- a/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C +++ b/src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,14 +37,75 @@ namespace laminarModels // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // template -tmp Maxwell::sigmaSource() const +Foam::PtrList +Maxwell::readModeCoefficients +( + const word& name, + const dimensionSet& dims +) const +{ + PtrList modeCoeffs(nModes_); + + if (nModes_ > 1) + { + if (this->coeffDict().found(name)) + { + IOWarningInFunction(this->coeffDict()) + << "For multi-mode entry '" << name << "' will be ignored." + << endl; + } + + forAll(modeCoefficients_, modei) + { + modeCoeffs.set + ( + modei, + new dimensioned + ( + name, + dims, + modeCoefficients_[modei].lookup(name) + ) + ); + } + } + else + { + if (modeCoefficients_.size() == 1) + { + IOWarningInFunction(this->coeffDict()) + << "For single mode 'modes' entry will be ignored." << endl; + } + + modeCoeffs.set + ( + 0, + new dimensioned + ( + name, + dims, + this->coeffDict_.lookup(name) + ) + ); + } + + return modeCoeffs; +} + + +template +tmp Maxwell::sigmaSource +( + const label modei, + volSymmTensorField& sigma +) const { return tmp ( new fvSymmTensorMatrix ( - sigma_, - dimVolume*this->rho_.dimensions()*sigma_.dimensions()/dimTime + sigma, + dimVolume*this->rho_.dimensions()*sigma.dimensions()/dimTime ) ); } @@ -77,6 +138,18 @@ Maxwell::Maxwell propertiesName ), + modeCoefficients_ + ( + this->coeffDict().found("modes") + ? PtrList + ( + this->coeffDict().lookup("modes") + ) + : PtrList() + ), + + nModes_(modeCoefficients_.size() ? modeCoefficients_.size() : 1), + nuM_ ( dimensioned @@ -87,15 +160,7 @@ Maxwell::Maxwell ) ), - lambda_ - ( - dimensioned - ( - "lambda", - dimTime, - this->coeffDict_.lookup("lambda") - ) - ), + lambdas_(readModeCoefficients("lambda", dimTime)), sigma_ ( @@ -110,6 +175,65 @@ Maxwell::Maxwell this->mesh_ ) { + if (nModes_ > 1) + { + sigmas_.setSize(nModes_); + + forAll(sigmas_, modei) + { + IOobject header + ( + IOobject::groupName("sigma" + name(modei), alphaRhoPhi.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::NO_READ + ); + + // Check if mode field exists and can be read + if (header.typeHeaderOk(true)) + { + Info<< " Reading mode stress field " + << header.name() << endl; + + sigmas_.set + ( + modei, + new volSymmTensorField + ( + IOobject + ( + header.name(), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ) + ); + } + else + { + sigmas_.set + ( + modei, + new volSymmTensorField + ( + IOobject + ( + header.name(), + this->runTime_.timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sigma_ + ) + ); + } + } + } + if (type == typeName) { this->printCoeffs(type); @@ -124,8 +248,14 @@ bool Maxwell::read() { if (laminarModel::read()) { + if (modeCoefficients_.size()) + { + this->coeffDict().lookup("modes") >> modeCoefficients_; + } + nuM_.readIfPresent(this->coeffDict()); - lambda_.readIfPresent(this->coeffDict()); + + lambdas_ = readModeCoefficients("lambda", dimTime); return true; } @@ -135,6 +265,7 @@ bool Maxwell::read() } } + template tmp Maxwell::R() const @@ -142,6 +273,7 @@ Maxwell::R() const return sigma_; } + template tmp Maxwell::devRhoReff() const @@ -205,7 +337,6 @@ void Maxwell::correct() const rhoField& rho = this->rho_; const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_; const volVectorField& U = this->U_; - volSymmTensorField& sigma = this->sigma_; fv::options& fvOptions(fv::options::New(this->mesh_)); laminarModel::correct(); @@ -213,40 +344,67 @@ void Maxwell::correct() tmp tgradU(fvc::grad(U)); const volTensorField& gradU = tgradU(); - uniformDimensionedScalarField rLambda - ( - IOobject + forAll(lambdas_, modei) + { + volSymmTensorField& sigma = nModes_ == 1 ? sigma_ : sigmas_[modei]; + + uniformDimensionedScalarField rLambda ( - IOobject::groupName("rLambda", this->alphaRhoPhi_.group()), - this->runTime_.constant(), - this->mesh_ - ), - 1.0/(lambda_) - ); + IOobject + ( + IOobject::groupName + ( + "rLambda" + + (nModes_ == 1 ? word::null : name(modei)), + this->alphaRhoPhi_.group() + ), + this->runTime_.constant(), + this->mesh_ + ), + 1/lambdas_[modei] + ); - // Note sigma is positive on lhs of momentum eqn - volSymmTensorField P - ( - twoSymm(sigma & gradU) - - nuM_*rLambda*twoSymm(gradU) - ); + // Note sigma is positive on lhs of momentum eqn + const volSymmTensorField P + ( + twoSymm(sigma & gradU) + - nuM_*rLambda*twoSymm(gradU) + ); - // Viscoelastic stress equation - tmp sigmaEqn - ( - fvm::ddt(alpha, rho, sigma) - + fvm::div(alphaRhoPhi, sigma) - + fvm::Sp(alpha*rho*rLambda, sigma) - == - alpha*rho*P - + sigmaSource() - + fvOptions(alpha, rho, sigma) - ); + // Viscoelastic stress equation + fvSymmTensorMatrix sigmaEqn + ( + fvm::ddt(alpha, rho, sigma) + + fvm::div + ( + alphaRhoPhi, + sigma, + "div(" + alphaRhoPhi.name() + ',' + sigma_.name() + ')' + ) + + fvm::Sp(alpha*rho*rLambda, sigma) + == + alpha*rho*P + + sigmaSource(modei, sigma) + + fvOptions(alpha, rho, sigma) + ); - sigmaEqn.ref().relax(); - fvOptions.constrain(sigmaEqn.ref()); - solve(sigmaEqn); - fvOptions.correct(sigma_); + sigmaEqn.relax(); + fvOptions.constrain(sigmaEqn); + sigmaEqn.solve("sigma"); + fvOptions.correct(sigma); + } + + if (sigmas_.size()) + { + volSymmTensorField sigmaSum("sigmaSum", sigmas_[0]); + + for (label modei = 1; modei modeCoefficients_; + + label nModes_; + dimensionedScalar nuM_; - dimensionedScalar lambda_; + + PtrList lambdas_; // Fields + //- Single or mode sum viscoelastic stress volSymmTensorField sigma_; + //- Mode viscoelastic stresses + PtrList sigmas_; + // Protected Member Functions + PtrList readModeCoefficients + ( + const word& name, + const dimensionSet& dims + ) const; + //- Return the turbulence viscosity tmp nu0() const { return this->nu() + nuM_; } - virtual tmp sigmaSource() const; + virtual tmp sigmaSource + ( + const label modei, + volSymmTensorField& sigma + ) const; public: diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarContraction/constant/turbulenceProperties b/tutorials/incompressible/pimpleFoam/laminar/planarContraction/constant/turbulenceProperties index 843b011bda..5d734b6821 100644 --- a/tutorials/incompressible/pimpleFoam/laminar/planarContraction/constant/turbulenceProperties +++ b/tutorials/incompressible/pimpleFoam/laminar/planarContraction/constant/turbulenceProperties @@ -23,14 +23,44 @@ laminar MaxwellCoeffs { nuM 0.002; + + // Single mode coefficient lambda 0.03; + + // Example 2-mode specification + // modes + // ( + // { + // lambda 0.01; + // } + + // { + // lambda 0.04; + // } + // ); } GiesekusCoeffs { nuM 0.002; + + // Single mode coefficients lambda 0.03; alphaG 0.1; + + // Example 2-mode specification + // modes + // ( + // { + // lambda 0.01; + // alphaG 0.05; + // } + + // { + // lambda 0.04; + // alphaG 0.2; + // } + // ); } printCoeffs on;