Cylindrical coordinates

[Åsmund Hjulstad]

I am trying to solve fluid flow between concentric cylinders, driven by the rotation of the outer cylinder. Flow is incompressible, stationary and axi-symmetric, so I wanted to solve it in 1D.

This might be me not understanding FEM, only, but perhaps someone has an easy pointer on how to do cylindrical coordinates.

The PDE  

   - mu * div (grad( u ) )  == 0

multiply with test function and integrate

   - mu * Integral [ div (grad(u)) * v r dr dtheta ]   == 0    (I multiply with r as it is polar coordinates)

It does not vary with theta, so I get a 2*pi factor that goes away (as rhs==0)

Integration by parts gives 

   - mu * Integral [ (u - r u') * v'  dr  ] == 0

I cannot get it to match an exact solution, though. 

Help most appreciated. 

[Anders Logg]

In cylindrical coordinates with d/dtheta = d/dz = 0, you have

-div(grad(u)) = -1/r d/dr (r du/dr)

Integrate against test function v and measure 2pi*r*dr (factor 2pi goes away):

-int 1/r*d/dr(r du/dr) v r dr = -int d/dr(r du/dr) v dr = int u' v' r dr where u' = du/dr and v' = v/dr.

You also get a boundary term at r = R which you need to handle differently depending on what boundary conditions you have. Note that the boundary term at r = 0 goes away.

Note that your physical domain (a disc with radius R) is different from the computational domain (the square [0, 2pi] x [0, R] or just [0, R]).

--

Anders

--
You received this message because you are subscribed to the Google Groups "fenics-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fenics-suppor...@googlegroups.com.
To post to this group, send email to fenics-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fenics-support/46f4499b-d047-48f5-8206-79a094686e0e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Åsmund Hjulstad]

Thanks. Still doesn't match exact solution, though, even when increasing to a very fine grid it gives me a relative error of 7e-4.

I wonder why it is correct to consider the flow velocity (u) as only a scalar, when it is in fact a vector {0, u[r], 0}   in (r,theta,z)-coordinates. 

My domain is between two concentric cylinders, so the computational domain is [r_inner, r_outer]. Boundary condition is u==0 on the inner surface and constant rotation on the outside. u[r_outer]==\omega r_outer, implemented in my fenics code using DirichletBC.

- Åsmund.

tirsdag 7. februar 2017 08.19.49 UTC+1 skrev Anders Logg følgende:

In cylindrical coordinates with d/dtheta = d/dz = 0, you have

-div(grad(u)) = -1/r d/dr (r du/dr)

Integrate against test function v and measure 2pi*r*dr (factor 2pi goes away):

-int 1/r*d/dr(r du/dr) v r dr = -int d/dr(r du/dr) v dr = int u' v' r dr where u' = du/dr and v' = v/dr.

You also get a boundary term at r = R which you need to handle differently depending on what boundary conditions you have. Note that the boundary term at r = 0 goes away.

Note that your physical domain (a disc with radius R) is different from the computational domain (the square [0, 2pi] x [0, R] or just [0, R]).

--

Anders

mån 6 feb. 2017 kl 22:13 skrev Åsmund Hjulstad <asmund....@gmail.com>:

I am trying to solve fluid flow between concentric cylinders, driven by the rotation of the outer cylinder. Flow is incompressible, stationary and axi-symmetric, so I wanted to solve it in 1D.

This might be me not understanding FEM, only, but perhaps someone has an easy pointer on how to do cylindrical coordinates.

The PDE  

   - mu * div (grad( u ) )  == 0

multiply with test function and integrate

   - mu * Integral [ div (grad(u)) * v r dr dtheta ]   == 0    (I multiply with r as it is polar coordinates)

It does not vary with theta, so I get a 2*pi factor that goes away (as rhs==0)

Integration by parts gives 

   - mu * Integral [ (u - r u') * v'  dr  ] == 0

I cannot get it to match an exact solution, though. 

Help most appreciated. 

--
You received this message because you are subscribed to the Google Groups "fenics-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fenics-suppor...@googlegroups.com.

To post to this group, send email to fenics...@googlegroups.com.

[Kristian]

If you use mshr to generate the mesh, you should know that the mesh is a polygon with 16 (I think) sides as default, but this can be increased.

To post to this group, send email to fenics-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fenics-support/f71f6934-d8de-43dc-8ddd-b89ffc72975c%40googlegroups.com.

[Jack Hale]

Probably time to post some code, maybe FEniCS Q&A might be better forum?

Although it's still a work in progress, we've been trying to solve
some curved shell problems in FEniCS using concepts from differential
geometry to map from the computational domain into the physical
domain, see e.g. Ciarlet
http://www6.cityu.edu.hk/rcms/publications/ln9.pdf

We are currently implementing the differential geometry terms as
Expressions but it should also be possible to derive all of the
differential geometry ingredients using UFL directly from the UFL
expression of the chart/map.

Perhaps the following code can help you, it solves Naghdi's non-linear
shell equations on a 2d manifold embedded in 3d space with cylindrical
coordinates:

https://bitbucket.org/unilucompmech/fenics-shells/src/b519f8815e9b669a8d70b6f3d125de6fc77f2f1a/demo/naghdi/demo_naghdi-arnoldbrezzi-batheexample.py?at=master&fileviewer=file-view-default

Jack S. Hale
Research Scientist
University of Luxembourg
Jack Hale

> To post to this group, send email to fenics-...@googlegroups.com.

> To view this discussion on the web visit

> https://groups.google.com/d/msgid/fenics-support/f71f6934-d8de-43dc-8ddd-b89ffc72975c%40googlegroups.com.

[Åsmund Hjulstad]

The current domain is just  1D

mesh = IntervalMesh(1000,ri,ro)

I just want to make sure I get the rheology descriptions right, and am doing this 1D test where there are exact solutions. The exact solution is derived using omega(r) instead of u(r), but that shouldn't cause any problems when accounted for properly.

Anyways, thanks for help. I will post running code in FEniCS Q&A, as Jack Hale suggested just now.

- Åsmund 

To unsubscribe from this group and stop receiving emails from it, send an email to fenics-support+unsubscribe@googlegroups.com.
To post to this group, send email to fenics-support@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/fenics-support/f71f6934-d8de-43dc-8ddd-b89ffc72975c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

mvh,
Åsmund Hjulstad

[Åsmund Hjulstad]

Just to followup with what I did to make it work, or at least get my solutions to match. (just in case someone searches for this)

I consider the flow velocity as a vector field in the derivations.   

- Div ( Grad ( uvec ) ) = 0   

Multiply with vvec, integrate by parts, u is constant on boundary.  The : (colon) is elementwise multiplication and then summation of all elements.

Int ( Grad ( uvec ) : Grad ( vvec )  dx )  = 0    

Then I replace Grad ( uvec ) with the expressions for cylindrical coordinates and uvec = {0, u(r), 0}. which is

0 -(u[r]/r)     0

u'[r] 0 0

0 0 0

and replace dx with r dr *2pi * 1   (and the 2*pi*1 goes away), which then gives me 

a_lin = mu_const * (u*v/r**2 + Dx(u,0)*Dx(v,0))*r*dx   # from my fenics script.

At the point I am interested in, it matches the exact solution to the first 8 digits. 

To post to this group, send email to fenics-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/fenics-support/f71f6934-d8de-43dc-8ddd-b89ffc72975c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

mvh,
Åsmund Hjulstad

--

mvh,
Åsmund Hjulstad