next up previous main_module_page.gifCFD_module_list.gif
Next: 2. Parabolic Equations Up: 1. Introduction Previous: 1.3 Finite differences

1.4 Finite differences in polar coordinates


Example 5   Consider $\nabla^2 U = 0$ in polar coordinates, i.e.,

\begin{displaymath}\frac{\partial^2 U}{\partial r^2} + \frac{1}{r} \frac{\partia...
...}
+ \frac{1}{r^2} \frac{\partial^2 U}{\partial \theta^2} = 0
\end{displaymath}

(See Figure 1.7.)

  
Figure 1.7: Mesh in polar coordinates
\begin{figure}
\begin{center}
\includegraphics[scale=0.9]{fd_fig_4_5.eps}
\end{center}
\end{figure}

Then we have

\begin{displaymath}\frac{u_{p-1,q} - 2 u_{p,q} + u_{p+1,q}}{h^2}
+ \frac{1}{p h...
...p,q} + u_{p,q+1}}{\left( \Delta \theta \right)^2} \right] = 0
\end{displaymath}

etc.


Example 6   In cylindrical coordinates the PDE $\frac{\partial U}{\partial t} = \nabla^2U$ can be written

 \begin{displaymath}\frac{\partial U}{\partial t} = \frac{\partial^2 U}{\partial ...
...2 U}{\partial \theta^2} + \frac{\partial^2 U}{\partial z^2}
\end{displaymath} (1.3)


  
Figure 1.8: Cylindrical grid system
\begin{figure}
\begin{center}
\includegraphics[scale=0.9]{DBI_ex1_fig1.eps}
\end{center}
\end{figure}

We can set up a grid system as shown with $\Delta \theta = 2 \pi / N$, $\Delta r = 1/M$ and represent derivatives by either central or forward differences. How do we deal with the point $r=0$?


Solution
If the solution is independent of $z$ then equation (1.3) becomes

 \begin{displaymath}\frac{\partial U}{\partial t} = \frac{\partial^2 U}{\partial ...
... r} +
\frac{1}{r^2} \frac{\partial^2 U}{\partial \theta^2}
\end{displaymath} (1.4)

and in Cartesian coordinates it may be written

 \begin{displaymath}\frac{\partial U}{\partial t} = \frac{\partial^2 U}{\partial ...
...al y^2}
\left( + \frac{\partial^2 U}{\partial z^2} \right)
\end{displaymath} (1.5)

We write

\begin{displaymath}\Delta \theta = \frac{2 \pi}{N}, \quad
\Delta r = \frac{1}{M}, \quad
\Delta t
\end{displaymath}

and

\begin{displaymath}U \left( r, \theta, t \right) \equiv U \left( p \Delta r, q \Delta \theta, n \Delta t \right)
\equiv U_{p,q,n}
\end{displaymath}

At a general point

\begin{alignat*}{2}
\frac{\partial U}{\partial t} & = \frac{U_{p,q,n+1} - U_{p,...
...)^2} \hspace{1cm} &
& \hbox{at $r$ , $\theta$ , $t$\space etc.}
\end{alignat*}
For nonzero values of $r$ there are no difficulties using the above formula but how do we deal with $\frac{\partial^2 U}{\partial r^2}$, $\frac{1}{r} \frac{\partial U}{\partial r}$ and $\frac{\partial^2 U}{\partial \theta^2}$ at $r=0$ ?
e.g.

\begin{displaymath}\frac{1}{r} \frac{\partial U}{\partial r} = \hbox{?} \quad \hbox{at} \quad r=0
\end{displaymath}

Therefore at $r = 0$ we use the PDE in Cartesian form, i.e. (1.5)

\begin{displaymath}\frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\part...
...artial^2 U}{\partial z^2} \right) \quad \hbox{for $\nabla^2$}
\end{displaymath}


  
Figure 1.9: 2D grid system
\begin{figure}
\begin{center}
\includegraphics[scale=0.9]{DBI_soln1_fig2.eps}
\end{center}
\end{figure}

Thus in 2D at $r=0$ we have

 \begin{displaymath}\begin{array}{r@{}l}
\nabla^2 U & {}= \displaystyle\frac{u_1...
...- u_0 \right] \bigg/ \left( \Delta r \right)^2
\end{array}
\end{displaymath} (1.6)

However, the $x$ and $y$ axes are two perpendicular axes and the generalisation of equation (1.6) gives

\begin{displaymath}\nabla^2 U = \frac{4 \left[ u_m - u_0 \right]}{\left( \Delta r \right)^2}
\end{displaymath}

where $u_m$ is the mean value of $u$ along $r = \Delta r$, i.e., on $q=1$.
When $\frac{\partial}{\partial \theta} \equiv 0$, equation (1.4) becomes

\begin{displaymath}\frac{\partial U}{\partial t} = \frac{\partial^2 U}{\partial r^2} + \frac{1}{r} \frac{\partial U}{\partial r}
\end{displaymath}

and the question arises as to how do we deal with the term $\frac{1}{r} \frac{\partial U}{\partial r}$ at $r=0$ ?
The answer is that we use l'Hopital's rule, i.e.

\begin{displaymath}\frac{1}{r} \frac{\partial U}{\partial r} = \frac{\frac{\part...
...1} = \frac{\partial^2 U}{\partial r^2}
\quad \hbox{at $r=0$}
\end{displaymath}

(N.B. $\frac{\partial U}{\partial r} = 0$ at $r=0$ if problem is symmetrical)
Therefore, at $r=0$

\begin{displaymath}\frac{\partial U}{\partial t} = 2 \frac{\partial^2 U}{\partial r^2}
\hspace{2cm} U(r,t) = U_{p,n}
\end{displaymath}

and

\begin{displaymath}\frac{\partial^2 U}{\partial r^2} = \frac{u_{-1,n} - 2u_{0,n} + u_{1,n}}{\left( \Delta r \right)^2}
\end{displaymath}


\begin{displaymath}\frac{\partial U}{\partial r} = 0 \quad \Longrightarrow \quad
\frac{u_{-1,n} - u_{1,n}}{2 \Delta r} = 0
\end{displaymath}

Therefore

\begin{displaymath}\boxed{
\frac{\partial^2 U}{\partial r^2} = \frac{2 \left( u...
...t)}{\left( \Delta r \right)^2}
\quad \hbox{at} \quad r=0
}
\end{displaymath}


[You can now do Question 1 of the Examples]


next up previous main_module_page.gifCFD_module_list.gif
Next: 2. Parabolic Equations Up: 1. Introduction Previous: 1.3 Finite differences
Maintained by: Dr S D Harris
2003-06-16