CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   fvc::div() strange behaviour (http://www.cfd-online.com/Forums/openfoam-solving/72435-fvc-div-strange-behaviour.html)

 ivan_cozza February 5, 2010 09:34

fvc::div() strange behaviour

Dear Foamers,
I noticed a strange behaviour in fvc::div() operator:

I have a fluctuating velocity field ut, evaluated as follows:

ut = fvc::curl(psi), where psi is a vector-stream function. Obviously, ut is divergence-free, as div(ut) = div(curl(psi)) = 0 analitically.

In order to verify this identity I calculated fvc::div(ut), but it's nothing but zero! I found high values of divergence.

But if I do this test as follows:

Now, divUt is zero, consistently with the analytical expression.

Someone can explain me this discrepancy?

Thanks, Ivan

 hjasak February 6, 2010 05:23

Nice, but you are forgetting the discretisation error. Remember, field approximation of any FVM field is piecewise linear (with discontinuities). Thus, even if you have analyical values of stream function in cell centres, the field will not be perfect (eg who says it's linear?).

First, you use curl which will interpolate from cell centres to faces for the Gauss Theorem = discretisation error. This one you can actually avoid by analytically picking up the stream function values at Face (rather than cell) centres.

Second, the velocity vector is numerically calculated in the cell centres again, and for div() you again need it on faces: same story.

In conclusion: you are completely right in the analytical world, but in a discrete representation you will pick up discretisation error along the way: this is what you are seeing.

Enjoy,

Hrv

 ivan_cozza February 6, 2010 07:09

Hrvoje,
thank you for your reply, very good idea to calculate the stream-function on cell faces, I will implement it soon.
I have taken into account the discretization error in my test, so I didn't expect to see div(ut) = 0, but something like div(ut) = "something small related to numerics" (and this "something small" is another time questionable, how much small this div must be?).
The very strange feature that I noticed is that if I calculate the same thing with two equivalent methods I find two very different results, like:

div(ut) = 1000