2D Finite Difference formulation in polar coordinates.

[maistral]

So I have this PDE:

d²T/dr² + 1/r dT/dr + d²T/dθ² = 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?

 

[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 occurrence 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.

 

[maistral]

Hi. Thanks for replying.

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 T₀, by using an energy balance approach. I got something like T₀ = (2T₁ + 2T₃) / 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 (T₁-T₀)/2Δr or (T₃-T₀)/2Δr is 0 right? Which do I replace? T₁ or T₃?

 

[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.

 

[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