[KTG]

Hi fellow foamers,


I am doing some k-epsilon steady sate simulations (simplefoam), which I don't have much experience with. I have a geometry that forces flow through a narrow pipe, which causes high gradients in velocity and pressure. Everything is converging, but the pressure is going super slow after a certain point.






I am trying to get the initial residuals below 1e-06. All quantities except the pressure get there reasonably quickly. I have to use insane relaxation factors though- like 0.00000000000001. I am guessing this is because of the pressure gradient, but I am also thinking I can't push that too much further without hitting numerical issues?


Any ideas are much appreciated. I am using:


Grad schemes:

cellLimited Gauss linear 1
Div schemes:
Gauss linearUpwindV cellMDLimited Gauss linear 1


Thanks!

[KTG]

looks like my residual plot didn't show up - its attached.

[KTG]

I am still having a lot of difficulty here. I noticed that the pressure is checkerboarding in my simulation, I have attached a picture. In the picture, a high speed flow is driven downward through the small tube pictured. It opens into a larger chamber, so there is a big pressure gradient where that checkerboard is occurring. This is strange, since the pressure is slowly converging and I am using bounded schemes. I have tried using an alternate turbulence model and refining the mesh, both approaches didn't really change much. I tried changing all the divergence schemes to upwind, and the problem is still there!

[KTG]

Here is another update on a Friday night. I probed the pressure through the small tube (backwards) and attached a figure. Pressure should be flowing from the high value (right side) to the zero value (left side). However, when it enters the tube it jumps down to some insane negative number! The strange oscillations at the boundaries of the small tube did not surprise me, but how is it possible that the velocity is converging so well? Its looks correct visually, and the residuals and continuity errors are really small. The flow is completely defying the pressure gradient and it doesn't make any sense.

[Tobermory]

Other things to check out / play with:

1. your mesh design/quality/resolution (is it a tet mesh, or hex?)
2. fvSchemes for your laplacian

You should not need to have such low relaxation factors - that's probably a sign that there is trouble elsewhere.

[Tobermory]

Originally Posted by KTG

Its looks correct visually, and the residuals and continuity errors are really small. The flow is completely defying the pressure gradient and it doesn't make any sense.

Maybe a result of you screwing down the relaxation factors so low? That effectively prevents the solution rebalancing ... If the model size is not too large, feel free to upload and maybe I can take a look and see if I can spot anything obvious.

[KTG]

Ahh yes, I think you are probably right about the relaxation factors allowing for that strange behavior. I wasn't getting any convergence without them, but had started with a potential flow solution, which was probably allowing the velocity to converge like that. I just doubled the mesh resolution and put all the relaxation factors at .5, and now the pressure is converging faster than anything else. I think the insane pressure gradients / velocities at the corners are causing the problem. I'll try messing with the mesh a little more and report back.



On a related note though, why does the yplus function object not seem to work for these types of simulations? I am getting about this for whatever mesh / velocity combination I try:

    patch walls y+ : min = 4.42377e-09, max = 0.284155, average = 0.0236263

I actually went up from the max of about .1 when the mesh was coarser and the velocity higher. Is there some setting I don't know about for this?


Thanks!

[KTG]

Well, I think I am getting closer with this, but am not there yet.

• I was able to get a convergence at the 10^-6 level using all upwind schemes. Could probably get additional orders of magnitude on the residuals.
• Tried going to div(phi,U) Gauss linearUpwindV cellMDLimited Gauss linear 1; and leaving everything else at upwind. The residuals hit a steady state that was not very good for pressure and velocity, though k and epsilon made it to around 10^-5.
• The yplus issue I mentioned earlier was due to the solution being poorly converged. In my final solution the yplus was in the high 20s.
• In the converged solution, epsilon in the small tube is very high - like 10^7. Does this make sense? In solutions that don't converge, epsilon blows up early on - usually within 5 iterations depending on how much relaxation there is.

Anyone have any ideas on how to get this to converge with numerical schemes more accurate than just plain upwinding? I attached a picture of the epsilon - velocity is about 86 m/s in the small tube. Thanks for reading!

[Tobermory]

I'm glad you're making progress. On the epsilon issue, remember that in a shear layer we can write , where and are turbulent velocity and length scales. So, if your velocity is high and lengthscale is small, then yes, epsilon will be large.

What did you do for your Laplacian and snGrad schemes, by the way? Did you look at reducing the limiter? This can help stabilise turbulence diffusion in some cases, which otherwise can cause a solution to blow up early on. Also - what are the metrics on your grid? How clean is the grid? It's very common that the root of all problems is the grid (followed closely by the boundaries!).

Good luck.

[KTG]

About the epsilon: yes, lets assume it is correct for now. I was thinking maybe there was a yplus issue, but have tried a number of mesh sizes and refinement doesn't seem to help. It occurred to me after getting convergence with the upwind schemes that I could just set the initial conditions to k and epsilon to something much closer to the final result. This helped a lot with avoiding blow ups. Its generally not a problem now, more I just get residuals that hit a sort of steady state, like the ones attached.



I have tried these schemes, but they don't seem to make much of a difference:

laplacianSchemes
{
//    default         Gauss linear corrected;
    default         Gauss linear limited corrected 0.5;
}
snGradSchemes
{
    default         limited .5;
}

In terms of the limiters, I tried .5 and .3, but I didn't really change much.




Here are some stats on my latest mesh:

Mesh stats
    points:           1269247
    faces:            3367790
    internal faces:   3171472
    cells:            1055701
    faces per cell:   6.19424
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     852998
    prisms:        65594
    wedges:        0
    pyramids:      0
    tet wedges:    250
    tetrahedra:    0
    polyhedra:     136859
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   12701
            5   7091
            6   24046
            7   26477
            8   11479
            9   29448
           10   9
           11   2
           12   21554
           13   2
           15   4040
           18   10

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points                  Surface topology
                  outlet      867      973  ok (non-closed singly connected)
                   inlet      867      973  ok (non-closed singly connected)
                   walls   193703   207155  ok (non-closed singly connected)
            tube_extract      881     2517  ok (non-closed singly connected)

Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking geometry...
    Overall domain bounding box (0.00766238 0.00763992 0.0287251) (0.0279774 0.05 0.07)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.01809e-16 3.31283e-16 -3.81013e-16) OK.
    Max cell openness = 4.49581e-16 OK.
    Max aspect ratio = 26.8693 OK.
    Minimum face area = 2.44844e-10. Maximum face area = 2.61156e-07.  Face area magnitudes OK.
    Min volume = 1.11115e-14. Max volume = 1.29617e-10.  Total volume = 7.18031e-06.  Cell volumes OK.
    Mesh non-orthogonality Max: 60.5062 average: 9.93102
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.7303 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

I have tried playing with the refinement a lot, but maybe you see something in the mesh quality I should try changing? I am actually not sure I am being too conservative with the mesh? I only have about 12 cells across the small channel. The wall spacing seemed to be about right, but maybe I am wrong... 12 cells across for a freestream of 86 m/s seems very small to me, but I never really use these turbulence models...



For the boundaries, I am using the classic "wind tunnel" type conditions (tube_extract is just a section of the walls where I was trying additional refinement):

U:

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 -.527 0); //uniform (0 -.08 0);
    }

/*
    {
        type   flowRateInletVelocity;
        volumetricFlowRate constant 0.00016667;
    }
*/
    outlet
    {
        type            zeroGradient;
    }

    "(walls|tube_extract)"
    {
        type            noSlip;
    }

}

p:

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

"(walls|tube_extract)"
    {
        type            zeroGradient;
    }

}

nut:

boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
"(walls|tube_extract)"
{
        type            nutkWallFunction;
        value           uniform 0;
    }

}

k:

internalField   uniform 100;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform .000438;
    }
    outlet
    {
        type            zeroGradient;
    }
"(walls|tube_extract)"    {
        type            kqRWallFunction;
        value           $internalField;
    }

}

epsilon

internalField   uniform 10000000;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform .00015;
    }
    outlet
    {
        type            zeroGradient;
    }
"(walls|tube_extract)"    {
        type            epsilonWallFunction;
        value           $internalField;
    }

}

If that generates any new ideas let me know - thanks!