This notebook contains course material from CBE20255 by Jeffrey Kantor (jeff at nd.edu); the content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.
This Jupyter notebook demonstrates the application of the general material balance equations to a variety of example problems.
The general principle of material balances is stated in words as
Accumulation = Inflow - Outflow + Generation - Consumption
An equation like this can be written for any conserved quantity, whether it is chemical species that chemical or bio engineers work with, or moeny in the case of finance, or populations for social scientists. Whenever terms on the right-hand side of the equation do not yield a zero result we are left with an unsteady-state balance, which are some of the most fascinating and challenging problems in engineering practice.
The first example concerns population growth. Looking birth and death rates in various countries around the globe, one comes across the following data maintained by the World Bank.
For Afghanistan (the first country on the World Bank charts), the 2011 population was estimated to be 29,105,480 with a birth rate of 37.0 births per year per 1,000 people, and a death rate as 8.0 per year per 1,000. Assuming there is no in-migration or out-migration, can we predict the population in 2014?
Writing out the balance equation, we get
$\underbrace{\mbox{Accumulation}}_{\frac{dP}{dt}} = \underbrace{\mbox{Inflow}}_{=0} - \underbrace{\mbox{Outflow}}_{=0} + \underbrace{\mbox{Generation}}_{births=\frac{37}{1000}P} - \underbrace{\mbox{Consumption}}_{deaths=\frac{8}{1000}P}$
the algebraic combination of the two terms on the left leaves us an equation for the rate of change of population
$\frac{dP}{dt} = 0.029\,P$
with an initial condition $P(2011) = 29,105,480$.
We'll use the shorthand notation $t_0 = 2011$ and $P(t_0) = 29,105,480$. Then separating the variables
$\int_{P(t_0)}^{P(t)}\frac{1}{P}\,dP = 0.029\,\int_{t_0}^t dt $
Doing the integrals (if this is the first time you've seen this, be sure you do this by hand!) leads to
$P(t) = P(t_0)\,e^{0.029(t-t_0)}$
This solution demonstrates the exponential growth of populations in circumstances where growth rate is constant and there are no other factors to consider.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
P_0 = 29105480
t_0 = 2011
t_1 = 2014
t = np.linspace(t_0,t_1)
P = P_0*np.exp(0.029*(t-t_0))
plt.plot(t,P)
plt.xlabel('Year')
plt.ylabel('Count')
plt.title("Proj. Population in Afganistan {:4d} to {:4d}".format(t_0,t_1))
<matplotlib.text.Text at 0x11172b400>
Another way of generating solutons to unsteady-state problems is by direct, numerical solution of the differential equations. This requires us to set up a function to numerically evaluate the right hand-side of the differential equation. That function, along with initial conditions and a list of time values for which we want to the know the solution, are passed to a solver that does the hard work.
def Inflow(t) : return 0
def Outflow(t): return 0
def Births(P) : return (37.0/1000)*P
def Deaths(P) : return (8.0/1000)*P
def Accumulation(P,t) : return Inflow(t) - Outflow(t) + Births(P) - Deaths(P)
P_0 = 29105480
t_0 = 2011
t_1 = 2014
t = np.linspace(t_0,t_1)
The SciPy
package provides a number of tools for the numerical solution of differential equations. Of these, perhaps the easiest to use is odeint
which is demonstrated here. odeint
needs to be imported from scipy.integrate
, and then given three pieces of information:
P = odeint(Accumulation,P_0,t)
plt.plot(t,P)
plt.xlabel('Year')
plt.ylabel('Population')
plt.title('Afghanistan')
<matplotlib.text.Text at 0x1118c46a0>
We have a tank partially filled with an initial volume V(t0) = 1000 liters of water that is contaminated by a toxic species X at an concentration CX(t0) = 100 ppm by volume. We wish to dilute this to 20 ppm by volume before sending it for further cleanup. We know adding 4000 liters of water would accomplish the goal, but for process monitoring purposes we'd like to compute concentration as a function of time assuming water enters at a steady rate of q(t) = 1 liter/second.
The total amount of X in the tank is the product of concentration times volume, that is CXV. The material balance for X is
For component X: $\underbrace{\mbox{Accumulation}_X}_{\frac{d(C_XV)}{dt}} = \underbrace{\mbox{Inflow}_X}_{=0} - \underbrace{\mbox{Outflow}_X}_{=0} + \underbrace{\mbox{Generation}_X}_{=0} - \underbrace{\mbox{Consumption}_X}_{=0}$
This is a two component system, so we need two material balances. We'll do the second balance on overall volume on the assumption that density is a constant this dilute solution. The overall balance is
For total volume: $\underbrace{\mbox{Accumulation}}_{\frac{dV}{dt}} = \underbrace{\mbox{Inflow}}_{=q} - \underbrace{\mbox{Outflow}}_{=0} + \underbrace{\mbox{Generation}}_{=0} - \underbrace{\mbox{Consumption}}_{=0}$
This gives us a pair of differential equations
$\begin{align*} \frac{d(C_XV)}{dt} & = 0 \\ \frac{dV}{dt} & = q \end{align*}$
What we're looking for is CX. But what we have are differential equations for product (CXV) and volume V. How can we get a differential equation for CX?
Using the chain rule
$\frac{d(C_XV)}{dt} = C_X\,\underbrace{\frac{dV}{dt}}_{=q} + V\,\frac{dC_X}{dt} = 0$
Substituting in the second equation leaves us with a pair of differential equations
$\begin{align*} \frac{dC_X}{dt} & = -\frac{q}{V} C_X \\ \frac{dV}{dt} & = q \end{align*}$
(Be sure you go through the algebra and understand both how and why we did this.)
# parameter values
q = 1 # liters/second
# initial conditions
C_0 = 100 # ppm
V_0 = 1000 # liters
# time grid in seconds
t = np.linspace(0,4000)
# compute list of derivative values
def deriv(X,t):
C,V = X
return [-q*C/V,q]
# solve with odeint
soln = odeint(deriv,[C_0,V_0],t)
# display the result
plt.plot(t,soln[:,0])
plt.xlabel('Time [seconds]')
plt.ylabel('[ppm/v]')
plt.title('Concentration of X')
<matplotlib.text.Text at 0x111a15630>
$\begin{align*} \frac{dH}{dt} & = r H \left(1-\frac{H}{k}\right) - \frac{a H L}{c + H}\\ \frac{dL}{dt} & = \frac{b a H L}{c + H} - d L \end{align*}$
a = 3.2
b = 0.6
c = 50
d = 0.56
k = 125
r = 1.6
def deriv(x,t):
H,L = x
dH = r*H*(1-H/k) - a*H*L/(c+H)
dL = b*a*H*L/(c+H) - d*L
return [dH,dL]
t = np.linspace(0,150,500)
ic = [20,30]
soln = odeint(deriv,ic,t)
plt.plot(t,soln)
plt.xlabel('Time [yr]')
plt.ylabel('Population');
plt.title('Lynx-Hare Population Model')
<matplotlib.text.Text at 0x111b2ecc0>
H = soln[:,0]
L = soln[:,1]
plt.plot(H,L)
plt.xlabel('Hare Population')
plt.ylabel('Lynx Population')
<matplotlib.text.Text at 0x111ca3978>