[florian_krause]

Hi grandgo,
if its a division by zero, then its definitely not in

yPlus = y.y() * uTauAvg / RASModel->nu();

or in nuEff what you are using. nu is non-zero as well as nuEff = nu + nu_t. By the way, why do you use nuEff to calculate yPlus?

So division by zero happens then in

uPlus = U / uTauAvg;

which makes sense if I look at the output vaina74 posted: avgUGradWall: 0. Clearly, the velocity gradient at the wall

U.boundaryField()[patchi].snGrad()

gives you zero.

Now, my question is, why is this? I dont know your case, but is it possible by any chance, that this is at time 0 and you initialize your internal velocity field as (0 0 0)? just a guess...

Best,
Florian

[grandgo]

hi florian,

By the way, why do you use nuEff to calculate yPlus?

i'm using nuEff, because i need the tool for a LES case.

nu is non-zero as well as nuEff = nu + nu_t

are you sure? i commented the "uplus line" out and it didnt work as well. same error message.

Now, my question is, why is this? I dont know your case, but is it possible by any chance, that this is at time 0 and you initialize your internal velocity field as (0 0 0)? just a guess...

i changed the internal velocity field a little bit from (0 0 0) to (0.001 0 0), but it's no use. the error occurs not at time t = 0 anyway....


best regards
grandgo

[florian_krause]

Hi grandgo,

my understanding is, kinematic viscosity nu shouldn't be zero, since its a physical property, therefore nuEff is not zero. I also think that you should divide by nu and not nuEff even if its a LES case, see

http://www.cfd-online.com/Wiki/Dimen...e_%28y_plus%29

and the source code of the yPlusLES utility.

It really puzzles me where should be a division by zero if you comment out

uPlus = U / uTauAvg;

since you only divide by nu apart from that.

Is this the only case where you get the error message?

Could you give me a small case where you get the same error message?

Best,
Florian

[grandgo]

Originally Posted by florian_krause

Hi grandgo,

my understanding is, kinematic viscosity nu shouldn't be zero, since its a physical property, therefore nuEff is not zero. I also think that you should divide by nu and not nuEff even if its a LES case, see

http://www.cfd-online.com/Wiki/Dimen...e_%28y_plus%29

and the source code of the yPlusLES utility.

It really puzzles me where should be a division by zero if you comment out

uPlus = U / uTauAvg;

since you only divide by nu apart from that.

Is this the only case where you get the error message?

Could you give me a small case where you get the same error message?

Best,
Florian

hi florian,

i don't know why nuEff instead of nu is used in yPlusLES utility. i guess it makes no difference.

and it's not the only case, where i get an error message. i picked the pitzDaily case in incompressible/pisoFoam/les and the same problem there.


best,
grandgo

[florian_krause]

Thats weird! I just run the pitzDaily case a few second until t=0.017 (random) and used my own plusPostLES utility-it works!

Here are the relevant code snippets of my plusPostLES:

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

    #include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
#   include "createMesh.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
        fvMesh::readUpdateState state = mesh.readUpdate();

        wallDist y(mesh, true);
        // Wall distance
        if (timeI == 0 || state != fvMesh::UNCHANGED)
        {
            Info<< "Calculating wall distance\n" << endl;
            Info<< "Writing wall distance to field "
                << y.name() << nl << endl;
            y.write();
        }

...

#       include "createPhi.H"

        singlePhaseTransportModel laminarTransport(U, phi);

        autoPtr<incompressible::LESModel> sgsModel
        (
            incompressible::LESModel::New(U, phi, laminarTransport)
        );

        volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();
        volScalarField nuEff = sgsModel->nuEff();

        const fvPatchList& patches = mesh.boundary();
        
        dimensionedScalar utauAvg("utauAvg", dimVelocity, 0);

        scalar nPatch = 0;
                
        forAll(patches, patchi)
        {
            const fvPatch& currPatch = patches[patchi];

            if (typeid(currPatch) == typeid(wallFvPatch))
            {
                yPlus.boundaryField()[patchi] =
                    d[patchi]
                   *sqrt
                    (
                        sgsModel->nu().boundaryField()[patchi]
                       *mag(UMean.boundaryField()[patchi].snGrad())
                    )
                   /sgsModel->nu().boundaryField()[patchi];
                const scalarField& Yp = yPlus.boundaryField()[patchi];

                uTau.boundaryField()[patchi] =
                  sqrt
                  (
            sgsModel->nu()
                   *mag(UMean.boundaryField()[patchi].snGrad())
                   );

                const fvPatchScalarField& utauWall = 
                  uTau.boundaryField()[patchi];
                
                dimensionedScalar utauTmp("utauTmp", dimVelocity, average(utauWall));
                
                utauAvg += utauTmp;

                nPatch ++;

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " y+ : min: " << min(Yp) << " max: " << max(Yp)
                    << " average: " << average(Yp) 
                    << " avgUMeanGradWall: " <<  average(mag(UMean.boundaryField()[patchi].snGrad())) << nl << endl;
            }
        }

        utauAvg /= nPatch; 
        
        Info << "Avg. shear velocity is: "
             << utauAvg << "; # of walls: " << nPatch << "." << endl;
       

        yStar = y.y() * utauAvg / sgsModel->nu();
        uStar = U / utauAvg;

Maybe some more info. I am using OF-1.7.x and didn't set any additional FPE options in my bashrc file.

Sorry, at the moment there is nothing else coming to my mind what might cause the error.

Regards,
Florian

[grandgo]

Originally Posted by florian_krause

Thats weird! I just run the pitzDaily case a few second until t=0.017 (random) and used my own plusPostLES utility-it works!

Here are the relevant code snippets of my plusPostLES:

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

    #include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
#   include "createMesh.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
        fvMesh::readUpdateState state = mesh.readUpdate();

        wallDist y(mesh, true);
        // Wall distance
        if (timeI == 0 || state != fvMesh::UNCHANGED)
        {
            Info<< "Calculating wall distance\n" << endl;
            Info<< "Writing wall distance to field "
                << y.name() << nl << endl;
            y.write();
        }

...

#       include "createPhi.H"

        singlePhaseTransportModel laminarTransport(U, phi);

        autoPtr<incompressible::LESModel> sgsModel
        (
            incompressible::LESModel::New(U, phi, laminarTransport)
        );

        volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();
        volScalarField nuEff = sgsModel->nuEff();

        const fvPatchList& patches = mesh.boundary();
        
        dimensionedScalar utauAvg("utauAvg", dimVelocity, 0);

        scalar nPatch = 0;
                
        forAll(patches, patchi)
        {
            const fvPatch& currPatch = patches[patchi];

            if (typeid(currPatch) == typeid(wallFvPatch))
            {
                yPlus.boundaryField()[patchi] =
                    d[patchi]
                   *sqrt
                    (
                        sgsModel->nu().boundaryField()[patchi]
                       *mag(UMean.boundaryField()[patchi].snGrad())
                    )
                   /sgsModel->nu().boundaryField()[patchi];
                const scalarField& Yp = yPlus.boundaryField()[patchi];

                uTau.boundaryField()[patchi] =
                  sqrt
                  (
            sgsModel->nu()
                   *mag(UMean.boundaryField()[patchi].snGrad())
                   );

                const fvPatchScalarField& utauWall = 
                  uTau.boundaryField()[patchi];
                
                dimensionedScalar utauTmp("utauTmp", dimVelocity, average(utauWall));
                
                utauAvg += utauTmp;

                nPatch ++;

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " y+ : min: " << min(Yp) << " max: " << max(Yp)
                    << " average: " << average(Yp) 
                    << " avgUMeanGradWall: " <<  average(mag(UMean.boundaryField()[patchi].snGrad())) << nl << endl;
            }
        }

        utauAvg /= nPatch; 
        
        Info << "Avg. shear velocity is: "
             << utauAvg << "; # of walls: " << nPatch << "." << endl;
       

        yStar = y.y() * utauAvg / sgsModel->nu();
        uStar = U / utauAvg;

Maybe some more info. I am using OF-1.7.x and didn't set any additional FPE options in my bashrc file.

Sorry, at the moment there is nothing else coming to my mind what might cause the error.

Regards,
Florian

hi florian,

why are you using UMean now? and what's Umean?

and whats the difference between utauavg and utauavg.value() ?

[florian_krause]

why are you using UMean now? and what's Umean?

sorry for the confusion with UMean in the code snippet... I used the utility for a fully developed turbulent pipe flow and time-averaged the velocity field. With that I calculated the mean (time-averaged) velocity gradient at the wall... I only wanted to show you the code of the utility I used for the pitzDaily case.

and whats the difference between utauavg and utauavg.value() ?

I don't understand what you mean by this! Are you talking about this one

utauAvg += utauTmp;

Its just a dimensionScalar with the dimensions of velocity, see defintion of utauTmp. Its not a class with a member value().

Regards,
Florian

[grandgo]

Originally Posted by florian_krause

sorry for the confusion with UMean in the code snippet... I used the utility for a fully developed turbulent pipe flow and time-averaged the velocity field. With that I calculated the mean (time-averaged) velocity gradient at the wall... I only wanted to show you the code of the utility I used for the pitzDaily case.

I don't understand what you mean by this! Are you talking about this one

utauAvg += utauTmp;

Its just a dimensionScalar with the dimensions of velocity, see defintion of utauTmp. Its not a class with a member value().

Regards,
Florian

no, in the original code you user this code

Info << "  avg. friction velocity uTau is: "
             << uTauAvg.value() << " (averaged over " << nPatch << " wall(s))" << nl <<endl;

i'm just wondering...

i don't know what the problem could be. i'm a newbie in c++ programming

the utility would fit perfectly to what i need.
i just go on trying and testing

best
grandgo

[florian_krause]

Its just a dimensionScalar with the dimensions of velocity, see defintion of utauTmp. Its not a class with a member value().

OK now I should say sorry for the non-sense...

It is a class, would have to check the documentation for the exact definition... but think it comes from the generic dimensioned Type class dimensioned<Type> with member function value(). That means, if you only use uTauAvg it will output name, value and dimension.

Even I am not a C++ expert...
sorry again!

[grandgo]

hi florian,

i'm struggling with the modification

the thing is: i'm not changing much at all and it still doesn't work.

here is what i changed in the createFields.H file:

        autoPtr<incompressible::LESModel> sgsModel
        (
            incompressible::LESModel::New(U, phi, laminarTransport)
        );

and in the .C file i substituted RASModel->nu() with sgsModel->nu().
that's all.

the openfoam version is 1.7.1.
i used the utility with the tutorial case: /incompressible/pisoFoam/les/pitzDaily, too. no chance.

i dont understand this.

EDIT: i think, my problem is history.

these two lines in the modified utility caused the problems in my case:

1) yPlus = y.y() * uTauAvg / nuEff;
2) uPlus = U / uTauAvg;

the problem with the first line was, that i used nuEff instead of sgsModel->nu(). and nuEff did become zero (nuEff = nu + nuSgs, with negative values of nuSgs possible). so i changed to sgsModel->nu().

the second line caused problems, because uTauAvg was zero at time = 0. as florian assumed, my 0/U internelField value is (0 0 0). so the solution is either to change the internalField value, or, as i did, change the code of the utility a little bit:

timeSelector::addOptions(false, true);

now the 0 directory is excluded in the calculation and everything is fine.

best
grandgo

[jona]

Hi *,

for converting the code of Florian to OF 2.0, "RASModel->nu()" had to be explicitely treated as "volScalarField&".

The corrected code is attached.

Have fun,

Jona

[fakekarma]

Hi jona,

I've tried to build the utility as follow:

:~/OpenFOAM/OpenFOAM-2.0.1/applications/utilities/postProcessing/wall/plusPostRANS$ wmake libso
wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file plusPostRANS.C
SOURCE=plusPostRANS.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-100 -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/meshTools/lnInclude -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/transportModels -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/turbulenceModels -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/turbulenceModels/incompressible/RAS/RASModel -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/turbulenceModels/incompressible/LES/LESModel -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/turbulenceModels/LES/LESdeltas/lnInclude -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/OpenFOAM/lnInclude -I/home/ag-peinke/OpenFOAM/OpenFOAM-2.1.0/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/plusPostRANS.o
'libNULL.so' is up to date.

so it seems to be no error, but I should now use it it, since it doesn't work simply typing plusPostRANS in the case folder. Should I do something more???

Thanks for your help,

Cheers,

Elia

[wouter]

hello,

type wmake, not wmake libso

hope this works
Wouter

[fakekarma]

Hallo Wouter,
with the command "wmake" works as it should for compiling.
Thanks very much.

Elia

[chegdan]

Hello All,

First of all, thanks to Florian for making this, I have used some other yPlus utilities floating around here and this one was particularly interesting. I was attempting to get a uPlus and yPlus for something other than a pipe flow and I have been reading this thread thinking about the problem. Have made a modification to the plusPostRANSUtility utility provided (for OF 2.1.x) and I nearly get the same answer on Florian's test case (running pisoFoam) with my utility and his. However, it is (edit was) dreadfully slow since it searches for the nearest uTau over all wall faces for each cell.

Right now it calculates a new y value (we will call it cellToFaceDistance), compares it to the current already calculated y value created using wallDist, and if there is a user specified fraction difference (i.e. (cellToFaceDistance-y)/cellToFaceDistance) than it is close enough and grabs that uTau for that cell and stops searching.

I have uploaded the code to github and it can be downloaded with

git clone git://github.com/chegdan/yPlusUplus.git

It is currently made for OF 2.1.x

EDIT: However, after thinking a bit...I'm not quite sure if this sort of u+ vs y+ is useful in cases where there is boundary layer separation around objects. For something where the streamlines are nice close to the object, it makes sense. For something with separation, it may not be best to non-dimensionalize the velocity field with respect to the friction velocity but the mean velocity. Either way, it was a good bit of code and if anyone can find it useful, enjoy.

Cheers.
Dan

[The King]

Hi! Did you found the error. I get the same, nothing found yet...

[eysteinn]

Originally Posted by florian_krause

Dear all,

This issue was already discussed in some length in various threads. Therefore I finally reviewed and attached my plusPostRANS utility.

It calculates y+ and u+ values when using one of the available low-Re RANS turbulence models. I hope the code as well as the output will be self-explaining

I also attached a case. It is a straight pipe (wedge) with periodic inlet/outlet boundary conditions. The Reynolds number based on the mean axial velocity is Re=5300. You can run the case using pisoFoam.

I would appreciate any feedback, it might be still improved at some point!

Best regards,
Florian

Hi Florian,

I have one question regarding your example case.
I noticed that you use zero gradient for epsilon and nut at the wall. Can you ellaborate on why that is?
(I somehow expected epsilon = 1e-9 or something small as you are solving for epsilon tilda and calculated for nut.)

Best regards,
Eysteinn

[owayz]

Hi Jona and Florian,
Thanks for the utility I am sure it will help a lot of people.
I am also interested in calculating the y+ for my case. But I have compressible simulation.
I want to know two things:
1- If its a transient simulation shouldn't we use the mean values of U, rho and other variables?
2- How could I extend this utility to calculate the same numbers of compressible flow?

Regards,
Ali

[owayz]

Hi,
I have made some small changes in the utility, developed by Florian/Jona. I have tried to make it work with compressible simulations.
Please check it, if you see its working as it is supposed to work.
Regards,
Awais Ali
P.S: For my case it reports a very low y+ values on the wall (unexpectedly). But the value that it writes seems to be Ok.

[owayz]

Hi All,
I have written this utility (again by modifying the utility previously developed) to calculate y+ and u+ values for compressible LES simulation.
Please let me know if you find this utility working. For me its working fine. I have used on my own case. The yPlus values are as they were expected to be.

Regards,
Awais