The Kitchin Research Group

If there is a pressure drop in a plug flow reactor, ¹ there are two equations needed to determine the exit conversion: one for the conversion, and one from the pressure drop.

Here is how to integrate these equations numerically in python.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

kprime = 0.0266
Fa0 = 1.08
alpha = 0.0166
epsilon = -0.15

def dFdW(F, W):
    'set of ODEs to integrate'
    X = F[0]
    y = F[1]
    dXdW = kprime / Fa0 * (1-X) / (1 + epsilon*X) * y
    dydW = -alpha * (1 + epsilon * X) / (2 * y)
    return [dXdW, dydW]

Wspan = np.linspace(0,60)
X0 = 0.0
y0 = 1.0
F0 = [X0, y0]
sol = odeint(dFdW, F0, Wspan)

# now plot the results
plt.plot(Wspan, sol[:,0], label='Conversion')
plt.plot(Wspan, sol[:,1], 'g--', label='y=$P/P_0$')
plt.legend(loc='best')
plt.xlabel('Catalyst weight (lb_m)')
plt.savefig('images/2013-01-08-pdrop.png')

Here is the resulting figure.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

One approach to solving BVPs is to use the shooting method. The reason we cannot use an initial value solver for a BVP is that there is not enough information at the initial value to start. In the shooting method, we take the function value at the initial point, and guess what the function derivatives are so that we can do an integration. If our guess was good, then the solution will go through the known second boundary point. If not, we guess again, until we get the answer we need. In this example we repeat the pressure driven flow example, but illustrate the shooting method.

In the pressure driven flow of a fluid with viscosity \(\mu\) between two stationary plates separated by distance \(d\) and driven by a pressure drop \(\Delta P/\Delta x\), the governing equations on the velocity \(u\) of the fluid are (assuming flow in the x-direction with the velocity varying only in the y-direction):

$$\frac{\Delta P}{\Delta x} = \mu \frac{d^2u}{dy^2}$$

with boundary conditions \(u(y=0) = 0\) and \(u(y=d) = 0\), i.e. the no-slip condition at the edges of the plate.

we convert this second order BVP to a system of ODEs by letting \(u_1 = u\), \(u_2 = u_1'\) and then \(u_2' = u_1''\). This leads to:

\(\frac{d u_1}{dy} = u_2\)

\(\frac{d u_2}{dy} = \frac{1}{\mu}\frac{\Delta P}{\Delta x}\)

with boundary conditions \(u_1(y=0) = 0\) and \(u_1(y=d) = 0\).

for this problem we let the plate separation be d=0.1, the viscosity \(\mu = 1\), and \(\frac{\Delta P}{\Delta x} = -100\).

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

d = 0.1 # plate thickness

def odefun(U, y):
    u1, u2 = U
    mu = 1
    Pdrop = -100
    du1dy = u2
    du2dy = 1.0 / mu * Pdrop
    return [du1dy, du2dy]

u1_0 = 0 # known
u2_0 = 1 # guessed

dspan = np.linspace(0, d)

U = odeint(odefun, [u1_0, u2_0], dspan)

plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.savefig('images/bvp-shooting-1.png')

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

d = 0.1 # plate thickness

def odefun(U, y):
    u1, u2 = U
    mu = 1
    Pdrop = -100
    du1dy = u2
    du2dy = 1.0 / mu * Pdrop
    return [du1dy, du2dy]

u1_0 = 0 # known
u2_0 = 10 # guessed

dspan = np.linspace(0, d)

U = odeint(odefun, [u1_0, u2_0], dspan)

plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.savefig('images/bvp-shooting-2.png')

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

d = 0.1 # plate thickness
Pdrop = -100
mu = 1

def odefun(U, y):
    u1, u2 = U
    du1dy = u2
    du2dy = 1.0 / mu * Pdrop
    return [du1dy, du2dy]

u1_0 = 0 # known
dspan = np.linspace(0, d)

def objective(u2_0):
    dspan = np.linspace(0, d)
    U = odeint(odefun, [u1_0, u2_0], dspan)
    u1 = U[:,0]
    return u1[-1]

u2_0, = fsolve(objective, 1.0)

# now solve with optimal u2_0
U = odeint(odefun, [u1_0, u2_0], dspan)

plt.plot(dspan, U[:,0], label='Numerical solution')
plt.plot([d],[0], 'ro')

# plot an analytical solution
u = -(Pdrop) * d**2 / 2 / mu * (dspan / d - (dspan / d)**2)
plt.plot(dspan, u, 'r--', label='Analytical solution')


plt.xlabel('d')
plt.ylabel('$u_1$')
plt.legend(loc='best')
plt.savefig('images/bvp-shooting-3.png')

from pycse import BVP_sh
import matplotlib.pyplot as plt

d = 0.1 # plate thickness
Pdrop = -100
mu = 1

def odefun(U, y):
    u1, u2 = U
    du1dy = u2
    du2dy = 1.0 / mu * Pdrop
    return [du1dy, du2dy]

x1 = 0.0; alpha = 0.0
x2 = 0.1; beta = 0.0
init = 2.0 # initial guess of slope at x=0

X,Y = BVP_sh(odefun, x1, x2, alpha, beta, init)
plt.plot(X, Y[:,0])
plt.ylim([0, 0.14])

# plot an analytical solution
u = -(Pdrop) * d**2 / 2 / mu * (X / d - (X / d)**2)
plt.plot(X, u, 'r--', label='Analytical solution')
plt.savefig('images/bvp-shooting-4.png')
plt.show()

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

Adapted from http://www.physics.arizona.edu/~restrepo/475B/Notes/sourcehtml/node24.html

We want to solve a linear boundary value problem of the form: y'' = p(x)y' + q(x)y + r(x) with boundary conditions y(x1) = alpha and y(x2) = beta.

For this example, we solve the plane poiseuille flow problem using a finite difference approach. An advantage of the approach we use here is we do not have to rewrite the second order ODE as a set of coupled first order ODEs, nor do we have to provide guesses for the solution. We do, however, have to discretize the derivatives and formulate a linear algebra problem.

we want to solve u'' = 1/mu*DPDX with u(0)=0 and u(0.1)=0. for this problem we let the plate separation be d=0.1, the viscosity \(\mu = 1\), and \(\frac{\Delta P}{\Delta x} = -100\).

The idea behind the finite difference method is to approximate the derivatives by finite differences on a grid. See here for details. By discretizing the ODE, we arrive at a set of linear algebra equations of the form \(A y = b\), where \(A\) and \(b\) are defined as follows.

\[A = \left [ \begin{array}{ccccc} % 2 + h^2 q_1 & -1 + \frac{h}{2} p_1 & 0 & 0 & 0 \\ -1 - \frac{h}{2} p_2 & 2 + h^2 q_2 & -1 + \frac{h}{2} p_2 & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & -1 - \frac{h}{2} p_{N-1} & 2 + h^2 q_{N-1} & -1 + \frac{h}{2} p_{N-1} \\ 0 & 0 & 0 & -1 - \frac{h}{2} p_N & 2 + h^2 q_N \end{array} \right ] \]

\[ y = \left [ \begin{array}{c} y_i \\ \vdots \\ y_N \end{array} \right ] \]

\[ b = \left [ \begin{array}{c} -h^2 r_1 + ( 1 + \frac{h}{2} p_1) \alpha \\ -h^2 r_2 \\ \vdots \\ -h^2 r_{N-1} \\ -h^2 r_N + (1 - \frac{h}{2} p_N) \beta \end{array} \right] \]

import numpy as np

# we use the notation for y'' = p(x)y' + q(x)y + r(x)
def p(x): return 0
def q(x): return 0
def r(x): return -100

#we use the notation y(x1) = alpha and y(x2) = beta

x1 = 0; alpha = 0.0
x2 = 0.1; beta = 0.0

npoints = 100

# compute interval width
h = (x2-x1)/npoints;

# preallocate and shape the b vector and A-matrix
b = np.zeros((npoints - 1, 1));
A = np.zeros((npoints - 1, npoints - 1));
X = np.zeros((npoints - 1, 1));

#now we populate the A-matrix and b vector elements
for i in range(npoints - 1):
    X[i,0] = x1 + (i + 1) * h

    # get the value of the BVP Odes at this x
    pi = p(X[i])
    qi = q(X[i])
    ri = r(X[i])

    if i == 0:
        # first boundary condition
        b[i] = -h**2 * ri + (1 + h / 2 * pi)*alpha; 
    elif i == npoints - 1:
        # second boundary condition
        b[i] = -h**2 * ri + (1 - h / 2 * pi)*beta; 
    else:
        b[i] = -h**2 * ri # intermediate points
    
    for j in range(npoints - 1):
        if j == i: # the diagonal
            A[i,j] = 2 + h**2 * qi
        elif j == i - 1: # left of the diagonal
            A[i,j] = -1 - h / 2 * pi
        elif j == i + 1: # right of the diagonal
            A[i,j] = -1 + h / 2 * pi
        else:
            A[i,j] = 0 # off the tri-diagonal
 
# solve the equations A*y = b for Y
Y = np.linalg.solve(A,b)

x = np.hstack([x1, X[:,0], x2])
y = np.hstack([alpha, Y[:,0], beta])

import matplotlib.pyplot as plt

plt.plot(x, y)

mu = 1
d = 0.1
x = np.linspace(0,0.1);
Pdrop = -100 # this is DeltaP/Deltax
u = -(Pdrop) * d**2 / 2.0 / mu * (x / d - (x / d)**2)
plt.plot(x,u,'r--')

plt.xlabel('distance between plates')
plt.ylabel('fluid velocity')
plt.legend(('finite difference', 'analytical soln'))
plt.savefig('images/pp-bvp-fd.png')
plt.show()

You can see excellent agreement here between the numerical and analytical solution.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post A heat exchanger must handle 2.5 L/s of water through a smooth pipe with length of 100 m. The pressure drop cannot exceed 103 kPa at 25 degC. Compute the minimum pipe diameter required for this application.

Adapted from problem 8.8 in Problem solving in chemical and Biochemical Engineering with Polymath, Excel, and Matlab. page 303.

We need to estimate the Fanning friction factor for these conditions so we can estimate the frictional losses that result in a pressure drop for a uniform, circular pipe. The frictional forces are given by \(F_f = 2f_F \frac{\Delta L v^2}{D}\), and the corresponding pressure drop is given by \(\Delta P = \rho F_f\). In these equations, \(\rho\) is the fluid density, \(v\) is the fluid velocity, \(D\) is the pipe diameter, and \(f_F\) is the Fanning friction factor. The average fluid velocity is given by \(v = \frac{q}{\pi D^2/4}\).

For laminar flow, we estimate \(f_F = 16/Re\), which is a linear equation, and for turbulent flow (\(Re > 2100\)) we have the implicit equation \(\frac{1}{\sqrt{f_F}}=4.0 \log(Re \sqrt{f_F})-0.4\). Of course, we define \(Re = \frac{D v\rho}{\mu}\) where \(\mu\) is the viscosity of the fluid.

It is known that \(\rho(T) = 46.048 + 9.418 T -0.0329 T^2 +4.882\times10^{-5}-2.895\times10^{-8}T^4\) and \(\mu = \exp\left({-10.547 + \frac{541.69}{T-144.53}}\right)\) where \(\rho\) is in kg/m^3 and \(\mu\) is in kg/(m*s).

The aim is to find \(D\) that solves: \(\Delta p = \rho 2 f_F \frac{\Delta L v^2}{D}\). This is a nonlinear equation in \(D\), since D affects the fluid velocity, the Re, and the Fanning friction factor. Here is the solution

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

T = 25 + 273.15
Q = 2.5e-3       # m^3/s
deltaP = 103000  # Pa
deltaL = 100     # m

#Note these correlations expect dimensionless T, where the magnitude
# of T is in K

def rho(T):
    return 46.048 + 9.418 * T -0.0329 * T**2 +4.882e-5 * T**3 - 2.895e-8 * T**4

def mu(T):
    return np.exp(-10.547 + 541.69 / (T - 144.53))

def fanning_friction_factor_(Re):
    if Re < 2100:
        raise Exception('Flow is probably not turbulent, so this correlation is not appropriate.')
    # solve the Nikuradse correlation to get the friction factor
    def fz(f): return 1.0/np.sqrt(f) - (4.0*np.log10(Re*np.sqrt(f))-0.4)
    sol, = fsolve(fz, 0.01)
    return sol

fanning_friction_factor = np.vectorize(fanning_friction_factor_)

Re = np.linspace(2200, 9000)
f = fanning_friction_factor(Re)

plt.plot(Re, f)
plt.xlabel('Re')
plt.ylabel('fanning friction factor')
# You can see why we use 0.01 as an initial guess for solving for the
# Fanning friction factor; it falls in the middle of ranges possible
# for these Re numbers.
plt.savefig('images/pipe-diameter-1.png')

def objective(D):
    v = Q / (np.pi * D**2 / 4)
    Re = D * v * rho(T) / mu(T)

    fF = fanning_friction_factor(Re)

    return deltaP - 2 * fF * rho(T) * deltaL * v**2 / D
    
D, = fsolve(objective, 0.04)

print('The minimum pipe diameter is {0} m\n'.format(D))

The minimum pipe diameter is 0.0389653369531 m

Any pipe diameter smaller than that value will result in a larger pressure drop at the same volumetric flow rate, or a smaller volumetric flowrate at the same pressure drop. Either way, it will not meet the design specification.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source