[nightdust]

Dear Foamers,

I am rather new to OpenFOAM (and C++) and set out to implement a new viscosity model: a piece-wise power law. Meaning that, depending on the strain rate, different power law parameters K, n will be chosen.

I have used the built-in power law as a starting point. I did some forum searches on the forum which informed my current code solution. At the moment, I still retained the original parameters k_ and n_ in my viscosity model, but later they will be supplanted by the scalarLists kList_ and nList_ .

I also created a small test case (square pipe), but if I try to run the simulation with my new own viscosity model, it crashes. (See attached pictures #1 and #2). I am able to use the original power law (that's why I retained the parameters k_ and n_) and print out my calculated values of nu, based on the local strain rate.

What makes me wonder is that the values of nu for the inlet, outlet an walls are still equal to the initial value I gave my volScalarField in the viscosity model. I understand that I only manipulated the values at the cell centres, and the values at the boundary (inlet, outlet, walls) are defined at the face centres. I thought that they were calculated based on the values of the cell centres.

Maybe there is a problem with a zero? I am not knowledgeable enough. But if I simulate the test case using the old power law calculation and print out my calculated new values, they seem to be calculated fine at every time step (see picture #3).


Since the error message is in German, I searched for its English equivalent. It should read: Segmentation error (core dumped).



What am I missing here? Thanks for your help!


My pieceWisePowerLaw.C file:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

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

#include "pieceWisePowerLaw.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace viscosityModels
{
    defineTypeNameAndDebug(pieceWisePowerLaw, 0);

    addToRunTimeSelectionTable
    (
        viscosityModel,
        pieceWisePowerLaw,
        dictionary
    );
}
}


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

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::pieceWisePowerLaw::calcNu() const
{
    // declare the scalars for K, n which will be used for the current cell
    scalar kToUse;
    scalar nToUse;

    volScalarField myNu
    (
        IOobject
        (
            "myNu",
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        U_.mesh(),
        dimensionedScalar
        (
            "myNuInit",
            dimensionSet(0,2,-1,0,0,0,0),
            10.0
        )
    );
    
    // loop over all cells
    forAll(myNu, cellI)
    {   
        // use the K and n parameters of the last interval by default
        kToUse = kList_[kList_.size()-1];
        nToUse = nList_[nList_.size()-1];       
       
        int j = 0;
        bool intervalFound = false;
        while(!intervalFound && j < kList_.size())
        {               
            // if the strain rate of the current cell centre is lower (or equal) than the kth threshold value,
            // choose K, n of the kth power law interval
            if(strainRate()()[cellI] <= gammaList_[j])
            {           
                kToUse = kList_[j];
                nToUse = nList_[j];
                // if the correct interval has been found, this variable will also cause the outer while loop
                //to terminate
                intervalFound = true;
            }                           

            j++;
        }

        // actual calculation of the viscosity (like power-law, but K, n are chosen based on the strain rate
        // of the current cell centre)
        myNu[cellI] = max
            (
                nuMin_.value(),
                min
                (
                    nuMax_.value(),               
                    kToUse*pow
                    (
                        max
                        (
                            strainRate()()[cellI],                   
                            scalar(VSMALL)
                        ),
                    nToUse - 1.0
                    )
                )
            );
    }
    
    Info << "myNu: " << myNu << nl;
    
    //return myNu;
    
    // the original return (power law)
    return max
    (
    nuMin_,
    min
    (
        nuMax_,
        k_*pow
        (
        max
        (
            dimensionedScalar("one", dimTime, 1.0)*strainRate(),
            dimensionedScalar("VSMALL", dimless, VSMALL)
        ),
        n_.value() - scalar(1.0)
        )
    )
    );
}


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

Foam::viscosityModels::pieceWisePowerLaw::pieceWisePowerLaw
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    viscosityModel(name, viscosityProperties, U, phi),
    pieceWisePowerLawCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")),
    k_(pieceWisePowerLawCoeffs_.lookup("k")),
    n_(pieceWisePowerLawCoeffs_.lookup("n")),
    nuMin_(pieceWisePowerLawCoeffs_.lookup("nuMin")),
    nuMax_(pieceWisePowerLawCoeffs_.lookup("nuMax")),
    gammaList_(pieceWisePowerLawCoeffs_.lookup("gammaList")),
    kList_(pieceWisePowerLawCoeffs_.lookup("kList")),
    nList_(pieceWisePowerLawCoeffs_.lookup("nList")),
    nu_
    (
        IOobject
        (
            name,
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        calcNu()
    )
{}


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

bool Foam::viscosityModels::pieceWisePowerLaw::read
(
    const dictionary& viscosityProperties
)
{
    viscosityModel::read(viscosityProperties);

    pieceWisePowerLawCoeffs_ = viscosityProperties.subDict(typeName + "Coeffs");

    pieceWisePowerLawCoeffs_.lookup("k") >> k_;
    pieceWisePowerLawCoeffs_.lookup("n") >> n_;
    pieceWisePowerLawCoeffs_.lookup("nuMin") >> nuMin_;
    pieceWisePowerLawCoeffs_.lookup("nuMax") >> nuMax_;
    pieceWisePowerLawCoeffs_.lookup("gammaList") >> gammaList_;
    pieceWisePowerLawCoeffs_.lookup("kList") >> kList_;
    pieceWisePowerLawCoeffs_.lookup("nList") >> nList_;

    return true;
}


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

My pieceWisePowerLaw.H file:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::viscosityModels::pieceWisePowerLaw

Description
     Standard power-law non-Newtonian viscosity model.

SourceFiles
    pieceWisePowerLaw.C

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

#ifndef pieceWisePowerLaw_H
#define pieceWisePowerLaw_H

#include "viscosityModel.H"
#include "dimensionedScalar.H"
#include "volFields.H"

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

namespace Foam
{
namespace viscosityModels
{

/*---------------------------------------------------------------------------*\
                           Class pieceWisePowerLaw Declaration
\*---------------------------------------------------------------------------*/

class pieceWisePowerLaw
:
    public viscosityModel
{
    // Private data

        dictionary pieceWisePowerLawCoeffs_;

    // these two variables will later be removed
        dimensionedScalar k_;
        dimensionedScalar n_;
        
    dimensionedScalar nuMin_;
        dimensionedScalar nuMax_;
    scalarList gammaList_;
    scalarList kList_;
    scalarList nList_;

        volScalarField nu_;


    // Private Member Functions

        //- Calculate and return the laminar viscosity
        tmp<volScalarField> calcNu() const;


public:

    //- Runtime type information
    TypeName("pieceWisePowerLaw");


    // Constructors

        //- Construct from components
        pieceWisePowerLaw
        (
            const word& name,
            const dictionary& viscosityProperties,
            const volVectorField& U,
            const surfaceScalarField& phi
        );


    // Destructor

        virtual ~pieceWisePowerLaw()
        {}


    // Member Functions

        //- Return the laminar viscosity
        virtual const volScalarField& nu() const
        {
            return nu_;
        }

        //- Correct the laminar viscosity
        virtual void correct()
        {
            nu_ = calcNu();
        }

        //- Read transportProperties dictionary
        virtual bool read(const dictionary& viscosityProperties);
};


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

} // End namespace viscosityModels
} // End namespace Foam

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

#endif

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