Four spacetime dimensional simulation of rheological waves in solids and the merits of thermodynamics
AArticle
Four Spacetime Dimensional Simulation of RheologicalWaves in Solids and the Merits of Thermodynamics
Áron Pozsár , Mátyás Szücs , Róbert Kovács * and Tamás Fülöp Department of Energy Engineering, Faculty of Mechanical Engineering, BME, 1521 Budapest, Hungary;[email protected] (Á.P.); [email protected] (M.S.); [email protected] (T.F.) Montavid Thermodynamic Research Group, 1112 Budapest, Hungary Department of Theoretical Physics, Wigner Research Centre for Physics, Institute for Particle and NuclearPhysics, 1525 Budapest, Hungary * Correspondence: [email protected]
Abstract:
The recent results attained from a thermodynamically conceived numerical scheme appliedon wave propagation in viscoelastic/rheological solids are generalized here, both in the sense that thescheme is extended to four spacetime dimensions and in the aspect of the virtues of a thermodynamicalapproach. Regarding the scheme, the arrangement of which quantity is represented where in discretizedspacetime, including the question of appropriately realizing the boundary conditions, is nontrivial. Inparallel, placing the problem in the thermodynamical framework proves to be beneficial in regards tomonitoring and controlling numerical artefacts—instability, dissipation error, and dispersion error. This,in addition to the observed preciseness, speed, and resource-friendliness, makes the thermodynamicallyextended symplectic approach that is presented here advantageous above commercial finite elementsoftware solutions.
Keywords: symplectic numerical methods; rheology; solids; waves; spacetime; thermodynamics
1. Introduction
Solids may be less “solid” than expected. Beyond elastic behaviour, they may exhibit damped anddelayed response. This viscoelastic/rheological reaction may not be simply explained by a viscosity-relatedadditional stress (the Kelvin–Voigt model of rheology), but the time derivative of stress may also be neededin 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 describingboth creep (declining increase of strain during constant stress) and relaxation (declining decrease of stressduring constant strain), as well as the simplest one, via which it is possible to interpret that the dynamicelasticity coefficients of rocks are different from, and larger than, their static counterpart [1–4]. Related tothe latter aspect, high-frequency waves have a larger propagation speed in PTZ media than low-frequencyones [4], which makes this model relevant for, e.g., seismic phenomena and acoustic rock mechanicalmeasurement 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 oftransients and waves. Incorporating such ‘fast’ effects is expected to only be realizable by means ofnumerical calculations in most practical situations.In many of the practical applications, such a numerical calculation must be performed many timeswith different material coefficients, for example, as part of a fitting procedure where the experimental dataare to be fitted. Hence, the numerical scheme should be fast, resource-friendly, yet reliable and preciseenough. a r X i v : . [ phy s i c s . c l a ss - ph ] D ec of 29 In addition, numerical calculations face three frequent challenges: instability (exponential blow-upof 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 commercialfinite element softwares for wave phenomena in PTZ models. What we found—already for the Hookeancase (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 thementioned kinds.Now, if a numerical scheme exhibits dissipation error for conservative systems, then it is expected tobehave similarly for nonconservative ones, so that one cannot separate the real dissipation of mechanicalenergy from the dissipation artefact of the scheme. Additionally, similarly, a real wavy behaviour cannotbe distinguished from the dispersion/wavy artefact.These have motivated us to develop our own numerical scheme, which performs better [8]. Similarlyto that the PTZ model can be obtained in a thermodynamical approach as an internal variable extension ofHooke elasticity [9], our starting point was a symplectic scheme for Hooke elasticity. Symplectic numericalschemes (see, e.g., [10]) provide much better large-time approximations, thanks to the fact that a symplecticnumerical integrator of a conservative system is actually the exact integrator of a ‘nearby’ (coinciding inthe zeroth order of the time step) conservative system.In recent years, numerous works have been born in order to develop extensions of symplectic schemesto nonconservative systems [11–20]. We also took this path, and devised such an extension, on the exampleof the PTZ model, in one space dimension [8], with the novelties that some of the discretized field valuesreside with half space and time steps with respect to some other field values. This made the symplecticEuler method—an originally order-one accurate scheme—accurate to a second order, and the spatialaccuracy 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 modelallows 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 generalnonconservative 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 forstability, 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 inone 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 tothree space dimensions (3D), and the other is exploiting the whole thermodynamical theory around thePTZ 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 vectorialand tensorial degrees of freedom.The 3D scheme is designed in order to keep the nice, second-order, behaviour of the discretizationin both the spatial and in temporal direction. Achieving this is not so trivial—different components ofvectors 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 equationsbecomes closed. We succeed in finding a rule for this that is general enough to hold for both stress anddisplacment 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 perspectivehas helped us a lot. Specifically, on one side, thermodynamical balances in their differential form arefour-divergences from the spacetime aspect, and they have an integral counterpart, which, via Gauss’ of 29 theorem, helps one to find where to represent which flux-type quantity. Additionally, knowing that, fromthe spacetime point of view, velocity is a timelike four-vector, [21,22] gives the information that velocityshould be shifted not only spatially, but also temporally. Oppositely, stress is a spacelike tensor, so notemporal shift is needed.
Figure 1.
An excitation pulse, generated at the left endpoint of a finite-size one space dimensionalHookean sample, regularly arrives at the right endpoint. First row: the results from the scheme thatwas introduced in [8] for two different pulse lengths. Second and third rows: the corresponding resultsobtained by the finite element software COMSOL [8], at various settings: left column second row: BackwardDifferentiation Formula (BDF) Maximum order 2, third row: BDF Maximum order 5; right column secondrow: Dormand–Prince (DP) 5, third row: Runge–Kutta (RK) 34. In each finite element solution, dissipationerror (decrease of the amplitude) and dispersion error (artificial oscillations) are both observable, evenduring the first three bounces. Meanwhile, the pulses in the first row keep their shape even after manybounces [8].
In parallel, thermodynamics is not only important from the aspects of balances. Namely, commercialfinite element softwares only focus on the set of equations to solve, i.e., on the minimally necessaryequations to follow the minimally necessary quantities. However, knowing, from the spacetimeperspective, that momentum and energy form a four-quantity (also in continuum theory on Galilean of 29 spacetime) [22], in addition to the customarily taken balance of momentum, the balance of energyis 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 wherethermomechanical coupling and heat conduction are neglected) is that, if, say, temperature is followedvia a separate discretized evolution equation, then the conservation of total energy—at the discretizedlevel—is not built-in, but rather is a property that will hold only approximately. Subsequently, checkinghow well this conservation holds along the numerical solution can provide a diagnostic tool. Thus, onemay check the degree of dissipation error (i.e., the degree of violation of total energy conservation) and ofdispersion error (spurious oscillations on total energy that should be a constant). This idea is demonstratedbelow, 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 numericallevel. Accordingly, entropy, and the related entropy production, may serve as an aid for reliable numericalcalculations. Certain efforts in this direction have already been made [11,12]. Here, we introduce anotherway of utilizing this general idea.Specifically, we focus on entropy production. At the continuum level, it must be positive definiteaccording to the second law of thermodynamics. However, when discretized, this property may alsobecome challenged. Naturally, if an explicitly positive definite expression is discretized, then it remainspositive definite. However, alternative forms—which only turn out to be positive definite when the furtherthermodynamical equations also hold—are not ab ovo positive definite and, correspondingly, may fail inbeing/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 rheologicalmodels, like the PTZ one, are derived from the internal variable approach [9]. Here, we discretize such anexpression of entropy production and show that its value becoming negative can forecast a loss of stabilityand 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 respectto an inertial reference system), due to which we do not have to differentiate between Eulerian andLagrangian position or make a distinction between spatial spacetime vectors, covectors, tensors, etc, aswell as material manifold related ones, mass density (cid:36) can also be treated as constant, and the relationshipbetween the symmetric strain tensor ε to the velocity field v is ∂ ε ∂ t = (cid:16) → ∇ ⊗ v + v ⊗ ← ∇ (cid:17) , (1)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 accordingto (cid:36) ∂ v ∂ t = σ · ← ∇ (2)With the deviatoric and spherical parts of tensors, σ sph = ( tr σ ) , σ dev = σ − σ sph , ε sph = ( tr ε ) , ε dev = ε − ε sph (3) of 29 ( denoting the unit tensor), the Hooke elasticity can be expressed as σ dev = E dev ε dev , σ sph = E sph ε sph , E dev = G , E sph = K , (4)and its PTZ generalization {among its vast literature, see [4,9,24] for its treatment in the irreversiblethermodynamical 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 , (5)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 (6)and the coefficient combinationsˆ I dev = ˆ E dev − τ dev E dev , ˆ I sph = ˆ E sph − τ sph E sph , (7)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 . (8)Taking (also for simplicity) a constant ‘isobaric’ specific heat c σ as well as neglected thermal expansionand heat conduction, the internal variable approach puts the following thermodynamical backgroundbehind the PTZ model: after eliminating the internal variable, its specific total energy e total = e kinetic + e thermal + e elastic + e rheological , (9) e kinetic = v , e elastic = E dev (cid:36) tr (cid:16) ε dev2 (cid:17) + E sph (cid:36) tr (cid:16) ε sph2 (cid:17) , e thermal = c σ T , e rheological = τ dev (cid:36) ˆ I dev tr (cid:16) ˆ σ dev2 (cid:17) + τ sph (cid:36) ˆ I sph tr (cid:16) ˆ σ sph2 (cid:17) with absolute temperature T , accompanied with specific entropy s and entropy production rate density π s s = c σ ln TT ref , (10) π s = T (cid:40) I dev tr (cid:34) ˆ σ dev (cid:32) ˆ E dev ∂ ε ∂ t dev − τ dev ∂ σ ∂ t dev (cid:33)(cid:35) + I sph tr (cid:34) ˆ σ sph (cid:32) ˆ E sph ∂ ε ∂ t sph − τ sph ∂ σ ∂ t dev (cid:33)(cid:35)(cid:41) (11) = T (cid:26) I dev tr (cid:16) ˆ σ dev2 (cid:17) + I sph tr (cid:16) ˆ σ sph2 (cid:17)(cid:27) , (12)for which the specific internal energy part e total − e kinetic fulfils the balance (cid:36) ∂ ( e total − e kinetic ) ∂ t = tr (cid:18) σ ∂ ε ∂ t (cid:19) (13) of 29 and specific entropy the balance (cid:36) ∂ s ∂ t = π s , (14)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 > π s ≥ e total that couples T and ε , and—correspondingly—there is no ε dependent term in s , due to neglectedthermal expansion.)Superficially, it seems redundant to also provide π s in the equivalent form (11). However, it is just thisnot-automatically-positive-definite form that will prove to be beneficial in the diagnostics of the numericalsolution.From either balance (13) or (14), the time derivative of temperature can also be expressed: ∂ T ∂ t = T (cid:36) c σ π s . (17)As a simple analysis of the PTZ model, for ‘slow’ processes, which is to be understood with respect tothe time scales τ dev , ˆ τ dev = ˆ E dev / E dev , τ sph , ˆ τ sph = ˆ E sph / E sph , (18)a rule-of-thumb approximation is to neglect the time derivative terms (to only keep the lowest timederivative term for each quantity) in (5). The result is nothing but the Hooke model (4), for which thelongitudinal and transversal wave propagation speeds are c longitudinal = (cid:115) E dev + E sph (cid:36) , c transversal = (cid:115) E dev (cid:36) . (19)Now, as opposed to this ‘static’ limit, let us consider the limit of ‘fast’ processes: then, it is the timederivative terms (the highest time derivative term for each quantity) that we keep. The result is the timederivative 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 , (20)where the inequalities follow from (15). Accordingly, the wave propagation speedsˆ c longitudinal = (cid:115) E dev ∞ + ˆ E sph ∞ (cid:36) > c longitudinal , ˆ c transversal = (cid:115) ˆ E dev ∞ (cid:36) > c transversal (21) of 29 follow. This, on one side, illustrates how the PTZ model can interpret that the dynamic elasticity coefficientsof rocks are larger than their static counterpart [1–4]. On the other side, the nontrivial—frequencydependent, therefore, dispersive—wave propagation indicates that the numerical solution of PTZ wavepropagation problems should contain the minimal possible amount of dispersion error, in order to giveaccount of the dispersive property of the continuum model itself. In parallel, the dissipative nature of thePTZ model requires the minimal possible amount of dissipation error to reliably describe the decrease ofwave amplitudes.
3. The Numerical Scheme
We take a Cartesian grid with spacings ∆ x , ∆ y , ∆ z , and time step ∆ t . Corresponding to the continuumformula (2), we introduce the finite difference discretization (cid:36) ( v x ) j + / l + / , m , n − ( v x ) j − / l + / , m , n ∆ t = ( σ xx ) jl + m , n − ( σ xx ) jl , m , n ∆ x + ( σ xy ) jl + / , m + / , n − ( σ xy ) jl + / , m − / , n ∆ y + ( σ xz ) jl + / , m , n + / − ( σ xz ) jl + / , m , n − / ∆ z , (22) (cid:36) ( v y ) j + / l , m + / , n − ( v y ) j − / l , m + / , n ∆ t = ( σ yx ) jl + / , m + / , n − ( σ yx ) jl − / , m + / , n ∆ x + ( σ yy ) jl , m + n − ( σ yy ) jl , m , n ∆ y + ( σ yz ) jl , m + / , n + / − ( σ yz ) jl , m + / , n − / ∆ z , (23) (cid:36) ( v z ) j + / l , m , n + / − ( v z ) j − / l , m , n + / ∆ t = ( σ zx ) jl + / , m , n + / − ( σ zx ) jl − / , m , n + / ∆ x + ( σ zy ) jl , m + / , n + / − ( σ zy ) jl , m − / , n + / ∆ y + ( σ zz ) jl , m , n + − ( σ zz ) jl , m , n ∆ z , (24)where the time index j refers to a value at t j = j · ∆ t , j + / to a value at t j + / = ( j + / ) · ∆ t , the spaceindex 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 shiftedin time by half; diagonal stress (and strain) components reside at integer positions, off-diagonal ones areshifted in the two directions that match with the two Cartesian indices; and, velocity components are onlyshifted in the direction matching with their Cartesian index (see Figure 2). From these formulae, the j + / indexed velocities can be expressed explicitly (as functions of earlier quantities). Figure 2.
Spatial arrangement of the discretized quantities (two-dimensional projection). The circles standfor diagonal tensor components, squares for offdiagonal ones, and triangles for vector components, differentcomponents with differently oriented triangles. Void quantities are prescribed by boundary condition (inthe case stress boundary conditions are considered, like here.) of 29
This same pattern—distribution of quantities—is used for the discretization of (2): ( ε xx ) j + l , m , n − ( ε xx ) jl , m , n ∆ t = ( v x ) j + / l + / , m , n − ( v x ) j + / l − / , m , n ∆ x , (25) (cid:0) ε yy (cid:1) j + l , m , n − (cid:0) ε yy (cid:1) jl , m , n ∆ t = (cid:0) v y (cid:1) j + / l , m + / , n − (cid:0) v y (cid:1) j + / l , m − / , n ∆ x , (26) ( ε zz ) j + l , m , n − ( ε zz ) jl , m , n ∆ t = ( v z ) j + / l , m , n + / − ( v z ) j + / l , m , n − / ∆ x , (27) ( ε xy ) j + l + / , m + / , n − ( ε xy ) jl + / , m + / , n ∆ t = (cid:40) ( v x ) j + / l + / , m + n − ( v x ) j + / l + / , m , n ∆ y + ( v y ) j + / l + m + / , n − ( v y ) j + / l , m + / , n ∆ x (cid:41) , (28) ( ε xz ) j + l + / , m , n + / − ( ε xz ) jl + / , m , n + / ∆ t = (cid:40) ( v x ) j + / l + / , m , n + − ( v x ) j + / l + / , m , n ∆ z + ( v z ) j + / l + m , n + / − ( v z ) j + / l , m , n + / ∆ x (cid:41) , (29) ( ε yz ) j + l , m + / , n + / − ( ε yz ) jl , m + / , n + / ∆ t = (cid:40) ( v y ) j + / l , m + / , n + − ( v y ) j + / l , m + / , n ∆ z + ( v z ) j + / l , m + n + / − ( v z ) j + / l , m , n + / ∆ y (cid:41) , (30)from which formulae the j + α (cid:16) σ dev pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) +( − α ) (cid:16) σ dev pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) + ( σ dev pq ) j + l (cid:48) , m (cid:48) , n (cid:48) − ( σ dev pq ) jl (cid:48) , m (cid:48) , n (cid:48) ∆ t = E dev (cid:20) α (cid:16) ε dev pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) + ( − α ) (cid:16) ε dev pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) (cid:21) + ˆ E dev ( ε dev pq ) j + l (cid:48) , m (cid:48) , n (cid:48) − ( ε dev pq ) jl (cid:48) , m (cid:48) , n (cid:48) ∆ t , (31) α (cid:16) σ sph pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) +( − α ) (cid:16) σ sph pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) + (cid:16) σ sph pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) − (cid:16) σ sph pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) ∆ t = E sph (cid:20) α (cid:16) ε sph pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) + ( − α ) (cid:16) ε sph pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) (cid:21) + ˆ E sph (cid:16) ε sph pq (cid:17) j + l (cid:48) , m (cid:48) , n (cid:48) − (cid:16) ε sph pq (cid:17) jl (cid:48) , m (cid:48) , n (cid:48) ∆ t , p , q = x , y , z , l (cid:48) , m (cid:48) , n (cid:48) = integers or half-integers depending on p , q , (32)where α = j + of 29 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 actuallythe 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 eachother, so second-order accuracy is achieved, while the conventional interpretation of the symplectic Eulermethod is first-order only. In parallel, since mechanical energy (the Hamiltonian) is a velocity dependentterm 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 programcode.For the aspects of thermodynamics, we also discretize (17) [here using the form (12)], explicitlyexpressing the j + / indexed temperature values from T j + / l , m , n − T j − / l , m , n ∆ t = (cid:36) c σ (cid:26) I dev tr (cid:16) ˆ σ dev2 (cid:17) jl , m , n + I sph tr (cid:16) ˆ σ sph2 (cid:17) jl , m , n (cid:27) , (33)where the notation (6) is utilized, and the traces are to be expanded in Cartesian components and theterms containing offdiagonal components—that reside at half space-shifted locations in two indices—areaveraged 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, beingthus treated as a plane-strain problem. Initially, the beam is in relaxed/equilibrium state (zero stress, strainand velocity, and homogeneous temperature). Subsequently, on one of its sides, a single normal stresspulse is applied, with profile σ yy ( t , x , 0, z ) = σ b (cid:110) (cid:104) − cos (cid:16) π t τ b (cid:17)(cid:105) · (cid:104) − cos (cid:16) π x − X /2 W b (cid:17)(cid:105)(cid:111) if 0 ≤ t ≤ τ b and X − W b ≤ x ≤ X + W b ,0 otherwise (34)(see Figure 3), where τ b is the duration of the pulse, σ b is its amplitude, W b is its spatial width, and X is thewidth of the beam. On the other sides, normal stress is constantly zero (free surfaces). Figure 3.
The spatial distribution of normal stress as the boundary condition on one side of a squarecross-sectioned beam that is infinitely long in the z direction. The excitation is a single cosine-shaped ‘bump’in time, too. A 50 ×
50 grid is considered in the x – y plane, for 240 time steps, where the time step is the largest atwhich stability is maintained. Notably, stability investigation is fairly involved for this problem and itrequires a separate whole study. Further settings (in appropriate units) are (cid:36) = E dev = E sph = τ b = σ b = W b / X = Figures 4 and 5 show snapshots of a stress component distribution and a velocity component.
Figure 4.
Distribution of a stress component (left column) and of a velocity component (right column) atvarious instants, in the Hooke case. From top to bottom: snapshots at instants ( ) τ b , τ b , ( ) τ b , 2 τ b ,respectively. 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 timedirection (for kinetic energy).
Figure 5.
Continuation of Figure 4: distribution of a stress component (left column) and of a velocitycomponent (right column) at various instants, in the Hooke case. From top to bottom: snapshots at instants ( ) τ b , 3 τ b , ( ) τ b , 4 τ b , respectively. Figure 6.
Mechanical energy types as functions of time, for the Hooke case.
In a PTZ medium, the solution of the analogous problem is similarly good. Figures 7–10present snapshots, where Figures 9–10 display two further quantities: (cid:114) tr (cid:16) ˆ σ dev2 (cid:17) (essentially theHuber–Mises–Hencky or von Mises equivalent stress) and temperature. Dissipation is nicely indicatedvia temperature. Settings (in addition to the above ones for the Hooke case) are ˆ E dev = τ dev = c σ = T = α = / . Figure 7.
Distribution of a stress component (left column) and of a velocity component (right column) atvarious instants, in the Poynting–Thomson–Zener (PTZ) case. From top to bottom: snapshots at instants ( ) τ b , τ b , ( ) τ b , 2 τ b , respectively. Figure 8.
Continuation of Figure 7: distribution of a stress component (left column) and of a velocitycomponent (right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( ) τ b , 3 τ b , ( ) τ b , 4 τ b , respectively. Figure 9.
Cont.
Figure 9.
Snapshots of the distribution of the stress invariant (cid:113) tr ( ˆ σ dev ) (left column) and of temperature(right column) at various instants, in the PTZ case. From top to bottom: snapshots at instants ( ) τ b , τ b , ( ) τ b , 2 τ b , respectively. Figure 10.
Cont.
Figure 10.
Continuation of Figure 9: Snapshots of the distribution of the stress invariant (cid:113) tr ( ˆ σ dev ) (leftcolumn) and of temperature (right column) at various instants, in the PTZ case. From top to bottom:snapshots at instants ( ) τ b , 3 τ b , ( ) τ b , 4 τ b , respectively. The diagnostic role of the various energies, and especially their sum, is a great help again for checkingwhether the simulation performs acceptably. Figure 11 illustrates how the scheme that is introduced abovebehaves in this respect.
Figure 11.
Total energy and the various energy types as functions of time, for the PTZ case.
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, andthen one of its sides being pressed by a cosine ‘bump’ in time, as well as in both spatial directions: σ yy ( t , x , 0, z ) = σ b (cid:110) (cid:104) − cos (cid:16) π t τ b (cid:17)(cid:105) · (cid:104) − cos (cid:16) π x − X /2 W b (cid:17)(cid:105) · (cid:104) − cos (cid:16) π z − X /2 W b (cid:17)(cid:105)(cid:111) if 0 ≤ t ≤ τ b and X − W b ≤ x , z ≤ X + W b ,0 otherwise (35)(see Figure 12), where the notations are analogous to the ones of the previous case: τ b is the duration of thepulse, σ b is its amplitude, W b is its spatial width, and X is the length of the edges of the cube. On the othersides, normal stress is constantly zero (free surfaces), like in the previous example. Figure 12.
The spatial distribution of normal stress as the boundary condition on one side of a cube. Theexcitation is also a single cosine-shaped ‘bump’ in time.
Here, we only present the PTZ model, on a 25 × ×
25 grid, with all other settings being the same asin the previous example.The solution proves similarly satisfatory, as in the former case. Figures 13 and 14 show snapshots ofthe distribution of the von Mises stress invariant on two mid-planes of the cubes.
Figure 13.
Cont.
Figure 13.
Snapshots of the distribution of the stress invariant (cid:113) tr ( ˆ σ dev ) at the mid-plane parallel to theexcited side (left column) and at a mid-plane orthogonal to it (right column) at various instants, in the PTZcase. From top to bottom: snapshots at instants ( ) τ b , τ b , ( ) τ b , 2 τ b , respectively. Figure 14.
Cont.
Figure 14.
Continuation of Figure 13: snapshots of the distribution of the stress invariant (cid:113) tr ( ˆ σ dev ) atthe 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 ( ) τ b , 3 τ b , ( ) τ b , 8 τ b ,respectively. Next, in this example, we demonstrate the usefulness of another thermodynamical quantity, entropyproduction 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 = T I ˆ σ · (cid:18) ˆ E ∂ ε ∂ t − τ ∂ σ ∂ t (cid:19) (36)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 jn I ˆ σ jn · (cid:32) ˆ E ε j + n − ε jn ∆ t − τ σ j + n − σ jn ∆ t (cid:33) , (37)1 T jn I ˆ σ jn · (cid:32) ˆ E ε j + n − ε j − n ∆ t − τ σ j + n − σ j − n ∆ t (cid:33) , (38)1 T jn I ˆ σ j + n · (cid:32) ˆ E ε j + n − ε jn ∆ t − τ σ j + n − σ jn ∆ t (cid:33) , (39)1 T jn I ˆ σ j + n · (cid:32) ˆ E ε j + n − ε j − n ∆ t − τ σ j + n − σ j − n ∆ t (cid:33) , (40)where T jn denotes the time average (cid:16) T j − / n + T j + / n (cid:17) /2.These four versions are integrated in space and plotted in Figure 15, the left column for a stablesetting and the right column for an unstable one. The energies are also displayed. Only 25 space cells havebeen chosen to enhance artefacts. Figure 15.
Entropy production rate according to the four discretized versions (37)–(40), respectively, for atime step resulting in a stable run (left column) {namely, when the rheological Courant number ˆ c ∆ t / ∆ x is 1, where ˆ c = (cid:113) ˆ E / ( τ(cid:36) ) is the high-frequency rheological wave propagation speed} and for a choiceleading to an unstable outcome (right column) {rheological Courant number ˆ c ∆ t / ∆ x = τ E / ˆ E = X is such that ˆ E / E = X / c with the low-frequency/elastic wavepropagation speed c = (cid:112) E / (cid:36) . See [8] for further details. Visibly, certain versions become negative when instability gets exposed. Moreover, some becomenegative 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 thepressed cube: the outcomes can be seen in Figure 16. An intentionally rude 10 × ×
10 grid is taken, thetime unit is divided to 125 time steps in order to lead to a stable solution and to 100 time steps producingan unstable one. Further settings are (cid:36) = E dev = E sph =
5, ˆ E dev = τ dev = c σ = T = τ b = σ b = W b / X = α = / . Figure 16.
Entropy production rate according to the four discretized versions (37)–(40) generalized to threespace dimensions, respectively, computed for a PTZ cube, with a time step resulting in a stable run (leftcolumn) and for a choice leading to an unstable outcome (right column).
One can see that an appropriately chosen discretization diagnoses instability.
Such numerical comparisons between stable and unstable parameter domains, which are probablyenhanced by some analytical considerations, can help in the future in deciding which discretized entropyproduction 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 finiteelement solutions was presented. Repeating it for the present three space dimensional scheme wouldbe informative. We found it advisable to start investigating the higher dimensional wave propagationdescribing the possibilities of finite element softwares in a much simpler setting than the one treatedabove, based on the (discouraging) experience that is gained with the one-dimensional case, and since arheological model like PTZ in the deviatoric–spherical formulation is difficult to translate to the capabilitiesof classic finite element softwares.Namely, we consider the wave equation in two spatial dimensions, for a single scalar degree offreedom u = u ( t , x , y ) , with constant (unit) coefficients and no source term in the equation: ∂ u ∂ t = ∂ u ∂ x + ∂ u ∂ y . (41)A rectangular sample is taken with ‘flux’ (Neumann) boundary condition (normal spatial derivativecomponent is prescribed): on one of the edges, one cosine-type pulse of the kind (34) is applied, while theother edges are kept ‘free’ (zero normal derivative). Initially, both u and ∂ u / ∂ t are set to zero (‘relaxedinitial 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 degreeof freedom after two units of time, obtained with various available settings of COMSOL v5.3a.
Figure 17.
Cont.
Figure 17. u ( x , y ) with various COMSOL settings (cf. Figure 1). Left column, from top to bottom: BDFMax. order 5 & min. order 2, BDF Max. order 3 & min. order 2, BDF Max. order 1 & min. order 1, BDFMax. order 5 & min. order 3 Strict time stepping, RK34 Default. Right column, from top to bottom: DP5PI-Smooth, DP5 PI-Standard, DP5 PI-Quick, Generalized- α (GA), GA Predictor Constant. Not shown (failedlike 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). One can see that the solution depends on the settings very seriously. There are large-scale differencesand fine-structured irregularities. Moreover, one has no means of validation which outcome is correct towhat extent.When taking into account that our full problem has 16 coupled degrees of freedom in three spatialdimensions and with further time derivatives (in the PTZ model), it does not seem to be reasonable to tryto represent it via any commercial finite element software as long as such a much simpler problem cannotbe confidently treated.
7. Discussion
The numerical scheme presented here, due to its symplectic root, second-order accuracy, and theequation-friendly and spacetime geometry friendly arrangement of discretized quantities, has been foundto provide reliable results in a fast and resource-friendly way. Being a finite difference scheme, it is notvery flexible to simulate arbitrary shaped samples, but, already, the extension of the Cartesian formulaeto cylindrical and spherical geometries promises useful applications, including the various wave-basedmeasurement methods that are used in rock mechanics (see, e.g., [26]), many of which rely on simple andeasily treatable sample shapes. Fitting a rheological model on experimental data may require many runsso good finite difference schemes find their applicability.The investigation of stability and dissipative and dispersion error is expected to be much moreinvolved 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 ofthe 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 buildingconnection to a finite element—spacetime finite element—approach, along which way objects of generalshape could also be treated. Notably, the current finite element paradigm has deficiencies and probablyneeds 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 controllingnumerical 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 beutilized.
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 NationalResearch, 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 theMinistry for Innovation and Technology. This paper was supported by the János Bolyai Research Scholarship of theHungarian 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átramountain range.
Class. Quantum Gravity , , 114001, doi:10.1088/1361-6382/aa69e3.2. Ván, P.; Barnaföldi, G.G.; Bulik, T.; Biró, T.; Czellár, S.; Cie´slar, 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. , ,1693–1734, doi:10.1140/epjst/e2019-900153-1.3. Davarpanah, S.M.; Ván P.; Vásárhelyi, B. Investigation of relationship between dynamic and static deformationmoduli of rocks. Geomech. Geophys. Geo-Energy Geo-Resour. , , 29, doi:10.1007/s40948-020-00155-z.
4. Fülöp, T. Wave propagation in rocks—Investigating the effect of rheology.
Period. Polytech. Civ. Eng. ,to appear, doi:10.3311/PPci.16096.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; Article number:ISRM-EUROCK-2020-017.6. Fülöp, T.; Szücs, M. Analytical solution method for rheological problems of solids. arXiv , arXiv:1810.06350.7. Rieth, Á.; Kovács R.; Fülöp, T. Implicit numerical schemes for generalized heat conduction equations.
Int. J. HeatMass Transf. , , 1177–1182, doi:10.1016/j.ijheatmasstransfer.2018.06.067.8. Fülöp, T.; Kovács R.; Szücs, M.; Fawaier, M. Thermodynamical extension of a symplectic numerical schemewith half space and time shifts demonstrated on rheological waves in solids. Entropy , , 155,doi:10.3390/e22020155.9. Asszonyi, Cs.; Fülöp, T.; Ván, P. Distinguished rheological models for solids in the framework of athermodynamical internal variable theory. Contin. Mech. Thermodyn. , , 971–986.10. Hairer, E.; Lubich C.; Wanner, G. Geometric Numerical Integration , 2nd ed.; Springer: Berlin/Heidelberg, Germany,2006.11. Zinner, C.P.; Öttinger, H.C. Numerical stability with help from entropy: Solving a set of 13 moment equations forshock tube problem.
J. Non-Equilib. Thermodyn. , , 43–69, doi:10.1515/jnet-2018-0038.12. Shang, X.; Öttinger, H.C. Structure-preserving integrators for dissipative systems based on reversible-irreversiblesplitting. arXiv , arXiv:1804.05114.13. Portillo, D.; García Orden, J.C.; Romero, I. Energy-Entropy-Momentum integration schemes for general discretenon-smooth dissipative problems in thermomechanics. Int. J. Numer. Methods Eng. , , 776–802, doi:10.1002/nme.5532.14. Vermeeren, M.; Bravetti, A.; Seri, M. Contact variational integrators. J. Phys. A Math. Theor. , , 445206.15. Gay-Balmaz, F.; Yoshimura, H. Variational discretization of the nonequilibrium thermodynamics of simplesystems. Nonlinearity , , 1673.16. Couéraud, B.; Gay-Balmaz, F. Variational discretization of thermodynamical simple systems on Lie groups. Discret. Cont. Dyn. Syst. S , , doi: 10.3934/dcdss.2020064.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. , , 1841–1858.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. , , 2235–2248.19. Berezovski, A.; Ván, P. Internal Variables in Thermoelasticity ; Springer: Cham, Switzerland, 2017.20. Janeˇcka, A.; Málek, J.; Pr ˚uša, V.; Tierra, G. Numerical scheme for simulation of transient flows of non-Newtonianfluids characterised by a non-monotone relation between the symmetric part of the velocity gradient and theCauchy stress tensor.
Acta Mech. , , 729–747.21. Fülöp, T.; Ván, P. Kinematic quantities of finite elastic and plastic deformation. Math. Methods Appl. Sci. , ,1825–1841.22. Ván, P. Galilean relativistic fluid mechanics. Contin. Mech. Thermodyn. , , 585–610,doi:10.1007/s00161-016-0545-7.23. Matolcsi T. Ordinary Thermodynamics: Nonequilibrium Homogeneous Processes ; Akadémiai Kiadó (Publishing Houseof the Hungarian Academy of Sciences): Budapest, Hungary, 2004.24. Asszonyi, Cs.; Csatár, A.; Fülöp, T. Elastic, thermal expansion, plastic and rheological processes—Theory andexperiment.
Period. Civ. Eng. , , 591–601, doi:10.3311/PPci.8628.25. Szücs, M.; Fülöp, T. Kluitenberg–Verhás rheology of solids in the GENERIC framework. J. Non-Equilib. Thermodyn. , , 247–259, doi:10.1515/jnet-2018-0074.26. Malhotra, V.M.; Carino, N.J. (Eds.) Handbook on Nondestructive Testing of Concrete , 2nd ed.; CRC Press: Boca Raton,FL, USA, 2003; doi:10.1201/9781420040050.
27. Szabó, B.; Kossa, A. Characterization of impacts of elastic-plastic spheres.
Period. Polytech. Mech. Eng. ,64