# Not getting correct numerical solution for Advection-Diffusion-Reaction eqn

Objective: I am trying to numerically solve $C(x,y,t)$ from the following advection-diffusion-reaction equation in 2D space (x,y) and time. I will be testing my numerical solution with an approximate 1D analytical solution (given on 2nd page in this file). However, the issue is that I am not able to get correct numerical solution.

\begin{align} \text{ADR Equation: }\frac{\partial C}{\partial t} + \nabla\left(v C - D\nabla{C} \right)= \alpha C \end{align}

Method: I discretized the above ADR equation in 2D using finite-difference implicit scheme and I get the following equation in difference form.

\begin{align} p_1C^{n+1}_{i,j-1}+p_2C^{n+1}_{i-1,j}+p_3C^{n+1}_{i,j}+p_4C^{n+1}_{i+1,j}+p_5C^{n+1}_{i,j+1} = C^{n}_{i,j} \end{align}

I want to solve this as a system of equations using $A^{n+1} C^{n+1}=C^{n}$

Variables: $A$ is a matrix formed from variables $p_1, p_2, p_3, p_4, p_5$, which are constants in time but vary along $(x,y)$ grid. Superscript $n$ represents time index and subscripts $i,j$ represent space along $x,y$.

I.C: $$\textrm{C(x,y,0) = C_i, assume C_i as 1}$$

B.C: $$\frac{\partial C}{\partial x} = \frac{\partial C}{\partial y} = 0 \text{ at boundaries}$$. $$C(0,y=1,t) = \begin{cases} C_0\quad & \text{if 0<t\leq t_0} \\ 0 \quad & \text{if t>t_0} \end{cases}$$ $\textrm{assume$C_0$and$t_0$as 5000 and 2 days, respectively}$

Numerical Solution: I estimated $C(x,y,t)$ for $t=0$ to $t=100$ days, where dimensions of $x$ and $y$ are shown below in results.

$$\textrm{Fig-1: Snapshots of C(x,y) at various t}$$

In Fig-2, I am showing $C(x,t)$ for $0<t\leq t_0$ and $t=0$ to $t=100$ days:

$$\textrm{Fig-2: Snapshots of C(x,t) obtained after averaging C(x,y,t) over Y-axis}$$

Analytical Solution: Here I am showing $C(x,t)$ for $0<t\leq t_0$ and $t=0$ to $t=100$ days.

$$\textrm{Fig-3: Snapshots of C(x,t). I am trying to match this with Fig-2}$$

Issue: I have been trying to figure out if there is a problem in my numerical solution developed as per the implicit scheme discussed above; however, I have spent several hours and I couldn't figure it out. Please let me know if I missed providing anything.

• Have you checked if your code behaves correctly for simpler special cases first, e.g., constant coefficient heat equation? – Kirill Jul 19 '15 at 22:29
• I've checked and it gives mixed indications on its output. The concentrations fall off very quickly to zero. – user5510 Jul 19 '15 at 22:50
• You should compare it to known exact solutions. Checking whether a solution approaches zero is not conclusive. – Kirill Jul 19 '15 at 22:57
• @Pupil: Your first link has a permissions issue. – Geoff Oxberry Jul 20 '15 at 7:22
• The problem was the way I was forming the transmissibility matrix $A$ in $A\times X=B$ – user5510 Aug 10 '15 at 1:44