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

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.

• How did you derive the pde with eliminated [CurlyPhi]-coordinate? May 19 at 13:10
• @UlrichNeumann It's radially symmetric May 19 at 13:21
• @BioPhysicist That's already clear. May 19 at 15:37
• @UlrichNeumann Then I guess I don't understand your question. $c$, $\mathbf v$, and $S$ have no $\theta$ dependence, so the pde has no $\theta$ in it. May 19 at 15:41

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. • "It follows that with div -r/|r| = -1/r there is absorption by if its not compensated by the area factor r" Isn't the flux on the boundary $v\cdot c-\alpha\cdot\partial c/\partial r$? So if $\mathbf v=0$ and $\partial c/\partial r=0$ that should be no flux. Even if not, isn't the $1/r$ taken care of on the flux integral across the border as $\text d\ell=r\text d\theta$? May 19 at 21:04
• You observe from the picture, that with your parameter, the particles collect at 0.3 where diffusion force, driving force , constant local production 0.75 and deletion -0.75c . seem to reach its equilibrium. With drift term 0 the area of constant production prevails at maximal radius by area. May 20 at 4:41
• Sorry, I'm still not following. If I have set up no transfer across the boundary, and there is production proportional to $1-c$, it should just approach $c=1$ everywhere. Is it an artifact of the step functions? May 20 at 4:58
• I think your approach should be formualted in cartesian coordinates, first. Boundary conditions on circles pose no problem for numeric integration. The problematic term is the production constant in $dr$. In $2d$ ist should be $r dr$ May 20 at 5:15
• The flux being $v_\bot\cdot c-\alpha\cdot\partial c/\partial n$ on the boundary is independent of the coordinates. The production doesn't cause a flux across the boundary. May 20 at 5:52