The Navier-Stokes equation derived in the previous post represents a springboard for the study of Newtonian fluids in a variety of scenarios. Because of the complexity of the equations of fluid mechanics, even this ‘approximation’ represents a degree of complexity that has yet to be completely understood. In fact, questions about the structure and existence of solutions to the Navier-Stokes equations are considered to be so important that the Clay Mathematics Institute is offering a Millennium Prize for proofs associated with this system as a step towards understanding turbulence.

As a result, the next set of posts will focus on subsets of these equations with the examples often being even more simplified. The hope is by looking at simple scenarios a degree of intuition may be developed. This approach is largely influenced by the fine book Elementary Fluid Mechanics by Acheson.

The first block of these will be concerned with the original jewel in the mathematical modeling of fluids in the 18th and 19th centuries – an ideal fluid.

An ideal fluid is a drastic simplification to the fluid flow model in which internal transports of energy or momentum between fluid elements are ignored. Thus the internal energy of a fluid element remains constant during the flow and its kinetic energy and momentum change only due to the external forces. In addition, the flow is assumed to be incompressible ($$\nabla \cdot \vec V = 0$$), meaning that the density of any fluid element remains constant during the flow and, as a further simplification, each fluid element will be assigned the same value of density. Note that these are two distinct assumptions about the density, a point Acheson goes to great lengths to emphasize in his discussion on pages 6, 23-4, and, in particular, 356, where he states quite clearly that:

Despite these seemingly drastic simplifications, there are numerous ordinary, terrestrial applications of the theory. The most common one would be the flow of water in a variety of ‘household’ circumstances. It should be noted that the assumption of incompressibility and a constant, homogeneous density is actually a better assumption than that the fluid has no viscosity (i.e. is inviscid). Fluid compressibility needs to be more carefully considered only when the flow speeds near the speed of sound in the fluid. Nonetheless, the ideal fluid was the jewel of yesteryear mostly because these assumptions actually meant one could solve a problem, even if the problem only partially matched what was physically observed. This gap between what can be modeled and what is observed will become a common theme in the posts to follow.

To get the equations for an ideal fluid, start with the Navier-Stokes equations

\[ \rho \frac{D \vec V}{D t} = – \nabla P + \mu^2 \nabla \vec V + \frac{1}{3} \mu \nabla (\nabla \cdot \vec V) + \rho \vec g \; , \]

(where the material derivative $$D/Dt = \partial_t + \vec V \cdot \nabla$$), and set the viscosity $$\mu$$ equal to zero, yielding Euler’s equations

\[ \frac{D \vec V}{dt} = -\frac{1}{\rho} \nabla P + \vec g \; .\]

In addition to these equations, we always have the continuity equation

\[ \frac{D \rho}{d t} + \rho \nabla \cdot \vec V = 0 \; ,\]

which, given the two incompressibility assumptions, simplifies to

\[ \nabla \cdot \vec V = 0 \; \]

and

\[ \rho = \mathrm{constant} \; .\]

Taken together, these are a set of 4 equations in the 4 unknowns $$\vec V$$ and $$P$$ (with $$\rho$$ assumed to be known for the particular situation).

Even in this simple form, the equations can be quite hard to solve as can be more readily seen if the equations are put into component form (after expanding the material derivative):

\[ \partial_t V_x + V_x \partial_x V_x + V_y \partial_y V_x + V_z \partial_z V_x = – \frac{1}{\rho} \partial_x P \; ,\]

\[ \partial_t V_y + V_x \partial_x V_y + V_y \partial_y V_y + V_z \partial_z V_y = – \frac{1}{\rho} \partial_y P \; ,\]

\[ \partial_t V_z + V_x \partial_x V_z + V_y \partial_y V_z + V_z \partial_z V_z = – \frac{1}{\rho} \partial_z P \; ,\]

and

\[ \partial_x V_x + \partial_y V_y + \partial_z V_z = 0 \; .\]

Obviously, these equations are not particularly illuminating, even for the great masters of the 18th and 19th centuries. It is remarkable then, that there is a very convenient simplification that was discovered by Daniel Bernoulli, called Bernoulli’s principle, which he published in 1738.

While Bernoulli didn’t have access to modern vector calculus (making his discovery all the more remarkable), we will not be too proud to employ it since it will make the derivation clearer and will serve as a springboard to subsequent discoveries.

Return to Euler’s equation

\[ \partial_t \vec V + (\vec V \cdot \nabla) \vec V = -\frac{1}{\rho} \nabla P + \vec g \; .\]

The most troublesome term is the nonlinear $$(\vec V \cdot \nabla) \vec V$$ term on the left-hand side. A convenient vector identity can be employed to separate this nonlinear term into solenoidal and irrotational pieces

\[ (\vec V \cdot \nabla) \vec V = (\nabla \times \vec V) \times \vec V + \nabla \left( \frac{1}{2} \vec V^2 \right) \].

We know define the vorticity vector (more about this in the next post) as

\[ \vec \omega = \nabla \times \vec V \; .\]

Euler’s equation now becomes (moving the irrotational term to the right-hand side)

\[ \partial_t \vec V + \vec \omega \times \vec V = -\nabla \left( \frac{1}{2} \vec V^2 \right) -\frac{1}{\rho} \nabla P + \vec g \; .\]

Since the gravitational force can also be derived from the gradient of a potential $$\vec g = -\nabla \chi$$, Euler’s equation can finally be written as

\[ \partial_t \vec V + \vec \omega \times \vec V = -\nabla H \; , \]

where

\[ H = \frac{P}{\rho} + \chi + \frac{1}{2} \vec V^2 \; \]

Now if the flow is steady, the time variation can be set to zero. Taking the dot product of both sides with $$\vec V$$ and noting that $$\vec V \cdot (\vec \omega \times \vec V) = \omega \cdot (\vec V \times \vec V) = 0$$ yields

\[ \vec V \cdot (\nabla H) = (\vec V \cdot \nabla) H = 0 \; .\]

The physical interpretation of this equation is that the value of H is conserved, for a given fluid element, as it flows. This, then, is Bernoulli’s principle which states that changes of a fluid element’s flow speed are balanced by corresponding changes in pressure or gravitational potential energy, such that the initial value of H remains constant.
Bernoulli’s principle precisely describes how high flow from a faucet, a garden hose, or fountain can reach when pushed upward from an orifice when it emerges with a particular pressure and speed. It represents the simplest form of energy conservation for a fluid.

A simple example of the Venturi tube will illustrate some of the applications. Consider the pipe shown below (assume cylindrical symmetry).

Let the inlet area and speed be $$A_1$$ and $$V_1$$, respectively, and the area and speed at the middle of the pipe be $$A_2$$ and $$V_2$$. The question we can answer is what is the difference between the inlet pressure and the pressure at the center. Following any stream line (thin blue lines), we know that the mass of the fluid element is conserved and so the total inlet mass must be conserved. The continuity equation takes the form

\[ \rho_1 A_1 V_1 = \rho_2 A_2 V_2 \; , \]

for mass per unit time.

Using Bernoulli’s equation we also get

\[ \frac{P_1}{\rho_1} + \frac{1}{2} V_1^2 = \frac{P_2}{\rho_2} + \frac{1}{2} V_2^2 \; .\]

Now, for an ideal fluid, the density is constant, and we can solve for $$V_2$$ in terms of the known areas and inlet speed. Substituting this result in Bernoulli’s gives a pressure difference of

\[ P_1 – P_2 = \frac{1}{2} \rho V_1^2 \left[ \left( \frac{A_1}{A_2} \right)^2 – 1 \right] \; \]

Clearly showing that the pressure at the center of the tube is lower than at the inlet with the difference depending on the inlet speed and the ratio of the two cross-sectional areas. This tube is used in the Venturi meter to determine/monitor flow speeds.

Next post will examine the role that the vorticity plays in ideal fluid flow.

This month’s column begins the process of narrowing down the general equations of fluid mechanics, developed over the last many columns, into a set of specific equations.  In particular, to make progress in describing specific situations, we must be able to relate the stress tensor to the state of the fluid.  This relation is analogous to the generalized Hooke’s law used in the study of elastic deformations.

While there are an amazingly large number of types of fluids out there in the world, this discussion and those that follow will be focused only on the simplest of these: Newtonian fluids for which there is a linear relationship between the applied shearing stress and the velocity gradients that result.  To further simplify matters, we shall also look only at the form of the Cauchy momentum equation for a Newtonian fluid, which is known by the famous name of the Navier-Stokes equation. The treatment here parallels the one found in An Introduction to Fluid Dynamics, by G.K. Batchelor.

There are interesting differences between a moving fluid and one at rest.  A fluid at rest necessarily exerts a normal stress independent of direction; no shear stresses are supported and the stress must to isotropic, leading to a very simple expression for the stress tensor

\[ T_{ij} = -p \, \delta_{ij} \; . \]

The parameter $$p$$ is the static fluid pressure, which can only be a function of position (e.g. the change in pressure with depth in a swimming pool).

Batchelor notes that there should be no expectation that this simple relation holds for a moving fluid, since, once in motion, the fluid can support shearing stress and will, in general, depart from isotropy.  Nonetheless, it helps to have a scalar quantity, similar to static pressure, that characterizes the “squeezing” of the fluid. This pressure is defined by averaging the normal component of stress over the surface of a small sphere centered on the position $$\vec x$$

\[ P = -\frac{1}{3} T_{ii} = – \frac{1}{3} T_{ij} \delta_{ij} = -\frac{1}{4 \pi} T_{ij} \int n_i n_j d \Omega(\vec n) \; , \]

where $$\vec n$$ is the normal to the surface of the sphere.  Batchelor emphasizes that $$P$$ is a mechanical definition with only a loose connection with thermodynamic notion of pressure, the latter requiring some notion of equilibrium, which is absent when the fluid is in motion.  Nonetheless, $$P$$ is a useful quantity because it is experimentally accessible from measurements of force and momentum (or, more accurately, from typical mechanical measurements that allow the inference of force and momentum).

Since $$P$$ is associated with the isotropic portion of the stress, it is mathematically convenient to decompose the stress tensor as

\[ T_{ij} = -P \delta_{ij} + d_{ij} \; \]

where $$d_{ij}$$ is a traceless tensor, known as the deviatoric stress, whose non-zero nature is due solely to the fluid motion.  It represents the internal transport of momentum, which can manifest itself as friction.

It is at this point, where the assumption that the fluid is Newtonian is brought to bear by expressing the deviatoric stress as being linearly proportional to the velocity gradient by

\[ d_{ij} = A_{ijk \ell} \frac{\partial V_k}{\partial x_\ell} \; \]

Where the tensor $$A_{ijk \ell}$$ depends on the local state of the fluid but not directly on the velocity distribution.  Since the deviatoric stress is symmetric in its indices we can conclude that $${\mathbf A}$$ is symmetric in its first two indices.

Next, we can decompose the velocity gradient as

\[ \frac{\partial V_k}{\partial x_{\ell}} = e_{k \ell} – \frac{1}{2} \epsilon_{k \ell m} \omega_m \; , \]

where $$e_{k \ell} = 1/2 (\partial_k V_\ell + \partial_\ell V_k)$$ is a symmetric tensor and $$\omega_m$$ is the vorticity defined, in terms of the anti-symmetric portion of the velocity gradient, by

\[ \omega_m = \epsilon_{m n p} \frac{1}{2} \left(\partial_n V_p – \partial_p V_n \right) \; .\]

The penultimate simplifying assumption in this long cavalcade of reductions, is to confine the discussion to gases and simple liquids, in which the ‘shape’ of the molecules is nearly spherical.  This assumption allows us to conclude that the deviatoric stress is independent of the orientation of the fluid element, thus imposing rotational symmetry on $$A_{i j k \ell}$$. This assumption should be viewed in contrast to a fluid whose molecular constituents are long chains and that, as a result, exhibits a preferred (i.e. easy) axis.

Imposition of rotational symmetry on $${\mathbf A}$$ means that it can only have the following general form

\[ A_{i j k \ell} = \mu \delta_{ik} \delta_{j \ell} + \mu’ \delta_{i\ell} \delta_{jk} + \mu^{\prime\prime} \delta_{ij}\delta_{kl} \; .\]

The ‘mu’ parameters describe the viscosity that arises from various shearing motions.  Because of the symmetries imposed on $${\mathbf d}$$ and upon $${\mathbf A}$$, they are not all independent.  Since $${\mathbf A}$$ is symmetric in its first two indices $$ \mu = \mu’$$ and

\[ A_{i j k \ell} = \mu ( \delta_{ik} \delta_{j \ell} + \delta_{i \ell} \delta_{jk} )+ \mu^{\prime\prime} \delta_{ij} \delta_{k\ell} \\ = \mu (\delta_{i \ell} \delta_{jk} + \delta_{ik} \delta_{j\ell} ) + \mu^{\prime\prime} \delta_{ij}\delta_{ \ell k} = A_{ij \ell k} \; ,\]

Clearly showing that $${\mathbf A}$$ is now symmetric in its third and fourth indices as well.  Note that this symmetry was made manifest by swapping the first and second terms and using the symmetry of the Kronecker delta in the term.

Armed with this knowledge of $${\mathbf A}$$, we return to the relation between the deviatoric stress and the velocity gradient to find

\[ d_{ij} = A_{ijk \ell} \frac{\partial u_k}{\partial x_{\ell}} = A_{ijk \ell} e_{k \ell} = \left[ \mu(\delta_{ik} \delta_{j\ell} + \delta_{i\ell} \delta_{jk}) + \mu^{\prime\prime} \delta_{ij} \delta_{k\ell} \right] e_{k\ell} \\ = 2 \mu e_{ij} + \mu^{\prime\prime} Tr({\mathbf e} \delta_{ij}) \; . \]

Note that the second equality follows from the symmetry of $${\mathbf A}$$ in its last two indices, the third equality by substituting in the specific form of $${\mathbf A}$$ and the final one by multiplying the contractions out and defining $$e_{kk} \equiv Tr({\mathbf e})$$.  The final step is to note that, since $${\mathbf d}$$ is traceless

\[ Tr({\mathbf d}) = 2 \mu Tr({\mathbf e}) + 3 \mu^{\prime\prime} Tr({\mathbf e}) = 0 \; \]

implying that $$\mu^{\prime\prime} = \, – 2/3 \mu$$.

After all this work, we now have a form of the deviatoric stress $${\mathbf d}$$ in terms of the symmetric part $${\mathbf e}$$ of the velocity gradient given by

\[ d_{ij} = 2 \mu \left( e_{ij} – \frac{1}{3} Tr({\mathbf e}) \delta_{ij} \right) \; . \]

This expression was the result of many hands in the mid-1800s, but often bears the name of Saint-Venant who derived essentially as was done above in 1843.

Suppose that the only non-zero velocity gradient component is $$\partial u_x/\partial y$$.  All components of deviatoric stress is zero except for the tangential stresses $$d_{xy} = d_{yx} = \mu \partial u_x/\partial y$$, which justifies identifying $$\mu$$ with the elementary definition of viscosity.

The Cauchy momentum equation (in conservative form) is now

\[ \frac{\partial \rho \vec V}{\partial t} + \nabla \cdot (\rho \vec V \vec V) = \nabla \cdot \left[ -p {\mathbf 1} + 2 \mu \left( e_{ij} – \frac{1}{3} Tr({\mathbf e}) \delta_{ij} \right) \right] + \vec f_b \; . \]

Finally, assume that the viscosity is a constant (that is we will initially confine ourselves to those cases where this is a very good approximation), to get the Navier-Stokes equations

\[ \frac{\partial \rho \vec V}{\partial t} + \nabla \cdot (\rho \vec V \vec V) = – \nabla p + \mu \nabla^2 \vec V + \frac{1}{3} \mu \nabla (\nabla \cdot \vec V) + \rho \vec g \; , \]

where we used $$Tr({\mathbf e}) = \nabla \cdot \vec V$$ and we specialized to gravity as the only body force: $$\vec f_b = \rho \vec g$$.  These cosmetic changes are design to facilitate comparison with the usual way the Navier-Stokes equations are written.

Next column will begin looking at information that can be interrogated from these equations and the solutions in simple situations.

This month’s column is a departure from the norm but with good reason.  The usual fare consists of the development and subsequent mathematical exploration of some model of a physical phenomenon.  However, this installment serves more as a survey covering a broad swath of possible fluid phenomena rather than a detail exposition of each type of fluid behavior.  The reason for this quite simple: the possible types of fluid behavior is far more enormous than most of us ever realize.  It’s so vast that it seems that the terminology and classification has not kept pace with the discoveries as evidence by the tangled naming convention and the wide-spread confusion that exists on the internet (and often within related Wikipedia articles).

As school children, each of us is taught that there are three phases of matter: solid, liquid, and gas.  We are further taught that while solids have definite volume and shape, liquids only have definite volume, and gases do not enjoy certainty in either property.

If one pursues a STEM career in college, then a bit more information is added to these explanations and one learns that fluids are the collective name that we give for substances that ‘flow’ under stress.  Sometimes one even finds out that the key distinction between solids and fluids is that solids deform from their natural shape to a new shape under the application of a stress (typically leaving the volume unchanged or nearly so) but then return to their natural shape when the load is removed.  The amount of the deformation is determined by the magnitude of the load and, for ideal elastic solids, the relationship is linear, meaning that a double of the stress results in a doubling of the elastic strain.  In contrast, fluids do not have a natural shape (i.e. independent of their container) and so when a stress is applied, the material begins to flow with no terminal shape set by the magnitude of the load.  Removing the application of the load doesn’t mean that the fluid will return to its initial shape but rather that the flow will continue until internal frictional forces bring the flow to a halt.  These frictional forces are termed viscous forces and the degree to which they act in a fluid determines the fluid’s viscosity.

Now none of the above is untrue, in so far as it goes, but it tends to leave the idea that there is a sharp distinction between phases of materials.  And it is that is exactly what common sense and daily experience tell us is not true.

The world is literally littered with common materials that sometimes behave like a solid and sometimes like a fluid.  Consider the oobleck craze that swept across the internet some years back, which is nicely presented in the following video.

;

The two main ingredients are powdered cornstarch and water.  By themselves, each of these materials flow under stress and neither has a definite shape even if they both have definite volume.  But combine them together and suddenly one has a material that sometime behaves as a liquid, when the applied stress is low, and other times like a solid, when the applied stress is high.

This behavior falls under the general heading of shear-thickening, the property whereby a fluid’s viscosity increases sharply as the applied stress increases.  An old school example of shear-thickening is Silly Putty, which is actually a viscoelastic material that behaves like an elastic solid on short time scales (it can holds its shape even while hammered and bounces nicely) but which flows slowly under low stress.

On the opposite side of the spectrum of flow behaviors is common household paint which flows far more easily when subject to the high stress as the brush or roller rub against the wall, allowing the paint to evenly wet the surface, but which becomes thicker under lower stress, allowing it to stay in place until it dries.  Many common kitchen items behave like this with ketchup

and whipped cream

being tasty examples.  Note how in the ketchup video, it flowed quite nicely when subjected to the high stress that resulted by pushing on the sides of the bottle but how it resisted gravity and crept at a much slower rate once it was in the glass.

Shear-thinning also can lead to some visually arresting behaviors like the Kaye effect in shampoo

Of course, not every material behaves like oobleck or Silly Putty or ketchup or shampoo.  Two of the most ubiquitous fluids, air and water, have constant viscosities (or nearly so) and their flow rate (i.e. strain rate) is proportional to the applied stress.

A general classification scheme exists that lumps all materials whose strain rate is proportional to the applied stress as Newtonian fluids and ‘all the others’ as non-Newtonian fluids.  And while the two most important fluids in our lives (at least external to our bodies) are Newtonian fluids, it seems pedagogically stunted to group the remaining fluids under the simple heading of non-Newtonian.  There are two broad categories under the non-Newtonian heading that were already discussed above:  shear-thickening and shear-thinning.  Graphically the differences between shear-constant (the name I argue should be attached to Newtonian fluids) and shear-thickening and -thinning are easily represented in the following diagram

The following table summarizes the various subcategories under shear-thinning and -thickening.

It is interesting to note that even this categorization is rough and the lines between the groupings blur.  For example, the time-dependency exhibited by rheopectic and thixotropic materials may simply be masking a thermodynamic dependence.  No one is really sure.  That is why the very common material of ketchup is listed in both shear-thinning sub-categories.  And as the discussion about ketchup on rheothing.com shows the lines between them are fuzzy.

Next month’s column will look at the simplest fluid model, the Newtonian fluid.  The advantage of this simplicity is that we can construct an explicit model of the behavior.  The disadvantage is that much of the interesting phenomena that make the world fun to play in will have to wait.

This post is the first in what should be a fairly long sequence of articles concerning fluid mechanics.  The aim of this series is to present an intuitive feel for the phenomena of fluid mechanics from a physical point-of-view and a clean way of understanding the mathematical models used to describe them.  This last point deserves some comments especially.

One need only consult the numerous textbooks that abound in the marketplace of ideas to walk away with a profound sense of bewilderment at all the possible ways to express the fundamental, governing equations of fluid mechanics.  One’s confusion further increases by both a lack of standardization of terminology and by the failure to keep a sharp line between essential principles, such as the material derivative or the conservation of energy, and those relationships tailored for a specific application, such as the connection between applied stress and the strain rate for a Newtonian fluid.

This post attempts to improve this situation by deriving the governing equations, with no simplifications or specializations, using a minimalist approach consisting of the following ingredients:

1. The definition of the material derivative to move between the particle-flow (Lagrangian) and field (Euler) points-of-view (see Continuum Mechanics 1 – Material Derivative for details)
2. The divergence theorem
3. The Reynolds transport theorem connecting the time rate-of-change of a volume of fluid to its flow velocity (see Continuum Mechanics 3 – Reynolds Transport Theorem for details)
4. Surface forces defined in terms of the stress tensor
5. The conservation of mass
6. Newton’s second law
7. The First Law of Thermodynamics

This presentation reflects the influence of many sources, but the primary ones are Mechanics by Symon, The Governing Equations of Fluid Mechanics by Anderson found in Chapter 2 of Computational Fluid Mechanics, edited by Wendt, and An Introduction to Fluid Mechanics by Batchelor.  All of these are fine resources but each has a particular drawback that this synthesis is designed to address.

The governing equations of fluid mechanics in the particle-flow picture are particularly simple.  In this picture, we imagine looking as a small fluid element surrounded by either other fluid elements or boundaries with other materials.  By tracking the interactions of this element with its environment, one can easily write down the governing equations in familiar terms.  Transformation to the field picture is then facilitated by applying a combination of the items 1-3 in the list above.

Broadly speaking, there are two ways to combine items 1-3 to arrive at the field picture: the ordinary and flux-conservative forms (with many variations in each).  As discussed in detail in Anderson, the flux-conservative form, written generally as

\[ \frac{\partial \alpha}{\partial t} + \nabla \cdot \vec j_{\alpha} = 0 \; , \]

for a physical quantity $$\alpha$$ and corresponding current $$\vec j_{\alpha}$$ is the preferred form as it leads to better behaved equations, especially when dealing with shocks.   The ‘fork in the road’ between these two ways will be obvious below and will center on either explicitly using mass conservation in the equations (ordinary) or ignoring it (flux-conservative).

Conservation of Mass = Continuity Equation

The conservation of mass in the particle-flow picture is simply

\[ \frac{d}{dt} m = 0 \; ,\]

which states that, no matter how the boundary of the element deforms as it moves, the amount of matter contained within stays constant.

Defining the particle mass in the usual way, $$m = \rho {\mathcal Vol}$$, in terms of the mass density and the volume, gives

\[ \frac{d}{dt} \left( \rho {\mathcal Vol} \right) = \frac{d \rho}{dt} {\mathcal Vol} + \rho \frac{d}{dt} {\mathcal Vol} \; . \]

Next using Reynold’s transport theorem to eliminate the time derivative of the volume gives

\[ \frac{d \rho}{dt} {\mathcal Vol} + \rho \nabla \cdot {\vec V} \, {\mathcal Vol} = 0 \; . \]

The next step consists of using the form of the material derivative in the field picture to arrive at

\[ \left( \frac{\partial \rho}{\partial t} + {\vec V}\cdot \nabla \rho \right) {\mathcal Vol} + \rho \nabla \cdot {\vec V} \, {\mathcal Vol} = 0 \; .\]

Combining the second and the third term and then dividing out the $${\mathcal Vol}$$ gives the flux-conservative form of the continuity equation

\[ \frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho {\vec V}\right) = 0 \; . \]

Newton’s Second Law = Cauchy Momentum Equation

Newton’s second law in the particle-flow picture looks identical to how it is framed in elementary classical mechanics.

\[ \frac{d}{dt} \left( m \vec V \right) = \vec F_{total}\; , \]

where the right-hand side is the sum over all the forces exerted on the element.

Break the right-hand side into surface forces due to exchanges of momentum only at the boundary of the fluid element and body forces (e.g. gravity) where momentum exchange takes place between a field and every part of the element:

\[ \vec F_{total} = \sum_{faces} {\vec f}_s + {\vec f}_b {\mathcal Vol} \; .\]

Expressing the surface force on a given face in terms of the stress tensor as $${\vec f}_s = {\mathbf T} \cdot {\hat n} \, {\mathcal Area} $$, simplifies the first term to

\[ \sum_{faces} {\vec f}_s = \sum_{faces} {\mathbf T} \cdot {\hat n} \, {\mathcal Area} \; . \]

Applying the divergence theorem gives the final form of the resultant of the surface forces

\[ \sum_{faces} {\vec f}_s = \nabla \cdot {\mathbf T} \, {\mathcal Vol} \; .\]

Now examine the left-hand side first, which we will call $${\dot {\vec p }}$$.

\[ {\dot {\vec p }} = \frac{d}{dt} \left(m {\vec V} \right) = \frac{d m}{dt} {\vec V} + m \frac{d}{dt} {\vec V} \; . \]

There is a natural temptation to set the first term to zero and then to follow the same strategy as above.  Suppose this is done then

\[ \rho \, {\mathcal Vol} \frac{d}{dt} {\vec V} = \nabla \cdot {\mathbf T} \, {\mathcal Vol} + {\vec f}_b {\mathcal Vol} \; .\]

Dividing out the $${\mathcal Vol}$$ and expanding the material derivative gives

\[ \rho \left( \frac{\partial \vec V}{\partial t} + {\vec V} \cdot \nabla {\vec V} \right) = \nabla \cdot {\mathbf T} + {\vec f}_b \; . \]

While this is a valid fluid equation, it is not in flux-conservative form and, as cited above, experience has shown that the best equations for computational modeling are in flux-conservative form.  So we retain the $$(dm/dt) {\vec V}$$ and perform similar steps as to the above to get the left-hand side into the form

\[ \frac{d}{dt} {\vec p} = \frac{dm}{dt} {\vec V} + m \frac{d {\vec V}}{dt} = \left( \frac{\partial \rho}{\partial t} + \nabla \cdot {\rho {\vec V}} \right) {\vec V} {\mathcal Vol} \\ + \rho \, {\mathcal Vol} \left( \frac{\partial {\vec V}}{\partial t} + {\vec V} \cdot \nabla {\vec V} \right) \; . \]

Factoring out $${\mathcal Vol}$$ and combining terms yields

\[ \frac{d}{dt} {\vec p} = {\mathcal Vol} \left( \frac{\partial (\rho {\vec V})}{\partial t} + \nabla \cdot ({\rho {\vec V}{\vec V}}) \right) \; . \]

Substituting in this form and dividing out the volume gives the Cauchy momentum equation in flux form

\[ \frac{\partial (\rho {\vec V})}{\partial t} + \nabla \cdot ({\rho {\vec V}{\vec V}}) = \nabla \cdot {\mathbf T} + {\vec f}_b \; .\]

First Law of Thermodynamics = The Energy Equation

The final equation is the first law of thermodynamics, which in the particle-flow picture looks like

\[ \frac{d}{dt} E = \frac{d {\mathcal W}}{dt} + \frac{ d{\mathcal Q}}{dt} \; , \]

where $$E$$ is the total energy of the element, $$W$$ is the work done to the fluid element, and $$Q$$ is the heat that flows into it.  As in treatment of the Cauchy momentum equation, we need to distinguish between work due to surface forces and work due to body forces.  We must also distinguish between heat produced within the fluid element and heat that flows in from the boundary.

Because there are so many possible variations, we will content ourselves with one of the possible forms.  Additional posts will explore other forms of this equation and additional assumptions.

To start, let’s look at the left-hand side.  The total energy is the sum element’s kinetic energy and its internal potential energy

\[ E = \frac{1}{2} m {\vec V}^2 + U \; .\]

Defining the total energy through the relation $$E = m e_t$$ with the specific total energy $$e_t = V^2/2 + u$$, we can express the left-hand side as

\[ \frac{d}{dt} E = \frac{d}{dt} \left( m e_t \right) = \frac{dm}{dt} e_t + m \frac{d e_t}{dt} \; .\]

As in the case of the momentum equation, keeping the time-derivative of mass explicitly around leads to the flux-conservative form.  This is done with the usual application of the Reynold’s transport theorem and the definition of the material derivative to arrive at

\[ \left[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec V) \right] e_t {\mathcal Vol} + \rho {\mathcal Vol} \left[ \frac{\partial e_t}{\partial t} + {\vec V} \cdot \nabla (e_t) \right] \; . \]

These two terms combine nicely to yield

\[ \left[ \frac{\partial \rho e_t}{\partial t} + \nabla \cdot (\rho e_t \vec V) \right] {\mathcal Vol} \; , \]

for the left-hand side.

Over on the right-hand side the time rate of work done is given by

\[ \frac{d W}{dt} = \vec F_{total} \cdot \vec V = \vec f_{body} \cdot \vec V \, {\mathcal Vol} + \sum_{faces} \vec f_{s} \cdot \vec V \; .\]

The last term can be put into a more useful form by again expressing the surface forces in terms of the stress tensor to get that the work done by these forces is then

\[ \sum_{faces} \vec f_s \cdot \vec V = \sum_{faces} ({\mathbf T} \cdot \vec V ) \cdot {\hat n} \, {\mathcal Area} \; .\]

Finally apply the divergence theorem to get

\[ \sum_{faces} \vec f_s \cdot \vec V = \nabla \cdot ({\mathbf T} \cdot \vec V ) {\mathcal Vol} \; ,\]

which, when substituted back, gives a final form for the rate of work performed of

\[ \frac{d W}{dt} = \left[ \vec f_{body} \cdot \vec V + \nabla \cdot ({\mathbf T} \cdot \vec V ) \right] {\mathcal Vol} \; .\]

The rate of change in the heat in the element is given by

\[ \frac{dQ}{dt} = \rho \dot q {\mathcal Vol} + \sum_{faces} \vec h \cdot \hat n {\mathcal Area} \; , \]

where the first term accounts for the volumetric heating occurring within the element and the second term represents the heat flux $$\vec h$$ entering through the faces.

Applying the divergence theorem to the last term gives a final form for the rate of change of the heat of

\[ \frac{dQ}{dt} = \left[ \rho \dot q + \nabla \cdot \vec h \right] {\mathcal Vol}\; . \]

Equating all of these pieces and dividing out the volume yields the general form of the energy equation in flux-conservative form

\[ \frac{\partial \rho e_t}{\partial t} + \nabla \cdot (\rho e_t \vec V) = \vec f_{body} \cdot \vec V + \nabla \cdot ({\mathbf T} \cdot \vec V ) + \rho \dot q + \nabla \cdot \vec h \; . \]

Next month’s columns will begin to look at specific physical solutions to this equations in simple settings.

This installment was intended to be a realigning of the recent series of articles on continuum mechanics to a set of columns specialized to fluid mechanics.  However, there is a nice derivation of the inner workings of the stress tensor that I developed while working on the new post on fluids that is worthy of some daylight.  In some sense, this post (and the one that follows) adds the rigor that was missing from the general discussion presented in the column entitled Continuum Mechanics 5 – The Stress Tensor.  So, this month’s offering is a complete rederivation of a coordinate-free representation of the stress tensor.

As discussed in that earlier post, the force on a surface is naturally a function of the surface normal.  What was taken on authority was that the mapping should be linear (i.e., tensorial) and symmetric in order that there be no unphysical translational and angular accelerations, respectively.  This post will prove that the mapping should be linear, and next month’s post will provide a rigorous derivation of the symmetry properties.

All that is needed to construct a mathematical derivation are some basic physical principles and elementary mathematics.  The starting point for the derivation is Cauchy’s tetrahedron,

a portion of the continuum where three sides of the figure (denoted by the superscripts $$(1), (2), (3)$$)  coincide with coordinate planes (i.e., $$x$$-$$y$$, $$y$$-$$z$$, and $$z$$-$$x$$ planes) and the fourth side (no superscript) closes the shape.

The physical principles that we will employ are Newton’s laws, the fact that the surface force across an interface is proportional to the area of the interface, and the fact that the volume scales as the cube of some characteristic length $$L$$ while the area scales as the square.

Newton’s second law applied to the tetrahedron yields

\[ m_{T} \, {\vec a}_{T} = {\vec f}^{(1)} + {\vec f}^{(2)} + {\vec f}^{(3)} + {\vec f} \equiv {\vec f}_{total} \; ,\]

where $$m_{T}$$ is the mass of the tetrahedron and $${\vec a}_{T}$$ is the acceleration. Since the mass of the tetrahedron is the product of the mass density times the volume, $$\rho {\mathcal V}_T$$, it scales as $$\alpha^3$$ as each linear dimension is scaled by a factor $$\alpha$$.  Each of the forces on the right-hand side is proportional to an area of the face (e.g., $${\vec f}^{(1)} \propto \delta S^{(1)}$$), it scales as $$\alpha^2$$ as each linear dimension is scaled by the same factor $$\alpha$$.

The entire equation, when subjected to this scaling, becomes

\[ m_{\mathcal T} \, \alpha^3 \, {\vec a}_{T} = \alpha^2 {\vec f}_{total} \; .\]

Dividing through by $$\alpha$$, gives the acceleration as

\[ {\vec a}_{T} = \frac{1}{\alpha \, m_{T}} {\vec f}_{total} \; .\]

In order to prevent an infinite acceleration as $$\alpha \rightarrow 0$$, we require that the forces exactly balance to zero.

\[ {\vec f}_{total} = {\vec f}^{(1)} + {\vec f}^{(2)} + {\vec f}^{(3)} + {\vec f} = 0 \; .\]

Now we want to characterize the forces.  Since we know that each force is proportional to the area and depends on the normal we can write the mapping generally as

\[ {\vec f}^{(i)} = {\mathbf T}\left({\hat n}^{(i)}\right) \delta S^{(i)} \; , \]

for $$i = 1, 2, 3$$ and ‘blank’ (that is, no superscript for the oblique face).

From Newton’s third law, we know that the force on one side of an interface must be equal and opposite to the force on the other side.  This relationship requires that

\[ {\mathbf T}({\hat n}) = – {\mathbf T}({\hat n}) \; . \]

By their orientation, the faces that coincide with a coordinate plane have unit normals that relate to the basis vectors as

\[ {\hat n}_1 = – {\hat e}_1 \; , \]

\[ {\hat n}_2 = – {\hat e}_2 \; , \]

and

\[ {\hat n}_3 = – {\hat e}_3 \; . \]

From basic geometry, the areas of each of these planes are related to the area of the fourth plane by

\[ \delta S_1 = {\hat n} \cdot {\hat e}_x \, \delta S \; , \]

\[ \delta S_2 = {\hat n} \cdot {\hat e}_y \, \delta S \; , \]

and

\[ \delta S_3 = {\hat n} \cdot {\hat e}_z \, \delta S \; . \]

Since the set $$\{{\hat e}_x, {\hat e}_x, {\hat e}_x,\}$$ spans the space, the normal to the oblique face can be decomposed as

\[ {\hat n} = ({\hat e}_x \cdot {\hat n}) {\hat e}_x + ({\hat e}_y \cdot {\hat n}) {\hat e}_y + ({\hat e}_z \cdot {\hat n}) {\hat e}_z \equiv a \, {\hat e}_x + b \, {\hat e}_y + c \, {\hat e}_z \; . \]

Returning to the force equilibrium, and expressing the forces in terms of the stress and areas, and the areas of the faces aligned with the coordinate planes in terms of the area of the oblique face, yields

\[ {\mathbf T}({\hat n}) \, \delta S = – {\mathbf T} \left( {\hat n}^{(1)} \right) \, ({\hat e}_x \cdot {\hat n} ) \, \delta S \\ – {\mathbf T} \left( {\hat n}^{(2)} \right) ({\hat e}_y \cdot {\hat n} ) \, \delta S \\ – {\mathbf T} \left({\hat n}^{(3)} \right) ({\hat e}_z \cdot {\hat n} ) \, \delta S \; . \]

Divide out the common factor of $$\delta S$$, replace the dot products of the normal with the basic vectors by the components, and exploit Newton’s third law to absorb the minus signs on the right-hand side

\[ {\mathbf T}({\hat n}) = a {\mathbf T} \left({\hat e}_{x}\right) + b {\mathbf T} \left({\hat e}_{y}\right) + c {\mathbf T} \left({\hat e}_{z}\right) \; . \]

Finally, substitute in the form of the normal into the left-hand side to get

\[ {\mathbf T}\left( a {\hat e}_x + b {\hat e}_y + c {\hat e}_z \right) = a {\mathbf T} \left({\hat e}_{x}\right) + b {\mathbf T} \left({\hat e}_{y}\right) + c {\mathbf T} \left({\hat e}_{z}\right) \; . \]

The foregoing result is the definition of a tensorial relationship, and the first half of the proof is complete.

This month’s post marks the end of the study of continuum mechanics, or rather it marks the end of the name.  The past nine posts have focused on the basic material common to elasticity and fluids and then moved into a brief study of elastic motion.  The material here will serve as the final post on elasticity, cleaning up some loose ends and providing some rigor.  Next month, the focus will turn to fluid mechanics and the series will be retitled and renumbered accordingly.

The three topics covered here provide some additional details and rigor to the discussions that have come before with a particular focus on the elastic equations of motion that are used to analyze the wave motion discussed in post 9.

Rigorous derivation of the strain tensor

When I was first exposed to elasticity, it was from Chapter 3 of Mathematical Methods, 3^rd edition by Arfken; and, in fact, I used that text as a guide for the material in post 7.  However, a brief consultation with more specialized texts will often yield a strain tensor that is a bit more complex than the one presented in either of those two sources.  A more complete derivation, which is patterned after the derivation here, supplemented with my own insights, is as follows.

Let the position of any point at any stage of the deformation be defined by $${\vec r}_t = \vec \phi (\vec r, t)$$,  where $$\vec r$$ is the initial position of the point and $$t$$ measures the extent of deformation.  The function $$\vec \phi(\vec r,t)$$ is called the flow mapping.  The parameter $$t$$ parameterizes the deformation, labeling the successive stages.  It doesn’t necessarily need to be regarded as the time but it is often convenient to treat it as such and it will used that way below.  The requirement that $$\vec r$$ is the original undeformed position imposes the obvious boundary: $$\vec \phi (\vec r, 0) = \vec r$$.

Now focus on two points $${\mathcal P}$$ and $${\mathcal Q}$$ before and after the derivation.   Their positions before the deformation are:

\[ \vec {\mathcal P}(0) = \vec \phi(\vec r,0) = \vec r \; \]

and

\[ \vec {\mathcal Q}(0) = \vec \phi(\vec r + d \vec r,0) = \vec r + d \vec r \; . \]

After the deformation, the same two points are now located at

\[ \vec {\mathcal P}(t) = \vec \phi(\vec r, t) \; \]

and

\[ \vec {\mathcal Q}(t) = \vec \phi(\vec r + d \vec r,t) = \vec \phi(\vec r,t) + \frac{\partial \vec \phi(\vec r, \epsilon) }{\partial \vec r} d \vec r \; , \]

where the assumptions that $$d \vec r$$ is small in magnitude and that $$\vec \phi$$ is differentiable have been made.

Now measure the distance between the points $${\mathcal P}$$ and $${\mathcal Q}$$.  Before the deformation, that distance squared is given by

\[ ds^2 = \left\Vert \vec {\mathcal P}(0) – \vec {\mathcal Q}(0) \right \Vert^2 = d \vec r \cdot d \vec r = d r_i d r_i \; . \]

After the deformation, the distance squared between those points, now moved to their new position, is given by

\[ ds_{t}^2 = \left\Vert \vec {\mathcal P}(t) – \vec {\mathcal Q}(t) \right \Vert^2 = \frac{\partial \phi_i}{\partial r_j} d r_j \frac{\partial \phi_i}{\partial r_k} d r_k \; . \]

Calculate the difference between these two lengths

\[ ds_{t}^2 – ds^2 = \frac{\partial \phi_i}{\partial r_j} d r_j \frac{\partial \phi_i}{\partial r_k} d r_k – d r_j d r_k \delta_{jk} = \left( \frac{\partial \phi_i}{\partial r_j} \frac{\partial \phi_i}{\partial r_k} – \delta_{jk} \right) dr_j dr_k\; ,\]

where the last term has been modified a bit by the inclusion of a Kronecker delta.

We now define the Green-Lagrange strain tensor by the relationship

\[ 2 \epsilon_{jk} = \frac{\partial \phi_i}{\partial r_j} \frac{\partial \phi_i}{\partial r_k} – \delta_{jk} \; .\]

The strain vector $$\vec u$$ measures the difference in position between where a particle originally was located and where it is located at time $$t$$

\[ \vec u(t) = \vec \phi(\vec r,t) – \vec r \; . \]

Spatial derivatives of $$\vec \phi$$ are related to spatial derivatives of the deformation vector $$u$$ via

\[ \frac{\partial \phi_i(\vec r,t)}{\partial r_j} = \delta_{ij} + \frac{\partial u_i}{\partial r_j} \; . \]

Substituting this expression into the Green-Lagrange strain tensor gives

\[ 2 \epsilon_{jk} = \left( \delta_{ij} + \frac{\partial u_i}{\partial r_j} \right) \left( \delta_{ik} + \frac{\partial u_i}{\partial r_k} \right) – \delta_{jk} \; . \]

Expanding and simplifying yields

\[ \epsilon_{jk} = \frac{1}{2} \left( \frac{\partial u_k}{\partial r_j} + \frac{\partial u_j}{\partial r_k} + \frac{\partial u_i}{\partial r_j}\frac{\partial u_i}{\partial r_k} \right) \; , \]

which is the same as the earlier derivation except for the nonlinear term that can often be dropped.

Rigorous derivation of the equations of motion

The same notions can be used to derive the equations of motion that led to the wave equation analyzed in the last post.  Consider a cube of the continuum centered at $$\vec r$$.

The force on any face is given by $$f_i = T_{ij} n_j Area$$.  Start first with the forces on the two faces perpendicular to the $$x$$-axis.  The center of the rightmost face is located at $$\vec r + (\delta x/2) \hat \imath$$ and it is at this point that the stress is evaluated.  Likewise, the center of the leftmost face lies at $$\vec r – (\delta x/2) \hat \imath$$.  For notational convenience, functions evaluated on the rightmost face will be denoted with the argument $$x_{+}$$ and those evaluated on the leftmost face with $$x_{-}$$.  In this notation, the components of the forces on these faces are:

\[ \left[ \begin{array}{c} f_x \\ f_y \\ f_z \end{array} \right] = \left[ \begin{array}{c} T_{xx}(x_+) – T_{xx}(x_-) \\ T_{yx}(x_+) – T_{yx}(x_-) \\ T_{zx}(x_+) – T_{zx}(x_-) \end{array} \right] dy dz = \left[ \begin{array}{c} T_{xx,x}(\vec r) \\ T_{yx,x}(\vec r) \\ T_{zx,x}(\vec r) \end{array} \right] dx dy dz \; , \]

where the last equalities result by taking a Taylor’s expansion of the finite difference and keeping terms of only first order ($$T_{ij,k} = \partial T_{ij} / \partial r_k$$).

Similar manipulations give the expressions

\[ \left[ \begin{array}{c} f_x \\ f_y \\ f_z \end{array} \right] = \left[ \begin{array}{c} T_{xy}(y_{+}) – T_{xy}(y_-) \\ T_{yy}(y_+) – T_{yy}(y_-) \\ T_{zy}(y_+) – T_{zy}(y_-) \end{array} \right] dx dz = \left[ \begin{array}{c} T_{xy,y}(\vec r) \\ T_{yy,y}(\vec r) \\ T_{zy,y}(\vec r) \end{array} \right] dx dy dz \; \]

and

\[ \left[ \begin{array}{c} f_x \\ f_y \\ f_z \end{array} \right] = \left[ \begin{array}{c} T_{xz}(z_+) – T_{xz}(z_-) \\ T_{yz}(z_+) – T_{yz}(z_-) \\ T_{zz}(z_+) – T_{zz}(z_-) \end{array} \right] dx dy = \left[ \begin{array}{c} T_{xz,z}(\vec r) \\ T_{yz,z}(\vec r) \\ T_{zz,z}(\vec r) \end{array} \right] dx dy dz \; , \]

for the components of the force generated on the $$y$$- and $$z$$-faces, respectively.

Combining these results and expressing the sums in index notation gives the final, compact expression for the force as

\[ f_i = T_{ij,j} dx dy dz \; .\]

Note that the mass of the continuum contained in the cube is $$\rho \, dx dy dz$$, where $$\rho$$ is the mass density.  Newton’s second law for the cube then relates the acceleration of the cube, which is expressed in terms of the second time-derivative of the flow mapping to the cube’s mass and the forces acting on the cube as

\[ (\rho \, dx dy dz) \frac{\partial^2 \phi(\vec r,t)}{\partial t^2} = f_i \; . \]

Using the relation between the flow mapping and the deformation vector and between the force and the stress tensor gives

\[ \rho \, \partial^2_t u_i = T_{ij,j} \; , \]

which is precisely the starting expression for the last post.

Relation between the Lame constants and the elastic moduli

This final section derives the relationships to express the elastic moduli derived in post 7 to the Lame constants, which were introduced as a shorthand for the wave equation in an isotropic medium.  The Lame constants took the form

\[ \lambda = \frac{\nu E}{(1+\nu)(1-2\nu)} \; \]

and

\[ \mu \frac{E}{2(1+\nu)} \; ,\]

where $$E$$ is Young’s modulus and $$\nu$$ is Poisson’s ratio.

It is a matter of some simple algebra to express Young’s modulus and Poisson’s ration in terms of $$\lambda$$ and $$\mu$$ as follows.  Manipulating the second equation gives  $$E = 2 (1+\nu) \mu$$.  Substituting this into the first equation gives $$\lambda = 2 \nu \mu / (1-2\nu)$$ from which one gets via simplification

\[ \nu = \frac{\lambda}{2 (\mu + \lambda) } \; . \]

Back-substituting gives

\[ E = \mu \frac{2\mu + 3 \lambda}{\mu + \lambda} \; \]

for Young’s modulus.

To relate the shear and bulk moduli to the Lame constants, return to the generalized Hooke’s law

\[ T_{ij} = 2 \mu \epsilon_{ij} + \lambda Tr(\epsilon) \delta_{ij} \; . \]

In the case of an isolated shear along the $$x$$-direction that grows along the $$y$$-direction, the strain is expressed as $$\partial_y u_x$$ and the only non-zero terms in the strain tensor are $$ \epsilon_{xy} = \epsilon_{yx} = 1/2 \, \partial_y u_x$$.  Thus,

\[ T_{xy} = 2 \mu \epsilon_{xy} = 2 \mu \frac{1}{2} \partial_y u_x = \mu \partial_y u_x \; .\]

Since the shear modulus is defined (see post 8) as the ratio of stress $$T_{xy}$$ and the strain $$\partial_y u_x$$, it is clear that $$G = \mu$$.

Finally, the bulk modulus applies when the continuum is in isotropic (i.e. hydrostatic) stress equilibrium.  Isotropic stress means that the strain tensor is diagonal with all its on-diagonal components equal: $$T_{xx} = T_{yy} = T_{zz} \equiv T_{iso}$$.  The strain tensor is also diagonal with $$\epsilon_{xx} = \epsilon_{yy} = \epsilon_{zz} \equiv \epsilon_{iso}$$.  The application of $$T_{iso}$$ represents a pressure that uniformly compresses the continuum.  The measure of that compression is given by the fractional change in volume.  Consider a small cube of initial length $$L$$ on a side so that applying $$T_{iso}$$ causes an isotropic strain of $$\Delta L$$ in all directions.  The volume before the strain is $$L^3$$ and the volume after is, to first order, $$L^3 – 3 L^2 \Delta L$$.  The fractional change in the volume is then

\[ \frac{\Delta V}{V} = -\frac{3 L^2 \Delta L}{L^3} = -3 \frac{\Delta L}{L} = -3 \epsilon_{iso} \; . \]

The bulk modulus is defined as the ratio

\[ K = -\frac{\Delta P}{\Delta V/V} = \frac{T_{iso}}{3 \epsilon_{iso}} \; .\]

From the generalize Hooke’s law

\[ T_{iso} = 2 \mu \epsilon_{iso} + \lambda 3 \epsilon_{iso} = (2 \mu + 3 \lambda) \epsilon_{iso} \; . \]

Comparing the two expressions gives

\[ K = \lambda + \frac{2}{3} \mu \; . \]

So, that’s it for now on elastic solids and that corner of continuum mechanics.  Next month, I turn towards the related subdiscipline of fluid mechanics.

Last month’s post set the stage for one of the more interesting and important aspects of elastic motion – the propagation of elastic deformation waves in a material.  The familiar travelling and standing waves seen on stringed instruments are the prototype example but the richness and complexity is increased when the allowed deformations take place in more dimensions.  The complex relationship of the generalized Hooke’s law relating the strain to the stress tensor given by

\[ T_{ij} = c_{ijkl} \epsilon_{kl} \; ,\]

leads to many more propagation modes compared to the one-dimensional textbook case.  Since every solid material can support elastic waves, understanding these modes becomes important in the mechanical design of all manner of objects, such as bridges, support columns, cars, planes, skyscrapers, and so on.  This analysis also serves as a cornerstone in seismology for understanding one of the most powerful and catastrophic events on the planet, the earthquake.

As discussed in the last post, there are 81 components in the stiffness tensor $$c_{ijkl}$$ but there are never more than 21 independent components due to various symmetries of the physical relationships connecting the stiffness tensor to the stress and strain tensors.  This is a substantial reduction but is still far too complicated for an initial foray into elastic waves so the first step is to look at the waves that arise in an isotropic medium.

To perform this analysis requires three steps.  The first step simplifies the generalized Hooke’s law in terms of the elastic moduli discussed in detail in the previous post.   The second step adapts Newton’s law for the elastic medium yielding the wave equation.  The last step is the decomposition of the wave equation into its independent components.

Step 1 – Simplifying Generalized Hooke’s Law

To simplify the generalized Hooke’s law for an isotropic material begin by focusing in on the strains that result along the $$x_1$$-axis as a result of the application of various stresses.  Applying a small tensile stress along the $$x_1$$-axis gives

\[ E \epsilon_{11} = T_{11} \; , \]

where $$E$$ is Young’s modulus.

Application of the same small tensile stress along the $$x_2$$- and $$x_3$$-axes yields

\[ -\frac{E}{\nu} \epsilon_{11} = T_{22} \; \]

and

\[ -\frac{E}{\nu} \epsilon_{11} = T_{33} \; , \]

where $$\nu$$ is Poisson’s ratio.  Note that no indices are needed on this Poisson’s ratio since the material is isotropic.

Putting these pieces together gives

\[ E \epsilon_{11} = (1+\nu) T_{11} – \nu \left(T_{11} + T_{22} + T_{33} \right) = (1+\nu) T_{11} + \nu Tr({\mathbf T}) \; ,\]

where $$Tr({\mathbf T})$$ is the trace of the stress tensor.

Arfken presents a rigorous argument for generalizing the above equation but it is obvious that it is consistent with the form

\[ E \, \epsilon_{ij} = (1+\nu) T_{ij} – \nu Tr({\mathbf T}) \delta_{ij} \; . \]

Taking the trace of both sides and performing some straightforward algebra gives

\[ T_{ij} = \left( \frac{E}{1+\nu} \right) \epsilon_{ij} + \left(\frac{\nu}{1+\nu}\right) \left(\frac{E}{1-2\nu}\right) Tr({\mathbf \epsilon}) \delta_{ij} \; . \]

It is traditional to recast this equation slightly by defining the Lame coefficients

\[ \lambda = \frac{\nu E}{ (1+\nu)(1-2\nu) } \; ,\]

and

\[ \mu = \frac{E}{2 (1+\nu)} \; .\]

The isotropic stress-strain relationship now takes on the compact form

\[T_{ij} = 2 \mu \epsilon_{ij} + \lambda Tr({\mathbf \epsilon}) \delta_{ij} \; , \]

which is the desired relation.

Step 2 – Adapting Newton’s Law

The next step is to adapt Newton’s law to a piece of the continuum.  The standard way of doing this is to look at a small cube of the continuum whose sides are spanned by the components $$u_i$$, similar defined as in the earlier post in which the strain tensor is first derived.  The resulting equation is then (see https://www.geophysik.uni-muenchen.de/~igel/Lectures/Sedi/sedi_weq.pdf)

\[ \rho \partial_t^2 u_i = T_{ij,j} \; .\]

Substituting in for $$T_{ij}$$ with the relationship from step 1 yields

\[ \rho \partial_t^2 u_i = \partial_j \left( 2 \mu \epsilon_{ij} + \lambda Tr({\mathbf \epsilon}) \delta_{ij} \right) \; .\]

By definition of the strain tensor

\[2 \epsilon_{ij} = \partial_i u_j + \partial_j u_i \; \]

and

\[ Tr({\mathbf \epsilon}) = \partial_k u_k \; . \]

Substituting these definitions in yields

\[ \rho \partial_t^2 u_i = \lambda \partial_i \partial_k u_k + \mu \partial_i \partial_j u_j + \mu \partial^2_j u_i \; .\]

In vector terms

\[ \rho \partial_t^2 {\vec u} = (\lambda + \mu) \nabla (\nabla \cdot \vec u) + \mu \nabla^2 \vec u \; . \]

Step 3 – Decomposing the Wave Equation

As it stands the, wave equation mixes the various independent modes.  This observation isn’t at all obvious but it will become apparent with the application of the vector identity

\[ \nabla^2 \vec u = \nabla(\nabla \cdot \vec u) – \nabla \times (\nabla \times \vec u) \; .\]

Substituting this relation into the wave equation and collecting terms gives a new form of the wave equation

\[ \rho \partial_t^2 \vec u = (\lambda + 2 \mu) \nabla (\nabla \cdot \vec u) – \mu \nabla \times (\nabla \times \vec u) \; .\]

The importance of this form is that there are now two uncoupled modes: the ‘divergence’ mode and the ‘curl’ mode.  To see this, take the divergence of both sides and exchange all partial derivatives with abandon.  Since the divergence of a curl is identically zero the result is

\[ \rho \partial_t^2 \nabla \cdot \vec u = (\lambda + 2 \mu) \nabla^2 (\nabla \cdot \vec u) \; .\]

From this equation, we can see that a compression/expansion wave propagates with a speed given by

\[ v_{P} = \sqrt{ \frac{\lambda + 2 \mu}{\rho} } \; . \]

The ‘P’ label stands for primary or pressure wave.

Likewise, taking the curl gives

\[ \rho \partial_t^2 \nabla \times \vec u = \mu \nabla^2 (\nabla \times \vec u) \; .\]

From this equation, we can see that a shear wave propagates with a speed given by

\[ v_{S} = \sqrt{ \frac{\mu}{\rho} } \; . \]

The ‘S’ label stands for secondary or shear wave.

As the links to Wikipedia show, this simple theory has been used to analyze earthquake propagation.

The last post introduced the idea of the deformation field $${\vec u}(\vec r)$$ and the corresponding strain tensor $$2\epsilon_{ij} = \partial_j u_i + \partial_i u_j$$, which characterizes the amount that any two points in the continuum move away from each other.

While formally correct, the mathematical form of the strain tensor tends to obscure any physical intuition as to what different components are all about.  Some small hint for how to interpret the strain tensor was touched upon in the closing of the last post.  This post provides a more in-depth coverage by introducing the elastic moduli.

The elastic moduli generalize the concept of a spring concept from elementary mechanics.  Roughly speaking, these numbers describe the resistance that the material has to one of three basic actions:  1) Young’s modulus (usually called $$E$$) – being stretched or compressed along a particular direction, 2) bulk modulus (usually called $$K$$) – being expanded or contracted in all directions, and 3) shear modulus (usually called $$G$$) – being sheared.  The primary distinction between 1) and 2) is that the density of the material need not change during a uniaxial deformation but must necessarily change for a full expansion or contraction.  Detail for each of these follows.

Youngs Modulus

Consider a uniform rod of material of a rectangular cross-section $$A$$ with a natural length (i.e. stress-free) $$L$$.  When subjected to a tension force $$F$$, the rod expands by $$\Delta L$$.

Measurements of the extent of the stretch as a function of the stress gives Young’s modulus as

\[ E = \frac{F/A}{\Delta L/L} = \frac{T_{xx}}{\epsilon_{xx}} \; , \]

where the last equality follows by defining the long axis of the bar as the $$x$$-axis and $$T_{xx}$$ and $$\epsilon_{xx}$$ are the corresponding components of the stress and strain tensors.

Interestingly, the story of uniaxial stress isn’t quite complete.  As anyone who has every play with rubber bands knows, when an elastic material lengthens in one direction it typically shrinks in the transverse directions; the band gets thinner in the two directions perpendicular to the ‘stretch’ direction.

Generally, the transverse directions do not shrink by the same amount since the degree of contraction depends on material composition and this difference is manifest in $$\Delta L’ \neq \Delta L”$$.  Poisson’s ratios are defined as

\[ \nu_{zx} = -\frac{\Delta L’}{\Delta L} \; \]

and

\[ \nu_{zy} = -\frac{\Delta L”}{\Delta L} \; ,\]

where the notation $$\nu_{ij}$$ measures the ratio of change along the $$j$$-axis as a function of the strain caused along the $$i$$-axis by the corresponding stress.

When $$\Delta L’ = \Delta L”$$ the material has an axis of symmetry about which there is an axial uniformity.  If Young’s modulus is the same in all direction, then the material is called isotropic and no indices are needed as there is only one Poisson’s ratio needed.  The assumption of isotropy works for many simple materials but certain elastic man-made materials have very directionally-dependent Poisson ratios.   The Wikipedia article on Poisson’s ratio lists six separate ratios $$\nu_{ij}$$ for the fabric Nomex.

Note that as the material undergoes its elongation and retraction in the direction of the stress and in the perpendicular directions, respectively, that the density can remain constant or that it change.  The exact determination depends on the material as does just how much the two perpendicular directions differ in the amount they contract.

Bulk Modulus

Consider a small unit of material suddenly placed inside a high-pressure environment.  An example would be the submergence of a metal cube at the bottom of the ocean.  Compression stress assaults the cube from all sides seeking to change the volume and, hence, the density of the material.

The bulk modulus measures how resistant the material is to this change and, defined as the ratio of the pressure $$P$$ to the fractional change in the volume $$\Delta V / V$$, becomes

\[ K = – \frac{\Delta P}{\Delta V/V} \approx -V \frac{d P}{d V} ; .\]

Shear Modulus

Consider a rod of material that is placed between two plates, the bottom one fixed and the top one free to move along the $$x$$-direction.  As it moves, the top plate exerts a stress $$T_{xy}$$ on the top face of area $$A$$, which deforms the material by varying degrees as a function of distance $$y$$ from the fixed plate.

If the height of the rod is $$L$$ and the amount of deformation at the top is $$\Delta x$$ then the shear modulus is defined as

\[ G_{xy} = \frac{F/A}{\Delta x/L} = \frac{T_{xy}}{\partial_y u_{x}} = \frac{T_{xy}}{2 \epsilon_{xy}} \; . \]

Similar expressions exist for the other possible combinations of directions.  An isotropic material needs only one shear modulus to describe it.

General Relation of Strain to Stress

The definitions of the three basic moduli relate strain to stress in an intuitive and geometric way but the tensor nature of the stress and strain are absent and any connection between the moduli are obscure at this point.  The final step is to synthesize them into a general relationship between the stress tensor $$T_{ij}$$ and the strain relationship $$2 \epsilon_{ij}$$ and to discover how the moduli connect to each other.

The full tensor relationship is given by the generalize Hooke’s law

\[ T_{ij} = c_{ijkl} \epsilon_{kl} \; , \]

where the fourth-rank stiffness tensor $$c_{ijkl}$$ plays the role of the spring constant in the simply form of Hooke’s law $$F = -k x$$.   There are 81 different components of the stiffness tensor but most of them are not independent.  The symmetric nature of the stress and strain tensors ($$T_{ij} = T_{ji}$$ and $$\epsilon_{ij} = \epsilon_{ji}$$) imposes the relationships

\[ c_{ijkl} = c_{jikl} = c_{ijlk} = c_{jilk} \]

reducing the number of independent components to 36.  A further reduction down to 21 independent components comes from noting that the stiffness matrix is also symmetric upon the interchange of the first and second pair of indices

\[ c_{ijkl} = c_{klij} \; .\]

This symmetry in the stiffness tensor is inherited from the potential energy for an elastic body

\[ PE = \frac{1}{2} c_{ijkl} \epsilon_{ij} \epsilon_{kl} \; , \]

which is clearly invariant under the exchange of $$\epsilon_{ij}$$ and $$\epsilon_{kl}$$.

Next month’s column focuses on isotropic materials, exploring the specific relationships between the moduli in this case, and by looking at how elastic waves propagate.

As part of the series on Continuum Mechanics, this post is the first of several focusing on elastic behavior of materials. Referring to the previous post, elastic materials are ones which, in contrast with fluids, deform under the application of a stress but return to their equilibrium shape upon the removal of the stress.

Before characterizing the strain mathematically, consider first some physical expectations from the simple bending of a beam under the load. To focus the discussion, look at the three specific points in orange, yellow, and black.

Prior to the deformation, the separation between the orange and yellow points is entirely along the horizontal, and the separation between the yellow and black points is entirely along the vertical. Denote the relative position of the orange point relative to the yellow as $$\delta {\vec r}_0$$ and then examine what happens to that relative position under a load. When the beam bends, the material along the top edge stretches while the material along the bottom compresses. This is seen by the fact that the new relative position, $$\delta {\vec r}_1$$, from orange to yellow is now shorter in length and has components in both directions (its original length and direction being shown in red). Likewise, the relative position between yellow and black changes its orientation substantially but negligibly changes its length since both points are along the edge and the bending stress barely compresses the material at the edge.

The deformation of the relative position from the orange and yellow points measures how much change $$\delta {\vec r}_0$$ suffers under the load of the bending stress. Care needs to be exercised when examining the change in orientation since a rigid body rotation does not constitute deformation even though it rotates the entire body. Thus, the measure of deformation should look something like

\[ deformation = \delta {\vec r}_1 – \delta {\vec r}_2 – rigid\;\;body\;\;rotations \; ,\]

with the additional understanding that this deformation will vary from point to point within the material, that is, the deformation is a field.

The mathematical characterization of deformation presented here closely follows Mathematical Methods for Physicists by Arfken, Third Edition.

Like the bending beam discussed above, Arfken focuses on two points at $${\mathcal P}_0$$ and $${\mathcal Q}_0$$.

Assume $${\mathcal P}_0$$ is located at $${\vec r}$$ and that the relative position of $${\mathcal Q}_0$$ is $$\delta {\vec r}$$. Under the applied stress, $${\mathcal P}_0$$ transports by $${\vec u}(\vec r)$$ to its new location at $${\mathcal P}_1$$. Likewise, $${\mathcal Q}_0$$ transports by $${\vec u}(\vec r + \delta {\vec r})$$ to its new location $$\delta {\vec r}’ = {\vec u} + \delta {\vec u}$$ relative to $${\mathcal P}_1$$. Because it describes the transport of each point in the material, $${\vec u}({\vec r})$$ is called the transport field.

Since a rigid-body motion would have transported both $${\mathcal P}_0$$ and $${\mathcal Q}_0$$ by a common $${\vec u}$$, independent of position, the deformation is clearly related to any transport difference, $$\delta {\vec u}$$, between the two.

The transport difference

\[ \delta {\vec u}({\vec r}) = {\vec u}({\vec r}+\delta {\vec r}) – {\vec u}({\vec r}) \; \]

expands, to first order, as

\[ \delta {\vec u} = \nabla_{\vec r} {\vec u} \cdot {\vec r} \; .\]

The gradient of the transport field is a second-rank tensor, called the transport tensor, that maps an initial relative position to a final transport difference. In components, its representation is

\[ \nabla {\vec u} = \frac{\partial u_i}{\partial r_j} \delta r_j \equiv \partial_j u_i \delta r_j \; . \]

It is important to realize that $$\nabla {\vec u}$$ is not yet the strain tensor being hunted. There is a possibility that there is a rotation hidden within.

The usual way to deal with this is to decompose the transport tensor into symmetric and anti-symmetric pieces defined through

\[ \delta u_i = \eta_{ij} \delta r_j + \xi_{ij} \delta r_j \; , \]

where

\[ \eta_{ij} = \frac{1}{2} \left( \partial_j u_i + \partial_i u_j \right) \, \]

and

\[ \xi_{ij} = \frac{1}{2} \left( \partial_j u_i – \partial_i u_j \right) \, . \]

Arfken notes that the anti-symmetric piece, which is related to the curl of the transport field, contains the rigid-body rotations. To see this, consider the following expression

\[ [ (\nabla \times {\vec u}) \times \delta {\vec r} ]_i = \epsilon_{ijk} \epsilon_{j\ell m} \partial_{\ell} u_m \delta r_k \; , \]

where $$[{\vec r}]_j \equiv r_j$$ is the jth component of the vector $${\vec r}$$.

Performing the usual index gymnastics yields
\[ [ (\nabla \times {\vec u}) \times \delta {\vec r} ]_i = \partial_j u_i \delta r_j – \partial_i u_j \delta r_j = 2 \xi_{ij} \delta r_j \; . \]

So the anti-symmetric portion of the transport difference represents a rotation of $${\mathcal Q}$$ about $${\mathcal P}$$ along the $$\nabla \times {\vec u}$$ axis by $$|\nabla \times {\vec u}|$$ radians.

The symmetric piece of the transport difference, $$\eta_{ij}$$, represents the actual strain that develops at $${\vec r}$$.

The diagonal components of $$\eta_{ij}$$ represent stretches or compressions and the off-diagonal components correspond to shears. To see these assertions consider the original orange and yellow points separated from each other by one unit of distance along the horizontal, which, for convenience, will be called the x-axis: ($$\delta {\vec x} = {\hat \imath} $$).

After the deformation, the transport difference is given by

\[ \delta {\vec u} = \eta_{11} {\hat \imath} + \eta_{12} {\hat \jmath} + \eta_{13} {\hat k} \; .\]

The position of yellow relative to orange is given by

\[ \delta {\vec y} = \delta {\vec x} + \delta {\vec u} = (1 + \eta_{11} ) {\hat \imath} + \eta_{12} {\hat \jmath} + \eta_{13} {\hat k} \; .\]

Projected into two dimensions, the before and after configurations look as below.

Note that the change in the relative position has components along both the $$x$$- and $$y$$-axes, with the later growing larger as the distance along the $$x$$-axes from the orange point grows. Since the length of the relative position changed, there is a scaling involved. Also, the growth along the $$x$$-axis of the vertical component represents a shear in contrast to a rotation.

With the basic notions of the strain tensor $$\eta_{ij}$$ well established, the next post will cover the mapping of this tensor to the various moduli of elasticity.

Up to this point in these posts on continuum mechanics, everything has been generally applicable.  The material derivative, the Reynolds and Flux transport theorems, and the nature of the stress tensor are applicable to any type of continuum.  Unfortunately, the party needs to end at some point and this post can be considered as the ‘last call’ that signals that end.

To understand why the generalities must stop, we need to discuss something mathematical and something physical.

On the mathematical side, we ultimately need to remind ourselves of a fact that is easy to lose sight of amidst all this machinery.  Namely, we are seeking to solve for the behavior of the continuum as a function of time.  We are seeking to generalize Newton’s laws to the continuum.

In their conventional setting, Newton’s laws amount to the following statement:

\[ m \frac{d^2 {\bar q}(t)}{d t^2} = {\bar f}({\bar q},{\dot {\bar q}},t) \; , \]

where $$\dot {\bar q}$$ stands for the time-derivative of $${\bar q}(t)$$ as a function of time.  The key point here is that the force on the right-hand side, $$\bar f$$ can be determined by knowing the current state itself (state being defined as the phase-space ordered pair of $$({\bar q},{\bar{\dot q}})$$ and time).  If this relationship didn’t exist, then literally there would be no way to evolve the system forward in time.  Some entity somewhere would need to supply the rule for how to compute the force at each step.

When generalized to a continuum setting, these underlying mathematical requires still hold sway but in a modified form.  To give some sense how this modified form might look, consider the wave equation,

\[ \frac{\partial^2 A(x,t)}{\partial t^2} = v^2 \frac{\partial^2 A(x,t)}{\partial x^2} \;, \]

as it is usually presented for the amplitude of a vibrating string $$A(x,t)$$.  This equation is usually derived in the continuum limit by looking at a small portion of the string and analyzing the tensions pulling on each side (see, e.g., Mechanics by Symon).  The right-hand side of the wave equation is the force on each infinitesimal length of string located at $$x$$.  To see this, note that another path to the wave equation uses a discrete, lumped-coupled mass-spring system.  This system is then transitioned into a continuum using the appropriate limit as the masses and the spacing between them both go to zero such that ratio of mass to spacing length goes to the mass density.  The article on the wave equation in Wikipedia, has an accessible presentation of this latter derivation using discrete masses and springs obeying Hooke’s law.

The point that is again worth emphasizing is that the force on the continuum is known in terms of the continuum itself.  Here the state of the system is characterized in terms of the spatial variation of $$A(x,t)$$ through its second partial derivative with respect to $$x$$.  This characterization looks odd at first glance but when expressed in terms of finite differences,

\[ \frac{\partial^2 A(x,t)}{\partial x^2} \approx \frac{A(x+h,t) – 2 A(x,t) + A(x-h,t)}{h^2} \; \]

it merely expresses the ‘force’ in terms of the state of the string at points $$x-h, x, x+h$$ and should be considered as no more odd than expressing the force as the combination

\[ \frac{\vec r}{\sqrt{x^2 + y^2 + z^2}} \; \]

as is done for central force laws.

Obviously, what we need to find for continuum mechanics is a relationship that specifies the stress in terms of some function of the state of the continuum.  We should prepare ourselves for the fact that the function will almost always depend on some spatial and, perhaps, temporal derivative.  We actually got a foretaste of that in the final result of the Reynolds transport theorem

\[ \frac{d}{dt} Vol(t) = Vol(t) div(\vec V) \; , \]

which relates the kinematic description of a small volume to the spatial variation of the flow velocity.

The final step is to sketch out the physical expectations for how that relationship might look functionally for various materials.  It is here that the fork in the road comes for which individual materials behave remarkedly different.

In the case of a simple elastic material, the body deforms with an application of a force but returns to its equilibrium shape with the removal of the load.  The picture looks something like this:

In continuum mechanics, the deformation of the material is described by comparing the initial configuration of the body (usually called the reference configuration) to the configuration after application of the load (usually called the current configuration).  The strain of the body (i.e. the response) finds expression in the relative movement of a set of points in the reference configuration to the same set in the current configuration.  Removing rigid-body motion leaves what is called the strain.  The outcome of this operation is to provide a direct functional connection $$\sigma = f(\epsilon)$$ between the observed strain, $$\epsilon$$, and the applied stress, $$\sigma$$, which is what is needed to have mathematical closure.   This relationship generalizes the usual Hooke’s law encountered in basic mechanics, which states that the force due to a spring is given by $$\vec F = -k \vec x$$, but the generalization, although still linear, is a bit more complicated due to tensorial nature of the stress and strain.

In contrast to elastic materials, whose deformation only persists as long as the stress is applied, fluids (liquids and gases) flow under the application of a stress.  Their deformations continue to change after the load is removed.

In these cases, the key observation comes from Newtonians laws.  The stress being an outside force, can only accelerate the material during its application.  The key relationship that relates the stress to an observed property of the body is of the form $$\sigma = f(\dot \epsilon)$$, which says that by observing the rate-of-change of the strain one can deduce the stress-applied.

To be complete, the two relationships detailed above only scratch the surface of the various types of materials and the resulting functional relationships between strain in a body and the stress applied.  There are hosts of other behaviors and often the situation is muddled due to the fact that the functional relationship is determined ‘the other way around’.  New materials are subjected to a variety of known stresses in order to observe the corresponding strains.  This approach results in an implicit relationship between strain and stress.  But regardless of whether the relationship is explicit or implicit, the mathematical closure of the theory requires that the force on the material be expressible as a function of the state of the material and, perhaps time, just as it is done for ordinary particle mechanics.

The coming posts will explore some of these relationships explicitly but no attempt to be complete and comprehensive will be made – nature is just too varied to believe that can be done.