Problem Set 2#

CBE 60258, University of Notre Dame. © Prof. Alexander Dowling, 2023

You may not distribution homework solutions without written permissions from Prof. Alexander Dowling.

Your Name and Email:

Import Libraries, Define Necessary Functions#

import numpy as np
import matplotlib.pyplot as plt

'''
This cell installs Pyomo and Ipopt on Google Colab. To run this notebook
locally (e.g., with anaconda), you first need to install Pyomo and Ipopt.
'''

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/IDAES/idaes-pse/main/scripts/colab_helper.py"
    import colab_helper
    colab_helper.install_idaes()
    colab_helper.install_ipopt()

import pyomo.environ as pyo

def Jacobian(f,x,delta = 1.0e-7):
    '''Approximate Jacobian using forward finite difference

    Args:
        f: vector-valued function
        x: point to build approximation J(x) around
        delta: finite difference step size

    Returns:
        J: square Jacobian matrix (approximation)

    '''
    # Determine size
    N = x.size
    
    #Evaluate function f at x
    fx = f(x) #only need to evaluate this once
    
    # Make sure f is square (no. inputs = no. outputs)
    try:
        assert N == fx.size, "Your problem is not square!"
    except AssertionError:
        print("x = ",x)
        print("fx = ",fx)
    
    
    # Allocate empty matrix
    J = np.zeros((N,N))

    idelta = 1.0/delta #division is slower than multiplication
    x_perturbed = x.copy() #copy x to add delta
    
    # Loop over elements of x and columns of J
    for i in range(N):
        # Perturb (apply step) to element i of x
        x_perturbed[i] += delta
        
        # Approximate column in Jacobian
        col = (f(x_perturbed) - fx) * idelta
        
        # Reset element of x
        x_perturbed[i] = x[i]
        
        # Save results
        J[:,i] = col
    # end for loop
    return J

def newton_system(f,x0,exact_Jac=None,delta=1E-7,epsilon=1.0e-6, LOUD=False):
    """Find the root of the function f via exact or inexact Newton-Raphson method
    Args:
        f: function to find root of
        x0: initial guess
        exact_Jac: function to calculate J. If None, use finite difference
        delta: finite difference tolerance. Only used if J is not specified
        epsilon: tolerance
        
    Returns:
        estimate of root
    """
        
    x = x0
    if (LOUD):
        print("x0 =",x0)
    iterations = 0
    fx = f(x)
    while (np.linalg.norm(fx) > epsilon):
        if exact_Jac is None:
            # Use finite difference
            J = Jacobian(f,x,delta)
        else:
            J = exact_Jac(x)
        
        RHS = -fx;
        
        # solve linear system
        # We could have used GaussElimPivotSolve here instead
        delta_x = np.linalg.solve(J,RHS)
        
        # Check if GE returned any NaN values
        if np.isnan(delta_x).any():
            print("Gaussian Elimination Failed!")
            print("J = \n",J,"\n")
            print("J is rank",np.linalg.matrix_rank(J),"\n")
            print("RHS = ",RHS)
        x = x + delta_x
        fx = f(x)
        if (LOUD):
            print("\nIteration",iterations+1,": x =",x,"\n norm(f(x)) =",np.linalg.norm(fx),"\t cond(J(x)) =",np.linalg.cond(J))
        iterations += 1
    print("\nIt took",iterations,"iterations")
    return x #return estimate of root

1. Solubility of Carbon Dioxide in Aqueous NaOH#

The solubility of \(\mathrm{CO_2}\) in water is a strong function of the pH, due to the conversion of dissolved \(\mathrm{CO_2}\) to carbonic acid (\(\mathrm{H_2CO_3}\)) and its dissociation into \(\mathrm{HCO_3^-}\) and \(\mathrm{CO_3^{2-}}\).

In this problem, you will predict the equilibrium concentration of all dissolved species in a 0.01 M NaOH solution. Assuming that the NaOH completely dissociates (it is a very strong base), we want to calculate the pH and species concentrations for a solution in equilibrium with 0.2 atm \(\mathrm{CO_2}\) gas (i.e. \(\mathrm{P_{CO_2}}\) = 0.2 atm). This would roughly be the concentration of \(\mathrm{CO_2}\) in a stoichiometric burning of coal in air with all oxygen consumed. This chemistry is important for carbon capture from a power plant with a base solution.

We have the following equilibrium relationships:

\[\frac{\lbrack CO_2 \rbrack}{P_{CO_2} } = K_{abs}\]

\[\frac{\lbrack H_2CO_3 \rbrack}{\lbrack CO_2 \rbrack } = K_{eq}\]

\[\frac{\lbrack HCO_3^- \rbrack \times \lbrack H^+ \rbrack}{\lbrack H_2CO_3 \rbrack } = K_{a_1}\]

\[\frac{\lbrack CO_3^{2-} \rbrack \times \lbrack H^+ \rbrack}{\lbrack HCO_3^- \rbrack } = K_{a_2}\]

\[\lbrack H^+ \rbrack \times \lbrack OH^- \rbrack = K_{w}\]

The equilibrium constants have the following values:

Constant	Value	Units
\(K_{abs}\)	\(3.36 \times 10^{-2}\)	M/atm
\(K_{eq}\)	\(1.7 \times 10^{-3}\)	-
\(K_{a_1}\)	\(2.5 \times 10^{-4}\)	M
\(K_{a_2}\)	\(4.8 \times 10^{-11}\)	M
\(K_{w}\)	\(1.00 \times 10^{-14}\)	M\(^2\)

Further, the electroneutrality condition for the system is given by the equation:

(5)#\[\begin{equation} \lbrack Na^+ \rbrack + \lbrack H^+ \rbrack = \lbrack OH^- \rbrack + \lbrack HCO_3^- \rbrack + 2\lbrack CO_3^{2-} \rbrack \end{equation}\]

1a. Naïve Model#

Write down the system of nonlinear equations in canonical form to predict the equilibrium concentrations. Clearly define all the variables that you use. In 2c you will refine the model to improve condition, so do not transform the the variable in this part.

Hint: Start by performing a degree of freedom analysis.

1b. Numeric Solution with Python#

Numerically solve your model from 2a. We will see in this problem there are many different answers depending on the order of equations and the initial point. This is because the problem is nonlinear and exacerbated by the fact our naïve model is poorly scaled.

# lists of species 
species = ["H", "HCO3", "CO3"]
other = ["CO2", "H2CO3", "OH"]

# Provided specification
K_abs = 3.36e-2 # M/atm
K_eq = 1.7e-3
K_a1 = 2.5e-4 # M
K_a2 = 4.8e-11 #M
K_w = 1e-14 # M^2

# CO2 pressure
P_co2 = 0.2 # atm

# Molarity of solution is Na+ concentration
Na_con = 0.01 # M

# define the set of equations
def my_equations_1(x):
    ''' Nonlinear system of equations in conancial form c(x) = 0
    
    Arg:
        x: vector of variables
        
    Returns:
        r: residual, 
    
    '''
    # allocate an empty array
    r = np.empty((3))
  
    # equations
    
    # Add your solution here
  
    return r

Specify an initial point and solve your model in Python.

# Add your solution here

Specify a different initial point and solve your model in Python. Hint: You might get a different answer with a different initial point.

# Add your solution here

1c. Refined Model#

Reformulate your model from a to a model that is better conditioned.

1d. Numeric Solution with Python#

Numerically solve your model from c. First, define the refined model in canonical form.

# lists of species
species = ["H", "HCO3", "CO3"]
other = ["CO2", "H2CO3", "OH"]

# Provided specification
K_abs = 3.36e-2 # M/atm
K_eq = 1.7e-3
K_a1 = 2.5e-4 # M
K_a2 = 4.8e-11 #M
K_w = 1e-14 # M^2

# CO2 pressure
P_co2 = 0.2 # atm

# Molarity of solution is Na+ concentration
Na_con = 0.01 # M

# define the set of equations
def my_equations2(x):
    ''' Nonlinear system of equations in conancial form c(x) = 0
    Arg:
        x: vector of variables
    Returns:
        r: residual
    '''
    # allocate an empty array
    r = np.empty((3))
    
    # equations
    
    # Add your solution here
    
    return r

Next, specify one initial point and solve the model.

# Add your solution here

Now, try another initial point.

# Add your solution here

Record your answers for this problem in the template provided in Canvas and turn in via Gradescope.

1e. Discussion#

Why is transforming the equations so important for this problem? Write 2 to 4 sentences. For full credit, please discuss i) condition number, ii) scaling, and iii) convergence tolerance for Newton’s method.

Your Answer:

2. Numeric Integration: Reaction Kinetics#

The species A undergoes a dimerization reaction:

(6)#\[\begin{equation} 2A\rightarrow B \end{equation}\]

With the rate:

(7)#\[\begin{equation} r_1 = k_1 \cdot C_A^2 \end{equation}\]

The reaction continues, however, with the further reaction:

(8)#\[\begin{equation} B+A\rightarrow C \end{equation}\]

With the rate:

(9)#\[\begin{equation} r_2 = k_2 \cdot C_A \cdot C_B \end{equation}\]

We would like to maximize the product B, so we need to stop the reaction when this is at a maximum. We will consider \(k_1=2\) and \(k_2=4\) L/(mol-s) with the initial conditions \(C_A=1\) and \(C_B=0\) mol/L.

As seniors, you will learn how to apply a mass balance to convert reaction rate expressions into differential equations:

(10)#\[\begin{equation} \frac{dC_A}{dt}=-2\cdot r_1 - r_2 = -2\cdot k_1 \cdot C_A^2-k_2 \cdot C_A \cdot C_B \end{equation}\]

(11)#\[\begin{equation} \frac{dC_B}{dt}=r_1 - r_2 = k_1 \cdot C_A^2-k_2 \cdot C_A \cdot C_B \end{equation}\]

Do the following:

Calculate the concentrations \(C_A\) and \(C_B\) as a function of time.
Plot the concentrations from \(t=0\) to \(t=3\).
Estimate the maximum concentration of B and when it occurs. Mark this point on your plot. Save this value in the float maxCB.

Hint: Because this is not system of linear differential equations, we should use scipy.integrate.solve_ivp. Check your notes from Class 12 and/or explore the Scipy documentation for examples: https://docs.scipy.org/doc/scipy-1.1.0/reference/generated/scipy.integrate.solve_ivp.html

#Solve ODES
def f(t, y):
    ''' RHS of differential equation for reaction kinetics
    Args:
        t: time
        y: values for differential questions, [CA, CB]
    Returns:
        dydt: first derivation of y w.r.t. t
    
    '''
    # unpack current values of y
    CA, CB = y      
    
    # define kinetic parameters
    k1 = 2
    k2 = 4
    
    #Store equations as a list in the variable dydt
    
    # Add your solution here
    
    return dydt

# Initial values
CA0 = 1.0     
CB0 = 0.0 

# Bundle initial conditions for ODE solver
y0 = [CA0, CB0]

t = np.arange(0, 3, 0.01)
tspan = [np.min(t), np.max(t)]

# Call the ODE solver
# Note: The solution profile changes slightly if we change the integration method.
soln = integrate.solve_ivp(f, tspan, y0, t_eval=t, method='RK23')

# You may wish to uncomment the following line to take a closer look.
# print(soln)

CA = soln.y[0]
CB = soln.y[1]

#Find maximum concentration of B
# Add your solution here

print("maxCB =",maxCB)

#Plot
# Add your solution here

# Removed autograder test. You may delete this cell.

3. Numeric Integration: Cannonball#

An iron cannonball of mass 1kg is fired upwards at a velocity of 200m/s. In addition to gravity, there is a drag due to air resistance given by \(F_{D} = -0.00292 U|U|\) where \(U\) is \(dh/dt\) (all units SI). Note that we have to use \(U|U|\) rather than \(U^2\) because we have to keep track of the sign.

Do the following:

Use Newton’s laws to construct a pair of coupled first order differential equations to track the velocity and height of the cannon ball. Write this system on paper. Turn this in via Gradescope.
Integrate the model in Python. Hint: How should you integrate this numerically?
Plot the height and velocity as a function of time. Determine the maximum elevation of the cannonball and mark it on the graph. Save this value in the float maxh. Then print the value.

For a token amount of extra credit, prepare your plot with velocity on the left vertical axis and elevation on the right vertical axis. See this example for starter code: https://matplotlib.org/gallery/api/two_scales.html

# Add your solution here

# Removed autograder test. You may delete this cell.

4. Numeric Integration of Partial Differential Equations with Pyomo#

%matplotlib inline

# Import plotting libraries
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D 

# Import Pyomo
import pyomo.environ as pyo

# Import Pyomo numeric integration features
from pyomo.dae import DerivativeVar, ContinuousSet

During your time at Notre Dame, you will likely want (or at least need) to solve a partial differential equation (PDE) system. In this problem, we will practice using Pyomo to numerically integrate a simple and common PDE. (Special thanks to Prof. Kantor for this problem.)

Transport of heat in a solid is described by the familiar thermal diffusion model:

\[ \begin{align*} \rho C_p\frac{\partial T}{\partial t} & = \nabla\cdot(k\nabla T) \end{align*} \]

We will assume the thermal conductivity \(k\) is a constant, and define thermal diffusivity in the conventional way

\[ \begin{align*} \alpha & = \frac{k}{\rho C_p} \end{align*} \]

We will further assume symmetry with respect to all spatial coordinates except \(x\) where \(x\) extends from \(-X\) to \(+X\). The boundary conditions are

\[\begin{split} \begin{align*} T(t,X) & = T_{\infty} & \forall t > 0 \\ \nabla T(t,0) & = 0 & \forall t \geq 0 \end{align*} \end{split}\]

where we have assumed symmetry with respect to \(r\) and uniform initial conditions \(T(0, x) = T_0\) for all \(0 \leq r \leq X\).

4a. Rescaling and Dimensionless Model#

We would like a dimensionless model for two reasons: first, we only need to solve the dimensionless model once, i.e., it becomes independent of input data. Second, the dimensionless models are often scaled better for numerical solutions.

Let’s consider the following proposed scaling procedure:

\[\begin{split} \begin{align*} T' & = \frac{T - T_0}{T_\infty - T_0} \\ x' & = \frac{r}{X} \\ t' & = t \frac{\alpha}{X^2} \end{align*} \end{split}\]

Show this scaling procedure gives the following dimensionless system:

\[ \begin{align*} \frac{\partial T'}{\partial t'} & = \nabla^2 T' \end{align*} \]

with auxiliary conditions

\[\begin{split} \begin{align*} T'(0, x') & = 0 & \forall 0 \leq x' \leq 1\\ T'(t', 1) & = 1 & \forall t' > 0\\ \nabla T'(t', 0) & = 0 & \forall t' \geq 0 \\ \end{align*} \end{split}\]

Turn in your work (pencil and paper) via Gradescope. Important: Here the prime \('\) indicates the scaled variables and coordinates. It does not indicate a derivative. Thus \(T'\) is scaled temperature, NOT the derivative of temperature (which begs the question of “with respect to what?”).

4b. Numeric Integration via Pyomo#

For simplicity, let’s consider planar coordinates. For a slab geometry, we want to numerical integrate the following PDE:

\[ \begin{align*} \frac{\partial T'}{\partial t'} & = \frac{\partial^2 T'}{\partial x'^2} \end{align*} \]

with auxiliary conditions

\[\begin{split} \begin{align*} T'(0, x') & = 0 & \forall 0 \leq x' \leq 1 \\ T'(t', 1) & = 1 & \forall t' > 0\\ \frac{\partial T'}{\partial x'} (t', 0) & = 0 & \forall t' \geq 0 \\ \end{align*} \end{split}\]

Complete the following Pyomo code to integrate this PDE.

# Create Pyomo model
m = pyo.ConcreteModel()

# Define sets for spatial and temporal domains
m.x = ContinuousSet(bounds=(0,1))
m.t = ContinuousSet(bounds=(0,2))

# Define scaled temperature indexed by time and space
m.T = pyo.Var(m.t, m.x)

# Define variables for the derivates
m.dTdt   = DerivativeVar(m.T, wrt=m.t)
m.dTdx   = DerivativeVar(m.T, wrt=m.x)
m.d2Tdx2 = DerivativeVar(m.T, wrt=(m.x, m.x))

# Define PDE equation
def pde(m, t, x):
    if t == 0:
        return pyo.Constraint.Skip
    elif x == 0 or x == 1:
        return pyo.Constraint.Skip
    # Add your solution here

m.pde = pyo.Constraint(m.t, m.x, rule=pde)

# Define first auxilary condition
def initial_condition(m, x):
    if x == 0 or x == 1:
        return pyo.Constraint.Skip
    # Add your solution here

m.ic = pyo.Constraint(m.x, rule = initial_condition)

# Define second auxilary condition
def boundary_condition1(m, t):
    # Add your solution here

m.bc1 = pyo.Constraint(m.t, rule = boundary_condition1)

# Define third auxilary condition
@m.Constraint(m.t)
def boundary_condition2(m, t):
    # Add your solution here  

m.bc2 = pyo.Constraint(m.t, rule=boundary_condition2)

# Define dummy objective
m.obj = pyo.Objective(expr=1)

# Discretize spatial coordinate with forward finite difference and 50 elements
pyo.TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, scheme='FORWARD', wrt=m.x)

# Discretize time coordinate with forward finite difference and 50 elements
pyo.TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, scheme='FORWARD', wrt=m.t)
pyo.SolverFactory('ipopt').solve(m, tee=True).write()

4c. Visualize Solution#

# Extract indices
x = sorted(m.x)
t = sorted(m.t)

# Create numpy arrays to hold the solution
xgrid = np.zeros((len(t), len(x)))
# Hint: define tgrid and Tgrid the same way
# Add your solution here

# Loop over time
for i in range(0, len(t)):
    # Loop over space
    for j in range(0, len(x)):
        # Copy values
        xgrid[i,j] = x[j]
        tgrid[i,j] = t[i]
        # Hint: how to access values from Pyomo variable m.T?
        # Add your solution here

# Create a 3D wireframe plot of the solution
# Hint: consult the matplotlib documentation
# https://matplotlib.org/stable/gallery/mplot3d/wire3d.html

# Add your solution here

Write a few sentences to describe the PDE solution. Is it what you expect based on your prior knowledge of this system? Each person brings different prior knwoledge to this class, you everyone should have a distinct answer. In other words, there is no “right answer”. Instead, this is helping you practice interpreting results based on your knowledge which is a critical skill in graduate school.

Discussion:

Submission Instructions and Tips#

Answer discussion questions in this notebook.
When asked to store a solution in a specific variable, please also print that variable.
Turn in this notebook via Gradescope.
Also turn in written (pencil and paper) work via Gradescope.
Even if you are not required to turn in pseudocode, you should always write pseudocode. It takes only a few minutes and can save you hours of time.
We are not using the autograder for CBE 60258, so please skip those instructions.

Data and Computing for Chemical Engineers

Problem Set 2

Contents