ISRN Applied MathematicsVolume 2011 (2011), Article ID 803562, 16 pagesdoi:10.5402/2011/803562

Research Article

Numerical Study of Characteristic Equations of Thermoelastic Waves Propagating in a Uniaxial Prestressed Isotropic Plate

Shih-Ming Hsu and Ching-Chung Yin

Department of Mechanical Engineering, National Chiao Tung University, Hsinchu 30010, Taiwan

Received 19 April 2011; Accepted 15 May 2011

Academic Editors: L.-Y. Li and M. Tabata

Copyright © 2011 Shih-Ming Hsu and Ching-Chung Yin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Thermoelastic waves propagating in an isotropic thin plate exerted by a uniaxial tensile stress are represented in this work. Characteristic equation of guided thermoelastic waves is formulated based on the theory of acoustoelasticity and classical thermoelasticity. Curve-tracing method for complex root finding is used to determine the attenuation, which is the imaginary part of the complex-value wavenumber. It is found that each plate mode of thermoelastic wave propagating in an isotropic plate with or without prestress has a minimum attenuation at a specific frequency except the A₀ mode. These modes are called the Lamé modes, which are the volume resonances in the thickness direction and propagate along the plate with the least energy dissipation. Frequency spectra of the phase velocity dispersion and attenuation of thermoelastic waves propagating along various orientations in the uniaxial prestressed thin plate have further been discussed.

1. Introduction

Determination of residual stresses in products is a major issue in most manufacturing industries. Both laser-induced ultrasound (LIU) [1, 2] and photoacoustics (PA) [3, 4] are the special techniques in “photothermal (PT) science” and have received intensive attention in nondestructive measurement of the residual stresses. The former determines the phase velocities from the far-field response generated by the pulsed laser, but the latter acquires the near-field response excited by the intensity-modulated CW laser. Both methods are correlative to the responses induced by “thermal acoustic waves.” Their dynamical behaviors are described based on the theory of thermoelasticity.

Research on the theory of dynamical thermoelasticity which includes the displacement and temperature coupled fields was studied and established well in a book [5]. However, for describing the thermal transmission that propagates by means of wave with finite speed, several modified theories for the hyperbolic-type energy equation have been proposed. Differing from the classical theory of thermoelasticity, these theories are called the “generalized thermoelasticity” and classified as the followings: L-S theories [6] with one relaxation time, G-L theory [7] with two relaxation times, and three types of G-N theory [8–10] involved with the energy dissipation.

On the other hand, the thermoelastic waves or laser-induced ultrasonic waves have become rather important in recent decades. Generalized thermoelastic wave propagation in the homogeneous transversely isotropic and anisotropic media has been analyzed [11, 12]. Afterward, based on four assumptions of the classical L-S, G-L, and G-N theories, Sharma et al. [13] investigated the thermoelastic wave propagation in a homogeneous isotropic plate, and showed the wavenumber spectra of the symmetric and antisymmetric modes for the insulated and isothermal boundary conditions. Verma and Hasebe [14, 15] and Verma [16] studied the wave propagation in the infinite homogeneous isotropic and orthotropic plates as well as the anisotropic layered media using the generalized thermoelasticity with relaxation times. Salnikov and Scott [17] developed the asymptotic models of long-wave low-frequency and short-wave approximations to analyze the dispersion relations of thermoelastic waves in an infinite homogeneous isotropic plate subject to either isothermal or thermally insulated traction-free boundary conditions.

This work presents an emerging method that unifies the advantages of LIU and PA without loss of generality for anisotropic inspection, which is of significance in the residual stress measurement. The classical theory of thermoelasticity [5] and the acoustoelasticity [18] are employed to model the thermoelastic waves propagating in a copper foil under distinct uniaxial stress state. The characteristic equations of phase velocity dispersion and attenuation are derived, and the spectra for each plate mode are determined numerically. Except for the fundamental antisymmetric (A₀) mode, each thermoelastic plate mode has a unique characteristic in its attenuation spectrum, in which a minimum value occurs at a specific frequency. guided thermoelastic waves induced at these specific frequencies can propagate farther away because of smaller attenuation.

2. Theoretical Formulation

2.1. Constitutive Relations

The thermoelastic effect in a stressed flat plate is formulated within the framework of the natural, initial, and final states, shown in Figure 1, originally proposed in [18] for the formulation of acoustoelasticity. Assume that there is no temperature change between the natural and initial states, that is, . The stress component and the increment of entropy satisfy the following constitutive relations: where . The physical field variables and indicate the strain changes and temperature rise caused by the external disturbance applied to the initial state. The material constants , , and denote the effective elastic constants, thermoelastic coupling coefficients, and thermal constant influenced by the initial strains , respectively. They are defined by     where and denotes the cubic dilatation. The terms , , , , , , and are the second- and third-order elastic constants, thermoelastic coupling coefficients, thermal constant, heat capacity, mass density, and temperature measured in the natural state. The relation between the initial strains and the initial (residual) stresses is given by .

Figure 1: Coordinate systems of the natural, initial or prestressed, and final states.

2.2. Governing Equations

The elastic wave propagation in a medium under residual stress must satisfy the equations of motion in the initial state of the form where is the mass density in the initial state and is defined by . The term is the body force applied to the initial state. Further, the balance of entropy and Fourier heat transfer equation in the initial state are of the form where , , and are the temperature, distributed body heat source, and surface heat flux in the initial state, respectively. is the effective thermal conductivity influenced by the initial strains and can be written as where is the thermal conductivity measured in the natural state. Substituting (2.1a), (2.1b) into (2.3), (2.4a), (2.4b) subsequently yields the partial differential equations of thermoelastic waves as follows:

2.3. Thin Isotropic Plate Subjected to a Uniaxial Prestress in the X₁-Direction

In this study, an isotropic thin plate in the natural state is considered. Referring the schematic diagram shown in Figure 2, the residual stresses are assumed to be homogeneous. Only normal stress is applied in the -direction, and shear stress components vanish, that is, (). According to the relation , the initial strains and the associated cubic dilatation of the isotropic plate are given as follows: The existing effective material constants , , and () in the initial state appeared in (2.2a), (2.2b), and (2.5) are derived as follows: According to (2.7) and (2.8), the special relations , , , , and are obtained. Rewriting (2.6) leads to the thermoelastic governing equations involving the effect of uniaxial prestress as follows: where , , and .

Figure 2: Schematic diagram of a copper foil prestressed in the -direction. The unit vector denotes the propagating direction in the -plane of guided thermoelastic waves.

2.4. Thermoelastic Waves Propagating along Any Direction from the -Axis

The schematic diagram of a thin plate exerted by normal stress in the -direction with the Cartesian coordinates is shown in Figure 2, in which the -axis is the thickness direction and the - and -axes both extend to infinity. This plate occupies the region , where is the thickness. In the absence of both and in (2.9), all field quantities are assumed to have a plane wave harmonic function based on the partial wave expansion (PWE) method. Then, the solutions of () and are represented as where and indicate the amplitudes of and , respectively. The term () denotes the wave vector of guided wave in the -plane and denotes the angular wavenumber in the thickness (-) direction. The components of wave vector are defined by and , where denotes the included angle of the arrowhead and the -axis as shown in Figure 2. The parameters and are the angular wavenumber and angular frequency defined as times the wavenumber and the frequency , which are the reciprocal of the wavelength and period, respectively.

Substitution of (2.10) into (2.9) is followed by the thermoelastic Christoffel equations for the coupled P, SV, SH, and thermal waves in the matrix form: where the unknown terms are defined as follows: According to the existence of a nontrivial solution in (2.11), the determinate vanishes and the eigenvalues () must satisfy the quartic equation of as follows: where the coefficients , and can be obtained and sorted by employing the software MAPLE for symbolic computation. The quartic equation (2.13) can be directly solved by quartic formula and eight complex roots () are subsequently obtained. For the definiteness of a single-value function, a constraint for root is selected as to avoid the exponential increasing. The symbol “+” corresponds to the downgoing waves propagating along the positive -direction. On the other hand, the symbol “−” denotes those upgoing waves traveling toward the negative -direction. The eigenvector components with respect to satisfy the proportional relation, where denotes the unknown amplitudes and the proportional factors are the determinants of corresponding submatrices (or minors) in (2.11) of the form Applying (2.1a) and (2.4b), the solutions of traction (, ) and heat influx () along the positive -direction can be represented in the form of where and indicate the amplitudes of and , respectively. Following the above procedure in a similar manner, the components with respect to also satisfy the proportional relation, where the proportional factors are given by Hence, applying the matrix technique [19] used in analysis for the layered structure, a combination of (2.10) and (2.16) yields where and . The submatrices and are the matrices with elements and as given in (2.15) and (2.18), respectively. The diagonal matrix has a rank 4 with plane wave forms . The uncertain vector with elements can be determined by the boundary conditions.

2.5. Characteristic Equation of Guided Thermoelastic Waves

Thermoelastic waves propagating in a thin plate are dispersive because of the geometric constraints on the upper and bottom boundaries. In addition, they are dissipative due to transformation between strain and thermal energies. Assume that the upper and bottom surfaces () are both traction-free and adiabatic conditions, that is, (). Applying the above boundary conditions to (2.19), the characteristic equations for the symmetric and antisymmetric modes of thermoelastic waves, indicated by the superscript “S” and “A”, are derived in the form where the elements of vector () are defined as follows:

In this study, due to energy dissipation of the propagating guided wave in the plate, the imaginary part of the complex-value wave number plays an important role. It is defined as , where and are the attenuation per unit distance (Np/mm) and per wavelength (Np/ν), respectively. The characteristic equation given in (2.19) can be represented as or . Moreover, because of the difficulty in obtaining an exact zero solution to numerically. Instead, minimizing is the common method for finding an approximate solution. Lowe [19] developed an effective method (also called the curve tracing algorithm) to find the complex roots of the characteristic equation for acoustic guided waves. Figures 3(a) and 3(b) show the dispersion and attenuation spectra for real wavenumber and imaginary wavenumber with respect to frequency . There exist many roots of with respect to one point on the -axis. Besides, the MATLAB v6.5 subroutine fminbnd [20] which is an algorithm by the golden section search and parabolic interpolation and the subroutine amoeba [21] based on the downhill simplex method are frequently used in the numerical computation.

(a)

(b)

Figure 3: Frequency spectra of (a) the dispersion (real wavenumber ) and (b) the attenuation (imaginary wavenumber ) for the symmetric modes (red lines) and antisymmetric modes (black lines) for the thermoelastic waves propagating in an isotropic copper foil under the natural state (without any prestress).

3. Numerical Results and Discussion

In this study, a thin copper foil is considered as an example. The material properties [22, 23] used in numerical computation are expressed in the units of “”, “”, “”, and “” and listed in the first column of Table 1. Furthermore, the effective material properties under two initial states with two uniaxial prestresses and applied in the -direction are assembled in the second and third columns of Table 1, respectively. Obviously, according to the relations , , and in the -plane, the material properties in two distinct initial states reveal the characteristic of nearly transversely isotropic. Material properties in the plane perpendicular to the direction of applied prestress are quasi-isotropic.

Table 1: Material properties of the copper foil under the natural state and two initial states with two uniaxial prestresses and applied in the -direction.

Natural state
Prestress
Prestress

Thickness(mm)

.1

Mass density (mg/mm³)

Elastic constants(mg/mmμs²)

Temperature(k°K)

Thermal constant(mg/mms²k°K²)

Thermoelastic coupling coeff. (mg/mms²k°K)

Thermal conductivity (mg/mms²k°K)

is the heat capacity (mg/s²k°K), and is the thermal expansion coefficient (1/k°K).

The speeds of the longitudinal (L0) and transverse (S0) waves, and Lamé mode (Lame0), and the Rayleigh wave (R0) in the natural state are, respectively, defined by      where Poison’s ratio for the copper foil. The formulas given in (3.1c) and (3.1d) can be referred to the books by Graff [24] and Achenbach [25], respectively. In the absence of the thermoelastic coupling effect, the speeds of the longitudinal (L) and transverse (S) waves along the X ₁-, X ₂-, and X ₃-directions in the initial (prestressed) state are defined as follows: where the superscript “” denotes the polarized direction of plane wave. Next, their corresponding differences of speeds between the natural and initial states are defined by The speeds of the Lamé modes along the - and -directions and their corresponding differences are given by and they are related to the speeds of transverse (S) waves along the -direction according to the thickness resonance behavior of the Lamé mode. Equation (3.9) can provid an exact solution the Lamé modes because of the nearly isotropic property in the -plane for whereas (3.8) only provided the trend of speed change for the Lamé modes owing to the orthorhombic property in the -plane. Substituting the material properties listed in Table 1 into (3.1a)–(3.1d) to (3.9), the resultant data for the natural states and two initial states, that is, two prestresses: and applyied in the -direction, are assembled in Table 2.

Table 2: Collected data of the wave speeds under the natural state and two initial states with two uniaxial prestresses and applied in the -direction.

Wave speed (mm/μs)
Longitudinal wave
Shear wave
Corresponding differences of wave speed due to (3.5)–(3.7)
Lamé mode

Natural state

—
—

Prestress

Prestress

According to the results via the complex root finding, the frequency spectra of the dispersion (real wavenumber ) and attenuation (imaginary wavenumber ) for the symmetric modes (red lines) and antisymmetric modes (black lines) of a thermoelastic wave propagating in the copper foil in the natural state (without any prestress) are shown in Figures 3(a) and 3(b). The frequency spectra of the phase velocity and semilogarithmic plot of attenuation are also shown in Figures 4(a) and 4(b), where the phase velocity is defined as . In Figure 4(a), the four additional blue dashed lines, labeled by "," " ," "," and "" indicate the phase velocities of values 4.590, 2.106, 2.978, and 1.972 mm/s due to (3.1a)–(3.1d). As frequency increases, both phase velocities of the A₀ and S₀ modes converge to a constant value corresponding to the Rayleigh wave speed 1.972 mm/s.

(a)

(b)

Figure 4: Frequency spectra of (a) the phase velocity and (b) the semilogarithmic plot of attenuation for the symmetric modes (red lines) and antisymmetric modes (black lines) for the thermoelastic waves propagating in an isotropic copper foil under the natural state (without any prestress).

In Figure 4(b), the attenuation spectrum of each mode, except the A₀ mode, has a close-to-zero minimum at a specific frequency, which is called the “Lamé mode” [24]. The Lamé modes travel at a specific wave speed and can be indicated using a dimensionless parameter for the symmetric modes Sn and for the antisymmetric modes Am . Reflecting on the “”-labeled blue dashed line shown in Figure 4(a), the specific frequencies, that is, for the Sn modes and for the Am modes, can be also observed. The Lamé modes represent the volume resonances in the thickness direction and propagate along the plate with the close-to-zero attenuation. This means that the thermoelastic waves can propagate farther away without energy dissipation. Moreover, the attenuations of the A₀ and S₀ modes merge together at the frequency range higher than 40 MHz. The convergent value is about  mm⁻¹ at 80 MHz.

Let the copper foil be exerted by a tensile stress in the -direction. Two cases of applied stresses and are considered. Figures 5(a) and 5(b) show the frequency spectra of phase velocity dispersion and semilogarithmic attenuation for the symmetric modes (red lines) and antisymmetric modes (black lines) of thermoelastic waves propagating along the -direction. Figures 6(a) and 6(b) show those of thermoelastic waves propagating along the -direction. Comparison of Figures 5(a) and 6(a) indicates that the shift of each dispersion curve depends on the bulk wave speed differences between the natural and stressed states. According to Table 1, the changes of bulk wave speeds and in the -direction due to the tensile prestress are −0.044 and 0.024 mm/s, respectively. As shown in Figure 5(a), the resulting phase velocities of plate waves increase over a broad frequency range but decrease in the vicinity of longitudinal wave speed 4.590 mm/s. Similarly, the phase velocities of the plate waves along the -direction decrease in accordance with the changes of bulk wave speeds −0.005 and −0.050 mm/s. Figure 6(a) shows that the dispersion curves of higher modes change distinguishably in the region of higher frequency and wavenumber.

(a)

(b)

Figure 5: Frequency spectra of (a) the phase velocity and (b) the semilogarithmic plot of attenuation for the symmetric modes (red lines) and antisymmetric modes (black lines) for the thermoelastic waves propagating along the -direction in a copper foil under the natural state (without any prestress) and two initial states with two uniaxial prestresses and applied in the -direction.

(a)

(b)

Figure 6: Frequency spectra of (a) the phase velocity and (b) the semilogarithmic plot of attenuation for the symmetric modes (red lines) and antisymmetric modes (black lines) of a thermoelastic wave propagating along the -direction in the copper foil under the natural state (without any prestress) and two initial states with two uniaxial prestresses and applied in the -direction.

On the other hand, as shown in Figures 5(b) and 6(b), the attenuation spectrum of each mode has a minimum at a specific frequency. These specific modes represent the volume resonances in the thickness direction and propagate along the plate with the least energy dissipation. Figure 5(b) shows that the minimum attenuation of each mode, except the A₀ mode, increases as the tensile prestress in the same direction as wave propagation increases. This phenomenon is caused by the uniaxial prestress and leads to the orthorhombic symmetry in the -plane, for example, the values , , , , , , , and shown in Table 1. Next, Figure 6(b) shows that the specific frequencies of the Lamé modes are reduced if a tensile prestress is applied in the orientation perpendicular to the direction of wave propagation. The reductions of these specific frequencies can result from the changes of isotropic property in the -plane, that is, the specific relations , , and shown in Table 1. Therefore, the feature of close-to-zero attenuation of the Lamé modes is directly correlated with the isotropic material property in the sagittal plane of propagating thermoelastic guided wave.

In the previous illustrations, the horizontally polarized motion (SH wave) has been decoupled from the thermoelastic waves propagating along the - and -directions. Figures 7(a) and 7(b) show the frequency spectra of phase velocity dispersion and semilogarithmic attenuation of thermoelastic waves propagating along the direction inclined at 45° to the -axis in a copper foil which is tensilely prestressed by in the -direction. Owing to that the off-axis traveling wave motions are neither polarized in the sagittal plane nor horizontally polarized, and the P, SV, SH, and thermal waves are coupled. As shown in Figure 7(b), the “even” numbered antisymmetric modes (solid black lines), except the A₀ mode, and the “odd” numbered symmetric modes (solid red lines) possess the Lamé modes. These modes have the minimum attenuation at some specific frequencies, which is similar to the feature previously mentioned. However, the Lamé mode vanishes from the “odd” numbered antisymmetric modes (dash-dotted black lines) and the “even” numbered symmetric modes (dash-dotted red lines) due to the participation of the horizontally polarized motions, that is, SH waves. This phenomenon represents that the energy of thermoelastic guided wave will dissipate into the region out of the sagittal plane during propagating along the orientation between the - and -directions. The special regions neighboring the intersecting dispersion curves represent coupling motion of different modes. Therefore, it is of difficulty to interpret the mode-converted response generated at those regions using LIU or PA technique.

(a)

(b)

Figure 7: Frequency spectra of (a) the phase velocity and (b) the semilogarithmic plot of attenuation for the symmetric modes (red lines) and antisymmetric modes (black lines) for the thermoelastic waves propagating at 45° orientation to the -axis in a copper foil under the initial state with uniaxial prestress applied in the -direction.

4. Conclusion

In this paper, a copper foil exerted by a uniaxial tensile prestress in the -direction is considered as an example. Applying the curve-tracing method for the root finding of complex wavenumber, the numerical evidence indicates that the response of thermoelastic waves can characterize the uniaxial residual stresses in the plate-like structures through the frequency spectra of phase velocity dispersion and wavenumber attenuation, especially in the direction of wave propagation parallel or perpendicular to the loading direction. Except for the A₀ mode, the attenuation spectra of thermoelastic waves have steep descents at the specific frequencies where their unique minima occur. The attenuation increases with increasing tensile prestress in the same direction as wave propagation. If the prestress orientation is perpendicular to the direction of thermoelastic wave propagation, the reductions in these specific frequencies of Lamé modes are proportional to the magnitudes of applied stress. Along the perpendicular direction, the phase velocities apparently decrease as the prestress increases. Furthermore, the isotropic material property in the sagittal plane of propagating thermoelastic guided wave can affect the appearance of close-to-zero attenuation.

Acknowledgment

This research was partially supported by National Science Council under Grant no. NSC 98-2221-E-009-007.

References

D. A. Hutchins, “Ultrasonic generation by pulsed lasers,” in Physical Acoustics: Principles and Methods, W. P. Mason and R. N. Thurston, Eds., vol. 18, pp. 21–123, Academic Press, New York, NY, USA, 1988.

C. B. Scruby and L. E. Drain, Laser Ultrasonics: Techniques and Applications, Adam Hilger, New York, NY, USA, 1990.

J. A. Sell, Ed., Photothermal Investigations of Solids and Fluids, Academic Press, London, UK, 1989.

D. Almond and P. Patel, Photoacoustic and Photothermal Science and Techniques, Chapman and Hall, London, UK, 1996.

W. Nowacki, Dynamic Problems of Thermoelasticity, P. H. Francis and R. B. Hetnarski, Eds, PWN-Polish Scientific, Warszawa, Poland, 2nd edition, 1986.

H. W. Lord and Y. Shulman, “A generalized dynamical theory of thermoelasticity,” Journal of the Mechanics and Physics of Solids, vol. 15, no. 5, pp. 299–309, 1967.

A. E. Green and K. A. Lindsay, “Thermoelasticity,” Journal of Elasticity, vol. 2, no. 1, pp. 1–7, 1972.

A. E. Green and P. M. Naghdi, “A re-examination of the basic postulates of thermomechanics,” Proceedings of the Royal Society. London. Series A, vol. 432, no. 1885, pp. 171–194, 1991.

A. E. Green and P. M. Naghdi, “On undamped heat waves in an elastic solid,” Journal of Thermal Stresses, vol. 15, no. 2, pp. 253–264, 1992.

A. E. Green and P. M. Naghdi, “Thermoelasticity without energy dissipation,” Journal of Elasticity, vol. 31, no. 3, pp. 189–208, 1993.

H. Singh and J. N. Sharma, “Generalized thermoelastic waves in transversely isotropic media,” Journal of the Acoustical Society of America, vol. 77, no. 3, pp. 1046–1053, 1985.

J. N. Sharma and H. Singh, “Generalized thermoelastic waves in anisotropic media,” Journal of the Acoustical Society of America, vol. 85, no. 4, pp. 1407–1413, 1989.

J. N. Sharma, D. Singh, and R. Kumar, “Generalized thermoelastic waves in homogeneous isotropic plates,” Journal of the Acoustical Society of America, vol. 108, no. 2, pp. 848–851, 2000.

K. L. Verma and H. Hasebe, “Dispersion of thermoelastic waves in a plate with and without energy dissipation,” International Journal of Thermophysics, vol. 22, no. 3, pp. 957–978, 2001.

K. L. Verma and N. Hasebe, “On the flexural and extensional thermoelastic waves in orthotropic plates with two thermal relaxation times,” Journal of Applied Mathematics, no. 1, pp. 69–83, 2004.

K. L. Verma, “On the propagation of waves in layered anisotropic media in generalized thermoelasticity,” International Journal of Engineering Science, vol. 40, no. 18, pp. 2077–2096, 2002.

V. Salnikov and N. Scott, “Asymptotic and numerical analysis of thermoelastic waves in an isotropic plate,” Journal of Applied Mathematics, vol. 71, no. 3, pp. 409–433, 2006.

Y. -H. Pao, W. Sachse, and H. Fukuoka, “Acoustoelasticity and ultrasonic measurements of residual stresses,” in Physical Acoustics: Principles and Methods, W. P. Mason and R. N. Thurston, Eds., vol. 17, pp. 61–143, Academic Press, New York, NY, USA, 1984.

M. J. S. Lowe, “Matrix techniques for modeling ultrasonic waves in multilayered media,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 42, no. 4, pp. 525–542, 1995.

MATLAB Function Reference Volume 2: F-O, The MathWorks, Inc, 2004

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York, NY, USA, 3rd edition, 2007.

J. R. Davis, Ed., ASM Metal Handbook: Desk Edition, ASM International Handbook Committee, 2nd edition, 1998.

R. A. H. Hamilton and J. E. Parrott, “The third-order elastic constants of quasi-isotropic materials,” Journal of Physics C, vol. 1, no. 4, pp. 829–834, 1968.

K. F. Graff, Wave Motion in Elastic Solids, Oxford University Press, New York, NY, USA, 1975.

J. D. Achenbach, Wave Propagation in Elastic Solids, Amsterdam, The Netherlands, North-Holland, 1973.