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