2
$\begingroup$

Governing equations in polar or cylindrical coordinates often have terms with $\frac{1}{r}$ involved. At $r = 0$, such terms blow up to become a "singularity." The Cartesian version of such governing equations will not have these "singularities" at $(x,y)=(0,0)$. Therefore, the singularity at $r = 0$ must be numerical.

How do you handle the singularity at $r = 0$ in polar or cylindrical coordinates via e.g., finite difference?

To break this down:

  • How are the radial coordinates discretized and differentiated?
  • Is there is a transformation of variables of $(r,\theta)$ involved?
  • How does any of this prevent blow-up at $r = 0$?

Let's assume the general case where there is no (anti-) symmetry to exploit and the governing equations rely on $r$ and $\theta$ explicitly, so the polar problem is at least two-dimensional. The governing partial differential equations could be the heat equation, Navier-Stokes equations, or Laplace's equation.

Some related answers are:

Neumann Boundary Condition at r=0 in Polar Coordinates (Numerical BCs)

Finite difference methods in cylindrical and spherical co-ordinate systems

$\endgroup$
9
  • 1
    $\begingroup$ Using some small nonzero $r_{min}$ instead of $r=0$ is an easy solution acceptable for many applications. $\endgroup$ Commented Aug 2, 2022 at 1:03
  • $\begingroup$ So, the radial mesh would be something like {$r_{min}$, $r_{min} +\Delta r$, $+ ...$, $+ r_{max}$} or would it be something like {$-r_{max}$, $+...$, $-r_{min}$, $r_{min}$, $+ ...$, $+ r_{max}$}? If it is the former, does $r_{min}$ require a "boundary" condition? What would that be? If it is the latter, is there a transformation of variables involved (I'm thinking of needing to change the definition of $\theta$)? $\endgroup$
    – Steve M
    Commented Aug 2, 2022 at 4:06
  • $\begingroup$ A related question on Math.SE $\endgroup$
    – Ruslan
    Commented Aug 2, 2022 at 12:54
  • $\begingroup$ One could use the grid, e.g., $\epsilon, \Delta r, 2 \Delta r, ..., r_{max}$. At $r_{min}=\epsilon$ the boundary condition would enforce Dirichlet or Neumann (or more complex) conditions that you'd use with $r_{min}$=0. $\endgroup$ Commented Aug 2, 2022 at 14:27
  • 1
    $\begingroup$ The problem you're having is called the "coordinate singularity" and is, among other places, addressed in this paper doi.org/10.1006/jcph.1999.6382 $\endgroup$
    – user20857
    Commented Aug 2, 2022 at 17:03

1 Answer 1

5
$\begingroup$

Let me show this specifically for the finite element discretization of the Laplace equation: $$ -\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) = f(r,\theta). $$ The finite element method is based on the weak formulation, which you obtain by multiplying the equation by a test function $\varphi(r,\theta)$ and integrating over the domain, then integrate by parts. You'd think that that leads to $$ \int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, d\theta dr, $$ followed by integration by parts.

But you don't have to do it that way -- the above approach simply assumes that $r$ and $\theta$ are independent variables. Instead, you use the proper area element, which is $2\pi r \, dr \, d\theta$ and you get the following form instead: $$ \int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, 2\pi r \, d\theta dr \\ = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, 2\pi r \, d\theta dr. $$ The $2\pi$ cancel on the two sides, and so what you have left is this: $$ \int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, r \, d\theta dr. $$ After integration by parts, you then have the following (omitting boundary terms for simplicity): $$ \int_0^R \int_0^{2\pi} r \frac{\partial \varphi(r,\theta)}{\partial r} \frac{\partial u(r,\theta)}{\partial r} + \frac{1}{r} \frac{\partial \varphi(r,\theta)}{\partial \theta} \frac{\partial u(r,\theta)}{\partial \theta} \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, r \, d\theta dr. $$ As you can see, at least from in front of the the $r$ derivatives, the singular weight has disappeared. If you choose appropriate quadrature rules, then the quadrature points for the second term will never lie at $r=0$, and so the still singular weight will not matter in the $\theta$-derivative term.

$\endgroup$
10
  • $\begingroup$ How will the $1/r$ singular weight be made irrelevant after integration? To generalize, the method to treat the singularity is to somehow remove the $r = 0$ meshpoint from consideration? $\endgroup$
    – Steve M
    Commented Aug 2, 2022 at 4:29
  • $\begingroup$ I've never seen things done this way before. Could you point me to what a "proper" area element is? It seems to be what you'd get when you transform $dx dy$ into polar. $\endgroup$
    – NNN
    Commented Aug 2, 2022 at 14:12
  • $\begingroup$ @Nachiket Yes, exactly. It's just a coordinate transform. The integral is over a domain in $R^2$, and the integral should be the same whether you express it in $x-y$ coordinates or $r-\theta$ coordinates. The area element that includes the factor $2\pi r$ does that. $\endgroup$ Commented Aug 3, 2022 at 2:40
  • $\begingroup$ @SteveM We can't compute integrals in practice. We have to use quadrature. The question is what quadrature formula you choose. $\endgroup$ Commented Aug 3, 2022 at 2:41
  • 1
    $\begingroup$ @Nachiket You can, but the question is whether you are able to integrate the other part by parts in a convenient way so that the resulting bilinear form becomes symmetric and positive definite. You will have noticed that that is the case with the $2\pi r$ weight but that that would not have been the case when using no weight at all. $\endgroup$ Commented Aug 5, 2022 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.