1-D Diffusion-Reaction Simulation
By Taylor Geisler
The Problem
In this simulation we will investigate the diffusion of two species A and B across a length L. One side of the mixture is maintained as pure A. The other side is pure B. If we leave these two species to diffuse and react until they reach steady state, what will be the final concentration of each species be across the length of the diffusion system?
First, we write the mass balance equation describing each species. The first term on the right hand side is a measure of the change of concentration due to diffusion. This term is derived from Fick’s first law. The second term \( S_{i} \) is the change in concentration due to reaction.
\(
\frac{\partial C_{i}}{\partial t}=-\frac{\partial}{\partial z}\left(-D_{i}\frac{\partial C_{i}}{\partial z}\right)+S_{i}
\)
At steady state \(\frac{\partial C_{i}}{\partial t}=0 \). Assume that the diffusion coefficient \( D_{i} \) is constant throughout the mixture. The last equation then becomes
\(
0=D_{i}\frac{\partial^{2}C_{i}}{\partial z^{2}}+S_{i}
\)
Using a finite difference approximation of the second derivative at location j gives
\(
0=D_{i}\left(\frac{C_{i\, j-1}-2C_{i\, j}+C_{i\, j+1}}{\Delta z^{2}}\right)+S_{i\, j}
\)
We will investigate the reaction \(A+B\overset{k}{\rightarrow}C \) with \( r_{A}=r_{B}=-r_{C}=-kC_{A}C_{B}\). Plugging the rate of reaction into the equation gives
\(
D_{i}\left(\frac{C_{i\, j-1}-2C_{i\, j}+C_{i\, j+1}}{\Delta z^{2}}\right)=kC_{A}C_{B}
\)
The solution will be easier if this nonlinear equation is linearized.
\(
kC_{A}C_{B}\approx k(C_{A}C_{B}^{*}+C_{A}^{*}C_{B}-C_{A}^{*}C_{B}^{*})
\)
Then for species A and B
\(
D_{A}\left(\frac{C_{A\, j-1}-2C_{A\, j}+C_{A\, j+1}}{\Delta z^{2}}\right)=k(C_{Aj}C_{B}^{*}+C_{A}^{*}C_{Bj}-C_{A}^{*}C_{B}^{*})
\)
\(
D_{B}\left(\frac{C_{B\, j-1}-2C_{B\, j}+C_{B\, j+1}}{\Delta z^{2}}\right)=k(C_{Aj}C_{B}^{*}+C_{A}^{*}C_{Bj}-C_{A}^{*}C_{B}^{*})
\)
We now have an equation for both species at each location j across the length of the system. Rearranging the equations into a form with all of the like terms grouped together gives
\(
\frac{D_{A}}{\Delta z^{2}}C_{Aj-1}+\left(\frac{-D_{A}2}{\Delta z^{2}}-kC_{B}^{*}\right)C_{Aj}+\frac{D_{A}}{\Delta z^{2}}C_{Aj+1}-kC_{A}^{*}C_{Bj}=-kC_{A}^{*}C_{B}^{*}
\)
\(
\frac{D_{B}}{\Delta z^{2}}C_{Bj-1}+\left(\frac{-D_{B}2}{\Delta z^{2}}-kC_{A}^{*}\right)C_{Bj}+\frac{D_{B}}{\Delta z^{2}}C_{Bj+1}-kC_{B}^{*}C_{Aj}=-kC_{A}^{*}C_{B}^{*}
\)
These two equations, evaluated at each location j across the system, can be represented by a system of linear equations. If we decide to sample at N locations across our system, we will have 2N linear equations to solve. Before a solution can be reached an educated guess of \(C^{*}\) must be made for each species at each location j. Once \(C^{*}\) is guessed, the linear system can be solved, giving values of \(C_{A}\) and \(C_{B}\) at each location j.
If we didn’t guess the exact solution, the solution will be different than what we originally guessed. We will now update the guess with the new solution values and solve the system of equations again. With each iteration the guesses converge on the actual solution. When our solution gives us the same concentration that we guessed, within a certain tolerance, we know that we have converged on the steady state solution.
The concentration of C was determined by subtracting the remaining concentrations of A and B that had not reacted at each location from the total concentration.
The Simulation
The solutions to these equations are plotted above. Experimenting with the constants in these equations gives interesting results. The diffusion limited case occurs when the reaction constant k is very low or the length is very small. The reaction limiting case is also very interesting. It occurs when k is very high or when the length of the system is large. Experiment with different values to better understand the effects that rates of reaction, diffusion coefficients, concentration, and scale have on diffusion.