The Elasticity Tensor

how to calculate and plot the elasticity properties

MTEX offers a very simple way to compute elasticity properties of materials. This includes Young's modulus, linear compressibility, Christoffel tensor, and elastic wave velocities.

On this page ...
Import an Elasticity Tensor
Young's Modulus
Linear Compressibility
Christoffel Tensor
Elastic Wave Velocity

Import an Elasticity Tensor

Let us start by importing the elastic stiffness tensor of an Olivine crystal in reference orientation from a file.

fname = fullfile(mtexDataPath,'tensor','Olivine1997PC.GPa');

cs = crystalSymmetry('mmm',[4.7646 10.2296 5.9942],'mineral','Olivin');

C = stiffnessTensor.load(fname,cs)
 
C = stiffnessTensor  
  unit   : GPa              
  rank   : 4 (3 x 3 x 3 x 3)
  mineral: Olivin (mmm)     
 
  tensor in Voigt matrix representation:
 320.5  68.2  71.6     0     0     0
  68.2 196.5  76.8     0     0     0
  71.6  76.8 233.5     0     0     0
     0     0     0    64     0     0
     0     0     0     0    77     0
     0     0     0     0     0  78.7

Young's Modulus

Young's modulus is also known as the tensile modulus and measures the stiffness of elastic materials. It is computed for a specific direction d by the command YoungsModulus.

d = vector3d.X
E = C.YoungsModulus(vector3d.X)
 
d = vector3d  
 size: 1 x 1
  x y z
  1 0 0
E =
  286.9284

If the direction is ommited Youngs modulus is returned as a spherical function.

% compute Young's modulus as a directional dependend function
E = C.YoungsModulus

% which can be evaluated at any direction
E.eval(d)

% or plot it
setMTEXpref('defaultColorMap',blue2redColorMap);
plot(C.YoungsModulus,'complete','upper')
 
E = S2FunHarmonic  
 mineral: Olivin (mmm)
 bandwidth: 8
ans =
  286.1748

Linear Compressibility

The linear compressibility is the deformation of an arbitrarily shaped specimen caused by an increase in hydrostatic pressure and can be described by a second rank tensor. Similarly as the Youngs modulus it can be computed by the command linearCompressibility for specific directions d or as a spherical function

% compute as a spherical function
beta = linearCompressibility(C)

% plot it
plot(beta,'complete','upper')

% evaluate the function at a specific direction
beta.eval(d)
 
beta = S2FunHarmonic  
 mineral: Olivin (mmm)
 bandwidth: 2
 antipodal: true
ans =
    0.0018

Christoffel Tensor

The Christoffel Tensor is symmetric because of the symmetry of the elastic constants. The eigenvalues of the 3x3 Christoffel tensor are three positive values of the wave moduli which corresponds to \rho Vp^2 , \rho Vs1^2 and \rho Vs2^2 of the plane waves propagating in the direction n. The three eigenvectors of this tensor are then the polarization directions of the three waves. Because the Christoffel tensor is symmetric, the polarization vectors are perpendicular to each other.

% It is computed for a specific direction x by the
% command <tensor.ChristoffelTensor.html ChristoffelTensor>.

T = ChristoffelTensor(C,vector3d.X)
 
T = ChristoffelTensor  
  rank   : 2 (3 x 3)   
  mineral: Olivin (mmm)
 
 320.5     0     0
     0  78.7     0
     0     0    77

Elastic Wave Velocity

The Christoffel tensor is the basis for computing the direction dependent wave velocities of the p, s1, and s2 wave, as well as of the polarization directions. Therefore, we need also the density of the material, e.g.,

rho = 3.355
rho =
    3.3550

which we can write directly into the ellastic stiffness tensor

C = addOption(C,'density',rho)
 
C = stiffnessTensor  
  unit   : GPa              
  density: 3.355            
  rank   : 4 (3 x 3 x 3 x 3)
  mineral: Olivin (mmm)     
 
  tensor in Voigt matrix representation:
 320.5  68.2  71.6     0     0     0
  68.2 196.5  76.8     0     0     0
  71.6  76.8 233.5     0     0     0
     0     0     0    64     0     0
     0     0     0     0    77     0
     0     0     0     0     0  78.7

Then the velocities are computed by the command velocity

[vp,vs1,vs2,pp,ps1,ps2] = velocity(C)
 
vp = S2FunTri  
 
 vertices: 1 x 18338
 values:   1 x 18338
 
vs1 = S2FunTri  
 
 vertices: 1 x 18338
 values:   1 x 18338
 
vs2 = S2FunTri  
 
 vertices: 1 x 18338
 values:   1 x 18338
pp = 
  S2AxisFieldTri with properties:

          tri: [1×1 S2Triangulation]
       values: [1×18338 vector3d]
     vertices: [1×18338 S2Grid]
    antipodal: 0
ps1 = 
  S2AxisFieldTri with properties:

          tri: [1×1 S2Triangulation]
       values: [1×18338 vector3d]
     vertices: [1×18338 S2Grid]
    antipodal: 0
ps2 = 
  S2AxisFieldTri with properties:

          tri: [1×1 S2Triangulation]
       values: [1×18338 vector3d]
     vertices: [1×18338 S2Grid]
    antipodal: 0

In order to visualize these quantities, there are several possibilities. Let us first plot the direction dependent wave speed of the p-wave

plot(vp,'complete','upper')
ans = 
  Contour with properties:

    LineColor: [0 0 0]
    LineStyle: '-'
    LineWidth: 0.5000
         Fill: 'on'
    LevelList: [1×10 double]
        XData: [91×361 double]
        YData: [91×361 double]
        ZData: [91×361 double]

  Use GET to show all properties

Next, we plot on the top of this plot the p-wave polarization direction.

hold on
plot(pp)
hold off

We may even compute with these spherical functions as width ordinary values. E.g. to visualize the speed difference between the s1 and s2 waves we do.

plot(vs1-vs2,'complete','upper')

hold on
plot(ps1)
hold off
ans = 
  Contour with properties:

    LineColor: [0 0 0]
    LineStyle: '-'
    LineWidth: 0.5000
         Fill: 'on'
    LevelList: [1×10 double]
        XData: [91×361 double]
        YData: [91×361 double]
        ZData: [91×361 double]

  Use GET to show all properties

When projected to a plane the different wave speeds

planeNormal = vector3d.X;

% options for sections
optSec = {'color','interp','linewidth',6,'doNotDraw'};

% options for quiver
optQuiver = {'linewidth',2,'autoScaleFactor',0.35,'doNotDraw'};
optQuiverProp = {'color','k','linewidth',2,'autoScaleFactor',0.25,'doNotDraw'};
prop = S2VectorFieldHarmonic.normal; % the propagation direction

% wave velocyties
close all
plotSection(vp,planeNormal,optSec{:},'DisplayName','Vp')
hold on
plotSection(vs1,planeNormal,optSec{:},'DisplayName','Vs1')
plotSection(vs2,planeNormal,optSec{:},'DisplayName','Vs2')

% polarisation directions
quiverSection(vp,pp,planeNormal,'color','c',optQuiver{:},'DisplayName','pp')
quiverSection(vs1,ps1,planeNormal,'color','g',optQuiver{:},'DisplayName','ps1')
quiverSection(vs2,ps2,planeNormal,'color','m',optQuiver{:},'DisplayName','ps2')

% plot propagation directions as reference
quiverSection(vp,prop,planeNormal,optQuiverProp{:},'DisplayName','x')
quiverSection(vs1,prop,planeNormal,optQuiverProp{:})
quiverSection(vs2,prop,planeNormal,optQuiverProp{:})
hold off

axis off tight
legend('Vp','Vs1','Vs2','pp','ps1','ps2','x','Location','eastOutSide')
mtexTitle('Phase velocity surface (km/s)')

mtexColorMap blue2red
mtexColorbar('Title','(km/s)','location','southOutSide')

Similarly, we can visualize the slowness surfaces (s/km)

% plot slowness surfaces
plotSection(1./vp,planeNormal,optSec{:},'DisplayName','Vp')
hold on
plotSection(1./vs1,planeNormal,optSec{:},'DisplayName','Vs1')
plotSection(1./vs2,planeNormal,optSec{:},'DisplayName','Vs2')

% polarisation directions
quiverSection(1./vp,pp,planeNormal,'color','c',optQuiver{:},'DisplayName','pp')
quiverSection(1./vs1,ps1,planeNormal,'color','g',optQuiver{:},'DisplayName','ps1')
quiverSection(1./vs2,ps2,planeNormal,'color','m',optQuiver{:},'DisplayName','ps2')

% plot propagation directions as reference
quiverSection(1./vp,prop,planeNormal,optQuiverProp{:},'DisplayName','x')
quiverSection(1./vs1,prop,planeNormal,optQuiverProp{:})
quiverSection(1./vs2,prop,planeNormal,optQuiverProp{:})
hold off
axis off tight
legend('Vp','Vs1','Vs2','pp','ps1','ps2','x','Location','eastOutSide')
mtexTitle('Slowness surface (km/s)')

mtexColorMap blue2red
mtexColorbar('Title','(s/km)','location','southOutSide')

set back default colormap

setMTEXpref('defaultColorMap',WhiteJetColorMap)