CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Numerical error or case error? Flow in a 3D pipe

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 6, 2010, 11:50
Default Numerical error or case error? Flow in a 3D pipe
  #1
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Hi everyone! i´ve trying to find out what i am doing wrong, but still cant figure out. Hope someone can help, since im stucked and work is not advancing.
The thing is, i am runnig laminar, incompressible, steady flow through a 3D pipe. I´m setting as bounday conditions a pressure at the inlet and a pressure at the outlet, in order to compare with Poiseuille´s analytical solution.

Poiseuille predicts:

Umax = ((p2 - p1) * R^2)/(4 * mu* L)

Since my flow is incompressible (simpleFoam), i have to enter nu instead of mu, that in my case is nu = 3.047e-6. So, p2 and p1 i have to enter them as "kinematic pressures", that is p/rho. So, the equation could be also:

Umax = ((p2' - p1') * R^2)/(4* nu * L),

being p1' and p2' the respective kinematic pressures.

For may case (L = 0.04 m , R = 0.0015, (p2' - p1') = 0.026), the equation predicts Umax = 0.119 m/s.

The thing is that when i run the case, the flow developes but Umax of the parabola is 0.84 m/s.!!!!!

I have actually run the case for different values of L, R and nu, but the numerical resulta always underestimates the analytical solution.

I attach the images of my results, and the case, in case anyone wants to chek the configuration or even run the case and take a look. I could not include the mesh in "constant" since it was too heavy and couldnt upload it. It is a simple 3D pipe with L = 0.04 and R = 0.0015)

Thanks very much nad hope someone could help me!!

Regard,

Fernando
Attached Images
File Type: png p.png (14.0 KB, 98 views)
File Type: png U.png (10.4 KB, 102 views)
File Type: png U_out.png (2.9 KB, 99 views)
Attached Files
File Type: gz pipe-3D.tar.gz (2.3 KB, 41 views)
fsalvucci is offline   Reply With Quote

Old   August 6, 2010, 14:47
Default
  #2
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
Have you tried refining your mesh?

You're using upwind on velocity, try using linear or GammaV, or anything higher order.

Not too sure about your boundary conditions either.
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 6, 2010, 19:05
Default
  #3
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 699
Rep Power: 12
msrinath80 is on a distinguished road
Quote:
Originally Posted by fsalvucci View Post
Hi everyone! i´ve trying to find out what i am doing wrong, but still cant figure out. Hope someone can help, since im stucked and work is not advancing.
The thing is, i am runnig laminar, incompressible, steady flow through a 3D pipe. I´m setting as bounday conditions a pressure at the inlet and a pressure at the outlet, in order to compare with Poiseuille´s analytical solution.

Poiseuille predicts:

Umax = ((p2 - p1) * R^2)/(4 * mu* L)

Since my flow is incompressible (simpleFoam), i have to enter nu instead of mu, that in my case is nu = 3.047e-6. So, p2 and p1 i have to enter them as "kinematic pressures", that is p/rho. So, the equation could be also:

Umax = ((p2' - p1') * R^2)/(4* nu * L),

being p1' and p2' the respective kinematic pressures.

For may case (L = 0.04 m , R = 0.0015, (p2' - p1') = 0.026), the equation predicts Umax = 0.119 m/s.

The thing is that when i run the case, the flow developes but Umax of the parabola is 0.84 m/s.!!!!!

I have actually run the case for different values of L, R and nu, but the numerical resulta always underestimates the analytical solution.

I attach the images of my results, and the case, in case anyone wants to chek the configuration or even run the case and take a look. I could not include the mesh in "constant" since it was too heavy and couldnt upload it. It is a simple 3D pipe with L = 0.04 and R = 0.0015)

Thanks very much nad hope someone could help me!!

Regard,

Fernando
Paraview is reporting 0.084 m/s NOT 0.84 m/s. Perform some additional refinement near the walls and tighten the pressure tolerances and all will be well
msrinath80 is offline   Reply With Quote

Old   August 8, 2010, 20:46
Default
  #4
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
In Poiseuille flow in a cylindrical pipe, the velocity is

U(r) = 2*Umean(1-(r/R)^2)

being Umax = 2*Umean, and being Umean the velocity you specify at a uniform inlet. No need to go through pressure drops :-)

Simply set the inlet at Umean, and you will get a perfectly parabolic profile if the pipe is long enough. Keep in mind that you must be at least at 10 diameters from the inlet and the outlet to obtain a fully developed flow.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   August 9, 2010, 10:53
Default
  #5
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Hi everyone for the replies. Let´s see:

Quote:
Originally Posted by l_r_mcglashan View Post
Have you tried refining your mesh?

You're using upwind on velocity, try using linear or GammaV, or anything higher order.

Not too sure about your boundary conditions either.

You mean in fvSchemes, in the divSchemes? You mean, instead of div(phi,U) Gauss upwind; i should try with linear?
fsalvucci is offline   Reply With Quote

Old   August 9, 2010, 10:55
Default
  #6
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
And,

Quote:
Originally Posted by msrinath80 View Post
Paraview is reporting 0.084 m/s NOT 0.84 m/s. Perform some additional refinement near the walls and tighten the pressure tolerances and all will be well
Yes, i typed it wrong. When you say "tighten the pressure tolerances", you mean, in fvSolution, in

solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.01;
}


i should change the tolerance for, let´s say, 1e-07 ?

Thanks a lot!!!
fsalvucci is offline   Reply With Quote

Old   August 18, 2010, 16:50
Default
  #7
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Hi everyone! Still with this problem.

I refined the mesh a lot. Still the Umax value is underestimated.

I changed in the fvSchemes, the scheme for divergence to:
div(phi,U) Gauss linear;
The simulation doesnt converge and U values go to 10e+8 m/s.

I tightened the p tolerance to 10e-7, still the Umax value is underestimated.

Anyone has a clue on what´s happening? If you want to check my case settings, i have attached it in a previous post.

Thanks!
fsalvucci is offline   Reply With Quote

Old   August 18, 2010, 18:51
Default
  #8
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Please make the case available, or it is impossible to see what is wrong.

P.S. I would use div(phi, U) Gauss limitedLinearV 1;
Your mesh is probably not fine enough to satisfy the stability criterion required by the linear scheme.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; August 18, 2010 at 18:52. Reason: Added P.S.
alberto is offline   Reply With Quote

Old   August 19, 2010, 04:54
Default
  #9
Senior Member
 
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 14
l_r_mcglashan will become famous soon enough
When I first started with CFD I did a laminar pipe flow (don't we all). Provided the case is set up properly, the two things to watch out for are mesh density and to make sure you are using a second order discretisation scheme on velocity. I've attached a graph showing the importance of mesh density.
Attached Files
File Type: gz convergence1.eps.gz (7.2 KB, 50 views)
__________________
Laurence R. McGlashan :: Website
l_r_mcglashan is offline   Reply With Quote

Old   August 19, 2010, 08:52
Default
  #10
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Quote:
Originally Posted by alberto View Post
Please make the case available, or it is impossible to see what is wrong.

P.S. I would use div(phi, U) Gauss limitedLinearV 1;
Your mesh is probably not fine enough to satisfy the stability criterion required by the linear scheme.

Best,
Look up, at the begginig of the post, i attach the pipe-3D file with the case. The mesh is too big to be attached.

I will try with that scheme, and see what happens. Thanks a lot!
fsalvucci is offline   Reply With Quote

Old   August 19, 2010, 10:28
Default
  #11
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by fsalvucci View Post
Look up, at the begginig of the post, i attach the pipe-3D file with the case. The mesh is too big to be attached.
- In the case there is no mesh => impossible to check it
- Pressure at outlet should be fixedValue.

Best,
Alberto
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   August 19, 2010, 10:55
Default
  #12
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Quote:
Originally Posted by alberto View Post
- - Pressure at outlet should be fixedValue.
It is actually fixedValue 0;

(outlet patch is salida, sorry for the spanish)
fsalvucci is offline   Reply With Quote

Old   August 19, 2010, 13:26
Default
  #13
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Quote:
Originally Posted by l_r_mcglashan View Post
When I first started with CFD I did a laminar pipe flow (don't we all). Provided the case is set up properly, the two things to watch out for are mesh density and to make sure you are using a second order discretisation scheme on velocity. I've attached a graph showing the importance of mesh density.
Ok, lets assume the case is correctly set up, and that mesh density is quite good (the mesh consists of 1,300,000 elements, the geometry is a pipe of 0.08 m lenght and 0.0015 m radius), the thing would be to choose a correct second order scheme for U, right?

The div(phi,U) scheme i was using was upwind. You mean that i should use another scheme? for example, the one that alberto suggests:

div(phi, U) Gauss limitedLinearV 1;

I´m running the simulation now with that scheme, ill se what happens.
What other suggestion canyou make me? Another scheme? Anything else?
Thanks a lot!
fsalvucci is offline   Reply With Quote

Old   August 19, 2010, 17:19
Default
  #14
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by fsalvucci View Post
It is actually fixedValue 0;

(outlet patch is salida, sorry for the spanish)
Sorry, I confused p and U files.

I do not think the setup is correct, since you solve for an incompressible flow and you force the pressure at both inlet and outlet.

Also, if you check the residuals your case does not converge.

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).

http://dl.dropbox.com/u/659842/pipe-3D_UInlet.tar.gz

If you do not want to use this setup, then your case is analogous to a periodic pipe with a fixed forcing term in the momentum equation, where inlet and outlet are replaced by cyclic conditions, and your UEqn is modified to include the forcing term.

Best,
Alberto
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   August 20, 2010, 08:03
Default
  #15
Senior Member
 
Phoevos
Join Date: Mar 2009
Posts: 103
Rep Power: 8
fivos is on a distinguished road
Dear Fernando,

I took a look into OpenFOAM simpleFoam solver for a similar case to the one you describe. I set up a 3d pipe with D = 0.001 m, L =0.01 m and set a pressure difference of 0.032 m2/s2 (kinematic pressure, real pressure 32 Pa assuming water with ρ= 1000kg/m3).

Expected Umean= (R^2 * Δp ) / (8 * μ * L ) = 0.1 m/sec (for water dynamic viscosity μ=1E-3 Pa.s) and expected Reynolds 100. Expected Umax = 2 * Umean = 0.2 m/sec

For the mesh : element size 0.00005 m (~20 elements on diameter), see also attached figure (mesh.png).

Turbulence is turned off in the RASproperties, but a laminar model is used as a RASModel.

Take a look at the calculated results below >

Calculated max velocity 0.198m/sec
Calculated vol. flow through inlet/outlet : 7.584e-8 m3/sec (see in the velocity.png at the table below the results window the integration of velocity U over the inlet surface, the third component). This means that there is an average velocity of 0.09656m/sec (close to the theoretical 0.1m/sec).

Those results were obtained after 1000 iterations.

A better mesh and more iterations might give better results.
Attached Images
File Type: jpg mesh.jpg (38.0 KB, 77 views)
File Type: jpg pressure.jpg (13.4 KB, 48 views)
File Type: jpg velocity.jpg (66.0 KB, 59 views)

Last edited by fivos; August 20, 2010 at 08:32.
fivos is offline   Reply With Quote

Old   August 20, 2010, 08:05
Default
  #16
Senior Member
 
Phoevos
Join Date: Mar 2009
Posts: 103
Rep Power: 8
fivos is on a distinguished road
Boundary conditions >

p file

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

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

internalField uniform 0;

boundaryField
{
outlet
{
type fixedValue;
value uniform 0 ;
}
Inlet
{
type fixedValue;
value uniform 0.032 ;
}
wall.2
{
type zeroGradient;
}
}

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

U file

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

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

internalField uniform (0 0 0);

boundaryField
{
outlet
{
type zeroGradient;
}
Inlet
{
type zeroGradient;
}
wall.2
{
type fixedValue;
value uniform (0 0 0);
}
}

// ************************************************** *********************** //
fivos is offline   Reply With Quote

Old   August 20, 2010, 08:08
Default
  #17
Senior Member
 
Phoevos
Join Date: Mar 2009
Posts: 103
Rep Power: 8
fivos is on a distinguished road
Finally
- controlDict >

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

application simpleFoam;

startFrom latestTime;

startTime 0;

stopAt endTime;

endTime 1000;

deltaT 1;

writeControl timeStep;

writeInterval 200;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression uncompressed;

timeFormat general;

timePrecision 6;

runTimeModifiable yes;


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

- fvSchemes >

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

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss linearUpwind Gauss linear;
div(phi,nuTilda) Gauss linearUpwind Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(U) linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}


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


- and fvSolution >

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

solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.1;
smoother GaussSeidel;
nPreSweeps 1;
nPostSweeps 2;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;
}

U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-08;
relTol 0.1;
}

nuTilda
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-08;
relTol 0.1;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}

relaxationFactors
{
default 0;
p 0.3;
U 0.7;
nuTilda 0.7;
}


// ************************************************** *********************** //
fivos is offline   Reply With Quote

Old   August 20, 2010, 08:15
Default
  #18
Senior Member
 
Phoevos
Join Date: Mar 2009
Posts: 103
Rep Power: 8
fivos is on a distinguished road
Also the residuals for the final steps >>

Time = 999

smoothSolver: Solving for Ux, Initial residual = 0.00059155, Final residual = 1.14422e-05, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.000577205, Final residual = 1.06019e-05, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 5.82043e-06, Final residual = 1.12108e-07, No Iterations 4
GAMG: Solving for p, Initial residual = 1.25485e-07, Final residual = 6.79484e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 9.51563e-09, Final residual = 9.51563e-09, No Iterations 0
time step continuity errors : sum local = 3.59221e-07, global = -2.41345e-09, cumulative = -1.71276e-05
ExecutionTime = 302.97 s ClockTime = 304 s

Time = 1000

smoothSolver: Solving for Ux, Initial residual = 0.000587961, Final residual = 1.13731e-05, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.000573677, Final residual = 1.05372e-05, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 5.78361e-06, Final residual = 1.11398e-07, No Iterations 4
GAMG: Solving for p, Initial residual = 1.24755e-07, Final residual = 6.75283e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 9.4657e-09, Final residual = 9.4657e-09, No Iterations 0
time step continuity errors : sum local = 3.57335e-07, global = -2.44649e-09, cumulative = -1.71301e-05
ExecutionTime = 303.7 s ClockTime = 304 s

End
fivos is offline   Reply With Quote

Old   August 20, 2010, 08:56
Default
  #19
Senior Member
 
Phoevos
Join Date: Mar 2009
Posts: 103
Rep Power: 8
fivos is on a distinguished road
A final post :

I scaled my geometry x 3 for all dimensions, so I would get R = 0.0015m and L = 0.03m.
I tried again with Δp'=0.026m2/s2 kinematic pressure, or 26Pa pressure. Expected mean velocity : 0.0015^2*26/(8 * 3.046e-3*0.03) = 0.08m/sec and max velocity 0.16.

Results calculated by OpenFoam are again OK (controlDict, fvShemes and fvSolution the same as before), see below.

Final residuals >

Time = 999

smoothSolver: Solving for Ux, Initial residual = 0.000158666, Final residual = 3.91358e-06, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.000153262, Final residual = 3.60992e-06, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 1.24232e-06, Final residual = 3.12179e-08, No Iterations 4
GAMG: Solving for p, Initial residual = 2.72346e-08, Final residual = 4.86485e-09, No Iterations 2
GAMG: Solving for p, Initial residual = 5.10714e-09, Final residual = 5.10714e-09, No Iterations 0
time step continuity errors : sum local = 5.6672e-08, global = -4.12545e-10, cumulative = -0.00149578
ExecutionTime = 354.72 s ClockTime = 355 s

Time = 1000

smoothSolver: Solving for Ux, Initial residual = 0.000157434, Final residual = 3.88319e-06, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.000152073, Final residual = 3.58192e-06, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 1.23268e-06, Final residual = 3.09756e-08, No Iterations 4
GAMG: Solving for p, Initial residual = 2.70237e-08, Final residual = 9.97178e-09, No Iterations 1
GAMG: Solving for p, Initial residual = 1.0082e-08, Final residual = 4.81098e-09, No Iterations 1
time step continuity errors : sum local = 5.33856e-08, global = -4.26169e-10, cumulative = -0.00149578
ExecutionTime = 355.75 s ClockTime = 356 s

End
Attached Images
File Type: jpg velocity2.jpg (64.9 KB, 46 views)
fivos is offline   Reply With Quote

Old   August 20, 2010, 09:04
Default
  #20
New Member
 
Fernando
Join Date: Feb 2010
Posts: 28
Rep Power: 7
fsalvucci is on a distinguished road
Well, i would like to thank you all very much for tour corncern and your replies.
I ´m running now my case with 1500 iterations instead of 800 as i did before, with the div scheme suggested (limitedLinearV 1), and i´ll see what happens. If it still doenst reach the steady state analytical solution, i will try with the setup fivos is suggesting.

So, ill write again in a couple of hours telling you what has happended.

Thank you all!!!!
fsalvucci is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Direct Numerical Simulation- for pipe flow RameshK Main CFD Forum 2 December 25, 2009 14:41
pipe in pipe heat exchanger JohannV FLUENT 3 December 3, 2009 03:53
My Revised "Time Vs Energy" Article For Review Abhi Main CFD Forum 2 July 9, 2002 09:08
Terrible Mistake In Fluid Dynamics History Abhi Main CFD Forum 12 July 8, 2002 09:11
fluid flow fundas ram Main CFD Forum 5 June 17, 2000 21:31


All times are GMT -4. The time now is 07:27.