March 2012
M	T	W	T	F	S	S
« Feb		Apr »
	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

Archive for March, 2012

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.

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 »

Classical Mechanics Euler Angles (lecture for phy354 taught by Prof E. Poppitz)

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

Pictorially.

We want to look at some of the trig behind expressing general rotations. We can perform a general rotation by a sequence of successive rotations. One such sequence is a rotation around the $z,x,z$ axes in sequence. Application of a rotation of angle $\phi$ takes us from our original (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig1}) frame to that of (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig2}). A second rotation around the (new) $x$ axis by angle $\theta$ takes us to (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig3}), and finally a rotation of $\psi$ around the (new) $z$ axis, takes us to (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig4}).

A composite image of all of these rotations taken together can be found in figure (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig5}).

Relating the two pairs of coordinate systems.

Let’s look at this algebraically instead, using figure (\ref{fig:classicalMechanicsEulerAngles:classicalMechanicsEulerAnglesFig6}) as a guide.

Step 1. Rotation of $\phi$ around $z$

$\begin{aligned}\begin{bmatrix}x' \\ y' \\ z'\end{bmatrix}=\begin{bmatrix}\cos\phi & \sin\phi & 0 \\ -\sin\phi & \cos\phi & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}x \\ y \\ z\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.1)$

Step 2. Rotation around $x'$ .

$\begin{aligned}\begin{bmatrix}x'' \\ y'' \\ z''\end{bmatrix}=\begin{bmatrix}1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta & 0 \\ 0 & -\sin\theta & \cos\theta & 0 \end{bmatrix}\begin{bmatrix}x' \\ y' \\ z'\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.2)$

Step 3. Rotation around $z''$ .

So, our full rotation is the composition of the rotation matrices

$\begin{aligned}\begin{bmatrix}x''' \\ y''' \\ z'''\end{bmatrix}=\begin{bmatrix}\cos\psi & \sin\psi & 0 \\ -\sin\psi & \cos\psi & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta & 0 \\ 0 & -\sin\theta & \cos\theta & 0 \end{bmatrix}\begin{bmatrix}\cos\phi & \sin\phi & 0 \\ -\sin\phi & \cos\phi & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}x \\ y \\ z\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.4)$

Let’s introduce some notation and write this as

$\begin{aligned}B_z(\alpha)=\begin{bmatrix}\cos\alpha & \sin\alpha & 0 \\ -\sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.5)$

$\begin{aligned}B_x(\theta) =\begin{bmatrix}1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta & 0 \\ 0 & -\sin\theta & \cos\theta & 0 \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.6)$

so that we have the mapping

$\begin{aligned}\mathbf{r} \rightarrow B_z(\psi) B_x(\theta) B_z(\phi) \mathbf{r}\end{aligned} \hspace{\stretch{1}}(2.7)$

Now let’s write

$\begin{aligned}\mathbf{r} = \begin{bmatrix}x_1 \\ x_2 \\ x_3\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.8)$

We’ll call

$\begin{aligned}A(\psi, \theta, \phi) = B_z(\psi) B_x(\theta) B_z(\phi) \end{aligned} \hspace{\stretch{1}}(2.9)$

so that

$\begin{aligned}x_i' = \sum_{j = 1}^3 A_{ij} x_j.\end{aligned} \hspace{\stretch{1}}(2.10)$

We’ll drop the explicit summation sign, so that the summation over repeated indexes are implied

$\begin{aligned}x_i' = A_{ij} x_j.\end{aligned} \hspace{\stretch{1}}(2.11)$

This matrix $A(\psi, \theta, \phi)$ is in fact a general parameterization of the $3 \times 3$ special orthogonal matrices. The set of three angles $\theta$ , $\phi$ , $\psi$ parameterizes all rotations in 3dd space. Transformations that preserve $\mathbf{a} \cdot \mathbf{b}$ and have unit determinant.

In symbols we must have

$\begin{aligned}A^\text{T} A = 1\end{aligned} \hspace{\stretch{1}}(2.12)$

$\begin{aligned}\det A = +1.\end{aligned} \hspace{\stretch{1}}(2.13)$

Having solved this auxiliary problem, we now want to compute the angular velocity.

We want to know how to express the coordinates of a point that is fixed in the body. i.e. We are fixing $x_i'$ and now looking for $x_i$ .

The coordinates of a point that has $x'$ , $y'$ and $z'$ in a body-fixed frame, in the fixed frame are $x$ , $y$ , $z$ . That is given by just inverting the matrix

$\begin{aligned}\begin{bmatrix}x_1 \\ x_2 \\ x_3\end{bmatrix}&=A^{-1}(\psi, \theta, \phi)\begin{bmatrix}x_1' \\ x_2' \\ x_3'\end{bmatrix} \\ &=B_z^{-1}(\phi)B_x^{-1}(\theta)B_z^{-1}(\psi)\begin{bmatrix}x_1' \\ x_2' \\ x_3'\end{bmatrix} \\ &=B_z^{\text{T}}(\phi)B_x^{\text{T}}(\theta)B_z^{\text{T}}(\psi)\begin{bmatrix}x_1' \\ x_2' \\ x_3'\end{bmatrix}\end{aligned}$

Here we’ve used the fact that $B_x$ and $B_z$ are orthogonal, so that their inverses are just their transposes.

We have finally

$\begin{aligned}\begin{bmatrix}x_1 \\ x_2 \\ x_3\end{bmatrix}=B_z(-\phi)B_x(-\theta)B_z(-\psi)\begin{bmatrix}x_1' \\ x_2' \\ x_3'\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(2.14)$

If we assume that $\psi$ , $\theta$ and $\phi$ are functions of time, and compute $d \mathbf{r}/dt$ . Starting with

$\begin{aligned}x_i = [A^{-1}(\phi, \theta, \psi)]_{ij} x_j',\end{aligned} \hspace{\stretch{1}}(2.15)$

$\begin{aligned}\Delta x_i = \left([A^{-1}(\phi + \Delta \phi, \theta + \Delta \theta, \psi + \Delta \psi)]_{ij} -[A^{-1}(\phi, \theta, \psi)]_{ij} \right)x_j'.\end{aligned} \hspace{\stretch{1}}(2.16)$

For small changes, we can Taylor expand and retain only the first order terms. Doing that and dividing by $\Delta t$ we have

$\begin{aligned}\frac{d x_i}{dt} =\left(\frac{\partial {}}{\partial {\psi}} A^{-1}_{ij} \dot{\psi}+\frac{\partial {}}{\partial {\theta}} A^{-1}_{ij} \dot{\theta}+\frac{\partial {}}{\partial {\phi}} A^{-1}_{ij} \dot{\phi}\right) x_j'\end{aligned} \hspace{\stretch{1}}(2.17)$

Now, we use

$\begin{aligned}x_j' = A_{jl} x_l,\end{aligned} \hspace{\stretch{1}}(2.18)$

so that we have

$\begin{aligned}\frac{d x_i}{dt} =\left(\left(\frac{\partial {}}{\partial {\psi}} A^{-1}_{ij} \right) A_{jl} \dot{\psi}+\left(\frac{\partial {}}{\partial {\theta}} A^{-1}_{ij} \right) A_{jl} \dot{\theta}+\left(\frac{\partial {}}{\partial {\phi}} A^{-1}_{ij} \right) A_{jl} \dot{\phi}\right) x_l.\end{aligned} \hspace{\stretch{1}}(2.19)$

We are looking for a relation of the form

$\begin{aligned}\frac{d\mathbf{r}}{dt} = \boldsymbol{\Omega} \times \mathbf{r}.\end{aligned} \hspace{\stretch{1}}(2.20)$

We can write this as

$\begin{aligned}\begin{bmatrix}v_x \\ v_y \\ v_z\end{bmatrix} =\left(\dot{\theta} \frac{\partial {A^{-1}}}{\partial {\theta}} A+\dot{\phi} \frac{\partial {A^{-1}}}{\partial {\phi}} A+\dot{\psi} \frac{\partial {A^{-1}}}{\partial {\psi}} A\right)\begin{bmatrix}x \\ y \\ z\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(2.21)$

Actually doing this calculation is asked of us in HW6. The final answer is

$\begin{aligned}\frac{d x_i}{dt} = \left(\dot{\phi} \epsilon_{ijk} j_n^\phi+\dot{\theta} \epsilon_{ijk} j_n^\theta+\dot{\psi} \epsilon_{ijk} j_n^\psi\right) x_j.\end{aligned} \hspace{\stretch{1}}(2.22)$

Here $\epsilon_{ijk}$ is the usual fully antisymmetric tensor with properties

$\begin{aligned}\epsilon_{ijk} =\left\{\begin{array}{l l}0 & \quad \mbox{when any of the indexes are equal.} \\ 1 & \quad \mbox{for any of$ latex ijk = 123, 231, 312$ (cyclic permutations of $123$ .} \\ -1 & \quad \mbox{for any of $ijk = 213, 132, 321$ .} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.23)$

[end] This and anything else should be ignored.# email qana907geri@post.wordpress.com -s autopost -a …

Posted in Math and Physics Learning. | Tagged: angular velocity, Euler angles, general 3D rotation, phy354 | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 18: Scaling to setup to solve the boundary layer problem. Taught by Prof. K. Das.

Posted by peeterjoot on March 21, 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.

Motivation.

We’ve been talking about impulsively started flow and the Stokes boundary problem.

We’ll now move on to a similar problem, that of fluid flow over a solid body.

Fluid flow over a solid body.

Consider figure (\ref{fig:continuumL18:continuumL18Fig1}), where we have an illustration of flow over a solid object with a boundary layer of thickness $\delta$ .

Flow over a solid object with boundary layer.

We have a couple scales to consider.

Velocity scale $U$ .
Length scale in the $y$ direction $\delta$ .
Length scale in the $x$ direction $L$ .

where

$\begin{aligned}L \gg \delta.\end{aligned} \hspace{\stretch{1}}(3.1)$

As always, we start with the Navier-Stokes equation, restricting ourselves to the steady state ${\partial {\mathbf{u}}}/{\partial {t}} = 0$ case. In coordinates, for incompressible flows, we have our usual $x$ momentum, $y$ momentum, and continuity equations

\begin{subequations}

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

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

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

\end{subequations}

Let’s look at the scaling of these equations, starting with the continuity equation 3.2c. This is roughly

$\begin{aligned}\frac{\partial {u}}{\partial {x}} &\sim \frac{U}{L} \\ \frac{\partial {v}}{\partial {y}} &\sim \frac{v}{\delta}\end{aligned} \hspace{\stretch{1}}(3.3)$

We require that these have to be of the same order of magnitude.

FIXME: why? This doesn’t make sense to me since in horizontal flow we had $v = 0$ , and the two components of the divergence are obviously of different scales.

If these are of the same scale we have

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

so that

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

$\begin{aligned}v \ll U\end{aligned} \hspace{\stretch{1}}(3.7)$

Looking at the viscous terms

$\begin{aligned}\nu \frac{\partial^2 {{u}}}{\partial {{x}}^2} &\sim \frac{\nu U}{L^2} \\ \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2} &\sim \frac{\nu U}{\delta^2}\end{aligned} \hspace{\stretch{1}}(3.8)$

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

So we can neglect the $x$ component of the Laplacian in our $x$ momentum equation 3.2a.

How about the inertial terms

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

Since these are of the same order (in the boundary regions) we cannot neglect either. We also cannot neglect the pressure gradient, since this is what induces the flow.

For the $y$ momentum equation we have

$\begin{aligned}\nu \frac{\partial^2 {{v}}}{\partial {{x}}^2} &\sim \nu \frac{\delta U}{L} \frac{1}{{L^2}} \sim \nu \frac{\delta U}{L^3} \ll \frac{\nu U}{\delta^2} \\ \nu \frac{\partial^2 {{v}}}{\partial {{y}}^2} &\sim \nu \frac{\delta U}{L} \frac{1}{{\delta^2}} \sim \nu \frac{U}{\delta L} \ll \frac{\nu U}{\delta^2}\end{aligned} \hspace{\stretch{1}}(3.13)$

We can neglect all the Laplacian terms in the $y$ momentum equation.

Question: Why compare the magnitude of the viscous terms for the $y$ momentum to the magnitude of the same terms in the $x$ momentum equation, and not to the LHS of the $y$ momentum equation.

Answer. That’s a valid point, but our equations are coupled, and contributions from one feed into the other.

We aren’t done yet. For the inertial terms in the $y$ momentum equation we have

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

Note that

$\begin{aligned}\frac{\delta U^2}{L^2} = \frac{\delta}{L} \left( \frac{U^2}{L} \right) \ll \frac{U^2}{L}\end{aligned} \hspace{\stretch{1}}(3.17)$

We see that both of the $y$ momentum inertial terms can be neglected in comparison to the $x$ momentum equations.

Putting all of this together, our equations of motion for the boundary flow are now reduced to

\begin{subequations}

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

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

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

\end{subequations}

Let’s try to see what we can do with our pressure term. Let’s start with Navier-Stokes in vector form

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

Writing the body force as a potential

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

so that we have

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

Using the vector identity

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

we can write

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

If we consider the non-viscous region of the flow (far from the boundary layer), we can kill the Laplacian term. Again, considering only the steady state, and assuming that we have irrotational flow ( $\boldsymbol{\nabla} \times \mathbf{u} = 0$ ) in the non-viscous region, we have

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

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

This is the Bernoulli equation.

Using this we can get an idea of the magnitude of the pressure term. We have

$\begin{aligned}- \frac{\partial {}}{\partial {x}} \left( \frac{p}{\rho} \right) &\sim \frac{\partial {}}{\partial {x}} \left( \frac{1}{{2}} u^2 \right) \\ &\sim \frac{2}{2} u \frac{\partial {u}}{\partial {x}} \\ &\sim U \frac{dU}{dx}\end{aligned}$

Our equations of motion are finally reduced to

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

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

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

\end{subequations}

Aside

Observe that if we operate on 3.23 with a divergence operation, then we don’t have to assume non-visous flow but in that case we can only say that (for irrotational flow) we have

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

Posted in Math and Physics Learning. | Tagged: Bernoulli equation, boundary layer, irrotational flow, Navier-Stokes, PHY454, viscous region | Leave a Comment »

Inclined flow without constant height assumption. Setting up the problem, but not actually solving it (reworking vorticity equations)

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

Pressure and vorticity equations with the non-linear term retained. (Reworked slightly)

In section 40-2 [1], the identity

$\begin{aligned}( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} = (\boldsymbol{\nabla} \times \mathbf{u}) \times \mathbf{u}+ \frac{1}{{2}} \boldsymbol{\nabla} ( \mathbf{u} \cdot \mathbf{u} ),\end{aligned} \hspace{\stretch{1}}(4.37)$

is used to put the vorticity equation into a form with one additional portion expressed as a gradient. This is a superior way to handle the inertial term because the curl of that gradient is then killed off.

I’ve expressed the curl as a wedge product, and not a cross product (either works since they related by a constant duality transformation). With the wedge product the identity 4.37 has different signs. That is

$\begin{aligned}\mathbf{u} \cdot (\boldsymbol{\nabla} \wedge \mathbf{u})&=(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}-\boldsymbol{\nabla}' (\mathbf{u}' \cdot \mathbf{u}) \\ &=(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}-\frac{1}{{2}} \boldsymbol{\nabla} (\mathbf{u} \cdot \mathbf{u}),\end{aligned}$

Here I’ve used the Hestenes overdot notation [2] to mark the operational range of the gradient $\boldsymbol{\nabla}$ (i.e. indicating that the gradient acts only on one of the $\mathbf{u}$ terms initially). That gives us

$\begin{aligned}( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} = \mathbf{u} \cdot (\boldsymbol{\nabla} \wedge \mathbf{u})+\frac{1}{{2}} \boldsymbol{\nabla} (\mathbf{u} \cdot \mathbf{u}).\end{aligned} \hspace{\stretch{1}}(4.38)$

Navier-Stokes (for incompressible flows) now takes the form

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

where in our problem we have killed the time dependence and have

$\begin{aligned}\phi = -g ( x \sin\alpha, -y \cos\alpha).\end{aligned} \hspace{\stretch{1}}(4.40)$

The divergence of $\mathbf{u} \cdot (\boldsymbol{\nabla} \wedge \mathbf{u})$ unfortunately isn’t zero. For exposition purposes, let’s write this out explicitly as a function of the vorticity components

$\begin{aligned}\Omega_{rk} = \partial_r u_k - \partial_k u_r.\end{aligned} \hspace{\stretch{1}}(4.41)$

Expanding out that divergence we have

$\begin{aligned}\boldsymbol{\nabla} \cdot (\mathbf{u} \cdot (\boldsymbol{\nabla} \wedge \mathbf{u})) &=\boldsymbol{\nabla} \cdot \left( u_r \partial_s u_t \mathbf{e}_r \cdot (\mathbf{e}_s \wedge \mathbf{e}_t) \right) \\ &=\boldsymbol{\nabla} \cdot \left( u_r \partial_s u_t (\delta_{rs} \mathbf{e}_t-\delta_{rt} \mathbf{e}_s)\right)\\ &=\partial_k\left( u_r \partial_s u_k \delta_{rs} - u_r \partial_k u_t \delta_{rt} )\right) \\ &=\partial_k \left(u_r\left( \partial_r u_k - \partial_k u_r \right)\right) \end{aligned}$

$\begin{aligned}\boldsymbol{\nabla} \cdot (\mathbf{u} \cdot (\boldsymbol{\nabla} \wedge \mathbf{u})) =\partial_k \left(u_r\Omega_{rk}\right).\end{aligned} \hspace{\stretch{1}}(4.42)$

Let’s also, for exposition, expand out the curl of this remaining non-linear term in coordinates. Being a bit smarter this time, we can avoid expressing $\Omega$ in terms of $\boldsymbol{\nabla}$ and $\mathbf{u}$ and leave it as a bivector explicitly. We have

$\begin{aligned}\boldsymbol{\nabla} \wedge (\mathbf{u} \cdot \boldsymbol{\Omega})&=\frac{1}{{2}} \boldsymbol{\nabla} \wedge ( u_m \Omega_{ab} \mathbf{e}_m \cdot (\mathbf{e}_a \wedge \mathbf{e}_b) ) \\ &=\frac{1}{{2}} \boldsymbol{\nabla} \wedge ( u_m \Omega_{ab} (\delta_{ma} \mathbf{e}_b - \delta_{mb} \mathbf{e}_a) )\\ &=\frac{1}{{2}} \boldsymbol{\nabla} \wedge ( u_a \Omega_{ab} \mathbf{e}_b-u_b \Omega_{ab} \mathbf{e}_a)\\ &=\boldsymbol{\nabla} \wedge ( u_a \Omega_{ab} \mathbf{e}_b)\end{aligned}$

This is

$\begin{aligned}\boldsymbol{\nabla} \wedge (\mathbf{u} \cdot \boldsymbol{\Omega}) = \partial_r (u_a \Omega_{ak}) \mathbf{e}_r \wedge \mathbf{e}_k.\end{aligned} \hspace{\stretch{1}}(4.43)$

The Navier-Stokes equations are now recast in terms of vorticity as

\begin{subequations}

$\begin{aligned}\boldsymbol{\nabla} \cdot (\mathbf{u} \wedge \boldsymbol{\Omega})=-\boldsymbol{\nabla}^2 \left( \frac{p}{\rho} + \frac{1}{2} \mathbf{u}^2 + \phi \right) \end{aligned} \hspace{\stretch{1}}(4.44a)$

$\begin{aligned}\frac{\partial {\boldsymbol{\Omega}}}{\partial {t}} + \boldsymbol{\nabla} \wedge (\mathbf{u} \cdot \boldsymbol{\Omega}) = \nu \boldsymbol{\nabla}^2 \boldsymbol{\Omega}\end{aligned} \hspace{\stretch{1}}(4.44b)$

$\begin{aligned}\boldsymbol{\Omega} = \boldsymbol{\nabla} \wedge \mathbf{u}\end{aligned} \hspace{\stretch{1}}(4.44c)$

$\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0,\end{aligned} \hspace{\stretch{1}}(4.44d)$

\end{subequations}

Having restated things with the $\boldsymbol{\nabla} \mathbf{u}^2$ term moved to the RHS, let’s also now write out 4.44a and 4.44b in coordinate form (we want this for the 2D case). This is

\begin{subequations}

$\begin{aligned}\partial_b \left( u_a \Omega_{a b} \right)=-\boldsymbol{\nabla}^2 \left( \frac{p}{\rho} + \frac{1}{2} \mathbf{u}^2 + \phi \right) \end{aligned} \hspace{\stretch{1}}(4.45a)$

$\begin{aligned}\frac{\partial {\Omega_{mn}}}{\partial {t}} +\partial_m (u_a \Omega_{an}) -\partial_n (u_a \Omega_{am}) =\nu \boldsymbol{\nabla}^2 \Omega_{mn}\end{aligned} \hspace{\stretch{1}}(4.45b)$

\end{subequations}

For our problem where we have only $u_1$ and $u_2$ components, and any $\partial_3$ operations are zero, we find

$\begin{aligned}\partial_b \left( u_a \Omega_{a b} \right)&=\partial_2 (u_1 \Omega_{12} ) + \partial_1( u_2 \Omega_{21} ) \\ &=(\partial_2 u_1 - \partial_1 u_2 )\Omega_12+ u_1 \partial_2 \Omega_{12} - u_2 \partial_1 \Omega_{12} \\ &=-\Omega_{12}^2 + (u_1 \partial_2 - u_2 \partial_1 ) \Omega_{12} \\ &=-\Omega_{12}^2 - i \cdot (\mathbf{u} \wedge \boldsymbol{\nabla}) \Omega_{12} \\ \end{aligned}$

and

$\begin{aligned}\partial_1 (u_a \Omega_{a 2}) - \partial_2 (u_a \Omega_{a 1})&=\partial_1 (u_1 \Omega_{1 2}) - \partial_2 (u_2 \Omega_{2 1}) \\ &=\partial_1 (u_1 \Omega_{1 2}) + \partial_2 (u_2 \Omega_{1 2}) \\ &=(\not{{\partial_1 u_1 + \partial_2 u_2}} ) \Omega_{1 2}+ (u_1 \partial_1 + u_2 \partial_2) \Omega_{1 2} \\ &=(\mathbf{u} \cdot \boldsymbol{\nabla}) \Omega_{1 2}\end{aligned}$

So we have

\begin{subequations}

$\begin{aligned}\Omega_{12}^2 + i \cdot (\mathbf{u} \wedge \boldsymbol{\nabla}) \Omega_{12}=\boldsymbol{\nabla}^2 \left( \frac{p}{\rho} + \frac{1}{2} \mathbf{u}^2 + \phi \right) \end{aligned} \hspace{\stretch{1}}(4.46a)$

$\begin{aligned}\frac{\partial {\Omega_{12}}}{\partial {t}} +(\mathbf{u} \cdot \boldsymbol{\nabla}) \Omega_{1 2}=\nu \boldsymbol{\nabla}^2 \Omega_{1 2}.\end{aligned} \hspace{\stretch{1}}(4.46b)$

\end{subequations}

Here I’ve used $i = \mathbf{e}_1 \mathbf{e}_2$ again, so that the pair of differential operators on the LHS of the respective equations above are

$\begin{aligned}i \cdot (\mathbf{u} \wedge \boldsymbol{\nabla}) &= -u_1 \partial_2 + u_2 \partial_1 \\ \mathbf{u} \cdot \boldsymbol{\nabla} &= u_1 \partial_1 + u_2 \partial_2.\end{aligned} \hspace{\stretch{1}}(4.47)$

Note that since $\boldsymbol{\nabla}^2 \phi$ could equal zero (as in our problem) we will likely have additional work to ensure that any solution that we find to this set of equations is also still a solution to our original first order Navier-Stokes equation.

Now what?

The strategy that I’d thought to attempt to tackle this problem, when I’d left it like 4.36 was something along the following lines

\begin{itemize}
\item First ignore the non-linear terms. Find solutions for the homogeneous vorticity and pressure Laplacian equations that satisfy our boundary value conditions, and use that to find a first solution for $h(x)$ .
\item Use this to solve for $u$ and $v$ from the vorticity.
\end{itemize}

However, after reworking it using the identity found in Feynman’s dry water chapter, I think it’s best not to try to solve it yet, and study some more first. I have a feeling that there are likely more such techniques that have been developed that will be useful to know before I try to plow my way through things.

Regardless, it’s interesting to see just how tricky the equations of motion become when one doesn’t make unrealistic assumptions. I have a feeling that to actually attempt this specific problem, I may very well need a computer and numerical techniques.

References

[1] R.P. Feynman, R.B. Leighton, and M.L. Sands. Feynman lectures on physics.[Lectures on physics], chapter The flow of dry water. Addison-Wesley Publishing Company. Reading, Massachusetts, 1963.

[2] D. Hestenes. New Foundations for Classical Mechanics. Kluwer Academic Publishers, 1999.

Posted in Math and Physics Learning. | Tagged: geometric algebra, Navier-Stokes, PHY454, vorticity | Leave a Comment »

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

Posted by peeterjoot on March 18, 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.

In an informal discussion after class, it was claimed that the steady state flow down a plane would have constant height, unless you bring surface tension effects into the mix. Part of that statement just doesn’t make sense to me. Consider the forces acting on the fluid in the figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig1})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.3\textheight]{inclinedFlowWithoutConstantHeightAssumptionFig1}
\caption{Gravitational force components acting on fluid flowing down a plane.}

\end{figure}

In the inclined reference frame we have a component of the force acting downwards (in the negative y-axis direction), and have a component directed down the x-axis. Wouldn’t this act to both push the fluid down the plane and push part of the fluid downwards? I’d expect this to introduce some vorticity as depicted.

While we are just about to start covered surface tension, perhaps this is just allowing the surface to vary, and then solving the Navier-Stokes equations that result. Let’s try setting up the Navier-Stokes equation for steady state viscous flow down a plane without any assumption that the height is constant and see how far we can get.

The equations of motion.

We’ll use the same coordinates as before, with the directions given as in figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig2}). However, this time, we let the height $h(x)$ of the fluid at any distance $x$ down the plane from the initial point vary.

\begin{figure}[htp]
\centering
\def\svgwidth{0.7\columnwidth}
\caption{Diagram of coordinates for inclined flow problem.}

\end{figure}

For viscous incompressible flow down the plane our equations of motion are

\begin{subequations}

$\begin{aligned}\rho \frac{\partial {u}}{\partial {t}} + \rho (u \partial_x + v \partial_y) u = -\partial_x p + \mu \left(\partial_{xx} + \partial_{yy}\right) u + \rho g \sin\alpha \end{aligned} \hspace{\stretch{1}}(2.1a)$

$\begin{aligned}\rho \frac{\partial {v}}{\partial {t}} + \rho (u \partial_x + v \partial_y) v = -\partial_y p + \mu \left(\partial_{xx} + \partial_{yy}\right) v - \rho g \cos\alpha \end{aligned} \hspace{\stretch{1}}(2.1b)$

$\begin{aligned}0 = -\partial_z p \end{aligned} \hspace{\stretch{1}}(2.1c)$

$\begin{aligned}0 = \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(2.1d)$

\end{subequations}

Now, can we kill the time dependent term? Even allowing for $u$ to vary with $x$ and introducing a non-horizontal flow component, I think that we can. If the flow at $x = 0$ is constant, not varying at all with time, I think it makes sense that we will have no time dependence in the flow for $x \ne 0$ . So, I believe that our starting point is as above, but with the time derivatives killed off. That is

\begin{subequations}

$\begin{aligned}\rho (u \partial_x + v \partial_y) u = -\partial_x p + \mu \left(\partial_{xx} + \partial_{yy}\right) u + \rho g \sin\alpha \\ \end{aligned} \hspace{\stretch{1}}(2.2a)$

$\begin{aligned}\rho (u \partial_x + v \partial_y) v = -\partial_y p + \mu \left(\partial_{xx} + \partial_{yy}\right) v - \rho g \cos\alpha \\ \end{aligned} \hspace{\stretch{1}}(2.2b)$

$\begin{aligned}0 = -\partial_z p \\ \end{aligned} \hspace{\stretch{1}}(2.2c)$

$\begin{aligned}0 = \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(2.2d)$

\end{subequations}

These don’t look particularily easy to solve, and we haven’t even set up the boundary value constraints yet. Let’s do that next.

The boundary value constraints.

One of out constraints is the no-slip condition for the velocity components at the base of the slope

$\begin{aligned}\boxed{u(x, 0) = v(x, 0) = 0.}\end{aligned} \hspace{\stretch{1}}(3.3)$

We should also have a zero tangential component to the traction vector at the interface. We need to consider some geometry, and refer to figure (\ref{fig:inclinedFlowWithoutConstantHeightAssumption:inclinedFlowWithoutConstantHeightAssumptionFig3})

\begin{figure}[htp]
\centering
\def\svgwidth{0.7\columnwidth}
\caption{Differential vector element.}

\end{figure}

A position vector on the surface has the value

$\begin{aligned}\mathbf{r} = x \hat{\mathbf{x}} + h \hat{\mathbf{y}}\end{aligned} \hspace{\stretch{1}}(3.4)$

so that a differential element on the surface, tangential to the surface is proportional to

$\begin{aligned}d\mathbf{r} = \left( \hat{\mathbf{x}} + \frac{d{{h}}}{dx} \hat{\mathbf{y}} \right) dx\end{aligned} \hspace{\stretch{1}}(3.5)$

so our unit tangent vector in the direction depicted in the figure is

$\begin{aligned}\hat{\boldsymbol{\tau}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( 1, h' \right).\end{aligned} \hspace{\stretch{1}}(3.6)$

The outwards facing normal has a value, up to a factor of plus or minus one, of

$\begin{aligned}\hat{\mathbf{n}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( h', -1 \right).\end{aligned} \hspace{\stretch{1}}(3.7)$

We can fix the orientation by considering the unit bivector

$\begin{aligned}\hat{\boldsymbol{\tau}} \wedge \hat{\mathbf{n}} &= \frac{1}{{1 + (h')^2}} \left( 1, h' \right) \wedge\left( h', -1 \right) \\ &=\begin{vmatrix}1 & h' \\ h' & -1\end{vmatrix}\mathbf{e}_1 \mathbf{e}_2 \\ &=( -1 - (h')^2 )\mathbf{e}_1 \mathbf{e}_2.\end{aligned}$

So we really want the other orientation

$\begin{aligned}\hat{\mathbf{n}} = \frac{1}{{\sqrt{1 + (h')^2}}} \left( -h', 1 \right).\end{aligned} \hspace{\stretch{1}}(3.8)$

Our traction vector relative to the normal $\hat{\mathbf{n}}$ is

$\begin{aligned}\mathbf{t} &= \mathbf{e}_i T_{ij} n_j \\ &= \mathbf{e}_i \left( -p \delta_{ij}+ \mu e_{ij}\right) n_j \\ &=-p \hat{\mathbf{n}} + \mu \mathbf{e}_i e_{ij} n_j\end{aligned}$

So the component in the tangential direction is

$\begin{aligned}\mathbf{t} \cdot \hat{\boldsymbol{\tau}} &=-p \not{{\hat{\mathbf{n}} \cdot \hat{\boldsymbol{\tau}}}} + \mu \mathbf{e}_i e_{ij} n_j tau_i \\ &=\frac{\mu}{1 + (h')^2}\begin{bmatrix}1 & h' \end{bmatrix}\begin{bmatrix}e_{11} & e_{12} \\ e_{21} & e_{22}\end{bmatrix}\begin{bmatrix}-h' \\ 1\end{bmatrix} \\ &=\frac{\mu}{1 + (h')^2}\begin{bmatrix}1 & h' \end{bmatrix}\begin{bmatrix}-h' e_{11} + e_{12} \\ -h' e_{21} + e_{22}\end{bmatrix} \\ &=\frac{\mu}{1 + (h')^2}\left(-h' e_{11} + e_{12} + h'( -h' e_{21} + e_{22} )\right) \\ &=\frac{\mu}{1 + (h')^2}\left(h' (e_{22} - e_{11} )+ e_{12} (1 - (h')^2 )\right) \end{aligned}$

Our strain tensor components, for a general 2D flow, are

$\begin{aligned}e_{11} &= \frac{\partial {u}}{\partial {x}} \\ e_{22} &= \frac{\partial {v}}{\partial {y}} \\ e_{12} &= \frac{1}{{2}} \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right).\end{aligned} \hspace{\stretch{1}}(3.9)$

So, a zero tangential traction vector component at the interface requires

$\begin{aligned}\boxed{0 = h' \left( {\left.{{\left( \frac{\partial {v}}{\partial {y}} - \frac{\partial {u}}{\partial {x}} \right) + \frac{1}{{2}} \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)}}\right\vert}_{{y = h}}\right)(1 - (h')^2 ).}\end{aligned} \hspace{\stretch{1}}(3.12)$

What is the normal component of the traction vector at the interface? We can calculate

$\begin{aligned}\mathbf{t} \cdot \hat{\mathbf{n}} &=-p \hat{\mathbf{n}} \cdot \hat{\mathbf{n}} + \mu \mathbf{e}_i e_{ij} n_j n_i \\ &=-p+\frac{\mu}{1 + (h')^2}\begin{bmatrix}-h' & 1\end{bmatrix}\begin{bmatrix}-h' e_{11} + e_{12} \\ -h' e_{21} + e_{22}\end{bmatrix} \\ &=-p+\frac{\mu}{1 + (h')^2}\left(-h'(-h' e_{11} + e_{12}) -h' e_{21} + e_{22}\right) \\ &=-p+\frac{\mu}{1 + (h')^2}\left(-2 h' e_{12} + (h')^2 e_{11} + e_{22}\right)\end{aligned}$

So this component of the traction vector is

$\begin{aligned}\mathbf{t} \cdot \hat{\mathbf{n}} =-p+\frac{\mu}{1 + (h')^2}\left(- h' \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)+ (h')^2 \frac{\partial {u}}{\partial {x}} + \frac{\partial {v}}{\partial {y}} \right)\end{aligned} \hspace{\stretch{1}}(3.13)$

For the purely recilinear flow, with $h' = 0$ and $v = 0$ , as a a consequence of NS and our assumptions, all but the pressure portion of this component of the traction vector was zero. The force balance equation for the interface was therefore just a matching of the pressure with the external (ie: air) pressure.

In this more general case we have the same thing, but the non-pressure portions of the traction vector are all zero in the medium. Outside of the fluid (in the air say), we have assumed no motion, so this force balance condition becomes

$\begin{aligned}{\left.{{\mathbf{t} \cdot \hat{\mathbf{n}}}}\right\vert}_{{\text{fluid}}}={\left.{{\mathbf{t} \cdot \hat{\mathbf{n}}}}\right\vert}_{{\text{air}}}.\end{aligned} \hspace{\stretch{1}}(3.14)$

Again assuming no motion of the air, with air pressure $p_A$ , this is

$\begin{aligned}\boxed{-p(x, h)+{\left.{{\frac{\mu}{1 + (h')^2}\left(- h' \left( \frac{\partial {v}}{\partial {x}} +\frac{\partial {u}}{\partial {y}}\right)+ (h')^2 \frac{\partial {u}}{\partial {x}} + \frac{\partial {v}}{\partial {y}} \right)}}\right\vert}_{{y = h}}= -p_A}\end{aligned} \hspace{\stretch{1}}(3.15)$

Observe that for the horizontal flow problem, where $h' = 0$ and $v = 0$ , this would have been nothing more than a requirement that $p(h) = p_A$ , but now that we introduce downwards motion and allow the height to vary, the pressure matching condition becomes a much more complex beastie.

Laplacian of Pressure and Vorticity.

Supposing that we are neglecting the non-linear term of the Navier-Stokes equation. For incompressible steady state flow, without any external forces, we would then have

\begin{subequations}

$\begin{aligned}0 = -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} \end{aligned} \hspace{\stretch{1}}(4.16a)$

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

\end{subequations}

How do we actually solve this beastie?

Separation of variables?

Considering this in 2D, assuming no z-dependence, with $\mathbf{u} = \mathbf{u}(x, y) = (u, v)$ we have

$\begin{aligned}0 &= -\partial_x p + \mu (\partial_{xx} + \partial_{yy} )u \\ 0 &= -\partial_y p + \mu (\partial_{xx} + \partial_{yy} )v \\ 0 &= \partial_x u + \partial_y v.\end{aligned} \hspace{\stretch{1}}(4.17)$

Attempting separation of variables seems like something reasonable to try. With

$\begin{aligned}u &= X(x) Y(y) \\ v &= R(x) S(y) \\ p &= P(x) Q(y)\end{aligned} \hspace{\stretch{1}}(4.20)$

we get

$\begin{aligned}0 &= -P' Q + \mu (X'' Y + X Y'') \\ 0 &= -P Q' + \mu (R'' S + R S'') \\ 0 &= X' Y + R S'\end{aligned} \hspace{\stretch{1}}(4.23)$

I couldn’t find a way to substitute any of these into the other that would allow me to separate them, but perhaps I wasn’t clever enough.

In terms of vorticity?

The idea of substuting the zero divergence equation $\boldsymbol{\nabla} \cdot \mathbf{u}$ will clearly lead to something a bit simpler. Treating the Laplacian as a geometric (clifford) product of two gradients we have

$\begin{aligned}0 &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u} \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \boldsymbol{\nabla} \mathbf{u} ) \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \not{{\boldsymbol{\nabla} \cdot \mathbf{u}}} + \boldsymbol{\nabla} \wedge \mathbf{u} ) \\ &= -\boldsymbol{\nabla} p + \mu \boldsymbol{\nabla} ( \boldsymbol{\nabla} \wedge \mathbf{u} ) \\ &= \boldsymbol{\nabla} \left( -p + \mu \boldsymbol{\nabla} \wedge \mathbf{u} \right)\end{aligned}$

Writing out the vorticity (bivector) in component form, and writing $i = \mathbf{e}_1 \wedge \mathbf{e}_2 = \mathbf{e}_1 \mathbf{e}_2$ for the 2D pseudoscalar, we have

$\begin{aligned}\boldsymbol{\nabla} \wedge \mathbf{u} &= (\mathbf{e}_1 \partial_x + \mathbf{e}_2 \partial_y) \wedge (\mathbf{e}_1 u + \mathbf{e}_2 v) \\ &= \mathbf{e}_1 \mathbf{e}_2 (\partial_x v - \partial_y u ) \\ &= i (\partial_x v - \partial_y u )\end{aligned}$

It seems natural to write

$\begin{aligned}\Theta = \partial_x v - \partial_y u,\end{aligned} \hspace{\stretch{1}}(4.26)$

so that Navier-Stokes takes the form

$\begin{aligned}0 = \boldsymbol{\nabla} (-p + i \mu \Theta ).\end{aligned} \hspace{\stretch{1}}(4.27)$

Operating on this from the left with another gradient we find that this combination of pressure and vorticity must satisfy the following multivector Laplacian equation

$\begin{aligned}0 = \boldsymbol{\nabla}^2 (-p + i \mu \Theta ).\end{aligned} \hspace{\stretch{1}}(4.28)$

However, note that $\boldsymbol{\nabla}^2$ is a scalar operator, and this zero identity has both scalar and pure bivector components. Both must separately equal zero

\begin{subequations}

$\begin{aligned}0 = \boldsymbol{\nabla}^2 p \end{aligned} \hspace{\stretch{1}}(4.29a)$

$\begin{aligned}0 = \boldsymbol{\nabla}^2 \Theta.\end{aligned} \hspace{\stretch{1}}(4.29b)$

\end{subequations}

Note that we can obtain 4.29 much more directly, if we know that’s what we want to do. Just operate on 4.16a with the gradient from the left right off the bat. We find

$\begin{aligned}0 &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^3 \mathbf{u} \\ &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \mathbf{u}) \\ &= -\boldsymbol{\nabla}^2 p + \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \wedge \mathbf{u}) \\ \end{aligned}$

Again, we have a multivector equation scalar and bivector parts, that must separately equal zero. With the magnitude of the vorticity $\Theta$ given by 4.26, we once again obtain 4.29. This can be done in plain old vector algebra as well by operating on the left not by the gradient, but separately with a divergence and curl operator.

Pressure and vorticity equations with the non-linear term retained.

If we add back in our body force, and assume that it is constant (i.e. gravity), then this this will get killed with the application of the gradient. We’ll still end up with one Laplacian for pressure, and one for vorticity. That’s not the case for the inertial $(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}$ term of Navier-Stokes. Let’s take the divergence and curl of this and see how we have to modify the Laplacian equations above.

Starting with the divergence, with summation implied over repeated indexes, we have

$\begin{aligned}\boldsymbol{\nabla} \cdot ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u})&=\partial_k ( \mathbf{u} \cdot \boldsymbol{\nabla} u_k ) \\ &=\partial_k ( u_m \partial_m u_k ) \\ &=(\partial_k u_m) (\partial_m u_k ) + u_m \partial_m \partial_k u_k \\ &=\sum_k (\partial_k u_k)^2+\sum_{k \ne m} (\partial_k u_m) (\partial_m u_k ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) ( \boldsymbol{\nabla} \cdot \mathbf{u} )\end{aligned}$

So we have

$\begin{aligned}\boldsymbol{\nabla} \cdot ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u})=\sum_k (\partial_k u_k)^2+2 \sum_{k < m} (\partial_k u_m) (\partial_m u_k ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) ( \boldsymbol{\nabla} \cdot \mathbf{u} )\end{aligned} \hspace{\stretch{1}}(4.30)$

We are working with the $\boldsymbol{\nabla} \cdot \mathbf{u} = 0$ incompressibility assumption so we kill off the last term. Our velocity ends up introducing a non-homogeneous forcing term to the Laplacian pressure equation and we’ve got something trickier to solve

$\begin{aligned}\boxed{\rho \sum_k (\partial_k u_k)^2+2 \rho \sum_{k < m} (\partial_k u_m) (\partial_m u_k ) =-\boldsymbol{\nabla}^2 p.}\end{aligned} \hspace{\stretch{1}}(4.31)$

Now let’s see how our vorticity Laplacian needs to be modified. Taking the curl of the implusive term we have

$\begin{aligned}\boldsymbol{\nabla} \wedge ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}) &=\mathbf{e}_k \partial_k \wedge ( u_m \partial_m u_n \mathbf{e}_n ) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \partial_k ( u_m \partial_m u_n ) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \left( (\partial_k u_m) (\partial_m u_n ) +u_m \partial_m \partial_k u_n \right) \\ &=(\mathbf{e}_k \wedge \mathbf{e}_n) \left( (\partial_k u_m) (\partial_m u_n ) + (\mathbf{u} \cdot \boldsymbol{\nabla}) \partial_k u_n\right) \end{aligned}$

So we have

$\begin{aligned}\boldsymbol{\nabla} \wedge ((\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}) =(\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u})+ (\mathbf{u} \cdot \boldsymbol{\nabla}) (\boldsymbol{\nabla} \wedge \mathbf{u})\end{aligned} \hspace{\stretch{1}}(4.32)$

Putting things back together, our vorticity equation is

$\begin{aligned}\rho (\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u})+ \rho (\mathbf{u} \cdot \boldsymbol{\nabla}) (\boldsymbol{\nabla} \wedge \mathbf{u})= \mu \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \wedge \mathbf{u})\end{aligned} \hspace{\stretch{1}}(4.33)$

Or, with

$\begin{aligned}\boldsymbol{\omega} = \boldsymbol{\nabla} \wedge \mathbf{u},\end{aligned} \hspace{\stretch{1}}(4.34)$

this is

$\begin{aligned}\boxed{(\boldsymbol{\nabla} u_m) \wedge (\partial_m \mathbf{u}) + (\mathbf{u} \cdot \boldsymbol{\nabla}) \boldsymbol{\omega} = \nu \boldsymbol{\nabla}^2 \boldsymbol{\omega}.}\end{aligned} \hspace{\stretch{1}}(4.35)$

It is this and 4.31 that we really have to solve. Before moving on, let’s write out all the non-boundary condition equations in coordinate form for the 2D case that we are interested in here. We have

\begin{subequations}

$\begin{aligned}\rho \left( \left(\frac{\partial {u}}{\partial {x}}\right)^2 +\left(\frac{\partial {u}}{\partial {y}}\right)^2 + 2 \frac{\partial {u}}{\partial {y}} \frac{\partial {v}}{\partial {x}} \right) = -\boldsymbol{\nabla}^2 p\end{aligned} \hspace{\stretch{1}}(4.36a)$

$\begin{aligned}2 \left( \frac{\partial {u}}{\partial {x}} \frac{\partial {v}}{\partial {y}}+\frac{\partial {v}}{\partial {x}} \frac{\partial {v}}{\partial {y}}\right)+ \left( u \frac{\partial {}}{\partial {x}} + v \frac{\partial {}}{\partial {y}} \right) \Theta=\nu \boldsymbol{\nabla}^2 \Theta\end{aligned} \hspace{\stretch{1}}(4.36b)$

$\begin{aligned}\Theta = \frac{\partial {v}}{\partial {x}} - \frac{\partial {u}}{\partial {y}}\end{aligned} \hspace{\stretch{1}}(4.36c)$

\end{subequations}

Our solution has to satisfy these equations, as well as still satisfying the original Navier-Stokes system 2.2 that includes our gravitational term, and also has to satisfy both of our boundary value constraints 3.12, 3.15, and have $u(x, 0) = v(x, 0) = 0$ . Wow, what a mess! And this is all just the steady state problem. Imagine adding time into the mix too!

Now what?

The strategy that I’d thought to attempt to tackle a problem with is something along the following lines

\begin{itemize}
\item First ignore the non-linear terms. Find solutions for the homogenous vorticity and pressure Laplacian equations that satisfy our boundary value conditions, and use that to find a first solution for $h(x)$ .
\item Use this to solve for $u$ and $v$ from the vorticity.
\end{itemize}

However, after having gotten this far, just setting up the problem for solution, I think that it’s perhaps best not to try to solve it yet, and study some more first. I have a feeling that there are likely a number of techniques that have been developed that will likely be useful to know before I try to plow my way through things. It’s interesting to see just how tricky the equations of motion become when one doesn’t make unrealistic assumptions. I have a feeling that to actually attempt this specific problem, I may very well need a computer and numerical techniques.

Posted in Math and Physics Learning. | Tagged: boundary value constraints, Navier-Stokes, PHY454, traction vector, vorticity vector | Leave a Comment »

PHY454H1S Continuum mechanics midterm reflection.

Posted by peeterjoot on March 17, 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.

I didn’t manage my time well enough on the midterm to complete it (and also missed one easy part of the second question). For later review purposes, here is either what I answered, or what I think I should have answered for these questions.

Problem 1.

$\mathbf{P}$ -waves, $\mathbf{S}$ -waves, and Love-waves.

\begin{itemize}
\item Show that in $\mathbf{P}$ -waves the divergence of the displacement vector represents a measure of the relative change in the volume of the body.

Answer.
The $\mathbf{P}$ -wave equation was a result of operating on the displacement equation with the divergence operator

$\begin{aligned}\boldsymbol{\nabla} \cdot \left( \rho \frac{\partial^2 {\mathbf{e}}}{\partial {{t}}^2} = (\lambda + \mu) \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{e}) + \mu \boldsymbol{\nabla}^2 \mathbf{e}\right)\end{aligned} \hspace{\stretch{1}}(2.1)$

we obtain

$\begin{aligned}\frac{\partial^2 {{}}}{\partial {{t}}^2} \left( \boldsymbol{\nabla} \cdot \mathbf{e} \right) = \frac{\lambda + 2 \mu}{\rho} \boldsymbol{\nabla}^2 (\boldsymbol{\nabla} \cdot \mathbf{e}).\end{aligned} \hspace{\stretch{1}}(2.2)$

We have a wave equation where the “waving” quantity is $\Theta = \boldsymbol{\nabla} \cdot \mathbf{e}$ . Explicitly

$\begin{aligned}\Theta &= \boldsymbol{\nabla} \cdot \mathbf{e} \\ &= \frac{\partial {e_1}}{\partial {x}}+\frac{\partial {e_2}}{\partial {y}}+\frac{\partial {e_3}}{\partial {z}}\end{aligned}$

Recall that, in a coordinate basis for which the strain $e_{ij}$ is diagonal we have

$\begin{aligned}dx' &= \sqrt{1 + 2 e_{11}} dx \\ dy' &= \sqrt{1 + 2 e_{22}} dy \\ dz' &= \sqrt{1 + 2 e_{33}} dz.\end{aligned} \hspace{\stretch{1}}(2.3)$

Expanding in Taylor series to $O(1)$ we have for $i = 1, 2, 3$ (no sum)

$\begin{aligned}dx_i' \approx (1 + e_{ii}) dx_i.\end{aligned} \hspace{\stretch{1}}(2.6)$

so the displaced volume is

$\begin{aligned}dV' &= dx_1dx_2dx_3(1 + e_{11})(1 + e_{22})(1 + e_{33}) \\ &=dx_1dx_2dx_3( 1 + e_{11} + e_{22} + e_{33} + O(e_{kk}^2) )\end{aligned}$

Since

$\begin{aligned}e_{11} &= \frac{1}{{2}} \left( \frac{\partial {e_1}}{\partial {x}} +\frac{\partial {e_1}}{\partial {x}} \right) = \frac{\partial {e_1}}{\partial {x}} \\ e_{22} &= \frac{1}{{2}} \left( \frac{\partial {e_2}}{\partial {y}} +\frac{\partial {e_2}}{\partial {y}} \right) = \frac{\partial {e_2}}{\partial {y}} \\ e_{33} &= \frac{1}{{2}} \left( \frac{\partial {e_3}}{\partial {z}} +\frac{\partial {e_3}}{\partial {z}} \right) = \frac{\partial {e_3}}{\partial {z}}\end{aligned} \hspace{\stretch{1}}(2.7)$

We have

$\begin{aligned}dV' = (1 + \boldsymbol{\nabla} \cdot \mathbf{e}) dV,\end{aligned} \hspace{\stretch{1}}(2.10)$

$\begin{aligned}\frac{dV' - dV}{dV} = \boldsymbol{\nabla} \cdot \mathbf{e}\end{aligned} \hspace{\stretch{1}}(2.11)$

The relative change in volume can therefore be expressed as the divergence of $\mathbf{e}$ , the displacement vector, and it is this relative volume change that is “waving” in the $\mathbf{P}$ -wave equation as illustrated in the following (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig1}) sample 1D compression wave

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFig1}
\caption{A 1D compression wave.}
\end{figure}

\item Between a $\mathbf{P}$ -wave and an $\mathbf{S}$ -wave which one is longitudinal and which one is transverse?

Answer.
$\mathbf{P}$ -waves are longitudinal.

$\mathbf{S}$ -waves are transverse.

\item Whose speed is higher?

Answer.
From the formula sheet we have

$\begin{aligned}\left( \frac{c_L}{c_T} \right)^2 &= \frac{ \lambda + 2 \mu}{\rho} \frac{\rho}{\mu} \\ &= \frac{\lambda}{\mu} + 2 \\ &> 1\end{aligned}$

so $\mathbf{P}$ -waves travel faster than $\mathbf{S}$ -waves.

\item Is Love wave a body wave or a surface wave?

Answer.

Love waves are surface waves, traveling in a medium that can slide on top of another surface. These are characterized by vorticity rotating backwards compared to the direction of propagation as shown in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig2})

\begin{figure}[htp]
\centering
\def\svgwidth{0.6\columnwidth}
\caption{Love wave illustrated.}
\end{figure}
\end{itemize}

(b) constitutive relation, Newtonian fluids, and no-slip conditions.

\begin{itemize}
\item In continuum mechanics what do you mean by \textit{constitutive relation}?
Answer. The constitutive relation is the stress-strain relation, generally

$\begin{aligned}\sigma_{ij} = c_{abij} e_{ab}\end{aligned} \hspace{\stretch{1}}(2.12)$

for isotropic solids we model this as

$\begin{aligned}\sigma_{ij} = \lambda e_{kk} \delta_{ij} + 2 \mu e_{ij}\end{aligned} \hspace{\stretch{1}}(2.13)$

and for Newtonian fluids

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

\item What is the definition of a Non-Newtonian fluid?
Answer.

A non-Newtonian fluid would be one with a more general constitutive relationship.

Grading note. I lost a mark here. I think the answer that was being looked for (as in [1]) was that a Newtonian fluid is one with a linear stress strain relationship, and a non-Newtonian fluid would be one with a non-linear relationship. According to [2] an example of a non-Newtonian material that we are all familiar with is Silly Putty. This linearity is also how a Newtonian fluid was defined in the notes, but I didn’t remember that (this isn’t really something we use since we assume all fluids and materials are Newtonian in any calculations that we do).

\item What do you mean by \emph{no-slip} boundary condition at a fluid-fluid interface?

Exam time management note. Somehow in my misguided attempt to be complete, I missed this question amongst the rest of my verbosity).

Answer.
The no slip boundary condition is just one of velocity matching. At a non-moving boundary, the no-slip condition means that we’ll require the fluid to also have no velocity (ie. at that interface the fluid isn’t slipping over the surface). Between two fluids, this is a requirement that the velocities of both fluids match at that point (and all the rest of the points along the region of the interaction.)

\item Write down the continuity equation for an incompressible fluid.
Answer.

An incompressible fluid has

$\begin{aligned}\frac{d\rho}{dt} = 0,\end{aligned} \hspace{\stretch{1}}(2.15)$

but since we also have

$\begin{aligned}0 &=\frac{d\rho}{dt} \\ &= - \rho (\boldsymbol{\nabla} \cdot \mathbf{u}) \\ &= \frac{\partial {\rho}}{\partial {t}} + (\mathbf{u} \cdot \boldsymbol{\nabla}) \rho \\ &= 0.\end{aligned}$

A consequence is that $\boldsymbol{\nabla} \cdot \mathbf{u} = 0$ for an incompressible fluid. Let’s recall where this statement comes from. Looking at mass conservation, the rate that mass leaves a volume can be expressed as

$\begin{aligned}\frac{dm}{dt}&= \int \frac{d\rho}{dt} dV \\ &= -\int_{\partial V} \rho \mathbf{u} \cdot d\mathbf{A} \\ &= -\int_V \boldsymbol{\nabla} \cdot (\rho \mathbf{u}) dV\end{aligned}$

(the minus sign here signifying that the mass is leaving the volume through the surface, and that we are using an outwards facing normal on the volume.)

If the surface bounding the volume doesn’t change with time (ie. ${\partial {V}}/{\partial {t}} = 0$ ) we can write

$\begin{aligned}\frac{\partial {}}{\partial {t}} \int \rho dV = -\int \boldsymbol{\nabla} \cdot (\rho \mathbf{u}) dV,\end{aligned} \hspace{\stretch{1}}(2.16)$

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

so that in differential form we have

$\begin{aligned}0 = \frac{\partial {\rho}}{\partial {t}} + \boldsymbol{\nabla} \cdot (\rho \mathbf{u}).\end{aligned} \hspace{\stretch{1}}(2.18)$

Expanding the divergence by chain rule we have

$\begin{aligned}\frac{\partial {\rho}}{\partial {t}} +\mathbf{u} \cdot \boldsymbol{\nabla} \rho = -\rho \boldsymbol{\nabla} \cdot \mathbf{u},\end{aligned} \hspace{\stretch{1}}(2.19)$

but this is just

$\begin{aligned}\frac{d\rho}{dt} = -\rho \boldsymbol{\nabla} \cdot \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(2.20)$

So, for an incompressible fluid (one for which $d\rho/dt =0$ ), we must also have $\boldsymbol{\nabla} \cdot \mathbf{u} = 0$ .

\end{itemize}

Problem 2.

Statement.

Consider steady simple shearing flow $\mathbf{u} = \hat{\mathbf{x}} u(y)$ as shown in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFigQ1}) with imposed constant pressure gradient ( $G = -dp/dx$ ), $G$ being a positive number, of a single layer fluid with viscosity $\mu$ .

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFigQ1}
\caption{Shearing flow with pressure gradient and one moving boundary.}
\end{figure}

The boundary conditions are no-slip at the lower plate ( $y = h$ ). The top plate is moving with a velocity $-U$ at $y = h$ and fluid is sticking to it, so $u(h) = -U$ , $U$ being a positive number. Using the Navier-Stokes equation.

\begin{itemize}
\item Derive the velocity profile of the fluid.

Answer

Our equations of motion are

\begin{subequations}

$\begin{aligned}0 = \boldsymbol{\nabla} \cdot \mathbf{u}\end{aligned} \hspace{\stretch{1}}(3.21a)$

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

\end{subequations}

Here, we’ve used the steady state condition and are neglecting gravity, and kill off our mass compression term with the incompressibility assumption. In component form, what we have left is

$\begin{aligned}0 &= \partial_x u \\ u \not{{\partial_x u}} &= -\partial_x p + \mu \boldsymbol{\nabla}^2 u \\ 0 &= -\partial_y p \\ 0 &= -\partial_z p\end{aligned} \hspace{\stretch{1}}(3.22)$

with $\partial_y p = \partial_z p = 0$ , we must have

$\begin{aligned}\frac{\partial {p}}{\partial {x}} = \frac{dp}{dx} = -G,\end{aligned} \hspace{\stretch{1}}(3.26)$

which leaves us with just

$\begin{aligned}0 &= G + \mu \boldsymbol{\nabla}^2 u(y) \\ &= G + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} \\ &= G + \mu \frac{d^2 u}{dy^2}.\end{aligned}$

Having dropped the partials we really just want to integrate our very simple ODE a couple times

$\begin{aligned}u'' = -\frac{G}{\mu}.\end{aligned} \hspace{\stretch{1}}(3.27)$

Integrate once

$\begin{aligned}u' = -\frac{G}{\mu} y + \frac{A}{h},\end{aligned} \hspace{\stretch{1}}(3.28)$

and once more to find the velocity

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

Let’s incorporate an additional constant into $B'$

$\begin{aligned}B' = \frac{G}{2 \mu} h^2 + B\end{aligned} \hspace{\stretch{1}}(3.30)$

so that we have

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

(I didn’t do use $B'$ this way on the exam, nor did I include the factor of $1/h$ in the first integration constant, but both of these should simplify the algebra since we’ll be evaluating the boundary value conditions at $y = \pm h$ .)

$\begin{aligned}u = \frac{G}{2 \mu} (h^2 - y^2) + \frac{A}{h} y + B\end{aligned} \hspace{\stretch{1}}(3.32)$

Applying the velocity matching conditions we have for the lower and upper plates respectively

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

Adding these we find

$\begin{aligned}B = -\frac{U}{2}\end{aligned} \hspace{\stretch{1}}(3.35)$

and subtracting find

$\begin{aligned}A = -\frac{U}{2}.\end{aligned} \hspace{\stretch{1}}(3.36)$

Our velocity is

$\begin{aligned}u = \frac{G}{2 \mu} (h^2 - y^2) - \frac{U}{2 h} y -\frac{U}{2}\end{aligned} \hspace{\stretch{1}}(3.37)$

or rearranged a bit

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

\item Draw the velocity profile with the direction of the flow of the fluid when $U = 0$ , $G \ne 0$ .

Answer

With $U = 0$ our velocity has a simple parabolic profile with a max of $\frac{G}{2 \mu} (h^2 - y^2)$ at $y = 0$

$\begin{aligned}u(y) = \frac{G}{2 \mu} (h^2 - y^2).\end{aligned} \hspace{\stretch{1}}(3.39)$

This is plotted in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig3})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFig3}
\caption{Parabolic velocity profile.}
\end{figure}

\item Draw the velocity profile with the direction of the flow of the fluid when $G = 0$ , $U \ne 0$ .

Answer

With $G = 0$ , we have a plain old shear flow

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

This is linear with minimum velocity $u = 0$ at $y = -h$ , and a maximum of $-U$ at $y = h$ . This is plotted in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig4})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFig4}
\caption{Shear flow.}
\end{figure}

\item Using linear superposition draw the velocity profile of the fluid with the direction of flow qualitatively when $U \ne 0$ , $G \ne 0$ . (i) low $U$ , (ii) large $U$ .
Exam time management note. Somehow I missed this question when I wrote the exam … I figured this out right at the end when I’d run out of time by being too verbose elsewhere. I’m really not very good at writing exams in tight time constraints anymore.

Answer

For low $U$ we’ll let the parabolic dominate, and can graphically add these two as in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig5})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFig5}
\caption{Superposition of shear and parabolic flow (low $U$ )}
\end{figure}
For high $U$ , we’ll let the shear flow dominate, and have plotted this in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig6})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumMidtermReflectionFig6}
\caption{Superposition of shear and parabolic flow (high $U$ )}
\end{figure}

\item Calculate the maximum speed when $U \ne 0$ , $G \ne 0$ .
Answer

Since our acceleration is

$\begin{aligned}\frac{du}{dy} = -\frac{G}{\mu} y - \frac{U}{2 h}\end{aligned} \hspace{\stretch{1}}(3.41)$

our extreme values occur at

$\begin{aligned}y_m = -\frac{U \mu}{2 h G}.\end{aligned} \hspace{\stretch{1}}(3.42)$

At this point, our velocity is

$\begin{aligned}u(y_m) &= \frac{G}{2 \mu} \left(h^2 - \left( \frac{U \mu}{2 h G} \right)^2\right) - \frac{U}{2} \left( 1 -\frac{U \mu}{2 h^2 G}\right) \\ &=\frac{G h^2}{2 \mu} -\frac{U}{2}+ \frac{U^2 \mu}{4 h^2 G} \left(1 -\frac{1}{{2}}\right)\end{aligned}$

or just

$\begin{aligned}u_{\text{max}} = \frac{G h^2}{2 \mu} -\frac{U}{2} + \frac{U^2 \mu}{8 h^2 G}.\end{aligned} \hspace{\stretch{1}}(3.43)$

\item Calculate the flux (the volume flow rate) when $U \ne 0$ , $G \ne 0$ .
Answer

An element of our volume flux is

$\begin{aligned}\frac{dV}{dt} = dy dz \mathbf{u} \cdot \hat{\mathbf{x}}\end{aligned} \hspace{\stretch{1}}(3.44)$

Looking at the volume flux through the width $\Delta z$ is then

$\begin{aligned}\text{Flux} &= \int_0^{\Delta z} dz \int_{-h}^h dy u(y) \\ &= \Delta z \int_{-h}^h dy \frac{G}{2 \mu} (h^2 - y^2) - \frac{U}{2} \left( 1 + \frac{y}{h} \right) \\ &= \Delta z \int_{-h}^h dy \frac{G}{2 \mu} \left(h^2 y - \frac{1}{{3}} y^3 \right) - \frac{U}{2} \left( y + \frac{y^2}{2 h} \right) \\ &= \Delta z \left( \frac{2 G h^3}{3 \mu} - U h \right)\end{aligned}$

\item Calculate the mean speed when $U \ne 0$ , $G \ne 0$ .
Answer

Exam time management note. I squandered too much time on other stuff and didn’t get to this part of the problem (which was unfortunately worth a lot). This is how I think it should have been answered.

We’ve done most of the work above, and just have to divide the flux by $2 h \Delta z$ . That is

$\begin{aligned}\left\langle{{u}}\right\rangle = \frac{G h^2}{3 \mu} - \frac{U}{2}.\end{aligned} \hspace{\stretch{1}}(3.45)$

\item Calculate the tangential force (per unit width) $F_x$ on the strip $0 \le x \le L$ of the wall $y = -h$ when $U \ne 0$ , $G \ne 0$ .
Answer

Our traction vector is

$\begin{aligned}T_1 &= \sigma_{1j} n_j \\ &= \left( -p \delta_{1j} + 2 \mu e_{1j} \right) \delta_{2j} \\ &= 2 \mu e_{12} \\ &= \mu \left( \frac{\partial {u}}{\partial {y}}+\not{{\frac{\partial {v}}{\partial {x}}}}\right)\end{aligned}$

So the $\hat{\mathbf{x}}$ directed component of the traction vector is just

$\begin{aligned}T_1 = \mu \frac{\partial {u}}{\partial {y}}.\end{aligned} \hspace{\stretch{1}}(3.46)$

We’ve calculated that derivative above in 3.41, so we have

$\begin{aligned}T_1 &= \mu \left( -\frac{G}{\mu} y - \frac{U}{2 h} \right) \\ &= - G y - \frac{U \mu}{2 h} \end{aligned}$

so at $y = -h$ we have

$\begin{aligned}T_1(-h) = G h - \frac{U \mu}{2 h}.\end{aligned} \hspace{\stretch{1}}(3.47)$

To see the contribution of this force on the lower wall over an interval of length $L$ we integrate, but this amounts to just multiplying by the length of the segment of the wall

$\begin{aligned}\int_0^L T_1(-h) dx = \left( G h - \frac{U \mu}{2 h} \right) L.\end{aligned} \hspace{\stretch{1}}(3.48)$

\end{itemize}

References

[1] Wikipedia. Newtonian fluid — wikipedia, the free encyclopedia [online]. 2011. [Online; accessed 17-March-2012]. http://en.wikipedia.org/w/index.php?title=Newtonian_fluid&oldid=460346447.

[2] Wikipedia. Non-newtonian fluid — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 17-March-2012]. http://en.wikipedia.org/w/index.php?title=Non-Newtonian_fluid&oldid=479813690.

Posted in Math and Physics Learning. | Tagged: constituative relation, continuity equation, incompressible fluid, Love wave, Navier-Stokes, Newtonian fluid, no slip condition, p-wave, PHY454, pressure gradient, s-wave, shear flow, strain, stress, superposition, traction vector, volume flux | Leave a Comment »

Geometric Algebra. The very quickest introduction.

Posted by peeterjoot on March 17, 2012

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

Motivation.

An attempt to make a relatively concise introduction to Geometric (or Clifford) Algebra. Much more complete introductions to the subject can be found in [1], [2], and [3].

Axioms

We have a couple basic principles upon which the algebra is based

Vectors can be multiplied.
The square of a vector is the (squared) length of that vector (with appropriate generalizations for non-Euclidean metrics).
Vector products are associative (but not necessarily commutative).

That’s really all there is to it, and the rest, paraphrasing Feynman, can be figured out by anybody sufficiently clever.

By example. The 2D case.

Consider a 2D Euclidean space, and the product of two vectors $\mathbf{a}$ and $\mathbf{b}$ in that space. Utilizing a standard orthonormal basis $\{\mathbf{e}_1, \mathbf{e}_2\}$ we can write

$\begin{aligned}\mathbf{a} &= \mathbf{e}_1 x_1 + \mathbf{e}_2 x_2 \\ \mathbf{b} &= \mathbf{e}_1 y_1 + \mathbf{e}_2 y_2,\end{aligned} \hspace{\stretch{1}}(3.1)$

and let’s write out the product of these two vectors $\mathbf{a} \mathbf{b}$ , not yet knowing what we will end up with. That is

$\begin{aligned}\mathbf{a} \mathbf{b} &= (\mathbf{e}_1 x_1 + \mathbf{e}_2 x_2 )( \mathbf{e}_1 y_1 + \mathbf{e}_2 y_2 ) \\ &= \mathbf{e}_1^2 x_1 y_1 + \mathbf{e}_2^2 x_2 y_2+ \mathbf{e}_1 \mathbf{e}_2 x_1 y_2 + \mathbf{e}_2 \mathbf{e}_1 x_2 y_1\end{aligned}$

From axiom 2 we have $\mathbf{e}_1^2 = \mathbf{e}_2^2 = 1$ , so we have

$\begin{aligned}\mathbf{a} \mathbf{b} = x_1 y_1 + x_2 y_2 + \mathbf{e}_1 \mathbf{e}_2 x_1 y_2 + \mathbf{e}_2 \mathbf{e}_1 x_2 y_1.\end{aligned} \hspace{\stretch{1}}(3.3)$

We’ve multiplied two vectors and ended up with a scalar component (and recognize that this part of the vector product is the dot product), and a component that is a “something else”. We’ll call this something else a bivector, and see that it is characterized by a product of non-colinear vectors. These products $\mathbf{e}_1 \mathbf{e}_2$ and $\mathbf{e}_2 \mathbf{e}_1$ are in fact related, and we can see that by looking at the case of $\mathbf{b} = \mathbf{a}$ . For that we have

$\begin{aligned}\mathbf{a}^2 &=x_1 x_1 + x_2 x_2 + \mathbf{e}_1 \mathbf{e}_2 x_1 x_2 + \mathbf{e}_2 \mathbf{e}_1 x_2 x_1 \\ &={\left\lvert{\mathbf{a}}\right\rvert}^2 +x_1 x_2 ( \mathbf{e}_1 \mathbf{e}_2 + \mathbf{e}_2 \mathbf{e}_1 )\end{aligned}$

Since axiom (2) requires our vectors square to equal its (squared) length, we must then have

$\begin{aligned}\mathbf{e}_1 \mathbf{e}_2 + \mathbf{e}_2 \mathbf{e}_1 = 0,\end{aligned} \hspace{\stretch{1}}(3.4)$

$\begin{aligned}\mathbf{e}_2 \mathbf{e}_1 = -\mathbf{e}_1 \mathbf{e}_2.\end{aligned} \hspace{\stretch{1}}(3.5)$

We see that Euclidean orthonormal vectors anticommute. What we can see with some additional study is that any colinear vectors commute, and in Euclidean spaces (of any dimension) vectors that are normal to each other anticommute (this can also be taken as a definition of normal).

We can now return to our product of two vectors 3.3 and simplify it slightly

$\begin{aligned}\mathbf{a} \mathbf{b} = x_1 y_1 + x_2 y_2 + \mathbf{e}_1 \mathbf{e}_2 (x_1 y_2 - x_2 y_1).\end{aligned} \hspace{\stretch{1}}(3.6)$

The product of two vectors in 2D is seen here to have one scalar component, and one bivector component (an irreducible product of two normal vectors). Observe the symmetric and antisymmetric split of the scalar and bivector components above. This symmetry and antisymmetry can be made explicit, introducing dot and wedge product notation respectively

$\begin{aligned}\mathbf{a} \cdot \mathbf{b} &= \frac{1}{{2}}( \mathbf{a} \mathbf{b} + \mathbf{b} \mathbf{a}) = x_1 y_1 + x_2 y_2 \\ \mathbf{a} \wedge \mathbf{b} &= \frac{1}{{2}}( \mathbf{a} \mathbf{b} - \mathbf{b} \mathbf{a}) = \mathbf{e}_1 \mathbf{e}_2 (x_1 y_y - x_2 y_1).\end{aligned} \hspace{\stretch{1}}(3.7)$

so that the vector product can be written as

$\begin{aligned}\mathbf{a} \mathbf{b} = \mathbf{a} \cdot \mathbf{b} + \mathbf{a} \wedge \mathbf{b}.\end{aligned} \hspace{\stretch{1}}(3.9)$

Pseudoscalar

In many contexts it is useful to introduce an ordered product of all the unit vectors for the space is called the pseudoscalar. In our 2D case this is

$\begin{aligned}i = \mathbf{e}_1 \mathbf{e}_2,\end{aligned} \hspace{\stretch{1}}(4.10)$

a quantity that we find behaves like the complex imaginary. That can be shown by considering its square

$\begin{aligned}(\mathbf{e}_1 \mathbf{e}_2)^2&=(\mathbf{e}_1 \mathbf{e}_2)(\mathbf{e}_1 \mathbf{e}_2) \\ &=\mathbf{e}_1 (\mathbf{e}_2 \mathbf{e}_1) \mathbf{e}_2 \\ &=-\mathbf{e}_1 (\mathbf{e}_1 \mathbf{e}_2) \mathbf{e}_2 \\ &=-(\mathbf{e}_1 \mathbf{e}_1) (\mathbf{e}_2 \mathbf{e}_2) \\ &=-1^2 \\ &= -1\end{aligned}$

Here the anticommutation of normal vectors property has been used, as well as (for the first time) the associative multiplication axiom.

In a 3D context, you’ll see the pseudoscalar in many places (expressing the normals to planes for example). It also shows up in a number of fundamental relationships. For example, if one writes

$\begin{aligned}I = \mathbf{e}_1 \mathbf{e}_2 \mathbf{e}_3\end{aligned} \hspace{\stretch{1}}(4.11)$

for the 3D pseudoscalar, then it’s also possible to show

$\begin{aligned}\mathbf{a} \mathbf{b} = \mathbf{a} \cdot \mathbf{b} + I (\mathbf{a} \times \mathbf{b})\end{aligned} \hspace{\stretch{1}}(4.12)$

something that will be familiar to the student of QM, where we see this in the context of Pauli matrices. The Pauli matrices also encode a Clifford algebraic structure, but we do not need an explicit matrix representation to do so.

Rotations

Very much like complex numbers we can utilize exponentials to perform rotations. Rotating in a sense from $\mathbf{e}_1$ to $\mathbf{e}_2$ , can be expressed as

$\begin{aligned}\mathbf{a} e^{i \theta}&=(\mathbf{e}_1 x_1 + \mathbf{e}_2 x_2) (\cos\theta + \mathbf{e}_1 \mathbf{e}_2 \sin\theta) \\ &=\mathbf{e}_1 (x_1 \cos\theta - x_2 \sin\theta)+\mathbf{e}_2 (x_2 \cos\theta + x_1 \sin\theta)\end{aligned}$

More generally, even in N dimensional Euclidean spaces, if $\mathbf{a}$ is a vector in a plane, and $\hat{\mathbf{u}}$ and $\hat{\mathbf{v}}$ are perpendicular unit vectors in that plane, then the rotation through angle $\theta$ is given by

$\begin{aligned}\mathbf{a} \rightarrow \mathbf{a} e^{\hat{\mathbf{u}} \hat{\mathbf{v}} \theta}.\end{aligned} \hspace{\stretch{1}}(5.13)$

This is illustrated in figure (1).

Plane rotation.

Notice that we have expressed the rotation here without utilizing a normal direction for the plane. The sense of the rotation is encoded by the bivector $\hat{\mathbf{u}} \hat{\mathbf{v}}$ that describes the plane and the orientation of the rotation (or by duality the direction of the normal in a 3D space). By avoiding a requirement to encode the rotation using a normal to the plane we have an method of expressing the rotation that works not only in 3D spaces, but also in 2D and greater than 3D spaces, something that isn’t possible when we restrict ourselves to traditional vector algebra (where quantities like the cross product can’t be defined in a 2D or 4D space, despite the fact that things they may represent, like torque are planar phenomena that do not have any intrinsic requirement for a normal that falls out of the plane.).

When $\mathbf{a}$ does not lie in the plane spanned by the vectors $\hat{\mathbf{u}}$ and $\hat{\mathbf{v}}$ , as in figure (2), we must express the rotations differently. A rotation then takes the form

$\begin{aligned}\mathbf{a} \rightarrow e^{-\hat{\mathbf{u}} \hat{\mathbf{v}} \theta/2} \mathbf{a} e^{\hat{\mathbf{u}} \hat{\mathbf{v}} \theta/2}.\end{aligned} \hspace{\stretch{1}}(5.14)$

3D rotation.

In the 2D case, and when the vector lies in the plane this reduces to the one sided complex exponential operator used above. We see these types of paired half angle rotations in QM, and they are also used extensively in computer graphics under the guise of quaternions.

References

[1] L. Dorst, D. Fontijne, and S. Mann. Geometric Algebra for Computer Science. Morgan Kaufmann, San Francisco, 2007.

[2] C. Doran and A.N. Lasenby. Geometric algebra for physicists. Cambridge University Press New York, Cambridge, UK, 1st edition, 2003.

[3] D. Hestenes. New Foundations for Classical Mechanics. Kluwer Academic Publishers, 1999.

Posted in Math and Physics Learning. | Tagged: anticommute, clifford algebra, commutative law, commute, cross product, distributive law, dot product, euclidean metrics, geometric algebra, geometric product, orthonormal basis, pauli matrix, pseudoscalar, quaternion, rotation, vector product, vector products, wedge product | 2 Comments »

« Previous Entries

	NVMCached: An NVM-ba… on Intel memory ordering, fence i…
	Adam C Scott on avoiding gdb signal noise…
	Ken on Scotiabank iTrade RESP ……
	Alan Ball on Oops. Fixing a drill hole in P…
	Peeter Joot's B… on Stokes theorem in Geometric…

Math, physics, perl, and programming obscurity.

Categories

Archives

Recent Posts

Meta

Recent Comments

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

Subscribe

Archive for March, 2012

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.

A mathematical description.

Problem Q3 (revisited).

Exact solutions.

Pictorially.

Relating the two pairs of coordinate systems.

Disclaimer.

Motivation.

Fluid flow over a solid body.

Aside

Pressure and vorticity equations with the non-linear term retained. (Reworked slightly)

Now what?

References

Motivation.

The equations of motion.

The boundary value constraints.

Laplacian of Pressure and Vorticity.

Separation of variables?

In terms of vorticity?

Pressure and vorticity equations with the non-linear term retained.

Now what?

Motivation.

Problem 1.

-waves, -waves, and Love-waves.

(b) constitutive relation, Newtonian fluids, and no-slip conditions.

Problem 2.

Statement.

References

Motivation.

Axioms

By example. The 2D case.

Pseudoscalar

Rotations

References

$\mathbf{P}$ -waves, $\mathbf{S}$ -waves, and Love-waves.