Eigen decomposition for 3D stress tensor

by Vince Cronin

Begun September 13, 2010;  last revised September 8, 2017
Unfinished Draft;  Subject to Revision

Introduction

If we have a cubic free body that is arbitrarily oriented in a stress field and is in a state of static equilibrium, and we know the tractions (normal and shear stresses) acting on each of its infinitesimal faces, we can use eigen decomposition of the traction tensor to define the orientation and magnitude of the corresponding principal stresses.  The orientation of the orthogonal basis vectors that are parallel to the principal axes are the eigenvectors of the symmetric traction tensor, and the magnitudes of the principal stresses ((Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])) are the corresponding eigenvalues.

Imagine that we have a symmetrical 3-D matrix representing the traction tensor as follows:

s = (({{10, 3, 6}, {3, 4, 7}, {6, 7, 2}}))

We will use the Mathematica function Eigenvectors to define the eigenvectors of matrix s:

   Eigenvectors[s];
   
Similarly, we can use the Mathematica function Eigenvalues to define the eigenvectors of matrix s

       Eigenvalues[s];

Instructions

This code requires six input data :  the stress acting normal to the X plane ((Subscript[σ, xx])), the stress acting normal to the Y plane ((Subscript[σ, yy])), the stress acting normal to the Z plane ((Subscript[σ, zz])), the shear stress acting on the X plane in the Y direction ((Subscript[σ, xy])), the shear stress acting on the X plane in the Z direction ((Subscript[σ, xz])), and the shear stress acting on the Y plane in the Z direction ((Subscript[σ, yz])).  In accord with the convention generally used by geologists, a compressive normal stress has a positive sign, tensional normal stress is negative, and a shear stress that tends to make a body rotate in an anticlockwise direction is a positive shear stress.

Input (unsigned/positive values only)

In the input section below, type the stress values on the appropriate lines of code (which are shaded light red) after the "=" sign and before the ";" mark.

The variable (Subscript[σ, xx]) is the normal stress acting on the X plane -- the plane that is perpendicular to the X coordinate axis of the free body diagram.  Input the unsigned (positive) magnitude.

(Subscript[σ, xx])=2;

The variable (Subscript[σ, yy]) is the normal stress acting on the Y plane -- the plane that is perpendicular to the Y coordinate axis of the free body diagram.  Input the unsigned (positive) magnitude.

(Subscript[σ, yy])=1.5;

The variable (Subscript[σ, zz]) is the normal stress acting on the Z plane -- the plane that is perpendicular to the Z coordinate axis of the free body diagram.  Input the unsigned (positive) magnitude.

(Subscript[σ, zz])=1.7;

The variable (Subscript[σ, xy]) is the shear stress acting on the X plane in the Y direction.  Input the unsigned (positive) magnitude.

(Subscript[σ, xy])=0.7;

The variable (Subscript[σ, xy]) is the shear stress acting on the X plane in the Y direction.  Input the unsigned (positive) magnitude.

(Subscript[σ, xz])=0.3;

The variable (Subscript[σ, xy]) is the shear stress acting on the X plane in the Y direction.  Input the unsigned (positive) magnitude.

(Subscript[σ, yz])=0.5;

Computation

(Subscript[σ, yx])=(Subscript[σ, xy]);

(Subscript[σ, zx])=(Subscript[σ, xz]);

(Subscript[σ, zy])=(Subscript[σ, yz]);

s=(({{(Subscript[σ, xx]), (Subscript[σ, xy]), (Subscript[σ, xz])}, {(Subscript[σ, yx]), (Subscript[σ, yy]), (Subscript[σ, yz])}, {(Subscript[σ, zx]), (Subscript[σ, zy]), (Subscript[σ, zz])}}));

The eigenvectors that result from the next step are unit vectors that are colinear with the principal axes of the ellipse.  Consequently, the three eigenvectors ({result1[[1, 1]], result1[[1, 2]], result1[[1, 3]]}, {result1[[2, 1]], result1[[2, 2]], result1[[2, 3]]} and {result1[[3, 1]], result1[[3, 2]], result1[[3, 3]]}) are perpendicular to each other.  The eigenvalues are the magnitudes of the principal stresses.

principalAxes=Eigenvectors[s];

axisLengths=Sort[Eigenvalues[s], Greater];

The variable "a" is the length of the major axis of the principal ellipse;  "b" is the length of the intermediate axis;  "c" is the length of the minor axis.

(Subscript[σ, 1])=axisLengths[[1]];

(Subscript[σ, 2])=axisLengths[[2]];

(Subscript[σ, 3])=axisLengths[[3]];

Results

The eigenvectors are

principalAxes

The magnitude of the greatest principal stress ((Subscript[σ, 1])) is

(Subscript[σ, 1])

The corresponding axial vector is

(Subscript[vectorσ, 1])={((Subscript[σ, 1]) principalAxes[[1,1]]),((Subscript[σ, 1]) principalAxes[[1,2]]),((Subscript[σ, 1]) principalAxes[[1,3]])}

The magnitude of the intermediate principal stress ((Subscript[σ, 2])) is

(Subscript[σ, 2])

The corresponding axial vector is

(Subscript[vectorσ, 2])={((Subscript[σ, 2]) principalAxes[[2,1]]),((Subscript[σ, 2]) principalAxes[[2,2]]),((Subscript[σ, 2]) principalAxes[[2,3]])}

The magnitude of the least principal stress ((Subscript[σ, 3])) in is

(Subscript[σ, 3])

The corresponding axial vector is

(Subscript[vectorσ, 3])={((Subscript[σ, 3]) principalAxes[[3,1]]),((Subscript[σ, 3]) principalAxes[[3,2]]),((Subscript[σ, 3]) principalAxes[[3,3]])}

Graphics

The parametric equations of an ellipsoid can be written as
x = a Cos[u] Sin[v]
y = b Sin[u] Sin[v]
z = c Cos[v]
for u in[0, 2 π] and v in[0, π].

test1=Table[{(Subscript[σ, 1]) Cos[u]Sin[v],(Subscript[σ, 2]) Sin[u]Sin[v],(Subscript[σ, 3]) Cos[v]},{u,0,2π,(π/20)},{v,0,π,(π/20)}];

test2=Flatten[test1,1];

ListPointPlot3D[test2,AspectRatio->1,PlotRange->{{-(Subscript[σ, 1]),(Subscript[σ, 1])},{-(Subscript[σ, 1]),(Subscript[σ, 1])},{-(Subscript[σ, 1]),(Subscript[σ, 1])}}]

test4=ParametricPlot3D[{(Subscript[σ, 1]) Cos[u]Sin[v],(Subscript[σ, 2]) Sin[u]Sin[v],(Subscript[σ, 3]) Cos[v]},{u,0,2π},{v,0,π},ColorFunction->"GrayTones",PlotStyle->Specularity[White,50],AxesLabel->{"Greatest Principal Stress","Intermediate Principal Stress", "Least Principal Stress"},Mesh->None,AspectRatio->1,PlotRange->{{-(Subscript[σ, 1]),(Subscript[σ, 1])},{-(Subscript[σ, 1]),(Subscript[σ, 1])},{-(Subscript[σ, 1]),(Subscript[σ, 1])}}]

References

Ferguson, John, 1994, Introduction to Linear Algebra in Geology :  London, Chapman & Hall, 203 p., ISBN 0 - 412 - 49350 - 0

Fossen, Haakon, 2010, Structural Geology :  New York, Cambridge University Press, 463 p., ISBN 978 - 0 - 521 - 51664 - 8

Oertel, G., 1996, Stress and Deformation -- A Handbook on Tensors in Geology :  New York, Oxford University Press, 292 p., ISBN 0 - 19 - 509503 - 0

Pollard, D.D., and Fletcher, R.C., 2005, Fundamentals of Structural Geology :  New York, Cambridge University Press, 500 p., ISBN 0 - 521 - 83927 - 0;  see especially Chapter 6.

The graphics code to produce the ellipsoid came from Eric W.Weisstein,2009,http://mathworld.wolfram.com/notebooks/Surfaces/Ellipsoid.nb.
For more information,see Eric'sMathWorldentryhttp://mathworld.wolfram.com/Ellipsoid.html.

Wolfram MathWorld has a web resource on eigenvectors and eigen decomposition available at 
http : // mathworld.wolfram.com/Eigenvector.html

The corresponding resource on ellipsoids is available at
http://mathworld.wolfram.com/Ellipsoid.html

Color schemes used in Mathematica are displayed at
http://reference.wolfram.com/mathematica/guide/ColorSchemes.html

Other Mathematica resources are available for download at 
Emmanuel Amiot, "Eigenvectors by Hand" from The Wolfram Demonstrations Project  http://demonstrations.wolfram.com/EigenvectorsByHand/
Yaroslav Bulatov, "Linear Transformation with Given Eigenvectors" from The Wolfram Demonstrations Project  http://demonstrations.wolfram.com/LinearTransformationWithGivenEigenvectors/