Content under Creative Commons Attriubtion license CC-BY 4.0, (c)2014 Ian A. Carr, for Prof. L. A. Barba's MAE 6286 course at The George Washington University.

Shoot First, Ask Questions Later

First, congratulations on making it to the 5th module of #numericalmooc! It's been quite a journey. There are two objectives in this notebook. One is to introduce a useful way of solving boundary value problems (BVP) called the shooting method. The second is to start applying our knowledge of numerics to the complex world of fluid motion. Let's get started!

The Blasius Boundary Layer

Computational fluid dynamics is a very complex field. We are going to get started with a relatively simple problem. As the title suggests, it's called the Blasius Boundary Layer.

In case some of you are a little rusty in your fluid mechanics knowledge, you can check out the wikipedia page on boundary layers to get an idea. Put simply, a boundary layer is a thin layer of fluid just above a surface. The effect of the fluid's viscosity is greatest within this thin layer. Boundary layers can get very complex depending on the flow conditions and geometry of the surface, but today we are going to look at a very simple one.

Figure 1 shows fluid moving over a flat plate which continues infinitely in the positive x direction. The fluid moves with a unidirectional velocity we're calling $u_{\infty}$. A student named Heinrich Blasius came up with an interesting solution to this problem in 1908, and that's what we are going to explore today.

Figure 1. Sketch of a boundary layer over a flat plate.

This flow field can be described with a single equation, a 3rd order, nonlinear ODE, but a single equation no less. To understand how this boils down to a single equation we have to go through a bit of derivation first. For those of you unfamiliar with fluids this is all the more critical in understanding this problem.

For this derivation we assume the flow is two-dimensional, steady, incompressible, and has zero pressure gradient. This means we can get rid of a lot of terms in the continuity equation and the equation of fluid motion (Navier-Stokes). Specifically the continuity equation reduces to

\begin{equation} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0, \end{equation}

and the equation of motion reduces to

\begin{equation} u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu. \frac{\partial^2 u}{\partial y^2}. \end{equation}

With some boundary conditions due to the no-slip condition at the surface of the plate

\begin{equation} u = v = 0 \; \textrm{at} \; y = 0,\\ u = U, \; \frac{\partial u}{\partial y} = 0\; \textrm{at} \; y = \infty. \end{equation}

From here Blasius reasoned the velocity profile should look similar for all values of x when compared to a non-dimentional distance from the wall. This means when the non-dimensional velocity is plotted against non-dimentional distance, the profile is the same for all scales. The obvious choice for non-dimentionalizing the distance from the wall is the boundary layer thickness, $\delta$ in Figure 1. Blasius decided to call this non-dimentional distance $\eta$, termed the similarity variable.

Just from these assumptions we can see the solution would take the form

\begin{equation} \frac{u}{U} = g(\eta) \; \textrm{where} \; \eta \propto \frac{y}{\delta}. \end{equation}

Inspired by the Stokes solution, Blasius defined $\eta$ as

\begin{equation} \eta = y \sqrt{\frac{U}{\nu x}}. \end{equation}

You still with me? Good. Now we have to introduce what's called a stream function. If you're unfamiliar with this concept it's a set of equations which describe the stream lines in a given flow. In this case the stream function, $\Psi$, is

\begin{equation} u = \frac{\partial \Psi}{\partial y} \; \textrm{and} \; v = -\frac{\partial \Psi}{\partial x} \end{equation}

As said, we are working with non-dimentional values due to the similarity assumption, so we want to make this stream function dimensionless. To do this we put $u$ and $v$ into the equation of motion, which gives

\begin{equation} f(\eta) = \frac{\Psi}{\sqrt{\nu x U}}. \end{equation}

Now that we have the stream function, $\Psi$, and the similarity variable, $\eta$, defined we can pull out the velocity components. Sparing some algebra, the velocity components are given by

\begin{equation} u = U \frac{df}{d \eta}, \end{equation}\begin{equation} v = \frac{1}{2} \sqrt{\frac{\nu U}{x}} \left[ \eta \frac{df}{d \eta} -f \right]. \end{equation}

Now we want to fit these velocity components into the equation of motion and to do that we need to differentiate them. These come out to

\begin{equation} \frac{\partial u}{\partial x} = -\frac{U}{2x} \eta \frac{d^2 f}{d \eta^2}, \end{equation}\begin{equation} \frac{\partial u}{\partial y} = U \sqrt{U/\nu x} \frac{d^2 f}{d \eta^2}, \end{equation}

and for that second derivative,

\begin{equation} \frac{\partial^2 u}{\partial y^2} = \frac{U^2}{\nu x} \frac{d^3 f}{d \eta^3}. \end{equation}

Finally, we plug those into the equation of motion to find

\begin{equation} 2 \frac{d^3 f}{d \eta^3} + f \frac{d^2 f}{d \eta^2} = 0. \end{equation}

And that's it! That's the Blasius Boundary Layer Solution, the single 3rd order non-linear ODE we need to solve. Not too bad right? We need to make sure we have the boundary conditions in terms of $f$ and $\eta$.

\begin{equation} f = \frac{df}{d \eta} = 0 \; \textrm{at} \; \eta = 0,\\ \frac{df}{d \eta} \rightarrow 1 \; \textrm{at} \; \eta \rightarrow \infty. \end{equation}

Now we need to talk about how we will go about solving this equation.

The Shooting Method

The shooting method provides us with the simplest way to solve the Blasius Boundary Layer. So to get the idea, suppose we have a (BVP) with Dirichlet boundary conditions

\begin{equation} y'' = f(x,y,y'), \; y(a) = \alpha, \; y(b) = \beta . \end{equation}

We then rearrange the BVP to the form

\begin{equation} y' = z, \\ z' = f(x,y,z), \\ y(a) = \alpha, \\ y(b) = \beta. \end{equation}

We can then turn the BVP above into an Initial Value Problem (IVP) by replacing the boundary condition $y(b)$ with

\begin{equation} z(a) = \theta, \end{equation}

where $\theta$ is some number. We can then solve the IVP with the method of our choosing (Euler, Runge-Kutta, etc.) to obtain the solution at $y(b)$. If, when we solve the IVP, $y(b) = \beta$, we've done it! We've solved the BVP! Mostly likely, though, $y(b) \neq \beta$ on the first try and we will have to try other values of $\theta$ before we find the right one.

So what would it take to apply the shooting method to the Blasius boundary layer?

Converting to an IVP

As said, the first step in implimenting the shooting method is turning the BVP into an IVP. To reiterate, we want to rewrite Blasius Boundary Layer Solution in such a way that it only involves the first derivative and can therefore be solved as an IVP. In this case the three equations are

\begin{equation} \frac{df}{d \eta} = f_1 \; \textrm{with } f(0) = 0, \end{equation}\begin{equation} \frac{df_1}{d \eta} = f_2 \; \textrm{with } f_1(0) = 0,\: f_1(\infty)=1, \end{equation}\begin{equation} \frac{df_2}{d \eta} = \frac{1}{2} f_2*f \; \textrm{with } f_2(0) = ?. \end{equation}

The term we have to 'shoot' for is $f_2(0)$, and we can use the knowledge that $f_1(\infty) = 1$ to do so.

Euler's Method

So now that we have the problem laid out, we need to pick a way of solving the ODE. Let's use the simplest one for now, Euler's method. You'll remember we worked with this method way back in the first module on phugoid motion. Applying Euler's method to the IVP equations gives

\begin{equation} f^{i+1} = f^i + f^i *d \eta \end{equation}\begin{equation} f^{i+1}_1 = f_1^i + f_2^i * d \eta \end{equation}\begin{equation} f^{i+1}_2 = f_2^i - \frac{1}{2} * f^i * f_2^i*d \eta \end{equation}

And to make sure, the subscripts indicate the level of separation when turning the BVP into an IVP and the superscripts indicate the iterative scheme step.

Finally, let's start coding this up.

Now, as said, we want to guess the initial value of $f_2$ such that our $f_1(\infty) = 1$ boundary condition is met. We could do this manually b changing the first value, which would look something like this.

But obviously, this is too much work and not what we're going to do. Instead let's just choose an arbitrarily high value for $\eta$, say 5. If we break it up into sufficentialy small steps we should be able to run through them until we find the one which matches our boundary condition. We then solve with each initial guess for $f_2(0)$ using Euler's method and print the final values to see how we did.

Also, let's take a look at how many guesses it took to get to the right value.

The last value of $f_1$, which is representative of it approaching $\infty$, is pretty close to 1, which is what we expect. The correct initial value for $f_2$ was about 0.33, which agrees with the literature. Also, let's make a mental note of the number of guesses it took to find the right one.

Now, let's look at how the non-dimentionalized velocity in the x direction changes as we move away from the wall.

Looks pretty good right? If you go to the Blasius Boundary Layer wikipedia page you can see our solution lines up pretty well with their figure.

We can do the same with the y-component of velocity.

Again, this figure matches the known solution.

Let's get back to the 66 guesses it took to get to the right initial value of $f_2$. We got to the solution by working our way up from 0 in small increments. That's not exactly a very clever way of doing things. Let's look into another method that might reduce the number of guesses we need to take.

The Newton-Raphson Method

We want to improve the way we "aim" before we use the shooting method. An intuative way of doing this would be to compare each shot to the target value and adjust accordingly. This is exactly what the Newton-Raphson method does.

Given a reasonable guess, the Newton-Raphson Method uses the derivative of the function to determine the next guess. This method can then be used iteratively to narrow in on the correct answer using the following equation.

\begin{equation} x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \end{equation}

As you can see this equation requires we know the derivative, $f'(x_n)$, of the function we are working with. In the case of the Blasius Boundary Layer equation, we need to take the derivative of the 3 IVPs determined earlier. This results in

\begin{equation} \frac{df_3}{d \eta} = f_4 \; \textrm{with} \; f_3(0) = 0, \end{equation}\begin{equation} \frac{df_4}{d \eta} = f_5, \; \textrm{with } f_4(0) = 0, \end{equation}\begin{equation} \frac{df_5}{d \eta} = \frac{1}{2}(f_3*f_2 + f*f_5)\; \textrm{with}\; f_5(0) = 1. \end{equation}

As before, we discretize these equations and apply Euler's method to form

\begin{equation} f_3^{i+1} = f_3^i + f_3^i *d \eta, \end{equation}\begin{equation} f_4^{i+1} = f_4^i + f_5^i * d \eta, \end{equation}\begin{equation} f_5^{i+1} = f_5^i - \frac{1}{2}(f_3^i* f_2^i + f^i * f_5^i) * d \eta . \end{equation}

We then use these equations in the Newton-Raphson method to narrow in on the result. The adaptation of the general Newton-Raphson equation for this problem is

\begin{equation} f_2(0) \leftarrow f_2(0) - \frac{f_2(\infty)-1}{f_5(\infty)}. \end{equation}

Okay, so now that we have all that on the table, let's code it up and see how we've improved.

The speed of the method depends on the quality of the initial guess. With most applications it's possible to make an educated guess for this value. Let's start with 0.2 for now.

Now let's loop through both the original IVPs and the derivatives. Because this method doesn't simply grow the way the simple shooting method did, we have to set up a condition for convergence. In the code we can just call this the error bound. We will set it to some arbitrarily small value.

Look at that!

We converge to the same values for $f_1(\infty)$ and for the correct initial condition for $f_2$ but with an order of magnitude fewer guesses!

As said, the speed of the method depends of the quality of the initial guess. Try changing the value of the initial guess to see how it effects the number of required guesses before convergence.

In the coming notebooks we will look at other interesting applications of our knowledge to fluids and other engineering problems.

Dig Deeper

We've shown the Newton-Raphson method greatly reduces the number of guesses required to converge on the answer, but the computation has a few more steps. Use the timeit functionality to see how the computational time compares.

Is it faster or slower? With a more complex problem do you think taking the time to exploy the Newton-Raphson Method would be advantageous?

References

• Pozrikidis, C. Introduction to Theoretical and Computational Fluid Dynamics. New York: Oxford UP, 1997. Print.

• Ben-Israel, Adi. "A Newton-Raphson Method for the Solution of Systems of Equations." Journal of Mathematical Analysis and Applications 15.2 (1966): 243-52. Web.

• Oteh, Uche. Mechanics of Fluids. Bloomington, IN: AuthorHouse, 2008. Print.