[fsalvucci]

Well, finally, i tried with the setup fivos recomended, with a really dense mesh, even with the bounday layer more refined. Umax of my analytical solution is 0.1199 m/s, and i´m obtaining after 4000 iterations 0.1089 m/s.
I dont know, i dont want to bother anymore, but still i would prefer a more exact solution and dont know what else to do.

So, thanks to you all for your help!!!!

Sincerely,

Fernando

[alberto]

Why do you keep fixing the pressure at both your boundary conditions?

Best,

[Locky3827]

Hi,
I have run laminar pipe flows with both Newtonian and several varieties of non-Newtonian fluids (with a corrected shear rate) and always got the correct velocities and pressure gradients when compared with the analytical solutions. I used a fixed pressure at one end of the domain only and fixed the inlet velocity.
Alberto has given you good advice.
You should not need an extra fine mesh or a large number of iterations unless something is wrong.
Regards,
Lachlan

[fivos]

@Fernando : I will take a look later for your case, I suppose you are still trying with the same conditions (L = 0.04 m , R = 0.0015, (p2' - p1') = 0.026, nu = 3.047e-6). I will try it once more with exactly the values you are using and see what happens.

Just a note : a (fine) mesh is not good if it is skewed. Have you checked your skewness? What elements are you using? It is harder to get convergence on a Pyramids/Prisms/Wedges mesh than a Hex mesh. I would suggest taking an other look in your mesh. Try a much coarser mesh (which would take a few time to run), but with good quality.

@ alberto : sometimes you want/have to use pressure inlet/outlet boundary conditions, because you don't have any idea of the developing flow. Of course the specific case is rather simple, since there are analytic relations. But I guess that Fernando is just benchmarking and experimenting on something bigger he might try later on.

By the way, why setting pressure at inlet / outlet is not correct? It is true that fixing static pressure at inlet and outlet is not the best thing to do (it is better to specify total pressure at inlet, static pressure at outlet, since you explicitly define the total potential of the flow), but it should work in that simple case.

Anyway I can't believe that this setup cannot work. I have already tried it with ~1% deviation.

[Locky3827]

From a pipeline point of view normally you have a flowrate in mind and you want to find the pressure gradient as this determines the load on pumps etc for a complete pipe system. Alternatively you can specify the pressure at the inlet and outlet and then calculate the flowrate and hence velocity.
Regards,
Lachlan

[fivos]

Well I have rerun the pipe simulation for L = 0.04 m , R = 0.0015, (p2' - p1') = 0.026, nu = 3.047e-6, for two different meshes, one with refined near wall region and another without refinement (mesh sizes ~200000 and ~100000 cells respectively).

Analytical values Umean = 0.06m/sec, Umax = 0.12m/sec

Results are OK after 1000 iterations with low residuals (below 10^-6) :

Without near wall refinement (pipe1.jpg):
Umax = ~0.11898 m/sec (deviation 0.85%)
From vol. flow Umean = 0.060053....m/sec (deviation near 0)
Velocity profile is not very smooth at edges of the pipe (see graph)

With near wall refinement (pipe2.jpg)
Umax = ~0.11918 m/sec (deviation 0.68%)
From vol. flow Umean = 0.0595757....m/sec (deviation near 0.7%)
Near wall velocities are now OK!!

This proves that the problem is not due to OpenFoam, or the boundary conditions but probably due to your mesh, Fernando (since you used my fvSolution,fvSchemes, p, U files).

Use a less refined mesh (to save time), but make sure it has good quality and avoid prisms/wedges/pyramids.

No need for excessive iterations too.

[fsalvucci]

Well, i made the simulation with the classical boundary conditions, v at inlet v=0.25 m/s and pressure at outlet 0. Flow should develope in a parabolic manner with Umax = 2*Umean, being U mean the flat velocity imposed at the inlet, so Umax should be 0.5 m/s. The result is a parabolic profile with Umax = 0.46... m/s, that is 8% of error.

Then, the problem should be the mesh. I´ve been using two meshes:
* a non-structured mesh composed of tetras for the uniform refined mesh
* a non-structured mesh composed of tetras and prisms (wedges) for the refinig of the boundary layer.

This is because as you suposed, i´m doing this as a benchmark, what i have to do is simultaions like these but in complex geometries, like curved pipes with branches. So, this geometries y have to reconstruct them from images, and then generate a surface in a CAD program and them mesh it. And im using Netgen, which meshes with unestructured elements. For the refinement of the boundary zone, i use engrid.

Is there any meshing tool you can suggest me, that, form a NURBS surface form the CAD program can mesh with hexas???

As usual, thanks!!! The help of all of you is being very valuable.

[fivos]

Yes the problem is the mesh. Keep in mind that using pyramids/wedges/prisms increases numerical viscosity, since the flow is not (generally speaking) aligned with the flow. That' s why such meshes are harder to converge, comparing to hex meshes, and require high order schemes.

Regarding the mesh generator, I am using Gambit 2.3.16. It has the ability to create structured and unstructured meshes on 2d and 3d, can create boundary layers, it is able to create geometry from primitives (box, sphere, cylinder, etc) and you can import geometry from CAD programs (e.g. solidworks). Unfortunately this program is not free, it is (or was, since it has been replaced by ANSYS design modeler) the mesh generator for Fluent. Generally I can say that it is a very powerful geometry and mesh generator (even comparing to its succesor, ANSYS design modeler).

A free mesh generator that can create structured meshes is Discretiser (here : http://www.discretizer.org/ ). I have not used it, but It appears to be an impressive project. I might give it a try sometime.

Take also a look here, for mesh generators generally : http://www-users.informatik.rwth-aachen.de/~roberts/software.html

Good luck with your work.

[alberto]

Originally Posted by fivos

@ alberto : sometimes you want/have to use pressure inlet/outlet boundary conditions, because you don't have any idea of the developing flow. Of course the specific case is rather simple, since there are analytic relations. But I guess that Fernando is just benchmarking and experimenting on something bigger he might try later on.

That was clear. My doubt was on the setup of the problem from a mathematical point of view, since he is using an incompressible solver.

[alberto]

Originally Posted by fsalvucci

Is there any meshing tool you can suggest me, that, form a NURBS surface form the CAD program can mesh with hexas???

If the geometry is very complicated, it is often not worth trying to mesh it with hexahedra. What you are doing on the mesh generation is correct: use hex cells on boundaries and tetrahedral cells in the core of the flow.

What you can do to improve the stability and the convergence is to properly set up the numerics.

A few suggestions:

- Gradient scheme: cellLimited leastSquares 1;
- Convection scheme:
div(phi,U) linearUpwindV cellLimited leastSquares 1;
For turbulence quantities: linearUpwind cellLimited leastSquares 1;

You might need to limit also the laplacian, replacing "corrected" with "limited 0.33" in case your grid is not very good.

Finally, if you have non-orthogonal cells (run checkMesh to see), perform more corrector steps on the pressure equation. As a rule of thumb the last corrector should have a small (nearly zero) number of iterations.

Short question: is your problem to determine the flow in a pipeline, given a pressure loss/gradient? If so, the less troublesome way to do this, even if it appears to be more complicated, is to map inlet and outlet (or impose periodicity if possible), and add a forcing term to the momentum equation.

Best,

[fsalvucci]

Ok, thanks a lot again with all the sugests, with the mesh and with de schemes. I´ll take it into account and perform a simulation with those schemes.

When i check the mesh, about non-orgthogonality and skewness:

Mesh non-orthogonality Max: 65.206 average: 11.4348
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.717552 OK.

So, i should set up some nonorthogonal correctors, right? like, how many?

My final purpose will be to run the simulations in the complex geometries specifying as inlet condition a velocity profile, actually, i have already performed them. but, i needed to show first that the simulation procedure and technique are correct and that for a case with analytical solution, everything is good.

So,for this last point, i will generate a mesh of a pipe with hexas, and for the simulations with complex geometries, i will compare results with an hexa mesh (if i can really get them with for example discretizer) and with the tetras and prisms mesh with these schemes alberto is suggesting.

Thanks a lot, really, you helped me a lot. Ill write in a few days telling you what happended.

Fernando

[alberto]

A useful script is attached to generate a mesh in a pipe (I just adapted them from an old post, the file is NOT my work).

The pipe section is split in a central square and four curvilinear quadrilaterals.

Simply edit the file setting:

L - Length
D - Diameter

NPS - Nodes along the square edge
NPD - Nodes from square to perimeter
NPY - Nodes along pipe length

Generate the blockMeshDict with:

m4 cylinderMesh.m4 > blockMeshDict

Save the blockMeshDict in your polyMesh directory and run blockMesh on your case.

The BC names are "inlet", "outlet", "walls".

P.S. I was insisting on not fixing the pressure at both sides because p is used to impose the continuity in incompressible solvers

Best,
Alberto

[alberto]

Here's a version with mesh grading, prepared by a student I am helping with OF (Antonello Baccoli).

P.S. You need to install the m4 interpreter on your workstation to use these script. It is usually already available or packaged for all the major Linux distributions.

Best,

[NorbertB]

Hi

@fivos : You should never put fixedValue on both inlet and outlet, this could lead to divergency issue , especially for an incompressible case.
Usually inlet is zeroGradient (or similar) for pressure and you specify your velocity.
Then, outlet is fixedValue for pressure and zeroGradient for velocity. Note that this last condition could "kill swirls" if your duct/pipe/geometry is not long enough. And this leads to errors.

[fivos]

Setting the pressure at the inlet and outlet represents a pressure inlet and outlet respectively, isn't it? I understand that pressure in an incompressible solver is used to enforce incompressibility (div(u)=0), but how are you going to model some cases where you have a known pressure difference and you want the velocity field? Actually I have came across a case where I should calculate the velocity field at a nozzle for a given total pressure. In that case I fixed the total pressure at the inlet and static pressure at the outlet, obtaining the velocity field (velocity field was zeroGradient at inlets and pressureInletOutletVelocity at outlets) and mass flow through the nozzle. How such a case should be modeled? Should I try different mass flows through the nozzle until I find the one giving the desired pressure gradient?

I don't have much experience in OpenFoam, I more familiar with Fluent/CFX. In these programs you are allowed to set pressure at a pressure inlet. It is true that setting pressure is not a very robust way to solve your problem, actually setting static pressures at inlet and outlet is indeed considered unreliable in the CFX manual, setting a total pressure at the inlet and static pressure at the outlet is advised to do instead. Since the way of solving the NS equations is more or less the same for all of these programs ( at least I guess so), how it is possible for them to model the pressure inlet?

I apologise for my ignorance, hope you can answer some of my questions.

[l_r_mcglashan]

I've attached a solver for using periodic boundaries where you have a fixed mass flow rate. Use cyclic boundary conditions on variables in the 0/ folder.

Have a look anyway.

[maddalena]

Hi everybody,
I am studying pipe flow as well. At the moment I want to investigate the effect of using different meshes on the same pipe, in order to assess what the overall numerical error will be when I will switch to a more complex geometry.
I followed this thread on the schemes setup; actually, these are the standard schemes I am used to apply. As for the BC, I also applied this:

Originally Posted by alberto

A proper setup would be to specify the velocity at the inlet and the pressure at the outlet, as in the case I run, which is equivalent to yours, and gives the exact result for the radial velocity profile (the mesh is much less refined than yours).

This is 0/U:

inlet           
    {
        type            fixedValue;
        value         uniform (6 0 0);  // Obtained as known massFlow/Area
    }

On a second step, I would like to impose the mass flow in place of the velocity at the inlet:

    inlet           
    {
        type            flowRateInletVelocity;
        flowRate        0.3205;                    // Volumetric/mass flow rate [m3/s or kg/s]
        value           $internalField;         // placeholder
    }

However, this leads to an unstable solution, which oscillates and gives unrealistic velocity distribution.
My questions are:

• what is the reason why the fixedValue U gives a stable solution and the massFlow does not, even if in the two cases the velocity is fixed?
• Is it ok to use this setup if I want to fix the mass flow?

Thank you for any help.

maddalena

[stainboy]

Dear colleagues,

could you please help me with validation of 3D Poiseuille flow. I'm quite depressed, I've spend hundred of hours to validate my simulation against analytical solution but I can't find where I'm making an error. So let's start from the beggining.

I would like to perform simulation of fluid which have kinematic viscosity nu=3.33-e6 (blood). Pipe have radius R = 0.1m and length L = 2m.

Under this link you can find stl mesh: Mesh

To generate FVM I'm using netgen. Mesh is generated without any problems and check mesh does not comply on the parameters. Check mesh is here

My approach is to use the simpleFoam algorithm and given pressure drop between the inlet and outlet. According to Poiseuille equation
V_max = (DP_kin * R^2)/4*nu*L. I would like to have V_max = 1m/s therefore DP = V_max *4 * nu * L / R^2 = 1 * 4 * 3.33e-6 * 2 / 10e-2 = 53.28e-2. I'm putting this number as value of the fixedValue bc type at the inlet patch in 0/p file, rest of the patches in this file are defined as zeroGradient and fixedValue = 0 respectively for wall and outlet. When it comes to 0/U definition I'm using as follows:

• wall - fixedValue uniform (0 0 0)
• inlet - pressureDirectedInletVelocity; inletDirection uniform ( 0 0 1 ); value uniform ( 0 0 0 );
• outlet - zeroGradient

In transportProperties I'm defining nu:

transportModel  Newtonian;
nu              nu [0 2 -1 0 0 0 0] 3.33e-06;

When solving case there is no problem with convergence at all but I was very suprised when I saw results in paraview.
Chart pressure is fine. The difference between the inlet and the outlet is the same as I determined in the boundary conditions file, but the speed is far from the expected results.

Here is the package with plots: Plots

Could you please help me with validation of my case? Here is the fully configured case, just run simpleFoam: Case

I'm using wrong equation or making some stupid mistake which I cannot find.

In addition I was performing many simulations with different meshes from 30K till 2M elements with and without boundary layer and did not obtaind correct results.

Thank you in advance for your help.

Here are some additional info related to solver configuration
fvSchemes

ddtSchemes
{
	default steadyState;
}

gradSchemes
{
	default Gauss linear;
}

divSchemes
{
	default Gauss linear;
	div(phi,U) Gauss upwind;
}

laplacianSchemes
{
	default Gauss linear corrected;
}

interpolationSchemes
{
	default linear;
}

snGradSchemes
{
	default corrected;
}

fluxRequired
{
	default no;
	p;
}

fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.7.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
	p
	{
		solver GAMG;
		agglomerator faceAreaPair;
		smoother GaussSeidel;
		nCellsInCoarsestLevel 10;
		mergeLevels 1;
		nPreSweeps 0;
		nPostSweeps 1;
		cacheAgglomeration on;
		tolerance 1e-9;
		relTol 0.001;
		
	}
	
	U
	{
		solver PBiCG;
		preconditioner
		{
			preconditioner DILU;
		}
		tolerance 1e-08;
		relTol 0.0001;
	}
}

SIMPLE
{
	nNonOrthogonalCorrectors 2;
	pRefCell 0;
	pRefValue 0;
	residualControl
	{
		p 1e-05;
		U 1e-05;
	}
}

relaxationFactors
{
	p 0.2;
	U 0.5;
}

potentialFlow
{
	nNonOrthogonalCorrectors 10;
}

[Locky3827]

Surely the flow would be turbulent in a 0.1 m pipe at this speed?

[stainboy]

Ah yes, in this configuration Re would be around 30K. I will check the case for Re = 1.

Thank you for suggestions.