turbulenceModels/laminar: Maxwell, Giesekus: Added multi-mode support

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.

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<BasicTurbulenceModel>::Giesekus
         type
     ),
 
-    alphaG_("alphaG", dimless, this->coeffDict_)
+    alphaGs_(this->readModeCoefficients("alphaG", dimless))
 {
     if (type == typeName)
     {
@@ -76,7 +76,7 @@ bool Giesekus<BasicTurbulenceModel>::read()
 {
     if (Maxwell<BasicTurbulenceModel>::read())
     {
-        alphaG_.read(this->coeffDict());
+        alphaGs_ = this->readModeCoefficients("alphaG", dimless);
 
         return true;
     }
@@ -89,12 +89,16 @@ bool Giesekus<BasicTurbulenceModel>::read()
 
 template<class BasicTurbulenceModel>
 tmp<fvSymmTensorMatrix>
-Giesekus<BasicTurbulenceModel>::sigmaSource() const
+Giesekus<BasicTurbulenceModel>::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
     );
 }
 

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<dimensionedScalar> alphaGs_;
 
 
     // Protected Member Functions
 
-        virtual tmp<fvSymmTensorMatrix> sigmaSource() const;
+        virtual tmp<fvSymmTensorMatrix> sigmaSource
+        (
+            const label modei,
+            volSymmTensorField& sigma
+        ) const;
 
 
 public:

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<class BasicTurbulenceModel>
-tmp<fvSymmTensorMatrix> Maxwell<BasicTurbulenceModel>::sigmaSource() const
+Foam::PtrList<Foam::dimensionedScalar>
+Maxwell<BasicTurbulenceModel>::readModeCoefficients
+(
+    const word& name,
+    const dimensionSet& dims
+) const
+{
+    PtrList<dimensionedScalar> 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<scalar>
+                (
+                    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<scalar>
+            (
+                name,
+                dims,
+                this->coeffDict_.lookup(name)
+            )
+        );
+    }
+
+    return modeCoeffs;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvSymmTensorMatrix> Maxwell<BasicTurbulenceModel>::sigmaSource
+(
+    const label modei,
+    volSymmTensorField& sigma
+) const
 {
     return tmp<fvSymmTensorMatrix>
     (
         new fvSymmTensorMatrix
         (
-            sigma_,
+            sigma,
-            dimVolume*this->rho_.dimensions()*sigma_.dimensions()/dimTime
+            dimVolume*this->rho_.dimensions()*sigma.dimensions()/dimTime
         )
     );
 }
@@ -77,6 +138,18 @@ Maxwell<BasicTurbulenceModel>::Maxwell
         propertiesName
     ),
 
+    modeCoefficients_
+    (
+        this->coeffDict().found("modes")
+      ? PtrList<dictionary>
+        (
+            this->coeffDict().lookup("modes")
+        )
+      : PtrList<dictionary>()
+    ),
+
+    nModes_(modeCoefficients_.size() ? modeCoefficients_.size() : 1),
+
     nuM_
     (
         dimensioned<scalar>
@@ -87,15 +160,7 @@ Maxwell<BasicTurbulenceModel>::Maxwell
         )
     ),
 
-    lambda_
+    lambdas_(readModeCoefficients("lambda", dimTime)),
-    (
-        dimensioned<scalar>
-        (
-            "lambda",
-            dimTime,
-            this->coeffDict_.lookup("lambda")
-        )
-    ),
 
     sigma_
     (
@@ -110,6 +175,65 @@ Maxwell<BasicTurbulenceModel>::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<volSymmTensorField>(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<BasicTurbulenceModel>::read()
 {
     if (laminarModel<BasicTurbulenceModel>::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<BasicTurbulenceModel>::read()
     }
 }
 
+
 template<class BasicTurbulenceModel>
 tmp<Foam::volSymmTensorField>
 Maxwell<BasicTurbulenceModel>::R() const
@@ -142,6 +273,7 @@ Maxwell<BasicTurbulenceModel>::R() const
     return sigma_;
 }
 
+
 template<class BasicTurbulenceModel>
 tmp<Foam::volSymmTensorField>
 Maxwell<BasicTurbulenceModel>::devRhoReff() const
@@ -205,7 +337,6 @@ void Maxwell<BasicTurbulenceModel>::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<BasicTurbulenceModel>::correct();
@@ -213,40 +344,67 @@ void Maxwell<BasicTurbulenceModel>::correct()
     tmp<volTensorField> tgradU(fvc::grad(U));
     const volTensorField& gradU = tgradU();
 
-    uniformDimensionedScalarField rLambda
+    forAll(lambdas_, modei)
-    (
+    {
-        IOobject
+        volSymmTensorField& sigma = nModes_ == 1 ? sigma_ : sigmas_[modei];
+
+        uniformDimensionedScalarField rLambda
         (
-            IOobject::groupName("rLambda", this->alphaRhoPhi_.group()),
+            IOobject
-            this->runTime_.constant(),
+            (
-            this->mesh_
+                IOobject::groupName
-        ),
+                (
-        1.0/(lambda_)
+                    "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
+        // Note sigma is positive on lhs of momentum eqn
-    volSymmTensorField P
+        const volSymmTensorField P
-    (
+        (
-        twoSymm(sigma & gradU)
+            twoSymm(sigma & gradU)
-      - nuM_*rLambda*twoSymm(gradU)
+          - nuM_*rLambda*twoSymm(gradU)
-    );
+        );
 
-    // Viscoelastic stress equation
+        // Viscoelastic stress equation
-    tmp<fvSymmTensorMatrix> sigmaEqn
+        fvSymmTensorMatrix sigmaEqn
-    (
+        (
-        fvm::ddt(alpha, rho, sigma)
+            fvm::ddt(alpha, rho, sigma)
-      + fvm::div(alphaRhoPhi, sigma)
+          + fvm::div
-      + fvm::Sp(alpha*rho*rLambda, sigma)
+            (
-     ==
+                alphaRhoPhi,
-        alpha*rho*P
+                sigma,
-      + sigmaSource()
+                "div(" + alphaRhoPhi.name() + ',' + sigma_.name() + ')'
-      + fvOptions(alpha, rho, sigma)
+            )
-    );
+          + fvm::Sp(alpha*rho*rLambda, sigma)
+         ==
+            alpha*rho*P
+          + sigmaSource(modei, sigma)
+          + fvOptions(alpha, rho, sigma)
+        );
 
-    sigmaEqn.ref().relax();
+        sigmaEqn.relax();
-    fvOptions.constrain(sigmaEqn.ref());
+        fvOptions.constrain(sigmaEqn);
-    solve(sigmaEqn);
+        sigmaEqn.solve("sigma");
-    fvOptions.correct(sigma_);
+        fvOptions.correct(sigma);
+    }
+
+    if (sigmas_.size())
+    {
+        volSymmTensorField sigmaSum("sigmaSum", sigmas_[0]);
+
+        for (label modei = 1; modei<sigmas_.size(); modei++)
+        {
+            sigmaSum += sigmas_[modei];
+        }
+
+        sigma_ == sigmaSum;
+    }
 }
 
 

src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H

@@ -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
@@ -25,9 +25,11 @@ Class
     Foam::laminarModels::Maxwell
 
 Description
-    Maxwell model for viscoelasticity using the upper-convected time
+    Generalised Maxwell model for viscoelasticity using the upper-convected time
-    derivative of the stress tensor.
+    derivative of the stress tensor with support for multiple modes.
+
     See http://en.wikipedia.org/wiki/Upper-convected_Maxwell_model
+        http://en.wikipedia.org/wiki/Generalized_Maxwell_model
 
     The model includes an additional viscosity (nu) from the transport
     model from which it is instantiated, which makes it equivalent to
@@ -37,6 +39,13 @@ Description
 
     Reference:
     \verbatim
+        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.
+
         Amoreira, L. J., & Oliveira, P. J. (2010).
         Comparison of different formulations for the numerical calculation
         of unsteady incompressible viscoelastic fluid flow.
@@ -76,24 +85,43 @@ protected:
 
         // Model coefficients
 
+            PtrList<dictionary> modeCoefficients_;
+
+            label nModes_;
+
             dimensionedScalar nuM_;
-            dimensionedScalar lambda_;
+
+            PtrList<dimensionedScalar> lambdas_;
 
 
         // Fields
 
+            //- Single or mode sum viscoelastic stress
             volSymmTensorField sigma_;
 
+            //- Mode viscoelastic stresses
+            PtrList<volSymmTensorField> sigmas_;
+
 
     // Protected Member Functions
 
+        PtrList<dimensionedScalar> readModeCoefficients
+        (
+            const word& name,
+            const dimensionSet& dims
+        ) const;
+
         //- Return the turbulence viscosity
         tmp<volScalarField> nu0() const
         {
             return this->nu() + nuM_;
         }
 
-        virtual tmp<fvSymmTensorMatrix> sigmaSource() const;
+        virtual tmp<fvSymmTensorMatrix> sigmaSource
+        (
+            const label modei,
+            volSymmTensorField& sigma
+        ) const;
 
 
 public:

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;