From a7eb3505362003826f3f31428acf59d755446003 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 11 Mar 2020 23:24:08 +0000 Subject: [PATCH] turbulenceModels/laminar: Maxwell, Giesekus: Added multi-mode support MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit By specifying a list of coefficients in turbulenceProperties, e.g. for the generalised Maxwell model: modes ( { lambda 0.01; } { lambda 0.04; } ); of for the generalised Giesekus model: modes ( { lambda 0.01; alphaG 0.05; } { lambda 0.04; alphaG 0.2; } ); Visco-elasticity stress tensors (sigma0, sigma1...) are solved for each mode and summed to create the effective stress of the complex fluid: Any number of modes can be specified and if only one mode is required the 'modes' entry is not read and the coefficients are obtained as before. The mode sigma? fields are read if present otherwise are constructed and initialised from the sigma field but all of the mode sigma? fields are written for restart and the sigma field contains the sum. References: http://en.wikipedia.org/wiki/Generalized_Maxwell_model Wiechert, E. (1889). Ueber elastische Nachwirkung. (Doctoral dissertation, Hartungsche buchdr.). Wiechert, E. (1893). Gesetze der elastischen Nachwirkung für constante Temperatur. Annalen der Physik, 286(11), 546-570. --- .../laminar/Giesekus/Giesekus.C | 14 +- .../laminar/Giesekus/Giesekus.H | 15 +- .../laminar/Maxwell/Maxwell.C | 248 ++++++++++++++---- .../laminar/Maxwell/Maxwell.H | 38 ++- .../constant/turbulenceProperties | 30 +++ 5 files changed, 285 insertions(+), 60 deletions(-) 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;