Numerical modeling of thermoelastic generation of ultrasound by laser irradiation in the coupled thermoelasticity.

Laser-generation of ultrasound is investigated in the coupled dynamical thermoelasticity in the presented paper. The coupled heat conduction and wave equations are solved using finite differences. It is shown that the application of staggered grids in combination with explicit integration of the wave equation facilitates the decoupling of the solution and enables the application of a combination of implicit and explicit numerical integration techniques. The presented solution is applied to model the generation of ultrasound by a laser source in isotropic and transversely isotropic materials. The influence of the coupling of the generalized thermoelasticity is investigated and it will be shown, that for ultra high frequency waves (i.e. 100GHz) generated by laser pulses with duration in the picosecond range, the thermal feedback becomes considerable leading to a strong attenuation of the longitudinal bulk wave. Moreover, the coupling leads to dispersion influencing the wave velocities at low frequencies. The numerical simulations are compared to theoretical results available in the literature. Wave fields generated by a line focused laser source are presented by the numerical model for isotropic and for transversely isotropic materials.

Introduction

Experimental application of laser-generated ultrasound in solids has a long and successful history. It has been extensively used in experimental acoustics in the past decades since it enables contactless excitation of different elastic waves, such as bulk, surface or guided waves. In search of better understanding and optimization of the generation mechanism, the modeling by a thermoelastic source has attracted great attention from the very beginning of the experimental application. Thus, the generation process is described by the coupled heat and wave equations, whereby various coupling terms are taken into account. One of the earliest work [1] has introduced a point-source representation of the thermal expansion as the so-called surface center of expansion. The corresponding elastic problem was solved by a spatial–temporal integral transform. Similar solutions have been pursued also by later authors applying an equivalent dipole loading in a pure elastic problem [2] or investigating the influence of the thermal diffusion on the generated bulk longitudinal wave [3]. These solutions have been extended to layered media [4] and to line sources in isotropic and transversely isotropic media [5,6]. For such analytic solutions, however, certain simplifications are required, such as a step-like temperature rise with a point or line spatial distribution for the surface heating. In general cases with more realistic distributions of the heat flux the inverse integral transforms remain a great challenge and consequently numerical inversion techniques were pursued. This approach was applied to investigate epicentral waveforms in [3,7,8], the generation process by line focused sources [9] or the influence of the hyperbolic term in the heat conduction equation [10]. Simplifications in the numerical inversion techniques are, however, controversial [11]. An alternative way to solve this coupled problem is shown in [12,13] using the reciprocity theorem with applications to a half space and plates.

Within laser ultrasonics the coupling of the heat and wave equations is due to thermal expansion leading to a semi-coupled problem. In the generalized thermoelasticity, however, also the thermal feedback of the propagating stress pulses is included, hence, the problem becomes fully coupled and both thermal and elastic waves arise [14]. Moreover, the coupling of the two fields influences the elastic wave velocities leading to additional dispersion and strong attenuation of the waves in isotropic [15,16] or in anisotropic media [17,18]. Although, the influence of the coupling on the wave velocity remains limited [15] the attenuation of pulses increases with frequency and it will be shown, that for laser pulses with duration of picoseconds or less the excited waves posses ultra high frequencies (∼100 GHz) whereby this attenuation becomes considerable.

Numerical methods, such as finite elements or finite differences have been rarely pursued to model laser-generation, though they offer a greater flexibility in their applications. A combination of an implicit-explicit finite difference technique has been applied in the one dimensional case to investigate the generation of bulk longitudinal waves [19] using the so-called magic time step [20]. The laser-generation of surface waves and bulk waves in coated and uncoated half spaces have been investigated using the finite element method in [21,22]. This solution applies a direct integration technique, though, the transient response could also be obtained by a combination of an Eigenvalue decomposition of the homogeneous equation and a harmonic analysis with an inverse Fourier-transform [23]. Numerical solutions are also capable to introduce full coupling as it was shown in [24] where a time domain finite element method was applied to study second sound effects due to heat propagation.

In the presented work laser-generation of ultrasound is investigated in the coupled dynamical thermoelasticity. In the presented numerical solution with staggered grids the two sets of equations decouple and therefore a combination of implicit and explicit finite difference techniques are applied. The hyperbolic heat conduction equation is solved by the implicit Wilson Θ integration technique and the wave equations by the explicit Euler method. Due to the combination of implicit and explicit techniques the same spatial and temporal discretisation are used for both sets of equations whereby the numerical efficiency remains high through the application of the explicit integration for the wave equation. The influence of coupling is investigated for line-focused pulsed laser sources with various pulse rise times. It will be shown, that for ultra high frequency waves (∼100 GHz) the thermal feedback becomes considerable through strong attenuation. Moreover, the coupling increases the wave velocities for lower frequencies by introducing dispersion. Wave fields generated by the line focused laser source are presented by the numerical model for isotropic and also for transversely isotropic materials.

Governing equations

The governing equations of the generalized thermoelasticity in the linear elastic case are given by the coupled heat conduction and wave equations describing the coupling of displacement and temperature fields. In the presented work a 2D case with plane-strain boundary conditions is investigated, which is a good approximation of a line-focused laser source [9]. The hyperbolic heat conduction equation is given as:where T denotes the temperatures, K the thermal conductivity, c the specific heat of the material at constant deformation, τ the material relaxation time, T0 the reference temperature, q the external heat flux, ρ the density, the normal strains and the thermal moduli. The wave equations are given as:where ux, u denotes the displacement components and are the stress components. In the presented work we consider an orthotropic material law with nine independent elastic constants and Hooke’s law is replaced by the Duhamel–Neumann equations [25, p. 421–422] for the linear elasticity, which considers the effect of thermal expansion:

According to the assumed orthotropy, the linear thermal expansion is described by three independent constants with three independent thermal moduli [26, p. 25–26]. The thermal moduli relate the stiffness matrix and the linear thermal expansion coefficients as:The heat conduction equation in Eq. (1) is given with its hyperbolic form, which assumes a finite heat propagation speed c, commonly approximated as five to ten times the longitudinal wave velocity c [9,10]. Here, we will assume that ce = 5c which leads to the following thermal relaxation time [10]:Although the hyperbolic term in the presented heat equation becomes significant only for laser heating with duration in the fs–ps range, the validity of this equation for such ultrashort heating has been questioned [27]. During the heating process the free electrons absorb the radiated energy leading to a subsequent heating of the metal lattice through collisions. This non-equilibrium heating process is described by the two-temperature model [27] which assigns different temperatures to the electrons and to the lattice. According to previous results [19,27] this effect becomes significant for heating shorter than 1 ps (assuming room temperatures). To avoid implications by this effect we will restrict our investigations to cases where the laser heating is significantly longer than 1 ps; hence, the shortest pulse duration considered in the following is 20 ps.

Explicit–implicit finite difference approximation of the coupled equations

The spatial discretisation of the two sets of equations in Eqs. (1)–(6) is done with finite differences using staggered grids. Hence, the temperature, displacement and stress components are discretised on different grids which are shifted by a half cell in one or in both directions, according to Fig. 1. The grids of the temperatures and normal stresses are identical and the grids of the displacements and shear stresses are shifted by a half grid dimension ().

Fig. 1

Staggered grid for the discretisation of the coupled thermoelastic equations. The temperatures T and the normal stresses σ, σ are discretised on an identical grid; the grids of the displacements u, u and shear stresses σ are shifted by a half cell dimensions .

The presented form of the coupled equations with stresses as additional variables (heterogenous solution [28,29]) is extensively applied in numerical solutions [19,30,31] since it requires only the first spatial derivatives to be discretised in the wave equation. The heterogenous, staggered grid formulism leads to further advantages in the investigated case due to decoupling of the discretised solution, as will be discussed later. Further details of discretisation with finite differences on staggered grids are not presented here, they can be found in [30-32], for example. Stress free boundary conditions to the wave equations and adiabatic boundary condition to the heat conduction equation are applied [19] and due to the symmetry of the problem only one half of the area is modeled using symmetric boundary conditions:

The spatially discretised equations could be written in the following matrix form:where M, C, K are the coefficient matrices. The vectors and f are given as:where q denotes the heat flux and are the first order finite differential operators in x and y directions, respectively. The nodal temperatures T, displacements and stresses are collected in vectors. The coefficient matrices M, C, K are given as follows:where E denotes the identity matrix and B is a coefficient matrix given in Appendix B. Due to the heterogeneous formulism in Eqs. (2)–(6) the discretised heat and wave equations are effectively decoupled and could be treated separately. Using a homogeneous solution, i.e. by inserting Eqs. (4)–(6) into Eqs. (2) and (3), off-diagonal elements would appear in K leading to a fully coupled set of equations. In the current uncoupled case the wave equation is discretised by the conditionally stable, explicit Euler-method; therefore, the displacements in the (n+1)th time step are explicitly solved as a function of the displacements u and stresses in previous time steps [30-32], as given in Appendix A.

Eq. (10) reduces therefore to the heat equation:where the coefficient matrix B is given in Appendix B. The strains appear on the right side of the equation as a ”thermal load” since they are calculated from the explicitly (and independently) integrated wave equation.

A conditionally stable, explicit, time domain discretisation of Eq. (15) would require a much smaller time step than the one used in the wave equation. Unconditionally stable, implicit integration techniques, such as the Wilson Θ method [33, p. 777] are available, and they enable the use of the same temporal and spatial discretisation as in the wave equation. Eq. (15) is therefore integrated by the Wilson Θ method [33, p. 777], described in Appendix C.

During the numerical solution, Eq. (15) is evaluated by the Wilson Θ method (Appendix C) and the resulting temperature field is applied in the wave equation as excitation. The wave equation is than solved explicitly (Appendix A). Due to the two coupling terms (thermal expansion and thermal feedback) both sets of equations are updated in every time step.

Numerical examples

The discretised equations are applied to simulate the generation process by a pulsed laser with a pulse rise time of 20 ps. In the simulations an isotropic aluminum alloy1 and a transversely isotropic zinc2 are investigated within an area of 6 μm × 6 μm3. In the simulation for zinc the free surface of the half space coincides with the isotropic plane of the crystal; hence the modeled plane is perpendicular to it. The reference temperature will be assumed as 293°K. Zinc has a lower reflectivity than aluminum; for a better comparison, however, we assume that the same amount of heat is absorbed at the surface of the material in both cases and we assume a reflectivity of 0.95 (5% of the pulse energy is absorbed) for both materials.

Excitation – surface heat flux

During the irradiation of a metallic surface with a short laser pulse the electromagnetic waves penetrate a thin skin close to the surface (in the range of ∼10 nm [36]) which could be modeled as a volume heating or approximated as a surface flux [9]. In the current work the half width of the applied laser source (∼1 μm) is much larger than the penetration depth and the surface flux approximation will be applied. A suitable spatio-temporal pulse shape is given as a Gaussian pulse [9,23]:where υ and R denote the temporal and spatial half widths of the pulse, respectively, and A is the magnitude of the pulse. To obtain realistic results we consider 50 μJ pulse energy focused to a line with a length of L = 10 mm and a width of 2R = 2 μm; in combination with the reflectivity this leads to an amplitude of A = 125 Jm-2.

In the presented heat and wave equations in Eqs. (1)–(6) the simplest orthotropic material law was assumed with thermal and material properties independent of the temperature. Temperature dependent thermophysical parameters using linear temperature terms were incorporated in [21] for aluminum using the same reflectivity of 5% as in the current work. The investigation has shown that with incident irradiation energy densities in the range of 10 kJm−2 (3.5 mJ pulse energy with a spot radius of 300 μm) this effect is negligible; for higher energy densities in the range of few hundred kJm−2, however, it becomes significant and the generated waveforms are strongly distorted. In our case the energy density is 2.5 kJm−2, hence we neglect the temperature dependence of the thermophysical parameters.

Validation, convergence

The numerical technique was validated by evaluating the dispersion relations of the generated bulk longitudinal waves along the epicentral axis in the uncoupled case. The spatio-temporal Fourier transform of the simulated data is shown in Fig. 2. The corresponding theoretical wave velocity for the uncoupled case in aluminum is given as and the quasilongitudinal wave velocity in zinc as . The theoretical dispersion relations are marked with circles and show an almost perfect agreement with the simulation. In the isotropic case the theoretical and simulated results agree well up to 100 GHz, though, for zinc the curve starts to deviate at ∼100 GHz as a consequence of the numerical dispersion.

Fig. 2

Validation of the simulations utilizing the generated bulk longitudinal wave for aluminum (a) and the quasi longitudinal wave in zinc (b). The circles denote the theoretical dispersion relation.

To evaluate the frequency content of the generated bulk longitudinal wave, time domain signals on the epicentral axis are captured, shown in Fig. 3a. The sharp, positive peaks in the signals are the generated bulk waves (or quasi longitudinal wave in zinc), often described as the precursor. The generated wave in zinc has a larger amplitude as the result of the higher linear thermal expansion coefficient of zinc along the c axis of the crystal. The frequency contents of these pulses are evaluated by a temporal Fourier transform (Fig. 3b) after multiplication with a Hanning window. The excited frequency spectrum in Fig. 3b reaches above 100 GHz in both cases, whereby the corresponding wavelengths are Λ ∼ 50–60 nm. The full width at half maximum is, however, much higher for aluminum (∼35 GHz vs. ∼20 GHz in zinc). Spatial counterparts of the time domain signals are shown in Fig. 3c as epicentral wave forms (displacements on the epicentral axis) for both materials.

Fig. 3

(a) Time domain signals on the epicentral axis with the quasi longitudinal wave (precursor) for zinc (solid line) and longitudinal bulk wave for aluminum (dotted line), generated by a 20 ps laser pulse. (b) Fourier-transform of the signals, with frequency contents up to 100 GHz. (c) Waveforms on the epicentral axis for both materials. (d) Comparison of three discretisations for aluminum to test the convergence of the technique.

For numerical simulations it is convenient to investigate the convergence or in case of transient solutions the influence of numerical dispersion [20], which is crucial for broadband pulses. The simulation in the isotropic case was repeated with three different discretisations for the same simulated area. For the original simulation the vertical dimension of the cell (the propagation direction of the bulk wave) was chosen as Δy = 4.43 nm. For comparison, simulations with Δy = 7.7 nm and Δy = 12.9 nm were done (keeping the aspect ratio of the cell constant) and the displacements on the epicentral axis were compared for aluminum in Fig. 3d at t = 495 ps. In the original case the discretisation at 100 GHz (Λ ∼ 60 nm) corresponds to ; for the additional simulations this ratio is reduced to ∼8 and ∼5, respectively. The presented displacements on the epicentral axis in Fig. 3d shows strong numerical dispersion for the poorer discretisations, visible due to the broadening of the pulse and the oscillation after the main pulse. Hence, for the investigations of the thermal feedback the originally chosen discretisation is sufficient but also necessary. This ratio is also in good correlation with the values suggested in literature [30-32].

Generated waveforms

First, we would like to visualize the generated waves in the near field, in the close vicinity of the source. The temperatures and the generated wave fields are presented in Fig. 4 at three different times. The upper row shows the results for aluminum, the lower row for zinc. At ∼0.06 ns after excitation the bulk longitudinal wave (quasilongitudinal mode for zinc) generated by the very short light pulse is already visible for both materials. The surface has already a large out-of-plane displacement due to the surface heating. At ∼0.16 ns after excitation the bulk longitudinal wave is well separated and a small compressed region becomes visible between the expanding surface and the propagating bulk wave, interpreted as the beginning of the surface wave. At ∼0.64 ns the generated Rayleigh waves become visible for both materials, however, they are still not well-separated from the region where surface heating and diffusion dominate; here, the vertically expanding area below the surface heating is the result of the presence of the compressive, unipolar Rayleigh waves. Their shape is clearly different in zinc as the result of the anisotropy. A comprehensive description of the generation process by a line focused laser in isotropic medium is given in [9], therefore, in the following we will focus on the effect of the thermal feedback.

Fig. 4

Snapshots of the generation process through temperatures (left) and magnitudes of displacements (right) for aluminum (upper row) and zinc (lower row) for three different times. At 0.06 and 0.16 ns after excitation the bulk wave is already visible. The generated Rayleigh wave becomes visible only after ∼0.64 ns.

Influence of the thermal feedback

Coupling, or thermal feedback, is the result of volume changes [37, p. 394], occured by normal strains. We expect the largest effect on the generated bulk longitudinal or quasilongitudinal waves along the epicentral axis and focus our investigations on these waves. A direct comparison of the displacements on the epicentral axis in Fig. 5a and b for simulations with and without thermal feedback shows a broadening of the pulses with a decrease in amplitude. Although, the effect is considerable for the isotropic aluminum, it becomes stronger for zinc. This attenuation is predicted in [37] as an effect of the thermal feedback. The propagating stress wave results in a small thermal disturbance along its propagation converting the mechanical energy into heat, hence the observed attenuation. Attenuation of the pulses on the epicentral axes are shown in Fig. 5c and d. In the simulation a pulse width of 2 μm was applied, but this width is strongly widening for larger propagation depths (Fig. 4c) due to leakage into the bulk. Hence, geometrical attenuation appears also for the uncoupled cases. For the quasilongitudinal wave in zinc the influence of the thermal feedback becomes particularly strong as shown in Fig. 5d.

Fig. 5

(a)–(b) Influence of the thermal feedback in the isotropic case (aluminum) and in the transversely isotropic case (zinc). (c)–(d) Variation of the waveforms on the epicentral axis at different times emphasizing the attenuation. (e) Time domain signals on the epicentral axis in different depths. (f) Fourier transform of the time domain signals in (e). (g) Influence of the thermal feedback on the generated Rayleigh waves. (h) Theoretical variation of the coupled longitudinal wave velocity and its attenuation over frequency.

The same trend is shown for the time domain signals on the epicentral axis in Fig. 5e. The Fourier transforms of these signals in Fig. 5f reveal, that the attenuation is strongly frequency dependent, visible for zinc, where the spectrum strongly narrows after 5.5 μm propagation depth and frequencies over ∼40 GHz disappear. For the surface waves (Fig. 5g) the influence is reduced, since the excited frequencies are considerably lower.

A further feature, visible in Fig. 5a and b, is the difference between the wave velocities in the coupled and uncoupled solutions. This effect is predicted by theoretical results in the coupled thermoelasticity. In the isotropic case the dispersion relation for the coupled longitudinal wave is given as [37, p. 394]:where ε denotes the thermoelastic coupling constant and C11 is an element of the stiffness tensor. Although, the presented equation is based on the parabolic heat equation, it approximates well the heat conduction for the investigated frequencies. This non-linear equation could be solved for the unknown phase velocity c by setting . The roots are complex and choosing the root for the elastic wave the dispersion relation could be evaluated (Fig. 5h). The phase velocity in the coupled thermoelasticity is therefore higher than the pure elastic wave velocity (cl = 6119 ms-1) and the waves are strongly attenuated, according to the negative imaginary part of the wave velocity in Fig. 5h. It has been shown [14,37], that the phase velocity is convergent for low and for high frequencies:

The influence of the coupled thermoelasticity remains therefore strong for both low and high frequencies. At low frequencies the phase velocity is increased with low attenuation and at high frequencies the phase velocity converges to the compressional wave velocity but it is strongly attenuated. The difference in the wave velocities in Fig. 5a is ∼55 ms-1, calculated from the different arrival times of the pulses, which is a good medium value for the theoretical variation of ∼20–80 ms-1 in Fig. 5h.

The thermal feedback in Eq. (1) is based on Kelvin’s formula [25, p. 459] and governed by the strain velocity and the strain acceleration, in our case:whereby T denotes the absolute temperature as the sum of the reference temperature and the disturbance . Assuming that we obtain after integration:The generated temperatures by the mechanical feedback depend upon the propagating strains and strain velocities which are a function of the frequency and therefore pulse rise time. This influence was investigated for different pulse rise times (20 ps, 50 ps, 100 ps, 250 ps and 500 ps) keeping the pulse energy constant. The results are shown in Fig. 6. Kelvin’s formula in Eq. (20) states that the generated thermal feedback is governed by the strains. Hence, for short pulses with very quick rise times the effect is stronger, as shown in Fig. 6a. A better comparison is obtained by eliminating the amplitude difference through a normalization as shown in Fig. 6b. With increasing pulse width the attenuation decreases, as predicted by the theory. Although, the attenuation from the thermal feedback disappears above ∼100 ps pulse rise time, a difference of ∼5% between the amplitudes of the coupled and uncoupled solutions remains such as the small time shifts between their arrival times.

Fig. 6

Influence of the pulse rise time on the thermal feedback. (a) The attenuation is decreasing with increasing pulse rise time and above approximately 100 ps becomes negligible. (b) Normalized amplitudes of the generated pulses; the deviation in the wave velocity remains. (c) Comparison of the influence of the parabolic and the hyperbolic terms; the term from the classical theory dominates. (d) Temperatures generated by the thermal feedback with different pulse rise times.

Although, the thermal feedback becomes significant for the 20 ps pulse rise time, the hyperbolic term remains negligible in the backward coupling, as shown in Fig. 6c by a comparison of the two terms () of the thermal feedback in Eq. (20). Here, the second term is magnified by a factor of 50 for a better visibility. The generated temperatures by the thermal feedback are presented in Fig. 6d; they remain below 4 K and strongly decrease with increasing pulse rise time.

Surface waves and far field

Finally, the generated far field wave forms will be presented for the isotropic and for the transversely isotropic materials generated by the line focused laser source. Here, a larger area 30 μm × 30 μm was used in the simulations which is discretised by 2418 × 2418 cells.4 Numerical results represented by the vertical surface displacements are shown in Fig. 7a and b. The generated Rayleigh waves are unipolar with no attenuation, as expected. The small, attenuating wave in front of the Rayleigh wave is usually assigned to the longitudinal bulk wave, or ”surface skimming longitudinal wave” [9,38]. The magnitude of the displacements are shown in Fig. 7c at t = 4.48 ns (aluminum) and t = 5.8 ns (zinc) after excitation. Besides the generated bulk longitudinal (L) and shear waves (S) also surface waves (R) are visible in aluminum. The longitudinal bulk wave shows an almost cylindrical wave front in the isotropic case with strongly varying pulse widths (and frequency content), as a consequence of the finite width of the excitation. The wave with a skewed angle to the surface is identified as the leaky-Rayleigh wave (LR) which is predicted for isotropic materials with a high Poisson’s ratio [39].

Fig. 7

Out-of plane surface displacements with Rayleigh and leaky-Rayleigh waves for aluminum (a) and for zinc (b). (c) Sections of the generated wave fields in the models for aluminum (left) and for zinc (right). Besides the generated bulk longitudinal (L), quasi longitudinal (QL) and shear waves (S) also Rayleigh (R) and leaky-Rayleigh (LR) waves are visible.

The wave forms in zinc are different due to the anisotropy. In the transversely isotropic case in the investigated plane there exist three different wave types [34]: quasi longitudinal (QL), quasi transverse (QT) and pure transverse (PT) modes. In Fig. 7c the QL and the QT modes could be detected. At the surface Rayleigh wave is visible, whereby its shape strongly diverges from the isotropic one. Similarly to the isotropic case, the longitudinal bulk wave is disappearing in close vicinity of the surface therefore, the observed surface wave is identified as a leaky-Rayleigh wave.

Conclusion

Laser-generation of ultrasound in the coupled thermoelasticity has been investigated in the presented paper by means of numerical techniques. A numerical model, based on the finite difference method has been introduced. As a result of the application of spatially staggered grids the two sets of discretised equations are decoupled, facilitating the use of a combination of explicit and implicit temporal integration techniques. Here, the heat conduction equation is integrated implicitly and the wave equation by an explicit technique. The number of the coupled linear equations is therefore, only one-third of the total degrees of freedom of the system. The model has been applied to investigate the influence of the coupling (i.e. thermal feedback) on the propagating stress pulses generated by picosecond laser pulses.

It has been also shown, that the coupling influences the wave velocities, by introducing dispersion into an originally non-dispersive system and the attenuation due to the coupling becomes considerable for ultra-high frequency pulses (∼100 GHz). These result are in accordance with the theoretical results, which predict dispersion and increased attenuation for elastic waves as a consequence of the coupling. By varying the pulse rise time it has been shown, that for bulk waves generated by pulses with rise time longer than ∼100 ps the attenuation due to the coupling becomes negligible. Although, the difference in the wave velocities between the coupled and uncoupled cases remains also for longer pulses with lower frequencies, in accordance with the theoretical predictions.

In addition, the model has been applied to investigate the generated waveforms by line-focused pulsed laser sources in isotropic and transversely isotropic media. Snapshots of the generated wave fields were presented revealing the presence of bulk shear and longitudinal and surface waves. The presented results show the applicability of numerical methods to solve and investigate the coupled thermoelasticity with various coupling terms.