CFD Online Discussion Forums - Sudden contration

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Sudden contration (https://www.cfd-online.com/Forums/openfoam-solving/207714-sudden-contration.html)

Sudden contration

Hello

I'm trying to learn to use OpenFOAM. For that purpose I tried to simulate a sudden contration of a circular pipe. I go from 80 mm to 40 mm see image below (clipped trough the half), where the inlet is the 80 mm.

Attachment 65803

When i run it runs fine, but when i check the results in paraview I get a pressure drop of 50. Converting it to bar using the density of water of 998, I get a pressure drop of 0.5 bar. When I try to calculate the pressure drop in a manual way (fx. using http://www.pressure-drop.mobi/0303.html) i get a pressure drop of 0.13 bar.

Changing the inlet velocity I always get roughly a factor of 4 between the two methods. I suspect that I'm doing something completely wrong in my simulation, but I don't know what it is.

I know one of the difficult things with CFD is to choose sensible parameters, so maybe that is where i'm going wrong. Any help is much appreciated.

Iøve created teh mesh with salome and used ideasUnvToFoam to create the mesh in OpenFOAM, you can see my mesh here (I know it is a bit coarse, but I thought it would still give a sensible result):

Attachment 65805

I've tried to setup the simpleFoam solver with the following parameters:

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volVectorField;

    object      U;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 1 -1 0 0 0 0];



internalField   uniform (0 0 0);



boundaryField

{

    channel_closed_solid_inlet

    {

        type            fixedValue;

        value           uniform (0 2 0);

    }



    channel_closed_solid_outlet

    {

        type            zeroGradient;

    }



    channel_closed_solid_walls

    {

        type            fixedValue;

                value                        uniform (0 0 0);

    }

}



// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    object      p;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -2 0 0 0 0];



internalField   uniform 0;



boundaryField

{

    channel_closed_solid_inlet

    {

        type            zeroGradient;

    }



    channel_closed_solid_outlet

    {

        type            fixedValue;

        value           uniform 0;

    }



    channel_closed_solid_walls

    {

        type            zeroGradient;

    }

}



// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    location    "0";

    object      k;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -2 0 0 0 0];



internalField   uniform 0.024576;



boundaryField

{

    channel_closed_solid_inlet

    {

        type            fixedValue;

        value           uniform 0.024576;

    }

    channel_closed_solid_outlet

    {

        type            zeroGradient;

    }

    channel_closed_solid_walls

    {

        type            kqRWallFunction;

        value           uniform 0.024576;

    }

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    location    "0";

    object      epsilon;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -3 0 0 0 0];



internalField   uniform 0.113;



boundaryField

{

    channel_closed_solid_inlet

    {

        type            fixedValue;

        value           uniform 0.113;

    }

    channel_closed_solid_outlet

    {

        type            zeroGradient;

    }

    channel_closed_solid_walls

    {

        type            epsilonWallFunction;

        value           uniform 0.113;

    }

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    location    "0";

    object      nut;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -1 0 0 0 0];



internalField   uniform 0;



boundaryField

{

    channel_closed_solid_inlet

    {

        type            calculated;

        value           uniform 0;

    }

    channel_closed_solid_outlet

    {

        type            calculated;

        value           uniform 0;

    }

    channel_closed_solid_walls

    {

        type            nutkWallFunction;

        value           uniform 0;

    }

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "constant";

    object      turbulenceProperties;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



simulationType RAS;



RAS

{

    // Tested with kEpsilon, realizableKE, kOmega, kOmegaSST, v2f,

    // ShihQuadraticKE, LienCubicKE.

    RASModel        kEpsilon;



    turbulence      on;



    printCoeffs     on;

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "constant";

    object      transportProperties;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



transportModel  Newtonian;



nu              nu [0 2 -1 0 0 0 0] 1e-6;



CrossPowerLawCoeffs

{

        nu0                nu0 [0 2 -1 0 0 0 0] 1e-6;

        nuInf        nuInf [0 2 -1 0 0 0 0 ] 1e-6;

        m                m [0 0 1 0 0 0 0] 1;

        n                n [0 0 0 0 0 0 0] 1;

}



BirdCarreuCoeffs

{

        nu0                nu0 [0 2 -1 0 0 0 0] 1e-6;

        nuInf        nuInf [0 2 -1 0 0 0 0 ] 1e-6;

        k                k [0 0 1 0 0 0 0] 0;

        n                n [0 0 0 0 0 0 0] 1;

}



// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "system";

    object      controlDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



application     simpleFoam;



startFrom       latestTime;



startTime       0;



stopAt          endTime;



endTime         2000;



deltaT          1;



writeControl    timeStep;



writeInterval   100;



purgeWrite      0;



writeFormat     ascii;



writePrecision  6;



writeCompression off;



timeFormat      general;



timePrecision   6;



runTimeModifiable true;





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "system";

    object      fvSchemes;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



ddtSchemes

{

    default         steadyState;

}



gradSchemes

{

    default         Gauss linear;

}



divSchemes

{

    default                 none;

    div(phi,U)              Gauss upwind;

    div(phi,k)              Gauss upwind;

    div(phi,epsilon)         Gauss upwind;

    div(phi,R)                     Gauss upwind;

    div(R)                                Gauss linear;

        div(phi,nuTilda)        Gauss upwind;

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

}



laplacianSchemes

{

    default         Gauss linear corrected;

}



interpolationSchemes

{

    default         linear;

}



snGradSchemes

{

    default         corrected;

}



wallDist

{

    method meshWave;

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "system";

    object      fvSolution;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



solvers

{

    p

    {

        solver          GAMG;

        tolerance       1e-06;

        relTol          1e-06;

        smoother        GaussSeidel;

                nPreSweeps                0;

                nPostSweeps                2;

                cacheAgglomeration on;

                nCellsInCoarsestLevel 10;

                mergeLevelse        1;

    }



    "(U|k|epsilon|R|nuTilda)"

    {

        solver          smoothSolver;

        smoother        symGaussSeidel;

        tolerance       1e-06;

        relTol          1e-06;

    }

}



SIMPLE

{

    nNonOrthogonalCorrectors 0;

}



relaxationFactors

{

    fields

        {

                p        0.3;

        }

                

        equations

        {

                U                 0.3;

                k                 0.3;

                epsilon 0.3;

                R                 0.3;

                nuTilda 0.3;

        }

}





// ************************************************************************* //

Check whether your y+ values at the walls are in the appropriate range. Being off can cause a faulty prediction of the pressure loss. Alternatively, you could try running the case using the kOmega model.

Furthermore, the mesh is not all that beautiful.

Thanks for the suggestions. I checked the yPlus values (didn't know about them, thanks for that), and they are way to high. Also I agree that the mesh is not good.

I use salome to mesh, and I'm also new using that. I tried to make the mesh better. What do you think about this new mesh? Looks better?

Attachment 65836

If a model allows manual meshing you should go that way. In you case it is straightforward.

If you look closely at the meshes the have strange regions both.

But I don't believe that your problem depends on that.

~

Which Re number did you use? May be there is a problem with the turbulence calculation. Pressure drops differ widely between laminar and turbulent flow.

The Reynolds number is very high. I think in the region of 200000 - 300000. You don't mean that it should be set somewhere in OpenFOAM do you? Or do you mean when I calculated the k and epsilon values?

Not sure what you mean by manual mesh? How do I do that?

Should I actually calculate different k and epsilon for the pipe walls with the large diameter and the small diameter?

> Reynolds number is very high...200000 - 300000

I recommend starting with lower Re numbers and looking what happens.

> what you mean by manual mesh? How do I do that?

blockMesh.

The reason I used Salome to do the mesh, was because I'm investigation how a OpenFOAM workflow can be implemented at our company, and in general the geometries will be very complex.

I tried going down to a speed of 0.1 m/s, and a Re below 20000.

I still get a factor of 4 between the manual calculation and OpenFOAM.

I think my convergence and yPlus is OK?

Code:

smoothSolver:  Solving for Ux, Initial residual = 1.22094e-07, Final residual = 1.22094e-07, No Iterations 0

smoothSolver:  Solving for Uy, Initial residual = 1.6853e-08, Final residual = 1.6853e-08, No Iterations 0

smoothSolver:  Solving for Uz, Initial residual = 1.25779e-07, Final residual = 1.25779e-07, No Iterations 0

GAMG:  Solving for p, Initial residual = 8.73275e-07, Final residual = 8.73275e-07, No Iterations 0

GAMG:  Solving for p, Initial residual = 8.73275e-07, Final residual = 8.73275e-07, No Iterations 0

time step continuity errors : sum local = 2.53288e-06, global = 5.52101e-08, cumulative = -6.2634e-06

smoothSolver:  Solving for epsilon, Initial residual = 7.9535e-07, Final residual = 7.9535e-07, No Iterations 0

smoothSolver:  Solving for k, Initial residual = 9.04923e-07, Final residual = 9.04923e-07, No Iterations 0

ExecutionTime = 264.54 s  ClockTime = 265 s

Code:

patch channel_closed_solid_walls_large y+ : min = 0.482836, max = 14.916, average = 4.53306

    patch channel_closed_solid_walls_small y+ : min = 6.79218, max = 48.9457, average = 23.8934

> ecause I'm investigation how a OpenFOAM workflow can be implemented at our company

I don't think that it works that way. OF (and other cfd software) is not some kind of advanced pocket computer. You need experience how to work with it. This includes a strategy to ensure the correctness of your results.

Yeah, I agree. This is a preliminary study to see how one can work with OpenFOAM. I think it is possible to study the workflow. That is how will we take a given component from design, to mesh, to solving, and if this is a viable way of working (or we need to invest in expensive commercial programs). I think my conclusion is that it is definetely possible to use OpenFOAM, and I like the fact that it is open source and c++, because as a physicists and c++ programmer, I see the possibility to program special cases important for our company in the future.

If we see that that is possible, and we find that it is worth investing money and time in it, we will start getting help from the outside, in form of courses and support, to gain more experience in these kind of simulations. I'm fully aware that it is not hard to get OpenFOAM up and running, but it is very difficult to make correct assumptions and interpretations of the results. Something we have to invest a lot of time in to be able to do.

My question here was more if someone could see why I didn't get my results to match, my simple calculations made by empiric formulas (formulas we have seen match experiments very well). If it is not possible to do so, and we need more training and experience to do it correctly, we will take that into considerations. Justed hoped that someone maybe had experience with simulating a contraction, and knew how to simulate it, or knew if precise results could be expected.

It could be that with CFD you can not expect to get results more precise than within a factor of 5, and that my simulation was succesfull (hope not :)). If that is the case that is also fine, but then we need to take that into considerations.

Hello,

did you try the have prism layers close to the wall? Usually tetrahedra are too diffusive to accurately solve the wall shear stress close to the wall. Besides this they introduce an additional diffusion due to the fact that the flow is not normal to the cell face.

So I strongly advise to use a high quality mesh close to the wall to accurately resolve the transport phenomena which are leading to the two recirculation regions in front upwind and downwind of the contraction.

In addition it is best practice (if the computational costs allows it) to have a y+ ~ 1 at the wall in order to avoid the usage of wall function which are known to not correctly reproduce separated flows.

Best

Michael

Hello

If I don't use wallfunctions what should i use for type in the k and epsilon files?

> It could be that with CFD you can not expect to get results more precise than within a factor of 5

You may expect much better results.

Start with a simple geometry, a simple physics (low to medium Re) and an straightforward mesh. Look what happens and change one variable after the other.

You try
- without experience
- using a mesh of doubtful quality
- consisting of tetraeder elements (instead of heaxaeder)
- to calculate an example withe high Reynolds number.

Too much complications for the start.

Quote:

Originally Posted by mAlletto (Post 708691)

What he said... And to reach that goal while still starting from Salome. I would recommend that you take a look at the most excellent youtube series content posted by Tobi (https://www.youtube.com/user/HolzmannCFD/videos). The manifold videos by Tobi basically outlines the workflow that you are looking into. In short it will teach you how to triangulate your geometry and then make a nice hex-mesh by using SnappyHexMesh (cfMesh could also be used of course)

Thanks for the great link. I was actually looking into using hex mesh. I guess the work flow will then be to prepare stl files in salome, and then use snappyHexMesh to create the mesh.

Hi Karlsen,
thank you very much for mentioning my screencast in that way. I feel honored. An additional information to the screencasts:

Quote:

The training videos offer a broad and intensive insight into the OpenFOAM® framework. The videos will guide you through all different kinds of OpenFOAM® related topics such as »Design preparation«, »Surface Triangulation«, »Meshing with snappyHexMesh«, »Setting up the Numerical Case for OpenFOAM®«, »Boundary Conditions« as well as »Advanced Boundary Conditions«. During the videos, plenty of tricks and hacks are given which one allows understanding the toolbox better while user mistakes can be significantly reduced. Learning to use the command line will improve your Linux level, increase your speed drastically and allow you to use the powerful scripting languages bash / csh / python in combination with OpenFOAM® such as combining OpenFOAM® and DAKOTA®.

More than 4.000 people all over the world evaluated the screencast already. Based on this lecture, these people could run their projects in OpenFOAM® and increased their skills drastically.

The screencast was recorded for OpenFOAM-4.x. Based on ongoing developments, code maintenance, sustainability and performance improvements, classes and functions can be named differently in newer releases. If you have any problems related to the videos and the newer releases, do not hesitate to write Holzmann CFD your inquiry (advancedMaterial@Holzmann-cfd.de).

Addition hints to the problem you are solving:

The boundary layer has to establish first, using a fixedValue will not give you a boundary layer. Thus, this layer has to be established first (I am not sure if this will influence the results drastically). However, it is worth to know that your inlet boundary condition can influence the phenomenon you are interested in; a rule of thumb is, e.g., for the pipe diameter of 100 mm, the pipe inlet - before the contraction - should be at least 300 mm away.
Related to the hint mentioned above, a parabolic inlet profile would be recommended.
Change the inlet values for k and epsilon to turbulentKinetic... and mixingLength....
The outlet should have a characteristic length to avoid influence to the internal field. E.g., if the outlet is 10 mm behind the contraction, it will influence the whole behavior of the system based on the fact, that the outlet pipe (in reality) will not stop here.
During my master study, we made such tests too. However, after the contraction, we had the expansion immediately and recorded the pressure drop before and behind the contraction. The comparison to the analytical solution was 99.99%. Therefore, I expect that there is something wrong in your set-up.
For Tetraeder meshes, the numerical schemes could be changed.
Furthermore, you should be aware of your first order scheme.
Summing up, there are still plenty of things you can consider improving your simulation. However, a discrepancy of 2.5 bar is a bit too high.
What comes into my mind, please be aware of the pressure you set. It is not the real pressure. The differences $\Delta p$ are fine but absolute values are different.

Good luck & as Karlsen mentioned, hex-dominant meshes are much better in a numerical sense. E.g., second order linear scheme and tet meshes for pure convection: https://www.youtube.com/watch?v=C0CcN7l37Fg
Additionally, as Gerhard said, the turbulence model can influence your solution (I am not a turbulence guy; therefore I cannot give you any feedback).

A few hints for open source software in general (if one is interested - valid for everybody) - the above content is not related to the topic
================================================== =======================================
Furthermore, I want to make one statement to OpenFOAM and Open Source software that many people do not or don´t want to recognize.
I published my material for more than five years for free till mid of 2018. Therefore I know the following circumstances:

Based on the continuous increase of my content, I was not able to spend more and more of my private time keeping the material up to date.
Sure, the community was always happy about my work and especially the screencasts and my book were successful, which was an honor to me.
However, the time I spend (several thousand hours) for preparing, re-organizing and re-building the material was not in relation to what I got back.
Summing up, I got 19 voluntary donations; most of them were between 5€ to 10€ - that is totally fine and very welcomed. But 19 people who give value to my free material compared to the statistics - the book was downloaded more than 20.000 times; screencasts were watched more than 6.000 times; tutorials were downloaded (summing all 40 single cases) more than 50.000 times - it is the infinite element, no one can see.
Moreover, considering my fix costs for server stuff and so one ... ... ... - yep, I provided material for free and I pay for that :D

In other words, doing a CFD analysis with 19 cells will not make sense at all. However, if you have 76.000 cells, you can do pretty cool things, don´t you? Okay, some people think they need 600.000.000 cells for their problems :D

It sounds a bit of advertising to donate to OpenFOAM, to my work or other open source software but it is not. I want to make you sensitive, that it is not naturally that open source software is free forever.

My wish for the future would be: If people are using open source software, they should give something back regarding:

own developments
bug reports
contributing to that software
or donating an amount of money you think is worth

An excellent example that people thought a natural free material is forever free is my material - it´s closed now (most of the stuff); You don´t believe how many emails I got, questioning when it is available again :)

================================================== =======================================

Enough of my talk. I hope that the hints at the beginning help to resolve the problem. Another option is that you post your case that we can check out the stuff in more detail. Additionally, the equation and reference to the equation would be beautiful to have; I expect that you use simple Bernoulli equation?

Further information for OpenFOAM

wiki.openfoam.com
openfoamwiki.net
József Nagy YouTube Channel
The blog of our Chine colleague
Chalmers CFD courses
OpenFOAM Workshop material
PFAU Meetings

There are many more, but these sources exist for a long time and
new things are added from time to time.

Thank you so much for taking your time to help me with such a detailed answer.

Just a couple of questions:

The boundary layer has to establish first, using a fixedValue will not give you a boundary layer. Thus, this layer has to be established first (I am not sure if this will influence the results drastically). However, it is worth to know that your inlet boundary condition can influence the phenomenon you are interested in; a rule of thumb is, e.g., for the pipe diameter of 100 mm, the pipe inlet - before the contraction - should be at least 300 mm away.
Related to the hint mentioned above, a parabolic inlet profile would be recommended.

I will try to extend the length of the pipe before and after the contraction to give time to make the flow develop. Not sure about the fixedValue thing, should I change the boundary condition, or will the extended pipe length with the fixedValue boundary condition, develop a boundary layer?

Change the inlet values for k and epsilon to turbulentKinetic... and mixingLength....

I will try this.

During my master study, we made such tests too. However, after the contraction, we had the expansion immediately and recorded the pressure drop before and behind the contraction. The comparison to the analytical solution was 99.99%. Therefore, I expect that there is something wrong in your set-up.

Great news :)

Furthermore, you should be aware of your first order scheme.

Not sure what you mean by this?

Summing up, there are still plenty of things you can consider improving your simulation. However, a discrepancy of 2.5 bar is a bit too high.
What comes into my mind, please be aware of the pressure you set. It is not the real pressure. The differences \delta p are fine but absolute values are different.

As far as I understand it is the pressure divided by density? Or do you mean that the pressures are not the absolute pressure, because I set my reference at the outlet to zero?

Quote:

Originally Posted by MrAndersDk (Post 708694)

Hello

If I don't use wallfunctions what should i use for type in the k and epsilon files?

you can find a good discussion about wall function in this thread:

https://www.cfd-online.com/Forums/op...ion-usage.html

Quote:

Originally Posted by MrAndersDk (Post 708791)

Furthermore, you should be aware of your first order scheme.

Not sure what you mean by this?

To calculate different quantities such as the pressure, or fluxes at the faces, we need discretization schemes (system/fvSchemes). You can adjust them according to your needs.

Quote:

Summing up, there are still plenty of things you can consider improving your simulation. However, a discrepancy of 2.5 bar is a bit too high.
What comes into my mind, please be aware of the pressure you set. It is not the real pressure. The differences $\Delta p$ are , but absolute values are different.

As far as I understand it is the pressure divided by density? Or do you mean that the pressures are not the absolute pressure because I set my reference at the outlet to zero?

Correct. The pressure difference $\delta p$ can be compared directly, but if you compare the pressure value with another pressure (Pa), you have to multiply it with the density.

Thanks for all the suggestions.

I've tried to use the suggestions. First of all I've tried learning some snappyHexMesh to genererate a better hexagonal mesh, and elongated the pipe length. I now have the following mesh:

Attachment 65932

Attachment 65933

Attachment 65934

As you can see the grid is much finer especially at the walls. However, there are still some strange cell's, don't know if that is a problem. The finer mesh is generating more of the expected turbulence in the corners and just after the contraction. However, i still get a too high pressure drop over the contraction compared to the empirical formula.

Attachment 65935

As you can see the pressure drop is between 3.8 to 3.2 coresponding to 32 mbar and 38 mbar for water. The empirical formula suggest 8 mbar for a flow of 0.5 m/s.

I've also tried using the suggestied boundary conditions

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    location    "0";

    object      epsilon;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -3 0 0 0 0];



internalField   uniform 0.00031185580835531597;



boundaryField

{

    inlet

    {

        type            turbulentMixingLengthDissipationRateInlet;

        mixingLength    0.05;

                value           1;

    }

    outlet

    {

        type            zeroGradient;

    }

    wall_large

    {

                //type                          zeroGradient;

        type            epsilonWallFunction;

        value           uniform 0.0000031185580835531597;

    }

        wall_small

    {

                //type                          zeroGradient;

        type            epsilonWallFunction;

        value           uniform 0.000008088538952456775;

    }

}





// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |

|  \\    /   O peration     | Version:  v1806                                 |

|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       volScalarField;

    location    "0";

    object      k;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



dimensions      [0 2 -2 0 0 0 0];



internalField   uniform 0.000040603056516117657;



boundaryField

{

    inlet

    {

        type            turbulentIntensityKineticEnergyInlet;

        intensity           0.05;

                value               1;

    }

    outlet

    {

        type            zeroGradient;

    }

    wall_large

    {

                //type                          zeroGradient;

        type            kqRWallFunction;

        value           uniform 0.000040603056516117657;

    }

        wall_small

    {

                //type                          zeroGradient;

        type            kqRWallFunction;

        value           uniform 0.00004828544369982472;

    }

}





// ************************************************************************* //

maybe the wall functions are still wrong? Also tried the zero gradient condition on the walls, but that didn't seem to help.

So all in all I like the new grid, and the ability to see turbulence effects with flow lines, but the pressure drop is still just as much off.

More suggestions how to improve my simulation?

Please provide your Test case as well as the empirical equation and the limits of the equation.

You should avoid excesive cell volume jumps close to the wall. In my experience introducing some rounding at the edges of the contraction reduces the pressure loss. What are realistic radii achieved in the production?

Hello

my files are here. I've tried many different grids and schemes but all seems to reproduce roughly the same result. It could be I should try a small rounding of the transition and corners to make the geometry more realistic.

Attachment 65952

The empirical relation comes from the book Pipe Flow - A practical and Comprehensive guide by Rennels and Hudson. It is given by

$dp = \frac{1}{2} \rho v^2 K_2$

where

is the velocity in the small part of the pipe.

is the loss constant, which for a diameter ratio of 0.5 is roughly 0.5. This formula also agrees well with this webpage calculator:

http://www.pressure-drop.mobi/0303.html

I've estimated k and epsilon from the following page

http://ichrome.com/blogs/archives/342

and also tried other formulas. I've also tried running the simulation without wall functions and used the zero gradient boundary condition.

Well, after your equation it is clear that you can satisfy your results based on Bernoulli. However, I will check the case today in the evening and will report back.

Edit. The geometry is missing.

Here is a link to all the files:

https://drive.google.com/file/d/11-w...ew?usp=sharing

My problem is that K_2 should be around 0.5, but my simulation corresponds to more like K_2 = 1.5.

I really appreciate you taking your time to help

Hi all,

I should do other things than investigating in that case:). However, I made it, and the result is as follow:

Diameter inlet 0.08 m
Diameter outlet 0.04 m
Velocity inlet 0.5 m/s
Velocity outlet (unknown)

Based on mass conservation, we can simply say $\dot V_1 = \dot V_2$ (density is constant). Therefore, the velocity at the outlet can be estimated to be around 2 m/s (~ Re=80000). So the flow field is very turbulent.

Also, we can take the Bernoulli equation that states:

$p_1 + \frac{1}{2} \rho |\textbf{U}_1|^2 = p_2 + \frac{1}{2} \rho |\textbf{U}_2|^2$

We know:

velocity inlet (1) is 0.05 m/s
velocity outlet (2) is 2 m/s
pressure outlet (2) is 0 [Pa/density] (simulation / incompressible)
Keep in mind there is no pressure of 0 Pa. We can assume to set 1000hPa at the outlet which does not influence the numerics

Therefore, the outlet pressure can be calculated to be:

$p_2 = p_1 + \frac{1}{2} \rho |\textbf{U}_1|^2 - \frac{1}{2} \rho |\textbf{U}_2|^2$

This is equal to 98125 Pa. The difference is 1875 Pa (= 18.75 mbar). In that calculation, we do not have any dissipation effects caused by vortexes/turbulences. Thus, as we already know that the Reynolds number is high, dissipation effects take place, and the pressure drop will rise. Therefore, the $\Delta p$ should be higher than 18 mbar and not lower.

I changed your mesh completely and made the simulation with almost 80% fewer cells. Based on the high Reynolds number, the steady-state simulation cannot be reached (fluctuations). In your fine mesh, it has to be much worse. However, the average value of the pressure at the outlet is:

Code:

Executing functionObjects

surfaceFieldValue patchAverage(name=inlet,p) write:

    areaAverage(inlet) of p = 3.22968

The inlet patch is set to 0 so the difference in pressure $\Delta p$ is 3.22968 Pa / kg / m^3. After multiplying with the density of the fluid, we get around 3229 Pa of pressure difference which corresponds to 32 mbar. This correlates to the above-given remarks; at least that it has to be higher. The equation you showed is correct but can you specify why K should be 0.5? In my opinion, you can use the equation only from inlet to the compression (pipe 1) - pressure loss inside that pipe - and after that inside the second pipe. As you should know, the K values include the length and diameter of the pipe as well as the roughness of the pipe (friction).

Assuming the pipe to be very smooth, we might use the simple formulation of Blasius:

$\lambda = \frac{0.3164}{Re^{0.25}}$

The following calculation is related to the case attached.

Pipe #1 (from the inlet to compression)

Length is 0.3 m
Diameter 0.08 m
Velocity 0.5 m/s
Re = 40000
Lambda = 0.022372

$K_1 = \lambda_1 \frac{l_1}{d_1} = 0.0223728 \frac{0.3}{0.08} = 0.084$

The pressure lose within the first pipe section is:

$\Delta p_1 = \frac{1}{2}\rho v_1^2 K_1 = 10.48 Pa$

Pipe #2 (from compression to outlet)

Length is 0.5 m
Diameter 0.04 m
Velocity 2 m/s
Re = 80000
Lambda = 0.0188132

$K_2 = \lambda_2 \frac{l_2}{d_2} = 0.0188132 \frac{0.5}{0.04} = 0.23516$

The pressure loose within the second pipe section is:

$\Delta p_2 = \frac{1}{2}\rho v_2^2 K_2 = 470.33 Pa$

Thus, the pressure loose based on the roughness inside the pipe is around 480 Pa.

The 480 Pa are based on the pipe system and its length. As already stated, based on Bernoulli we will get an automatic pressure drop which is around 1875 Pa.So we end up in total with 23 mbar pressure loose. I guess - but I am not sure - that the K value, respectively the $\lambda$ value takes care of turbulence dissipation inside the pipe too (becasue the Reynolds number is included). So as you can see from the K values calculated based on Blasius, we get complete different ranges for the pressure drop in the two pipe sections. I would trust that the pressure drop is around 30 mbar here.

I already spend 3 hours for that topic now and to make it a fundamental work, you should:

Choose a roughness of your pipe
Estimate the lambda value for each pipe section
Estimate the K value of each pipe section
Use Bernoulli to get the pressure drop based on continuity
Add the estimated pipe-pressure drops
Recalculate with OpenFOAM and appropriate nut wall functions (roughness)
Compare the results
Realize that it fits well but not 100 %

Keep in mind that the pressure drop in the CFD simulation is more complicated as one might think at the beginning. Friction at the pipe wall (nut wall function), dissipation due to turbulence and recirculation, and bla blub :) even though this example is very simple, it demonstrates the complexity of fluid dynamics and the knowledge one should have to be really good within that topic.

Hopefully your question is solved now.
By the way. As you might realize, that the mesh sizes are entirely different (mine has around 20.000 cells) the results are more or less in the same area ~ 30 mbar. There would be other things I would like to say, such as the turbulence resolution, RANS, ... but well - I go to bed now. By the way, I got 30 mbar pressure drop via laminar calculation too.

Case Download OpenFOAM Foundation 6.x

Tobias you are too kind, I've never expected you to do all that. I'm amazed by this communitys, will to help me, even though I know almost nothing, and asks stupid questions. Not just you Tobias but all who have contributed to my question. I could understand from your first post, that you are doing all of this for free, even though it is actually your profession, and you didn't have to. I've made a small contribution to your cause, because people like you are invaluable for open source projects to survive.

Now back to my question :). I've actually also ran the simulation on a much coarser grid and did it laminar, and got the same result roughly as you point out. The reason I kept refining my grid was because of my ignorance, and that people pointed out it could be my resolution of the grid around the wall that was not good enough.

Now I actually see that I was actually so blinded by all this CFD, that I forgot the most fundamental phsyics, and lost sight at was I was actually calculating in OpenFOAM.

So just to recap, and test my understanding (please correct me if I'm wrong, or confirm my understanding).

If we call the first segment (large pipe) for 1 and the last segment for 2, then the following equation must be satisfied:

$p_1 + \frac{1}{2} \rho v_1^2 - \Delta p = p_2 + \frac{1}{2} \rho v_2^2$

where $\Delta p$ is the static pressure loss due to the contraction (neglecting the pipe friction loss, because it is expected to be minor compared to the contraction). It is this $\Delta p$ my formula is actually calculating, and not

as I stupidly have done. So isolating $\Delta p$ in my formula above i get

$\Delta p = p_1 - p_2 + \frac{1}{2} \rho (v_1^2 - v_2^2)$

inserting the values from the simulation for p_1 and p_2 and the velocities v_1 = 0.5 and v_2 = 2.0 I get

$\Delta p = 30 mbar - 0 mbar + \frac{1}{2} \rho (0.5^2 - 2.0^2) = 11 mbar$

This is acceptable close to the 8 mbar estimated from the empiric formulas, taking into account that such empiric formulas is also very uncertain. Maybe one could easily argue that the 3 mbar, that the simulation overestimates $\Delta p$ stems from the surface friction which my formula doesn't take into account.

So to recap further, I'm an idiot :). However, even though I did something stupid, all of your input has pushed me further, and my understanding of OpenFOAM has somewhat increased.

So how should I proceed to become better at CFD (not at OpenFOAM. I think I'm actually convinced that I can learn OpenFOAM by videos, tutorials and this forum. Even though I find it difficult to understand how to make qualified decisions on grids, schemes and so on, this problem has convinced me that my focus probably not should be so much on OpenFOAMs quirky details). Is it possible to get better at CFD without courses or a full education? I've found books on general CFD by Wilcox, Munson, Anderson, etc. is that the way to go?

I have a PhD in other areas of physics (yeah I know, who would have thought so...), and have taking math and c++ on a high level, so I'm quite proficient in learning new complex skills, and think I have a solid foundation. But I'm unsure of the right path to becoming proficient in CFD. Any advice or guidance in this area, is much appreciated.

by the way. I'm also going to bed now, but I will dissect your solution Tobias tomorrow, to steal any ideas of how you set it up.

Hi,
first of all, thank you for the kind words and your feedback & recap. Based on my lack of knowledge and almost rusty fluid dynamics lecture, I remember some equation you mention (Bernoulli + dp), but I cannot find my material right now. However, a short search gave me the following link (http://www.thermopedia.com/content/659/) where a comparison/equation is provided to calculate the pressure drop (it seems to be a bit different to yours) for cases we are looking at. One thing is clear at all: The sudden jump in terms of the pipe diameter should give a pressure drop based on viscous shearing/friction $\Delta p$ .

Summing up. I refreshed my knowledge of the topic and might do some further investigations with the roughness of a pipe, respectively the nut wall functions - roughness (if there is any time). However, as I realized, there are still plenty of knowledge gaps I should reconsider to investigate some day.

To your question of CFD information. I read a few books such as Bird et al. (Transport Phenomena), Wilcox (Turbulence), Ferziger et al. (Numerical ...), Moukalled et al., the Ph.D. thesis of Jasak Hrvjoe or Philipp Cardiff and many more. Personally, some literature was harder to read, others more straightforward. A summary of different books/papers/thesis is given in my book (still free on research gate). However, considering my book, there are still a lot of things to improve, and it might be somehow written in a 'different' way - I don´t know.
The good thing nowadays is, that there is plenty of information out for CFD in general and OpenFOAM. E.g., the wiki.openfoam.com site or the courses of Chalmers :)

As I am not intelligent and smart as other genius people out there, I always have to derive the equations on my own. After that, commonly I get the point - but not always :), and in some cases, I struggle to figure it out at all. Therefore, picking up your statement 'I am an idiot' does not hold and is disproved.

To get an overview of the influence of numerical schemes you can check out my voluntary part of my website voluntary.holzmann-cfd.de. However, this is just an overview :)

thanks for the reply.

I think the equation you linked to are the same as mine, rewritting for mass flow instead of velocity, and expressing v_1 through v_2 using the area ratios (constant volume flow).

Okay, I was not sure because of the more complex formulation. As you said, K is around 0.5, I expected it is the S value of the equation, and thus, the C value is somehow missing. However, a simple calculation of the formula will give the answer.

I actually think they write that S is the area ratio of the two pipe section. C is K or at least closely related to it.