Conservation of a physical quantity when using Neumann boundary conditions applied to the advection-diffusion equation

Question

I don't understand the different behaviour of the advection-diffusion equation when I apply different boundary conditions. My motivation is the simulation of a real physical quantity (particle density) under diffusion and advection. Particle density should be conserved in the interior unless it flows out from the edges. By this logic, if I enforce Neumann boundary conditions the ends of the system such as $\frac{\partial \phi}{\partial x}=0$ (on the left and the right sides) then the system should be "closed" i.e. if the flux at the boundary is zero then no particles can escape.

For all the simulations below, I have applied the Crank-Nicolson discretization to the advection-diffusion equation and all simulation have $\frac{\partial \phi}{\partial x}=0$ boundary conditions. However, for the first and last rows of the matrix (the boundary condition rows) I allow $\beta$ to be changed independently of the interior value. This allows the end points to be fully implicit.

Below I discuss 4 different configurations, only one of them is what I expected. At the end I discuss my implementation.

Diffusion only limit

Here the advection terms are turned off by setting the velocity to zero.

Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points

The quantity is not conserved as can be seen by the pulse area reducing.

Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries

By using fully implicit equation on the boundaries I achieve what I expect: no particles escape. You can see this by the area being conserved as the particle diffuse. Why should the choice of $\beta$ at the boundary points influence the physics of the situation? Is this a bug or expected?

Diffusion and advection

When the advection term is included, the value of $\beta$ at the boundaries does not seem to influence the solution. However, for all cases when the boundaries seem to be "open" i.e. particles can escape the boundaries. Why is this the case?

Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points

Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries

Implementation of the advection-diffusion equation

Starting with the advection-diffusion equation,

$ \frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2} + \boldsymbol{v}\frac{\partial \phi}{\partial x} $

Writing using Crank-Nicolson gives,

$ \frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right] $

Note that $\beta$=0.5 for Crank-Nicolson, $\beta$=1 for fully implicit, and, $\beta$=0 for fully explicit.

To simplify the notation let's make the substitution,

$ s = D\frac{\Delta t}{(\Delta x)^2} \\ r = \boldsymbol{v}\frac{\Delta t}{2 \Delta x} $

and move the known value $\phi_{j}^{n}$ of the time derivative to the right-hand side,

$ \phi_{j}^{n+1} = \phi_{j}^{n} + s \left( 1-\beta \right) \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + s \beta \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) + r \left( 1 - \beta \right) \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + r \beta \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) $

Factoring the $\phi$ terms gives,

$ \underbrace{\beta(r - s)\phi_{j-1}^{n+1} + (1 + 2s\beta)\phi_{j}^{n+1} -\beta(s + r)\phi_{j+1}^{n+1}}_{\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}}} = \underbrace{ (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}}_{\boldsymbol{M\cdot}\boldsymbol{\phi^n}} $

which we can write in matrix form as $\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}} = \boldsymbol{M}\cdot\boldsymbol{\phi^{n}}$ where,

$ \boldsymbol{A} = \left( \begin{matrix} 1+2s\beta & -\beta(s + r) & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & \beta(r-s) & 1+2s\beta \\ \end{matrix} \right) $

$ \boldsymbol{M} = \left( \begin{matrix} 1-2s(1-\beta) & (1 - \beta)(s + r) & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & (1 - \beta)(s - r) & 1-2s(1-\beta) \\ \end{matrix} \right) $

Applying Neumann boundary conditions

NB is working through the derivation again I think I have spotted the error. I assumed a fully implicit scheme ($\beta$=1) when writing the finite difference of the boundary condition. If you assume a Crank-Niscolson scheme here the complexity become too great and I could not solve the resulting equations to eliminate the nodes which are outside the domain. However, it would appear possible, there are two equation with two unknowns, but I couldn't manage it. This probably explains the difference between the first and second plots above. I think we can conclude that only the plots with $\beta$=0.5 at the boundary points are valid.

Assuming the flux at the left-hand side is known (assuming a fully implicit form),

$ \frac{\partial\phi_1^{n+1}}{\partial x} = \sigma_L $

Writing this as a centred-difference gives,

$ \frac{\partial\phi_1^{n+1}}{\partial x} \approx \frac{\phi_2^{n+1} - \phi_0^{n+1}}{2\Delta x} = \sigma_L $

therefore, $ \phi_0^{n+1} = \phi_{2}^{n+1} - 2 \Delta x\sigma_L $

Note that this introduces a node $\phi_0^{n+1}$ which is outside the domain of the problem. This node can be eliminated by using a second equation. We can write the $j=1$ node as,

$ \beta(r - s)\phi_0^{n+1} + (1+2s\beta)\phi_1^{n+1} - \beta(s+r)\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} $

Substituting in the value of $\phi_0^{n+1}$ found from the boundary condition gives the following result for the $j$=1 row,

$ (1+2s\beta)\phi_1^{n+1} - 2s\beta\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} + 2\beta(r-s)\Delta x\sigma_L $

Performing the same procedure for the final row (at $j$=$J$) yields,

$ -2s\beta\phi_{J-1}^{n+1} + (1+2s\beta)\phi_J^{n+1} = (1-\beta)(s - r)\phi_{J-1}^{n} + (1 - 2s(1-\beta))\phi_{J}^{n} + 2\beta(s+r)\Delta x\sigma_R $

Finally making the boundary rows implicit (setting $\beta$=1) gives,

$ (1+2s)\phi_1^{n+1} - 2s\phi_2^{n+1} = \phi_{j-1}^{n} + 1\phi_{j}^{n} + 2(r-s)\Delta x\sigma_L $

$ -2s\phi_{J-1}^{n+1} + (1+2s)\phi_J^{n+1} = \phi_{J}^{n} + 2(s+r)\Delta x\sigma_R $

Therefore with Neumann boundary conditions we can write the matrix equation, $\boldsymbol{A}\cdot\phi^{n+1} = \boldsymbol{M}\cdot\phi^{n} + \boldsymbol{b_N}$,

where,

$ \boldsymbol{A} = \left( \begin{matrix} 1+2s & -2s & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & -2s & 1+2s \\ \end{matrix} \right) $

$ \boldsymbol{M} = \left( \begin{matrix} 1 & 0 & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & 0 & 1 \\ \end{matrix} \right) $

$ \boldsymbol{b_N} = \left( \begin{matrix} 2 (r - s) \Delta x \sigma_L & 0 & \ldots & 0 & 2 (s + r) \Delta x \sigma_R \end{matrix} \right)^{T} $

My current understanding

• I think the difference between the first and second plots is explained by noting the error outlined above.

• Regarding the conservation of the physical quantity. I believe the cause is that, as pointed out here, the advection equation in the form I have written it doesn't allow propagation in the reverse direction so the wave just passes through even with zero-flux boundary conditions. My initial intuition regarding conservation only applied when advection term is zero (this is solution in plot number 2 where the area is conserved).

• Even with Neumann zero-flux boundary conditions $\frac{\partial \phi}{\partial x} = 0$ the mass can still leave the system, this is because the correct boundary conditions in this case are Robin boundary conditions in which the total flux is specified $j = D\frac{\partial \phi}{\partial x} + \boldsymbol{v}\phi = 0$. Moreover the Neunmann condition specifies that mass cannot leave the domain via diffusion, it says nothing about advection. In essence what we have hear are closed boundary conditions to diffusion and open boundary conditions to advection. For more information see the answer here, Implementation of gradient zero boundary conditon in advection-diffusion equation.

Would you agree?

• A tutorial on how to implement Robin boundary conditions.

Answer 1

I think that one of your problems is that (as you observed in your comments) Neumann conditions are not the conditions you are looking for, in the sense that they do not imply the conservation of your quantity. To find the correct condition, rewrite your PDE as

$$ \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\partial x} + v \phi \right) + S(x,t) .$$

Now, the term that appears in parentheses, $ D\frac{\partial \phi}{\partial x} + v \phi = 0 $ is the total flux and this is the quantity that you must put to zero on the boundaries to conserve $\phi$. (I have added $S(x,t)$ for the sake of generality and for your comments.) The boundary conditions that you have to impose are then (supposing your space domain is $(-10,10)$)

$$ D\frac{\partial \phi}{\partial x}(-10) + v \phi(-10) = 0 $$

for the left side and

$$ D\frac{\partial \phi}{\partial x}(10) + v \phi(10) = 0 $$

for the right side. These are the so-called Robin boundary condition (note that Wikipedia explicitly says these are the insulating conditions for advection-diffusion equations).

If you set up these boundary conditions, you get the conservation properties that you were looking for. Indeed, integrating over the space domain, we have

$$ \int \frac{\partial \phi}{\partial t} dx = \int \frac{\partial}{\partial x} \left( D \frac{\partial \phi}{\partial x} + v \phi \right) dx + \int S(x,t) dx$$

Using integration by parts on the right hand side, we have

$$ \int \frac{\partial \phi}{\partial t} dx = \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(10) - \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(-10) + \int S(x,t) dx$$

Now, the two central terms vanish thanks to the boundary conditions. Integrating in time, we obtain

$$ \int_0^T \int \frac{\partial \phi}{\partial t} dx dt = \int_0^T \int S(x,t) dx dt$$

and if we are allowed to switch the first two integrals,

$$ \int \phi(x,T) dx - \int \phi(x,0) dx = \int_0^T \int S(x,t) dx$$

This shows that the domain is insulated from the exterior. In particular, if $S=0$, we get the conservation of $\phi$.