[ArathoN]

After some times I retried to simulate a channel with cyclic BCs with the solver pimpleFoam and using RAS models (i want only to compute the mean velocity profile).

The problem is i can't get the pressure variable to converge it will always stay at around 0.4-0.6 and I really tried everything.

First some info on the geometry and physical properties:

H=0.025m
Lx=5*H  (i tried longer Lx but nothing changed)
nu=1.5e-5
Re(H,Uc)=15000
Re(H,Ubar)=13500
I=7%  (from experimental data that i have)
with Uc=velocity at centerline, Ub=Ubar and H = width of the channel.

I estimated y__firstcell with yplus=1,4,30 and considering different cell-to-cell expansion (keeping it under 1.15), in the file attached yPlus is around the unit. I tried increasing the yplus until yplus=40 but this didn't change the weird behaviour of the pressure and Uy residuals.


Now the initialization of teh turbulent flow are

K=0.48

For epsilon this is a little bit tricky, now i've read all the code about the solver and the turbulence correction and i noticed that in the "turbulenceModels" epsilon has different definition depending if it is related to a cell adjacent to teh wall (via teh wallfunction, there is a lone where it will update G and epsilon) or not (e.g in kEpsilon.C there is teh update of nut values and from this relation you can define epsilon). So there are two possible definition one based on the turbulence length (o mixing length) and the other on nut/nu.

epsilon(l)=22  with l=10%H
epsilon(nut/nu)=139   with nut/nu=10

The mesh should be fine:

Overall domain bounding box (0 0 0) (0.125 0.025 0.025)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-1.15504e-18 7.86934e-17 2.92744e-16) OK.
    Max cell openness = 2.51663e-16 OK.
    Max aspect ratio = 22.206 OK.
    Minimum face area = 7.0364e-08. Maximum face area = 3.3667e-05.  Face area magnitudes OK.
    Min volume = 1.7591e-09. Max volume = 4.20837e-08.  Total volume = 7.8125e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.42118e-06 OK.
    Coupled point location match (average 8.25619e-18) OK.

What i tried to modify:
-schemes and fvsolution parameters
-initial conditions
-yplus and cell-to-cell expansion (or R)
-lenght of the channel
-Bcs Wall function or not
Unfortunately this didn't help.

Here the archive with all the files, I atatched also one graph of teh residual profile (done with pyFoam). I hope someone could figure out the problem.

As a side note even with the LES simulation i couldn't get the right redisuals, the p and uy residuals were the problems. I used the same flow properties with a finer mesh.

EDIT: Maybe i've found the error, if i choose in controldict to run the simulation with fixed deltat the courant number will explode after some steps. I don't know how to solve it, I tried changing the timeStep with no success. What velocity U and spacing do you consider when defining courant number? I usually try to get the lowest timestep so i use the highest velocity and the lowest spacing grade (I

EDIT2: I tried also a simulation with no turbulence but again the residual floats around the unit.


OT: I have 3 questions:
1- what relation do you use to define epsilon?
2- what do you put in the label pRefCell? I've chosen first 1001 but it gave me an error then i switched to 101, I know that it is teh identification on one cell of the domain where you initialize the pressure (so you'll not have an ill-conditioned problem) but it the choice of the cell important?
3- Reading teh code I noticed that the utility yPlusRAS only calls the calculation of yPlus done in the WF. Moreover studying the WF there two different definitions of the yPlus(utau,y,nu) and yStar(k). However i didn't find nothing RELIABLE about the comparison of the two definitions and on how they interact with each other (i know how to define them from the wall equations)

These question are only my left doubts after going deep into the code and i couldn't answer them.

[RodriguezFatz]

Hi.
I am pretty sure we can get this running!

To your second question: pRef is not used for initialization as I understand it, but for setting a reference point in the pressure field for the complete simulation. Thus, the cell you choose is set to a fixed pressure (as you declare it with "pRefValue = 0"). This is needed for incompressible simulations, if there is no absolute pressure information at the boundaries. Anyway, it never changes the results of the simulation, because it is just an offset in the pressure field. You can choose any cell and any value.

Ok, now to solve the problems.
1) The mesh you sent in the zip-file: What y+ do you expect here? I see you use RAS-model kepsilon, so since this is a high-Re model I guess you have y+ >> ?
2) How can I create the mesh? I only import files from ICEM, so I don't know how to use the blockMeshDict.
3) Why do you use a transient solver? Why don't you use simpleFoam?
4) Is this 2D?

[ArathoN]

Thank you for your help.

Originally Posted by RodriguezFatz

1) The mesh you sent in the zip-file: What y+ do you expect here? I see you use RAS-model kepsilon, so since this is a high-Re model I guess you have y+ >> ?

In that mesh yPlus is around the unit, but i considered yplus=1,2,4,30,50. I wanted to see if teh yPlus was the problem or not, unfortunately this didn't change the results. I know that kEpsilon if a High-Re models but i thought that it should give results even in low-Re meshes even though they will be quite wrong. I'll test some other turbulent models.

Originally Posted by RodriguezFatz

2) How can I create the mesh? I only import files from ICEM, so I don't know how to use the blockMeshDict.

The mesh is created by the utility blockMesh so type in the terminal:

blockMesh

Originally Posted by RodriguezFatz

3) Why do you use a transient solver? Why don't you use simpleFoam?

I'm using pimpleFoam because in this forum they said it can work with the cyclic BCs but after i read the code of pimpleFoam and simpleFoam I found out that both of them should support this BC and the pressureGradientExplicitSource argument. Infact they both call the fvOptions.correct(U) to correct the velocity field.
I tried one simulation with simpleFoam but i had the same results.

Originally Posted by RodriguezFatz

4) Is this 2D?

yeah, I only want the mean velocity profile and not the perturbations fields.


I think i'm doing something really stupid but can't figure out where, because even with the les simulation i had the same problem, the residuals of (p,Uy) were always too high around 0.5-1.

[HanSolo123]

Originally Posted by ArathoN

yeah, I only want the mean velocity profile and not the perturbations fields.

I think this is wrong. To get the mean velocity profile, you need to average the (turbulent) fluctuations for a long time. To get those fluctiations, you have to recover the turbulent elements in the flow. This can only be done in a 3D simulation, not 2D!

BTW: your k, which is needed to solve the NSE, is build from the 3 velocity components. But in your case, there are only 2. When you have a look at DNS data or in those LES threads, you can easily see that the 3rd component has a huge impact. When you leave it, you will make a mistake.
So my advice: go for 3D.

[ArathoN]

Originally Posted by HanSolo123

I think this is wrong. To get the mean velocity profile, you need to average the (turbulent) fluctuations for a long time. To get those fluctiations, you have to recover the turbulent elements in the flow. This can only be done in a 3D simulation, not 2D!

BTW: your k, which is needed to solve the NSE, is build from the 3 velocity components. But in your case, there are only 2. When you have a look at DNS data or in those LES threads, you can easily see that the 3rd component has a huge impact. When you leave it, you will make a mistake.
So my advice: go for 3D.

I know that a 3D simulation is way better than a 2D, however i can't stay days waiting for a solution and I only need a fully developed turbulent velocity profile to use it as an entry for other cases.

Tbh the RAS simulation, given the simplicity of the case, will always give with the right condition and mesh a mean velocity profile similar to the experimental data. However the perturbation will be far away from the real data especially in the regions near the mixing layer because they are totally modeled with and assumption that in these region isn't so valid (hyp. boussinesq).

About the averaging, the RAS equation are based on the mean value of the variable and its perturbation for every time step, so in my opinion the velocity field represent the mean part of that variable. Let me explain properly in the turbulence equation (eg. pEqn or Ueqn) you are using the mean value of the velocity (see every textbook about RANS eg. turbulent flow of Pope) so the output will be the mean value. But this can be wrong.

[RodriguezFatz]

This case will converge to steady state, so if you use a RAS model, 2d is absolutely sufficient. It will give exactly the same result as 3d with symmetry b.c. in the third dimension. You are right, arathon. Also, since it will converge to steady state, there is no need to "average the fluctuations for a long time" since there are no fluctuations. All fluctuations are captured by the RAS model.

About the turbulence model: You use pretty wrong boundary conditions. First you must decide whether you want to have an y+ of 1 or let's say 100.
For the first case you need a low-Re model with low-Re boundary conditions. Let's say you go for the LaunderSharmaKE. For the second case you need a high-Re model with high-Re b.c. As I read in this forum the standard kepsilon is a bad choice since it is pretty unstable. Better use the realizableKE.

Then decide on the solver. I think you want simpleFoam. You need to change the fvSolution file so it knows what to do when SIMPLE is started. What your fvSolution is also missing is the under-relaxation. You won't get simple stable without under-relaxation. I didn't get PIMPLE stable without it, too.

What do I need to change in the blockMeshDict to get a high- and a low-Re mesh? If you tell me, I will get you both cases running...

[ArathoN]

Originally Posted by RodriguezFatz

This case will converge to steady state, so if you use a RAS model, 2d is absolutely sufficient. It will give exactly the same result as 3d with symmetry b.c. in the third dimension. You are right, arathon. Also, since it will converge to steady state, there is no need to "average the fluctuations for a long time" since there are no fluctuations. All fluctuations are captured by the RAS model.

About the turbulence model: You use pretty wrong boundary conditions. First you must decide whether you want to have an y+ of 1 or let's say 100.
For the first case you need a low-Re model with low-Re boundary conditions. Let's say you go for the LaunderSharmaKE. For the second case you need a high-Re model with high-Re b.c. As I read in this forum the standard kepsilon is a bad choice since it is pretty unstable. Better use the realizableKE.

Then decide on the solver. I think you want simpleFoam. You need to change the fvSolution file so it knows what to do when SIMPLE is started. What your fvSolution is also missing is the under-relaxation. You won't get simple stable without under-relaxation. I didn't get PIMPLE stable without it, too.

What do I need to change in the blockMeshDict to get a high- and a low-Re mesh? If you tell me, I will get you both cases running...

Thank you very much, I noticed that kEpsilon is totally bad in OF with respect to Fluent but wrongly i thought that it will give converge even with wrong solution. I tried with the same BCs but using the kOmegaSST (with the appropriate fix to make it work on low-re meshes) and it converged pretty fast. I used a different Re from experimental data so I'll validate the solution today and i'll update.

I've confirmed that in OF 2.3 simplefoam supports cyclic BCs so i'll switch to it. I know about the under-relaxation, reading the OF code of solver and turbulence models was really a good journey of discovering how the concept learn on uni lectures are applied.

Last about the yplus, i evaluate the y of the first cell center from wall using the desired yplus (if you want i can write all the process). Then i calculate the spacing deltay of the first cell and using this utility i define the number of grid points and the grading.
In block/hex you have the definition of the block then the grid point in (x,y,z) afterward the grading (in OF it's the total expansion ratio deltalast/deltafirst).

[RodriguezFatz]

Originally Posted by ArathoN

Last about the yplus, i evaluate the y of the first cell center from wall using the desired yplus (if you want i can write all the process). Then i calculate the spacing deltay of the first cell and using this utility i define the number of grid points and the grading.
In block/hex you have the definition of the block then the grid point in (x,y,z) afterward the grading (in OF it's the total expansion ratio deltalast/deltafirst).

Sorry, what I meant was: What line and what value do I have to change with the current setting of Ubar=8.1m/s to get 1) a low-re and 2) a high-re grid?

[ArathoN]

Originally Posted by RodriguezFatz

Sorry, what I meant was: What line and what value do I have to change with the current setting of Ubar=8.1m/s to get 1) a low-re and 2) a high-re grid?

for yplus=1, here some different setup with different cell-to-cell expansion (max 1.15):

	Ny 	Ryb	        Ryt
	55	10.37	0.0964
	40	16.25	0.0615
	30	23.92	0.0418
	25	30.28	0.0330

for yplus=100,30

Y+=30
	Ny 	Ryb     	Ryt
	7	1.116	0.8961
	7	1	        1.0000

Y+=100
	Ny 	Ryb	        Ryt
	2	1.232	0.8116883117

Ny=Nyb=Nyt.

You need to change the value of the variables at the top of blockMeshDict file and then type the command

blockMesh

EDIT: I tried using simpleFoam with kOmegaSST but omega after some timesteps will explode and now the same happens with pimpleFoam.

[RodriguezFatz]

These channels are crap... just give me some more time. I don't know why this "just runs" in Fluent and I always have stability problems in OpenFoam like hell...

[ArathoN]

Originally Posted by RodriguezFatz

These channels are crap... just give me some more time. I don't know why this "just runs" in Fluent and I always have stability problems in OpenFoam like hell...

I confirmed that channel flow with cyclic plus pressure gradient (enablem on fvOptions) doesn't work with steadyState condition.

If i choose simplefoam Ubar and gradP (and residual of turbulence variable) will rise until the solver will crash, the same if you use pimpleFoam with ddtschemes as steadyState. But if i choose the ddtscheme as Euler the simulation will not crash.

I'm really curious on why this happens, with steady state the Ubar and gradp will go aroud e10 at the first step. I tried to change schemes and bcs with no success, I even changed deltaT.


I'm looking at the code:

Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve
        << ", pressure gradient = " << gradP << endl;

with

dGradP_ = (mag(Ubar_) - magUbarAve)/rAUave;
scalar gradP = gradP0_ + dGradP_;

and

 
magUbarAve += (flowDir_ & U[cellI])*volCell; (for all the cells)
magUbarAve /= V_;

With simplefoam teh outputted values are (Ubar=8.1):

Pressure gradient source: uncorrected Ubar = 2.27933, pressure gradient = 56090

however at the initial time step all the internal flow should have a velocity almost equal to Ubar

magUbarAve = Ubar

because the velocity U[cellI] initialized as equal to Ubar in /0/U is parallel to Ubar then sum(volCel)/V should be alomost equal to teh unit therefore dGradP had to be small.

IMO the probem lies on the instabilities caused by the absence of the temporal derivative term that is one of the contributor to the convergence and stability of the solver (idea made after i've read "Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows " by Hrvoje Jasak). For this reason in simplefoam you need the underelaxaion factors to mitigate the absence of the temporal derivate contribution. Then this maybe affect the rAU matrix, so that the instabilities will propagate to the variable dGradP, infact there is no correction introduced for such a situation.

However i can't explain how the Ubar value is like that inthe first timestep.

Edit: maybe the root of the problem:

magUbarAve =/= Ubar

is caused bythe underelaxation done to UEqn this maybe give a velocity field lower based on the URF. I'll try using URF for U as 1(so no UR) and see what happens. I'll update as soon as teh simulation will finish. the problem.

Edit 2: changing URF doesn't change the value of magUbarAve.

[RodriguezFatz]

Originally Posted by ArathoN

I confirmed that channel flow with cyclic plus pressure gradient (enablem on fvOptions) doesn't work with steadyState condition.

If i choose simplefoam Ubar and gradP (and residual of turbulence variable) will rise until the solver will crash, the same if you use pimpleFoam with ddtschemes as steadyState. But if i choose the ddtscheme as Euler the simulation will not crash.

I'm really curious on why this happens, with steady state the Ubar and gradp will go aroud e10 at the first step. I tried to change schemes and bcs with no success, I even changed deltaT.

OK, i think this time you are wrong. I did this before - with a channel and with a short pipe - and it works. It's just so freaking unstable / dependent in a lot oft settings at the beginning...

[ArathoN]

Originally Posted by RodriguezFatz

OK, i think this time you are wrong. I did this before - with a channel and with a short pipe - and it works. It's just so freaking unstable / dependent in a lot oft settings at the beginning...

Only if i change the dttschemes the simulation will crash, it's like he doesn't like any steady solver.

I even though maybe the pressure predictor is the cause and i've set it as 1 (so pisoFoam) but the simulation is running normally. I changed schemes, BCs, turbulence model and mesh properties but the only thing that make the simulation crash is the temporal derivative scheme.

From what i gathered in the previous post with steady solver it computes wrongly the Ubar averaged over the domain and this leads to the crash. Unfortunately i can't understand why because the evaluation of volCell and V() is the same with steady or unsteady solvers as for U[CellI].

[RodriguezFatz]

You did something wrong with simplefoam. I will send you the low-re case tomorrow. I don't know if the high re case works because it just such few cells in y direction.

[ArathoN]

Originally Posted by RodriguezFatz

You did something wrong with simplefoam. I will send you the low-re case tomorrow. I don't know if the high re case works because it just such few cells in y direction.

ok thank you, if you want i can send you my case with kOmegaSST.

By the way how do you choose the grading of the grid after fixing yplus?

I know that it should stay under 1.2 (cell.to-cell expansion), the aspect ratio should better be under 20 (read it in some book advised in openfoamwiki) and for wall resolved we need at least 10 pints into the boundary layer. Are there any other consideration you account off?


EDIT: i plotted p with paraFoam and i get such results, aren't they quite strange? The pressure should always decay, why the is a rise in p? In this case Lx=20H shouls i lower it to filter this.

[RodriguezFatz]

ArathoN,

Unfortunately, I didn't manage to get it converged.
I get a picture quite similar to the residual plot of your first post. Surprisingly, the results look good and match pretty good with Fluent's results anyway. Now, I find y-velocity values (absolute) of about +/- 1e-14 in my domain. I don't understand why these y-velocities of nearly zero can produce such a high residual. And why the pressure residual remains that high.

Anyway, for the LaunderSharmaKE model, you can go with these b.c. (low-Re mesh!!!):
k:

internalField   uniform 1;
boundaryField
{
    "(bottomWall|topWall)"
    {

        type            fixedValue;
    value         uniform 1.0e-12;
    }

    "(inlet_top|outlet_top|inlet_bot|outlet_bot)"
    {
        type            cyclic;
    }

    frontAndBack
    {
        type            empty;
    }
    
}

epsilon

internalField   uniform 0.1;
boundaryField
{
    "(bottomWall|topWall)"
    {
        
    type fixedValue;
        value           uniform 1e-8;

    }
    
    "(inlet_top|outlet_top|inlet_bot|outlet_bot)"
    {
        type            cyclic;
    }

    frontAndBack
    {
        type            empty;
    }

}

nut

internalField   uniform 1e-3;
boundaryField
{
    "(bottomWall|topWall)"
    {
        type            nutLowReWallFunction;
        value           uniform 0;
    }

    "(inlet_top|outlet_top|inlet_bot|outlet_bot)"
    {
        type            cyclic;
    }

    frontAndBack
    {
        type            empty;
    }
}

p and U just as You already do it.

fvSolution:

solvers
{
       "(p|pFinal)"
     {
    solver           GAMG;
        tolerance        1e-99;
        relTol           1e-5;
        smoother         DICGaussSeidel;
        nPreSweeps       0;
        nPostSweeps      1;
    nFinestSweeps    2;
    scaleCorrection  true;
    directSolveCoarsestLevel false;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 500;
        mergeLevels      1;
        maxIter         15;
     }


    "(U|k|epsilon|omega|UFinal|kFinal|epsilonFinal|omegaFinal)"    
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-99;
        relTol           1.0e-6;
    maxIter         15;
        };
}

SIMPLE
{
        nNonOrthogonalCorrectors 0;
    residualControl
    {
            p               1e-10;
            U               1e-10;
        k        1e-10;
        epsilon        1e-10;
        omega        1e-10;
    }
    pRefCell            1;
    pRefValue           0;
}

relaxationFactors
{
    fields
    {
        p                 0.3;
    }
    equations
    {

        "(U)"             0.7;
        "(k)"             0.5;
        "(epsilon)"       0.5;
        "(omega)"         0.8;
    }
}

and fvSchemes:

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
   
           
           default Gauss linear;
        grad(sqrt(k)) faceMDLimited leastSquares 1;
}

divSchemes
{
    div((nuEff*dev(T(grad(U))))) Gauss linear;
    div((nuEff*T(grad(U)))) Gauss linear;

    div(phi,U)  bounded Gauss linearUpwind grad(linUpw); 
    div(phi,k)  bounded Gauss linearUpwind grad(linUpw); 
    div(phi,epsilon) bounded Gauss linearUpwind grad(linUpw); 
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p               ;
}

I checked all different divSchemes and gradient with/without limiters. I don't think this is the reason for the residual. I always managed to get these cases to convergence eventually. Maybe someone else can help here...

[ArathoN]

Finally I was able to obtain a converged solution with kOmagaSST and the u-profile is in accordance with the experimental data. I tried the laudershamaKE and while it didn't converge, res(p) ~0.1 & res(uy) ~0.01, the u profile was good and quite similar to komegasst.

By the way I reduced the tolerances on fvsolution yours were too restrictive and I noticed that after e-9 the solution will not change.

Finally when I see the flow in parafoam I noticed a weird behavior of p/uy at the inlet for komegasst and at the upper wall for laundersharma , there are pressure/uy spikes for no reason typical of a periodic condition on paper. See the images maybe this is the cause of the residual profile and maybe the pressure gradient option is still a bit bugged.

I couldn't upload the images with the Android app, here a link to imgur http://imgur.com/a/t9lJy

[RodriguezFatz]

ArathoN, you can find a low-Re version of the kOmegaSST model here (post #5):
http://www.cfd-online.com/Forums/ope...tml#post491414
I implemented the same version that is implemented in Fluent, when you check the box "low-Re-blabla" in Fluent. The regular SST model in openFoam does not include the damping functions near the wall. Thus, even if you set low-Re b.c. in openFoam, you will not get matching profiles necessarily.

[ArathoN]

I'll try it and report, i noticed that the laundersharmaKE gives better evaluation of wallshearstress hence better utau however the uprofile Is a bit worser.
OT: can you share with me the pdf of menter i couldn't find it.

[ArathoN]

Originally Posted by RodriguezFatz

ArathoN, you can find a low-Re version of the kOmegaSST model here (post #5):
http://www.cfd-online.com/Forums/ope...tml#post491414
I implemented the same version that is implemented in Fluent, when you check the box "low-Re-blabla" in Fluent. The regular SST model in openFoam does not include the damping functions near the wall. Thus, even if you set low-Re b.c. in openFoam, you will not get matching profiles necessarily.

I tried the low-re version and it improved a little bit the data, however this would work only for the channel in a backward facing step it worsened the solution.

I confirm that the pressure is the main problem when using cyclic condition with the simpleFoam solver, if you try to visualize the pressure field in parafoam after simulation you'll see that it is non-sense and unphysical (random high and low pressure spots). I tried to use a different Reynolds number and the same thing happened the pressure and Uy residual are a mess however the U-profile is quite good.