Simulation of Thermoelastic Wave Propagation by Means of a Composite Wave-Propagation Algorithm

A modiﬁcation of the wave-propagation algorithm is used as a tool for determining contact quantities in a ﬁnite-volume scheme for the numerical simulation of two-dimensional thermoelastic wave propagation in inhomogeneous media. The modiﬁcation is needed to provide the satisfaction of the thermodynamic consistency conditions between adjacent elements. It appears that the algorithm is thermodynamically consistent except for the limiter functions. Therefore, a composite scheme is used where the Godunov step is applied after each three second-order Lax–Wendroff steps. Elimination of source terms is made following the method of balancing source terms after independent solution of the heat conduction equation.


INTRODUCTION
The question of conceiving an accurate numerical scheme for thermoelasticity was raised recently in conjunction with the problem of the simulation of the progress of phase-transition fronts in crystalline substances.The latter is a free boundary problem which involves rapid localized changes in the field solution.It is with this problem in mind that the authors have developed a numerical scheme dealing first with the case of materially inhomogeneous thermoelastic conductors with smooth or abrupt property variations but no phase changes.
It should be noted that there are many possibilities for the numerical solution of hyperbolic systems of equations.However, certain assumptions about the smoothness of solution are typically used to approximate derivatives in standard finite-difference methods.These approximations are not valid near discontinuities in the material parameters.Therefore, standard methods often fail completely if the parameters vary drastically on the grid size.By contrast, the recently developed wave-propagation algorithm [1] has been found quite natural for the modeling of wave propagation in rapidly varying heterogeneous media [2].This algorithm combines high resolution with multidimensional wave propagation.Special limiter functions are applied to reduce spurious oscillations near discontinuities.As a result, sharp resolution of shocks along with nearly second-order accuracy of smooth solutions is obtained.
However, the application of limiters depends on the considered problem [2] and still seems to be more of an art than a science.Moreover, the fulfilling of thermodynamic consistency conditions cannot be immediately controlled by using nonlinear wave limiters.The thermodynamic consistency conditions are extremely important for the modeling of phase transitions in solids.It seems, therefore, that the recently proposed composite schemes [3] are more convenient for our purposes, because of the use of filters that are consistent with differential equations.In the present paper, necessary modifications are introduced for the conservation of the thermodynamic consistency of the numerical scheme.
The system of equations for thermoelastic wave propagation, given in the second section, is represented in an integral form which is appropriate for a finite-volume scheme.Contact quantities that appear in this formulation are determined by means of the two-dimensional wave-propagation algorithm [1] as it is briefly described in the fourth section in the example of elastic waves.A thermodynamically consistent composite scheme is obtained by application of the Godunov step after each three second-order Lax-Wendroff steps.In the fifth section, the elimination of source terms in the equations of thermoelasticity is made following the method of balancing source terms [4] after independent solution of the heat conduction equation.Results of computation for certain test problems show the efficiency and physical consistency of the algorithm.

BASIC EQUATIONS OF THERMOELASTICITY
We shall consider the classical thermoelasticity of heat conductors.When the geometrical nonlinearities are neglected, the main two equations of thermoelasticity are the local balance of momentum at each regular material point (in the absence of body force) [5], [6], and the heat conduction equation (in the absence of body force and anelasticity), where ρ 0 = ρ0 (x), S is the entropy, T is the absolute temperature, q is the heat flux vector, and v and σ are the velocity and Cauchy stress tensor, respectively.Here σ, S, and q are given by with a free energy per unit volume given by if the u i are the Cartesian components of the elastic displacement.The indicated explicit dependence on the point x means that the body is materially inhomogeneous in general.For the sake of example, we more particularly consider linear isotropic thermoelasticity for which where the three different contributions that relate to elastic energy, thermal energy, and thermoelastic interaction energy, respectively, are given by Here T 0 is a spatially uniform reference temperature and only small deviations from it are envisaged.Simultaneously, we assume the Fourier law of heat conduction, The more usual dilatation coefficient α is related to the thermoelastic coefficient m, and the Lamé coefficients λ and μ by m = −α(3λ + 2μ).We can then rewrite the relevant bulk equations of inhomogeneous linear isotropic thermoelasticity as the following three equations of which the second one is none other than the time derivative of the Duhamel-Neumann thermoelastic constitutive equation, From the mathematical viewpoint, the problem is to find a solution to the system (6)-( 8) of equations with corresponding initial and boundary conditions.This is conceivably difficult in general; therefore an efficient and robust numerical method is required.Obviously, numerical schemes deal with discrete elements representing the continuous body.The formulation ( 6)-( 8) paves the way for a reformulation of the governing equations in terms of parameters related to the discrete elements.It happens that the latter can be viewed from both numerical and thermodynamic viewpoints.

FINITE-VOLUME METHOD AND THE THERMODYNAMICS OF DISCRETE SYSTEMS
It is salient to remind the reader of the notion of discrete systems in thermodynamics (Schottky [7], Muschik [8]).In such thermodynamics, the thermodynamic state space is extended by means of contact quantities in order to describe nonequilibrium states.In this perspective a discrete system is a domain G of R 3 , which is separated from its environment G * by a partition ∂G.The interaction between G and G * is described by contact quantities.In a Schottky system per se, this interaction consists of heat, work, and mass exchanges.For instance, considering heat exchange Q, the contact temperature is defined by the inequality for vanishing work and mass exchange rates.Here T * is the thermostatic temperature of the equilibrium environment.From Eq. ( 9) it follows that Q and the term in brackets always have the same sign.If we now suppose that exactly one equilibrium environment exists for each arbitrary discrete system for which the net heat exchange between them vanishes, then the defining inequality ( 9) determines the contact temperature of the system as the thermostatic temperature T * of the system's environment for which this net exchange vanishes.The dynamic pressure, p, and the dynamic chemical potential, μ, are defined analogously by The contact quantities so defined provide a complete thermodynamic description of nonequilibrium states of a separated discrete system.Note, however, that the values of the defined contact quantities differ from the values of usual bulk parameters of the case of local equilibrium.For interacting elements, which are the case in our study, the values of bulk and contact quantities of adjacent elements are additionally connected by thermodynamic consistency conditions [9].
In the required extension of the concepts of the thermodynamics of discrete systems to the thermoelastic case, we divide the body in a finite number of identical elements.This is also the strategy of the numerical method of finite volumes.The state of each element is then identified with the thermodynamic state of a discrete system associated with that element; each element is assumed to be in local equilibrium.In thermoelasticity, in addition to and the defining inequality (9), which governs heat exchange, we must define a contact dynamic stress tensor i j since the state space includes the deformation.Analogously to (9), which holds for εij = 0, we have thus an inequality that reminds us of the Hill-Mandel principle of maximal dissipation [10].Here σ * i j is the Cauchy stress tensor in the environment.Now it remains to make the connection between the bulk quantities that appear in Eqs. ( 6)-( 8) and the contact quantities.This is achieved as follows.The contact stress tensor i j is defined at the boundary ∂ V of the volume element V and V i denoting, by duality, the contact deformation velocity at that boundary of unit outward normal n i .The integration over the finite-volume element of Eqs. ( 6)-( 8) and the definition of the strain rate yields the following set of integral forms when contact quantities are substituted for bulk ones, where and source terms resulting from material inhomogeneities (labeled "inh") and thermoelastic couplings (labeled "te") are given by Equation (13) follows from the last of (3) by taking the time derivative and integrating over volume V .As m here does not depend on time, the contribution ϕ te i j could possibly be rewritten as a time derivative and grouped with the left-hand side of Eq. ( 14), leaving only terms originating from material inhomogeneities as source terms.However, we keep the above formalism in order to distinguish between the purely elastic and thermoelastic cases.The reason for this is that in "classical" thermoelasticity, the thermoelastic coupling in Eq. ( 8) is usually considered small and hence negligible, leaving from (8) an equation determining temperature independently of other fields (cf.[5], [6]), even in inhomogeneous materials, that is in that "thermal stress" approximation (15) is reduced to Equations ( 12)-(15)-or (16)-are in a form suitable for a numerical approach by the method of finite volumes.It remains, however, to relate bulk and contact quantities (e.g., q i and i j ).For that purpose, noting that volume elements are cells, thermodynamic consistency conditions have been established [9] (subscripts 1 and 2 denote adjacent cells), where E is the internal energy per unit volume, and the energy of interaction E int is undetermined yet.These conditions provide the control for updating rules for the state of adjacent cells.

TWO-DIMENSIONAL ELASTIC WAVES
The connection between bulk and contact quantities in the elastic case will be established by means of the wave-propagation algorithm [1].The general form for the integral equations ( 12)-( 16) in the considered case is A finite-volume scheme corresponding to these equations in two space dimensions can be represented in the form (the source term disappears for each locally homogeneous elementary volume) where AQ ± and B Q ± are the contact quantities in the horizontal and vertical directions, respectively.At the same time, the system of equations of linear elasticity can be presented in the form of conservation laws.In an inhomogeneous medium in two space dimensions this system of equations is considered in a matrix form The wave-propagation algorithm for the solution of these equations is based on solving Riemann problems at the interface between grid cells [1] .The fluctuations arising from Riemann problems in the xand y-directions, respectively, are determined exactly w ( p+) .
Here α ( p) i j and β ( p) i j are eigenvalues of matrices A i j and B i j , where coefficients depend on the state of the cell (i j); W ( p) and w ( p) are horizontal and vertical waves corresponding to the local Riemann problem; superscripts "+" and "−" denote positive and negative eigenvalues.
The construction of the wave propagation algorithm begins with establishing the firstorder Godunov scheme in terms of these fluctuations [1]: It is easy to see that the connection between fluctuations and contact quantities can be established in a very simple way: Using this connection, we can exploit almost all advantages of the wave-propagation algorithm, namely, high-resolution and multidimensional motion.
In the thermoelastic case, the thermodynamic derivative in the conditions (17) can be expressed in terms of stress components as where i and j are fixed.It is supposed that the introduced contact quantities (capital letters) are connected with the energy of interaction in an analogous way: Therefore, the thermodynamic consistency conditions for two cells numbered 1 and 2 and separated by a coordinate line are σ (1)  i j + (1)  i j − T (1) ∂σ (1)   i j It should be noted that values of the contact stresses, determined by means of (22), satisfy the thermodynamic consistency conditions (25) in their isothermal form.The Godunov's method is extended to a high-resolution method by adding correction terms [1].The form of the extended method is where up i j is the update for a first-order upwind Godunov method.The second-order correction terms take the form [1] Fk w ( p)   and provide second-order accuracy.The obtained algorithm is a variant of the well-known Lax-Wendroff method [11].The thermodynamic consistency conditions (25) remain satisfied in this case.The multidimensional motion is accomplished by splitting each fluctuation A * q i j and B * q i j into two transverse fluctuations [1], which will be called B + A * q i j (the up-going transverse fluctuation), B − A * q i j (the down-going transverse fluctuation), A + B * q i j (the right-going transverse fluctuation), and A − B * q i j (the left-going transverse fluctuation).
The total effect of vertical transverse propagation at the interface between cells (i − 1 j) and (i j) is determined as Analogously, the total effect of horizontal transverse propagation at the interface between cells (i j − 1) and (i j) is The introduction of the transverse fluctuations improves the stability limit up to Courant number 1 and fulfills the thermodynamic consistency conditions.
The whole algorithm is implemented in the form where up i j is the update for the first-order upwind Godunov method, and trans i j represents the effect of transverse fluctuations, It is well known that the Lax-Wendroff scheme produces oscillations behind discontinuities [3].The usual way to reduce spurious oscillations is to introduce limiter functions to modify the second-order corrections near discontinuities.Here, instead of limiters, the recent idea of using filters that are consistent with differential equations [3] is applied to provide the satisfaction of thermodynamic consistency conditions.
The corresponding composite scheme is obtained by application of the Godunov step after each three second-order Lax-Wendroff steps.Obviously, the thermodynamic consistency conditions remain satisfied at each step.

Boundary Conditions
Boundary conditions are specified in terms of bulk quantities for boundary elements.For stress-free boundaries, the bulk stresses are set equal to zero.The contact quantities are not specified at the outer boundaries.

Numerical Tests
We have performed a number of simulations to estimate the correctness and capabilities of the algorithm.They include, in particular, the following: • Elastic wavefronts from point source at the boundary of a homogeneous medium; • Elastic wave propagation in a layered medium; • Elastic waves in a medium with periodically and randomly distributed inclusions.
All calculations were performed by means of the composite Lax-Wendroff-Godunov scheme with the Courant number equal to 1. Waves were excited by the prescription of nonzero normal component of the stress tensor within a finite aperture at the boundary for a short time.
As the first example, a short-time symmetrical excitation in the middle of a boundary of a homogeneous rectangular domain by sinusoidal normal stress is considered.All other boundaries are stress-free.Physical parameters for aluminium were used in the calculation: c p = 6420 m/s, c s = 3040 m/s, ρ 0 = 2700 kg/m 3 .
Figures 1 and 2 represent snapshots of wavefronts after 50 and 250 time steps, respectively, in terms of total displacement.It seems that the proposed algorithm gives a good correspondence with the almost classical situation.Figure 3 shows a snapshot of the propagation of elastic wave, which is excited at a part of the bottom boundary with linear retardation, in a layered medium (corresponding density distribution is shown in Fig. 4).In the inhomogeneous case, we need to prescribe additionally the physical parameters for copper: c p = 4560 m/s, c s = 2600 m/s, ρ 0 = 8960 kg/m 3 .The picture is asymmetric because of asymmetric loading at the boundary.The same medium is used for the performing of standard experiment with refining of the spatial grid.The results of calculations are presented in Fig. 5 in terms of normal stress distribution along the centerline.The density distribution is also shown schematically.One can see that refining the grid gives more detailed stress distribution but does not change the qualitative results.Figures 6 and 7   wave in a medium with periodically distributed inclusions of another material.For randomly distributed inclusions, as shown in Fig. 8, we have to start with a more asymmetric picture, but then some regularity appears with the larger part of the signal finally propagating quite symmetrically and practically aligned with the initial source, as can be seen in Fig. 9.This is typical of the "isotropic averaging effect" of random distributions.

THERMOELASTIC WAVES
In the thermoelastic case, the equation for the temperature in two dimensions, is solved by means of a simple algorithm for the heat conductivity equation in an inhomogeneous medium [12], where B = λ + 2/3μ is the bulk modulus.According to the method of balancing source terms [4], the modified bulk velocities in both xand y-directions, are used for the calculation of normal dynamic stresses and corresponding contact velocities, which allows us to eliminate the source terms.The computational derivative for the temperature is calculated as When the contact temperatures are defined as where superscripts "r," "l," "u," and "d" denote left, right, up, and down boundaries, respectively, we can fulfill the thermodynamic consistency conditions (25).The first-order Godunov method as well as transverse propagation, second-order correction, and composition are then applied as above.
Results of simulation for thermoelastic wave propagation in inhomogeneous media are presented in Figs.10-13.Figure 11   with an interface between two distinct, but otherwise spatially homogeneous, thermoelastic materials (copper and aluminium).The step-wise distribution of material properties is shown in Fig. 10.The wave was excited by a purely thermal shock at a part of the bottom boundary, as shown by a narrow black rectangle in Fig. 10; all other boundaries are stress free.Figure 11 shows the first time the wave hits the interface.Figure 13 shows snapshots of the mechanical trace (normal stress) of a thermoelastic wave inside a medium with laterally continuously varying properties (as in Fig. 12).The initial excitation is the same as for the previous case.The wave front curves as it propagates, because of lateral inhomogeneity.Animated pictures can be obtained via e-mail by request.

CONCLUDING REMARKS
The concept of contact quantities which was introduced for the thermodynamic description of the nonequilibrium state of elements in the thermodynamics of discrete systems [8] is successfully used for the construction of a finite-volume numerical scheme for thermoelatic wave propagation in inhomogeneous media.The contact quantities appeared naturally in the integral balance laws for thermoelasticity; the laws serve as a basis for the formulation of the finite-volume algorithm.At the same time, these quantities can be easily placed in correspondence with the fluctuations that appear in the formulation of the wave-propagation algorithm [1].This allows us to transform the wave-propagation algorithm into a thermodynamically consistent scheme using the idea of composition proposed in [3].The modification obtained conserves practically all advantages of the wave-propagation algorithm, namely, high resolution, multidimensional motion, and stability up to a Courant number equal to unity.
We assumed that the energy of interaction is connected with contact quantities in the same way as the usual elastic energy and bulk quantities are connected, because only elastic interactions are considered.The thermodynamic consistency conditions introduced for adjacent elements provide the equivalence of the thermodynamic description at each level of refinement of the mesh.In the present paper, the conditions are fulfilled automatically for the composite scheme used, because no entropy production at interfaces appears in the case of linear thermoelastic waves.Within the composite wave-propagation algorithm, every discontinuity in parameters is taken into account by solving the Riemann problem at each interface between elements.The reflection and transmission of waves at each interface are handled automatically for the inhomogeneous media considered (without entropy production at these interfaces).In fact, it is shown that the algorithm is thermodynamically consistent.
The thermodynamic consistency conditions also allow extension of the proposed algorithm to inhomogeneities such as moving phase transition fronts where entropy production takes place.In such a case-with dissipative phase transition fronts-the conditions mentioned above cannot be fulfilled automatically (especially in their combination that provides the criterion of progress of the fronts), and they are then used as additional relations at the interface.In such more complex situations, both the entropy balance equation and the expression for the driving force, which follows from the balance equation for the pseudomomentum, must be used.However, this extension is presently beyond the scope of the paper.It corresponds to work still in progress.

FIG. 1 .
FIG. 1. Contour plot of elastic wavefronts from a point source in homogeneous medium, Courant number 1.0, 50 time steps.

FIG. 2 .
FIG. 2. Contour plot of elastic wavefronts from a point source in homogeneous medium, Courant number 1.0, 250 time steps.
FIG. 10.Step-wise distribution of material properties inside the computational domain.