# Finite Difference advection-reaction-diffusion in spherical coordinates, problem with Diffusion

I have a problem with the following equation:

$$V_r\frac{\partial C_A}{\partial r} + \frac{V_\theta}{r}\frac{\partial C_A}{\partial \theta} = \frac{2}{Pe_A} \left[ \frac{\partial ^2 C_A}{\partial r ^2} + \frac{2}{r} \frac{\partial C_A}{\partial r} - k C_A\right]$$

I'm trying to solve this with finite difference technique. My discretization :

$$V_r \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} + \frac{V_\theta}{r} \frac{A_{i,j+1}-A_{i,j-1}}{2\Delta \theta} - \frac{2}{Pe} \left [ \frac{A_{i+1,j}-2A_{i,j}+A_{i-1,j}}{\Delta r ^2}+\frac{2}{r} \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} - kA_{i,j}\right] = 0$$

I first tried it with this, but because I used central-difference in doesn't work for $Pe(cell)>2$, I had strange oscillations at the boundary. So I thought I use a upwind scheme now.

$V_r$ and $V_\theta$ are given through polynoms as a function ($r, \theta$).

So for $\theta<90 °$ I Use second order upwind for the convection term $V_r > 0$.

And for $\theta >90 °$ in the other direction. $V_\theta$ is always $<0$, so I only use this form for this.

My problem is, for high Peclet number and low reaction constants, the concentration is rising up where it shouldn't. This happens in a region, where the convection term gets stronger.

What am I supposed to do ? Is there a difference about using upwind in polar coordinates ?

Edit :

I want go more in detail. This is my coordinate system :

The problem appears at the zero degress line. I also have seen that the problem is strongly dependent on $V_\theta$. If I multiply $V_\theta$ with 0.9 for example there is no problem. The concentration is not rising. I've looked at the speed vectors in more detail. The problem is vanishing when the speed vectors near the zero degree line are horizontal or show away from it. I have Neumann boundary condition there ( also at 180°). ($\frac{\partial C_A}{\partial \theta}=0$). Now I also have seen that for the regions where this "pseudo-source" exist the boundary condition isn't met perfectly. There is a gradient. So I am thinking maybe this acts like a "pseudo-source ". I also have drichilet boundary condition for r=1 and r-->infinite. ( My calculation begins at r=1).

My grid number is the following ( for example I have three radial and five angular steps ). My notation is c(i,j), i is the radial direction and j the angular direction. So say for i=1, it is c(1,j) with j=0...4 from phi=0...180°. Let's have a look for zero degress with constants l_i : $c10*l1+...c12*l2+c14*l4=0$.

I need no ghoistpoints here. I just sad c12=c14 and rearranged the equation to $c10*l1...+c12(l2+l4)$. For 180 ° I need ghoistpoints because of second order upwind scheme I have c(i,j-2) which is outside my domain, but at 180 ° there is no problem.

Should I implement Neumann on another way ? Or do you have any other suggestions ?

• I typically run into these problems because i want to go too fast with my implementation. Sometimes its better to build the simplest version of a code and add complexity. In your case i would advise: 1. write a code for a cartesian grid, 2. just consider diffusion, forgot about the advection and reaction term for now. When you have that working you should have a linear concentration profile. Then add complexity by moving to a spherical grid and then add a reaction term and finally the convection term. For me this is much less prone to mistakes, theoretical and practical ones. – nluigi Sep 26 '16 at 9:24