Thermodynamically Extended Symplectic Numerical Scheme with Half Space and Time Shift Applied for Rheological Waves in Solids

On the example of the Poynting–Thomson–Zener rheological model for solids, which exhibits both dissipation and wave propagation – with nonlinear dispersion relation –, we introduce and investigate a finite difference numerical scheme for continuum thermodynamical problems. The key element is the positioning of the discretized quantities with shifts by half space and time steps with respect to each other. The arrangement is chosen according to the spacetime properties of the quantities and of the equations governing them. Numerical stability, dissipative error and dispersive error are analysed in detail. With the best settings found, the scheme is capable of making precise and fast predictions.


Introduction
Numerical solution methods for dissipative problems are important and are a nontrivial topic. Already for reversible systems, the difference between a symplectic and nonsymplectic finite difference method is striking: the former can offer reliable prediction that stays near the exact solution even at extremely large time scales while the latter may provide a solution that drifts away from the exact solution steadily. For dissipative systems, the situation is harder. Methods that were born with reversibility in mind may apparently fail for a nonreversible problem. For example, a finite element software is able to provide, at the expense of large run-time, quantitatively and even qualitatively wrong outcome while a simple finite difference scheme solves the same problem fast and precisely [1].
Thermodynamics also modifies the way of thinking concerning numerical modelling. Even if quantities known from mechanics form a closed system of equations to solve numerically, monitoring temperature (or other thermodynamical quantities) for a nonreversible system can give insight on the processes and phenomena, for example, pointing out the presence of viscoelasticity/rheology, and displaying when plastic changes start [10]. In addition, temperature can also react, in the form of thermal expansion and heat conduction, even in situations where one is not prepared for this 'surprise' [11].
Furthermore, in a sense, thermodynamics is a stability theory. Therefore, how thermodynamics ensures asymptotic stability for systems may give new ideas on how stability and suppression of errors arXiv:1908.07975v1 [physics.class-ph] 21 Aug 2019 can be achieved for numerical methods. A conceptually closer relationship is desirable between these two areas.
Along these lines, here, we present a study where a new numerical scheme is suggested and applied for a continuum thermodynamical model. The scheme proves to be an extension of a symplectic method. In parallel, our finite difference scheme introduces a shifted arrangement of quantities by half space and time steps with respect to each other, according to the spacetime nature of the involved quantities and the nature of equations governing them, since balances, kinematic equations, and Onsagerian equations all have their own distinguished discretized realization. This also makes the scheme one order higher precise as the original symplectic scheme.
The continuum system we take as the subject of our investigation is important on its own -it is the Poynting-Thomson-Zener rheological model for solids. This model exhibits dissipation and wave propagation (actually: dispersive wave propagation) both, and is thus ideal for testing various aspects and difficulties. Meanwhile, its predictions are relevant for many solids, typically ones with complicated microor mesoscopic structure like rocks [12][13][14], plastics [10], asphalt etc. This non-Newtonian rheological model can explain why slow and fast measurements and processes give different results.
Solutions in the force equilibrial and space independent limit have proved successful in explaining experimental results [10]. Space dependent -but still force equilibrial -analytical solutions can model opening of a tunnel, gradual loading of thick-walled tubes and spherical tanks, and other problems [15]. The next level is to leave the force equilibrial approximation, partly in order to cover and extend the force equilibrial results but also to be utilized for evaluating measurements that include wave propagation as well. The present work is, in this sense, the next step in this direction.

Properties of the continuum model
The system we consider is a homogeneous solid with Poynting-Thomson-Zener rheology, in small-strain approximation 1 , in one space dimension (1D). Notably, the numerical scheme we introduce in the following section can be generalized to 2D and 3D with no difficulty 2 -the present 1D treatment is to keep the technical details at a minimum so we can focus on the key ideas.
The set of equations we discuss is, accordingly, ∂v ∂t where is mass density (constant in the small-strain approximation), (1) tells how the spatial derivative of stress σ determines the time derivative of the velocity field v (volumetric force density being omitted for simplicity), (2) is the kinematic relationship between the strain field ε 3 and v, and the rheological relationship (3) contains, in addition to Young's modulus E, two positive coefficientsÊ, τ.
1 Hence, there is no need to distinguish between Lagrangian and Eulerian variables, and between material manifold vectors/covectors/tensors/. . . and spatial spacetime ones. 2 The results of our ongoing research on 2D and 3D are to be communicated later. 3 In the present context, ε can be used as the thermodynamical state variable for elasticity, but not in general, see [17,18].
The Poynting-Thomson-Zener model is a subfamily within the Kluitenberg-Verhás model family, which family can be obtained via a nonequilibrium thermodynamical internal variable approach [16]. The Poynting-Thomson-Zener model looks particularly simple, after eliminating the internal variable, both in specific energy e total and in specific entropy s: On the other side, for processes much faster then the two time scales, keeping the highest time derivatives leads to that is, for stress and strain changes (e.g., for deviations from initial values), the system effectively behaves like a Hooke one, with 'dynamic' Young's modulus The corresponding effective wave equation possesses the wave speed For a more rigorous and closer investigation of these aspects, the dispersion relation can be derived. Namely, on the line −∞ < x < ∞, any (not too pathological) field can be given as a continuous linear combination of e ikx space dependences, where the 'wave number' k is any real parameter. If such a (Fourier) decomposition is done at, say, t = 0, then the subsequent time dependence of one such mode may be particularly simple: (15) with some appropriate ω -complex, in general -; the factor i in the first component is introduced in order to be in tune with later convenience. A space and time dependence e −iωt e ikx = e Im ωt e −i Re ωt e ikx = e −(− Im ω)t e ik(x− Re ω k t) expresses travelling with constant velocity Re ω k and exponential decrease (for dissipative systems like ours, Im ω < 0). In general, it depends on the number of fields and on the order of time derivatives how many ω's are possible. In our case, the relationship between compatible ω and k -the dispersion relationis straightforward to derive: In the limit |ω| → 0 (limit of slow processes), we find while in the opposite limit |ω| → ∞ (limit of fast processes), the result is Both results confirm the findings above [(11) and (14), respectively]. This is a point where we can see the importance of the Poynting-Thomson-Zener model. Namely, when measuring Young's modulus (or, in 3D, the two elasticity coefficients) of a solid, the speed of uniaxial loading, or the frequency of sound in an acoustic measurement, may influence the outcome and an adequate/sufficient interpretation may come in terms of a Poynting-Thomson-Zener model. Indeed, in rock mechanics, dynamic elastic moduli are long known to be larger than their static counterparts (a new and comprehensive study on this, [19], is in preparation), in accord with the thermodynamics-originated inequality in (14) (or its 3D version).

The numerical scheme
The classic attitude to finite difference schemes is that all quantities are registered at the same discrete positions and at the same discrete instants. An argument against this practice is that, when dividing a sample into finite pieces, some physical quantities have a meaning related to the bulk, the centre of a piece, while others have a physical role related to the boundaries of a unit. For example, (specific) extensive and density quantities would naturally live at a centre, while currents/fluxes are boundary related by their physical nature/role. When one has a full -at the general level, 4D -spacetime perspective 6 then it turns out that quantities may "wish" to be shifted with respect to each other by a half in time as well. This latter aspect is less straightforward to visualize but the structure of the equations -for example, the structure of balanceshelps us to reveal what intends to be shifted with respect to what. This is what we realize for the present system. Discrete space and time values are chosen as and discrete values of stress are prescribed to these spatial and temporal coordinates: Then, investigating (1), we decide to put velocity values half-shifted with respect to stress values both in space and time: 7 v j n at t j − ∆t 2 , x n + ∆x 2 , and discretize (1) as Next, studying (2) suggests analogously to have strain values half-shifted with respect to velocity values both in time and space. Therefore, strain is to reside at the same spacetime location as stress: 6 Traditional physical quantities are usually time-and spacelike components of four-vectors, four-covectors, four-cotensors etc. and (2) is discretized as Finally, for the Hooke model, (10) is discretized plainly as we can recognize the steps of the symplectic Euler method [20] (with the Hamiltonian corresponding to e kinetic + e elastic ). Now, a symplectic method is highly favourable because of its extremely good large-time behaviour, including preservation of energy/etc. conservation. While (27) coincides with the symplectic Euler method computationally, the present interpretation of the quantities is different, because of the space and time shifts. One advantageous consequence is that, due to the reflection symmetries (see Figure 1), our scheme makes second order precise predictions (understood in powers of ∆t and ∆x), while the symplectic Euler method makes only first order precise ones [20]. Indeed, setting and assuming that precisely, the error of the prediction for v after Taylor series expansion, cancellations, and the use of (1).
Analogously, with second order preciseness of prediction ε j+1 n can be proved. In case of the Poynting-Thomson-Zener model, we need to discretize (3). Here, both σ and its derivative, and both ε and its derivative, appear. Hence, shifting does not directly help us. This is what one can expect for dissipative, irreversible, relaxation type, equations in general. However, an interpolation-like solution is possible: where α = 1/2 is expected to provide second order precise prediction, and other seminal values are α = 1 (the explicite case, which is expected to be stiff) and α = 0 (the fully implicite case). For generic α, (32) looks implicit. However, actually, thermodynamics has brought in an ordinary differential equation type extension to the Hooke continuum, not a partial differential equation type one, and a linear one, in fact. Thus (32) can be rewritten in explicit form, assuming Second order preciseness of (33) for α = 1/2 is then straightforward to verify, in complete analogy to the two previous proofs.

Stability
One may specify a space step ∆x according to the given need, adjusted to the desirable spatial resolution. In parallel, the time step ∆t is reasonably chosen as considerably smaller than the involved time scales (e.g., τ andτ of our example system). Now, a finite difference scheme may prove to be unstable for the taken ∆x and ∆t, making numerical errors (which are generated unavoidably because of floating-point round-off) increase essentially exponentially and ruining usefulness of what we have done. Therefore, first, a stability analysis is recommended, to explore the region of good pairs of ∆x, ∆t for the given scheme and system.
We continue with this step for our scheme and system, performing a von Neumann investigation [21], where the idea is similar to the derivation of the dispersion relation. There, we study time evolution of continuum Fourier modes [see (15)], while here examine whether errors, expanded in modes with e ikx n space dependence, increase or not, during an iteration by one time step. For such linear situations as ours -when the iteration step means a multiplication by a matrix -, such a mode may simply get a growth factor ξ (that is k dependent but space independent), in other words, the iteration matrix (frequently called 'transfer matrix') has these modes as eigenvectors with the corresponding eigenvalues ξ. Then, |ξ| < 1 (for all k) ensures stability. Furthermore, |ξ| = 1 means stability if the algebraic multiplicity of ξ -its multiplicity as a root of the characteristic polynomial of the transfer matrix -equals its geometric multiplicity -the number of linearly independent eigenvectors (eigenmodes), i.e., the dimension of the eigensubspace ( [22], page 186, Theorem 4.13; [23], page 381, Proposition 2).
With boundary conditions specified, one can say more. 8 Boundary conditions may allow only certain combinations of e ikx n as eigenmodes of the transfer matrix. Consequently, this type of analysis is more involved and is, therefore, usually omitted. As a general rule-of-thumb, one can expect that |ξ| > 1 for some e ikx n indicates instability also for modes obeying the boundary conditions, while |ξ| ≤ 1 for all e ikx n suggests stability for all modes allowed by the boundary conditions 9 .

Hooke case
In the Hooke case, the 'plane wave modes' for the two bookkept quantities v, ε can, for later convenience, be written as the condition on k related to that k outside such a 'Brillouin zone' makes the description redundant.
Realizing the iteration steps (27) as matrix products leads, for the amplitudes introduced in (35), to For space dependences (35), 8 All systems require boundary or asymptotic conditions. We also specify some in the forthcoming section on applications. 9 Namely, the problem of differing multiplicities for |ξ| = 1 can be wiped out by the boundary conditions. in other words, to the eigenvalue problem Let us introduce the notation for the Courant number of our scheme for the Hooke system. Comparing the characteristic polynomial of T, with its form written via its roots, reveals, on one side, that, in order to have both |ξ + | ≤ 1 and |ξ + | ≤ 1, both magnitudes have to be 1 (since their product is 1), which, on the other side, also implies as both C and S are non-negative. If CS < 1 then the two roots, are complex, with unit modulus, and are the complex conjugate of one another. Especially simple -and principally distinguished, as we see in the next sections -is the case C = 1: then with the remarkable property that arg ξ ± are linearly depending on k -so to say, both branches of the discrete dispersion relation are linear.
In parallel, if CS = 1 then the two roots coincide, ξ ± = −1. The algebraic multiplicity 2 is accompanied with geometric multiplicity 1: only the multiples of are eigenvectors. If C = 1 then this affects only one mode, S = 1, k = π k , and if that mode is prohibited by the boundary conditions then the choice C = 1 ensures a stable scheme.
With C > 1, CS ≤ 1 would be violated by a whole interval for k [recall (37)], which may not be cured by boundary conditions so the best candidate (largest ∆t for a fixed ∆x, or the smallest possible ∆x for fixed ∆t) to have stability is C = 1.

Poynting-Thomson-Zener case
For the Poynting-Thomson-Zener system, the von Neumann stability analysis of our discretization studies the modes gives The characteristic polynomial is noŵ Three roots are considerably more difficult to directly analyse. One alternative is to use Jury's criteria [24] for whether the roots are within the unit circle of the complex plane, and another possibility is to apply the Möbius transformation ξ = η+1 η−1 on (52) and utilize the Routh-Hurwitz criteria whether the mapped roots are within the left half plane. The two approaches provide the same result. Nevertheless, one criterion provided by one of these two methods may not directly be one criterion of the other method. It is only the combined result (the intersection of the conditions) that agrees. Accordingly, it can be beneficial to perform both investigations because a simple condition provided by one of the routes may be labouring to recognize as consequence of the conditions directly offered by the other route.
Jury's criteria, for our case, are as follows. First,P(1) > 0 gives Second, (−1) 3P (−1) > 0 yields which, in light of (54), reduces to Third, the matrices a 3 a 2 0 a 3 ± 0 a 0 a 0 a 1 have to be positive innerwise 10 . The '+' branch leads to which is weaker than (56), because there the rhs is larger by (1 − α) + τ ∆t C 2 S 2 [and cf. (54)]. Meanwhile, the '−' branch induces conditionτ which we have already met in (9), as the thermodynamical requirement (6) at the continuum level, and which also induces, via (57), which is stronger than (54). This also allows to rearrange (56) and exploit it as Conditions (58)-(60) summarize the obtained stability requirements; the first referring to the constants of the continuum model only, the second relating α and ∆t of the scheme, and the third limiting ∆x (through C) in light of α and ∆t.
If, instead of Jury's criteria, one follows the Routh-Hurwitz path, on the Möbius transformed polynomialQ then, having b 3 > 0, roots lie in the left half plane if all corner subdeterminants of As expected, these conditions prove to be equivalent to the ones obtained via Jury's criteria -we omit the details to avoid redundant repetition.

Kelvin-Voigt model
Although the focus of the present paper is on the hyperbolic-like case corresponding to τ > 0, the above calculations are valid for τ = 0, the Kelvin-Voigt subfamily as well. As a brief analysis of this case, (58) is trivially satisfied withτ > 0. (59) gives the nontrivial condition α < 1/2. 11 Finally, (60) gives which looks like some mixture of a stability condition for a scheme for a parabolic problem like Fourier heat conduction and of a condition for a simple reversible wave propagation.

Beyond Kelvin-Voigt
When τ > 0 thenĈ [recall (14)] becomes important. The most interesting case is α = 1/2, where the scheme gives second order precise predictions: (59) holds trivially, and (60) can be rewritten asĈ With boundary conditions also present, we may extend this condition tô Considering the two other potentially interesting cases as well: If α = 1 then (59) induces ∆t < 2τ, which is not a harsh requirement since the time step must usually be much smaller then the time scales of the system in order to obtain a physically acceptable numerical solution. In parallel,Ĉ is limited from above by a number smaller than 1. On the other side, when α = 0 then (59) is automatically true again, and nowĈ is limited from above by a number larger than 1. Since we may need ∆t τ for a satisfactory solution, this O ∆t τ gain over 1 is not considerable.

Hooke case
It is worth looking back to the Hooke limit of (60): τ =τ = 0 (with whatever α) tells C < 1. One can see that the |ξ| < 1 stability requirement gives conservative results and does not tell us how far the obtained inequalities are from equalities.

Numerical results
The calculations communicated here are carried out with zero v, ε, σ as initial conditions, and with stress boundary conditions: on one end of the sample, a cosine shaped pulse is applied, while the other end is free (stress is zero). With τ b denoting the temporal width of the pulse, the excitation is, accordingly, For making our quantities dimensionless -suitable for numerical calculations -the following units are used: the length of the sample, c (so a Hookean wave arrives at the other end during unit time), σ b , and, for temperature, c σ . Henceforth, with respect to this time unit, 0.2 is used for τ b , 1.25 for τ, and 5 forτ, implyingĉ = 2c ≡ 2. Temperature is calculated according to the discretized form of (8), with the natural choice that temperature values reside at the same place as stress and strain but half-shifted in time (T j n at t j − ∆t 2 , x n ). When plotting, say, elastic energy of the whole sample at time t j , a simple E 2 ∑ n ε , are counted with weight 1 2 . Second, kinetic energy and thermal energy, both being based on quantities half-shifted in time, are calculated as a time average, their value at t j taken as the average of their value at t j − ∆t 2 and t j + ∆t 2 .

Hookean wave propagation
For the Hooke system, our scheme is symplectic, with very reliable long-time behaviour. This is well visible in Figure 3: the shape is nicely preserved, no numerical artefacts are visible in the spacetime picture, and the sum of elastic and kinetic energy is conserved.

Poynting-Thomson-Zener wave propagation
For the Poynting-Thomson-Zener system, we find that the principally optimal choice of α = 1/2 does outperform α = 0 (withĈ = 1). Figure 4 shows such a comparison: α = 1/2 produces a reliable signal shape quite independently of space resolution, while α = 0 needs more than N = 1000 space cells to reach the same reliability. Actually, α = 1/2 offers that realibility already at N = 50, and even N = 25 'does a decent job', as depicted in Figure 5.
With α = 1/2, the spacetime picture and total energy conservation are not less satisfactory, as visible in Figure 6.  The physical explanation of the signal shape (Figures 4-5) is that the fastest modes propagate with speedĉ (recall Section 2), transporting the front of the signal, while slow modes travel with c <ĉ, gradually falling behind, and forming a little-by-little thickening tail.
In parallel, the spacetime picture shows that this tail effect is less relevant than the overall decrease of the signal, due to dissipation.
Finally, concerning the energy results, the remarkable fact is that all ingredients v, ε, σ, T are calculated via discretized time integration, therefore, total energy conservation is not built-in but is a test of the quality of the whole numerical approach.

Hooke case
The Hooke system might appear as a simple introductory task for numerics. This is actually far from true. Already the Hooke case displays both dissipation error and dispersion error if not treated with appropriate care [25]. While the greatest danger, instability, is about exponential exploding of error, dissipation error is 'the opposite': when the signal decreases in time, losing energy due to numerical artefact only. This type of error is related to |ξ| < 1 modes, which indicates that one should try to stay on the unit circle with ξ. On the other side, in addition to the modulus of ξ, its argument can also cause trouble: if arg ξ is not linear in k then dispersion error is induced, which is observable as unphysical waves generated numerically around signal fronts. These errors are present even in a symplectic scheme like ours, as illustrated in Figure 7. More insight is provided by Figure 8.

Poynting-Thomson-Zener case
In case of a dissipative system like the Poynting-Thomson-Zener one, it is hard to detect the dissipative error, i.e., to distinguish it from the physical dissipation. The dispersion error remains visible, as Figure 9 shows. Usually, one would need to set ∆t to be much smaller than τ,τ (and τ b ) to obtain a physically acceptable approximation. Rewriting the coefficients (53) as in the limit ∆t τ → 0, the characteristic polynomial reduces to with roots satisfying ξ 0 = |ξ + | = |ξ − | = 1, excluding thus dissipation error. Especially simple and distinguished is the caseĈ = 1, when the roots are providing dispersion relations linear in k and, hence, getting rid of dispersion error as well.
With slightly nonzero ∆t τ , these nice properties get detuned but only up to O ∆t τ , as shown in Figures 10-12 (prepared at dimensionless time step value 0.01; the detuning appears weaker for α = 1/2 than for α = 0).

Discussion
Choosing a good finite difference numerical scheme for a continuum thermodynamical problem is not easy. A good starting point can be a symplectic scheme for the reversible part, as done here, too. Another advantage is provided by a shifted arrangement of quantities by half space and time steps, suited to balances, to the kinematic equations, to the Onsagerian equations etc.
Even with all such preparations, instability is a key property to ensure. And when all these are settled, dissipative and dispersive errors can invalidate our calculation, which may not be recognized when the continuum system is dissipative and when it allows wave propagation as well.
In the future, the study provided here can be supplemented by comparison with analytical and finite element solutions.
Another logical continuation of this line of research is extension of the scheme to 2D and 3D spacethis is actually work in progress [25].  Figures 10-11, the roots are not exactly on the unit circle -here, ∆t dependence of |ξ 0 | and |ξ ± | is displayed, at a neutral value k∆x = π/4, forĈ = 1 and α = 1/2. Concerning the thermodynamical system to be investigated, the whole Kluitenberg-Verhás family -which the present Poynting-Thomson-Zener model is a subclass of -could be studied. The presence of second derivative of strain, and actually already the Kelvin-Voigt subfamily, brings in the aspect of parabolic characteristics so useful implications may be gained for other thermodynamical areas like non-Fourier heat conduction as well.
In conclusion, reliable numerical methods for thermodynamical systems, which avoid all the various pitfalls, are an important direction for future research.