View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002292OpenFOAMBugpublic2016-10-14 14:312016-10-20 09:11
Reporterqiqi Assigned Tohenry  
PriorityhighSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version15.04
Fixed in Versiondev 
Summary0002292: SpalartAllmaras in pisoFoam not restarting correctly
Descriptionturbulence->validate(); in pisoFoam.C changes nut, thereby making the flow different from the case without save-and-restart.
Steps To ReproduceRun a simulation for 10 steps, save it, and continue for another 10 steps. Compare the outcome of the second 10 steps with the outcome of running 20 steps from the beginning. The difference will be in the 5th digit. But for chaotic turbulent flows, this difference will amplify over time, making simulations with and without restart totally different.
TagsNo tags attached.

Activities

henry

2016-10-14 14:42

manager   ~0007008

Chaotic turbulent flows are indeed chaotic and only the averaged statistics are reproducible. If you average the results for a long time, restart and restart the averaging do the averaged statistics agree between the two runs?

qiqi

2016-10-14 14:49

reporter   ~0007009

That's correct. I found this bug, which does not occur for other LES models, when experimenting with numerical methods for sensitivity analysis. Being able to restart and continue on exactly the same trajectory is important for me, and likely for other numerical methods, e.g., for computing Lyapunov exponents of the flow dynamics.

henry

2016-10-14 15:00

manager   ~0007010

All turbulence->validate() does is call correctNut which updates the nut field in the same manner as it does after the calculation of nuTilda. This call is necessary to ensure the nut field is physical and consistent with nuTilda at the beginning of the run, particularly when starting for user-input fields in which nut may be specified as 0. There is no difference in this procedure between SpalartAllmaras and any of the other models.

qiqi

2016-10-15 02:23

reporter   ~0007011

You are right. turbulence->validate() should not have changed the nut that was computed in the last time step. But somehow it does when using SpalartAllmarasDES.

henry

2016-10-15 07:29

manager   ~0007012

Here is the correctNut code from the SpalartAllmarasDES model:

template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut
(
    const volScalarField& fv1
)
{
    this->nut_ = nuTilda_*fv1;
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
{
    correctNut(fv1(this->chi()));
}

I don't see anything wrong with this.

henry

2016-10-15 07:44

manager   ~0007013

Incidentally SpalartAllmarasDDES and SpalartAllmarasIDDES are derived from SpalartAllmarasDES and inherit the same correctNut code.

wyldckat

2016-10-18 13:47

updater   ~0007027

I suspect that the original diagnosis isn't fully correct. The re-calculation of the "nut" field is very likely only a symptom of the real problem.
Furthermore, the original report refers to "SpalartAllmaras" and not "SpalartAllmarasDES", which increases my suspicion.

@qiqi: Can you please provide a test case that replicates the problem that you're seeing?
I ask this because the most likely suspect is actually the boundary conditions at the inlet(s), given that if it's using a random inlet source for flow speed and/or pressure values, that would justify the problem.

If this is the case, either the inlet values need to be pre-planned for each time step or a re-programmable random number generator would have to be implemented. The random generators in C++11 do have this feature, namely to get a seed for at a particular time step, to be re-use in the next run.

qiqi

2016-10-18 15:11

reporter   ~0007028

The testcase is encapsulated in

https://github.com/qiqi/fds/blob/master/tests/test_pisofoam4_restart_io.py

The case is in https://github.com/qiqi/fds/tree/master/tests/data/pisofoam_restart

Note that all the boundary conditions are "clean" -- fixedValue, zeroGradient.

The test passes only with a modified version of pisoFoam, the only difference between which and the shipped version is the validate(); line.

wyldckat

2016-10-19 00:27

updater   ~0007031

Many thanks for the test case! With it I've been able to reproduce the same problem, but I haven't fully figured out yet the real source of the problem. The narrowest I've gotten when trailing back was that this expression seems to be at the genesis of the problem:

   this->nut_ = nuTilda_*fv1;


I was able to do a quick test after running pisoFoam once for 1s and then using the following hack in pisoFoam:

      turbulence->validate();
      runTime++;
      runTime.writeNow();
      return 0;

The "nut" field was considerably different, although still in similar ranges of values.


OK, I believe I've got it. I see what's wrong. In "SpalartAllmarasDES<BasicTurbulenceModel>::correct()", "fv1" is calculated with the "nuTilda" from before "nuTildaEqn" was solved. "correctNut(fv1)" at the end of the method uses the outdated fv1 value, resulting in the discrepancy. Removing "fv1" from the call should re-establish consistency.

Attached is the file "SpalartAllmarasDES.C" that replaces "src/TurbulenceModels/turbulenceModels/LES/SpalartAllmarasDES/SpalartAllmarasDES.C".

@qiqi: I haven't enough time today to test this, so any chance you can test this fix? It should work as intended, but it's better to make sure ;)

wyldckat

2016-10-19 00:28

updater  

SpalartAllmarasDES.C (10,873 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "SpalartAllmarasDES.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace LESModels
{

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const
{
    return nuTilda_/this->nu();
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1
(
    const volScalarField& chi
) const
{
    const volScalarField chi3("chi3", pow3(chi));
    return chi3/(chi3 + pow3(Cv1_));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2
(
    const volScalarField& chi,
    const volScalarField& fv1
) const
{
    return 1.0 - chi/(1.0 + chi*fv1);
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::S
(
    const volTensorField& gradU
) const
{
    return sqrt(2.0)*mag(symm(gradU));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega
(
    const volTensorField& gradU
) const
{
    return sqrt(2.0)*mag(skew(gradU));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
(
    const volScalarField& chi,
    const volScalarField& fv1,
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    return
    (
        max
        (
            Omega
          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
            Cs_*Omega
        )
    );
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r
(
    const volScalarField& nur,
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    tmp<volScalarField> tr
    (
        min
        (
            nur
           /(
                max
                (
                    Omega,
                    dimensionedScalar("SMALL", Omega.dimensions(), SMALL)
                )
                *sqr(kappa_*dTilda)
            ),
            scalar(10)
        )
    );
    tr.ref().boundaryFieldRef() == 0.0;

    return tr;
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw
(
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    const volScalarField r(this->r(nuTilda_, Omega, dTilda));
    const volScalarField g(r + Cw2_*(pow6(r) - r));

    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
(
    const volScalarField& chi,
    const volScalarField& fv1,
    const volTensorField& gradU
) const
{
    tmp<volScalarField> tdTilda(CDES_*this->delta());
    min(tdTilda.ref().ref(), tdTilda(), y_);
    return tdTilda;
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut
(
    const volScalarField& fv1
)
{
    this->nut_ = nuTilda_*fv1;
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
{
    correctNut(fv1(this->chi()));
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    LESeddyViscosity<BasicTurbulenceModel>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    sigmaNut_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmaNut",
            this->coeffDict_,
            0.66666
        )
    ),
    kappa_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "kappa",
            this->coeffDict_,
            0.41
        )
    ),
    Cb1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb1",
            this->coeffDict_,
            0.1355
        )
    ),
    Cb2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb2",
            this->coeffDict_,
            0.622
        )
    ),
    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
    Cw2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw2",
            this->coeffDict_,
            0.3
        )
    ),
    Cw3_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw3",
            this->coeffDict_,
            2.0
        )
    ),
    Cv1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cv1",
            this->coeffDict_,
            7.1
        )
    ),
    Cs_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cs",
            this->coeffDict_,
            0.3
        )
    ),
    CDES_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "CDES",
            this->coeffDict_,
            0.65
        )
    ),
    ck_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "ck",
            this->coeffDict_,
            0.07
        )
    ),

    nuTilda_
    (
        IOobject
        (
            "nuTilda",
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh_
    ),

    y_(wallDist::New(this->mesh_).y())
{
    if (type == typeName)
    {
        this->printCoeffs(type);
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
{
    if (LESeddyViscosity<BasicTurbulenceModel>::read())
    {
        sigmaNut_.readIfPresent(this->coeffDict());
        kappa_.readIfPresent(*this);

        Cb1_.readIfPresent(this->coeffDict());
        Cb2_.readIfPresent(this->coeffDict());
        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
        Cw2_.readIfPresent(this->coeffDict());
        Cw3_.readIfPresent(this->coeffDict());
        Cv1_.readIfPresent(this->coeffDict());
        Cs_.readIfPresent(this->coeffDict());

        CDES_.readIfPresent(this->coeffDict());
        ck_.readIfPresent(this->coeffDict());

        return true;
    }
    else
    {
        return false;
    }
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::
DnuTildaEff() const
{
    return tmp<volScalarField>
    (
        new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
    );
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::k() const
{
    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));
    return sqr(this->nut()/ck_/dTilda(chi, fv1, fvc::grad(this->U_)));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
{
    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));

    tmp<volScalarField> tLESRegion
    (
        new volScalarField
        (
            IOobject
            (
                "DES::LESRegion",
                this->mesh_.time().timeName(),
                this->mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            neg(dTilda(chi, fv1, fvc::grad(this->U_)) - y_)
        )
    );

    return tLESRegion;
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    LESeddyViscosity<BasicTurbulenceModel>::correct();

    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));

    tmp<volTensorField> tgradU = fvc::grad(U);
    const volScalarField Omega(this->Omega(tgradU()));
    const volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
    const volScalarField Stilda(this->Stilda(chi, fv1, Omega, dTilda));

    tmp<fvScalarMatrix> nuTildaEqn
    (
        fvm::ddt(alpha, rho, nuTilda_)
      + fvm::div(alphaRhoPhi, nuTilda_)
      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
     ==
        Cb1_*alpha*rho*Stilda*nuTilda_
      - fvm::Sp
        (
            Cw1_*alpha*rho*fw(Stilda, dTilda)*nuTilda_/sqr(dTilda),
            nuTilda_
        )
      + fvOptions(alpha, rho, nuTilda_)
    );

    nuTildaEqn.ref().relax();
    fvOptions.constrain(nuTildaEqn.ref());
    solve(nuTildaEqn);
    fvOptions.correct(nuTilda_);
    bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
    nuTilda_.correctBoundaryConditions();

    correctNut();
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace LESModels
} // End namespace Foam

// ************************************************************************* //
SpalartAllmarasDES.C (10,873 bytes)   

wyldckat

2016-10-19 00:28

updater  

SpalartAllmarasDES-2.C (10,873 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "SpalartAllmarasDES.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace LESModels
{

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const
{
    return nuTilda_/this->nu();
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1
(
    const volScalarField& chi
) const
{
    const volScalarField chi3("chi3", pow3(chi));
    return chi3/(chi3 + pow3(Cv1_));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2
(
    const volScalarField& chi,
    const volScalarField& fv1
) const
{
    return 1.0 - chi/(1.0 + chi*fv1);
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::S
(
    const volTensorField& gradU
) const
{
    return sqrt(2.0)*mag(symm(gradU));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega
(
    const volTensorField& gradU
) const
{
    return sqrt(2.0)*mag(skew(gradU));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
(
    const volScalarField& chi,
    const volScalarField& fv1,
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    return
    (
        max
        (
            Omega
          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
            Cs_*Omega
        )
    );
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r
(
    const volScalarField& nur,
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    tmp<volScalarField> tr
    (
        min
        (
            nur
           /(
                max
                (
                    Omega,
                    dimensionedScalar("SMALL", Omega.dimensions(), SMALL)
                )
                *sqr(kappa_*dTilda)
            ),
            scalar(10)
        )
    );
    tr.ref().boundaryFieldRef() == 0.0;

    return tr;
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw
(
    const volScalarField& Omega,
    const volScalarField& dTilda
) const
{
    const volScalarField r(this->r(nuTilda_, Omega, dTilda));
    const volScalarField g(r + Cw2_*(pow6(r) - r));

    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
(
    const volScalarField& chi,
    const volScalarField& fv1,
    const volTensorField& gradU
) const
{
    tmp<volScalarField> tdTilda(CDES_*this->delta());
    min(tdTilda.ref().ref(), tdTilda(), y_);
    return tdTilda;
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut
(
    const volScalarField& fv1
)
{
    this->nut_ = nuTilda_*fv1;
    this->nut_.correctBoundaryConditions();
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
{
    correctNut(fv1(this->chi()));
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    LESeddyViscosity<BasicTurbulenceModel>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    sigmaNut_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmaNut",
            this->coeffDict_,
            0.66666
        )
    ),
    kappa_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "kappa",
            this->coeffDict_,
            0.41
        )
    ),
    Cb1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb1",
            this->coeffDict_,
            0.1355
        )
    ),
    Cb2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cb2",
            this->coeffDict_,
            0.622
        )
    ),
    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
    Cw2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw2",
            this->coeffDict_,
            0.3
        )
    ),
    Cw3_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cw3",
            this->coeffDict_,
            2.0
        )
    ),
    Cv1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cv1",
            this->coeffDict_,
            7.1
        )
    ),
    Cs_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cs",
            this->coeffDict_,
            0.3
        )
    ),
    CDES_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "CDES",
            this->coeffDict_,
            0.65
        )
    ),
    ck_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "ck",
            this->coeffDict_,
            0.07
        )
    ),

    nuTilda_
    (
        IOobject
        (
            "nuTilda",
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh_
    ),

    y_(wallDist::New(this->mesh_).y())
{
    if (type == typeName)
    {
        this->printCoeffs(type);
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
{
    if (LESeddyViscosity<BasicTurbulenceModel>::read())
    {
        sigmaNut_.readIfPresent(this->coeffDict());
        kappa_.readIfPresent(*this);

        Cb1_.readIfPresent(this->coeffDict());
        Cb2_.readIfPresent(this->coeffDict());
        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
        Cw2_.readIfPresent(this->coeffDict());
        Cw3_.readIfPresent(this->coeffDict());
        Cv1_.readIfPresent(this->coeffDict());
        Cs_.readIfPresent(this->coeffDict());

        CDES_.readIfPresent(this->coeffDict());
        ck_.readIfPresent(this->coeffDict());

        return true;
    }
    else
    {
        return false;
    }
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::
DnuTildaEff() const
{
    return tmp<volScalarField>
    (
        new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
    );
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::k() const
{
    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));
    return sqr(this->nut()/ck_/dTilda(chi, fv1, fvc::grad(this->U_)));
}


template<class BasicTurbulenceModel>
tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
{
    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));

    tmp<volScalarField> tLESRegion
    (
        new volScalarField
        (
            IOobject
            (
                "DES::LESRegion",
                this->mesh_.time().timeName(),
                this->mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            neg(dTilda(chi, fv1, fvc::grad(this->U_)) - y_)
        )
    );

    return tLESRegion;
}


template<class BasicTurbulenceModel>
void SpalartAllmarasDES<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    LESeddyViscosity<BasicTurbulenceModel>::correct();

    const volScalarField chi(this->chi());
    const volScalarField fv1(this->fv1(chi));

    tmp<volTensorField> tgradU = fvc::grad(U);
    const volScalarField Omega(this->Omega(tgradU()));
    const volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
    const volScalarField Stilda(this->Stilda(chi, fv1, Omega, dTilda));

    tmp<fvScalarMatrix> nuTildaEqn
    (
        fvm::ddt(alpha, rho, nuTilda_)
      + fvm::div(alphaRhoPhi, nuTilda_)
      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
     ==
        Cb1_*alpha*rho*Stilda*nuTilda_
      - fvm::Sp
        (
            Cw1_*alpha*rho*fw(Stilda, dTilda)*nuTilda_/sqr(dTilda),
            nuTilda_
        )
      + fvOptions(alpha, rho, nuTilda_)
    );

    nuTildaEqn.ref().relax();
    fvOptions.constrain(nuTildaEqn.ref());
    solve(nuTildaEqn);
    fvOptions.correct(nuTilda_);
    bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
    nuTilda_.correctBoundaryConditions();

    correctNut();
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace LESModels
} // End namespace Foam

// ************************************************************************* //
SpalartAllmarasDES-2.C (10,873 bytes)   

wyldckat

2016-10-19 00:30

updater   ~0007032

OK... I'm not use to yet to the possibility to attach files along with the comment... please ignore the file "SpalartAllmarasDES-2.C", which is a duplicate of the previous upload "SpalartAllmarasDES.C".

wyldckat

2016-10-19 10:47

updater   ~0007036

@qiqi: I forgot to mention that after replacing the file, please run the "Allwmake" script in the folder "src/TurbulenceModels", because it needs to update at least a couple of libraries that depend on this template class.

wyldckat

2016-10-19 23:20

updater   ~0007049

@Henry: I've confirmed that the provided test case works as intended when using the attached file "SpalartAllmarasDES.C", for replacing "src/TurbulenceModels/turbulenceModels/LES/SpalartAllmarasDES/SpalartAllmarasDES.C"

The change was simply this:

  - correctNut(fv1);
  + correctNut();

so that the final calculated "nut" field within the method "correct()" is based on the latest "nuTilda" and not based on the one from before solving the equation.


I've also done a quick look through the other sibling template classes within LES and the only ones that use "correctNut()" with additional arguments are "dynamicKEqn" and "dynamicLagrangian", but those two rely on calculations done based on "U", which doesn't change within the "correct()" method.

qiqi

2016-10-20 02:10

reporter   ~0007051

Thanks for fixing this issue. Great job Henry.

qiqi

2016-10-20 02:15

reporter   ~0007052

I should be thanking @wyldckat. Thanks wyldckat. Great job.

henry

2016-10-20 09:11

manager   ~0007056

Resolved in OpenFOAM-dev by commit 667dc10af61341690286ad2420b996937f94464f

Issue History

Date Modified Username Field Change
2016-10-14 14:31 qiqi New Issue
2016-10-14 14:42 henry Note Added: 0007008
2016-10-14 14:49 qiqi Note Added: 0007009
2016-10-14 15:00 henry Note Added: 0007010
2016-10-15 02:23 qiqi Note Added: 0007011
2016-10-15 07:29 henry Note Added: 0007012
2016-10-15 07:44 henry Note Added: 0007013
2016-10-18 13:47 wyldckat Note Added: 0007027
2016-10-18 15:11 qiqi Note Added: 0007028
2016-10-19 00:27 wyldckat Note Added: 0007031
2016-10-19 00:28 wyldckat File Added: SpalartAllmarasDES.C
2016-10-19 00:28 wyldckat File Added: SpalartAllmarasDES-2.C
2016-10-19 00:30 wyldckat Note Added: 0007032
2016-10-19 10:47 wyldckat Note Added: 0007036
2016-10-19 23:13 wyldckat Assigned To => henry
2016-10-19 23:13 wyldckat Status new => assigned
2016-10-19 23:20 wyldckat Note Added: 0007049
2016-10-20 02:10 qiqi Note Added: 0007051
2016-10-20 02:15 qiqi Note Added: 0007052
2016-10-20 09:11 henry Status assigned => resolved
2016-10-20 09:11 henry Resolution open => fixed
2016-10-20 09:11 henry Fixed in Version => dev
2016-10-20 09:11 henry Note Added: 0007056