Next Article in Journal
Lottery and Auction on Quantum Blockchain
Previous Article in Journal
Modeling Spectral Properties in Stationary Processes of Varying Dimensions with Applications to Brain Local Field Potential Signals
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Four Spacetime Dimensional Simulation of Rheological Waves in Solids and the Merits of Thermodynamics

1
Department of Energy Engineering, Faculty of Mechanical Engineering, BME, 1521 Budapest, Hungary
2
Montavid Thermodynamic Research Group, 1112 Budapest, Hungary
3
Wigner Research Centre for Physics, Department of Theoretical Physics, Institute for Particle and Nuclear Physics, 1525 Budapest, Hungary
*
Author to whom correspondence should be addressed.
Entropy 2020, 22(12), 1376; https://doi.org/10.3390/e22121376
Submission received: 1 November 2020 / Revised: 25 November 2020 / Accepted: 30 November 2020 / Published: 5 December 2020
(This article belongs to the Section Thermodynamics)

Abstract

:
The recent results attained from a thermodynamically conceived numerical scheme applied on wave propagation in viscoelastic/rheological solids are generalized here, both in the sense that the scheme is extended to four spacetime dimensions and in the aspect of the virtues of a thermodynamical approach. Regarding the scheme, the arrangement of which quantity is represented where in discretized spacetime, including the question of appropriately realizing the boundary conditions, is nontrivial. In parallel, placing the problem in the thermodynamical framework proves to be beneficial in regards to monitoring and controlling numerical artefacts—instability, dissipation error, and dispersion error. This, in addition to the observed preciseness, speed, and resource-friendliness, makes the thermodynamically extended symplectic approach that is presented here advantageous above commercial finite element software solutions.

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
ε t = 1 2 v + v ,
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
ϱ v t = σ ·
With the deviatoric and spherical parts of tensors,
σ sph = 1 3 tr σ 1 , σ dev = σ σ sph , ε sph = 1 3 tr ε 1 , ε dev = ε ε sph
(1 denoting the unit tensor), the Hooke elasticity can be expressed as
σ dev = E dev ε dev , σ sph = E sph ε sph , E dev = 2 G , E sph = 3 K ,
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
σ dev + τ dev σ dev t = E dev ε dev + E ^ dev ε dev t , σ sph + τ sph σ sph t = E sph ε sph + E ^ sph ε sph t ,
the coefficients will be treated as constants hereafter.
In order to make the subsequent formulae more intelligible, we introduce
σ ^ dev = σ dev E dev ε dev , σ ^ sph = σ sph E sph ε sph
and the coefficient combinations
I ^ dev = E ^ dev τ dev E dev , I ^ sph = E ^ sph τ sph E sph ,
with the aid of which (5) gets simplified to
σ ^ dev + τ dev σ ^ dev t = I ^ dev ε dev t , σ ^ sph + τ sph σ ^ sph t = I ^ sph ε sph t .
Taking (also for simplicity) a constant ‘isobaric’ specific heat c σ 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
e total = e kinetic + e thermal + e elastic + e rheological , e kinetic = 1 2 v 2 , e elastic = E dev 2 ϱ tr ε dev 2 + E sph 2 ϱ tr ε sph 2 , e thermal = c σ T , e rheological = τ dev 2 ϱ I ^ dev tr σ ^ dev σ dev 2 + τ sph 2 ϱ I ^ sph tr σ ^ sph σ sph 2
with absolute temperature T, accompanied with specific entropy s and entropy production rate density π s
s = c σ ln T T ref ,
(11) π s = 1 T 1 I ^ dev tr σ ^ dev E ^ dev ε t dev τ dev σ t dev + 1 I ^ sph tr σ ^ sph E ^ sph ε t sph τ sph σ t dev (12) = 1 T 1 I ^ dev tr σ ^ dev σ dev 2 + 1 I ^ sph tr σ ^ sph σ sph 2 ,
for which the specific internal energy part e total e kinetic fulfils the balance
ϱ e total e kinetic t = tr σ ε t
and specific entropy the balance
ϱ s t = π s ,
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,
I ^ dev > 0 , I ^ sph > 0
and
π s 0
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 e total 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 π s 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:
T t = T ϱ c σ π s .
As a simple analysis of the PTZ model, for ‘slow’ processes, which is to be understood with respect to the time scales
τ dev , τ ^ dev = E ^ dev / E dev , τ sph , τ ^ sph = E ^ sph / E sph ,
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
c longitudinal = 2 E dev + E sph 3 ϱ , c transversal = E dev 2 ϱ .
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:
σ dev = E dev ε dev , σ sph = E sph ε sph , E dev = E ^ dev / τ dev > E dev , E sph = E ^ sph / τ sph > E sph ,
where the inequalities follow from (15). Accordingly, the wave propagation speeds
c ^ longitudinal = 2 E ^ dev + E ^ sph 3 ϱ > c longitudinal , c ^ transversal = E ^ dev 2 ϱ > c transversal
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 Δ x , Δ y , Δ z , and time step Δ t . Corresponding to the continuum formula (2), we introduce the finite difference discretization
ϱ v x l + 1 2 , m , n j + 1 2 v x l + 1 2 , m , n j 1 2 Δ t = σ x x l + 1 , m , n j σ x x l , m , n j Δ x + σ x y l + 1 2 , m + 1 2 , n j σ x y l + 1 2 , m 1 2 , n j Δ y + σ x z l + 1 2 , m , n + 1 2 j σ x z l + 1 2 , m , n 1 2 j Δ z ,
ϱ v y l , m + 1 2 , n j + 1 2 v y l , m + 1 2 , n j 1 2 Δ t = σ y x l + 1 2 , m + 1 2 , n j σ y x l 1 2 , m + 1 2 , n j Δ x + σ y y l , m + 1 , n j σ y y l , m , n j Δ y + σ y z l , m + 1 2 , n + 1 2 j σ y z l , m + 1 2 , n 1 2 j Δ z ,
ϱ v z l , m , n + 1 2 j + 1 2 v z l , m , n + 1 2 j 1 2 Δ t = σ z x l + 1 2 , m , n + 1 2 j σ z x l 1 2 , m , n + 1 2 j Δ x + σ z y l , m + 1 2 , n + 1 2 j σ z y l , m 1 2 , n + 1 2 j Δ y + σ z z l , m , n + 1 j σ z z l , m , n j Δ z ,
where the time index j refers to a value at t j = j · Δ t , j + 1 2 to a value at t j + 1 2 = j + 1 2 · Δ t , the space index l refers to a value at x l = l · Δ x , 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 j + 1 2 indexed velocities can be expressed explicitly (as functions of earlier quantities).
This same pattern—distribution of quantities—is used for the discretization of (2):
ε x x l , m , n j + 1 ε x x l , m , n j Δ t = v x l + 1 2 , m , n j + 1 2 v x l 1 2 , m , n j + 1 2 Δ x ,
ε y y l , m , n j + 1 ε y y l , m , n j Δ t = v y l , m + 1 2 , n j + 1 2 v y l , m 1 2 , n j + 1 2 Δ x ,
ε z z l , m , n j + 1 ε z z l , m , n j Δ t = v z l , m , n + 1 2 j + 1 2 v z l , m , n 1 2 j + 1 2 Δ x ,
ε x y l + 1 2 , m + 1 2 , n j + 1 ε x y l + 1 2 , m + 1 2 , n j Δ t = 1 2 v x l + 1 2 , m + 1 , n j + 1 2 v x l + 1 2 , m , n j + 1 2 Δ y + v y l + 1 , m + 1 2 , n j + 1 2 v y l , m + 1 2 , n j + 1 2 Δ x ,
ε x z l + 1 2 , m , n + 1 2 j + 1 ε x z l + 1 2 , m , n + 1 2 j Δ t = 1 2 v x l + 1 2 , m , n + 1 j + 1 2 v x l + 1 2 , m , n j + 1 2 Δ z + v z l + 1 , m , n + 1 2 j + 1 2 v z l , m , n + 1 2 j + 1 2 Δ x ,
ε y z l , m + 1 2 , n + 1 2 j + 1 ε y z l , m + 1 2 , n + 1 2 j Δ t = 1 2 v y l , m + 1 2 , n + 1 j + 1 2 v y l , m + 1 2 , n j + 1 2 Δ z + v z l , m + 1 , n + 1 2 j + 1 2 v z l , m , n + 1 2 j + 1 2 Δ y ,
from which formulae the j + 1 indexed strains can be explicitly expressed, and for the discretized version of (5),
α σ p q dev l , m , n j + ( 1 α ) σ p q dev l , m , n j + 1 + σ p q dev l , m , n j + 1 σ p q dev l , m , n j Δ t = E dev α ε p q dev l , m , n j + ( 1 α ) ε p q dev l , m , n j + 1 + E ^ dev ε p q dev l , m , n j + 1 ε p q dev l , m , n j Δ t ,
α σ p q sph l , m , n j + ( 1 α ) σ p q sph l , m , n j + 1 + σ p q sph l , m , n j + 1 σ p q sph l , m , n j Δ t = E sph α ε p q sph l , m , n j + ( 1 α ) ε p q sph l , m , n j + 1 + E ^ sph ε p q sph l , m , n j + 1 ε p q sph l , m , n j Δ t , p , q = x , y , z , l , m , n = integers or half - integers depending on p , q ,
where α = 1 / 2 ensures second-order accuracy of the whole scheme (the proof is analogous to the one in [8]), from which— together with (3)—the j + 1 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 1 from old 1 and old 2 , new 2 from new 1 and old 2 ’). The interpretation is the improvement: here, new 1 and new 2 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 j + 1 2 indexed temperature values from
T l , m , n j + 1 2 T l , m , n j 1 2 Δ t = 1 ϱ c σ 1 I ^ dev tr σ ^ dev σ dev 2 l , m , n j + 1 I ^ sph tr σ ^ sph σ sph 2 l , m , n j ,
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 l , m , n , first neighbours only.

4. Solution for a Rectangular Beam, and the Role of Total Energy

The first example on which we demonstrate the scheme is a square cross-sectioned long beam, being thus treated as a plane-strain problem. Initially, the beam is in relaxed/equilibrium state (zero stress, strain and velocity, and homogeneous temperature). Subsequently, on one of its sides, a single normal stress pulse is applied, with profile
σ y y t , x , 0 , z = σ b 1 2 1 cos 2 π t τ b · 1 2 1 cos 2 π x X / 2 W b if 0 t τ b and X W b 2 x X + W b 2 , 0 otherwise
(see Figure 3), where τ b is the duration of the pulse, σ b is its amplitude, W b is its spatial width, and X is the width of the beam. On the other sides, normal stress is constantly zero (free surfaces).
A 50 × 50 grid is considered in the xy plane, for 240 time steps, where the time step is the largest at which stability is maintained. Notably, stability investigation is fairly involved for this problem and it requires a separate whole study. Further settings (in appropriate units) are ϱ = 1 , E dev = 4 , E sph = 10 , τ b = 0.3 , σ b = 5 , W b / X = 0.6 .

4.1. Hooke Case

Figure 4 and Figure 5 show snapshots of a stress component distribution and a velocity component.
In a movie format, it is more spectacular how reliably the simulation performs.
Furthermore, it is not only the eye that could judge the reliability: with the help of thermodynamics, energy—in the Hooke case, mechanical energy— proves to be a useful diagnostic tool:
  • if it explodes then there is instability;
  • if it deviates from a constant then there is dissipation error; and,
  • if it is wavy/oscillating then there is dispersion error.
The scheme presented here also functions satisfactorily in this aspect, as displayed in Figure 6. For energy, we perform summation, over the integer centred discrete cells, of the energy terms discretized along the above lines, including that averages, like for (33), are taken wherever necessary, also in the time direction (for kinetic energy).

4.2. PTZ Case

In a PTZ medium, the solution of the analogous problem is similarly good. Figure 7, Figure 8, Figure 9 and Figure 10 present snapshots, where Figure 9 and Figure 10 display two further quantities: tr σ ^ dev σ dev 2 (essentially the Huber–Mises–Hencky or von Mises equivalent stress) and temperature. Dissipation is nicely indicated via temperature. Settings (in addition to the above ones for the Hooke case) are E ^ dev = 1.4 , τ dev = 0.2 , c σ = 0.001 , T 0 = 0.1 , α = 1 2 .
The diagnostic role of the various energies, and especially their sum, is a great help again for checking whether the simulation performs acceptably. Figure 11 illustrates how the scheme that is introduced above behaves in this respect.

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:
σ y y t , x , 0 , z = σ b { 1 2 1 cos 2 π t τ b · 1 2 1 cos 2 π x X / 2 W b σ b { 1 2 1 cos 2 π t τ b · 1 2 1 cos 2 π z X / 2 W b } if 0 t τ b and X W b 2 x , z X + W b 2 , 0 otherwise
(see Figure 12), where the notations are analogous to the ones of the previous case: τ b is the duration of the pulse, σ b is its amplitude, W b 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 25 × 25 × 25 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
π s = 1 T 1 I ^ σ ^ · E ^ ε t τ σ t
Let us introduce four different discretizations of this product, embodying the patterns
  • old · ( new old ) ,
  • old · ( new older ) ,
  • new · ( new old ) , and
  • new · ( new older ) :
1 T n j 1 I ^ σ ^ n j · E ^ ε n j + 1 ε n j Δ t τ σ n j + 1 σ n j Δ t ,
1 T n j 1 I ^ σ ^ n j · E ^ ε n j + 1 ε n j 1 2 Δ t τ σ n j + 1 σ n j 1 2 Δ t ,
1 T n j 1 I ^ σ ^ n j + 1 · E ^ ε n j + 1 ε n j Δ t τ σ n j + 1 σ n j Δ t ,
1 T n j 1 I ^ σ ^ n j + 1 · E ^ ε n j + 1 ε n j 1 2 Δ t τ σ n j + 1 σ n j 1 2 Δ t ,
where T n j denotes the time average T n j 1 2 + T n j + 1 2 / 2 .
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 10 × 10 × 10 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 ϱ = 1 , E dev = 3 , E sph = 5 , E ^ dev = 20 , τ dev = 0.391 , c σ = 0.001 , T 0 = 0.1 , τ b = 0.25 , and σ b = 3 , W b / X = 0.6 , α = 1 2 .
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.

6. Two-Dimensional Wave Propagation According to the Finite Element Software COMSOL

In [8], a comparison of solutions via the one space dimensional scheme with corresponding finite element solutions was presented. Repeating it for the present three space dimensional scheme would be informative. We found it advisable to start investigating the higher dimensional wave propagation describing the possibilities of finite element softwares in a much simpler setting than the one treated above, based on the (discouraging) experience that is gained with the one-dimensional case, and since a rheological model like PTZ in the deviatoric–spherical formulation is difficult to translate to the capabilities of classic finite element softwares.
Namely, we consider the wave equation in two spatial dimensions, for a single scalar degree of freedom u = u ( t , x , y ) , with constant (unit) coefficients and no source term in the equation:
2 u t 2 = 2 u x 2 + 2 u y 2 .
A rectangular sample is taken with ‘flux’ (Neumann) boundary condition (normal spatial derivative component is prescribed): on one of the edges, one cosine-type pulse of the kind (34) is applied, while the other edges are kept ‘free’ (zero normal derivative). Initially, both u and u / t are set to zero (‘relaxed initial state’).
To stay similar to the numerical calculations that are presented here, a 50 × 50 spatial grid is taken. The software that we use is COMSOL v5.3a, where this wave equation problem is a built-in possibility. The pulse duration is 0.3 time unit. Figure 17 displays the resulting spatial distribution of the scalar degree of freedom after two units of time, obtained with various available settings of COMSOL v5.3a.
One can see that the solution depends on the settings very seriously. There are large-scale differences and fine-structured irregularities. Moreover, one has no means of validation which outcome is correct to what extent.
When taking into account that our full problem has 16 coupled degrees of freedom in three spatial dimensions and with further time derivatives (in the PTZ model), it does not seem to be reasonable to try to represent it via any commercial finite element software as long as such a much simpler problem cannot be confidently treated.

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.

Author Contributions

Á.P.: numerical simulations, figures, analysis. M.S.: analysis, numerical simulations, figures. R.K.: numerical simulations, figures, manuscript text, funding. T.F.: conceptual idea, numerical simulations, figures, manuscript text. All authors have read and agreed to the published version of the manuscript.

Funding

The research reported in this paper and carried out at BME has been supported by the grants National Research, Development and Innovation Office-NKFIH KH 130378, and K124366(124508), and by the NRDI Fund (TKP2020 NC, Grant No. BME-NC) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry for Innovation and Technology. This paper was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

Acknowledgments

The authors are thankful to Tamás Kovács for the discussions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Barnaföldi, G.G.; Bulik, T.; Cieslar, M.; Dávid, E.; Dobróka, M.; Fenyvesi, E.; Gondek-Rosinska, D.; Gráczer, Z.; Hamar, G.; Huba, G.; et al. First report of long term measurements of the MGGL laboratory in the Mátra mountain range. Class. Quantum Gravity 2017, 34, 114001. [Google Scholar] [CrossRef]
  2. Ván, P.; Barnaföldi, G.G.; Bulik, T.; Biró, T.; Czellár, S.; Cieślar, M.; Czanik, C.; Dávid, E.; Debreceni, E.; Denys, M.; et al. Long term measurements from the Mátra Gravitational and Geophysical Laboratory. Eur. Phys. J. 2019, 228, 1693–1734. [Google Scholar] [CrossRef] [Green Version]
  3. Davarpanah, S.M.; Ván, P.; Vásárhelyi, B. Investigation of relationship between dynamic and static deformation moduli of rocks. Geomech. Geophys. Geo-Energy Geo-Resour. 2020, 6, 29. [Google Scholar] [CrossRef] [Green Version]
  4. Fülöp, T. Wave propagation in rocks—Investigating the effect of rheology. Period. Polytech. Civ. Eng. 2020. to appear. [Google Scholar] [CrossRef]
  5. Fülöp, T.; Szücs, M. A solution method for determining rheological time dependence around tunnels. In Proceedings of the EUROCK2020, Trondheim, Norway, 12–14 October 2020. [Google Scholar]
  6. Fülöp, T.; Szücs, M. Analytical solution method for rheological problems of solids. arXiv 2018, arXiv:1810.06350. [Google Scholar]
  7. Rieth, Á.; Kovács, R.; Fülöp, T. Implicit numerical schemes for generalized heat conduction equations. Int. J. Heat Mass Transf. 2018, 126, 1177–1182. [Google Scholar] [CrossRef] [Green Version]
  8. Fülöp, T.; Kovács, R.; Szücs, M.; Fawaier, M. Thermodynamical extension of a symplectic numerical scheme with half space and time shifts demonstrated on rheological waves in solids. Entropy 2020, 22, 155. [Google Scholar] [CrossRef] [Green Version]
  9. Asszonyi, C.; Fülöp, T.; Ván, P. Distinguished rheological models for solids in the framework of a thermodynamical internal variable theory. Contin. Mech. Thermodyn. 2015, 27, 971–986. [Google Scholar] [CrossRef] [Green Version]
  10. Hairer, E.; Lubich, C.; Wanner, G. Geometric Numerical Integration, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  11. Zinner, C.P.; Öttinger, H.C. Numerical stability with help from entropy: Solving a set of 13 moment equations for shock tube problem. J. Non-Equilib. Thermodyn. 2019, 44, 43–69. [Google Scholar] [CrossRef] [Green Version]
  12. Shang, X.; Öttinger, H.C. Structure-preserving integrators for dissipative systems based on reversible-irreversible splitting. arXiv 2018, arXiv:1804.05114. [Google Scholar] [CrossRef] [Green Version]
  13. Portillo, D.; García Orden, J.C.; Romero, I. Energy-Entropy-Momentum integration schemes for general discrete non-smooth dissipative problems in thermomechanics. Int. J. Numer. Methods Eng. 2017, 112, 776–802. [Google Scholar] [CrossRef]
  14. Vermeeren, M.; Bravetti, A.; Seri, M. Contact variational integrators. J. Phys. A Math. Theor. 2019, 52, 445206. [Google Scholar] [CrossRef] [Green Version]
  15. Gay-Balmaz, F.; Yoshimura, H. Variational discretization of the nonequilibrium thermodynamics of simple systems. Nonlinearity 2018, 31, 1673. [Google Scholar] [CrossRef] [Green Version]
  16. Couéraud, B.; Gay-Balmaz, F. Variational discretization of thermodynamical simple systems on Lie groups. Discret. Cont. Dyn. Syst. S 2020, 13. [Google Scholar] [CrossRef] [Green Version]
  17. Romero, I. Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part I: Monolithic integrators and their application to finite strain thermoelasticity. Comput. Methods Appl. Mech. Eng. 2010, 199, 1841–1858. [Google Scholar] [CrossRef]
  18. Romero, I. Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part II: Fractional step methods. Comput. Methods Appl. Mech. Eng. 2010, 199, 2235–2248. [Google Scholar] [CrossRef]
  19. Berezovski, A.; Ván, P. Internal Variables in Thermoelasticity; Springer: Cham, Switzerland, 2017. [Google Scholar]
  20. Janečka, A.; Málek, J.; Průša, V.; Tierra, G. Numerical scheme for simulation of transient flows of non-Newtonian fluids characterised by a non-monotone relation between the symmetric part of the velocity gradient and the Cauchy stress tensor. Acta Mech. 2019, 230, 729–747. [Google Scholar] [CrossRef] [Green Version]
  21. Fülöp, T.; Ván, P. Kinematic quantities of finite elastic and plastic deformation. Math. Methods Appl. Sci. 2012, 35, 1825–1841. [Google Scholar] [CrossRef] [Green Version]
  22. Ván, P. Galilean relativistic fluid mechanics. Contin. Mech. Thermodyn. 2017, 29, 585–610. [Google Scholar] [CrossRef] [Green Version]
  23. Matolcsi, T. Ordinary Thermodynamics: Nonequilibrium Homogeneous Processes; Akadémiai Kiadó (Publishing House of the Hungarian Academy of Sciences): Budapest, Hungary, 2004. [Google Scholar]
  24. Asszonyi, C.; Csatár, A.; Fülöp, T. Elastic, thermal expansion, plastic and rheological processes—Theory and experiment. Period. Civ. Eng. 2016, 60, 591–601. [Google Scholar] [CrossRef] [Green Version]
  25. Szücs, M.; Fülöp, T. Kluitenberg–Verhás rheology of solids in the GENERIC framework. J. Non-Equilib. Thermodyn. 2019, 44, 247–259. [Google Scholar] [CrossRef] [Green Version]
  26. Malhotra, V.M.; Carino, N.J. (Eds.) Handbook on Nondestructive Testing of Concrete, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2003. [Google Scholar] [CrossRef]
  27. Szabó, B.; Kossa, A. Characterization of impacts of elastic-plastic spheres. Period. Polytech. Mech. Eng. 2020, 64, 165–171. [Google Scholar] [CrossRef] [Green Version]
Figure 1. An excitation pulse, generated at the left endpoint of a finite-size one space dimensional Hookean sample, regularly arrives at the right endpoint. First row: the results from the scheme that was introduced in [8] for two different pulse lengths. Second and third rows: the corresponding results obtained by the finite element software COMSOL [8], at various settings: left column second row: Backward Differentiation Formula (BDF) Maximum order 2, third row: BDF Maximum order 5; right column second row: Dormand–Prince (DP) 5, third row: Runge–Kutta (RK) 34. In each finite element solution, dissipation error (decrease of the amplitude) and dispersion error (artificial oscillations) are both observable, even during the first three bounces. Meanwhile, the pulses in the first row keep their shape even after many bounces [8].
Figure 1. An excitation pulse, generated at the left endpoint of a finite-size one space dimensional Hookean sample, regularly arrives at the right endpoint. First row: the results from the scheme that was introduced in [8] for two different pulse lengths. Second and third rows: the corresponding results obtained by the finite element software COMSOL [8], at various settings: left column second row: Backward Differentiation Formula (BDF) Maximum order 2, third row: BDF Maximum order 5; right column second row: Dormand–Prince (DP) 5, third row: Runge–Kutta (RK) 34. In each finite element solution, dissipation error (decrease of the amplitude) and dispersion error (artificial oscillations) are both observable, even during the first three bounces. Meanwhile, the pulses in the first row keep their shape even after many bounces [8].
Entropy 22 01376 g001
Figure 2. Spatial arrangement of the discretized quantities (two-dimensional projection). The circles stand for diagonal tensor components, squares for offdiagonal ones, and triangles for vector components, different components with differently oriented triangles. Void quantities are prescribed by boundary condition (in the case stress boundary conditions are considered, like here.)
Figure 2. Spatial arrangement of the discretized quantities (two-dimensional projection). The circles stand for diagonal tensor components, squares for offdiagonal ones, and triangles for vector components, different components with differently oriented triangles. Void quantities are prescribed by boundary condition (in the case stress boundary conditions are considered, like here.)
Entropy 22 01376 g002
Figure 3. The spatial distribution of normal stress as the boundary condition on one side of a square cross-sectioned beam that is infinitely long in the z direction. The excitation is a single cosine-shaped ‘bump’ in time, too.
Figure 3. The spatial distribution of normal stress as the boundary condition on one side of a square cross-sectioned beam that is infinitely long in the z direction. The excitation is a single cosine-shaped ‘bump’ in time, too.
Entropy 22 01376 g003
Figure 4. Distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Hooke case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Figure 4. Distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Hooke case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Entropy 22 01376 g004
Figure 5. Continuation of Figure 4: distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Hooke case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Figure 5. Continuation of Figure 4: distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Hooke case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Entropy 22 01376 g005aEntropy 22 01376 g005b
Figure 6. Mechanical energy types as functions of time, for the Hooke case.
Figure 6. Mechanical energy types as functions of time, for the Hooke case.
Entropy 22 01376 g006
Figure 7. Distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Poynting–Thomson–Zener (PTZ) case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Figure 7. Distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the Poynting–Thomson–Zener (PTZ) case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Entropy 22 01376 g007
Figure 8. Continuation of Figure 7: distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Figure 8. Continuation of Figure 7: distribution of a stress component (left column) and of a velocity component (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Entropy 22 01376 g008aEntropy 22 01376 g008b
Figure 9. Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) (left column) and of temperature (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Figure 9. Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) (left column) and of temperature (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Entropy 22 01376 g009aEntropy 22 01376 g009b
Figure 10. Continuation of Figure 9: Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) (left column) and of temperature (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Figure 10. Continuation of Figure 9: Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) (left column) and of temperature (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 4 τ b , respectively.
Entropy 22 01376 g010aEntropy 22 01376 g010b
Figure 11. Total energy and the various energy types as functions of time, for the PTZ case.
Figure 11. Total energy and the various energy types as functions of time, for the PTZ case.
Entropy 22 01376 g011
Figure 12. The spatial distribution of normal stress as the boundary condition on one side of a cube. The excitation is also a single cosine-shaped ‘bump’ in time.
Figure 12. The spatial distribution of normal stress as the boundary condition on one side of a cube. The excitation is also a single cosine-shaped ‘bump’ in time.
Entropy 22 01376 g012
Figure 13. Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) at the mid-plane parallel to the excited side (left column) and at a mid-plane orthogonal to it (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Figure 13. Snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) at the mid-plane parallel to the excited side (left column) and at a mid-plane orthogonal to it (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 1 / 2 ) τ b , τ b , ( 3 / 2 ) τ b , 2 τ b , respectively.
Entropy 22 01376 g013aEntropy 22 01376 g013b
Figure 14. Continuation of Figure 13: snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) at the mid-plane parallel to the excited side (left column) and at a mid-plane orthogonal to it (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 8 τ b , respectively.
Figure 14. Continuation of Figure 13: snapshots of the distribution of the stress invariant tr ( σ ^ dev 2 ) at the mid-plane parallel to the excited side (left column) and at a mid-plane orthogonal to it (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( 5 / 2 ) τ b , 3 τ b , ( 7 / 2 ) τ b , 8 τ b , respectively.
Entropy 22 01376 g014
Figure 15. Entropy production rate according to the four discretized versions (37)–(40), respectively, for a time step resulting in a stable run (left column) (namely, when the rheological Courant number c ^ Δ t / Δ x is 1, where c ^ = E ^ / ( τ ϱ ) is the high-frequency rheological wave propagation speed) and for a choice leading to an unstable outcome (right column) (rheological Courant number c ^ Δ t / Δ x = 1.003 ). PTZ model with τ E / E ^ = 0.25 ; the sample length X is such that E ^ / E = 5 X / c with the low-frequency/elastic wave propagation speed c = E / ϱ . See [8] for further details.
Figure 15. Entropy production rate according to the four discretized versions (37)–(40), respectively, for a time step resulting in a stable run (left column) (namely, when the rheological Courant number c ^ Δ t / Δ x is 1, where c ^ = E ^ / ( τ ϱ ) is the high-frequency rheological wave propagation speed) and for a choice leading to an unstable outcome (right column) (rheological Courant number c ^ Δ t / Δ x = 1.003 ). PTZ model with τ E / E ^ = 0.25 ; the sample length X is such that E ^ / E = 5 X / c with the low-frequency/elastic wave propagation speed c = E / ϱ . See [8] for further details.
Entropy 22 01376 g015aEntropy 22 01376 g015b
Figure 16. Entropy production rate according to the four discretized versions (37)–(40) generalized to three space dimensions, respectively, computed for a PTZ cube, with a time step resulting in a stable run (left column) and for a choice leading to an unstable outcome (right column).
Figure 16. Entropy production rate according to the four discretized versions (37)–(40) generalized to three space dimensions, respectively, computed for a PTZ cube, with a time step resulting in a stable run (left column) and for a choice leading to an unstable outcome (right column).
Entropy 22 01376 g016
Figure 17. u ( 2 , x , y ) with various COMSOL settings (cf. Figure 1). Left column, from top to bottom: BDF Max. order 5 & min. order 2, BDF Max. order 3 & min. order 2, BDF Max. order 1 & min. order 1, BDF Max. order 5 & min. order 3 Strict time stepping, RK34 Default. Right column, from top to bottom: DP5 PI-Smooth, DP5 PI-Standard, DP5 PI-Quick, Generalized- α (GA), GA Predictor Constant. Not shown (failed like DP5 PI-Standard or GA Predictor Constant): DP5 Pi-Off, RK34 Manual, Cash-Karp 5 Free or Strict (Manual with a small time step gives a result similar to that of RK Default).
Figure 17. u ( 2 , x , y ) with various COMSOL settings (cf. Figure 1). Left column, from top to bottom: BDF Max. order 5 & min. order 2, BDF Max. order 3 & min. order 2, BDF Max. order 1 & min. order 1, BDF Max. order 5 & min. order 3 Strict time stepping, RK34 Default. Right column, from top to bottom: DP5 PI-Smooth, DP5 PI-Standard, DP5 PI-Quick, Generalized- α (GA), GA Predictor Constant. Not shown (failed like DP5 PI-Standard or GA Predictor Constant): DP5 Pi-Off, RK34 Manual, Cash-Karp 5 Free or Strict (Manual with a small time step gives a result similar to that of RK Default).
Entropy 22 01376 g017aEntropy 22 01376 g017b
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Pozsár, Á.; Szücs, M.; Kovács, R.; Fülöp, T. Four Spacetime Dimensional Simulation of Rheological Waves in Solids and the Merits of Thermodynamics. Entropy 2020, 22, 1376. https://doi.org/10.3390/e22121376

AMA Style

Pozsár Á, Szücs M, Kovács R, Fülöp T. Four Spacetime Dimensional Simulation of Rheological Waves in Solids and the Merits of Thermodynamics. Entropy. 2020; 22(12):1376. https://doi.org/10.3390/e22121376

Chicago/Turabian Style

Pozsár, Áron, Mátyás Szücs, Róbert Kovács, and Tamás Fülöp. 2020. "Four Spacetime Dimensional Simulation of Rheological Waves in Solids and the Merits of Thermodynamics" Entropy 22, no. 12: 1376. https://doi.org/10.3390/e22121376

APA Style

Pozsár, Á., Szücs, M., Kovács, R., & Fülöp, T. (2020). Four Spacetime Dimensional Simulation of Rheological Waves in Solids and the Merits of Thermodynamics. Entropy, 22(12), 1376. https://doi.org/10.3390/e22121376

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop