ESM 5344			Spring 2008, Semester Project			Douglas A. Harold

Initialization

ClearAll[θ,ϕ,ρ,c,ν,λ,v1,v2,v3,vs1,vs2,vs3,α1,α2,α3,νLine,α1Line,α2Line,α3Line,X1,X2,X3,Xmax,k,Sol,inc](*Clears all previous asignments*)

The following creates Figure1

Figure1:=Graphics3D[{
{Blue,Thick,Line[{{0,0,0},{1,0,0}}]},
{Blue,Thick,Line[{{0,0,0},{0,1,0}}]},
{Blue,Thick,Line[{{0,0,0},{0,0,1}}]},
{Green,Thick,Line[{{0,0,0},{1,1,1}}]},
{Green,Thin,Line[{{0,0,0},{1,1,0}}]},
{LightBlue,Opacity[.7],Polygon[{{0,0,0},{1,1,1},{1,1,0}}]},
Text[Style["X1",Medium,Bold,Blue],{1,0,0}1/2,Background->Automatic],
Text[Style["X2",Medium,Bold,Blue],{0,1,0}1/2,Background->Automatic],
Text[Style["X3",Medium,Bold,Blue],{0,0,1}1/2,Background->Automatic],
Text[Style["ν",Large,Bold,Green],{1,1,1}1/2,Background->Automatic],
Text[Style["ϕ",Large,Bold,Black],{0.25,.125,0},Background->Automatic],
Text[Style["θ",Large,Bold,Black],{0.25,0.25,0.125},Background->Automatic]
}]

The following establishes the angle increment between successive values of θ & ϕ

inc=5°;(*May be changed if desired to improve resolution or reduce processing time*)

The following creates empty tables for use elsewhere.

c=Table[0,{i,6},{j,6}];(*All values are initially set to zero*)
λ = Table[0, {i, 3}, {j, 3}];(*All values are initially set to zero*)

v1=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)
v2=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)
v3=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)
α1=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)
α2=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)
α3=Table[0,{i,((180°)/inc)},{j,((360°)/inc)}];(*All values are initially set to zero*)

The following creates the dynamic Figure2

Figure2:=Manipulate[
{
X1={Blue,Thick,Line[{{0,0,0},{Xmax,0,0}}]};
X2={Blue,Thick,Line[{{0,0,0},{0,Xmax,0}}]};
X3={Blue,Thick,Line[{{0,0,0},{0,0,Xmax}}]}; 
(*X1, X2 & X3 are the axes defining the principal crystalLographic planes*)

νLine={Green,Thick,Line[{{0,0,0},v1[[(ϕ/inc),(θ/inc)]]}]};
(*Defines the Poynting vector*)

α1Line={Red,Thick,Line[{v1[[(ϕ/inc),(θ/inc)]],(v1[[(ϕ/inc),(θ/inc)]]+0.1Xmax α1[[(ϕ/inc),(θ/inc)]]),(v1[[(ϕ/inc),(θ/inc)]]-0.1Xmax α1[[(ϕ/inc),(θ/inc)]])}]}; 
(*Defines the direction of vibration for the largest velocity*)

α2Line={Red,Thick,Line[{v2[[(ϕ/inc),(θ/inc)]],(v2[[(ϕ/inc),(θ/inc)]]+0.1Xmax α2[[(ϕ/inc),(θ/inc)]]),(v2[[(ϕ/inc),(θ/inc)]]-0.1Xmax α2[[(ϕ/inc),(θ/inc)]])}]};
(*Defines the direction of vibration for the intermediate velocity*)

α3Line={Red,Thick,Line[{v3[[(ϕ/inc),(θ/inc)]],(v3[[(ϕ/inc),(θ/inc)]]+0.1Xmax α3[[(ϕ/inc),(θ/inc)]]),(v3[[(ϕ/inc),(θ/inc)]]-0.1Xmax α3[[(ϕ/inc),(θ/inc)]])}]};
(*Defines the direction of vibration for the minimum velocity*)

If[VelocitySurface, 
(*Determines if the Velocity Surface checkbox is checked*)

(*This section plots the entire velocity surfaces when the checkbox is checked*)
Graphics3D[{X1,X2,X3,νLine,α1Line,α2Line,α3Line,
{Lighter[Green,0.5],Point[vs1]},
{Lighter[Gray,0.25],Point[vs2]},
Point[vs3],
Text[Style["X1",Medium,Bold,Blue],{1,0,0}Xmax/2,Background->Automatic],
Text[Style["X2",Medium,Bold,Blue],{0,1,0}Xmax/2,Background->Automatic],
Text[Style["X3",Medium,Bold,Blue],{0,0,1}Xmax/2,Background->Automatic]},
Axes->True,AxesLabel->{x1,x2,x3},ImageSize->k,
PlotRange->{{-Xmax,Xmax},{-Xmax,Xmax},{-Xmax,Xmax}}],

(*This section plots the velocity surfaces in the plane defined by ϕ when the checkbox is not checked*)
Graphics3D[{X1,X2,X3,νLine,α1Line,α2Line,α3Line,
Point[v1[[ϕ/inc]]],
Point[v2[[ϕ/inc]]],
Point[v3[[ϕ/inc]]],
Text[Style["X1",Medium,Bold,Blue],{1,0,0}Xmax/2,Background->Automatic],
Text[Style["X2",Medium,Bold,Blue],{0,1,0}Xmax/2,Background->Automatic],
Text[Style["X3",Medium,Bold,Blue],{0,0,1}Xmax/2,Background->Automatic]},
Axes->True,AxesLabel->{x1,x2,x3},ImageSize->k,
PlotRange->{{-Xmax,Xmax},{-Xmax,Xmax},{-Xmax,Xmax}}]
] (*END IF*)
},
{{ϕ,inc},inc,180°,inc,Appearance->"Labeled" },
{{θ,inc},inc,360°,inc,Appearance->"Labeled" },
Delimiter,
{{Xmax,10000,"Zoom"},20000,1000,500},
{{k,700,"Image Size"},100,1000,1},
Delimiter,
{VelocitySurface,{False,True},Checkbox}
]

The following routine makes the stiffness matrix symmetric.

Programc:=
For[i=1,i<7,i++,
For[j=1,j<7,j++,
c[[j,i]]=c[[i,j]]
]]

The following defines the λ matrix.

λ[[1,1]]:=c[[1,1]] ((ν[[1]])^2)+c[[6,6]] ((ν[[2]])^2)+c[[5,5]] ((ν[[3]])^2)+
2 (c[[1,6]] ν[[1]] ν[[2]]+c[[1,5]] ν[[1]] ν[[3]]+c[[5,6]] ν[[2]] ν[[3]])

λ[[2,2]]:=c[[6,6]] ((ν[[1]])^2)+c[[2,2]] ((ν[[2]])^2)+c[[4,4]] ((ν[[3]])^2)+
2 (c[[2,6]] ν[[1]] ν[[2]]+c[[4,6]] ν[[1]] ν[[3]]+c[[2,4]] ν[[2]] ν[[3]])

λ[[3,3]]:=c[[5,5]] ((ν[[1]])^2)+c[[4,4]] ((ν[[2]])^2)+c[[3,3]] ((ν[[3]])^2)+
2 (c[[5,4]] ν[[1]] ν[[2]]+c[[5,3]] ν[[1]] ν[[3]]+c[[4,3]] ν[[2]] ν[[3]])

λ[[1,2]]:=c[[1,6]] ((ν[[1]])^2)+c[[6,2]] ((ν[[2]])^2)+c[[5,4]] ((ν[[3]])^2)+
 (c[[1,2]]+c[[6,6]]) ν[[1]] ν[[2]]+
(c[[1,4]]+c[[5,6]]) ν[[1]] ν[[3]]+
(c[[6,4]]+c[[5,2]]) ν[[2]] ν[[3]]

λ[[1,3]]:=c[[1,5]] ((ν[[1]])^2)+c[[6,4]] ((ν[[2]])^2)+c[[5,3]] ((ν[[3]])^2)+
 (c[[1,4]]+c[[6,5]]) ν[[1]] ν[[2]]+
(c[[1,3]]+c[[5,5]]) ν[[1]] ν[[3]]+
(c[[6,3]]+c[[5,4]]) ν[[2]] ν[[3]]

λ[[2,3]]:=c[[6,5]] ((ν[[1]])^2)+c[[2,4]] ((ν[[2]])^2)+c[[4,3]] ((ν[[3]])^2)+
 (c[[6,4]]+c[[2,5]]) ν[[1]] ν[[2]]+
(c[[6,3]]+c[[4,5]]) ν[[1]] ν[[3]]+
(c[[2,3]]+c[[4,4]]) ν[[2]] ν[[3]]

The following routine makes the λ matrix symmetric.

Programλ:=
For[i=1,i<4,i++,
For[j=1,j<4,j++,
λ[[j,i]]=λ[[i,j]]
]]

The following reorganizes the velocities into a form that is more convenient for plotting the velocity surfaces.

vs1:=Flatten[v1,1];
vs2:=Flatten[v2,1];
vs3:=Flatten[v3,1];

The following program solves the eigenvalue problem for each value of ν

Program1:=
For[{i=1,ϕ=-inc},i≤((180°)/inc),i++,
{ϕ=ϕ+inc,
For[{j=1,θ=-inc},j≤((360°)/inc),j++,
{θ=θ+inc,

Sol=Eigensystem[λ],
(*Find the eigenvalues & eigenvectors for each successive direction ν*)

v1[[i,j]]=ν(Sqrt[((Sol[[1,1]])/ρ)]),
(*Largest velocity*)

v2[[i,j]]=ν(Sqrt[((Sol[[1,2]])/ρ)]),
(*Intemediate velocity*)

v3[[i,j]]=ν(Sqrt[((Sol[[1,3]])/ρ)]),
(*Smallest velocity*)

α1[[i,j]]=Sol[[2,1]], 
(*Vibration direction corresponding to the largest wave velocity*)

α2[[i,j]]=Sol[[2,2]],
(*Vibration direction corresponding to the intermediate wave velocity*)

α3[[i,j]]=Sol[[2,3]]    
(*Vibration direction corresponding to the smallest wave velocity*)

}]
}];

Phase Velocity Surfaces

Please note that if asked one should "Enable Dynamic" and then evaluate the notebook before scrolling down. Furthermore, the notebook should be reevaluated whenever material properties are changed. Lastly, note that in order to clean up this section and make it easier to read, many programs and functions have been moved into the Initialization section above. That section may be expanded by clicking on the grey triangle on the left.

SET UP

Identify Variables

Definition of variables:

θ 		Angle up from the X1 - X2 plane
ϕ  		Angle about the X3 axis from the X1 - X3 plane

ρ  		Mass density of the material
c  		Stiffness Tensor written in symmetric 6 X 6 matrix form
λ  		3 X 3 symmetric matrix determined by c & ν, as follows
			(Subscript[λ, ij]) = (Subscript[c, ijkl]) (Subscript[ν, k])(Subscript[ν, l])

ν, νLine  	Poynting Vector of unit length

v1, vs1  		Largest velocity in the ν direction
v2, vs2  		Intermediate velocity in the ν direction
v3, vs3  		Smallest velocity in the ν direction

α1, α1Line  	Vibration direction for the largest velocity in the ν direction
α2, α2Lline  	Vibration direction for the intermediate velocity in the ν direction
α3, α3Line  	Vibration direction for the smallest velocity in the ν direction

X1. X2, X3  	Axes defining the principal crystallographic planes

Xmax, k  	variables used for scaling the output
Sol  		Variable used for temporarily storing the eigensystem solution

inc  		Increment angle between successive values of θ or ϕ

Introduction

1. Set up the 6 x 6 stiffness matrix

Although the stiffness tensor (Subscript[C, ijkl]) is really a (4^th)order tensor with as many as 81 independent constants, symmetry on both the stress and strain tensors reduce the number of independent constants to 36 (which may be written in a 6 x 6 matrix. Furthermore, if a strain energy function exists for the material then this matrix is symmetric and the number of independent elastic constants is reduces to 21. This will be assumed for this solution.

Define the experimental values for the material below. These may be changed as desired, but one should use base units or the equivalent. For example, one could use Pa, but should avoid using kPa, etc. The stiffness matrix is a 6 x 6 symmetric matrix. Therefore, one should only define values on or above the main diagonal. Thus, as many as 21 independent constants may defined. Note the syntax, (Subscript[c, ij]) = c[[i,j]] and that i & j each vary from 1 to 6.

Mass density ρ of the material ((kg/(m^3))):

ρ=2.02((10^3));

Stiffness constants (Pa):

c[[1,1]]=4.92((10^10));
c[[2,2]]=2.44((10^10));
c[[3,3]]=3.54((10^10));
c[[4,4]]=1.05((10^10));
c[[5,5]]=1.22((10^10));
c[[6,6]]=2.82((10^10));
c[[1,2]]=2.48((10^10));
c[[1,3]]=2.45((10^10));
c[[2,3]]=1.45((10^10));

The following program makes the stiffness matrix symmetric:

Programc;

Here is what the stiffness matrix now looks like:

MatrixForm[c]

2. Define the ν vector

ν:={Cos[θ]Cos[ϕ],Cos[θ]Sin[ϕ],Sin[θ]};

The variables θ & ϕ define the direction of the unit vector ν as shown in Figure1.

Figure1

The phase velocities and vibration directions are solutions to the eigenvalue problem defined by Christoffel's Equation for Linear Elastic Homogeneous Anisotropic (LEHA) material.

	((Subscript[C, ijkl])(Subscript[ν, j])(Subscript[ν, l])-((ρv)^2)(Subscript[δ, ik]))(Subscript[α, k])=0

It is convenient to make the following definition:

	(Subscript[λ, ij]) = (Subscript[c, ijkl]) (Subscript[ν, k])(Subscript[ν, l])

Then it is obvious that the velocities in the direction defined by ν is found by computing the eivenvalue of the 3 x 3 matrix:

	(Subscript[λ, ij])-((ρv)^2)(Subscript[δ, ij])

The eigenvalues are ((ρv)^2) ; therefore, the velocities are found as follows:

	v=(Sqrt[(Eigenvalue/ρ)])
	
Furthermore, for each eigenvalue there is a corresponding eigenvector. This vector (Subscript[α, k]) defines the direction of vibration for a wave traveling in the direction ν at the wave speed v. As expected from linear algebra, there are three wave velocities and three corresponding vibration directions for each direction ν. Usually these velocities are distinct, but are not necessarily so. Whenever degeneracies in the eigenvalue problem occur (i.e. two or more velocities are equal for a given direction ν), the vibration directions are found by assuming that they vary smoothly through the degeneracy. This effect can be observed in the dynamic simulation below (Figure2).

Additionally, as is obvious in Figure2, the greatest wave speed is usually associated with a longitudinal (or quasilongitudinal) wave, while the two slower waves are usually transverse (or quasitransverse).

SOLUTION

Computation

Now fill in the values for the six independent values (those on & above the main diagonal) according to the following definition:

	(Subscript[λ, ij]) = (Subscript[c, ijkl]) (Subscript[ν, k])(Subscript[ν, l])

The following program makes λ symmetric:

Programλ;

This program finds the solution to the eigenvalue problem for each value of θ & ϕ.

Program1;

Graphics

Figure2 is a dynamic plot of the results found earlier. Note that the vibration directions are shown by the red lines for the ν direction shown in green. The plot is much more responsive when the velocity surfaces are not plotted. Thus, one should set the zoom level and the ν direction (by adjusting θ & ϕ) before showing the velocity surfaces. To rotate the graph, use the mouse to click and drag the image.

Figure2