Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘traction vector’

Continuum mechanics fluids review.

Posted by peeterjoot on April 25, 2012

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


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

Vector displacements.

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

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

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


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

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


Allowing us to write

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

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


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

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




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

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


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

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

Relative change in volume

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

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

Conservation of mass

Utilizing Green’s theorem we argued that

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

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

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

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

Constitutive relation

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

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

Conservation of momentum (Navier-Stokes)

As in elasticity, our momentum conservation equation had the form

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

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

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

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

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


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

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


No slip condition

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

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

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


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

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


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

Traction vector matching at an interface.

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

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

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

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

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

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

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

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


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

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

Worked problems from class.

A number of problems were tackled in class

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

Channel and shear flow

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

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

We find that this reduces to

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

with solution

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

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

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

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

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

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

we find

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

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

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

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

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

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

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

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

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

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


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

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

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

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

Mass conservation through apertures.

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

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

and therefore for incompressible fluids

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

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

Curve for tap discharge.

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

Figure 1: Tap flow measurement.

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

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

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

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

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

Mass conservation gives us

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


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

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

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

Figure 2: Comparison of measured stream radii and calculated.

Bernoulli equation.

With the body force specified in gradient for

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

and utilizing the vector identity

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

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

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


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

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

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

Surface tension.

Surfaces, normals and tangents.

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

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

Computing the gradient we find

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

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


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

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


Laplace pressure.

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

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

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

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

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

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

Surface tension gradients.

Considering the tangential component of the traction vector difference we find

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

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

Surface tension for a spherical bubble.

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

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

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

Figure 3: Spherical bubble in liquid.

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

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

A sample problem. The meniscus curve.

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

Figure 4: Curvature of fluid against a wall.

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

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

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

For fluid at rest, Navier-Stokes takes the form

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

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

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


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

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

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

We have then

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

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

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

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

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

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

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

The text introduces the capillary constant

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

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

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

we can integrate to find

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

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

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

Integrating this with Mathematica I get

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

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

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

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

Non-dimensionality and scaling.

With the variable transformations

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

we can put Navier-Stokes in dimensionless form

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

Here R is Reynold’s number

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

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

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

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

Eulerian and Lagrangian.

We defined

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

Boundary layers.

Impulsive flow.

We looked at the time dependent unidirectional flow where

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

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

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

and were able to show that

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

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

FIXME: really need to plot 16.80.

Oscillatory flow.

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

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

with an assumed solution of the form

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

we found


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

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


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

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

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


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

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

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


With boundary conditions

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

With a similarity variable

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

and stream functions

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


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

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

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

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

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

Singular perturbation theory.

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

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

where the inverse of Reynold’s number

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

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

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

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

for which the exact solution was found to be

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

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

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

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

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

This gives us

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

for which we find

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


We characterized stability in terms of displacements writing

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

and defining

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

Thermal stability: Rayleigh-Benard problem.

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

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

and introducing an equation for the base state

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

we found

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

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

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

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

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

With an assumption that density change and temperature are linearly related

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

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

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

We also applied our perturbation to the energy balance equation

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

We determined that the base state temperature obeyed

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

with solution

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

This and application of the perturbation gave us

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

We used this to non-dimensionalize with

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

And found (primes dropped)


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

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


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

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

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


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

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

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

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

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

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

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

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

which is the boundary for thermal stability or instability.

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


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

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

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

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

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

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

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

Posted by peeterjoot on April 8, 2012

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


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

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

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

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

In either case we have

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

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

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

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

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

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

Verifying the relationship.

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

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

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

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

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

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

so we find

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

Cylindrical strain tensor.

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

Our gradient in cylindrical coordinates is the familiar

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

and our cylindrical velocity is

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

Our curl is then

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

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

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

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

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

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

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

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

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

Putting things together we have

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

For our stress tensor

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

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


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

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

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


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

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

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

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

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

Putting things together we have

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

For our stress tensor

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

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


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

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

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


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

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

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

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

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

Putting things together we have

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

For our stress tensor

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

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


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

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

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




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

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

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

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

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

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


Spherical strain tensor.

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

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

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

and our gradient mirrors this structure

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

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


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

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

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


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

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

Finally, our velocity in spherical coordinates is just

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

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

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

So we have

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

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

The directional derivative portion of our strain is

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

The other portion of our strain tensor is

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

Putting these together we find

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

Which gives

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

For our stress tensor

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

we can now read off our components by taking dot products


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

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

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


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

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

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

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

So we have

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

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

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

Putting these together we find

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

Which gives

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

For our stress tensor

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

we can now read off our components by taking dot products


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

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

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


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

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

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

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

So we have

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

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

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

Putting these together we find

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

Which gives

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

For our stress tensor

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

we can now read off our components by taking dot products


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

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

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


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



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

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

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

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

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

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



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

[2] Peeter Joot. Continuum mechanics., chapter {Introduction and strain tensor.}

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

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


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

\caption{Gravitational force components acting on fluid flowing down a plane.}


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.

\caption{Diagram of coordinates for inclined flow problem.}


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


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


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


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

\caption{Differential vector element.}


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


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


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


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

\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.

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


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.

\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.

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}


\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

\caption{A 1D compression wave.}

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

\mathbf{P}-waves are longitudinal.

\mathbf{S}-waves are transverse.

\item Whose speed is higher?

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?


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

\caption{Love wave illustrated.}

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

\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?

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

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.

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.


Problem 2.


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.

\caption{Shearing flow with pressure gradient and one moving boundary.}

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.

\item Derive the velocity profile of the fluid.


Our equations of motion are


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


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.


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})
\caption{Parabolic velocity profile.}

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


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})
\caption{Shear flow.}

\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.


For low U we’ll let the parabolic dominate, and can graphically add these two as in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig5})
\caption{Superposition of shear and parabolic flow (low U)}
For high U, we’ll let the shear flow dominate, and have plotted this in figure (\ref{fig:continuumMidtermReflection:continuumMidtermReflectionFig6})
\caption{Superposition of shear and parabolic flow (high U)}

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

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.

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.

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.

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)



[1] Wikipedia. Newtonian fluid — wikipedia, the free encyclopedia [online]. 2011. [Online; accessed 17-March-2012].

[2] Wikipedia. Non-newtonian fluid — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 17-March-2012].

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

PHY454H1S Continuum mechanics. Problem Set 2. Rectilinear flow with pressure gradient. Two layer flow driven by moving surface.

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


This problem set is as yet ungraded.

Problem Q1.


Imagine a steady rectilinear blood flow of the form \mathbf{u} = u(y) \hat{\mathbf{y}} through a two dimensional artery. It is driven by a constant pressure gradient G = -dp/dx maintained by an external `heart’. The top and bottom walls of the artery are 2h distance apart and the fluid satisfies no-slip boundary conditions at the walls. Assuming that the fluid is Newtonian,

\item Show that the Navier-Stokes equation reduces to

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

where \mu is the viscosity of the blood.

\item Show that the velocity profile of the fluid inside the artery is a parabolic profile.
\item What is the maximum speed of the fluid? Draw the velocity profile to show where the maximum speed occurs inside the artery.
\item If due to smoking etc., the viscosity of the blood gets doubled, then what should be the new pressure gradient to be maintained by the `heart’ to keep the liquid flux through the artery at the same level as the non-smoking one?

Solution. Part 1. Navier-Stokes equation for the system.

The Navier-Stokes equation, for an incompressible unidirectional fluid \mathbf{u} = (u, 0, 0), assuming that there is no z dependence, takes the form


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

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

\begin{aligned}0 = - \frac{\partial {p}}{\partial {z}} \end{aligned} \hspace{\stretch{1}}(2.2c)

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


With a steady state assumption we kill the {\partial {u}}/{\partial {t}} term, and 2.2d kills of the x-component of the Laplacian and our non-linear inertial term on the LHS, leaving just


\begin{aligned}0 = - \frac{\partial {p}}{\partial {x}} + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} \end{aligned} \hspace{\stretch{1}}(2.3a)

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

\begin{aligned}0 = - \frac{\partial {p}}{\partial {z}}.\end{aligned} \hspace{\stretch{1}}(2.3c)


With \frac{\partial {p}}{\partial {z}} = \frac{\partial {p}}{\partial {y}} = 0, we have \frac{\partial {p}}{\partial {x}} = dp/dx = -G, so 2.3a is reduced to

\begin{aligned}0 = G + \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2}.\end{aligned} \hspace{\stretch{1}}(2.4)

Finally, since we have an assumption of no z-dependence ({\partial {u}}/{\partial {z}} = 0) and from the incompressibility assumption 2.2d ({\partial {u}}/{\partial {x}} = 0), we have

\begin{aligned}\frac{\partial^2 {{u}}}{\partial {{y}}^2} = \frac{d^2 u}{dy^2} = -\frac{G}{\mu},\end{aligned} \hspace{\stretch{1}}(2.5)

as desired.

Solution. Part 2. Velocity profile.

Integrating 2.5 twice, we have

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

Application of the no-slip boundary value condition u(\pm h) = 0, we have

\begin{aligned}0 &= -\frac{G}{2 \mu} h^2 + A h + B \\ 0 &= -\frac{G}{2 \mu} h^2 - A h + B \end{aligned} \hspace{\stretch{1}}(2.7)

Adding and subtracting these, we find


\begin{aligned}A = 0 \end{aligned} \hspace{\stretch{1}}(2.9a)

\begin{aligned}B = \frac{G}{2 \mu} h^2,\end{aligned} \hspace{\stretch{1}}(2.9b)


so the velocity is given by the parabolic function

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

Solution. Part 3. Maximum speed of the flow.

It is clear that the maximum speed of the fluid is found at y = 0

\begin{aligned}u(0) = \frac{G h^2}{2 \mu}\end{aligned} \hspace{\stretch{1}}(2.11)

The velocity profile for this flow is drawn in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig1}).

\caption{Velocity profile for 1D constant pressure gradient steady state flow.}

Solution. Part 4. Effects of viscosity doubling.

With our velocity being dependent on the G/\mu ratio, it is clear, even without calculating the flux, that we will need twice the pressure gradient if the viscosity is doubled to maintain the same flux through the artery and veins. To demonstrate this more thoroughly we can calculate this mass flux. For an element of mass leaving a portion of the conduit, bounded by the plane normal to \hat{\mathbf{x}} we have

\begin{aligned}\frac{dm}{dt} &= \rho \frac{dV}{dt} \\ &= \rho dz dy \mathbf{u} \cdot \hat{\mathbf{x}}\end{aligned}

Integrating this over a width \Delta z, our flux through the plane is

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

Doubling the blood viscosity for our smoker, our respective fluxes are

\begin{aligned}\text{Flux}_{\text{smoker}} &= \Delta z \frac{2 G_{\text{smoker}} h^3}{3 (2 \mu)}  \\ \text{Flux}_{\text{non-smoker}} &= \Delta z \frac{2 G h^3}{3 \mu}.\end{aligned} \hspace{\stretch{1}}(2.12)

Demanding equality before and after smoking we find

\begin{aligned}G_{\text{smoker}} = 2 G.\end{aligned} \hspace{\stretch{1}}(2.14)

where G is the magnitude of the pressure gradient before the bad habits kicked in. The smoker’s poor little heart (soon to be a big overworked and weak heart) has to generate pressure gradients that are twice as big to get the same quantity of blood distributed through the body.

Problem Q2.


Consider steady simple shearing flow with no imposed pressure gradient (G = 0) of a two layer fluid with viscosity

\begin{aligned}\mu =\left\{\begin{array}{l l}\mu^{(1)} & \quad \mbox{latex -h < y < 0,$} \\ \mu^{(2)} & \quad \mbox{0 < y < h.}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.15)$

The boundary conditions are no-slip at the lower plate (y = -h) and at y = 0. The top plate is moving with a velocity -U at y = h and fluid is sticking to it. using the continuity of tangential (shear) stress at the interface (y = 0)

\item Derive the velocity profile of the two fluids.
\item Calculate the maximum speed.
\item Calculate the mean speed.
\item Calculate the flux (the volume flow rate.)
\item Calculate the tangential force (per unit width) F_x on the strip 0 \le x \le L of the wall y = -h.
\item Calculate the tangential force (per unit width) F_x^0 on the strip 0 \le x \le L at the interface y = 0 by the top fluid on the lower fluid.

Part 1. Derive the velocity profile of the two fluids.

With coordinates referring to figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig2}), our steady flow for layers 1 and 2 has the form

\caption{Two layer flow induced by moving wall.}


\begin{aligned}0 = - \frac{\partial {p}}{\partial {x}} + \mu^{(i)} \frac{\partial^2 {{u^{(i)}}}}{\partial {{y}}^2}\end{aligned} \hspace{\stretch{1}}(3.16a)

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

\begin{aligned}0 = - \frac{\partial {p}}{\partial {z}},\end{aligned} \hspace{\stretch{1}}(3.16c)


as we found in Q1. Only the boundary value conditions and the driving pressure are different here. In this problem and the next, we have constant pressure gradients dp/dx = -G to deal with, so we really have just the pair of equations

\begin{aligned}0 = G + \mu^{(i)} \frac{d^2 u^{(i)}(y)}{dy^2},\end{aligned} \hspace{\stretch{1}}(3.17)

to solve. For this Q2 problem we have G = 0, so the algebra to match our boundary value constraints becomes a bit easier. Our boundary value constraints are


\begin{aligned}u^{(1)}(-h) = 0 \\ \end{aligned} \hspace{\stretch{1}}(3.18a)

\begin{aligned}u^{(2)}(h) = -U \\ \end{aligned} \hspace{\stretch{1}}(3.18b)

\begin{aligned}u^{(1)}(0) = u^{(2)}(0),\end{aligned} \hspace{\stretch{1}}(3.18c)


plus one more to match the tangential components of the traction vector with respect to the normal \hat{\mathbf{n}} = (0, 1, 0). The components of that traction vector are

\begin{aligned}t_i &= \left( -p \delta_{ij} + 2 \mu e_{ij} \right) n_j \\ &= \left( -p \delta_{ij} + 2 \mu e_{ij} \right) \delta_{2j} \\ &= -p \delta_{i2} + 2 \mu e_{i2},\end{aligned}

but we are only interested in the horizontal component t_1 which is

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

So the matching the tangential components of the traction vector at the interface gives us our last boundary value constraint

\begin{aligned}{\left.{{\mu^{(1)} \frac{\partial {y}}{\partial {u}}}}\right\vert}_{{y = 0}} = {\left.{{\mu^{(2)} \frac{\partial {y}}{\partial {u}}}}\right\vert}_{{y = 0}},\end{aligned} \hspace{\stretch{1}}(3.19)

and we are ready to do our remaining bits of algebra. We wish to solve the pair of equations

\begin{aligned}u^{(1)} &= A^{(1)} y + B^{(1)} \\ u^{(2)} &= A^{(2)} y + B^{(2)},\end{aligned} \hspace{\stretch{1}}(3.20)

for the four integration constants A^{(i)} and B^{(i)} using our boundary value constraints. The linear system to solve is

\begin{aligned}0 &= -A^{(1)} h + B^{(1)} \\ -U &= A^{(2)} h + B^{(2)} \\ B^{(1)} &= B^{(2)} \\ \mu^{(1)} A^{(1)} &= \mu^{(2)} A^{(2)}.\end{aligned} \hspace{\stretch{1}}(3.22)

With B = B^{(i)}, we have

\begin{aligned}0 &= -A^{(1)} h + B \\ -U &= \frac{\mu^{(1)}}{\mu^{(2)}} A^{(1)} h + B \end{aligned} \hspace{\stretch{1}}(3.26)

Subtracting these to solve for A^{(1)} we find

\begin{aligned}-U = h A^{(1)} \left( \frac{\mu^{(1)}}{\mu^{(2)}} + 1 \right).\end{aligned} \hspace{\stretch{1}}(3.28)

This gives us everything we need

\begin{aligned}A^{(1)} &= -\frac{U \mu^{(2)}}{h(\mu^{(1)} + \mu^{(2)})} \\ A^{(2)} &= -\frac{U \mu^{(1)}}{h(\mu^{(1)} + \mu^{(2)})} \\ B^{(1)} = B^{(2)} &= -\frac{U \mu^{(2)}}{\mu^{(1)} + \mu^{(2)}},\end{aligned} \hspace{\stretch{1}}(3.29)

Referring back to 3.20 our velocities are

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

Checking, we see at a glance we see that we have u^{(2)}(h) = -U, u^{(1)}(-h) = 0, u^{(1)}(0) = u^{(2)}(0), and {\left.{{\mu^{(1)} du^{(1)}/dy}}\right\vert}_{{y=0}} = \mu^{(2)} {\left.{{du^{(2)}/dy}}\right\vert}_{{y=0}} as desired.

As an example, let’s add some numbers. With mercury and water in layers {(1)} and {(2)} respectively, we have

\begin{aligned}\mu^{(1)} &= 0.001526 \quad \text{Pa-s} \\ \mu^{(2)} &= 0.00089 \quad \text{Pa-s} \end{aligned} \hspace{\stretch{1}}(3.33)

so that our velocity is

\begin{aligned}u(y) = \left\{\begin{array}{l l}-1.37 U \left(1 + \frac{y}{h}\right) & \quad \mbox{latex y \in [-h, 0]$} \\ -1.37 U \left(1 + 1.7 \frac{y}{h}\right) & \quad \mbox{y \in [0, h]}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.35)$

This is plotted with h = U = 1 in figure (\ref{fig:continuumProblemSet2:continuumProblemSet2Fig3})
\caption{Two layer shearing flow with water over mercury.}

Part 2. Calculate the maximum speed.

With u^{(1)}(-h) = 0, and u^{(1)} linearly decreasing, then u^{(2)} linearly decreasing further from the value at y = 0, it is clear that the maximum speed, no matter the viscosities of the fluids, is on the upper moving interface. This maximum takes the value {\left\lvert{u^{(2)}(h)}\right\rvert} = U.

Part 3. Calculate the mean speed.

As linear functions the average speeds of the respective fluids fall on the midpoints y = \pm h/2. These are

\begin{aligned}\left\langle{{u^{(1)}}}\right\rangle &= -\frac{U \mu^{(2)} }{ 2 (\mu^{(1)} + \mu^{(2)})} \\ \left\langle{{u^{(2)}}}\right\rangle &= -\frac{U \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \left( 1 + \frac{\mu^{(1)}}{2 \mu^{(2)}} \right) \end{aligned} \hspace{\stretch{1}}(3.36)

Averaging these two gives us the overall average, so we find

\begin{aligned}\left\langle{{u(y)}}\right\rangle = -\frac{U }{ 4 (\mu^{(1)} + \mu^{(2)})} \left( 3 \mu^{(2)} + \mu^{(1)} \right) \end{aligned} \hspace{\stretch{1}}(3.38)

Part 4. Calculate the flux

We can calculate the volume flux, much like the mass flux (although the mass flux seems a more sensible quantity to calculate). Looking at the rate of change of an element of fluid passing through the y-z plane we have

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

Integrating over the total height, for a width \Delta z we have

\begin{aligned}\text{Volume Flux} &= \Delta z \int_{-h}^h u(y) dy \\ &= \Delta z 2h \left\langle{{u}}\right\rangle \\ \end{aligned}

So our volume flux through a width \Delta z is

\begin{aligned}\text{Volume Flux} =-\frac{2 h U \Delta z}{ 4 (\mu^{(1)} + \mu^{(2)})} \left( 3 \mu^{(2)} + \mu^{(1)} \right).\end{aligned} \hspace{\stretch{1}}(3.40)

Part 5. Tangential force on the wall.

From 3.32 the tangential components of our traction vectors are

\begin{aligned}t^{(1)} &= \mu^{(1)} \frac{d u^{(1)}}{dy} \\ &= -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{1}{{h}}\end{aligned}


\begin{aligned}t^{(2)} &= \mu^{(2)} \frac{d u^{(2)}}{dy} \\ &= -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{1}{{h}}\end{aligned}

We see that the tangential component of the traction vector is a constant throughout both fluids. Allowing this force to act on a length L of the lower wall, our force per unit width over that strip is just

\begin{aligned}F = -\frac{U \mu^{(1)} \mu^{(2)} }{(\mu^{(1)} + \mu^{(2)})} \frac{L}{h}.\end{aligned} \hspace{\stretch{1}}(3.41)

The negative value here makes sense since it is acting to push the fluid backwards in the direction of the upper wall motion.

Part 6. Tangential force on the lower fluid.

Due to constant nature of the tangential component of the traction vector shown above, the force per unit width of the upper fluid acting on the lower fluid, is also given by 3.41.

Problem Q3.


If on top of the problem described above a constant pressure gradient G = -dp/dx is applied between the boundaries y = \pm h, describe qualitatively what type of flow profile you would expect in the steady state. Draw the velocity profiles for two cases (i) \mu^{(1)} > \mu^{(2)} (ii) \mu^{(1)} < \mu^{(2)}. Explain your result.


We showed earlier that the Navier-Stokes equations for this Q3 case, where G is non-zero were given by 3.17, which restated is

\begin{aligned}0 = G + \mu^{(i)} \frac{d^2 u^{(i)}(y)}{dy^2}.\end{aligned} \hspace{\stretch{1}}(4.42)

Our solutions will now necessarily be parabolic, of the form

\begin{aligned}u^{(i)}(y) = -\frac{G}{2 \mu^{(i)}} y^2 + A^{(i)} y + B^{(i)},\end{aligned} \hspace{\stretch{1}}(4.43)

with the tangential traction vector components given by

\begin{aligned}t^{(i)} = -G y + A^{(i)} \mu^{(i)}\end{aligned} \hspace{\stretch{1}}(4.44)

The boundary value constants become a bit messier to solve for, and should we wish to do so we’d have to solve the system

\begin{aligned}0 &= -\frac{G}{2 \mu^{(1)}} h^2 - A^{(1)} h + B^{(1)} \\ -U &= -\frac{G}{2 \mu^{(2)}} h^2 + A^{(2)} h + B^{(2)} \\ B^{(1)} &= B^{(2)} \\ A^{(1)} \mu^{(1)} &= A^{(2)} \mu^{(2)} \end{aligned} \hspace{\stretch{1}}(4.45)

Without actually solving this system we should expect that our solution will have the form of our pure shear flow, with parabolas superimposed on these linear flows. 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}).

\caption{Superposition of constant pressure gradient and shear flow solutions (\mu^{(1)} > \mu^{(2)}).}

\caption{Superposition of constant pressure gradient and shear flow solutions (\mu^{(1)} < \mu^{(2)}).}

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, we get

\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}}(4.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.

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

PHY454H1S Continuum Mechanics. Lecture 15: More on surface tension and Reynold’s number. Taught by Prof. K. Das.

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


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

Review. Surface tension clarifications

For a surface like figure (\ref{fig:continuumL15:continuumL15Fig1})
\caption{Vapor liquid interface.}

we have a discontinuous jump in density. We will have to consider three boundary value constraints

\item Mass balance. This is the continuity equation.
\item Momentum balance. This is the Navier-Stokes equation.
\item Energy balance. This is the heat equation.

We have not yet discussed the heat equation, but this is required for non-isothermal problems.

The boundary condition at the interface is given by the stress balance. Denoting the difference in the traction vector at the interface by

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

Here the gradient is in the tangential direction of the surface as in figure (\ref{fig:continuumL15:continuumL15Fig2})
\caption{normal and tangent vectors on a curve.}

In the normal direction

\begin{aligned}[\mathbf{t}]_1^2 \cdot \hat{\mathbf{n}}&= (\mathbf{t}_2 - \mathbf{t}_1) \cdot \hat{\mathbf{n}} \\ &= -\frac{\sigma}{2 R} \end{aligned}

With the traction vector having the value

\begin{aligned}\mathbf{t} &= \mathbf{e}_i T_{ij} n_j \\ &= \mathbf{e}_i \left( -p \delta_{ij} + \mu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right)\right)n_j\end{aligned}

We have in the normal direction

\begin{aligned}\mathbf{t} \cdot \mathbf{n} =n_i \left( -p \delta_{ij} + \mu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right)\right) n_j\end{aligned} \hspace{\stretch{1}}(2.2)

With \mathbf{u} = 0 on the surface, and n_i \delta_{ij} n_j = n_j n_j = 1 we have

\begin{aligned}\mathbf{t} \cdot \mathbf{n} = -p\end{aligned} \hspace{\stretch{1}}(2.3)

Returning to (\mathbf{t}_2 - \mathbf{t}_1) \cdot \hat{\mathbf{n}} we have

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

This is the Laplace pressure. Note that the sign of the difference is significant, since it effects the direction of the curvature. This is depicted pictorially in figure (\ref{fig:continuumL15:continuumL15Fig3})
\caption{pressure and curvature relationships}

Question In [1] the curvature term is written

\begin{aligned}\frac{1}{{2 R}} \rightarrow \frac{1}{{R_1}} + \frac{1}{{R_2}}.\end{aligned} \hspace{\stretch{1}}(2.5)

Why the difference? Answer: This is to account for non-spherical surfaces, with curvature in two directions. Illustrating by example, imagine a surface like as in figure (\ref{fig:continuumL15:continuumL15Figq})
\caption{Example of non-spherical curvature.}

Followup required to truly understand things

While, this review clarifies things, we still don’t really know how the surface tension term \sigma is defined. Nor have we been given any sort of derivation of 2.1, from which the end result follows.

I’m assuming that \sigma is a property of the two fluids at the interface, so if you have, for example, oil and vinegar in a bottle, we have surface tension and curvature that’s probably related to how the two of these interact. If there is still a mixing or settling process occurring, I’d imagine that this could even vary from point to point on the surface (imagine adding soap to a surface where stuff can float until the soap mixes in enough that things start sinking in the radius of influence of the soap).

Surface tension gradients.

Now consider the tangential component of the traction vector

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

So we see that for a static fluid, we must have

\begin{aligned}\boldsymbol{\nabla}_I \sigma = 0\end{aligned} \hspace{\stretch{1}}(2.7)

For a static interface there cannot be any surface tension gradient. This becomes very important when considering stability issues. We can have surface tension induced flow called capillary, or mandarin (?) flow.

Reynold’s number.

In Navier-Stokes after making non-dimensionalization changes of the form

\begin{aligned}x \rightarrow L x'\end{aligned} \hspace{\stretch{1}}(3.8)

the control parameter is like Reynold’s number.


\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}\end{aligned} \hspace{\stretch{1}}(3.9)

We call the term

\begin{aligned}\rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u}\end{aligned} \hspace{\stretch{1}}(3.10)

the inertial term. It is non-zero only when something is being “carried along with the velocity”. Consider a volume fixed in space and one that is moving along with the fluid as in figure (\ref{fig:continuumL15:continuumL15Fig4})
\caption{Moving and fixed frame control volumes in a fluid.}

All of our viscosity dependence shows up in the Laplacian term, so we can roughly characterize the Reynold’s number as the ratio

\begin{aligned}\text{Reynold's number}&\rightarrow\frac{{\left\lvert{\text{effect of inertia}}\right\rvert}}{{\left\lvert{\text{effect of viscosity}}\right\rvert}}  \\ &=\frac{ {\left\lvert{\rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} }\right\rvert} }{{\left\lvert{ \mu \boldsymbol{\nabla}^2 \mathbf{u}}\right\rvert}} \\ &\sim\frac{ \rho U^2/L }{\mu U/L^2} \\ &\sim\frac{ \rho U L }{\mu }\end{aligned}

In figures (\ref{fig:continuumL15:continuumL15Fig5}), (\ref{fig:continuumL15:continuumL15Fig6}) we have two illustrations of viscous and non-viscous regions the first with a moving probe pushing its way through a surface, and the second with a wing set at an angle of attack that generates some turbulence. Both are illustrations of the viscous and inviscous regions for the two flows. Both of these are characterized by the Reynold’s number in some way not really specified in class. One of the points of mentioning this is that when we are in an essentially inviscous region, we can neglect the viscosity (\mu \boldsymbol{\nabla}^2 \mathbf{u}) term of the flow.


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

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

PHY454H1S Continuum Mechanics. Lecture 14: Non-dimensionality and scaling. Taught by Prof. K. Das.

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


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

Review. Surfaces

We are considering a surface as depicted in (\ref{fig:continuumL14:continuumL14Fig13})

\caption{Variable surface geometries}

With the surface height given by

\begin{aligned}z = h(x, t),\end{aligned} \hspace{\stretch{1}}(2.1)

where this describes the interface. Taking the difference

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

we define a surface. We considered a small displacement as in (\ref{fig:continuumL14:continuumL14Fig14}).

\caption{A vector differential element}

Recall that if \phi is a constant, then \boldsymbol{\nabla} \phi is a normal to the surface. We showed this by considering the differential

\begin{aligned}0 &= d\phi \\ &= \frac{\partial {\phi}}{\partial {x}} dx+\frac{\partial {\phi}}{\partial {y}} dy+\frac{\partial {\phi}}{\partial {z}} dz \\ &=(\boldsymbol{\nabla} \phi) \cdot d\mathbf{r}.\end{aligned}

We can construct the unit normal by scaling. For our 1D example we have

\begin{aligned}\hat{\mathbf{n}} &= \frac{\boldsymbol{\nabla} \phi}{{\left\lvert{\boldsymbol{\nabla} \phi}\right\rvert}} \\ &= \frac{1}{{{\left\lvert{\boldsymbol{\nabla} \phi}\right\rvert}}} \left(\frac{\partial {\phi}}{\partial {x}},\frac{\partial {\phi}}{\partial {y}}\right) \end{aligned}

so that our unit normal is

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

A unit tangent can also be constructed by inspection

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

Traction vector at the interface.

Recall that our stress tensor has the form

\begin{aligned}T_{ij} = - p \delta_{ij} + \rho \nu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right)\end{aligned} \hspace{\stretch{1}}(3.5)

(here we are switching notations for the stress since we will be using \sigma for surface tension in this section)

The traction vector components are

\begin{aligned}t_i = T_{ij} n_j =- p n_i + \rho \nu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right) n_j\end{aligned} \hspace{\stretch{1}}(3.6)

Considering a control volume as illustrated in we can arrive at what we call the jump stress balance equation

figure (\ref{fig:continuumL14:continuumL14fig3})
\caption{Control volume for liquid air interface}

\begin{aligned}[\mathbf{T} \hat{\mathbf{n}}]^2_1 = \frac{2 \sigma}{R} \hat{\mathbf{n}} - \boldsymbol{\nabla}_I \sigma\end{aligned} \hspace{\stretch{1}}(3.7)


\begin{aligned}\sigma &= \text{surface tension} \\ R &= \text{radius of curvature} \\ \boldsymbol{\nabla}_I &= \text{gradient along the interface}\end{aligned} \hspace{\stretch{1}}(3.8)

and the suffix 2 and prefix 1 indicates that we are considering the interface between fluids labelled 1 and 2 (liquid and air respectively in the diagram).

For a derivation see Prof after class?

Force balance along the normal direction gives

\begin{aligned}\hat{\mathbf{n}} [\mathbf{T} \hat{\mathbf{n}}]^2_1 = \hat{\mathbf{n}} \cdot \frac{2 \sigma}{R} \hat{\mathbf{n}} - \not{{\hat{\mathbf{n}} \cdot (\boldsymbol{\nabla}_I \sigma)}}\end{aligned} \hspace{\stretch{1}}(3.11)

If you do this calculation, you will get

\begin{aligned}[-p]^2_1 = \frac{ 2 \sigma}{R}\end{aligned} \hspace{\stretch{1}}(3.12)

I think this was called the Laplace equation?

Question: How was \sigma defined? A: Energy per unit area.

Figure (\ref{fig:continuumL14:continuumL14fig4}) was given as part of an explaination of surface tension and curvature, but I missed part of that discussion. Perhaps this is elaborated on in the class notes?

\caption{Molecular gas and liquid interactions at a surface.}

Non dimensionalization and scaling


By scaling we mean how much detail do you want to look at in the analysis. Consider the figure (\ref{fig:continuumL14:continuumL14fig5a}) where we imagine that we zoom in on something that appears smooth from a distance, but perhaps grandular close up as in figure (\ref{fig:continuumL14:continuumL14fig5b}). Picking the length scale to be used in this case can be very important.

\caption{Coarse scaling example.}

\caption{Fine grain scaling example (a zoom).}

FIXME: prof wrote:

\begin{aligned}y = A x^\alpha\end{aligned} \hspace{\stretch{1}}(4.13)

what was that about?

Rescaling by characteristic length and velocity.

Suppose that a fluid is flowing with

\item a characteristic velocity U, with dimensions [U] \sim L T^{-1}
\item a characteristic length scale L

Considering the dimensions of the terms in the NS equation

\begin{aligned}[\rho] = M L^{-3}\end{aligned} \hspace{\stretch{1}}(4.14)

\begin{aligned}[p] = M L T^{-2} L^{-2} = M L^{-1} T^{-2}\end{aligned} \hspace{\stretch{1}}(4.15)

\begin{aligned}[t] = T = \frac{L}{U}\end{aligned} \hspace{\stretch{1}}(4.16)


\begin{aligned}[p] = [\rho U^2] = M L^{-3} L^2 T^{-2}  = M L^{-1} T^{-2}\end{aligned} \hspace{\stretch{1}}(4.17)

Now let’s alter the NS equation using some scaling to put it into a dimensionless form

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

\begin{aligned}\frac{\partial {\mathbf{u}}}{\partial {t}} \rightarrow  \frac{\partial {(U \mathbf{u}')}}{\partial {\left(\frac{L}{U} t'\right)}} = \frac{U^2}{L} \frac{\partial {\mathbf{u}'}}{\partial {t'}}\end{aligned} \hspace{\stretch{1}}(4.19)

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{x}} \frac{\partial {}}{\partial {x}}+\hat{\mathbf{y}} \frac{\partial {}}{\partial {y}}+\hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\rightarrow \hat{\mathbf{x}} \frac{\partial {}}{\partial {L x'}}\hat{\mathbf{y}} \frac{\partial {}}{\partial {L y'}}\hat{\mathbf{z}} \frac{\partial {}}{\partial {L z'}}\end{aligned} \hspace{\stretch{1}}(4.20)

so that

\begin{aligned}\boldsymbol{\nabla} \rightarrow \frac{1}{{L}} \boldsymbol{\nabla}'\end{aligned} \hspace{\stretch{1}}(4.21)

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} \rightarrow \left( U \mathbf{u}' \cdot \frac{1}{{L}} \boldsymbol{\nabla}' \right) U \mathbf{u}' = \frac{U^2}{L} (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}'\end{aligned} \hspace{\stretch{1}}(4.22)

\begin{aligned}\frac{1}{{\rho}} \boldsymbol{\nabla} p \rightarrow \frac{1}{{L}} \frac{\boldsymbol{\nabla}' (\not{{\rho}} U^2) }{\not{{\rho}}} p' = \frac{U^2}{L} \boldsymbol{\nabla}' p'\end{aligned} \hspace{\stretch{1}}(4.23)

\begin{aligned}\nu \boldsymbol{\nabla}^2 \mathbf{u} \rightarrow \frac{\nu}{L^2} \boldsymbol{\nabla}' U \mathbf{u}' = \frac{\nu U}{L^2} \boldsymbol{\nabla}' \mathbf{u}'\end{aligned} \hspace{\stretch{1}}(4.24)

Putting everything together, NS takes the form

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


\begin{aligned}\frac{\partial {\mathbf{u}'}}{\partial {t'}} + (\mathbf{u}' \cdot \boldsymbol{\nabla}') \mathbf{u}' = \boldsymbol{\nabla}' p' + \frac{\nu}{U L} \boldsymbol{\nabla}' \mathbf{u}'\end{aligned} \hspace{\stretch{1}}(4.26)

Introducing the Reynold’s number

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

We have NS in dimensionless form

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

The implications of this will be discussed further in the next lecture.

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

PHY454H1S Continuum Mechanics. Lecture 11: Worked examples of Navier-Stokes solutions. Taught by Prof. K. Das.

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


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

Navier-Stokes equation.

The Navier-Stokes equation (our fluids equivalent to Newton’s second law) was found to be

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

In this course we’ll focus on the incompressible case where we have

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

We watched a video of the rocking tank as in figure (\ref{fig:continuumL11:continuumL11fig1}). The boundary condition that accounted for the matching of the die marker is that we have no slipping at the interface. Writing \boldsymbol{\tau} for the tangent to the interface then this condition at the interface is described mathematically by the conditions

\caption{Rocking tank velocity matching.}

\begin{aligned}\mathbf{u}_A \cdot \boldsymbol{\tau} &= \mathbf{u}_B \cdot \boldsymbol{\tau} \\ \mathbf{u}_A \cdot \hat{\mathbf{n}} &= \mathbf{u}_B \cdot \hat{\mathbf{n}}.\end{aligned} \hspace{\stretch{1}}(2.3)

Referring to figure (\ref{fig:continuumL11:continuumL11fig2}) where the tangents and normals are depicted an example representation of the normal and tangent vectors for the fluids are

\caption{Normals and tangents at interface for 2D system}

\begin{aligned}\boldsymbol{\tau} &= \begin{bmatrix}1 \\ 0\end{bmatrix} \\ \hat{\mathbf{n}} &= \begin{bmatrix}0 \\ 1\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(2.5)

For the traction vector

\begin{aligned}T_i = \sigma_{ij} n_j,\end{aligned} \hspace{\stretch{1}}(2.7)

we also have at the interface we must have matching of

\begin{aligned}\boldsymbol{\tau} \cdot \mathbf{T}.\end{aligned} \hspace{\stretch{1}}(2.8)

More explicitly, in coordinates this is

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

Steady incompressible rectilinear (unidirectional) flow.

In this case we can fix our axis so that

\begin{aligned}\mathbf{u} = \hat{\mathbf{x}} u(x, y, z, t),\end{aligned} \hspace{\stretch{1}}(3.10)

where the velocity components in the other directions

\begin{aligned}v &= 0 \\ w &= 0\end{aligned} \hspace{\stretch{1}}(3.11)

are both zero. Symbolically, the steady state condition is

\begin{aligned}\frac{\partial {\mathbf{u}}}{\partial {t}} = 0.\end{aligned} \hspace{\stretch{1}}(3.13)

We start with the incompressibility condition, which written explicitly, is

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


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

This implies

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

so our velocity can only be function of the y and z coordinates only

\begin{aligned}u = u(y, z).\end{aligned} \hspace{\stretch{1}}(3.17)

The non-linear term of the Navier-Stokes equation takes the form

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\left(u \frac{\partial {}}{\partial {x}}+\not{{v}} \frac{\partial {}}{\partial {y}}+\not{{w}} \frac{\partial {}}{\partial {z}} \right) (\hat{\mathbf{x}} u( y, z) + \hat{\mathbf{y}} \not{{v}} + \hat{\mathbf{z}} \not{{w}} ) \\ &= \hat{\mathbf{x}} u \not{{\frac{\partial {u}}{\partial {x}}}} \\ &= 0.\end{aligned}

With incompressibility and u = v = 0 conditions killing this term, and the steady state condition 3.13 killing the \rho {\partial {\mathbf{u}}}/{\partial {t}} term, the Navier-Stokes equation for this incompressible unidirectional steady state flow (in the absence of body forces) is reduced to

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

In coordinates this is


\begin{aligned}\frac{\partial {p}}{\partial {x}} = \mu \left( \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{\partial^2 {{u}}}{\partial {{z}}^2} \right) \end{aligned} \hspace{\stretch{1}}(3.19a)

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

\begin{aligned}\frac{\partial {p}}{\partial {z}} = 0.\end{aligned} \hspace{\stretch{1}}(3.19c)


Operating on the first with an x partial we find

\begin{aligned}\frac{\partial^2 {{p}}}{\partial {{x}}^2} = \mu \left( \frac{\partial^2 {{}}}{\partial {{y}}^2} \not{{\frac{\partial {u}}{\partial {x}}}} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \not{{ \frac{\partial {u}}{\partial {x}} }} \right) = 0\end{aligned} \hspace{\stretch{1}}(3.20)

Since we have

\begin{aligned}\frac{\partial^2 {{p}}}{\partial {{x}}^2} = 0\end{aligned} \hspace{\stretch{1}}(3.21)

we also have

\begin{aligned}\frac{d^2 p}{dx^2} = 0,\end{aligned} \hspace{\stretch{1}}(3.22)

so our pressure must be linear with position

\begin{aligned}p = A x + B,\end{aligned} \hspace{\stretch{1}}(3.23)

as illustrated in figure (\ref{fig:continuumL11:continuumL11fig3})

\caption{Pressure gradient in 1D system.}

\begin{aligned}p =\left\{\begin{array}{l l}p_0 & \quad \mbox{latex x = 0$} \\ p_L & \quad \mbox{x = L} \end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.24)$

we have

\begin{aligned}p = \frac{p_L - p_0}{L} x + p_0\end{aligned} \hspace{\stretch{1}}(3.25)


\begin{aligned}\frac{dp}{dx} = \frac{p_L - p_0}{L} = \text{constant} \equiv -G\end{aligned} \hspace{\stretch{1}}(3.26)

Example: Shearing flow.

The flows of this sort don’t have to be trivial. For example, even with constant pressure (p_0 = p_L) as in figure (\ref{fig:continuumL11:continuumL11fig4}) we can have a “shearing flow” where the fluids at the top surface are not necessarily moving at the same rates as the fluid below that surface. We have fluid flow in the x direction only, and our velocity is a function only of the y coordinate.

\caption{Velocity variation with height in shearing flow.}

\begin{aligned}\mathbf{u} &= \hat{\mathbf{x}} u(y) \\ G &= 0 \\ u(0) &= 0 \\ u(h) &= U.\end{aligned} \hspace{\stretch{1}}(3.27)

For such a flow 3.19a simplifies to

\begin{aligned}\frac{d^2 u}{dy^2} = 0\end{aligned} \hspace{\stretch{1}}(3.31)

with solution

\begin{aligned}u = \frac{U}{h} y + u(0) = \frac{U}{h} y.\end{aligned} \hspace{\stretch{1}}(3.32)

Example: Channel flow

\begin{aligned}\mathbf{u} &= \hat{\mathbf{x}} u(y) \\ G &= - \frac{dp}{dx} \ne 0\end{aligned} \hspace{\stretch{1}}(3.33)

This time our simplified Navier-Stokes equation 3.19a is reduced to something slightly more complicated

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

with solution

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

The boundary value conditions with the coordinate system in use illustrated in figure (\ref{fig:continuumL11:continuumL11fig5}) require the velocity to be zero at the interface (the pipe walls preventing flow in the interior of the pipe)

\caption{1D Channel flow coordinate system setup.}

\begin{aligned}u(\pm h) = -\frac{G}{2 \mu} h^2 \pm A h + B = 0\end{aligned} \hspace{\stretch{1}}(3.37)

One solution, immediately evident is,

\begin{aligned}A &= 0 \\ B &= \frac{G}{2 \mu} h^2,\end{aligned} \hspace{\stretch{1}}(3.38)

so our solution becomes

\begin{aligned}u = \frac{G}{2 \mu} \Bigl( h^2 - y^2 \Bigr),\end{aligned} \hspace{\stretch{1}}(3.40)

a parabolic velocity flow. This is illustrated graphically in figure (\ref{fig:continuumL11:continuumL11fig6}).

\caption{Parabolic velocity distribution.}

It’s clear that this is maximized by y = 0, but we can also see this by computing

\begin{aligned}\frac{du}{dy} = \frac{G}{\mu} y = 0.\end{aligned} \hspace{\stretch{1}}(3.41)

This maximum is

\begin{aligned}u_{\text{max}} = \frac{G}{2\mu} h^2\end{aligned} \hspace{\stretch{1}}(3.42)

The flux, or flow rate is

\begin{aligned}Q &= \iint_S \mathbf{u} \cdot \hat{\mathbf{x}} ds \\ &= \int_0^1 dz \int_{-h}^h dy u(y) \\ &=\frac{2 G h^3}{3}\end{aligned}

Let’s now compute the strain (e_{ij}) and the stress (\sigma_{ij} = -p \delta_{ij} + 2 \mu e_{ij})

\begin{aligned}e_{12} &= e_{21} = \frac{1}{{2}} \left( \frac{\partial {u}}{\partial {y}} \right) = - \frac{G y}{2 \mu} \\ e_{11} &= \frac{\partial {u}}{\partial {x}} = 0 \\ e_{22} &= \frac{\partial {v}}{\partial {y}} = 0\end{aligned} \hspace{\stretch{1}}(3.43)


\begin{aligned}\sigma_{12} = 2 \mu e_{12} = -G y\end{aligned} \hspace{\stretch{1}}(3.46)

This can be used to compute the forces on the inner surfaces of the tube. As illustrated in figure (\ref{fig:continuumL11:continuumL11fig7}), our normals at \pm h are \mp \hat{\mathbf{y}} respectively. The traction vector in the y direction is at y = h is

\caption{Normals in 1D channel flow system}

\begin{aligned}T_i = \sigma_{i 2} {\left.{{n_2}}\right\vert}_{{y = h}} = G h,\end{aligned} \hspace{\stretch{1}}(3.47)

so that

\begin{aligned}\mathbf{T} = \hat{\mathbf{x}} G h\end{aligned} \hspace{\stretch{1}}(3.48)

(here the x directionality comes from the i = 1 index of the stress tensor).

\begin{aligned}F_{x_L} = \hat{\mathbf{x}} \cdot \mathbf{T} = G h\end{aligned} \hspace{\stretch{1}}(3.49)

The total force is then

\begin{aligned}\int_0^L F_x dx = + G h L\end{aligned} \hspace{\stretch{1}}(3.50)

FIXME: ask in class. This is the tangential force at the boundary of the wall. What is it a force on? If it is tangential, how can it act on the wall? It could act on an impediment placed right up next to the wall, but if that’s the case, why are we integrating from x = 0 to x = L?

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

PHY456H1S Continuum mechanics. Problem Set 1. Stress, Strain, Traction vector. Force free equilibrium.

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


This problem set is as yet ungraded.

Problem Q1.


For the stress tensor

\begin{aligned}\sigma =\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}\text{M Pa}\end{aligned} \hspace{\stretch{1}}(2.1)

Find the corresponding strain tensor, assuming an isotropic solid with Young’s modulus E = 200 \times 10^9 \text{N}/\text{m}^2 and Poisson’s ration \nu = 0.35.


We need to express the relation between stress and strain in terms of Young’s modulus and Poisson’s ratio. In terms of Lam\’e parameters our model for the relations between stress and strain for an isotropic solid was given as

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

Computing the trace

\begin{aligned}\sigma_{kk} = (3 \lambda + 2 \mu) e_{kk},\end{aligned} \hspace{\stretch{1}}(2.3)

allows us to invert the relationship

\begin{aligned}2 \mu e_{ij} = \sigma_{ij} - \lambda \frac{\sigma_{kk}}{3 \lambda + 2 \mu} \delta_{ij}.\end{aligned} \hspace{\stretch{1}}(2.4)

In terms of Poisson’s ratio \nu and Young’s modulus E, our Lam\’e parameters were found to be

\begin{aligned}\lambda &= \frac{ E \nu }{(1 - 2 \nu)(1 + \nu)} \\ \mu &= \frac{E}{2(1 + \nu)},\end{aligned} \hspace{\stretch{1}}(2.5)


\begin{aligned}3 \lambda + 2 \mu&= \frac{ 3 E \nu }{(1 - 2 \nu)(1 + \nu)} + \frac{E}{1 + \nu} \\ &= \frac{E}{1 + \nu} \left( \frac{3 \nu}{1 - 2 \nu} + 1\right) \\ &= \frac{E}{1 + \nu} \frac{1 + \nu}{1 - 2 \nu} \\ &= \frac{E}{1 - 2 \nu}.\end{aligned}

Our stress strain model for the relationship for an isotropic solid becomes
we find

\begin{aligned}\frac{E}{1 + \nu} e_{ij}&=\sigma_{ij}-\frac{ E \nu }{(1 - 2 \nu)(1 + \nu)} \frac{1 - 2 \nu}{E}\sigma_{kk} \delta_{ij} \\ &=\sigma_{ij}-\frac{ \nu }{1 + \nu}\sigma_{kk} \delta_{ij} \\ \end{aligned}


\begin{aligned}e_{ij}=\frac{1}{{E}}\left((1 + \nu)\sigma_{ij}-\nu\sigma_{kk} \delta_{ij}\right).\end{aligned} \hspace{\stretch{1}}(2.7)

As a sanity check note that this matches (5.12) of [1], although they use a notation of \sigma instead of \nu for Poisson’s ratio. We are now ready to tackle the problem. First we need the trace of the stress tensor

\begin{aligned}\sigma_{kk} = (6 + 1 + 3) \text{M Pa} = 10 \text{M Pa},\end{aligned} \hspace{\stretch{1}}(2.8)

\begin{aligned}e_{ij}&=\frac{1}{{E}}\left((1 + \nu)\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}-10 \nu\begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}\right)\text{M Pa} \\ &=\frac{1}{{E}}\left(\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}+ 0.35\begin{bmatrix}-4 & 0 & 2 \\ 0 & -9 & 1 \\ 2 & 1 & -7\end{bmatrix}\right)\text{M Pa} \\ &=\frac{1}{{2 \times 10^{5}}}\left(\begin{bmatrix}6 & 0 & 2 \\ 0 & 1 & 1 \\ 2 & 1 & 3\end{bmatrix}+ 0.35\begin{bmatrix}-4 & 0 & 2 \\ 0 & -9 & 1 \\ 2 & 1 & -7\end{bmatrix}\right)\end{aligned}

Expanding out the last bits of arithmetic the strain tensor is found to have the form

\begin{aligned}e_{ij}=\begin{bmatrix} 23 & 0 & 13.5 \\  0 & -10.75 & 6.75 \\  13.5 & 6.75 & 2.75\end{bmatrix} 10^{-6}.\end{aligned} \hspace{\stretch{1}}(2.9)

Note that this is dimensionless, unlike the stress.

Problem Q2.


Small displacement field in a material is given by

\begin{aligned}e_1 &= 2 x_1 x_2 \\ e_2 &= x_3^2 \\ e_3 &= x_1^2 - x_3\end{aligned} \hspace{\stretch{1}}(3.10)


\item the infinitesimal strain tensor e_{ij},
\item the principal strains and the corresponding principal axes at (x_1, x_2, x_3) = (1, 2, 4),
\item Is the body under compression or expansion?

Solution. infinitesimal strain tensor e_{ij}

Diving right in, we have

\begin{aligned}e_{11}&= \frac{\partial {e_1}}{\partial {x_1}} \\ &= \frac{\partial {}}{\partial {x_1}}2 x_1 x_2 \\ &= 2 x_2\end{aligned}

\begin{aligned}e_{22}&= \frac{\partial {e_2}}{\partial {x_2}} \\ &= \frac{\partial {}}{\partial {x_2}} x_3^2 \\ &= 0\end{aligned}

\begin{aligned}e_{33}&= \frac{\partial {e_3}}{\partial {x_3}} \\ &= \frac{\partial {}}{\partial {x_3}} ( x_1^2 - x_3 ) \\ &= -1\end{aligned}

\begin{aligned}e_{12}&=\frac{1}{{2}} \left(\frac{\partial {e_2}}{\partial {x_1}}+\frac{\partial {e_1}}{\partial {x_2}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_1}} x_3^2 }}+\frac{\partial {}}{\partial {x_2}} 2 x_1 x_2\right) \\ &=x_1\end{aligned}

\begin{aligned}e_{23}&=\frac{1}{{2}} \left(\frac{\partial {e_3}}{\partial {x_2}}+\frac{\partial {e_2}}{\partial {x_3}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_2}} (x_1^2 - x_3 )}}+\frac{\partial {}}{\partial {x_3}} x_3^2\right) \\ &=x_3\end{aligned}

\begin{aligned}e_{31}&=\frac{1}{{2}} \left(\frac{\partial {e_1}}{\partial {x_3}}+\frac{\partial {e_3}}{\partial {x_1}}\right) \\ &=\frac{1}{{2}}\left(\not{{\frac{\partial {}}{\partial {x_3}} 2 x_1 x_2 }}+\frac{\partial {}}{\partial {x_1}} (x_1^2 - x_3 )\right) \\ &=x_1\end{aligned}

In matrix form we have

\begin{aligned}\mathbf{e} =\begin{bmatrix}2 x_2 & x_1 & x_1 \\ x_1 & 0 & x_3 \\ x_1 & x_3 & -1 \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.13)

Solution. principle strains and axes

At the point (1, 2, 4) the strain tensor has the value

\begin{aligned}\mathbf{e} =\begin{bmatrix}4 & 1 & 1 \\ 1 & 0 & 4 \\ 1 & 4 & -1\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(3.14)

We wish to diagonalize this, solving the characteristic equation for the eigenvalues \lambda

\begin{aligned}0 &=\begin{vmatrix}4 -\lambda & 1 & 1 \\ 1 & -\lambda & 4 \\ 1 & 4 & -1 -\lambda\end{vmatrix} \\ &=(4 -\lambda )\begin{vmatrix} -\lambda & 4 \\  4 & -1 -\lambda\end{vmatrix}-\begin{vmatrix}1 & 1 \\ 4 & -1 -\lambda\end{vmatrix}+\begin{vmatrix}1 & 1 \\ -\lambda & 4 \\ \end{vmatrix} \\ &=(4 - \lambda)(\lambda^2 + \lambda - 16)-(-1 -\lambda - 4)+(4 + \lambda) \\ \end{aligned}

We find the characteristic equation to be

\begin{aligned}0 = -\lambda^3 + 3 \lambda^2 + 22\lambda - 55.\end{aligned} \hspace{\stretch{1}}(3.15)

This doesn’t appear to lend itself easily to manual solution (there are no obvious roots to factor out). As expected, since the matrix is symmetric, a plot (\ref{fig:continuumL8:continuumProblemSet1Q2fig1}) shows that all our roots are real

\caption{Q2. Characteristic equation.}

Numerically, we determine these roots to be

\begin{aligned}\{5.19684, -4.53206, 2.33522\}\end{aligned} \hspace{\stretch{1}}(3.16)

with the corresponding basis (orthonormal eigenvectors), the principle axes are

\begin{aligned}\left\{\hat{\mathbf{p}}_1,\hat{\mathbf{p}}_2,\hat{\mathbf{p}}_3\right\}=\left\{\begin{bmatrix}0.76291 \\ 0.480082 \\ 0.433001\end{bmatrix},\begin{bmatrix}-0.010606 \\ -0.660372 \\ 0.750863\end{bmatrix},\begin{bmatrix}-0.646418 \\ 0.577433 \\ 0.498713\end{bmatrix}\right\}.\end{aligned} \hspace{\stretch{1}}(3.17)

Solution. Is body under compression or expansion?

To consider this question, suppose that as in the previous part, we determine a basis for which our strain tensor e_{ij} = p_i \delta_{ij} is diagonal with respect to that basis at a given point \mathbf{x}_0. We can then simplify the form of the stress tensor at that point in the object

\begin{aligned}\sigma_{ij}&=\frac{E}{1 + \nu} \left(e_{ij} + \frac{\nu}{1 - 2 \nu} e_{mm} \delta_{ij}\right) \\ &=\frac{E}{1 + \nu} \left(p_i + \frac{\nu}{1 - 2 \nu} e_{mm}\right)\delta_{ij}.\end{aligned}

We see that the stress tensor at this point is also necessarily diagonal if the strain is diagonal in that basis (with the implicit assumption here that we are talking about an isotropic material). Noting that the Poisson ratio is bounded according to

\begin{aligned}-1 \le \nu \le \frac{1}{{2}},\end{aligned} \hspace{\stretch{1}}(3.18)

so if our trace is positive (as it is in this problem for all points x_2 > 1/2), then any positive principle strain value will result in a positive stress along that direction). For example at the point (1,2,4) of the previous part of this problem (for which x_2 > 1/2), we have

\begin{aligned}\sigma_{ij}=\frac{E}{1 + \nu}\begin{bmatrix}5.19684+ \frac{3 \nu}{1 - 2 \nu}  & 0 & 0 \\ 0 & -4.53206+ \frac{3 \nu}{1 - 2 \nu}  & 0 \\ 0 & 0 & 2.33522+ \frac{3 \nu}{1 - 2 \nu}\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(3.19)

We see that at this point the (1,1) and (3,3) components of stress is positive (expansion in those directions) regardless of the material, and provided that

\begin{aligned}\frac{3 \nu}{1 - 2 \nu} > 4.53206\end{aligned} \hspace{\stretch{1}}(3.20)

(i.e. \nu > 0.375664) the material is under expansion in all directions. For \nu < 0.375664 the material at that point is expanding in the \hat{\mathbf{p}}_1 and \hat{\mathbf{p}}_3 directions, but under compression in the \hat{\mathbf{p}}_2 directions.

(save to disk and run with either Mathematica or the free Wolfram CDF player ( ) )

For a Mathematica notebook that visualizes this part of this problem see This animates the stress tensor associated with the problem, for different points (x,y,z) and values of Poisson’s ratio \nu, with Mathematica manipulate sliders available to alter these (as well as a zoom control to scale the graphic, keeping the orientation and scale fixed with any variation of the other parameters). This generalizes the solution of the problem (assuming I got it right for the specific (1,2,4) point of the problem). The vectors are the orthonormal eigenvectors of the tensor, scaled by the magnitude of the eigenvectors of the stress tensor (also diagonal in the basis of the diagonalized strain tensor at the point in question). For those directions that are under expansive stress, I’ve colored the vectors blue, and for compressive directions, I’ve colored the vectors red.

This requires either a Mathematica client or the free Wolfram CDF player, either of which can run the notebook after it is saved to your computer’s hard drive.

Problem Q3.


The stress tensor at a point has components given by

\begin{aligned}\sigma =\begin{bmatrix}1 & -2 & 2 \\ -2 & 3 & 1 \\ 2 & 1 & -1\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(4.21)

Find the traction vector across an area normal to the unit vector

\begin{aligned}\hat{\mathbf{n}} = ( \sqrt{2} \mathbf{e}_1 - \mathbf{e}_2 + \mathbf{e}_3)/2\end{aligned} \hspace{\stretch{1}}(4.22)

Can you construct a tangent vector \boldsymbol{\tau} on this plane by inspection? What are the components of the force per unit area along the normal \hat{\mathbf{n}} and tangent \boldsymbol{\tau} on that surface? (hint: projection of the traction vector.)


The traction vector, the force per unit volume that holds a body in equilibrium, in coordinate form was

\begin{aligned}P_i = \sigma_{ik} n_k\end{aligned} \hspace{\stretch{1}}(4.23)

where n_k was the coordinates of the normal to the surface with area df_k. In matrix form, this is just

\begin{aligned}\mathbf{P} = \sigma \hat{\mathbf{n}},\end{aligned} \hspace{\stretch{1}}(4.24)

so our traction vector for this stress tensor and surface normal is just

\begin{aligned}\mathbf{P} &=\frac{1}{{2}}\begin{bmatrix}1 & -2 & 2 \\ -2 & 3 & 1 \\ 2 & 1 & -1\end{bmatrix}\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{2}}\begin{bmatrix}\sqrt{2} + 2 + 2 \\ -2\sqrt{2} - 3 + 1 \\ 2\sqrt{2} - 1 -1\end{bmatrix} \\ &=\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}\end{aligned}

We also want a vector in the plane, and can pick

\begin{aligned}\boldsymbol{\tau} = \frac{1}{{\sqrt{2}}}\begin{bmatrix}0 \\ 1 \\ 1\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(4.25)


\begin{aligned}\boldsymbol{\tau}' = \begin{bmatrix}\frac{1}{{\sqrt{2}}} \\ \frac{1}{{2}} \\ -\frac{1}{{2}}\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(4.26)

It’s clear that either of these is normal to \hat{\mathbf{n}} (the first can also be computed by normalizing \hat{\mathbf{n}} \times \mathbf{e}_1, and the second with one round of Gram-Schmidt). However, neither of these vectors in the plane are particularly interesting since they are completely arbitrary. Let’s instead compute the projection and rejection of the traction vector with respect to the normal. We find for the projection

\begin{aligned}(\mathbf{P} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}&=\frac{1}{{4}}\left(\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}\cdot \begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \right)\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix}  \\ &=\frac{1}{{4}}\left( 1 + 2\sqrt{2}+\sqrt{2} +1 +\sqrt{2} - 1\right)\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix}  \\ &=\frac{1}{{2}}\left( 1 + 4\sqrt{2}\right)\hat{\mathbf{n}}\end{aligned}

Our rejection, the component of the traction vector in the plane, is

\begin{aligned}(\mathbf{P} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}} &=\mathbf{P} - (\mathbf{P} \cdot \hat{\mathbf{n}})\hat{\mathbf{n}} \\ &=\frac{1}{{2}}\begin{bmatrix}\sqrt{2}/2 + 2 \\ -\sqrt{2} -1 \\ \sqrt{2} - 1\end{bmatrix}-\frac{1}{{4}}(1 + r \sqrt{2})\begin{bmatrix}\sqrt{2} \\ -1 \\ 1\end{bmatrix} \\ &=\frac{1}{{4}}\begin{bmatrix}\sqrt{2} \\ -3 \\ -5\end{bmatrix}\end{aligned}

This gives us a another vector perpendicular to the normal \hat{\mathbf{n}}

\begin{aligned}\hat{\boldsymbol{\tau}} = \frac{1}{{6}}\begin{bmatrix}\sqrt{2} \\ -3 \\ -5\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(4.27)

Wrapping up, we find the decomposition of the traction vector in the direction of the normal and its projection onto the plane to be

\begin{aligned}\mathbf{P} = \frac{1}{{2}}(1 + 4\sqrt{2}) \hat{\mathbf{n}}+\frac{3}{2} \hat{\boldsymbol{\tau}}.\end{aligned} \hspace{\stretch{1}}(4.28)

The components we can read off by inspection.

Problem Q4.


The stress tensor of a body is given by

\begin{aligned}\sigma =\begin{bmatrix}A \cos x & y^2 & C x \\ y^2 & B \sin y & z \\ C x & z & z^3\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(5.29)

Determine the constant A, B, and C if the body is in equilibrium.


In the absence of external forces our equilibrium condition was

\begin{aligned}\partial_k \sigma_{ik} = 0.\end{aligned} \hspace{\stretch{1}}(5.30)

In matrix form we wish to operate (to the left) with the gradient coordinate vector

\begin{aligned}0 &= \sigma \stackrel{ \leftarrow }{\boldsymbol{\nabla}} \\ &=\begin{bmatrix}A \cos x & y^2 & C x \\ y^2 & B \sin y & z \\ C x & z & z^3\end{bmatrix}\begin{bmatrix}\stackrel{ \leftarrow }{\partial}_x \\ \stackrel{ \leftarrow }{\partial}_y \\ \stackrel{ \leftarrow }{\partial}_z \\ \end{bmatrix} \\ &=\begin{bmatrix}\partial_x (A \cos x) + \partial_y(y^2) + \not{{\partial_z(C x)}} \\ \not{{\partial_x (y^2)}} + \partial_y(B \sin y) + \partial_z(z) \\ \partial_x (C x) + \not{{\partial_y(z)}} + \partial_z(z^3)\end{bmatrix} \\ &=\begin{bmatrix}-A \sin x + 2 y \\ B \cos y + 1 \\ C + 3 z^2 \end{bmatrix} \\ \end{aligned}

So, our conditions for equilibrium will be satisfied when we have

\begin{aligned}A &= \frac{2 y }{\sin x} \\ B &= -\frac{1}{\cos y} \\ C &= -3 z^2,\end{aligned} \hspace{\stretch{1}}(5.31)

provided y \ne 0, and y \ne \pi/2 + n\pi for integer n. If equilibrium is to hold along the y = 0 plane, then we must either also have A = 0 or also impose the restriction x = m \pi (for integer m).

A couple other mathematica notebooks

Some of the hand calculations done in this problem set I’ve confirmed using Mathematica. Those notebooks are available here


These all require either a Mathematica client or the free Wolfram CDF player. Note that I haven’t figured out a way to get a browser based CDF player to play these without explicit download.


[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.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , | 10 Comments »

PHY454H1S Continuum Mechanics. Lecture 5: Constitutive relationship. Taught by Prof. K. Das.

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


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

Review: Cauchy Tetrahedron.

Referring to figure (\ref{fig:continuumL5:continuumL5fig1})
\caption{Cauchy tetrahedron direction cosines.}

recall that we can decompose our force into components that refer to our direction cosines n_i = \cos\phi_i

\begin{aligned}f_1 &= \sigma_{11} n_1 + \sigma_{12} n_2 + \sigma_{13} n_3 \\ f_2 &= \sigma_{21} n_1 + \sigma_{22} n_2 + \sigma_{23} n_3 \\ f_3 &= \sigma_{31} n_1 + \sigma_{32} n_2 + \sigma_{33} n_3\end{aligned} \hspace{\stretch{1}}(2.1)

Or in tensor form

\begin{aligned}f_i = \sigma_{ij} n_j.\end{aligned} \hspace{\stretch{1}}(2.4)

We call this the traction vector and denote it in vector form as

\begin{aligned}\mathbf{T} = \boldsymbol{\sigma} \cdot \hat{\mathbf{n}}=\begin{bmatrix}\sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33}\end{bmatrix}\begin{bmatrix}n_1 \\ n_2 \\ n_3\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.5)

Constitutive relation.

Reading: section 2, section 4 and section 5 from the text [1].

We can find the relationship between stress and strain, both analytically and experimentally, and call this the Constitutive relation. We prefer to deal with ranges of distortion that are small enough that we can make a linear approximation for this relation. In general such a linear relationship takes the form

\begin{aligned}\sigma_{ij} = c_{ijkl} e_{kl}.\end{aligned} \hspace{\stretch{1}}(3.6)

Consider the number of components that we are talking about for various rank tensors

\begin{aligned}\begin{array}{l l}\mbox{latex 0^\text{th}$ rank tensor} & \mbox{3^0 = 1 components} \\ \mbox{1^\text{st} rank tensor} & \mbox{3^1 = 3 components} \\ \mbox{2^\text{nd} rank tensor} & \mbox{3^2 = 9 components} \\ \mbox{3^\text{rd} rank tensor} & \mbox{3^3 = 81 components}\end{array}\end{aligned} \hspace{\stretch{1}}(3.7)$

We have a lot of components, even for a linear relation between stress and strain. For isotropic materials we model the constitutive relation instead as

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

For such a modeling of the material the (measured) values \lambda and \mu (shear modulus or modulus of rigidity) are called the Lam\’e parameters.

It will be useful to compute the trace of the stress tensor in the form of the constitutive relation for the isotropic model. We find

\begin{aligned}\sigma_{ii}&= \lambda e_{kk} \delta_{ii} + 2 \mu e_{ii} \\ &= 3 \lambda e_{kk} + 2 \mu e_{jj},\end{aligned}


\begin{aligned}\sigma_{ii} = (3 \lambda + 2 \mu) e_{kk}.\end{aligned} \hspace{\stretch{1}}(3.9)

We can now also invert this, to find the trace of the strain tensor in terms of the stress tensor

\begin{aligned}e_{ii} = \frac{\sigma_{kk}}{3 \lambda + 2 \mu}\end{aligned} \hspace{\stretch{1}}(3.10)

Substituting back into our original relationship 3.8, and find

\begin{aligned}\sigma_{ij} = \lambda \frac{\sigma_{kk}}{3 \lambda + 2 \mu} \delta_{ij} + 2 \mu e_{ij},\end{aligned} \hspace{\stretch{1}}(3.12)

which finally provides an inverted expression with the strain tensor expressed in terms of the stress tensor

\begin{aligned}\boxed{2 \mu e_{ij} =\sigma_{ij} - \lambda \frac{\sigma_{kk}}{3 \lambda + 2 \mu} \delta_{ij}.}\end{aligned} \hspace{\stretch{1}}(3.12)

Special cases.

Hydrostatic compression

Hydrostatic compression is when we have no shear stress, only normal components of the stress matrix \sigma_{ij} is nonzero. Strictly speaking we define Hydrostatic compression as

\begin{aligned}\sigma_{ij} = -p \delta_{ij},\end{aligned} \hspace{\stretch{1}}(3.13)

i.e. not only diagonal, but with all the components of the stress tensor equal.

We can write the trace of the stress tensor as

\begin{aligned}\sigma_{ii} = - 3 p = (3 \lambda + 2 \mu) e_{kk}.\end{aligned} \hspace{\stretch{1}}(3.14)

Now, from our discussion of the strain tensor e_{ij} recall that we found in the limit

\begin{aligned}dV' = (1 + e_{ii}) dV,\end{aligned} \hspace{\stretch{1}}(3.15)

allowing us to express the change in volume relative to the original volume in terms of the strain trace

\begin{aligned}e_{ii} = \frac{dV' - dV}{dV}.\end{aligned} \hspace{\stretch{1}}(3.16)

Writing that relative volume difference as \Delta V/V we find

\begin{aligned}- 3 p = (3 \lambda + 2 \mu) \frac{\Delta V}{V},\end{aligned} \hspace{\stretch{1}}(3.17)


\begin{aligned}- \frac{ p V}{\Delta V} = \left( \lambda + \frac{2}{3} \mu \right) = K,\end{aligned} \hspace{\stretch{1}}(3.18)

where K is called the Bulk modulus.

Uniaxial stress

Again illustrated in the plane as in figure (\ref{fig:continuumL5:continuumL5fig2})
\caption{Uniaxial stress.}

Expanding out 3.12 we have for the 1,1 element of the strain tensor

\begin{aligned}\boldsymbol{\sigma} =\begin{bmatrix}\sigma_{11} & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.19)

\begin{aligned}2 \mu e_{11}&= \sigma_{11} - \frac{\lambda ( \sigma_{11} + \not{{\sigma_{22}}} ) }{3 \lambda + 2 \mu} \\ &= \sigma_{11} \frac{3 \lambda + 2 \mu - \lambda }{3 \lambda + 2 \mu} \\ &= 2 \sigma_{11} \frac{\lambda + \mu }{3 \lambda + 2 \mu}\end{aligned}


\begin{aligned}\frac{\sigma_{11}}{e_{11}} = \frac{\mu(3 \lambda + 2 \mu)}{\lambda + \mu } = E\end{aligned} \hspace{\stretch{1}}(3.20)

where E is Young’s modulus. Young’s modulus in the text (5.3) is given in terms of the bulk modulus K. Using \lambda = K - 2\mu/3 we find

\begin{aligned}E &=\frac{\mu(3 \lambda + 2 \mu)}{\lambda + \mu } \\ &=\frac{\mu(3 (K - 2\mu/3)+ 2 \mu)}{K - 2\mu/3 + \mu } \\ &=\frac{3 K \mu}{ K + \mu/3 } \end{aligned}

\begin{aligned}\boxed{E =\frac{\mu(3 \lambda + 2 \mu)}{\lambda + \mu } =\frac{9 K \mu}{ 3 K + \mu } }\end{aligned} \hspace{\stretch{1}}(3.21)

FIXME: figure (\ref{fig:continuumL5:continuumL5fig3}) reference?

\caption{stress associated with Young’s modulus}

We define Poisson’s ratio \nu as the quantity

\begin{aligned}\frac{e_{22}}{e_{11}} = \frac{e_{33}}{e_{11}} = - \nu.\end{aligned} \hspace{\stretch{1}}(3.22)

Note that we are still talking about uniaxial stress here. Referring back to 3.12 we have

\begin{aligned}2 \mu e_{2 2}&= \sigma_{2 2} - \lambda \frac{\sigma_{k k}}{3 \lambda + 2 \mu} \delta_{2 2} \\ &= \sigma_{2 2} - \lambda \frac{\sigma_{k k}}{3 \lambda + 2 \mu} \\ &= - \frac{\lambda \sigma_{11}}{3 \lambda + 2 \mu}\end{aligned}

Recall (3.20) that we had

\begin{aligned}\sigma_{11} = \frac{\mu (3 \lambda + 2 \mu)}{\lambda + \mu} e_{11}.\end{aligned} \hspace{\stretch{1}}(3.23)

Inserting this gives us

\begin{aligned}2 \mu e_{22} = - \frac{\lambda}{\not{{3 \lambda + 2 \mu}}} \frac{ \mu (\not{{3 \lambda + 2\mu}})}{\lambda + \mu} e_{11}\end{aligned}


\begin{aligned}\boxed{\nu = -\frac{e_{22}}{e_{11}} = \frac{\lambda}{2 (\lambda + \mu)}.}\end{aligned} \hspace{\stretch{1}}(3.24)

We can also relate the Poisson’s ratio \nu to the shear modulus \mu

\begin{aligned}\mu = \frac{E}{2(1 + \nu)}\end{aligned} \hspace{\stretch{1}}(3.25)

\begin{aligned}\lambda = \frac{E \nu}{(1 - 2 \nu)(1 + \mu)}\end{aligned} \hspace{\stretch{1}}(3.26)

\begin{aligned}e_{11} &= \frac{1}{{E}}\left( \sigma_{11} - \nu(\sigma_{22} + \sigma_{33}) \right) \\ e_{22} &= \frac{1}{{E}}\left( \sigma_{22} - \nu(\sigma_{11} + \sigma_{33}) \right) \\ e_{33} &= \frac{1}{{E}}\left( \sigma_{33} - \nu(\sigma_{11} + \sigma_{22}) \right)\end{aligned} \hspace{\stretch{1}}(3.27)

These ones are (5.14) in the text, and are easy enough to verify (not done here).

Appendix. Computing the relation between Poisson’s ratio and shear modulus.

Young’s modulus is given in 3.21 (equation (43) in the Professor’s notes) as

\begin{aligned}E = \frac{\mu(3 \lambda + 2 \mu)}{\lambda + \mu },\end{aligned} \hspace{\stretch{1}}(3.30)

and for Poisson’s ratio 3.24 (equation (46) in the Professor’s notes) we have

\begin{aligned}\nu = -\frac{e_{22}}{e_{11}} = \frac{\lambda}{2 (\lambda + \mu)}.\end{aligned} \hspace{\stretch{1}}(3.31)

Let’s derive the other stated relationships (equation (47) in the Professor’s notes). I get

\begin{aligned}2 (\lambda + \mu) \nu = \lambda \\ \implies \\ \lambda ( 2 \nu - 1 ) = - 2\mu\nu\end{aligned}


\begin{aligned}\lambda = \frac{ 2 \mu \nu} { 1 - 2 \nu }\end{aligned}

For substitution into the Young’s modulus equation calculate

\begin{aligned}\lambda + \mu &= \frac{ 2 \mu \nu} { 1 - 2 \nu } + \mu \\ &= \mu \left( \frac{ 2 \nu} { 1 - 2 \nu } + 1 \right)  \\ &= \mu \frac{ 2 \nu + 1 - 2 \nu} { 1 - 2 \nu }  \\ &= \frac{ \mu} { 1 - 2 \nu }  \\ \end{aligned}


\begin{aligned}3 \lambda + 2 \mu &= 3 \frac{ \mu} { 1 - 2 \nu } - \mu \\ &= \mu \frac{ 3 - (1 - 2 \nu)} { 1 - 2 \nu } \\ &= \mu \frac{ 2 + 2 \nu} { 1 - 2 \nu } \\ &= 2 \mu \frac{ 1 + \nu} { 1 - 2 \nu } \\ \end{aligned}

Putting these together we find

\begin{aligned}E &= \frac{\mu(3 \lambda + 2 \mu)}{\lambda + \mu } \\ &= \mu 2 \mu \frac{ 1 + \nu} { 1 - 2 \nu } \frac{ 1 - 2 \nu}{\mu} \\ &= 2 \mu ( 1 + \nu ) \\ \end{aligned}

Rearranging we have

\begin{aligned}\mu = \frac{E}{2 (1 + \nu)}.\end{aligned} \hspace{\stretch{1}}(3.32)

This matches (5.9) in the text (where \sigma is used instead of \nu).

We also find

\begin{aligned}\lambda &= \frac{ 2 \mu \nu} { 1 - 2 \nu } \\ &= \frac{ \nu} { 1 - 2 \nu } \frac{E }{1 + \nu}.\end{aligned}


[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.

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