May 2024
M	T	W	T	F	S	S
	1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30	31

Posts Tagged ‘PHY454’

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.)]

Motivation.

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

Statement.

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{subequations}

$\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)$

\end{subequations}

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{subequations}

$\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)$

\end{subequations}

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}$

and

$\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}$

and

$\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{subequations}

$\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)$

\end{subequations}

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{subequations}

$\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)$

\end{subequations}

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{subequations}

$\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)$

\end{subequations}

so we have

\begin{subequations}

$\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)$

\end{subequations}

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)$

References

[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]. http://en.wikipedia.org/w/index.php?title=Taylor\%E2\%80\%93Couette_flow&oldid=483281707.

Posted in Math and Physics Learning. | Tagged: Couette flow, cylindrical coordinates, friction, Navier-Stokes, PHY454, torque | Leave a Comment »

Putting the stress tensor (and traction vector) into explicit vector form.

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

Motivation.

Exersize 6.1 from [1] is to show that the traction vector can be written in vector form (a rather curious thing to have to say) as

$\begin{aligned}\mathbf{t} = -p \hat{\mathbf{n}} + \mu ( 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla})\mathbf{u} + \hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})).\end{aligned} \hspace{\stretch{1}}(1.1)$

Note that the text uses a wedge symbol for the cross product, and I’ve switched to standard notation. I’ve done so because the use of a Geometric-Algebra wedge product also can be used to express this relationship, in which case we would write

$\begin{aligned}\mathbf{t} = -p \hat{\mathbf{n}} + \mu ( 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u} + (\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{n}}).\end{aligned} \hspace{\stretch{1}}(1.2)$

In either case we have

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{n}}=\hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})=\boldsymbol{\nabla}' (\hat{\mathbf{n}} \cdot \mathbf{u}') - (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u}\end{aligned} \hspace{\stretch{1}}(1.3)$

(where the primes indicate the scope of the gradient, showing here that we are operating only on $\mathbf{u}$ , and not $\hat{\mathbf{n}}$ ).

After computing this, lets also compute the stress tensor in cylindrical and spherical coordinates (a portion of that is also problem 6.10), something that this allows us to do fairly easily without having to deal with the second order terms that we encountered doing this by computing the difference of squared displacements.

We’ll work primarily with just the strain tensor portion of the traction vector expressions above, calculating

$\begin{aligned}2 {\mathbf{e}}_{\hat{\mathbf{n}}}=2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla})\mathbf{u} + \hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})=2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla})\mathbf{u} + (\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{n}}.\end{aligned} \hspace{\stretch{1}}(1.4)$

We’ll see that this gives us a nice way to interpret these tensor relationships. The interpretation was less clear when we computed this from the second order difference method, but here we see that we are just looking at the components of the force in each of the respective directions, dependent on which way our normal is specified.

Verifying the relationship.

Let’s start with the the plain old cross product version

$\begin{aligned}(\hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u}) + 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u})_i&=n_a (\boldsymbol{\nabla} \times \mathbf{u})_b \epsilon_{a b i} + 2 n_a \partial_a u_i \\ &=n_a \partial_r u_s \epsilon_{r s b} \epsilon_{a b i} + 2 n_a \partial_a u_i \\ &=n_a \partial_r u_s \delta_{ia}^{[rs]} + 2 n_a \partial_a u_i \\ &=n_a ( \partial_i u_a -\partial_a u_i ) + 2 n_a \partial_a u_i \\ &=n_a \partial_i u_a + n_a \partial_a u_i \\ &=n_a (\partial_i u_a + \partial_a u_i) \\ &=\sigma_{i a } n_a\end{aligned}$

We can also put the double cross product in wedge product form

$\begin{aligned}\hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})&=-I \hat{\mathbf{n}} \wedge (\boldsymbol{\nabla} \times \mathbf{u}) \\ &=-\frac{I}{2}\left(\hat{\mathbf{n}} (\boldsymbol{\nabla} \times \mathbf{u})- (\boldsymbol{\nabla} \times \mathbf{u}) \hat{\mathbf{n}}\right) \\ &=-\frac{I}{2}\left(-I \hat{\mathbf{n}} (\boldsymbol{\nabla} \wedge \mathbf{u})+ I (\boldsymbol{\nabla} \wedge \mathbf{u}) \hat{\mathbf{n}}\right) \\ &=-\frac{I^2}{2}\left(- \hat{\mathbf{n}} (\boldsymbol{\nabla} \wedge \mathbf{u})+ (\boldsymbol{\nabla} \wedge \mathbf{u}) \hat{\mathbf{n}}\right) \\ &=(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{n}}\end{aligned}$

Equivalently (and easier) we can just expand the dot product of the wedge and the vector using the relationship

$\begin{aligned}\mathbf{a} \cdot (\mathbf{c} \wedge \mathbf{d} \wedge \mathbf{e} \wedge \cdots )=(\mathbf{a} \cdot \mathbf{c}) (\mathbf{d} \wedge \mathbf{e} \wedge \cdots ) - (\mathbf{a} \cdot \mathbf{d}) (\mathbf{c} \wedge \mathbf{e} \wedge \cdots ) +\end{aligned} \hspace{\stretch{1}}(2.5)$

so we find

$\begin{aligned}((\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{n}} + 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u})_i&=(\boldsymbol{\nabla}' (\mathbf{u}' \cdot \hat{\mathbf{n}})-(\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u}+ 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u})_i \\ &=\partial_i u_a n_a+n_a \partial_a u_i \\ &=\sigma_{ia} n_a.\end{aligned}$

Cylindrical strain tensor.

Let’s now compute the strain tensor (and implicitly the traction vector) in cylindrical coordinates.

Our gradient in cylindrical coordinates is the familiar

$\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \hat{\boldsymbol{\phi}} \frac{1}{{r }}\frac{\partial {}}{\partial {\phi}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}},\end{aligned} \hspace{\stretch{1}}(3.6)$

and our cylindrical velocity is

$\begin{aligned}\mathbf{u} = \hat{\mathbf{r}} u_r + \hat{\boldsymbol{\phi}} u_\phi + \hat{\mathbf{z}} u_z.\end{aligned} \hspace{\stretch{1}}(3.7)$

Our curl is then

$\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u}&=\left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \hat{\boldsymbol{\phi}} \frac{1}{{r }}\frac{\partial {}}{\partial {\phi}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right)\wedge\left(\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\phi}} u_\phi + \hat{\mathbf{z}} u_z\right) \\ &=\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\phi}}\left(\partial_r u_\phi -\frac{1}{{r}} \partial_\phi u_r\right)+\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{z}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)+\hat{\mathbf{z}} \wedge \hat{\mathbf{r}}\left(\partial_z u_r - \partial_r u_z\right)+\frac{1}{{r}} \hat{\boldsymbol{\phi}} \wedge \left((\partial_\phi \hat{\mathbf{r}}) u_r+(\partial_\phi \hat{\boldsymbol{\phi}}) u_\phi\right)\end{aligned}$

Since $\partial_\phi \hat{\mathbf{r}} = \hat{\boldsymbol{\theta}}$ and $\partial_\phi \hat{\boldsymbol{\phi}} = -\hat{\mathbf{r}}$ , we have only one cross term and our curl is

$\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u}=\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\phi}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)+\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{z}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)+\hat{\mathbf{z}} \wedge \hat{\mathbf{r}}\left(\partial_z u_r - \partial_r u_z\right).\end{aligned} \hspace{\stretch{1}}(3.8)$

We can now move on to compute the directional derivatives and complete the strain calculation in cylindrical coordinates. Let’s consider this computation of the stress for normals in each direction in term.

With $\hat{\mathbf{n}} = \hat{\mathbf{r}}$ .

Our directional derivative component for a $\hat{\mathbf{r}}$ normal direction doesn’t have any cross terms

$\begin{aligned}2 (\hat{\mathbf{r}} \cdot \boldsymbol{\nabla}) \mathbf{u}&=2 \partial_r\left(\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\phi}} u_\phi + \hat{\mathbf{z}} u_z\right) \\ &=2\left(\hat{\mathbf{r}} \partial_r u_r + \hat{\boldsymbol{\phi}} \partial_r u_\phi + \hat{\mathbf{z}} \partial_r u_z\right).\end{aligned}$

Projecting our curl bivector onto the $\hat{\mathbf{r}}$ direction we have

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

Putting things together we have

$\begin{aligned}2 \mathbf{e}_{\hat{\mathbf{r}}}&=2\left(\hat{\mathbf{r}} \partial_r u_r + \hat{\boldsymbol{\phi}} \partial_r u_\phi + \hat{\mathbf{z}} \partial_r u_z\right)-\hat{\boldsymbol{\phi}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)+\hat{\mathbf{z}}\left(\partial_z u_r - \partial_r u_z\right) \\ &=\hat{\mathbf{r}}\left(2 \partial_r u_r\right)+\hat{\boldsymbol{\phi}}\left(2 \partial_r u_\phi-\partial_r u_\phi+\frac{1}{{r}} \partial_\phi u_r- \frac{u_\phi}{r}\right)+\hat{\mathbf{z}}\left(2 \partial_r u_z+\partial_z u_r - \partial_r u_z\right).\end{aligned}$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\mathbf{r}}} = - p \hat{\mathbf{r}} + 2 \mu e_{\hat{\mathbf{r}}},\end{aligned} \hspace{\stretch{1}}(3.9)$

we can now read off our components by taking dot products to yield

\begin{subequations}

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

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

$\begin{aligned}\sigma_{r z}=\mu \left( \frac{\partial {u_z}}{\partial {r}}+\frac{\partial {u_r}}{\partial {z}}\right).\end{aligned} \hspace{\stretch{1}}(3.10c)$

\end{subequations}

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\phi}}$ .

Our directional derivative component for a $\hat{\boldsymbol{\phi}}$ normal direction will have some cross terms since both $\hat{\mathbf{r}}$ and $\hat{\boldsymbol{\phi}}$ are functions of $\phi$

$\begin{aligned}2 (\hat{\boldsymbol{\phi}} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\frac{2}{r}\partial_\phi\left(\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\phi}} u_\phi + \hat{\mathbf{z}} u_z\right) \\ &=\frac{2}{r}\left(\hat{\mathbf{r}} \partial_\phi u_r + \hat{\boldsymbol{\phi}} \partial_\phi u_\phi + \hat{\mathbf{z}} \partial_\phi u_z+(\partial_\phi \hat{\mathbf{r}}) u_r + (\partial_\phi \hat{\boldsymbol{\phi}}) u_\phi\right) \\ &=\frac{2}{r}\left(\hat{\mathbf{r}} (\partial_\phi u_r - u_\phi) + \hat{\boldsymbol{\phi}} (\partial_\phi u_\phi + u_r )+ \hat{\mathbf{z}} \partial_\phi u_z\right) \\ \end{aligned}$

Projecting our curl bivector onto the $\hat{\boldsymbol{\phi}}$ direction we have

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\phi}}&=(\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\phi}}) \cdot \hat{\boldsymbol{\phi}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)+(\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{z}}) \cdot \hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)+(\hat{\mathbf{z}} \wedge \hat{\mathbf{r}}) \cdot \hat{\boldsymbol{\phi}}\left(\partial_z u_r - \partial_r u_z\right) \\ &=\hat{\mathbf{r}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)-\hat{\mathbf{z}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)\end{aligned}$

Putting things together we have

$\begin{aligned}2 \mathbf{e}_{\hat{\boldsymbol{\phi}}}&=\frac{2}{r}\left(\hat{\mathbf{r}} (\partial_\phi u_r - u_\phi) + \hat{\boldsymbol{\phi}} (\partial_\phi u_\phi + u_r )+ \hat{\mathbf{z}} \partial_\phi u_z\right)+\hat{\mathbf{r}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)-\hat{\mathbf{z}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right) \\ &=\hat{\mathbf{r}}\left(\frac{1}{r}\partial_\phi u_r-\frac{u_\phi}{r}+\partial_r u_\phi\right)+\frac{2}{r} \hat{\boldsymbol{\phi}}\left(\partial_\phi u_\phi + u_r\right)+\hat{\mathbf{z}}\left(\frac{1}{r} \partial_\phi u_z + \partial_z u_\phi\right).\end{aligned}$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\boldsymbol{\phi}}} = - p \hat{\boldsymbol{\phi}} + 2 \mu e_{\hat{\boldsymbol{\phi}}},\end{aligned} \hspace{\stretch{1}}(3.11)$

we can now read off our components by taking dot products to yield

\begin{subequations}

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

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

$\begin{aligned}\sigma_{\phi r}=\mu \left(\frac{1}{r}\frac{\partial {u_r}}{\partial {\phi}}-\frac{u_\phi}{r}+\frac{\partial {u_\phi}}{\partial {r}}\right).\end{aligned} \hspace{\stretch{1}}(3.12c)$

\end{subequations}

With $\hat{\mathbf{n}} = \hat{\mathbf{z}}$ .

Like the $\hat{\mathbf{r}}$ normal direction, our directional derivative component for a $\hat{\mathbf{z}}$ normal direction will not have any cross terms

$\begin{aligned}2 (\hat{\mathbf{z}} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\partial_z\left(\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\phi}} u_\phi + \hat{\mathbf{z}} u_z\right) \\ &=\hat{\mathbf{r}} \partial_z u_r + \hat{\boldsymbol{\phi}} \partial_z u_\phi + \hat{\mathbf{z}} \partial_z u_z\end{aligned}$

Projecting our curl bivector onto the $\hat{\mathbf{z}}$ direction we have

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\phi}}&=(\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\phi}}) \cdot \hat{\mathbf{z}}\left(\partial_r u_\phi-\frac{1}{{r}} \partial_\phi u_r+ \frac{u_\phi}{r}\right)+(\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{z}}) \cdot \hat{\mathbf{z}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)+(\hat{\mathbf{z}} \wedge \hat{\mathbf{r}}) \cdot \hat{\mathbf{z}}\left(\partial_z u_r - \partial_r u_z\right) \\ &=\hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)-\hat{\mathbf{r}}\left(\partial_z u_r - \partial_r u_z\right)\end{aligned}$

Putting things together we have

$\begin{aligned}2 \mathbf{e}_{\hat{\mathbf{z}}}&=2 \hat{\mathbf{r}} \partial_z u_r + 2 \hat{\boldsymbol{\phi}} \partial_z u_\phi + 2 \hat{\mathbf{z}} \partial_z u_z+\hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)-\hat{\mathbf{r}}\left(\partial_z u_r - \partial_r u_z\right) \\ &=\hat{\mathbf{r}}\left(2 \partial_z u_r -\partial_z u_r + \partial_r u_z\right)+\hat{\boldsymbol{\phi}}\left(2 \partial_z u_\phi +\frac{1}{{r}} \partial_\phi u_z- \partial_z u_\phi\right)+\hat{\mathbf{z}}\left(2 \partial_z u_z\right) \\ &=\hat{\mathbf{r}}\left(\partial_z u_r + \partial_r u_z\right)+\hat{\boldsymbol{\phi}}\left(\partial_z u_\phi +\frac{1}{{r}} \partial_\phi u_z\right)+\hat{\mathbf{z}}\left(2 \partial_z u_z\right).\end{aligned}$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\mathbf{z}}} = - p \hat{\mathbf{z}} + 2 \mu e_{\hat{\mathbf{z}}},\end{aligned} \hspace{\stretch{1}}(3.13)$

we can now read off our components by taking dot products to yield

\begin{subequations}

$\begin{aligned}\sigma_{z z}=-p + 2 \mu \frac{\partial {u_z}}{\partial {z}}\end{aligned} \hspace{\stretch{1}}(3.14a)$

$\begin{aligned}\sigma_{z r}=\mu \left(\frac{\partial {u_r}}{\partial {z}}+ \frac{\partial {u_z}}{\partial {r}}\right)\end{aligned} \hspace{\stretch{1}}(3.14b)$

$\begin{aligned}\sigma_{z \phi}=\mu \left(\frac{\partial {u_\phi}}{\partial {z}}+\frac{1}{{r}} \frac{\partial {u_z}}{\partial {\phi}}\right).\end{aligned} \hspace{\stretch{1}}(3.14c)$

\end{subequations}

Summary.

\begin{subequations}

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

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

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

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

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

\end{subequations}

Spherical strain tensor.

Having done a first order cylindrical derivation of the strain tensor, let’s also do the spherical case for completeness. Would this have much utility in fluids? Perhaps for flow over a spherical barrier?

We need the gradient in spherical coordinates. Recall that our spherical coordinate velocity was

$\begin{aligned}\frac{d\mathbf{r}}{dt} = \hat{\mathbf{r}} \dot{r} + \hat{\boldsymbol{\theta}} (r \dot{\theta}) + \hat{\boldsymbol{\phi}} ( r \sin\theta \dot{\phi} ),\end{aligned} \hspace{\stretch{1}}(4.16)$

and our gradient mirrors this structure

$\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \hat{\boldsymbol{\theta}} \frac{1}{{r }}\frac{\partial {}}{\partial {\theta}} + \hat{\boldsymbol{\phi}} \frac{1}{{r \sin\theta}} \frac{\partial {}}{\partial {\phi}}.\end{aligned} \hspace{\stretch{1}}(4.17)$

We also previously calculated \inbookref{phy454:continuumL2}{eqn:continuumL2:1010} the unit vector differentials

\begin{subequations}

$\begin{aligned}d\hat{\mathbf{r}} = \hat{\boldsymbol{\phi}} \sin\theta d\phi + \hat{\boldsymbol{\theta}} d\theta\end{aligned} \hspace{\stretch{1}}(4.18a)$

$\begin{aligned}d\hat{\boldsymbol{\theta}} = \hat{\boldsymbol{\phi}} \cos\theta d\phi - \hat{\mathbf{r}} d\theta\end{aligned} \hspace{\stretch{1}}(4.18b)$

$\begin{aligned}d\hat{\boldsymbol{\phi}} = -(\hat{\mathbf{r}} \sin\theta + \hat{\boldsymbol{\theta}} \cos\theta) d\phi,\end{aligned} \hspace{\stretch{1}}(4.18c)$

\end{subequations}

and can use those to read off the partials of all the unit vectors

$\begin{aligned}\frac{\partial \hat{\mathbf{r}}}{\partial \{r,\theta, \phi\}} &= \{0, \hat{\boldsymbol{\theta}}, \hat{\boldsymbol{\phi}} \sin\theta \} \\ \frac{\partial \hat{\boldsymbol{\theta}}}{\partial \{r,\theta, \phi\}} &= \{0, -\hat{\mathbf{r}}, \hat{\boldsymbol{\phi}} \cos\theta \} \\ \frac{\partial \hat{\boldsymbol{\phi}}}{\partial \{r,\theta, \phi\}} &= \{0, 0, -\hat{\mathbf{r}} \sin\theta -\hat{\boldsymbol{\theta}} \cos\theta \}.\end{aligned} \hspace{\stretch{1}}(4.19)$

Finally, our velocity in spherical coordinates is just

$\begin{aligned}\mathbf{u} = \hat{\mathbf{r}} u_r + \hat{\boldsymbol{\theta}} u_\theta + \hat{\boldsymbol{\phi}} u_\phi,\end{aligned} \hspace{\stretch{1}}(4.22)$

from which we can now compute the curl, and the directional derivative. Starting with the curl we have

$\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u}&=\left( \hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \hat{\boldsymbol{\theta}} \frac{1}{{r }}\frac{\partial {}}{\partial {\theta}} + \hat{\boldsymbol{\phi}} \frac{1}{{r \sin\theta}} \frac{\partial {}}{\partial {\phi}} \right) \wedge\left( \hat{\mathbf{r}} u_r + \hat{\boldsymbol{\theta}} u_\theta + \hat{\boldsymbol{\phi}} u_\phi \right) \\ &=\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\theta}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r\right)\\ & +\hat{\boldsymbol{\theta}} \wedge \hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta\right)\\ & +\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi\right)\\ & +\frac{1}{{r}} \hat{\boldsymbol{\theta}} \wedge \left(u_\theta \underbrace{\partial_\theta \hat{\boldsymbol{\theta}}}_{-\hat{\mathbf{r}}}+u_\phi \underbrace{\partial_\theta \hat{\boldsymbol{\phi}}}_{0}\right)\\ & +\frac{1}{{r \sin\theta}} \hat{\boldsymbol{\phi}} \wedge \left(u_\theta \underbrace{\partial_\phi \hat{\boldsymbol{\theta}}}_{\hat{\boldsymbol{\phi}} \cos\theta}+u_\phi \underbrace{\partial_\phi \hat{\boldsymbol{\phi}}}_{-\hat{\mathbf{r}} \sin\theta - \hat{\boldsymbol{\theta}} \cos\theta}\right).\end{aligned}$

So we have

$\begin{aligned}\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u}&=\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\theta}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)\\ & +\hat{\boldsymbol{\theta}} \wedge \hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)\\ & +\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.23)$

With $\hat{\mathbf{n}} = \hat{\mathbf{r}}$ .

The directional derivative portion of our strain is

$\begin{aligned}2 (\hat{\mathbf{r}} \cdot \boldsymbol{\nabla}) \mathbf{u}&=2 \partial_r (\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\theta}} u_\theta + \hat{\boldsymbol{\phi}} u_\phi ) \\ &=2 (\hat{\mathbf{r}} \partial_r u_r + \hat{\boldsymbol{\theta}} \partial_r u_\theta + \hat{\boldsymbol{\phi}} \partial_r u_\phi ).\end{aligned}$

The other portion of our strain tensor is

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{r}}&=(\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\theta}}) \cdot \hat{\mathbf{r}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)\\ & +(\hat{\boldsymbol{\theta}} \wedge \hat{\boldsymbol{\phi}}) \cdot \hat{\mathbf{r}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)\\ & +(\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{r}}) \cdot \hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right) \\ &=-\hat{\boldsymbol{\theta}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)\\ & +\hat{\boldsymbol{\phi}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right).\end{aligned}$

Putting these together we find

$\begin{aligned}2 {\mathbf{e}}_{\hat{\mathbf{r}}}&=2 (\hat{\mathbf{r}} \cdot \boldsymbol{\nabla})\mathbf{u} + (\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\mathbf{r}} \\ &=2 (\hat{\mathbf{r}} \partial_r u_r + \hat{\boldsymbol{\theta}} \partial_r u_\theta + \hat{\boldsymbol{\phi}} \partial_r u_\phi )-\hat{\boldsymbol{\theta}}\left(\partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)+\hat{\boldsymbol{\phi}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right) \\ &=\hat{\mathbf{r}}\left(2 \partial_r u_r\right)+\hat{\boldsymbol{\theta}}\left(2 \partial_r u_\theta-\partial_r u_\theta + \frac{1}{{r}} \partial_\theta u_r - \frac{u_\theta}{r}\right)+\hat{\boldsymbol{\phi}}\left(2 \partial_r u_\phi+ \frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right).\end{aligned}$

Which gives

$\begin{aligned}2 {\mathbf{e}}_{\hat{\mathbf{r}}}=\hat{\mathbf{r}}\left(2 \partial_r u_r\right)+\hat{\boldsymbol{\theta}}\left(\partial_r u_\theta+ \frac{1}{{r}} \partial_\theta u_r - \frac{u_\theta}{r}\right)+\hat{\boldsymbol{\phi}}\left(\partial_r u_\phi+ \frac{1}{{r \sin\theta}} \partial_\phi u_r- \frac{u_\phi}{r}\right)\end{aligned} \hspace{\stretch{1}}(4.24)$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\mathbf{r}}} = - p \hat{\mathbf{r}} + 2 \mu e_{\hat{\mathbf{r}}},\end{aligned} \hspace{\stretch{1}}(4.25)$

we can now read off our components by taking dot products

\begin{subequations}

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

$\begin{aligned}\sigma_{r \theta}=\mu \left(\frac{\partial {u_\theta}}{\partial {r}}+ \frac{1}{{r}} \frac{\partial {u_r}}{\partial {\theta}} - \frac{u_\theta}{r}\right)\end{aligned} \hspace{\stretch{1}}(4.26b)$

$\begin{aligned}\sigma_{r \phi}=\mu \left(\frac{\partial {u_\phi}}{\partial {r}}+ \frac{1}{{r \sin\theta}} \frac{\partial {u_r}}{\partial {\phi}}- \frac{u_\phi}{r}\right).\end{aligned} \hspace{\stretch{1}}(4.26c)$

\end{subequations}

This is consistent with (15.20) from [3] (after adjusting for minor notational differences).

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\theta}}$ .

Now let’s do the $\hat{\boldsymbol{\theta}}$ direction. The directional derivative portion of our strain will be a bit more work to compute because we have $\theta$ variation of the unit vectors

$\begin{aligned}(\hat{\boldsymbol{\theta}} \cdot \boldsymbol{\nabla}) \mathbf{u} &= \frac{1}{r} \partial_\theta (\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\theta}} u_\theta + \hat{\boldsymbol{\phi}} u_\phi ) \\ &= \frac{1}{r} \left( \hat{\mathbf{r}} \partial_\theta u_r + \hat{\boldsymbol{\theta}} \partial_\theta u_\theta + \hat{\boldsymbol{\phi}} \partial_\theta u_\phi \right)+\frac{1}{r} \left( (\partial_\theta \hat{\mathbf{r}}) u_r + (\partial_\theta \hat{\boldsymbol{\theta}}) u_\theta + (\partial_\theta \hat{\boldsymbol{\phi}}) u_\phi \right) \\ &= \frac{1}{r}\left(\hat{\mathbf{r}} \partial_\theta u_r + \hat{\boldsymbol{\theta}} \partial_\theta u_\theta + \hat{\boldsymbol{\phi}} \partial_\theta u_\phi \right)+\frac{1}{r} \left( \hat{\boldsymbol{\theta}} u_r - \hat{\mathbf{r}} u_\theta \right).\end{aligned}$

So we have

$\begin{aligned}2 (\hat{\boldsymbol{\theta}} \cdot \boldsymbol{\nabla}) \mathbf{u}=\frac{2}{r} \hat{\mathbf{r}} (\partial_\theta u_r- u_\theta)+ \frac{2}{r} \hat{\boldsymbol{\theta}} (\partial_\theta u_\theta+ u_r) + \frac{2}{r} \hat{\boldsymbol{\phi}} \partial_\theta u_\phi,\end{aligned} \hspace{\stretch{1}}(4.27)$

and can move on to projecting our curl bivector onto the $\hat{\boldsymbol{\theta}}$ direction. That portion of our strain tensor is

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\theta}}&=(\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\theta}}) \cdot \hat{\boldsymbol{\theta}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)\\ & +(\hat{\boldsymbol{\theta}} \wedge \hat{\boldsymbol{\phi}}) \cdot \hat{\boldsymbol{\theta}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)\\ & +(\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{r}}) \cdot \hat{\boldsymbol{\theta}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right) \\ &=\hat{\mathbf{r}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)-\hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right).\end{aligned}$

Putting these together we find

$\begin{aligned}2 {\mathbf{e}}_{\hat{\boldsymbol{\theta}}}&=2 (\hat{\boldsymbol{\theta}} \cdot \boldsymbol{\nabla})\mathbf{u} + (\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\theta}} \\ &= \frac{2}{r} \hat{\mathbf{r}} (\partial_\theta u_r - u_\theta )+ \frac{2}{r} \hat{\boldsymbol{\theta}} (\partial_\theta u_\theta + u_r )+ \frac{2}{r} \hat{\boldsymbol{\phi}} \partial_\theta u_\phi \\ &+\hat{\mathbf{r}}\left(\partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)-\hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta + \frac{u_\phi \cot\theta}{r}\right).\end{aligned}$

Which gives

$\begin{aligned}2 {\mathbf{e}}_{\hat{\boldsymbol{\theta}}}=\hat{\mathbf{r}} \left( \frac{1}{r} \partial_\theta u_r + \partial_r u_\theta- \frac{u_\theta}{r}\right)+\hat{\boldsymbol{\theta}} \left( \frac{2}{r} \partial_\theta u_\theta+ \frac{2}{r} u_r\right)+\hat{\boldsymbol{\phi}} \left(\frac{1}{r} \partial_\theta u_\phi+ \frac{1}{{r \sin\theta}} \partial_\phi u_\theta- \frac{u_\phi \cot\theta}{r}\right).\end{aligned} \hspace{\stretch{1}}(4.28)$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\boldsymbol{\theta}}} = - p \hat{\boldsymbol{\theta}} + 2 \mu e_{\hat{\boldsymbol{\theta}}},\end{aligned} \hspace{\stretch{1}}(4.29)$

we can now read off our components by taking dot products

\begin{subequations}

$\begin{aligned}\sigma_{\theta \theta}=-p+\mu \left( \frac{2}{r} \frac{\partial {u_\theta}}{\partial {\theta}}+ \frac{2}{r} u_r\right)\end{aligned} \hspace{\stretch{1}}(4.30a)$

$\begin{aligned}\sigma_{\theta \phi}=\mu \left(\frac{1}{r} \frac{\partial {u_\phi}}{\partial {\theta}}+ \frac{1}{{r \sin\theta}} \frac{\partial {u_\theta}}{\partial {\phi}}- \frac{u_\phi \cot\theta}{r}\right)\end{aligned} \hspace{\stretch{1}}(4.30b)$

$\begin{aligned}\sigma_{\theta r}= \mu \left(\frac{1}{r} \frac{\partial {u_r}}{\partial {\theta}} + \frac{\partial {u_\theta}}{\partial {r}}- \frac{u_\theta}{r}\right).\end{aligned} \hspace{\stretch{1}}(4.30c)$

\end{subequations}

This again is consistent with (15.20) from [3].

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\phi}}$ .

Finally, let’s do the $\hat{\boldsymbol{\phi}}$ direction. This directional derivative portion of our strain will also be a bit more work to compute because we have $\hat{\boldsymbol{\phi}}$ variation of the unit vectors

$\begin{aligned}(\hat{\boldsymbol{\phi}} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\frac{1}{r \sin\theta} \partial_\phi (\hat{\mathbf{r}} u_r + \hat{\boldsymbol{\theta}} u_\theta + \hat{\boldsymbol{\phi}} u_\phi) \\ &=\frac{1}{r \sin\theta}(\hat{\mathbf{r}} \partial_\phi u_r+\hat{\boldsymbol{\theta}} \partial_\phi u_\theta+\hat{\boldsymbol{\phi}} \partial_\phi u_\phi+(\partial_\phi \hat{\mathbf{r}} )u_r+(\partial_\phi \hat{\boldsymbol{\theta}} )u_\theta+(\partial_\phi \hat{\boldsymbol{\phi}} )u_\phi) \\ &=\frac{1}{r \sin\theta}(\hat{\mathbf{r}} \partial_\phi u_r+\hat{\boldsymbol{\theta}} \partial_\phi u_\theta+\hat{\boldsymbol{\phi}} \partial_\phi u_\phi+\hat{\boldsymbol{\phi}} \sin\thetau_r+\hat{\boldsymbol{\phi}} \cos\thetau_\theta-(\hat{\mathbf{r}} \sin\theta+ \hat{\boldsymbol{\theta}} \cos\theta)u_\phi)\end{aligned}$

So we have

$\begin{aligned}2 (\hat{\boldsymbol{\phi}} \cdot \boldsymbol{\nabla}) \mathbf{u}=2 \hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \frac{u_\phi}{r}\right)+2 \hat{\boldsymbol{\theta}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_\theta-\frac{1}{{r}} \cot\theta u_\phi\right)+2 \hat{\boldsymbol{\phi}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_\phi+ \frac{1}{{r}} u_r+ \frac{1}{{r}} \cot\theta u_\theta\right),\end{aligned} \hspace{\stretch{1}}(4.31)$

and can move on to projecting our curl bivector onto the $\hat{\boldsymbol{\phi}}$ direction. That portion of our strain tensor is

$\begin{aligned}(\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\phi}}&=(\hat{\mathbf{r}} \wedge \hat{\boldsymbol{\theta}}) \cdot \hat{\boldsymbol{\phi}}\left( \partial_r u_\theta - \frac{1}{{r}} \partial_\theta u_r + \frac{u_\theta}{r}\right)\\ & +(\hat{\boldsymbol{\theta}} \wedge \hat{\boldsymbol{\phi}}) \cdot \hat{\boldsymbol{\phi}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)\\ & +(\hat{\boldsymbol{\phi}} \wedge \hat{\mathbf{r}}) \cdot \hat{\boldsymbol{\phi}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right) \\ &=\hat{\boldsymbol{\theta}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)\\ &-\hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right).\end{aligned}$

Putting these together we find

$\begin{aligned}2 {\mathbf{e}}_{\hat{\boldsymbol{\theta}}}&=2 (\hat{\boldsymbol{\phi}} \cdot \boldsymbol{\nabla})\mathbf{u} + (\boldsymbol{\nabla} \wedge \mathbf{u}) \cdot \hat{\boldsymbol{\phi}} \\ &=2 \hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \frac{u_\phi}{r}\right)+2 \hat{\boldsymbol{\theta}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_\theta-\frac{1}{{r}} \cot\theta u_\phi\right)+2 \hat{\boldsymbol{\phi}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_\phi+ \frac{1}{{r}} u_r+ \frac{1}{{r}} \cot\theta u_\theta\right) \\ &+\hat{\boldsymbol{\theta}}\left(\frac{1}{{r}} \partial_\theta u_\phi - \frac{1}{{r \sin\theta}} \partial_\phi u_\theta+ \frac{u_\phi \cot\theta}{r}\right)-\hat{\mathbf{r}}\left(\frac{1}{{r \sin\theta}} \partial_\phi u_r - \partial_r u_\phi- \frac{u_\phi}{r}\right).\end{aligned}$

Which gives

$\begin{aligned}2 {\mathbf{e}}_{\hat{\boldsymbol{\phi}}}=\hat{\mathbf{r}} \left( \frac{ \partial_\phi u_r }{r \sin\theta}- \frac{u_\phi}{r}+ \partial_r u_\phi\right)+\hat{\boldsymbol{\theta}} \left(\frac{\partial_\phi u_\theta}{r \sin\theta}- \frac{u_\phi \cot\theta}{r}+\frac{\partial_\theta u_\phi}{r}\right)+2 \hat{\boldsymbol{\phi}} \left(\frac{\partial_\phi u_\phi}{r \sin\theta}+ \frac{u_r}{r}+ \frac{\cot\theta u_\theta}{r}\right).\end{aligned} \hspace{\stretch{1}}(4.32)$

For our stress tensor

$\begin{aligned}\boldsymbol{\sigma}_{\hat{\boldsymbol{\phi}}} = - p \hat{\boldsymbol{\phi}} + 2 \mu e_{\hat{\boldsymbol{\phi}}},\end{aligned} \hspace{\stretch{1}}(4.33)$

we can now read off our components by taking dot products

\begin{subequations}

$\begin{aligned}\sigma_{\phi \phi}=-p+2 \mu \left(\frac{1}{{r \sin\theta}} \frac{\partial {u_\phi}}{\partial {\phi}}+ \frac{u_r}{r}+ \frac{\cot\theta u_\theta}{r}\right)\end{aligned} \hspace{\stretch{1}}(4.34a)$

$\begin{aligned}\sigma_{\phi r}=\mu \left( \frac{1}{r \sin\theta} \frac{\partial {u_r}}{\partial {\phi}}- \frac{u_\phi}{r}+ \frac{\partial {u_\phi}}{\partial {r}}\right)\end{aligned} \hspace{\stretch{1}}(4.34b)$

$\begin{aligned}\sigma_{\phi \theta}= \mu \left(\frac{1}{r \sin\theta} \frac{\partial {u_\theta}}{\partial {\phi}}- \frac{u_\phi \cot\theta}{r}+\frac{1}{{r}} \frac{\partial {u_\phi}}{\partial {\theta}}\right).\end{aligned} \hspace{\stretch{1}}(4.34c)$

\end{subequations}

This again is consistent with (15.20) from [3].

Summary

\begin{subequations}

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

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

$\begin{aligned}\sigma_{\theta \phi}= \mu \left(\frac{1}{r \sin\theta} \frac{\partial {u_\theta}}{\partial {\phi}}- \frac{u_\phi \cot\theta}{r}+\frac{1}{{r}} \frac{\partial {u_\phi}}{\partial {\theta}}\right).\end{aligned} \hspace{\stretch{1}}(4.35e)$

\end{subequations}

References

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

[2] Peeter Joot. Continuum mechanics., chapter {Introduction and strain tensor.} http://sites.google.com/site/peeterjoot2/math2012/phy454.pdf.

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

Posted in Math and Physics Learning. | Tagged: consitutative relationship, curl, geometric algebra, PHY454, strain tensor, stress tensor, traction vector, wedge product | Leave a Comment »

Compilation of class notes for phy454h1s, continuum mechanics (including all lectures)

Posted by peeterjoot on April 6, 2012

As mentioned before the midterm, I’ve got a compilation of class notes here:

https://sites.google.com/site/peeterjoot2/math2012/phy454.pdf

This has been updated now with the full set of lecture notes from class. Included in this document are:

my lecture notes, often auxmented by whatever I felt required to make them understandable to myself
some solved problems
some stuff I was playing with on the side
links to the various mathematica notebooks (some of which have some cool stuff, and can be viewed after download with the free wolfram CDF player).

All of the following, previously posted individually are included, and will no longer be maintained independently

Apr 4, 2012 Thermal stability.

Mar 31, 2012 Channel flow with step pressure gradient

Mar 29, 2012 Stability.

Mar 28, 2012 Asymptotic solutions of ill conditioned equations.

Mar 23, 2012 Boundary layers and stream functions. The Blasius problem.

Mar 23, 2012 Problem Set 3. Velocity scaling, non-dimensionalisation, boundary layers.

Mar 21, 2012 Scaling to setup to solve the boundary layer problem.

Mar 16, 2012 Impulsive flow. Boundary layers. Oscillatory driven flow.

Mar 14, 2012 Midterm reflection.

Mar 14, 2012 Boundary Layers. Time dependent flow.

Mar 9, 2012 More on surface tension and Reynold’s number.

Mar 7, 2012 Non-dimensionality and scaling.

Mar 6, 2012 Inclined flow without constant height assumption. Setting up the problem, but not actually solving it.

Mar 2, 2012 Problem Set 2. Rectilinear flow with pressure gradient. Two layer flow driven by moving surface.

Mar 2, 2012 Hydrostatics. Surface normals and tangent vectors.

Mar 1, 2012 Steady state inclined flow down a plane of two layers of incompressible viscous fluids.

Mar 1, 2012 Steady state inclined flow down a plane of two layers of incompressible viscous equal density fluids.

Feb 17, 2012 Flow in a pipe. Gravity driven flow of a film.

Feb 15, 2012 Worked examples of unidirectional Navier-Stokes solutions.

Feb 10, 2012 Navier-Stokes equation.

Feb 8, 2012 Newtonian fluids. Mass conservation. Constitutive relation. Incompressible fluids.

Feb 8, 2012 Problem Set 1. Stress, Strain, Traction vector. Force free equilibrium.

Feb 3, 2012 Phasor description of elastic waves. Fluid dynamics.

Feb 1, 2012 P-waves and S-waves.

Jan 27, 2012 Compatibility condition and elastostatics.

Jan 25, 2012 Constitutive relationship.

Jan 20, 2012 Strain tensor components.

Jan 18, 2012 Strain tensor review. Stress tensor.

Jan 13, 2012 Introduction and strain tensor.

Jan 11, 2012 Overview.

Posted in Math and Physics Learning. | Tagged: continuum mechanics, fluid mechanics, Mathematica, PHY454 | 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.)]

Disclaimer.

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

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)$
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)$
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)$
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{subequations}

$\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)$

\end{subequations}

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)$

where

$\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)$

where

$\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)$

Writing

$\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)$

where

$\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{subequations}

$\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)$

\end{subequations}

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

$\text{Real}(\sigma) > 0$ ( $\Delta T > \Delta T_e$ ), we have an instable system.
$\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{subequations}

$\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)$

\end{subequations}

Multimedia presentations.

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.
A simulated mushroom cloud occurring with one fluid seeping into another. This looks it matches what we find under Rayleigh-Taylor instability in [3].
plume, motion up through a denser fluid.
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].
Jet of water injected into a rotating tub on a turntable. Jet forms and surfaces.

References

[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]. http://en.wikipedia.org/w/index.php?title=Kelvin\%E2\%80\%93Helmholtz_instability&oldid=484301421.

[3] Wikipedia. Rayleigh-taylor instability — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-April-2012]. http://en.wikipedia.org/w/index.php?title=Rayleigh\%E2\%80\%93Taylor_instability&oldid=483569989.

[4] Wikipedia. Plateau-rayleigh instability — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-April-2012]. http://en.wikipedia.org/w/index.php?title=Plateau\%E2\%80\%93Rayleigh_instability&oldid=478499841.

Posted in Math and Physics Learning. | Tagged: Navier-Stokes, PHY454, Prandtl number, Rayleigh number, thermal stability | Leave a Comment »

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

Posted by peeterjoot on April 4, 2012

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

Disclaimer.

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

Review. Laminar boundary layer theory.

In the boundary layer we found

Continuity equation
$\begin{aligned}\frac{\partial u}{\partial x} + \frac{\partial {v}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.1)$
$x$ momentum equation
$\begin{aligned}u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = - \frac{1}{{\rho}} \frac{\partial {p}}{\partial x} + \nu \frac{\partial^2 u}{\partial y^2}\end{aligned} \hspace{\stretch{1}}(2.2)$
$y$ momentum equation
$\begin{aligned}\frac{\partial {p}}{\partial y} = 0\end{aligned} \hspace{\stretch{1}}(2.3)$

In the inviscid region $p$ is a constant in $y$

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

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

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

Starting with

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

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

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

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

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

so that our boundary layer equations are

\begin{subequations}

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

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

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

\end{subequations}

With boundary conditions

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

Fluid flow over a flat plate (Blasius problem).

Define a similarity variable $\eta$

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

Suppose we want

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

Since we have

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

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

We can make the transformation

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

We can introduce stream functions

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

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

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

Now introduce a similarity variable

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

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

Also note that this last step requires:

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

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

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

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

but that doesn’t obviously give us 3.20?

Ah. We argued that

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

and that the larger of the viscous terms was

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

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

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

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

Our boundary conditions are

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

Deriving the equation of motion.

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

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

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

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

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

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

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

so that

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

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

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

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

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

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

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

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

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

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

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

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

and

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

and

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

Our $x$ component of Navier-Stokes now takes the form

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

or (assuming $g \ne 0$ )

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

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

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

It’s argued that we expect

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

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

\begin{subequations}

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

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

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

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

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

\end{subequations}

and boundary value conditions

\begin{subequations}

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

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

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

\end{subequations}

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

Numeric solution.

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

Figure 2: Boundary layer solution to flow over plate.

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

Singular perturbation theory.

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

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

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

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

finally

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

If we have

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

and

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

then

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

when $\text{Re} >> 1$ .

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

Our EOM is

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

with

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

as $x, y \rightarrow L$

performing a non-dimensionalization we have

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

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

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

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

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

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

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

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

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

This is illustrated in figure (3).

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

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

When $x \gg \epsilon$ we have approximately

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

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

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

References

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

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

Posted in Math and Physics Learning. | Tagged: Blasius problem, Laminar boundary layers, Navier-Stokes, PHY454, Singular perturbation theory | Leave a Comment »

Channel flow with step pressure gradient

Posted by peeterjoot on April 1, 2012

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

Motivation and statement.

This is problem 2.5 from [1].

Viscous fluid is at rest in a two-dimensional channel between stationary rigid walls with $y = \pm h$ . For $t \ge 0$ a constant pressure gradient $P = -dp/dx$ is imposed. Show that $u(y, t)$ satifies

$\begin{aligned}\frac{\partial {u}}{\partial {t}} = \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{P}{\rho},\end{aligned} \hspace{\stretch{1}}(1.1)$

and give suitable initial and boundary conditions. Find $u(y, t)$ in form of a Fourier series, and show that the flow approximates to steady channel flow when $t \gg h^2/\nu$ .

Solution

With only horizontal components to the flow, the Navier-Stokes equations for incompressible flow are

\begin{subequations}

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

$\begin{aligned}\frac{\partial {u}}{\partial {x}} = 0\end{aligned} \hspace{\stretch{1}}(2.2b)$

\end{subequations}

Substitution of 2.2b into 2.2a gives us

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

Our equation to solve is therefore

$\begin{aligned}\boxed{\frac{\partial {u}}{\partial {t}} = \Theta(t) \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.}\end{aligned} \hspace{\stretch{1}}(2.4)$

This equation, for $t < 0$ , allows for solutions

$\begin{aligned}u = A y + B\end{aligned} \hspace{\stretch{1}}(2.5)$

but the problem states that the fluid is at rest initially, so we don’t really have to solve anything (i.e. $A = B = 0$ ).

The no-slip conditions introduce boundary value conditions $u(\pm h, t) = 0$ .

For $t \ge 0$ we have

$\begin{aligned}\frac{\partial {u}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.6)$

If we attempt separation of variables with $u(y, t) = Y(y) T(t)$ , our equation takes the form

$\begin{aligned}T' Y = \frac{P}{\rho} + \nu T Y''.\end{aligned} \hspace{\stretch{1}}(2.7)$

We see that the non-homogenous term prevents successful application of separation of variables. Let’s modify our problem by attempting to recast our equation into a homogenous form by adding a particular solution for the steady state flow problem. That problem was the solution of

$\begin{aligned}\frac{\partial^2 {{}}}{\partial {{y}}^2} u_s(y, 0) = -\frac{P}{\rho \nu}\end{aligned} \hspace{\stretch{1}}(2.8)$

which has solution

$\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right) + A y + B.\end{aligned} \hspace{\stretch{1}}(2.9)$

The freedom to incoporate an $h^2$ constant into the equation as an integration constant has been employed, knowing that it will kill the $y^2$ contributions at $y = \pm h$ to make the boundary condition matching easier. Our no-slip conditions give us

$\begin{aligned}0 &= A h + B \\ 0 &= -A h + B.\end{aligned} \hspace{\stretch{1}}(2.10)$

Adding this we have $2 B = 0$ , and subtracting gives us $2 A h = 0$ , so a specific solution that matches our required boundary value (and initial value) conditions is just the steady state channel flow solution we are familiar with

$\begin{aligned}u_s(y, 0) = \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.12)$

Let’s now assume that our general solution has the form

$\begin{aligned}u(y, t) = u_H(y, t) + u_s(y, 0).\end{aligned} \hspace{\stretch{1}}(2.13)$

Applying the Navier-Stokes equation to this gives us

$\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \frac{P}{\rho} + \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2} + \nu \frac{\partial^2 {{u_s}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.14)$

But from 2.8, we see that all we have left is a homogenous problem in $u_H$

$\begin{aligned}\frac{\partial {u_H}}{\partial {t}} = \nu \frac{\partial^2 {{u_H}}}{\partial {{y}}^2},\end{aligned} \hspace{\stretch{1}}(2.15)$

where our boundary value conditions are now given by

$\begin{aligned}0 &= u_H(\pm h, t) + u_s(\pm h) \\ &= u_H(\pm h, t),\end{aligned}$

and

$\begin{aligned}0 &= u(y, 0) \\ &= u_H(y, 0) + \frac{P}{2 \rho \nu} \left( h^2 - y^2 \right),\end{aligned}$

or
\begin{subequations}

$\begin{aligned}u_H(\pm h, t) = 0\end{aligned} \hspace{\stretch{1}}(2.16a)$

$\begin{aligned}u_H(y, 0) = -\frac{P}{2 \rho \nu} \left( h^2 - y^2 \right).\end{aligned} \hspace{\stretch{1}}(2.16b)$

\end{subequations}

Now we can apply separation of variables with $u_H = T(t) Y(y)$ , yielding

$\begin{aligned}T' Y = \nu T Y'',\end{aligned} \hspace{\stretch{1}}(2.17)$

$\begin{aligned}\frac{T'}{T} = \nu \frac{Y''}{Y} = \text{constant} = - \nu \alpha^2.\end{aligned} \hspace{\stretch{1}}(2.18)$

Here a positive constant $\nu \alpha^2$ has been used assuming that we want a solution that is damped with time.

Our solutions are

$\begin{aligned}T &\propto e^{- \nu \alpha^2 t} \\ Y &= A \sin \alpha y + B \cos\alpha y,\end{aligned} \hspace{\stretch{1}}(2.19)$

$\begin{aligned}u_H(y, t) = \sum_\alpha e^{-\alpha^2 \nu t} \left( A_\alpha \sin \alpha y + B_\alpha \cos\alpha y \right).\end{aligned} \hspace{\stretch{1}}(2.21)$

We have constraints on $\alpha$ due to our boundary value conditions. For our sin terms to be solutions we require

$\begin{aligned}\sin (\alpha (\pm h)) = \sin n \pi\end{aligned} \hspace{\stretch{1}}(2.22)$

and for our cosine terms to be solutions we require

$\begin{aligned}\cos (\alpha (\pm h)) = \cos \left( \frac{\pi}{2} + n \pi \right)\end{aligned} \hspace{\stretch{1}}(2.23)$

$\begin{aligned}\alpha &= \frac{2 n \pi}{2 h} \\ \alpha &= \frac{2 n + 1 \pi}{2 h}\end{aligned} \hspace{\stretch{1}}(2.24)$

respectively.

Our homogeneous solution therefore takes the form

$\begin{aligned}u_H(y, t) = C_0 + \sum_{m > 0} C_m e^{ -(m \pi/2h)^2 \nu t } \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{$ latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{ $m$ odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.26)$

Our undetermined constants should be provided by the boundary value constraint at $t = 0$ 2.16b, leaving us to solve the Fourier problem

$\begin{aligned}-\frac{P}{2 \mu} \left( h^2 - y^2 \right)=\sum_{m \ge 0} C_m \left\{\begin{array}{l l}\sin \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{$ latex m$ even} \\ \cos \left( \frac{ m \pi y }{2 h} \right) & \quad \mbox{ $m$ odd} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.27)$

Multiplying by a sine and integrating will clearly give zero (even times odd function over a symmetric interval). Let’s see if there’s any scaling required to select out the $C_m$ term

$\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{2h}{\pi} \int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) \pi dy/2h \\ &=\frac{2h}{\pi} \int_{-\pi/2}^{\pi/2}\cos m x \cos n x dx \\ &=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( \cos( (m - n) \pi/2 ) +\cos( (m + n) \pi/2 ) \right) dx\end{aligned}$

Note that since $m$ and $n$ must be odd, $m \pm n = 2 c$ for some integer $c$ , so this integral is zero unless $m = n$ (consider $m = 2 a + 1, n = 2 b + 1$ ). For the $m = n$ term we have

$\begin{aligned}\int_{-h}^h \cos \left( \frac{ m \pi y }{2 h} \right) \cos \left( \frac{ n \pi y }{2 h} \right) dy&=\frac{h}{\pi} \int_{-\pi/2}^{\pi/2}\left( 1+\cos( m \pi ) \right) dx \\ &=h\end{aligned}$

Therefore, our constants $C_m$ (for odd $m$ ) are given by

$\begin{aligned}C_m &= -\frac{P h }{2 \mu} \int_{-h}^h\left( 1 - \left( \frac{y}{h}\right)^2 \right)\cos \left( \frac{ m \pi y }{2 h} \right) dy \\ &=-\frac{P h^2 }{2 \mu} \int_{-1}^1\left( 1 - x^2 \right)\cos \left( \frac{ m \pi x }{2} \right) x \\ \end{aligned}$

With $m = 2 n + 1$ , we have

$\begin{aligned}C_{2 n + 1} = -\frac{16 P h^2 (-1)^n}{\mu \pi^3 (2 n + 1)^3}.\end{aligned} \hspace{\stretch{1}}(2.28)$

For that calculation see: phy454/channelFlowWithStepPressureGradient.cdf.

Our complete solution is

$\begin{aligned}u(y, t) =\frac{P h^2}{2 \mu} \left( 1 - \left( \frac{y}{h} \right)^2 \right)- \frac{16 P h^2 }{\mu \pi^3 }\sum_{n = 0}^\infty \frac{(-1)^n}{(2 n + 1)^3}e^{ -((2 n + 1) \pi/2h)^2 \nu t } \cos \left( \frac{ (2 n + 1) \pi y }{2 h} \right) .\end{aligned} \hspace{\stretch{1}}(2.29)$

The largest of the damped exponentials above is the $n = 0$ term which is

$\begin{aligned}e^{ - \pi^2 \nu t /h^2 },\end{aligned} \hspace{\stretch{1}}(2.30)$

so if $\nu t >> h^2$ these terms all die off, leaving us with just the steady state.

Rather remarkably, this Fourier series is actually a very good fit even after only a single term. Using the viscosity and density of water, $h = 1 \text{cm}$ , and $P = 3 \times \mu_{\text{water}} \times (2 \text{cm}/{s})/ h^2$ (parameterizing the pressure gradient by the average velocity it will induce), a plot of the parabola that we are fitting to and the difference of that from the first Fourier term is shown in figure (1).

Figure 1: Parabolic channel flow steady state, and difference from first Fourier term.

The higher order corrections are even smaller. Even the first order deviations from the parabola that we are fitting to is a correction on the scale of $1/100$ of the height of the parabola. This is illustrated in figure (2) where the magnitude of the first 5 deviations from the steady state are plotted.

Figure 2: Difference from the steady state for the first five Fourier terms.

An animation of the time evolution above can be found in figure (3). If this animation is unavailable, it can also be found at http://youtu.be/0vZuv9HBtmo.

Figure 3: Time evolution of channel flow velocity profile after turning on a constant pressure gradient.

It’s also interesting to look at the very earliest part of the time evolution. Observe in figure (4) (or http://youtu.be/dDkx8iLwOew) the oscillatory phenomina. Could some of that be due to not running with enough Fourier terms in this early part of the evolution when more terms are probably significant?

Figure 4: Early time evolution of channel flow velocity profile after turning on a constant pressure gradient.

References

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

Posted in Math and Physics Learning. | Tagged: animation, channel flow, Mathematica, Navier-Stokes, PHY454, time evolution | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 20: Asymptotic solutions of ill conditioned equations. Taught by Prof. K. Das.

Posted by peeterjoot on March 30, 2012

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

Disclaimer.

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

Two ill conditioned LDEs.

We’ll consider two cases, both ones that we can solve exactly

With $u(0) = 1$ , and letting $\epsilon \rightarrow 0$ , we’ll look at solutions of the ill conditioned LDE
$\begin{aligned}\epsilon \frac{du}{dy} + u = y\end{aligned} \hspace{\stretch{1}}(2.1)$
With $u(0) = 0$ , $u(1) = 2$ , and $0 < \epsilon \ll 1$ we’ll look at the second order ill conditioned LDE
$\begin{aligned}\epsilon \frac{d^2u}{dy^2} + \frac{du}{dy} = 1\end{aligned} \hspace{\stretch{1}}(2.2)$

The first order LDE.

Exact solution.

Our homogeneous equation is

$\begin{aligned}\epsilon \frac{du}{dy} + u = 0,\end{aligned} \hspace{\stretch{1}}(2.3)$

with solution

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

Looking for a solution of the form

$\begin{aligned}u = A(y) e^{-y/\epsilon},\end{aligned} \hspace{\stretch{1}}(2.5)$

we find

$\begin{aligned}\epsilon A' e^{-y/\epsilon} = y,\end{aligned} \hspace{\stretch{1}}(2.6)$

and integrate to find

$\begin{aligned}A(y) = (y - \epsilon) e^{x/\epsilon} + C.\end{aligned} \hspace{\stretch{1}}(2.7)$

Application of the boundary value constraints give us

$\begin{aligned}u = ( 1 + \epsilon ) e^{-y/\epsilon} + y - \epsilon.\end{aligned} \hspace{\stretch{1}}(2.8)$

This is plotted in figure (1).

Figure 1: Plot of exact solution to simple first order ill conditioned LDE.

Limiting cases.

We want to consider the limiting case where

$\begin{aligned}0 < \epsilon \ll 1,\end{aligned} \hspace{\stretch{1}}(2.9)$

and we let $\epsilon \rightarrow 0$ . If $y = O(1)$ , then we have

$\begin{aligned}u \approx 1 \times 0 + y,\end{aligned} \hspace{\stretch{1}}(2.10)$

or just

$\begin{aligned}u = y.\end{aligned} \hspace{\stretch{1}}(2.11)$

However, if $y = O(\epsilon)$ then we have to be more careful constructing an approximation. When $y$ is very small, but $\epsilon$ is also of the same order of smallness we have

$\begin{aligned}e^{-y/\epsilon} \ne 0.\end{aligned} \hspace{\stretch{1}}(2.12)$

If $\epsilon \rightarrow 0$ and $y \rightarrow O(\epsilon)$

$\begin{aligned}e^{-y/\epsilon} \rightarrow e^{-\epsilon O(1) /\epsilon} \rightarrow e^{-O(1)}\end{aligned} \hspace{\stretch{1}}(2.13)$

$\begin{aligned}u \approx e^{-y/\epsilon}\end{aligned} \hspace{\stretch{1}}(2.14)$

Approximate solution in the inner region.

When $y = O(\epsilon)$ define a new scale

$\begin{aligned}Y = \frac{y}{\epsilon}\end{aligned} \hspace{\stretch{1}}(2.15)$

$\begin{aligned}y = O(1)\end{aligned} \hspace{\stretch{1}}(2.16)$

so that our LDE takes the form

$\begin{aligned}\frac{du}{dY} + u = \epsilon Y\end{aligned} \hspace{\stretch{1}}(2.17)$

When $\epsilon \rightarrow 0$ we have

$\begin{aligned}\frac{du}{dY} + u \approx 0.\end{aligned} \hspace{\stretch{1}}(2.18)$

We have solution

$\begin{aligned}\ln u = -Y + \ln C,\end{aligned} \hspace{\stretch{1}}(2.19)$

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

Question: Couldn’t we just Laplace transform.
Answer given: We’d still get into trouble when we take $\epsilon \rightarrow 0$ . My comment: I don’t think that’s strictly true. In an example like this where we have an exact solution, a Laplace transform technique should also yield that solution. I think the real trouble will come when we attempt to incorporate the non-linear inertial terms of the Navier-Stokes equation.

Second order example.

Exact solution.

We saw above in the first order system that our specific solution was polynomial. While that was found by the method of variation of parameters, it seems obvious in retrospect. Let’s start by looking for such a solution, starting with a first order polynomial

$\begin{aligned}u = A y + B.\end{aligned} \hspace{\stretch{1}}(2.21)$

Application of our LDE operator on this produces

$\begin{aligned}A &= 1 \\ B &= 0.\end{aligned}$

Now let’s move on to find a solution to the homogeneous equation

$\begin{aligned}\epsilon \frac{d^2u}{dy^2} + \frac{du}{dy} = 0\end{aligned} \hspace{\stretch{1}}(2.22)$

As usual, we look for the characteristic equation by assuming a solution of the form $u = e^{m y}$ . This gives us

$\begin{aligned}\epsilon m^2 + m = (\epsilon m + 1) m = 0,\end{aligned} \hspace{\stretch{1}}(2.23)$

with roots

$\begin{aligned}m = 0, -1/\epsilon.\end{aligned} \hspace{\stretch{1}}(2.24)$

So our homogeneous equation has the form

$\begin{aligned}u(y) = A e^{-y/\epsilon} + B,\end{aligned} \hspace{\stretch{1}}(2.25)$

and our full solution is

$\begin{aligned}u(y) = A e^{-y/\epsilon} + B + y\end{aligned} \hspace{\stretch{1}}(2.26)$

with the constants $A$ and $B$ to be determined from our boundary value conditions. We find

$\begin{aligned}0 &= u(0) = A + B + 0 \\ 2 &= u(1) = A e^{-1/\epsilon} + B + 1.\end{aligned}$

We’ve got $B = -A$ and by subtracting

$\begin{aligned}A ( e^{-1/\epsilon} -1 ) = 1.\end{aligned} \hspace{\stretch{1}}(2.27)$

So the exact solution is

$\begin{aligned}u = y + \frac{e^{-y/\epsilon} - 1}{e^{1/\epsilon} - 1}.\end{aligned} \hspace{\stretch{1}}(2.28)$

This is plotted in figure (2).

Figure 2: Plot of our ill conditioned second order LDE

Solution in the regular region.

For small $\epsilon$ relative to $y$ our LDE is approximately

$\begin{aligned}\frac{du}{dy} = 1,\end{aligned} \hspace{\stretch{1}}(2.29)$

which has solution

$\begin{aligned}u = y + C\end{aligned} \hspace{\stretch{1}}(2.30)$

Our $u(1) = 2$ boundary value constraint gives us

$\begin{aligned}C = 1.\end{aligned} \hspace{\stretch{1}}(2.31)$

Our solution in the regular region where $\epsilon \rightarrow 0$ and $y = O(1)$ is therefore just

$\begin{aligned}u = y + 1.\end{aligned} \hspace{\stretch{1}}(2.32)$

Solution in the ill conditioned region.

Now let’s consider the inner region. We’ll see below that when $y = O(\epsilon)$ , and we allow both $\epsilon$ and $y$ tend to zero independently, we have approximately

$\begin{aligned}u \sim 1 - e^{-y/\epsilon}.\end{aligned} \hspace{\stretch{1}}(2.33)$

We’ll now show this. We start with a helpful change of variables as we did in the first order case

$\begin{aligned}Y = \frac{y}{\epsilon}.\end{aligned} \hspace{\stretch{1}}(2.34)$

When $y = O(\epsilon)$ and $Y = O(1)$ we have

$\begin{aligned}\not{{\epsilon}} \frac{d^2 u}{\not{{\epsilon}} \epsilon dY^2} + \frac{1}{{\epsilon}} \frac{du}{dY} = 1\end{aligned} \hspace{\stretch{1}}(2.35)$

$\begin{aligned}\frac{d^2 u}{dY^2} + \frac{du}{dY} = \epsilon.\end{aligned} \hspace{\stretch{1}}(2.36)$

This puts the LDE into a non ill conditioned form, and allows us to let $\epsilon \rightarrow 0$ . We have approximately

$\begin{aligned}\frac{d^2 u}{dY^2} + \frac{du}{dY} = 0.\end{aligned} \hspace{\stretch{1}}(2.37)$

We’ve solved this in our exact solution work above (in a slightly more general form), and thus in this case we have just

$\begin{aligned}u = A + B e^{-Y}\end{aligned} \hspace{\stretch{1}}(2.38)$

at $Y = 0$ we have

$\begin{aligned}u = A + B = 0.\end{aligned} \hspace{\stretch{1}}(2.39)$

so that

$\begin{aligned}B = -A.\end{aligned} \hspace{\stretch{1}}(2.40)$

and we find for the inner region

$\begin{aligned}u = A (1 - e^{-Y}) = A( 1 - e^{-y/\epsilon} ) \sim 1 - e^{-y/\epsilon}.\end{aligned} \hspace{\stretch{1}}(2.41)$

Taking these independent solutions for the inner and outer regions and putting them together into a coherent form (called matched asymptotic expansion) is a rich and tricky field. For info on that we’ve been referred to [1].

References

[1] EJ Hinch. Perturbation methods, volume 6. Cambridge Univ Pr, 1991.

Posted in Math and Physics Learning. | Tagged: Ill conditioned LDEs, Navier-Stokes, PHY454, regular perturbation, singular perturbation. | Leave a Comment »

PHY454H1S Continuum mechanics. Problem Set 3. Velocity scaling, non-dimensionalisation, boundary layers.

Posted by peeterjoot on March 30, 2012

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

Disclaimer.

This problem set is as yet ungraded.

Problem Q1.

Background.

In fluid convection problems one can make several choices for characteristic velocity scales. Some choices are given below for example:

$U_1 = g \alpha d^2 \nabla T/\nu$
$U_2 = \nu/d$
$U_3 = \sqrt{ g \alpha d \nabla T }$
$U_4 = \kappa/ d$

where $g$ is the acceleration due to gravity, $\alpha = ({\partial {V}}/{\partial {T}})/V$ is the coefficient of volume expansion, $d$ length scale associated with the problem, $\nabla T$ is the applied temperature difference, $\nu$ is the kinematic viscosity and $\kappa$ is the thermal diffusivity.

For $U_1$ :
Observing that

$\begin{aligned}\left[ \frac{\partial {\mathbf{u}}}{\partial {t}} \right] = [ \nu \boldsymbol{\nabla}^2 \mathbf{u} ]\end{aligned} \hspace{\stretch{1}}(2.1)$

we must have

$\begin{aligned}[\nu] = \frac{1}{{[t] [\boldsymbol{\nabla}^2]}} = \frac{1}{{T}} L^2\end{aligned} \hspace{\stretch{1}}(2.2)$

We also find

$\begin{aligned}[\alpha] = \frac{1}{{[V]}} \left[ \frac{\partial {V}}{\partial {T}} \right] = \frac{1}{{[K]}},\end{aligned} \hspace{\stretch{1}}(2.3)$

so that

$\begin{aligned}[U_1] = \frac{L}{T \not{{T}}} \not{{\frac{1}{{K}}}} L^2 \not{{K}} \frac{\not{{T}}}{L^2} = \frac{L}{T}\end{aligned} \hspace{\stretch{1}}(2.4)$
For $U_2$ :
$\begin{aligned}[U_2] = \frac{L^2}{T} \frac{1}{{L}} = \frac{L}{T}.\end{aligned} \hspace{\stretch{1}}(2.5)$
For $U_3$
$\begin{aligned}[U_3] = \sqrt{ \frac{L}{T^2} \frac{1}{{K}} L K } = \frac{L}{T}\end{aligned} \hspace{\stretch{1}}(2.6)$
For $U_4$
According to http://scienceworld.wolfram.com/physics/ThermalDiffusivity.html, the thermal diffusivity is defined by

$\begin{aligned}\frac{\partial {T}}{\partial {t}} = \kappa \boldsymbol{\nabla}^2 T\end{aligned} \hspace{\stretch{1}}(2.7)$

so that

$\begin{aligned}[\kappa] = \frac{1}{{[t][\boldsymbol{\nabla}^2]}} = \frac{L^2}{T}\end{aligned} \hspace{\stretch{1}}(2.8)$

That gives us

$\begin{aligned}[U_4] = \frac{L^2}{T} \frac{1}{{L}} = \frac{L}{T}.\end{aligned} \hspace{\stretch{1}}(2.9)$

We’ve verified that all of these have dimensions of velocity.

Part 1. Statement. Check whether the dimensions match in each case above.

Solution.

Part 2. Statement. Pure liquid.

For pure liquid, say pure water at room temperature, one has the following estimates in cgs units:

$\begin{aligned}\alpha &\sim 10^{-4} \\ \kappa &\sim 10^{-3} \\ \nu &\sim 10^{-2}\end{aligned}$

For a $d \sim 1 \text{cm}$ layer depth and a ten degree temperature drop convective velocities have been experimentally measured of about $10^{-2}$ .

With $g \sim 10^{-3}$ , calculate the values of $U_1$ , $U_2$ , $U_3$ , and $U_4$ . Which ones of the characteristic velocities $(U_1$ , $U_2$ , $U_3, U_4)$ do you think are suitable for nondimensionalising the velocity in Navier-Stokes/Energy equation describing the water convection problem?

We have

$\begin{aligned}U_1 &\sim 10^{-3} 10^{-4} (1)^2 10^1 \frac{1}{{10^{-2}}} = 10^{-4} \\ U_2 &\sim 10^{-2}/1 = 10^{-2} \\ U_3 &\sim \sqrt{ 10^{-3} 10^{-4} (1) 10^1 } = 10^{-3} \\ U_4 &\sim 10^{-3}/1 = 10^{-3}\end{aligned} \hspace{\stretch{1}}(2.10)$

Use of $U_2 = \nu/d$ gives the closest match to the measured characteristic velocity of $10^{-2}$ .

Solution.

Part 3. Statement. Mantle convection.

For mantle convection, we have

$\begin{aligned}\alpha &\sim 10^{-5} \\ \nu &\sim 10^{21} \\ \kappa &\sim 10^{-2} \\ d &\sim 10^8 \\ \nabla T &\sim 10^3,\end{aligned}$

and the actual convective mantle velocity is $10^{-8}$ . Which of the characteristic velocities should we use to nondimensionalise Navier-Stokes/Energy equations describing mantle convection?

Solution

Let’s compute the characteristic velocities again with the mantle numbers

$\begin{aligned}U_1 &\sim \frac{10^{-3} 10^{-5} 10^{16} 10^3 }{10^{21}} = 10^{-10} \\ U_2 &\sim \frac{10^{21}}{10^8} = 10^{13} \\ U_3 &\sim \sqrt{ 10^{-3} 10^{-5} 10^8 10^3 } \sim 10^1 \\ U_4 &\sim \frac{10^{-2}}{10^8} = 10^{-10}\end{aligned} \hspace{\stretch{1}}(2.14)$

Both $U_1$ and $U_4$ come close to the actual convective mantle velocity of $10^{-8}$ . Use of $U_1$ to nondimensionalise is probably best, since it has more degrees of freedom, and includes the gravity term that is probably important for such large masses.

Problem Q2.

Statement

Nondimensionalise N-S equation

$\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 g \hat{\mathbf{z}}\end{aligned} \hspace{\stretch{1}}(3.18)$

where $\hat{\mathbf{z}}$ is the unit vector in the $z$ direction. You may scale:

velocity with the characteristic velocity $U$ ,
time with $R/U$ , where $R$ is the characteristic length scale,
pressure with $\rho U^2$ ,

Reynolds number $\text{Re} = R U \rho/ \mu$ and Froude number $\text{Fr} = g R/U$ .

Solution

Let’s start by dividing by $g \rho$ , to make all terms (most obviously the $\hat{\mathbf{z}}$ term) dimensionless.

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

Our suggested replacements are

$\begin{aligned}\mathbf{u} &= U \mathbf{u}' \\ \frac{\partial {}}{\partial {t}} &= \frac{U}{R} \frac{\partial {}}{\partial {t'}} \\ p &= \rho U^2 p' \\ \boldsymbol{\nabla} &= \frac{1}{{R}} \boldsymbol{\nabla}'.\end{aligned} \hspace{\stretch{1}}(3.20)$

Plugging these in we have

$\begin{aligned}\frac{U^2}{g R} \frac{\partial {\mathbf{u}'}}{\partial {t'}} + \frac{U^2}{g R} (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = - \frac{\not{{\rho}} U^2}{g \not{{\rho}} R} \boldsymbol{\nabla}' p' + \frac{\mu U}{g \rho R^2} {\boldsymbol{\nabla}'}^2 \mathbf{u}' + \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(3.24)$

Making a $\text{Fr} = gR/U$ replacement, using the Froude number, we have

$\begin{aligned}\frac{U}{\text{Fr}} \frac{\partial {\mathbf{u}'}}{\partial {t'}} + \frac{U}{\text{Fr}} (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = - \frac{U}{\text{Fr}} \boldsymbol{\nabla}' p' + \frac{\mu }{\text{Fr} \rho R} {\boldsymbol{\nabla}'}^2 \mathbf{u}' + \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(3.25)$

Scaling by $\text{Fr}/U$ we tidy things up a bit, and also allow for insertion of the Reynold’s number

$\begin{aligned}\frac{\partial {\mathbf{u}'}}{\partial {t'}} + (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = - \boldsymbol{\nabla}' p' + \frac{1}{\text{Re}} {\boldsymbol{\nabla}'}^2 \mathbf{u}' + \frac{\text{Fr}}{U}\hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(3.26)$

Observe that the dimensions of Froude’s number is that of velocity

$\begin{aligned}[\text{Fr}] = [g] T = \frac{L}{T},\end{aligned} \hspace{\stretch{1}}(3.27)$

so that the end result is dimensionless as desired. We also see that Froude’s number, characterizes the significance of the body force for fluid flow at the characteristic velocity. This is consistent with [1] where it was stated that the Froude number is used to determine the resistance of a partially submerged object moving through water, and permits the comparison of objects of different sizes (complete with pictures of canoes of various sizes that Froude built for such study).

Problem Q3.

Statement

In case of Stokes’ boundary layer problem (see class note) calculate shear stress on the plate $y = 0$ . What is the phase difference between the velocity of the plate $U(t) = U_0 \cos\omega t$ and the shear stress on the plate?

Solution

We found in class that the velocity of the fluid was given by

$\begin{aligned}u(y, t) = U_0 e^{-\lambda y} \cos(\lambda y - \omega t)\end{aligned} \hspace{\stretch{1}}(4.28)$

where

$\begin{aligned}\lambda = \sqrt{\frac{\omega}{2 \nu}}\end{aligned} \hspace{\stretch{1}}(4.29)$

Calculating our shear stress we find

$\begin{aligned}\mu \frac{\partial {u}}{\partial {y}} &= U_0 \lambda \mu e^{-\lambda y}\left(-1- \sin(\lambda y - \omega t)\right)\end{aligned}$

and on the plate ( $y = 0$ ) this is just

$\begin{aligned}{\left.{{\mu \frac{\partial {u}}{\partial {y}}}}\right\vert}_{{y = 0}} = U_0 \lambda \mu (-1 + \sin(\omega t)).\end{aligned} \hspace{\stretch{1}}(4.30)$

We’ve got a constant term, plus one that is sinusoidal. Observing that

$\begin{aligned}\cos x &= \text{Real}( e^{ix} ) \\ \sin x &= \text{Real}( -i e^{ix} ) = \text{Real}( e^{i (x - \pi/2)} ),\end{aligned} \hspace{\stretch{1}}(4.31)$

The phase difference between the non-constant portion of the shear stress at the plate, and the plate velocity $U(t) = U_0 \cos\omega t$ is just $-\pi/2$ . The shear stress at the plate lags the driving velocity by 90 degrees.

References

[1] Wikipedia. Froude number — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 27-March-2012]. http://en.wikipedia.org/w/index.php?title=Froude_number&oldid=479498080.

Posted in Math and Physics Learning. | Tagged: boundary layers, characteristic length, characteristic velocity, Froude's number, Navier-Stokes, nondimensionalisation, PHY454, Reynold's number, velocity scaling | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 21: Stability. Taught by Prof. K. Das.

Posted by peeterjoot on March 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.)]

Disclaimer.

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

Stability. Some graphical illustrations.

What do we mean by stability? A configuration is stable if after a small disturbance it returns to it’s original position. A couple systems to consider are shown in figures (\ref{fig:continuumL21:continuumL21Fig1a}), (\ref{fig:continuumL21:continuumL21Fig1b}) and (\ref{fig:continuumL21:continuumL21Fig1c}).

Figure 1: Stable well configuration.

Figure 2: Instable peak configuration.

Figure 3: Stable tabletop configuration.

We can examine how a displacement $\delta x$ changes with time after making it. In a stable configuration without friction we will induce an oscillation as plotted in figure (\ref{fig:continuumL21:continuumL21Fig2a}) for the parabolic configuration. With friction we’ll have a damping effect. This is plotted for the parabolic well in figure (\ref{fig:continuumL21:continuumL21Fig2b}).

Figure 4: Displacement time evolution in undamped well system.

Figure 5: Displacement time evolution in damped well system.

For the inverted parabola our displacement takes the form of figure (\ref{fig:continuumL21:continuumL21Fig3})

Figure 6: Time evolution of displacement in instable parabolic configuration.

For the ball on the table, assuming some friction that stops the ball, fairly quickly, we’ll have a displacement as illustrated in figure (\ref{fig:continuumL21:continuumL21Fig3b})

Figure 7: Time evolution of displacement in tabletop configuration.

Characterizing stability

Let’s suppose that our displacement can be described in exponential form

$\begin{aligned}\delta x \sim e^{\sigma t}\end{aligned} \hspace{\stretch{1}}(3.1)$

where $\sigma$ is the \textit{growth rate of perturbation}, and is in general a complex number of the form

$\begin{aligned}\sigma = \sigma_\text{R} + i \sigma_\text{I}\end{aligned} \hspace{\stretch{1}}(3.2)$

Case I. Oscillatory unstability

A system of the form

$\begin{aligned}\sigma_{\text{R}} &= 0 \\ \sigma_{\text{I}} &> 0\end{aligned}$

\textit{oscillatory unstable}. An example of this is the undamped parabolic system illustrated above.

Case II. Marginal unstability.

$\begin{aligned}\delta x &\sim e^{\sigma_{\text{R}} t} e^{i \sigma_{\text{I}} t} \\ &\sim e^{\sigma_{\text{R}} t} \left( \cos \sigma_{\text{I}} t + i \sin \sigma_{\text{I}} t \right)\end{aligned}$

We’ll call systems of the form

$\begin{aligned}\sigma_{\text{I}} &= 0 \\ \sigma_{\text{R}} &> 0\end{aligned}$

\textit{marginally unstable}. We can have unstable systems with $\sigma_{\text{I}} \ne 0$ but still $\sigma_{\text{R}} > 0$ , but these are less common.

Case III. Neutral stability.

$\begin{aligned}\sigma &= 0 \\ \sigma_{\text{R}} &= \sigma_{\text{I}} = 0\end{aligned}$

An example of this was the billiard table example where the ball moved to a new location on the table after being bumped slightly.

A mathematical description.

For a discussion of stability in fluids we’ll not only have to incorporate the Navier-Stokes equation as we’ve done, but will also have to bring in the heat equation. Unfortunately that isn’t in the scope of this course to derive. Let’s consider as system heated on a bottom plate, and consider the fluid and convection due to heating. This system is illustrated in figure (\ref{fig:continuumL21:continuumL21Fig4})

Figure 8: Fluid in cavity heated on the bottom plate.

We start with Navier-Stokes as normal

$\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 \hat{\mathbf{z}} g.\end{aligned} \hspace{\stretch{1}}(4.3)$

For steady state with $\mathbf{u} = 0$ initially (our base state), we’ll call the following the equation of the base state

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

We’ll allow perturbations of each of our variables

$\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}$

After perturbation NS takes the form

$\begin{aligned}(\rho_s + \delta \rho )\frac{\partial {(0 + \delta \mathbf{u})}}{\partial {t}} + (\rho_s + \delta \rho) ((0 + \delta \mathbf{u}) \cdot \boldsymbol{\nabla}) (0 + \delta \mathbf{u}) = - \boldsymbol{\nabla} (p_s + \delta p) + \mu \boldsymbol{\nabla}^2 (0 + \delta \mathbf{u}) - (\rho_s + \delta \rho) \hat{\mathbf{z}} g\end{aligned} \hspace{\stretch{1}}(4.5)$

Retaining only terms that are of first order of smallness.

$\begin{aligned}\rho_s \frac{\partial {\delta \mathbf{u}}}{\partial {t}} = - \boldsymbol{\nabla} p_s - \boldsymbol{\nabla} \delta p + \mu \boldsymbol{\nabla}^2 \delta \mathbf{u} - \rho_s \hat{\mathbf{z}} g - \delta \rho \hat{\mathbf{z}} g\end{aligned} \hspace{\stretch{1}}(4.6)$

applying our equation of base state 4.4, we have

$\begin{aligned}\rho_s \frac{\partial {\delta \mathbf{u}}}{\partial {t}} = \not{{\rho_s \hat{\mathbf{z}} g}} - \boldsymbol{\nabla} \delta p + \mu \boldsymbol{\nabla}^2 \delta \mathbf{u} - \not{{\rho_s \hat{\mathbf{z}} g}} - \delta \rho \hat{\mathbf{z}} g,\end{aligned} \hspace{\stretch{1}}(4.7)$

$\begin{aligned}\rho_s \frac{\partial {\delta \mathbf{u}}}{\partial {t}} = - \boldsymbol{\nabla} \delta p + \mu \boldsymbol{\nabla}^2 \delta \mathbf{u} - \delta \rho \hat{\mathbf{z}} g.\end{aligned} \hspace{\stretch{1}}(4.8)$

we can write

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

Applying the divergence operation on both sides, and using $\boldsymbol{\nabla} \cdot \mathbf{u} = 0$ so that $\boldsymbol{\nabla} \cdot \delta \mathbf{u} = 0$ we have

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

$\begin{aligned}\frac{1}{{\rho_s}} \boldsymbol{\nabla}^2 \delta p = - (\hat{\mathbf{z}} \cdot \boldsymbol{\nabla} ) \frac{\delta \rho}{\rho_s} g.\end{aligned} \hspace{\stretch{1}}(4.11)$

Assuming that $\rho_s$ is constant (actually that’s already been done above), we can cancel it, leaving

$\begin{aligned}\boldsymbol{\nabla}^2 \delta p = - (\hat{\mathbf{z}} \cdot \boldsymbol{\nabla} ) g \delta \rho = -g \frac{\partial {}}{\partial {z}} \delta \rho.\end{aligned} \hspace{\stretch{1}}(4.12)$

operating once more with ${\partial {}}/{\partial {z}}$ we have

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

Going back to 4.9 and taking only the $z$ component we have

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

$\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \nu \boldsymbol{\nabla}^2 \right) \boldsymbol{\nabla}^2 \delta w &= -\frac{1}{{\rho_s}} \frac{\partial { \boldsymbol{\nabla}^2 \delta p}}{\partial {z}} - \frac{g}{\rho_s} \boldsymbol{\nabla}^2 \delta \rho \\ &= -\frac{g}{\rho_s} \frac{\partial^2 {{\delta \rho}}}{\partial {{z}}^2} - \frac{g}{\rho_s} \boldsymbol{\nabla}^2 \delta \rho \\ &= -\frac{g}{\rho_s} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right)\delta \rho \\ &=g \alpha \left( \frac{\partial^2 {{}}}{\partial {{x}}^2}+\frac{\partial^2 {{}}}{\partial {{y}}^2}\right)\delta T\end{aligned}$

in the last step we use the following assumed relation for temperature

$\begin{aligned}\delta \rho = - \rho_s \alpha \delta T.\end{aligned} \hspace{\stretch{1}}(4.15)$

Here $\alpha$ is the coefficient of thermal expansion. This is just a statement that expansion and temperature are related (as we heat something, it expands), with the ratio of the density change relative to the original being linearly related to the change in temperature.

We have finally

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

Solving this is the Rayleigh-Benard instability problem.

While this is a fourth order differential equation, it’s still the same sort of problem logically as we’ve been working on. Our boundary value conditions at $z = 0$ are

$\begin{aligned}u, v, w, \delta u, \delta v, \delta w = 0.\end{aligned} \hspace{\stretch{1}}(4.17)$

Also relevant will be a similar equation relating temperature and fluid flow rate

$\begin{aligned}\left( \frac{\partial {}}{\partial {t}} - \kappa \boldsymbol{\nabla}^2 \right) \delta T = \Delta T \frac{\delta w}{d},\end{aligned} \hspace{\stretch{1}}(4.18)$

which we will cover in the next (and final) lecture of the course.

Posted in Math and Physics Learning. | Tagged: instability, Navier-Stokes, PHY454, Rayleigh-Benard problem. | Leave a Comment »

PHY454H1S Continuum mechanics. Problem Set 2. Two layer flow driven by moving surface. (revisited)

Posted by peeterjoot on March 26, 2012

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

Problem Q3 (revisited).

I’d produced the following sketches. For a higher viscosity bottom layer $\mu^{(1)} > \mu^{(2)}$ , this should look something like figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig4}) whereas for the higher viscosity on the top, these would be roughly flipped as in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig5}).

This superposition can be justified since we have no $(\mathbf{u} \cdot \boldsymbol{\nabla})\mathbf{u}$ term in the Navier-Stokes equations for these systems.

Exact solutions.

The figures above are kind of rough. It’s not actually hard to solve the system above. After some simplification, I find with Mathematica the following solution

$\begin{aligned}u^{(1)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(1)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(1)}} \\ u^{(2)}(y) &= -\frac{\mu^{(2)} U -G h^2}{\mu^{(1)}+\mu^{(2)}}-\frac{y \left(G h^2 (\mu^{(2)} -\mu^{(1)}) +2 \mu^{(1)} \mu^{(2)} U \right)}{2 h \mu^{(2)} \left(\mu^{(1)}+\mu^{(2)}\right)}-\frac{G y^2}{2 \mu^{(2)}}.\end{aligned} \hspace{\stretch{1}}(3.49)$

Should we wish a more exact plot for any specific values of the viscosities, we could plot exactly with software the vector field described by these velocities.

I suppose it is cheating to use Mathematica and then say that the solution is easy? To make amends for being lazy with my algebra, let’s show that it is easy to do manually too. I’ll do the same problem manually, but generalize it slightly. We can do this easily if we just be a bit marter with our integration constants. Let’s solve the problem for the upper and lower walls moving with velocity $V_2$ and $V_1$ respectively, and let the heights from the interface be $h_2$ and $h_1$ respectively.

We have the same set of differential equations to solve, but now let’s write our solution with the undetermined coefficients expressed as

$\begin{aligned}u^{(2)} &= -\frac{G}{2 \mu^{(2)}}\left( h_2^2 - y^2 \right) + \frac{A_2}{\mu^{(2)}} (y - h_2) + B_2 \\ u^{(1)} &= -\frac{G}{2 \mu^{(1)}}\left( h_1^2 - y^2 \right) + \frac{A_1}{\mu^{(1)}} (y + h_1) + B_1.\end{aligned} \hspace{\stretch{1}}(3.55)$

Now it’s super easy to match the boundary conditions at $y = -h_1$ and $y = h_2$ (the lower and upper walls respectively). Clearly the integration constants $B_1,B_2$ are just the velocities. Matching the tangential component of the traction vectors at $y = 0$ we have

$\begin{aligned}A_2 = A_1\end{aligned} \hspace{\stretch{1}}(3.55)$

and matching velocities at $y = 0$ gives us

$\begin{aligned}-\frac{G}{2 \mu^{(2)}}h_2^2 - \frac{A_2}{\mu^{(2)}} h_2 + V_2 = -\frac{G}{2 \mu^{(1)}}h_1^2 + \frac{A_1}{\mu^{(1)}} h_1 + V_1.\end{aligned} \hspace{\stretch{1}}(3.55)$

This gives us

$\begin{aligned}u^{(2)} &= \frac{G h_2^2}{2 \mu^{(2)}}\left(\frac{y^2}{h_2^2} - 1 \right) + A \frac{ h_2}{\mu^{(2)}} \left( \frac{y}{h_2} - 1 \right) + V_2 \\ u^{(1)} &= \frac{G h_1^2}{2 \mu^{(1)}}\left(\frac{y^2}{h_1^2} - 1 \right) + A \frac{ h_1}{\mu^{(1)}} \left( \frac{y}{h_1} + 1 \right) + V_1 \\ A &=\frac{V_2 - V_1+ \frac{G h_1^2}{2 \mu^{(1)}}-\frac{G h_2^2}{2 \mu^{(2)}}}{\frac{h_1}{\mu^{(1)}}+\frac{h_2}{\mu^{(2)}}}.\end{aligned} \hspace{\stretch{1}}(3.55)$

Plotting this with sliders or animation in Mathematica is a fun way to explore visualizing this. The results vary widely depending on the various parameters. Such a Mathematica notebook can be found in https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/problemSetIIQ3PlotWithManipulate.cdf. A pair of cached animations with variation of the pressure gradient for $v_1 = 0$ , $h_1 = h_2$ , showing the superposition of the shear and channel flow solutions can be found here

With $\mu^{(1)} > \mu^{(2)}$ fluids with densities of mercury and water in the lower (1) and upper (2) layers respectively.
With $\mu^{(2)} > \mu^{(1)}$ , this is plotted fluids with densities of water and mercury in the lower (1) and upper (2) layers respectively.

Posted in Math and Physics Learning. | Tagged: animation, channel flow, Manipulate, Mathematica, Navier-Stokes, PHY454, shear flow, superpositiion | Leave a Comment »

« Previous Entries

	Determining the alig… on C structure alignment pad…
	Manas shetty on Cartesian to spherical change…
	peeterjoot on Derivative recurrence relation…
	Daniel Pires on Derivative recurrence relation…
	peeterjoot on Curious problem using the vari…

Math, physics, perl, and programming obscurity.

Categories

Archives

Recent Posts

Meta

Recent Comments

People not reading this blog: 7,179,979,522 minus:

Subscribe

Posts Tagged ‘PHY454’

Motivation.

Statement.

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

Solution: Part 2. Fixing the constants.

Solution: Part 3. Friction and torque.

References

Motivation.

Verifying the relationship.

Cylindrical strain tensor.

With .

With .

With .

Summary.

Spherical strain tensor.

With .

With .

With .

Summary

References

Disclaimer.

Review. Rayleigh Benard Problem

Application of the perturbation to the energy equation.

Non-dimensionalisation of the thermal velocity equation.

Non-dimensionalization of the energy equation

Normal mode analysis.

Back to our coupled equations.

Multimedia presentations.

References

Disclaimer.

Review. Laminar boundary layer theory.

Fluid flow over a flat plate (Blasius problem).

Deriving the equation of motion.

Numeric solution.

Singular perturbation theory.

References

Motivation and statement.

Solution

References

Disclaimer.

Two ill conditioned LDEs.

The first order LDE.

Exact solution.

Limiting cases.

Approximate solution in the inner region.

Second order example.

Exact solution.

Solution in the regular region.

Solution in the ill conditioned region.

References

Disclaimer.

Problem Q1.

Background.

Part 1. Statement. Check whether the dimensions match in each case above.

Solution.

Part 2. Statement. Pure liquid.

Solution.

Part 3. Statement. Mantle convection.

Solution

Problem Q2.

Statement

Solution

Problem Q3.

Statement

Solution

References

Disclaimer.

Stability. Some graphical illustrations.

Characterizing stability

Case I. Oscillatory unstability

Case II. Marginal unstability.

Case III. Neutral stability.

With $\hat{\mathbf{n}} = \hat{\mathbf{r}}$ .

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\phi}}$ .

With $\hat{\mathbf{n}} = \hat{\mathbf{z}}$ .

With $\hat{\mathbf{n}} = \hat{\mathbf{r}}$ .

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\theta}}$ .

With $\hat{\mathbf{n}} = \hat{\boldsymbol{\phi}}$ .