CFD Online Discussion Forums - Interpolation Scheme on grad(U)

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.
Code:
#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