# Numerically solving the Advection-diffusion equation with no-flux boundary condition leads to violation of mass conservation

I am trying to solve numerically the advection-diffusion equation of the following form

$$\frac{\partial C}{\partial t}=\alpha\frac{\partial^2 C}{\partial x^2}+\beta \frac{\partial C}{\partial x}$$

with no flux BCs namely

$$\alpha\frac{\partial C}{\partial x}+\beta C=0$$

at $$x=0$$ and $$x=L$$.

I solve this as one would normally solve the diffusion equation, namely by using finite differences and the forward Euler algorithm.

$$C_N(t+dt)=C_N(t)+dt\left(\alpha\frac{C_{N+1}(t)-2C_N(t)+C_{N-1}(t)}{dx^2}+\beta\frac{C_{N+1}(t)-C_{N-1}(t)}{2dx}\right)$$

The BC's I implement via the ghost point method i.e. $$N=n+1$$ and $$N=0$$ are excluded and eliminated via the boundary conditions. For $$N=1$$ we get

$$\alpha\frac{C_{2}-C_{0}}{2dx}+\beta C_1=0$$

or

$$C_{0}=C_{2}+\frac{2dx\beta}{\alpha}C_1$$

and for $$N=n$$

we get similarly

$$C_{n+1}=C_{n-1}-\frac{2dx\beta}{\alpha}C_n$$

which gives us the time-evolution equations at the boundary as

$$C_1(t+dt)=C_1(t)+dt\left(\alpha\frac{2C_{2}(t)-2C_1(t)+\frac{2dx\beta}{\alpha}C_1(t)}{dx^2}-\frac{\beta^2}{\alpha}C_1\right)$$

and

$$C_n(t+dt)=C_n(t)+dt\left(\alpha\frac{2C_{n-1}(t)-2C_n(t)-\frac{2dx\beta}{\alpha}C_n(t)}{dx^2}-\frac{\beta^2}{\alpha}C_n\right)$$.

Now when I try to solve this numerically mass is not conserved, but the mass non-conservation can also easily be seen by considering the system initially being of homogenous concentration $$C_0$$. During the first iteration, all the interior points remain fixed while at the boundary we get

$$C_n(t+dt)=C_0+dt\left(\alpha\frac{-\frac{2dx\beta}{\alpha}C_0}{dx^2}-\frac{\beta^2}{\alpha}C_0\right)$$

and

$$C_1(t+dt)=C_0+dt\left(\alpha\frac{\frac{2dx\beta}{\alpha}C_0}{dx^2}-\frac{\beta^2}{\alpha}C_0\right)$$.

Now this clearly does not conserve the mass of the system as it should per the no-flux boundary conditions. What is causing this?

Edit: By mass here I mean the integral of the concentration

$$\int_0^L C(x) dx$$

which is calculated using the trapezoidal rule.

Edit 2: Here is how the system currently evolves over time

The code can also be found here

# 1. Import libraries
import numpy as np
import matplotlib.pyplot as plt

# 2. Set up parameters
nx = 101
dx = 0.01
nt = 20000
alpha = 0.1
beta=-1
dt = 0.0001

# 3. Initial conditions

q = np.ones(nx)
qn= q.copy()

# 4. Solving the PDE
for n in range(nt):

qn=q.copy()

q[1] = qn[1] + alpha * dt / dx ** 2 * (qn[2]-2*qn[1]+qn[2]+2*dx*(beta/alpha)*qn[1]) +\
( beta*dt/(2*dx))*(qn[2] - (qn[2]+2*dx*(beta/alpha)*qn[1]))
for i in range(2, nx-2):

q[i] = qn[i] + alpha * dt / dx ** 2 * (qn[i+1]-2*qn[i]+qn[i-1]) +\
( beta*dt/(2*dx))*(qn[i+1] - qn[i-1])

q[nx-2] = qn[nx-2] + alpha * dt / dx ** 2 * (qn[nx-3]-2*qn[nx-2]+(qn[nx-3]-(beta/alpha)*2*dx*qn[nx-2])) +\
beta*dt/(2*dx)*(-qn[nx-3] +(qn[nx-3]-(beta/alpha)*2*dx*qn[nx-2]) )

# 5. Plot results
plt.plot(np.linspace(0, 1, 101), q)
plt.xlim(dx, 1-dx)
plt.show()

• Neither central finite difference nor forward Euler time stepping should be expected to exactly conserve mass. Can you please add more details as to how you are computing the mass at each step and what it looks like over time for your current implementation? Commented May 10 at 22:29
• @whpowell96 Thank you. I added some clarifications to my post as well as an plot of the "mass" which by here I mean the integral of the concentration. As you can see from the attached image the "mass" gets completely drained out of the system. Commented May 10 at 22:54
• This could be due to numerical diffusion from the central differences. Could you upload some plots of the solution at different times? Commented May 10 at 23:38
• @whpowell96 Sure, I uploaded some. Commented May 11 at 0:46

For $$\alpha, \beta = \mathrm{const}$$ you can rewrite your equation

$$\frac{\partial C}{\partial t}=\alpha\frac{\partial^2 C}{\partial x^2}+\beta \frac{\partial C}{\partial x}, \qquad (1),$$

to

$$\frac{\partial C}{\partial t}=\frac{\partial\left(\alpha \frac{\partial C}{\partial x} +\beta C\right)}{\partial x},\qquad (2).$$

With $$(2)$$ you can see that a conservative numerical method requires a consistent approximation of flux and boundary condition

$$f=-\left(\alpha\frac{\partial C}{\partial x}+\beta C\right),$$

in order to preserve

$$\frac{\partial C}{\partial t}+\frac{\partial f}{\partial x} = 0.$$

• Could you elaborate on what you mean by a "consistent approximation"? Commented May 11 at 12:26