Implementing and testing boundary conditions in polar coordinates for convection-diffusion equation

Question

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.

The equation I am implementing involves a radial velocity field that depends on the dependent variable $c$ that I have forced to drop to $0$ at the boundaries using step functions, and a source term that goes to $0$ as $c$ approaches $1$: $$\mathbf v=(1-c)\cdot\Theta(9.5-r)\cdot\Theta(r-.0105)\cdot\hat r$$ $$S=0.75(1-c)$$

And the differential equation I am solving is

$$\frac{\partial c}{\partial t}=10\nabla^2c-10\nabla\cdot\left(\mathbf vc\right)+S$$

I impose no diffusive flux on the boundaries (which is an annulus with radii of $0.01$ and $10$) and initial condition of $c=0$. Since $\mathbf v=0$ and $\partial c/\partial r=0$ on the boundaries, I believe there should be no net flux across the boundaries.

My code implementation is as follows:

ndEQ = D[c[r, t], t] == 10.*Laplacian[c[r, t], {r, \[Theta]}, "Polar"] - 
10.*Div[{(1 - c[r, t])*UnitStep[9.5 - r]*UnitStep[r - .0105], 0}*
   c[r, t], {r, \[Theta]}, "Polar"] + .75*(1 - c[r, t]);

initialCondition = c[r, 0.] == 0.;

bcOuter = (D[c[r, t], r] /. r -> 10.) == 0.; 

bcInner = (D[c[r, t], r] /. r -> .01) == 0.; 

sol = NDSolveValue[{ndEQ, initialCondition, bcInner, bcOuter}, 
   c, {r, .01, 10.}, {t, 0, 5.}];

With no flux across the boundaries and a non-zero source term, I believe $c$ should actually (uniformly?) approach $1$ over time, but this is not the case:

Plot[Table[sol[r, t], {t, 0., 5., .5}], {r, .01, 10.},PlotRange -> {0, 1.1}, AxesLabel -> {"r", "c"}]

Am I missing something with how the boundary conditions are enforced using polar coordinates? I ran an analogous set of equations on a 1D space and the system behaved as I would expect (uniform increase to $1$), so I believe there is something I am missing in the polar coordinate system.

Answer 1

I replaced your step functions by HeavisidePi, the window function. Then the Div term produces DiracDelta (hope that Div in "Polar" handles distributions correctly: yes, checked ) and played around with prefactors of Pi and Delta.

As it seems, the Div term is absorbing with a limit at 0.3 .

      dEQ = D[ c[r, t], t] == 
      (10.*Laplacian[  c[r, t], {r, \[Theta]}, "Polar"] - 
       10.*Div[{(1 -  c[r, t]) HeavisidePi[(r/9.151 - 1/2) ], 0}* c[r, t], 
       {r, \[Theta]}, "Polar"] + 
       .75*(1 - c[r, t]));

I just noticed that in fact the vector field is radial constant. It follows that with div -r/|r| = -1/r there is absorption by if its not compensated by the area factor r. Finally, the question is if the area integral is constant in time.

I must admit, that I forgot the central point.

The $1-d$ diffuion equation with drift term $n/r$ is the radial part of the Bessel process, were the 1/r drift is the effect of increasing volume with radius.

n=1 is the statistical critical value of recurrence. For n=1 in free space, the particles visit any circle infinitly often, for n=2 (spherical) they never come back.

For reflecting walls the effect is that for n=0 the limit distribution is 1d-flat, for n=1 it is 2d-flat, for n=3 its a gaussian with reflections etc.

The 2d-polar Laplacian produces the correct $1/r \partial_r$ drift term due to the increasing volume with r of the 2d Wiener process.

But the polar $\nabla$, acting on a constant radial vector does, what one may call "change of the dimension of space continously" for the underlying Wiener process.

So expect, qualititatively crowding at the inner or outer boundary depending on the constant of the drift term $1/r$.

Speaking as a physicist: its probably not possible on the cell geometry to introduce a radial drift field, that changes the Bessel process dimension from 2 to 11.