In [1]:
import numpy
from matplotlib import pyplot
%matplotlib inline
size = 7
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}} $$and
$$ \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$.
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)$.
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.
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$.
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])
Out[4]:
With some trial and error we see that $H(0)\approx0.33$ satifies our far field condition on $G$.
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.
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({z})={G}".format(z=delB,G=G[i-1])
Out[5]:
$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!
Find where the 99.9% thickness is by changing $i$ in the code above.
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)]
Out[6]:
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$
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])
print(row_format.format(name,delf,*delta12(f,eta[1]*delf)))
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))
print(head_format.format("del","del_1","del_2","H"))
delL = delAndPrint(eta,"Linear")
delK = delAndPrint(karman,"Karman")
delP = delAndPrint(pohlF,"Pohlhausen")
print(row_format.format("Blasius",delB,*delta12(G,dz)))
Other than the somewhat arbitrary $\delta$ value, these compare very well with the Blasius solution.
Finally, lets compare visually:
In [10]:
pyplot.figure(figsize=(size,size))
pyplot.axis([0,1,0,7])
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')
pyplot.show()
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)
css_styling()
Out[11]: