<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  Change $\lambda$ below to see its effect on the profile shape.
Step2:  Quiz 1
What value of $\lambda$ denotes separated flow?
$\lambda$<-12
$\lambda$=0
$\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.
Step3:  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?
$u_e = 0$
$u_e' = 0$
$\delta/x << 1$
$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
Step4:  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?
$g_1 = 0$
$g_1 = Re_\delta$
$ \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.
Step5:  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)} $$
Step6:  Lets plot the functions of  $\lambda$ to get a feel for how the boundary layer will develop.
Step7:  Quiz 4
What will happen if $\lambda>\lambda_0$?
Flat plate boundary layer flow.
The boundary layer will shrink.
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
Step8:  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$
Step9:  Looks good, only 1% error. 
Bonus
Step10:  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$
Step11:  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
Step12:  Let plot the boundary layer thickness on the circle compared to the exact solution for a laminar flat plate boundary layer from Blausius.
Step13:  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$?
Increasing $\nu$ delays separation
Decreasing $\nu$ delays separation
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
Step14:  Please ignore the cell below. It just loads our style for the notebooks.
<ASSISTANT_TASK:>
Python Code:
import numpy
def pohlF(eta): return 2*eta-2*eta**3+eta**4
def pohlG(eta): return eta/6*(1-eta)**3
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')
Explanation: 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
End of explanation
pohlPlot(lam=7)
Explanation: Change $\lambda$ below to see its effect on the profile shape.
End of explanation
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.
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')
Explanation: Quiz 1
What value of $\lambda$ denotes separated flow?
$\lambda$<-12
$\lambda$=0
$\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.
End of explanation
def g_1(lam): return df_0(lam)-lam*(disp_ratio(lam)+2*mom_ratio(lam))
Explanation: 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?
$u_e = 0$
$u_e' = 0$
$\delta/x << 1$
$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]$$
End of explanation
from scipy.optimize import bisect
lam0 = bisect(g_1,-12,12)         # use bisect method to find root between -12...12
print 'lambda_0 = ',lam0
Explanation: 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?
$g_1 = 0$
$g_1 = Re_\delta$
$ \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.
End of explanation
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'
Explanation: 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)} $$
End of explanation
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)
Explanation: Lets plot the functions of  $\lambda$ to get a feel for how the boundary layer will develop.
End of explanation
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
Explanation: Quiz 4
What will happen if $\lambda>\lambda_0$?
Flat plate boundary layer flow.
The boundary layer will shrink.
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?
$\psi_{i+1}$ is unknown
$g$ is unknown
$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:
End of explanation
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]))
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)
Explanation: 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$
End of explanation
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
Explanation: 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
End of explanation
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
Explanation: 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$
End of explanation
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!
Explanation: 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:
End of explanation
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)
Explanation: Let plot the boundary layer thickness on the circle compared to the exact solution for a laminar flat plate boundary layer from Blausius.
End of explanation
# your code here
Explanation: 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$?
Increasing $\nu$ delays separation
Decreasing $\nu$ delays separation
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? 
Use the flat plate estimate, I'm sure that will be fine...
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}$.
End of explanation
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom.css', 'r').read()
    return HTML(styles)
css_styling()
Explanation: Please ignore the cell below. It just loads our style for the notebooks.
End of explanation
<END_TASK>