## Solving the Blasius equation

| categories: bvp | tags: fluids | View Comments

In fluid mechanics the Blasius equation comes up (http://en.wikipedia.org/wiki/Blasius_boundary_layer) to describe the boundary layer that forms near a flat plate with fluid moving by it. The nonlinear differential equation is:

\begin{eqnarray} f''' + \frac{1}{2} f f'' &=& 0 \\ f(0) &=& 0 \\ f'(0) &=& 0 \\ f'(\infty) &=& 1 \end{eqnarray}

This is a nonlinear, boundary value problem. The point of solving this equation is to get the value of $$f''(0)$$ to evaluate the shear stress at the plate.

We have to convert this to a system of first-order differential equations. Let $$f_1 = f$$, $$f_2 = f_1'$$ and $$f_3 = f_2'$$. This leads to:

\begin{eqnarray} f_1' = f_2 \\ f_2' = f_3 \\ f_3' = -\frac{1}{2} f_1 f_3 \\ f_1(0) = 0 \\ f_2(0) = 0 \\ f_2(\infty) = 1 \end{eqnarray}

It is not possible to specify a boundary condition at $$\infty$$ numerically, so we will have to use a large number, and verify it is “large enough”. From the solution, we evaluate the derivatives at $$\eta=0$$, and we have $$f''(0) = f_3(0)$$.

We have to provide initial guesses for f_1, f_2 and f_3. This is the hardest part about this problem. We know that f_1 starts at zero, and is flat there (f'(0)=0), but at large eta, it has a constant slope of one. We will guess a simple line of slope = 1 for f_1. That is correct at large eta, and is zero at η=0. If the slope of the function is constant at large $$\eta$$, then the values of higher derivatives must tend to zero. We choose an exponential decay as a guess.

Finally, we let a solver iteratively find a solution for us, and find the answer we want. The solver is in the pycse module.

import numpy as np
from pycse import bvp

def odefun(F, x):
f1, f2, f3 = F
return [f2,
f3,
-0.5 * f1 * f3]

def bcfun(Fa, Fb):
return [Fa[0],        # f1(0) =  0
Fa[1],        # f2(0) = 0
1.0 - Fb[1]]  # f2(inf) = 1

eta = np.linspace(0, 6, 100)
f1init = eta
f2init = np.exp(-eta)
f3init = np.exp(-eta)

Finit = np.vstack([f1init, f2init, f3init])

sol = bvp(odefun, bcfun, eta, Finit)

print "f''(0) = f_3(0) = {0}".format(sol[2, 0])

import matplotlib.pyplot as plt
plt.plot(eta, sol[0])
plt.xlabel('$\eta$')
plt.ylabel('$f(\eta)$')
plt.savefig('images/blasius.png')

f''(0) = f_3(0) = 0.332491109552


org-mode source

## Another look at nonlinear BVPs

| categories: bvp | tags: | View Comments

Boundary value problems may have more than one solution. Let us consider the BVP:

\begin{eqnarray} y'' + |y| &=& 0 \\ y(0) &=& 0 \\ y(4) &=& -2 \end{eqnarray}

We will see this equation has two answers, depending on your initial guess. We convert this to the following set of coupled equations:

\begin{eqnarray} y_1' &=& y_2 \\ y_2' &=& -|y_1| \\ y_1(0)&=& 0\\ y_1(4) &=& -2 \end{eqnarray}

This BVP is nonlinear because of the absolute value. We will have to guess solutions to get started. We will guess two different solutions, both of which will be constant values. We will use pycse.bvp to solve the equation.

import numpy as np
from pycse import bvp
import matplotlib.pyplot as plt

def odefun(Y, x):
y1, y2 = Y
dy1dx = y2
dy2dx = -np.abs(y1)
return [dy1dx, dy2dx]

def bcfun(Ya, Yb):
y1a, y2a = Ya
y1b, y2b = Yb

return [y1a, -2 - y1b]

x = np.linspace(0, 4, 100)

y1 = 1.0 * np.ones(x.shape)
y2 = 0.0 * np.ones(x.shape)

Yinit = np.vstack([y1, y2])

sol = bvp(odefun, bcfun, x, Yinit)

plt.plot(x, sol[0])

# another initial guess
y1 = -1.0 * np.ones(x.shape)
y2 = 0.0 * np.ones(x.shape)

Yinit = np.vstack([y1, y2])

sol = bvp(odefun, bcfun, x, Yinit)

plt.plot(x, sol[0])
plt.legend(['guess 1', 'guess 2'])
plt.savefig('images/bvp-another-nonlin-1.png')
plt.show()


This example shows that a nonlinear BVP may have different solutions, and which one you get depends on the guess you make for the solution. This is analogous to solving nonlinear algebraic equations (which is what is done in solving this problem!).

org-mode source

## Boundary value problem in heat conduction

| categories: bvp | tags: heat transfer | View Comments

For steady state heat conduction the temperature distribution in one-dimension is governed by the Laplace equation:

$$\nabla^2 T = 0$$

with boundary conditions that at $$T(x=a) = T_A$$ and $$T(x=L) = T_B$$.

The analytical solution is not difficult here: $$T = T_A-\frac{T_A-T_B}{L}x$$, but we will solve this by finite differences.

For this problem, lets consider a slab that is defined by x=0 to x=L, with $$T(x=0) = 100$$, and $$T(x=L) = 200$$. We want to find the function T(x) inside the slab.

We approximate the second derivative by finite differences as

$$f''(x) \approx \frac{f(x-h) - 2 f(x) + f(x+h)}{h^2}$$

Since the second derivative in this case is equal to zero, we have at each discretized node $$0 = T_{i-1} - 2 T_i + T_{i+1}$$. We know the values of $$T_{x=0} = \alpha$$ and $$T_{x=L} = \beta$$.

$A = \left [ \begin{array}{ccccc} % -2 & 1 & 0 & 0 & 0 \\ 1 & -2& 1 & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 1 & -2 \end{array} \right ]$

$x = \left [ \begin{array}{c} T_1 \\ \vdots \\ T_N \end{array} \right ]$

$b = \left [ \begin{array}{c} -T(x=0) \\ 0 \\ \vdots \\ 0 \\ -T(x=L) \end{array} \right]$

These are linear equations in the unknowns $$x$$ that we can easily solve. Here, we evaluate the solution.

import numpy as np

#we use the notation T(x1) = alpha and T(x2) = beta
x1 = 0; alpha = 100
x2 = 5; beta = 200

npoints = 100

# preallocate and shape the b vector and A-matrix
b = np.zeros((npoints, 1));
b[0] = -alpha
b[-1] = -beta

A = np.zeros((npoints, npoints));

#now we populate the A-matrix and b vector elements
for i in range(npoints ):
for j in range(npoints):
if j == i: # the diagonal
A[i,j] = -2
elif j == i - 1: # left of the diagonal
A[i,j] = 1
elif j == i + 1: # right of the diagonal
A[i,j] = 1

# solve the equations A*y = b for Y
Y = np.linalg.solve(A,b)

x = np.linspace(x1, x2, npoints + 2)
y = np.hstack([alpha, Y[:,0], beta])

import matplotlib.pyplot as plt

plt.plot(x, y)

plt.plot(x, alpha + (beta - alpha)/(x2 - x1) * x, 'r--')

plt.xlabel('X')
plt.ylabel('T(X)')
plt.legend(('finite difference', 'analytical soln'), loc='best')
plt.savefig('images/bvp-heat-conduction-1d.png')


org-mode source

## Plane Poiseuille flow - BVP solve by shooting method

| categories: bvp | tags: fluids | View Comments

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$$.

## 1 First guess

We need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0) and see if the solution goes through the u_2(d)=0 boundary value.

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')


Here we have undershot the boundary condition. Let us try a larger guess.

## 2 Second guess

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')


Now we have clearly overshot. Let us now make a function that will iterate for us to find the right value.

## 3 Let fsolve do the work

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')


You can see the agreement is excellent!

This also seems like a useful bit of code to not have to reinvent regularly, so it has been added to pycse as BVP_sh. Here is an example usage.

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()


org-mode source

## Plane poiseuelle flow solved by finite difference

| categories: bvp | tags: fluids | View Comments

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.