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

non zero divergence for incompressible flow!

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

Like Tree5Likes
  • 5 Post By alberto

Reply
 
LinkBack Thread Tools Display Modes
Old   June 17, 2010, 16:16
Default non zero divergence for incompressible flow!
  #1
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi all,

I'm simulating wake vortices in ground effect. I use icoFoam (incompressible) solver. But when I compute the divergence of U, I find non zero divergence (~-0.5 to 0.5).

I compute the divergence with : foamCalc div U and my own utility (using Gauss linear scheme). Both way I got the same results (with a difference less then 1%).

Can somebody explain me why I have such high divergence (I'm not expecting to have div(U) = 0 but more like div(U) ~= 0.01)

Thank you,

Pascal
Pascal_doran is offline   Reply With Quote

Old   June 18, 2010, 08:06
Default
  #2
Senior Member
 
Pavan
Join Date: May 2009
Location: Melbourne
Posts: 101
Rep Power: 8
rieuk is on a distinguished road
Hey Pascal,

Have you checked for convergence at each time-step (like looking at the initial residuals)? Maybe you need to increase your nCorr (number of PISO corrector steps)? Just guessing...
rieuk is offline   Reply With Quote

Old   June 18, 2010, 10:36
Default
  #3
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
You must check div(phi), using the phi coming out from the pressure equation, not div(U).

Actually the code does that for you when it prints the local and global continuity error. Search for the file continuityErrs.H to check.

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   June 18, 2010, 18:46
Default
  #4
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Thank you for your help

But can you explain me why it is better to check div(phi) than to check div(U) both should be near zero? And by the way what is a good value for div(U) and div(phi) for an incompressible flow (I know that in theory it should be div(U) = 0 ...)

Pascal
Pascal_doran is offline   Reply With Quote

Old   June 18, 2010, 19:16
Default
  #5
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
Hi Pascal,

Quote:
Originally Posted by Pascal_doran View Post
Thank you for your help
But can you explain me why it is better to check div(phi) than to check div(U) both should be near zero?
It is a question of consistency with how you discretise and solve the equations.

Short answer: you impose div(phi) = 0, meaning div(U_f) = 0, being U_f the velocity on the cell faces, not div(U) = 0, being U the velocity in cell centers.

Long answer: If you go through the derivation, you can think to semi-discretize the momentum equation as

A*U = H - grad(p)/rho

so that

U = (H/A) - grad(p)/(A*rho)

At this point you use the continuity equation

div(phi) = div(Uf) = 0

The continuity equation tells you the divergence of the flux is zero, and the flux is computed at cell faces, not at cell centres.

By interpolating the predicted U, obtained from the semi-discrete momentum equation, you get (S is the surface area vector, and |S| its norm)

phi = Uf = (H/A)_f - snGrad(p)/(rho*A)_f |S|

and, replacing this in the continuity equation

div((H/A)_f) - laplacian(1/(rho*A)_f |S|, p) = 0

Note that in OpenFOAM phi is computed as

phi = fvc::interpolate(U) + ddtPhiCorr(rUA, U, phi)

where U = H/A, and rUA = 1/A. The term ddtPhiCorr originates from the collocated grid arrangement.

Quote:
And by the way what is a good value for div(U) and div(phi) for an incompressible flow (I know that in theory it should be div(U) = 0 ...)
Very small. Take a look at icoFoam tutorials for the cavity case, which return something like:

time step continuity errors : sum local = 8.06059e-09, global = 9.28097e-19, cumulative = 9.00754e-18

Best,
gregor, sharonyue, kosan and 2 others like this.
__________________
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   June 18, 2010, 19:52
Default
  #6
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Thank you Alberto!

I really appreciate the short and long answer! It really helps me.

Regards,

Pascal
Pascal_doran is offline   Reply With Quote

Old   June 18, 2010, 23:00
Default
  #7
Senior Member
 
Pavan
Join Date: May 2009
Location: Melbourne
Posts: 101
Rep Power: 8
rieuk is on a distinguished road
But U_f and U should the same velocity distributions (or almost since U_f is just the interpolated U at cell surfaces) so why should the divergence be different for each of them?
rieuk is offline   Reply With Quote

Old   June 19, 2010, 11:00
Default
  #8
Member
 
Patricio Bohorquez
Join Date: Mar 2009
Location: Jaén, Spain
Posts: 94
Rep Power: 8
pbohorquez is on a distinguished road
Hi Pavan:

icoFoam solves for the integral form of Navier-Stokes equations. Therefore, you should compute the mass imbalance (in the computational cell) as the integral of div(U) in the cell volume V_P. Subsequently, you should apply Gauss theorem to transform the volumetric integral of div(U) as the sum of phi on the face. Do you know what does div(phi) mean?
pbohorquez is offline   Reply With Quote

Old   June 19, 2010, 21:42
Default
  #9
Senior Member
 
Pavan
Join Date: May 2009
Location: Melbourne
Posts: 101
Rep Power: 8
rieuk is on a distinguished road
Thanks Patricio, I forgot about the way it's calculating div in the discrete space.
rieuk is offline   Reply With Quote

Old   June 20, 2010, 15:31
Default
  #10
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
Maybe we should put this in the wiki FAQ
__________________
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   September 20, 2010, 16:17
Default
  #11
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi Alberto,

Do you have an idea how to compute the vorticity based on the flux phi like this :
curl(phi)

When I tried to compiled I got this error message :
Code:
zVorticityPhi.C:82: error: no matching function for call to ‘curl(Foam::surfaceScalarField&)’
make: *** [Make/linux64GccDPOpt/zVorticityPhi.o] Error 1
Here the main part of my code:

Code:
        IOobject phiheader
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );

    if (phiheader.headerOk())
    {
        Info<< "    Lecture de phi" << endl;
        surfaceScalarField phi(phiheader, mesh);
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
        Info<< "    Calcul de vorticity" << endl;
        volVectorField vorticity
        (
            IOobject
            (
                "vorticity",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ
            ),
            fvc::curl(phi)
        );
    }
Thank you,
Pascal
Pascal_doran is offline   Reply With Quote

Old   September 20, 2010, 16:38
Default
  #12
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
Hello, phi is a scalarField, so curl(phi) is not defined. You have to use a vectorField to compute the curl.

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   September 20, 2010, 17:18
Default
  #13
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Thanks for your reply,
So are you saying that I can't compute vorticity directly from the flux? Because phi must be a surfaceScalar and the vorticity must a volVector. What are you suggesting me to do since div(phi) is more accurate than div(U) I guess that curl based on phi would be more accurate than the curl based on U? What do you think?

Pascal
Pascal_doran is offline   Reply With Quote

Old   September 20, 2010, 17:23
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
What do you use the vorticity for?
__________________
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   September 20, 2010, 17:27
Default
  #15
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
I use the vorticity for tracking the position of wake vortices and for stability analysis.
Pascal_doran is offline   Reply With Quote

Old   September 20, 2010, 17:36
Default
  #16
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
What I said is phi is a scalar quantity (it is the U_f \cdot surface), while the curl operation is only defined for vectors.

The U in cell centres, which is what you visualize in paraview is not "inaccurate". It does not satisfy the continuity equation strictly, since the continuity constraint is applied to the flux.

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   September 21, 2010, 05:09
Default
  #17
Member
 
Patricio Bohorquez
Join Date: Mar 2009
Location: Jaén, Spain
Posts: 94
Rep Power: 8
pbohorquez is on a distinguished road
Quote:
Originally Posted by Pascal_doran View Post
Thanks for your reply,
So are you saying that I can't compute vorticity directly from the flux? Because phi must be a surfaceScalar and the vorticity must a volVector. What are you suggesting me to do since div(phi) is more accurate than div(U) I guess that curl based on phi would be more accurate than the curl based on U? What do you think?

Pascal
Hi Pascal,

that is a good question. But how can you relate phi with curl(U)? May be I am wrong, but I think that Gauss theorem cannot be applied to curl operator in order to express a volume integral as surface integrals. So I have difficulties to figure out a way of implementing a "conservative" discretization of curl. I think that it is usually treated as a source term.

To visualize vortex shedding you can just employ the "vorticity" command implemented by default in OF, and select the component to visualize in paraView. In the presence of non-orthogonal cells you may find jumps at element boundaries, as discussed in Tomboulides and Orszag (JFM, 2000, 416:45-73), so take care of them. I will post a nice picture showing vortex sheding in the near future.

Best wishes,
Patricio
pbohorquez is offline   Reply With Quote

Old   September 21, 2010, 10:22
Default
  #18
Member
 
Pascal
Join Date: Jun 2009
Location: Montreal
Posts: 65
Rep Power: 8
Pascal_doran is on a distinguished road
Hi Patricio,

I think you're totally right. I was just wondering what was the most efficient way to compute the vorticity. In my case the mesh is orthogonal so I will keep using the vorticity utility.

Thanks

Pascal
Pascal_doran is offline   Reply With Quote

Reply

Tags
incompressible divergence

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
Flow meter Design CD adapco Group Marketing CD-adapco 3 June 21, 2011 08:33
multiphase flow, quick divergence of contuinity eq violet FLUENT 7 February 3, 2009 18:54
flow past a missle: how to solve divergence? xiaofish FLUENT 0 September 9, 2007 22:53
reversed flow at velocity inlet / mass flow inlet ib FLUENT 1 March 26, 2007 13:11
transform navier-stokes eq. to euler-eq. pxyz Main CFD Forum 37 July 7, 2006 08:42


All times are GMT -4. The time now is 14:23.