I'm trying to find steady solutions of Navier-Stokes in an axisymmetric problem at low Reynolds numbers. We're using P2 elements for velocity and P1 for pressure
As the Navier-Stokes equations have $\frac{1}{r^2}$ terms in them we multiply the governing equations by $r$. Instead of using $dx$ at the end of the variational from we use $r dx$ as this gives us the axisymmetric finite elements.
This introduces a lot of $r^2$ and $r$ terms in the weak form. I was originally getting these by:
r = project(Expression('x[1] '),P1)
r2 = r*r
This worked well and we got what appeared to be reasonable solutions.
Whilst debugging a completely different problem I thought that actually seeing as $r^2$ is quadratic it should be approximated with P2 elements. I then changed it to:
r2 =project(Expression('x[1]*x[1] '),P2)
I thought this would make little difference but it causes very large oscillations along the line $r=0$. These oscillations are confined to a very small region near $r=0$ when using a very fine mesh.
My questions are therefore: Which implementation of $r^2$ is technically correct? How does FEniCS treat r=r*r (I would have though that if I bumped up the quadrature enough I would start to see the same oscillations I got when using P2 elements for r2 but I didn't observe this when I tried)? Is this a problem common with axisymmetric problems and are there any tricks round it?
Thanks for any help in advance!
