1. Introduction
Solids may be less “solid” than expected. Beyond elastic behaviour, they may exhibit damped and delayed response. This viscoelastic/rheological reaction may not be simply explained by a viscosity-related additional stress (the Kelvin–Voigt model of rheology), but the time derivative of stress may also be needed in the description, with the simplest such model being the so-called standard or Poynting–Thomson–Zener (PTZ) one [see its details below]. Namely, the PTZ model is the simplest model that enables describing both creep (declining increase of strain during constant stress) and relaxation (declining decrease of stress during constant strain), as well as the simplest one, via which it is possible to interpret that the dynamic elasticity coefficients of rocks are different from, and larger than, their static counterpart [
1,
2,
3,
4]. Related to the latter aspect, high-frequency waves have a larger propagation speed in PTZ media than low-frequency ones [
4], which makes this model relevant for, e.g., seismic phenomena and acoustic rock mechanical measurement methods.
Analytical solutions for problems in PTZ and more complex rheological solid media exist (see, e.g., [
5,
6]), but mostly in the force-equilibrial/quasistatic approximation, which cannot give account of transients and waves. Incorporating such ‘fast’ effects is expected to only be realizable by means of numerical calculations in most practical situations.
In many of the practical applications, such a numerical calculation must be performed many times with different material coefficients, for example, as part of a fitting procedure where the experimental data are to be fitted. Hence, the numerical scheme should be fast, resource-friendly, yet reliable and precise enough.
In addition, numerical calculations face three frequent challenges: instability (exponential blow-up of the solution), dissipation error (artificial decrease of amplitudes and energies), and dispersion error (artificial oscillations near fast changes). A good scheme keeps these artefacts under control.
Being driven by (primarily rock mechanical) applications in scope, we have tried to use commercial finite element softwares for wave phenomena in PTZ models. What we found—already for the Hookean case (but also for non-Fourier heat conduction [
7])—was disappointing: the solutions ran very slowly, with large memory and CPU demand, and they were burdened by considerable numerical artefacts of the mentioned kinds.
Now, if a numerical scheme exhibits dissipation error for conservative systems, then it is expected to behave similarly for nonconservative ones, so that one cannot separate the real dissipation of mechanical energy from the dissipation artefact of the scheme. Additionally, similarly, a real wavy behaviour cannot be distinguished from the dispersion/wavy artefact.
These have motivated us to develop our own numerical scheme, which performs better [
8]. Similarly to that the PTZ model can be obtained in a thermodynamical approach as an internal variable extension of Hooke elasticity [
9], our starting point was a symplectic scheme for Hooke elasticity. Symplectic numerical schemes (see, e.g., [
10]) provide much better large-time approximations, thanks to the fact that a symplectic numerical integrator of a conservative system is actually the exact integrator of a ‘nearby’ (coinciding in the zeroth order of the time step) conservative system.
In recent years, numerous works have been born in order to develop extensions of symplectic schemes to nonconservative systems [
11,
12,
13,
14,
15,
16,
17,
18,
19,
20]. We also took this path, and devised such an extension, on the example of the PTZ model, in one space dimension [
8], with the novelties that some of the discretized field values reside with half space and time steps with respect to some other field values. This made the symplectic Euler method—an originally order-one accurate scheme—accurate to a second order, and the spatial accuracy was also second order.
Indeed, our scheme has performed well: produced, in a much faster and resource-friendlier way, much more artefact-free solutions, as demonstrated in
Figure 1. We note that, although the PTZ model allows for an exact integrator in the nonconservative part of the model, along analogous lines, as in [
12], we have refrained from using it, since, in the future, we also wish to use the same scheme for more general nonconservative systems, so we intended to test robustness in the dissipative aspect.
In addition to comparison to finite-element solutions, in [
8], we analytically derived the criteria for stability, and showed how dissipation and dispersion error can be kept small.
The first study done while using our scheme was numerically measuring wave propagation speed in one space dimensional samples, and finding good agreement with the corresponding analytical result [
4].
The next step is reported here, with two novelties. The first is the generalization of the scheme to three space dimensions (3D), and the other is exploiting the whole thermodynamical theory around the PTZ model for diagnostics regarding the credibility of the numerical solution, since stability, dissipation, and dispersion error are much harder to investigate in the case of a 3D model, with its numerous vectorial and tensorial degrees of freedom.
The 3D scheme is designed in order to keep the nice, second-order, behaviour of the discretization in both the spatial and in temporal direction. Achieving this is not so trivial—different components of vectors and tensors are placed at different discretized positions in order to fulfil the aim. In parallel, the boundaries also pose a challenge: quantities must be placed in such a way that the set of equations becomes closed. We succeed in finding a rule for this that is general enough to hold for both stress and displacment boundary conditions, where these two may differ at different sides of the 3D sample.
In finding the arrangement of discretized quantities suggested here, the spacetime perspective has helped us a lot. Specifically, on one side, thermodynamical balances in their differential form are four-divergences from the spacetime aspect, and they have an integral counterpart, which, via Gauss’ theorem, helps one to find where to represent which flux-type quantity. Additionally, knowing that, from the spacetime point of view, velocity is a timelike four-vector, [
21,
22] gives the information that velocity should be shifted not only spatially, but also temporally. Oppositely, stress is a spacelike tensor, so no temporal shift is needed.
In parallel, thermodynamics is not only important from the aspects of balances. Namely, commercial finite element softwares only focus on the set of equations to solve, i.e., on the minimally necessary equations to follow the minimally necessary quantities. However, knowing, from the spacetime perspective, that momentum and energy form a four-quantity (also in continuum theory on Galilean spacetime) [
22], in addition to the customarily taken balance of momentum, the balance of energy is also present. This enables one to follow, in addition to the mechanically considered quantities, internal energy—or, if practice favours so, temperature. The point in doing so (even in situations where thermomechanical coupling and heat conduction are neglected) is that, if, say, temperature is followed via a separate discretized evolution equation, then the conservation of total energy—at the discretized level—is not built-in, but rather is a property that will hold only approximately. Subsequently, checking how well this conservation holds along the numerical solution can provide a diagnostic tool. Thus, one may check the degree of dissipation error (i.e., the degree of violation of total energy conservation) and of dispersion error (spurious oscillations on total energy that should be a constant). This idea is demonstrated below, on the example of the PTZ model (also equipped with the thermodynamical constituents).
Furthermore, thermodynamics also provides entropy, which is known to serve as a Lyapunov function, ensuring asymptotic stability (see, e.g., [
23]). Now, stability also becomes challenged at the numerical level. Accordingly, entropy, and the related entropy production, may serve as an aid for reliable numerical calculations. Certain efforts in this direction have already been made [
11,
12]. Here, we introduce another way of utilizing this general idea.
Specifically, we focus on entropy production. At the continuum level, it must be positive definite according to the second law of thermodynamics. However, when discretized, this property may also become challenged. Naturally, if an explicitly positive definite expression is discretized, then it remains positive definite. However, alternative forms—which only turn out to be positive definite when the further thermodynamical equations also hold—are not
ab ovo positive definite and, correspondingly, may fail in being/remaining so along a numerical solution. Such forms are provided in a natural way, for instance, when the balance of entropy is connected to the balance of internal energy, such as when rheological models, like the PTZ one, are derived from the internal variable approach [
9]. Here, we discretize such an expression of entropy production and show that its value becoming negative can forecast a loss of stability and blowing-up of the solution.
2. The Continuum PTZ Model and the Thermodynamics Behind
We consider a homogeneous and isotropic solid, in the small-deformation approximation (with respect to an inertial reference system), due to which we do not have to differentiate between Eulerian and Lagrangian position or make a distinction between spatial spacetime vectors, covectors, tensors, etc, as well as material manifold related ones, mass density
can also be treated as constant, and the relationship between the symmetric strain tensor
to the velocity field
v is
with
and
denoting the spatial derivative operation acting to the right and left, respectively.
The stress tensor
is also assumed to be symmetric, and it governs the time evolution of
v according to
With the deviatoric and spherical parts of tensors,
(1 denoting the unit tensor), the Hooke elasticity can be expressed as
and its PTZ generalization (among its vast literature, see [
4,
9,
24] for its treatment in the irreversible thermodynamical internal variable approach and [
25], for its presentation in the GENERIC (General Equation for the Non- Equilibrium Reversible–Irreversible Coupling) framework), is
the coefficients will be treated as constants hereafter.
In order to make the subsequent formulae more intelligible, we introduce
and the coefficient combinations
with the aid of which (
5) gets simplified to
Taking (also for simplicity) a constant ‘isobaric’ specific heat
as well as neglected thermal expansion and heat conduction, the internal variable approach puts the following thermodynamical background behind the PTZ model: after eliminating the internal variable, its specific total energy
with absolute temperature
T, accompanied with specific entropy
s and entropy production rate density
for which the specific internal energy part
fulfils the balance
and specific entropy the balance
as can be found along the lines of [
9] (including its Appendix B), and it is straightforward to check. Because of the second law of thermodynamics,
and
follow for the PTZ model [
9], where (
16) is already apparent from the form (12). (Recall that heat conduction is neglected, so there are no heat and entropy flux terms in the balances. In parallel, there is no term in
that couples
T and
, and—correspondingly—there is no
dependent term in
s, due to neglected thermal expansion.)
Superficially, it seems redundant to also provide in the equivalent form (11). However, it is just this not-automatically-positive-definite form that will prove to be beneficial in the diagnostics of the numerical solution.
From either balance (
13) or (
14), the time derivative of temperature can also be expressed:
As a simple analysis of the PTZ model, for ‘slow’ processes, which is to be understood with respect to the time scales
a rule-of-thumb approximation is to neglect the time derivative terms (to only keep the lowest time derivative term for each quantity) in (
5). The result is nothing but the Hooke model (
4), for which the longitudinal and transversal wave propagation speeds are
Now, as opposed to this ‘static’ limit, let us consider the limit of ‘fast’ processes: then, it is the time derivative terms (the highest time derivative term for each quantity) that we keep. The result is the time derivative of an effective/‘dynamic’ Hooke model:
where the inequalities follow from (
15). Accordingly, the wave propagation speeds
follow. This, on one side, illustrates how the PTZ model can interpret that the dynamic elasticity coefficients of rocks are larger than their static counterpart [
1,
2,
3,
4]. On the other side, the nontrivial—frequency dependent, therefore, dispersive—wave propagation indicates that the numerical solution of PTZ wave propagation problems should contain the minimal possible amount of dispersion error, in order to give account of the dispersive property of the continuum model itself. In parallel, the dissipative nature of the PTZ model requires the minimal possible amount of dissipation error to reliably describe the decrease of wave amplitudes.
3. The Numerical Scheme
We take a Cartesian grid with spacings
,
,
, and time step
. Corresponding to the continuum formula (
2), we introduce the finite difference discretization
where the time index
j refers to a value at
,
to a value at
, the space index
l refers to a value at
,
m is the space index in the
y direction, and
n in the
z direction. Accordingly, stress (and strain) values reside at integer time instants, while the velocity ones are shifted in time by half; diagonal stress (and strain) components reside at integer positions, off-diagonal ones are shifted in the two directions that match with the two Cartesian indices; and, velocity components are only shifted in the direction matching with their Cartesian index (see
Figure 2). From these formulae, the
indexed velocities can be expressed explicitly (as functions of earlier quantities).
This same pattern—distribution of quantities—is used for the discretization of (
2):
from which formulae the
indexed strains can be explicitly expressed, and for the discretized version of (
5),
where
ensures second-order accuracy of the whole scheme (the proof is analogous to the one in [
8]), from which— together with (
3)—the
indexed stresses can be explicitly expressed, except for stress boundary locations, where we express strain (and know stress from the boundary condition).
Actually, regarding boundary conditions, the rule we found for both stress boundary condition and velocity (or displacement) boundary condition is that, if a quantity is missing for determining another boundary quantity, then that missing quantity is to be added outside the boundary. This also works for mixed boundary conditions, with different ones meeting at edges of a rectangular sample, for example. In what follows, we present stress boundary condition examples (relevant, e.g., for a wide class of rock mechanical applications).
The pattern of which quantity to reside where—at integer or half-integer space and time indices—could also be conjectured from the structure of the equations, but, as said in the Introduction, the spacetime viewpoint helps a lot to find this arrangement a geometrically—spacetime geometrically— natural one.
In the reversible special case of the Hooke system, this scheme is symplectic. It is actually the symplectic Euler method (in words: ‘new from old and old, new from new and old’). The interpretation is the improvement: here, new and new are shifted in time with respect to each other, so second-order accuracy is achieved, while the conventional interpretation of the symplectic Euler method is first-order only. In parallel, since mechanical energy (the Hamiltonian) is a velocity dependent term plus a strain dependent term (stress becomes a simple linear function of strain), our scheme is explicit. This also remains true at the PTZ level, so one can expect—and find, actually—a fast-running program code.
For the aspects of thermodynamics, we also discretize (
17) [here using the form (12)], explicitly expressing the
indexed temperature values from
where the notation (
6) is utilized, and the traces are to be expanded in Cartesian components and the terms containing offdiagonal components—that reside at half space-shifted locations in two indices—are averaged around the location
, first neighbours only.
5. Solution for a Cube, and the Role of Entropy Production
In the second example treated, a cube is considered, initially relaxed, as in the previous example, and then one of its sides being pressed by a cosine ‘bump’ in time, as well as in both spatial directions:
(see
Figure 12), where the notations are analogous to the ones of the previous case:
is the duration of the pulse,
is its amplitude,
is its spatial width, and
X is the length of the edges of the cube. On the other sides, normal stress is constantly zero (free surfaces), like in the previous example.
Here, we only present the PTZ model, on a grid, with all other settings being the same as in the previous example.
The solution proves similarly satisfatory, as in the former case.
Figure 13 and
Figure 14 show snapshots of the distribution of the von Mises stress invariant on two mid-planes of the cubes.
Next, in this example, we demonstrate the usefulness of another thermodynamical quantity, entropy production rate density.
It is instructive to start with showing how the idea works in the simpler, one space dimensional, setting (the rod discussed in [
8]). The one space dimensional analogue of (11) is
Let us introduce four different discretizations of this product, embodying the patterns
where
denotes the time average
.
These four versions are integrated in space and plotted in
Figure 15, the left column for a stable setting and the right column for an unstable one. The energies are also displayed. Only 25 space cells have been chosen to enhance artefacts.
Visibly, certain versions become negative when instability gets exposed. Moreover, some become negative even before that, showing, at an early stage, that there is a problem to come.
Next, let us see how the three space dimensional generalizations behave for the problem of the pressed cube: the outcomes can be seen in
Figure 16. An intentionally rude
grid is taken, the time unit is divided to 125 time steps in order to lead to a stable solution and to 100 time steps producing an unstable one. Further settings are
,
,
,
,
,
,
,
, and
,
,
.
One can see that an appropriately chosen discretization diagnoses instability.
Such numerical comparisons between stable and unstable parameter domains, which are probably enhanced by some analytical considerations, can help in the future in deciding which discretized entropy production rate is useful for diagnosing what.
7. Discussion
The numerical scheme presented here, due to its symplectic root, second-order accuracy, and the equation-friendly and spacetime geometry friendly arrangement of discretized quantities, has been found to provide reliable results in a fast and resource-friendly way. Being a finite difference scheme, it is not very flexible to simulate arbitrary shaped samples, but, already, the extension of the Cartesian formulae to cylindrical and spherical geometries promises useful applications, including the various wave-based measurement methods that are used in rock mechanics (see, e.g., [
26]), many of which rely on simple and easily treatable sample shapes. Fitting a rheological model on experimental data may require many runs so good finite difference schemes find their applicability.
The investigation of stability and dissipative and dispersion error is expected to be much more involved than in the corresponding one space dimensional situation, where the analysis was done in [
8]. Nevertheless, it is an important task for the future, because the outcomes support efficient applications of the scheme.
It is an interesting challenge to apply the presented scheme for other dissipative situations (like [
27], just to mention one example).
A systematic and general framework could be obtained, which is also beneficial for other purposes, if the spacetime background is strengthened further, by using four-quantities, four-equations on them, and formulating discretization in a fully four-geometrical way. It would, for example, help in building connection to a finite element—spacetime finite element—approach, along which way objects of general shape could also be treated. Notably, the current finite element paradigm has deficiencies and probably needs to be renewed, as indicated by results found here and earlier works [
7,
8].
In addition, the use the thermodynamical full description of a system for monitoring and controlling numerical artefacts during a computer simulation is a promising perspective. The steps made here: recognizing the usefulness of total energy and its various parts, and of entropy production rate density, are hoped to contribute to a future routine in numerical environments, where thermodynamics could be utilized.