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...
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
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)
Using the Pohlhausen profile, the various factors in the momentum integral equation are defined as
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]:
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)$!
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.
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'$$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
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)}} $$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]:
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:
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.
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)
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]:
The circle solution is completely different due to the external pressure gradients.
This is in good agreement with Hoerner's Fluid-Dynamic Drag, which states that theoretical laminar separation occurs at ~$110^o$ for this case.
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$.
This is not true of a turbulent boundary layer.
How can you compute the total friction drag coefficient $C_F$ on the circle?
Hint: numpy.trapz
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]: