# A 2D Finite Difference formulation in polar coordinates.

1. May 2, 2017

### maistral

So I have this PDE:

d2T/dr2 + 1/r dT/dr + d2T/dθ2 = 0.

How do I implement dT/dr || [r = 0] = 0? Also, what should I do about 1/r?

This is actually the first time I am going to attack FDF in polar/cylindrical coordinates. I can finite-difference the base equation fairly decently; I am just having a hard time in implementing the derivative boundary condition at r = 0. Can someone give me an idea what to do?

2. May 2, 2017

### bigfooted

In practice, no computations are made on the line r=0. I assume you know how to discretize and how to obtain your finite difference matrix-vector system.
Suppose you discretize and you have N+1 nodes from j=0..N in the radial direction where the nodes with index j=0 correspond to r=0. This leads to a matrix-vector system of size (N+1). The Neumann boundary condition (first order) can be written as $\frac{T_1 - T_0}{\Delta R_1} = 0$ or: $T_0 = T_1$

The solution on the computational nodes at r=0 are known when $T_1$ is known. So you replace every occurence of $T_0$ with $T_1$ and you delete the first row and column of the discretization matrix and you solve the size N matrix system. When the solution $T_1 .. T_N$ is found, you simply add $T_0 = T_1$ to your solution vector.

3. May 3, 2017

### maistral

I managed to do what you suggested in a 1-D system (in radial directions only). I wanted to stencil my problem in a 2-D system (radial and angular directions). This is my problem:

I know the number of points is... I'm doing this to kill off certain doubts.

I stencil around T0, by using an energy balance approach. I got something like T0 = (2T1 + 2T3) / 4. Is this correct? This is actually the problem I am facing.

If I implement the boundary condition that dT/dr [r = 0] = 0, that would mean either (T1-T0)/2Δr or (T3-T0)/2Δr is 0 right? Which do I replace? T1 or T3?

Last edited: May 3, 2017
4. May 3, 2017

### pasmith

Physically $r = 0$ is an interior point of the domain. You can't impose a condition on the value of $T$ or $\partial T/\partial r$ there.

5. May 3, 2017

### bigfooted

If he is simulating only a quarter of a circle, then r=0 is on the boundary of the domain (it is not so clear in the text, but the problem description says $\Delta \phi = \pi / 2$ ).

Then to implement the Neumann boundary condition:
on the horizontal line, T1=T0
on the 45-degree line, T2=T0
on the vertical line, T3=T0
So T0=T1=T2=T3