Non-convex dissipation potentials in multiscale non-equilibrium thermodynamics

TL;DR

This paper explores how non-convex dissipation potentials influence stability and phase transitions in non-equilibrium thermodynamics, providing insights into complex fluid behaviors and supercritical boiling.

Contribution

It introduces the concept of CR-stability related to dissipation potential convexity and develops evolution equations for perturbations in a thermodynamically consistent framework.

Findings

Loss of convexity leads to dissipative phase transitions.

CR-stability analysis explains complex fluid behaviors.

Non-convex dissipation potentials relate to experimental phenomena.

Abstract

Reformulating constitutive relation in terms of gradient dynamics (being derivative of a dissipation potential) brings additional information on stability, metastability and instability of the dynamics with respect to perturbations of the constitutive relation, called CR-stability. CR-instability is connected to the loss of convexity of the dissipation potential, which makes the Legendre-conjugate dissipation potential multivalued and causes dissipative phase transitions, that are not induced by non-convexity of free energy, but by non-convexity of the dissipation potential. CR-stability of the constitutive relation with respect to perturbations is then manifested by constructing evolution equations for the perturbations in a thermodynamically sound way. As a result, interesting experimental observations of behavior of complex fluids under shear flow and supercritical boiling curve can be explained.

Figures

Equations

$$\mathbb{D}=\left[\alpha {\left(1+\beta {\left|{\mathbb{T}}_{\delta }\right|}^2\right)}^n+\gamma \right]{\mathbb{T}}_{\delta },$$

$$\mathfrak{h}({\mathbb{T}}_{\delta },\mathbb{D})=0,$$

$$\mathbb{D}=\frac{\partial {\Xi }^{\ast }}{\partial {\mathbb{T}}_{\delta }},$$

$${\mathbb{T}}_{\delta }=\frac{\partial \Xi }{\partial \mathbb{D}},$$

$${\Xi }^{\ast }\left({\mathbb{T}}_{\delta }\right)=\{\begin{array}{ll}\frac{\alpha }{2\beta (n+1)}\left[(1+\beta {\left|{\mathbb{T}}_{\delta }\right|}^2{)}^{n+1}-1\right]+\frac{\gamma }{2}{\left|{\mathbb{T}}_{\delta }\right|}^2,& n\ne -1,\\ \frac{\alpha }{2\beta }\ln \left(1+\beta {\left|{\mathbb{T}}_{\delta }\right|}^2\right)+\frac{\gamma }{2}{\left|{\mathbb{T}}_{\delta }\right|}^2,& n=-1,\end{array}$$

$$⟨\mathbb{D},{\mathbb{T}}_{\delta }⟩=\underset{\ge 0}{\underbrace{\Xi (\mathbb{D})}}+\underset{\ge 0}{\underbrace{{\Xi }^{\ast }({\mathbb{T}}_{\delta })}}\ge 0,$$

$$⟨\mathbb{D},{\Xi }_{\mathbb{D}}⟩\ge 0\text{mbox}or⟨{\mathbb{T}}_{\delta },{\Xi }_{{\mathbb{T}}_{\delta }}^{\ast }⟩\ge 0.$$

$$0\le {\mathbb{T}}_{\delta }:{\Xi }_{{\mathbb{T}}_{\delta }}^{\ast }=\left|{\mathbb{T}}_{\delta }\right|\frac{{\mathbb{T}}_{\delta }}{\left|{\mathbb{T}}_{\delta }\right|}:{\Xi }_{|{\mathbb{T}}_{\delta }|}^{\ast }\frac{\partial \left|{\mathbb{T}}_{\delta }\right|}{\partial {\mathbb{T}}_{\delta }}=\left|{\mathbb{T}}_{\delta }\right|\frac{{\mathbb{T}}_{\delta }}{\left|{\mathbb{T}}_{\delta }\right|}:{\Xi }_{|{\mathbb{T}}_{\delta }|}^{\ast }\frac{{\mathbb{T}}_{\delta }}{\left|{\mathbb{T}}_{\delta }\right|}=\left|{\mathbb{T}}_{\delta }\right|{\Xi }_{|{\mathbb{T}}_{\delta }|}^{\ast },$$

$$\frac{\partial \mathbf{u}}{\partial t}=\underset{=\text{mbox}reversibleterms}{\underbrace{-\nabla p-\mathrm{div}\left(\frac{\mathbf{u}\otimes \mathbf{u}}{\rho }\right)}}+\mathrm{div}({\mathbb{T}}_{\delta }),$$

$${{\Xi }_{{\mathbb{T}}_{\delta }}^{\ast }|}_{{\tilde{\mathbb{T}}}_{\delta }(\mathbb{D})}=\mathbb{D},$$

$$\frac{\partial{}}{\partial{{{{\mathbb{T}}}_{\delta}}}}\left(-\Xi^{*}({{{\mathbb{T}}}_{\delta}})+\langle{{{\mathbb{T}}}_{\delta}},\mathbb{D}\rangle\right)\Big{|}_{\tilde{{\mathbb{T}}}_{\delta}(\mathbb{D})}=\mathbb{0},$$

$$\Xi(\mathbb{D})=_{\mathrm{def}}\left(-\Xi^{*}({{{\mathbb{T}}}_{\delta}})+\langle{{{\mathbb{T}}}_{\delta}},\mathbb{D}\rangle\right)\Big{|}_{\tilde{{\mathbb{T}}}_{\delta}(\mathbb{D})}.$$

$${\Xi }_{\mathbb{D}}={\tilde{\mathbb{T}}}_{\delta }(\mathbb{D}),$$

$$\frac{\partial{}}{\partial{\mathbb{D}}}\left(-\Xi(\mathbb{D})+\langle\mathbb{D},{{{\mathbb{T}}}_{\delta}}\rangle\right)\Big{|}_{\tilde{\mathbb{D}}({{{\mathbb{T}}}_{\delta}})}=\mathbb{0},$$

$${\Xi }^{\ast }({\mathbb{T}}_{\delta })=-\Xi \left(\tilde{\mathbb{D}}({\mathbb{T}}_{\delta })\right)+⟨\tilde{\mathbb{D}}({\mathbb{T}}_{\delta }),{\mathbb{T}}_{\delta }⟩,$$

$${\int }_{|{\mathbb{T}}_{\delta }|(B)}^{|{\mathbb{T}}_{\delta }|(E)}|\mathbb{D}(T)|\mathrm{d}T={\int }_{|{\mathbb{T}}_{\delta }|(B)}^{|{\mathbb{T}}_{\delta }|(E)}{\Xi }_{|{\mathbb{T}}_{\delta }|(T)}^{\ast }\mathrm{d}T={\Xi }^{\ast }(|{\mathbb{T}}_{\delta }(E)|)-{\Xi }^{\ast }(|{\mathbb{T}}_{\delta }(B)|),$$

$$\frac{{\Xi }^{\ast }({\mathbb{T}}_{\delta }(E))-{\Xi }^{\ast }({\mathbb{T}}_{\delta }(B))}{|{\mathbb{T}}_{\delta }(E)|-|{\mathbb{T}}_{\delta }(B)|}={\Xi }_{|{\mathbb{T}}_{\delta }(E)|}^{\ast }={\Xi }_{|{\mathbb{T}}_{\delta }(B)|}^{\ast }.$$

$${\int }_{|{\mathbb{T}}_{\delta }|(B)}^{|{\mathbb{T}}_{\delta }|(E)}|\mathbb{D}(T)|\mathrm{d}T=(|{\mathbb{T}}_{\delta }|(E)-|{\mathbb{T}}_{\delta }|(B))\underset{=|\mathbb{D}(|{\mathbb{T}}_{\delta }|(B)|)|}{\underbrace{|\mathbb{D}(|{\mathbb{T}}_{\delta }(E)|)|}},$$

$$\frac{\partial \mathbf{u}}{\partial t}={\left(\frac{\partial \mathbf{u}}{\partial t}\right)}_{rev}+{{\Xi }_{{\mathbf{u}}^{\ast }}|}_{{\mathbf{u}}^{\ast }={\Psi }_{\mathbf{u}}^{(CH\to E)}}.$$

$${\Psi }^{(CH\to E)}=\int \mathrm{d}\mathbf{r}S(\rho ,\mathbf{u},e)-\frac{1}{T_0}E(\rho ,\mathbf{u},e)+\frac{{\mu }_0}{T_0}M(\rho ,\mathbf{u},e),$$

$$S=\int \mathrm{d}\mathbf{r}s\left(\rho ,e-\frac{{\mathbf{u}}^2}{2\rho }\right),$$

$${\mathbf{u}}^{\ast }={\left(\frac{\partial {\Psi }^{(CH\to E)}}{\partial \mathbf{u}}\right)}_{\rho ,s}=S_{\mathbf{u}}=-\frac{1}{T}\frac{\mathbf{u}}{\rho }.$$

$$\Xi =\Xi \left(\Gamma ({\mathbf{u}}^{\ast })\right),$$

$$\Gamma (\bullet )=-\frac{1}{2}\left(\nabla \bullet +(\nabla \bullet {)}^{\top }\right),$$

$${\Xi }_{{\mathbf{u}}^{\ast }}=⟨{\Xi }_{\mathbb{D}},\frac{\partial \mathbb{D}}{\partial {\mathbf{u}}^{\ast }}⟩=⟨{\Xi }_{\mathbb{D}},\frac{\partial \Gamma ({\mathbf{u}}^{\ast })}{\partial {\mathbf{u}}^{\ast }}⟩={\Gamma }^{\top }\left({\Xi }_{\mathbb{D}}\right),$$

$$\int \mathrm{d}\mathbf{r}\Gamma (\mathbf{a})\cdot \mathbf{b}=\int \mathrm{d}\mathbf{r}\mathbf{a}\cdot {\Gamma }^{\top }(\mathbf{b}),\forall \text{mbox}symmetrictensorfields\mathbf{a},\mathbf{b}.$$

$$\frac{\partial \mathbf{u}}{\partial t}={\left(\frac{\partial \mathbf{u}}{\partial t}\right)}_{rev}+{\Gamma }^{\top }({\mathbb{T}}_{\delta }),$$

$$\frac{\partial \mathbf{u}}{\partial t}$$

$$\frac{\partial {\mathbb{T}}_{\delta }^{\ast }}{\partial t}$$

$${\mathbb{T}}_{\delta }={\Psi }_{{\mathbb{T}}_{\delta }^{\ast }}^{(EH\to E)},\text{ or }\frac{\partial }{\partial {\mathbb{T}}_{\delta }^{\ast }}\left(-{\Psi }^{(EH\to E)}+\int \mathrm{d}\mathbf{r}{\mathbb{T}}_{\delta }:{\mathbb{T}}_{\delta }^{\ast }\right)=0,$$

$${\Psi }^{(EH\to E)}{=}_{\mathrm{def}}{S}^{(EH\to E)}-\frac{1}{T_0}{E}^{(EH\to E)}+\frac{{\mu }_0}{T_0}{M}^{(EH\to E)}.$$

$$S^{(EH)}=S^{(CH)}-\int \mathrm{d}\mathbf{r}\frac{1}{2}\beta {\left|{\mathbb{T}}_{\delta }^{\ast }\right|}^2=\int \mathrm{d}\mathbf{r}s\left(\rho ,e-\frac{{\mathbf{u}}^2}{2\rho }\right)-\frac{1}{2}\beta {\left|{\mathbb{T}}_{\delta }^{\ast }\right|}^2.$$

$${\Psi }^{(EH\to E)}={\Psi }^{(CH\to E)}-\int \mathrm{d}\mathbf{r}\frac{1}{2}\beta {\left|{\mathbb{T}}_{\delta }^{\ast }\right|}^2,$$

Full text

††journal: Continuum Mechanics and Thermodynamics

1 Introduction

Consider two concentric cylinders, the outer rotating with respect to the inner one. The gap between the cylinders is filled with a non-Newtonian fluid. In such an experiment either the force required to rotate the outer cylinder (shear stress) or speed of the rotation (shear rate) can be controlled, the other being measured, see Fig. 1.111We assume that the measured shear stress and shear rate corresponds to the actual properties of the fluid within the gap. An interesting experimental observation is that when varying shear stress, shear rate behaves continuously, while when varying shear rate, shear stress exhibits a jump boltenhagen.p.hu.y.ea:observation .

We ask the following two questions:

How to describe the relation between shear stress and shear rate by means of a constitutive relation compatible with non-equilibrium thermodynamics? 2. 2.

What is the reason that when controlling the shear stress, the response (shear rate) is continuous, while when controlling the shear rate, the response (shear stress) is discontinuous?

The first question is answered in Sec. 2, where an appropriate model is fitted to the experimental data and then reformulated in terms of gradient dynamics. The second question is partially answered in Sec. 3, where the multivalued conjugate gradient structure is identified. Subsequently, in Sec. 4, the dynamics is lifted to an extended state space. This enables us to discuss various regimes of the dynamics and its stability, metastability and instability with respect to perturbations of the constitutive relation—the CR-stability, and to complete answering the second question in Sec. 5.1, which is the essential result of this paper. Moreover, we propose that the phenomenon of critical heat flux in Sec. 8 could be approached by the method of non-convex dissipation potentials.

Novelty of this paper lies in the following points. A dissipation potential is identified that yields constitutive relation (1). The dissipation potential is not convex everywhere (it is convex in the vicinity of equilibrium), and CR-stability of the dynamics generated by this dissipation potential (with respect to perturbations of the dissipation potential) is discussed in the context of mesoscopic multiscale thermodynamics. The conjugate dissipation potential, which is obtained by a form of Legendre transformation, is multivalued, which, to our best knowledge, is mentioned in the thermodynamic literature for the first time whatsoever. Moreover, a thermodynamic reason is identified which favors seeking constitutive relation in the form $\mathbb{D}({{{\mathbb{T}}}_{\delta}})$, i.e. shear rate as a function of shear stress. Finally, the experimental behavior presented in Fig. 1 is explained. A model based on a non-convex dissipation potential is also proposed for the critical heat flux in Sec 8.

2 Constitutive relation

The experimental dependence of shear rate on shear stress plotted in Fig. 1 can be described by a model of non-Newtonian fluid introduced by Le Roux and Rajagopal roux.c.rajagopal.kr:shear

[TABLE]

where $\alpha$ and $\beta$ are positive numbers, $\gamma$ is a non-negative number, $n$ is a real number, $\mathbb{D}=_{\mathrm{def}}\frac{1}{2}\left(\nabla\mathbf{v}+\nabla\mathbf{v}^{\top}\right)$ is the symmetric part of the velocity gradient and ${{{\mathbb{T}}}_{\delta}}=_{\mathrm{def}}{\mathbb{T}}-\left(\frac{1}{3}\operatorname{Tr}{\mathbb{T}}\right)\mathbb{I}$ is the deviatoric part of the Cauchy stress ${\mathbb{T}}$. Setting $\gamma=0$, constitutive relation (1) reduces to the stress power-law model studied by Málek et al. malek.j.prusa.ea:generalizations . The fit is plotted in Fig. 1.

Note that only the deviatoric part of the stress tensor is taken into account and that the fluid is considered to be effectively incompressible. It is possible that abandoning such a constraint could lead to other interesting behavior, but the purpose of this paper is achieved while respecting this constraint.

Constitutive relation (1) can be seen as a special case of general implicit constitutive relation rajagopal.kr:on*3 ; rajagopal.kr:on*4

[TABLE]

where ${\mathfrak{h}}$ is a tensorial function. It is however not clear whether the fully implicit case brings any advantage and whether it is physically substantiated. On the other hand, the special case $\mathbb{D}({{{\mathbb{T}}}_{\delta}})$ has a sound thermodynamic meaning and practical advantages, see Sec. 4 and perlacova.t.prusa.v:tensorial .

Constitutive relation (1) can be also regarded as a consequence of gradient dynamics222Gradient dynamics is a cornerstone of the GENERIC framework GO ; OG , where it generates the dissipative evolution. It has strong connection to the principle of large deviations Mielke-Peletier , and it can be related to other techniques associated with the maximization of entropy production SEA ; Rajagopal2004 , where for instance quantum dissipation is addressed. However, we prefer using gradient dynamics, since it seems to be more straightforward in the case of non-quadratic dissipation potentials janecka.a.pavelka.m:gradient . Note that an alternative non-potential version of GENERIC is preferred by others HCO ; Hutter2013 . Although the latter version is more general as it allows for antisymmetric dissipative coupling, it is not clear whether such a coupling is necessary and correct Miroslav-WhyGeneric . We prefer the potential version of GENERIC because reversibility and irreversibility are then clearly distinguished PRE2014 and because we can then perform Legendre transformations, which are crucial for the manuscript at hand.

[TABLE]

where the dissipation potential $\Xi$ depends on the shear rate, and it is conjugate to $\Xi^{*}$ by means of the Legendre transformation, see Sec. 3 for more details. The former case, Eq. (3a), is covered by dissipation potential

[TABLE]

whose derivative gives Eq. (1) in the sense of Eq. (3a). This dissipation potential is plotted in Fig. 2. The latter case, Eq. (3b), is covered by the dissipation potential $\Xi$ specified in Sec. 3.

The important feature of the dissipation potential $\Xi^{*}$, Eq. (4), is that it is convex in the vicinity of zero while non-convex for higher shear stresses. Convexity is then regained for even higher shear stresses. By analogy with equilibrium thermodynamics Callen , it can be expected that the convex parts generate CR-stable or CR-metastable evolution while the non-convex parts lead to CR-instability. This is indeed so as discussed in Sec. 4 and Sec. 5.1.

2.1 Second law of thermodynamics

On the second law of thermodynamics and convexity: Assume for a moment that the dissipation potential $\Xi(\mathbb{D})$ is smooth and convex. The potential is zero at the origin an thus is positive everywhere except in the thermodynamic equilibrium. Due to the convexity, the origin is also the only point where derivative of the dissipation potential is zero, ${{{\mathbb{T}}}_{\delta}}=\Xi_{\mathbb{D}}|_{\mathbb{D}=0}=\mathbb{0}$. The conjugate dissipation potential $\Xi^{*}({{{\mathbb{T}}}_{\delta}})$, which is formed by Legendre transformation, is also convex and zero at the origin. Moreover, it holds that

[TABLE]

see e.g. p. 267 in Roubicek . Entropy production is the product of thermodynamic forces and fluxes as usually in dGM or in the GENERIC evolution equation for momentum (14), and this last inequality represents the positive sign of entropy production. Convexity of the dissipation potential thus implies compatibility with the second law of thermodynamics.

On the second law and non-convexity: Convexity of the dissipation potential is a sufficient condition for compatibility with the second law. However, it is not a necessary condition. The second law is fulfilled if and only if the dissipation potential generates non-negative entropy production,

[TABLE]

Practically always in non-equilibrium thermodynamics this condition is required to hold locally (not only after integration over the whole system). Moreover, let us consider only the case where the dissipation potential depends only on a norm of the thermodynamic force or flux (as always in this paper). The local version of condition (6) is then

[TABLE]

which is satisfied for instance if the function $\Xi^{*}(|{{{\mathbb{T}}}_{\delta}}|)$ is monotonous, see Fig. 2. In this sense, a non-convex dissipation potential which is monotonous as a function of $\left|{{{\mathbb{T}}}_{\delta}}\right|$ fulfills the second law of thermodynamics.

A generalization to dissipation potentials that are not necessarily simple functions of a norm of the forces or fluxes was found for instance in Mira-Milan (Eqs. (28)-(30) therein).

It should be also borne in mind that condition (6) expresses non-negativity of entropy production only if the thermodynamic potential whose derivatives generate the thermodynamic forces is convex. Therefore, we always assume that entropy is a concave function of the state variables and that the thermodynamic potential is convex.

We believe that the experimental observations can be understood also without non-convex dissipation potentials, but on more detailed levels of description. The non-convex dissipation potential should be a result of a reduction from a more detailed level of description equipped with a convex dissipation potential. It is not clear, however, what is the more detailed level of description. Perhaps a further extension in the sense of internal variables (that can be often seen as a CR-extension) could generate the same results using a convex dissipation potential. But unless such extension is known, let us proceed with a non-convex dissipation potential.

2.2 Approach to equilibrium and existence

Convexity of the dissipation potential near equilibrium, which can be interpreted as approximately quadratic behavior, is moreover in tight relation with stability of the equilibrium state. Apart from the convexity of the dissipation potential, convexity of the energy plays a crucial role in stability analysis. When these two requirements are satisfied, the equilibrium state is stable with respect to small fluctuations, as can be seen for example from analysis of the spectrum of the linearized GENERIC near equilibrium, see PRE2014 , or as discussed in the paper by Matolcsi Matolcsi-Lyapunov . Thermodynamic equilibrium is stable if the dissipation potential is convex and the entropy concave near equilibrium.

On the other hand, although convexity of the dissipation potential (also monotonicity of the derivative) implies increase of the entropy by condition (6), it does not make up for an exact mathematical proof of the existence of solutions and of the approach to equilibrium from arbitrarily far-from-equilibrium states. It is, however, an important and essential ingredient from both the physical and the mathematical points of view. This can be well illustrated on the case of the Boltzmann equation where the exact mathematics has been done.

Boltzmann equation is an illustration of GENERIC with the non-quadratic but convex dissipation potential, see GO . The proof of the global existence of solutions has been done by Di Perna and Lions Lions , for which the latter author received the Fields Medal. The Boltzmann $H$-theorem is not sufficient but plays an essential role in the proof. The approach to the equilibrium state was also shown by Desvilletes and Villani Desvilettes-Villani , where again the Boltzmann $H$-theorem does not suffice but plays an essential role.

Summing up, the increase of the entropy as implied by monotonicity of the dissipation potential does not, indeed, imply immediately all the mathematical rigor in the proofs of the existence of solutions and the approach to the equilibrium from any far-from-equilibrium state, but it is an essential ingredient in such proofs (that have been done only for a few mesoscopic time evolution equations so far).

2.3 Shear and vorticity banding

Let us now return to Fig. 1. It is possible to interpret the $\mathbb{D}({{{\mathbb{T}}}_{\delta}})$ dependence in Fig. 1 (open triangles) in terms of gradient dynamics. For each ${{{\mathbb{T}}}_{\delta}}$ one has a unique value of $\mathbb{D}$, which is the slope of the function $\Xi^{*}$ at the particular value of ${{{\mathbb{T}}}_{\delta}}$. In the next section, we try to answer how to interpret the dependence ${{{\mathbb{T}}}_{\delta}}(\mathbb{D})$.

Before going on with the details on the conjugate dissipation potential, let us briefly comment on the relation of this manuscript to the gradient and vorticity banding phenomena, excellently reviewed in Fielding ; Olmsted . Shear banding is usually experimentally observed in the cylindrical Couette flow (as in this paper), where there might develop instabilities in form of (visible) bands. If the bands develop in the radial direction, the phenomenon is referred to as gradient banding while if they develop in the direction of the rotational axis, the phenomenon is referred to as vorticity banding.

Gradient banding usually leads to stress–strain curves with multiple shear rates corresponding to a single shear stress (”horizontal” S-curve) while vorticity banding leads to curves where multiple shear stresses correspond to a single shear rate (”vertical” S-curve). Although such distinction between the two phenomena is typical, it is not always true Olmsted . When crossing the region of unstable shear rates in the case of gradient banding, the region of higher shear rates is usually characterized by lower effective viscosity, i.e. the fluid is shear thinning. On the other hand, vorticity banding often corresponds to increased viscosity—shear thickening fluids. However, there have been reported experiments where such distinction was not valid Olmsted .

Since it is difficult to characterize gradient/vorticity banding in the terms of horizontal/vertical S-curves or shear thinning/thickening fluids, we shall discuss the S-curves themselves. There are several models leading to horizontal S-curves for the stress–strain dependency (multiple shear rates for a single shear stress). For instance, the Johnson–Segalman model Johnson-Segalman or rather its weakly non-local variant Olmsted-Radulescu , thermodynamics with internal variables Verhas ; Asszonyi2015 , viscoelastic models undergoing scission of polymeric chains Miroslav-micelles or viscoelastic models with relaxation of internal structure Miroslav-suspensions ; Germann-Beris indicate the horizontal S-shape of the stress–strain relation. The literature is rather vast, and we do not attempt to make the list of references complete. Interested reader could consult, for example, recent special issue of the Journal of Rheology Fielding-special . Most of the models usually exhibit shear thinning, and it is believed that shear thickening should be captured by models with non-linear dependence of the constitutive relation on the shear stress Radulescu , which is a motivation for the constitutive relation (1). This paper thus aims at modeling of the vertical S-shaped stress–strain curve with shear thickening, typically observed in vorticity banding.

3 Conjugate representation

The purpose of this section is to present details on the Legendre duality between $\Xi$ and $\Xi^{*}$. A consequence of the non-convexity of $\Xi^{*}$ is that $\Xi$ is multivalued, i.e. there are multiple shear stresses corresponding to certain shear rates. The shear rates that correspond to multiple stresses then constitute the CR-metastable and CR-unstable regions of the dynamics, which are discussed in Sec. 5.1. In fact, Sec. 5.1 contains the essence of this paper, and reading it before the following more technical parts could be advisable. Note also that Legendre transformation can be regarded as constrained extremization of the generating potential, which often turns out to be a practical method, see e.g. Junker2014 .

The balance of linear momentum reads

[TABLE]

where $\mathbf{u}$ is the density of momentum (per volume). The irreversible part of this evolution equation is in the form of divergence of the shear stress ${{{\mathbb{T}}}_{\delta}}$ (thermodynamic flux). Shear rate $\mathbb{D}$ is thus interpreted as a thermodynamic force. The reversible terms, see Appendix B for the definition of reversibility and irreversibility, are the gradient of the pressure and the divergence of convective terms, and their Hamiltonian nature is discussed in App. C.

Only the field of momentum density plays the role of state variable, since we assume that the fluid is incompressible ($\operatorname{div}\mathbf{v}=0$) and isothermal. When carrying out the CR-extension below (Sec. 4.2), also the field of conjugate stress tensor becomes a state variable.

[TABLE]

[TABLE]

The graph $\tilde{{\mathbb{T}}}_{\delta}(\mathbb{D})$ or $\tilde{\mathbb{D}}({{{\mathbb{T}}}_{\delta}})$ is plotted in Fig. 4, where also the lever rule is indicated. Of course, on one hand, we have

[TABLE]

while on the other hand from Fig. 2 it is clear that

[TABLE]

Combining these two results leads to

[TABLE]

which means that the shaded areas in Fig. 4 are equal—the Maxwell lever rule known from equilibrium thermodynamics Callen .

The Maxwell lever rule plays an analogical role in equilibrium phase transitions. The hysteresis is observed if fluctuations are sufficiently small. If they are large enough, the system evolves along the straight Maxwell line. The same is true in the case of the non-convex dissipation potential.

The Legendre duality between $\Xi$ as a function of forces and $\Xi^{*}$ as a function of fluxes is often advantageous because one of the forms can be more accessible. The former choice is used for example in chemical kinetics Grchemkin while the latter form for example used in plasticity modeling Mielke2003 .

Stability of constitutive relation: At many places within this paper, we discuss CR-stability of constitutive relations generated by non-convex dissipation potentials. By CR-stability of constitutive relations we mean whether the constitutive relation is CR-stable with respect to perturbations. If the dissipation potential is convex everywhere, the dual (Legendre transformed) dissipation potential is convex as well, and for each value of the thermodynamic force there is only one value of the flux. On the other hand, if the potential is non-convex, the dual dissipation potential is multi-valued and for some values of thermodynamic forces (or fluxes) there are several possible values of fluxes (or forces). The graph of the constitutive relation then has several branches as in Fig. 3. In such regions where multiple admissible fluxes (or forces) are present, perturbation of the flux can cause evolution of the system towards another branch of the graph of the constitutive relation. The constitutive relation is said to be CR-unstable if any perturbation leads to selection of a different branch, CR-metastable if only large enough perturbation switches the system to a different branch, and CR-stable if the current branch of the constitutive relation is not changed by any perturbation of the flux.

4 Extended hydrodynamics (CR-extension)

In this section the set of hydrodynamic state variables555The state variables are $(\rho,\mathbf{u},s)$ in the energetic representation or $(\rho,\mathbf{u},e)$ in the entropic representation. is expanded by promoting shear stress to an independent state variable, which is the reason for the name extended hydrodynamics. The idea of extending the set of state variables is ubiquitous in non-equilibrium thermodynamics, for example Extended Irreversible Thermodynamics (EIT) Jou-EIT or thermodynamics with internal variables Van-Berezovski . For simplicity, only the isothermal setting is considered in this paper. The isothermal hydrodynamic evolution is thus enriched to a form of extended hydrodynamics. The extra state variable has its own evolution equation, which will be shown to be the derivative of a multiscale thermodynamic Lyapunov (MTL) function with respect to the conjugate of the extra state variable (Sec. 4.5). The extra state variable then evolves towards the maxima of the Lyapunov function. The figures in Sec. 5.1 show how the Lyapunov function depends on the shear rate and towards what values the shear stress evolves. This will give a graphical representation of the CR-stability analysis of the system. Motivation for reading this rather technical section can thus be found in Sec. 5.1.

Dissipation potential $\Xi^{*}$ is convex near equilibrium and for high shear stresses, but it loses convexity in the intermediate range of the shear stress. Convexity leads to existence of a Lyapunov functional, see e.g. GO , and Braun–Le Chatelier principle Braun-LeChat , that guarantees stability of the evolution. The former paper is crucial for understanding everything that follows. Instability can be thus anticipated in the non-convex regions. Let us now have a look at a non-equilibrium thermodynamic description of stability, metastability and instability. The analysis is formulated within a recent framework developed by Grmela MG-CR , the CR-thermodynamics (CR stands for constitutive relation), which is a form of mesoscopic multiscale thermodynamics.666We use a slightly different notation compared with the original paper MG-CR . Instead of thermodynamic potentials, we use the negatives of the potentials, and we call them MTL functions. The reason is that we can then always use the same form of Legendre transformation. Another difference is that the extra variable is denoted by ${\mathbb{T}}_{\delta}^{*}$ instead of ${{{\mathbb{T}}}_{\delta}}$ so that the usual meaning of flux is kept.

4.1 GENERIC structure of hydrodynamics

Evolution of momentum density is governed by Eq. (8), which can be rewritten as (within the GENERIC framework)

[TABLE]

The reversible part is Hamiltonian, see in Appendix C, and it expresses the Euler equations of compressible hydrodynamics. The MTL function (the negative of the thermodynamic potential from MG-CR ) driving classical hydrodynamics (CH) to equilibrium (E) is

[TABLE]

where $e$ is total energy density (per volume). Total entropy is

[TABLE]

where $s(\rho,\varepsilon)$ denotes local equilibrium entropy density (which depends on density and internal energy density $\varepsilon=e-\frac{\mathbf{u}^{2}}{2\rho}$). Total energy is $E=\int\mathrm{d}\mathbf{r}\,e(\mathbf{r})$, and $M=\int\mathrm{d}\mathbf{r}\,\rho$ is the total mass. Constant $T_{0}$ is the temperature an isolated system would have after relaxation to equilibrium and constant $\mu_{0}$ the equilibrium chemical potential. Local temperature is defined as the derivative of entropy with respect to the total energy density, i.e. $T^{-1}(\mathbf{r})=_{\mathrm{def}}S_{e(\mathbf{r})}$. Derivative of the MTL function with respect to momentum is then

[TABLE]

Assuming isothermal conditions, the conjugate momentum, $\mathbf{u}^{*}$, is proportional to velocity $\mathbf{v}=\mathbf{u}/\rho$.

The dissipation potential is often in the form

[TABLE]

$\Gamma$ being a linear operator. In our case

[TABLE]

which leads to shear rate $\Gamma(\mathbf{u}^{*})_{\mathbf{u}^{*}=\Psi_{\mathbf{u}}^{(CH\rightarrow E)}}=\frac{1}{T_{0}}\mathbb{D}$. The constant prefactor $T_{0}$ is going to be ignored for the sake of simplicity. In a slightly overloaded notation, we can then write

[TABLE]

where the adjoint operator $\Gamma^{\top}$ is defined by

[TABLE]

This also specifies the scalar product $\langle\bullet,\bullet\rangle$ of two functions as the spatial integral of product of the functions ($L^{2}$ scalar product). In our case, we have $\Gamma^{\top}(\bullet)=\operatorname{div}(\bullet)$. Eventually, the evolution equation for momentum density becomes

[TABLE]

where ${{{\mathbb{T}}}_{\delta}}=\Xi_{\mathbb{D}}$. Note that the shear stress ${{{\mathbb{T}}}_{\delta}}$ is a thermodynamic flux because it is conjugate to the force ($\mathbb{D}$), which comes from differentiation of the MTL function (or entropy). Equation (22), which is the standard evolution equation for momentum density in fluid dynamics, is thus fully compatible with the GENERIC expression (14).

4.2 Promoting conjugate stress to a new state variable

The goal is now to promote the flux to an independent state variable. But ${{{\mathbb{T}}}_{\delta}}$ can not be promoted itself because derivative of a functional with respect to ${{{\mathbb{T}}}_{\delta}}$ would then have to be on right hand side of evolution equation for $\mathbf{u}$ instead of ${{{\mathbb{T}}}_{\delta}}$ itself. Rather, a conjugate flux ${\mathbb{T}}_{\delta}^{*}$ is the new variable, that is coupled with momentum density. Since ${{{\mathbb{T}}}_{\delta}}$ is interpreted the shear stress, ${\mathbb{T}}_{\delta}^{*}$ is interpreted as conjugate shear stress. As usual in non-equilibrium thermodynamics RedExt , the dissipation can be moved onto a higher (more detailed) level of description, i.e. into the evolution equation for the extra state variable. Dissipative evolution of momentum density can be then recovered by letting the extra variable relax. The reversible coupling between the original state variables and the extra state variable is required to be antisymmetric in order to fulfill the Onsager–Casimir reciprocal relations Casimir1945 ; PRE2014 ; HCO and thus we write (still following MG-CR ) evolution equations for extended hydrodynamics (EH) as follows

[TABLE]

To close the evolution equations, the conjugates are evaluated as the corresponding derivatives of a new extended MTL function $\Psi^{(EH\to E)}$, $\mathbf{u}^{*}=\Psi^{(EH\rightarrow E)}_{\mathbf{u}}$ and ${\mathbb{T}}_{\delta}^{**}={{{\mathbb{T}}}_{\delta}}=\Psi^{(EH\rightarrow E)}_{{\mathbb{T}}_{\delta}^{*}}$. Note that ${\mathbb{T}}_{\delta}^{**}$ was identified with ${{{\mathbb{T}}}_{\delta}}$ because the generalized Legendre transformation is idempotent (second application of the transformation restores the original functional), i.e. ${\mathbb{T}}_{\delta}^{**}={{{\mathbb{T}}}_{\delta}}$. The MTL function $\Psi^{(EH\rightarrow E)}$, which drives evolution of the extended hydrodynamics towards equilibrium, has yet to be specified in order to define the conjugate stress tensor completely. Evolution equations (23) represent evolution on the level of extended hydrodynamics.

Parity with respect to the time-reversal PRE2014 of ${{{\mathbb{T}}}_{\delta}}$ is now even, which is in contrast to the odd parity of ${{{\mathbb{T}}}_{\delta}}$ before promoting ${\mathbb{T}}_{\delta}^{*}$ to a state variable. The physical reason is that we have now a more detailed description of the shear stress, which can be regarded for example as the stress formulated by kinetic theory (App. D). Similarly, ${\mathbb{T}}_{\delta}^{*}$ and ${\mathbb{T}}_{\delta}^{**}$ are even with respect to time reversal, since MTL functions and dissipation potentials, with respect to which the fluxes are conjugate, are always even. Therefore, evolution equation (23a) is now reversible (in contrast to Eq. (8)) while equation (23b) contains a reversible term (first) and an irreversible term (second). By moving to a more detailed level of description, parity of the quantity promoted to an extra state variable may change because the quantity becomes independent and is no longer enslaved by the constitutive relation valid before the promotion. As a consequence, irreversible terms may become reversible.

4.3 Specification of the MTL function

The conjugate stress will be now specified more precisely. The relation between ${{{\mathbb{T}}}_{\delta}}$ and ${\mathbb{T}}_{\delta}^{*}$ is a Legendre transformation,

[TABLE]

which is carried out with respect to the MTL function

[TABLE]

This MTL function drives the extended hydrodynamics towards equilibrium. Note that all quantities with $(EH\rightarrow E)$ superscript depend on the hydrodynamic fields $\mathbf{x}^{(CH)}=_{\mathrm{def}}(\rho,\mathbf{u},e)$ and the conjugate stress tensor ${\mathbb{T}}_{\delta}^{*}$. Quantities $E^{(EH\rightarrow E)}$, $S^{(EH\rightarrow E)}$ and $M^{(EH\rightarrow E)}$, which have so far been left undetermined, are the building blocks of the corresponding MTL function, and they can reflect the particular nature of the physical system. Alternatively, the extra evolution equation can be interpreted as a numerical scheme leading to constitutive relation (1).

Let us consider at least one possible choice of $\Psi^{(EH\rightarrow E)}$. Within the kinetic theory (App. D), kinetic energy of particles is equal to the sum of kinetic energy of the overall motion $\mathbf{u}^{2}/(2\rho)$ and pressure, and the deviatoric part of stress tensor thus does not contribute to energy (when neglecting long-range pair interactions). On the other hand, additional knowledge expressed as the extra state variable reduces entropy, and the entropy on the $(EH)$ level is thus lower than the entropy on the level of classical hydrodynamics $(CH)$. Similarly as in Extended Irreversible Thermodynamics Jou-EIT , the first approximation of the $(EH)$ entropy is

[TABLE]

Coefficient $\beta$ is positive so that the entropy is concave. Energy and mass on the $(EH)$ level are the same as on the $(CH)$ level, thus the MTL function driving evolution from $(EH)$ to equilibrium is

[TABLE]

from which it follows that

[TABLE]

In summary, for the particular choice of the entropy (26), the conjugate stress tensor ${\mathbb{T}}_{\delta}^{*}$ is proportional to the standard shear stress tensor ${{{\mathbb{T}}}_{\delta}}$ and the evolution equations (23) for the constitutive relation (1) and the operator $\Gamma$ defined by (19) read

[TABLE]

The second equation can be apparently interpreted as the time derivative regularization of the original constitutive relation.

4.4 Duality between ${{{\mathbb{T}}}_{\delta}}$ and ${\mathbb{T}}_{\delta}^{*}$

The inverse Legendre transformation from ${\mathbb{T}}_{\delta}^{*}$ to ${{{\mathbb{T}}}_{\delta}}$ follows. The conjugate MTL function $\Psi^{*,(EH\rightarrow E)}$ can be constructed as

[TABLE]

where ${\mathbb{T}}_{\delta}^{*}({{{\mathbb{T}}}_{\delta}})$ is the solution to (24), i.e. Eq. (28) for the particular choices of entropy and energy. Conversely, the Legendre transformation

[TABLE]

leads back to the original MTL function,

[TABLE]

where ${{{\mathbb{T}}}_{\delta}}({\mathbb{T}}_{\delta}^{*})$ are the solutions to Eq. (31). The Legendre duality between ${{{\mathbb{T}}}_{\delta}}$ and ${\mathbb{T}}_{\delta}^{*}$ has thus been established.

4.5 Fast reducing evolution

Now we are coming to the highlight of the thermodynamic extension. Let us now identify the evolution of ${\mathbb{T}}_{\delta}^{*}$ that leads back to classical hydrodynamics, i.e. the reducing (fast) evolution leading to the reduced (hydrodynamic) evolution of variables $\mathbf{x}$. Equation (23b) can be rewritten as

[TABLE]

where

[TABLE]

with $\mathbb{D}=\Gamma(\mathbf{u}^{*})$, is the MTL function driving evolution from extended hydrodynamics to classical hydrodynamics. This MTL function will be referred to as the reducing MTL function. Equation (33) then means that conjugate shear stress ${\mathbb{T}}_{\delta}^{*}$ evolves in the direction of steepest ascent of $\Psi^{(EH\rightarrow CH)}(\mathbb{D},{{{\mathbb{T}}}_{\delta}})$. MTL function (34) generates Legendre transformation from $\Xi^{*}$ to $\Xi$ as

[TABLE]

which is equivalent to Eq. (9b). The constitutive relation is thus recovered after the extra variable (${\mathbb{T}}_{\delta}^{*}$) has relaxed to a state enslaved by the slow variable ($\mathbf{u}$).

The reducing evolution (33) can be alternatively rewritten as

[TABLE]

where Eq. (24) was used, and compatibility of the reducing evolution (33) with the second law of thermodynamics is guaranteed by

[TABLE]

for $\Psi^{(EH\rightarrow E)}$ concave. Note that only derivatives with respect to the fast variable are taken into account because only evolution of the fast variable is considered in the reducing evolution. The MTL function $\Psi^{(EH\rightarrow E)}$ has to be thus chosen concave, and the reducing evolution is compatible with the second law of thermodynamics.

On the other hand, dissipation potentials $\Xi$ and $\Xi^{*}$ were not convex and thus

[TABLE]

does not need to have definite sign and evolution on the hydrodynamic level can have apparent negative entropy production. Indeed, if shear rate and shear stress were spatially homogeneous with only for example $xy$ component nonzero, the slope of $\Xi^{*}$ could be negative for positive $\mathrm{D}_{xy}$. The way out of such a problem is the lift to extended hydrodynamics, where entropy production keeps its positivity. In other words, apparent negative entropy production would indicate evolution of instability (phase transition), and can be regarded as an invitation to a higher (more detailed) level of description.

After having demonstrated how the constitutive relation can be interpreted as a consequence of a more detailed evolution, study of behavior of the more detailed evolution gives information about CR-stability of the evolution with respect to the constitutive relation itself. The constitutive relation is stable if any perturbation decreases in time. It is metastable if only small perturbations decrease in time, and it is unstable if no perturbation decreases in time. The relation between CR-stability and the usual mathematical concept of stability is discussed in 5.2.

4.6 Alternative forms of the extension

Evolution equations (23) are relatively simple and useful for qualitative analysis of CR-stability carried out in Sec. 5.1, but they are not objective due to the presence of partial time derivative of the conjugate stress instead of an objective derivative. Moreover, Jacobi identity is not fulfilled for the system of equations (checked by program Kroeger2010 ). Hence, the evolution equations have to be regarded rather as a toy model useful for explaining construction of the extension and the CR-stability implications.

The extension can be made more physically sound by choosing a better operator $\Gamma$. Such a choice, however, must be accompanied by an another choice of the quantity to be promoted to an extra state variable. For example, the choice

[TABLE]

with $\mathbf{c}$ being the new variable (conformation tensor) leads to the adjoint operator

[TABLE]

It should be stressed that this extension is considered in the energetic representation, where entropy density is among the state variables instead of the total energy density. Let us refer to this extension as to the conformation tensor (CT) level.

The evolution equation for the extra variable is then

[TABLE]

which is analogical to Eq. (23). The reversible part is the upper-convected derivative of $\mathbf{c}$ (with an extra compressibility term, which is sometimes referred as to the Truesdell rate). Indeed, in the energetic representation total entropy depends only on the entropy density and derivative of MTL function (25) becomes

[TABLE]

and Eq. (41) then becomes

[TABLE]

where the left hand side is the upper-convective derivative of $\mathbf{c}$. Note that velocity is given by $\mathbf{u}/\rho$. Such an extended evolution also fulfills the Jacobi identity, is thus Hamiltonian, see e.g. Gn1 , and is also objective. Objectivity seems to be in tight relation with validity of Jacobi identity. Evolution equation (41) is the analogy of Eq. (23) and this whole paragraph is the analogy of Sec. 4.2.

Conjugate dissipation potential $\Xi^{*}$ is defined through Legendre transformation

[TABLE]

which is analogical to the transformations carried out in Sec. 3.

The relation between $\mathbf{c}$ and $\mathbf{c}^{*}$ can be obtained as in Sec. 4.4, i.e. through Legendre transformation with respect to a MTL function $\Psi^{(CT\rightarrow E)}$. The MTL function is constructed from entropy, energy and total mass (as in Eq. (25)), see e.g. Sarti-Marrucci ; Gn1 . A particular choice of entropy could be

[TABLE]

where $H$, $Q$ and $\beta$ are some material-dependent parameters and the last term expresses weakly non-local effects as in Vitek-stress-diffusion . The weakly non-local effects could play a role when the conformation tensor should exhibit steep spatial variations, as is usual in shear banding. Entropy (45) is just one possible choice made here for being explicit. However, the concrete form of entropy depends on the class of materials considered.

Finally, Eq. (41) can be rewritten as

[TABLE]

where the reducing MTL function is

[TABLE]

This is completely analogical to Sec. 4.5, and analysis of the second law of thermodynamics could proceed as in that section.

This alternative extension brings this paper closer to the common modeling of shear banding, where viscoelastic models are usually employed (see the discussion in Sec. 2). Similarly, by choosing a different operator $\Gamma$, one can derive evolution equations for Reynolds stress or weakly non-local vorticity, that can be found in Miroslav-turbulence ; PhysicaD-2015 . In the following section, however, we stick to the extension where conjugate stress tensor is the extra state variable for simplicity.

5 Interpretation of CR-stability

5.1 Thermodynamic interpretation

The reducing MTL function $\Psi^{(EH\rightarrow CH)}$ is plotted in Fig. 5 for several values of shear rates $\mathbb{D}$.

The curves plotted in the figure correspond to the following different regimes of the reducing evolution, indicated also in Fig. 3.

1. I

For low shear rates $\left|\mathbb{D}\right|$ (Regime I when comparing to Fig. 3) there is only one local maximum, which is near the origin ($\left|{{{\mathbb{T}}}_{\delta}}\right|=0\,\mathrm{Pa}$), as the curve for $\left|\mathbb{D}\right|=0.05\,\mathrm{s^{-1}}$ indicates. Bearing in mind that $\Psi^{(EH\rightarrow E)}$ has to be chosen concave, Eq. (36) tells that shear stress evolves in the direction of derivative of the reducing MTL function. In other words, if the derivative is positive, ${{{\mathbb{T}}}_{\delta}}$ increases, and if it is negative, ${{{\mathbb{T}}}_{\delta}}$ decreases. Thus, if ${{{\mathbb{T}}}_{\delta}}$ is initially to the left of the local maximum, where the derivative is positive, it tends to the local maximum. If it is to the right of the maximum, where the derivative is negative, it again tends to the maximum. ${{{\mathbb{T}}}_{\delta}}$ is thus attracted to the local (and global) maximum. The evolution towards this maximum is thus CR-stable. 2. II

For slightly higher shear rates (Regime II), another local maximum appears in higher shear stresses (curve for $\left|\mathbb{D}\right|=0.12\,\mathrm{s^{-1}}$). Similarly as in Regime I (low shear rates), this new local maximum attracts ${{{\mathbb{T}}}_{\delta}}$ from its neighborhood. However, since the first maximum is still higher, settling of ${{{\mathbb{T}}}_{\delta}}$ in the second maximum can be seen as a CR-metastable state. After a sufficient fluctuation, ${{{\mathbb{T}}}_{\delta}}$ will end up in the first (higher) maximum. Moreover, a local minimum appears between the local maxima. The local minimum is an CR-unstable state, since ${{{\mathbb{T}}}_{\delta}}$ is attracted to the local maxima regardless how close to the local minimum it is. That can be seen from signs of derivative of the reducing MTL function to the left and to the right of the local minimum. 3. III

For even higher shear rates (Regime III and curve $\left|\mathbb{D}\right|=0.2\,\mathrm{s^{-1}}$), the second maximum becomes higher than the first one, and the first maximum thus becomes CR-metastable while the second becomes CR-stable. The local minimum is still located between the maxima, and although it is a stationary point ($\dot{{\mathbb{T}}}_{\delta}=\mathbb{0}$ in the local minimum), such a state is CR-unstable, since any arbitrarily small fluctuation in the stress would make the stress go towards one of the maxima. Note that the value of the shear rate between Regimes II and III, when the two local maxima have the same height, corresponds to the shear rate selected by the Maxwell lever rule (see Sec. 2). 4. IV

Finally, (Regime IV) the first local maximum disappears (as well as the local minimum) for even higher shear rates (curve $\left|\mathbb{D}\right|=0.35\,\mathrm{s^{-1}}$), and the second maximum becomes the only CR-stable state.

The described dynamics of the stress tensor can be seen also in Fig. 3. Regime I (low shear rates) corresponds to the part of curve $A-B$ left of point $D$, where only one ${{{\mathbb{T}}}_{\delta}}=\Xi_{\mathbb{D}}$ is possible. Regime II (one CR-stable state, one CR-unstable state and one CR-metastable state) is indicated in that figure as well. The upper solid curve is the CR-stable state, the lower solid curve is the CR-metastable state and the dashed curve is the CR-unstable state. Regime III (one CR-metastable state, one CR-unstable state and one CR-stable state) is described analogically. Finally, Regime IV has only one state, which is CR-stable. The reducing evolution (36) is thus able to distinguish between CR-stable, CR-metastable and CR-instable states in Fig. 3, and the analysis is similar to the treatment of phase transitions in equilibrium thermodynamics Landau5 ; Callen . We are thus coming to a kind of dissipative phase transitions.

Let us now return to the experimental results in Fig. 3 and to the second question raised in the introduction. When slowly raising the shear rate, we start in Regime I, where shear stress is determined uniquely. After passing to Regime II, an another possible state appears, but shear stress evolves continuously due to carefulness of the experimentalists (smallness of fluctuations). Passing also through Regime III, point $C$ (boundary between Regimes III and IV) is reached, where the CR-metastable state becomes CR-unstable and shear stress jumps to the CR-stable state. That is why shear stress behaves discontinuously.

Moreover, it could be expected that when going back from high shear rates to the small ones, shear stress would fall as low as to point $D$ in Fig. 1, from which it would jump to the CR-stable state on branch $A$. One should thus observe hysteresis indicated in that figure. We are not aware of such a backward experiment having been carried out, but modeling of hysteresis by means of introducing a new evolution equation is usual in literature kubin.lp.poirier.jp:relaxation .

On the other hand, when slowly varying ${{{\mathbb{T}}}_{\delta}}$, there is always a unique value of $\mathbb{D}$, determined by slope of dissipation potential $\Xi^{*}$ in Fig. 2, and the slope varies continuously. That is why no jump in shear rate is observed in the experiment.

5.2 Stability of the regularized constitutive relation

The same dynamical behavior can be concluded directly from the stability analysis of the regularized constitutive relation (29b). Its one-dimensional counterpart with a general conjugate dissipation potential $\Xi^{*}$ reads

[TABLE]

In the shear-rate-controlled experiment we assume that $\left|\mathbb{D}\right|=\left|\mathbb{D}\right|(t)$ is given. Equation (48) is then a non-autonomous ordinary differential equation with stationary points $\left|\mathbb{D}\right|-\Xi^{*}_{\left|{{{\mathbb{T}}}_{\delta}}\right|}\left(\left|{{{\mathbb{T}}}_{\delta}}\right|\right)$. It can equivalently written as an extended autonomous system

[TABLE]

From this system, we can deduce the phase field in the $\left|\mathbb{D}\right|$–$\left|{{{\mathbb{T}}}_{\delta}}\right|$ plane. In Figure 6, we show the phase field for the particular choice of the dual dissipation potential (4).

We can easily see that the solutions are quickly (depending on the parameter $\beta$) attracted towards the increasing parts of the constitutive relation and repelled from the decreasing part. In other words, the decreasing part is unstable in the usual mathematical sense while the increasing parts are asymptotically stable. This leads to the jumps (and possibly hysteresis) in the values of the stress $\left|{{{\mathbb{T}}}_{\delta}}\right|$ when changing $\left|\mathbb{D}\right|$.

On the other hand, if we prescribe $\left|{{{\mathbb{T}}}_{\delta}}\right|=\left|{{{\mathbb{T}}}_{\delta}}\right|(t)$ (shear stress controlled experiment) and we further assume that its time derivative is small enough, (48) reduces back to the original constitutive relation. In this case, there are no jumps.

At least in this case the CR-stability analysis leads to the same conclusion as the usual stability analysis for ordinary differential equations provided the evolution of the extra state variable (here conjugate stress) is assumed much faster than evolution of the remaining variables.

5.3 Linearized hydrodynamic stability

The CR-instability of the decreasing part of the constitutive relation can be also verified by means of the linearized hydrodynamic stability. Perturbing a known simple shear base flow, the aim is to study the temporal evolution of the disturbance. Linearizing the governing equations for the disturbance yields a generalized eigenvalue problem. Since constitutive relation (1) does not have to be invertible, the resulting problem might be in the form of a system of equations, unlike its classical counterpart—the single Orr–Sommerfeld equation, see Orr orr.w:stability*1 ; orr.w:stability and Sommerfeld sommerfeld.a:beitrag . Using a variant of a pseudospectral collocation method, e.g. Canuto and Hussaini canuto.c.hussaini.my.ea:spectral*1 ; canuto.c.hussaini.my.ea:spectral , we are able to solve the problem and numerically confirm that once the flow is governed by a non-monotone constitutive relation, the perturbation grows in time and the flow is unstable. The derivation of the problem and description of the numerical method is rather lengthy and it will be subject of a future publication.

The hydrodynamic stability is usually studied for a Poiseuille or Couette flow of a fluid in some archetypal geometry like a two-dimensional channel or a pipe with prescribed boundary conditions. In the classical setting, this yields a single differential eigenvalue equation, while in our case, it leads to a system of equations. In either case, it has to be solved numerically and the spectral collocation method has proven effective, see Orszag . This is different from, for example, Ván Van-stability who could use the Routh–Hurwitz criteria, since he studied the stability of the equilibrium and he considered the perturbation in the form of an unbounded plane wave, resulting in a third-order polynomial. We cannot make use of the mentioned criteria as we want to study stability of non-equilibrium flows with imposed perturbations also satisfying the boundary conditions, and consequently we do not obtain just a polynomial. Anyway, determining the stability of the equilibrium can be also done from the kinetic energy of the disturbance.

6 Braun–Le Chatelier principle

Braun–Le Chatelier principle in equilibrium thermodynamics is the statement that a perturbation of a thermodynamic force acting on a system triggers such internal processes inside the system that the effect of the perturbation is reduced after the processes have relaxed LeChatelier . The principle can be derived from convexity of the energy functional Landau5 , and it is equivalent with stability of the system. Analogical derivation starting from convexity of a dissipation potential leads to an analogical result, Braun–Le Chatelier principle for dissipative thermodynamics Braun-LeChat ,

[TABLE]

where $J_{int}$ is a rate of change of an internal variable and $X_{int}$ is the thermodynamic force corresponding on the internal variable. In other words, apparent viscosity (the derivative of stress with respect to strain) is lower after the internal variable has relaxed ($J_{int}=0$) than just after the perturbation (change of shear rate).

Therefore, Braun–Le Chatelier principle, which is a consequence of relaxation of an internal variable and convexity of a dissipation potential, leads to shear thinning. This is in agreement with that shear thinning is usually observed in the gradient banding experiments, which are well described by a viscoelastic model. The internal variable dynamics of which leads to the viscoelastic behavior relaxes and causes shear thinning in agreement with the Braun–Le Chatelier principle.

On the other hand, shear thickening is more difficult to explain, which supports the use of non-linear or even non-convex dissipation potentials.

7 Non-smooth dissipation potentials

This section stands a little apart from the rest of the paper, as it addresses yield stress behavior instead of S-curves. It is included because the authors believe that discussing non-smoothness of dissipation potentials is a natural continuation of the discussion of non-convexity.

Inspired by Rajagopal2004 , dissipation potential exhibiting yield stress (no shear rate for sufficiently low shear stress) is for example

[TABLE]

where $\alpha$ and $\beta$ are positive constants, see Fig. 8. The yield stress (given by the coefficient $\alpha$) is caused by the fact that for low shear stresses (low slopes of the dissipation potential) there is no corresponding $\mathbb{D}$. More exactly, the only shear rate corresponding to low shear stresses is zero.

Dissipation potential (51) is not smooth at $\left|\mathbb{D}\right|=0$, where it has no derivative. Therefore, Legendre transformation to $\Xi$ has to be generalized by means of subdifferentials, see e.g. Roubicek ,

[TABLE]

where $\partial\Xi$ is the set of slopes of all hyperplanes touching $\Xi$ at a point $\mathbb{D}$. When constructing the subdifferential we assume that $\Xi$ is convex and continuous in a neighborhood of the point of non-smoothness. Continuity then guarantees that the set is non-empty (a consequence of Hahn–Banach theorem as commented in Roubicek ). Non-smoothness in non-convex regions can be handled as in Remark 5.8 in Roubicek . When $\Xi$ is smooth, i.e. it has standard functional derivative, inclusion (52) becomes equality and the Legendre transformation restores its usual form.

Inclusion (52) can be interpreted as the following problem: For a given ${{{\mathbb{T}}}_{\delta}}$ find all values of $\mathbb{D}$ such that the inclusion is fulfilled. Solution to the problem is a graph $\mathbb{D}({{{\mathbb{T}}}_{\delta}})$. It can be seen in Fig. 8 that for sufficiently high shear stresses, there is always a non-zero $\mathbb{D}$ for which $\Xi_{\mathbb{D}}={{{\mathbb{T}}}_{\delta}}$. However, only $\partial\Xi$ at $\left|\mathbb{D}\right|=0$ contains the low shear stresses, which means that when the shear stress is sufficiently low, only $\left|\mathbb{D}\right|=0$ is the solution to the problem. In other words, there is no shear rate for sufficiently low shear stresses, but when shear stress crosses a threshold, shear rates starts to increase. We will return to this point graphically later on.

The conjugate dissipation potential is defined as

[TABLE]

and it is depicted in Fig. 8. Since the conjugate dissipation potential is smooth, the backward Legendre transformation is carried out with derivatives as usually.

Let us now have a look at these dissipation potentials from the perspective of multiscale mesoscopic thermodynamics as in Sec. 4. The reducing MTL function is again $-\Xi^{*}({{{\mathbb{T}}}_{\delta}})+\int\mathrm{d}\mathbf{r}\,{{{\mathbb{T}}}_{\delta}}:\mathbb{D}$ and it is plotted in Fig. 9.

The shear stress again tends to the maximum of the MTL function for the given shear rate. When $\left|\mathbb{D}\right|=0$, norm of shear stress can end up anywhere between 0 and 1 as the function is flat in that region. When shear rate is higher, the maximum is located to the right of $\left|{{{\mathbb{T}}}_{\delta}}\right|=1$. In other words, if there is non-zero shear rate, shear stress has to be higher than the threshold $\left|{{{\mathbb{T}}}_{\delta}}\right|=1$. The material thus exhibits yield stress behavior and its irreversible evolution is CR-stable.

8 Critical heat flux

In 1934, Nukiyama nukiyama.s:maximum experimentally showed that by controlling the heat flux from a hot platinum wire to boiling water, the temperature of the wire jumps between the heating branch and the cooling branch and exhibits a hysteresis loop. Then in 1937, Drew and Mueller drew.tb.mueller.ac:boiling were able to control the temperature and observed that the heat flux initially grows with the temperature, then decreases until it finally grows again. The boiling curve, heat flux versus temperature difference (thermodynamic force) dependence, is then non-monotonous, see also Stosic Stosic and the book by Bergman and co-authors bergman.tl:fundamentals . In Fig. 10, we have plotted the experimental data by Nukiyama nukiyama.s:maximum and fitted them with a relation of the form (1).

As the flux-force relation can be described with an analogue of (1), it can be likewise generated by a non-convex dissipation potential. Such a formulation would then lead to implications on the CR-stability of the constitutive relation and the experimentally observed hysteresis. For instance, the CR-unstable part is inaccessible when controlling the heat flux while careful temperature variation makes it available. This is analogous to the stress/shear rate behavior described in this paper. We hope to address this issue more explicitly in the near future.

9 Conclusion

Constitutive relation (1), that fits well the considered experimental data (Fig. 1), can be formulated by means of dissipation potential $\Xi^{*}({{{\mathbb{T}}}_{\delta}})$, defined in Eq. (4). This dissipation potential is convex near equilibrium, but loses convexity for higher shear stresses. Due to the loss of convexity the Legendre-conjugate dissipation potential $\Xi(\mathbb{D})$ is multivalued. This reflects the experimental observation that varying shear stress leads to continuous evolution of shear rate (value of $\mathbb{D}=\Xi^{*}_{{{{\mathbb{T}}}_{\delta}}}$ is unique) while varying shear rate leads to a jump (and possibly hysteresis) in shear stress (value of ${{{\mathbb{T}}}_{\delta}}=\Xi_{\mathbb{D}}$ is not unique). The phenomenon of critical heat flux can be addressed analogically.

By extending classical hydrodynamics to extended hydrodynamics as in Sec. 4, one can show which parts of the multivalued dissipation potential $\Xi$ are CR-stable (with respect to perturbations of the constitutive relation), which are CR-metastable and which CR-unstable. It is then possible to identify between which branches of the dissipation potential the shear stress will jump, see Sec. 5.1. The findings are compatible with the experimental data.

Acknowledgment

We are grateful to Vít Průša, Václav Klika and Onřej Souček for their invaluable comments when discussing the manuscript and for pointing out further research directions (e.g. hysteresis modeling in rheology). We also want to express our gratitude to Gian Paolo Beretta for drawing our attention to the Critical heat flux.

Most ideas in this work are either inspired or brought up by Miroslav Grmela, founder of our thermodynamic tribe.

This work was supported by Czech Science Foundation, project no. 17-15498Y. This work has been supported by Charles University Research program No. UNCE/SCI/023.

Adam Janečka acknowledges the support of Project No. LL1202 in the programme ERC-CZ funded by the Ministry of Education, Youth and Sports of the Czech Republic and the support of Project 260 449/2017 “Student research in the field of physics didactics and mathematical and computer modelling”.

Appendix A Generalizations of the Legendre transformation

For a smooth, strictly convex (with positive second derivative) real function $f(x)$, the Legendre transform (see Callen ) is defined as

[TABLE]

where $x$ satisfies $f^{\prime}(x)=u$. Since the function $f$ is strictly convex, its derivative is monotonic and thus invertible, i.e. $x=\left(f^{\prime}\right)^{-1}(u)$. Then we can write

[TABLE]

The Legendre transform of a convex function is concave which can be immediately seen from

[TABLE]

Zia et al. zia.rkp.ea:making advocated for a definition with the opposite sign

[TABLE]

which is a symmetric representation of the Legendre transform, $f^{*}_{L}(u)+f(x)=ux$, and converses convexity. Both versions are invertible and the bi-conjugate is the original function $f^{**}=f$. In the second version, the formula for the inverse transform is identical to the primal transform.

The generalization of the Legendre transform to arbitrary functions is the Fenchel transform (also called the Fenchel–Legendre transform or the Young–Fenchel conjugate)

[TABLE]

The Fenchel transform is not invertible, the inversion of the conjugate yields the convex hull of the original function

[TABLE]

thus some information is lost in the process.

The Fenchel transform can be defined in the same manner for functionals, see bauschke.hh.lucet.y:what ; Roubicek . Let $V$ be a real Hilbert space and $f:V\to{\mathbb{R}}$, then the Fenchel transform $f^{*}_{F}$ at $u\in V^{*}$ is

[TABLE]

where $\langle\bullet,\bullet\rangle$ is the corresponding inner product. The conjugate $f^{*}_{F}$ is always convex and lower semicontinuous and $f^{**}_{F}=f$ if and only if $f$ is convex and lower semicontinuous.

For arbitrary differentiable functions, dorst.l.boomgaard.r:analytical ; dorst.l.boomgaard.r:morphological introduced the invertible multivalued slope transform

[TABLE]

where $\operatorname{stat}_{x}f(x)$ is the set of all stationary values of $f(x)$ with respect to $x$

[TABLE]

Then the inverse slope transform is

[TABLE]

Throughout the paper, we use this general definition of the Legendre-like conjugation for it is invertible and no information is lost. For functions with consecutively variating convex and concave parts, each with an invertible derivative (as in our case), the slope transform corresponds to application of the classical Legendre transform to each part. Then, the number of different functions in the multivalued Legendre transform agrees to the convex/concave parts of the original function, see maragos.p:slope . For functionals in general, we believe that the slope transform should be invertible as well, unfortunately we could not find any relevant mathematical reference.

Appendix B Reversibility and irreversibility

Definition of reversible dynamics is based on the time-reversal transformation (TRT), see e.g. PRE2014 or HCO . TRT inverts momenta of all particles, $\mathbf{p}\stackrel{{\scriptstyle TRT}}{{\rightarrow}}-\mathbf{p}$ and the direction of time increment $\mathrm{d}t\stackrel{{\scriptstyle TRT}}{{\rightarrow}}-\mathrm{d}t^{\prime}$. If an evolution equation is transformed by TRT to itself, it is reversible. If, on the other hand, some terms on the right hand side transform differently than the left hand side (time derivative of a state variable), these terms generate irreversible evolution.

Let us illustrate the concept of reversibility on the Hamilton canonical equations, which are transformed by TRT as follows.

[TABLE]

After TRT, the Hamilton canonical equations are equivalent with the original equations, which means that they are reversible.

On the other hand, Hamilton canonical equations with damping, e.g.

[TABLE]

are not completely reversible because the damping term in the momentum equation transforms differently than the left hand side of the equation. The equations consist of both reversible and irreversible evolution.

Note also that the equations obtained from a variational principle (extremization of an action) are not automatically reversible although Hamilton canonical equations (without damping) are. Such equations can be found for instance in textbook Landau5 , where a dissipation function is introduced into the principle of least action.

On the other hand, Hamiltonian kinematics derived by projection from a more detailed reversible Hamiltonian kinematics is surely reversible, see PhysicaD-2015 , and this is the case the Boltzmann kinematics (kinetic theory) and fluid mechanics. When talking about reversible and irreversible terms in this paper, we mean the reversibility with respect to TRT and we refer to the cited literature in case of any doubts.

Appendix C Reversible part of momentum balance

The reversible part of the momentum balance (evolution equation for momentum density) is generated by the Poisson bracket of classical hydrodynamics (see e.g. Arnold ),

[TABLE]

Supplied with energy, e.g.

[TABLE]

the reversible evolution of an arbitrary functional of the hydrodynamic variables $A(\rho,\mathbf{u},s)$ is

[TABLE]

By rewriting time-derivative of the functional as

[TABLE]

and comparing with Eq. (68), one can read the reversible part of evolution equations for the hydrodynamic fields,

[TABLE]

where $p=-\varepsilon(\rho,s)+\rho\varepsilon_{\rho}+s\varepsilon_{s}$ is the pressure. The reversible evolution of hydrodynamic fields with energy (67) is of course the Euler equations. The reversible terms from Eq. (8) are on the right hand side of Eq. (70b).

Strictly speaking, Eq. (8) is in the entropic representation (in variables $(\rho,\mathbf{u},e)$) while Eq. (70b) is in the energetic representation (variables $(\rho,\mathbf{u},s)$) and one should thus perform transformation from the latter representation to the former, see e.g. HCO .

Appendix D Kinetic theory

Let us recall a few fundamental findings from classical kinetic theory, well reviewed in dGM . The state variable of kinetic theory is the one-particle distribution function $f(\mathbf{r},\mathbf{v})$, which depends on the position and the velocity. To express energy of interactions between particles exactly, one would need to use a two-particle distribution function, which is inaccessible in the one-particle kinetic theory. Therefore, let us stay in the region of ideal or nearly ideal (diluted) gas, where total energy is

[TABLE]

$m$ being mass of one particle. Density, momentum density and entropy density then are

[TABLE]

Evolution equation for the one-particle distribution function in absence of external force is the Boltzmann equation

[TABLE]

where the last term represents the collision integral, which generates irreversible evolution of $f$. Evolution equation for momentum density then becomes

[TABLE]

Note that the collision integral does not contribute to balance of momentum due to conservation of momentum in collisions. The term on the right-hand side generates reversible evolution of $\mathbf{u}$ on the level of kinetic theory, since it is even with respect to time-reversal, see e.g. PRE2014 , and it can be regarded as divergence of stress tensor

[TABLE]

The last term in Eq. (75) constitutes the shear stress, and the standard expression for Newtonian fluids can be recovered near equilibrium for example by means of Champan–Enskog expansion, see e.g. dGM ; MG-CR . The important observation here is that on the kinetic level of description the shear stress (or extra stress) is even with respect to time-reversal while on the hydrodynamic level, where it is proportional to shear rate, it is odd.