Solving a one-component 1D diffusion-reaction equation - Multiphase Reactor Modeling using PyMRM

In this exercise you will build an implicit solver for a diffusion equation. A one-component stationary diffusion equation with constant coefficients has the form

D\ \frac{\partial^{2}c}{\partial x^{2}} = 0

(1)

This equation needs to be solved on a domain of length $L$ with on the left and right boundaries a boundary condition of the form

a\frac{\partial c}{\partial n} + bc = d,

(2)

where $n$ indicates the outward pointing normal direction, i.e. $- x$ for the left boundary and $x$ for the right boundary.

Questions:

Assume a spatial discretization of $N$ cells with points $j$ located in cell centers. Write down the spatial discretization formula for the diffusion terms for points $j = 1,\ldots,N - 2$ . (Note that cells $j = 0$ and $j = N - 1$ , are neighboring a boundary).
Write down the spatial discretization formula’s for points $j = 0$ and $j = N - 1$ by implementing a Dirichlet boundary condition ( $a = 0$ , $b = 1$ ).
Write down the resulting formula’s in matrix-vector form.
Implement the matrix-vector equation in Python using SciPy. Define a sparse matrix using the SciPy sparse array, e.g., by means of: scipy.sparse.diags_array.
Verify the implementation for Dirichlet BCs on both sides.
Extend the boundary condition implementation to general mixed boundary conditions (any value of $a$ and $b$ ). Write down the resulting formula’s in matrix-vector form.
Verify the implementation of the mixed boundary conditions by choosing values of $a$ , $b$ and $d$ on both sides and compare the result with the expected analytical solution.

Dispersion with a first order reaction obeys

D\ \frac{\partial^{2}c}{\partial x^{2}} - kc = 0

(3)

Questions:

Include the first order reaction term to the matrix-vector equation.
Verify the implementation for varying values of $k$ and general mixed boundary conditions by comparison with the expected analytical solution.

The unsteady diffusion (first order) reaction equation is

\frac{\partial c}{\partial t} = D\ \frac{\partial^{2}c}{\partial x^{2}} - kc

(4)

Questions:

Add an implementation of the accumulation term using Euler-backward time discretization.
Solve the equations in time and verify the accumulation term by comparison with analytical solutions for (simple) initial and boundary conditions.