Implementing code for zero flux condition in Advection-Diffusion equation

Question

I would like to know how to implement a zero flux condition for the avdection-diffusion equation defined by:

Analysing the above we can realise that zero flux condition is satisfied when:

.

So, I write a code using finite difference scheme:

import numpy as np
import matplotlib.pyplot as plt

nx = 101
dx = 0.01
nt = 200
c = -1.
D = .1
dt = 0.0001

u = np.zeros(nx)
u[0:50] = 2
un = u.copy()

for n in range(nt):
    for i in range(1, nx-1):
    un = u.copy()
    u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])-\
           c*dt/dx*(un[i] - un[i-1])

    u[0] =  u[1]*(c*dx+D)/D #BC zero flux at left side

plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

Where

u[1]*(c*dx+D)/D

represent the zero flux condition which result from:

However, the result does not satisfy the zero flux condition, so the mass is not conservative along the time.

Anybody can help me to detect my mistakes?

Thanks in advance,

Answer 1

I have figured out the answer to my question, and I will share it here in case someone find it helpful.

I wanted to establish zero-flux conservative boundary conditions to the advection difussion equation, this represents the Robin Boundary Conditions. Given it is almost impossible to write equation in SO, you can review an excellent explanation here: https://scicomp.stackexchange.com/questions/5434/conservation-of-a-physical-quantity-when-using-neumann-boundary-conditions-appli

The below code include the RBC to the advection diffusion equation, which solve my problem.

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

# 2. Set up parameters
nx = 101
dx = 0.01
nt = 7000
c = .5
D = .1
dt = 0.0001

# 3. Initial conditions
u = np.zeros(nx)
u[10:35] = 4
un = u.copy()

# 4. Solving the PDE
for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] + D * dt / dx ** 2 * (un[i+1]-2*un[i]+un[i-1]) -\
           c*dt/dx*(un[i] - un[i-1])
    u[0] = u[1] - dx * c / D * u[1]    # Robin Boundary Condition LEFT
    u[-1] = u[-2] + dx * c / D * u[-1] # Robin Boundary Condition RIGHT

# 5. Plot results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()