Buoyancy and Stability: An Introduction

Ever since people set out to sea in ships, the issues of buoyancy and stability have been of importance. In spite of this, the treatment it receives in textbooks is often lacking. Following is an overview of the subject; basic understanding of the principles is essential in performing the experiment and interpreting the results.


Buoyancy is ultimately what makes things float, such as the buoy in Figure 1. This is true whether the material the boat is made of is lighter than water (like the balsa wood rafts Thor Heyerdahl and his crew crossed the Pacific with in 1947) or heavier than water. The latter would include objects from the buoy shown to the ships of the U.S. Navy.

Figure 1
Figure 1

The basic concept is very simple: for anything placed in a fluid medium, the upward force the medium exerts on the body is equal to the weight of the fluid the body displaces. This is not only true of bodies placed in water; it is also true of those in air. The difference is that, for those in air, the weight of the air displaced is usually not enough to “float” the aircraft. A notable exception are dirigibles such as the “Goodyear blimp,” which is filled with helium, a gas lighter than air. Another lighter-than-air gas used is hydrogen. This is very combustible, as everyone was reminded of when the Hindenburg caught fire in New Jersey in 1938.

Most buoyancy applications are marine ones, and it is those we will concentrate on in this experiment. We will also concentrate on rectangular forms and flat-bottomed vessels, which simplifies the math somewhat. However, these principles can be extended to just about any floating craft.

Using a flat-bottomed craft also makes it easier to understand why displacing a fluid works. Consider first the following: how the force of the fluid on the flat hull of a craft varies with depth1:

Figure 2: Illustrating Water Pressure Increasing in Proportion to Draught

For a fluid at rest, the hydrostatic pressure increases linearly with depth, thus


Figure 3 Inadequate Freeboard
Figure 3 Inadequate Freeboard

where p is the hydrostatic pressure, γ is the unit weight of the water, and D is the depth from the water’s surface to the bottommost point of the vessel, usually called the draught. This distance from the water line to the top of the rectangle (the gunwale) is called the freeboard; the results of inadequate freeboard can be seen in Figure 3.

In any case, for a vessel of beam (width) W and a length L the volume it displaces is given by the equation


Combining and rearranging these two equations,


For the boat to float, it has to be in static equilibrium, and so the downward force of the weight of the boat Wboat must equal the upward force Fbuoyant. Therefore,


So we’ve established a relationship between the weight of the boat and the volume of water it displaces. The “far right” hand side only applies to boats with a flat bottom and straight sides.

What this means is that there are three ways we can weigh an existing boat:

  1. We can simply weigh it on a scale. For small boats this isn’t too difficult; larger ones can be tricky.  We can then estimate how far it will sink into the water.

  2. We can measure the freeboard, then obtain D and, knowing L, W and the unit weight of water, we can compute the weight of the boat. This works easily for rectangular boats; for real boats, you have to determine the relationship between the actual waterline and the displacement, then see where the actual waterline ends up.

  3. We can use an overflow method, which is okay for small experiments (like Archimedes used) but not so hot on a larger scale.  But this illustrates our concept.

Procedure for determining volume of water displacement2:

Figure 4 Displacement


Buoyancy is a fairly straightforward concept, although it may be a little hard to grasp up front. Stability—the ability of the ship to resist overturning—is a little more difficult, although it’s obviously important, as the following diagram of a ship with waves coming at the beam shows3.

Figure 5: Showing the Transverse Movements of an Easy-Rolling Vessel Among Waves, and Also of a Raft

Let’s define (or recall) a couple of terms.

Centre of Gravity: this is easy, mathematically this is the centroid of the mass or weight of the ship. An illustration of this is below.

Figure 6: Centre of Gravity of a Ship

Centre of Buoyancy: this is a little trickier, this is the centroid of the cross-sectional area of the ship under the water line, as shown below.

Figure 7: Centre of Buoyancy of a Box-Shaped Vessel

As you can see, for a box-shaped vessel which is not listing (i.e., leaning at an angle) or has no squat (i.e., not angled along the length of the boat) the centre of buoyancy is located halfway down the draught of the vessel, halfway across the beam, and dead amidships.

The centre of gravity and the centre of buoyancy are not necessarily at the same place; in fact, they are usually different. That difference determines both the stability of the ship and, literally, how it rolls.

We know that motor vehicles with high centres of gravity (such as off-road vehicles) are more prone to turn over in use than those with lower centres of gravity. Ships are the same; we need to have a way to decide how stable a ship is and whether there is a point that a ship becomes unconditionally stable or unconditionally unstable.

As long as a ship is upright, and both the centre of gravity and the centre of buoyancy are in the centre of the ship in all respects, it is theoretically possible for a ship never to turn over. As a practical matter this is impossible; even very large ships like cruise ships, which use their size to resist roll in most wave situations, are going to roll some. Below is a diagram which shows the centre of gravity and the centre of buoyancy for a ship which is upright and which is inclined 14º.

Figure 8: Vessel Floating Upright, and Inclined 14 Degrees

We need to look at this carefully and note the following:

  • The point G is the centre of gravity of the ship.

  • The point B or B’ is the centre of buoyancy of the ship. In the course of inclination the centre of buoyancy will change because the shape of the cross-section under the waterline changes; this is fairly simple to calculate for rectangular ships and more complicated for curved hull shapes.

  • The point M is the metastatic point of the ship. The distance GM is called the metastatic height of the ship.

  • If point G is below point B or B’, the ship is unconditionally stable; it will not turn over unless G and B’ is changed by taking on water, shifting cargo in the ship, etc.

  • If point G is below point M, the ship is conditionally stable, and if point G is above point M, the ship is unconditionally unstable.

The reason for this last point is simple: the ship above is rolling in a clockwise direction. The resisting moment of the buoyancy, calculated by (GZ)(Wbuoyant) is counter-clockwise, as the buoyant force is upward. This is true as long as G is below M. If G moves upward above M, then the now driving moment (GZ)(Wbuoyant) turns clockwise, the same direction as the rolling of the ship, and the ship will generally turn over4.

Thus the location of M, abstract as it may seem, becomes a critical part of the design of a ship. But how is it done? There are two methods we will discuss here of determining the metastatic height of a ship.

Determining Metastatic Height

Theoretical Method

This method uses the following formula to determine the location of the metacentre:


For a rectangular vessel, the moment of inertia is the same as we used in mechanics of materials, i.e., LW3/12, and is applied as follows:

The displacement volume was given earlier. We then compute the distance between the metacentre M and the centre of buoyancy B as follows:

Note carefully that this is NOT the metacentric height GM; it is then necessary to subtract the distance from the centre of buoyancy to the centre of gravity from this result to obtain GM. This is done as follows:

Figure 9: Computing GM From the Height of the Metacentre Above the Centre of Buoyancy

It’s worth noting here that the location of point M is independent of the centre of gravity and dependent upon the geometry of the ship and its volume under the water line (or total weight.)

Timing the Roll

This method is sort of an “old salt’s” rule of thumb method. First, let’s define the roll time. The roll time is the time it takes for a ship to start from rest at an angle of roll (port or starboard,) roll to the opposite side, and return to the original orientation. This can be approximated by the equation5


where tr = roll time of ship, seconds
GM = metastatic height of ship, meters or feet
W = beam of ship, meters or feet
C = constant based on units of GM and B
= 0.44 for units of feet
= 0.80 for units of meters

Solving for metastatic height,


This is significant for another reason: another rule of thumb used in yacht design is that the roll time is seconds should be between 1 and 1.1 times the beam W of the boat in meters6. Yachts with shorter roll times tend to “check” or quickly come back to centre when in rough seas; this can be a hard experience for passengers and make a complete mess of stowed cargo. Yachts with longer roll times will come over like they’re about to capsize, and then slowly roll back around. The usual result of this is seasickness and a miserable ride.

1Walton, T. (1899) Know Your Own Ship. Charles Griffin and Company, London, England. Much of the material that follows on buoyancy and stability comes from this work.

2Keep in mind that the unit weight of sea water is greater than fresh. Why is this so?

3Seamanship pointer: if you’re in a boat and are facing high waves, wake, etc., best way to take them is to point the blow into the direction the waves are coming from, not take the waves on the beam.

4Whether it actually does gets into issues of freeboard, rate of roll, etc., which are beyond this presentation or experiment.

5Nudelman, N. (1992) Yacht Design Course. Lesson 6: Stability. Westlawn Institute of Marine Technology, Stamford, CT.

6Gerr, D. (1992) The Nature of Boats. International Marine Publishing, Camden, ME. Much of the material on this method of determining GM is taken from this source.

When Catholic Academia Bails on Philosophy, We’re All in Trouble

Which is what some of it, at least, has done:

A similar crisis has shaken the philosophical estate within the church. Before 1970 philosophy enjoyed an enviable prominence in the curriculum of Catholic colleges. This Neo-Scholastic philosophy was certainly structured around the perennial questions—Does God exist? What is virtue?—but it was an odd, manual Thomism in which students never actually read Aquinas. A smug catechetical certitude seemed to lurk behind the paint-by-numbers proofs and the gleeful one-paragraph refutations of modern “adversaries.”

That world has disappeared; its chastened replacement in the Catholic academy bears the stamp of marginality: minimal curricular presence, hyper-specialization, incoherence among the squabbling philosophical factions.

This cultural recession of philosophy has encouraged some Catholics to abandon philosophy as a central component of the church’s discourse. The issue has become especially neuralgic in the dispute over the formation of clergy. But the project of a nonphilosophical Catholicism is fraught with peril.

In some ways, the greatest blessing God bestowed on me in my Christian formation was to dodge education at one of these institutions, which enabled me to take in Aquinas and the like on the side.  Not only did it avoid the problems there, but inculcation in philosophically structured Christianity has helped me to avoid some of the sillier–and more dangerous–trends in Charismatic and Pentecostal thought.

Catholic theology in particular is pretty much toast outside of a philosophical framework.  And that throws away one of the major advantages that Catholicism has.  That’s a pity, the rest of us need the discipline.

The blunt truth is that Evangelical theology is an oxymoron precisely because it rejects any philosophical framework.  The Bible, however, was written in the flow of human history and experience, where we live, and was intended to address that experience directly.  Sooner or later, however, the question of “why?” will come up, and without a philosophical framework that question is unanswerable.

Today Pentecostal and Charismatic theology is at a crossroads because it inherited Evangelical theology’s basic thought structure without its limiting assumptions.  The result is that some in the academy (and elsewhere) are about to take the leap outside of Christianity without knowing it, and we all know where that ends.

There are problems with philosophy too; I tackled that issue some time back here.  Some of the problems we have now could be avoided by jettisoning much of modern philosophy altogether.  But throwing out the baby with the bath water isn’t the answer, for Catholics or anyone else.

Getting Closer on the St. Andrew’s Sex Scandal?

It’s not the closest thing, but…

Bishop Audrey Scanlan of the Episcopal Diocese of Central Pennsylvania yesterday removed the Rev. Howard White from the priesthood.

White, 75, was among several adults who sexually abused students at St. George’s School in Middleton, Rhode Island in the 1970s and 80s, according to a report released recently by independent investigators on behalf of the school.

In the wake of media reports about sexual abuse at St. George’s, several people from other Episcopal dioceses in which White had worked said that he sexually abused them when they were young.

It’s interesting to note that, according to this, “(Headmaster) Zane terminated White in 1974 after he said White admitted ‘sexually abusing a sophomore boy and attempting to sexually abuse at least two and likely three others.'”  The wheels of justice turn very slowly in this “enlightened” Episcopal Church, since someone was put on notice 42 years ago.  It puts Episcopal Bishop Porter Taylor’s statement that White “…had been identified by former students of St. George’s School in Rhode Island as having engaged in sexual misconduct in the early 1970s while he served on the staff at that school” in a different light.  He should have used the plaintiff attorney’s favourite phrase “by his own admission,” but he didn’t.

So what does this have to do with the scandals at the other end of the East Coast?  As I noted in my last piece on the subject, the link is the Rev. George Andrews II, who went from St. George’s to St. Andrew’s as Headmaster in the late 1980’s.  Will he and other Episcopal reverends be implicated?  We will see.

And as far as that “other” sex scandal we have in the Presidential election, I’ll stand by my previous position that, not so far in the future, such things won’t be scandalous any more, not in our society.  As was the case with the Roman Catholic Church, the left will milk the scandal cow as long as it lives before they butcher it, and Evangelical leadership is just plain clueless.

Mohr’s Circle Analysis Using Linear Algebra and Numerical Methods


Mohr’s Circle-or more generally the stress equilibrium in solids-is a well known method to analyse the stress state of a two- or three-dimensional solid. Most engineers are exposed to its derivation using ordinary algebra, especially as it relates to the determination of principal stresses and invariants for comparison with failure criteria. In this presentation, an approach using linear algebra will be set forth, which is interesting not only from the standpoint of the stress state but from linear algebra considerations, as the problem makes an excellent illustration of linear algebra concepts from a real geometric system. A numerical implementation of the solution using Danilevsky’s Method will be detailed.

1. Introduction

The analysis of the stress state of a solid at a given point is a basic part of mechanics of materials.Although such analysis is generally associated with the theory of elasticity, in fact it is based on static equilibrium, and is also valid for the plastic region as well. In this way it is used in non-linear finite element analysis, among other applications. The usual objective of such an analysis is to determine the principal stresses at a point, which in turn can be compared to failure criteria to determine local deformation.In this approach the governing equations will be cast in a linear algebra form and the problem solved in this way, as opposed to other types of solutions given in various textbooks. Doing it in this way can have three results:

  1. It allows the abstract concepts of linear algebra to be illustrated well in a physical problem.
  2. It allows the physics of the determination of principal stresses to be seen in a different way with the mathematics.
  3. It opens up the problem to numerical solution, as opposed to the complicated closed-form solutions usually encountered, of the invariants, principal stresses or direction cosines.

2. Two-Dimensional Analysis

2.1. Eigenvalues, Eigenvectors and Principal Stresses

The simplest way to illustrate this is to use two-dimensional analysis. Even with this simplest case, the algebra can become very difficult very quickly, and the concepts themselves obscured. Consider first the stress state shown in Figure 1, with the notation which will be used in the rest of the article.


Figure 1: Stress State Coordinates (modified from Verruijt and van Bars [8])

The theory behind this figure is described in many mechanics of materials textbooks; for this case the presentation of Jaeger and Cook [4] was used. The element is in static equilibrium along both axes. In order for the summation of moments to be zero,

\tau_{xy}=\tau_{yx} (1)

The angles are done a little differently than usual; this is to allow an easier transition when three-dimensional analysis is considered. The direction cosines based on these angles are as follows:

l=cos\alpha\ (2)

m=cos\beta\ (3)

Now consider the components p_{x},\,p_{y} of the stress vector p on their respective axes. Putting these into matrix form, they are computed as follows:

\left[\begin{array}{cc}  \sigma_{{x}} & \tau_{{\it xy}}\\  \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}}  \end{array}\right]\left[\begin{array}{c}  l\\  \noalign{\medskip}m  \end{array}\right]=\left[\begin{array}{c}  p_{{x}}\\  \noalign{\medskip}p_{{y}}  \end{array}\right] (4)

which is a classic Ax = b type of problem. At this point there are some things about the matrix in Equation 4 that need to be noted (DeRusso et al. [2]):

  1. It is square.
  2. It is real and symmetric. Because of these two properties:
    1. The eigenvalues (and thus the principal stresses, as will be shown) are real. Since for two-dimensional space the equations for the principal stresses are quadratic, this is not a “given” from pure algebra alone.
    2. The eigenvectors form an orthogonal set, which is important in the diagonalization process.
  3. The sum of the diagonal entries of the matrix, referred to as the trace, is equal to the sum of the values of the eigenvalues. As will be seen, this means that, as we rotate the coordinate axes, the trace remains invariant.

At this point there are two related questions that need to be asked. The first is whether static equilibrium will hold if the coordinate axes are rotated. The obvious answer is “yes,” otherwise there would be no real static equilibrium. The stresses will obviously change in the process of rotation. These values can be found using a rotation matrix and multiplying the original matrix as follows (Strang [7]):

\left[\begin{array}{cc}  cos\alpha & -sin\alpha\\  sin\alpha & cos\alpha  \end{array}\right]\left[\begin{array}{cc}  \sigma_{{x}} & \tau_{{\it xy}}\\  \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}}  \end{array}\right]=\left[\begin{array}{cc}  \sigma'_{{x}} & \tau'_{{\it xy}}\\  \noalign{\medskip}\tau'_{{\it xy}} & \sigma'_{{y}}  \end{array}\right] (5)

The rotation matrix is normally associated with Wallace Givens, who taught at the University of Tennessee at Knoxville. The primed values represent the stresses in the rotated coordinate system. We can rewrite the rotation matrix as follows, to correspond with the notation given above:^

\left[\begin{array}{cc}  l & -m\\  m & l  \end{array}\right]\left[\begin{array}{cc}  \sigma_{{x}} & \tau_{{\it xy}}\\  \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}}  \end{array}\right]=\left[\begin{array}{cc}  \sigma'_{{x}} & \tau'_{{\it xy}}\\  \noalign{\medskip}\tau'_{{\it xy}} & \sigma'_{{y}}  \end{array}\right]\ (6)

2.1 Eigenvalues, Eigenvectors and Principal Stresses

The second question is this: is there an angle (or set of direction cosines) where the shear stresses would go away, leaving only normal stresses? Because of the properties of the matrix, the answer to this is also “yes,” and involves the process of diagonalizing the matrix. The diagonalized matrix (the matrix with only non-zero values along the diagonal) can be found if the eigenvalues of the matrix can be found, i.e., if the following equation can be solved for \lambda :

\left[\begin{array}{cc}  \sigma_{{x}}-\lambda & \tau_{{\it xy}}\\  \noalign{\medskip}\tau_{{\it xy}} & \sigma_{{y}}-\lambda  \end{array}\right]=0\ (7)

To accomplish this, we take the determinant of the left hand side of Equation 7, namely

{\lambda}^{2}-\left(\sigma_{{x}}+\sigma_{{y}}\right)\lambda-\left({\tau_{{\it xy}}}^{2}-\sigma_{x}\sigma_{y}\right)=0 (8)

If we define

J_{1}=\sigma_{{x}}+\sigma_{{y}} (9)


J{}_{2}={\tau_{{\it xy}}}^{2}-\sigma_{x}\sigma_{y} (10)

Equation 8 can be rewritten as

\lambda^{2}-J_{1}\lambda-J_{2}=0 (11)

The quantities I_{1},\,I_{2} are referred to as the invariants; they do not change as the axes are rotated. The first invariant is the trace of the system, which was predicted to be invariant earlier. These are very important,

especially in finite element analysis, where failure criteria are frequently computed relative to the invariantsand not the principal stresses in any combination. This is discussed at length in Owen and Hinton [6].

This is the characteristic polynomial of the matrix of Equation 7. Most people generally don’t associatem atrices and polynomials, but every matrix has an associated characteristic polynomial, and conversely a polynomial can have one or more corresponding matrices.The solution of Equation 8 produces the two eigenvalues of the matrix. The first eigenvalue of the matrix is

\lambda_{1}=1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}  (12)

and the second

\lambda_{1}=1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}} (13)

We can also determine the eigenvectors from this. Without going into the process of determining these,the first eigenvector is

\bar{x}_{1}=\left[\begin{array}{c}  {\frac{-1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}}{\tau_{{\it xy}}}}\\  1  \end{array}\right]\   (14)

and the second is

\bar{x}_{2}=\left[\begin{array}{c}  {\frac{-1/2\,\sigma_{{y}}+1/2\,\sigma_{{x}}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}}{\tau_{{\it xy}}}}\\  1  \end{array}\right]   (15)

One thing we can do with these eigenvectors is to normalize them, i.e., have it so that the sum of the squares of the entries is unity. Doing this yields

\bar{x}_{1}=\left[\begin{array}{c}  {\frac{-1/2\,\sigma_{{y}}{\tau_{{\it xy}}}^{-1}+1/2\,\sigma_{{x}}{\tau_{{\it xy}}}^{-1}+1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}{\tau_{{\it xy}}}^{-1}}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}\\  {\frac{1}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}  \end{array}\right] (16)


\bar{x}_{2}=\left[\begin{array}{c}  {\frac{-1/2\,\sigma_{{y}}{\tau_{{\it xy}}}^{-1}+1/2\,\sigma_{{x}}{\tau_{{\it xy}}}^{-1}-1/2\,\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}{\tau_{{\it xy}}}^{-1}}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}\\  {\frac{1}{\sqrt{1/2\,{\frac{{\sigma_{{y}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-{\frac{\sigma_{{x}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{y}}}{{\tau_{{\it xy}}}^{2}}}+1/2\,{\frac{{\sigma_{{x}}}^{2}}{{\tau_{{\it xy}}}^{2}}}-1/2\,{\frac{\sqrt{{\sigma_{{y}}}^{2}-2\,\sigma_{{x}}\sigma_{{y}}+{\sigma_{{x}}}^{2}+4\,{\tau_{{\it xy}}}^{2}}\sigma_{{x}}}{{\tau_{{\it xy}}}^{2}}}+2}}}  \end{array}\right] (17)

The eigenvectors, complicated as the algebra is, are useful in that we can diagonalize the matrix as follows:

S^{-1}AS=\Lambda (18)


S=\left[\begin{array}{cc}  \bar{x}_{1_{1}} & \bar{x}_{2_{1}}\\  \bar{x}_{1_{2}} & \bar{x}_{2_{2}}  \end{array}\right] (19)

This is referred to as a similarity or collinearity transformation. Similar matrices are matrices with the same eigenvalues. The matrix \Lambda , although simpler in form than A , has the same eigenvalues as the “original”matrix.

Either the original or the normalized forms of the eigenvectors can be used; the result will be the same as the scalar multiples will cancel in the matrix inversion. The normalization was done to illustrate the
relationship between the eigenvectors, the diagonalization matrix, and the Givens rotation matrix, since for the last sin^{2}\alpha+cos^{2}\alpha=1, an automatically normalized relationship. It is thus possible to extract the angle of the principal stresses from the eigenvectors.

Inverting the result of Equation 19 and multiplying through Equation 18 yields at last

\Lambda=\left[\begin{array}{cc}  \lambda_{1} & 0\\  0 & \lambda_{2}  \end{array}\right]=\left[\begin{array}{cc}  \sigma_{1} & 0\\  0 & \sigma_{3}  \end{array}\right] (20)

The two eigenvalues from Equations 12 and 13 are the principal stresses at the stress point in question.(The use of different subscripts for the eigenvalues and principal stresses comes from too many years dealing with Mohr-Coulomb failure theory.) One practical result of Equations 9 and 20 and the underlying theory is that, for any coordinate orientation,

\sigma_{x}+\sigma_{y}=\sigma_{1}+\sigma_{3} (21)
This is a very handy check when working problems such as this, especially since the algebra is so involved.

2.2. Two-Dimensional Example

To see all of this “in action” consider the soil stress states shown in Figure 2.


Figure 2: Stress State Example (from Navy [5])

Consider the stress state “A.” The stresses are expressed as a ratio of a pressure P at the surface; furthermore, following geotechnical engineering convention, compressive stresses are positive and the ordinate is the “z” axis.. For this study the convention of Figure 1 will be adopted and thus \sigma_{x}=-0.77,\,\sigma_{y}=-0.98,\,\tau_{xy}=0.05   .By direct substitution (and carrying the results to precision unjustified in geotechnical engineering) we obtain the following:

  • \lambda_{1}=\sigma_{1}=-.7587029665P (Equations 12 and 20.)
  • \lambda_{2}=\sigma_{3}=-.9912970335P (Equations 12 and 20.) Note that the principal stresses are reversed; this is because the sign convention is reversed, and thus what was formerly the smaller of the stresses is now the larger. Also, the axis of the first principal stress changes because the first principal stress itself has changed.
  • S=\left[\begin{array}{cc}  .9754128670 & -.220385433\\  .2203854365 & .9754128502  \end{array}\right] (Equation 19.) Note that the normalized values of the eigenvectors (Equations 16 and 17) are used. Also note the similarity between this matrix and the Givens rotation matrix of Equation 5.  The diagonalization process is simply a rotation process, as the physics of the problem suggest. (The diagonal terms should be equal to each other and the absolute values ofthe off-diagonal terms should be also; they are not because of numerical errors in Maple where they were computed.
  • \alpha=12.73167230^{\circ};\,\beta=77.26832747^{\circ} (Equations 2, 3 and 6.)
  • In both cases the trace of A and \Lambda is the same, namely 1.75P (Equation 21.)

The results are thus the same. The precision is obviously greater than the “slide rule era” Navy [5], but the results illustrate that numerical errors can creep in, even with digital computation.

3. Three-Dimensional Analysis

3.1. Regular Principal Stresses

It is easy to see that, although the concepts are relatively simple, the algebra is involved.It is for this reason that the two-dimensional state was illustrated. With three dimensional stresses, the principles are the same, but the algebra is even more difficult.To begin, one more direction cosine must be defined,

n = cos\gamma (22)

where \gamma is of course the angle of the direction of the rotated coordinate system relative to the original one. The rotated system does not have to be the principal axis system, although that one is of most interest.
For three dimensions, Equation 4 should be written as follows:

\left[\begin{array}{ccc}  \sigma_{x} & \tau_{xy} & \tau_{xz}\\  \tau_{xy} & \sigma_{y} & \tau_{yz}\\  \tau_{xz} & \tau_{yz} & \sigma_{z}  \end{array}\right]\left[\begin{array}{c}  l\\  m\\  n  \end{array}\right]=\left[\begin{array}{c}  p_{x}\\  p_{y}\\  p_{z}  \end{array}\right] (23)

All of the moment summations of the shear stresses-which result in the symmetry of the matrix-have been included. The eigenvalues of the matrix are thus the solution of

\left[\begin{array}{ccc}  \sigma_{x}-\lambda & \tau_{xy} & \tau_{xz}\\  \tau_{xy} & \sigma_{y}-\lambda & \tau_{yz}\\  \tau_{xz} & \tau_{yz} & \sigma_{z}-\lambda  \end{array}\right]=0 (24)

The characteristic equation of this matrix-and thus the equation that solves for the eigenvalues-is

\lambda^{3}-J_{1}\lambda^{2}-J_{2}\lambda-J_{3}=0 (25)


J_1=\sigma_{x}+\sigma_{y}+\sigma_{z} (26)

J_2=\tau_{xy}^{2}+\tau_{xz}^{2}+\tau_{yz}^{2}-\left(\sigma_{x}\sigma_{y}+\sigma_{x}\sigma_{z}+\sigma_{y}\sigma_{z}\right) (27)

J_3=\sigma_{x}\sigma_{y}\sigma_{z}+2\tau_{xy}\tau_{xz}\tau_{yz}-\left(\sigma_{x}\tau_{yz}^{2}+\sigma_{y}\tau_{xz}^{2}+\sigma_{z}\tau_{xy}^{2}\right) (28)

It is easy to show the following:

  • Equation 26 reduces to Equation 9 if \sigma_{z}=0 .
  • Equation 27 reduces to Equation 10 if additionally \tau_{xz}=\tau_{yz}=0 .
  • For all three conditions, Equation 25 reduces to 11.

From the previous derivation, we can thus expect

 \lambda_{1}=\sigma_{1} (29)
 \lambda_{2}=\sigma_{2} (30)
\lambda_{3}=\sigma_{3} (31)
Solving a cubic equation such as this can be accomplished using the Cardan Formula. Generally, it is possible that we will end up with complex roots, but because of the properties of the matrix we are promised real roots. The eigenvalues are as follows: