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 formulation and solution of material balances for an ethylene oxide flowsheet in Python using the symbolic algebra package Sympy.
While material balances for chemical processes are readily expressed as linear equations, extracting the matrix of coefficients for the linear equations can be tedious and error-prone. Fortunately, symbolic calculations can be used to solve material balance problems typical of introductory chemical engineering courses. This is demonstrated below using an example from the textbook * Introduction to Chemical Processes: Principles, Analysis, and Synthesis* by Regina Murphy.
"TASNEE 001" by Secl - Own work. Licensed under CC BY 3.0 via Wikimedia Commons.
The problem is to analyze the capability of an existing process for the production of ethylene oxide by the reaction of ethylene and oxygen
$$2\,C_2H_4 + O_2 \longrightarrow 2\,C_2H_4O$$The target production is 1.7 × 106 kgmol/year of an ethylene oxide product with 98 mol% purity. The feedrate of ethylene is specified as 196 kgmol/hr, and of oxygen is specified as 84.5 kgmol/hr. The reactor has a nominal single pass conversion of 6% ethylene. The separator recovers 98% of the ethylene and 99.5% of the oxygen and 98% of the ethylene for recycle, and 97% of the ethylene oxide for the product stream.
The problem is to find the nominal product purity and production rates, and to examine the sensitivity of product purity and production to the equipment specifications.
Sympy is a library written in pure python for symbolic algebra. The solution strategy is to represent the stream variables and reaction extents as symbolic variables, then express the material balances and process specifications as symbolic equations, and finally to use the Sympy solver to find a nominal solution and to perform parametric analyses. The following cell initializes sympy
and necessary printing and display functions.
%matplotlib inline
import sympy as sym
sym.init_printing()
from IPython.display import display
The stream variables are systematically created using the Sympy var
function, and added to the list streams
.
stream_table = []
stream_table.append(sym.var('nE1 nO1 nEO1')) # Stream 1
stream_table.append(sym.var('nE2 nO2 nEO2')) # Stream 2
stream_table.append(sym.var('nE3 nO3 nEO3')) # Stream 3
stream_table.append(sym.var('nE4 nO4 nEO4')) # Stream 4
stream_table.append(sym.var('nE5 nO5 nEO5')) # Stream 5
display(stream_table)
This problem includes only a single reaction in a single reactor. A corresponding extent of reaction variable is created, and a master list of all process variables is created.
extents = [sym.var('X')]
display(extents)
Create list of all variables
variables = []
for x in extents:
variables.append(x)
for s in stream_table:
for v in s:
variables.append(v)
display(variables)
Material balances are written for each process unit using the symbolic stream and extent of reaction variables created above. Each material balance is expressed as the net rate of accumulation that will be set to zero to specify steady-state operation. The material balances are gathered into a list for each process unit, then the lists are concantenated to create a list of all material balances.
mixer = [
sym.Eq(0, nE1 + nE5 - nE2),
sym.Eq(0, nO1 + nO5 - nO2),
sym.Eq(0, nEO5 - nEO2)]
reactor = [
sym.Eq(0, nE2 - nE3 - 2*X),
sym.Eq(0, nO2 - nO3 - X),
sym.Eq(0, nEO2 - nEO3 + 2*X)]
separator = [
sym.Eq(0, nE3 - nE4 - nE5),
sym.Eq(0, nO3 - nO4 - nO5),
sym.Eq(0, nEO3 - nEO4 - nEO5)]
material_balances = mixer + reactor + separator
for eqn in material_balances:
display(eqn)
Process specifications are written as equalities using the Sympy Eq
function.
feed_spec = [
sym.Eq(nE1, 196.0),
sym.Eq(nO1, 84.5),
sym.Eq(nEO1, 0.0)]
reactor_spec = [
sym.Eq(nE2 - nE3, 0.06*nE2)]
separator_spec = [
sym.Eq(nE5, 0.98*nE3),
sym.Eq(nO5, 0.995*nO3),
sym.Eq(nEO4, 0.97*nEO3)]
specifications = feed_spec + reactor_spec + separator_spec
for eqn in specifications:
display(eqn)
A simple degree of freedom analysis is to compare the number of variables to the number of equations.
nVars = 0
for s in stream_table:
for v in s:
nVars += 1
print("Stream: {0:2d} Variable: {1:5s}".format(nVars,v.name))
print("\n%d Extents of Reaction\n" % len(extents))
for v in extents:
print("Extent: ", v.name)
print("\n%d Variables = %d Stream Variables + %d Extents of Reaction \n" \
% (len(variables),nVars,len(extents)))
Stream: 1 Variable: nE1 Stream: 2 Variable: nO1 Stream: 3 Variable: nEO1 Stream: 4 Variable: nE2 Stream: 5 Variable: nO2 Stream: 6 Variable: nEO2 Stream: 7 Variable: nE3 Stream: 8 Variable: nO3 Stream: 9 Variable: nEO3 Stream: 10 Variable: nE4 Stream: 11 Variable: nO4 Stream: 12 Variable: nEO4 Stream: 13 Variable: nE5 Stream: 14 Variable: nO5 Stream: 15 Variable: nEO5 1 Extents of Reaction Extent: X 16 Variables = 15 Stream Variables + 1 Extents of Reaction
equations = material_balances + specifications
print("\n%d Equations = %d Material Balances + %d Specifications" \
% (len(equations),len(material_balances),len(specifications)))
print("\n%d Material Balances\n" % len(material_balances))
for mb in material_balances:
print(mb)
print("\n%d Specifications\n" % len(specifications))
for spec in specifications:
print(spec)
16 Equations = 9 Material Balances + 7 Specifications 9 Material Balances Eq(0, nE1 - nE2 + nE5) Eq(0, nO1 - nO2 + nO5) Eq(0, -nEO2 + nEO5) Eq(0, -2*X + nE2 - nE3) Eq(0, -X + nO2 - nO3) Eq(0, 2*X + nEO2 - nEO3) Eq(0, nE3 - nE4 - nE5) Eq(0, nO3 - nO4 - nO5) Eq(0, nEO3 - nEO4 - nEO5) 7 Specifications Eq(nE1, 196.000000000000) Eq(nO1, 84.5000000000000) Eq(nEO1, 0.0) Eq(nE2 - nE3, 0.06*nE2) Eq(nE5, 0.98*nE3) Eq(nO5, 0.995*nO3) Eq(nEO4, 0.97*nEO3)
soln = sym.solve(material_balances + specifications)
for k in soln.keys():
print("Variable {0:4s}: {1:8.2f}".format(str(k),round(soln[k],2)))
Variable nO3 : 1976.14 Variable nEO3: 153.85 Variable nEO1: 0.00 Variable nE3 : 2338.07 Variable X : 74.62 Variable nO4 : 9.88 Variable nEO5: 4.62 Variable nO2 : 2050.76 Variable nE1 : 196.00 Variable nEO2: 4.62 Variable nEO4: 149.24 Variable nO5 : 1966.26 Variable nE2 : 2487.31 Variable nO1 : 84.50 Variable nE4 : 46.76 Variable nE5 : 2291.31
# display solution for each variable, rounded to 1 decimal place
for v in variables:
display(sym.Eq(v,round(soln[v],1)))
purity = soln[nEO4]/(soln[nEO4]+soln[nE4]+soln[nO4])
production = 24*350*(soln[nEO4] + soln[nE4] + soln[nO4])/1000000
print("Annual Production %4.2f million kgmol/year at %5.3f purity." \
% (production,purity))
Annual Production 1.73 million kgmol/year at 0.725 purity.
The problem asked for an analysis of the sensitivity of the problem results to changes in unit performance. This is implemented by restating the specifications where a key parameter is replaced by a symbolic variable, and the process the performance plotted as a function of the parameter.
feed_spec = [
sym.Eq(nE1, 196.0),
sym.Eq(nO1, 84.5),
sym.Eq(nEO1, 0.0)]
fconv = sym.var('fconv')
reactor_spec = [
sym.Eq(nE2 - nE3, fconv*nE2)]
separator_spec = [
sym.Eq(nE5, 0.98*nE3),
sym.Eq(nO5, 0.995*nO3),
sym.Eq(nEO4, 0.97*nEO3)]
specifications = feed_spec + reactor_spec + separator_spec
for s in specifications:
display(s)
Recycle calculations are introduce a strong dependence of flow rates on parameters such as fraction conversion in the reactor and fractional recovery in separation units. To see this, here we solve for the flowrate of E2 as a function of fractional conversion of ethylene in the reactor.
From the material balances
E2 = E1 + E5
E3 = (1-fconv) E2
E5 = 0.98 E3
Take a moment and solve these by hand.
soln = sym.solve(material_balances + specifications, exclude=[fconv])
soln
display(sym.Eq(nE2,soln[nE2]))
sym.plot(soln[nE2],(fconv,0,1),xlabel='fconv',ylabel='kg-mol/hour',title='E2 flowrate');
for v in variables:
display(sym.Eq(v, soln[v]))
purity = soln[nEO4]/(soln[nEO4]+soln[nE4]+soln[nO4])
display(sym.simplify(purity))
sym.plot(purity,(fconv,0,.2),xlabel='fconv',ylabel='fraction EO',title='Purity');
production = 24*350*(soln[nEO4] + soln[nE4] + soln[nO4])/1000000
display(sym.simplify(production))
sym.plot(production,(fconv,0,0.2),xlabel='fconv',ylabel='kg-mol/year',title='Production');
Will improving the recovery of Ethylene Oxide to the product stream improve the product purity? Make a prediction, then check against the results of the following calculations.
feed_spec = [
sym.Eq(nE1,196.0),
sym.Eq(nO1, 84.5)]
reactor_spec = [
sym.Eq(nE2 - nE3, 0.06*nE2)]
frcvr = sym.var('frcvr')
separator_spec = [
sym.Eq(nE5, 0.98*nE3),
sym.Eq(nO5, 0.995*nO3),
sym.Eq(nEO4, frcvr*nEO3)]
specifications = feed_spec + reactor_spec + separator_spec
for s in specifications:
display(s)
soln = sym.solve(material_balances + specifications, exclude=[frcvr])
purity = soln[nEO4]/(soln[nEO4]+soln[nE4]+soln[nO4])
display(purity)
sym.plot(purity,(frcvr,0.9,1.00),xlabel='frcvr',ylabel='fraction EO',title='Purity');
production = 24*350*(soln[nEO4] + soln[nE4] + soln[nO4])/1000000.0
display(production)
sym.plot(production,(frcvr,0.9,1.00),xlabel='frcvr',ylabel='kg-mol/year',title='Production');
Increasing the fraction of ethylene recovered for recycle should improve product purity. Let's see what happens when we raise it.
feed_spec = [
sym.Eq(nE1,196.0),
sym.Eq(nO1, 84.5)]
reactor_spec = [
sym.Eq(nE2 - nE3, 0.06*nE2)]
frcvr = sym.var('frcvr')
separator_spec = [
sym.Eq(nE5, frcvr*nE3),
sym.Eq(nO5, 0.995*nO3),
sym.Eq(nEO4, 0.97*nEO3)]
specifications = feed_spec + reactor_spec + separator_spec
for s in specifications:
display(s)
soln = sym.solve(material_balances + specifications,exclude=[frcvr])
purity = soln[nEO4]/(soln[nEO4]+soln[nE4]+soln[nO4])
display(sym.simplify(purity))
sym.plot(purity,(frcvr,0.9,1),xlabel='frcvr',ylabel='fraction EO',title='Purity');
production = 24*350*(soln[nEO4] + soln[nE4] + soln[nO4])/1000000.0
display(sym.simplify(production))
sym.plot(production,(frcvr,0.9,1.00),xlabel='frcvr',ylabel='kg-mol/year',title='Production');
As we've discovered, the oxygen feedrate is not sufficient to consume all of the ethylene. Let's explore what happens if we change the oxygen feedrate.
oxyfeed = sym.var('oxyfeed')
feed_spec = [
sym.Eq(nE1,196.0),
sym.Eq(nO1,oxyfeed)]
reactor_spec = [
sym.Eq(nE2 - nE3, 0.06*nE2)]
separator_spec = [
sym.Eq(nE5, 0.98*nE3),
sym.Eq(nO5, 0.995*nO3),
sym.Eq(nEO4, 0.97*nEO3)]
specifications = feed_spec + reactor_spec + separator_spec
for s in specifications:
display(s)
soln = sym.solve(material_balances + specifications, exlude=[oxyfeed])
purity = soln[nEO4]/(soln[nEO4]+soln[nE4]+soln[nO4])
display(sym.simplify(purity))
sym.plot(purity,(oxyfeed,74.5,196/2),xlabel='O1',ylabel='fraction EO',title='Purity');
production = 24*350*(soln[nEO4] + soln[nE4] + soln[nO4])/1000000
display(sym.simplify(production))
sym.plot(production,(oxyfeed,74.5,196/2),xlabel='O1',ylabel='kg-mol/year',title='Production');