fvc::div() strange behaviour
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:
volTensorField gradUt = fvc::grad(ut);
volScalarField divUt = gradUt.component(0) + gradUt.component(4) + gradUt.component(8); (Div is the trace of Grad, so...)
Now, divUt is zero, consistently with the analytical expression.
Someone can explain me this discrepancy?
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.
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
tr(grad(ut)) = 1x10^-9
but div(ut) = tr(grad(ut)).
I wonder if I'm missing something...
|All times are GMT -4. The time now is 23:23.|