I am modeling a system using the convection-diffusion equation on a 2D, radially symmetric space. I wanted to do some sanity checks to make sure I am coding it correctly. I set up a situation where I believe there should be no net flux across the boundaries, but it seems like the solutions I am getting does have flux. I asked about my Mathematica code implementation on Mathematica Stack Exchange, but it seems like I am missing some important mathematical details that is causing me to overlook something.
For the equation $$\frac{\partial c}{\partial t}=\alpha\nabla^2c-\nabla\cdot(\mathbf vc)+S$$
I thought the flux across a boundary (independent of coordinates) was given by $(\mathbf vc-\alpha\nabla c)\cdot\hat n$, where $\hat n$ is the unit vector normal to the boundary surface. In polar coordinates with a radial velocity $\mathbf v=v_r\hat r$ field and a region that is an annulus, this could be then written as $v_rc-\alpha\cdot\partial c/\partial r$ on the inner and outer circles.
Therefore, if we use a function $v_r$ that goes to $0$ on the boundaries, and if we impose no flux boundary conditions $\partial c/\partial r=0$, shouldn't we have no net flux into or out of the annular region? I would think so, but my simulations as well as this answer on my Mathematica SE post seems to suggest otherwise.
What else needs to be true of this radially symmetric system so that there is actually no flux at the boundaries? Or does the choice of $v_r$ within the region cause unintended sources to pop up?
If $S = 0$, then mathematically, the solution can also be of the form $c(p) = e^{<k,p>+\omega t}$ with $\omega = \alpha (t \nabla^2 \omega + (k + t \nabla \omega).(k + t \nabla \omega)) - v (k_1+k_2+k_3 + t \nabla . \omega ) - \nabla.v$.
So if $v = constant$ (may be $v = 0$), function $c(p)$ will be uniform across a plane $<k,p> = constant$ and may not be uniform on a surface of sphere and may not die out.
Please add more details on actual $v,S$ so that it will be useful for others to answer.
Cheers man !