## Solving the Blasius equation

Posted March 11, 2013 at 10:44 AM | categories: bvp | tags: fluids | View Comments

Updated March 11, 2013 at 04:34 PM

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

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