Peeter Joot's Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘strain’

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 »

Continuum mechanics elasticity review.

Posted by peeterjoot on April 23, 2012

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

Motivation.

Review of key ideas and equations from the theory of elasticity portion of the class.

Strain Tensor

Identifying a point in a solid with coordinates x_i and the coordinates of that portion of the solid after displacement, we formed the difference as a measure of the displacement

\begin{aligned}u_i = x_i' - x_i.\end{aligned} \hspace{\stretch{1}}(2.1)

With du_i = {\partial {u_i}}/{\partial {x_j}} dx_j, we computed the difference in length (squared) for an element of the displaced solid and found

\begin{aligned}dx_k' dx_k' - dx_k dx_k = \left( \frac{\partial {u_j}}{\partial {x_i}} + \frac{\partial {u_i}}{\partial {x_j}} + \frac{\partial {u_k}}{\partial {x_i}} \frac{\partial {u_k}}{\partial {x_j}} \right) dx_i dx_j,\end{aligned} \hspace{\stretch{1}}(2.2)

or defining the \textit{strain tensor} e_{ij}, we have

\begin{aligned}(d\mathbf{x}')^2 - (d\mathbf{x})^2= 2 e_{ij} dx_i dx_j\end{aligned} \hspace{\stretch{1}}(2.3a)

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

In this course we use only the linear terms and write

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

Unresolved: Relating displacement and position by strain

In [1] it is pointed out that this strain tensor simply relates the displacement vector coordinates u_i to the coordinates at the point at which it is measured

\begin{aligned}u_i = e_{ij} x_j.\end{aligned} \hspace{\stretch{1}}(2.5)

When we get to fluid dynamics we perform a linear expansion of du_i and find something similar

\begin{aligned}dx_i' - dx_i = du_i = \frac{\partial {u_i}}{\partial {x_k}} dx_k = e_{ij} dx_k + \omega_{ij} dx_k\end{aligned} \hspace{\stretch{1}}(2.6)

where

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

Except for the antisymmetric term, note the structural similarity of 2.5 and 2.6. Why is it that we neglect the vorticity tensor in statics?

Diagonal strain representation.

In a basis for which the strain tensor is diagonal, it was pointed out that we can write our difference in squared displacement as (for k = 1, 2, 3, no summation convention)

\begin{aligned}(dx_k')^2 - (dx_k)^2 = 2 e_{kk} dx_k dx_k\end{aligned} \hspace{\stretch{1}}(2.8)

from which we can rearrange, take roots, and apply a first order Taylor expansion to find (again no summation convention)

\begin{aligned}dx_k' \approx (1 + e_{kk}) dx_k.\end{aligned} \hspace{\stretch{1}}(2.9)

An approximation of the displaced volume was then found in terms of the strain tensor trace (summation convention back again)

\begin{aligned}dV' \approx (1 + e_{kk}) dV,\end{aligned} \hspace{\stretch{1}}(2.10)

allowing us to identify this trace as a relative difference in displaced volume

\begin{aligned}e_{kk} \approx \frac{dV' - dV}{dV}.\end{aligned} \hspace{\stretch{1}}(2.11)

Strain in cylindrical coordinates.

Useful in many practice problems are the cylindrical coordinate representation of the strain tensor

\begin{aligned}2 e_{rr} &= \frac{\partial {u_r}}{\partial {r}}  \\ 2 e_{\phi\phi} &= \frac{1}{{r}} \frac{\partial {u_\phi}}{\partial {\phi}} +\frac{1}{{r}} u_r  \\ 2 e_{zz} &= \frac{\partial {u_z}}{\partial {z}}  \\ 2 e_{zr} &= \frac{\partial {u_r}}{\partial {z}} + \frac{\partial {u_z}}{\partial {r}} \\ 2 e_{r\phi} &= \frac{\partial {u_\phi}}{\partial {r}} - \frac{1}{{r}} u_\phi + \frac{1}{{r}} \frac{\partial {u_r}}{\partial {\phi}} \\ 2 e_{\phi z} &= \frac{\partial {u_\phi}}{\partial {z}} +\frac{1}{{r}} \frac{\partial {u_z}}{\partial {\phi}}.\end{aligned} \hspace{\stretch{1}}(2.12)

This can be found in [2]. It was not derived there or in class, but is not too hard, even using the second order methods we used for the Cartesian form of the tensor.

An easier way to do this derivation (and understand what the coordinates represent) follows from the relation found in section 6 of [3]

\begin{aligned}2 \mathbf{e}_i e_{ij} n_j = 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u} + \hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u}),\end{aligned} \hspace{\stretch{1}}(2.18)

where \hat{\mathbf{n}} is the normal to the surface at which we are measuring a force applied to the solid (our Cauchy tetrahedron).

The cylindrical tensor coordinates of 2.12 follow from
2.18 nicely taking \hat{\mathbf{n}} = \hat{\mathbf{r}}, \hat{\boldsymbol{\phi}}, \hat{\mathbf{z}} in turn.

Compatibility condition.

For a 2D strain tensor we found an interrelationship between the components of the strain tensor

\begin{aligned}2 \frac{\partial^2 e_{12}}{\partial x_1 \partial x_2} =\frac{\partial^2 {{e_{22}}}}{\partial {{x_1}}^2} +\frac{\partial^2 {{e_{11}}}}{\partial {{x_2}}^2},\end{aligned} \hspace{\stretch{1}}(2.19)

and called this the compatibility condition. It was claimed, but not demonstrated that this is what is required to ensure a deformation maintained a coherent solid geometry.

I wasn’t able to find any references to this compatibility condition in any of the texts I have, but found [4], [5], and [6]. It’s not terribly surprising to see Christoffel symbol and differential forms references on those pages, since one can imagine that we’d wish to look at the mappings of all the points in the object as it undergoes the transformation from the original to the deformed state.

Even with just three points in a plane, say \mathbf{a}, \mathbf{b}, \mathbf{c}, the general deformation of an object doesn’t seem like it’s the easiest thing to describe. We can imagine that these have trajectories in the deformation process \mathbf{a} = \mathbf{a}(\alpha, \mathbf{b} = \mathbf{b}(\beta), \mathbf{c} = \mathbf{c}(\gamma), with \mathbf{a}', \mathbf{b}', \mathbf{c}' at the end points of the trajectories. We’d want to look at displacement vectors \mathbf{u}_a, \mathbf{u}_b, \mathbf{u}_c along each of these trajectories, and then see how they must be related. Doing that carefully must result in this compatibility condition.

Stress tensor.

By sought and found a representation of the force per unit area acting on a body by expressing the components of that force as a set of divergence relations

\begin{aligned}f_i = \partial_k \sigma_{i k},\end{aligned} \hspace{\stretch{1}}(3.20)

and call the associated tensor \sigma_{ij} the \textit{stress}.

Unlike the strain, we don’t have any expectation that this tensor is symmetric, and identify the diagonal components (no sum) \sigma_{i i} as quantifying the amount of compressive or contractive force per unit area, whereas the cross terms of the stress tensor introduce shearing deformations in the solid.

With force balance arguments (the Cauchy tetrahedron) we found that the force per unit area on the solid, for a surface with unit normal pointing into the solid, was

\begin{aligned}\mathbf{t} = \mathbf{e}_i t_i = \mathbf{e}_i \sigma_{ij} n_j.\end{aligned} \hspace{\stretch{1}}(3.21)

Constitutive relation.

In the scope of this course we considered only Newtonian materials, those for which the stress and strain tensors are linearly related

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

and further restricted our attention to isotropic materials, which can be shown to have the form

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

where \lambda and \mu are the Lame parameters and \mu is called the shear modulus (and viscosity in the context of fluids).

By computing the trace of the stress \sigma_{ii} we can invert this to find

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

Uniform hydrostatic compression.

With only normal components of the stress (no shear), and the stress having the same value in all directions, we find

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

and identify this combination -3 \lambda - 2 \mu as the pressure, linearly relating the stress and strain tensors

\begin{aligned}\sigma_{ij} = -p e_{ij}.\end{aligned} \hspace{\stretch{1}}(3.26)

With e_{ii} = (dV' - dV)/dV = \Delta V/V, we formed the Bulk modulus K with the value

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

Uniaxial stress. Young’s modulus. Poisson’s ratio.

For the special case with only one non-zero stress component (we used \sigma_{11}) we were able to compute Young’s modulus E, the ratio between stress and strain in that direction

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

Just because only one component of the stress is non-zero, does not mean that we have no deformation in any other directions. Introducing Poisson’s ratio \nu in terms of the ratio of the strains relative to the strain in the direction of the force we write and then subsequently found

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

We were also able to find

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

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

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

Displacement propagation

It was argued that the equation relating the time evolution of a one of the vector displacement coordinates was given by

\begin{aligned}\rho \frac{\partial^2 {{u_i}}}{\partial {{t}}^2} = \frac{\partial {\sigma_{ij}}}{\partial {x_j}} + f_i,\end{aligned} \hspace{\stretch{1}}(4.35)

where the divergence term {\partial {\sigma_{ij}}}/{\partial {x_j}} is the internal force per unit volume on the object and f_i is the external force. Employing the constitutive relation we showed that this can be expanded as

\begin{aligned}\rho \frac{\partial^2 {{u_i}}}{\partial {{t}}^2} = (\lambda + \mu) \frac{\partial^2 u_k}{\partial x_i \partial x_k}+ \mu\frac{\partial^2 u_i}{\partial x_j^2},\end{aligned} \hspace{\stretch{1}}(4.36)

or in vector form

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

P-waves

Operating on 4.37 with the divergence operator, and writing \Theta = \boldsymbol{\nabla} \cdot \mathbf{u}, a quantity that was our relative change in volume in the diagonal strain basis, we were able to find this divergence obeys a wave equation

\begin{aligned}\frac{\partial^2 {{\Theta}}}{\partial {{t}}^2} = \frac{\lambda + 2 \mu}{\rho} \boldsymbol{\nabla}^2 \Theta.\end{aligned} \hspace{\stretch{1}}(4.38)

We called these P-waves.

S-waves

Similarly, operating on 4.37 with the curl operator, and writing \boldsymbol{\omega} = \boldsymbol{\nabla} \times \mathbf{u}, we were able to find this curl also obeys a wave equation

\begin{aligned}\rho \frac{\partial^2 {{\boldsymbol{\omega}}}}{\partial {{t}}^2} = \mu \boldsymbol{\nabla}^2 \boldsymbol{\omega}.\end{aligned} \hspace{\stretch{1}}(4.39)

These we called S-waves. We also noted that the (transverse) compression waves (P-waves) with speed C_T = \sqrt{\mu/\rho}, traveled faster than the (longitudinal) vorticity (S) waves with speed C_L = \sqrt{(\lambda + 2 \mu)/\rho} since \lambda > 0 and \mu > 0, and

\begin{aligned}\frac{C_L}{C_T} = \sqrt{\frac{ \lambda + 2 \mu}{\mu}} = \sqrt{ \frac{\lambda}{\mu} + 2}.\end{aligned} \hspace{\stretch{1}}(4.40)

Scalar and vector potential representation.

Assuming a vector displacement representation with gradient and curl components

\begin{aligned}\mathbf{u} = \boldsymbol{\nabla} \phi + \boldsymbol{\nabla} \times \mathbf{H},\end{aligned} \hspace{\stretch{1}}(4.41)

We found that the displacement time evolution equation split nicely into curl free and divergence free terms

\begin{aligned}\boldsymbol{\nabla}\left(\rho \frac{\partial^2 {{\phi}}}{\partial {{t}}^2} - (\lambda + 2\mu) \boldsymbol{\nabla}^2 \phi\right)+\boldsymbol{\nabla} \times\left(\rho \frac{\partial^2 {\mathbf{H}}}{\partial {{t}}^2} - \mu \boldsymbol{\nabla}^2 \mathbf{H}\right)= 0.\end{aligned} \hspace{\stretch{1}}(4.42)

When neglecting boundary value effects this could be written as a pair of independent equations

\begin{aligned}\rho \frac{\partial^2 {{\phi}}}{\partial {{t}}^2} - (\lambda + 2\mu) \boldsymbol{\nabla}^2 \phi = 0\end{aligned} \hspace{\stretch{1}}(4.43a)

\begin{aligned}\rho \frac{\partial^2 {\mathbf{H}}}{\partial {{t}}^2} - \mu \boldsymbol{\nabla}^2 \mathbf{H}= 0.\end{aligned} \hspace{\stretch{1}}(4.43b)

This are the irrotational (curl free) P-wave and solenoidal (divergence free) S-wave equations respectively.

Phasor description.

It was mentioned that we could assume a phasor representation for our potentials, writing

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

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

finding

\begin{aligned}\mathbf{u} = i \mathbf{k} \phi + i \mathbf{k} \times \mathbf{H}.\end{aligned} \hspace{\stretch{1}}(4.45)

We did nothing with neither the potential nor the phasor theory for solid displacement time evolution, and presumably won’t on the exam either.

Some wave types

Some time was spent on non-qualitative descriptions and review of descriptions for solutions to the time evolution equations we did not attempt

  1. P-waves [7]. Irrotational, non volume preserving body wave.
  2. S-waves [8]. Divergence free body wave. Shearing forces are present and volume is preserved (slower than S-waves)
  3. Rayleigh wave [9]. A surface wave that propagates near the surface of a body without penetrating into it.
  4. Love wave [10]. A polarized shear surface wave with the shear displacements moving perpendicular to the direction of propagation.

For reasons that aren’t clear both the midterm and last years final ask us to spew this sort of stuff (instead of actually trying to do something analytic associated with them).

References

[1] R.P. Feynman, R.B. Leighton, and M.L. Sands. Feynman lectures on physics.[Lectures on physics], chapter Elastic Materials. Addison-Wesley Publishing Company. Reading, Massachusetts, 1963.

[2] L.D. Landau, EM Lifshitz, JB Sykes, WH Reid, and E.H. Dill. Theory of Elasticity: Vol. 7 of Course of Theoretical Physics. 1960.

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

[4] Wikipedia. Compatibility (mechanics) — wikipedia, the free encyclopedia [online]. 2011. [Online; accessed 23-April-2012]. http://en.wikipedia.org/w/index.php?title=Compatibility_(mechanics)&oldid=463812965.

[5] Wikipedia. Infinitesimal strain theory — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 23-April-2012]. http://en.wikipedia.org/w/index.php?title=Infinitesimal_strain_theory&oldid=478640283.

[6] Wikipedia. Saint-venant’s compatibility condition — wikipedia, the free encyclopedia [online]. 2011. [Online; accessed 23-April-2012]. http://en.wikipedia.org/w/index.php?title=Saint-Venant\%27s_compatibility_condition&oldid=436103127.

[7] Wikipedia. P-wave — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 1-February-2012]. http://en.wikipedia.org/w/index.php?title=P-wave&oldid=474119033.

[8] Wikipedia. S-wave — wikipedia, the free encyclopedia [online]. 2011. [Online; accessed 1-February-2012]. http://en.wikipedia.org/w/index.php?title=S-wave&oldid=468110825.

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

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

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

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 »

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

Disclaimer.

This problem set is as yet ungraded.

Problem Q1.

Statement

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.

Solution

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)

and

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

or

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

Statement

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)

Find

\begin{enumerate}
\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?
\end{enumerate}

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

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumProblemSet1Q2fig1}
\caption{Q2. Characteristic equation.}
\end{figure}

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 ( http://www.wolfram.com/cdf-player/ ) )

For a Mathematica notebook that visualizes this part of this problem see https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2animated.cdf. 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.

Statement

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

Solution

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)

or

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

Statement

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.

Solution

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

\begin{itemize}
\item
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q1.cdf
\item
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2.cdf
\item
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q2animated.cdf
\item
https://raw.github.com/peeterjoot/physicsplay/master/notes/phy454/continuumProblemSet1Q3.cdf
\end{itemize}

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.

References

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

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

 
Follow

Get every new post delivered to your Inbox.