This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Difference between grad(u) and epsilon = sym(grad(u)) in Stokes demos

+1 vote

Looking through the FEniCS 1.6.0 demos I've noticed that the Stokes and Navier-Stokes examples use the velocity gradient grad(u) rather than the strain rate tensor epsilon(u) = sym(grad(u)) = 1/2*(grad(u) + transpose(grad(u))). At first I thought this was because for these problems the two operators are equivalent. But, I recently compared the results of demo_stokes-taylor-hood.py which uses grad(u) and those using epsilon(u) as defined above, i.e.

a = (inner(sym(grad(u)), grad(v)) - div(v)*p + q*div(u))*dx

instead of

a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx

Here is a summary of the results (min and max magnitudes):

            |      u      |      p
            |-------------|-----------
            |  min   max  |  min   max
------------|-------------|-----------
grad(u)     |    0   2.1  |  -66   160
epsilon(u)  |    0   2.2  |  -33    78

Could someone please explain why the FEniCS demos use grad(u) instead of epsilon(u), why the results are different, and which results are correct?

asked Mar 8, 2016 by benzwick FEniCS Novice (350 points)

1 Answer

+1 vote

Both formulations are consistent in case of Dirchlet conditions and constant viscosity. In case of Neumann conditions the outflow behavior differs, but you cannot say which one is "correct", because it depends on what you (physically) assume to happen past the boundary (traction is zero vs. gradient is zero, in normal direction). In the discrete setting however, even for pure Dirchlet conditions and constant viscosity, the operators slightly differ, which can lead to different results at finite mesh sizes. In your case I assume you are dealing with nonempty outflow boundaries, which should explain the rather large deviations. Try to visualize the velocity vectors at the outflow boundary to see the difference.

Update: I just saw that your scaling is not correct. It should be 2*sym(grad(u)), otherwise this scales down your pressure (if the rhs is zero), and this is exactly what you observe. What I have written above still applies.

answered Mar 8, 2016 by Christian Waluga FEniCS Expert (12,310 points)
edited Mar 9, 2016 by Christian Waluga
...