How to handle boundary conditions in Crank-Nicolson solution of IVP-BVP?

Question

I'm trying to solve the PDE for $c(r,t)$ $$c_t=(1/r)(rJ)_r$$ using Crank-Nicolson, and I'm having difficulty with the boundary conditions. $J$ is the flux, the initial condition is $c(0,r)=c_{init}$, and the flux is 0 at both boundaries $r=r_0$ and $r=r_1$. The flux is $J=rDc_r-s\omega^2r^2c$; $D$ and $s$ may be functions of $c$, but are often constants, as they are treated in this example.

So for the Crank-Nicolson scheme I use these approximations: $$c_t=(c_{n+1,j}-c_{n,j})/\Delta t$$ $$c_{rr}=\left({1\over 2}\right)\left[{c_{n,j+1}-2c_{n,j}+c_{n,j-1}\over (\Delta r)^2}+{c_{n+1,j+1}-2c_{n+1,j}+c_{n+1,j-1}\over (\Delta r)^2}\right]$$ $$c_r=\left({1\over 2}\right)\left[{c_{n,j+1}-c_{n,j-1}\over \Delta r}+{c_{n+1,j+1}-c_{n+1,j-1}\over \Delta r}\right]$$ $$c=\left({1\over 2}\right)\left[c_{n,j}+c_{n+1,j}\right]$$ where $n=$ the time step ($t=t_0+n\Delta t$), and $j=$ the position ($r=r_0+j\Delta r$, $j=0,...,J$).

For the boundary conditions, I want $J(r_0,t)=0$ and $J(r_1,t)=0$ for all $t$.

My first attempt was to use "ghost points" at $j=-1$ and $j=J+1$. I used the centered difference at $j=0$ and $j=J$, respectively, for the derivative $c_r$, and the values $c_{n,0}$ and $c_{n,J}$: $$c_r={c_{n,1}-c_{n,-1}\over 2\Delta r}$$ $$c_r={c_{n,J+1}-c_{n,J-1}\over 2\Delta r}$$ Substituting these approximations into the expression for the flux $J$ yields $c_{n,-1}$ in terms of $c_{n,0}$ and $c_{n,1}$, and $c_{n,J+1}$ in terms of $c_{n,J-1}$ and $c_{n,J}$.

With the ghost points eliminated, the $J+1 \times J+1$ tridiagonal matrices ${\bf A}$ and ${\bf B}$ are used to solve for the vector ${\bf c}_{n+1}$ from ${\bf c}_n$ for each time step.

Yet, when the numerical solution is calculated, instead of the expected exponential curve, becoming large near $r=r_1$, $c$ becomes tiny everywhere, and the values over time make it apparent that there is flux out at both boundaries. That is the issue I am grappling with.

I thought about using a Crank-Nicolson average to approximate $c_r$ at the boundaries, but then I get four additional unknowns, and I don't know how to deal with this.

Any suggestions?

Answer 1

Let me give an answer that is a general comment on prescribed zero flux for advection-diffusion (or convection-diffusion) PDE that is an important topic and it might be (but not necessary) the problem in your situation. Zero flux boundary condition is very non-standard in this case, so you might check if there is no misinterpretation of your problem.

Your flux $J$ can be split into two parts: $$ J=J_{diff} + J_{adv} = - r D c_r + s \omega^2 r^2 c $$ Let, formally, $s \omega^2 = 0$ then $J$ is purely diffusive flux and zero flux boundary conditions at $r_0$ and $r_1$ make sense as it is parabolic PDE. The situation is different if you consider for a moment that $D=0$ and you have only advective flux with "speed" $v(r)=s \omega^2 r^2$ and your PDE would be hyperbolic. For this type of problems one prescribes boundary conditions only for "inflow" boundaries, when the flow enters the domain, in your case when $v(r_0)>0$ or $v(r_1)<0$. In the opposite case for so-called "outflow" boundaries (consider the opposite inequalities) NO boundary conditions are required as the flow simply leaves the domain.

Let me consider your constants $s$ and $\omega$ being positive. Typical boundary conditions of advection for inflow are Dirichlet boundary conditions, i.e. in your case $c(r_0,t) = c_{given}(t)$. No boundary conditions are required for the advective flux at $r=r_1$.

Now if all constants are non-zero, a standard treatment for advection-diffusion flux at boundary is that you consider them separately. For instance prescribe the zero diffusive flux everywhere and prescribe the concentration in the advective flux only at inflow boundary.

If you prescribed in your code the zero advective-diffusive flux everywhere it still make sense for the inflow, because it is equivaklent to zero diffusive flux and zero concentration at advective flux. But for the outflow it can lead to non-physical behavior, because in fact what you are prescribing is that the diffusive flux is equal minus advective flux. This is especially serious if you have advection dominated problem when so called Peclet number is large, i.e. $D<s \omega^2 r_1^2$.

If this still does not solve your problem, give us a few pictures of time development for your numerical solution, just one or two time steps from constant initial condition.

P.S. In general the zero flux B.C. in Crank-Nicolson method are considered for each time point separately, so the ghost values at $n$ and $n+1$ are eliminated separately.

Answer 2

The main reason my solution did not work was that

$$c_r=\left({1\over 2}\right)\left[{c_{n,j+1}-c_{n,j-1}\over \Delta r}+{c_{n+1,j+1}-c_{n+1,j-1}\over \Delta r}\right]$$ is incorrect. For the centered-difference approximation for $c_r$, the denominator should be $2\Delta r$, not $\Delta r$. There is an additional factor of $1\over 2$ from the Crank-Nicolson stencil.

Making this correction fixed most, but not all, of the mass loss. There is still some problem, most likely in the application of the boundary conditions.