I am trying to solve a problem with 2D Convection-Diffusion equation with U = Concentration ($mg/m^{2}$) using Implicit Upwind Finite Difference Method like this

$$ \frac{\partial U}{\partial t} + v_{y}\frac{\partial U}{\partial y} = D(\frac{\partial^{2} U}{\partial x^{2}}+ \frac{\partial^{2} U}{\partial y^{2}}) $$

Set $\sigma_D = D\frac{\Delta t}{\Delta x^{2}} = D\frac{\Delta t}{\Delta y^{2}}$ and $\sigma_{C} = \frac{v_{y}\Delta t}{\Delta y}$

I have dicretized the equation above using Implicit Upwind scheme as I have vy = -0.01 m/s going down towards the ground:

$$ U_{i,j}^{m+1}(1+4\sigma_D+\sigma_{C})-\sigma_{D}(U_{i+1,j}^{m+1}+U_{i-1,j}^{m+1} + U_{i,j+1}^{m+1})-(\sigma_{C}+\sigma_{D})U_{i,j-1}^{m+1} = U_{i,j}^{m} $$ with domain [x,y] = [10km,4km]. Dirichlet BCs are at the top, left and right of the domain and Robin BC is at the bottom (modelled as ground):

$$ D\frac{\partial U}{\partial y} - v_{y}U = 0 $$

Set $\sigma_{R} = \frac{v_{y}\Delta y}{D}$

Dicretized Robin BC using downwind method (since the wind is going down):

$$ v_{y}U_{i,1}^{m+1}-D(\frac{U_{i,1}^{m+1} - U_{i,0}^{m+1}}{\Delta y}) $$

Rearrange the equation above gives:

$$ U_{i,1}^{m+1}(\sigma_{R}-1)+U_{i,0}^{m+1}=0 $$

Initial Conditon: 1000kg of chemical is dropped in the middle of domain [5km,2km] ( U = Concentration [mg/m2] ):

$$ U_{init} = \frac{1000E6}{\Delta x \Delta y} $$

This is my sample code to run with this method:

Lx = 10E3
Ly = 4E3
dx, dy = 20, 20
nx = int(Lx/dx + 1)
ny = int(Ly/dy + 1)
D = 0.5 
x = np.linspace(0.0, Lx, nx)
y = np.linspace(0.0, Ly, ny)
mid_x = int(nx/2)
mid_y = int(ny/2)
mass = 1000E6 #mg
vy = -0.01
sigma_D = 1
dt = sigma_D * min(dx, dy)**2 / D
sigma_c = vy*dt/dy
day = 4
time = 3600*24*day
nt = int(time/dt)
def IC(nx, ny, mass, dx, dy, mid_x, mid_y):
    Dis = np.zeros((ny, nx))
    Con = mass / (dy*dx)
    Dis[mid_y, mid_x] = Con
    return Dis, Con
U0, C = IC(nx, ny, mass, dx, dy, mid_x, mid_y)

def Implicit(U0, nt, dx, dy, D, vx, vy, frn, dt, nx, ny):
    sigma_D = D*dt/(dx)**2
    sigma_c = vy*dt/dy
    sigma_R = vy*dy/D
    AA = csr_matrix((nx*ny, nx*ny)).tolil(copy = True)
    #Boundary Condition
    for j in range(ny):
        for i in range(nx):
            AA[i + j*nx, i + j*nx] = 1
    #Inner nodes 
    for j in range(1, ny - 1):
        for i in range(1, nx - 1):
            #n + 1 side
            AA[i + j*nx, i + j*nx]     = 1 + 4*sigma_D + sigma_c
            AA[i + j*nx, i+1 + j*nx]   = -sigma_D
            AA[i + j*nx, i-1 + j*nx]   = -sigma_D
            AA[i + j*nx, i + (j+1)*nx] = -sigma_D 
            AA[i + j*nx, i + (j-1)*nx] =  -(sigma_D +sigma_c)
    #Robin BC
    for j in range(0, 1):
        for i in range(nx):
            AA[i + j*nx, i + j*nx] = -1
            AA[i + j*nx, i + (j+1)*nx] = -sigma_R + 1
    Matrix_1 = AA.tocsr()
    for n in tqdm(range(nt)):
        U0[0] = 0
        U0 = csr_matrix(U0.reshape(ny*nx, 1))
        U1 = scipy.sparse.linalg.bicgstab(Matrix_1, U0.todense())[0]

        U0 = U1.copy().reshape((ny, nx)) 
    return U0

This is what I got from the code after 4 days (3600*24*4 seconds):
Concentration Plot in 3D Concentration Plot

Total concentration over time
Total Concentration

I am not sure why my total concentration keeps losing over time. Theoretically, it should stay the same over time since when it hits the ground, it is stuck there.

I believe my boundary condition for Robin BC might be wrong but I am not sure where I went wrong. Any suggestions?

  • $\begingroup$ I believe both your boundary conditions are not conservative. Neumann boundary conditions would be. Robin boundary condition models a heat exchange with the surroundings $\endgroup$
    – VorKir
    Oct 21, 2019 at 0:11
  • $\begingroup$ You can compute the concentration that leaves the domain, by approximating the fluxes $U \vec{v}\cdot \vec{n} + \vec{\nabla U}\cdot \vec{n}$, with $\vec{n}$ the unit vector normal to the boundary, pointing to the outside. Integrate these fluxes along the boundaries to find how much U you are loosing per unit time. For example, the flux of concentration that leaves through the upper boundary (// to $\vec{y}$) is: $\phi = \sum_{j=0}^{ny} \Delta y \left( \frac{U_{dirichlet} - U_{ij}}{\Delta x / 2} \right)$. You must force convective and diffusive boundary fluxes to 0, to ensure conservation. $\endgroup$
    – Laurent90
    Nov 21, 2020 at 10:21

1 Answer 1


I think you have some confusion about convection-diffusion equation (CDE) as well as assuming a non-physical velocity profile here. Let's rewrite your CDE as:

$$\frac{\partial \phi}{\partial t} + \nabla \cdot (-D \nabla \phi + \mathbf{u} \phi) = 0$$

Where $\phi$ is concentration, $D$ is diffusion coefficient, and $\mathbf{u}$ is velocity profile of your fluid. You said you have three Dirichlet boundary conditions at the top, left, and right, so they are formulated as:

top: $\phi(x,y_{top}) = \phi_{top}$

left: $\phi(x_{left},y) = \phi_{left}$

right: $\phi(x_{right},y) = \phi_{right}$

and finally, as far as I understand from your description, you want to have zero flux ($\mathbf{J} = -D \nabla \phi + \mathbf{u} \phi$) at the bottom. But you need to know that putting a constant velocity in $y$ direction, is somewhat unphysical. Why? Cause you know that at the walls (left and right), your velocity must be zero, but constant velocity does not satisfy these conditions. I suggest you replace constant velocity in $y$ direction with this formula:

$$\mathbf{u} = u_{max} (x-x_{left}) (x - x_{right}) \mathbf{j}$$

You want to make sure this new velocity profile is incompressible:

$$\nabla \cdot \mathbf{u} = \frac{\partial u_{x}}{\partial x} + \frac{\partial u_{y}}{\partial y} = 0$$

Finally, you want to have outflow boundary condition at the bottom. The most popular assumption at the outflow is: diffusive flux ($-D \nabla \phi$) is negligible at the outlet in comparison to convective flux ($\phi \mathbf{u}$), So:

bottom: $-D \frac{\partial \phi}{\partial y}|_{y = y_{bottom}} = 0$

But, you need to know that the thing which will remain constant is the total mass in the region, so:

$$\frac{d\Phi}{d t} = \frac{d}{dt} \int_{\Omega} \phi d^{2} \mathbf{r} = 0$$

So, in order to check the correctiveness of your implementation, you could find $\Phi$ the total mass in the computational domain ($\Omega$) and you should not see any drift or increase in total mass.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.