Eigen decomposition for 2D stress tensor

by Vince Cronin

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

Introduction

If we have a cubic free body that is  oriented such that the X and Y axes are in the (Subscript[σ, 1])- (Subscript[σ, 3]) plane, 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[σ, 3])) are the corresponding eigenvalues.

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

s = (({{10, 3}, {3, 4}}))

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];

Wolfram MathWorld has a web resource on eigenvectors and eigen decomposition available at 
http : // mathworld.wolfram.com/Eigenvector.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/

Instructions

This code requires three input data :  the stress acting normal to the X plane ((Subscript[σ, xx])), the stress acting normal to the Y plane ((Subscript[σ, yy])), and the shear stress acting on the X plane in the Y direction ((Subscript[σ, xy])).  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

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.  By the structural geologists' sign convention, a compressive normal stress is a positive value, and a tensional normal stress is a negative value.

(Subscript[σ, xx])=19;

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.  By the structural geologists' sign convention, a compressive normal stress is a positive value, and a tensional normal stress is a negative value.

(Subscript[σ, yy])=10;

The variable (Subscript[σ, xy]) is the shear stress acting on the X plane in the Y direction.  By the structural geologists' sign convention, a shear stress has a positive value if it tends to rotate the body in an anticlockwise direction around the origin of the coordinate system.

(Subscript[σ, xy])=-5;

Computation

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

s=(({{(Subscript[σ, xx]), (Subscript[σ, xy])}, {(Subscript[σ, yx]), (Subscript[σ, yy])}}));

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

result1=Eigenvectors[s];

The angle between the X axis and the nearest principal stress axis is called "theta" in the code below.

theta=If[((Subscript[σ, xy])≠0),(-1*(ArcTan[result1[[1,2]]/result1[[1,1]]])),0];

The eigenvalues are the magnitudes of the principal stresses/strains (depending on whether the input data are stresses or strains)

axisLengths=Eigenvalues[N[s]];

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

a=If[(axisLengths[[1]]>axisLengths[[2]]),axisLengths[[1]],axisLengths[[2]]];

b=If[(axisLengths[[1]]<axisLengths[[2]]),axisLengths[[1]],axisLengths[[2]]];

rotatedEllipse=Table[{(a Cos[t] Cos[theta])+(b Sin[t]*(-Sin[theta])),(a Cos[t] Sin[theta])+(b Sin[t] Cos[theta])},{t,0,2π,π/36}];

The first plot file (plot1) plots an ellipse with semimajor axis length of "a", semiminor axis length "b", and major axis inclined "—theta" degrees (i.e., "angle" radians) relative to the X axis of the free body.

plot1=ListLinePlot[%, AspectRatio ->1, AxesLabel->{"X axis of free body","Y axis of free body"},PlotRange->{{-a,a},{-a,a}}, PlotStyle->{Thick,Blue}];

The second plot file (plot2) creates a line along the major axis of the ellipse.

plot2 = ListLinePlot[{{a Cos[theta], a Sin[theta]}, {-a Cos[theta], -a Sin[theta]}}, AspectRatio->1,PlotRange->{{-a,a},{-a,a}}, PlotStyle->{Red}];

The third plot file (plot3) creates a line along the minor axis of the ellipse.

plot3 = ListLinePlot[{{b Cos[theta + (π/2)], b Sin[theta + (π/2)]}, {-b Cos[theta + (π/2)], -b Sin[theta + (π/2)]}}, AspectRatio->1,PlotRange->{{-a,a},{-a,a}}, PlotStyle->{Red}];

The fourth plot file (plot4) merges the previous three plot files into a single file for display.

plot4=Show[plot1,plot2,plot3];

Results

The magnitude of the greatest principal stress (a) in megaPascals is

a

The magnitude of the least principal stress (b) in megaPascals is

b

The angle (theta) between the X axis of the free body and the greatest principal stress axis (longer red line in plot below), expressed in degrees, is

N[theta*(180/π)]

where a positive value of theta means that the greatest principal stress axis is found with an anticlockwise rotation of theta degrees from the X axis.

The plot of the stress ellipse (blue) relative to the coordinate axes that are established in the free body (plot4) is below.  The axes are labeled in megaPascals.

plot4

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.

Wolfram MathWorld has a web resource on eigenvectors and eigen decomposition available at 
http : // mathworld.wolfram.com/Eigenvector.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/