On the identification of two internal tensorial variables and a heat transport equation with inertial, thermal viscosity and vorticity terms

1 Introduction

Due to the recent interest in phonon hydrodynamics [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13] and, in particular, in heat vortices [14], [15], [16], [17], [18], [19], several new extensions of the Guyer–Krumhansl equation have been developed in the context of thermal transport [20], [21]. Thermodynamic foundations for generalized equations for heat transport has been given by non-equilibrium thermodynamics [22], [23], [24], [25], [26], [27], [28], [29], [30] with internal variables (see [24] and [31], [32], [33], [34]). Vortices in the motion of phonons were described on the level of kinetic theory [18], as well as in phonon hydrodynamics [19], using the framework of General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) approach (see [35], [36], [37], [38], [39]). This is in contrast with other approaches where the shear viscosity of the phonon fluid is taken into account [5], [6], [7], but nonlinear terms propagating the vortices is missing. To better represent the phonon vortices, it is convenient to enlarge such frameworks by taking into account some rotational viscosity of the phonon fluid, analogous to that occurring in complex fluids [40], [41], [42].

In a previous paper [14], two nonlocal tensorial variables (the symmetric part Q ^s and the antisymmetric part Q ^a of a tensorial variable Q) were added, and a subsequent approximation leads to a generalized Guyer–Krumhansl equation. In this paper, we further analyze this approach.

In particular, we identify a posteriori the two internal variables by means of an appropriate asymptotic method, in order to obtain a heat transfer equation more general than the one derived in [14], that contains additional contributions describing special viscous and vortical motions of the phonons inside the considered rigid heat-conducting media.

The existence of vortices does not imply the need of a rotational viscosity; vortices arise in the usual Navier–Stokes equation with only shear viscosity (or with shear viscosity and bulk viscosity). In fact, the concept of thermal viscosity is not new: it is a part of the model of phonon hydrodynamics, used for more than 20 years in the literature. It does not come from dimensional considerations, but from the mathematical structure of the equation for the heat flux, in situations in which it is analogous to the equation for velocity in viscous fluids. The essential physical aspect of rotational viscosity is the irreversible transfer of some part of the macroscopic or mesoscopic rotational energy of vortices to microscopic rotational energy of the constituent particles [40], [41], [42].

In the case of complex fluids, this refers to elongated molecules which may rotate with respect to some of their axes, and which take organized macroscopic rotational energy into disorganized molecular rotations. In the case of phonon hydrodynamics, such microscopic rotations could arise in the case that the basic constituents of the crystal are not single particles, but couples of particles which do not only vibrate around their respective equilibrium positions but also rotate around their equilibrium orientation. A part of the organized rotational momentum of the phonon fluid would then go into disorganized rotational motion of the couples constituting the crystal. Thus, it would be found not in any crystal, but in particular crystals having such microscopic constitutions.

In this sense, we are enlarging the concept of phonon hydrodynamics to a more complex interpretations than the simplest phonon hydrodynamics. Some attempts in this direction have been done by analogy of some rheological models describing non-Newtonian fluids, as power-law fluids [11] or non-linear viscoelastic fluids [12], but considering the symmetric part Q ^s of a tensorial internal variable Q. Here, instead, we also consider the antisymmetric part Q ^a of such a tensor.

The presence of non-local effects in generalized equations for heat transport is of special interest in systems in which l (the mean-free path of the heat carriers, phonons) is comparable to (or bigger than) L, the size of the system (for instance, it could be the radius of a cylindrical channel), i.e., for Knudsen number l L ≥ 1 . This may be valid for systems small enough, as the nanosystems, having their size comparable to the mean free path, or when (see [14]) “the mean free path is long enough to become comparable to the size of the system, as for low-temperature phonons (superfluid helium-4 [43]), or in rarefied systems, as for instance it occurs in aerodynamics at low environmental pressure”. In nanosystems, generalizations of the Guyer–Krumhansl equation for the heat transfer are used [10].

In this paper we use the Gyarmati approach (see [26], [27]), where the specific entropy s depends on the equilibrium state variables but also presents homogeneous quadratic terms depending on vectorial, tensorial variables, and current densities. Then, we show an analogous model formulated within the GENERIC framework.

The obtained results may be relevant in several technological fields. In [18], phonon vortices at heavy impurities in two-dimensional materials were investigated theoretically and discovered. Using the density-functional-theory, calculations for two configurations of Si (Silicon) impurities in Graphene, Si–C ₃ and Si–C ₄ and of other materials were carried out and it was shown that the vortex formation is due to the vibrational displacements of large-impurities heavy atoms at a characteristic frequency. Experimental measurements were performed using special monochromated electron microscopes, and it was discovered that vortices present circular patterns of vibrating atoms, reflecting the local symmetry and connecting to a phonon wave. Latest generations of microscopes could make more accurate investigations. It is expected that phonon vortices at defects influence the thermal conductivity and other thermal properties. In [19], a von Kármán vortex street was predicted for higher heat fluxes in the presence of impurities.

The paper is organized as follows. In Section 2 we recall the model for rigid media formulated in [14] while providing further clarifications. The model has two non-local macroscopic internal variables, the symmetric part and the antisymmetric part of a second order tensor Q, Q ^s and Q ^a . The evolution equations of Q ^s and Q ^(a) and a generalized Guyer–Krumhansl equation, describing kinematic and viscous motions of phonons with a rotational viscosity are presented. In Section 3, we derive other forms of the generalized heat equation obtained in Section 2, bringing in closer to the Guyer–Krumhansl equation. Section 4 shows similarities between our model, hydrodynamics of micropolar fluids, and hydrodynamics of phonons. In Section 5, we identify the two internal tensorial variables, Q ^s and Q ^(a), by asymptotic expansion of the evolution equations, keeping the zeroth and first orders. The small parameters with respect to which the equations are expanded express the ratios of the relaxation times of each variable and the typical timescale. Furthermore, the stability conditions of equilibrium solutions of this obtained equation are studied. Finally, Section 6 contains a geometric analogue of the models, formulated within the GENERIC approach, which suggests an alternative approach for the formulation of Extended Irreversible Thermodynamics and leads to a non-Clausius relation between the entropy flux and heat flux.

2 Recalling the model

The aim of the model presented in [14] is to describe diffusive, ballistic, viscous, and vortical motions of phonons. On the macroscopic scale it is possible to investigate all these behaviours within the framework of non-equilibrium thermodynamics with internal variables. In this Section, we recall the model, that leads to a generalized Guyer–Krumhansl equation, while providing further insight.

The thermodynamic state space C contains, besides the internal energy u and the heat flux J ^(q), two non-local macroscopic internal variables, Q ^s and Q ^a , which express the symmetric part and the antisymmetric part of a second order tensor Q, C = C u , J ( q ) , Q s , Q a .

Balance laws within the thermodynamic state space are described by the following equations:

1. Internal energy balance equation:

   (1) ρ d u d t + ∇ ⋅ J ( q ) = 0 ,

   where ρ is the mass density field of the considered medium (supposed to be constant), d/dt is the material time derivative, u is the internal energy density, and J ^(q) is the heat flux.

2. Entropy inequality:

   (2) σ ( s ) = ρ d s d t + ∇ ⋅ J ( s ) ⩾ 0 ,

   where the fields of specific entropy s, the entropy flux density J ^(s), and the entropy production density σ ^(s) are constitutive quantities, functions of the independent variables of the thermodynamics state space.

As neither the heat flux nor the internal variables are a priori conserved, we employ the Gyarmati approach [26], [27] to find the remaining evolution equations. The specific entropy s is not function only on the local equilibrium state variables, but also depends on J ^(q), Q ^s , and Q ^a . To ensure stability, entropy is a concave function of the state variables, and thus near equilibrium it can be approximated by a quadratic form,

where s ^(eq)(u) is the equilibrium entropy function, α ₁, α 2 s and α 2 a are inductivity constants (positive for the sake of thermodynamic stability), and the symbol “:” denotes double contraction, for example Q s : Q s = Q i j s Q i j s   ( i , j = 1,2,3 ) .

Derivatives of the entropy with respect to the state variables are then

The entropy flux J ^(s) is assumed having the form

where β 1 s and β 2 s are constant material coefficients, that conveniently represent the deviation of J ^(s) from the local equilibrium. Quantity Q i ( a ) = − ϵ ijk Q j k a is the axial vector associated with antisymmetric tensor Q ^a . The symbol “⋅” indicates one contraction, for instance Q s ⋅ J ( q ) = Q i j s J j ( q ) , and we have used the formula Q ^a ⋅ J ^(q) = Q ^(a) × J ^(q), with the symbol “×” indicating the vector product, see Appendix A for more details.

The evolution equations for the independent variables J ^(q), Q ^s and Q ^(a) are then found from the entropy production, which is obtained by taking the time derivative of entropy,

where the upper dot “^⋅” denotes the total time derivative. The entropy production σ ^s , see Equation (2), is then

In (7) we have taken into account that ∂ ∂ x i ( Q i j J j ( q ) ) = Q i j , i J j ( q ) + Q i j J j , i ( q ) = ( Q j i ) , i T J j ( q ) + Q i j ( J i , j ( q ) ) T = J j ( q ) ( Q j i ) , i T + Q i j ( J i , j ( q ) ) T , where the symbol “ ^T ” means transposition and a comma in lower indices denotes the partial spatial derivate, for instance J i , j ( q ) = ∂ J i ( q ) ∂ x j = ∇ J ( q ) i j . Furthermore, we have decomposed the gradient of heat flux ∇J ^(q) in its symmetric part ∇ J ( q ) s and antisymmetric part ∇ J ( q ) a . Note also that we consider only the case of a rigid heat conductor at rest, so that the substantial time derivative is equal to the partial temporal derivative ∂ ∂ t .

Using the expression ∂ ∂ x i ( ϵ ijk Q j ( a ) J k ( q ) ) = ϵ ijk Q j , i ( a ) J k ( q ) + ϵ ijk Q j ( a ) J k , i ( q ) = J k ( q ) ϵ kij Q j , i ( a ) − Q j ( a ) ϵ jik J k , i ( q ) , the formula for entropy production can be then also rewritten as

where the following notation is used, β 1 a = β 1 ( a ) , and α 2 ( a ) = 2 α 2 a , as well as identities   Q i j a Q ̇ i j a = ϵ ijk Q k ( a ) ϵ ijl Q ̇ l ( a ) = ( δ j j δ k l − δ l j δ k j ) Q k ( a ) Q ̇ l ( a ) = 3 Q ( a ) ⋅ Q ̇ ( a ) − Q ( a ) ⋅ Q ̇ ( a ) = 2 Q ( a ) ⋅ Q ̇ ( a ) .

The entropy production σ ^s (8) is given by the sum of products of thermodynamic fluxes and forces. Assuming that the media under consideration are perfectly isotropic and thus have the symmetry property of being invariant under orthogonal transformations, we obtain the following linear relations among the chosen thermodynamic fluxes J ^(q), Q ^s , Q ^(a) and the corresponding thermodynamic affinities:

where λ ₁ (heat conductivity), λ 2 s , and λ 2 ( a ) are positive scalar coefficients. Note that J ^(q), Q ^s , and Q ^(a) are a polar vector, a second order tensor, and an axial vector, respectively, that do not couple for the Curie principle [23].

Then, in the isotropic case the Onsager matrix is diagonal. In Section 6, we show that this diagonality is partly a consequence of the Hamiltonian structure of the model (taking into account only the reversible evolution).

The two equations (9b) and (9c) are coupled with the evolution equation for the heat flux (9a). If we assume that these two equations relax quickly to their quasi-stationary states where Q ̇ s = 0 = Q ̇ ( a ) , the remaining equation for the heat flux turns to a generalized Guyer–Krumhansl equation

with

which was identified in [14].

In (10), coefficients l 1 2 and l 2 2 represent the square of the mean free path of phonons related to the exchange of linear momentum and internal angular momentum, respectively. Coefficient l 1 2 is multiplied by a term describing weakly non-local effects and, in particular, the viscous motions of phonons. Coefficient l 2 2 is multiplied by a non-local term related to the rotational motion of the heat flux, describing rotational viscous motion of the heat carriers in crystals of elongated or polarized constituents. These terms are analogous to the terms that in the Navier–Stokes equation for viscous fluids describe the shear and rotational viscous effects of fluid particles.

When coefficient l 2 2 can be neglected, Equation (10) reduces to the Guyer–Krumhansl heat equation,

describing both diffusive motion of phonons and ballistic motion of phonons subjected to collisions with the walls of the solid.

3 Enlarged Guyer–Krumhansl description

In this Section we show how Equation (10) is connected with various levels of description of beyond-Fourier heat conduction. Heat equation (10) was found by assuming that Q ̇ s and Q ̇ ( a ) were negligible in their evolution equations (9b) and (9c). In other words, we assumed that the tensorial quantities Q ^s and Q ^(a) relax quickly to their respective quasi-equilibria, characterized by Q s = − λ 2 s β 1 s ( ∇ J ( q ) ) s and Q ( a ) = λ 2 ( a ) β 1 ( a ) ∇ × J ( q ) . When we plug these constitutive relations into the remaining evolution equation for the heat flux (9a), we obtain

which is another form generalized Guyer–Krumhansl equation (10).

In case of β 1 ( a ) = 0 , when the entropy does not depend on the skew-symmetric tensor Q ^a , Equation (13) simplifies to the Guyer–Krumhansl equation,

see Equation (65) in [20].

Moreover, if we assume that the entropy is independent of the symmetric tensor Q ^s as well, the equation further simplifies to the Cattaneo equation [44],

describing at the macroscopic scale thermal signals (second-sound phenomena) with relaxation time τ and finite velocity of propagation.

Finally, if the dependence of entropy on the heat flux is negligible as well, the equation restores the usual Fourier heat conduction J ^(q) = −λ ₁∇T, describing thermal signals with relaxation time τ null, valid at large space scales and at low frequencies.

4 Hydrodynamic analogies

Equation (10) may be of special interest in the context of phonon hydrodynamics. It contains the usual Guyer–Krumhansl equation (the third term on the right-hand side), that is related to the viscous behaviour of phonons, and the fourth term on the right-hand side, that describes the rotational viscous behaviour arising from rotational motions of the phonon fluid (heat vortices). This term may be physically relevant when the particles constituting the crystal are not single particles but dipoles (see Figure 1). In this case, the rotational viscosity describes the transfer from organized macroscopic rotational motion of the heat vortex to microscopic disorganized individual rotations of the dipoles.

Figure 1:

The dots represent particles. In (a), the crystal is made by single particles; in (b), it is made by dipoles. In (a), the particles vibrate; in (b), they vibrate and rotate. The crystal may have any kind of geometry; in particular, we may assume, for the sake of simplicity, that it has a cubic crystallization.

Equation (10) can be also written in the form

with ν ₁ and ν ₂ playing the role of shear and kinematic viscosities, given by ν 1 = l 1 2 τ and ν 2 = l 2 2 τ . J ^(q) is given by J ^(q) = (ρc _v T)v, with v being the average velocity of phonon flow and c _v the specific heat per unit mass (see for instance [19]).

From the point of view of hydrodynamic analogies, Equation (16) is similar to the Navier–Stokes equation for micropolar fluids, having elongated, not spherical molecules,

where v is the velocity field, ∇p the pressure field, μ ₀, μ ₁ and μ ₂ are the bulk, the shear and the rotational viscosities, and ω is the local spin or intrinsic angular momentum (not considered in this paper). The rotational degrees of freedom, that vibrate with respect to their equilibrium position, contribute to the macroscopic rotational energy (see [45], [46], [47], [48]). In (17), ω appears as related to the balance equation of the angular momentum, we consider ω = 0. If the liquid flows in a porous environment, an additional term − α v ρ on the right-hand side of Navier–Stokes equation (17) could be present (a Darcy term), α being a porosity coefficient. Such Darcy term would play a role analogous to the first term on the right side-hand side of equation (16), J ( q ) τ .

Moreover, we can have in a solid a lattice formed by elongated (not spherical) particles (as for instance dipoles), free to rotate around their center of mass, heat transfer will incorporate vibrational energy and rotational energy. The interaction of the heat vortices with the rotational degrees of freedom of the components of the lattice, will transfer collective rotational motion of the heat flux to individual disordered random motion. This transfer contributes to dissipation and is phenomenologically described by the rotational viscosity.

From the point of view of the analogy between the hydrodynamics of phonons and the hydrodynamics of Navier Stokes fluids, vortices are found from the Navier–Stokes equation for isotropic fluids using the Reynolds method, that is used to study turbulence. This method leads to finding a term called Reynolds stress tensor, that describes vortices in the isotropic fluids (see [40]).

In this paper the two internal variables Q ^s and Q ^(a) describe thermal phenomena inside the medium from the macroscopic point of view, as for instance the introduction of the “anisotropy-grain tensor” describes particular physical fields responsible for the internal structure and/or geometry of anisotropic media in ferroelastic crystals. In Section 5, we add some remarks about the use of the internal variable as a powerful tool to model the behaviour of media.

5 Identification of the two internal variables Q ^s and Q ^a and a heat equation with special thermal viscosity terms

In Sections 2 and 3, the two internal variables Q ^s and Q ^a were eliminated to obtain the generalized heat transport equations (10) and (13). In the current Section, we determine the expressions for the two internal variables in more detail.

Internal variables can be used to describe complex phenomena continuous media [24], [27], [31], [32], [33], [34]. In several papers, the macroscopic internal variables are introduced by their definitions to describe the internal structure of a medium. For example, the anisotropy-grain tensor describes particular physical fields responsible for the internal structure and/or geometry in ferroelastic crystals [49]. It is assumed that the medium as a whole is homogeneous and isotropic, where there exist particular physical fields responsible for its internal anisotropy described by internal variables. There are other example of internal variables also in rheology and other physical fields [14]. In [50] in a model for media presenting several microscopic phenomena of dielectric polarization, n hidden variables influencing the dielectric polarization are introduced. By means of a particular demonstration, they are interpreted as n partial specific polarization polar vectors, which are parts of a total specific polarization and describe the dielectric relaxation of polar isotropic media [51]. Using these specific partial polarizations as internal variables, a generalized Debye equation for the dielectric relaxation is obtained for isotropic media. Moreover, in [52], a model for isotropic media presenting several microscopic phenomena of magnetization showed that it is possible to describe these microscopic phenomena by introducing n partial specific magnetizations as internal variables. Using these specific partial magnetizations, a generalized Snoek equation for magnetic relaxation phenomena was obtained by eliminating the internal variables. This is the case of complex media where different kinds of molecules, contributing to the total magnetization, have different magnetic susceptibilities and relaxation times. These physical situations arise for example in magnetic resonance in medicine, biology and other science fields, where different molecules present particular magnetic relaxation phenomena. In other papers the internal variables are eliminated by an appropriate method from the formalism as in [13]. In Reference [31], Muschik defines the internal variable by seven concepts, among which: the internal variables need a microscopic/molecular model for their physical interpretation. In this Section, we use asymptotic expansion to identify the two internal variables Q ^s and Q ^a . It seems that this asymptotic method to determine internal variables has not been used in literature so far for the purpose of internal variable identification. It allows to find a generalized Guyer–Krumhansl-like heat equation that describes more complex heat phenomena.