Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘Navier-Stokes’

Steady state velocity profile of stirred cup of non-bottomless coffee (continued).

Posted by peeterjoot on April 30, 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.)]

Solving for the flow below the stirring point.

Previously, we found the functional form for an azimuthal flow that has z-axis dependence. Attempting apply boundary conditions that including a cylindrical mixing device got us into trouble, since our no-slip condition would simultaneously require zero and non-zero velocity at the point of contact of the mixing interface and the bottom of the vessel (i.e. where the stir stick contacts the bottom of the cup). We can avoid this issue by constraining the mixing to only occur above the bottom of the cup, and then look at the flow that this induces below the mixing point.

Let’s attempt to solve the boundary value problem. We’ll position the mixing cylinder at a height d above the bottom of the cup, and set the radius of that inner cylinder to s as before, with the mixing occurring at an angular velocity of \Omega. We now want apply these boundary value constrains to the velocity function we’ve found for the steady state problem. Provided z < d, this will have the form

\begin{aligned}u(r, z, 0) = s \Omega \sum C_i J_1(\lambda_i r/R) \sinh(\lambda_i z/R)\end{aligned} \hspace{\stretch{1}}(1.1)

where \lambda_i are the zeros of the order one Bessel function J_1. Observe that our boundary value conditions of u(R, z, 0) = 0 and u(r, 0, 0) = 0 are automatically satisfied. Note that because of the u(R, z, 0) = 0 equality here, this form of solution is no good if we are mixing at the edge of the cup, so we require s \ne R.

Our remaining boundary value condition u(s, d, 0) = s \Omega, means that we have to solve

\begin{aligned}1 = \sum C_i J_1(\lambda_i s/R) \sinh(\lambda_i d/R)\end{aligned} \hspace{\stretch{1}}(1.2)

Having had the same sort of problem in our steady state bottomless coffee problem, we know how to solve this

\begin{aligned}C_i \sinh(\lambda_i d/R) = \frac{\int_0^1 w J_1 (\lambda_i w) dw}{\int_0^1 w J_1^2 (\lambda_i w) dw}.\end{aligned} \hspace{\stretch{1}}(1.3)

So we find that below the point of stirring (z = d), our steady state solution for the velocity should be

\begin{aligned}\mathbf{u}(r, z, 0) = s \Omega \hat{\boldsymbol{\phi}} \sum_{i=1}^\infty\frac{\int_0^1 w J_1 (\lambda_i w) dw}{\int_0^1 w J_1^2 (\lambda_i w) dw}J_1(\lambda_i r/R) \frac{\sinh(\lambda_i z/R)}{\sinh(\lambda_i d/R) }\end{aligned} \hspace{\stretch{1}}(1.4)

For a cup size of R = 5 \text{cm}, stir radius of s = 3 \text{cm}, stir depth d = 2 \text{cm}, and angular velocity \Omega = 2 \pi \text{rad}/\text{s}, we find

\begin{aligned}\begin{aligned}u(r, z, 0) &=27.1119 J_{1}(7.66341 r) \sinh (7.66341 z)-3.42819 J_{1}(14.0312 r) \sinh (14.0312 z) \\ &+4.97807 J_{1}(20.3469 r) \sinh (20.3469 z)-1.53542 J_{1}(26.6474 r) \sinh (26.6474 z) + \cdots\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.5)

We can plot this velocity field as in figure (1), but it’s hard to see the radial dependence.

Figure 1: Vector field plot of the velocity field below the stir depth.

The radial dependence of the magnitude of the velocity can be seen in figure (2), which plots u(r, d, 0). We see the zero velocity at the edges of the cup as expected, and once we hit the stir height, matching of stir velocity. An animation showing the variation of the radial velocity profile at various depths up to the stir height is available at, and also below.

Figure 2: Radial velocity dependence at the stir height.

Observing this animation we see that the velocity is dominated by the first term in the Bessel series, essentially just scaled by the hyperbolic sine that multiplies it.

Posted in Math and Physics Learning. | Tagged: , , , | Leave a Comment »

Steady state velocity profile of stirred cup of non-bottomless coffee.

Posted by peeterjoot on April 29, 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.)]


Having tackled the problem of the spin down of a bottomless cup of coffee, lets try setting up the harder problem with a non-infinite cup. It turns out that this is a lot harder. Even the steady state solution requires a superposition of Bessel functions (whereas we only had that in the bottomless problem when we tried to find the time evolution after ceasing the stirring).

I can find the form of the solution, but don’t know how to actually apply the boundary value constraints. I also find the no-slip constraints themselves become problematic, because they lead to an inconsistency at the point of contact of a moving and a static interface. Before continuing, I need to find out how to deal with that inconsistency.

Navier-Stokes for the problem.

Working in cylindrical coordinates is the only sensible option. Let’s look first at the steady state problem, with the cup being stirred at at constant rate with the unrealistic rotating cylinder stir stick used previously. We’ll assume an azimuthal velocity profile

\begin{aligned}\mathbf{u} = u(r, z, t) \hat{\boldsymbol{\phi}}.\end{aligned} \hspace{\stretch{1}}(2.1)

Let’s first verify that this satisfies our incompressible condition

\begin{aligned}\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u}&=\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z \right) \cdot ( u \hat{\boldsymbol{\phi}} ) \\ &=\not{{\hat{\mathbf{r}} \cdot \hat{\boldsymbol{\phi}}}} \partial_r u + \frac{\hat{\boldsymbol{\phi}}}{r} \cdot (u \partial_\phi \hat{\boldsymbol{\phi}} + \hat{\boldsymbol{\phi}} \not{{\partial_\phi u}})+ \not{{\hat{\mathbf{z}} \cdot \hat{\boldsymbol{\phi}}}} \partial_z u.\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.2)

The only term that survives is the \partial_\phi \hat{\boldsymbol{\phi}} = -\hat{\mathbf{r}} but that’s perpendicular to \hat{\boldsymbol{\phi}} so we have 0 = \boldsymbol{\nabla} \cdot \mathbf{u} as desired. This leaves us with

\begin{aligned}\frac{u}{r} \partial_\phi ( u \hat{\boldsymbol{\phi}} )= -\frac{1}{{\rho}} \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z \right) p+ \nu \left( \frac{1}{{r}} \partial_r ( r \partial_r ) + \frac{1}{r^2} \partial_{\phi \phi} + \partial_{zz} \right) u \hat{\boldsymbol{\phi}} + \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.3)

Splitting into \hat{\boldsymbol{\phi}}, \hat{\mathbf{r}}, \hat{\mathbf{z}} coordinates respectively we have

\begin{aligned}0 = -\frac{1}{{r \rho}} \partial_\phi p + \nu \left( \frac{1}{{r}} \partial_r ( r \partial_r u ) - \frac{u}{r^2} + \partial_{zz} u \right) \end{aligned} \hspace{\stretch{1}}(2.4a)

\begin{aligned}- \frac{u^2}{r} = -\frac{1}{{\rho}} \partial_r p \end{aligned} \hspace{\stretch{1}}(2.4b)

\begin{aligned}0 = -\frac{1}{{\rho}} \partial_z p - g.\end{aligned} \hspace{\stretch{1}}(2.4c)

Demanding a symmetrical solution kills off the pressure term in 2.4a, and we can proceed with separation of variables to find the allowable solutions for the velocity before imposing our boundary value conditions. With u = R(r) Z(z) we have

\begin{aligned}\frac{1}{{r R}} ( R' + r R'' ) - \frac{1}{r^2} = -\frac{Z''}{Z} = -\frac{1}{{a^2}} = \frac{1}{{b^2}}\end{aligned} \hspace{\stretch{1}}(2.5)

Should we pick a negative or a positive constant here (ie: trig or hyperbolic functions for Z)? Allowing for both temporarily, we find for R

\begin{aligned}0 = r R' + r^2 R'' + R \left( -1 + \frac{r^2}{a^2} \right) \end{aligned} \hspace{\stretch{1}}(2.6a)

\begin{aligned}0 = r R' + r^2 R'' + R \left( -1 - \frac{r^2}{b^2} \right).\end{aligned} \hspace{\stretch{1}}(2.6b)

Using Mathematica (phy454/coffeeCupWithBottom.cdf) to look up the solutions, we find that this was a Bessel equation with solutions

\begin{aligned}R(r) \in \text{span} \{ J_1(r/a), Y_1(r/a) \} \end{aligned} \hspace{\stretch{1}}(2.7a)

\begin{aligned}R(r) \in \text{span} \{ J_1(i r/b), Y_1(i r/b) \}.\end{aligned} \hspace{\stretch{1}}(2.7b)

However, the second set are not real valued for r > 0. This means we want the hyperbolic solutions for Z. Because we also have a boundary value constrain of u = 0 on the bottom of the cup, we can pick only hyperbolic sine solutions. As before, we’ll “stir” the coffee with a cylinder at radius s (with cup radius R), our solution has the form

\begin{aligned}u(r, z, 0) = \left\{\begin{array}{l l}s \Omega\sum C_a J_1(r/a) \sinh(z/a)& \quad \mbox{latex r < s$} \\ s \Omega\sum \left( D_a J_1(r/a) + E_a Y_1(r/a) \right) \sinh(z/a) & \quad \mbox{r \in [s, R]} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.8)$

Observe that, regardless of a we have u(0, z, 0) = 0 because J_1(0) = 0. We also can’t scale a as \lambda_i/s (where \lambda_i are the zeros of J_1) when the stirring isn’t at the edge since we need u(s, z, 0) = \Omega s. That suggests that we probably want to write our system as

\begin{aligned}u(r, z, 0) = \left\{\begin{array}{l l}s \Omega\sum C_i J_1(\lambda_i r/R) \sinh(\lambda_i z/R)& \quad \mbox{latex r < s$} \\ s \Omega\sum \left( D_i J_1(\lambda_i r/R) + E_i Y_1(\lambda_i r/R) \right) \sinh(\lambda_i z/R) & \quad \mbox{r \in [s, R]} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.9)$

The boundary value constraints of u(r = s) = \Omega s and u(r = R) = 0 require the solution of

\begin{aligned}1 = \sum C_i J_1(\lambda_i s/R) \sinh(\lambda_i z/R) \end{aligned} \hspace{\stretch{1}}(2.10a)

\begin{aligned}1 = \sum \left( D_i J_1(\lambda_i s/R) + E_i Y_1(\lambda_i s/R) \right) \sinh(\lambda_i z/R)\end{aligned} \hspace{\stretch{1}}(2.10b)

\begin{aligned}0 = \sum E_i Y_1(\lambda_i) \sinh(\lambda_i z/R).\end{aligned} \hspace{\stretch{1}}(2.10c)

We know how to solve a system like 2.10a if we did not have the \sinh term in the mix. Can we look for a numerical (least squares?) solution to this more general problem. Will it be possible to find something that’s a good fit regardless of the value of z?

We have other troubles too. We can’t simultaneously satisfy the boundary value condition of u(s, z) = \Omega s and also u(s, 0) = 0. Should we attempt something like a least squares solution and start going near z = 0 the small \sinh values near there will start forcing the C_i‘s arbitrarily high.

We have to somehow modify the no-slip condition so that when we have a moving interface in contact with a static interface we somehow deal with the fact that the no-slip constraints cannot simultaneously require the velocity to match both the moving and static interfaces at that point of contact.

Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

Surface for spinning bucket of water.

Posted by peeterjoot on April 29, 2012

[Click here for a PDF of this post with nicer formatting]


Here’s a problem from the 2009 phy1530 final that was appropriate for exam prep for this course too. It also serves as a nice example of how to determine a surface as a function of pressure, something I want to do for the non-bottomless coffee problem to be attempted.


An undergraduate student is assigned a problem about an ideal fluid rotating at a constant angular velocity \Omega under gravity g. The velocity field is \mathbf{u} = (-\Omega y, \Omega x, 0). Here, x and y are horizontal and z points up. The student is supposed to find the surfaces of constant pressure, and hence the shape of the free surface of water in a rotating bucket. The free surface corresponds to the surface for which p = p_0, where p_0 is the atmospheric pressure. Surface tension is neglected.

  • On their homework assignment, the student writes:

    “By Bernoulli’s equation:

    \begin{aligned}B = \frac{p}{\rho} + \frac{1}{{2}}u^2 + g z\end{aligned} \hspace{\stretch{1}}(2.1)

    where B is a constant. So the constant pressure surface at p = p_0 is

    \begin{aligned}z = \left( \frac{B}{g} - \frac{p_0}{\rho g}\right)- \frac{\Omega^2 }{2 g} \left( x^2 + y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.2)

    But this seems to show that the surface of the water in a rotating bucket is \textit{highest in the middle}! What is wrong with the student’s argument?

  • Write down the Euler equations in component form and integrate them directly to find the pressure p, and hence obtain the correct parabolic shape for the free surface.

Solution. Part 1.

Let’s recall how we derived Bernoulli’s theorem. We started with Navier-Stokes and used the identity

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} = \boldsymbol{\nabla} \frac{1}{{2}} \mathbf{u}^2 + (\boldsymbol{\nabla} \times \mathbf{u}) \times \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(3.3)

Navier-Stokes for a steady state incompressible flow, with external body force per unit volume \rho \mathbf{g} = -\rho \boldsymbol{\nabla} \chi take the form

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

For the non-viscous (“dry-water”) case where we take \mu = \nu \rho = 0, and treat the density \rho as a constant we find

\begin{aligned}\mathbf{u} \times (\boldsymbol{\nabla} \times \mathbf{u})=\boldsymbol{\nabla} \left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} + \chi\right).\end{aligned} \hspace{\stretch{1}}(3.5)

Observe that we only arrive at Bernoulli’s theorem if the flow is also irrotational (as well as incompressible and non-viscous), as we require an irrotational flow where \boldsymbol{\nabla} \cdot \mathbf{u} = 0 to claim that the gradient on the RHS is zero.

In this problem we do not have an irrotational flow, which can be demonstrated by direct calculation. We have

\begin{aligned}\begin{aligned}\boldsymbol{\nabla} \times \mathbf{u}&=\Omega\begin{vmatrix}\hat{\mathbf{x}} & \hat{\mathbf{y}} & \hat{\mathbf{z}} \\ \partial_x & \partial_y & 0 \\ -y & x & 0\end{vmatrix} \\ &=2 \hat{\mathbf{z}} \Omega \\ &\ne 0\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.6)

In fact we have

\begin{aligned}\begin{aligned}\mathbf{u} \times (\boldsymbol{\nabla} \times \mathbf{u})&=2 \Omega^2\begin{vmatrix}\hat{\mathbf{x}} & \hat{\mathbf{y}} & \hat{\mathbf{z}} \\ -y & x & 0  \\ 0 & 0 & 1\end{vmatrix} \\ &=2 \Omega^2 (\hat{\mathbf{x}} + \hat{\mathbf{y}})\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.7)

The closest we can get to Bernoulli’s theorem for this problem is

\begin{aligned}2 \Omega^2 (\hat{\mathbf{x}} + \hat{\mathbf{y}})= \boldsymbol{\nabla} \left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} + g z\right).\end{aligned} \hspace{\stretch{1}}(3.8)

We can say that the directional derivatives in directions perpendicular to \hat{\mathbf{x}} + \hat{\mathbf{y}} are zero, and that

\begin{aligned}\begin{aligned}2 \Omega^2 &= (\partial_x + \partial_y) \left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} + g z\right) \\ &= (\partial_x + \partial_y) \left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} \right) \\ \end{aligned}\end{aligned} \hspace{\stretch{1}}(3.9)

Perhaps those could be used to solve for the surface, but we no longer have something that is obviously integrable.

Because \mathbf{u} \cdot (\mathbf{u} \times (\boldsymbol{\nabla} \times \mathbf{u})) = 0, we can also say that

\begin{aligned}\begin{aligned}0 &= \mathbf{u} \cdot \boldsymbol{\nabla}\left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} + g z\right) \\ &= \Omega ( y \partial_x - x \partial_y )\left( \frac{1}{{2}} \mathbf{u}^2 + \frac{p}{\rho} \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.10)

Perhaps this could also be used to find the surface?

Solution. Part 2.

We want to write down the steady state, incompressible, non-viscous Navier-Stokes equations. The first of these is trivially satisfied by our assumed solution

\begin{aligned}0 = \boldsymbol{\nabla} \cdot \mathbf{u} = \partial_x (-\Omega y) + \partial_y(\Omega x).\end{aligned} \hspace{\stretch{1}}(4.11)

For the inertial term we’ve got

\begin{aligned}\begin{aligned}\mathbf{u} \cdot \boldsymbol{\nabla} \mathbf{u}&=\Omega^2 (-y \partial_x + x \partial_y (-y, x, 0) \\ &=\Omega^2 (-x, -y, 0),\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.12)

Leaving us with

\begin{aligned}\begin{aligned}-\Omega^2 x &= -\frac{1}{{\rho}} \partial_x p \\ -\Omega^2 y &= -\frac{1}{{\rho}} \partial_y p \\           0 &= -\frac{1}{{\rho}} \partial_z p - g\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.13)

Integrating these, we seek seek simultaneous solutions to

\begin{aligned}\begin{aligned}p &= \frac{1}{{2}} \rho \Omega^2 x^2 + f(y,z) \\ p &= \frac{1}{{2}} \rho \Omega^2 y^2 + g(x,z) \\ p &= h(x, y) - \rho g z.\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.14)

It’s clear that one solution would be

\begin{aligned}p = p_0 + \frac{1}{{2}} \rho \Omega^2 (x^2 + y^2) - \rho g z.\end{aligned} \hspace{\stretch{1}}(4.15)

where p_0 is some constant to be determined, dependent on where we set our origin. Putting the origin of the coordinate system at the lowest point in the parabolic profile (x, y, z) = (0, 0, 0), we have p(0, 0, 0) = p_0, which fixes p_0 as the atmospheric pressure. If the radius of the bucket is R, the max height h of the surface above that point is also found on this surface of constant pressure

\begin{aligned}p_0 = p_0 + \frac{1}{{2}} \rho \Omega^2 R^2 - \rho g h,\end{aligned} \hspace{\stretch{1}}(4.16)


\begin{aligned}h = \frac{\Omega^2 R^2 }{2 g}.\end{aligned} \hspace{\stretch{1}}(4.17)

Appendix. Proof of vector identities used.

\begin{aligned}\begin{aligned}\left( \boldsymbol{\nabla} \frac{1}{{2}} \mathbf{u}^2 + (\boldsymbol{\nabla} \times \mathbf{u}) \times \mathbf{u} \right)_i&=\partial_i \frac{1}{{2}} u_j u_j + \partial_a u_b \epsilon_{a b r} u_s \epsilon_{r s i} \\ &=u_j \partial_i u_j + \partial_a u_b u_s \delta^{[a b]}_{s i} \\ &=u_j \partial_i u_j + u_s \partial_s u_i - u_s \partial_i u_s \\ &= (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u})_i \end{aligned}\end{aligned} \hspace{\stretch{1}}(5.18)

Also observe that our claim that \mathbf{u} \cdot (\mathbf{u} \times (\boldsymbol{\nabla} \mathbf{u})) = 0 follows easily after expansion in coordinates

\begin{aligned}\mathbf{u} \cdot (\mathbf{u} \times (\boldsymbol{\nabla} \times \mathbf{u}) )=u_i u_s ( \partial_s u_i - \partial_i u_s ).\end{aligned} \hspace{\stretch{1}}(5.19)

We’ve got a symmetric and antisymmetric factor in the summation, so the end result is zero.

Posted in Math and Physics Learning. | Tagged: , , , | Leave a Comment »

Continuum mechanics fluids review.

Posted by peeterjoot on April 25, 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.)]


Review of key ideas and equations from the fluid dynamics portion of the class.

Vector displacements.

Those portions of the theory of elasticity that we did cover have the appearance of providing some logical context for the derivation of the Navier-Stokes equation. Our starting point is almost identical, but we now look at displacements that vary with time, forming

\begin{aligned}d\mathbf{x}' = d\mathbf{x} + d\mathbf{u} \delta t.\end{aligned} \hspace{\stretch{1}}(2.1)

We compute a first order Taylor expansion of this differential, defining a symmetric strain and antisymmetric vorticity tensor


\begin{aligned}e_{ij} = \frac{1}{{2}} \left(\frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_j}}{\partial {x_i}} \right).\end{aligned} \hspace{\stretch{1}}(2.2a)

\begin{aligned}\omega_{ij} = \frac{1}{{2}} \left(\frac{\partial {u_i}}{\partial {x_j}}-\frac{\partial {u_j}}{\partial {x_i}} \right)\end{aligned} \hspace{\stretch{1}}(2.2b)


Allowing us to write

\begin{aligned}dx_i' = dx_i + e_{ij} dx_j \delta t + \omega_{ij} dx_j \delta t.\end{aligned} \hspace{\stretch{1}}(2.3)

We introduced vector and dual vector forms of the vorticity tensor with


\begin{aligned}\Omega_k = \frac{1}{{2}} \partial_i u_j \epsilon_{i j k}\end{aligned} \hspace{\stretch{1}}(2.4a)

\begin{aligned}\omega_{i j} = -\Omega_k \epsilon_{i j k},\end{aligned} \hspace{\stretch{1}}(2.4b)




\begin{aligned}\boldsymbol{\omega} = \boldsymbol{\nabla} \times \mathbf{u}\end{aligned} \hspace{\stretch{1}}(2.5a)

\begin{aligned}\boldsymbol{\Omega} = \frac{1}{{2}} (\boldsymbol{\omega})_a \mathbf{e}_a.\end{aligned} \hspace{\stretch{1}}(2.5b)


We were then able to put our displacement differential into a partial vector form

\begin{aligned}d\mathbf{x}' = d\mathbf{x} + \left( \mathbf{e}_i (e_{ij} \mathbf{e}_j) \cdot d\mathbf{x} + \boldsymbol{\Omega} \times d\mathbf{x} \right) \delta t.\end{aligned} \hspace{\stretch{1}}(2.6)

Relative change in volume

We are able to identity the divergence of the displacement as the relative change in volume per unit time in terms of the strain tensor trace (in the basis for which the strain is diagonal at a given point)

\begin{aligned}\frac{dV' - dV}{dV \delta t} = \boldsymbol{\nabla} \cdot \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(3.7)

Conservation of mass

Utilizing Green’s theorem we argued that

\begin{aligned}\int \left( \frac{\partial {\rho}}{\partial {t}} + \boldsymbol{\nabla} \cdot (\rho \mathbf{u}) \right) dV = 0.\end{aligned} \hspace{\stretch{1}}(4.8)

We were able to relate this to rate of change of density, computing

\begin{aligned}\frac{d\rho}{dt} = \frac{\partial {\rho}}{\partial {t}} + \mathbf{u} \cdot \boldsymbol{\nabla} \rho =- \rho \boldsymbol{\nabla} \cdot \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(4.9)

An important consequence of this is that for incompressible fluids (the only types of fluids considered in this course) the divergence of the displacement \boldsymbol{\nabla} \cdot \mathbf{u} = 0.

Constitutive relation

We consider only Newtonian fluids, for which the stress is linearly related to the strain. We will model fluids as disjoint sets of hydrostatic materials for which the constitutive relation was previously found to be

\begin{aligned}\sigma_{ij} = - p \delta_{ij} + 2 \mu e_{ij}.\end{aligned} \hspace{\stretch{1}}(5.10)

Conservation of momentum (Navier-Stokes)

As in elasticity, our momentum conservation equation had the form

\begin{aligned}\rho \frac{du_i}{dt} = \frac{\partial {\sigma_{ij}}}{\partial {x_j}} + \rho f_i,\end{aligned} \hspace{\stretch{1}}(6.11)

where f_i are the components of the external (body) forces per unit volume acting on the fluid.

Utilizing the constitutive relation and explicitly evaluating the stress tensor divergence {\partial {\sigma_{ij}}}/{\partial {x_j}} we find

\begin{aligned}\rho \frac{d\mathbf{u}}{dt} =\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u}+ \mu \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{u}) + \rho \mathbf{f}.\end{aligned} \hspace{\stretch{1}}(6.12)

Since we treat only incompressible fluids in this course we can decompose this into a pair of equations


\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u}+ \rho \mathbf{f}.\end{aligned} \hspace{\stretch{1}}(6.13a)

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0\end{aligned} \hspace{\stretch{1}}(6.13b)


No slip condition

We’ll find in general that we have to solve for our boundary value conditions. One of the important constraints that we have to do so will be a requirement (experimentally motivated) that our velocities match at an interface. This was illustrated with a rocker tank video in class.

This is the no-slip condition, and includes a requirement that the fluid velocity at the boundary of a non-moving surface is zero, and that the fluid velocity on the boundary of a moving surface matches the rate of the surface itself.

For fluids A and B separated at an interface with unit normal \hat{\mathbf{n}} and unit tangent \hat{\boldsymbol{\tau}} we wrote the no-slip condition as


\begin{aligned}\mathbf{u}_A \cdot \hat{\boldsymbol{\tau}} = \mathbf{u}_B \cdot \hat{\boldsymbol{\tau}}\end{aligned} \hspace{\stretch{1}}(7.14a)

\begin{aligned}\mathbf{u}_A \cdot \hat{\mathbf{n}} = \mathbf{u}_B \cdot \hat{\mathbf{n}}.\end{aligned} \hspace{\stretch{1}}(7.14b)


For the problems we attempt, it will often be enough to consider only the tangential component of the velocity.

Traction vector matching at an interface.

As well as matching velocities, we have a force balance requirement at any interface. This will be expressed in terms of the traction vector

\begin{aligned}\mathbf{T} = \mathbf{e}_i \sigma_{ij} n_j = \boldsymbol{\sigma} \cdot \hat{\mathbf{n}}\end{aligned} \hspace{\stretch{1}}(8.15)

where \hat{\mathbf{n}} = n_j \mathbf{e}_i is the normal pointing from the interface into the fluid (so the traction vector represents the force of the interface on the fluid). When that interface is another fluid, we are able to calculate the force of one fluid on the other.

In addition the the constraints provided by the no-slip condition, we’ll often have to constrain our solutions according to the equality of the tangential components of the traction vector

\begin{aligned}{\left.{{\tau_i (\sigma_{ij} n_j)}}\right\vert}_{{A}} ={\left.{{\tau_i (\sigma_{ij} n_j)}}\right\vert}_{{B}},\end{aligned} \hspace{\stretch{1}}(8.16)

We’ll sometimes also have to consider, especially when solving for the pressure, the force balance for the normal component of the traction vector at the interface too

\begin{aligned}{\left.{{n_i (\sigma_{ij} n_j)}}\right\vert}_{{A}} ={\left.{{n_i (\sigma_{ij} n_j)}}\right\vert}_{{B}}.\end{aligned} \hspace{\stretch{1}}(8.17)

As well as having a messy non-linear PDE to start with, our boundary value constraints can be very complicated, making the subject rich and tricky.


A number of problems we did asked for the flux rate. A slightly more sensible physical quantity is the mass flux, which adds the density into the mix

\begin{aligned}\int \frac{dm}{dt} = \rho \int \frac{dV}{dt} = \rho \int (\mathbf{u} \cdot \hat{\mathbf{n}}) dA\end{aligned} \hspace{\stretch{1}}(8.18)

Worked problems from class.

A number of problems were tackled in class

  1. channel flow with external pressure gradient
  2. shear flow
  3. pipe (Poiseuille) flow
  4. steady state gravity driven film flow down a slope.

Channel and shear flow

An example that shows many of the features of the above problems is rectilinear flow problem with a pressure gradient and shearing surface. As a review let’s consider fluid flowing between surfaces at z = \pm h, the lower surface moving at velocity v and pressure gradient dp/dx = -G we find that Navier-Stokes for an assumed flow of \mathbf{u} = u(z) \hat{\mathbf{x}} takes the form

\begin{aligned}0 &= \partial_x u + \partial_y (0) + \partial_z (0) \\ u \not{{\partial_x u}} &= - \partial_x p + \mu \partial_{zz} u \\ 0 &= -\partial_y p \\ 0 &= -\partial_z p\end{aligned} \hspace{\stretch{1}}(9.19)

We find that this reduces to

\begin{aligned}\frac{d^2 u}{dz^2} = -\frac{G}{\mu}\end{aligned} \hspace{\stretch{1}}(9.23)

with solution

\begin{aligned}u(z) = \frac{G}{2\mu}(h^2 - z^2) + A (z + h) + B.\end{aligned} \hspace{\stretch{1}}(9.24)

Application of the no-slip velocity matching constraint gives us in short order

\begin{aligned}u(z) = \frac{G}{2\mu}(h^2 - z^2) + v \left( 1 - \frac{1}{{2h}} (z + h) \right).\end{aligned} \hspace{\stretch{1}}(9.25)

With v = 0 this is the channel flow solution, and with G = 0 this is the shearing flow solution.

Having solved for the velocity at any height, we can also solve for the mass or volume flux through a slice of the channel. For the mass flux \rho Q per unit time (given volume flux Q)

\begin{aligned}\int \frac{dm}{dt} =\int \rho \frac{dV}{dt} =\rho (\Delta A) \int \mathbf{u} \cdot \hat{\boldsymbol{\tau}},\end{aligned} \hspace{\stretch{1}}(9.26)

we find

\begin{aligned}\rho Q =\rho (\Delta y) \left( \frac{2 G h^3}{3 \mu} + h v\right).\end{aligned} \hspace{\stretch{1}}(9.27)

We can also calculate the force of the boundaries on the fluid. For example, the force per unit volume of the boundary at z = \pm h on the fluid is found by calculating the tangential component of the traction vector taken with normal \hat{\mathbf{n}} = \mp \hat{\mathbf{z}}. That tangent vector is found to be

\begin{aligned}\boldsymbol{\sigma} \cdot (\pm \hat{\mathbf{n}}) = -p \hat{\mathbf{z}} \pm 2 \mu \mathbf{e}_i e_{ij} \delta_{j 3} = - p \hat{\mathbf{z}} \pm \hat{\mathbf{x}} \mu \frac{\partial {u}}{\partial {z}}.\end{aligned} \hspace{\stretch{1}}(9.28)

The tangential component is the \hat{\mathbf{x}} component evaluated at z = \pm h, so for the lower and upper interfaces we have

\begin{aligned}{\left.{{(\boldsymbol{\sigma} \cdot \hat{\mathbf{n}}) \cdot \hat{\mathbf{x}}}}\right\vert}_{{z = -h}} &= -G (-h) - \frac{v \mu}{2 h} \\ {\left.{{(\boldsymbol{\sigma} \cdot -\hat{\mathbf{n}}) \cdot \hat{\mathbf{x}}}}\right\vert}_{{z = +h}} &= -G (+h) + \frac{v \mu}{2 h},\end{aligned} \hspace{\stretch{1}}(9.29)

so the force per unit area that the boundary applies to the fluid is

\begin{aligned}\text{force per unit length of lower interface on fluid} &= L \left( G h - \frac{v \mu}{2 h} \right) \\ \text{force per unit length of upper interface on fluid} &= L \left( -G h + \frac{v \mu}{2 h} \right).\end{aligned} \hspace{\stretch{1}}(9.31)

Does the sign of the velocity term make sense? Let’s consider the case where we have a zero pressure gradient and look at the lower interface. This is the force of the interface on the fluid, so the force of the fluid on the interface would have the opposite sign

\begin{aligned}\frac{v \mu}{2 h}.\end{aligned} \hspace{\stretch{1}}(9.33)

This does seem reasonable. Our fluid flowing along with a positive velocity is imparting a force on what it is flowing over in the same direction.


We covered hydrostatics as a separate topic, where it was argued that the pressure p in a fluid, given atmospheric pressure p_a and height from the surface was

\begin{aligned}p = p_a + \rho g h.\end{aligned} \hspace{\stretch{1}}(10.34)

As noted below in the surface tension problem, this is also a consequence of Navier-Stokes for \mathbf{u} = 0 (following from 0 = -\boldsymbol{\nabla} p + \rho \mathbf{g}).

We noted that replacing the a mass of water with something of equal density would not change the non-dynamics of the situation. We then went on to define Buoyancy force, the difference in weight of the equivalent volume of fluid and the weight of the object.

Mass conservation through apertures.

It was noted that mass conservation provides a relationship between the flow rates through apertures in a closed pipe, since we must have

\begin{aligned}\rho_1 A_1 v_1 = \rho_2 A_2 v_2,\end{aligned} \hspace{\stretch{1}}(11.35)

and therefore for incompressible fluids

\begin{aligned}A_1 v_1 = A_2 v_2.\end{aligned} \hspace{\stretch{1}}(11.36)

So if A_1 > A_2 we must have v_1 < v_2.

Curve for tap discharge.

We can use this to get a rough idea what the curve for water coming out a tap would be. Suppose we measure the volume flux, putting a measuring cup under the tap, and timing how long it takes to fill up. We then measure the radii at different points. This can be done from a photo as in figure (1).

Figure 1: Tap flow measurement.

After making the measurement, we can get an idea of the velocity between two points given a velocity estimate at a point higher in the discharge. For a plain old falling mass, our final velocity at a point measured from where the velocity was originally measured can be found from Newton’s law

\begin{aligned}\Delta v = g t\end{aligned} \hspace{\stretch{1}}(11.37)

\begin{aligned}\Delta z = \frac{1}{{2}} g t^2 + v_0 t\end{aligned} \hspace{\stretch{1}}(11.38)

Solving for v_f = v_0 + \Delta v, we find

\begin{aligned}v_f = v_0 \sqrt{ 1 + \frac{2 g \Delta z}{v_0^2} }.\end{aligned} \hspace{\stretch{1}}(11.39)

Mass conservation gives us

\begin{aligned}v_0 \pi R^2 = v_f \pi r^2\end{aligned} \hspace{\stretch{1}}(11.40)


\begin{aligned}r(\Delta z) = R \sqrt{ \frac{v_0}{v_f} } = R \left( 1 + \frac{2 g \Delta z}{v_0^2} \right)^{-1/4}.\end{aligned} \hspace{\stretch{1}}(11.41)

For the image above I measured a flow rate of about 250 ml in 10 seconds. With that, plus the measured radii at 0 and 6 \text{cm}, I calculated that the average fluid velocity was 0.9 \text{m}/\text{s}, vs a free fall rate increase of 1.3 \text{m}/\text{s}. Not the best match in the world, but that’s to be expected since the velocity has been considered uniform throughout the stream profile, which would not actually be the case. A proper treatment would also have to treat viscosity and surface tension.

In figure (2) is a plot of the measured radial distance compared to what was computed with 11.41. The blue line is the measured width of the stream as measured, the red is a polynomial curve fitted to the raw data, and the green is the computed curve above.

Figure 2: Comparison of measured stream radii and calculated.

Bernoulli equation.

With the body force specified in gradient for

\begin{aligned}\mathbf{g} = - \boldsymbol{\nabla} \chi\end{aligned} \hspace{\stretch{1}}(12.42)

and utilizing the vector identity

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = \boldsymbol{\nabla} \times (\boldsymbol{\nabla} \times \mathbf{u}) + \boldsymbol{\nabla} \left( \frac{1}{{2}} \mathbf{u}^2 \right),\end{aligned} \hspace{\stretch{1}}(12.43)

we are able to show that the steady state, irrotational, non-viscous Navier-Stokes equation takes the form

\begin{aligned}\boldsymbol{\nabla} \left( \frac{p}{\rho} + \chi + \frac{1}{{2}} \mathbf{u}^2 \right) = 0.\end{aligned} \hspace{\stretch{1}}(12.44)


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

This is the Bernoulli equation, and the constants introduce the concept of streamline.

FIXME: I think this could probably be used to get a better idea what the tap stream radius is, than the method used above. Consider the streamline along the outermost surface. That way you don’t have to assume that the flow is at the average velocity rate uniformly throughout the stream. Try this later.

Surface tension.

Surfaces, normals and tangents.

We reviewed basic surface theory, noting that we can parameterize a surface as in the following example

\begin{aligned}\phi = z - h(x, t) = 0.\end{aligned} \hspace{\stretch{1}}(13.46)

Computing the gradient we find

\begin{aligned}\boldsymbol{\nabla} \phi = -\hat{\mathbf{x}} \frac{\partial {h}}{\partial {x}} \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(13.47)

Recalling that the gradient is normal to the surface we can compute the unit normal and unit tangent vectors


\begin{aligned}\hat{\mathbf{n}} =\frac{1}{{\sqrt{ 1 + \left( \frac{\partial {h}}{\partial {x}} \right)^2 }}}\left( -\frac{\partial {h}}{\partial {x}}, 1 \right)\end{aligned} \hspace{\stretch{1}}(13.48a)

\begin{aligned}\hat{\boldsymbol{\tau}} =\frac{1}{{\sqrt{ 1 + \left( \frac{\partial {h}}{\partial {x}} \right)^2 }}}\left( 1, \frac{\partial {h}}{\partial {x}} \right)\end{aligned} \hspace{\stretch{1}}(13.48b)


Laplace pressure.

We covered some aspects of this topic in class. [1] covers this topic in typical fairly hard to comprehend) detail, but there’s lots of valuable info there. section 2.4.9-2.4.10 of [2] has small section that’s a bit easier to understand, with less detail. Recommended in that text is the “Surface Tension in Fluid Mechanics” movie which can be found on youtube in three parts \youtubehref{DkEhPltiqmo}, \youtubehref{yiixltf\_HKw}, \youtubehref{5d6efCcwkWs}, which is very interesting and entertaining to watch.

It was argued in class that the traction vector differences at the surfaces between a pair of fluids have the form

\begin{aligned}\mathbf{t}_2 - \mathbf{t}_1 = -\frac{\sigma}{2 R} \hat{\mathbf{n}} - \boldsymbol{\nabla}_I \sigma\end{aligned} \hspace{\stretch{1}}(13.49)

where \boldsymbol{\nabla}_I = \boldsymbol{\nabla} - \hat{\mathbf{n}} (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) is the tangential (interfacial) gradient, \sigma is the surface tension, a force per unit length value, and R is the radius of curvature.

In static equilibrium where \mathbf{t} = -p \hat{\mathbf{n}} (since \boldsymbol{\sigma} = 0 if \mathbf{u} = 0), then dotting with \hat{\mathbf{n}} we must then have

\begin{aligned}p_2 - p_1 = \frac{\sigma}{2 R} \end{aligned} \hspace{\stretch{1}}(13.50)

Surface tension gradients.

Considering the tangential component of the traction vector difference we find

\begin{aligned}(\mathbf{t}_2 - \mathbf{t}_1) \cdot \hat{\boldsymbol{\tau}} = - \hat{\boldsymbol{\tau}} \cdot \boldsymbol{\nabla}_I \sigma\end{aligned} \hspace{\stretch{1}}(13.51)

If the fluid is static (for example, has none of the creep that we see in the film) then we must have \boldsymbol{\nabla}_I \sigma = 0. It’s these gradients that are responsible for capillary flow and other related surface tension driven motion (lots of great examples of that in the film).

Surface tension for a spherical bubble.

In the film above it is pointed out that the surface tension equation we were shown in class

\begin{aligned}\Delta p = \frac{2 \sigma}{R},\end{aligned} \hspace{\stretch{1}}(13.52)

is only for spherical objects that have a single radius of curvature. This formula can in fact be derived with a simple physical argument, stating that the force generated by the surface tension \sigma along the equator of a bubble (as in figure (3)), in a fluid would be balanced by the difference in pressure times the area of that equatorial cross section. That is

Figure 3: Spherical bubble in liquid.

\begin{aligned}\sigma 2 \pi R = \Delta p \pi R^2\end{aligned} \hspace{\stretch{1}}(13.53)

Observe that we obtain 13.52 after dividing through by the area.

A sample problem. The meniscus curve.

To get a better feeling for this, let’s look to a worked problem. The most obvious one to try to attempt is the shape of a meniscus of water against a wall. This problem is worked in [1], but it is worth some extra notes. As in the text we’ll work with z axis up, and the fluid up against a wall at x = 0 as illustrated in figure (4).

Figure 4: Curvature of fluid against a wall.

The starting point is a variation of what we have in class

\begin{aligned}p_1 - p_2 = \sigma \left( \frac{1}{{R_1}} + \frac{1}{{R_2}} \right),\end{aligned} \hspace{\stretch{1}}(13.54)

where p_2 is the atmospheric pressure, p_1 is the fluid pressure, and the (signed!) radius of curvatures positive if pointing into medium 1 (the fluid).

For fluid at rest, Navier-Stokes takes the form

\begin{aligned}0 = -\boldsymbol{\nabla} p_1 + \rho \mathbf{g}.\end{aligned} \hspace{\stretch{1}}(13.55)

With \mathbf{g} = -g \hat{\mathbf{z}} we have

\begin{aligned}0 = -\frac{\partial {p_1}}{\partial {z}} - \rho g,\end{aligned} \hspace{\stretch{1}}(13.56)


\begin{aligned}p_1 = \text{constant} - \rho g z.\end{aligned} \hspace{\stretch{1}}(13.57)

We have p_2 = p_a, the atmospheric pressure, so our pressure difference is

\begin{aligned}p_1 - p_2 = \text{constant} - \rho g z.\end{aligned} \hspace{\stretch{1}}(13.58)

We have then

\begin{aligned}\text{constant} -\frac{\rho g z}{\sigma} = \frac{1}{{R_1}} + \frac{1}{{R_2}}.\end{aligned} \hspace{\stretch{1}}(13.59)

One of our axis of curvature directions is directly along the y axis so that curvature is zero 1/R_1 = 0. We can fix the constant by noting that at x = \infty, z = 0, we have no curvature 1/R_2 = 0. This gives

\begin{aligned}\text{constant} -0 = 0 + 0.\end{aligned} \hspace{\stretch{1}}(13.60)

That leaves just the second curvature to determine. For a curve z = z(x) our absolute curvature, according to [3] is

\begin{aligned}{\left\lvert{\frac{1}{{R_2}}}\right\rvert} = \frac{{\left\lvert{z''}\right\rvert}}{(1 + (z')^2)^{3/2}}.\end{aligned} \hspace{\stretch{1}}(13.61)

Now we have to fix the sign. I didn’t recall any sort of notion of a signed radius of curvature, but there’s a blurb about it on the curvature article above, including a nice illustration of signed radius of curvatures can be found in this wikipedia radius of curvature figure for a Lemniscate. Following that definition for a curve such as z(x) = (1-x)^2 we’d have a positive curvature, but the text explicitly points out that the curvatures are will be set positive if pointing into the medium. For us to point the normal into the medium as in the figure, we have to invert the sign, so our equation to solve for z is given by

\begin{aligned}-\frac{\rho g z}{\sigma} = -\frac{z''}{(1 + (z')^2)^{3/2}}.\end{aligned} \hspace{\stretch{1}}(13.62)

The text introduces the capillary constant

\begin{aligned}a = \sqrt{2 \sigma/ g \rho}.\end{aligned} \hspace{\stretch{1}}(13.63)

Using that capillary constant a to tidy up a bit and multiplying by a z' integrating factor we have

\begin{aligned}-\frac{2 z z'}{a^2} = -\frac{z'' z'}{(1 + (z')^2)^{3/2}},\end{aligned} \hspace{\stretch{1}}(13.64)

we can integrate to find

\begin{aligned}A - \frac{z^2}{a^2} = \frac{1}{(1 + (z')^2)^{1/2}}.\end{aligned} \hspace{\stretch{1}}(13.65)

Again for x = \infty we have z = 0, z' = 0, so A = 1. Rearranging we have

\begin{aligned}\int dx = \int dz \left( \frac{1}{{(1 - z^2/a^2)^2}} - 1 \right)^{-1/2}.\end{aligned} \hspace{\stretch{1}}(13.66)

Integrating this with Mathematica I get

\begin{aligned}x - x_0 =\sqrt{2 a^2-z^2} \sgn(a-z)+ \frac{a}{\sqrt{2}} \ln \left(\frac{a \left(2 a-\sqrt{4 a^2-2 z^2} \sgn(a-z)\right)}{z}\right).\end{aligned} \hspace{\stretch{1}}(13.67)

It looks like the constant would have to be fixed numerically. We require at x = 0

\begin{aligned}z'(0) = \frac{-\cos\theta}{\sin\theta} = -\cot \theta,\end{aligned} \hspace{\stretch{1}}(13.68)

but we don’t have an explicit function for z.

Non-dimensionality and scaling.

With the variable transformations

\begin{aligned}\mathbf{u} &\rightarrow U \mathbf{u}' \\ p &\rightarrow \rho U^2 p' \\ t &\rightarrow \frac{L}{U} t' \\ \boldsymbol{\nabla} &\rightarrow \frac{1}{{L}} \boldsymbol{\nabla}'\end{aligned} \hspace{\stretch{1}}(14.69)

we can put Navier-Stokes in dimensionless form

\begin{aligned}\frac{\partial {\mathbf{u}'}}{\partial {t'}} + (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = \boldsymbol{\nabla}' p' + \frac{1}{{R}} \boldsymbol{\nabla}' \mathbf{u}'.\end{aligned} \hspace{\stretch{1}}(14.73)

Here R is Reynold’s number

\begin{aligned}R = \frac{L U}{\nu}\end{aligned} \hspace{\stretch{1}}(14.74)

A relatively high or low Reynold’s number will effect whether viscous or inertial effects dominate

\begin{aligned}R \sim \frac{{\left\lvert{\text{effect of inertia}}\right\rvert}}{{\left\lvert{\text{effect of viscosity}}\right\rvert}} \sim \frac{ {\left\lvert{\rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} }\right\rvert} }{{\left\lvert{ \mu \boldsymbol{\nabla}^2 \mathbf{u}}\right\rvert}}\end{aligned} \hspace{\stretch{1}}(14.75)

The importance of examining where one of these effects can dominate was clear in the Blassius problem, where doing so allowed for an analytic solution that would not have been possible otherwise.

Eulerian and Lagrangian.

We defined

  1. Lagrangian: the observer is moving with the fluid.
  2. Eulerian: the observer is fixed in space, watching the fluid.

Boundary layers.

Impulsive flow.

We looked at the time dependent unidirectional flow where

\begin{aligned}\mathbf{u} &= u(y, t) \hat{\mathbf{x}} \\ \frac{\partial {u}}{\partial {t}} &= \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2} \\ u(0, t) &=\left\{\begin{array}{l l}0 & \quad \mbox{for latex t < 0$} \\ U & \quad \mbox{for t \ge 0}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(16.76)$

and utilized a similarity variable \eta with u = U f(\eta)

\begin{aligned}\eta = \frac{y}{2 \sqrt{\nu t}}\end{aligned} \hspace{\stretch{1}}(16.79)

and were able to show that

\begin{aligned}u(y, t) = U_0 \left(1 - \text{erf}\left(\frac{y}{2 \sqrt{\nu t}}\right) \right).\end{aligned} \hspace{\stretch{1}}(16.80)

The aim of this appears to be as an illustration that the boundary layer thickness \delta grows with \sqrt{\nu t}.

FIXME: really need to plot 16.80.

Oscillatory flow.

Another worked problem in the boundary layer topic was the Stokes boundary layer problem with a driving interface of the form

\begin{aligned}U(t) = U_0 e^{i \Omega t}\end{aligned} \hspace{\stretch{1}}(16.81)

with an assumed solution of the form

\begin{aligned}u(y, t) = f(y) e^{i \Omega t},\end{aligned} \hspace{\stretch{1}}(16.82)

we found


\begin{aligned}u(y, t) =U_0 e^{-\lambda y} \cos\left( -i (\lambda y - \Omega t) \right).\end{aligned} \hspace{\stretch{1}}(16.83a)

\begin{aligned}\lambda = \sqrt{\frac{\Omega}{2 \nu}}\end{aligned} \hspace{\stretch{1}}(16.83b)


This was a bit more obvious as a boundary layer illustration since we see the exponential drop off with every distance multiple of \sqrt{\frac{2 \nu}{\Omega}}.

Blassius problem (boundary layer thickness in flow over plate).

We examined the scaling off all the terms in the Navier-Stokes equations given a velocity scale U, vertical length scale \delta and horizontal length scale L. This, and the application of Bernoulli’s theorem allowed us to make construct an approximation for Navier-Stokes in the boundary layer


\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}}(16.84a)

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

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


With boundary conditions

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

With a similarity variable

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

and stream functions

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


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

we were able to show that our velocity dependence was given by the solutions of

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

This was done much more clearly in [4] and I worked this problem myself with a hybrid approach (non-dimensionalising as done in class).

FIXME: the end result of this is a plot (a nice one can be found in [5]). That plot ends up being one that’s done in terms of the similarity variable \eta. It’s not clear to me how this translates into an actual velocity profile. Should plot these out myself to get a feel for things.

Singular perturbation theory.

The non-dimensional form of Navier-Stokes had the form

\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}}(17.93)

where the inverse of Reynold’s number

\begin{aligned}\text{Re} = \frac{U L}{\nu}\end{aligned} \hspace{\stretch{1}}(17.94)

can potentially get very small. That introduces an ill-conditioning into the problems that can make life more interesting.

We looked at a couple of simple LDE systems that had this sort of ill conditioning. One of them was

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

for which the exact solution was found to be

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

The rough idea is that we can look in the domain where x \sim \epsilon and far from that. In this example, with x far from the origin we have roughly

\begin{aligned}\epsilon \times 1 + u = x \approx 0 + u\end{aligned} \hspace{\stretch{1}}(17.97)

so we have an asymptotic solution close to u = x. Closer to the origin where x \sim O(\epsilon) we can introduce a rescaling x = \epsilon y to find

\begin{aligned}\epsilon \frac{1}{{\epsilon}} \frac{du}{dy} + u = \epsilon y.\end{aligned} \hspace{\stretch{1}}(17.98)

This gives us

\begin{aligned}\frac{du}{dy} + u \approx 0,\end{aligned} \hspace{\stretch{1}}(17.99)

for which we find

\begin{aligned}u \propto e^{-y} = e^{-x/\epsilon}.\end{aligned} \hspace{\stretch{1}}(17.100)


We characterized stability in terms of displacements writing

\begin{aligned}\delta x = e^{(\sigma_R + i \sigma_I) t}\end{aligned} \hspace{\stretch{1}}(18.101)

and defining

  • Oscillatory unstability. \sigma_{\text{R}} = 0, \sigma_{\text{I}} > 0.
  • Marginal unstability. \sigma_{\text{I}} = 0, \sigma_{\text{R}} > 0.
  • Neutral stability. \sigma_{\text{I}} = 0, \sigma_{\text{R}} = 0.

Thermal stability: Rayleigh-Benard problem.

We considered the Rayleigh-Benard problem, looking at thermal effects in a cavity. Assuming perturbations of the form

\begin{aligned}\mathbf{u} &= \mathbf{u}_\text{base} + \delta \mathbf{u} = 0 + \delta \mathbf{u} \\ p &= p_s + \delta p \\ \rho &= \rho_s + \delta \rho\end{aligned}

and introducing an equation for the base state

\begin{aligned}\boldsymbol{\nabla} p_s = -\rho_s \hat{\mathbf{z}} g\end{aligned} \hspace{\stretch{1}}(18.102)

we found

\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \nu \boldsymbol{\nabla}^2 \right) \delta \mathbf{u} = -\frac{1}{{\rho_s}} \boldsymbol{\nabla} \delta p - \frac{\delta \rho}{\rho_s} \hat{\mathbf{z}} g\end{aligned} \hspace{\stretch{1}}(18.103)

Operating on this with {\partial {}}/{\partial {z}} \boldsymbol{\nabla} \cdot () we find

\begin{aligned}\boldsymbol{\nabla}^2 \frac{\partial {\delta p}}{\partial {z}} = -g \frac{\partial^2 {{\delta \rho}}}{\partial {{z}}^2},\end{aligned} \hspace{\stretch{1}}(18.104)

from which we apply back to 18.103 and take just the z component to find

\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \nu \boldsymbol{\nabla}^2 \right) \delta w = -\frac{1}{{\rho_s}} \frac{\partial {\delta p}}{\partial {z}} - \frac{\delta \rho}{\rho_s} g\end{aligned} \hspace{\stretch{1}}(18.105)

With an assumption that density change and temperature are linearly related

\begin{aligned}\delta \rho = - \rho_s \alpha \delta T,\end{aligned} \hspace{\stretch{1}}(18.106)

and operating with the Laplacian we end up with a relation that follows from the momentum balance equation

\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \nu \boldsymbol{\nabla}^2 \right) \boldsymbol{\nabla}^2 \delta w=g \alpha \left(\frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right)\delta T.\end{aligned} \hspace{\stretch{1}}(18.107)

We also applied our perturbation to the energy balance equation

\begin{aligned}\frac{\partial {T}}{\partial {t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) T = \kappa \boldsymbol{\nabla}^2 T\end{aligned} \hspace{\stretch{1}}(18.108)

We determined that the base state temperature obeyed

\begin{aligned}\kappa \frac{\partial^2 {{}}}{\partial {{z}}^2} T_s = 0,\end{aligned} \hspace{\stretch{1}}(18.109)

with solution

\begin{aligned}T_s = T_0 - \frac{\Delta T}{d} z.\end{aligned} \hspace{\stretch{1}}(18.110)

This and application of the perturbation gave us

\begin{aligned}\frac{\partial {\delta T}}{\partial {t}} + \delta \mathbf{u} \cdot \boldsymbol{\nabla} T_s = \kappa \boldsymbol{\nabla}^2 \delta T.\end{aligned} \hspace{\stretch{1}}(18.111)

We used this to non-dimensionalize with

\begin{aligned}\begin{array}{l l}x,y,z & \quad \mbox{with latex d$} \\ t & \quad \mbox{with d^2/\nu} \\ \delta w & \quad \mbox{with \kappa/d} \\ \delta T & \quad \mbox{with \Delta T}\end{array}\end{aligned} \hspace{\stretch{1}}(18.112)$

And found (primes dropped)


\begin{aligned}\boldsymbol{\nabla}^2 \left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta w = \mathcal{R} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2} +\frac{\partial^2 {{}}}{\partial {{y}}^2} \right) \delta T\end{aligned} \hspace{\stretch{1}}(18.113a)

\begin{aligned}\left( \text{P}_r \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta T = \delta w,\end{aligned} \hspace{\stretch{1}}(18.113b)


where we’ve introduced the Rayleigh number and Prandtl number’s

\begin{aligned}\mathcal{R} = \frac{g \alpha \Delta T d^3}{\nu \kappa},\end{aligned} \hspace{\stretch{1}}(18.114a)

\begin{aligned}\text{P}_r = \frac{\nu}{\kappa}\end{aligned} \hspace{\stretch{1}}(18.114b)


We were able to construct some approximate solutions for a problem similar to these equations using an assumed solution form

\begin{aligned}\delta w &= w(z) e^{ i ( k_1 x + k_2 y) + \sigma t} \\ \delta T &= \Theta(z) e^{ i ( k_1 x + k_2 y) + \sigma t}\end{aligned} \hspace{\stretch{1}}(18.115)

Using these we are able to show that our PDEs are similar to that of

\begin{aligned}{\left.{{w = D^2 w = D^4 w}}\right\vert}_{{z = 0, 1}} = 0,\end{aligned} \hspace{\stretch{1}}(18.117)

where D = {\partial {}}/{\partial {z}}. Using the trig solutions that fall out of this we were able to find the constraint

\begin{aligned}0 =\begin{vmatrix}\left( n^2 \pi^2 + k^2 \right)^2 + \sigma \left( n^2 \pi^2 + k^2\right) & -\mathcal{R} k^2 \\ -1 & n^2 \pi^2 + k^2 + \text{P}_r \sigma\end{vmatrix},\end{aligned} \hspace{\stretch{1}}(18.118)

which for \sigma = 0, this gives us the critical value for the Rayleigh number

\begin{aligned}\mathcal{R} = \frac{(k^2 + n^2 \pi^2)^3}{k^2},\end{aligned} \hspace{\stretch{1}}(18.119)

which is the boundary for thermal stability or instability.

The end result was a lot of manipulation for which we didn’t do any sort of applied problems. It looks like a theory that requires a lot of study to do anything useful with, so my expectation is that it won’t be covered in detail on the exam. Having some problems to know why we spent two days on it in class would have been nice.


[1] L.D. Landau and E.M. Lifshitz. A Course in Theoretical Physics-Fluid Mechanics. Pergamon Press Ltd., 1987.

[2] S. Granger. Fluid Mechanics. Dover, New York, 1995.

[3] Wikipedia. Curvature — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 25-April-2012].

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

[5] Wikipedia. Blasius boundary layer — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 28-March-2012].

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , | Leave a Comment »

Spin down of coffee in a bottomless cup.

Posted by peeterjoot on April 25, 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.)]


Here’s a variation of a problem outlined in section 2 of [1], which looked at the time evolution of fluid with initial rotational motion, after the (cylindrical) rotation driver stops, later describing this as the spin down of a cup of tea. I’ll work the problem in more detail than in the text, and also make two refinements.

  • I drink coffee and not tea.
  • I stir my coffee in the interior of the cup and not on the outer edge.

Because of the second point I’ll model my stir stick as a rotating cylinder in the cup and not by somebody spinning the cup itself to stir the tea. This only changes the solution for the steady state part of the problem.


We’ll work in cylindrical coordinates following the conventions of figure (1).

Figure 1: Fluid flow in nested cylinders.


We’ll assume a solution that with velocity azimuthal in direction, and both pressure and velocity that are only radially dependent.

\begin{aligned}\mathbf{u} = u(r) \hat{\boldsymbol{\phi}}.\end{aligned} \hspace{\stretch{1}}(2.1)

\begin{aligned}p = p(r)\end{aligned} \hspace{\stretch{1}}(2.2)

Let’s first verify that this meets the non-compressible condition that eliminates the \mu \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{u}) term from Navier-Stokes

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u}&=\left(\hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z\right) \cdot \left(u \hat{\boldsymbol{\phi}}\right) \\ &=\hat{\boldsymbol{\phi}} \cdot\left(\hat{\mathbf{r}} \partial_r u + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi u + \hat{\mathbf{z}} \partial_z u\right) +u\left(\hat{\mathbf{r}} \cdot \partial_r \hat{\boldsymbol{\phi}} + \frac{\hat{\boldsymbol{\phi}}}{r} \cdot \partial_\phi \hat{\boldsymbol{\phi}} + \hat{\mathbf{z}} \cdot \partial_z \hat{\boldsymbol{\phi}}\right) \\ &=\hat{\boldsymbol{\phi}} \cdot \hat{\mathbf{r}} \partial_r u +u\frac{\hat{\boldsymbol{\phi}}}{r} \cdot \left(-\hat{\mathbf{r}}\right) \\ &= 0.\end{aligned}

Good. Now let’s express each of the terms of Navier-Stokes in cylindrical form. Our time dependence is

\begin{aligned}\rho \partial_t u(r, t) \hat{\boldsymbol{\phi}}=\rho \hat{\boldsymbol{\phi}} \partial_t u.\end{aligned} \hspace{\stretch{1}}(2.3)

Our inertial term is

\begin{aligned}\begin{aligned}\rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\frac{\rho u}{r} \partial_\phi (u \hat{\boldsymbol{\phi}}) \\ &=\frac{\rho u^2}{r} (-\hat{\mathbf{r}}).\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.4)

Our pressure term is

\begin{aligned}-\boldsymbol{\nabla} p=-\hat{\mathbf{r}} \partial_r p,\end{aligned} \hspace{\stretch{1}}(2.5)

and our Laplacian term is

\begin{aligned}\begin{aligned}\mu \boldsymbol{\nabla}^2 \mathbf{u}&=\mu \left( \frac{1}{{r}} \partial_r ( r \partial_r) + \frac{1}{{r^2}} \partial_{\phi\phi} + \partial_{z z}\right)u(r) \hat{\boldsymbol{\phi}} \\ &=\mu \left( \frac{\hat{\boldsymbol{\phi}}}{r} \partial_r ( r \partial_r u) + \frac{-\hat{\mathbf{r}} u}{r^2} \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.6)

Putting things together, we find that Navier-Stokes takes the form

\begin{aligned}\rho \hat{\boldsymbol{\phi}} \partial_t u+\frac{\rho u^2}{r} (-\hat{\mathbf{r}})=-\hat{\mathbf{r}} \partial_r p+\mu \left( \frac{\hat{\boldsymbol{\phi}}}{r} \partial_r ( r \partial_r u) + \frac{-\hat{\boldsymbol{\phi}} u}{r^2} \right),\end{aligned} \hspace{\stretch{1}}(2.7)

which nicely splits into an separate equations for the \hat{\boldsymbol{\phi}} and \hat{\mathbf{r}} directions respectively

\begin{aligned}\frac{1}{{\nu}} \partial_t u=\frac{1}{r} \partial_r ( r \partial_r u) - \frac{u}{r^2}\end{aligned} \hspace{\stretch{1}}(2.8a)

\begin{aligned}\frac{\rho u^2}{r}=\partial_r p.\end{aligned} \hspace{\stretch{1}}(2.8b)

Steady state solution

Before t = 0 we seek the steady state, the solution of

\begin{aligned}r \partial_r ( r \partial_r u) - u = 0.\end{aligned} \hspace{\stretch{1}}(2.9)

We’ve seen that

\begin{aligned}u(r) = A r + \frac{B}{r}\end{aligned} \hspace{\stretch{1}}(2.10)

is the general solution, and can now fit this to the boundary value constraints. For the interior portion of the cup we have

\begin{aligned}{\left.{{A r + \frac{B}{r}}}\right\vert}_{{r = 0}} = 0\end{aligned} \hspace{\stretch{1}}(2.11)

so B = 0 is required. For the interface of the “stir-stick” (moving fast enough that we can consider it having a cylindrical effect) at r = R_1 we have

\begin{aligned}A R_1 = \Omega R_1,\end{aligned} \hspace{\stretch{1}}(2.12)

so the interior portion of our steady state coffee velocity is just

\begin{aligned}\mathbf{u} = \Omega r \hat{\boldsymbol{\phi}}.\end{aligned} \hspace{\stretch{1}}(2.13)

Between the cup edge and the stir-stick we have to solve

\begin{aligned}A R_1 + \frac{B}{R_1} &= \Omega R_1 \\ A R_2 + \frac{B}{R_2} &= 0,\end{aligned} \hspace{\stretch{1}}(2.14)


\begin{aligned}A R_1^2 + B &= \Omega R_1^2 \\ A R_2^2 + B &= 0.\end{aligned} \hspace{\stretch{1}}(2.16)

Subtracting we find

\begin{aligned}A = -\frac{\Omega R_1^2}{R_2^2 - R_1^2}\end{aligned} \hspace{\stretch{1}}(2.18a)

\begin{aligned}B = \frac{\Omega R_1^2 R_2^2}{R_2^2 - R_1^2},\end{aligned} \hspace{\stretch{1}}(2.18b)

so our steady state coffee flow is

\begin{aligned}\mathbf{u} =\left\{\begin{array}{l l}\Omega r \hat{\boldsymbol{\phi}}& \quad \mbox{latex r \in [0, R_1]$} \\ \frac{\Omega R_1^2}{R_2^2 – R_1^2} \left( \frac{R_2^2}{r} -r \right)\hat{\boldsymbol{\phi}}& \quad \mbox{r \in [R_1, R_2]} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.19)$

Time evolution.

We can use a separation of variables technique with u(r, t) = R(r) T(t) to find

\begin{aligned}\frac{1}{{\nu}} \frac{T'}{T} = \frac{1}{{R}} \left( \frac{1}{r} \partial_r ( r \partial_r R) - \frac{R}{r^2}\right)= -\lambda^2,\end{aligned} \hspace{\stretch{1}}(2.20)

which gives us

\begin{aligned}T \propto e^{-\lambda^2 \nu t},\end{aligned} \hspace{\stretch{1}}(2.21)

and R specified by

\begin{aligned}0 = r^2 \frac{d^2 R}{dr^2} + r \frac{d R}{dr} + R \left( r^2 \lambda^2 - 1 \right).\end{aligned} \hspace{\stretch{1}}(2.22)

Checking [2] (9.1.1) we see that this can be put into the standard form of the Bessel equation if we eliminate the \lambda term. We can do that writing z = r \lambda, \mathcal{R}(z) = R(z/\lambda) and noting that r d/dr = z d/dz and r^2 d^2/dr^2 = z^2 d^2/dz^2, which gives us

\begin{aligned}0 = z^2 \frac{d^2 \mathcal{R}}{dr^2} + z \frac{d \mathcal{R}}{dr} + \mathcal{R} \left( z^2 - 1 \right).\end{aligned} \hspace{\stretch{1}}(2.23)

The solutions are

\begin{aligned}\mathcal{R}(z) = J_{\pm 1}(z), Y_{\pm 1}(z).\end{aligned} \hspace{\stretch{1}}(2.24)

From (9.1.5) of the handbook we see that the plus and minus variations are linearly dependent since J_{-1}(z) = -J_1(z) and Y_{-1}(z) = -Y_1(z), and from (9.1.8) that Y_1(z) is infinite at the origin, so our general solution has to be of the form

\begin{aligned}\mathbf{u}(r, t) = \hat{\boldsymbol{\phi}} \sum_\lambda c_\lambda e^{-\lambda^2 \nu t} J_{1}(r \lambda).\end{aligned} \hspace{\stretch{1}}(2.25)

In the text, I see that the transformation \lambda \rightarrow \lambda/a (where a was the radius of the cup) is made so that the Bessel function parameter was dimensionless. We can do that too but write

\begin{aligned}\mathbf{u}(r, t) = \hat{\boldsymbol{\phi}} \sum_\lambda c_\lambda e^{-\frac{\lambda^2}{R_2^2} \nu t} J_{1}\left(\lambda \frac{r}{R_2}\right).\end{aligned} \hspace{\stretch{1}}(2.26)

Our boundary value constraint is that we require this to match 2.19 at t = 0. Let’s write R_2 = R, R_1 = a R, z = r/R, so that we are working in the unit circle with z \in [0, 1]. Our boundary problem can now be expressed as

\begin{aligned}\frac{1}{{\Omega R}} \sum_\lambda c_\lambda J_{1}\left(\lambda z\right)=\left\{\begin{array}{l l}z & \quad \mbox{latex z \in [0, a]$} \\ \frac{1}{\frac{R^2}{a^2} – 1} \left( \frac{1}{{z}} – z\right)& \quad \mbox{z \in [a, 1]} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.27)$

Let’s pull the \Omega R factor into c_\lambda and state the problem to be solved as

\begin{aligned}\mathbf{u}(r, t) = \Omega R \hat{\boldsymbol{\phi}} \sum_{i=1}^n c_i e^{-\frac{\lambda_i^2}{R^2} \nu t} J_{1}\left(\lambda_i \frac{r}{R}\right)\end{aligned} \hspace{\stretch{1}}(2.28a)

\begin{aligned}\sum_{i = 1}^n c_i J_{1}\left(\lambda_i z\right) = \phi(z)\end{aligned} \hspace{\stretch{1}}(2.28b)

\begin{aligned}\phi(z) = \left\{\begin{array}{l l}z & \quad \mbox{latex z \in [0, a]$} \\ \frac{a^2}{1 – a^2} \left( \frac{1}{{z}} – z\right)& \quad \mbox{z \in [a, 1]} \\ \end{array}\right..\end{aligned} \hspace{\stretch{1}}(2.28c)$

Looking at section 2.7 of [3] it appears the solutions for c_i can be obtained from

\begin{aligned}c_i = \frac{\int_0^1 r\phi(z) J_1(\lambda_i z) dz}{\int_0^1 r J_1^2(\lambda_i z) dz},\end{aligned} \hspace{\stretch{1}}(2.29)

where \lambda_i are the zeros of J_1.

To get a feel for these, a plot of the first few of these fitting functions is shown in figure (2).

Figure 2: First four zero crossing Bessel functions.


Using Mathematica in bottomlessCoffee.cdf, these coefficients were calculated for a = 0.6. The n = 1, 3, 5 approximations to the fitting function are plotted with a comparison to the steady state velocity profile in figure (3).

Figure 3: Bessel function fitting for the steady state velocity profile for n = 1, 3, 5.


As indicated in the text, the spin down is way too slow to match reality (this can be seen visually in the worksheet by animating it).


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

[2] M. Abramowitz and I.A. Stegun. {\em Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

[3] H. Sagan. Boundary and eigenvalue problems in mathematical physics. Dover Pubns, 1989.

Posted in Math and Physics Learning. | Tagged: , , , , , | Leave a Comment »

Flow between infinite moving inner cylinder and outer cylinder (3D visualizations).

Posted by peeterjoot on April 16, 2012

[Click here for a PDF of this post with nicer formatting.]

We previously looked at the problem of an infinite cylinder of radius R_1 is moving with velocity v parallel to its axis. It is places inside another cylinder of radius R_2. The axes of the two cylinders coincide. The fluid is incompressible, with viscosity \mu and density \rho, the flow is assumed to be stationary, and no external pressure gradient is applied, and found the solution to have the form

\begin{aligned}\boxed{w = \frac{G}{4 \mu} (R_2^2 - r^2) + \left(v - \frac{G}{4 \mu} (R_2^2 - R_1^2)\right)\frac{\ln r/R_2}{\ln R_1/R_2}.}\end{aligned} \hspace{\stretch{1}}(3.27)

This system with such simple geometry lends itself well to 3D visualization of the velocity profiles

An animation of this can be found at The Mathematica notebook that generated this is available at phy454/twoCylinders3D.cdf. The 3D visualization portion of that notebook has also been made available as part of the wolfram demonstrations project.

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »

Flow between infinite moving inner cylinder and outer cylinder.

Posted by peeterjoot on April 14, 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.)]


A problem from the 2009 phy1530 final that appears appropriate for phy454 exam prep.


An infinite cylinder of radius R_1 is moving with velocity v parallel to its axis. It is places inside another cylinder of radius R_2. The axes of the two cylinders coincide. The fluid is incompressible, with viscosity \mu and density \rho, the flow is assumed to be stationary, and no external pressure gradient is applied.

  • Find and sketch the velocity field of the fluid between the cylinders.
  • Find the friction force per unit length acting on each cylinder.
  • Find and sketch the pressure field of the liquid.
  • If an external pressure gradient is present, how do you think your answer will change? Sketch your expectation for the velocity and pressure in this case.


Part 1,3. Velocity and pressure

Let’s start with the illustration of figure (1) to fix coordinates.

Figure 1: Coordinates for flow between two cylinders.


We’ll assume that we can find a solution of the following form


\begin{aligned}\mathbf{u} = w(r) \hat{\mathbf{z}}\end{aligned} \hspace{\stretch{1}}(3.1a)

\begin{aligned}p = p(r).\end{aligned} \hspace{\stretch{1}}(3.1b)


We’ll also work in cylindrical coordinates where our gradient is

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi + \hat{\mathbf{z}} \partial_z.\end{aligned} \hspace{\stretch{1}}(3.2)

Let’s look at the various terms of the Navier-Stokes equation. Our non-linear term is

\begin{aligned}\mathbf{u} \cdot \boldsymbol{\nabla} \mathbf{u} = w \partial_z ( w(r) \hat{\mathbf{z}} ) = 0,\end{aligned} \hspace{\stretch{1}}(3.3)

Our Laplacian term is

\begin{aligned}\mu \boldsymbol{\nabla}^2 \mathbf{u} &= \mu \left( \frac{1}{{r}} \partial_r ( r \partial_r ) + \frac{1}{{r^2}} \partial_{\phi\phi} + \partial_{z z} \right) w(r) \hat{\mathbf{z}} \\ &=\frac{\mu}{r} ( r w' )' \hat{\mathbf{z}}.\end{aligned}

Putting the pieces together we have

\begin{aligned}0 = - \hat{\mathbf{r}} p' + \frac{\mu}{r} (r w')' \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(3.4)

Decomposing these into one equation for each component we have

\begin{aligned}p' = 0,\end{aligned} \hspace{\stretch{1}}(3.5)


\begin{aligned}(r w')' = 0.\end{aligned} \hspace{\stretch{1}}(3.6)

The pressure can be trivially solved

\begin{aligned}p(r) = \text{constant},\end{aligned} \hspace{\stretch{1}}(3.7)

and for our velocity equation we get

\begin{aligned}r w' = A,\end{aligned} \hspace{\stretch{1}}(3.8)

Short of satisfying our boundary value constraints our velocity is

\begin{aligned}w = A \ln r + B.\end{aligned} \hspace{\stretch{1}}(3.9)

Our boundary value conditions are given by

\begin{aligned}w(R_2) &= 0 \\ w(R_1) &= v,\end{aligned} \hspace{\stretch{1}}(3.10)

so our integration constants are given by

\begin{aligned}0 &= A \ln R_2 + B \\ v &= A \ln R_1 + B.\end{aligned} \hspace{\stretch{1}}(3.12)

Taking differences we’ve got

\begin{aligned}v = A \ln( R_1/R_2 ).\end{aligned} \hspace{\stretch{1}}(3.14)

So our constants are


\begin{aligned}A = \frac{v}{ \ln( R_1/R_2 ) }\end{aligned} \hspace{\stretch{1}}(3.15a)

\begin{aligned}B = -\frac{v \ln R_2}{ \ln( R_1/R_2 ) },\end{aligned} \hspace{\stretch{1}}(3.15b)



\begin{aligned}\boxed{w(r) = \frac{ v \ln (r/R_2) }{ \ln( R_1/R_2 ) }.}\end{aligned} \hspace{\stretch{1}}(3.16)

A plot of this function can be found in figure (2), and the Mathematica notebook that generated this plot can be found in phy454/twoCylinders.cdf. That notebook has some slider controls that can be used interactively.

Figure 2: Velocity plot due to inner cylinder dragging fluid along with it.


Part 2. Friction force per unit length.

For the frictional force per unit area on the fluid by the inner cylinder we have

\begin{aligned}(\boldsymbol{\sigma} \cdot \hat{\mathbf{r}} ) \cdot \hat{\mathbf{z}}&=-p \hat{\mathbf{r}} \cdot \hat{\mathbf{z}} +2 \mu \frac{1}{{2}} \left( \frac{\partial {u_z}}{\partial {r}}+\not{{\frac{\partial {u_r}}{\partial {z}}}}\right) \\ &=\mu v \frac{\ln r }{\ln (R_1/R_2) }\end{aligned}

So the forces on the inner and outer cylinders for a strip of width \Delta z is


\begin{aligned}\text{frictional force on inner cylinder} = -2 \pi R_1 \Delta z \mu v \hat{\mathbf{z}} \frac{\ln R_1 }{\ln (R_1/R_2) }\end{aligned} \hspace{\stretch{1}}(3.17a)

\begin{aligned}\text{frictional force on inner cylinder} =2 \pi R_2 \Delta z \mu v \hat{\mathbf{z}} \frac{\ln R_2 }{\ln (R_1/R_2) }\end{aligned} \hspace{\stretch{1}}(3.17b)


Part 4. With external pressure gradient.

With an external pressure gradient imposed we expect a superposition of a parabolic flow profile with what we’ve calculated above. With

\begin{aligned}G = - \frac{dp}{dz},\end{aligned} \hspace{\stretch{1}}(3.18)

our Navier-Stokes equation will now take the form

\begin{aligned}0 = - \hat{\mathbf{r}} p' - (-G \hat{\mathbf{z}}) + \frac{\mu}{r} (r w')' \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(3.19)

We want to solve the LDE

\begin{aligned}- \frac{ G r}{\mu} =(r w')' = r w'' + w'\end{aligned} \hspace{\stretch{1}}(3.20)

The homogeneous portion of this equation

\begin{aligned}(r w')' = 0,\end{aligned} \hspace{\stretch{1}}(3.21)

we have already solved finding w = C \ln r + D. It looks reasonable to try a polynomial solution for the specific solution. Let’s try a second order polynomial


\begin{aligned}w = A r^2 + B r \end{aligned} \hspace{\stretch{1}}(3.22a)

\begin{aligned}w' = 2 A r + B\end{aligned} \hspace{\stretch{1}}(3.22b)

\begin{aligned}w'' = 2 A.\end{aligned} \hspace{\stretch{1}}(3.22c)


We need

\begin{aligned}-\frac{G r}{\mu} = 2 A r + 2 A r + B,\end{aligned} \hspace{\stretch{1}}(3.23)

So B = 0 and 4 A = -G/\mu, and our general solution has the form

\begin{aligned}w = -\frac{G}{4 \mu} r^2 + C \ln r + D.\end{aligned} \hspace{\stretch{1}}(3.24)

requiring just the boundary condition fitting. Let’s tweak the constants slightly, writing

\begin{aligned}w = \frac{G}{4 \mu} (R_2^2 - r^2) + C \ln r/R_2 + D,\end{aligned} \hspace{\stretch{1}}(3.25)

so that D = 0 falls out of the w(R_2) = 0 constraint. Our last integration constant is then determined by the solution of

\begin{aligned}v = \frac{G}{4 \mu} (R_2^2 - R_1^2) + C \ln R_1/R_2.\end{aligned} \hspace{\stretch{1}}(3.26)


\begin{aligned}\boxed{w = \frac{G}{4 \mu} (R_2^2 - r^2) + \left(v - \frac{G}{4 \mu} (R_2^2 - R_1^2)\right)\frac{\ln r/R_2}{\ln R_1/R_2}.}\end{aligned} \hspace{\stretch{1}}(3.27)

A plot of this, with a pressure gradient small enough that we still see the logarithmic profile is shown in figure (3). An animation of this with different values for R_1, v, and G/4\mu is available on, but the Mathematica notebook above can also be used.

Figure 3: Pressure gradient added.


Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

Couette flow of a viscous incompressible fluid (visualizing it).

Posted by peeterjoot on April 12, 2012

[Click here for a PDF of this post(and previous).]

In Couette flow of a viscous incompressible fluid, I looked at solving this problem. As a followup, here are some visualizations of the velocity field as the outer cylinder angular velocity is changed.

The Mathematica workbook that was used to generate the figures is also available for more direct play, and includes slider controls for the cylinder radii and angular velocities.

Posted in Math and Physics Learning. | Tagged: , , , , , , | Leave a Comment »

Couette flow of a viscous incompressible fluid.

Posted by peeterjoot on April 9, 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.)]


A problem from this years phy1530 problem set 2 that appears appropriate for phy454 exam prep.


Consider the steady flow between two long cylinders of radii R_1 and R_2, R_1 > R_1, rotating about their axes with angular velocities \Omega_1, \Omega_2. Look for a solution of the form, where \hat{\boldsymbol{\phi}} is a unit vector along the azimuthal direction:


\begin{aligned}\mathbf{u} = v(r) \hat{\boldsymbol{\phi}}\end{aligned} \hspace{\stretch{1}}(2.1a)

\begin{aligned}p = p(r).\end{aligned} \hspace{\stretch{1}}(2.1b)


  • Write out the Navier-Stokes equations and find differential equations for v(r) and p(r). You should find that these equations have relatively simple solutions, i.e.,

    \begin{aligned}v(r) = a r + \frac{b}{r}.\end{aligned} \hspace{\stretch{1}}(2.2)

  • Fix gthe constants a and b from the boundary conditions. Determine the pressure p(r).
  • Compute the friction forces that the fluid exerts on the cylinders, and compute the torque on each cylinder. Show that the total torque on the fluid is zero (as must be the case).

This is also a problem that I recall was outlined in section 2 from [1]. Some of the instabilities that are mentioned in the text are nicely illustrated in [2].

We illustrate our system in figure (1).

Figure 1: Coutette flow configuration


Solution: Part 1. Navier-Stokes and resulting differential equations.

Navier-Stokes for steady state incompressible flow has the form


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

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0.\end{aligned} \hspace{\stretch{1}}(3.3b)


where the gradient has the form

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi.\end{aligned} \hspace{\stretch{1}}(3.4)

Let’s first verify that the incompressible condition 3.3b is satisfied for the presumed form of the solution we seek. We have

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} &=\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) \cdot (v(r) \hat{\boldsymbol{\phi}}(\phi) ) \\ &=(\hat{\mathbf{r}} \cdot \hat{\boldsymbol{\phi}}) v' + \frac{\hat{\boldsymbol{\phi}}^2}{r} \partial_\phi v(r)+ \frac{v(r) \hat{\boldsymbol{\phi}}}{r} \cdot \partial_\phi \hat{\boldsymbol{\phi}} \\ &= \frac{v(r) \hat{\boldsymbol{\phi}}}{r} \cdot (-\hat{\mathbf{r}}) \\ &= 0\end{aligned}

Good. Now let’s write out the terms of the momentum conservation equation 3.3a. We’ve got

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\frac{ v}{r} \partial_\phi ( v \hat{\boldsymbol{\phi}} ) \\ &=-\frac{ v^2 \hat{\mathbf{r}}}{r},\end{aligned}


\begin{aligned}-\frac{1}{{\rho}} \boldsymbol{\nabla} p&=-\frac{1}{{\rho}} \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) p(r) \\ &=-\frac{\hat{\mathbf{r}} p'}{\rho},\end{aligned}


\begin{aligned}\nu \boldsymbol{\nabla}^2 \mathbf{u}&=\nu \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) \cdot\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) (v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \partial_{rr} + \frac{1}{{r^2}} \partial_{\phi\phi}+ \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \cdot (\hat{\mathbf{r}} \partial_r)\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \partial_{rr} + \frac{1}{{r^2}} \partial_{\phi\phi}+ \frac{1}{r} \partial_r\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \frac{1}{{r}} \partial_{r} (r \partial_r) + \frac{1}{{r^2}} \partial_{\phi\phi}\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \frac{1}{{r}} (r v')' - \frac{v}{r^2} \right)\hat{\boldsymbol{\phi}}\end{aligned}

Equating \hat{\mathbf{r}} and \hat{\boldsymbol{\phi}} components we have two equations to solve


\begin{aligned}r (r v')' - v = 0\end{aligned} \hspace{\stretch{1}}(3.5a)

\begin{aligned}p' = \frac{\rho v^2}{r}.\end{aligned} \hspace{\stretch{1}}(3.5b)


Expanding out our velocity equation we have

\begin{aligned}r^2 v'' + r v' - v = 0,\end{aligned} \hspace{\stretch{1}}(3.6)

for which we’ve been told to expect that 2.2 is a solution (and it has the two integration constants we require for a solution to a homogeneous equation of this form). Let’s verify that we’ve computed the correct differential equation for the problem by trying this solution

\begin{aligned}r^2 v'' + r v' - v &=r^2 \left( a -\frac{b}{r^2} \right)' + r \left( a -\frac{b}{r^2} \right) - a r - \frac{b}{r} \\ &=r^2 \frac{2 b}{r^3} + \not{{a r}} - \frac{b}{r} - \not{{a r}} - \frac{b}{r} \\ &=\frac{2 b}{r} - \frac{2 b}{r} \\ &= 0.\end{aligned}

Given the velocity, we can now determine the pressure up to a constant

\begin{aligned}p' &= \frac{\rho}{r} \left( a r + \frac{b}{r} \right)^2 \\ &= \frac{\rho}{r} \left( a^2 r^2 + \frac{b^2}{r^2} + 2 a b \right) \\ &= \rho \left( a^2 r + \frac{b^2}{r^3} + 2 \frac{a b}{r} \right)\end{aligned}


\begin{aligned}p_r -p_0= \rho \left( \frac{1}{{2}} a^2 r^2 - \frac{b^2}{2 r^2} + 2 a b \ln r \right)\end{aligned} \hspace{\stretch{1}}(3.7)

Solution: Part 2. Fixing the constants.

To determine our integration constants we recall that velocity associated with a radial position \mathbf{x} = r \hat{\mathbf{r}} in cylindrical coordinates takes the form

\begin{aligned}\frac{\mathbf{x}}{dt} = \dot{r} \hat{\mathbf{r}} + r \hat{\boldsymbol{\phi}} \dot{\phi},\end{aligned} \hspace{\stretch{1}}(4.8)

where \dot{\phi} is the angular velocity. The cylinder walls therefore have the velocity

\begin{aligned}v = r \dot{\phi},\end{aligned} \hspace{\stretch{1}}(4.9)

so our boundary conditions (given a no-slip assumption for the fluids) are

\begin{aligned}v(R_1) &= R_1 \Omega_1 \\ v(R_2) &= R_2 \Omega_2.\end{aligned} \hspace{\stretch{1}}(4.10)

This gives us a pair of equations to solve for a and b

\begin{aligned}R_1 \Omega_1 &= a R_1 + \frac{b}{R_1} \\ R_2 \Omega_2 &= a R_2 + \frac{b}{R_2}.\end{aligned} \hspace{\stretch{1}}(4.12)

Multipling each by R_1 and R_2 respectively gives us

\begin{aligned}b = R_1^2 (\Omega_1 - a) = R_2^2 (\Omega_2 - a).\end{aligned} \hspace{\stretch{1}}(4.14)

Rearranging for a we find

\begin{aligned}R_1^2 \Omega_1 - R_2^2 \Omega_2 = (R_1^2 - R_2^2) a,\end{aligned} \hspace{\stretch{1}}(4.15)


\begin{aligned}a = \frac{ R_2^2 \Omega_2 - R_1^2 \Omega_1}{R_2^2 - R_1^2}.\end{aligned} \hspace{\stretch{1}}(4.16)

For b we have

\begin{aligned}b &= R_1^2 (\Omega_1 - a) \\ &=\frac{R_1^2 }{R_2^2 - R_1^2}(\Omega_1 ( R_2^2 - \not{{R_1^2}}) - R_2^2 \Omega_2 + \not{{R_1^2 \Omega_1}}),\end{aligned}


\begin{aligned}b = \frac{R_1^2 R_2^2}{R_2^2 - R_1^2} (\Omega_1 -\Omega_2).\end{aligned} \hspace{\stretch{1}}(4.17)

This gives us


\begin{aligned}v(r) = \frac{1}{{R_2^2 - R_1^2}}\left(\left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right) r+\frac{R_1^2 R_2^2}{r} (\Omega_1 -\Omega_2)\right)\end{aligned} \hspace{\stretch{1}}(4.18a)

\begin{aligned}\begin{aligned}p_r -&p_0= \frac{\rho }{(R_2^2 - R_1^2)^2} \times \\ &\left( \frac{1}{{2}} \left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right)^2r^2 -\frac{R_1^4 R_2^4}{2 r^2} (\Omega_1 - \Omega_2)^2+ 2 \left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right) R_1^4 R_2^4 (\Omega_1 - \Omega_2)^2 \ln r\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.18b)


Solution: Part 3. Friction and torque.

We can expand out the identity for the traction vector

\begin{aligned}\mathbf{t} = \mathbf{e}_i \sigma_{ij} n_j= -p \hat{\mathbf{n}} + \mu \left( 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u} + \hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})\right),\end{aligned} \hspace{\stretch{1}}(5.19)

in cylindrical coordinates and find


\begin{aligned}\sigma_{rr}=-p + 2 \mu \not{{\frac{\partial {u_r}}{\partial {r}}}}\end{aligned} \hspace{\stretch{1}}(5.20a)

\begin{aligned}\sigma_{\phi \phi}=-p + 2 \mu \left(\frac{1}{{r}}\not{{\frac{\partial {u_\phi}}{\partial {\phi}}}} + \not{{\frac{u_r}{r}}}\right)\end{aligned} \hspace{\stretch{1}}(5.20b)

\begin{aligned}\sigma_{z z}=-p + 2 \mu \not{{\frac{\partial {u_z}}{\partial {z}}}}\end{aligned} \hspace{\stretch{1}}(5.20c)

\begin{aligned}\sigma_{r \phi}=\mu \left( \frac{\partial {u_\phi}}{\partial {r}}+\frac{1}{{r}} \not{{\frac{\partial {u_r}}{\partial {\phi}}}}- \frac{u_\phi}{r}\right)\end{aligned} \hspace{\stretch{1}}(5.20d)

\begin{aligned}\sigma_{\phi z}=\mu \left(\frac{1}{r} \not{{\frac{\partial {u_z}}{\partial {\phi}}}} + \not{{\frac{\partial {u_\phi}}{\partial {z}}}}\right)\end{aligned} \hspace{\stretch{1}}(5.20e)

\begin{aligned}\sigma_{z r}=\mu \left(\not{{\frac{\partial {u_r}}{\partial {z}}}}+ \not{{\frac{\partial {u_z}}{\partial {r}}}}\right),\end{aligned} \hspace{\stretch{1}}(5.20f)


so we have


\begin{aligned}\sigma_{rr} = \sigma_{\phi \phi} = \sigma_{z z} = -p \end{aligned} \hspace{\stretch{1}}(5.21a)

\begin{aligned}\sigma_{\phi z} = \sigma_{z r} = 0\end{aligned} \hspace{\stretch{1}}(5.21b)

\begin{aligned}\sigma_{r \phi} = \mu \left( \frac{\partial {u_\phi}}{\partial {r}} - \frac{u_\phi}{r} \right)\end{aligned} \hspace{\stretch{1}}(5.21c)


We want to expand the last of these

\begin{aligned}\sigma_{r \phi} &= \mu \left( \frac{\partial {u_\phi}}{\partial {r}} - \frac{u_\phi}{r} \right) \\ &= \mu \left( a r + \frac{b}{r}\right)' \\ &= \mu \left( a - \frac{b}{r^2}\right).\end{aligned}

So the traction vector, our force per unit area on the fluid at the inner surface (where the normal is \hat{\mathbf{r}}), is

\begin{aligned}\mathbf{t}_1 = -p \hat{\mathbf{r}} + \mu \left( a - \frac{b}{r^2} \right) \hat{\boldsymbol{\phi}},\end{aligned} \hspace{\stretch{1}}(5.22)

and our torque per unit area from the inner cylinder on the fluid is thus

\begin{aligned}\boldsymbol{\tau}_1 = r \hat{\mathbf{r}} \times \mathbf{t} = r \mu \left( a - \frac{b}{r^2} \right) \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(5.23)

Observing that our stress tensors flip sign for an inwards normal, our torque per unit area from the outer cylinder is

\begin{aligned}\boldsymbol{\tau}_2 = r \hat{\mathbf{r}} \times (-\mathbf{t}_1) = -r \mu \left( a - \frac{b}{r^2} \right) \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(5.24)

For the complete torque on the fluid due to a strip of width \Delta z the magnitudes of the total torque from each cylinder are respectively

\begin{aligned}\tau_1 = 2 \pi r^2 \Delta z \mu \left( a - \frac{b}{r^2} \right)\end{aligned} \hspace{\stretch{1}}(5.25)

\begin{aligned}\tau_2 = - 2 \pi r^2 \Delta z \mu \left( a - \frac{b}{r^2} \right)\end{aligned} \hspace{\stretch{1}}(5.26)

As expected these torques on the fluids sum to zero

\begin{aligned}\tau_2 + \tau_1 = 0.\end{aligned} \hspace{\stretch{1}}(5.27)

Evaluating these at R_1 and R_2 respectively gives us the torques on the fluid by the cylinders, so inverting these provides the torques on the cylinders by the fluid

\begin{aligned}\text{torque on inner cylinder (1) by the fluid} = -(2 \pi R_1^2) \Delta z \mu \left( a - \frac{b}{R_1^2} \right).\end{aligned} \hspace{\stretch{1}}(5.28)

For the outer cylinder the total torque on a strip of width \Delta z is

\begin{aligned}\text{torque on outer cylinder (2) by the fluid} = (2 \pi R_2^2) \Delta z \mu \left( a - \frac{b}{R_2^2} \right).\end{aligned} \hspace{\stretch{1}}(5.29)


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

[2] Wikipedia. Taylor-couette flow — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 10-April-2012].\%E2\%80\%93Couette_flow&oldid=483281707.

Posted in Math and Physics Learning. | Tagged: , , , , , | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 22: Thermal stability. Taught by Prof. K. Das.

Posted by peeterjoot on April 5, 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.)]


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

Review. Rayleigh Benard Problem

Reading: section 9.3 from [1].

Illustrated in figure (1) is the heated channel we’ve been discussing

Channel with heat applied to the base


We’ll take initial conditions

\begin{aligned}\frac{\partial {\mathbf{u}}}{\partial {t}} &= 0 \\ \frac{\partial {T}}{\partial {t}} &= 0\end{aligned} \hspace{\stretch{1}}(2.1)

\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = - \boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} + \rho \mathbf{g}.\end{aligned} \hspace{\stretch{1}}(2.3)

Our energy equation is

\begin{aligned}\frac{\partial {T}}{\partial {t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) T = \kappa \boldsymbol{\nabla}^2 T.\end{aligned} \hspace{\stretch{1}}(2.4)

We have this \mathbf{u} \cdot \boldsymbol{\nabla} term because our heat can be carried from one place to the other, due to the fluid motion. We’d not have this convective term for heat dissipation in solids because elements of a solid are not moving around in the bulk.

We’ll also use

\begin{aligned}\boldsymbol{\nabla} p_s = \rho_s \mathbf{g}\end{aligned} \hspace{\stretch{1}}(2.5)

In the steady (base) state we have

\begin{aligned}0 = \kappa \boldsymbol{\nabla}^2 T = \kappa \left( \frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}+\frac{\partial^2 {{}}}{\partial {{z}}^2} \right) T,\end{aligned} \hspace{\stretch{1}}(2.6)

but since we are only considering spatial variation with z we have

\begin{aligned}\kappa \frac{\partial^2 {{}}}{\partial {{z}}^2} T_s = 0,\end{aligned} \hspace{\stretch{1}}(2.7)

with solution

\begin{aligned}T_s = T_0 - \frac{\Delta T}{d} z.\end{aligned} \hspace{\stretch{1}}(2.8)

We found that after application of the perturbation

\begin{aligned}\mathbf{u} &\rightarrow 0 + \delta \mathbf{u} \\ p &\rightarrow p_s + \delta p \\ \rho &\rightarrow p_s + \delta \rho \\ T &\rightarrow p_s + \delta T\end{aligned} \hspace{\stretch{1}}(2.9)

to the base state equations, our perturbed Navier-Stokes equation was

\begin{aligned}\boldsymbol{\nabla}^2 \left( \frac{\partial {}}{\partial {t}} - \nu \boldsymbol{\nabla}^2 \right) \delta w = g \alpha \left(\frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right)\delta T.\end{aligned} \hspace{\stretch{1}}(2.13)

Application of the perturbation to the energy equation.

\begin{aligned}\frac{\partial {T_s + \delta T}}{\partial {t}} + (\delta \mathbf{u} \cdot \boldsymbol{\nabla}) (T_s + \delta T) = \kappa \boldsymbol{\nabla}^2 (T_s + \delta T)\end{aligned} \hspace{\stretch{1}}(3.14)

We’ve got

\begin{aligned}\frac{\partial {T_s}}{\partial {t}} = 0.\end{aligned} \hspace{\stretch{1}}(3.15)

Using this, and 2.7, and neglecting any terms of second order smallness we have

\begin{aligned}\boxed{\frac{\partial {\delta T}}{\partial {t}} + \delta \mathbf{u} \cdot \boldsymbol{\nabla} T_s = \kappa \boldsymbol{\nabla}^2 \delta T.}\end{aligned} \hspace{\stretch{1}}(3.16)

We’d like to solve this and 2.13 simultaneously.

Non-dimensionalisation of the thermal velocity equation.

We’d like to scale

\begin{aligned}\begin{array}{l l}x,y,z & \quad \mbox{with latex d$} \\ t & \quad \mbox{with d^2/\nu} \\ \delta w & \quad \mbox{with \kappa/d} \\ \delta T & \quad \mbox{with \Delta T}\end{array}\end{aligned} \hspace{\stretch{1}}(3.17)$

Sanity check of dimensions

  1. viscosity dimensions

    \begin{aligned}[\nu] \sim \left[ \frac{1}{\text{T}} \right]/ \left[\frac{1}{{\text{L}^2}} \right] \sim \frac{\text{L}^2}{\text{T}},\end{aligned} \hspace{\stretch{1}}(3.18)

  2. thermal conductivity dimensionsSince [ \kappa \boldsymbol{\nabla}^2 ] \sim 1/\text{T}, we have

    \begin{aligned}[ \kappa ] \sim \frac{\text{L}^2}{\text{T}}\end{aligned} \hspace{\stretch{1}}(3.19)

  3. time scaling

    \begin{aligned}\left[\frac{d^2}{\nu}\right] \sim \frac{\text{L}^2}{\text{L}^2 \text{T}^{-1}} \sim \text{T}.\end{aligned} \hspace{\stretch{1}}(3.20)

  4. velocity scaling

    \begin{aligned}[\kappa/d] \sim \frac{\text{L}^2}{\text{T}} \frac{1}{{\text{L}}} \sim [\delta w]\end{aligned} \hspace{\stretch{1}}(3.21)

Looks like everything checks out.

Let’s apply this rescaling to our perturbed velocity equation 2.13

\begin{aligned}{\boldsymbol{\nabla}'}^2 \left( \frac{\nu}{d^2} \frac{\partial {}}{\partial {t'}} - \frac{\nu}{d^2} {\boldsymbol{\nabla}'}^2 \right)\frac{\kappa}{d}\delta w'=g \alpha \Delta T\left( \frac{\partial^2 {{}}}{\partial {{x'}}^2}+\frac{\partial^2 {{}}}{\partial {{y'}}^2}\right) \delta T'.\end{aligned} \hspace{\stretch{1}}(3.24)

Introducing the \textit{Rayleigh number}

\begin{aligned}\mathcal{R} = \frac{g \alpha \Delta T d^3}{\nu \kappa},\end{aligned} \hspace{\stretch{1}}(3.23)

and dropping primes, we have

\begin{aligned}\boxed{\boldsymbol{\nabla}^2 \left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right)\delta w=\mathcal{R} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right) \delta T.}\end{aligned} \hspace{\stretch{1}}(3.24)

My class notes original had \mathcal{R} with the value

\begin{aligned}\frac{g \Delta T d^2}{\nu \alpha},\end{aligned}

but performing this non-dimensionalization shows that this was either quoted incorrectly, or typed wrong in the heat of the moment. A check against the text shows (equation (9.27)), shows that 3.23 is correct.

Non-dimensionalization of the energy equation

Rescaling our energy equation 3.16 we find

\begin{aligned}\frac{\nu}{d^2} \frac{\partial {\delta T'}}{\partial {t'}} + \frac{\kappa}{d^2} \delta \mathbf{u}' \cdot \boldsymbol{\nabla}' T_s' = \frac{\kappa}{d^2} {\boldsymbol{\nabla}'}^2 \delta T'.\end{aligned} \hspace{\stretch{1}}(3.25)

Introducing the \textit{Prandtl number}

\begin{aligned}\text{P}_r = \frac{\nu}{\kappa},\end{aligned} \hspace{\stretch{1}}(3.26)

and dropping primes our non-dimensionalized energy equation takes the form

\begin{aligned}\boxed{\left( \text{P}_r\frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta T = \delta w.}\end{aligned} \hspace{\stretch{1}}(3.27)

Normal mode analysis.

We’ve got a pair of nasty looking coupled equations 3.24, and 3.27. Repeated so that we can see them together


\begin{aligned}\boldsymbol{\nabla}^2 \left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta w = \mathcal{R} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2} +\frac{\partial^2 {{}}}{\partial {{y}}^2} \right) \delta T\end{aligned} \hspace{\stretch{1}}(3.28a)

\begin{aligned}\left( \text{P}_r \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta T = \delta w,\end{aligned} \hspace{\stretch{1}}(3.42)


it’s clear that we can decouple these by inserting 3.42 into 3.28a. Doing that gives us a beastly 6th order spatial equation for the perturbed temperature

\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \left( \text{P}_r \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \boldsymbol{\nabla}^2 \delta T = \mathcal{R} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2} +\frac{\partial^2 {{}}}{\partial {{y}}^2} \right) \delta T.\end{aligned} \hspace{\stretch{1}}(3.29)

It’s pointed out in the text we have all the x and y derivatives coming together we can apply separation of variables with

\begin{aligned}\delta T = \Theta(z) f(x, y) e^{\sigma t},\end{aligned} \hspace{\stretch{1}}(3.30)

provided we introduce some restrictions on the form of f(x, y). Here \sigma (if real) is the growth rate. Applying the Laplacian to this assumed solution we find

\begin{aligned}\boldsymbol{\nabla}^2 \delta T = e^{\sigma t} \left( \Theta''(z) f(x, y) +\Theta(z) \boldsymbol{\nabla}_t^2 f(x, y) \right),\end{aligned} \hspace{\stretch{1}}(3.31)


\begin{aligned}\boldsymbol{\nabla}_t^2 = \frac{\partial^2 {{}}}{\partial {{x}}^2} +\frac{\partial^2 {{}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(3.32)

For 3.31 to be separable we require a constant proportionality

\begin{aligned}\boldsymbol{\nabla}_t^2 f(x, y) \propto f(x, y),\end{aligned} \hspace{\stretch{1}}(3.33)


\begin{aligned}\boldsymbol{\nabla}_t^2 f(x, y) \pm k^2 f(x, y) = 0.\end{aligned} \hspace{\stretch{1}}(3.34)

Picking + k^2 so that we don’t have hyperbolic solutions, f must have the form

\begin{aligned}f(x, y) = e^{i (k_x x + k_y y)},\end{aligned} \hspace{\stretch{1}}(3.35)


\begin{aligned}k^2 = k_x^2 + k_y^2.\end{aligned} \hspace{\stretch{1}}(3.36)

Our separation of variables function now takes the form

\begin{aligned}\delta T = \Theta(z) e^{ i ( k_1 x + k_2 y) + \sigma t}.\end{aligned} \hspace{\stretch{1}}(3.37)


\begin{aligned}D = \frac{\partial}{\partial z},\end{aligned} \hspace{\stretch{1}}(3.38)

our beastly equation to solve is then given by

\begin{aligned}0 &= \left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}_t^2 - D^2 \right) \left( \text{P}_r \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}_t^2 -D^2 \right) (\boldsymbol{\nabla}_t^2 -D^2)\Theta(z) e^{ i ( k_1 x + k_2 y) + \sigma t}- \mathcal{R} \boldsymbol{\nabla}_t^2 \Theta(z) e^{ i ( k_1 x + k_2 y) + \sigma t}\\ &=\Bigl( \left( \sigma + k^2 - D^2 \right) \left( \text{P}_r \sigma + k^2 - D^2 \right) \left( -k^2 -D^2 \right) + \mathcal{R} k^2 \Bigr) \Theta(z) e^{ i ( k_1 x + k_2 y) + \sigma t}\end{aligned}

This is now an equation for only \Theta(z)

\begin{aligned}0 = \Bigl( \left( \sigma + k^2 - D^2 \right) \left( \text{P}_r \sigma + k^2 - D^2 \right) \left( -k^2 -D^2 \right) + \mathcal{R} k^2 \Bigr) \Theta(z).\end{aligned} \hspace{\stretch{1}}(3.39)

Conceptually we have just a plain old LDE, and should we decide to expand this out we have something of the form

\begin{aligned}0 = \Bigl( \alpha D^6 + \beta D^4 + \gamma D^2 + \zeta \Bigr) \Theta(z).\end{aligned} \hspace{\stretch{1}}(3.40)

Our standard toolbox method to solve this is to assume a solution \Theta(z) = e^{m z} and compute the characteristic equation. We’d have to solve

\begin{aligned}0 = \alpha m^6 + \beta m^4 + \gamma m^2 + \zeta.\end{aligned} \hspace{\stretch{1}}(3.41)

Let’s back up a bit instead. Looking back to 3.42, it’s clear that we’ll have the same separable form for our perturbed velocity since we have

\begin{aligned}\left( \text{P}_r \sigma + k^2 - D^2 \right) \Theta(z) e^{i (\mathbf{k} \cdot \mathbf{x}) } e^{\sigma t} = \delta w,\end{aligned} \hspace{\stretch{1}}(3.42)


\begin{aligned}\mathbf{k} = (k_x, k_y, 0).\end{aligned} \hspace{\stretch{1}}(3.43)

Assuming a solution of the form

\begin{aligned}\delta w = w(z) e^{ i ( k_1 x + k_2 y) + \sigma t},\end{aligned} \hspace{\stretch{1}}(3.44)

our velocity is then fully specified in terms of the temperature, since we have

\begin{aligned}w(z) = \left( \text{P}_r \sigma + k^2 - D^2 \right) \Theta(z).\end{aligned} \hspace{\stretch{1}}(3.45)

Back to our coupled equations.

Having gleamed an idea what the form of our solutions is, we can simplify our original coupled system, writing


\begin{aligned}(D^2 - k^2) ( \sigma - (D^2 - k^2) ) w(z) = -\mathcal{R} k^2 \Theta(z)\end{aligned} \hspace{\stretch{1}}(3.46a)

\begin{aligned}(D^2 - k^2 -\text{P}_r \sigma ) \Theta(z) = -w(z).\end{aligned} \hspace{\stretch{1}}(3.46b)


Considering the boundary conditions, if the heating is even then at z = t = 0 we can’t have any variation with x and y, so can only have \Theta(0) = 0. Thus at the boundary, from 3.46a, we have

\begin{aligned}0 = {\left.{{(D^2 - k^2) ( \sigma - (D^2 - k^2) ) w(z)}}\right\vert}_{{z = 0}}.\end{aligned} \hspace{\stretch{1}}(3.47)

From the continuity equation \boldsymbol{\nabla} \cdot \mathbf{u} = 0, the text argues that we also have {\partial {\delta w}}/{\partial {z}} = 0 on the boundary, so that on that plane we also have

\begin{aligned}0 = {\left.{{w(z) = D w(z)}}\right\vert}_{{z = 0}}.\end{aligned} \hspace{\stretch{1}}(3.48)

Expanding out 3.47 then gives us

\begin{aligned}0 = {\left.{{(D^4 - D^2 ( 2 k^2 + \sigma ) + k^2(\sigma + k^2))w}}\right\vert}_{{z = 0}}.\end{aligned} \hspace{\stretch{1}}(3.50)


\begin{aligned}0 = {\left.{{(D^4 - D^2 ( 2 k^2 + \sigma ))w}}\right\vert}_{{z = 0}}.\end{aligned} \hspace{\stretch{1}}(3.50)

These boundary value constraints 3.48, and 3.50, plus the coupled system equations 3.46 are the complete problem to solve. To get a feel for the solution of this system, consider the system with the following simpler set of boundary value constraints

\begin{aligned}{\left.{{w = D^2 w = D^4 w}}\right\vert}_{{z = 0, 1}} = 0,\end{aligned} \hspace{\stretch{1}}(3.51)

which in the text is described as the artificial problem of thermal instability for boundaries that are stress free (FIXME: it’s not clear to me what that means without some thought … return to this). For such a system on the boundaries z = 0, 1 (noting that we are still in dimensionless quantities), we have solutions

\begin{aligned}w(z) &= A_n \sin\left( n \pi z \right) \\ \Theta(z) &= B_n \sin\left( n \pi z \right).\end{aligned} \hspace{\stretch{1}}(3.52)

Note that we have

\begin{aligned}D^2 - k^2 = \left( n \pi \right)^2 - k^2.\end{aligned} \hspace{\stretch{1}}(3.54)

Inserting 3.52 into our system 3.46, we have

\begin{aligned}\left( \left( n \pi \right)^2 - k^2 \right)\left( \sigma - \left( \left( n \pi \right)^2 - k^2 \right)\right) A_n \sin\left( n \pi z \right)+ \mathcal{R} k^2 B_n \sin\left( n \pi z \right)&= 0 \\ \left( \left( n \pi \right)^2 - k^2 - \text{P}_r \sigma \right) B_n \sin\left( n \pi z \right)+ A_n \sin\left( n \pi z \right)&= 0\end{aligned} \hspace{\stretch{1}}(3.55)

For any A_n, B_n, we must then have

\begin{aligned}0 = \begin{vmatrix}\left( n^2 \pi^2 - k^2 \right)^2 - \sigma \left( n^2 \pi^2 - k^2\right) & -\mathcal{R} k^2 \\ 1 & n^2 \pi^2 - k^2 - \text{P}_r \sigma\end{vmatrix}.\end{aligned} \hspace{\stretch{1}}(3.57)

For \sigma = 0, this gives us the critical value for the Rayleigh number

\begin{aligned}\mathcal{R} = \frac{(k^2 - n^2 \pi^2)^3}{k^2},\end{aligned} \hspace{\stretch{1}}(3.58)

the value that separates our stable and unstable solutions. On the other hand for

  1. \text{Real}(\sigma) > 0 (\Delta T > \Delta T_e), we have an instable system.
  2. \text{Real}(\sigma) < 0 (\Delta T > \Delta T_e), we have a stable system.

FIXME: The text has a positive sign on the n^2 term above. He actually solves the quadratic for \sigma, but I don’t see how that would make a difference. Is there an error here, or a typo in the text?

This is illustrated in figure (2).

Critical Rayleigh's number.


The instability means that we’ll have instable flows as illustrated in figure (3)

Instability due to heating.


Solving for these critical points we find


\begin{aligned}k_e^2 = \frac{\pi^2}{2}\end{aligned} \hspace{\stretch{1}}(3.59a)

\begin{aligned}\mathcal{R}_e = \frac{27 \pi^4}{4}\end{aligned} \hspace{\stretch{1}}(3.59b)


Multimedia presentations.

  1. Kelvin-Helmholtz instability.Colored salt water underneath, with unsalty water on top. Apparatus tilted causing flow of one over the other. Instability of the interface.

    See [2] for a really cool animation of a simulation of this effect. It ends up looking very fractal. Also interesting is the picture of this observed for real in the atmosphere of Saturn.

  2. A simulated mushroom cloud occurring with one fluid seeping into another. This looks it matches what we find under Rayleigh-Taylor instability in [3].
  3. plume, motion up through a denser fluid.
  4. Plateau-Rayleigh instability. Drop pinching off. See instability in the fluid channel feeding the drop. A crude illustration of this can be found in figure (4).

    Crude illustration of instability leading to a drop pinching off.


    A better illustrations (and animations) can be found in [4].

  5. Jet of water injected into a rotating tub on a turntable. Jet forms and surfaces.


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

[2] Wikipedia. Kelvin-helmholtz instability — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-April-2012].\%E2\%80\%93Helmholtz_instability&oldid=484301421.

[3] Wikipedia. Rayleigh-taylor instability — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-April-2012].\%E2\%80\%93Taylor_instability&oldid=483569989.

[4] Wikipedia. Plateau-rayleigh instability — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-April-2012].\%E2\%80\%93Rayleigh_instability&oldid=478499841.

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »