EEJQSTEC WAVE PRGPAGATEON 3N AMSOTMPIC MEDM Thais far the Degree of Mn 3. mum STATE wNS‘iEliSi'l'Y JijERSI-SN hi". STAUME 19:3? 'w'r.1::>;'a‘ MW all—m.» LIBRARY Michigan State University ELASTIC WAVE PROPAGATION IN ANISOTROPIC MEDIA BY Juergen H. Staudte AN ABSTRACT OF A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Physics 1967 ABSTRACT ELASTIC WAVE PROPAGATION IN ANISOTROPIC MEDIA: by Juergen H. Staudte The propagation of plane elastic waves in an aniso- tropic material is investigated. Using a compact tensor nota- tion, the theory is developed such that the fundamental quantities necessary to describe wave propagation can be calculated. These include the velocities of the wavefronts of the acoustic beam, the particle displacement vectors, and the directions of the energy flow. A computer program, which numerically calculates these quantities, is included. A As the wave propagation properties of light are similar and in general, are better known, the optical theory is reviewed, using Maxwell's equations as a basis. Light propagation is studied by using such concepts as wave surfaces and normal velocity surfaces; the same surfaces are generated for elastic waves from the numerical results of the computer program. From this, the anomaly of cuspidal edges on the wave surfaces of elastic waves is discussed. The similarities between conical refraction for light and elastic waves are discussed. The pertinent quantities associated with elastic conical refraction are calculated for quartz. ELASTIC WAVE PROPAGATION IN ANISOTROPIC MEDIA BY Juergen H. Staudte A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Physics 1967 \J an \N. \\ \- \ '\ ACKNOWLEDGEMENT The author wishes to express his gratitude to Professor E. A. Hiedemann and Dr. B. D. Cook for their guid- ance in this work. Thanks are also due to Dr. F. Ingenito and Dr. W. R. Klein for many interesting and helpful discus- sions. The financial support of this research by the National Science Foundation is gratefully acknowledged. Acknowledgement is also due the Office of Naval Research for providing much of the equipment necessary for this research. ii CHAPTER II. TABLE OF CONTENTS INTRODUCTION . . . . . . . . . . . . . . THE THEORY OF THE PROPAGATION OF LIGHT IN ANISOTROPIC MEDIA . . . . . . . . . . III. WAVES IN IV. A. B. C. APPENDIX A. BIBLIOGRAPHY . . . . THE THEORY OF THE PROPAGATION OF ELASTIC ANISOTROPIC MEDIA . . . . . . . SPECIAL TOPICS . . . . . . . . . . . . . wave Surfaces with Cuspidal Edges Graphical Representation of Power Flow 0 O O O O O O O O O 0 O O Conical Refraction . . . . . . Mathematical Details for Elastic Wave Propagation . . . . . . Instructions for the Use of the Computer Program . . . . . . Computer Program . . . . . . iii 1h 25 25 28 29 35 LL3 A9 57 FIGURE 10 11 12 13 1h LIST OF FIGURES Diagram showing the relationship of vectors of an electromagnetic wave. . . . . . . . . . . . . . Diagram showing the propagation of characteristics PAGE of an electromagnetic wave in an anisotropic medium 9 The wave surface (inner) and normal velocity sur- face (outer) as defined by the ray vector and normal velocity vector . . . . . . . . . . . . Construction of the wave surface from the normal velocity and ray velocity . . . . . . . . . . . Construction of the envelopes defining plane waves from point sources using Huygens' principle . . Wave surfaces of aluminum.. . . . . . . . . . Photographs showing the wavefronts and power flow of quasilongitudinal (a) and quasishear (b) waves in quart z 0 O O O O O O O 0 O O O O O O O O 0 Schematic diagram of the displacement eigenvector and directions of power flow for an elastic wave prOpagating in the direction of the y-axis of '~ quartz . . . . . . . . . . . . . . . . . . . . Wave surfaces (a) and normal velocity surfaces (b) Of Zine. I O O O O O O C C C O O C O C O . wave surfaces (a) and normal velocity surfaces (b) of nickel . . . . . . . . . . . . . . . . . Graphical construction of the normal velocities and the direction of the propagation vector of plane waves which have power flows in the [110] direction of KBr. . . . . . . . . . . . . . . . Graphical representation of power flow of the quasilongitudinal mode in quartz . . . . . . . Conical refraction of light in biaxial crystals . Conical refraction in quartz . . . . . . . . . . iv 11 ll ll 18 20 21 26 26 27 3o 31 NOTATION CONVENTION A vector D linearly related to a vector E in the manner 'D, E... en. en E. Di 6;. ea," 613 EL D 3 £3. E31 E33 E3 where 1, 2 and 3 represent x, y and 2 respectively, can be written in the simpler form Di: 6;; E3 . Further simplification cert: be achieved by introducing the Einstein summation convention: when a small letter subscript occurs twice in the same term, summation with respect to that subscript is to be automatically understood. The above equation can, therefore, be written as D‘ ' 5;; ES . A quadratic term such as U: will mean U" U; and, therefore, summation will be understood. The subscriptsi, a , k and l will represent the coordinate components. The capital letter N used as a subscript will only be a label for the N'th root and, therefore, is excluded from the usual rules of tensor summation. All summation will be understood to be from 1 to 3. Any vector quantity such as E can have auxiliary quantities defined as follows: a) E the magnitude of 'E' ; A - b) E the unit vector in the direction of E ; c) El; the i'th component of the vector E . V (I. btx11lm‘q 1‘. (:1 (b (fi:£33!‘17{£?51 LIST OF SYMBOLS Dot above a letter denotes differentiation with respect to time. Unit vector along the i'th coordinate axis Spatial coordinate: x1 = x; x = y;:x = z. 2 3 A )hiiaé «a Xq,(:i Amplitude vector The i'th direction cosine of the N'th elastic displacement eigenvector Magnetic flux density vector Speed of light in a vacuum Elastic constants (rank four tensor) The i'th direction cosine of the N'th electric diaplacement eigenvector Electric displacement Electric field vector Magnetic field intensity vector Plane wave propagation vector A. 14. e; where 1". are the direction cosines denoting the direction of The N'th normal velocity vector Elastic Poynting vector The N'th ray vector Electric Poynting vector Time Particle displacement vector vi It"? 9» (v-‘c‘ 9 *3 8 Phase velocity of plane wave Defined as [/61 ’u' Potential energy per unit volume Kroenecker delta; 3;}- = Permittivity or strain tensor Components of diagonalized - tensor Permiability in a vacuum Stress tensor Angular frequency vii CHAPTER I INTRODUCTION The study of light and elastic wave propagation have a common history in their early development. With the publication of Huygens' treatise1 in 1690, the study of wave propagation began. Since light was thought to be prepagated through the aether, which was conceived to be an elastic solid, the elastic theory was devel- oped to a great extent. Not until 1817 was the transverse nature of light suggested by Young and, ten years later, incorporated into the theory of crystal optics by Fresnelz. From this point on, the longitudinal nature of wave pr0pagation was neglected for some time. Christoffel3 in 1877 wrote a memoir on elastic wave propagation in anisotropic media and developed a method for calculating phase velocities for three plane waves originating from an arbitrary disturbance in a plane. The technique for the determination of the characteristic displacement vectors associated with these plane waves was given by Lord Kelvin in his Baltimore Lectures)+ in 190%. In 1927, Love5 developed an expression to cal- culate the direction of the energy flow for each of the plane waves. wave behavior in crystals is usually analyzed by studying the propagation of a disturbance which originated at a point, or which is a plane wave. Both concepts are useful and correct, and they can be derived from each other. Huygens' principle considers only point sources; plane waves can be constructed from an infinite set of point sources. When the plane wave approach is used, the point source can be constructed from an infinite set of plane waves that travel in all directions. The surface of equal phase of a disturbance that is generated by a point source is called the wave surface. By using Huygens' construction, the direction and velocity of energy flow fo plane source can be determined. rA When a disturbance is a plane wave, the velocity normal to the plane of equal disturbance is called the normal velocity or the phase velocity of the plane wave. In practice, the point source approach is most useful in graphical techniques, and for con- structing models for better comprehension of the problem; but for analytical work, the plane wave techniques are more suitable because more powerful mathematical methods are available. The disturbances can be of quite different nature. Only two will be discussed in this paper, namely the electromag- netic wave in the form of light, and the focus of our discussion, the elastic wave. Light is studied in Chapter II for comparison and for the better understanding of the wave surfaces and wave vectors or elastic waves. Many different approaches for solving elastic wave propagation problems have been introduced over the years. Usually the different methods did not introduce anything new but they rather presented the methods of Kelvin and Love in a simpler form by taking special conditions into consideration. This paper will discuss in detail the earlier mentioned developments in tensor notation. The results of an experiment which renders visible elastic waves are used to illustrate the main features of elastic propagation in anistropic media. From some numerical results calculated using the computer program in Appendix C, several propagation phenomena of elastic waves are discussed. These include the existence of bicuspidal edges on wave surfaces. The appearance of bicuspidal edges imply that as many as five distinct signals from a point source may be observed. The possibility of this type of surface was first dis- covered by Musgrave6 about ten years ago. Surfaces calculated using the theory presented here are the same as those calculated by Musgrave - using a different technique. Conical refraction in optics is discussed and the basic principles are extended to elastic waves. A numerical example is presented for quartz. Also, a graphical method for presenting energy flow and characteristic particle displacement is suggested and demonstrated. CHAPTER II THE THEORY OF THE PROPAGATION OF LIGHT IN ANISOTROPIC MEDIA The fundamental properties of light proPagation in crystals can be derived from‘Maxwell's equations. For an infinite, homogenous, nondissipative material with no sources present, Maxwell's equations can be written as follows: 17xE§=-%J§ V. (1) V xfi: f5 , (2) and v.6‘VOfito . (3) The interaction of the electric field E with nonferrous materials is of the order of 105 larger than the interaction of the mag- netic field intensity F, . For this reason, only interactions of the electric field with the material are considered. The electric displacement 5 is linearly related to E by the dielec- tric or permittivity tensor in the following manner: 0 8 00 it Q. For convenience, we let the axes of the chosen coordin- ate system coincide with the principal axes of the ex; quadric, and then we can write D,= é1E1 , w) etc. The wave equation for an electromagnetic wave can now be constructed by substituting Eq. 2 into the curl of Eq. 1 which gives .0 VXVxE=-:-§i5=v(v-§)-V’E, (6) where [1, is the permeability(B=/J°H) and C. is the speed of light in a vacuum. A plane wave solution of Eq. 6 will be assumed; i.e., E=gexP[i(W«F-wt)]. (7) L( is the propagation vector whose components can be written as “La K1; where I; are direction cosines. It is also noted that the phase velocity 17' of the plane wave is fifi' . Equation 7 can, therefore, be rewritten as E‘éié; €XPLLK<£§X§’UL)]) (8) where 6%L is the coordinate unit vector. Substituting Eq. 8 into Eq. 6, we get ésx‘v‘D: -.- K‘1,;(I~E)- “11M: Ego, <9) Simplifying and rewriting Eq. 9 in component form and noting that Z’I‘ i , we obtain for one component ,L‘nn‘D,“ %+ AMI), (10) ct or D g 2' (IO?) ’ l/1£.-§/Aozyl~ I (11) and similar expressions for the other components. If we calcu- # late V'D‘O , we find £2011?) .1_1 == C) (12) U/EL — ”a D" ’ or on multiplying by/lo/Z'E, 1 .1; '_ C) (13) 'lfleilub"’1"' . Let us define Vi al/EL/lr Equation 11 and Eq. 12 can now be rewritten as D = 5'1. 1' (11+) ‘ I‘D V,‘-\7’- and all #1 ' =0 (15) V; " 177' Written out and simplified, Eq. 15 becomes 1 l ‘ " 16 A. “’1“ WW;- v‘) + UV:— v‘Xg— u‘) + £30,} 0%,; 03:0} ) 2. Equation 16 is a quadratic equation in 13 and therefore has 1 two roots which will be written as 17" where N = l, 2. A capital letter N used as a subscript will only be a label for the bl 'th root and, therefore, is excluded from the usual rules of tensor summation. These two velocities II" , known as the characteristic or eigenvelocities, are associated with two directions or eigenvectors into which an arbitrary electric displacement is resolved. These eigenvectors can be calculated by substituting tn: into Eq. 1h, which yields the components of the displacement vectors. The direction cosines d”; of the eigenvectors are now obtained from the normalized ratios of these components 16’- dNL 3 War-1 ) "5 where r", 2 , Du: z I . '“ l ._ (17) ]r =;.lZ!2.‘: 4&3E(\fl—”Uh) ' and ”2‘ D": 1. (VII' 35') ’ Kn : D", :——£3(V‘2- £— 3 DH: £.(v;-v£~) It can be shown that the two displacement eigenvectors are orthogonal by using Eq. 16 to prove that their scalar product is zero. With the information obtained up to this point, we can find the power flow using Poynting's theorem. The power flows along a vector 5 , known as the Poynting vector, where 3=%(EXH). (18) We are not interested in the magnitude of the power flow but just the direction. The electric vector, which is necessary to calculate 3; in Eq. 18, can be obtained from the electric dis- placement eigenvector by using Eq. 5. Two electric vectors is” are found and, therefore, two different gig will result. From Maxwell's equations, we know that D is perpendicular to K , and both if and 15 are perpendicular to F; . Equation 18 shows 3 to be perpendicular to E and H and, therefore,“5 , E ,R, and S are in the same plane (see Fig. 1). Figure 2 shows a three dimensional picture of a beam of light of arbitrary displacement A D A T E *5 A T— f’K 3 A 5 Figure 1. DEgram showing the relationship of vectors of an electromagnetic wave. For a plane wave with displacement and electric field E , the energy will flow along 5 and the wavefronts are normal to R . All vectors shown are in the plans. of the paper. The magnetic field ll (not shown) points perpendicularly out of the paper. I) D) MAE Figure 2. Diagram showing the propagation .7 characteristics of an electromagnetic wave in an anisotropic medium. The incoming light of arbitrary polarization upon entering the medium at B is re- solved into two beams of discrete polar- ization 5' and D; . The wave fronts A, B, C and C are parallel and denote the successive positions after increments of unit time. The normal velocity vectors n. , N1 show the magnitude of the vel- ocity in the direction of the ray vectors . , R1 indicate the ray velocities in the directions of the power flows. 10 entering an anisotropic solid at normal incidence to the surface. Illustrated in the figure are all previously discussed vectors for the case where the beam is split into two beams: one where 5 is parallel to E , and the other where it is not. If R- is varied and the subsequent end points of ii. and R; (as defined in Fig. 2) are plotted, the envelope of these points will produce two surfaces which are the wave surfaces. Using the same procedure for n. and N1. , two other surfaces are produced which are the normal velocity surfaces. From Fig. 2, we can see that h“ and ii; are perpendicular to their respective wavefronts. The end points of F7. and R; (which originate in the plane B) lie in the planes of their respective wavefronts. Each wavefront will be shown later to be tangent to a wave surface at some point. Figure 3 illustrates in two dimen- sions the construction of both surfaces by the method described above. The locus of endpoints of the vectors R. and n define the wave surface and normal velocity surface, respectively. For a particular R. and n , the line connecting the endpoints of the two vectors will lie in the plane of the wavefront. The construction and the relationship of the two types of surfaces to each other can be seen in Fig. 3. Having constructed the wave surfaces, the phase velocity and the ray direction can be obtained graphically from Fig. h. For a plane wave of direction 5F) (see Fig. ’4), the wavefront or line?“ must be constructed ‘1. in such a way as to be perpendicular to 0N and tangent to a point ll Figure 3. The wave surface (inner)and normal velocity surface (outer; as defined by the ray vector and normal velocity VGCCOI‘. 9&3“: Figure h. Construction of the wave surface from the normal velocity and ray velocity. \ M”. ’0 a o7 o " \ "_ \ Figure‘sz Construction of the envelopes defining plane waves from point sources using Huygens' principle. 12 on the wave surface which occurs at the only possible point of contact. This point of contact R will give us the direction of the ray 5? . The ray velocity is defined by the length OR ; the normal velocity '0'” is given by ON . The polarization or Ielectric displacement can be recovered by considering Fig. 1. The points R , 0 , and N define the plane in which the vectors 5 , E , R , and 3 lie. Therefore, we illustrate in Fig. 1L a simple construction that yields 5 and E . As two surfaces exist for electromagnetic waves, the previous considerations will hold for both, and therefore enable us to obtain the quantities for both resolved rays. If a point source of light were placed in a material, generally two sets of wavefronts would emanate from this point; each set having the shape of one of the wave surfaces. This is true, since such a point source may be considered to be composed of an infinite set of plane waves. Conversely, a plane source can be thought of as an infinite set of point sources which radiate coherently. The electric diaplacements of wave surfaces from each of the point sources will interfere and will only have nonvanish- ing amplitude in the envelope of all the wave surfaces. This construction is known as Huygens' construction and is illustrated in fig. 5. It shall be noted that quantitatively the use of the wave surfaces has many limitations. The graphical method which must be applied produces only approximate numerical results. 13 Moreover, these representations are limited to two dimensions. Nevertheless, they are very useful in explaining propagation phenomena qualitatively, in predicting anomalies, and in showing the overall picture. Their usefulness will become obvious when conical refraction is studied in a later chapter. In the next chapter the same approach as used in this chapter will be taken for obtaining some insight into the behavior of elastic waves; also, we will draw heavily on the results of this chapter. CHAPTER III THE THEORY OF THE PROPAGATION OF ELASTIC WAVES IN ANISOTROPIC MEDIA In order to study the behavior of elastic waves, we must determine essentially the same vectors, surfaces and veloc- ities as in the case of light. Again, we will construct the wave equation and the solution will be assumed to be a plane wave. We will let the stresses, strains and particle displacements be linearly related to each other and nonlinear effects, which are small, will be neglected. For a linear anisotropic material the following tensor relationship expresses Hooke's law: 513 " Gian. en. (19) where 513 are the components of stress tensor, 6&1. , the strain tensor and Ciikl are known as the elastic constants. The compon- ents of strain are related to the particle displacement l} by a-L 3.2L 2.9L 2 ELL glam-+3“). (O) Newton's second law for a continuous medium can be written as F1 = Paul (21) where F5 is the equilibrium density of the medium. This force F3 arises from stress gradients in the material, 1h 15 F:a : :5 Gris . 25 )QL (22) Combining Eq. 21 and Eq. 22, we get the following equation for a linear homogenous, continuous, anisotropic medium, 0- = 301:5 (23) o 6 25)(L . Substituting Eq. 19 and Eq. 20 into Eq. 23, we obtain the wave equation: " I 3 aUk 30.. .=-— c.. i a +..__ (an) For the above differential equation we will assume a plane wave solution giVen by U" AEXP[L(V-7-wt)], (25) The phase velocity of a plane wave is given by O:%and Ktkl, thus the displacement U can be written A O U = A5 65 exp[»k(x;l.,;-ut)], (26) Substituting Eq. 26 into Eq. 2h, we get (See Appendix A for details) (27) (Pill- 851 PowwAi .0) 16 where [-1 ‘-’ 1,11,. C'~ (28) Equation 27 has non-trival solutions only if the determinant of the coefficients of A5 is equal to zero: Fa _ 5&9on = 0° (29) Equation 29 is a cubic equation in P031 and, therefore, we can obtain three roots or velocities ”N . The three velocities are the phase velocities of three plane waves with the same propaga- tion direction I . In general, the wave of largest velocity is mostly longitudinal in nature and is commonly called the quasi- longitudinal wave. The other two waves are quasishear waves. In the study of light we have seen that only two trans- verse waves with different velocities can exist. Each eigen- velocity is associated with a displacement eigen vector which we calculate by substituting the eigenvelocities back into the equa- tion for is . The same technique can be used for elastic waves by substituting 31,: into Eq. 27. The displacement eigenvector components can be derived by solving for the ratios of A~5(See Appendix A). The normalized displacement eigenvector components are again direction cosines and are written as ONL' A plane wave of arbitrary diaplacement and direction I. will be resolved into three orthogonal displacement eigenvectors. The displacement 17 eigenvectors are completely uncoupled and therefore travel in- dependently with their corresponding eigenvelocities. For a particular plane wave, the amplitudes of the displacement eigen- vectors can be written as I - A . (30) AN " CN e}. ON)- where the CN‘S are the displacement amplitudes into which the initial displacement has been resolved. The three resolved plane waves can now be written by inserting Eq. 30 into Eq. 26 obtaining UN 2" CN an}, 6,; explLKUaia-vntfl. (31) As in light, we should suspect that the power flow of each of the three possible plane waves for a given I? is not along R- and neither parallel nor perpendicular to A" . The components of the power flow across a unit surface that is normal to the crystallographic axis are derived in Appendix A, and are given as A: ‘ ‘ (32) fi,‘ :=’ (713“Jj . Substituting Eq. 19, Eq. 20, and Eq. 31 into Eq. 32, we get (See Appendix A) .. J. a 3—3- . (33) Eu; ‘ 2C" ‘0'” C53 in. a", (1,a~k+ ham.) . 18 " P" P‘ I Normalizing under the condition that N! + ”a + ”3 a , we obtain the direction cosines of the energy flux or the components of the elastic Poynting vector. The angle of deviation 5" of the elastic Poynting vectors PM from the propagation vector K can be written as 605 5N=1LPNL (3“) where a,” in this equation are the direction cosines. The magnitude of the ray vector can be calculated from 17N (RAY) COS 5» The velocity of the disturbance which carries the energy is ‘0'" (35) called the ray velocity and, as in light, the wavesurface can be constructed by plotting all possible endpoints of the ray vector. In Fig. 6 an example is given of such surfaces for elastic waves in aluminum. , I’C’m Figure 6. Wave surfaces of aluminum. 19 For light the deviation angles of the power flow from the wave normals are at most a few degrees; but, for elastic waves, deviations can be as much as a radian. The deviation of power flow for a plane wave prOpagating in the direction o£:;;axis in quartz was numerically calculated to be about 25 degrees for all three waves. For this specific case, the results of an experiment9 which renders the sound beam.visible will be used to illustrate the properties of eigenvelocities, diSplacement eigenvectors, and the ray deviation. Figure 7 shows the direction of power flow and the wavefronts of a quasilongitudinal and a quasishear elastic wave. In order to fully interpret the photographs of Figure 7, a short review of the case of a plane wave propagating along the Y-axis in a sample of quartz is appropriate. Figure 8 shows two displacement eigenvectors and their associated power flow direc- tions as calculated from known elastic constants. The third dis- placement eigenvector is out of the page and its associated power flow is close to power flow of the shown quasishear wave. If a displacement is imposed along the Y-axis, as for example, by an X-cut quartz transducer which is mounted on a surface perpendicular to the Y-axis, two modes of displacement can be excited simultaneously. The transducer does not have a displacement component along the X-axis and therefore no motion is resolved in the direction the eigenvector which is directed along the X-axis. The quartz-air interface on the right side of 20 u 111 1 11”!"i. W;..E::; . ..11 111111111 1 11”,, u .‘1HH'I ..111111 (a) 1111111 1 .. . 111.111'11'111. 1111111111": WWI1MW 1"“‘1 “"1 11131111) Figure 7. Photographs showinggthe wavefronts and power flows of quasilongitudinal (affand quasishear (b) waves in quartz. Figure 8. 21 +Y Schematic diagram of the displacement eigenvector and directions of power flow for an elastic wave propagating in the direction 65 the y-axis in quartz. 22 Figure 7 reflects the two sound beams back along the same paths which they have traveled before reflection. By adjusting the driving frequency of the transducer, a stationary wave can be established. It is these two stationary waves that are rendered visible in Figure 7. In the lower left of Fig. 7a the outline of the x-cut quartz transducer, whose motion is along the y-axis, can be seen. This figure shows the quasilongitudinal mode with its wavefronts and power flow clearly visible. The motion of the x-cut trans- ducer is, of course, resolved into a longitudinal wave and a transverse or shear wave. The energy flow of the quasishear wave is downward and, therefore, is reflected about and eventually dissipated. Consequently, this shear wave has no chance to develop sufficient amplitude to be visible. The longitudinal beam is spreading and therefore fills up every possible space for this mode. The upper left and lower right triangular regions of the block contain no sound due to reflection from the upper and lower surfaces. If we move the transducer to the upper corner, a longitudinal resonance is now impossible because it would be reflected about as the quasishear wave is in Fig. 7a; however, the shear mode can be tuned to resonance. Figure 7b shows the quasishear mode in resonance. Due to the lower velocity of the quasishear mode, the separation of the wavefronts is less. In Chapter II, the electrical equivalent quantities were derived as in this chapter. From those quantities, the normal velocity surface, the wave surface, and the concept of 23 point sources and plane waves in Huygens' construction were devel- oped. Since the latter mentioned concepts are identical for elastic waves, no discussion is necessary but only a short comparison is appropriate. Of course, due to the greater complexity of elastic waves compared to electromagnetic waves, additional anomalies are found, and they will be discussed in the next chapter. The question whether Huygens' construction and its con- sequences are valid for elastic waves can, to some extent, be answered by comparing our numerical results with numerical results 6.7 given by Musgrave who used Huygens' construction as the founda- tion for his expressions. Musgrave transformed Eq. 27 into a form which can be interpreted as the analytical expression of the normal velocity surface. By studying Fig. 5, we can see that the wave surface can be graphically obtained from the normal velocity sur- faces by using Huygens' construction. Musgrave10 derived an analytic expression for the points on the wave surfaces by using the equations of the normal velocity surfaces and the equations required for constructing the tangent surfaces. Those points on the wave surface belong to a particular plane wave propagation vector and therefore give the ray velocities and ray directions directly. Musgrave calculated by his method the normal velocity 7 surfaces and the wave surfaces for some crystals of cubic and hexagonal6 symmetry. In.Appendix C, a computer program is described which numerically calculates all necessary quantities with the equations given in this chapter. The expressions given in this 2h chapter and the expressions given by Musgrave are mathematically of different form. Since it is not completely obvious by looking at the two forms that they will give the same numerical results, several surfaces previously calculated by Musgrave were recalcu- lated for comparison. As expccted, the results were identical and, therefore, both methods are compatible for calculating the direction of energy flow. The use of wave surfaces and Huygens' principle should give us valuable clues to the behavior of elastic waves. Fig. 6 shows such surfaces for aluminum. Its behavior is close to isotropic and therefore the surfaces are smooth and simple. More complicated surfaces are discussed in the next section. CHAPTER IV SPECIAL TOPICS A. Wave Surfaces With Cuspidal Edges A wave surface can be interpreted as the surface which describes the boundary between the disturbed and undisturbed region due to a wave emanating from a point source. For light, the wave surfaces are either of one sheet or, in general, at most of two sheets. As two wave surfaces are possible, in general, a flash from a point source would result in two flashes when observed through an anisotropic material far away from the source. From Fig. 6, we can see that elastic waves have three wave surfaces and, therefore, we would receive three pulses from a single pulse originating from a point. If we look caredully at the wave surfaces shown in Fig. 9 and 10, we see cuSpidal edges on the wave surfaces of one shear mode. Observation in the angular regions of the cuspidal edges will result in the reception of as many as five pulses rather than the normal three. This peculiar phenomenon usually occurs when the deviation of the power flow is large. A detailed study of cuspidal edges has been given by Musgrave8. We can use Huygens' construction to illustrate what happens when we consider plane waves rather than a point source. Using the techniques previously eXplained for Fig. 5, we have explored the bicuspidal edges on the wavesurface of KBr. In this example, the bicuspidal edge in the [001] plane is symmetric about the <110> direction. In Fig. 11, we have drawn tangent lines to the wave- 35 26 “u u 1 (a) f if if A 0;)” Figure 9. Wave surfaces (5) and normal velocity surfacesIIbI:of zinc. I—iCOT—u a Figure 10. Wave surfaces (a) and normal velocity surfaces (b) of nickel. Figure 11. 27 _ [ono] ’ \\ quasilongitudinal \ L \‘ [IOO] 1 2 3 ‘ I. xlO3 m/soc . [010] filflflbhoor . . .1100] ' 2 3 x '03 n/uc Graphical construction of the normal velocities and the direction of the prOpagation vectors of plane waves which have power flows in the [110] direction of KBr. 28 surfaces at every point where a surface crosses the 110 direction. The normal velocities and direction of the prOpagation vectors were then ascertained by constructing perpendicular intersectors which also passed through the origin. We find that five waves have their power flow in a particular 110 direction: namely, three wave (one longitudinal and two shear) with wave normals in the 110 direction, and two other shear waves whose respective propa- gation directions are approximately 100 and 800 from the 100 direction. Hence, the bicuspidal edges on wavefronts can be interpreted to mean that more than one plane wave have their power flowsin the same direction. Using the approximate values given above as a guide, better values for the propagation direction with power flow in the 110 direction were found by calculating the direction of power flow for several propagation directions about nearly 100 from the axis. These desired propagation directions were found to be 8.90 and 81.10 respectively. B. Graphical Representation of Power Flow It is now clear that wave surfaces are able to describe elastic wave behavior in a material qualitatively; however, quan- titatively they are not too useful except for approximate values. Another limitation is the difficulty of representing them in three dimensions. In order to pontray elastic wave behavior more accur- ately in three dimensions, several graphical techniques have been 29 devisedll. Figure 12 is an example of a graph which can be used to present the direction of the power flow or the displacement eigenvector for a particular mode in a specific material. Figure 12 shows the power flow of the quasilongitudinal mode in quartz. This graph gives a good overall view of energy flow of one mode in a specific material. The square grid represents the power flow direction of the longitudinal wave and the superimposed grid repre- sents the direction of the propagation vector ii . If ii is known and plotted on the superimposed grid, this point read from the square grid will give the direction of the power flow. The reverse procedure can be used if the direction of the power flow is known or if a particular beam direction is desired. C. Conical Refraction In optical biaxial crystals, light exhibits a peculiar phenomenon, known as conical refraction. An investigation of this phenomenon in light can be extended to elastic waves. By analyzing the wave surfaces of optical biaxial crystals, it is found that there are four directions where the intersection of the wave sur- faces resembles the top of a volcano. Figure 13a shows the cross- section of a hypothetical wave surface (this is a construct for explanation only). If a plane wave of single propagation direc- tion ”(0 were possible, the surface normal to 1‘0 ‘would touch the entire rim of the volcano-like wave surfaces. In the case 3O 3:33 23 mo 33m uoaon mo coaumucmmouqou Hmoanmwuo .NH 93w?“ cc 2 1.. I. 08 no. 31 W (b) Figure 13. Conical refraction of light in biaxial crystals. (a) shows the cross section of intersecting wave surfaces which1’ives rise to conical refraction. (b) shows the power flow of two slightly off axis wavevectors. where the propagation vectors deviate greatly from 12; , as F‘. and El; shown in Fig. 13, the energy would flow, according to Huygens' construction, to the points P. and a, for R, with normal velocities 17“ and V“ respectively and for R1 to the points P1 and Q3. with normal velocities 1),“ and ”a respectively. If we let X, and R1. approach 120 , the points P. and Q1 and P1 and Q. will be nearly at the same position and the normal velocities will approach each other. Any finite beam of light which propagates along “0 will consist of an infinite number of 32 — plane waves whose propagation vector is very nearly l‘o . Under ordinary conditions, a narrow pencil of light can be achieved at best. In Fig. 13b, the power flows of the two propagation vectors are shown. A pencil of light will produce an entire cone of light. The identical conditions are present in certain directions in anisotropic materials for transverse elastic waves. These con- ditions are met along trigonal axes such as in the z-direction of quartz or in the (111) directions in certain cubic crystals. For the above directions, the cone is circular, but in other special directions where conical refraction occurs the cross section of the cone is not necessarily circular. With the computer program in the appendix, calculations have been made for conical refraction occuring about the z-direction for quartz. The displacement eigen- vectors and the power flows have been computed and plotted in Fig. 1h for propagation vectors which are .01 degrees off the z-axis. Each R vector has two shear modes with their power flows nearly 16.80 from the z-axis. In Fig. 1h, the surface of a unit sphere is mapped into the plane of the paper with its pole (z-axis) pointing out of the paper. Points in this plane represent endpoints of vectors with their origin at the center of the sphere. The endpoints of the ik vectors all lie on the dashed circle and their corresponding power flow lie on the solid circle. Some A typical F< vectors are shown as solid and open dots on the dashed circle. The intersections of the line which passes through the 33 Figure 1h. Conical refraction in quartz. 311 A pole and the endpoint of a particular 1‘ vector with the solid circle gives the corresponding power flow vectors. 0n the end- points of the power flow vectors, the corresponding diSplacement eigenvectors all which lie in the xy-plane, are shown as arrows. The solid and dashed arrows correspond to solid and open dots, respectively. The direction of the power flow and the displace- ment vectors are very insensitive to small angular changes of the propagation vector R from the z-axis. Any finite source will have a distribution of a vectors about the z-axis and, there- fore, the power will flow to all points on the solid circle. In summation, I would like to point out that, in general, observed light propagation phenomena in anisotropic media are a result of the decomposition of a wave of arbitrary polarization into two eigenpolarizations which travel independently, with different normal velocities and ray directions. As we have seen, this.is.alae the case in sound, and, therefore, the method of analysis used in studying optical problems, such as ordinary refraction, conical refraction, reflection, circular and el- can liptical polarization, interference, etc., andAbe directly applied to elastic waves. APPENDIX A MATHEMATICAL DETAILS FOR ELASTIC WAVE PROPAGATION In this appendix we develop the details which are necessary in order to make proper use of the computer program in Appendix C; at the same time, the theory outlined in Chapter III is worked out in full. While going through the mathematical steps, reference to the computer program is made by giving the range of card numbers in double brackets as ((27-3h)). These cards will be the computer program's equivalent mathematical steps. Equation 2%, which is the wave equation of an anisotropic, linear homogeneous solid is written as follows: O. I A]. POUS=TC£5H 3 (43” +430) ( ) bx); ax]. a)". ' The solution assumed for the above wave equation is a plane wave U , (Eq. 25) where .. In . U = A; 8,; expthyQ " 17“]. (A2) Differentiating Eq. A2, we obtain bu- t A: “Us. “221.1,. (A3) 3X; 3X; ’ and u 1... _- . 1 1 U5 - U3 K '0'. (A11) 35 36 Substituting Eq. A3 and Eq. Ah into Eq. A1, we find 1 ’- ~-.'.. . ’- 1 (A5 P,” K Ui'zcisu“ “(1,141.3 X‘A;)fi’ ) where GL'Lthlgrvt). 1 Dividing by K 1.1, the above equation becomes P." 11“ = 2. C11“ 11(XLA.+1.A,), (A6) Since (3‘5“(1kAlJ is equivalent to C‘is‘lh‘I‘Ah) , Eq. A6 can be rewritten as = .L . (A7) P 1.7 A3 1 21(C53uLC15Lh)Ah Defining a new variable, namely . = . c.. (A8) PM. 1‘ A ‘1” 1 and noting that C “Cl , Eq. A6 can be rewritten as Liht iith (11(90’1553/‘111‘ 0. (A9) The above has non-trivial solutions only if the determinant of the coefficients is equal to zero: Pit. .. po 17" 61h” = O ) (A10) 37 '7.“ m” r l1, Fm +711 .9 1 5° C2 9 ,7] 0) II C) I3 [:3 [13‘ Pout The above is a cubic equation in P. '07" and therefore has three roots. The three roots are the eigenvelocities which are associated with three eigenvectors. The numerical calculations are made in the computer program from card #5 through 93. Substituting the eigenvelocities into Eq. A9, we obtain E. - Po“: [:1 r", AN! [:1 ’11-. 9,739: Ft) AM. ‘0. as [:5 F33. P00: ANS To calculate the three components of the three eigenvectors, we ) solve Eqs. All simultaneously for the ratio's of the AN]... 5. 38 One of the three possible forms of ratio's {3 b is written as f°11°ws= (wk-106)): Am -.-. I . ra A... t , P1 All}: «.- Es(r:.’gvu)” ELI-:3 . (A12) A m RU},- 8,0,3}- TL [:3 2 {i =Am = mm wan-mm AN‘ “1(113“ 9017;)‘C3n3 If two or all off-diagonal terms of r are zero, as is the case if R is along the crystal axis, the above ratios are indeter- minate. Therefore, we must test for that condition ((57)) and proceed to solve Eqs. All in a different way. From Eq. A10, we can see that if all off-diagonal terms are zero, it follows that Expat") Gk.) and the displacements are along the corresponding crystal axis ((132-138)). If only one off-diagonal term is non- zero, we obtain one po 17:; immediately with the corresponding A”; being along a coordinate axis. The remaining nonzero minor will be a quadratic equation in 9.1,; . Substituting the two roots of the quadratic equation into Eqs. A11 and nothing that the com- ponent of the pure displacement vector found initially will be set to zero, we obtain the ratios of the ANA; by solving Eqs. All simultaneously ((106-132) ) . In order to obtain the eigendisplacement direction 39 u cosines a”; , we must normalize by aNL: fi/(P;f((l39-l’+5)). The direction of power flow in an anisotropic material can be obtained by considering the energy change in a finite vol- ume per unit time. The energy change per unit time must be equal to the power flow out of the volume. If the components of power flow out of the volume are considered to be along the coordinate axes across a unit area, the elastic equivalent to the Poynting . vector can be obtained. If V is the potential energy per unit volume, the stress can be written as 3V (A13) 36“} - 61.3 ’ or am: “Eidezg’cuuiéu dc. (Am) Therefore, we find 5...; .. l .. \h/" ‘:l;¢.“ (igidalfizj ‘ ii (335 €55) ‘1 0 (A15) The total energy per unit volume is the sum kinetic energy 1- and potential energy W , or 6' (A16) 4- i . E =T+W=Ji Pot-’13. .L .. 1‘) The rate of change of energy in the unit volume can be obtained by differentiating (Eq. A16) with respect to time, E‘ Busus *z°71$€33+1°:)e~) (A17) The last two terms can be shown to be equal by use of Hooke's .0 law. From Eq. 22, P U; can be rewritten and substituted into 0 Eq. A16, resulting in - 'am- . 30- a . E: A U. + “—4. = .0 . (A18) 9X; ’5 "3X; aXL(G;!U5)' The total energy change in a finite volume can thus be written as ' 3 .. -— ' . ( ) Converting this integral into a surface integral by use of the divergence theorem, we get -‘ - '. (A20) P. - Ev ~ JACK»), db, s where: /5; is the component of surface area normal to the direction 1. Thus the elastic equivalent to the Poynting vector is P. = 01:3 0' , (A21) #1 Substituting Eq. 17 and Eq. 19 into Eq. A21, we obtain | 25‘) 2’th1 0 .g—c“ __L°_E.+_____. . (A22) PM 2 ‘JKL( 3x1. axh)u~5’ From Eq. 30, the incident plane wave has been resolved into diSplacement eigenvectors, where three plane waves were obtained of the form MC»: mi. 9.; ”P“ “1» (M1; " RN] . (A23) Taking the time and space derivatives of Eq. A23, we obtain . 1!. U~5=-,‘,un (Nam!- (Azu) and 3(1):. Liz :- u~C~1 (A25) Substituting Eq. AZH and Eq. A25 into Eq. A22, we find that . l. 2 l. 1. 1 PAN... 1 CN kN v" an; (ONE 1!. + a,“ ‘1); 1”. (A26) Time-averaging Eq. A26 over one cycle and noting that tif'fi?’ N we obtain ‘1- .L ’- ‘0 Puz‘zcn 1," (1.5", °~5(°~ulg* ain‘t). #2 Equation A27 gives us the elastic Poynting vector with the components of power flow across unit area in the direction of the crystal axes. The numerical calculations have been accomplished in the computer program by card number 1&6 through 156 and the normali- zation is accomplished by card numbers 157 through 161. z The normal velocity 17” is calculated from PD’D'N by cards 161 and 163. APPENDIX B INSTRUCTIONS FOR THE USE OF THE COMPUTER PROGRAM The program, written Fortran II, calculates: l)’ The directions of the displacement eigenvectors, in terms of a) direction cosines, b) angles in spherical coordinates. 2) The direction of the associated power flow, (ray direction) in terms of a) direction cosines b) angles in spherical coordinates. 3) The deviation of the ray direction from the plane wave normal. h) The normal and ray velocities. The input necessary for calculating these quantities is: l) The elastic constants of the material, 2) The density of the material, 3) The direction of propagation of the desired elastic plane wave. 3321.12 The input is achieved in the following sequence: Card Field Format Explanation l 1 5A8 Reads the program name 2-10 1-9 9F7.h Inputs the elastic constants. Table I lists the elastic constants in sequence as they are read. 1*3 Card 11 ll 11 11 11 Field 1 2 10 Format F10.3 F7.l 316 316 316 Ah Explanation Inputs density of material in kg/m3. Inputs the divisor by which the numbers of fields 3 through 9 are divided. The result of this division is the conversion of these numbers into degrees of angle. Reads the angles of the prOpagation vector in spherical coordinates. Field 3-6 determine the angle 9, with starting angle, angle limit, and angle steps listed, respectively. Angles g'in spherical coordinates are read in the same way ascthe angle 9. The angle 6 goes through all its values before another value is taken of the angle 9. The following three numbers (one in each field) define the reference axes for the angles 6 and fl. 1 2 3 ....¢ is rotated in the xy-plane measured from the x-axis toward the y-axis. 9 is referred from the z-axis. 1+5 Card Field Format Explanation 132 ........¢ is rotated in the xz-plane measured from the x-axis toward the z-axis. 9 is referred from the y-axis. 2 3 1 .....g is rotated in the yz- plane measured from the y-axis toward the z axis. 9 is referred from the x-axis. ll ll 16 The following code numbers give con- tinuation instructions: O....stop the program. 1....repeat card 11 with different data. 2....read in new set of cards as cards 1-11 with new data. OUTPUT For every set of calculations the direction of the propagation vector is specified in the heading. There are three major columns, each containing the calculated quantities for one mode with its associated eigen- velocity given in the last line. The abbreviations at the Cards 1 2 11 3 16 A 15 5 16 6 12 7 1A 8 15 9 1A 10 13 2 16 66 56 66 26 A6 56 A6 36 A6 TABLE I Fields (F7.h) 3 h 5 6 7 8 9 15 16 12 1A 15 1A 13 56 66 26 A6 56 A6 36 55 56 25 #5 55 1+5 35 56 66 26 A6 56 A6 36 25 26 22 2A 25 2A 23 A5 A6 2A AA A5 AA 3A 55 56 25 A5 55 A5 35 A5 A6 2A AA A5 AA 3A 35 36 23 3h 35 311 33 The numbers in the Table are the subscripts of the elastic constants, e.g., on card 5, field 6, 1+6 indicates the value of theelastic constant C116 should be inputted. (For details concerning four to two subscript conversion see Nye 11) beginning of each line indicate the following quantities: DVDC. Dvm. O 0 Displacement vector given in direction cosines showing the x,y and 2 components respectively. .Displacement vector given in spherical coordinates. ‘The first angle is measured from the x-axis and the second from thezaxis. 117 PVDG....Elastic Poynting vector given in direction cosines (same order as in DVDC). PVDG....Elastic Poynting vector given in sper- ical coordinates (same order as in DVDG). RVDV....Ray velocity and ray deviation in m/sec and degrees, respectively. PLWV....Plane wave phase velocity. As a sample problem, let us calculate the quantities which are necessary to plot Fig. 12. The material is quartz and, therefore, we name the program QUARTZ CRYSTAL as shown in sample data card 1. Nye11 lists the nonzero terms and equivalent values of the elastic constants for crystal class of symmetry 32 (such as alpha quartz). The following relationships between the constants were given: c11 = c22 °33 = c33 °AA = cAA 666 _ 1/2(611-612) c12 = c12 013 = c23 c1A = c56 ='°2A All other constants are zero. A8 12 The numerical values of the constants as given by Farnell are C11 = .869A c33 = 1.0680 °AA = .5762 c12 = .0696 c13 = .1560 °1A = .17A3 . These numerical values are shown in sample cards 2 - 10. Taking symmetry into account, we need only to vary the angle 9 through 90 degrees from the z-axis and the angle 6 through 120 degrees in the xy-plane in order to obtain the overall power flow behavior in a sample of quartz. The incremental angle was chosen to be 5 degrees. Sample data card 11 shows these instructions for 9 and 6. It is recommended that several features of the program be changed or deleted if other information is required. The refer- ences in the theory of Appendix A, the cross reference and the flow chart in Appendix C should enable such changes to be made without much difficulty. For greater precision in the calculated results, the program should be changed to double precision. By adding a plot routine to the program, wave surfaces, angles of deviation and many other parameters may be plotted. 101 READ 103, INP01,INF02,INP03,INPoA,INP05 APPENDIX C COMPUTER PROGRAM ODIMENSION 01(3,3.3),02(3.3.3).C3(3.3.3).D(3).G(3.3).V(3).S(3.3). 1DC(3.3).SN0RM(3).T(3.3).B(3),THETA1(a).FEE1(3).THETA2(3),FEE2(3). 2DIVIDE(3).RAYV(3),DEV(3) 10A DO 10 J=l,3 10 201 3200READ 609,Row,DIv,xNSTRT1,NLIM1.NSTEP1,NSTRT2,NLIM2,NSTEP2,IP,JP,KP, 331 300 111 READ 106, ((c1(J,R,L),L=1,3).K=1.3) READ 106, ((c2(J,R,L),L=1,3),K=1.3) READ 106, ((c3(J,K,L),L=1,3),K=1,3) DO 201 J=1,3 DO 201 K=1, 3 D0 201 L=l,3 C1(J,K,L)=C1(J,K,L)*l.Ell 02(J,K,L)=CZ(J,K,L)*1.E11 C3(J,K,L)=C3(J,K,L)*l.Ell lNOW PRINT 610 PRINT 128, INF01,INP02,INF03,INP0A,INP05 PRINT 388,IP,JP,KP NSTRT1=NSTRT1+1 NSTRT2=NSTRT2+1 NLIM1=NLD11 +1 NLIM2=NLIM2+1 D0 300 ID=NSTRTT,.NLIN1,NSTEP1 DO 300 JD=NSTRT2, NLIM2,NSTEP2 DEGl=ID-l DEG2=JD-l DEG1=DEGl/DIV DEG2=DEG2/DIV RADI=DEGI*.017A5329 RAD2=D262*.017A5329 D(IP)=COSF(RAD2)*SINF(RAD1) D(JP)=SINF(RAD1)*SINF(RAD2) D(KP)=COSF(RAD1) GO To 111 PRINT 7A9,DE62,IP,DE61,KP PRINT 750:((DC(I:J):J=1)3):I=113) PRINT 751.((FEE1(I).THEIA1(I)),I=1.3) PRINT 752.((S(K.L).K=1.3),L=1,3) PRINT 753,((PEE2(I),THETA2(1))I=1,3) PRINT 755,((RAYV(I),DEV(I)), I=l.3) PRINT 75A,v 1),v 2 ,v(3) IF (Now-1)361,320,101 DO 113 1:1,3 G(l,I)=O. G(2,I)=O. G(3,I)=O. ”9 l 10 15 20 25 3O 35 MO MS A8 50 D0 113 J=l,3 LL9 113 DC(l,J)=O. 50 D011A J=l,3 D011A I=1,3 D011A K=l.3 G(l.J)=Cl(I,J,K)*D(I)*D(K)+G(1,J) G(2,J)=C2(I,J,K)*D(I)*D(K)+G(2,J) 55 11“ G(3,J)=03(I,J,K)*D(I)*D(K)+G(3,J) IF (C(1,2)*G(1.3)+G(1,2)*G(2.3)+G(1.3)*G(2,3))107,302.107 107 B(1)=-G(1.1)-G(2,2)-G(3.3) BlOV3=B(1)/3.0 (03(2)=c(1,1)*c(2,2)+c(1,1)*c(3,3)+c(2,2)*c(3,3)-c(1,2)**2-G(1,3)**2 60 l-G(2,3)**2 61 03(3)=G(1, 3)**2*G(2.2)+G(1,2)**2*G(3. 3)+G(2.3)**2*G(1.1) 62 1 -G(1,1)*G(2,2)*G(3,3)-2-0*:G(1,2)*G(1.3)*G(2,3) AEF=B(2)-B(l)*BlOV3 BET=2.0*B10V3**3—B(2)*BlOV3+B(3) 65 BETOV2=BET/2.0 ALFOV3=ALF/3.0 CUAOV3=ALFOV3**3 SQBOV2=BETOV2**2 DEL=SQBOV2+CUAOV3 70 IP(DEL)7A0,720,730 720 GAM:SQRTF(-ALFOV3) IF(BET)722,722,721 721 v(1) = -2.0*GAM-BlOV3 V(2)=GAM-B10V3 75 V(3)=V(2) GO TO 750 722 V(1)=2.0*GAMrBlOV3 V(2)=-GAM-BlOV3 V(3)=V(2) 80 GO To 750 730 CUAOV3=~SQBOV2 GO TO 720 7A0 QUOT=SQBOV2/CUAOV3 ROOT=SQRTF(-QUOT) 85 11F