[Tibo99]

Hi everyone,

Since the flow I want to simulate show that dU/dz is non-linear, I want to implement an interpolation scheme on velocity U by applying the scheme on grad(U) in the sub dictionnary 'gradSchemes', file fvSchemes.

Everything compile successfully (no message of error), but the results (U, k ,epsilon) are not what I should observe (comparing with litterature).

Does anyone see something really bad that could cause an issue in the following code that I am endup with?

FYI, I used the interpolation scheme 'localMax' as a template.

Thank you very much for the feedback.

#ifndef UInterpolationScheme_H
#define UInterpolationScheme_H

#include "surfaceInterpolationScheme.H"
#include "volFields.H"
#include "surfaceFields.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                 Class UInterpolationScheme Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class UInterpolationScheme
:
    public surfaceInterpolationScheme<Type>
{
    // Private Data

    scalar z0_ = 0;

public:

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


    // Constructors

        //- Construct from Istream.
        UInterpolationScheme
        (
            const fvMesh& mesh,
            Istream& is
        )
        :
            surfaceInterpolationScheme<Type>(mesh)
        {
            z0_ = readScalar(is);
        }

        //- Construct from faceFlux and Istream
        UInterpolationScheme
        (
            const fvMesh& mesh,
            const surfaceScalarField& faceFlux,
            Istream& is
        )
        :
            surfaceInterpolationScheme<Type>(mesh)
        {}

    // Member Functions

        //- Return the interpolation weighting factors
        virtual tmp<surfaceScalarField> weights
        (
            const GeometricField<Type, fvPatchField, volMesh>&
        ) const
        {
            NotImplemented;

            return tmp<surfaceScalarField>(nullptr);
        }

        //- Return the face-interpolate of the given cell field
        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
        interpolate
        (
            const GeometricField<Type, fvPatchField, volMesh>& vf
        ) const
        {
            const fvMesh& mesh = vf.mesh();

            tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tvff
            (
                new GeometricField<Type, fvsPatchField, surfaceMesh>
                (
                    IOobject
                    (
                        "UInterpolationScheme::interpolate(" + vf.name() + ')',
                        mesh.time().timeName(),
                        mesh
                    ),
                    mesh,
                    vf.dimensions()
                )
            );
            GeometricField<Type, fvsPatchField, surfaceMesh>& vff = tvff.ref();

            // computing interpolation on the internal faces-------------------
            forAll(mesh.Cf(), facei)
            {
                // getting owner and neighbour cell indices first
                label owner_id = mesh.faceOwner()[facei];
                label neighbour_id = mesh.faceNeighbour()[facei];

                // finding which cell is to be P and N based on height
                scalar z_owner = mesh.C()[owner_id].z();
                scalar z_neighbour = mesh.C()[neighbour_id].z();

                label P_id, N_id;

                if (z_owner < z_neighbour)
                {
                    P_id = owner_id;
                    N_id = neighbour_id;
                }
                else
                {
                    N_id = owner_id;
                    P_id = neighbour_id;
                }

                // getting Z coordinates of P and N cells & n face
                scalar Z_N = mesh.C()[N_id].z()+z0_;
                scalar Z_P = mesh.C()[P_id].z()+z0_;
                scalar Z_n = mesh.Cf()[facei].z()+z0_;

                scalar numerator = Foam::log(Z_N/Z_n);
                scalar denominator = Foam::log(Z_N/Z_P);

                scalar alphaU(1.0);
                if(numerator != 0 || denominator != 0)
                    alphaU = numerator/denominator;

                vff[facei] = alphaU*vf[P_id] + (1-alphaU)*vf[N_id];

      

            }

            typename GeometricField<Type, fvsPatchField, surfaceMesh>::
                Boundary& bff = vff.boundaryFieldRef();

            forAll(bff, patchi)
            {
                const fvPatchField<Type>& pf = vf.boundaryField()[patchi];
                Field<Type>& pff = bff[patchi];

                if (pf.coupled())
                {
                    tmp<Field<Type>> tpif(pf.patchInternalField());
                    const Field<Type>& pif = tpif();

                    tmp<Field<Type>> tpnf(pf.patchNeighbourField());
                    const Field<Type>& pnf = tpnf();

                    forAll(pff, i)
                    {
                        pff[i] = max(pif[i], pnf[i]);
                    }
                }
                else
                {
                    pff = pf;
                }
            }

            return tvff;
        }

    // Member Operators

};


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

} // End namespace Foam

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

#endif

[HPE]

just randomly shooting - might be the source of problems? discrepancy in fvc::grad(U)

[Tibo99]

Thank for answering,

Actually, I came across this thread some weeks ago and I was suspecting that could be the issue.

But, in the k-epsilon turbulence model, just before resolving the transport equation, we can see these lines here:

tmp<volTensorField> tgradU = fvc::grad(U);
     const volScalarField::Internal GbyNu
     (
         this->type() + ":GbyNu",
         tgradU().v() && dev(twoSymm(tgradU().v()))
     );
     const volScalarField::Internal G(this->GName(), nut()*GbyNu);
     tgradU.clear();

Correct me if I'm wrong, but because of what its described in the thread that you share, they add these lines to make sure that the term are the right one used in the transport equation?

Also, I performed some testing by using my new interpolation scheme files and then instead of implementing the new formula, I was using the following formula:

vff[facei] = (vf[P_id] + vf[N_id])/2;

,and then I compare the results by using 'grad(U) Gauss linear;'. I got the same results.

What is your thought?

Thanks again for answering.