Boundary Layer Solver

This notebook will develop a numerical method for solving the boundary layer momentum integral equation using Pohlhausen velocity profiles.

Momentum integral equation

In the boundary layer portion of the course we derived the governing equations for a boundary layer using the concept of a velocity profile

$$u = u_e(x) f(\eta), \quad \eta=\frac y{\delta(x)}$$

where $u_e$ is the local free stream velocity and $\delta$ is the boundary layer thickness. Note that $x$ is the distance along the wall from the leading edge and $y$ is the distance from the wall.



This enables the development of the momentum integral equation

$$ \frac 12 c_f = \frac{u_e'}{u_e}(\delta_1+2\delta_2)+\delta_2' $$

which balances the local wall friction with the change in the boundary layer profile. The tick mark indicates a derivative, ie $u_e'=\frac{du_e}{dx}$.

The goal is to use the momentum equation to determine how the boundary layer develops, predicting the friction drag and the point of separation.

The velocity $u_e$ (and $u_e'$) is considered to be prescribed by the potential flow solution, but there are still too many unknowns. We need to choose a profile to develop this further...

Pohlhausen profile

The Pohlhausen profile is used to describe a laminar velocity profile exposed to external pressure gradients. The profile is defined as

$$ \frac u {u_e} = f(\eta) = P_F(\eta)+\lambda P_G(\eta) $$

where $\lambda$ is the shape factor, given by $$ \lambda = \frac {\delta^2}\nu u_e'$$

and the profile shapes are defined by

  • $ P_F = 2\eta-2\eta^3+\eta^4 $ is the flat plate profile
  • $ P_G = \frac\eta 6 (1-\eta)^3 $ is the modification for pressure gradients

These can be easly defined using a set of python functions


In [1]:
import numpy

def pohlF(eta): return 2*eta-2*eta**3+eta**4
def pohlG(eta): return eta/6*(1-eta)**3

In [2]:
from matplotlib import pyplot
%matplotlib inline

def pohlPlot(lam):
    pyplot.xlabel(r'$u/u_e$', fontsize=16)
    pyplot.axis([-0.1,1.1,0,1])
    pyplot.ylabel(r'$y/\delta$', fontsize=16)
    eta = numpy.linspace(0.0,1.0,100)
    pyplot.plot(pohlF(eta),eta, lw=1, c='k', label=r'$P_F$')
    pyplot.plot(pohlF(eta)+lam*pohlG(eta),eta, lw=2, c='g', label=r'$P_F+\lambda P_G$')
    pyplot.legend(loc='upper left')

Change $\lambda$ below to see its effect on the profile shape.


In [3]:
pohlPlot(lam=7)


Quiz 1

What value of $\lambda$ denotes separated flow?

  1. $\lambda$<-12
  2. $\lambda$=0
  3. $\lambda$>12

Using the Pohlhausen profile, the various factors in the momentum integral equation are defined as

  • $\frac{\delta_1}\delta = \int_0^1 (1-f) d\eta = \frac3{10}-\lambda\frac1{120}$
  • $\frac{\delta_2}\delta = \int_0^1 f(1-f) d\eta = \frac{37}{315}-\lambda\frac1{945}-\lambda^2\frac1{9072}$
  • $\frac 12 c_f Re_\delta =f'(0)= 2+\lambda\frac1{6}$

where $Re_\delta = \frac{u_e\delta}\nu$ is the local boundary layer Reynolds number.


In [4]:
def disp_ratio(lam): return 3./10.-lam/120.
def mom_ratio(lam): return 37./315.-lam/945.-lam**2/9072.
def df_0(lam): return 2+lam/6.

In [5]:
pyplot.xlabel(r'$\lambda$', fontsize=16)
lam = numpy.linspace(-12,12,100)
pyplot.plot(lam,disp_ratio(lam), lw=2, label=r'$\delta_1/\delta$')
pyplot.plot(lam,mom_ratio(lam), lw=2, label=r'$\delta_2/\delta$')
pyplot.plot(lam,df_0(lam)/10., lw=2, label=r'$c_f Re_\delta/20$')
pyplot.legend(loc='upper right')


Out[5]:
<matplotlib.legend.Legend at 0x1063c2950>

Note that these are all polynomial functions of $\lambda$. Since $u_e$ is given by potential flow and $\lambda = \frac {\delta^2}\nu u_e'$, the only unknown in the momentum equation is now $\delta(x)$!

Stagnation point condition

Now we need to write the momentum equation in terms of $\delta$ (and $\lambda$) and solve. This equation needs to be valid from the leading edge all the way to the point of separation.

For any body with finite thickness the boundary layer will begin at the stagnation point at the front of the body. However, describing the boundary layer at a stagnation point is somewhat tricky.

Quiz 2

Which relationships are true at a stagnation point?

  1. $u_e = 0$
  2. $u_e' = 0$
  3. $\delta/x << 1$
  4. $c_f$ is singular

That's no good - the momentum equation will be singular at the leading edge. We can avoid this problem by multiplying the whole equation by $Re_\delta$, leading to:

$$ \frac 12 c_f Re_\delta = \frac\delta\nu u_e' [\delta_1+2\delta_2]+Re_\delta \delta_2'$$

The first term on the RHS can be simplified by dividing the brackets by $\delta$ and mutiplying by $\delta$ outside out to produce the definition of $\lambda$. This lets us group the terms only dependant on $\lambda$ together to define

$$ g_1(\lambda) = \frac 12 c_f Re_\delta - \lambda \left[\frac{\delta_1}{\delta}+2\frac{\delta_2}\delta\right]$$

In [6]:
def g_1(lam): return df_0(lam)-lam*(disp_ratio(lam)+2*mom_ratio(lam))

Using this definition, the momentum equation is

$$ g_1(\lambda) = Re_\delta \delta_2'$$

Quiz 3

The equation above further simplifies at the stagnation point. Which is correct?

  1. $g_1 = 0$
  2. $g_1 = Re_\delta$
  3. $ \frac 12 c_f = 0$

Solving this equations will determine our initial condition $\lambda_0$. Using my vast google skills I found the bisect function in scipy.optimize which will solve for the root.


In [7]:
from scipy.optimize import bisect
lam0 = bisect(g_1,-12,12)         # use bisect method to find root between -12...12
print 'lambda_0 = ',lam0


lambda_0 =  7.05232310118

With the value of $\lambda_0$ determined, the initial condition $\delta_0$ is simply

$$ \delta_0 = \sqrt{\frac{\nu \lambda_0}{u_e'(x_0)}} $$

Pohlhausen momentum equation

The only thing left to do is write $\delta_2'$ in terms of $\delta'$. Using $F=\frac{\delta_2}\delta$ we have

$$ \delta_2' = \frac{d}{dx}(F\delta) $$

From the line plot above, we see that $F$ is nearly unchanged across the whole range of $\lambda$, so we will treat it as a constant. Therefore the complete Pohlhausen momentum equation is

$$ g_1 = Re_\delta F \delta'$$

Isolating the derivative, we have

$$ \delta'= \frac{g_1(\lambda)}{Re_\delta F(\lambda)} $$

In [8]:
def ddx_delta(Re_d,lam):
    if Re_d==0: return 0                     # Stagnation point condition
    return g_1(lam)/mom_ratio(lam)/Re_d      # delta'

Lets plot the functions of $\lambda$ to get a feel for how the boundary layer will develop.


In [9]:
pyplot.xlabel(r'$\lambda$', fontsize=16)
pyplot.ylabel(r'$g_1/F$', fontsize=16)
pyplot.plot(lam,ddx_delta(1,lam), lw=2)
pyplot.scatter(lam0,0, s=100, c='r')
pyplot.text(lam0,3, r'$\lambda_0$',fontsize=15)


Out[9]:
<matplotlib.text.Text at 0x107363710>
Quiz 4

What will happen if $\lambda>\lambda_0$?

  1. Flat plate boundary layer flow.
  2. The boundary layer will shrink.
  3. The Pohlausen equation will be singular.

Ordinary differential equations

The momentum equation above is an ordinary differentilal equation (ODE), having the form

$$ \psi' = g(\psi(x),x) $$

where the derivative is only a function of the variable $\psi$ and one indepentant variable $x$. All ODEs have an important feature in common:

Mathematics fundamental: ODEs
Systems' whose evolution depends only on their current state

This makes them easier to solve. If we integrate the ODE from $x_0$ to $x_1$ we have

$$ \psi(x_1) = \psi_1= \psi_0+\int_{x_0}^{x_1} g(\psi(x),x) dx $$

which means all we need to solve for $\psi_1$ is the initial condition $\psi_0$ and an estimate the RHS integral. And once we have $\psi_1$, we can get $\psi_2$, etc. In general we have

$$ \psi_{i+1}= \psi_i+\int_{x_i}^{x_{i+1}} g(\psi(x),x) dx \quad i=0,\ldots, N-1$$

This means the ODE can be solved by marching from $x=0$ to $x=L$. Compare this to the vortex panel method and its linear system of equations that needed to be solved simultaneously using matrices... This is easy.

Numerical integration

You've seen numerical ways to determine the area under a curve before, like the trapezioal rule

$$ \int_{x_i}^{x_{i+1}} f(x) dx \approx \frac12[f(x_1)+f(x_{i+1})] \Delta x$$

where $\Delta x=x_{i+1}-x_1$

Quiz 5

What is the important difference between the integral above and the ODE integral?

  1. $\psi_{i+1}$ is unknown
  2. $g$ is unknown
  3. $g$ is nonlinear

This means we have to split the numerical method into two steps. First we estimate the integral as $g(\psi_i,x_i)\Delta x$. This lets us predict an estimate of $\psi_{i+1}$

$$ \tilde\psi_{i+1}= \psi_i+ g(\psi_i,x_i)\Delta x $$

However, this one-sided estimate of the integral is very rough. In the next step we correct the prediction using the trapeziodal rule

$$ \psi_{i+1}= \psi_i+ \frac12[g(\psi_i,x_i)+g(\tilde\psi_{i+1},x_{i+1})]\Delta x$$

This is often called the predictor/corrector method, or Heun's method.

Let's code it up:


In [10]:
def heun(g,psi_i,i,dx,*args):
    g_i = g(psi_i,i,*args)                      # integrand at i
    tilde_psi = psi_i+g_i*dx                    # predicted estimate at i+1
    g_i_1 = g(tilde_psi,i+1,*args)              # integrand at i+1
    return psi_i+0.5*(g_i+g_i_1)*dx             # corrected estimate

In this code we've made the integrand g a function of psi and the index i. Note that g_i_1=$g_{i+1}$ and we've passed $i+1$ as the index. We've also left the option for additional arguments to be passed to g as *args which is required for the boundary layer ODE.

Before we get to that, lets test heun using $\psi'=\psi$ with $\psi_0=1$, since we know the solution is $\psi = e^x$


In [11]:
N = 20                                      # number of steps
x = numpy.linspace(0,numpy.pi,N)            # set up x array from 0..pi
psi = numpy.full_like(x,1.)                 # psi array with phi0=1
def g_test(psi,i): return psi               # define derivative function
for i in range(N-1):                        # march!
    psi[i+1] = heun(g_test,psi[i],i,(x[i+1]-x[i]))

In [12]:
pyplot.plot(x,psi)
pyplot.plot(x,numpy.exp(x))
print 'exp(pi) ~ ', psi[N-1],', error = ',1-psi[N-1]/numpy.exp(numpy.pi)


exp(pi) ~  22.8496648117 , error =  0.0125764524705

Looks good, only 1% error.

Bonus: What is the error if we don't do the correction step?

Boundary layer on a circle

Returning to the boundary layer ODE, we first define a function which can be integrated by heun


In [13]:
def g_pohl(delta_i,i,u_e,du_e,nu):
    Re_d = delta_i*u_e[i]/nu            # compute local Reynolds number
    lam = delta_i**2*du_e[i]/nu         # compute local lambda 
    return ddx_delta(Re_d,lam)          # get derivative

where u_e, du_e, and nu are the extra arguments, needed to compute $Re_\delta$ and $\lambda$.

Then we use this function and heun to march from the initial condition $\lambda_0,\delta_0$ along the boundary layer until we reach the point of separation at $\lambda<-12$


In [14]:
def march(x,u_e,du_e,nu):
    delta0 = numpy.sqrt(lam0*nu/du_e[0])                # set delta0
    delta = numpy.full_like(x,delta0)                   # delta array
    lam = numpy.full_like(x,lam0)                       # lambda array
    for i in range(len(x)-1):                           # march!
        delta[i+1] = heun(g_pohl,delta[i],i,x[i+1]-x[i],    # integrate BL using...
                          u_e,du_e,nu)                          # additional arguments
        lam[i+1] = delta[i+1]**2*du_e[i+1]/nu               # compute lambda
        if abs(lam[i+1])>12: break                          # check stop condition
    return delta,lam,i                                  # return with separation index

and we're done!


Let's test it on the flow around a circle.

In this case the boundary layer will march around the circle from $s=0,\ldots,R\pi$. Lets set the parameters $R=1$, $U_\infty=1$ and $Re_R=10^5$, such that $\nu=10^{-5}$. The tangential velocity around a circular cylinder using potential flow is simply

$$u_e = 2\sin(s)$$

Now that we've defined march we can set-up and solve for the boundary layer in just a few lines of code:


In [15]:
nu = 1e-4                                   # viscosity
N = 32                                      # number of steps
s = numpy.linspace(0,numpy.pi,N)            # distance goes from 0..pi
u_e = 2.*numpy.sin(s)                       # velocity
du_e = 2.*numpy.cos(s)                      # gradient
delta,lam,iSep = march(s,u_e,du_e,nu)       # solve!

Let plot the boundary layer thickness on the circle compared to the exact solution for a laminar flat plate boundary layer from Blausius.


In [19]:
pyplot.ylabel(r'$\delta/R$', fontsize=16)
pyplot.xlabel(r'$s/R$', fontsize=16)
pyplot.plot(s[:iSep+1],delta[:iSep+1],lw=2,label='Circle')
pyplot.plot(s,s*5/numpy.sqrt(s/nu),lw=2,label='Flat plate')
pyplot.legend(loc='upper left')
pyplot.scatter(s[iSep],delta[iSep], s=100, c='r')
pyplot.text(s[iSep]+0.1,delta[iSep],'separation between\n'
            +'%.2f' % s[iSep]+'<s<'+'%.2f' % s[iSep+1],fontsize=12)


Out[19]:
<matplotlib.text.Text at 0x107249f10>

The circle solution is completely different due to the external pressure gradients.

  • The boundary layer growth is stunted on the front body
  • $\delta$ increases rapidly as the flow approaches the midbody
  • The flow separates around $1.87$ radians $\approx 107^o$

This is in good agreement with Hoerner's Fluid-Dynamic Drag, which states that theoretical laminar separation occurs at ~$110^o$ for this case.

Quiz 6

How does the separation point depend on $\nu$?

  1. Increasing $\nu$ delays separation
  2. Decreasing $\nu$ delays separation
  3. Chaning $\nu$ has no effect on separation

We know analytically that $\delta$ scales as $\sqrt\nu$ (which you can double check using the code above), and therefore $\lambda=\frac{\delta^2}\nu u_e'$ doesn't depend on $\nu$ at all. Since the separation point is determined by $\lambda$, this is also independant of $\nu$.

Fluids fundamental: Separation Point
The point of laminar separation is independant of $Re$

This is not true of a turbulent boundary layer.

Quiz 7

How can you compute the total friction drag coefficient $C_F$ on the circle?

  1. Use the flat plate estimate, I'm sure that will be fine...
  2. Compute $\tau_w=\frac 12 c_f \rho u_e^2 $ and integrate numerically

Hint: numpy.trapz

Your turn

Determine $C_F=\frac {2F_F}{\rho U_\infty^2 S}$ , where $F_F = \int \tau_w s_x ds$ is the 2D friction drag and $S$ is the 2D surface area, and compare it to the flat plate solution: $1.33 Re^{-1/2}$.


In [17]:
# your code here

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


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


Out[18]:

In [18]: