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

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 A0 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.


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 0 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.

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, Θ 0 Θ i . The stress component T IJ and the increment of entropy Ξ satisfy the following constitutive relations: where I, J, K, L 1, 2, 3. The physical field variables u K,L and ΔΘ indicate the strain changes and temperature rise caused by the external disturbance applied to the initial state. The material constants c IJKL , λ IJ , and α denote the effective elastic constants, thermoelastic coupling coefficients, and thermal constant influenced by the initial strains u i K,L , respectively. They are defined by 3 denotes the cubic dilatation. The terms c IJKL , c IJKLMN , λ IJ , α, C E , ρ 0 , and Θ 0 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 u i K,L and the initial residual stresses T i IJ is given by T i IJ c IJKL u i K,L .

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 ρ i is the mass density in the initial state and is defined by ρ i 1 − e i NN ρ 0 . The term ρ i b I 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 Figure 2: Schematic diagram of a copper foil prestressed in the X 1 -direction. The unit vector n denotes the propagating direction in the X 1 X 2 -plane of guided thermoelastic waves.
where Θ i Θ 0 , ρ i h, and q I are the temperature, distributed body heat source, and surface heat flux in the initial state, respectively. k IJ is the effective thermal conductivity influenced by the initial strains u i K,L and can be written as where k IJ 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.6

Thin Isotropic Plate Subjected to a Uniaxial Prestress in the X 1 -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 T i IJ are assumed to be homogeneous. Only normal stress is applied in the X 1 -direction, and shear stress components vanish, that is, T i

2.8
According to 2.7 and 2.8 , the special relations c 22 c 33 , c 13 c 12 , c 55 c 66 , λ 2 λ 3 , and k 2 k 3 are obtained. Rewriting 2.6 leads to the thermoelastic governing equations involving the effect of uniaxial prestress T i 11 as follows: assumed to have a plane wave harmonic function e i ξ J X J ζX 3 −ωt based on the partial wave expansion PWE method. Then, the solutions of u I I 1, 2, 3 and ΔΘ are represented as where A I and A 4 indicate the amplitudes of u I and ΔΘ, respectively. The term ξ J J 1, 2 denotes the wave vector of guided wave in the X 1 X 2 -plane and ζ denotes the angular wavenumber in the thickness X 3 -direction. The components of wave vector ξ J are defined by ξ 1 ξ cos θ and ξ 2 ξ sin θ, where θ denotes the included angle of the arrowhead n and the X 1 -axis as shown in Figure 2. The parameters ξ and ω are the angular wavenumber and angular frequency defined as 2π times the wavenumber k and the frequency f, 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:

2.12
According to the existence of a nontrivial solution in 2.11 , the determinate vanishes and the eigenvalues ± ζ k k 1, 2, 3, 4 must satisfy the quartic equation of ζ 2 as follows: where the coefficients a 8 , a 6 , a 4 , a 2 , and a 0 can be obtained and sorted by employing the software MAPLE for symbolic computation. The quartic equation 2.13 can be directly ISRN Applied Mathematics 7 solved by quartic formula and eight complex roots ± ζ k k 1, 2, 3, 4 are subsequently obtained. For the definiteness of a single-value function, a constraint for root ζ k is selected as Im ζ k ≥ 0 to avoid the exponential increasing. The symbol " " corresponds to the downgoing waves propagating along the positive X 3 -direction. On the other hand, the symbol "−" denotes those upgoing waves traveling toward the negative X 3 -direction. The k with respect to ± ζ k satisfy the proportional relation, 14 where C ± k denotes the unknown amplitudes and the proportional factors p ± 1k , p ± 2k , p ± 3k , p ± 4k are the determinants of corresponding submatrices or minors in 2.11 of the form

2.15
Applying 2.1a and 2.4b , the solutions of traction t I ≡ T 3I , I 1, 2, 3 and heat influx q in ≡ − q 3 along the positive X 3 -direction can be represented in the form of where B I and B 4 indicate the amplitudes of T 3I and −q 3 , respectively. Following the above procedure in a similar manner, the components k with respect to ± ζ k also satisfy the proportional relation, where the proportional factors q ± 1k , q ± 2k , q ± 3k , q ± 4k are given by

2.18
Hence, applying the matrix technique 19 used in analysis for the layered structure, a combination of 2.10 and 2.16 yields The submatrices P ± and Q ± are the 4 × 4 matrices with elements p ± ik and q ± ik as given in 2.15 and 2.18 , respectively. The diagonal matrix D ± X 3 has a rank 4 with plane wave forms e ±iζ k X 3 . The uncertain 4 × 1 vector C ± with elements C ± k can be determined by the boundary conditions.

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 X 3 ±h/2 are both traction-free and adiabatic conditions, that is, T 3I q 3 0 I 1, 2, 3 . 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  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 k k r ik i k r 1 iγ ν /2π , where k i and γ ν are the attenuation per unit distance Np/mm and per wavelength Np/ν , respectively. The characteristic equation Ω S,A f, k given in 2.19 can be represented as Ω S,A f, k r , k i or Ω S,A f, k r , γ ν . Moreover, because of the difficulty in obtaining an exact zero solution to Ω S,A f, k r , k i 0 numerically. Instead, minimizing |Ω S,A f, k r , k i | 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 k r and imaginary wavenumber k i with respect to frequency f. There exist many roots of f, k i with respect to one point on the k r -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.

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 "mg", "mm", "μs", and "k • K" and listed in the first column of Table 1. Furthermore, the effective material properties under two initial states with two uniaxial prestresses 0.02c 44 and 0.04c 44 applied in the X 1direction are assembled in the second and third columns of Table 1, respectively. Obviously, according to the relations c 22 c 33 ≈ c 23 2c 44 , λ 2 λ 3 , and k 2 k 3 in the X 2 X 3 -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 T i 11 are quasi-isotropic.
.74 is the heat capacity mg/μs 2 ·k • K , and γ 0.0165 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 ν 0.367 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 1 -, X 2 -, and X 3 -directions in the initial prestressed state are defined as follows: where the superscript " IJK " 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 X 1 -and X 2 -directions and their corresponding differences are given by and they are related to the speeds of transverse S waves along the X 3 -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 X 2 X 3 -plane for whereas 3.8 only provided the trend of speed change for the Lamé modes owing to the orthorhombic property in the X 1 X 3 -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: 0.02c 44 and 0.04c 44 applyied in the X 1 -direction, are assembled in Table 2.
According to the results via the complex root finding, the frequency spectra of the dispersion real wavenumber k r and attenuation imaginary wavenumber k i 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 c ph and semilogarithmic plot of attenuation k i are also shown in Figures 4 a and 4 b , where the phase velocity c ph is defined as f/k r . In Figure 4 a , the four additional blue dashed lines, labeled by "c L0 ," " c S0 ," "c Lame0 ," and "c R0 " 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 0 and S 0 modes converge to a constant value corresponding to the Rayleigh wave speed c R0 1.972 mm/μs.
In Figure 4 b , the attenuation spectrum of each mode, except the A 0 mode, has a closeto-zero minimum at a specific frequency, which is called the "Lamé mode" 24 . The Lamé modes travel at a specific wave speed c Lame0 and can be indicated using a dimensionless parameter k r h n 1/2 for the symmetric modes S n n 0, 1, 2, . . . and k r h m for the antisymmetric modes A m m 1, 2, 3, . . . . Reflecting on the "c Lame0 "-labeled blue dashed line shown in Figure 4 a , the specific frequencies, that is, f c Lame0 /h · n 1/2 for the S n modes and f c Lame0 /h · m for the A m 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 0 and S 0 modes merge together at the frequency range higher than 40 MHz. The convergent value is about 0.5 × 10 −3 mm −1 at 80 MHz.
Let the copper foil be exerted by a tensile stress in the X 1 -direction. Two cases of applied stresses 0.02c 44 and 0.04c 44 are considered. Figures 5 a and 5    Attenuation (Np/mm) Figure 4: Frequency spectra of a the phase velocity c ph and b the semilogarithmic plot of attenuation k i 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 .

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 X 1 -direction. Figures 6 a and 6 b show those of thermoelastic waves propagating along the X 2 -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 Δc 100 L1 and Δc 001 S1 in the X 1 -direction due to the tensile prestress T i −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.
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 0 mode, increases as the tensile prestress in the same direction as wave propagation increases. This phenomenon is caused by the uniaxial prestress T i 11 and leads to the orthorhombic symmetry in the X 1 X 3 -plane, for example, the values c 11 , c 33 , c 13 , c 55 , λ 1 , λ 3 , k 1 , and k 3 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 X 2 X 3 -plane, that is, the specific relations c 22 c 33 ≈ c 23 2c 44 , λ 2 λ 3 , and k 2 k 3 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 X 1 -and X 2 -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 X 1 -axis in a copper foil which is tensilely prestressed by 0.02c 44 in the  Figure 6: Frequency spectra of a the phase velocity c ph and b the semilogarithmic plot of attenuation k i for the symmetric modes red lines and antisymmetric modes black lines of a thermoelastic wave propagating along the X 2 -direction in the copper foil under the natural state without any prestress and two initial states with two uniaxial prestresses 0.02c 44 and 0.04c 44 applied in the X 1 -direction.  A m modes (m is even) with lame mode A m modes (m is odd) without lame mode S n modes (n is even) without lame mode S n modes (n is odd) with lame mode  A m modes (m is even) with lame mode A m modes (m is odd) without lame mode S n modes (n is even) without lame mode S n modes (n is odd) with lame mode b Figure 7: Frequency spectra of a the phase velocity c ph and b the semilogarithmic plot of attenuation k i for the symmetric modes red lines and antisymmetric modes black lines for the thermoelastic waves propagating at 45 • orientation to the X 1 -axis in a copper foil under the initial state with uniaxial prestress 0.02c 44 applied in the X 1 -direction. X 1 -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 0 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 X 1 -and X 2 -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.

Conclusion
In this paper, a copper foil exerted by a uniaxial tensile prestress in the X 1 -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 0 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.