Thermomass Theory in the Framework of GENERIC

Thermomass theory was developed to deal with the non-Fourier heat conduction phenomena involving the influence of heat inertia. However, its structure, derived from an analogy to fluid mechanics, requires further mathematical verification. In this paper, General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) framework, which is a geometrical and mathematical structure in nonequilibrium thermodynamics, was employed to verify the thermomass theory. At first, the thermomass theory was introduced briefly; then, the GENERIC framework was applied in the thermomass gas system with state variables, thermomass gas density ρh and thermomass momentum mh, and the time evolution equations obtained from GENERIC framework were compared with those in thermomass theory. It was demonstrated that the equations generated by GENERIC theory were the same as the continuity and momentum equations in thermomass theory with proper potentials and eta-function. Thermomass theory gives a physical interpretation to the GENERIC theory in non-Fourier heat conduction phenomena. By combining these two theories, it was found that the Hamiltonian energy in reversible process and the dissipation potential in irreversible process could be unified into one formulation, i.e., the thermomass energy. Furthermore, via the framework of GENERIC, thermomass theory could be extended to involve more state variables, such as internal source term and distortion matrix term. Numerical simulations investigated the influences of the convective term and distortion matrix term in the equations. It was found that the convective term changed the shape of thermal energy distribution and enhanced the spreading behaviors of thermal energy. The distortion matrix implies the elasticity and viscosity of the thermomass gas.

In order to characterize non-Fourier heat conduction, various models had been proposed based on different assumptions and physical pictures. Hyperbolic heat conduction equations were preferred because they solved the paradox of infinite thermal disturbance speed [15] and, in recent studies, investigated the phonon momentum conservation process [16,17], which implied the inertial phenomena of heat flux [18,19]. Cattaneo-Vernotte equation [20,21], the first mathematical hyperbolic heat conduction equation, assumed a relaxation process between temperature gradient and heat flux, giving a telegraph equation when coupled with the energy conservation equation. However, This paper aims to apply the thermomass theory in the framework of GENERIC, demonstrating the mathematical reasonability of the governing equations of thermomass gas and the physical foundation of the heat conduction in the framework of GENERIC.

Thermomass Theory
The part will introduce the thermomass theory briefly. Heat is conserved in heat conduction, which indicates the mass nature of heat. Therefore, it is also called thermomass. The transport process of heat can be regarded as the flowing behaviors of thermomass gas in porous media. Thus, in this way, the heat conduction phenomena can be described by fluid mechanics. Guo et al. [36] noticed these characteristics and proposed the thermomass theory, which focused on the extra mass increased due to thermal vibration. The theory aims to develop the dynamics of thermomass gas in heat conduction. In fact, Grmela [55] shared the same idea and regarded the hyperbolic heat conduction in fluid as the transport of two component fluid in his paper in 1983. The first component was fluid particle and the second was called caloric component.
The relativistic mass of an object is defined based on the Einstein's mass-energy relation: where E is the mechanical energy and c is light speed. Similar to the mechanical energy, the relativistic mass of thermal energy in a dielectric solid is: which is called thermomass, as extra mass added to the rest mass. As for heat conduction in the dielectric solid, thermomass is transported, while the molecules stay in the original places. The dynamics of thermomass gas is based on two state variables, ρ h and u h , representing the local thermomass gas density and drift velocity of thermomass gas in per unit volume, respectively, which are defined as: where ρ is the material density and ρ h is the density of thermomass gas. C V is the specific heat. T and q are the local temperature and heat flux, respectively. Then, three conservation laws are established, namely the thermomass continuity equation, the thermomass momentum equation and the thermomass energy equation: where u h u h in Equation (7) is a second order tensor and u h ·u h is a scalar. The first equation corresponds to the heat conservation equation, which is: Entropy 2020, 22, 227 4 of 24 In the thermomass momentum Equation (7), the first term represents the temporal inertial force of thermomass while the second term represents the spatial inertial force. P h is the thermomass gas pressure defined as a quadratic function [33]: and f h is the resistance, usually assumed to have a linear relationship with the drift velocity of thermomass gas: The coefficient µ can be derived from the comparisons with Fourier's law. Of course, different resistance models lead to different time evolution equations. For example, Dong et al. [39] adopted Brinkman's equation for the resistance force to describe the boundary influences on thermomass propagation in microscale medium. When the temporal inertial force and spatial inertial force can be ignored, the thermal pressure is equal to the thermal resistance [36], namely: At that time, Equation (12) is actually Fourier's law, so that µ can be expressed in the form of thermal conductivity k as: Then, the thermomass momentum equation (Equation (7)) becomes: where qq is a second order tenser. The time evolution equation of heat flux q is derived with temperature T on the foundation of thermomass gas dynamics. This equation can reduce to the conventional Fourier's conduction law when all the inertia effects are omitted. It reduces to the hyperbolic heat conduction equation when the spatial inertia is neglected and it is also able to describe the non-Fourier heat conduction phenomena even in the steady state, when the temporal inertia is omitted. Therefore, Equation (14) is called the general heat conduction law in the references [35,36]. The conservation equation of thermomass energy can be derived from the coupling of thermomass continuity equation and momentum equation, while new quantities are defined as the thermomass kinetic energy (Equation (15)) and thermomass potential energy (Equation (16)): In the transport process where the drift velocity of thermomass gas is not so low that the kinetic energy of thermomass cannot be neglected, the kinetic energy and potential energy of thermomass convert into each other mutually cycle by cycle, and the system goes through the oscillations of heat flux and temperature. It is important to emphasize that there are two systems in a solid. The first one is the porous skeleton of rest mass of molecules and the other is the thermomass gas, which is full of the porous skeleton and can flow through it, if a temperature difference is applied to the solid. Heat conduction in a solid is in fact the flowing of thermomass gas in a porous medium, which can be described by the principles of fluid dynamics. The quantities in the porous skeleton system are those we are familiar with, such as density ρ and specific heat C V . In the thermomass gas system, the quantities are labeled by the subscript h, such as ρ h and u h . ρ h is the density of thermomass gas, which is the function of the density of rest mass skeleton, ρ, as shown in Equation (4). It is assumed that the propagation of thermomass gas is only influenced by the thermomass gas system, free from the skeleton system.

Thermomass Theory in the Framework of GENERIC
GENERIC theory divides the whole process into the reversible part and the irreversible part, controlled by two potentials, namely the Hamiltonian potential and the dissipation potential. Thermomass gas can be regarded as a kind of special fluid, and thus its transport process can be dealt with in the framework of GENERIC with the analogy to fluid mechanics. Because different state variables take effect in different cases, choosing proper state variables is significant, which makes the GENERIC theory easier to be extended. In the thermomass system, there are two elementary state variables to choose, thermomass gas density ρ h and thermomass gas flow m h , which represents the thermomass gas momentum, corresponding to the thermomass gas density ρ h and thermomass gas drift velocity u h in thermomass theory by: In this paper, the thermomass system based on the solid system is studied. Then, the Poisson bracket are derived via Hamiltonian mechanics [51]: where Einstein's summation convention is used for the subscript i, j and m h,i represents the i th component of the vector m h . Equation (18) is obtained from the analogy to the Poisson bracket in fluid mechanics [51]. In this equation, A is a function with respect to thermomass density and thermomass momentum A(ρ h , m h ) and represents some system properties. H is the Hamiltonian energy in thermomass system. Poisson bracket depicts the reversible process, including the propagation and the transition between thermomass kinetic energy and thermomass potential energy. The Poisson bracket is supposed to satisfy three basic principles, namely the linearity, These principles are testified for Equation (18) since its structure is similar to that in fluid mechanics. Then, the following equations are obtained via integration by parts: .
The time evolution equations of the reversible process can be derived by substituting Equation (19) into Equation (21): A so-called pressure potential is introduced to simplify the expression, as: The construction of such pressure potential is common in Hamiltonian mechanics and it implies the transition between kinetic energy and potential energy. Equation (23) is transformed into: by substituting Equation (24) into Equation (23). These Equations (22) and (23) describe a non-dissipative transport process of thermomass and thermomass momentum. They are supposed to correspond with the following hyperbolic heat conduction equations without dissipation: It has to be clarified that the non-dissipative process described by Equations (26) and (27) does not exist in reality. They are analyzed only as the reversible part of the whole heat transport process and what occurs in the process is the combination of the reversible part and the irreversible part. Through the comparisons between Equations (26) and (27) and Equations (22) and (23), it is found that a proper definition of Hamiltonian potential H makes these two equation sets become the same. A specific form of H can be chosen as: which is the thermomass energy in the thermomass theory, consisting of kinetic energy and potential energy. Moreover, ε is a constant with dimension m 5 kg −1 s −2 to make sure that the Hamiltonian energy H has the dimension J/m 3 , which actually means the energy density in the unit volume. Taking into consideration the relationships Equations (4), (5), (28), Equation (22) and Equation (23) are transformed into: and the parameter ε is defined as: When the variation of temperature is not large, linear approximation can be adopted and the variation of energy density e in the expression (31) is omitted. In this way, the coefficient ε 2 e/c 2 in the right of the Equation (30) gives a linear energy-dependent relationship. Of course, the formulations of Hamiltonian energy are not unique, a linear or a cubic relationship for thermomass potential energy fits the equations as well. However, it is demonstrated that the linear relationship: leads to a constant pressure potential: which is not suitable in the physical system. The cubic relation: requires that the thermal conductivity relies on the quadratic energy density e 2 , which is not considered in this paper: As for the irreversible processes, eta-function and dissipation potential are introduced and the time evolution equations are in entropy representation. Eta-function η, whose derivative with respect to energy density e(x,t) in local equilibrium state equals the inverse temperature: has the meaning of entropy. The subscript eq represents the equilibrium state. η e , abbreviation for the derivative of η with respect to e is called the conjugate variable of e. In the GENERIC theory, contact geometry and conjugate variables are widely adopted. In equilibrium state, the eta-function has the expression: in solid physics with constant D. According to the maximum entropy principle, the eta-function approaches its maximum during the process from nonequilibrium state to equilibrium state. η(ρ h , m h ) is the function of thermomass density and thermomass momentum and its variation is determined by these two state variables: The introduction of dissipation potential in GENERIC theory extends the irreversible process from linear gradient dynamics to nonlinear cases. The time evolution can be written in the form of inner product: and when A is the eta-function itself, it is supposed to satisfy that: Keeping that in mind, the derivative of A with respect to time could be seen as the sum of the time evolution of state variables (Equation (38)), Equation (39) could be transformed into: As for how to get a proper dissipation potential, there are many choices. It seems that the simplest form is the quadratic one: One thing to be noticed is that the dissipation potential might be not only the function of state variables, but also the gradient of state variables and so on. Therefore, the specific form of dissipation potential is pertinent to the material types, microstructures and other characteristic properties. Here, one typical form can be chosen as: The absolute value of η 0 is difficult to calculate in practice as it is related with the entropy in equilibrium state. However, luckily, only the derivatives of eta-function with respect to the state variables are necessary in the equation set. The dissipation potential is free from the conjugate variable η ρ h , because the thermomass is conserved, rather than being reduced during the heat conduction process. The time evolution can be derived by combining the reversible process and irreversible process: In the nonequilibrium heat conduction process, the thermomass is transferred and conserved while the thermomass momentum is reduced until a steady state is established.

Comparisons Between GENERIC Theory and Thermomass Theory
In the thermomass theory, the momentum equation is in the form of thermomass density and drift velocity of thermomass gas: When the thermomass theory is applied in the framework of GENERIC, the counterpart equation (Equation (47)) is expressed in the form of thermomass density and thermomass momentum. In order to obtain more obvious comparisons, the same variables are substituted in these equations. They are the energy density e(x,t) and heat flux q(x,t), which are more common in hyperbolic heat conduction equations. Substituting the relationships Equations (4) and (5) into Equations (46) and (47): are derived. The first equation (Equation (49)) is the heat conservation equation. The second equation (Equation (50)) is compared with the momentum equation in thermomass theory (Equation (7)), which is also rewritten in the form of e and q: Entropy 2020, 22, 227 9 of 24 The first term in the left of Equation (51) represents the temporal inertial effect of heat and the second term is the convective term, representing the spatial inertial effect. Hyperbolic heat conduction phenomena arise because of the temporal inertial force of thermomass. Then, the spatial inertial force influences the interactions between the heat fluxes in different directions. The third term in Equation (51) is the pressure gradient, characterized by the gradient of energy density. The last term is the resistance, leading to the decrease of heat flux. In the whole process, the heat is conserved and transferred, while the thermomass momentum is reduced.
By comparing Equation (50) and Equation (51), it is found that they have the same structure and similar terms. The difference is in the definition of the coefficients. In the thermomass momentum equation, the pressure term is derived from the thermal pressure in solid. The resistance is obtained from the analogy to the porous media, and it is assumed that the resistance force has a linear relation with the drift velocity of thermomass gas [33,36]. In the GENERIC equation, the third term, the pressure potential, is obtained from the reversible process, representing the combination of thermomass kinetic energy and thermomass potential energy. Furthermore, the fourth term comes from the dissipation in the irreversible process, leading to entropy production and thermomass energy dissipation [51].

Thermomass Energy
In the GENERIC theory, the reversible process and irreversible process are controlled by two individual potentials. In the reversible process, Hamiltonian potential energy H, which consists of thermomass kinetic energy and thermomass potential energy, takes effect, while in the irreversible process, dissipation potential governs and keeps damping monotonously. It was found in this paper that these two potentials have similar forms and they can be unified as one potential energy, namely the thermomass energy.
It was noticed that in the reversible process, thermomass energy is conserved, while in the irreversible process, resistance dissipates the thermomass energy. These behaviors are much like the controlling potentials in GENERIC theory. Therefore, the Hamiltonian potential and dissipation potential in the thermomass system with state variables ρ h and m h can be written as: The following relationships are required to be satisfied: where τ is the relaxation time. In this way, it can be answered what is dissipated in the process in heat conduction. It is the thermomass energy H. In the reversible process, thermomass energy is transmitted and conserved, with the transformation between thermomass kinetic energy and thermomass potential energy, while in the irreversible process, thermomass energy is dissipated monotonously. The behaviors of thermomass are much similar with those of entropy. However, the irreversible process leads to the increase of entropy production, while the thermomass energy is reduced due to its dissipation. The thermomass energy analysis has two advantages. Firstly, it provides a uniform formulation for the Hamiltonian potential and the dissipation potential. From nonequilibrium state to equilibrium state, the thermomass energy keeps being dissipated monotonously, which can be used as the Lyapunov function. Secondly, it gives a clearer and more intrinsic physical picture for the heat transport in GENERIC theory, where reversible process is described in energy representation while the irreversible one is in entropy representation.

Extension in the Framework of GENERIC
The reduction from microscale view to macroscale view ignores some details and focuses only on the variation of macroscopic patterns. Nevertheless, different levels of descriptions are required in different situations. For example, in an equilibrium system, only one state variable is necessary, namely the total thermal energy E. However, in a system described by Fourier's law, the distribution of local energy density e(x) is supposed to be considered. If more proper variables are added to the list, more details of the system are revealed. However, the word proper is significant, which requires to pick up the relatively more general variables and leave out the less significant ones and that the state variables are rational. For example, when both heat flux and the gradient of heat flux influence the heat conduction process, the heat flux is obviously more general, while the gradient of it might be only a good supplement.
One advantage of the GENERIC framework is that it is easy to extend because it provides the mathematical structure. Thus, flexible state variables can be chosen according to different systems. Of course, the most important thing is that these state variables are reasonable for the target system. One possible extension is the distortion matrix, in which case the propagation of heat suffers from viscosity, and heat fluxes in different directions have stronger interactions with each other.
At first, the concept of volume density of labels a is introduced [51,56]. In Lagranian mechanics, the coordinate is set to be with the thermomass particles and the particle trajectories are recorded. However, when these equations are transformed into Euler representation, the coordinates might be compressed or twisted. The volume density of labels a, which is a vector, implies where the particles come from. The gradient of a, called distortion matrix D, shows the torsion and deformation of the given volume, which has been applied in elasto-plastic solid successfully [56]. The distortion matrix D is usually used to describe the solid-like fluid or fluid-like solid, which has the flowing and elasticity characteristics at the same time. In heat conduction, even though the thermomass gas is regarded as a kind of fluid, the real transport process is complex and these kinds of fluids might have the solid characteristics. Therefore, it is reasonable to apply the theory for elasto-plastic solid to thermomass gas. Actually, the historical effect of heat flux with temperature [7] or the vibration interference [57], for example, can be seen as some kind of elasticity. It is thought that the distortion matrix is relevant with the initial state of thermal energy distribution and represents the viscous effect and elastic effect.
The time evolution equations with label field a are derived as the first step. The derivatives of label field a with respect to displacement and momentum are defined as: Thus, the Poisson bracket, consisting of thermomass density ρ h , thermomass momentum m h and label field a has the form of: Via integration by parts, the time evolution of these state variables yields: From the equations above, it is found that label field takes part in a pure convective process in reversible systems. The force potential π † differs from that in Equation (24)  .
However, the irreversible process might destroy the continuity of label field a and Equation (55) might not be always satisfied, which results in the failure of Jacobi identity. The last part in Equation (62) is added to the equation to make sure that it fulfills the Jacobi identity even for incompatible distortion fields. The time evolution equations derived are: Due to the complex distortion matrix, a universal form of force potential is difficult to find. A possible form of H is: In fact, Equation (65) implies that the heat flux is not only determined by the gradient along this direction, but also interacts with the gradients of variables along other directions. The physical meaning of the distortion matrix, the elasticity of thermomass, can be regarded as a nonlocal effect, because it takes into consideration the interactions between the thermomass particles and their surroundings. Its formulation is similar with that of Rational Extended Thermodynamics [47]. However, the closure problems are avoided and it closes the equation set by the application of dissipation potential, as shown in the following part.
The new variables also have their own part in the dissipation process, which is irreversible. When the phonon entropy is the function of label field a or distortion matrix field D, or even the gradient of the distortion matrix, which is a third-order tensor, the dissipation terms have different characteristics. Here, for simplicity, it is assumed that the eta-function is the function of the distortion matrix D. Then the time evolution Equation (64) is rewritten as: One possible example of Ξ and η is proposed. The dissipation potential is set to be a quadratic function, as: The irreversible process requires that the dissipation potential always keeps damping during the process of propagation, thus, the terms is supposed to give the equation the tendency to equilibrium. Analogous with the behavior of heat flux, the relationship between eta-function and distortion matrix field is set to be: In the above equations, β and ϕ are the coefficients related to the material thermal properties. The complete time evolution equation for D is then rewritten as: The state variables that are necessary to describe the system are changing in different problems and at different description levels. Taking more state variables does not always mean better, because while some patterns are negligible, more details can mean more meaningless efforts and disordered importance. In view of the current experimental measurement techniques, label field and distortion matrix field are both difficult to observe. However their effects are presented by the thermomass density ρ h and thermomass momentum m h . This article makes an effort to demonstrate that the heat conduction principles for thermomass system might be more complex than was thought and the theory can be developed to a more general stage.
Inner heat source absorbs and emits heat, providing a new term to the dissipation potential. Different from distortion matrix, inner heat source is not a state variable, but a given system parameter. In a system represented by the state variables ρ h and m h , the inner heat source has nothing to do with Hamiltonian mechanics and it influences the heat conduction process by changing the dissipation potential. The time evolution equations of thermomass system with inner heat source are supposed to be: according to the physical interpretation, where s h is the inner thermomass source, representing the mass density variation per unit volume per second. Therefore, to make new potential quantities compatible with Equations (71) and (72), the Hamiltonian potential and dissipation potential are changed into: The analyses of dissipation potential show that the thermomass source potential is written as the product of the derivative of eta-function with respect to the thermomass density ρ h and the inner thermomass source s h . If this inner source is a hot source s h > 0, it adds to the dissipation potential. Otherwise, it reduces the dissipation potential. Distortion matrix and inner heat sources are used as the examples to demonstrate how to extend the thermomass theory in the framework of GENERIC. The state variables are chosen according to the characteristics of the system. Inner heat sources gives a new term to the dissipation potential, and the distortion matrix allows for stronger interactions between heat fluxes in different directions, describing the viscous effect and elasticity of thermomass gas.

Influences of the Convective Term
One of the characteristics of Equation (50), the time evolution equation for the thermomass momentum, is the second term which implies the momentum convection effects. The influences of the convective term on the propagation behaviors are investigated in this section by numerical simulations. Here, the heat conduction problems are solved in one-dimensional regimes and two-dimensional regimes.
The transient heat conduction phenomenon in a large plate with finite length as shown in Figure 1, is studied, which can be regarded as a one-dimensional heat conduction problem. The left boundary is heated by a sine-shape heat pulse for the time from 0 to t pulse : where Y 1 is the amplitude of heat pulse. Before calculations, the state variables and coefficients are normalized and the numerical algorithm comes from computational fluid mechanics [58][59][60]. After normalization, Y 1 is chosen to be 1. More numerical details can be seen in the Appendix A. If the thermal conductivity is constant, the effects of the convective term can be highlighted. Thus, in Figure 2 the problem with constant thermal conductivity is analyzed. That is, the equation: is solved, where the second term denotes the convective effect. In the simulations, k and τ are both set to be 1 to make the calculations easier. The choices of these coefficients are important, of course. However, it does not influence the qualitative characteristics of the equations, as it has been adopted in the reference [59][60][61][62]. The propagation patterns of the energy density with and without the second term are displayed in Figure 2. The results for the equations without the second term are plotted in solid line while the dash line denotes those for the equations with convective term. It can be seen that the heat pulse at the left boundary produces a thermal wave in the medium and it keeps the original pattern during propagation while the peak value is dissipated as shown by the solid line. When the convective term is considered, the shape is changed and a sharp wave front comes into being. The convective term shows the drift behaviors of thermomass conveyed by momentum. Thus, it is an obvious deduction when the mass nature of heat is taken into consideration. Meanwhile, the existence of convective term speeds up the thermal wave and leads to higher dissipation rate for the peak value of the patterns. Entropy 2019, 21, x FOR PEER REVIEW 14 of 25 Figure 1. Schematic diagram of the numerical regime for one-dimensional heat conduction problem.
where Y1 is the amplitude of heat pulse. Before calculations, the state variables and coefficients are normalized and the numerical algorithm comes from computational fluid mechanics [58][59][60]. After normalization, Y1 is chosen to be 1. More numerical details can be seen in the Appendix A. If the thermal conductivity is constant, the effects of the convective term can be highlighted. Thus, in Figure 2 the problem with constant thermal conductivity is analyzed. That is, the equation: is solved, where the second term denotes the convective effect. In the simulations, k and τ are both set to be 1 to make the calculations easier. The choices of these coefficients are important, of course. However, it does not influence the qualitative characteristics of the equations, as it has been adopted in the reference [59][60][61][62]. The propagation patterns of the energy density with and without the second term are displayed in Figure 2. The results for the equations without the second term are plotted in solid line while the dash line denotes those for the equations with convective term. It can be seen that the heat pulse at the left boundary produces a thermal wave in the medium and it keeps the original pattern during propagation while the peak value is dissipated as shown by the solid line. When the convective term is considered, the shape is changed and a sharp wave front comes into being. The convective term shows the drift behaviors of thermomass conveyed by momentum. Thus, it is an obvious deduction when the mass nature of heat is taken into consideration. Meanwhile, the existence of convective term speeds up the thermal wave and leads to higher dissipation rate for the peak value of the patterns.  the second term are displayed in Figure 2. The results for the equations without the second term are plotted in solid line while the dash line denotes those for the equations with convective term. It can be seen that the heat pulse at the left boundary produces a thermal wave in the medium and it keeps the original pattern during propagation while the peak value is dissipated as shown by the solid line. When the convective term is considered, the shape is changed and a sharp wave front comes into being. The convective term shows the drift behaviors of thermomass conveyed by momentum. Thus, it is an obvious deduction when the mass nature of heat is taken into consideration. Meanwhile, the existence of convective term speeds up the thermal wave and leads to higher dissipation rate for the peak value of the patterns.  To show more characteristics of the momentum convection, a two-dimensional heat conduction problem is considered where heat transport occurs in a plate with finite length and width and infinite height. Then, the left boundary is heated by a partial heat pulse: where Y 2 is also set to be 1. If it is not clarified specifically, all the coefficients are set to be 1. The distribution shape of the heat flux is displayed in Figure 3. Then, the energy distribution along with time is shown in Figure 4 and different colors denotes different values of the energy density, which has a corresponding relationship with temperature. Figure 4a-c are the energy distribution propagation patterns e n for Equations (49) and (50) Figure 5. It was found that the convective term enhances the spreading behaviors of the thermal energy. In macroscale descriptions, the spreading enhancement caused by the convective term is similar with the viscous effect. However, they are different because the convective effect is reversible and the viscous effect is irreversible.
it spreads away just like water ripples. The majority of the energy propagates forward, namely along with the original direction of the heat flux, while the nonequilibrium energy distribution produces the heat flux in other directions. The difference of the energy distribution (ec-en) is plotted in Figure 5. It was found that the convective term enhances the spreading behaviors of the thermal energy. In macroscale descriptions, the spreading enhancement caused by the convective term is similar with the viscous effect. However, they are different because the convective effect is reversible and the viscous effect is irreversible.  produces the heat flux in other directions. The difference of the energy distribution (ec-en) is plotted in Figure 5. It was found that the convective term enhances the spreading behaviors of the thermal energy. In macroscale descriptions, the spreading enhancement caused by the convective term is similar with the viscous effect. However, they are different because the convective effect is reversible and the viscous effect is irreversible.

Influences of the Distortion Matrix
Distortion matrix in thermomass theory, analogous to that in solid mechanics, implies the concept of thermomass elasticity. For example:

Influences of the Distortion Matrix
Distortion matrix in thermomass theory, analogous to that in solid mechanics, implies the concept of thermomass elasticity. For example: means the compression or tension of the label field and: implies the torsion of the label field. Assuming that the Hamiltonian potential is the function of the distortion matrix: then the heat conduction is also influenced by the distortion matrix coefficients D xx and D yy , which reveals a tendency of recovering to the initial state. In the numerical simulations, the nondimensional coefficients χ xx , χ yy are set to be 0.5 and 0.5, and other coefficients are 1.
The energy distribution e D at time t = 0.3, which is derived from the equation set with Hamiltonian energy (Equation (76), is plotted in Figure 6a and the energy distribution e c which is another view of Figure 4e is displayed in Figure 6b. These two figures are compared to highlight the influences of the distortion matrix. Moreover, their difference, e D -e c , is shown in Figure 7. The distortion coefficients D xx and D yy speeds up the wave propagation and the wave shape are also altered. Meanwhile, the coefficients enhance the fluctuations of wave. In Figure 7, another small wave peak appears behind the first wave pattern, and in some sense, it shows the effect of elasticity. Considering that the distortion matrix field goes through an irreversible process at the same time, the introduction of these coefficients will add to the dissipation and strengthen the viscous effects. One obvious feature of the viscous effects is that the peak value of the energy density in Figure 6a is lower than that in Figure 6b. One obvious feature of the viscous effects is that the peak value of the energy density in Figure 6a is lower than that in Figure 6b.

Conclusions
The thermomass theory and GENERIC theory deal with the non-Fourier heat conduction phenomena in different ways. The core idea of the thermomass theory is to reveal the mass nature of heat during heat conduction based on Einstein's mass-energy relation. Therefore, heat transfer is actually the flow of thermomass gas in porous media, thus, the principle of fluid mechanics can be used to describe the heat conduction process in the media. The momentum conservation equation of thermomass gas is actually a wave equation with thermal inertia force and damping. The GENERIC theory constructs an abstract geometrical structure for various transport processes. Its core idea is to adopt the contact geometry to derive the time evolution equations of the state variables, by handling the reversible and irreversible process separately. The reversible process is expressed in energy representation governed by Hamiltonian energy, while the irreversible process is expressed in entropy representation dominated by dissipation potential.
Analyses show that the governing equations of heat conduction derived from GENERIC theory are the same as that derived from thermomass theory, if proper coefficients and potentials are

Conclusions
The thermomass theory and GENERIC theory deal with the non-Fourier heat conduction phenomena in different ways. The core idea of the thermomass theory is to reveal the mass nature of heat during heat conduction based on Einstein's mass-energy relation. Therefore, heat transfer is actually the flow of thermomass gas in porous media, thus, the principle of fluid mechanics can be used to describe the heat conduction process in the media. The momentum conservation equation of thermomass gas is actually a wave equation with thermal inertia force and damping. The GENERIC theory constructs an abstract geometrical structure for various transport processes. Its core idea is to adopt the contact geometry to derive the time evolution equations of the state variables, by handling the reversible and irreversible process separately. The reversible process is expressed in energy representation governed by Hamiltonian energy, while the irreversible process is expressed in entropy representation dominated by dissipation potential.
Analyses show that the governing equations of heat conduction derived from GENERIC theory are the same as that derived from thermomass theory, if proper coefficients and potentials are chosen. The reversible process in GENERIC theory corresponds to conversion process between the thermomass potential energy and the thermomass kinetic energy in thermomass theory, while the irreversible process in GENERIC theory, which is described by dissipation potential, corresponds to the dissipation process of thermomass energy due to the resistance in thermomass theory. GENERIC theory provides a concrete mathematical foundation for the thermomass theory, while thermomass theory gives GENERIC more physical pictures. All these indicate that both GENERIC theory and thermomass theory are not only compatible, but also supportive and complementary mutually.
The advantage of thermomass theory is that it not only reveals the physical meanings of the Hamiltonian energy and the dissipation potential in GENERIC, but also uniformly expresses these two kinds of energy as one unified quantity. Hamiltonian energy in the reversible process corresponds to the conversion process of potential energy and kinetic energy of thermomass during wave transport of heat (thermomass), and the dissipation potential in the irreversible process corresponds to the dissipation process of thermomass energy. One advantage of GENERIC theory is that it provides a flexible structure, which is easy to include more new state variables, such as the inner heat source and the distortion matrix field. Numerical methods are adopted to investigate the influences of the convective term and distortion matrix term on the wave propagation patterns. The convective term establishes a sharp wave front and enhances the spreading of thermal energy. Thus, it leads to more dispersive energy distribution. The distortion matrix field considers the viscosity and elasticity of thermomass gas and constructs a new series of heat conduction equations, as is shown in Section 5. The application of these new state variables enables thermomass theory to describe more complex heat conduction phenomena.

Conflicts of Interest:
The authors declare no conflict of interest.

Letters
Greek symbols A function of state variables

Appendix A. Details of the Numerical Simulations
In one-dimensional heat conduction problems, the governing equation set is: The second term in Equation (A2) implies the convective effect. Therefore, when the convective term is ignored, the time evolution equation (Equation (A2)) becomes: The temperature T has a linear relationship with the energy density e: where the density ρ and heat specific C V are constant. Therefore, the evolution of temperature is the same to that of thermal energy density. In a two-dimensional regime, the governing equation set becomes: The above equations are normalized by several standard quantities, to make the numerical calculations easier. The new quantities are defined as: Moreover, the following relationships are supposed to be satisfied for the standard quantities: To simplify the expression, the superscript * in the new equations is omitted. As for the one-dimensional heat conduction problem, where heat is transported in a plate with finite length, periodic boundary conditions are applied to the width and height. The nondimensional length is set to be 1. A heat flux boundary is applied at the left boundary and thermal perturbation propagates in the medium from the left to the right. The total grid number is 500 and the unit temporal step is ∆t = 1 × 10 −5 . Then, as for the two-dimensional heat conduction problem, the plate has finite length and finite width. The height boundary is set to be periodic. The grid is 200 × 200 and the temporal step is ∆t = 1 × 10 −5 . All the coefficients in the equations are set to be 1, if they are not clarified specifically. The initial conditions are set to be: As for the distortion matrix D ij , it is: The boundary condition for the left boundary is: in one-dimensional problems and: in two-dimensional problems. The other boundaries are set to be adiabatic. That is: The energy density at the boundary is derived from the energy conservation equation. The distortion matrix at the boundary are set to be: Finite difference method is adopted to solve the problem numerically. Staggered grids are used for heat flux and energy density and distortion fields. The energy and distortion fields are placed at the grid nodes while the heat fluxes are positioned between the nodes, as shown in Figure A1.
As for the time derivative, the forward difference scheme is applied: The superscript n and n + 1 represent the time levels. The convective term is solved with the first order upwind scheme. As: The spatial difference is written as: when m > 0, and it is: when m < 0. The upwind difference scheme is suitable for the wave propagation problem. (A26) in two-dimensional problems. The other boundaries are set to be adiabatic. That is: (A27) The energy density at the boundary is derived from the energy conservation equation. The distortion matrix at the boundary are set to be: Finite difference method is adopted to solve the problem numerically. Staggered grids are used for heat flux and energy density and distortion fields. The energy and distortion fields are placed at the grid nodes while the heat fluxes are positioned between the nodes, as shown in Figure A1. As for the time derivative, the forward difference scheme is applied: The superscript n and n + 1 represent the time levels. The convective term is solved with the first order upwind scheme. As: The spatial difference is written as: when m > 0, and it is: when m < 0. The upwind difference scheme is suitable for the wave propagation problem. Another difference algorithm is the central difference, which is suitable for the relationship between heat flux and temperature. For instance, in Figure A2, the derivative of temperature with respect to x coordinate is simulated as: where the heat flux q x,i+0.5,j is located between the nodes T i+1,j and T i,j . As for a given volume unit, such as Figure A2, sometimes the values which are necessary in the numerical algorithm are lacking.

(A36)
where the nearest two points are chosen. As for other values which are lacking at the given points, interpolation is widely adopted. Sometimes, the four-point interpolation method is used. For example, in Figure A2, there is no value for q y at the give point q x,i+0.5,j . Thus, the value q y,i+0.5,j is set to be: q y,i+0.5,j = q y,i+1, j+0.5 + q y,i,j+0.5 + q y,i+1,j−0.5 + q y,i,j−0.5

. (A37)
In this way, the equation set are solved explicitly to obtain the quantities in the time level n + 1 with the data in the time level n. The solution of boundary conditions is significant in this method. Here, extrapolation at the boundary are used to provide the lacking values.
, 1, +0.5 , , +0.5 , 1, -0.5 , , -0.5 , 0.5, + 4 y i j y i j y i j y i j y i j q q q q q + + + + + = . (A37) In this way, the equation set are solved explicitly to obtain the quantities in the time level n + 1 with the data in the time level n. The solution of boundary conditions is significant in this method. Here, extrapolation at the boundary are used to provide the lacking values. However, it is difficult to analyze the stability for this nonlinear numerical algorithm. Thus, we make the time steps as small as possible to make this solution stable. The grid size dependency and time step dependency have been testified (see Table A1).  However, it is difficult to analyze the stability for this nonlinear numerical algorithm. Thus, we make the time steps as small as possible to make this solution stable. The grid size dependency and time step dependency have been testified (see Table A1). The peak value of the energy density at time t = 0.1 is calculated to testify the size independence of the spatial grid and time step. The value when the grid size is 200 × 200 and time step is 1 × 10 −5 is adopted as the reference and the deviation percent at other conditions are derived. It can be seen that the current scheme is reasonable (see Table A2). Table A2. Time step dependency of the current numerical algorithm.