In [1]:
import numpy
from matplotlib import pyplot
%matplotlib inline
size = 7

The Blasius boundary layer solution

The Blasius Equation

The Prandtl equations on a flat plate are

$$ u \frac{\partial u}{\partial x}+ v \frac{\partial u}{\partial y} = \nu \frac{\partial^2 u}{\partial y^2} $$$$ \frac{\partial u}{\partial x}+ \frac{\partial v}{\partial y} = 0. $$

We don't have an equation for $v$, so we use the stream function $\Psi$ to define the velocities as

$$ u = \frac{\partial \Psi}{\partial y}, \quad v = -\frac{\partial \Psi}{\partial x} $$

and automatically satisfy the continuity equation.

From scaling arguments we know the boundary layer thickness $\delta$ scales as

$$ \frac\delta x \sim \frac 1 {\sqrt{Re}} $$


$$ \frac u U = g\left(\frac y\delta\right) $$

for a laminar flat plate boundary layer. Therefore we choose

$$ \Psi = \sqrt{U \nu x}\enspace F(z), \quad z=y\sqrt{\frac U{\nu x}} $$

such that

$$ u = UF'(z) $$

where $F'=\tfrac{\partial}{\partial z}F$. Substitution into the momentum equation and simplification gives the Blasius equation

$$ F''' +\frac 12 FF'' = 0 $$

with the boundary conditions $F(0) = F'(0) = 0$ and $F'(\infty)=1$.

Numerical Solution

There is no known solution for this equation, so we will resort to a numerical method.


If this were an initial-value first-order differential equation we could use the classic Runge-Kutta method

In [2]:
def RK4(fz,fp,dz):
    k1 = fp(fz)
    k2 = fp(fz+dz/2*k1)
    k3 = fp(fz+dz/2*k2)
    k4 = fp(fz+dz*k3)
    return fz+dz/6*(k1+2*k2+2*k3+k4)

where fz=$f(z)$, fp=$f'$ is a function defining the first derivative of $f$, and dz=$\Delta z$ is the step size. RK4 returns the value of $f(z+\Delta z)$.

First order-system

In order to apply the Runge-Kutta method to the Blasius equation, we have to replace the third-order ODE with three first-order equations $$ F'=G $$ $$ G'=H $$ $$ H'= F''' = -\frac 12 F H $$

The right hand side of these equations defines the first derivative of the system $f=[F,G,H]$:

In [3]:
def fp(f):
    return numpy.array([f[1],f[2],-0.5*f[0]*f[2]])

and this function can be used in the RK4 function to integrate the system.

Initial conditions

But first, we need to set the initial conditions. The boundary conditions given above are $F(0)=0$, $F'(0)=G(0)=0$ and $F'(\infty)=G(\infty)=1$. However, the Runge-Kutta method requires that we set $F''(0)=H(0)$, which is unknown. The simplest strategy is to guess and adjust $H(0)$ until $G(\infty)=1$.

Challenge question

Change H0 in the code below until you achieve the far-field condition

In [4]:
H0 = 0.332           # Initial Condition
dz = 0.05            # Step size
Z = 10               # Integration limit

n = int(Z/dz)        # Set-up arrays
F = numpy.zeros(n)
G = numpy.zeros(n)
H = numpy.zeros(n)
H[0] = H0

for i in range(n-1): # Integrate from z=0..Z
    fz = numpy.array([F[i],G[i],H[i]])
    fz2 = RK4(fz,fp,dz)
    F[i+1] = fz2[0]
    G[i+1] = fz2[1]
    H[i+1] = fz2[2]

"H(0)={H0:0.3f}, G({Z})={Ginf}".format(H0=H0,Z=Z,Ginf=G[n-1])

'H(0)=0.332, G(10)=0.999884880035'

With some trial and error we see that $H(0)\approx0.33$ satifies our far field condition on $G$.

Challenge question

Is $Z=10$ really close enough to $\infty$? Is $\Delta z=0.05$ small enough? Adjust the step size and integration limit to see how sensitive this solution is.

Blasius profile properties

Boundary layer thickness: $\delta$

Because of our boundary conditions, $G=u/U=1$ only at $z=\infty$, but $G$ is very close to 1 for finite values of $z$

In [5]:
i = 100
delB = (i-1)*dz


$G=u/U=0.99$ at $z=y\sqrt{\frac U{\nu x}}\approx 5$. This is typically used as the Blasius boundary layer thickness - but this is arbitrary!

Challenge question

Find where the 99.9% thickness is by changing $i$ in the code above.

Displacement and momentum thickness

The displacement thickness and momentum thickness values are well defined irregardless of $\delta$.

$$ \delta_1 = \int_0^\infty 1-\frac u U dy $$$$ \delta_2 = \int_0^\infty \frac u U \left(1-\frac u U\right) dy $$

In [6]:
def delta12(f,dz): # Print del1,del2
    del1 = numpy.trapz(1-f,dx=dz)
    del2 = numpy.trapz(f*(1-f),dx=dz)
    return (del1,del2,del1/del2)

["%0.3f" % j for j in delta12(G,dz)]

['1.722', '0.665', '2.590']

And we see $\delta_1 \sqrt{\frac U{\nu x}} \approx 1.7$,$\ $ $\delta_2\sqrt{\frac U{\nu x}} \approx 0.66$, and the shape factor $H = \delta_1/\delta_2 \approx 2.6$

Comparison to polynomial profiles

We next compare the Blasius profiles to the linear, Karman (quadratic), and Pohlhausen (quintic) profiles.

First we compute the boundary layer thickness factor $\delta \sqrt{\frac U{\nu x}} = \sqrt{\frac{2f'(0)}F}$ where $f(\eta)$ is the polynomial defining the profile, $\eta=y/\delta$, and $F=\delta_2/\delta$.

In [7]:
def delta(f,deta):
    fp = (f[1]-f[0])/deta # Finite difference estimate of f'(0)
    F = numpy.trapz(f*(1-f),dx=deta)
    return numpy.sqrt(2*fp/F)

Since $z=\eta\delta\sqrt{\frac U{\nu x}}$, we can then use the delta12 function defined above to compute the other thicknesses in terms of $z$.

In [8]:
row_format ="{:>10}"+"{:>8.3f}"*(4)
def delAndPrint(f,name):
    delf = delta(f,eta[1])
    return delf

Now we need only to define the profiles and compare:

In [9]:
eta = numpy.linspace(0.0,1.0,100)
karman = 2*eta-eta**2
pohlF = 2*eta-2*eta**3+eta**4

head_format =("{:>18}"+"{:>8}"*(3))
delL = delAndPrint(eta,"Linear")
delK = delAndPrint(karman,"Karman")
delP = delAndPrint(pohlF,"Pohlhausen")

               del   del_1   del_2       H
    Linear   3.464   1.732   0.577   3.000
    Karman   5.464   1.821   0.728   2.500
Pohlhausen   5.836   1.751   0.685   2.555
   Blasius   4.950   1.722   0.665   2.590

Other than the somewhat arbitrary $\delta$ value, these compare very well with the Blasius solution.

Finally, lets compare visually:

In [10]:
pyplot.xlabel('u/U', fontsize=16)
pyplot.ylabel('z', fontsize=16)
pyplot.plot(eta,eta*delL, lw=3, label='Linear')
pyplot.plot(karman,eta*delK, lw=3, label='Karman')
pyplot.plot(pohlF,eta*delP, lw=3, color='k', label='Pohlhausen')
pyplot.plot(G,numpy.linspace(0,Z,n), lw=3, label='Blasius')
pyplot.legend(loc='upper left')

We can see that both the second and forth order estimates are close to the Blasius profile, but that neither matches exactly. This indictes that while these profiles match the boundary conditions, they do not satisfy the Prandtl equations throughout the boundary layer.

Please ignore the cell below. It just loads our style for the notebooks.

In [11]:
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom.css', 'r').read()
    return HTML(styles)