1 Introduction

Various kinds of solids, including rocks, plastics, and biological samples, exhibit viscoelastic behavior. This induces dissipation and an irreversible source for temperature increase. Inhomogeneous temperature caused by this and by external thermal effects—like thermal contact with a hydraulically operated loading device, heat transfer and thermal radiation arriving from the environment etc.—initiates heat conduction. Thermal expansion couples heat conduction with the viscoelastic mechanics via the kinematic changes, and it is this resulting coupled problem that one has to solve when wave propagation and other time-dependent processes in such media are to be described with sufficient accuracy. This ‘sufficient accuracy’ may be a quite demanding requirement, as for example in the problem of determining Newtonian noise emerging from the gravitational effect of vibrations (mass rearrangements) in the surrounding rock medium around underground gravitational wave detectors like KAGRA and the planned third-generation Einstein Telescope.

In order to obtain solutions, one way is to utilize available commercial software. Results of such studies can be found in [1, 2] where considerably simpler problems, elastic propagation of a stress pulse in one and two spatial dimensions, respectively, have been solved via numerous different built-in solution methods of the finite-element software COMSOL. The results were disappointing: some methods produced unstable outcome irrespective of settings, some others led to considerable dispersion error (artificial wavy patterns), some produced remarkable dissipation error (artificial damping), and the more-or-less acceptable-looking solutions were considerably different so there appeared a need for some tools to quantify the correctness of the solution and thus to tell from the solutions which may be closest to the correct one. For a more complex problem with numerous fields and various couplings, this path is more discouraging.

Another option is to develop an own numerical method, in a framework that offers insight to the solution’s quality in various ways, enabling monitoring and controlling. It also should be reliable for large times (many wave oscillations), artificial dissipation should be small as it should not be confused with true—in our present case, viscoelastic and heat conduction related—dissipation, and dispersion error should also be small and distinguishable from the true oscillations.

Choosing this latter option, in [1], in one space dimension, large-time reliability of elastic wave propagation was ensured by the symplectic Euler method, reinterpreted as a space and time staggered finite-difference scheme, which increased the first-order accuracy of the symplectic Euler method to second-order. Notably, the reason why a symplectic method performs outstandingly even over many time steps is that it proves to be the exact integrator of a nearby other reversible (Hamiltonian) system (see, e.g., [3, p. 343]). Also in [1], this symplectic scheme was generalized for the Poynting–Thomson–Zener (PTZ) viscoelastic model (see Fig. 1). Stability analysis and investigation of the dispersion relation have revealed the best settings for having stability, no dispersion error, and minimal dissipation error. The generalization preserved second-order accuracy of the scheme both in space and in time, and the finite-difference scheme remained explicit, making simulations run fast, with low computational resource demand.

It is this numerical setup that has been used, in [4], for numerically identifying the PTZ signal propagation velocity, finding good agreement with the corresponding analytical results.

Fig. 1
figure 1

Left: rheological circuit representation of the classic (zero thermal expansion) PTZ model \(\sigma + \tau \frac{\partial \sigma }{\partial t} = E \epsilon + {\hat{E}} \frac{\partial \epsilon }{\partial t}\) between the stress-related quantity \(\sigma \) and the strain-related \(\epsilon \). This is the simplest model that exhibits both creep and stress relaxation. The second law of thermodynamics requires \({\hat{I}} \equiv {\hat{E}} - \tau E \ge 0\) or \({{\hat{E}} \ge \tau E}\) among the positive constant coefficients, inducing that the dynamic (high-frequency) modulus \(E_\infty = {\hat{E}} / \tau \) is never smaller than the static modulus E. Right: experimentally found relationship between dynamic and static bulk moduli for various types of rock (data taken from [5]). An analogous tendency is found for the shear moduli. Accordingly, many rocks can successfully be modeled by a PTZ model in the volumetric/spherical part and another PTZ model in the deviatoric part, in the spherical–deviatoric decomposition [6]

Full size image

Next, in [2], the scheme was extended to three spatial dimensions, for rectangular/cuboid geometries. The quantities provided by thermodynamics were shown beneficial in monitoring the solution:

  • preservation of total energy demonstrated that the energy-friendly (Hamiltonian-friendly) property of the underlying symplectic scheme of the reversible part of the system was preserved successfully during the generalization to the irreversible level;

  • the increase of temperature displayed where dissipation is the strongest; and

  • entropy production proved helpful in observing loss of stability early: its turning into negative indicated violation of the second law of thermodynamics and thus the emergence of instability.

Here, we present the inclusion of thermal expansion and heat conduction into the approach started in [1, 2]. The challenge is that temperature has to be raised to the level of a fully active field that has its own dynamics (heat conduction) and carries the coupling between this dynamics and the mechanical side via its role in the kinematic aspects through thermal expansion.

Can we still have an explicit finite-difference scheme, with second-order accuracy in both space and time, with low dispersion error and dissipation error? Additionally, can we further enhance the role of the thermodynamics-provided quantities in monitoring the quality of the solution?

In what follows, we find affirmative answers to these questions.

2 Continuum equations on the involved quantities

The thermodynamically consistent amalgamation of PTZ rheology, thermal expansion, and heat conduction is the following set of continuum equations (see detailed explanations below):

$$\begin{aligned} \varrho {\dot{v}}&= \frac{\partial \sigma }{\partial x} ,\end{aligned}$$
(1)
$$\begin{aligned} \sigma&= \sigma _{\textrm{el}}+ {\hat{\sigma }},\end{aligned}$$
(2)
$$\begin{aligned} L&= \frac{\partial v}{\partial x} ,\end{aligned}$$
(3)
$$\begin{aligned} {\hat{\sigma }}+ \tau \frac{\partial {\hat{\sigma }}}{\partial t}&= {\hat{I}} L ,\end{aligned}$$
(4)
$$\begin{aligned} e_{\textrm{tot}}&= \underbrace{\frac{v^2}{2}}_{e_{\textrm{kin}}} \mathrel + e_{\textrm{int}},\end{aligned}$$
(5)
$$\begin{aligned} e_{\textrm{int}}&= \underbrace{ c_{\textrm{ex}}\left( T - T_{\textrm{ex}}\right) }_{e_{\textrm{thermal}}} \mathrel + \underbrace{ \frac{E}{2\varrho } \,\! \epsilon \left( \epsilon + 2 \alpha T_{\textrm{ex}}\right) }_{e_{\textrm{el}}} \mathrel + \underbrace{ \frac{\tau }{2 \varrho {\hat{I}}} {\hat{\sigma }}^2}_{e_{\textrm{rheol}}}\end{aligned}$$
(6)
$$\begin{aligned} s&= c_{\textrm{ex}}\ln \frac{T}{T_{\textrm{ex}}} + \frac{E \alpha }{\varrho } \epsilon ,\end{aligned}$$
(7)
$$\begin{aligned} j_e&= - \lambda \frac{\partial T}{\partial x} ,\end{aligned}$$
(8)
$$\begin{aligned} j_s&= \frac{1}{T} j_e,\end{aligned}$$
(9)
$$\begin{aligned} \pi _s&= \frac{\partial \frac{1}{T}}{\partial x} j_e+ \frac{1}{{\hat{I}}} \frac{1}{T} \;\! {\hat{\sigma }}\left( {\hat{I}} L - \tau \frac{\partial {\hat{\sigma }}}{\partial t} \right) ,\end{aligned}$$
(10)
$$\begin{aligned} \varrho \frac{\partial e_{\textrm{int}}}{\partial t}&= - \frac{\partial j_e}{\partial x} + \sigma L ,\end{aligned}$$
(11)
$$\begin{aligned} \varrho \frac{\partial s}{\partial t}&= - \frac{\partial j_s}{\partial x} + \pi _s,\end{aligned}$$
(12)
$$\begin{aligned} \frac{\partial \epsilon }{\partial t}&= L ,\end{aligned}$$
(13)
$$\begin{aligned} \epsilon&= D + \alpha ( T - T_{\textrm{ex}}) ,\end{aligned}$$
(14)
$$\begin{aligned} \sigma _{\textrm{el}}&= E D . \end{aligned}$$
(15)

Here, we work in one spatial dimension (x), in the small-strain regime where, in leading order, density \(\varrho \) is constant and substantive derivatives can be replaced by partial time derivatives \(\frac{\partial }{\partial t}\) (see, e.g., [7]). The balance of momentum (1) determines the evolution of velocity v via the divergence of stress \({\sigma }\). Stress is the sum of an elastic (reversible) term \(\sigma _{\textrm{el}}\) and an irreversible one \({\hat{\sigma }}\), both explained below. From velocity, the velocity gradient L is defined, with the aid of which the Onsagerian time evolution of \({\hat{\sigma }}\) can be expressed as (4) ([6] applied to the PTZ case), in which \({\hat{I}}\) is positive (cf. the caption of Fig. 1).

Furthermore, prior to the balance (11) of internal energy \(e_{\textrm{int}}\) and the balance (12) of entropy s, the related constitutive ingredients need to be provided: first, specific total energy \(e_{\textrm{tot}}\) is the sum of specific kinetic energy \(e_{\textrm{kin}}\) and specific internal energy \(e_{\textrm{int}}\). The form (6) of \(e_{\textrm{int}}\) originates from [8], then it is simplified along the lines of [9] (page 7), namely, we restrict ourselves to a first-order expanded form in absolute temperature T around a chosen temperature value \(T_{\textrm{ex}}\) (expansion point) with constant specific heat capacity \(c_{\textrm{ex}}\).

Besides T, the second state variable of \(e_{\textrm{int}}\) is the geometric quantity \(\epsilon \), which equals strain when the process starts from unstressed and thermodynamically equilibrial state at \(T = T_{\textrm{ex}}\). In \(\epsilon \), we consider a second-order Taylor-expanded form, the quadratic term providing the customary Hookean elastic energy expression with Young’s modulus E, while the term linear in \(\epsilon \) gives rise to thermal expansion, with a constant thermal expansion coefficient \(\alpha \).

Finally, the third state variable is \({\hat{\sigma }}\) itself (in a viscoelastic model beyond the level of PTZ, a more general internal variable would stand here, see [6], Appendix B), and the constant coefficient in \(e_{\textrm{rheol}}\) is fixed by thermodynamical consistency with (4) and (11). Next, \(s( T, \epsilon )\) is also fixed by thermodynamical consistency (namely, the Gibbs relation with the reversible sum \(e_{\textrm{thermal}}+ e_{\textrm{el}}\)). Heat current density \(j_e\) obeys Fourier’s law (8) with heat conduction coefficient \(\lambda \), and entropy current density \(j_s\) is in the standard relationship (9) with \(j_e\). Entropy production rate density \(\pi _s\) consists of two terms, the first proving to be positive definite in light of (8) and the second because of (4), in accord with the second law of thermodynamics.

These have been the necessary ingredients for the balance of internal energy (11) and for that of entropy (12). Finally, the change rate of \(\epsilon \) is determined by L and, as done in [7] (cf. page 21), \(\epsilon \) is decomposed into an elasticity-related part denoted here by D and a thermal part \(\alpha ( T - T_{\textrm{ex}})\). D is the part that is related directly to elastic stress, through (15). We note that displacement u has not been used above but can be calculated from velocity as \(u_{t_0}^{} (t) = \int _{t_0}^{t} v \big ( {{\tilde{t}}} \big )\, \textrm{d}{{\tilde{t}}}\).

This is the set of equations that is the generalization of the one treated in [1, 4], when thermal expansion and heat conduction are added, both in their simplest form. Apparently, many quantities play a role here. In such a case, a frequent approach is to eliminate as many of the quantities as possible, reducing to a minimum the number of degrees of freedom to be followed. Below, we prefer to keep all these quantities as they provide flexibility in applying different boundary conditions required by different applications, and as they are shown below to prove useful for monitoring the quality of the solution. In case of such a coupled problem, any insight of this kind is highly valuable.

3 Discretization

The space and time staggered finite-difference scheme applied here (see Table 1) is also the generalization of the one used in [1, 4]. Accordingly, quantities reside either at integer multiples \(x_n = n \Delta x\) of the space step \(\Delta x\) or at half-integer ones \(x_{n-1/2} = (n-1/2) \Delta x\). Regarding time, not all quantities are classified as taken either at integer time steps \(t^j = j \Delta t\) or at half-integer ones \(t^{j-1/2} = (j - 1/2) \Delta t\): some are calculated both at integer time steps and at half-integer ones. In such cases, the ones displayed in grey in Table 1 are not predicted values but intermediate ones, similarly to the ones in Runge–Kutta-type schemes. They are introduced because Fourier heat conduction breaks the leapfrog pattern that was present in [1, 4] but we still aim at a second-order accurate scheme.

Table 1 The finite-difference scheme
Full size table

Regarding initial conditions, we assume an equilibrial state (each of v, \(\sigma _{\textrm{el}}\), \({\hat{\sigma }}\), L, \(\epsilon \), D, \({j_e}\), \(j_s\), \(\pi _s\) is zero) with homogeneous \(T = T_{\textrm{ex}}\) for \(t^0\) and \(t^{-1/2}\), and, for definiteness, the boundaries \(x_{-1/2}\), \(x_{N+1/2}\) are thermally adiabatic and a single velocity pulse

$$\begin{aligned} v ( t, x_{-1/2} ) = {\left\{ \begin{array}{ll} \frac{v_{\textrm{b}}}{2} \left[ 1 - \cos \left( 2 \pi \frac{t}{t_{\textrm{b}}} \right) \right] &{} \text {if} \quad 0 \le t \le t_{\textrm{b}},\\ 0 &{} \text {otherwise} \end{array}\right. } \end{aligned}$$
(16)

is applied at the left endpoint, where \(t_{\textrm{b}}\) and \(v_{\textrm{b}}\) denote the temporal width and the height of the pulse, respectively, while the right endpoint is fixed (so \(v = 0\)). Notably, these boundary conditions forbid heat-type energy change and, after \(t_{\textrm{b}}\), mechanical work is also forbidden at both boundaries so total energy should be conserved after the pulse, and even its increase during the pulse should be predicted as the time integral of mechanical power at the left endpoint.

The explicit scheme conceived is the following chain of steps. Above each equal sign, the corresponding continuum equation is displayed (if applicable: for certain “grey” intermediate quantities there is no such equation). As can be seen, the steps can be grouped into three blocks, the second and third performing the same steps with a shift of a half time step.

(17)
(18)
(19)

This scheme is second-order accurate both in space and time. This can be seen as follows. Each of the formulae in (17)–(19) is devised in such a way that it has a time instant and a space location around which each of the terms in that formula provides an accurate approximation of the corresponding continuum term not only at leading-order but also at next-to-leading order, when expanding into Taylor series in time and space both. For example, the third equation in the group (17) is centered around \(t^j\) and \(x_n\), and the expansion of the average on the right-hand side gives

(20)

Next, average is actually a special case of linear interpolation and a more general statement also holds: any linear interpolation or linear extrapolation is accurate not only at leading order but also at next-to-leading order: e.g., the interpolation/extrapolation of an f(t) from \(f (t - p\Delta t)\) and \(f(t - q\Delta t)\) (here, p, q may be negative) is

(21)

Such extrapolations are applied in the first equation in the group (19), extrapolating from and , and extrapolating from and (case \(p = 1/2\), \({q = 1}\)).

Further, the difference-ratio approximation of a space derivative is also accurate even at next-to-leading order around the midpoint: e.g.,

(22)

The analogous statement holds for time derivatives. Then, if a formula approximating a derivative with an error of \({\mathcal {O}} (\Delta t^2)\) is rewritten to give the time evolution of the field, the error of this approximation is \({\mathcal {O}} (\Delta t^3)\). For example,

$$\begin{aligned} \frac{\epsilon ^{j+1}_n - \epsilon ^j_n}{\Delta t} = L^{j+1/2}_n \quad \Longrightarrow \quad \epsilon ^{j+1}_n = \epsilon ^j_n + \Delta tL^{j+1/2}_n \end{aligned}$$
(23)

yields an approximation for the change of the \(\epsilon \) value that is accurate up to second order, i.e., up to \(\Delta t^2\).

This conclusion can also be demonstrated numerically, as we can see in the subsequent Section.

4 Simulation results

In the calculations, quantities have been made dimensionless (in notation: overtilde) with respect to the sample length X as length unit, the travel time along this length of the rheological wave with the maximal speed \({\hat{c}} = \sqrt{ \frac{ {\hat{E}}/\tau }{\varrho } }\) as the time unit, E (which fixes the mass unit), and \(c_{\textrm{ex}}\) (which fixes the temperature unit). Correspondingly, the dimensionless parameters and variables are

$$\begin{aligned} {\tilde{t}}= & {} \frac{{\hat{c}}}{X} t, \quad \Delta {\tilde{t}} = \frac{{\hat{c}}}{X} \Delta t, \quad {\tilde{\tau }} = \frac{{\hat{c}}}{X} \tau , \quad {\tilde{t}}_{\textrm{b}} = \frac{{\hat{c}}}{X} t_{\textrm{b}}, \end{aligned}$$
(24)
$$\begin{aligned} {\tilde{x}}= & {} \frac{1}{X} x, \quad \Delta {\tilde{x}} = \frac{1}{X} \Delta x, \quad {\tilde{X}} = \frac{X}{X} = 1, \quad \tilde{{\hat{c}}} = \frac{{\hat{c}}}{{\hat{c}}} = 1, \quad {\tilde{v}}_{\textrm{b}} = \frac{1}{{\hat{c}}} v_{\textrm{b}}, \quad \end{aligned}$$
(25)
$$\begin{aligned} {\tilde{\varrho }}= & {} \frac{{\hat{c}}^2}{E} \varrho , \quad {\tilde{E}} = \frac{E}{E} = 1, \quad \tilde{{\hat{E}}} = \frac{{\hat{c}}}{X E} {\hat{E}}, \quad \tilde{{\hat{I}}} = \frac{{\hat{c}}}{X E} {\hat{I}}, \quad \end{aligned}$$
(26)
$$\begin{aligned} {\tilde{\lambda }}= & {} \frac{{\hat{c}}}{E X c_{\textrm{ex}}} \lambda , \quad {\tilde{T}}_{\textrm{ex}} = \frac{c_{\textrm{ex}}}{{\hat{c}}^2} T_{\textrm{ex}}, \quad {\tilde{c}}_{\textrm{ex}} = \frac{c_{\textrm{ex}}}{c_{\textrm{ex}}} = 1, \quad {\tilde{\alpha }} = \frac{{\hat{c}}^2}{c_{\textrm{ex}}} \alpha . \end{aligned}$$
(27)

Further, the nondimensional counterparts of the various field variables are

$$\begin{aligned} {\tilde{D}}= & {} D, \quad {\tilde{\epsilon }} = \epsilon , \quad {\tilde{v}} = \frac{1}{{\hat{c}}} v, \quad {\tilde{L}} = \frac{X}{{\hat{c}}} L, \quad {\tilde{\sigma }}_{\textrm{el}} = \frac{1}{E} \sigma _{\textrm{el}}, \quad {\tilde{T}} = \frac{c_{\textrm{ex}}}{{\hat{c}}^2} T, \quad \end{aligned}$$
(28)
$$\begin{aligned} {\tilde{e}}_{\textrm{int}}= & {} \frac{1}{{\hat{c}}^2} e_{\textrm{int}}, \quad {\tilde{j_e}} = \frac{1}{E {\hat{c}}} j_e, \quad {\tilde{s}} = \frac{1}{c_{\textrm{ex}}} s, \quad {\tilde{j}}_{s} = \frac{{\hat{c}}}{E c_{\textrm{ex}}} j_s, \quad {\tilde{\pi }}_{s} = \frac{X {\hat{c}}}{E c_{\textrm{ex}}} \pi _s\end{aligned}$$
(29)

(and, for the other specific energy and stress quantities, accordingly). The dimensionless versions of (1)–(15) and (17)–(19) are the straightforward counterparts where a tilde is applied above all quantities.

During testing, we have scanned a large part of the space of nondimensional parameters and report here all the observed features we have found important to mention. Unless stated otherwise, the rheological Courant number \({\hat{C}} = {\hat{c}} \frac{\Delta t}{\Delta x}\) has been set to 1, \(N = 30\) space cells have been used, the simulated dimensionless time was 15, and \({\tilde{v}}_{\textrm{b}}=0.01\), \({\tilde{t}}_{\textrm{b}}=0.5\), \({\tilde{T}}_{\textrm{ex}}=1\) have been used. (Details on why \({\hat{C}} = 1\) appears the best choice follow later.) Five materials (Stripa granite [10,11,12], plastic PA6 [13, 14], and three artificial ones for probing certain extremities) have been considered, deliberately chosen to examine various circumstances induced by the rather different material properties. With the inclusion of the three artificial materials (Table 2), we can demonstrate all the remarkable effects we have found.

Table 2 Properties of the real materials used for the various simulations
Full size table
Fig. 2
figure 2

Numerical evidence demonstrating second-order accuracy of the numerical scheme. Displayed is the \(L^2\) error along the sample during one time step, shown for two quantities (presentation inspired by [15]). The error appears at the order of \(\Delta t^3\). Calculation done with material A (see Table 3)

Full size image
Table 3 Nondimensional material properties used for the various simulations
Full size table

First presented is the numerical demonstration of the second-order accuracy of the numerical scheme. As visible in Fig. 2, the simulation result agrees with the analytical one of the previous Section.

In the next group of results shown (see Fig. 3), artificial material B has been chosen, in order to make the effects of the coupled phenomena fully visible.

As can be seen, total energy is preserved reliably, and even its initial increase during the pulse agrees with the transferred energy (time integrated mechanical power at the boundary). All energies tend to their expected equilibrial asymptotic value. As a special case, elastic energy tends to a nonzero constant since dissipation raised temperature, which generated a nonzero \(\epsilon \) via thermal expansion. In parallel, entropy production rate is reassuringly positive, in consistency with the continuum expectation due to the second law of thermodynamics.

Fig. 3
figure 3

First row Space-time dependence of velocity (left) and of temperature (right) when both mechanical dissipation and heat conduction are strong. Other mechanical and kinematic quantities behave similarly to velocity. Temperature inhomogeneity created by dissipation (which is concentrated along the path of the traveling pulse) gets homogenized by heat conduction. Second row time dependence of the corresponding various energies (space integrals) and of entropy production rate (space integral of \(\pi _s\)). Results shown for material B (see Table 3)

Full size image

For comparison, these calculations have been repeated with heat conduction “switched off”; see Fig. 4 for the differences. The only obvious difference is that the distribution of temperature cannot get homogenized. The turning points create the most dominant amount of dissipation (viscoelasticity is most active where changes are the most considerable) and temperature increase remains localized there. In contrast, within the bulk, the many back-and-forth travels of the pulse generate dissipation approximately homogeneously.

Fig. 4
figure 4

Results analogous to those in Fig. 3 but with no heat conduction

Full size image

For the Stripa granite and for the PA6 plastic, heat conduction was found to produce a weak effect, somewhat similarly to the outcomes shown in Fig. 4. The other possible source for differences was the interplay among the viscoelastic creep and relaxation time scales, the travel time of the pulse along the sample, and the duration of the pulse. In many realistic cases, dissipation became visible only after many bounces. Naturally, the effect of viscoelasticity is visible even in such cases as the propagation speed of the pulse, \({\hat{c}} = \sqrt{ \frac{ {\hat{E}}/\tau }{\varrho } }\), is larger than the elastic/low-frequency value \(c = \sqrt{ \frac{E}{\varrho } }\).

Next follows the study of dispersion error. In Fig. 5, outcomes are shown for PA6 plastic, with the rheological Courant number lowered from 1 to 0.5, at room temperature (\({\tilde{T}}_{\textrm{ex}}=0.39\)). No problem is visible in the total energy, however, the various energy terms unveil the presence of dispersion error. Moreover, the dubious change in the maxima/minima of the energy terms raises concern regarding dissipation error.

Fig. 5
figure 5

With Courant number \({\hat{C}} = 0.5\), dispersion error is observable. Note the gradually emerging ghost images around the pulse in the spacetime diagrams, and the artificial oscillations in internal and kinetic energies. Result shown for material PA6 (see Table 3)

Full size image

This and similar results visually indicate that, for this solid mechanical model, rheological Courant number 1 is a superior choice compared to other choices. Actually, one can also have a more quantitative support for this conviction. Namely, a starting point is the previous works [1, 2] (studying the case of no heat conduction and no thermal expansion) where it is not only observed but also verified, via analyzing the dispersion relation, that Courant number 1 is the optimal value regarding dispersion error. Due to the parabolic nature of heat conduction, a plausible expectation is that the same may apply here as well. Here, analytical studies would be much more involved because of the more complicated model so we include a numerical study where the Courant number is scanned between 0.3 and 1 while all other parameters (including the space step) are kept fixed. (The time step is changed according to \(\Delta t= \frac{{\hat{C}}}{{\hat{c}}} \Delta x\).)

Fig. 6
figure 6

Left Courant number dependence of the shape of the propagating excitation. As \({\hat{C}}\) is gradually decreased from 1, the artificial oscillation behind the pulse gets larger and larger. Right The \(L^2\) deviation of the excitation shape from that of the \({\hat{C}} = 1\) one. Calculation performed with material PA6

Full size image

Now, looking back to Fig. 6, we can see that, with optimal numerical settings, the scheme is energy-friendly not only in the sense of total energy conservation but the various energy terms themselves are also reliably calculated; while for suboptimal numerical settings the various energy terms display the appearance of numerical artefacts.

Finally, instability is investigated. In our coupled problem, it is not only the Courant number that poses a constraint on the time and space steps but heat conduction also brings in its stability criterion.

Fig. 7
figure 7

Loss of stability, and how discretized entropy production rate becoming negative reveals the problem early, well before energies become problematic. Calculation with material C (see Table 3)

Full size image

It may be difficult to cope with both requirements. Figure 7 shows such a case where we drove the simulation somewhat above the heat conduction related stability criterion, which was found to be \(\frac{\lambda }{\varrho c_{\textrm{ex}}} \frac{\Delta t}{\Delta x^2} \le 0.333\) via experimenting. As a generalization of the idea introduced in [2], we discretize \(\pi _s\) in a deliberately non-positive definite form—although, at the continuum level, the equations would ensure its positive definiteness—:

(30)

This way we make it possible that the discretized \(\pi _s\) takes negative value. Since the second law of thermodynamics is related to Lyapunov asymptotic stability, its violation at the discrete level is expected to reveal numerical loss of stability. With an adequately chosen discretization, entropy production rate may signal instability much before the problem manifests itself in other quantities. In the example visible in Fig. 7, instability becomes apparent at dimensionless time \({{\tilde{t}}} \approx 2\) while the space integral \({\tilde{\Pi }}_s\) of \({\tilde{\pi }}_s\) becomes negative already at \({{\tilde{t}}} \approx 1\). The idea introduced in [2] works in the present generalized context, too: entropy production rate proves to be a good instability indicator.

5 Discussion

The proposed explicit finite-difference scheme has second-order accuracy in both space and time, has maintained the good energy-preserving property and large-time reliability of the symplectic scheme of the reversible core of the system, and can avoid dispersion and dissipation error. Regarding monitoring the solution, time dependence of the various energy terms proves to make dispersion error visible, and entropy production rate can be utilized to point out instability at an early stage. These make the scheme applicable for simulating viscoelastic wave propagation coupled to heat conduction via thermal expansion, where dissipation resulting from viscoelasticity and heat conduction could be miscalculated due to dissipation error, wave-based oscillations could be miscalculated due to dispersion error, and numerical settings might lead to unstable solution.

Despite the somewhat complex structure of the scheme (17)–(19), it may readily be generalized to three spatial dimensions, analogously to how [2] generalized the scheme set out in [1]. This would enable to treat rectangular/cuboid geometries. Next, along the lines of [16], extension to cylindrical setups also appears realizable. Then, based on the experience obtained here, one becomes able to simulate wave propagation in cylindrical rock samples for applications like the dilatational resonant frequency method, to evaluate measurements of the Anelastic Strain Recovery method with a higher precision, and to determine gravitational consequences of waves around underground tunnels of gravitational wave detectors. Other areas include simulations for plastics, biological samples, and other solid materials where viscoelasticity is remarkable and thermal effects also play a considerable role. Some of these (e.g., [17]) would require the large-strain generalization of the current framework, which may be carried out based on the spacetime-oriented approach [18]. The scheme communicated here is expected to be generalizable to incorporate beyond-Fourier heat conduction (see, e.g., the recent review [19]) as well, where calculation of the additional irreversibility-related degrees of freedom may be successfully addressed via the scheme-building strategy presented here – steps in this direction have been reported in [20, 21].