|February 5, 2010, 08:34||
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?
|February 6, 2010, 04:23||
Join Date: Mar 2009
Location: London, England
Posts: 1,715Rep Power: 19
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.
|February 6, 2010, 06:09||
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...
|Thread||Thread Starter||Forum||Replies||Last Post|
|Strange behaviour when using LienCubicKE and NonlinearKEShih||hani||OpenFOAM Running, Solving & CFD||20||March 6, 2013 10:06|
|private boundary condition in controlDict strange behaviour||markc||OpenFOAM Programming & Development||12||November 23, 2011 10:18|
|Strange behaviour on outlet boundary using LES||segersson||OpenFOAM Running, Solving & CFD||0||December 9, 2009 03:57|
|strange behaviour of the one-equation SGS model||cfdmarkus||OpenFOAM Running, Solving & CFD||18||August 14, 2009 05:33|
|Strange multicomponent source behaviour||Zitron||CFX||4||July 12, 2007 15:32|