08 - Equations to Solve Problems

08 - Equations to Solve Problems#

Solving systems of equations is a common task in engineering. In this section, we’ll review how to solve systems of equations using Python with multiple different approaches. The exampels will come from evaluating mass balances on a processing system.

Solving Linear Systems of Equations#

Separation process principles are critical to many industries. Lets say you have a flashing unit that is separating a mixture of benzene and toluene. The feed to the unit is 40 mol% benzene and 60 mol% toluene. The vapor leaving the unit is 80 mol% benzene and 20 mol% toluene. The liquid leaving the unit is 10 mol% benzene and 90 mol% toluene. The feed rate is 100 mol/hr. The vapor and liquid leaving the unit are in equilibrium. What are the flow rates of the vapor and liquid leaving the unit?

This problem can be solved by writing a system of equations and solving for the unknowns. In this case, we have two unknowns, the vapor and liquid flow rates. We can write two equations, one for the benzene balance and one for the total flows. One of the equations would be a species balance and the other equation a total mass balance. A picture of the system is shown below.

The inlet enters at the left and the vapor exists the top and the liquid exists the bottom. This is equivalent to a single stage of a distillation column. The equations are:

\[\begin{split} \begin{align} F &= L + V \\ z_i F &= x_i L + y_i V \end{align} \end{split}\]

which for this case \(z_i\) = 0.4, \(F\) = 100 mol/hr, \(x_i\) = 0.1, \(y_i\) = 0.8. We can solve this in many different ways. We’ll solve this system of equations multiple ways:

first we’ll do it using NumPy (linear algebra),
then we’ll do it using SciPy (fsolve),
then we’ll do it using SymPy (symbolic algebra),
then we’ll do it using Gauss Elimination (a numerical method),
and with Gekko (a software package for dynamic optimization).

Solve Linear System with Numpy#

The above equations can be solved numerically as shown below using the numpy package. First you’ll adjust the equations algebraically to be in the form Ax = b. Then you’ll use the numpy.linalg.solve function to solve the system of equations.

The above equation rearranged into the form Ax = b is:

\[\begin{split} \begin{align} L + V &= 100 \\ 0.1 L + 0.8 V &= 40 \end{align} \end{split}\]

The A matrix is:

\[\begin{split} \begin{bmatrix} 1 & 1 \\ 0.1 & 0.8 \end{bmatrix} \end{split}\]

and the b vector is:

\[\begin{split} \begin{bmatrix} 100 \\ 40 \end{bmatrix} \end{split}\]

import numpy as np
A = np.array([[1,1],[0.1,0.8]])
b = np.array([100,40])

x = np.linalg.solve(A,b)
print(f'The solution is {x} for the system of equations with the first value equal to L and the second value equal to V')

The solution is [57.14285714 42.85714286] for the system of equations with the first value equal to L and the second value equal to V

Solve Linear System with Scipy’s Fsolve#

The above equations can be solved numerically as shown below using the scipy package. Fsolve can be used to solve for multiple variables in both a non-linear and linear scenario given there are continuous functions and the number of equations = number of unknowns. Method used to find roots is in the documentation: fsolve is a wrapper around MINPACK’s hybrd and hybrj algorithms.

from scipy.optimize import fsolve
def eq (vars):
    L,V = vars #unpack the tuple
    # equation is set equal to zero (eq1, eq2 etc)
    eq1 = 100 - (L+V)
    eq2 = 40-0.1*L-0.8*V
    #you return what eq1 and eq2 equal given the L and V params
    # fsolve will go through an algorithm changing L and V values 
    # until both eq1 and eq2 will be equal to zero 
    return eq1,eq2

# this is a list of guess values
guesses = [1,1]


# eq(guess)
l,v = fsolve(eq,guesses)
print(f'L  = {l:0.2f} and V = {v:0.2f}')
## Code credit to Jacob Morisset

L  = 57.14 and V = 42.86

Remember thqt the notation for fsolve is:

fsolve(function name(with input only related to the guess value(s)), guess value(s)). Further examples in the next lecture.

Solve Linear System with Sympy#

The above equations can be solved symbolically as shown below using the sympy package. Typically used for symbolic solutions to equations. The sympy package can also be used to solve differential equations symbolically. The solver is for linear equations or nonlinear equations, see the documentation for further details: https://docs.sympy.org/latest/modules/solvers/solvers.html.

import sympy as sym
from IPython.display import display, Math
L = sym.symbols('L')
V = sym.symbols('V')
ans = sym.solve([100 - (L+V), 0.4*100-0.1*L-0.8*V], [L,V])
#or equiavlently you can use 
#ans = sym.linsolve([100 - (L+V), 0.4*100-0.1*L-0.8*V],[L,V])
display(Math(sym.latex(ans)))

\[\displaystyle \left\{ L : 57.1428571428571, \ V : 42.8571428571429\right\}\]

Solve Linear System with Gaussian Elimination#

Linear algebra techniques can be used to take the inverse of the A matrix and multiply it by the b vector to solve for the unknowns. One way to get the inverse of the A matrix is to use the Gauss Elimination method. This method is shown below. Similar to the linalg solution above, the system of equations are arranged in the form Ax = b.

import numpy as np

# Define the augmented matrix [A|b]
A = np.array([[1, 1],
                [0.1, 0.8]])

b = np.array([100, 40])

# Augmented matrix dimensions
n_rows, n_cols = A.shape

# Forward elimination
for k in range(n_rows - 1):
    for i in range(k + 1, n_rows):
        factor = A[i, k] / A[k, k]
        A[i, k:] -= factor * A[k, k:]
        b[i] -= factor * b[k]

# Back substitution
x = np.zeros(n_rows)
for i in range(n_rows - 1, -1, -1):
    x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

print("Solution:")
print(x)

Solution:
[57.14285714 42.85714286]

Solve Linear System with Gekko#

Gekko is a software package for dynamic optimization. It can also be used to solve systems of equations. Developed here at Brigham Young University by John Hedengren. Further information here: https://gekko.readthedocs.io/en/latest/ and here: https://www.mdpi.com/2227-9717/6/8/106.

#Solve equation with Gekko
from gekko import GEKKO
m = GEKKO() # Initialize gekko
L = m.Var(value=1) # Initialize variables
V = m.Var(value=1)
m.Equations([L+V==100,0.1*L+0.8*V==40]) # Equations
m.solve(disp=False) # Solve
print(f'The solution is {L.value[0]:0.2f} for L and {V.value[0]:0.2f} for V')

The solution is 57.14 for L and 42.86 for V

The above answers are consistent with the answers from the other methods.

Solving Non-Linear Systems of Equations#

Similar tools can be used to solve non-linear systems of equations. Scipy’s fsolve is commonly used. There are other tools including Gekko, and Newton’s method.

One application of solving a non-linear equatoin in chemical engineering is calculating the pressure loss from flow down a pipe using the Colbrook equation to determine the Darcy friction factor. The Darcy friction factor is used to calculate the pressure drop in a pipe. The Darcy friction factor is a function of the Reynolds number and the relative roughness of the pipe. The Reynolds number is a function of the fluid velocity, density, viscosity, and pipe diameter.

#Colebrook equation, Re is Reynolds number, D is pipe diameter, epsilon is the pipe roughness
def Colbrook(f,Re,D,epsil):
    return 1/np.sqrt(f) + 2*np.log10(2.51/(Re*np.sqrt(f))+epsil/(3.71*D)) 

The Colbrook equation is:

\[ \begin{equation} \frac{1}{\sqrt{f}} = -2.0 \log_{10} \left( \frac{\epsilon/D}{3.71} + \frac{2.51}{\text{Re} \sqrt{f}} \right) \end{equation} \]

where \(f\) is the darcy friction factor, \(\epsilon\) is the pipe roughness, \(D\) is the pipe diameter, and Re is the Reynolds number equal to the product of the fluid density, velocity, and diameter divided by the fluid viscosity. This is valid with a relative roughness of \(\epsilon/D < 0.05\) and \(10^4 < \text{Re} < 10^8\).

The pressure loss down the length of pipe is calculated by:

\[ \begin{equation} \Delta P = f \frac{L}{D} \frac{\rho v^2}{2} \end{equation} \]

where \(\Delta P\) is the pressure drop, \(f\) is the Darcy friction factor, \(L\) is the length of the pipe, \(D\) is the pipe diameter, \(\rho\) is the fluid density, and \(v\) is the fluid velocity.

#pressure loss down a length of pipe in terms of the friction factor
def pressure_loss(f,L,D,rho,v):
    return f*L/D*rho*v**2/2

Example problem: calculate the pressure drop down a 100 m long pipe with a 2 inch diameter with a fluid velocity of 1 m/s. The fluid has a density of 1000 kg/m\(^3\) and a viscosity of 0.001 kg/m-s. The pipe roughness is 0.0001 m. What is the pressure drop?

Non-linear solution with scipy fsolve#

#first calculate the Reynolds number
def Reynolds(rho,v,D,mu):
    return rho*v*D/mu
#specify pressure losses
def pressure_loss(f,L,D,rho,v):
    return f*L/D*rho*v**2/2
#set constants
rho = 1000 #kg/m^3
vel = 1 #m/s
pipeDia = 2*0.0254 #inch to m
mu = 0.001 #Pa*s
epsil = 0.0001 #m
pipeLen = 100 #m
#solve for the friction factor
ffact =  fsolve(lambda f: Colbrook(f,Reynolds(rho,vel,pipeDia,mu),pipeDia,epsil), 0.01)[0]
#calculate the pressure drop
ploss = pressure_loss(ffact,pipeLen,pipeDia,rho,vel) #Pa
print(f'The Reynolds number is {Reynolds(rho,vel,pipeDia,mu):.0f}')
print(f'The friction factor is {ffact:.4f}')
print(f'The pressure drop is {ploss/1000:.2f} Pa (or {ploss/1000*0.145:.2f} psi)')

The Reynolds number is 50800
The friction factor is 0.0264
The pressure drop is 25.97 Pa (or 3.76 psi)

Non-linear solution with Gekko#

Solve for the friction factor using the GEKKO package.

#Solve for the friction factor using Gekko
m = GEKKO() # Initialize the gekko model
f = m.Var(value=0.01) # Initialize variables
m.Equation(1/m.sqrt(f) + 2*m.log10(2.51/(Reynolds(rho,vel,pipeDia,mu)*m.sqrt(f))+epsil/(3.71*pipeDia))==0) # Equations
m.solve(disp=False) # Solve
print(f'The friction factor is {f.value[0]:.4f}')

The friction factor is 0.0264