[highpressuretube]

Dear all,

I am still having troubles simulating a throttle. I have gone from Salome tetraeder elements to blockmesh hex 2D, than blockmesh hex 3D. Meshcheck says ok, yplus goes up to 60 but I read for k-epsilon modellisation 30-60 should still be ok. But still:

The expected pressure drop disappears!
In order to reach reasonably fast convergence I started with low RE model and augmented the reynolds number by increasing the kinematic viscosity after reaching convergence until I got to the measured viscosity I have from experiments. (Thats where the bumps in the residual plot come from, each time viscosity is multiplicated by factor 10).

As long as RE is indicating laminar flow, behaviour is as expected. When it comes to highly turbulent flow, I don't understand what happens. Especially since measurements indicate an expected specific pressure drop of 3000 m^2/s^2 and not zero. :-( I am out of ideas. Please help me understand.

Best regards,
highpressuretube

[anishtain4]

What kind of turbulence model are you using? what is the wall model? Are you using turbulence model for the laminar cases too?
What do you mean by "Augmented the reynolds number by increasing the kinematic viscosity"? Increasing viscosity decreases the reynolds number!

[highpressuretube]

Sorry for the late answer - was on vacation :-) But unfortunately the problem didn't disappear yet.

So turbulence model is k-epsilon.
I use turbulence model for the laminar case too an the result are similar to laminar model only.

The confusion about augmenting viscosity was because I augmented the power of the negative exponent from 2.255x10E-2 to 2.255x10E-6 so of course I was decreasing viscosity. All plausible there.

What exactly do you mean by what is the wall model?

[highpressuretube]

I have plotted k and epsilon trying to find indications as to where is the error. Any ideas?

[arvindpj]

In simple foam, the pressure is divided by density. Use the same scales for all 4 cases. I do see pressure drop in the Re=30000 case as well.

[highpressuretube]

But how can I have zero pressure at the inlet?

I enter with 9,4314m/s at the inlet, which should correspond to about 4l/min because of inlet diameter 3mm, right?

At 4l/min we measured a pressure drop of around 24bars between simulated inlet and outlet. With density of 823g/cm^3, I calculate an expected pressure drop of just about 2950, don't I?

This is nothing close to what I see in the calculation, where inlet pressure, thus delta pressure is zero.

Regards,
highpressuretube

[arvindpj]

Originally Posted by highpressuretube

But how can I have zero pressure at the inlet?

I enter with 9,4314m/s at the inlet, which should correspond to about 4l/min because of inlet diameter 3mm, right?

At 4l/min we measured a pressure drop of around 24bars between simulated inlet and outlet. With density of 823g/cm^3, I calculate an expected pressure drop of just about 2950, don't I?

This is nothing close to what I see in the calculation, where inlet pressure, thus delta pressure is zero.

Regards,
highpressuretube

Also, in SimpleFoam pressure at the inlet is relative to the outlet. You can set the outlet pressure to fixedValue = 0 and see the DelP at the intlet. What are your inlet and outlet BCs?

[highpressuretube]

Well, certainly inlet pressure is relative, but the last picture with RE = 30.000 shows: dP = P_outlet - P_inlet = 0 - 0 = 0. So no pressure drop over throttle there.

Anyway my bcs:

p File:

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

internalField   uniform 0;

boundaryField
{
inlet
{
       type            zeroGradient;
}
outlet
{
        type            fixedValue; 
        value           uniform 0;
}
walls
{
        type            zeroGradient;
}
side1
{
        type            cyclic;
	neighbourPatch	side2;
}
side2
{
	type            cyclic;
	neighbourPatch	side1;
}
}

U File:

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

internalField   uniform (0 0 0);

boundaryField
{
inlet
{
        type            fixedValue; 
        value           uniform (9.4314 0.000000 0.000000);
}
outlet
{
	type            inletOutlet;
	inletValue	uniform (0. 0. 0.); 
        value           uniform (0. 0. 0.);
}
walls
{
        type            fixedValue; 
        value           uniform (0. 0. 0.);
}
side1
{
        type            cyclic;
	neighbourPatch	side2;
}
side2
{
	type            cyclic;
	neighbourPatch	side1;
}
}

Any rare assumptions there?

[arvindpj]

BCs look ok. You could try symmetry instead of cyclic, but this should not affect the results.

Final question, has the solution converged?

[highpressuretube]

I already tried symmetry instead of cyclic. As you suggest no difference there.

When it comes to convergence, looking at the residuals graph you might be right that the last jump to Re 30000 could run a little longer. I set it up with the default values from the pipe cyclic example.

What would I need to change to run longer? The final tol, the relative tol or both?

solvers
{
    p
    {
        solver          PCG;//GAMG;
        smoother        GaussSeidel;
	preconditioner	DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|k|epsilon)"
    {
        solver          PBiCG;//smoothSolver;
	preconditioner	DILU;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|epsilon)Final"
    {
        $U;
        tolerance       1e-05;
        relTol          0;
    }
    Phi
    {
        solver          GAMG;
        smoother        DIC;

        tolerance       1e-06;
        relTol          0.01;
    }
}

I am not sure that I understood correctly how tolerance and reltol are connected to the residuals graph because the latter never drops below 0.0001.

Thanks for all your input,
highpressuretube

[highpressuretube]

I looked closer into convergence and ran the first Re=30 for significantly longer.

Can anyone give me his opinion to the residuals graph. It looks locally diverging but it not getting any better by letting it run longer. So i would stop at 1250 iterations.

Besides that the calculation already run overnight so the calculation cost are already way to high to get there.

[gkarlsen]

Let me see if I understood your work correctly: You are running a simulation with constant velocity, and with lower and lower viscosity. You are observing lower and lower differential pressure. Qualitatively that would appear to make perfect sense, no?

The actual pressure gradient between inlet and outlet is much smaller than you expect/measure. My first guess would be that the dimensions of your geometry are messed up (e.g. the classic mm to m). Could that be it?

Kind regards,
Geir

[highpressuretube]

Qualitatively it makes perfect sense, apart from dP = 0 for RE = 30.000. It should become smaller but it looks like reversal in the picture and that together with the difference to the measurements got me confused.

I specified the blockmeshdict in mm with parameter convertToMeters 0.001 at the beginning and gave the flow in m/s. Should be ok there.

[arvindpj]

Could you share the blockMeshDict, let me try running this case and check.

[highpressuretube]

Of course, that would be too kind of you.

[arvindpj]

Originally Posted by highpressuretube

Of course, that would be too kind of you.

The blockMeshDict looks OK; except few cells exceeding the non-orthogonality criterion of 75 degs. You can look those up by checking the mesh quality by running:

checkMesh -constant -allTopology

And you can still take account of mesh non-orthogonality by increasing the nNonOrthogonalCorrectors to 2 in the fvSolution.

Also, instead of decreasing the kinematic viscosity to increase the Reynold number, for better stability, I used a ramped inlet b.c for velocity as below:

 
    inlet
    {
        type                uniformFixedValue;
        uniformValue        table
                            (
                                (0   (0 0 0))
                                (100 (1 0 0))
                            );
    }

As expected, I do see a pressure drop for Re=30000. Attached is the pressure variation along the axis.

If you are just studying the pressure drop, you could make use of coarser mesh with wall functions. This could save precious computational resources.
I am also attaching the pressure predictions for wall functions approach as well. With proper grid refinement study, you could get a better estimate of pressure drop with far less computational cost.

As a final thought, you could use cfMesh to generated good hex-dominant mesh in no time.

Cheers!

[philippose]

Hello Katrin,
A Good day to you.

I just looked through the issue you have had, and also ran a quick simulation setup using your blockmesh file.

Here are some comments and suggestions...:

1. For such wall bounded flows, I have had much better and consistent results (over the last 10 years) using the k-Omega SST turbulence model along with the nutUSpaldingWallFunction in the "nut" file for the walls.

I strongly suggest you move from k-Epsilon to k-Omega SST + nutUSpaldingWallFunction

2. Use and "inletOutlet" condition for the k and Omega boundary conditions for the outlet and set the "inletValue" to be equal to the internalField value (as initialisation)

3. In this particular case (non-orthogonal cells) use "2" non-orthogonal correctors

4. You could do with a much more coarser mesh in all the other areas of your mesh except for the region of the narrow pipe.


Using the settings given above along with "GAMG" for pressure and "smoothSolver" for all the other variables, I am getting very good results even by around 500 iterations (not converged yet though... may need to go up to around 800 - 1000 iterations)

delta-Pressure is around 3100 (at around 500 iterations... not yet converged), and the velocity and pressure fields look good.


Wish you a good day ahead.

Regards,
Philippose

[highpressuretube]

Hallo Phillipose,

thanks for the hints. 3100 delta pressure is exactly what I would expect from the measurements.

I will thus try to implement your input over xmas break and see where it takes me.

The dissappeared delta pressure seems indeed to be lack of convergence.

Why do I say that? Because I implemented Arvindpj's suggestions of using a table instead of ramping the viscosity and set the nonOrthogonalCorrectors to 2.
The result: delta pressure came back, though not the measured 3100 but I got rather stuck with more than 10.000.

I wonder how you reached 34Pa arvindpj? Could you share your solver choice and how you calculatet the 34Pa, too? Or did you set a totalPressure boundary?

Well lots of thoughts to try, thanks sooooo much for the input and I will get back when I have new results.

Have a very nice xmas and a happy new year. (I will be with my computer but offline for the break.)

highpressuretube

[arvindpj]

Attached are the set-ups.
Cheers: