Problem setup¶
One way of treating waste water is to run it through an axially dispersed plug-flow reactor. As the contaminated water flows through the reactor, the solution reacts with catalysts that line the sides of the reactor, reducing the concentration of the contaminant along the length of the tube.
The following differential equation describes the steady-state concentration of the contaminant that reacts with first-order kinetics in an axially dispersed plug-flow reactor:
$$ D \dfrac{d^2c}{dx^2} - U \dfrac{dc}{dx} - kc = 0,$$where $D$ is the dispersion coefficient [m$^2$/hr], $c$ is the concentration of the contaminant [mmol/L], $x$ is the distance along the tube [m], $U$ is the velocity of the flow [m/hr], and $k$ is the reaction rate [hr$^{-1}$].
The boundary conditions for this reactor are called Danckwerts boundary conditions, and are formulated as:
$$Uc_{in} = \left.Uc\right|_{x=0} - D\left.\dfrac{dc}{dx}\right|_{x = 0};$$$$\left.\dfrac{dc}{dx}\right|_{x = L} = 0,$$where $c_{in}$ is the concentration of the contaminant at the inflow [mmol/L] and $L$ is the length of the reactor [m].
The analytical solution of this ordinary differential equation with the Danckwerts boundary conditions is
$$ c(x) = \dfrac{U c_{in}\left(\lambda_2 e^{\lambda_2 L} e^{\lambda_1 x} - \lambda_1 e^{\lambda_1 L} e^{\lambda_2 x}\right)}{(U - D\lambda_1)\lambda_2 e^{\lambda_2 L} - (U - D\lambda_2)\lambda_1 e^{\lambda_1 L}},$$where
$$\lambda_1 = \dfrac{U}{2D}\left(1 + \sqrt{1 + \dfrac{4kD}{U^2}}\right), \quad \text{ and } \quad \lambda_2 = \dfrac{U}{2D}\left(1 - \sqrt{1 + \dfrac{4kD}{U^2}}\right).$$The governing differential equation and boundary conditions can be reformulated into a system of algebraic equations by replacing the derivatives with centered finite-difference approximations. Namely,
$$ \dfrac{d^2c}{dx^2} \approx \dfrac{c_{i-1} - 2c_i + c_{i+1}}{\Delta x^2};$$$$ \dfrac{dc}{dx} \approx \dfrac{c_{i+1} - c_{i-1}}{2\Delta x},$$where $c_{i}$ is the concentration at some point $i$ in the tube (corresponding to the location $x_i$) and $\Delta x = x_{i+1} - x_i$ is the step size between those points.
Suppose you are a waste water engineer trying to determine the effectiveness of a plug-flow reactor your city is planning on purchasing. The water your plant processes contains dangerous levels of nitrates that violate EPA requirements and pose a threat to your city's citizens. The new reactor must reduce the level of nitrates down to below the EPA threshold of $0.16128$ [mmol/L] before the water can be distributed to homes.
You want to model the concentration of nitrates as a function of the fluid's position in the plug-flow reactor ($x$) to ensure that the new equipment meets your city's needs. The parameters you have are: $D = 500$ [m$^2$/hr], $U = 100$ [m/hr], $k = 2$ [hr$^{-1}$], $L = 70$ [m], and $c_{in} = 0.64512$ [mmol/L]. To complete the model:
Rewrite, using LaTeX, the governing equation by substituting in the centered finite-difference formulas above. Collect terms so that your equation has the form $$$$$$ e_i c_{i-1} + f_i c_i + g_i c_{i+1} = 0, \quad i = 1, 2, \dots, n-2.$$$$$$ Here, the domain is divided into $n$ nodes, and $e_1 = e_2 = \dots = e_{n-2}$, $f_1 = f_2 = \dots = f_{n-2}$, and $g_1 = g_2 = \dots = g_{n-2}$ are coefficients of the interior nodes in terms of $D$, $U$, $k$, and $\Delta x$ (do not substitute yet). Report the expressions you derived for $e_i$, $f_i$, and $g_i$ in $i = 1, 2, \dots, n-2$.
Rewrite the Danckwerts boundary conditions by substituting in the centered finite-difference formulas. Use these new expressions along with the expressions from Problem 1 and collect terms to derive the equations:$$$$$$f_0 c_0 + g_0 c_1 = r_0;$$ $$$$ $$e_{n-1} c_{n-2} + f_{n-1} c_{n-1} = r_{n-1}.$$ $$$$ Again, $f_0$, $g_0$, $r_0$, $e_{n-1}$, $f_{n-1}$, and $r_{n-1}$ should be expressed in the given parameters ($D$, $U$, $k$, $c_{in}$, and $\Delta x$). Report the six coefficient expressions you derived.
Create a function in Python that uses the Thomas algorithm to solve a tridiagonal system of equations (such as this problem). Your function should take in three arrays that represent the values on the diagonals (i.e., $e_i$, $f_i$, and $g_i$ for $i = 0, 1, \dots n-1$) and an array that represents the values on the right-hand-side of the equations (i.e., $r_i$ for $i = 0, 1, \dots, n-1$). Your function should return an array that represents the concentration values at each node. Indicate where this function is coded on your report.
Using the expressions you derived above and the Thomas algorithm function, determine the concentration of the contaminant at each point along the length of the tube, where $x_0 = 0$ and $x_{n-1} = L$. Using a coarse mesh, divide the tube into five nodes ($n = 5$). Next, repeat the same problem, but with a fine mesh by dividing the tube into one hundred nodes ($n = 100$). Plot the concentration results from the coarse mesh and fine mesh (using scatter plots or line plots with markers) along with the analytical solution.
How similar are your three solutions (the analytical, the coarse mesh, and the fine mesh)? If your solutions are not similar, what would explain the differences and why? Briefly describe/quantify your error and explain any sources.
Will the plug-flow reactor you tested release water that meets the EPA nitrate requirements? Briefly explain why or why not.
Use Python to generate all plots and to perform all calculations. Report your expressions from problems 1 and 2 using LaTeX in your report.