DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN CERAMIC
MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classified information.
___________________________
Chad Matthew Spurlock
Certificate of Approval:
___________________ ___________________
Mary L. Hughes Chai Yoo
Assistant Professor Professor
Civil Engineering Civil Engineering
___________________ ___________________
G. Ed Ramey Stephen L. McFarland
Professor Dean
Civil Engineering Graduate School
DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN CERAMIC
MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS
Chad Matthew Spurlock
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
December 15, 2006
iii
DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN CERAMIC
MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS
Chad Matthew Spurlock
Permission is granted to Auburn University to make copies of this thesis at its discretion,
upon request of individuals or institutions and at their expense. The author reserves all
publication rights.
______________________________
Signature of the Author
______________________________
Date
iv
THESIS ABSTRACT
DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN CERAMIC
MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS
Chad Matthew Spurlock
Master of Science, December 15, 2006
(B.S., Mathematics, Mississippi State University, 2003)
110 Typed Pages
Directed by Mary L. Hughes
MATLAB programs were developed to compute the elastic constants and damage
of ceramic matrix composites using ultrasonic wave speed measurements. The matrix of
elastic stiffnesses, C, relates the material stresses and strains. The ultrasonic wave
velocities are related to the elastic constants through the Christoffel equation. The
immersion ultrasonic wave speed measurement method has been used extensively to
determine elastic constants of anisotropic media. The computer programs contained
herein are designed to recover the elastic stiffnesses and damage magnitudes for materials
of orthotropic symmetry (i.e. having nine independent elastic stiffnesses) from data
generated from the ultrasonic wave speed measurements. The elasticity matrix, C, is
recovered using three methods: minimization of the least-squares of the Christoffel
equation, minimization of the sum of squares of the deviations between experimental
velocities and the solution of the Christoffel equation, and use of the rotation of axes
v
equation for fourth-order tensors. Damage is defined in the macroscopic sense as the
normalized variation of the elastic stiffnesses under loading.
vi
Style manual or journal used: A Manual for Writers of Term Papers, Theses, and
Dissertations. Sixth Edition. Kate L. Turabian.__________________________________
Computer software used: MATLAB 7.0.4, Microsoft Word, Microsoft Excel,
Mathematica 4.1, Adobe Illustrator___________________________________________
vii
TABLE OF CONTENTS
LIST OF TABLES.............................................................................................................. x
LIST OF FIGURES ........................................................................................................... xi
CHAPTER I INTRODUCTION..........................................................................................1
1.1 Background........................................................................................................1
1.2 Objective............................................................................................................2
CHAPTER II LITERATURE REVIEW .............................................................................3
2.1 Composite Materials ..........................................................................................3
2.2 Specific Types of Ceramic Matrix Composite Materials; Applications............7
2.3 Mechanical Behavior .........................................................................................8
2.4 Types of Material Symmetry ...........................................................................11
2.5 Rotation of Axes ..............................................................................................13
2.6 Micromechanics and Macromechanics............................................................14
2.7 Testing Methods...............................................................................................14
2.8 Ultrasonics .......................................................................................................15
CHAPTER III DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN
CERAMIC MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS...........................................................................................................19
3.1 Ultrasonic Test Methods..................................................................................19
3.2 The Time Difference between the Group and Phase Velocities ......................22
viii
3.3 Calculating the Phase Velocity and Refraction Angle.....................................28
3.4 Wave Behavior.................................................................................................30
3.5 The Wave Vector .............................................................................................33
3.6 Solving for the Elastic Stiffnesses, C...............................................................40
3.6.1 Minimization of Velocities ...............................................................40
3.6.2 Minimization of the Christoffel Equation.........................................41
3.6.3 Rotation Method ...............................................................................42
3.7 Solving the Systems of Equations....................................................................43
3.7.1 The Levenberg-Marquardt Method...................................................44
3.7.2 Least Squares Solution......................................................................45
3.8 Damage ............................................................................................................45
CHAPTER IV DESCRIPTION OF MATLAB PROGRAMS ..........................................48
4.1 Introduction......................................................................................................48
4.2 Minimization of the Christoffel Equation........................................................50
4.3 Minimization of the Variation between Experimental and Calculated
Velocities ...............................................................................................................50
4.4 Least Squares Solution using Rotation of Axes...............................................50
CHAPTER V COMPARISON OF MATLAB PROGRAMS ...........................................51
5.1 Reconstruction of Velocity Data......................................................................51
5.2 Results..............................................................................................................58
CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS....................................60
6.1 Conclusions......................................................................................................60
6.2 Recommendations............................................................................................61
ix
REFERENCES ..................................................................................................................63
APPENDIX A Minimization of the Christoffel Equation .................................................65
APPENDIX B Minimization of the Variation between Experimental and Calculated
Velocities ...........................................................................................................................74
APPENDIX C Least Squares Solution using Rotation of Axes ........................................86
x
LIST OF TABLES
2.1 Properties of some high performance ceramics.............................................................7
2.2 Some important ceramic reinforcements .......................................................................7
2.3 Examples of ceramic/ceramic composites.....................................................................8
2.4 Notations of elastic stiffnesses of the triclinic, monoclinic, orthotropic, tetragonal, and
cubic symmetry systems ....................................................................................................12
2.5 Notations of elastic stiffnesses of the triclinic, trigonal, transversely isotropic, and
isotropic symmetry systems...............................................................................................13
5.1 Values of elastic stiffnesses for 1-D SiC-SiC..............................................................52
5.2 Simulated velocity and angle measurements for quasi-longitudinal waves ................55
5.3 Simulated velocity and angle measurements for fast transverse waves ......................56
5.4 Simulated velocity and angle measurements for slow transverse waves.....................57
5.5 Comparison of results of MATLAB programs and Aristegui and Baste data.............59
xi
LIST OF FIGURES
2.1 Types of reinforcement in composites...........................................................................3
2.2 Toughening mechanisms in composites ........................................................................5
3.1 General test setup for an ultrasonic immersion test.....................................................20
3.2 Rotation of sample in order to excite transverse waves...............................................20
3.3 Under-load ultrasonic device .......................................................................................22
3.4 Representation of the times for different acoustic paths..............................................24
3.5 Schematic of different acoustic paths for the group and phase velocity vectors .........25
3.6 Schematic diagram of waves excited in a non-principal plane....................................31
3.7 Schematic of the composite plate ................................................................................33
3.8 View of the incident plane ...........................................................................................38
4.1 Initial orientations of sample and ultrasonic transmitter and receiver.........................49
5.1 Unidirectional fiber composite showing the axes defining the three planes of
measurement ......................................................................................................................51
5.2 Reconstructed slowness curves and experimentally measured data points for a 1-D
SiC-SiC composite.............................................................................................................52
1
CHAPTER I
INTRODUCTION
1.1 Background
Composite materials hold great promise for future structural applications. They
are very attractive for their high strength-to-weight ratio but also exhibit many other
desirable characteristics, such as heat resistance and hardness. Composite materials are
composed of two components: the reinforcement and the matrix. The reinforcement
provides the load carrying capacity for the composite and in many cases is in the form of
extremely high strength fibers. The matrix serves as a bond for the reinforcement and
transfers stresses between fibers. Both the matrix and fibers may be metallic, polymeric,
or ceramic; ceramic matrix composites are a relatively recent development in the
composites field (Chawla 2003) and are the main composite material of interest in this
study. The use of composite materials for structural applications has increased greatly in
the last decade; however, for composite materials to be implemented in structural
applications with confidence, their behavior under loading must be studied more closely.
Damage develops in composite materials through a variety of mechanisms, such
as matrix cracking, fiber pullout, and fiber-matrix debonding. These mechanisms
increase the toughness of composite materials (the process by which this occurs will be
discussed in Chapter 2), which are often combinations of extremely strong but brittle
materials. The development of damage negatively affects the values of the elastic
2
constants of the composite, making a complete description of material properties over a
range of stresses crucial for design.
Ultrasonic testing methods have been used widely to determine the elastic
constants of materials (Aristegui and Baste 1997, Mouchtachi et al. 2004, Rokhlin and
Wang 1992). The advantages of these methods over traditional tensile and flexure tests
are that ultrasonic tests are nondestructive and that the full set of elastic constants can be
determined from a single sample, which is especially important for composites whose
properties may vary significantly due to manufacturing irregularities.
1.2 Objective
In this work, three MATLAB programs have been developed to compute the
amount of damage and the elastic stiffnesses for materials with orthotropic or higher
(fewer independent elastic constants) symmetry using data collected from a one-way
through-transmission ultrasonic test conducted while the specimen is simultaneously
being subjected to uniaxial tension. The user inputs time readings, applied stress levels,
and orientation of the sample; in return, the program calculates the elastic stiffnesses,
elastic compliances, damage values, and pertinent material constants. The values are
recovered using three methods: minimization of the least-squares of the Christoffel
equation (Achenbach 1973), minimization of the sum of squares of the deviations
between experimental velocities and the solutions of the Christoffel equation, and use of
the rotation of axes equation for fourth-order tensors.
3
CHAPTER II
LITERATURE REVIEW
2.1 Composite Materials
A composite material is composed of two or more distinct components, the
reinforcement and the matrix. Each constituent may be homogeneous, but the composite
material is macroscopically heterogeneous (Suhling 2004). The reinforcement can be
distributed throughout the matrix as either particles or fibers; several reinforcement
schemes are shown in Figure 2.1. Typically, the different components of a composite can
be distinguished visually (Jones 1999).
Figure 2.1 Types of reinforcement in composites (Warren 1992).
Particulate composite materials are composed of three-dimensional particles
suspended in a matrix. Particles with a principally two-dimensional geometry are called
?flakes? or ?platelets? (Jones 1999). Fibrous reinforcement may exist as either
4
continuous or discontinuous strands, and continuous fibers may be oriented in any
direction, although unidirectional fibers are the type most often studied for structural
applications (Suhling 2004). Discontinuous fibers, also called ?whiskers? or ?short
fibers?, are capable of higher strengths than long fibers (Jones 1999), but the strength of a
fiber is used only if it is aligned in the direction of tensile loading (Wachtman 1989).
The reinforcement and matrix constituents of a composite material can be
metallic, polymeric, or ceramic. Composites are commonly classified according to their
matrix material. Three general classifications of composite materials are polymer matrix
composites (PMCs), metal matrix composites (MMCs), and ceramic matrix composites
(CMCs). Of these, ceramic matrix composites are the most recent entrants into the field
of composites (Chawla 2003).
The main advantage of using composite materials is that they can be designed to
meet specific design needs. If well-designed, a composite exhibits qualities of each of its
constituents as well as unique qualities neither material possesses alone. Material
characteristics that may be improved by forming a composite are strength, stiffness,
corrosion resistance, wear resistance, strength-to-weight ratios, fatigue life, extreme
temperature response, thermal conductivity, acoustical properties, crack (fracture)
resistance, cost, and ease of fabrication. The form, structural arrangement, material
properties, and interactions of a composite?s constituents determine the behavior and
properties of that composite material (Suhling 2004).
As mentioned previously, ceramic matrix composites define a class of composite
materials. Ceramics are chemically stable, inorganic, crystalline, non-metallic
compounds (Warren 1992). Ceramic materials are desirable for their high strength and
5
hardness, low density, and high resistance to chemicals, wear, and heat, but they are
limited by their brittleness, or low toughness, which makes ceramic materials highly
sensitive to flaws. The primary goal of using ceramic materials in composite form is to
retain the attractive qualities while increasing the fracture toughness (Chawla 2003).
Improvements in processing have reduced the size and frequency of defects in ceramics,
but certain features of composite materials, such as boundaries between constituents, act
as flaws. Furthermore, defects formed during service weaken brittle materials (Warren
1992). Some specific toughening mechanisms achieved by introducing fiber
reinforcement to a ceramic material are crack deflection, fiber bridging, fiber pullout,
microcracking, and debonding at the fiber/matrix interface (Chawla 2003). These
processes are illustrated in Figure 2.2.
Figure 2.2 Toughening mechanisms in composites (Warren 1992).
6
For a fiber reinforced composite material, failure generally occurs by cracking
and subsequent failure of the matrix followed by failure of the fibers. In many instances,
a developing crack may be deflected around a reinforcing strand rather than pass through
the fiber. This type of crack deflection produces an apparent toughness increase by
reducing the stress intensity at the crack tip.
As cracks grow in the matrix, the stronger fibers remain intact and bridge the gaps
in the matrix. Fiber bridging allows the continued loading of the composite even after
matrix failure. The energy expended in pulling a fiber from the matrix can increase the
overall toughness of the composite material. Bridging fibers that have fractured close to
a crack will pull out from the matrix rather than fracture again (Warren 1992). The
greater the fiber pull-out length, the greater the energy absorbed by the composite.
A zone of microcracks may form ahead of a growing crack. This area of
microcracks contributes to an overall toughness increase by forming a zone of lower
elastic modulus and absorbing strain release energy. To be effective, microcrack zones
should be restricted to individual sites in the composite in order to avoid microcrack
linkage (Warren 1992).
Characteristics of the bond between the fibers and the matrix also determine the
mechanical behavior of composite materials. If the interfacial bond between the matrix
and the fiber is sufficiently strong, a crack will propagate through a fiber with little
obstruction rather than traveling around it (Wachtman 1989). Debonding of the fiber and
matrix is an energy-dissipating mechanism that complements processes like crack
deflection and fiber pullout, leading to an overall toughness increase in the composite.
7
2.2 Specific Types of Ceramic Matrix Composite Materials; Applications
Structural ceramics are high performance ceramics such as oxides, nitrides, and
carbides of silicon, aluminum, titanium, and zirconium (Chawla 2003). Tables 2.1, 2.2,
and 2.3 list properties of materials that are well-suited for ceramic matrix materials,
important ceramic reinforcements, and examples of ceramic fiber/ceramic matrix
composites, respectively. The most common woven CMCs use carbon (C) fibers or
silicon carbide (SiC) fibers. Carbon fibers greatly improve a composite?s toughness, but
their capacity is limited at high temperatures and by corrosive or oxidizing environments.
Silicon carbide fibers are more resistant to oxidation and can be used at higher
temperatures than carbon fibers (Wachtman 1989). The woven fiber-matrix
combinations that have been investigated the most are C-C, C-SiC, and SiC-SiC.
Table 2.1 Properties of some high performance ceramics (Chawla 2003).
Material Young's
Modulus E
(GPa)
Poisson's Ratio
?
Coefficient of
Thermal Expansion
? (10-6 /oK)
SiC 420 0.22 4
Al2O3 380 0.25 8
Cordierite 130 0.25 2
Mullite 215 0.25 4
Sodalime glass 70 0.23 9
Table 2.2 Some important ceramic reinforcements (Chawla 2003).
Particle Sic, TiC, Al2O3
Discontinuous Fibers
(a) Whiskers Sic, TiB2, Al2O3
(b) Short Fibers Glass, Al2O3, SiC, (Al2O3+SiO2),
vapor grown carbon fibers
Continuous Fibers
(a) Oxide Al2O3, (Al2O3+SiO2), ZrO2,
silica-based glasses
(b) Nonoxide B, C, SiC, Si3N4, BN
8
Table 2.3 Examples of ceramic/ceramic composites (Warren 1992).
Composite type Matrix-Reinforcement
Particulate Al2O3-ZrO2
Al2O3-TiC
Al2O3-SiC
SiC-TiB2
Si3N4-TiC
Si3N4-ZrO2
Platelets Al2O3-SiC
Si3N4-SiC
Short random fibers (Whiskers) Al2O3-SiC
Si3N4-SiC
Long, parallel fibers Glass-C
Glass-SiC
Cross-plied Glass-C
Glass-SiC
SiC-SiC
Woven C-C
SiC-SiC
CMCs are most valued for their high temperature strength and performance
capabilities as well as for their wear resistance (Warren 1992), and they are studied for
applications in spacecraft programs as well as for inclusion in fighter planes, missiles,
and rockets (Chawla 2003). Low density, high strength composites are very often used in
numerous aerospace applications in which high strength-to-weight ratios must be
achieved. Non-aerospace applications include cutting tool inserts, wear parts in
machinery, nozzles, valve seals, and bearings. Other uses of ceramic materials include
non-structural applications, such as thermal insulation (Warren 1992).
2.3 Mechanical Behavior
A composite material may be composed of homogeneous constituents, but it is
macroscopically heterogeneous. The material properties of a homogeneous and isotropic
material are independent of both position in the body and orientation at a point in the
9
body. Most composite materials are heterogeneous and anisotropic, or nonisotropic.
The material properties depend on position and orientation of the fibers in the body and
on orientation with respect to load application (Jones 1999).
Stress (?) and strain (? ) components at a point in a body can be represented by
second-rank tensors, as follows:
?
?
?
?
?
?
?
?
?
?
=
333231
232221
131211
???
???
???
? , (2.1)
?
?
?
?
?
?
?
?
?
?
=
333231
232221
131211
???
???
???
? . (2.2)
For composite materials, the directions 1, 2, and 3 represented by the subscripts
correspond to the Cartesian coordinate directions x, y, and z. From Hooke?s law, each
stress component is proportional to each strain component. In equation form, the
constitutive relationships between stress and strain are expressed as:
Cij ijkl kl? ?= , (2.3)
Sij ijkl kl? ?= , (2.4)
where i, j, k, l = 1, 2, or 3,
Cijkl is the matrix of elastic stiffnesses, and
Sijkl is the matrix of elastic compliances.
Both C and S are fourth-rank tensors, and, as the equations stand, each has eighty-one
terms.
10
Due to the symmetry of the stress and strain components ( ,ij ji ij ji? ? ? ?= = ), the
maximum number of independent stiffnesses and compliances is reduced to thirty-six.
Through energy arguments, the number is further reduced to twenty-one. Thus, in four
subscript notation for both Cijkl and Sijkl , ij = ji, kl = lk, and ijkl = klij. It is more
convenient to use contracted notation:
Cq qr r? ?= , (2.5)
Sq qr r? ?= , (2.6)
where q, r = 1, 2, 3, 4, 5, or 6 and correspond to ij, kl = 11, 22, 33, 23, 13, or 12.
In expanded form:
1 11 12 13 14 15 16 1
2 22 23 24 25 26 2
3 33 34 35 36 3
4 44 45 46 4
5 55 56 5
6 66 6
C C C C C C
C C C C C
C C C C
C C C
Sym C C
C
? ?
? ?
? ?
? ?
? ?
? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?? ? ? ?=
? ?? ? ? ?
? ?? ? ? ?
? ?? ? ? ?
? ?? ? ? ?? ?? ? ? ?
? ? ? ? ? ?
, (2.7)
??
?
?
?
??
?
?
?
?
??
?
?
?
??
?
?
?
?
??
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
=
??
?
?
?
??
?
?
?
?
??
?
?
?
??
?
?
?
?
6
5
4
3
2
1
66
5655
464544
36353433
2625242322
161514131211
6
5
4
3
2
1
?
?
?
?
?
?
?
?
?
?
?
?
S
SSSym
SSS
SSSS
SSSSS
SSSSSS
. (2.8)
The physical interpretation of the components of the compliance matrix is as follows:
square6 S11, S22, and S33 relate an extensional stress in one coordinate direction to an
extensional strain in the same coordinate direction, and Sqq = 1/Eq, where Eq is the
Young?s Modulus in the q direction.
11
square6 S12, S13, and S23 relate an extensional strain to a perpendicular extensional stress; this
relationship is called extension-extension coupling, and Sqr = -?rq / Er = -?qr / Eq, where
?qr is Poisson?s ratio in the qr plane.
square6 S15, S16, S24, S26, S34, and S35 relate an extensional strain to a shear stress in the same
plane, also called shear-extension coupling.
square6 S14, S25, and S36 relate an extensional strain to a shear stress in a perpendicular plane;
also called shear-extension coupling.
square6 S44, S55, and S66 relate a shear strain to a shear stress in the same plane, and Sqq = 1/Gq,
where Gq is the shear modulus in the q plane.
square6 S45, S46, and S56 relate a shear strain to a shear stress in a perpendicular plane, also
called shear-shear coupling (Hearmon 1961).
2.4 Types of Material Symmetry
As previously stated, the maximum number of independent elastic constants is
twenty-one. This number is reduced if the material has symmetry elements. There are
thirty-two classes of symmetry based on three types of symmetry: planes of reflection
symmetry, pure rotation axes, and axes of rotary inversion. These operations are not
necessarily independent, and only nine different symmetry systems result from the thirty-
two classes of symmetry (Hearmon 1961). The values for the elastic stiffnesses for
material symmetry systems are summarized in Tables 2.4 and 2.5. The first row gives the
name of the symmetry system, the second row gives the number of independent elastic
stiffnesses, and the remaining rows list the notations of the elastic stiffnesses.
12
The types of symmetry most common for composite materials are triclinic (no
symmetry elements), monoclinic (one plane of reflection), orthotropic (three mutually
orthogonal planes of reflection), transversely isotropic (one plane of isotropy), and
isotropic (all planes are planes of reflection). A plane of reflection converts a point to its
mirror image. For example, if the 1-2 plane is a plane of reflection (also called a plane of
material symmetry) the point (x1, x2, x3) is converted to (x1, x2, -x3). For a plane of
isotropy, all planes perpendicular to that plane are planes of reflection. The other
symmetry elements and symmetry systems are exhibited in certain types of crystals
(Suhling 2004).
Table 2.4 Notations of elastic stiffnesses of the triclinic, monoclinic, orthotropic,
tetragonal, and cubic symmetry systems (Hearmon 1961).
Triclinic Monoclinic Orthotropic
(Orthorhombic)
Tetragonal Cubic
21 13 13 13 9 7 6 3
C11 C11 C11 C11 C11 C11 C11 C11
C12 C12 C12 C12 C12 C12 C12 C12
C13 C13 C13 C13 C13 C13 C13 C12
C14 C14 0 0 0 0 0 0
C15 0 C15 0 0 0 0 0
C16 0 0 C16 0 C16 0 0
C22 C22 C22 C22 C22 C11 C11 C11
C23 C23 C23 C23 C23 C13 C13 C12
C24 C24 0 0 0 0 0 0
C25 0 C25 0 0 0 0 0
C26 0 0 C26 0 -C16 0 0
C33 C33 C33 C33 C33 C33 C33 C11
C34 C34 0 0 0 0 0 0
C35 0 C35 0 0 0 0 0
C36 0 0 C36 0 0 0 0
C44 C44 C44 C44 C44 C44 C44 C44
C45 0 0 C45 0 0 0 0
C46 0 C46 0 0 0 0 0
C55 C55 C55 C55 C55 C44 C44 C44
C56 C56 0 0 0 0 0 0
C66 C66 C66 C66 C66 C66 C66 C44
13
Table 2.5 Notations of elastic stiffnesses of the triclinic, trigonal, transversely isotropic,
and isotropic symmetry systems (Hearmon 1961).
Triclinic Trigonal Transversely
Isotropic
(Hexagonal)
Isotropic
21 7 6 5 2
C11 C11 C11 C11 C11
C12 C12 C12 C12 C12
C13 C13 C13 C13 C12
C14 C14 C14 0 0
C15 C15 0 0 0
C16 0 0 0 0
C22 C11 C11 C11 C11
C23 C13 C13 C13 C12
C24 -C14 -C14 0 0
C25 -C15 0 0 0
C26 0 0 0 0
C33 C33 C33 C33 C11
C34 0 0 0 0
C35 0 0 0 0
C36 0 0 0 0
C44 C44 C44 C44 ? (C11-C12)
C45 0 0 0 0
C46 -C15 0 0 0
C55 C44 C44 C44 ? (C11-C12)
C56 C14 0 0 0
C66 ? (C11-C12) ? (C11-C12) ? (C11-C12) ? (C11-C12)
2.5 Rotation of Axes
When composite materials are subjected to loads in service conditions, the
direction of loading does not always coincide with the principal material directions (1, 2,
3). In this case, it is necessary to express the stress-strain relations using a new set of
axes, (x?, y?, z?) (Suhling 2004). The stiffnesses and the compliances are both fourth-
rank tensors, and on transforming from one set of axes to another:
C' Cijkl im jn ko lp mnopa a a a= , (2.9)
S' Sijkl im jn ko lp mnopa a a a= , (2.10)
14
where a represents a matrix of direction cosines relating one set of axes to another.
As previously stated, the maximum number of independent elastic constants is
twenty-one. When transforming from one set of axes to another, there will be twenty-one
equations, each containing twenty-one terms. When the axes undergo a rotation, the
number of independent constants remains the same, but the number of terms in the
stiffness and compliance tensors in the new coordinate system (x?, y?, z?) may increase
(Hearmon 1961).
2.6 Micromechanics and Macromechanics
Two basic approaches are used to determine the properties of a composite
material: the micromechanical approach and the macromechanical approach. The goal
of the micromechanical approach is to determine the material properties of a composite in
terms of the properties of its constituents. Micromechanics can be used to predict the
composite?s stiffnesses and compliances and is used when designing a composite.
Macromechanics assumes that the composite is homogeneous, and its material behavior
is based on the average apparent mechanical properties of the composite. The
macromechanical behavior of a composite is experimentally determined using the
composite material as a whole, rather than by testing each constituent, and is used in
design of components utilizing composite materials (Jones 1999).
2.7 Testing Methods
Standard material test procedures, such as tensile, flexure, and torsion tests, can
be used to determine elastic constants of composite materials, but calculation methods
may need to be modified due to the anisotropy of the material (Hearmon 1961) and the
further anisotropy induced as the composite experiences damage. Additionally,
15
fabrication of test samples for the various tests is not always easy due to the hardness of
many composites and due to the difficulties inherent in manufacturing composites at
different orientations. Furthermore, because these tests are destructive, the same sample
cannot be used for all tests, introducing the possibility of manufacturing irregularities
among specimens. Non-destructive techniques, such as vibrational and ultrasonic
methods, offer the capability to determine all the elastic constants for the same sample
quickly (Rokhlin and Wang 1992). Ultrasonic wave methods have been used extensively
to measure the elastic stiffnesses of anisotropic materials. Test setups that use ultrasonic
pulses include single-through transmission, double through-transmission, and point-
source/point-receiver techniques.
2.8 Ultrasonics
By combining the stress-strain relationship,
klijklij ?? C= , (2.11)
and the definition of strain,
( )ijjiij uu ,,21 +=? , (2.12)
with the linear momentum balance equation,
iijij uf &&??? =+, , (2.13)
one can obtain the following relationship
( ) ikjlljkijkl uuuC &&?=+ ,,21 , (2.14)
where ? is the mass density,
fi are the body forces (assumed to be zero),
u is the displacement vector,
16
and the double dot notation denotes differentiation twice with respect to time.
Combining Equation 2.14 with the equation for a plane harmonic displacement wave:
( )( )tqxidAu ppmm ?= ?exp , (2.15)
gives:
( ) 0=? kikljijkl dqqC ?? , (2.16)
where A is an independent constant,
d is the unit vector defining direction of motion of the particle displaced by the
plane harmonic wave,
i is the imaginary number,
? is the real-valued angular frequency of the plane harmonic wave,
x is the position vector,
q is the slowness vector (defined below),
t is the time variable,
? is Kronecker Delta, and
m, p = 1, 2, or 3.
For a nontrivial solution:
.0det =? ikljijkl qqC ?? (2.17)
The components of the slowness vector are defined as:
j jq n V= , (2.18)
where nj are components of the vector of the direction of wave propagation, and V is the
phase velocity.
17
Equation 2.17 can be rewritten using what are commonly referred to as the
Christoffel stiffnesses:
ik ijkl j lC n n? = , (2.19)
so that:
2det 0ik ikV? ?? ? = , (2.20)
Equation 2.20 is known as the Christoffel equation. It can be expressed in expanded
form as:
0
2
332313
23
2
2212
1312
2
11
=
????
????
????
V
V
V
?
?
?
, (2.21)
where
32563115211623552266211111 C2C2C2CCC nnnnnnnnn +++++=? , (2.22)
32243146212623442222216622 C2C2C2CCC nnnnnnnnn +++++=? , (2.23)
32343135214523332244215533 C2C2C2CCC nnnnnnnnn +++++=? , (2.24)
( ) ( ) ( ) 32254631561421661223452226211612 CCCCCCCCC nnnnnnnnn ++++++++=? , (2.25)
( ) ( ) ( ) 32453631551321561423352246211513 CCCCCCCCC nnnnnnnnn ++++++++=? , (2.26)
( ) ( ) ( ) 32442331453621254623342224215623 CCCCCCCCC nnnnnnnnn ++++++++=? , (2.27)
.,, 233213311221 ?=??=??=?
All of the eigenvalues of ?ik are real and positive and their corresponding
eigenvectors are orthogonal. The physical interpretation of this is that for a given
direction of wave propagation there will be three phase velocities and the three
18
corresponding displacement vectors will be orthogonal. In an anisotropic case, the
displacements are neither truly longitudinal nor truly transverse in character (Achenbach
1973).
19
CHAPTER III
DETERMINATION OF ELASTIC CONSTANTS AND DAMAGE IN CERAMIC
MATRIX COMPOSITES USING ULTRASONIC WAVE SPEED
MEASUREMENTS
3.1 Ultrasonic Test Methods
Common ultrasonic pulse methods include point source/point receiver and
immersion methods. In the point source/point receiver setup, the specimen is in direct
contact with the transducers. The wave velocities can be measured only in a direction
normal to the face of the sample. Specimens are required to be cut in several orientations
to compute the full set of elastic constants. This method is suitable for crystal but is not
feasible for many composite materials. Immersion methods allow the measurement of
wave velocities over a range of directions. In the one-way transmission setup, the
specimen is placed between two ultrasonic transducers, one for transmission and one for
reception (Markham 1970). The double through-transmission method uses a single
transducer that functions as both transmitter and receiver. The signal passes through the
sample, reflects off a back plate reflector, and returns to the transmitter/receiver (Rokhlin
and Wang 1992). The work in this thesis is based on the one-way through-transmission
test procedure.
20
Figure 3.1 shows a general view of an ultrasonic immersion tank. For this
method, the specimen is placed on a turntable in a tank filled with a liquid, usually water.
Figure 3.2 shows a schematic of the specimen on the turntable between two ultrasonic
transducers. The turntable can be rotated in two directions corresponding to the angles ?
and ? defined in the figure (Markham 1970).
Figure 3.1 General test setup for an ultrasonic immersion test (Markham 1970).
Figure 3.2 Rotation of sample in order to excite transverse waves. N is the axis normal to
the specimen; i is the incidence angle (the angle between the incoming wave and N); ? is
the degree of tilt; ? is the degree of spin (Markham 1970).
21
Audoin and Baste (1994) combined the ultrasonic immersion technique with a
uniaxial tension test. Figure 3.3 shows a schematic of a composite material test sample
positioned between the transmitting and receiving transducers, as uniaxial tension is
being applied. The sample is loaded in the 3 direction, and the transducers can rotate in
two directions (about the 1 and 2 axes).
The full set of elastic constants can be determined for each level of tensile load
applied, and the effects of damage on the sample can be ascertained by noting the
variation in the values of the elastic constants as damage progresses (Audoin and Baste
1994). In order that true plane waves are produced in the sample, the wavelength of the
ultrasonic transducer should be large compared to the fiber diameter and small compared
to the dimension of the specimen (Markham 1970).
22
Figure 3.3 Under-load ultrasonic device. The sample is loaded along direction 3. The
transducers are moved by a combination of two rotations, ?1 and ?2 (Audoin and Baste
1994).
3.2 The Time Difference between the Group and Phase Velocities
The phase velocity is needed to compute the elastic constants, as indicated in the
Christoffel equation (Equation 2.20), but the time measurements used in early immersion
tests were actually computed for the group velocity. In anisotropic materials, these two
velocities generally deviate from each other. This deviation called into question the
validity of early immersion test results. However, Rokhlin and Wang (1992) show that
23
there is no time difference between the group velocity and the phase velocity for an
arbitrary angle of incidence. In measurements of this type, T0 is the time it takes for the
signal to travel from the transmitter to receiver through the reference medium, e.g. water,
without a test sample in place. Tp is the total time for the signal to reach the receiver
when passing through the specimen at a particular phase velocity. Tg is the total time it
takes for the signal to reach the receiver when passing through the specimen at a
particular group velocity. The times for the group and phase velocity vectors in the
sample are tg and tp, respectively. Referring to Figure 3.4, the times can be written as:
,43210 ttttT +++= (3.1)
,41 tttT gg ++= (3.2)
1 3 4.p pT t t t t= + + + (3.3)
In this figure, time progresses from the top of the figure to the bottom. Boxes showing
"T" and "R" represent the locations of the transmitting transducer and receiving
transducer, respectively, and h indicates the thickness of the sample. The left-most
vertical line indicates the times associated with Tg (the path taken by the wave associated
with the group velocity is shown by the path T-O-A-R), the vertical line near the middle
of the figure indicates the times associated with Tp (the path taken by the wave associated
with the phase velocity is shown by the path T-O-B-C-R), and the right-most vertical line
indicates the times associated with T0 (the reference path does not pass through the
sample). The times associated with the group and phase velocities can be written as the
appropriate path length divided by the corresponding velocity (V0, Vg, and Vp are the
24
reference velocity (velocity with no specimen present), group velocity, and phase
velocity, respectively).
Figure 3.4 Representation of the times for different acoustic paths. T0 is the reference
time; Tp is the time for the phase velocity vector; Tg is the time for the group velocity
vector.
Figure 3.5 gives an alternate view of the test paths shown in Figure 3.4. It should
be noted that the phase velocity vector lies in the incident plane (the plane defined by the
incoming wave and the axis normal to the sample), but the group velocity vector does not
necessarily lie in this plane. In Figure 3.5, OA is the path of the group velocity vector
through the sample, and OB is the path of the phase velocity vector through the sample.
The angle ? is the angle of deviation between the group and phase velocity vectors. The
25
angle ? is the in-incident-plane component of ?, ? is the out-of-plane component, and
??? coscoscos = .
Figure 3.5 Schematic of different acoustic paths for the group and phase velocity
vectors. h is the sample thickness; ? is the angle of deviation between the group and
phase velocity vectors; ? is the in-incident-plane component of ?; ? is the out-of-
plane component of ? ; i? is the incident angle; r? is the refraction angle (Rokhlin
and Wang 1992).
Considering Figures 3.4 and 3.5, the times specified can be rewritten as:
gg VOAt /= , (3.4)
pp VOBt /= , (3.5)
where )cos()cos( ???
r
hOA
+= , and
)cos( r
hOB
?= .
26
The quantity ( )st? represents the difference between the group velocity and the phase
velocity in the specimen:
pgs ttt ?=)(? . (3.6)
( )st? can be more conveniently expressed through the following development:
Snell?s Law (shown below as Equation 3.7) relates the angle of incidence ( i? ),
the angle of refraction ( r? ), and the velocities outside and inside an interface of two
substances ( pV and oV here). In this case, the angle of refraction corresponds to the
phase velocity, so:
o
i
p
r
VV
)sin()sin( ?? = . (3.7)
Using Snell?s Law and the relationship between the group and phase velocities:
( ) )cos()cos(cos ??? ggp VVV == , (3.8)
the time difference between the group and phase velocities in the specimen, ( )st? , can be
rewritten as:
( ) ( )
)cos()cos(
)sin()sin(
)cos()cos(
)sin()sin()cos()cos()cos()cos(
)sin(
)sin(
)cos(
1
)cos(
)cos(
)sin(
)sin(
)cos()cos()cos(
coscos
)cos()cos()cos()(
rro
i
rr
rrr
ro
i
rrro
i
rprp
rprg
s
V
h
V
h
V
h
V
h
V
h
V
h
V
ht
???
??
???
??????
?
?
???
?
?
?
????
??
?????
+=
???
?
???
?
+
+?=
???
?
???
? ?
+=
?+=
?+=
(3.9)
Solving for t3 (as indicated in Figure 3.4) gives:
27
oVBCt /3 = . (3.10)
BC can be written as the difference between the vertical components of the group
velocity path, OA , and the phase velocity path, OB (both shown in Figure 3.5):
cos( ) cos( )
cos( ) cos( )
cos( ) cos( ) sin( )sin( ) cos( ) cos( ) sin( )sin( )
cos( ) cos( ) cos( ) cos( )
sin( )sin( )cos( ) sin( )sin( )cos( )
cos(
r i r i
r r
r i r i r i r i
r r r r
r i r r i r
h hBC
h
h
? ? ? ? ?
? ? ?
? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ? ? ?
+ ? ?= ?
+
? ?+ += + ? ?
? ?+ +? ?
+ ? +=
( )
( )
2 2
)cos( )
sin( )cos( ) cos( )sin( ) sin( )cos( )
cos( )cos( )
sin( )sin( ) cos( ) cos( ) sin( )sin( )
cos( ) cos( )
sin( )sin( )cos ( ) sin( )sin( )sin ( )
cos(
r r
r r i r
r r
r i r r
r r
i r i r
r
h
h
? ? ?
? ? ? ? ? ?
? ? ?
? ? ? ? ? ?
? ? ?
? ? ? ? ? ?
?
? ?
? ?+? ?
? ?+ ?
? ?+
? ?=
? ??
? ?? ?+
? ?
+=
) cos( )
sin( )sin( )
cos( ) cos( )
r
i
r r
h
? ?
? ?
? ? ?
? ?
? ?+? ?
= +
(3.11)
Thus:
)cos()cos( )sin()sin(3
rro
i
V
ht
???
??
+= , (3.12)
and therefore:
stt )(3 ?= . (3.13)
This result verifies that there is no total time difference between the times associated with
the group and phase velocities, as shown below in Equation 3.14.
( )
( )
1 3 4 1 4
1 4 1 4 .
p p p s
p p g p g g
T t t t t t t t t
T t t t t t t t t T
?= + + + = + + +
= + + ? + = + + = (3.14)
28
3.3 Calculating the Phase Velocity and Refraction Angle
To calculate the phase velocity, knowledge of the signal time delay and the length
of the acoustic path is required. The time difference between the reference time, T0, and
Tp is denoted by t? and is given by:
ppo ttTTt ?=?=? 2 , (3.15)
where )cos( )cos(2
ro
ir
V
ht
?
?? ?= , and
)cos( rpp V
ht
?= .
Substituting these values into Equation 3.15 and using the Snell's Law relationship given
in Equation 3.7 yields:
p
r
o
i
rp
r
o
i
rpro
ir
o
i
rpro
irir
VV
VV
VVV
VVh
t
)cos()cos(
)cos(
1)(sin)cos(
)cos(
1
)cos(
)sin()sin()cos(
)cos(
1
)cos(
)sin()sin()cos()cos(
2
??
?
??
??
???
??
????
?=
?+=
?+=
?+=?
(3.16)
Using this result yields:
htVV
o
i
p
r ??= )cos()cos( ?? , (3.17)
and noting that Snell?s Law is expressed as:
o
i
p
r
VV
)sin()sin( ?? = , (3.18)
29
a formulation for the refraction angle, r? , can be developed using the incidence angle, i? ,
time difference, t? , reference velocity, V0, and sample thickness, h. Dividing Equation
3.18 by Equation 3.17 gives:
h
Vt
h
t
VV
V
V
o
i
i
o
i
p
o
ip
r ?
?
=
???
?
???
? ??= )cos(
)sin(
)cos(
)sin(
)tan(
?
?
?
?
? , (3.19)
so that:
??
?
?
?
?
??
?
?
?
?
??=
h
Vt o
i
i
r
)cos(
)sin(arctan
?
?? . (3.20)
(This expression for the refraction angle includes quantities that are all known or can be
measured.) The phase velocity is recovered by substituting this result into Snell?s Law:
)sin( )sin(
i
ro
p
VV
?
?= . (3.21)
It is also possible to square and add Equations 3.17 and 3.18:
2222 )sin()cos()sin()cos(
???
?
???
?+
???
?
???
? ??=
???
?
???
?+
???
?
???
?
o
i
o
i
p
r
p
r
Vh
t
VVV
???? , (3.22)
which simplifies to:
( ) ??
?
?
???
? ?+??=
2
2
22
cos211
h
t
hV
t
VV o
i
op
? . (3.23)
Consequently:
( )
2/1
2
2
2
cos21 ?
???
?
???
? ?+??=
h
t
hV
t
VV o
i
o
p
? , (3.24)
30
and the refraction angle, r? , is recovered by substituting this result into Snell?s Law:
( )??
?
?
???
?=
o
ip
r V
V ?? sinarcsin . (3.25)
3.4 Wave Behavior
By changing the angle of incidence of the ultrasonic pulse (via rotation of the
sample as shown in Figure 3.2 or via rotation of the transducers as shown in Figure 3.3),
the wave will split by mode conversion on entering the solid, as indicated in Figure 3.6.
If the incident direction of the wave coincides with an axis of symmetry, only the
longitudinal wave velocity is transmitted through the sample. If the incident wave is
introduced in a plane of symmetry, or principal plane, the incoming wave splits into two
components: one quasi-transverse and one quasi-longitudinal. If it is introduced in a
non-principal plane, the wave splits into two quasi-transverse modes and a longitudinal
mode; these are known as the fast shear wave (QT1), the slow shear wave (QT2), and a
quasi-longitudinal wave (QL). Each of the three waves travels at a different speed, and
each arrives at the receiver at a different time (Markham 1970).
31
Figure 3.6 Schematic diagram of waves excited in a non-principal plane. QL stand for
quasi-longitudinal mode; QT1 and QT2 stand for fast and slow quasi-transverse modes,
respectively; the incident plane is defined by the angle ?; ?i is the incident angle (the
angle between x1 and the incident wave) (Aristegui and Baste 1997).
Often, the composite materials to be tested are very thin in the 1 direction, and
direct measurements in the 2-3 plane are unavailable (Chu and Rokhlin 1992).
Measurements are limited to planes perpendicular to the plane of the sample.
Most ceramic matrix composites are of an orthotropic or higher order of
symmetry, i.e. they possess fewer than twenty-one independent elastic constants.
Orthotropic materials have three planes of symmetry and nine independent elastic
constants. For measurements in the 1-2 plane, which is a plane of symmetry, the phase
32
velocity vector, n, can be represented as (n1, n2, 0). The Christoffel equation (Equation
2.20) reduces to:
( )
( ) .0
00
0
0
22
244
2
155
22
222
2
166216612
216612
22
266
2
111
=
?+
?++
+?+
VnCnC
VnCnCnnCC
nnCCVnCnC
?
?
?
(3.26)
The 33? term, 22442155 nCnC + , corresponds to transverse waves polarized in the 1-3 plane,
which are not produced when scanning in the 1-2 plane.
For measurements in the 1-3 plane, the second accessible plane of symmetry, the
wave vector is (n1, 0, n3), and the Christoffel equation reduces to:
( )
( )
.0
0
00
0
22
333
2
155315513
22
344
2
166
315513
22
355
2
111
=
?++
?+
+?+
VnCnCnnCC
VnCnC
nnCCVnCnC
?
?
?
(3.27)
The 22? , 23442166 nCnC + , term corresponds to transverse waves polarized in the 1-2 plane,
which are not produced when scanning in the 1-3 plane.
As a result, only seven of the nine independent constants, C11, C22, C33, C55, C66,
C12, and C13, can be determined from measurements in the two accessible principal
planes for an orthotropic material. The remaining constants, C23 and C44, must be
recovered from measurements taken in a non-principal plane. Often this plane is referred
to as 1-2?, where 2? is at an angle between the 2 and 3 axes, usually 45 degrees. Figure
3.7 shows the incoming wave vectors in the three data planes (Mouchtachi et al. 2004).
Complications may arise for woven composites. In this case, if the 1-2? plane is
oriented at 45 degrees, the plane may also be a plane of symmetry. Only one transverse
wave mode will be produced. This affects the results of C23 and C44, two constants in the
33
2-3 plane which are not directly measurable. Measurements in a non-symmetry plane
may prove to be more useful (Aristegui and Baste 1997).
Figure 3.7 Schematic of the composite plate. Axis 1 is normal to the plate; the angles
defining the direction of the wave are ? in plane (1, 2), ? in plane (1, 2?) and ? in plane
(1, 3). Plane (1, 2?) is an arbitrary non-principal plane defined by angle ? (Mouchtachi
et al. 2004).
3.5 The Wave Vector
In order to calculate the elastic constants, knowledge of not only the velocity but
also the direction of propagation of the wave in the sample is required. Any rotation
relating two sets of axes may be described using three rotations (Weisstein). The
34
rotations may be written individually as rotation matrices B, C, and D, and the final
rotation matrix A is their product:
A = B C D, (3.28)
?
?
?
?
?
?
?
?
?
?
=
333231
232221
131211
A
aaa
aaa
aaa
. (3.29)
There are many combinations of rotations that may be utilized to reach the same final
position. For this study, a choice was made to use the ?x y z? (pitch-roll-yaw)
convention. D is a rotation around the z-axis, C is a rotation around the y-axis, and B is a
rotation around the x-axis (Weisstein). Not all of these rotations are required for
ultrasonic immersion tests, but they are possible configurations. In this case, tilt is the
rotation around the 3 axis, spin is the rotation around the 2 axis, and the in-plane (the
plane of the sample) rotation is around the 1 axis, as pictured in Figures 3.2 and 3.3. In
matrix form:
( ) ( )
( ) ( )
?
?
?
?
?
?
?
?
?
?
?=
100
0cossin
0sincos
tilt ??
??
, (3.30)
( ) ( )
( ) ( ) ??
?
?
?
?
?
?
?
? ?
=
??
??
cos0sin
010
sin0cos
spin , (3.31)
( ) ( )
( ) ( )??
?
?
?
?
?
?
?
?
?
=
??
??
cossin0
sincos0
001
rotation plane-in . (3.32)
35
(Though some of the variable nomenclature is the same, the angles indicated here are not
the same angles as those mentioned in Section 3.3 to describe determination of the phase
velocity and refraction angle).
The product of these three rotation matrices yields a matrix which relates the final
orientation of the sample to its starting position. If the (1, 2, 3) directions correspond to
the coordinate system of the sample in its final position, and the (x, y, z) directions
correspond to the coordinate system of the sample in its initial position, then:
( )( )( )
??
???
??
???=
??
???
??
????=
??
???
??
???
z
y
x
z
y
x
A
tiltspinplanein
3
2
1
, (3.33)
where
?
?
?
?
?
?
?
?
?
?
=
333231
232221
131211
A
aaa
aaa
aaa
and
11
12
13
21
22
23
31
32
33
cos cos ,
cos sin ,
sin ,
cos sin cos sin sin ,
cos cos sin sin sin ,
cos sin ,
sin sin cos cos sin ,
cos sin cos sin sin ,
cos cos .
a
a
a
a
a
a
a
a
a
? ?
? ?
?
? ? ? ? ?
? ? ? ? ?
? ?
? ? ? ? ?
? ? ? ? ?
? ?
=
=
= ?
= ? +
= +
=
= +
= ? +
=
36
The a11 term shown above is equal to the cosine of the angle between the 1 axis
(normal to the plane of the sample) and the x-axis, which is the direction of the incoming
wave. This angle is the incidence angle, i? , so:
??? coscoscos =i . (3.34)
The refracted wave lies in the incident plane. This plane is defined by two
vectors: the incoming wave and the plate normal, or the x-axis and the 1-axis. The cross
product of these two vectors gives a vector normal to the incident plane, as shown below:
( )???
?????
sincos,sin,0
sinsincoscoscos
0011 =
?
=?
kji
x . (3.35)
Now, a coordinate system (x?, y?, z?) is defined for the refracted wave. z? is the
direction of the refracted wave, x? lies in the incident plane, and y? is normal to the
incident plane. A rotation matrix relating these two coordinate systems is required. As
stated before, the vector ( )??? sincos,sin,0 is normal to the incident plane. Rotation
around the x axis allows y? to be defined as this vector. After this rotation, x? and z? will
be in the incident plane, and the coordinate system may be rotated about the y? axis until
z? is in the direction of wave propagation. The general form for a rotation matrix around
the x-axis is:
( ) ( )
( ) ( )
1 0 0
0 cos sin
0 sin cos
? ?
? ?
? ?
? ??? ?
(3.36)
37
The direction cosines for the vector ( )??? sincos,sin,0 must be calculated. This
is accomplished by dividing the dot product of the vector and the appropriate unit vector
by the length of the vector as follows:
( ) ( )( ) 0sincos,sin,0 0,0,1sincos,sin,0 =???? ??? (3.37)
( ) ( )( )
???
?
???
???
222 sincossin
sin
sincos,sin,0
0,1,0sincos,sin,0
+
=? (3.38)
( ) ( )( )
???
??
???
???
222 sincossin
sincos
sincos,sin,0
1,0,0sincos,sin,0
+
=? (3.39)
The length of the vector can be simplified in the following manner:
( )
.sin
cos1
coscos1
cos1cossinsincossin
2
22
222222
i
i
?
?
??
??????
=
?=
?=
?+=+
(3.40)
The rotation matrix for rotation around the x-axis then becomes:
1
1 0 0
sin cos sin0
sin sin
cos sin sin0
sin sin
i i
i i
r ? ? ?? ?
? ? ?
? ?
? ?
? ?
? ?
? ?=
? ?
? ?
? ??
? ?
? ?
. (3.41)
After this rotation, the y? axis is normal to the incident plane and the z? axis is in
the incident plane. An appropriate rotation around the y? axis will align the z? axis with
the refracted wave. Figure 3.8 gives the perspective shown by looking at the incident
plane with the y? axis out of the page:
38
Figure 3.8 View of the incident plane. Axis 1 is normal to the plate;
Axis 3 is the long axis of the plate; i? is the angle of incidence;
r? is the angle of refraction.
The angle of rotation necessary to align the z' axis with the refracted wave is
ri ??
pi ?+
2
3 .
The second rotation matrix is thus determined to be:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????? ?+?????? ?+
?????? ?+??????? ?+
=
riri
riri
r
??pi??pi
??pi??pi
2
3cos0
2
3sin
010
2
3sin0
2
3cos
2 (3.42)
The relationship between the (x, y, z) and (x?, y?, z?) coordinate systems is
defined as:
( )( )
??
???
??
???=
??
???
??
???
z
y
x
rr
z
y
x
12
'
'
'
(3.43)
39
The matrix A mentioned previously related the (1, 2, 3) and (x, y, z) coordinate systems
through Equation 3.33. This relationship can be used to develop the relationship between
the (x', y', z') coordinate system and the (1, 2, 3) coordinate systems as follows:
The inverse of a rotation matrix is equal to its transpose. Therefore:
( ) ( ) ( )
??
???
??
????=
??
???
??
???=
??
???
??
???
3
2
1
planeinspintilt
3
2
1
A TTTT
z
y
x
. (3.44)
Substituting this relationship into Equation 3.43 gives:
( )( )( ) ( ) ( )T T T2 1
' 1
' tilt spin in plane 2
' 3
' 1
' 2 .
' 3
x
y r r
z
x
y R
z
? ? ? ?? ? ? ?
= ?? ? ? ?
? ? ? ?? ? ? ?
? ? ? ?? ? ? ?
=? ? ? ?
? ? ? ?? ? ? ?
(3.45)
Equation 3.45, then, defines a rotation matrix, R, relating the coordinate system of the
sample with that of the refracted wave using input values of degrees of spin, tilt, and in-
plane rotation, and the calculated values of angle of incidence and angle of refraction.
The components of the matrix R are as follows:
( ) ( )( ) ???? ????????? 2
22
11
sincsccos
sincoscossincsccoscos
iir
iriirR
??
????= (3.46)
( )
( )( )
( ) ( )??????????
???????
??????
sinsinsincoscossincsccoscos
sinsincossincossin
sinsincsccoscos12
+??
+???
?=
iir
ir
iirR
(3.47)
( )
( )( )
( ) ( )??????????
???????
??????
sincossinsincossincsccoscos
sincoscossinsinsin
sincsccoscoscos13
+???
+??
?=
iir
ir
iirR
(3.48)
40
021 =R (3.46)
( )??????????? sinsinsincoscossincscsinsincsccos222 ++= iiR (3.49)
( )??????????? sincossinsincossincsccossincsccos223 +?+= iiR (3.50)
( ) ( )
( ) ????
?????????
2
22
31
sinsincsc
sinsincsccoscoscoscos
iri
iriirR
?+
?+??=
(3.51)
( )
( )( )
( )( )??????????
???????
??????
sinsinsincoscossinsincsccos
sinsincossincoscos
sinsinsincsccos32
+?+
+???
??=
iri
ir
iriR
(3.52)
( )
( )( )
( )( )??????????
???????
??????
sincossinsincossinsincsccos
sincoscossinsincos
sincossincsccos33
+??+
+??
??=
iri
ir
iriR
(3.53)
3.6 Solving for the Elastic Stiffnesses, C
3.6.1 Minimization of Velocities
The Christoffel equation (Equation 2.20) can be expanded and rewritten as a
cubic equation in 2V? . The solution of this equation gives three roots: the quasi-
longitudinal velocity (VL), the slow transverse mode velocity (VT1), and the fast
transverse mode velocity (VT2). The equation can be solved for the particular velocities
using Cardan?s solution of a cubic equation (Mouchtachi et al. 2004), given by:
;33cos322 ??? ??
?
??
?
? ?=
LV (3.54)
;332cos3221 ?pi?? ??
?
??
?
? +?=
TV (3.55)
;332cos3222 ?pi?? ??
?
??
?
? ??=
TV (3.56)
41
where 3/ 2arccos ,
2 3
b
a
? ?
? ??
? ?? =
? ?? ?
? ?? ?? ?? ?
( )
( )
( )
2
3
11 22 33
2 2 2
12 13 23 11 22 11 33 22 33
2 2 2
11 22 33 12 13 23 11 23 22 13 33 12
,3
2 ,3 3
,
, and
2 .
a
b
? ?
?? ??
?
?
?
= ?
? ?= ? + ? ?
? ?
= ? ? + ? + ?
= ? ? + ? + ? ? ? ? ? ? ? ? ? ?
= ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ?
The elastic constants are recovered by minimizing the sum of squares of the
deviations between the experimentally measured velocities and Cardan?s solution, as
follows:
( ) ,21min
1
2?
=
?
m
i
c
i
e
i VV (3.57)
where m is the number of measurements,
e
iV are the experimental velocities, and
c
iV are Cardan?s solutions (Rokhlin and Wang 1992).
3.6.2 Minimization of the Christoffel Equation
An alternative method to that of minimization of the deviations of the velocities
using Cardan?s solution is minimization of the least-squares of the Christoffel equation
for each experimental velocity. For N directions of propagation, the Christoffel equation
can be expanded to a set of N cubic equations in m unknowns, where m is the number of
elastic constants to be determined. Generally, there are many more measurements than
42
elastic constants, and the set of equations is an overdetermined one (Castagnede et al.
1989). The recovery of the elastic constants is accomplished by minimizing the
functional:
( ) ( )( ){ } ,,
1
2?
=
=
N
p
ijppij CnfCF ? (3.58)
where ,2Vp ?? =
n is the wave propagation vector, and
( )( )ijpp Cnf ,? is the Christoffel equation.
For this method and for the previous one, the wave propagation vector is needed.
This may be computed using the rotation matrix developed in Section 3.5. The z?
direction corresponds to the refracted wave. Therefore, the third row of the rotation
matrix, R, may be used to calculate the wave vector.
3.6.3 Rotation Method
A third method for determining the elastic constants, proposed by Mayer and
Heidemann (1959), uses the tensor transformation equation for fourth rank tensors:
C' Cijkl im jn ko lp mnopa a a a= , (3.59)
where a denotes a matrix of direction cosines relating one set of axes to another.
In this method, rather than solving the Christoffel equation, the coordinate system
is rotated such that z? is in the direction of wave propagation. If the wave direction is in
the z? direction, C?33 correspond to the longitudinal, or quasi-longitudinal, wave velocity.
For a wave propagating in the z? direction, C?55 corresponds to the transverse wave
polarized in the x? direction, and C?44 corresponds to the transverse wave polarized in the
43
y? direction. For an orthotropic material, the equations used to determine these three
values are:
23
2
33
2
3213
2
33
2
3112
2
32
2
31
66
2
32
2
3155
2
33
2
3144
2
33
2
32
33
4
3322
4
3211
4
3133
222
444
'
CaaCaaCaa
CaaCaaCaa
CaCaCaC
+++
+++
++=
(3.60)
( ) ( ) ( )
233332232213333123211232312221
66
2
3221312255
2
3321312344
2
33223223
33
2
33
2
2322
2
32
2
2211
2
31
2
2144
222
'
CaaaaCaaaaCaaaa
CaaaaCaaaaCaaaa
CaaCaaCaaC
+++
++++++
++=
(3.61)
( ) ( ) ( )
233332131213333113111232311211
66
2
3211311255
2
3311311344
2
33123213
33
2
33
2
1322
2
32
2
1211
2
31
2
1155
222
'
CaaaaCaaaaCaaaa
CaaaaCaaaaCaaaa
CaaCaaCaaC
+++
++++++
++=
(3.62)
This method produces a system, usually overdetermined, of linear equations,
rather than a system of higher order equations as produced by the previous two methods.
For all three methods, measurements in three planes are required to solve for the full set
of elastic constants for an orthotropic material.
3.7 Solving the Systems of Equations
The numerical minimization procedure used to determine the set of elastic
constants when optimizing the Christoffel equation or the deviations of the experimental
and calculated velocities has traditionally been Newton?s algorithm. This algorithm
converges rapidly, but for higher-dimensional problems is unreliable unless initialized
near the exact solution. Mouchtachi et al. use the Levenberg-Marquardt method (LM
method), which combines features of the gradient method and Newton?s algorithm. For
the same set of data, the LM method and Newton?s algorithm produce different sets of
elastic constants. The results of the LM method are more consistent with findings using
44
classical tensile tests (Mouchtachi et al. 2004). The Levenberg-Marquardt method is
included in MATLAB and is used in this study to optimize the Christoffel equation and
the deviations of velocities.
3.7.1 The Levenberg-Marquardt Method
The Levenberg-Marquardt method is a numerical scheme used for nonlinear least
squares minimization. It combines the advantages of Newton?s algorithm and the
gradient descent method. The gradient descent method is updated by adding the negative
of the scaled gradient at each iteration as shown below:
fxx ii ??=+ ?1 . (3.63)
Newton?s algorithm can be expressed as:
( )( ) ( )iiii xfxfxx ???= ?+ 121 , (3.64)
where ( )( )ixf2? is the Hessian matrix (H) evaluated at ix .
The gradient method is accurate but slow to converge. Newton?s algorithm converges
rapidly but is sensitive to starting location. The Levenberg-Marquardt method combines
the advantages of these methods and is expressed as:
[ ]( ) ( )iii xfHdiagHxx ?+?= ?+ 11 ? . (3.65)
The Levenberg-Marquardt method uses the following update procedure:
1. Execute the algorithm above.
2. Evaluate the error.
3. If the error has increased, retract the step, increase ? , and return to step 1.
4. If the error has decreased, accept the step, and decrease ? .
45
The Levenberg-Marquardt method performs well for systems with hundreds of
parameters and achieves convergence much faster than the gradient descent method
(Ranganathan 2004).
MATLAB solves nonlinear least-square problems using the LM method with the
lsqnonlin command with the ?LargeScale? option set to ?off.? The command:
x = lsqnonlin(fun, x0), (3.66)
minimizes the sum of squares of the function fun with initialization at x0.
3.7.2 Least Squares Solution
When solving for the set of elastic constants using the rotation method, a solution
to an overdetermined system of linear equations is required. Given a general set of
overdetermined set of equations written in matrix form:
Ax =b, (3.67)
MATLAB solves for the values of the vector x in the least squares sense using the left
division command:
x = A\ b. (3.68)
(Palm 2001).
3.8 Damage
The micromechanical description of damage uses the parameters of cracks,
characteristic length, and characteristic spacing of constituents, as well as other factors, to
predict the macroscopic behavior of the composite. In the past, this approach to
describing damage progression has required that all the damage mechanisms be known a
priori and that the geometry and distribution of cracks be simple and regular. This is
46
rarely true for real materials. A direct measure of damage that can be experimentally
measured is the change in the elastic stiffness constants as damage progresses.
Under any load, an elastic stiffness can be represented as:
0C C CC= ? , (3.69)
where C0 is the elastic stiffness of the uncracked material, and Cc is the loss of stiffness
due to damage. The resulting change in the stiffness tensor, ?ij, is selected as a variable
representing the current state of damage of the material, as shown in Equation 3.70
below:
0C Cij ij ij? = ? i, j = 1 to 6. (3.70)
In many previous models used to describe damage, a scalar damage parameter has
been defined to be zero for an undamaged state and to be one at failure. This type of
description is inadequate for describing the tensorial nature of damage progression.
The normalized damage variable used in this study, D, is defined as the relative
change in the elastic stiffnesses. The terms on the diagonal of the elasticity matrix are
equal to zero at failure:
C 0ii = . (3.71)
So the limit of the damage variable is:
lim 0Cii ii? = . (3.72)
For the off-diagonal terms, the maximum value corresponds to a zero value of a minor of
the elasticity tensor (Baste and Audoin 1991), expressed as:
2C C C 0ij ii jj ijM = ? = . (3.73)
Hence:
47
lim 0C C Cij ij ii jj? = ? , .6,...,1, jiji ?= (3.74)
In this equation, the sign ? depends on the change in Cij, which can increase or decrease.
The components of the normalized damage tensor, D, are then given by:
lim 0C1 , 1 to 6,Cii iiii
ii ii
D i??= = ? = (3.75)
( ) ( )
0
lim 0 0 0 0
C C ,
C sign(C C ) C 1 C 1
, 1,...,6 .
ij ij ij
ij
ij ij ij ij ii ii jj jj
D
D D
i j i j
?
?
?= =
+ ? ? ?
= ?
(3.76)
The components of the damage tensor are measurable, have an identifiable physical
meaning, and form a finite set of data (Baste and Audoin 1991).
48
CHAPTER IV
DESCRIPTION OF MATLAB PROGRAMS
4.1 Introduction
The programs included were developed using MATLAB 7.0.4. The user begins
by inputting the number of samples to be evaluated. For each specimen, the user inputs
preliminary data: density (g/cm3), thickness (mm), reference velocity (m/s), reference
time (s), and a list of stress levels (MPa) that the sample is subjected to. The user then
inputs the rotation angles of the sample, and the resulting longitudinal and transverse
wave time values corresponding to these stress levels, as lists. A list of data is input in
brackets with a space between each value. For a set of readings:
time readings = [ ]1 2 3 ... nt t t t ,
spin = [ ]1 2 3 ... n? ? ? ? ,
tilt = [ ]1 2 3 ... n? ? ? ? ,
in-plane rotation = [ ]n???? ...321 .
The ith entry of each list forms a group describing one orientation and time reading. For
measurements in the 1-2 plane, only the angle of tilt is used. In the 1-3 plane only the
angle of spin is used. For a test in a non-principal plane, 1-2?, the angle of in-plane
rotation is constant and is equal to the angle between the 2? and 2 axes. The angle of tilt
is then equal to the incidence angle.
49
Figure 4.1 shows the initial setup for the ultrasonic immersion test. The
ultrasonic transmitter, T, and the receiver, R, are positioned on the 1 axis. The 2 and 3
axes define the plane of the sample, with the 3 axis in the long direction.
Figure 4.1 Initial orientations of sample and ultrasonic transmitter and receiver.
When using the minimization of velocities solution method, one must distinguish
the character of the transverse waves. Usually, the transverse waves in the 1-3 plane are
considered to be fast transverse waves, and those produced in the 1-2 plane are
considered slow transverse waves. Waves polarized in the fiber direction are faster than
the waves polarized perpendicular to the fiber direction (Rokhlin and Wang 1992).
However, for the minimization of the Christoffel equation, only the time readings and
corresponding angle measures are necessary. For the solution of rotation equations
method, one must distinguish between the types of polarization in a non-principal plane.
The fast transverse waves are considered as having in-plane polarization, and the slow
50
transverse waves have out-of-plane polarization. For the minimization of the velocities
and Christoffel equation methods, the user is also prompted to input a list of estimated
initial values for the nine elastic stiffnesses, C11, C22, C33, C44, C55, C66, C12, C13, and C23,
for the Levenberg-Marquardt method.
After the angle measures and time readings are input for each stress level, the
programs output tables of the elastic stiffnesses, C, elastic compliances, S, material
constants, E, G, and ?, and damage, D, for each stress level. Plots of each elastic
stiffness, elastic compliance, and damage variable versus the applied stress levels are also
created.
4.2 Minimization of the Christoffel Equation
The MATLAB 7.0.4 program developed for the minimization of the Christoffel
equation is presented in Appendix A.
4.3 Minimization of the Variation between Experimental and Calculated Velocities
The MATLAB 7.0.4 program developed for the minimization of the variation
between the experimental and calculated velocities is presented in Appendix B.
4.4 Least Squares Solution using Rotation of Axes
The MATLAB 7.0.4 program developed for the rotation of axes equation is
presented in Appendix C.
51
CHAPTER V
COMPARISON OF MATLAB PROGRAMS
5.1 Reconstruction of Velocity Data
The three MATLAB programs presented in Appendices A, B, and C were tested
using a set of velocity values back-calculated from ultrasonic immersion tests conducted
on a 1-D (unidirectional fiber) SiC-SiC ceramic matrix composite using data reported by
Aristegui and Baste (1997). The velocities were calculated in three data planes: the 1-2
plane, the 1-3 plane, and the 1-45o plane. The 45o axis is in the 2-3 plane and is forty-five
degrees from the 2 and 3 axes.
Figure 5.1 Unidirectional fiber composite showing the axes
defining the three planes of measurement: 1-2, 1-3, and 1- 2?.
The mass density of the composite material was 2.5 g/cm3, the sample thickness
was 3 mm, and the values of the nine independent elastic stiffnesses, determined from
52
experimental ultrasonic measurements similar to those described in this report, and
reported by Aristegui and Baste, are summarized in Table 5.1.
Table 5.1 Values of elastic stiffnesses for 1-D SiC-SiC (GPa) (Aristegui and Baste 1997).
C11 C22 C33 C44 C55 C66 C12 C13 C23
76 134 396 81 37.4 24.6 29 35 98
Figure 5.2 shows the measured data points and the slowness curves reconstructed
from the nine elastic stiffnesses in three planes of measurement (Aristegui and Baste
1997). The graphs in Figure 5.2 are polar plots of the slowness, pV1 , versus the
refraction angle, r? , in the three data planes. The longitudinal wave velocities are closest
to the origin, and the slow transverse wave velocities are farthest from the origin. The
data in Tables 5.2, 5.3, and 5.4 roughly correspond to the data points in Figure 5.2.
Figure 5.2 Reconstructed slowness curves and experimentally measured data points for a
1-D SiC-SiC composite in the (a) 1-2, (b) 1-3, and (c) 1-45o planes (Aristegui and Baste
1997).
Because Aristegui and Baste reported only the set of elastic constants determined
from the experimental measurements, and not the actual longitudinal and transverse wave
53
time data, a set of experimental velocity measurements for the orthotropic material was
simulated to check the validity of the MATLAB programs using the known values of the
elastic stiffnesses.
For the matrix
?
?
?
?
?
?
?
?
?
?
???
???
???
332313
232212
131211
,
where 23552266211111 CCC nnn ++=? ,
2
344
2
222
2
16622 CCC nnn ++=? ,
2
333
2
244
2
15533 CCC nnn ++=? ,
( ) 21661212 CC nn+=? ,
( ) 31551313 CC nn+=? ,
( ) 32442323 CC nn+=? , and
1n , 2n , and 3n are the components of the wave vector, n ,
the three eigenvalues are equal to 2LV? , 21TV? , and 22TV? , where ? is the density of the
composite, LV is the velocity of the quasi-longitudinal wave, 1TV is the velocity of the
slow transverse wave, and 2TV is the velocity of the fast transverse wave. The maximum
eigenvalue is equal to 2LV? , the minimum eigenvalue is equal to 21TV? , and the middle
eigenvalue is equal to 22TV? . These velocities were calculated in three planes for the 1-D
SiC-SiC ceramic matrix composite every 2o over a range of refraction angles, r? , from 0o
to 80o. The incidence angle, i? , for each wave was calculated using Snell?s Law:
54
??
?
?
???
?=
p
o
ri V
V?? sinarcsin
where the reference velocity, oV , was taken equal to 1450 m/s.
For measurements in the 1-2, 1-3, and 1-45o planes, the incidence angle is equal
to the degree of tilt, degree of spin, and degree of tilt, respectively. In the 1-45o plane,
the in-plane rotation is equal to 45o for all velocities. The set of velocities and angles
used in the three MATLAB programs represents a typical range of measurements of data
for an ultrasonic immersion test, as suggested by the data points shown in Figure 5.2.
The velocities, input angles, and refraction angles input into the MATLAB programs are
summarized in Tables 5.2, 5.3, and 5.4.
55
Table 5.2 Simulated velocity and angle measurements for quasi-longitudinal waves.
Quasi-longitudinal wave
velocity, VL (km/s) Degree of tilt Degree of spin
Degree of in-
plane rotation Refraction angle
5.51362 0 0 0 0
5.513821 0.525852 0 0 2
5.514463 1.050984 0 0 4
5.515655 1.574645 0 0 6
5.517582 2.096013 0 0 8
5.520505 2.614164 0 0 10
5.524759 3.128038 0 0 12
5.530757 3.63641 0 0 14
5.538983 4.137865 0 0 16
5.549982 4.630782 0 0 18
5.564348 5.11334 0 0 20
5.5827 5.583537 0 0 22
5.605644 6.039249 0 0 24
5.633735 6.478313 0 0 26
5.667429 6.898643 0 0 28
5.51362 0 0 0 0
5.517954 0 0.525458 0 2
5.53143 0 1.04776 0 4
5.555493 0 1.563351 0 6
5.592631 0 2.067874 0 8
5.646376 0 2.555848 0 10
5.721062 0 3.020606 0 12
5.82108 0 3.45482 0 14
5.949656 0 3.851803 0 16
6.107671 0 4.207147 0 18
6.293336 0 4.519724 0 20
6.502908 0 4.791426 0 22
5.51362 0 0 45 0
5.515885 0.525655 0 45 2
5.5229 1.049379 0 45 4
5.535335 1.569046 0 45 6
5.554343 2.082135 0 45 8
5.5816 2.58553 0 45 10
5.619303 3.075358 0 45 12
5.670077 3.546944 0 45 14
5.736694 3.99502 0 45 16
5.821587 4.414295 0 45 18
5.926249 4.800325 0 45 20
6.050806 5.150363 0 45 22
56
Table 5.3 Simulated velocity and angle measurements for fast transverse waves.
Fast transverse wave
velocity, VT2 (km/s) Degree of tilt Degree of spin
Degree of in-
plane rotation Refraction angle
4.35504 0 4.620025 0 14
4.436949 0 5.168129 0 16
4.503251 0 5.710394 0 18
4.5527 0 6.253684 0 20
4.586332 0 6.801757 0 22
4.606323 0 7.356026 0 24
4.615083 0 7.916548 0 26
4.614779 0 8.482751 0 28
4.607199 0 9.053836 0 30
4.593754 0 9.628949 0 32
4.575552 0 10.20724 0 34
4.553468 0 10.78788 0 36
4.528203 0 11.37001 0 38
4.50033 0 11.95277 0 40
4.470334 0 12.53521 0 42
4.414405 6.450419 0 45 20
4.461921 6.992343 0 45 22
4.497133 7.535655 0 45 24
4.520775 8.082787 0 45 26
4.534332 8.63439 0 45 28
4.5396 9.189815 0 45 30
4.538427 9.747472 0 45 32
4.532615 10.30499 0 45 34
4.523924 10.85918 0 45 36
4.514128 11.40594 0 45 38
4.505067 11.94001 0 45 40
57
Table 5.4 Simulated velocity and angle measurements for slow transverse waves.
Slow transverse wave
velocity, VT1 (km/s) Degree of tilt Degree of spin
Degree of in-
plane rotation Refraction angle
3.637803 10.78521 0 0 28
3.670122 11.39321 0 0 30
3.696181 11.99847 0 0 32
3.715556 12.60484 0 0 34
3.728078 13.21545 0 0 36
3.733817 13.83269 0 0 38
3.733043 14.4582 0 0 40
3.726173 15.0929 0 0 42
3.713735 15.73712 0 0 44
3.696317 16.39062 0 0 46
3.674538 17.05266 0 0 48
3.649026 17.72208 0 0 50
3.620398 18.39734 0 0 52
3.589253 19.07647 0 0 54
3.556167 19.75713 0 0 56
3.52169 20.43658 0 0 58
3.48635 21.11166 0 0 60
3.450647 21.77881 0 0 62
3.41506 22.43399 0 0 64
3.380046 23.07279 0 0 66
3.346039 23.69034 0 0 68
3.313451 24.28144 0 0 70
3.282671 24.84055 0 0 72
3.254062 25.36195 0 0 74
3.227962 25.83982 0 0 76
3.204678 26.26836 0 0 78
3.184486 26.64204 0 0 80
4.465526 0 12.54893 0 42
4.438631 0 13.11634 0 44
4.405591 0 13.69506 0 46
4.371547 0 14.27014 0 48
4.336808 0 14.84025 0 50
4.301664 0 15.4039 0 52
4.266389 0 15.95947 0 54
4.231248 0 16.50516 0 56
4.196496 0 17.03902 0 58
3.939178 14.81517 0 45 44
3.953095 15.29891 0 45 46
3.957615 15.79969 0 45 48
3.953035 16.31932 0 45 50
3.940427 16.85631 0 45 52
3.921286 17.40689 0 45 54
3.897177 17.96616 0 45 56
3.869516 18.52896 0 45 58
3.839503 19.09025 0 45 60
3.808127 19.64523 0 45 62
3.776201 20.18926 0 45 64
58
5.2 Results
From Table 5.5 below, it is clear that the results of the elastic stiffness values
delivered from the MATLAB programs are an excellent match with those reported by
Aristegui and Baste in two cases. The minimization of the Christoffel equation and the
minimization of the variation of the experimental and computed velocities provide
superior results compared to those obtained by the method of solving a system of rotation
equations. Certainly this outcome suggests that the latter solution method is not as well-
suited to problems of this type as are the former two solution methods.
The solution method that utilizes minimization of the Christoffel equation shows
advantages compared to the other two methods presented. Specifically, for this method,
no knowledge of the direction of polarization is required. A misinterpretation of the type
of polarization may lead to erroneous results when minimizing the velocities or solving
the system of rotation equations (i.e., when using the programs containing the other
solution methods).
Furthermore, the rotation solution method assumes an axis system for the input
wave that defines the direction of the wave, z?, and two orthogonal axes, x? and y?.
However, the refracted waves are nearly, but not exactly, longitudinal or transverse. The
orthogonal coordinate system (x?, y?, z?) approximates these quasi-longitudinal and
quasi-transverse waves as truly longitudinal or transverse, and it is therefore not exact.
The effects of these approximations are evident in the poor solutions given by the results
from the third method, reported in Table 5.5.
59
Computation of the damage variable was not performed as part of this check, as
only one set of velocities (corresponding to the response at only one stress level) was
calculated, since only one set of data was available for comparison.
Table 5.5 Comparison of results of MATLAB programs and Aristegui and Baste data.
C11 C22 C33 C44 C55 C66 C12 C13 C23
Data from
Aristegui and
Baste, 1997
76 134 396 81 37.4 24.6 29 35 98
Minimization of
Christoffel
equation
76 134 396 81 37.4 24.6 29 35 98
Minimization of
experimental
velocities and
Cardan?s
solution
76 134 396 81 37.4 24.6 29 35 98
Solution of
system of
rotation
equations
75.2 128 253 88.8 53.6 24.5 33.1 66.2 116
60
CHAPTER VI
CONCLUSIONS AND RECOMMENDATIONS
6.1 Conclusions
Three MATLAB programs have been written that can compute the full
complement of nine elastic constants for an orthotropic composite material, given a set of
experimental data collected from a series of ultrasonic immersion tests performed on a
composite sample that is simultaneously subjected to uniaxial tension. The data
necessary as input to the programs are a set of transit times corresponding to tests at
varying stress levels, and at varying orientations of the ultrasonic transducers with respect
to the composite specimen. The three programs differ only by their numerical solution
method; they each compute the same set of elastic constants. Each program also
computes a value indicating the amount of progressive damage imparted to the composite
sample due to the increasing uniaxial tensile stress.
By comparison with the set of experimentally determined elastic constants
reported by Aristegui and Baste (1997), it was shown that the programs can accurately
determine the nine elastic constants pertinent for an orthotropic material. Unfortunately,
there was no experimental data available for comparison of damage values.
When using the solution scheme that utilizes the minimization of velocities, the
solutions given for longitudinal, slow transverse, and fast transverse waves in Equations
3.54, 3.55, and 3.56 correspond to the maximum, minimum, and median eigenvalues.
61
For measurements made in a non-principal plane, all three waves are produced, and
identifying each wave type is straightforward. For measurements made in a principal
plane, only one transverse wave is produced, the character of which is unknown. For a
unidirectional composite, the transverse waves in the 1-2 plane are typically slow
transverse waves, and the transverse waves in the 1-3 plane are typically fast transverse
waves. However, as seen in Figure 5.2 (b), in a principal plane, the character of the
transverse wave may change from a fast transverse wave to a slow transverse wave as
indicated by the crossover of the slowness curves. Using simulated velocity data (as was
done to check the accuracy of the programs for calculating elastic constants), it was
simple to correctly identify the angle at which the wave changed from the fast transverse
wave to the slow transverse wave in the 1-3 plane. For an actual ultrasonic immersion
test, if and when the transverse wave changes character is unknown and thus correct
characterization is difficult.
6.2 Recommendations
It is recommended that for future study, the solution method of minimizing the
Christoffel equation should be expanded to include all possible types of symmetry, i.e. up
to twenty-one independent elastic constants. This would not only increase the range of
materials suitable for testing, but would also allow for a complete description of the
damage variable. A comparison of the MATLAB programs with experimental damage
data should also be performed.
Additionally, for the study conducted herein, it was assumed the composite
material maintains its class of symmetry during loading. This condition is satisfied if the
direction of crack growth coincides with the planes of symmetry of the composite, but
62
this may not always be the case (Audoin and Baste 1994). Modification of the program
to account for a change in symmetry class with loading is recommended.
Furthermore, further investigation using the minimization of the Christoffel
equation solution method should include a study to investigate how the number of
velocity readings, the range of refraction angles, and the change in initial estimates using
the Levenberg-Marquardt method affects the accuracy of the method.
63
REFERENCES
Achenbach, J. D. 1973. Wave propagation in elastic solids. North-Holland Series in
Applied Mathematics and Mechanics 16. ed. H. A. Lauwerier and W. T. Koiter.
Amsterdam: North-Holland Publishing Company.
Aristegui, C., and S. Baste. 1997. Optimal recovery of the elasticity tensor of general
anisotropic materials from ultrasonic velocity data. Journal of the Acoustical
Society of America 101, no. 2: 813-833.
Audoin, B., and S. Baste. 1994. Ultrasonic evaluation of stiffness tensor changes and
associated anisotropic damage in a ceramics matrix composite. Journal of Applied
Mechanics 61, (June): 309-16.
Baste, S., and B. Audoin. 1991. On internal variables in anisotropic damage. European
Journal of Mechanics, A/Solids 10, no. 6: 587-606.
Castagnede, B., J. T. Jenkins, W. Sachse, and S. Baste. 1990. Optimal determination of
the elastic constants of composite materials from ultrasonic wave-speed
measurements. Journal of Applied Physics 67, no. 6: 2753-2761.
Chawla, K. K. 2003. Ceramic matrix composites. 2nd ed. Boston: Kluwer Academic.
Chu, Y. C., and S. I. Rokhlin. 1992. Determination of macro- and micromechanical and
interfacial elastic properties of composites from ultrasonic data. Journal of the
Acoustical Society of America 92, no. 2, pt. 1: 920-931.
Hearmon, R. F. S. 1961. An introduction to applied anisotropic elasticity. London:
Oxford UP.
Jones, R. M. 1999. Mechanics of composite materials. 2nd ed. New York: Brunner-
Routledge.
Markham, M. F. 1970. Measurement of the elastic constants of fibre composites by
ultrasonics. Composites, (March): 145-49.
Mayer, W. G., and E. A. Hiedemann. 1959. Ultrasonic determination of elastic constants
and structural irregularities in transparent single crystals. (measurements in
sapphire). Acta Crystallographica 12: 1-6.
64
Mouchtachi, A., R. El Guerjouma, J. C. Baboux, D. Rouby, and D. Bouami. 2004.
Optimal determination of the elastic constants of woven 2D SiC/SiC composite
materials. Journal of Physics D: Applied Physics, 37: 3323-3329.
Palm, W. J. 2001. Introduction to MATLAB 6 for engineers. New York: McGraw- Hill.
Ranganathan, A. 2004. The Levenberg-Marquardt algorithm. August 1, 2006.
.
Rokhlin, S. I., and W. Wang. 1992. Double through-transmission bulk wave method for
ultrasonic phase velocity measurement and determination of elastic constants of
composite materials. Journal of the Acoustical Society of America 91, no. 6:
3303-12.
Suhling, J. C. 2004. Lecture notes, ME 7360, Mechanics of Composite Materials,
Department of Mechanical Engineering, Auburn University, Auburn, AL.
Wachtman, J. B., ed. 1989. Structural ceramics. Treatise on Materials Science and
Technology. Vol. 29. San Diego: Academic Press.
Warren, R., ed. 1992. Ceramic-matrix composites. New York: Chapman and Hall.
Weisstein, E. W. Euler angles. From MathWorld?A Wolfram Web Resource. August 1,
2006. .
65
APPENDIX A
Minimization of the Christoffel Equation
%This program develops and solves a system of nonlinear equations
%based on the Christoffel equation of an orthotropic composite
%Solves for nine elastic stiffnesses using Levenberg-Marquardt method
%Valid for incident waves on 2-3 plane, plane of the sample
%3 is the long axis, for unidirectional fiber composites 3 is the axis of fiber direction
%1 is the out-of-plane axis
%degrees of spin, tilt, rotation only valid up to 90 degrees
N = input('Enter the number of test specimens ');
disp('Input a list of data as [# # # #]');
for i = 1:N
disp('For sample');disp(i)
density = input('Input density of sample (g/cm^3) ');
density=1000*density;%Convert g/cm^3 to kg/m^3
h = input('Input sample thickness (cm) ');
h=h/100;%Convert cm to m
To = input('Input reference time (sec) ');
Vo = input('Input velocity of incident wave (velocity in immersion medium)(m/s) ');
stress = input('Input a list of stress levels (MPA) ');
for j = 1:length(stress)
disp('For stress level');disp(stress(j));
%Enter velocities and angles for each type of wave
disp('For measurements in the 1-2 plane, only tilt will be used')
66
disp('For measurements in the 1-3 plane, only spin will be used')
disp('For the longitudinal (quasi-longitudinal) waves:');
psi1 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha1 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi1 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t1 = input('Input a list of time readings for longitudinal waves (quasi-longitudinal waves)(sec)
');
disp('For the transverse (shear) waves:');
psi2 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha2 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi2 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t2 = input('Input a list of time readings for the transverse (shear) waves (sec) ');
%Calculate the time difference for each wave
T1 = To-t1;
T2 = To-t2;
%Calculate incidence angle for each wave
theta1 =acos(cos(psi1).*cos(alpha1));
theta2 =acos(cos(psi2).*cos(alpha2));
%Calculate the refraction angle for each wave
r1 =atan(sin(theta1)./(cos(theta1)-(T1*Vo/h)));
r2 =atan(sin(theta2)./(cos(theta2)-(T2*Vo/h)));
%Calculate the phase velocity for each wave
V1 =Vo*sin(r1)./sin(theta1);
V2 =Vo*sin(r2)./sin(theta2);
67
%For an incidence angle of zero (wave on 1 axis), correct values of V, r, and theta
for k=1:length(t1)
if psi1(k)==0 & alpha1(k)==0
V1(k)=Vo*h/(h-Vo*T1(k));
r1(k)=0;
theta1(k)=0;
end
end
%Square velocities, multiply by density, combine into one list
lambda1 = density*V1.*V1;
lambda2 = density*V2.*V2;
lambda = [lambda1 lambda2];
%Combine refraction, tilt, spin, inplane, and incidence angles for the three waves
r=[r1 r2];
psi = [psi1 psi2];
phi=[phi1 phi2];
alpha = [alpha1 alpha2];
theta=[theta1 theta2];
%The wave propagation vector for each wave
a1=[-cos(alpha).*cos(r-theta).*cos(psi)+cos(psi).^2.*csc(theta).*sin(alpha).^2.*sin(r-
theta)+csc(theta).*sin(r-theta).*sin(psi).^2];
a2=[-cos(psi).*csc(theta).*sin(r-theta).*sin(phi).*sin(psi)-cos(r-theta).*(-
cos(phi).*sin(alpha)+cos(alpha).*sin(phi).*sin(psi))+cos(psi).*csc(theta).*sin(alpha).*sin(r-
theta).*(cos(alpha).*cos(phi)+sin(alpha).*sin(phi).*sin(psi))];
a3=[-cos(phi).*cos(psi).*csc(theta).*sin(r-theta).*sin(psi)-cos(r-
theta).*(sin(alpha).*sin(phi)+cos(alpha).*cos(phi).*sin(psi))+cos(psi).*csc(theta).*sin(alpha).*sin(r-
theta).*(-cos(alpha).*sin(phi)+cos(phi).*sin(alpha).*sin(psi))];
%Correct values for incidence angle=0
68
for q=1:length(r)
if alpha(q)==0 & psi(q)==0
a1(q)=1;a2(q)=0;a3(q)=0;
end
end
%Create function F
%Christoffel equation for orthotropic symmetry
%the (c) variables are %C11 C22 C33 C44 C55 C66 C12 C13 C23
F= @ (c)
(a1.*a3.*c(8)+a1.*a3.*c(5)).*((a2.*a3.*c(9)+a2.*a3.*c(4)).*(a1.*a2.*c(7)+a1.*a2.*c(6))-
(a1.*a3.*c(8)+a1.*a3.*c(5)).*(a2.^2.*c(2)+a3.^2.*c(4)+a1.^2.*c(6)-lambda))-
(a2.*a3.*c(9)+a2.*a3.*c(4)).*(-
1*(a1.*a3.*c(8)+a1.*a3.*c(5)).*(a1.*a2.*c(7)+a1.*a2.*c(6))+(a2.*a3.*c(9)+a2.*a3.*c(4)).*(a1.^2.*c1
1+a3.^2.*c(5)+a2.^2.*c(6)-lambda))+(-
1*(a1.*a2.*c(7)+a1.*a2.*c(6)).^2+(a2.^2.*c(2)+a3.^2.*c(4)+a1.^2.*c(6)-
lambda).*(a1.^2.*c(1)+a3.^2.*c(5)+a2.^2.*c(6)-lambda)).*(a3.^2.*c(3)+a2.^2.*c(4)+a1.^2.*c(5)-
lambda);
%For first stress level, user inputs initial values of C as a list
%For subsequent stresses, program uses previous solution x
if j==1
xo=input('input initial values for [C11 C22 C33 C44 C55 C66 C12 C13 C23] in Pa');
else xo=x;
end
%Sets the lsqnonlin command to run LM method
options = optimset('LargeScale','off');
%Nonlinear least squares optimization of Christoffel equation using LM method
69
x=lsqnonlin(F,xo)
%store values in C table
%each row is for a different stress level
C(j,:)=x;
end
%Convert Pa to GPa
C=C/1000000000;
%Put values from completed C table into 6x6xstress array
for m=1:length(stress)
C6(1,1,m)=C(m,1);
C6(2,2,m)=C(m,2);
C6(3,3,m)=C(m,3);
C6(4,4,m)=C(m,4);
C6(5,5,m)=C(m,5);
C6(6,6,m)=C(m,6);
C6(1,2,m)=C(m,7);
C6(1,3,m)=C(m,8);
C6(2,3,m)=C(m,9);
C6(2,1,m)=C(m,7);
C6(3,1,m)=C(m,8);
C6(3,2,m)=C(m,9);
%Solve for elastic compliances
%S = inverse C
S6(:,:,m) = inv(C6(:,:,m));
%Store S values from 6x6xstress array in table
S(m,1)=S6(1,1,m);
S(m,2)=S6(2,2,m);
S(m,3)=S6(3,3,m);
70
S(m,4)=S6(4,4,m);
S(m,5)=S6(5,5,m);
S(m,6)=S6(6,6,m);
S(m,7)=S6(1,2,m);
S(m,8)=S6(1,3,m);
S(m,9)=S6(2,3,m);
end
%Compute damage variable
for p=1:length(stress)
D6(:,:,p)=1-C6(:,:,p)./C6(:,:,1);
end
%Compute off-diagonal damage terms
for k=1:length(stress)
for i=1:6
for j=1:6
if i~=j
D6(i,j,k)=(C6(i,j,1)-C6(i,j,k))./(C6(i,j,1)+sign(C6(i,j,1)-C6(i,j,k)).*(C6(i,i,1).*(1-
D(i,i,k)).*C6(j,j,1).*(1-D(j,j,k))).^(1/2));
end
end
end
end
%Store damage terms in table D
for m=1:length(stress)
D(m,1)=D6(1,1,m);
D(m,2)=D6(2,2,m);
D(m,3)=D6(3,3,m);
D(m,4)=D6(4,4,m);
71
D(m,5)=D6(5,5,m);
D(m,6)=D6(6,6,m);
D(m,7)=D6(1,2,m);
D(m,8)=D6(1,3,m);
D(m,9)=D6(2,3,m);
end
%Compute Young's moduli, Shear moduli, Poisson's ratios
E1=1./S(:,1);
E2=1./S(:,2);
E3=1./S(:,3);
G23=1./S(:,4);
G13=1./S(:,5);
G12=1./S(:,6);
pr12 = -E1.*S(:,7);
pr21 = -E2.*S(:,7);
pr13 = -E1.*S(:,8);
pr31 = -E3.*S(:,8);
pr23 = -E2.*S(:,9);
pr32 = -E3.*S(:,9);
%Display results
disp(' Stress C11 C22 C33 C44 C55 C66 C12 C13 C23')
disp([stress' C])
disp(' Stress S11 S22 S33 S44 S55 S66 S12 S13 S23')
disp([stress' S])
disp(' Stress E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23
pr32')
disp([stress' E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23 pr32])
72
disp('For sample');disp(i);
plot(stress,C),legend('C11','C22','C33','C44','C55','C66','C12','C13','C23');
figure
plot(stress,S),legend('S11','S22','S33','S44','S55','S66','S12','S13','S23');
%Plot grid of elastic stiffnesses versus stress graphs
figure
subplot(331);plot(stress,C(:,1)),xlabel('Stress (MPa)'),ylabel('C11 (GPa)');
subplot(332);plot(stress,C(:,2)),xlabel('Stress (MPa)'),ylabel('C22 (GPa)');
subplot(333);plot(stress,C(:,3)),xlabel('Stress (MPa)'),ylabel('C33 (GPa)');
subplot(334);plot(stress,C(:,4)),xlabel('Stress (MPa)'),ylabel('C44 (GPa)');
subplot(335);plot(stress,C(:,5)),xlabel('Stress (MPa)'),ylabel('C55 (GPa)');
subplot(336);plot(stress,C(:,6)),xlabel('Stress (MPa)'),ylabel('C66 (GPa)');
subplot(337);plot(stress,C(:,7)),xlabel('Stress (MPa)'),ylabel('C12 (GPa)');
subplot(338);plot(stress,C(:,8)),xlabel('Stress (MPa)'),ylabel('C13 (GPa)');
subplot(339);plot(stress,C(:,9)),xlabel('Stress (MPa)'),ylabel('C23 (GPa)');
figure
%Plot grid of elastic compliances versus stress graphs
subplot(331);plot(stress,S(:,1)),xlabel('Stress (MPa)'),ylabel('S11 (GPa-1)');
subplot(332);plot(stress,S(:,2)),xlabel('Stress (MPa)'),ylabel('S22 (GPa-1)');
subplot(333);plot(stress,S(:,3)),xlabel('Stress (MPa)'),ylabel('S33 (GPa-1)');
subplot(334);plot(stress,S(:,4)),xlabel('Stress (MPa)'),ylabel('S44 (GPa-1)');
subplot(335);plot(stress,S(:,5)),xlabel('Stress (MPa)'),ylabel('S55 (GPa-1)');
subplot(336);plot(stress,S(:,6)),xlabel('Stress (MPa)'),ylabel('S66 (GPa-1)');
subplot(337);plot(stress,S(:,7)),xlabel('Stress (MPa)'),ylabel('S12 (GPa-1)');
subplot(338);plot(stress,S(:,8)),xlabel('Stress (MPa)'),ylabel('S13 (GPa-1)');
subplot(339);plot(stress,S(:,9)),xlabel('Stress (MPa)'),ylabel('S23 (GPa-1)');
73
figure
%Plot grid of damage versus stress graphs
subplot(331);plot(stress,D(:,1)),xlabel('Stress (MPa)'),ylabel('D11');
subplot(332);plot(stress,D(:,2)),xlabel('Stress (MPa)'),ylabel('D22');
subplot(333);plot(stress,D(:,3)),xlabel('Stress (MPa)'),ylabel('D33');
subplot(334);plot(stress,D(:,4)),xlabel('Stress (MPa)'),ylabel('D44');
subplot(335);plot(stress,D(:,5)),xlabel('Stress (MPa)'),ylabel('D55');
subplot(336);plot(stress,D(:,6)),xlabel('Stress (MPa)'),ylabel('D66');
subplot(337);plot(stress,D(:,7)),xlabel('Stress (MPa)'),ylabel('D12');
subplot(338);plot(stress,D(:,8)),xlabel('Stress (MPa)'),ylabel('D13');
subplot(339);plot(stress,D(:,9)),xlabel('Stress (MPa)'),ylabel('D23');
end
74
Appendix B
Minimization of the Variation between Experimental and Calculated Velocities
%This program develops and solves a system of nonlinear equations
%based on the variation of the experimental velocities and Cardan's solution of the cubic
Christoffel equation
%for an orthotropic composite
%Solves for nine elastic stiffnesses using Levenberg-Marquardt method
%Valid for incident waves on 2-3 plane, plane of the sample
%3 is the long axis, for unidirectional fiber composites 3 is the axis of fiber direction
%1 is the out-of-plane axis
%degrees of spin, tilt, rotation only valid up to 90 degrees
N = input('Enter number of test specimens ');
disp('Input a list of data as [# # # #]');
for i = 1:N
disp('For sample');disp(i)
density = input('Input density of sample (g/cm^3) ');
density=1000*density;%Convert g/cm^3 to kg/m^3
h = input('Input sample thickness (mm) ');
h=h/1000;%Convert cm to m
To = input('Input reference time (sec) ');
Vo = input('Input velocity of incident wave (velocity in immersion medium)(m/s) ');
stress = input('Input a list of stress levels (MPA) ');
75
for j = 1:length(stress)
disp('For stress level');disp(stress(j));
%Enter velocities and angles for each type of wave
disp('For measurements in the 1-2 plane, only tilt will be used')
disp('For measurements in the 1-3 plane, only spin will be used')
disp('For the longitudinal (quasi-longitudinal) waves:');
psi1 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha1 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi1 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t1 = input('Input a list of time readings for longitudinal waves (quasi-longitudinal waves)(sec)
');
disp('For the fast transverse (shear) waves:');
psi2 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha2 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi2 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t2 = input('Input a list of time readings for fast transverse (shear) waves (sec) ');
disp('For the slow transverse (shear) waves:');
psi3 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha3 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi3 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t3 = input('Input a list of time readings for slow transverse (shear) waves (sec) ');
%Calculate the time difference for each wave
T1 = To-t1;
76
T2 = To-t2;
T3 = To-t3;
%Calculate the incidence angle for each wave
theta1 =acos(cos(psi1).*cos(alpha1))
theta2 =acos(cos(psi2).*cos(alpha2))
theta3 =acos(cos(psi3).*cos(alpha3))
%Calculate the refraction angle for each wave
r1 =atan(sin(theta1)./(cos(theta1)-(T1*Vo/h)));
r2 =atan(sin(theta2)./(cos(theta2)-(T2*Vo/h)));
r3 =atan(sin(theta3)./(cos(theta3)-(T3*Vo/h)));
%Calculate the phase velocity for each wave
V1 =Vo*sin(r1)./sin(theta1);
V2 =Vo*sin(r2)./sin(theta2);
V3 =Vo*sin(r3)./sin(theta3);
%For an incidence angle of zero (wave on 1 axis), correct values of V, r, and theta
for k=1:length(t1)
if psi1(k)==0 & alpha1(k)==0
V1(k)=Vo*h/(h-Vo*T1(k));
r1(k)=0;
theta1(k)=0;
end
end
%Square velocities, multiply by density, combine into one list
lambda1 = density*V1.*V1;
lambda2 = density*V2.*V2;
lambda3 = density*V3.*V3;
lambda = [lambda1 lambda2 lambda3];
77
%longitudinal wave vectors
La1=[-cos(alpha1).*cos(r1-
theta1).*cos(psi1)+cos(psi1).^2.*csc(theta1).*sin(alpha1).^2.*sin(r1-theta1)+csc(theta1).*sin(r1-
theta1).*sin(psi1).^2];
La2=[-cos(psi1).*csc(theta1).*sin(r1-theta1).*sin(phi1).*sin(psi1)-cos(r1-theta1).*(-
cos(phi1).*sin(alpha1)+cos(alpha1).*sin(phi1).*sin(psi1))+cos(psi1).*csc(theta1).*sin(alpha1).*sin(
r1-theta1).*(cos(alpha1).*cos(phi1)+sin(alpha1).*sin(phi1).*sin(psi1))];
La3=[-cos(phi1).*cos(psi1).*csc(theta1).*sin(r1-theta1).*sin(psi1)-cos(r1-
theta1).*(sin(alpha1).*sin(phi1)+cos(alpha1).*cos(phi1).*sin(psi1))+cos(psi1).*csc(theta1).*sin(alp
ha1).*sin(r1-theta1).*(-cos(alpha1).*sin(phi1)+cos(phi1).*sin(alpha1).*sin(psi1))];
%fast shear wave vectors
Qa1=[-cos(alpha2).*cos(r2-
theta2).*cos(psi2)+cos(psi2).^2.*csc(theta2).*sin(alpha2).^2.*sin(r2-theta2)+csc(theta2).*sin(r2-
theta2).*sin(psi2).^2];
Qa2=[-cos(psi2).*csc(theta2).*sin(r2-theta2).*sin(phi2).*sin(psi2)-cos(r2-theta2).*(-
cos(phi2).*sin(alpha2)+cos(alpha2).*sin(phi2).*sin(psi2))+cos(psi2).*csc(theta2).*sin(alpha2).*sin(
r2-theta2).*(cos(alpha2).*cos(phi2)+sin(alpha2).*sin(phi2).*sin(psi2))];
Qa3=[-cos(phi2).*cos(psi2).*csc(theta2).*sin(r2-theta2).*sin(psi2)-cos(r2-
theta2).*(sin(alpha2).*sin(phi2)+cos(alpha2).*cos(phi2).*sin(psi2))+cos(psi2).*csc(theta2).*sin(alp
ha2).*sin(r2-theta2).*(-cos(alpha2).*sin(phi2)+cos(phi2).*sin(alpha2).*sin(psi2))];
%slow shear wave vectors
Ta1=[-cos(alpha3).*cos(r3-
theta3).*cos(psi3)+cos(psi3).^2.*csc(theta3).*sin(alpha3).^2.*sin(r3-theta3)+csc(theta3).*sin(r3-
theta3).*sin(psi3).^2];
78
Ta2=[-cos(psi3).*csc(theta3).*sin(r3-theta3).*sin(phi3).*sin(psi3)-cos(r3-theta3).*(-
cos(phi3).*sin(alpha3)+cos(alpha3).*sin(phi3).*sin(psi3))+cos(psi3).*csc(theta3).*sin(alpha3).*sin(
r3-theta3).*(cos(alpha3).*cos(phi3)+sin(alpha3).*sin(phi3).*sin(psi3))];
Ta3=[-cos(phi3).*cos(psi3).*csc(theta3).*sin(r3-theta3).*sin(psi3)-cos(r3-
theta3).*(sin(alpha3).*sin(phi3)+cos(alpha3).*cos(phi3).*sin(psi3))+cos(psi3).*csc(theta3).*sin(alp
ha3).*sin(r3-theta3).*(-cos(alpha3).*sin(phi3)+cos(phi3).*sin(alpha3).*sin(psi3))];
%Correct values for incidence angle=0
for q=1:length(t1)
if alpha1(q)==0 & psi1(q)==0
La1(q)=1;La2(q)=0;La3(q)=0;
end
end
%The coefficients (c) in order
%C11 C22 C33 C44 C55 C66 C12 C13 C23
%Compute Cardan's solution in terms of elastic stiffnesses (c) for longitudinal waves
LG11 = @ (c) c(1)*La1.^2 + c(6)*La2.^2 + c(5)*La3.^2;
LG22 = @ (c) c(6)*La1.^2 + c(2)*La2.^2 + c(4)*La3.^2;
LG33 = @ (c) c(5)*La1.^2 + c(4)*La2.^2 + c(3)*La3.^2;
LG12 = @ (c) La1.*La2.*(c(7)+c(6));
LG13 = @ (c) La1.*La3.*(c(8)+c(5));
LG23 = @ (c) La2.*La3.*(c(9)+c(4));
Ldelta = @ (c) -LG11(c)-LG22(c)-LG33(c);
Lbeta = @ (c) -LG12(c).^2-LG13(c).^2-
LG23(c).^2+LG11(c).*LG22(c)+LG11(c).*LG33(c)+LG22(c).*LG33(c);
Lgamma = @ (c) -LG11(c).*LG22(c).*LG33(c)-
2*LG12(c).*LG13(c).*LG23(c)+LG11(c).*LG23(c).^2+LG22(c).*LG13(c).^2+LG33(c).*LG12(c).^2;
79
La=@ (c) (Ldelta(c).^2)/3-Lbeta(c);
Lb=@ (c) Lgamma(c)-Ldelta(c).*Lbeta(c)/3+2*(Ldelta(c)/3).^3;
Lpsic = @(c) acos(-Lb(c)./(2*(La(c)/3).^(3/2)));
Vc1=@ (c) 2*cos(Lpsic(c)/3).*(La(c)/3).^(1/2)-Ldelta(c)/3;
%Compute Cardan's solution in terms of elastic stiffnesses (c) for fast shear waves
QG11 = @ (c) c(1)*Qa1.^2 + c(6)*Qa2.^2 + c(5)*Qa3.^2;
QG22 = @ (c) c(6)*Qa1.^2 + c(2)*Qa2.^2 + c(4)*Qa3.^2;
QG33 = @ (c) c(5)*Qa1.^2 + c(4)*Qa2.^2 + c(3)*Qa3.^2;
QG12 = @ (c) Qa1.*Qa2.*(c(7)+c(6));
QG13 = @ (c) Qa1.*Qa3.*(c(8)+c(5));
QG23 = @ (c) Qa2.*Qa3.*(c(9)+c(4));
Qdelta = @ (c) -QG11(c)-QG22(c)-QG33(c);
Qbeta = @ (c) -QG12(c).^2-QG13(c).^2-
QG23(c).^2+QG11(c).*QG22(c)+QG11(c).*QG33(c)+QG22(c).*QG33(c);
Qgamma = @ (c) -QG11(c).*QG22(c).*QG33(c)-
2*QG12(c).*QG13(c).*QG23(c)+QG11(c).*QG23(c).^2+QG22(c).*QG13(c).^2+QG33(c).*QG12(c)
.^2;
Qa=@ (c) (Qdelta(c).^2)/3-Qbeta(c);
Qb=@ (c) Qgamma(c)-Qdelta(c).*Qbeta(c)/3+2*(Qdelta(c)/3).^3;
80
Qpsic = @(c) acos(-Qb(c)./(2*(Qa(c)/3).^(3/2)));
Vc2=@ (c) 2*cos((Qpsic(c)-2*pi)/3).*(Qa(c)/3).^(1/2)-Qdelta(c)/3;
%Compute Cardan's solution in terms of elastic stiffnesses (c) for slow shear waves
TG11 = @ (c) c(1)*Ta1.^2 + c(6)*Ta2.^2 + c(5)*Ta3.^2;
TG22 = @ (c) c(6)*Ta1.^2 + c(2)*Ta2.^2 + c(4)*Ta3.^2;
TG33 = @ (c) c(5)*Ta1.^2 + c(4)*Ta2.^2 + c(3)*Ta3.^2;
TG12 = @ (c) Ta1.*Ta2.*(c(7)+c(6));
TG13 = @ (c) Ta1.*Ta3.*(c(8)+c(5));
TG23 = @ (c) Ta2.*Ta3.*(c(9)+c(4));
Tdelta = @ (c) -TG11(c)-TG22(c)-TG33(c);
Tbeta = @ (c) -TG12(c).^2-TG13(c).^2-
TG23(c).^2+TG11(c).*TG22(c)+TG11(c).*TG33(c)+TG22(c).*TG33(c);
Tgamma = @ (c) -TG11(c).*TG22(c).*TG33(c)-
2*TG12(c).*TG13(c).*TG23(c)+TG11(c).*TG23(c).^2+TG22(c).*TG13(c).^2+TG33(c).*TG12(c).^2
;
Ta=@ (c) (Tdelta(c).^2)/3-Tbeta(c);
Tb=@ (c) Tgamma(c)-Tdelta(c).*Tbeta(c)/3+2*(Tdelta(c)/3).^3;
Tpsic = @(c) acos(-Tb(c)./(2*(Ta(c)/3).^(3/2)));
Vc3=@ (c) 2*cos((Tpsic(c)+2*pi)/3).*(Ta(c)/3).^(1/2)-Tdelta(c)/3;
81
%Combine the computed velocities for each wave
Vc=@ (c)[Vc1(c) Vc2(c) Vc3(c)];
%Create function of difference between experimental and computed velocities
F=@ (c) lambda-Vc(c);
%For first stress level, user inputs initial values of C as a list
%For subsequent stresses, program uses previous solution x
if j==1
xo=input('input initial values for [C11 C22 C33 C44 C55 C66 C12 C13 C23]');
else xo=x;
end
%Sets the lsqnonlin command to run LM method
options = optimset('LargeScale','off');
%Minimize sum of squares of the deviations between experimental and Cardan's solution
velocities using LM method
x=lsqnonlin(F,xo);
%store values in C table
%each row is for a different stress level
C(j,:)=x;
end
%Convert Pa to GPa
C=C/1000000000;
%Put values from completed C table into 6x6xstress array
for m=1:length(stress)
C6(1,1,m)=C(m,1);
C6(2,2,m)=C(m,2);
C6(3,3,m)=C(m,3);
C6(4,4,m)=C(m,4);
C6(5,5,m)=C(m,5);
82
C6(6,6,m)=C(m,6);
C6(1,2,m)=C(m,7);
C6(1,3,m)=C(m,8);
C6(2,3,m)=C(m,9);
C6(2,1,m)=C(m,7);
C6(3,1,m)=C(m,8);
C6(3,2,m)=C(m,9);
%Solve for elastic compliances
%S = inverse C
S6(:,:,m) = inv(C6(:,:,m));
%Store S values from 6x6xstress array in table
S(m,1)=S6(1,1,m);
S(m,2)=S6(2,2,m);
S(m,3)=S6(3,3,m);
S(m,4)=S6(4,4,m);
S(m,5)=S6(5,5,m);
S(m,6)=S6(6,6,m);
S(m,7)=S6(1,2,m);
S(m,8)=S6(1,3,m);
S(m,9)=S6(2,3,m);
end
%Compute damage variable
for p=1:length(stress)
D6(:,:,p)=1-C6(:,:,p)./C6(:,:,1);
end
%Compute off-diagonal damage terms
for k=1:length(stress)
for i=1:6
83
for j=1:6
if i~=j
D6(i,j,k)=(C6(i,j,1)-C6(i,j,k))./(C6(i,j,1)+sign(C6(i,j,1)-C6(i,j,k)).*(C6(i,i,1).*(1-
D(i,i,k)).*C6(j,j,1).*(1-D(j,j,k))).^(1/2));
end
end
end
end
%Store damage terms in table D
for m=1:length(stress)
D(m,1)=D6(1,1,m);
D(m,2)=D6(2,2,m);
D(m,3)=D6(3,3,m);
D(m,4)=D6(4,4,m);
D(m,5)=D6(5,5,m);
D(m,6)=D6(6,6,m);
D(m,7)=D6(1,2,m);
D(m,8)=D6(1,3,m);
D(m,9)=D6(2,3,m);
end
%Compute Young's moduli, Shear moduli, Poisson's ratios
E1=1./S(:,1);
E2=1./S(:,2);
E3=1./S(:,3);
G23=1./S(:,4);
G13=1./S(:,5);
84
G12=1./S(:,6);
pr12 = -E1.*S(:,7);
pr21 = -E2.*S(:,7);
pr13 = -E1.*S(:,8);
pr31 = -E3.*S(:,8);
pr23 = -E2.*S(:,9);
pr32 = -E3.*S(:,9);
%Display results
disp(' Stress C11 C22 C33 C44 C55 C66 C12 C13 C23')
disp([stress' C])
disp(' Stress S11 S22 S33 S44 S55 S66 S12 S13 S23')
disp([stress' S])
disp(' Stress E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23
pr32')
disp([stress' E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23 pr32])
disp('For sample');disp(i);
plot(stress,C),legend('C11','C22','C33','C44','C55','C66','C12','C13','C23')
figure
plot(stress,S),legend('S11','S22','S33','S44','S55','S66','S12','S13','S23')
%Plot grid of elastic stiffnesses versus stress graphs
figure
subplot(331);plot(stress,C(:,1)),xlabel('Stress (MPa)'),ylabel('C11 (GPa)');
subplot(332);plot(stress,C(:,2)),xlabel('Stress (MPa)'),ylabel('C22 (GPa)');
subplot(333);plot(stress,C(:,3)),xlabel('Stress (MPa)'),ylabel('C33 (GPa)');
subplot(334);plot(stress,C(:,4)),xlabel('Stress (MPa)'),ylabel('C44 (GPa)');
85
subplot(335);plot(stress,C(:,5)),xlabel('Stress (MPa)'),ylabel('C55 (GPa)');
subplot(336);plot(stress,C(:,6)),xlabel('Stress (MPa)'),ylabel('C66 (GPa)');
subplot(337);plot(stress,C(:,7)),xlabel('Stress (MPa)'),ylabel('C12 (GPa)');
subplot(338);plot(stress,C(:,8)),xlabel('Stress (MPa)'),ylabel('C13 (GPa)');
subplot(339);plot(stress,C(:,9)),xlabel('Stress (MPa)'),ylabel('C23 (GPa)');
figure
%Plot grid of elastic compliances versus stress graphs
subplot(331);plot(stress,S(:,1)),xlabel('Stress (MPa)'),ylabel('S11 (GPa-1)');
subplot(332);plot(stress,S(:,2)),xlabel('Stress (MPa)'),ylabel('S22 (GPa-1)');
subplot(333);plot(stress,S(:,3)),xlabel('Stress (MPa)'),ylabel('S33 (GPa-1)');
subplot(334);plot(stress,S(:,4)),xlabel('Stress (MPa)'),ylabel('S44 (GPa-1)');
subplot(335);plot(stress,S(:,5)),xlabel('Stress (MPa)'),ylabel('S55 (GPa-1)');
subplot(336);plot(stress,S(:,6)),xlabel('Stress (MPa)'),ylabel('S66 (GPa-1)');
subplot(337);plot(stress,S(:,7)),xlabel('Stress (MPa)'),ylabel('S12 (GPa-1)');
subplot(338);plot(stress,S(:,8)),xlabel('Stress (MPa)'),ylabel('S13 (GPa-1)');
subplot(339);plot(stress,S(:,9)),xlabel('Stress (MPa)'),ylabel('S23 (GPa-1)');
figure
%Plot grid of damage versus stress graphs
subplot(331);plot(stress,D(:,1)),xlabel('Stress (MPa)'),ylabel('D11');
subplot(332);plot(stress,D(:,2)),xlabel('Stress (MPa)'),ylabel('D22');
subplot(333);plot(stress,D(:,3)),xlabel('Stress (MPa)'),ylabel('D33');
subplot(334);plot(stress,D(:,4)),xlabel('Stress (MPa)'),ylabel('D44');
subplot(335);plot(stress,D(:,5)),xlabel('Stress (MPa)'),ylabel('D55');
subplot(336);plot(stress,D(:,6)),xlabel('Stress (MPa)'),ylabel('D66');
subplot(337);plot(stress,D(:,7)),xlabel('Stress (MPa)'),ylabel('D12');
subplot(338);plot(stress,D(:,8)),xlabel('Stress (MPa)'),ylabel('D13');
subplot(339);plot(stress,D(:,9)),xlabel('Stress (MPa)'),ylabel('D23');
end
86
Appendix C
Least Squares Solution using Rotation of Axes
%This program uses the rotation of axes equation for fourth rank tensors
%to develop and solve a system of linear equations with the nine unknown elastic constants
%of an orthotropic composite
%Valid for incident waves on 2-3 plane, plane of the sample
%3 is the long axis, for unidirectional fiber composites 3 is the axis of fiber direction
%1 is the out-of-plane axis
%degrees of spin, tilt, rotation only valid up to 90 degrees
N = input('Enter number of test specimens ');
disp('Input a list of data as [# # # #]');
for i = 1:N
disp('For sample');disp(i)
density = input('Input density of sample (g/cm^3) ');
density=1000*density;%Convert g/cm^3 to kg/m^3
h = input('Input sample thickness (mm) ');
h=h/1000;%Convert cm to m
To = input('Input reference time (sec) ');
Vo = input('Input velocity of incident wave (velocity in immersion medium)(m/s) ');
stress = input('Input a list of stress levels (MPA) ');
for j = 1:length(stress)
disp('For stress level');disp(stress(j));
87
disp('For principal planes only');
disp('For measurements in the 1-2 plane, only tilt will be used');
disp('For measurements in the 1-3 plane, only spin will be used');
disp('For the longitudinal (quasi-longitudinal) waves:');
psi1p = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha1p = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
t1p = input('Input a list of time readings for longitudinal waves (quasi-longitudinal
waves)(sec) ');
disp('For the slow transverse (shear) waves:');
psi3p = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha3p = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
t3p = input('Input a list of time readings for transverse (shear) waves (sec) ');
disp('For non-principal planes only');
disp('For the longitudinal (quasi-longitudinal) waves:');
psi1np = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha1np = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi1np = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t1np = input('Input a list of time readings for longitudinal waves (quasi-longitudinal
waves)(sec) ');
disp('For the fast transverse (shear) waves:');
88
psi2 = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha2 = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi2 = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t2 = input('Input a list of time readings for transverse (shear) waves with polarization out of
the incident plane(sec) ');
disp('For the slow transverse (shear) waves:');
psi3np = input('Input a list of degrees of spin (around 2 axis)(rad) ');
alpha3np = input('Input a list of degrees of tilt (around 3 axis)(rad) ');
phi3np = input('Input a list of in-plane rotation angles (around 1 axis) (rad) ');
t3np = input('Input a list of time readings for transverse (shear) waves with polarization in the
incident plane(sec) ');
t1=[t1p t1np];
t3=[t3p t3np];
psi1=[psi1p psi1np];
alpha1=[alpha1p alpha1np];
phi1=phi1np;
psi3=[psi3p psi3np];
alpha3=[alpha3p alpha3np];
phi3=phi3np;
%Calculate the time difference for each wave
T1 = To-t1;
T2 = To-t2;
89
T3 = To-t3;
%Calculate the incidence angle, theta, for each wave
theta1 =acos(cos(psi1).*cos(alpha1));
theta2 =acos(cos(psi2).*cos(alpha2));
theta3 =acos(cos(psi3).*cos(alpha3));
%Calculate the refraction angle, r, for each wave
r1 =atan(sin(theta1)./(cos(theta1)-(T1*Vo/h)));
r2 =atan(sin(theta2)./(cos(theta2)-(T2*Vo/h)));
r3 =atan(sin(theta3)./(cos(theta3)-(T3*Vo/h)));
%Calculate phase velocity of each wave, m/s
V1 =Vo*sin(r1)./sin(theta1);
V2 =Vo*sin(r2)./sin(theta2);
V3 =Vo*sin(r3)./sin(theta3);
%For an incidence angle of zero (wave on 1 axis), correct values of V, r, and theta
for k=1:length(t1)
if psi1(k)==0 & alpha1(k)==0
V1(k)=Vo*h/(h-Vo*T1(k));
r1(k)=0;
theta1(k)=0;
end
end
%Combine velocities, square and multiply by density
V=[V1 V2 V3];
lambda = (density*V.*V)';
90
%List of rotation matrix values for longitudinal waves
La11=-cos(r1-theta1).*cos(psi1).^2.*csc(theta1).*sin(alpha1).^2-
cos(alpha1).*cos(psi1).*sin(r1-theta1)-cos(r1-theta1).*csc(theta1).*sin(psi1).^2;
La12=cos(r1-theta1).*cos(psi1).*csc(theta1).*sin(phi1).*sin(psi1)-sin(r1-theta1).*(-
cos(phi1).*sin(alpha1)+cos(alpha1).*sin(phi1).*sin(psi1))-cos(r1-
theta1).*cos(psi1).*csc(theta1).*sin(alpha1).*(cos(alpha1).*cos(phi1)+sin(alpha1).*sin(phi1).*sin(p
si1));
La13=cos(r1-theta1).*cos(phi1).*cos(psi1).*csc(theta1).*sin(psi1)-sin(r1-
theta1).*(sin(alpha1).*sin(phi1)+cos(alpha1).*cos(phi1).*sin(psi1))-cos(r1-
theta1).*cos(psi1).*csc(theta1).*sin(alpha1).*(-
cos(alpha1).*sin(phi1)+cos(phi1).*sin(alpha1).*sin(psi1));
La21=0;
La22=cos(psi1).^2.*csc(theta1).*sin(alpha1).*sin(phi1)+csc(theta1).*sin(psi1).*(cos(alpha1).*cos(
phi1)+sin(alpha1).*sin(phi1).*sin(psi1));
La23=cos(phi1).*cos(psi1).^2.*csc(theta1).*sin(alpha1)+csc(theta1).*sin(psi1).*(-
cos(alpha1).*sin(phi1)+cos(phi1).*sin(alpha1).*sin(psi1));
La31=-cos(alpha1).*cos(r1-
theta1).*cos(psi1)+cos(psi1).^2.*csc(theta1).*sin(alpha1).^2.*sin(r1-theta1)+csc(theta1).*sin(r1-
theta1).*sin(psi1).^2;
La32=-cos(psi1).*csc(theta1).*sin(r1-theta1).*sin(phi1).*sin(psi1)-cos(r1-theta1).*(-
cos(phi1).*sin(alpha1)+cos(alpha1).*sin(phi1).*sin(psi1))+cos(psi1).*csc(theta1).*sin(alpha1).*sin(
r1-theta1).*(cos(alpha1).*cos(phi1)+sin(alpha1).*sin(phi1).*sin(psi1));
La33=-cos(phi1).*cos(psi1).*csc(theta1).*sin(r1-theta1).*sin(psi1)-cos(r1-
theta1).*(sin(alpha1).*sin(phi1)+cos(alpha1).*cos(phi1).*sin(psi1))+cos(psi1).*csc(theta1).*sin(alp
ha1).*sin(r1-theta1).*(-cos(alpha1).*sin(phi1)+cos(phi1).*sin(alpha1).*sin(psi1));
%List of rotation matrix values for out-of-plane polarization waves
91
Qa11=-cos(r2-theta2).*cos(psi2).^2.*csc(theta2).*sin(alpha2).^2-
cos(alpha2).*cos(psi2).*sin(r2-theta2)-cos(r2-theta2).*csc(theta2).*sin(psi2).^2;
Qa12=cos(r2-theta2).*cos(psi2).*csc(theta2).*sin(phi2).*sin(psi2)-sin(r2-theta2).*(-
cos(phi2).*sin(alpha2)+cos(alpha2).*sin(phi2).*sin(psi2))-cos(r2-
theta2).*cos(psi2).*csc(theta2).*sin(alpha2).*(cos(alpha2).*cos(phi2)+sin(alpha2).*sin(phi2).*sin(p
si2));
Qa13=cos(r2-theta2).*cos(phi2).*cos(psi2).*csc(theta2).*sin(psi2)-sin(r2-
theta2).*(sin(alpha2).*sin(phi2)+cos(alpha2).*cos(phi2).*sin(psi2))-cos(r2-
theta2).*cos(psi2).*csc(theta2).*sin(alpha2).*(-
cos(alpha2).*sin(phi2)+cos(phi2).*sin(alpha2).*sin(psi2));
Qa21=0;
Qa22=cos(psi2).^2.*csc(theta2).*sin(alpha2).*sin(phi2)+csc(theta2).*sin(psi2).*(cos(alpha2).*cos(
phi2)+sin(alpha2).*sin(phi2).*sin(psi2));
Qa23=cos(phi2).*cos(psi2).^2.*csc(theta2).*sin(alpha2)+csc(theta2).*sin(psi2).*(-
cos(alpha2).*sin(phi2)+cos(phi2).*sin(alpha2).*sin(psi2));
Qa31=-cos(alpha2).*cos(r2-
theta2).*cos(psi2)+cos(psi2).^2.*csc(theta2).*sin(alpha2).^2.*sin(r2-theta2)+csc(theta2).*sin(r2-
theta2).*sin(psi2).^2;
Qa32=-cos(psi2).*csc(theta2).*sin(r2-theta2).*sin(phi2).*sin(psi2)-cos(r2-theta2).*(-
cos(phi2).*sin(alpha2)+cos(alpha2).*sin(phi2).*sin(psi2))+cos(psi2).*csc(theta2).*sin(alpha2).*sin(
r2-theta2).*(cos(alpha2).*cos(phi2)+sin(alpha2).*sin(phi2).*sin(psi2));
Qa33=-cos(phi2).*cos(psi2).*csc(theta2).*sin(r2-theta2).*sin(psi2)-cos(r2-
theta2).*(sin(alpha2).*sin(phi2)+cos(alpha2).*cos(phi2).*sin(psi2))+cos(psi2).*csc(theta2).*sin(alp
ha2).*sin(r2-theta2).*(-cos(alpha2).*sin(phi2)+cos(phi2).*sin(alpha2).*sin(psi2));
%List of rotation matrix values for in-plane polarization waves
Ta11=-cos(r3-theta3).*cos(psi3).^2.*csc(theta3).*sin(alpha3).^2-
cos(alpha3).*cos(psi3).*sin(r3-theta3)-cos(r3-theta3).*csc(theta3).*sin(psi3).^2;
92
Ta12=cos(r3-theta3).*cos(psi3).*csc(theta3).*sin(phi3).*sin(psi3)-sin(r3-theta3).*(-
cos(phi3).*sin(alpha3)+cos(alpha3).*sin(phi3).*sin(psi3))-cos(r3-
theta3).*cos(psi3).*csc(theta3).*sin(alpha3).*(cos(alpha3).*cos(phi3)+sin(alpha3).*sin(phi3).*sin(p
si3));
Ta13=cos(r3-theta3).*cos(phi3).*cos(psi3).*csc(theta3).*sin(psi3)-sin(r3-
theta3).*(sin(alpha3).*sin(phi3)+cos(alpha3).*cos(phi3).*sin(psi3))-cos(r3-
theta3).*cos(psi3).*csc(theta3).*sin(alpha3).*(-
cos(alpha3).*sin(phi3)+cos(phi3).*sin(alpha3).*sin(psi3));
Ta21=0;
Ta22=cos(psi3).^2.*csc(theta3).*sin(alpha3).*sin(phi3)+csc(theta3).*sin(psi3).*(cos(alpha3).*cos(
phi3)+sin(alpha3).*sin(phi3).*sin(psi3));
Ta23=cos(phi3).*cos(psi3).^2.*csc(theta3).*sin(alpha3)+csc(theta3).*sin(psi3).*(-
cos(alpha3).*sin(phi3)+cos(phi3).*sin(alpha3).*sin(psi3));
Ta31=-cos(alpha3).*cos(r3-
theta3).*cos(psi3)+cos(psi3).^2.*csc(theta3).*sin(alpha3).^2.*sin(r3-theta3)+csc(theta3).*sin(r3-
theta3).*sin(psi3).^2;
Ta32=-cos(psi3).*csc(theta3).*sin(r3-theta3).*sin(phi3).*sin(psi3)-cos(r3-theta3).*(-
cos(phi3).*sin(alpha3)+cos(alpha3).*sin(phi3).*sin(psi3))+cos(psi3).*csc(theta3).*sin(alpha3).*sin(
r3-theta3).*(cos(alpha3).*cos(phi3)+sin(alpha3).*sin(phi3).*sin(psi3));
Ta33=-cos(phi3).*cos(psi3).*csc(theta3).*sin(r3-theta3).*sin(psi3)-cos(r3-
theta3).*(sin(alpha3).*sin(phi3)+cos(alpha3).*cos(phi3).*sin(psi3))+cos(psi3).*csc(theta3).*sin(alp
ha3).*sin(r3-theta3).*(-cos(alpha3).*sin(phi3)+cos(phi3).*sin(alpha3).*sin(psi3));
%The coefficients in order
%C11 C22 C33 C44 C55 C66 C12 C13 C23
93
%C'33
A1=[La31.^4;
La32.^4;
La33.^4;
4.*La32.^2.*La33.^2;
4.*La31.^2.*La33.^2;
4.*La31.^2.*La32.^2;
2.*La31.^2.*La32.^2;
2.*La31.^2.*La33.^2;
2.*La32.^2.*La33.^2];
%C'44
A2=[Qa21.^2.*Qa31.^2;
Qa22.^2.*Qa32.^2;
Qa23.^2.*Qa33.^2;
(Qa23.*Qa32+Qa22.*Qa33).^2;
(Qa23.*Qa31+Qa21.*Qa33).^2;
(Qa22.*Qa31+Qa21.*Qa32).^2;
2.*Qa21.*Qa22.*Qa31.*Qa32;
2.*Qa21.*Qa23.*Qa31.*Qa33;
2.*Qa22.*Qa23.*Qa32.*Qa33];
%C'55
A3=[Ta11.^2.*Ta31.^2;
Ta12.^2.*Ta32.^2;
Ta13.^2.*Ta33.^2;
(Ta13.*Ta32+Ta12.*Ta33).^2;
94
(Ta13.*Ta31+Ta11.*Ta33).^2;
(Ta12.*Ta31+Ta11.*Ta32).^2;
2.*Ta11.*Ta12.*Ta31.*Ta32;
2.*Ta11.*Ta13.*Ta31.*Ta33;
2.*Ta12.*Ta13.*Ta32.*Ta33];
%Correct values in A1 for incidence angle of zero
for k=1:length(lambda)
if psi(k)==0 & alpha(k)==0
A1(:,k)=0;
A1(1,k)=1;
end
end
%Combines the transpose of A1,A2,A3 into one matrix A
A=[A1';A2';A3'];
%Solve for %C11 C22 C33 C44 C55 C66 C12 C13 C23 in the least squares sense
X=A\lambda;
%store values in C table
%each row is for a different stress level
C(j,:)=X';
end
%Convert Pa to GPa
C=C/1000000000;
%Put values from completed C table into 6x6xstress array
for m=1:length(stress)
C6(1,1,m)=C(m,1);
C6(2,2,m)=C(m,2);
95
C6(3,3,m)=C(m,3);
C6(4,4,m)=C(m,4);
C6(5,5,m)=C(m,5);
C6(6,6,m)=C(m,6);
C6(1,2,m)=C(m,7);
C6(1,3,m)=C(m,8);
C6(2,3,m)=C(m,9);
C6(2,1,m)=C(m,7);
C6(3,1,m)=C(m,8);
C6(3,2,m)=C(m,9);
%Solve for elastic compliances
%S = inverse C
S6(:,:,m) = inv(C6(:,:,m));
%Store S values from 6x6xstress array in table
S(m,1)=S6(1,1,m);
S(m,2)=S6(2,2,m);
S(m,3)=S6(3,3,m);
S(m,4)=S6(4,4,m);
S(m,5)=S6(5,5,m);
S(m,6)=S6(6,6,m);
S(m,7)=S6(1,2,m);
S(m,8)=S6(1,3,m);
S(m,9)=S6(2,3,m);
end
%Compute damage variable
for p=1:length(stress)
D6(:,:,p)=1-C6(:,:,p)./C6(:,:,1);
end
96
%Compute off-diagonal damage terms
for k=1:length(stress)
for i=1:6
for j=1:6
if i~=j
D6(i,j,k)=(C6(i,j,1)-C6(i,j,k))./(C6(i,j,1)+sign(C6(i,j,1)-C6(i,j,k)).*(C6(i,i,1).*(1-
D(i,i,k)).*C6(j,j,1).*(1-D(j,j,k))).^(1/2));
end
end
end
end
%Store damage terms in table D
for m=1:length(stress)
D(m,1)=D6(1,1,m);
D(m,2)=D6(2,2,m);
D(m,3)=D6(3,3,m);
D(m,4)=D6(4,4,m);
D(m,5)=D6(5,5,m);
D(m,6)=D6(6,6,m);
D(m,7)=D6(1,2,m);
D(m,8)=D6(1,3,m);
D(m,9)=D6(2,3,m);
end
%Compute Young's moduli, Shear moduli, Poisson's ratios
E1=1./S(:,1);
E2=1./S(:,2);
E3=1./S(:,3);
97
G23=1./S(:,4);
G13=1./S(:,5);
G12=1./S(:,6);
pr12 = -E1.*S(:,7);
pr21 = -E2.*S(:,7);
pr13 = -E1.*S(:,8);
pr31 = -E3.*S(:,8);
pr23 = -E2.*S(:,9);
pr32 = -E3.*S(:,9);
%Display results
disp(' Stress C11 C22 C33 C44 C55 C66 C12 C13 C23')
disp([stress' C])
disp(' Stress S11 S22 S33 S44 S55 S66 S12 S13 S23')
disp([stress' S])
disp(' Stress E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23
pr32')
disp([stress' E1 E2 E3 G23 G13 G12 pr12 pr21 pr13 pr31 pr23 pr32])
disp('For sample');disp(i)
plot(stress,C),legend('C11','C22','C33','C44','C55','C66','C12','C13','C23')
figure
plot(stress,S),legend('S11','S22','S33','S44','S55','S66','S12','S13','S23')
%Plot grid of elastic stiffnesses versus stress graphs
figure
98
subplot(331);plot(stress,C(:,1)),xlabel('Stress (MPa)'),ylabel('C11 (GPa)');
subplot(332);plot(stress,C(:,2)),xlabel('Stress (MPa)'),ylabel('C22 (GPa)');
subplot(333);plot(stress,C(:,3)),xlabel('Stress (MPa)'),ylabel('C33 (GPa)');
subplot(334);plot(stress,C(:,4)),xlabel('Stress (MPa)'),ylabel('C44 (GPa)');
subplot(335);plot(stress,C(:,5)),xlabel('Stress (MPa)'),ylabel('C55 (GPa)');
subplot(336);plot(stress,C(:,6)),xlabel('Stress (MPa)'),ylabel('C66 (GPa)');
subplot(337);plot(stress,C(:,7)),xlabel('Stress (MPa)'),ylabel('C12 (GPa)');
subplot(338);plot(stress,C(:,8)),xlabel('Stress (MPa)'),ylabel('C13 (GPa)');
subplot(339);plot(stress,C(:,9)),xlabel('Stress (MPa)'),ylabel('C23 (GPa)');
figure
%Plot grid of elastic compliances versus stress graphs
subplot(331);plot(stress,S(:,1)),xlabel('Stress (MPa)'),ylabel('S11 (GPa-1)');
subplot(332);plot(stress,S(:,2)),xlabel('Stress (MPa)'),ylabel('S22 (GPa-1)');
subplot(333);plot(stress,S(:,3)),xlabel('Stress (MPa)'),ylabel('S33 (GPa-1)');
subplot(334);plot(stress,S(:,4)),xlabel('Stress (MPa)'),ylabel('S44 (GPa-1)');
subplot(335);plot(stress,S(:,5)),xlabel('Stress (MPa)'),ylabel('S55 (GPa-1)');
subplot(336);plot(stress,S(:,6)),xlabel('Stress (MPa)'),ylabel('S66 (GPa-1)');
subplot(337);plot(stress,S(:,7)),xlabel('Stress (MPa)'),ylabel('S12 (GPa-1)');
subplot(338);plot(stress,S(:,8)),xlabel('Stress (MPa)'),ylabel('S13 (GPa-1)');
subplot(339);plot(stress,S(:,9)),xlabel('Stress (MPa)'),ylabel('S23 (GPa-1)');
figure
%Plot grid of damage versus stress graphs
subplot(331);plot(stress,D(:,1)),xlabel('Stress (MPa)'),ylabel('D11');
subplot(332);plot(stress,D(:,2)),xlabel('Stress (MPa)'),ylabel('D22');
subplot(333);plot(stress,D(:,3)),xlabel('Stress (MPa)'),ylabel('D33');
subplot(334);plot(stress,D(:,4)),xlabel('Stress (MPa)'),ylabel('D44');
subplot(335);plot(stress,D(:,5)),xlabel('Stress (MPa)'),ylabel('D55');
subplot(336);plot(stress,D(:,6)),xlabel('Stress (MPa)'),ylabel('D66');
99
subplot(337);plot(stress,D(:,7)),xlabel('Stress (MPa)'),ylabel('D12');
subplot(338);plot(stress,D(:,8)),xlabel('Stress (MPa)'),ylabel('D13');
subplot(339);plot(stress,D(:,9)),xlabel('Stress (MPa)'),ylabel('D23');
end