Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

PHY454H1S Continuum Mechanics. Lecture 19: Boundary layers and stream functions. The Blasius problem. Taught by Prof. K. Das.

Posted by peeterjoot on April 4, 2012

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Disclaimer.

Peeter’s lecture notes from class. May not be entirely coherent.

Review. Laminar boundary layer theory.

In the boundary layer we found

  • Continuity equation

    \begin{aligned}\frac{\partial u}{\partial x} + \frac{\partial {v}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.1)

  • x momentum equation

    \begin{aligned}u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = - \frac{1}{{\rho}} \frac{\partial {p}}{\partial x} + \nu \frac{\partial^2 u}{\partial y^2}\end{aligned} \hspace{\stretch{1}}(2.2)

  • y momentum equation

    \begin{aligned}\frac{\partial {p}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.3)

In the inviscid region p is a constant in y

\begin{aligned}\frac{\partial {p}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.4)

This will be approximately true in the boundary layer too as illustrated in figure (1).

Figure 1: Illustrating the boundary layer thickness and associated pressure variation.

 

Starting with

\begin{aligned}\frac{\partial {\mathbf{u}}}{\partial {t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = -\boldsymbol{\nabla} \left( \frac{p}{\rho} + \chi \right) + \nu \boldsymbol{\nabla}^2 \mathbf{u},\end{aligned} \hspace{\stretch{1}}(2.5)

we were able to show that inviscid irrotational incompressible flows are governed by Bernoulli’s equation

\begin{aligned}\frac{p}{\rho} + \chi + \frac{1}{{2}} \mathbf{u}^2 = \text{constant}\end{aligned} \hspace{\stretch{1}}(2.6)

In the absence of body forces (or constant potentials), we have

\begin{aligned}-\frac{\partial {}}{\partial x}\frac{p}{\rho} = \frac{\partial {}}{\partial x} \left(\frac{1}{{2}} \mathbf{u}^2 \right) \sim U \frac{\partial {U}}{\partial x}\end{aligned} \hspace{\stretch{1}}(2.7)

so that our boundary layer equations are

\begin{subequations}

\begin{aligned}u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = U \frac{dU}{dx} + \nu \frac{\partial^2 u}{\partial y^2}\end{aligned} \hspace{\stretch{1}}(2.8a)

\begin{aligned}\frac{\partial {p}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.8b)

\begin{aligned}\frac{\partial u}{\partial x} + \frac{\partial {v}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.8c)

\end{subequations}

With boundary conditions

\begin{aligned}U(x, 0) &= 0 \\ U(x, \infty) &= U(x) \\ V(x, 0) &= 0\end{aligned} \hspace{\stretch{1}}(2.9)

Fluid flow over a flat plate (Blasius problem).

Define a similarity variable \eta

\begin{aligned}\eta \sim \frac{y}{\sqrt{\nu t}}\end{aligned} \hspace{\stretch{1}}(3.12)

Suppose we want

\begin{aligned}\eta \sim f(y, x)\end{aligned} \hspace{\stretch{1}}(3.13)

Since we have

\begin{aligned}x \sim U t,\end{aligned} \hspace{\stretch{1}}(3.14)

or

\begin{aligned}t \sim \frac{x}{U}.\end{aligned} \hspace{\stretch{1}}(3.15)

We can make the transformation

\begin{aligned}\eta = \frac{y}{\sqrt{2 \frac{\nu x}{U}}}\end{aligned} \hspace{\stretch{1}}(3.16)

We can introduce stream functions

\begin{aligned}u &= \frac{\partial {\psi}}{\partial y} \\ v &= -\frac{\partial {\psi}}{\partial x}\end{aligned} \hspace{\stretch{1}}(3.17)

We can check that this satisfies the continuity equation since we have

\begin{aligned}\frac{\partial u}{\partial x} + \frac{\partial {v}}{\partial y} &=\frac{\partial^2 \psi}{\partial x \partial y}-\frac{\partial^2 \psi}{\partial y \partial x} \\ &= 0\end{aligned}

Now introduce a similarity variable

\begin{aligned}f(\eta) = \frac{\psi}{\delta U_0} = \frac{\psi}{\sqrt{2 \nu x U_0}}\end{aligned} \hspace{\stretch{1}}(3.19)

Note that we’ve also suddenly assumed that U = U_0 (a constant, which will also kill the U' term in the N-S equation). This isn’t really justified by anything we’ve done so far, but asking about this in class, it was stated that this is a restriction for this formulation of the Blasius problem.

Also note that this last step requires:

\begin{aligned}\delta = \sqrt{ 2 \nu x/U_0 }.\end{aligned} \hspace{\stretch{1}}(3.20)

This at least makes sense dimensionally since we have [\sqrt{\nu x/U}] = \sqrt{ (L^2/T) (L) (T/L)} = L, but where did this definition of \delta come from?

In [] it is mentioned that this is a result of the scaling argument. We did have some scaling arguments that included \delta in the expressions from last lecture, one of which was 2.7

\begin{aligned}\delta \sim \frac{v L}{U},\end{aligned} \hspace{\stretch{1}}(3.21)

but that doesn’t obviously give us 3.20?

Ah. We argued that

\begin{aligned}v \frac{\partial u}{\partial y} \sim \frac{U^2}{L}\end{aligned} \hspace{\stretch{1}}(3.22)

and that the larger of the viscous terms was

\begin{aligned}\nu \frac{\partial^2 u}{\partial y^2} \sim \frac{\nu U}{\delta^2}.\end{aligned} \hspace{\stretch{1}}(3.23)

If we require that these are the same order of magnitude, as argued in section 8.3 of [2], then we find 3.21.

Regardless, given this change of variables, we can apparently compute

\begin{aligned}f''' + f f'' = 0.\end{aligned} \hspace{\stretch{1}}(3.24)

Our boundary conditions are

\begin{aligned}\begin{array}{l l}f = f' = 0 & \quad \mbox{latex \eta = 0$} \\ f’ = 1 & \quad \mbox{\eta = \infty} \\ \end{array}\end{aligned} \hspace{\stretch{1}}(3.25)$

Deriving the equation of motion.

Attempting to derive 3.24 using the definitions above gets a bit messy. It’s messy enough that I mistakenly thought that we couldn’t possible arrive at a differential equation that has a plain old (non-derivative) f in it as in 3.24 above. The algebra involved in taking the derivatives got the better of me. This derivation is treated a different way in [2]. For the purpose of completeness (and because that derivation also leaves out some details), lets do this from start to finish with all the gory details, but following the outline provided in the text.

Instead of pre-determining the form of the similarity variable exactly, we can state it in terms of an unknown and to be determined function of position writing

\begin{aligned}\eta = \frac{y}{g(x)}.\end{aligned} \hspace{\stretch{1}}(3.26)

We still introduce stream our stream functions 3.17, but require that our horizontal velocity component is only a function of our similarity variable

\begin{aligned}u = U h(\eta),\end{aligned} \hspace{\stretch{1}}(3.27)

where h(\eta) is to be determined, and is scaled by our characteristic velocity U. Observe that, as above, we are assuming that U(x) = U, a constant (which also kills off the U U' term in the Navier-Stokes equation.) Given this form of u, we note that

\begin{aligned}u(\eta) = \frac{\partial {\psi}}{\partial y} = \frac{\partial {\psi}}{\partial {\eta}} \frac{\partial {\eta}}{\partial y} = \frac{1}{{g(x)}} \frac{\partial {\psi}}{\partial {\eta}},\end{aligned} \hspace{\stretch{1}}(3.28)

so that

\begin{aligned}\psi = U g(x) \int_0^\eta h(\eta') d\eta' + k(x).\end{aligned} \hspace{\stretch{1}}(3.29)

It’s argued in the text that we also want \psi to be a streamline, so that \psi(\eta = 0) = 0 implying that k(x) = 0. I don’t honestly follow the rational for that, but it’s certainly convenient to set k(x) = 0, so lets do so and see where things go. With

\begin{aligned}f(\eta) = \int_0^\eta h(\eta') d\eta'.\end{aligned} \hspace{\stretch{1}}(3.30)

Observe that f(0) is necessarily zero with this definition. We can now write

\begin{aligned}\psi(\eta, x) = U g(x) f(\eta).\end{aligned} \hspace{\stretch{1}}(3.31)

This is like what we had in class, with the exception that instead of a constant relating \psi and f(\eta) we also have a function of x. That’s exactly what we need so that we can end up with both f and derivatives of f in our Navier-Stokes equation.

Now let’s do the mechanical bits, computing all the derivatives. We can compute v to start with

\begin{aligned}v &= -\frac{\partial {\psi}}{\partial x} \\ &= - U \frac{\partial {}}{\partial x} \left( g(x) f(\eta) \right) \\ &= - U \left( g' f + g f' \frac{\partial {\eta}}{\partial x} \right) \\ &= - U \left( g' f - g f' \frac{y g'}{g^2} \right) \\ &= - U \left( g' f - f' \frac{y g'}{g} \right) \\ &= - U \left( g' f - f' g' \eta \right) \\ &= U g' \left( f' \eta - f \right)\end{aligned}

We had initially u = U h(\eta), but f'(\eta) = h(\eta), so we’ve now got both u and v specified in terms of f and g and their derivatives

\begin{aligned}u &= U f' \\ v &= U g' \left( f' \eta - f \right).\end{aligned} \hspace{\stretch{1}}(3.32)

We’ve got a bunch of the u derivatives that we have to compute

\begin{aligned}\frac{\partial u}{\partial x} &=U \frac{\partial { f' }}{\partial x} \\ &=U \frac{\partial { f' }}{\partial {\eta}} \frac{\partial {\eta}}{\partial x} \\ &=-U f'' \frac{y g'}{g^2} \\ &=-U f'' \eta \frac{g'}{g},\end{aligned}

and

\begin{aligned}\frac{\partial u}{\partial y} &=U \frac{\partial { f' }}{\partial y} \\ &=U f'' \frac{\partial {\eta}}{\partial y} \\ &=U f'' \frac{1}{{g}},\end{aligned}

and

\begin{aligned}\frac{\partial^2 u}{\partial y^2} &=U \frac{1}{{g}} \frac{\partial {f''}}{\partial y} \\ &=U \frac{1}{{g}} f''' \frac{\partial {\eta}}{\partial y} \\ &=U \frac{1}{{g^2}} f'''.\end{aligned}

Our x component of Navier-Stokes now takes the form

\begin{aligned}0 &=u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}-\nu \frac{\partial^2 u}{\partial y^2} \\ &=U f' \left( -U f'' \eta \frac{g'}{g}\right)+U g' \left( f' \eta - f \right)U f'' \frac{1}{{g}}- \nu U \frac{1}{{g^2}} f''' \\ &=\frac{U}{g^2}\left( \not{{-U f' g f'' \eta g' }}+\not{{U g g' f' \eta f''}}-U g g' f f'' - \nu f'''\right)\end{aligned}

or (assuming g \ne 0)

\begin{aligned}f''' + \frac{U g g' }{\nu} f f'' = 0.\end{aligned} \hspace{\stretch{1}}(3.34)

Now, it we wish U g g'/\nu = 1 (to make this equation as easy to solve as possible), we can integrate to find the required form of g(x). This gives

\begin{aligned}\frac{U g^2 }{2 \nu} = x + C.\end{aligned} \hspace{\stretch{1}}(3.35)

It’s argued that we expect

\begin{aligned}\frac{\partial u}{\partial y} = U f'' \frac{1}{{g}},\end{aligned} \hspace{\stretch{1}}(3.36)

to become singular at x = 0, so we should set C = 0. This leaves us with

\begin{subequations}

\begin{aligned}g(x) = \sqrt{\frac{2 x \nu}{U}}\end{aligned} \hspace{\stretch{1}}(3.37a)

\begin{aligned}\eta = \frac{y}{\sqrt{\frac{2 x \nu}{U}}}\end{aligned} \hspace{\stretch{1}}(3.37b)

\begin{aligned}f''' + f f'' = 0\end{aligned} \hspace{\stretch{1}}(3.37c)

\begin{aligned}u = U f'\end{aligned} \hspace{\stretch{1}}(3.37d)

\begin{aligned}v = (\eta f' - f) \frac{\nu}{g},\end{aligned} \hspace{\stretch{1}}(3.37e)

\end{subequations}

and boundary value conditions

\begin{subequations}

\begin{aligned}f(0) = 0\end{aligned} \hspace{\stretch{1}}(3.38a)

\begin{aligned}f'(0) = 0\end{aligned} \hspace{\stretch{1}}(3.38b)

\begin{aligned}f'(\infty) = 1.\end{aligned} \hspace{\stretch{1}}(3.38c)

\end{subequations}

where 3.38b follows from u(0) = 0 and 3.37d, and 3.38c follows from the fact that u tends to U.

Numeric solution.

We can solve this numerically and find solutions that look like figure (2).

Figure 2: Boundary layer solution to flow over plate.

 

This is the Blasius solution to the problem of fluid flow over a flat plate.

Singular perturbation theory.

In the boundary layer analysis we’ve assumed that our inertial term and viscous terms were of the same order of magnitude. Lets examine the validity of this assumption

\begin{aligned}u \frac{\partial u}{\partial x} \sim \nu \frac{\partial^2 u}{\partial y^2}\end{aligned} \hspace{\stretch{1}}(4.39)

or

\begin{aligned}\frac{U}{L} \sim \frac{\nu U}{\delta^2}\end{aligned} \hspace{\stretch{1}}(4.40)

or

\begin{aligned}\delta \sim \frac{\nu L}{U}\end{aligned} \hspace{\stretch{1}}(4.41)

finally

\begin{aligned}\frac{\delta}{L} \sim \frac{\nu}{U L} \sim (\text{Re})^{-1/2}.\end{aligned} \hspace{\stretch{1}}(4.42)

If we have

\begin{aligned}\frac{\delta}{L} \ll 1,\end{aligned} \hspace{\stretch{1}}(4.43)

and

\begin{aligned}\frac{\delta}{L} \ll \frac{1}{{\sqrt{\text{Re}}}}\end{aligned} \hspace{\stretch{1}}(4.44)

then

\begin{aligned}\frac{\delta}{L} \ll 1\end{aligned} \hspace{\stretch{1}}(4.45)

when \text{Re} >> 1.

(this is the whole reason that we were able to do the previous analysis).

Our EOM is

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = -\frac{\boldsymbol{\nabla} p}{\rho } + \nu \boldsymbol{\nabla}^2 \mathbf{u}\end{aligned} \hspace{\stretch{1}}(4.46)

with

\begin{aligned}\mathbf{u} &\rightarrow U \\ p &\rightarrow U^2 \rho\end{aligned} \hspace{\stretch{1}}(4.47)

as x, y \rightarrow L

performing a non-dimensionalization we have

\begin{aligned}(\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = -\boldsymbol{\nabla}' p' + \frac{\nu}{U L} {\boldsymbol{\nabla}'}^2 \mathbf{u}'\end{aligned} \hspace{\stretch{1}}(4.49)

or

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = -\boldsymbol{\nabla} p + \frac{1}{{\text{Re}}} {\boldsymbol{\nabla}}^2 \mathbf{u}\end{aligned} \hspace{\stretch{1}}(4.50)

to force \text{Re} \rightarrow \infty, we can write

\begin{aligned}\frac{1}{{\text{Re}}} = \epsilon.\end{aligned} \hspace{\stretch{1}}(4.51)

so that as \epsilon \rightarrow 0 we have \text{Re} \rightarrow \infty.

With a very small number modifying the highest degree partial term, we have a class of differential equations that doesn’t end up converging should we attempt a standard perturbation treatment. An example that is analogous is the differential equation

\begin{aligned}\epsilon \frac{du}{dx} + u = x,\end{aligned} \hspace{\stretch{1}}(4.52)

where \epsilon \ll 1 and u(0) = 1. The exact solution of this ill conditioned differential equation is

\begin{aligned}u = (1 + \epsilon) e^{-x/\epsilon} + x - \epsilon\end{aligned} \hspace{\stretch{1}}(4.53)

This is illustrated in figure (3).

Figure 3: Solution to the ill conditioned first order differential equation.

 

Study of this class of problems is called \textit{Singular perturbation theory}.

When x \gg \epsilon we have approximately

\begin{aligned}u \sim x,\end{aligned} \hspace{\stretch{1}}(4.54)

but when x \sim O(\epsilon), x/\epsilon \sim O(1) we have approximately

\begin{aligned}u \sim e^{-x/\epsilon}.\end{aligned} \hspace{\stretch{1}}(4.55)

References

[1] Wikipedia. Blasius boundary layer — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 28-March-2012]. http://en.wikipedia.org/w/index.php?title=Blasius_boundary_layer&oldid=480776115.

[2] D.J. Acheson. Elementary fluid dynamics. Oxford University Press, USA, 1990.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: