Boundary value problems

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

d = 0.1
mu = 1
deltaPdeltax = -100

The function defining the BVP has to return an array that has a row for each equation, and a column for each value in the grid.

def bvp(y, U):
    u1, u2 = U
    du1dy = u2
    du2dy = np.ones(y.shape) / mu * deltaPdeltax
    return [du1dy, du2dy]

The boundary condition function will get the whole numeric solution at each boundary. We want \(u1(a) = 0\) and \(u1(b)=0\).

def bc(Ua, Ub):
    u1a, u2a = Ua
    u1b, u2b = Ub
    return [u1a, u1b]

Next, we need an initial guess for u1 and u2 on a grid of points. You have to make some decisions here. You need a guess that is reasonably close, but not hard to construct. Here, we anticipate a solution that looks parabolic, and that goes through the points: (0, 0), (d, 0), and some point at (d / 2, ?), where ? represents the point of maximum velocity in middle. We can easily get this polynomial with np.polyfit.

?np.polyfit

We don't know what the maximum velocity is, so we make a guess, say 0.5. Then, we get the parameters, and apply them to an array of y values.

pars = np.polyfit([0, d / 2, d],  # x-points
                  [0, 0.5, 0],    # y-points
                  2)              # polynomial order
pars

array([-2.00000000e+02,  2.00000000e+01, -4.48803257e-16])

Now, we can define a Y grid and define the guess for the first U1.

Y = np.linspace(0, d)

U1 = np.polyval(pars, Y)

We also need a guess for U2, and in this case we know that \(u2 = u1'\), so we just use that.

U2 = np.gradient(U1, Y, edge_order=2)

U = np.array([U1, U2])
print(U.shape)

(2, 50)

You should always visualize the guess to make sure it does what you want. It is hard to make these!

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(Y, U[0], label='u1')
plt.gca().tick_params('y', colors='b')
plt.ylabel('u1')

plt.twinx()
plt.plot(Y, U[1], 'r', label='u2')
plt.gca().tick_params('y', colors='r')
plt.ylabel('u2')
plt.legend()

<Figure size 432x288 with 2 Axes>

Now, we are ready to solve the BVP.

from scipy.integrate import solve_bvp

sol = solve_bvp(bvp, bc, Y, U)
print(sol.message)
plt.plot(sol.x, sol.y[0])
plt.xlabel('y')
plt.ylabel('U')

The algorithm converged to the desired accuracy.

Text(0, 0.5, 'U')

<Figure size 432x288 with 1 Axes>

exercise Try using different guesses, e.g. lines, or triangle shapes, etc. What else looks like this shape? Half a cycle of a sin wave? A semi-circle?

exercise How can you show this is a solution to the BVP?

Another typical boundary value problem in chemical engineering is the concentration profile inside a catalyst particle. Here is the dimensionless equation for a second order reaction in a slab. Note here we have a boundary condition on the derivative at the origin. This kind of condition means either there is no flux at this position, or that the slab is symmetric about this position.

\(\frac{d^2c}{dx^2} = \Phi^2 c^2\)

with \(c'(0)\) = 0 and \(c(1) = 1\)

We again convert this to a system of first order differential equations like this:

Let c1 = c, c1' = c', and c2 = c1', so c2' = c1'' = c''

Then we have:

\(c1' = c2\)

\(c2' = \Phi^2 c1^2\)

with boundary conditions \(c1'(0) = 0\) and \(c1(0) = 1\).

We begin with the required functions:

Phi = 50

def bvp(x, C):
    c1, c2 = C
    dc1dx = c2
    dc2dx = Phi**2 * c1**2
    return [dc1dx, dc2dx]

def bc(Ca, Cb):
    c1a, c2a = Ca
    c1b, c2b = Cb

    # Now, evaluate the derivatives at the first boundary condition
    c1prime, c2prime = bvp(0, [c1a, c2a])
    return [c1prime,  # will all equal zero
            c1b - 1]  # c1(b) = 1

We need an initial guess. We make a naive one, that \(c(x) = 1\) in the slab, i.e. there is no reaction. As usual, we visualize the guess to be sure it does what we intended.

X = np.linspace(0, 1)

C1 = np.ones(X.shape)
C2 = np.gradient(C1, X)

plt.plot(X, C1)

[<matplotlib.lines.Line2D at 0x181eab7110>]

<Figure size 432x288 with 1 Axes>

Now we solve the system.

C = [C1, C2]
sol = solve_bvp(bvp, bc, X, C)
sol.message

'The algorithm converged to the desired accuracy.'

plt.plot(sol.x, sol.y[0])
plt.xlabel('x')
plt.ylabel('C')
plt.xlim([0, 1])
plt.ylim([0, 1])

(0, 1)

<Figure size 432x288 with 1 Axes>

You can see the solution looks nothing like our initial guess. In this case, a high thiele modulus means most of the reaction happens near the catalyst surface, and the interior of the slab has hardly any reactant in it. This solution is consistent with that.

The effectiveness factor for this system is defined by:

\(E = \int_0^1 c^2 dx\)

We can estimate this with the trapezoid or Simpson's method (remember that the solution is a vector of numbers).

c = sol.y[0]
print(np.trapz(c**2, sol.x))

from scipy.integrate import simps
print(simps(c**2, sol.x))

0.016528962860989847 0.01633462094482091

Or, we can use the dense_output of the solution with quad.

from scipy.integrate import quad

def integrand(x):
    c1, c2 = sol.sol(x)
    return c1**2

quad(integrand, 0, 1)

(0.01632998588357348, 9.76148224805893e-09)

excercise Repeat this example for different values of Φ.

exercise Try different kinds of guesses. Think of a guess that has the properties of the boundary conditions, e.g. c'(0) = 0, and c(1) = 1.

exercise Evaluate the quality of the solution based on the equations.