Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘vorticity vector’

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

Posted by Peeter Joot 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: , , , , | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 8: Phasor description of elastic waves. Fluid dynamics. Taught by Prof. K. Das.

Posted by Peeter Joot on February 6, 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.

Review. Elastic wave equation

Starting with

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

and applying a divergence operation we find

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

This is the P-wave equation. Applying a curl operation we find

\begin{aligned}\rho \frac{\partial^2 {{\boldsymbol{\omega}}}}{\partial {{t}}^2} &= C_T^2 \boldsymbol{\nabla}^2 \boldsymbol{\omega} \\ \boldsymbol{\omega} &= \boldsymbol{\nabla} \times \mathbf{e} \\ C_T^2 &= \frac{\lambda + 2\mu}{\rho}.\end{aligned} \hspace{\stretch{1}}(2.5)

This is the S-wave equation. We also found that

\begin{aligned}\frac{C_L}{C_T} > 1,\end{aligned} \hspace{\stretch{1}}(2.8)

and concluded that P waves are faster than S waves. What we haven’t shown is that the P waves are longitudinal, and that the S waves are transverse.

Assuming a gradient and curl description of our displacement

\begin{aligned}\mathbf{e} = \boldsymbol{\nabla} \phi + \boldsymbol{\nabla} \times \mathbf{H} = \mathbf{P} + \mathbf{S},\end{aligned} \hspace{\stretch{1}}(2.9)

we found

\begin{aligned}(\lambda + 2 \mu) \boldsymbol{\nabla}^2 \phi - \rho \frac{\partial^2 {{\phi}}}{\partial {{t}}^2} &= 0 \\ \mu \boldsymbol{\nabla}^2 \mathbf{H} - \rho \frac{\partial^2 {\mathbf{H}}}{\partial {{t}}^2} &= 0,\end{aligned} \hspace{\stretch{1}}(2.10)

allowing us to separately solve for the P and the S wave solutions respectively. Now, let’s introduce a phasor representation (again following section 22 of the text [1])

\begin{aligned}\phi &= A \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right) \\ \mathbf{H} &= \mathbf{B} \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right)\end{aligned} \hspace{\stretch{1}}(2.12)

Operating with the gradient we find

\begin{aligned}\mathbf{P}&= \boldsymbol{\nabla} \phi \\ &= \mathbf{e}_k \partial_k A \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right) \\ &= \mathbf{e}_k \partial_k A \exp\left( i ( k_m x_m - \omega t) \right) \\ &= \mathbf{e}_k i k_k A \exp\left( i ( k_m x_m - \omega t) \right) \\ &= i \mathbf{k} A \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right) \\ &= i \mathbf{k} \phi\end{aligned}

We can also write

\begin{aligned}\mathbf{P} = \mathbf{k} \phi'\end{aligned} \hspace{\stretch{1}}(2.14)

where \phi' is the derivative of \phi “with respect to its argument”. Here argument must mean the entire phase \mathbf{k} \cdot \mathbf{x} - \omega t.

\begin{aligned}\phi' = \frac{ d\phi( \mathbf{k} \cdot \mathbf{x} - \omega t )}{ d(\mathbf{k} \cdot \mathbf{x} - \omega t) } = i \phi\end{aligned} \hspace{\stretch{1}}(2.15)

Actually, argument is a good label here, since we can use the word in the complex number sense.

For the curl term we find

\begin{aligned}\mathbf{S}&= \boldsymbol{\nabla} \times \mathbf{H} \\ &= \mathbf{e}_a \partial_b H_c \epsilon_{a b c} \\ &= \mathbf{e}_a \partial_b \epsilon_{a b c} B_c \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right) \\ &= \mathbf{e}_a \partial_b \epsilon_{a b c} B_c \exp\left( i ( k_m x_m - \omega t) \right) \\ &= \mathbf{e}_a i k_b \epsilon_{a b c} B_c \exp\left( i ( \mathbf{k} \cdot \mathbf{x} - \omega t) \right) \\ &= i \mathbf{k} \times \mathbf{H}\end{aligned}

Again writing

\begin{aligned}\mathbf{H}' = \frac{ d\mathbf{H}( \mathbf{k} \cdot \mathbf{x} - \omega t )}{ d(\mathbf{k} \cdot \mathbf{x} - \omega t) } = i \mathbf{H}\end{aligned} \hspace{\stretch{1}}(2.16)

we can write the S wave as

\begin{aligned}\mathbf{S} = \mathbf{k} \times \mathbf{H}'\end{aligned} \hspace{\stretch{1}}(2.17)

Some waves illustrated.

The following wave types were noted, but not defined:

\begin{itemize}
\item Rayleigh wave. This is discussed in section 24 of the text (a wave that propagates near the surface of a body without penetrating into it). Wikipedia has an illustration of one possible mode of propagation [2].
\item Love wave. These aren’t discussed in the text, but wikipedia [3] describes them as polarized shear waves (where the figure indicates that the shear displacements are perpendicular to the direction of propagation).
\end{itemize}

Some illustrations from the class notes were also shown. Hopefully we’ll have some homework assignments where we do some problems to get a feel for how to apply the formalism.

Fluid dynamics.

In fluid dynamics we look at displacements with respect to time as illustrated in figure (\ref{fig:continuumL8:continuumL8fig1})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL8fig1}
\caption{Differential displacement.}
\end{figure}

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

In index notation

\begin{aligned}dx_i'&= dx_i + du_i \delta t \\ &= dx_i + \frac{\partial {u_i}}{\partial {x_j}} dx_j \delta t\end{aligned}

We define

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

a symmetric tensor. We also define

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

Effect of e_{ij} when diagonalized

\begin{aligned}e_{ij} ==\begin{bmatrix}e_{11} & 0 & 0 \\ 0 & e_{22} & 0 \\ 0 & 0 & e_{33}\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.21)

so that in this frame of reference we have

\begin{aligned}dx_1' &= ( 1 + e_{11} \delta t) dx_1 \\ dx_2' &= ( 1 + e_{22} \delta t) dx_2 \\ dx_3' &= ( 1 + e_{33} \delta t) dx_3\end{aligned} \hspace{\stretch{1}}(3.22)

Let’s find the matrix form of the antisymmetric tensor. We find

\begin{aligned}\omega_{11} = \omega_{22} = \omega_{33} = 0\end{aligned} \hspace{\stretch{1}}(3.25)

Introducing a vorticity vector

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

we find

\begin{aligned}\omega_{12} &= \frac{1}{{2}}\left( \frac{\partial {u_1}}{\partial {x_2}} -\frac{\partial {u_2}}{\partial {x_1}} \right) = - \frac{1}{{2}} (\boldsymbol{\nabla} \times \mathbf{u})_3 \\ \omega_{23} &= \frac{1}{{2}}\left( \frac{\partial {u_2}}{\partial {x_3}} -\frac{\partial {u_3}}{\partial {x_2}} \right) = - \frac{1}{{2}} (\boldsymbol{\nabla} \times \mathbf{u})_1 \\ \omega_{31} &= \frac{1}{{2}}\left( \frac{\partial {u_3}}{\partial {x_1}} -\frac{\partial {u_1}}{\partial {x_3}} \right) = - \frac{1}{{2}} (\boldsymbol{\nabla} \times \mathbf{u})_2\end{aligned} \hspace{\stretch{1}}(3.27)

Writing

\begin{aligned}\Omega_i = \frac{1}{{2}} \omega_i\end{aligned} \hspace{\stretch{1}}(3.30)

we find the matrix form of this antisymmetric tensor

\begin{aligned}\omega_{ij}=\begin{bmatrix}0 & -\Omega_3 & \Omega_2 \\ \Omega_3 & 0 & -\Omega_1 \\ -\Omega_2 & \Omega_1 & 0 \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.31)

\begin{aligned}dx_1'&= dx_1 + \left( \not{{\omega_{11}}} dx_1 + \omega_{12} dx_2 + \omega_{13} dx_3 \right) \delta t \\ &= dx_1 + \left( \omega_{12} dx_2 + \omega_{13} dx_3 \right) \delta t \\ &= dx_1 + \left( \Omega_2 dx_3 - \Omega_3 dx_2 \right) \delta t\end{aligned}

Doing this for all components we find

\begin{aligned}d\mathbf{x}' = d\mathbf{x} + (\boldsymbol{\Omega} \times d\mathbf{x}) \delta t.\end{aligned} \hspace{\stretch{1}}(3.32)

The tensor \omega_{ij} implies rotation of a control volume with an angular velocity \boldsymbol{\Omega} = \boldsymbol{\omega}/2 (half the vorticity vector).

In general we have

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

Making sense of things.

After this first fluid dynamics lecture I was left troubled. We’d just been barraged with a set of equations pulled out of a magic hat, with no notion of where they came from. Unlike the contiuum strain tensor, which was derived by considering differences in squared displacements, we have an antisymmetric term now. Why did we have no such term considering solids?

After a bit of thought I think I see where things are coming from. We have essentially looked at a first order decomposition of the displacement (per unit time) of a point in terms of symmetric and antisymmetric terms. This is really just a gradient evaluation, split into coordinates

\begin{aligned}x_i' &= x_i + (\boldsymbol{\nabla} u_i) \cdot d\mathbf{x} \delta t \\ &= x_i + \frac{\partial {u_i}}{\partial {x_j}} dx_j \delta t \\ &= x_i + \frac{1}{{2}}\left(\frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_i}}{\partial {x_j}} \right)dx_j \delta t +\frac{1}{{2}}\left(\frac{\partial {u_i}}{\partial {x_j}} -\frac{\partial {u_i}}{\partial {x_j}} \right)dx_j \delta t  \\ &=x_i + e_{ij} dx_j \delta t + \omega_{ij} dx_j \delta t\end{aligned}

Here, as in the solids case, we have

\begin{aligned}\mathbf{u} = \mathbf{x}' - \mathbf{x}\end{aligned} \hspace{\stretch{1}}(3.34)

References

[1] L.D. Landau, EM Lifshitz, JB Sykes, WH Reid, and E.H. Dill. Theory of elasticity: Vol. 7 of course of theoretical physics. 1960.

[2] Wikipedia. Rayleigh wave — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-February-2012]. http://en.wikipedia.org/w/index.php?title=Rayleigh_wave&oldid=473693354.

[3] Wikipedia. Love wave — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 4-February-2012]. http://en.wikipedia.org/w/index.php?title=Love_wave&oldid=474355253.

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

 
Follow

Get every new post delivered to your Inbox.

Join 36 other followers