Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘shear flow’

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

Motivation.

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

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

\end{subequations}

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

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

\end{subequations}

or

\begin{subequations}

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

\end{subequations}

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

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

\end{subequations}

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

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

\end{subequations}

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.

Flux

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.

Hydrostatics.

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)

or

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

or

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

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

\end{subequations}

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)

or

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

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

\end{subequations}

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

\begin{aligned}u \frac{\partial {u}}{\partial {x}} + v \frac{\partial {u}}{\partial {y}} = U \frac{dU}{dx} + \nu \frac{\partial^2 {{u}}}{\partial {{y}}^2}\end{aligned} \hspace{\stretch{1}}(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)

\end{subequations}

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)

and

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

Stability.

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

\begin{aligned}\boldsymbol{\nabla}^2 \left( \frac{\partial {}}{\partial {t}} - \boldsymbol{\nabla}^2 \right) \delta w = \mathcal{R} \left( \frac{\partial^2 {{}}}{\partial {{x}}^2} +\frac{\partial^2 {{}}}{\partial {{y}}^2} \right) \delta T\end{aligned} \hspace{\stretch{1}}(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)

\end{subequations}

where we’ve introduced the Rayleigh number and Prandtl number’s
\begin{subequations}

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

\end{subequations}

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.

References

[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]. http://en.wikipedia.org/w/index.php?title=Curvature&oldid=488021394.

[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]. http://en.wikipedia.org/w/index.php?title=Blasius_boundary_layer&oldid=480776115.

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

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

Posted by peeterjoot on March 26, 2012

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

Problem Q3 (revisited).

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

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

Exact solutions.

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

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

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

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

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

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

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

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

and matching velocities at y = 0 gives us

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

This gives us

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

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

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

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

PHY454H1S Continuum mechanics midterm reflection.

Posted by peeterjoot on March 17, 2012

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

Motivation.

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

Problem 1.

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

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

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

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

we obtain

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

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

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

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

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

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

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

so the displaced volume is

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

Since

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

We have

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

or

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

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

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

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

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

\mathbf{S}-waves are transverse.

\item Whose speed is higher?

Answer.
From the formula sheet we have

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

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

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

Answer.

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

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

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

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

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

for isotropic solids we model this as

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

and for Newtonian fluids

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

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

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

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

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

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

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

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

An incompressible fluid has

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

but since we also have

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

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

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

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

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

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

or

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

so that in differential form we have

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

Expanding the divergence by chain rule we have

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

but this is just

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

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

\end{itemize}

Problem 2.

Statement.

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

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

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

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

Answer

Our equations of motion are

\begin{subequations}

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

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

\end{subequations}

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

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

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

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

which leaves us with just

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

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

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

Integrate once

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

and once more to find the velocity

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

Let’s incorporate an additional constant into B'

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

so that we have

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

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

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

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

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

Adding these we find

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

and subtracting find

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

Our velocity is

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

or rearranged a bit

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

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

Answer

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

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

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

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

Answer

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

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

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

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

Answer

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

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

Since our acceleration is

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

our extreme values occur at

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

At this point, our velocity is

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

or just

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

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

An element of our volume flux is

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

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

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

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

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

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

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

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

Our traction vector is

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

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

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

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

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

so at y = -h we have

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

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

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

\end{itemize}

References

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

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

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

Compilation of class notes for phy454h1s, continuum mechanics (so far).

Posted by peeterjoot on February 17, 2012

Have collected all my pre-midterm continuum mechanics notes into a single document. The individual pdfs below are still available, but won’t be updated further.

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

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

Feb 10, 2012 Navier-Stokes equation.

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

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

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

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

Jan 27, 2012 Compatibility condition and elastostatics.

Jan 25, 2012 Constitutive relationship.

Jan 20, 2012 Strain tensor components.

Jan 18, 2012 Strain tensor review. Stress tensor.

Jan 13, 2012 Introduction and strain tensor.

Jan 11, 2012 Overview.

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

PHY454H1S Continuum Mechanics. Lecture 12: Flow in a pipe. Gravity driven flow of a film. Taught by Prof. K. Das.

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

Disclaimer.

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

Review. Steady rectilinear flow.

Steady:

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

Rectilinear is a unidirectional flow such as

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

\begin{enumerate}
\item
Utilizing an incompressibility assumption \boldsymbol{\nabla} \cdot \mathbf{u} = 0, so for this case we have

\begin{aligned}\frac{\partial {u}}{\partial {x}} = 0\end{aligned}

or

\begin{aligned}u = u(y, z)\end{aligned}

Note that Prof. Das called this a continuity requirement, and justified this label with the relation

\begin{aligned}\frac{d\rho}{dt} = \rho (\boldsymbol{\nabla} \cdot \mathbf{u}),\end{aligned} \hspace{\stretch{1}}(2.3)

which was a consequence of mass conservation. It’s still not clear to me why he would call this a continuity requirement.

\item Nonlinear term is zero. (\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = 0
\item p = p(x). Since \frac{d^2 p}{dx^2} = 0 we also have \frac{dp}{dx} = -G, a constant.

\item \mu \left( \frac{\partial^2 {{u}}}{\partial {{y}}^2} + \frac{\partial^2 {{u}}}{\partial {{z}}^2} \right) = G

\end{enumerate}

Solution by intuition.

Two examples that we have solved analytically are illustrated in figure (\ref{fig:continuumL12:continuumL12fig1}) and figure (\ref{fig:continuumL12:continuumL12fig2})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig1}
\caption{Simple shear flow}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig2}
\caption{Channel flow}
\end{figure}

Sometimes we can utilize solutions already found to understand the behaviour of more complex systems. Combining the two we can look at flow over a plate as in figure (\ref{fig:continuumL12:continuumL12fig3})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig3}
\caption{Flow on a plate}
\end{figure}

Example 2. Fluid in a container. If the surface tension is altered on one side, we induce a flow on the surface, leading to a circulation flow. This can be done for example, by introducing a heat source or addition of surfactant.

This is illustrated in figure (\ref{fig:continuumL12:continuumL12fig4})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig4}
\caption{Circulation flow induced by altering surface tension.}
\end{figure}

This sort of flow is hard to analyze, only first done by Steve Davis in the 1980′s. The point here is that we can use some level of intuition to guide our attempts at solution.

Flow down a pipe.

Reading: section 2 from [1].

Recall that the Navier-Stokes equation is

\begin{aligned}\boxed{\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}}(4.4)

We need to express this in cylindrical coordinates (r, \theta, z) as in figure (\ref{fig:continuumL12:continuumL12fig5})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig5}
\caption{Flow through a pipe.}
\end{figure}

Our gradient is

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

For our Laplacian we find

\begin{aligned}\boldsymbol{\nabla}^2 &= \left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right) \cdot\left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right) \\ &=\partial_{rr} + \frac{\hat{\boldsymbol{\theta}}}{r} \cdot (\partial_\theta \hat{\mathbf{r}}) \partial_r+ \frac{1}{{r}} \partial_\theta \left( \frac{1}{{r}} \partial_\theta \right)+ \partial_{zz} \\ &=\partial_{rr} + \frac{1}{{r}} \partial_r + \frac{1}{{r^2}} \partial_{\theta\theta} + \partial_{zz},\end{aligned}

which we can write as

\begin{aligned}\boldsymbol{\nabla}^2 = \frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2}.\end{aligned} \hspace{\stretch{1}}(4.6)

NS takes the form

\begin{aligned}\boxed{\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} &+ \rho \left(u_r \frac{\partial {}}{\partial {r}} + \frac{u_\theta}{r} \frac{\partial {}}{\partial {\theta}} + u_z \frac{\partial {}}{\partial {z}} \right) \mathbf{u} = \\ &- \left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right)p + \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)\mathbf{u} + \rho \mathbf{f}.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(4.7)

For steady state and incompressible fluids in the absence of body forces we have

\begin{aligned}\left(\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{\hat{\boldsymbol{\theta}}}{r} \frac{\partial {}}{\partial {\theta}} + \hat{\mathbf{z}} \frac{\partial {}}{\partial {z}}\right)p = \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)\mathbf{u},\end{aligned} \hspace{\stretch{1}}(4.8)

or, in coordinates

\begin{aligned}\frac{\partial {p}}{\partial {r}} &= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_r \\ \frac{1}{r} \frac{\partial {p}}{\partial {\theta}}&= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_\theta \\ \frac{\partial {p}}{\partial {z}}&= \mu \left(\frac{1}{{r}} \frac{\partial {}}{\partial {r}} \left( r \frac{\partial {}}{\partial {r}} \right) + \frac{1}{{r^2}} \frac{\partial^2 {{}}}{\partial {{\theta}}^2} + \frac{\partial^2 {{}}}{\partial {{z}}^2} \right)u_z\end{aligned} \hspace{\stretch{1}}(4.9)

With an assumption that we have no radial or circulatory flows (u_r = u_\theta = 0), and with u_z = w assumed to only have a radial dependence, our velocity is

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

and an assumption of linear pressure dependence

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

then NS takes the final simple form

\begin{aligned}\frac{1}{{r}} \frac{d}{dr} \left( r \frac{dw}{dr} \right) = - \frac{G}{\mu}.\end{aligned} \hspace{\stretch{1}}(4.14)

Solving this we have

\begin{aligned}r \frac{dw}{dr} = - \frac{G r^2}{2\mu} + A\end{aligned} \hspace{\stretch{1}}(4.15)

\begin{aligned}w = -\frac{G r^2}{4 \mu} + A \ln(r) + B\end{aligned} \hspace{\stretch{1}}(4.16)

Requiring finite solutions for r = 0 means that we must have A = 0. Also w(a) = 0, we have B = G a^2/4 \mu so we must have

\begin{aligned}w(r) = \frac{G}{4 \mu}( a^2 - r^2 )\end{aligned} \hspace{\stretch{1}}(4.17)

Example: Gravity driven flow of a liquid film

(This is one of our Professor’s favorite problems).

Coordinates as in figure (\ref{fig:continuumL12:continuumL12fig6})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig6}
\caption{Gravity driven flow down an inclined plane.}
\end{figure}

\begin{aligned}\mathbf{u} = \hat{\mathbf{x}} u(y)\end{aligned} \hspace{\stretch{1}}(5.18)

Boundary conditions

\begin{enumerate}
\item u(y = 0) = 0
\item Tangential stress at the air-liquid interface y = h is equal.

\begin{aligned}\boldsymbol{\tau} \cdot (\boldsymbol{\sigma}_l \cdot \hat{\mathbf{n}}) = \boldsymbol{\tau} \cdot (\boldsymbol{\sigma}_a \cdot \hat{\mathbf{n}}),\end{aligned} \hspace{\stretch{1}}(5.19)

\end{enumerate}

We write

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

and seek simultaneous solutions to the pair of stress tensor equations

\begin{aligned}\sigma_{ij}^l &= - p \delta_{ij} + \mu^l \left( \frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_j}}{\partial {x_i}}\right) \\ \sigma_{ij}^a &= - p \delta_{ij} + \mu^a \left( \frac{\partial {u_i}}{\partial {x_j}} +\frac{\partial {u_j}}{\partial {x_i}}\right).\end{aligned} \hspace{\stretch{1}}(5.23)

In general this requires an iterated approach, solving for one with an initial approximation of the other, then switching and tuning the numerical method carefully for convergence.

We expect that the flow of liquid will induce a flow of air at the interface, but may be able to make a one-sided approximation. Let’s see how far we get before we have to introduce any approximations and compute the traction vector for the liquid

\begin{aligned}\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}} &= \begin{bmatrix}-p & \mu^l {\partial {u}}/{\partial {y}} & 0 \\ \mu^l {\partial {u}}/{\partial {y}} & -p & 0 \\ 0 & 0 & 0\end{bmatrix}\begin{bmatrix}0 \\ 1 \\ 0\end{bmatrix} \\ &=\begin{bmatrix}\mu^l {\partial {u}}/{\partial {y}} \\ -p \\ 0\end{bmatrix}\end{aligned}

So

\begin{aligned}\boldsymbol{\tau} \cdot (\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}})=\begin{bmatrix}1 & 0 & 0\end{bmatrix}\begin{bmatrix}\mu^l {\partial {u}}/{\partial {y}} \\ -p \\ 0\end{bmatrix}=\mu^l \frac{\partial {u}}{\partial {y}}\end{aligned} \hspace{\stretch{1}}(5.25)

Our boundary value condition is therefore

\begin{aligned}{\left.{{\mu^l \frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} ={\left.{{\mu^a \frac{\partial {u^a}}{\partial {y}}}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.26)

When can we decouple this, treating only the liquid? Observe that we have

\begin{aligned}{\left.{{\frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} ={\left.{{\frac{\mu^a}{\mu^l} \frac{\partial {u^a}}{\partial {y}}}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.27)

so if

\begin{aligned}\frac{\mu_a}{\mu_l} \ll 1\end{aligned} \hspace{\stretch{1}}(5.28)

we can treat only the liquid portion of the problem, with a boundary value condition

\begin{aligned}{\left.{{\frac{\partial {u^l}}{\partial {y}}}}\right\vert}_{{y = h}} = 0.\end{aligned} \hspace{\stretch{1}}(5.29)

Let’s look at the component of the traction vector in the direction of the normal (liquid pressure acting on the air)

\begin{aligned}\hat{\mathbf{n}} \cdot (\boldsymbol{\sigma}^l \cdot \hat{\mathbf{n}}) = \hat{\mathbf{n}} \cdot (\boldsymbol{\sigma}^a \cdot \hat{\mathbf{n}}) \end{aligned} \hspace{\stretch{1}}(5.30)

or

\begin{aligned}\begin{bmatrix}0 & 1 & 0\end{bmatrix}\begin{bmatrix}\mu^l \frac{\partial {u}}{\partial {y}} \\ -p^l \\ 0\end{bmatrix}= -{\left.{{p^l}}\right\vert}_{{y = h}} = -{\left.{{p^a}}\right\vert}_{{y = h}}\end{aligned} \hspace{\stretch{1}}(5.31)

i.e. We have pressure matching at the interface.

Our body force is

\begin{aligned}\mathbf{f} = \begin{bmatrix}g \sin\alpha \\ -g \cos\alpha \\ 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(5.32)

Referring to the Navier-Stokes equation 4.4, we see that our only surviving parts are

\begin{subequations}

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

\begin{aligned}0 = -\frac{\partial {p}}{\partial {y}} - \rho g \cos\alpha \end{aligned} \hspace{\stretch{1}}(5.33b)

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

\end{subequations}

The last gives us p \ne p(z). Integrating the second we have

\begin{aligned}p = \rho g y \cos\alpha + p_1\end{aligned} \hspace{\stretch{1}}(5.34)

Since p = p_{\text{atm}} at y = h, we have

\begin{aligned}p_{\text{atm}} = \rho g h \cos\alpha + p_1\end{aligned} \hspace{\stretch{1}}(5.35)

Our first NS equation 5.33a becomes

\begin{aligned}0 = \mu \frac{\partial^2 {{u}}}{\partial {{y}}^2} + g \sin\alpha,\end{aligned} \hspace{\stretch{1}}(5.36)

or

\begin{aligned}\frac{\partial^2 {{u}}}{\partial {{y}}^2} = -\frac{g}{\mu} \sin\alpha\end{aligned} \hspace{\stretch{1}}(5.37)

Solving we have

\begin{aligned}u = - \rho g \frac{\sin\alpha}{2 \mu} y^2 + A y + B\end{aligned} \hspace{\stretch{1}}(5.38)

With

\begin{aligned}u(0) &= 0 \\ {\left.{{\frac{\partial {u}}{\partial {y}}}}\right\vert}_{{y = h}} &= 0\end{aligned} \hspace{\stretch{1}}(5.39)

\begin{aligned}u = \rho g \frac{\sin\alpha}{2 \mu} \left( 2 h y - y^2 \right) .\end{aligned} \hspace{\stretch{1}}(5.41)

This velocity distribution is illustrated figure (\ref{fig:continuumL12:continuumL12fig7}).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL12fig7}
\caption{Velocity streamlines for flow down a plane.}
\end{figure}

It’s important to note that in these problems we have to derive our boundary value conditions! They are not given.

In this discussion, the height h was assumed to be constant, with the tangential direction constant and parallel to the surface that the liquid is flowing on. It’s claimed in class that this is actually a consequence of surface tension only! That’s not at all intuitive, but will be covered when we learn about “stability conditions”.

Study note.

Memorizing the NS equation is required for midterm, but more complex stuff (like cylindrical forms of the strain tensor if required) will be given.

References

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

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

 
Follow

Get every new post delivered to your Inbox.