CFD Online Discussion Forums - OF div schemes greatly influence the result

Dear all:
I'm doing a simulation about wind aroud a flat roof. The model can be seen below.
[image]https://www.cfd-china.com/assets/uploads/files/1611984747439-wholemodel.png[/img]

In the reference paper, the author applied 'second order upwind' for space discretization in Fluent. And I have used the same model parameters (eg. domain size, minimum element size, boundary condition, turbulence model). simpleFOAM and realizable k-epsilon model are used. However I find that if I adopt different divergence schemes (namely, gauss upwind, linearUpwind and limitedLinear), the results will be different from each other and only the result of 'gauss upwind' agrees with that of reference paper.

From what I gather, in openfoam 'gauss linearUpwind' represents 'second order upwind' which was adopted in the reference paper. But when I use 'gauss linearUpwind' for the divergence schemes in OpenFOAM, the simulated result (namely, friction velocity on roof) differs greatly from that of reference paper. In comparison, the performance improves if I adopted 'gauss upwind' which represented first order upwind. The results of different div schemes can be seen below.
[image]https://www.cfd-china.com/assets/uploads/files/1611984885997-0f0345fa-d59b-4030-9578-0dcb60ab0d6f-image.png[/img]

To describe my question detailedly, some details of my fvSchems and fvSolutions are also posted.

Code:

ddtSchemes

{

    default         steadyState;

}



gradSchemes

{

    default         Gauss linear;



    limited         cellLimited leastSquares 1;

    //grad(U)         $limited;

    grad(k)         $limited;

    grad(epsilon)     $limited;

}



divSchemes

{

    default         none;



    div(phi,U)      Gauss limitedLinear 1;//bounded Gauss linearUpwind limited, Gauss upwind;



    turbulence      bounded Gauss linearUpwind limited;

    div(phi,k)      $turbulence;

    div(phi,epsilon) $turbulence;



    div((nuEff*dev2(T(grad(U))))) Gauss linear;

}



laplacianSchemes

{

    default         Gauss linear corrected;

}



interpolationSchemes

{

    default         linear;

}



snGradSchemes

{

    default         corrected;

}



wallDist

{

    method meshWave;

}

Code:

solvers

{

    p

    {

        //solver          PCG;

        //preconditioner  DIC;

        solver          GAMG;

        tolerance       1e-7;

        relTol          0.1;

        smoother        DICGaussSeidel;

    }



    "(U|k|epsilon|omega)"

    {

        solver          PBiCGStab;

        preconditioner  DILU;

        //solver          smoothSolver;

        //smoother        symGaussSeidel;

        tolerance       1e-7;

        relTol          0.001;

    }

}



SIMPLE

{

    nNonOrthogonalCorrectors 2;

    consistent      yes;



    residualControl

    {

        p               1e-7;

        U               1e-7;

        "(k|epsilon|omega|f|v2)" 1e-7;

    }

}



relaxationFactors

{

    fields

    {

        p               0.3;

    }

    equations

    {

        U               0.7;

        "(k|omega|epsilon).*" 0.7;

    }

}

Actually, in many papers, first order upwind scheme are not encouraged to use. But in my simulation, it seems that I have to choose 'gauss upwind' to get a 'right' result compared with reference papar of which the author has adopted second order upwind schemes. And higher-order div schemes lead to something wrong.
That makes me confused. Could anyone please give me some suggestions?
Thank you
samuel