Finite Amplitude Stability of Internal Steady Flows of the Giesekus Viscoelastic Rate-Type Fluid

Abstract

Using a Lyapunov type functional constructed on the basis of thermodynamical arguments, we investigate the finite amplitude stability of internal steady flows of viscoelastic fluids described by the Giesekus model. Using the functional, we derive bounds on the Reynolds and the Weissenberg number that guarantee the unconditional asymptotic stability of the corresponding steady internal flow, wherein the distance between the steady flow field and the perturbed flow field is measured with the help of the Bures–Wasserstein distance between positive definite matrices. The application of the theoretical results is documented in the finite amplitude stability analysis of Taylor–Couette flow.

1. Introduction

Flows of viscoelastic fluids exhibit the phenomenon dubbed “elastic turbulence” or “inertia-less turbulence”. The flows of viscoelastic fluids can become—unlike the flows of the standard viscous fluids—unstable or “turbulent” at very low values of the Reynolds number. This behaviour indicates that the instability or transition to “turbulence” is driven by a nonstandard mechanism. Namely, it is not driven by the nonlinearity due to the inertial term in balance of linear momentum, but it must be attributed to the nonlinearity in the governing equation for the “elastic” part of the Cauchy stress tensor. The key challenge is to identify the parameter values that prohibit the onset of “elastic” instability or that trigger the “elastic” instability, and to describe the transition scenarios leading from the laminar to the “turbulent” flow. This task requires one to perform some sort of nonlinear stability analysis, since a nonlinear interaction between the finite amplitude perturbations might be decisive.

The phenomenon of “elastic turbulence” has been thoroughly investigated both from the experimental as well as theoretical point of view (see reviews by Petrie and Denn [1], Larson [2], Shaqfeh [3], Morozov and van Saarloos [4] or Li et al. [5]). In particular, the experiments reported by Groisman and Steinberg [6] stimulated enormous research activity regarding the elastic turbulence. On the other hand, theoretical results mainly follow from direct numerical simulations based on various viscoelastic rate-type models (see Dubief et al. [7], Lieu et al. [8], Grilli et al. [9], Page and Zaki [10], Biancofiore et al. [11], Valente et al. [12], Lee and Zaki [13] and Plan et al. [14] for some recent contributions). The need to resort to sophisticated numerical simulations in order to get qualitative insight into the flow dynamics is not surprising.

The reason is that the instabilities in viscoelastic fluids are very likely of subcritical nature (see Meulenbroek et al. [15]). The subcritical nature of the instability implies, as remarked by Morozov and van Saarloos [4], that

[Linear stability] (if it exists) is not very relevant for the existence of dynamics of the patterns that typically arise before the instability of the base state occurs.

(Note that a similar issue arises even for the standard Navier–Stokes fluid (see, for example, Baggett and Trefethen [16] for a discussion of several low-dimensional models of subcritical transition).) This means that the linear stability analysis, that is stability analysis with respect to infinitesimal perturbations, is of limited applicability in the investigation of the transition scenarios, albeit it can still provide important insight into the problem. (See, for example, Beris et al. [17], Blonce [18], Grillet et al. [19] and Pourjafar and Sadeghy [20], Pourjafar and Sadeghy [21] for linear stability analysis of flows of viscoelastic fluids described by the Giesekus viscoelastic rate-type model.) Moreover, quoting again Morozov and van Saarloos [4]

[Subcritical instability] is governed by all kinds of nonlinear self-enhancing interactions and so there is almost never a simple approximation scheme that allows one to explore the infinite dimensional space of interactions in all details, and determine which direction corresponds to the smallest threshold [for instability]. Thus, in practice, one can explore such situations, in theoretical studies as well as in experiments, only for a given class of perturbations.

On top of that, even if the technique such as weakly nonlinear analysis is apparently successful, then, as Meulenbroek et al. [15] put it,

One should also keep in mind that our expansion is only carried out to lowest order in the nonlinearity, so one may wonder about the robustness of these results as long as higher order terms in the expansion are unknown.

In what follows, we want to address the lack of analytical results for the stability problem of flows of viscoelastic fluids subject to finite amplitude perturbations. In particular, using a Lyapunov type technique, we investigate the stability of internal steady flows of viscoelastic fluids described by the Giesekus model, and we derive bounds on the values of the Reynolds number and the Weissenberg number that guarantee the flow stability subject to any (finite) perturbation. The result provides a sufficient condition for stability, hence it can be seen as complementary result to the search for the smallest threshold for instability via approximation methods. The derived bounds are interesting not only on their own. What is perhaps equally interesting is the way the bounds are derived. The derivation heavily relies on the underlying thermodynamical arguments and the notion of energy, which is an approach that seems to be discouraged in the nonlinear stability analysis of viscoelastic fluids (see Doering et al. [22]).

The paper is organised as follows. In Section 2, we describe the Giesekus model, and we briefly comment on its thermodynamical underpinnings. In particular, we identify the energy storage mechanisms and the entropy production mechanisms that are implied by the evolution equations for the Giesekus model. Once the thermodynamical background is summarised, we formulate the governing equations for an internal steady flow (see Section 3), and we proceed with the stability analysis of this non-equilibrium steady state. The stability is analysed using a Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ constructed by the thermodynamically based method proposed by Bulíček et al. [23]. The functional used in the stability analysis of a steady flow $\widehat{\mathit{v}}$ in a domain $\Omega$ is constructed in Section 4, and it is given by the formula

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})={\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}[-\ln \det (\mathbb{I}+{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})+\mathrm{Tr}\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)]\mathrm{d}\mathrm{v},$$

(1)

where $\tilde{\mathit{v}}$ denotes the perturbation of the velocity field, $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ is related to the perturbation of the stress field, and $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ is related to the stress field in the steady flow. (See the corresponding sections for the notation.) The fact that Equation (1) can serve as a Lyapunov type functional is closely related to the proper choice of the distance function that characterises the proximity of the perturbation $\mathit{W}$ and the corresponding steady state $\widehat{\mathit{W}}$. In our case, the distance function is introduced using the Bures–Wasserstein distance

$${\mathrm{dist}}_{\mathbb{P}\left(d\right),\mathrm{BW}}(\mathbb{A},\mathbb{B}){=}_{\mathrm{def}}{\{\mathrm{Tr}\mathbb{A}+\mathrm{Tr}\mathbb{B}-2\mathrm{Tr}\left[{\left({\mathbb{A}}^{\frac{1}{2}}\mathbb{B}{\mathbb{A}}^{\frac{1}{2}}\right)}^{\frac{1}{2}}\right]\}}^{\frac{1}{2}},$$

(2)

(see Bhatia et al. [24]), which measures the distance between the symmetric positive semidefinite matrices $\mathbb{A}$ and $\mathbb{B}$. Once we generalise Equation (2) to the setting of spatially distributed fields of symmetric positive semidefinite matrices, we exploit the concept of Lyapunov functional, and we derive bounds on the Reynolds number and the Weissenberg number that guarantee the flow stability with respect to any perturbation (see Theorem 1). These bounds are universal for any flow geometry.

The bounds on the Reynolds number and the Weissenberg number are then explicitly evaluated in Section 6 in the case of Taylor–Couette type flow. Further, we also perform direct numerical simulations that allow us to quantitatively document some features of the perturbation dynamics. The results are commented in Section 7.

2. Giesekus Model

2.1. Governing Equations

The governing equations for the Giesekus fluid (see Giesekus [25]), in the absence of external force, read

$$\begin{array}{cc}\hfill \mathrm{div}\mathit{v}& =0,\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \rho \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}& =\mathrm{div}\mathbb{T},\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\nu }_1\stackrel{▿}{\overline{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}& =-\mu [\alpha {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}^2+(1-2\alpha ){\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-(1-\alpha )\mathbb{I}],\hfill \end{array}$$

(3c)

where $\mathit{v}$ denotes the velocity, $\rho$ denotes the density and ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is an extra tensorial quantity whose physical meaning is given below. Finally, the symbol $\mathbb{T}$ denotes the Cauchy stress tensor that is given by the formulae

$$\mathbb{T}=m\mathbb{I}+{\mathbb{T}}_{\delta },{\mathbb{T}}_{\delta }=2\nu \mathbb{D}+\mu {\left({\mathbb{B}}_{{\kappa }_{p\left(t\right)}}\right)}_{\delta },$$

(3d)

where m denotes the mean normal stress (pressure) and $\mathbb{D}{=}_{\mathrm{def}}\frac{1}{2}(\nabla \mathit{v}+{(\nabla \mathit{v})}^{\top })$ denotes the symmetric part of the velocity gradient. Symbols $\nu$, ${\nu }_1$, $\mu$ and $\alpha$, $\alpha \in (0,1)$ denote material parameters. Note that, if $\alpha =0$, then one recovers the standard Maxwell/Oldroyd-B models. The value $\alpha =0$ is however not covered in the presented stability analysis.

The remaining notation is the standard one: $\frac{\mathrm{d}}{\mathrm{d}t}{=}_{\mathrm{def}}\frac{\partial }{\partial t}+\mathit{v}•\nabla$ denotes the material time derivative, and

$$\stackrel{▿}{\mathbb{A}}{=}_{\mathrm{def}}\frac{\mathrm{d}\mathbb{A}}{\mathrm{d}t}-\mathbb{L}\mathbb{A}-\mathbb{A}{\mathbb{L}}^{\top },$$

(4)

denotes the upper convected derivative, where $\mathbb{L}{=}_{\mathrm{def}}\nabla \mathit{v}$, and the symbol ${\mathbb{A}}_{\delta }{=}_{\mathrm{def}}\mathbb{A}-\frac{1}{3}\left(\mathrm{Tr}\mathbb{A}\right)\mathbb{I}$ denotes the traceless part of the corresponding tensor. In virtue of the incompressibility constraint in Equation (3a), we have ${\mathbb{D}}_{\delta }=\mathbb{D}$. Note that, if one uses a simple substitution $\mathbb{S}{=}_{\mathrm{def}}\mu ({\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-\mathbb{I})$, and if one redefines the pressure, $p{=}_{\mathrm{def}}-m+\frac{1}{3}\left(\mathrm{Tr}\mathbb{S}\right)\mathbb{I}$, then Equation (3c) and (3d) transform to $\lambda \stackrel{▿}{\overline{\mathbb{S}}}+\mathbb{S}+\frac{\alpha \lambda }{{\nu }_1}{\mathbb{S}}^2=2{\nu }_1\mathbb{D}$ and $\mathbb{T}=-p\mathbb{I}+2\nu \mathbb{D}+\mathbb{S}$, where $\lambda {=}_{\mathrm{def}}\frac{{\nu }_1}{\mu }$. This is another frequently used form of the governing equations for the Giesekus fluid.

2.2. Thermodynamic Basis

The Giesekus model has been originally derived without any reference to thermodynamics. However, we want to design a Lyapunov type functional using concepts from non-equilibrium thermodynamics, hence we need to explore thermodynamical underpinnings of the model. The issue of finding a thermodynamic basis for viscoelastic rate-type models is claimed to be resolved by a plethora of theories for thermodynamics of complex fluids (see, for example, Leonov [26], Mattos [27], Wapperom and Hulsen [28], Dressler et al. [29] or Ellero et al. [30]). (Notably, the treatise by Dressler et al. [29] contains a rich bibliography, and describes the issue from the viewpoint of the GENERIC formalism (see Grmela and Öttinger [31], Öttinger and Grmela [32] and Pavelka et al. [33])). In the present analysis, we exploit the approach proposed by Rajagopal and Srinivasa [34] that is relatively simple and that provides one a purely phenomenologically based concept of visco-elastic response.

The fact that the approach by Rajagopal and Srinivasa [34] closely follows the phenomenological concept of visco-elastic material is best seen in the interpretation of the quantity ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ that appears in the formula for the Cauchy stress tensor. This quantity can be interpreted as the left Cauchy–Green tensor associated with the elastic part of the fluid response. Using the approach by Rajagopal and Srinivasa [34], the derivation of Maxwell/Oldroyd-B type models was discussed by Málek et al. [35] (see also Hron et al. [36]). More complex viscoelastic rate-type models that document the applicability of the approach in more involved settings are discussed in Málek et al. [37], Málek et al. [38] and in Dostalík et al. [39]. Following Dostalík et al. [39], we know that the Giesekus fluid is a fluid with the specific Helmholtz free energy $\psi$ in the form

$$\psi {=}_{\mathrm{def}}-c_{\mathrm{V}}^{\mathrm{iNSE}}\theta [\ln \left(\frac{\theta }{{\theta }_{\mathrm{ref}}}\right)-1]+\frac{\mu }{2\rho }(\mathrm{Tr}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-3-\ln \det {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}),$$

(5)

where $\theta$ denotes the absolute temperature, ${\theta }_{\mathrm{ref}}$ denotes a constant reference temperature, $c_{\mathrm{V}}^{\mathrm{iNSE}}$ is a positive material parameter (specific heat capacity at constant volume) and $\mu$ is another positive material parameter. The specific Helmholtz free energy describes the energy storage ability of the fluid, and the chosen ansatz is the same as for the standard Maxwell/Oldroyd-B fluid. This implies that the Giesekus fluid and Maxwell/Oldroyd-B fluids differ, from the perspective of the current approach, only in their entropy production mechanisms (see below).

Specifying the Helmholtz free energy as a function of $\theta$ and ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$, one can use the standard thermodynamical identities $\eta =-\frac{\partial \psi }{\partial \theta }$ and $e=\psi +\theta \eta$, and obtain explicit formulae for the specific entropy $\eta$, and the specific internal energy e

$$\eta =c_{\mathrm{V}}^{\mathrm{iNSE}}\ln \left(\frac{\theta }{{\theta }_{\mathrm{ref}}}\right),e=c_{\mathrm{V}}^{\mathrm{iNSE}}\theta +\frac{\mu }{2\rho }(\mathrm{Tr}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-3-\ln \det {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}).$$

(6a)

Note that adding the kinetic energy to the mechanical part of the internal energy e, that is to the term $\frac{\mu }{2\rho }(\mathrm{Tr}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-3-\ln \det {\mathbb{B}}_{{\kappa }_{p\left(t\right)}})$, we can define the specific mechanical energy via

$$e_{\mathrm{mech}}{=}_{\mathrm{def}}\frac{1}{2}{\left|\mathit{v}\right|}^2+\frac{\mu }{2\rho }(\mathrm{Tr}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-3-\ln \det {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}).$$

(6b)

Once the Helmholtz free energy, and consequently also the internal energy, is specified, one can derive the evolution equation for the entropy that has the structure $\rho \frac{\mathrm{d}\eta }{\mathrm{d}t}+\mathrm{div}{\mathit{j}}_{\eta }=\xi ,$ where ${\mathit{j}}_{\eta }$ denotes the entropy flux and $\xi$ stands for the entropy production. In the case of Giesekus fluid the entropy production is given by the formula $\xi =\frac{\zeta }{\theta }$, where

$$\zeta {=}_{\mathrm{def}}2\nu \mathbb{D}:\mathbb{D}+\frac{{\mu }^2}{2{\nu }_1}\mathrm{Tr}[\alpha {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}^2+(1-3\alpha ){\mathbb{B}}_{{\kappa }_{p\left(t\right)}}+(1-\alpha ){\mathbb{B}}_{{\kappa }_{p\left(t\right)}}^{-1}+(3\alpha -2)\mathbb{I}]+\kappa \frac{{|\nabla \theta |}^2}{\theta }.$$

(7)

(We use the notation $\mathbb{A}:\mathbb{A}{=}_{\mathrm{def}}\mathrm{Tr}\left(\mathbb{A}{\mathbb{A}}^{\top }\right)$ for the scalar product on the space of matrices, and $\left|\mathbb{A}\right|$ for the corresponding Frobenius norm.) Since ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is a symmetric positive definite matrix, it is easy to check that the entropy production in Equation (7) is a nonegative quantity, hence the second law of thermodynamics is satisfied.

The fact that ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is a symmetric positive definite matrix follows directly from the governing equations in Equation (3) via an argument similar to that of Boyaval et al. [40]. It is also a consequence of the fact that ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is in the approach by Rajagopal and Srinivasa [34] constructed as the left Cauchy–Green tensor, which means that ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ can be decomposed as ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}={\mathbb{F}}_{{\kappa }_{p\left(t\right)}}{\mathbb{F}}_{{\kappa }_{p\left(t\right)}}^{\top }$. We exploit the positivity of ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ quite frequently in our analysis.

Finally, we introduce three more important quantities that play crucial role in the construction of Lyapunov type functional via the method proposed by Bulíček et al. [23]. Namely, we introduce the net total energy $E_{\mathrm{tot}}$, the net mechanical energy $E_{\mathrm{mech}}$ and the net entropy S of the fluid occupying the domain $\Omega$,

$$E_{\mathrm{tot}}{=}_{\mathrm{def}}{\int }_{\Omega }\rho [\frac{1}{2}{\left|\mathit{v}\right|}^2+e]\mathrm{d}\mathrm{v},E_{\mathrm{mech}}{=}_{\mathrm{def}}{\int }_{\Omega }\rho e_{\mathrm{mech}}\mathrm{d}\mathrm{v},S{=}_{\mathrm{def}}{\int }_{\Omega }\rho \eta \mathrm{d}\mathrm{v}.$$

(8)

2.3. Scaling

The equations in Equation (3) governing the evolution of mechanical variables can be transformed to a dimensionless form by introducing the characteristic length $x_{\mathrm{char}}$, characteristic time $t_{\mathrm{char}}$. (Note that the tensor field ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ already is a dimensionless quantity.) Using the following relations between the original quantities and their dimensionless versions denoted by stars $x=x_{\mathrm{char}}{x}^{\star }$, $t=t_{\mathrm{char}}{t}^{\star }$, $\mathit{v}=\frac{x_{\mathrm{char}}}{t_{\mathrm{char}}}{\mathit{v}}^{\star }$, $m=\frac{\nu }{t_{\mathrm{char}}}{m}^{\star }$, we obtain

$$\begin{array}{cc}\hfill {\mathrm{div}}^{\star }{\mathit{v}}^{\star }& =0,\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}^{\star }}{\mathrm{d}{t}^{\star }}& ={\mathrm{div}}^{\star }{\mathbb{T}}^{\star },\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill {\stackrel{▿}{\overline{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}}^{\star }& =-\frac{1}{\mathrm{Wi}}[\alpha {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}^2+(1-2\alpha ){\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-(1-\alpha )\mathbb{I}],\hfill \end{array}$$

(9c)

where the dimensionless Cauchy stress tensor ${\mathbb{T}}^{\star }$ is given by

$${\mathbb{T}}^{\star }=\frac{1}{\mathrm{Re}}{m}^{\star }\mathbb{I}+{\left({\mathbb{T}}^{\star }\right)}_{\delta },{\left({\mathbb{T}}^{\star }\right)}_{\delta }=\frac{2}{\mathrm{Re}}{\left({\mathbb{D}}^{\star }\right)}_{\delta }+\Xi {\left({\mathbb{B}}_{{\kappa }_{p\left(t\right)}}\right)}_{\delta }.$$

(9d)

In Equation (9), we introduce three dimensionless numbers—the Reynolds number $\mathrm{Re}$, the Weissenberg number $\mathrm{Wi}$ and dimensionless shear modulus $\Xi$—via the formulae $\mathrm{Re}{=}_{\mathrm{def}}\frac{\rho x_{\mathrm{char}}^2}{\nu t_{\mathrm{char}}}$, $\mathrm{Wi}{=}_{\mathrm{def}}\frac{{\nu }_1}{\mu t_{\mathrm{char}}}$, and $\Xi {=}_{\mathrm{def}}\frac{\mu t_{\mathrm{char}}^2}{\rho x_{\mathrm{char}}^2}$. It remains to introduce a scaling factor for the net mechanical energy $E_{\mathrm{mech}}$, which is used for the construction of the Lyapunov type functional in Section 4. Using the scaling $E_{\mathrm{mech}}=\frac{\rho x_{\mathrm{char}}^5}{t_{\mathrm{char}}^2}{E}_{\mathrm{mech}}^{\star }$, we obtain

$$E_{\mathrm{mech}}^{\star }\left({\mathit{W}}^{\star }\right)={\int }_{{\Omega }^{\star }}[\frac{1}{2}{\left|{\mathit{v}}^{\star }\right|}^2+\frac{\Xi }{2}(\mathrm{Tr}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-3-\ln \det {\mathbb{B}}_{{\kappa }_{p\left(t\right)}})]\mathrm{d}{\mathrm{v}}^{\star }.$$

(10)

Hereafter, we omit the star denoting dimensionless quantities unless otherwise specified.

The scaling is chosen in such a way that, if $\mathrm{Wi}\to 0+$, then ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ approaches the identity tensor. Indeed, if $\mathrm{Wi}\to 0+$ then, Equation (9c) implies that

$$\mathbb{O}=\alpha {\mathbb{B}}_{{\kappa }_{p\left(t\right)}}^2+(1-2\alpha ){\mathbb{B}}_{{\kappa }_{p\left(t\right)}}-(1-\alpha )\mathbb{I},$$

(11)

and the solution of Equation (11) is ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}=\mathbb{I}$. Moreover, if ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}=\mathbb{I}$, then the second term in Equation (10), that is the elastic contribution to the mechanical energy, vanishes, and the mechanical energy of the fluid reduces to the kinetic energy only. Finally, if ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}=\mathbb{I}$, then the additional term in the Cauchy stress tensor in Equation (9d) vanishes. This means that for $\mathrm{Wi}\to 0+$ the governing equations in Equation (9) reduce to the standard incompressible Navier–Stokes equations.

2.4. Boundary Conditions

The governing equations in Equation (9) must be supplemented with boundary conditions for the velocity $\mathit{v}$. We are interested in internal flow problems, where one prescribes Dirichlet boundary conditions on a part of the flow domain $\Omega \subset {\mathbb{R}}^3$, and periodic boundary conditions on another part of the domain. Such a domain is usually called the periodic cell. (For example, in the case of flow in between two infinite concentric rotating cylinders, the Dirichlet boundary condition is prescribed on the surfaces of the cylinders, while the periodic boundary condition is prescribed in the direction of the axis of the cylinders.) On the parts of the boundary corresponding to the periodicity directions, say ${\Gamma }_1$, we therefore prescribe periodic boundary condition for $\mathit{v}$, while, on the remaining part of the boundary, say ${\Gamma }_2$, we prescribe the no-penetration and the no-slip boundary condition,

$$\begin{array}{cc}\hfill {\mathit{v}•\mathit{n}|}_{{\Gamma }_2}& =0,\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill {(\mathbb{I}-\mathit{n}\otimes \mathit{n})\mathit{v}|}_{{\Gamma }_2}& =\mathit{V},\hfill \end{array}$$

(12b)

where $\mathit{n}$ is the unit outward normal to the boundary of $\Omega$ and $\mathit{V}$ is a given velocity in the tangential direction to the boundary. This means that the fluid adheres to the boundary, and, moreover, if $\mathit{V}\ne \mathbf{0}$, then, in general, the energy is exchanged between the fluid and its surroundings. Indeed, the balance of the net total energy reads $\frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}={\int }_{\partial \Omega }\left(\mathbb{T}\mathit{v}\right)•\mathit{n}\mathrm{d}\mathrm{s}-{\int }_{\partial \Omega }{\mathit{j}}_q•\mathit{n}\mathrm{d}\mathrm{s},$ where $\partial \Omega$ denotes the boundary of the domain $\Omega$ and ${\mathit{j}}_q$ denotes the heat flux. Consequently, if $\mathit{v}\ne \mathbf{0}$ on the boundary, then the term ${\int }_{\partial \Omega }\left(\mathbb{T}\mathit{v}\right)•\mathit{n}\mathrm{d}\mathrm{s}$ does not, in general, vanish or is compensated by the second term on the right-hand side, and the net total energy might even change in time.

Concerning the boundary conditions for the perturbation $\tilde{\mathit{v}}$ with respect to the reference state $\widehat{\mathit{v}}$ (see below), we see that, if $\widehat{\mathit{v}}$ satisfies Equation (12), then the perturbed state $\mathit{v}=\widehat{\mathit{v}}+\tilde{\mathit{v}}$ also satisfies Equation (12) provided that

$${\tilde{\mathit{v}}|}_{{\Gamma }_2}=\mathbf{0}.$$

(13)

The periodic boundary condition on ${\Gamma }_1$ is preserved for the perturbation $\tilde{\mathit{v}}$.

In the following, we frequently use the identity

$${\int }_{\partial \Omega }\mathit{f}•\mathit{n}\mathrm{d}\mathrm{s}=0,$$

(14)

where $\mathit{f}:\partial \Omega \to {\mathbb{R}}^3$ is a smooth function such that $\mathit{f}$ fulfills the periodic boundary condition on ${\Gamma }_1$ and $\mathit{f}=\mathbf{0}$ on ${\Gamma }_2$. Note that the identity holds even if one part of the boundary, whether ${\Gamma }_1$ or ${\Gamma }_2$, is not present.

3. Base Flow—Non-Equilibrium Steady State

3.1. Notation for the Stability Analysis

We are interested in the evolution of the triplet $\mathit{W}{=}_{\mathrm{def}}[\mathit{v},m,{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}],$ which solves the evolution equations in Equation (3). We further use the notation $\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$ for the triplet corresponding to a non-equilibrium steady state solution, and $\tilde{\mathit{W}}=[\tilde{\mathit{v}},\tilde{m},\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$ for the perturbation with respect to the non-equilibrium steady state. This means that the triplet describing the complete perturbed state is given as a sum of the reference state $\widehat{\mathit{W}}$ and the perturbation $\tilde{\mathit{W}}$ with respect to the reference state

$$\begin{array}{cc}\hfill \mathit{W}& =\widehat{\mathit{W}}+\tilde{\mathit{W}},\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill [\mathit{v},m,{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}]& =[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]+[\tilde{\mathit{v}},\tilde{m},\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}].\hfill \end{array}$$

(15b)

(Note that sometimes we work only with the pair $\mathit{W}{=}_{\mathrm{def}}[\mathit{v},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}]$, since the pressure is insubstantial in our analysis.) The term non-equilibrium steady state is chosen in accordance with the practice in thermodynamics, and it means that the entropy is produced at the steady state $\widehat{\mathit{W}}$. In particular, the adjective non-equilibrium does not refer to the stability of the steady state.

3.2. Governing Equations in a Steady State

The steady state$\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$whose stability we want to investigate is a solution to the equations in Equation (9) where the partial time derivatives are identically equal to zero. In particular, we assume that the state described by the triplet $[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$ solves the system

$$\begin{array}{cc}\hfill \mathrm{div}\widehat{\mathit{v}}& =0,\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill (\widehat{\mathit{v}}•\nabla )\widehat{\mathit{v}}& =\mathrm{div}\mathbb{T}(\widehat{\mathit{W}}),\hfill \end{array}$$

(16b)

$$\begin{array}{cc}\hfill (\widehat{\mathit{v}}•\nabla )\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-\widehat{\mathbb{L}}\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{\mathbb{L}}}^{\top }& =-\frac{1}{\mathrm{Wi}}[\alpha {\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2+(1-2\alpha )\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-(1-\alpha )\mathbb{I}].\hfill \end{array}$$

(16c)

subject to the boundary conditions in Equation (12) on ${\Gamma }_2$, that is

$${\widehat{\mathit{v}}•\mathit{n}|}_{{\Gamma }_2}=0,{(\mathbb{I}-\mathit{n}\otimes \mathit{n})\widehat{\mathit{v}}|}_{{\Gamma }_2}=\mathit{V},$$

(17)

and the periodic boundary conditions on ${\Gamma }_1$. Here, the symbol $\mathbb{T}(\widehat{\mathit{W}})$ denotes the Cauchy stress tensor induced by the triplet $[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$, that is

$$\mathbb{T}(\widehat{\mathit{W}})=\frac{1}{\mathrm{Re}}\widehat{m}\mathbb{I}+\frac{2}{\mathrm{Re}}\widehat{\mathbb{D}}+\Xi {\left(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)}_{\delta },$$

(18)

where $\widehat{\mathbb{D}}=\frac{1}{2}(\widehat{\mathbb{L}}+{\widehat{\mathbb{L}}}^{\top })$, and $\widehat{\mathbb{L}}=\nabla \widehat{\mathit{v}}$.

Note that, if $\mathit{V}=\mathbf{0}$, that is if no mechanical energy is supplied to the fluid, then the system would admit an equilibrium solution

$$[\widehat{\mathit{v}},\widehat{m},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]=[\mathbf{0},c,\mathbb{I}],$$

(19)

where c is an arbitrary number. (This is the standard ambiguity in the identification of the pressure well known from the case of Navier–Stokes fluid.) Here, we use the adjective equilibrium to emphasise that such a steady state would lead to zero entropy production. Indeed, if ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}=\mathbb{I}$ and $\mathit{v}=\mathbf{0}$, then the (mechanical part) of the entropy production in Equation (7) vanishes.

On the other hand, if $\mathit{V}\ne \mathbf{0}$, then one must in general expect that the steady fields $\widehat{\mathit{v}}$ and $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ are spatially inhomogeneous, and consequently the entropy production in Equation (7) is positive. This means that the system produces the entropy; hence, it is, from the thermodynamical point of view, out of equilibrium. Consequently, as discussed above, we use the adjective non-equilibrium and we refer to the base flow as of non-equilibrium steady state.

3.3. Concept of Stability

Concerning the stability of the non-equilibrium steady state, we are interested in its asymptotic stability. If we have a non-equilibrium steady state $\widehat{\mathit{W}}$ that solves Equation (16), then we want to know whether the perturbation $\mathit{W}=\widehat{\mathit{W}}+\tilde{\mathit{W}}$ of the non-equilibrium steady state $\widehat{\mathit{W}}$ tends back to the non-equilibrium steady state $\widehat{\mathit{W}}$ as time goes to infinity. In our case, the evolution of the perturbed state $\mathit{W}$ is governed by the equations in Equation (9) that must be solved subject to the given boundary conditions in Equation (12) and subject to initial conditions

$${\mathit{v}|}_{t=0}=\widehat{\mathit{v}}+{\tilde{\mathit{v}}}_0,{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}|}_{t=0}=\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+{\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)}_0.$$

(20)

The non-equilibrium steady state $\widehat{\mathit{W}}$ is said to be asymptotically stable provided that the solutions that start close enough to the steady state not only remain close enough to the steady state but also eventually converge to the steady state, that is the triplet $\mathit{W}$ converges to $\widehat{\mathit{W}}$ as time goes to infinity,

$$\mathit{W}\stackrel{t\to +\infty }{\to }\widehat{\mathit{W}},$$

(21)

for all sufficiently small initial data ${\tilde{\mathit{v}}}_0$ and ${\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)}_0$ (see (Henry [41], Definition 4.1.2) or Yoshizawa [42] for the formal definition and further discussion).

Ideally, one would like to obtain stronger results. Namely, one would like to have an unconditional result that states that the non-equilibrium steady state is recovered as time goes to infinity regardless of the choice of initial perturbation. This behaviour is expected if one deals with non-equilibrium steady states that are driven by a small energy inflow that is by a small boundary velocity $\mathit{V}$, or, in other words, if one deals with non-equilibrium steady states that are not far away from the equilibrium steady state.

The key task in the stability analysis is the choice of a metric/norm on the state space to give a meaning to the statement in Equation (21). Namely, we need to answer the question as how to characterise the distance between $\widehat{\mathit{W}}$ and $\mathit{W}$, since Equation (21) means

$$\mathrm{dist}(\widehat{\mathit{W}},\mathit{W})\stackrel{t\to +\infty }{\to }0,$$

(22)

where $\mathrm{dist}(·,·)$ is a given metric that is not necessarily induced by a norm. Since ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is at a given spatial point $\mathit{x}\in \Omega$ a positive definite matrix, it seems reasonable to design the metric in such a way that it reflects this fact. This means that we have to rely on a metric on the set of positive definite matrices. There are several possible definitions of the metric on these sets (see Appendix A). If we use the Bures–Wasserstein distance,

$${\mathrm{dist}}_{\mathbb{P}\left(d\right),\mathrm{BW}}(\mathbb{A},\mathbb{B}){=}_{\mathrm{def}}{\{\mathrm{Tr}\mathbb{A}+\mathrm{Tr}\mathbb{B}-2\mathrm{Tr}\left[{\left({\mathbb{A}}^{\frac{1}{2}}\mathbb{B}{\mathbb{A}}^{\frac{1}{2}}\right)}^{\frac{1}{2}}\right]\}}^{\frac{1}{2}},$$

(23)

(see Equation (A1) in Appendix A), and if we generalise this concept to the spatially distributed tensor fields, then we can define the distance between $\widehat{\mathit{W}}$ and $\mathit{W}$ as

$$\mathrm{dist}(\widehat{\mathit{W}},\mathit{W}){=}_{\mathrm{def}}{({\parallel \widehat{\mathit{v}}-\mathit{v}\parallel }_{L^2(\Omega )}^2+{\left[{\mathrm{dist}}_{{\mathbb{P}}_{\Omega }\left(d\right),\mathrm{BW}}(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}})\right]}^2)}^{\frac{1}{2}},$$

(24)

(see Equation (A7) and Appendix A for a discussion of the notation and correctness of this definition). It turns out that this concept of distance nicely fits to the dynamical system we are interested in.

The term “stability” is used in many other contexts; hence, we briefly comment on these other notions of stability. In particular, we emphasise what is in the present work not meant by the stability. First, we are not interested in the stability in the sense of continuous dependence on initial data, which is the concept of stability investigated in Dafermos [43] and various subsequent works especially in the theory of hyperbolic systems (see Dafermos [44]). The stability in the sense of continuous dependence on initial data means (see, for example, Schaeffer and Cain [45]) that

[…] if the initial data for an initial value problem are altered slightly, then the perturbed solution diverges from the original solution no faster than at a controlled exponential rate.

Apparently, the asymptotic stability we are interested in is a more ambitious concept, since we want the perturbed solution to converge back to the original solution (non-equilibrium steady state). Second, we are not interested in the stability of the steady state subject to infinitesimal perturbations, that is in the linearised stability. We are interested in the evolution of finite amplitude perturbations.

Finally, we emphasise that in our analysis we work with perturbations that are solution to the governing equations in the classical sense. (All derivatives are understood as the classical derivatives, not as generalised derivatives such as distributional derivatives and so forth.) In particular, we do not consider the perturbations that solve the governing equations only in a weak sense, although it is an important issue worth further investigation. The reader interested in the discussion of the state-of-the-art rigorous mathematical theory of equations governing the motion of viscoelastic fluids is kindly referred to the work of Masmoudi [46] or Barrett and Süli [47].

4. Lyapunov Functional

4.1. Concept of Lyapunov Functional

Let us briefly recall the concept of Lyapunov functional (see Henry [41]). We consider a system of governing equations in the form

$$\frac{\mathrm{d}\mathit{X}}{\mathrm{d}t}=\mathit{F}\left(\mathit{X}\right),$$

(25)

where $\widehat{\mathit{X}}$ is a steady state, that is $\mathit{F}(\widehat{\mathit{X}})=\mathbf{0}$, and where ${\parallel ·\parallel }_{\mathrm{st}}$ denotes a norm on the underlying state space. We say that the functional $\mathcal{V}(\tilde{\mathit{X}}\parallel \widehat{\mathit{X}})$ is a strict Lyapunov functional of the steady state $\widehat{\mathit{X}}$ provided that:

• There exists a neighbourhood of $\widehat{\mathit{X}}$ such that the functional is bounded from below by a function f of the distance between the steady state $\widehat{\mathit{X}}$ and the perturbation $\mathit{X}$, that is

  $$\mathcal{V}(\tilde{\mathit{X}}\parallel \widehat{\mathit{X}})\ge f\left({\parallel \widehat{\mathit{X}}-\mathit{X}\parallel }_{\mathrm{st}}\right),$$

  (26a)

  where f is a continuous strictly increasing function such that $f\left(0\right)=0$ and $f\left(r\right)>0$ whenever $r>0$.
• The time derivative of $\mathcal{V}(\tilde{\mathit{X}}\parallel \widehat{\mathit{X}})$ is negative and bounded from above by a function g of the distance between the steady state $\widehat{\mathit{X}}$ and the perturbation $\mathit{X}$, that is

  $$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{V}(\tilde{\mathit{X}}\parallel \widehat{\mathit{X}})\le -g\left({\parallel \widehat{\mathit{X}}-\mathit{X}\parallel }_{\mathrm{st}}\right),$$

  (26b)

  where g is a continuous strictly increasing function such that $g\left(0\right)=0$ and $g\left(r\right)>0$ whenever $r>0$.

If the given system of governing equations admits a strict Lyapunov functional near the state $\widehat{\mathit{X}}$, then we know that the steady state $\widehat{\mathit{X}}$ is asymptotically stable (see (Henry [41], Theorem 4.1.4)). This means that the solution $\mathit{X}=\tilde{\mathit{X}}+\widehat{\mathit{X}}$ that starts in the neighbourhood of $\widehat{\mathit{X}}$ satisfies

$${\parallel \widehat{\mathit{X}}-\mathit{X}\parallel }_{\mathrm{st}}\stackrel{t\to +\infty }{\to }0.$$

(27)

While the concept of Lyapunov type functional is very simple, it is difficult to apply in a particular setting. The main difficulty is to find a Lyapunov type functional. (Note that in the infinite dimensional setting one can not easily exploit LaSalle’s invariance principle since it requires precompactness of the trajectories, which is a qualitative property that goes beyond our assumption regarding the existence of the classical solution. Consequently, having a Lyapunov type functional with Equation (26b) replaced by the mere negativity everywhere except at the equilibrium, $\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{V}(\tilde{\mathit{X}}\parallel \widehat{\mathit{X}})<0$ for $\tilde{\mathit{X}}\ne \mathbf{0}$, is not a viable option.) Fortunately, since we are interested in equations describing a physical system, we can try to search for the functional using physical concepts.

If we were dealing with the stability of a homogeneous equilibrium steady state in a thermodynamically isolated system, then a Lyapunov type functional could be constructed using the net entropy S and the net total energy $E_{\mathrm{tot}}$ functional. It can be shown that the appropriate Lyapunov type functional is in this setting reads

$${\mathcal{V}}_{\mathrm{eq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})=-[S(\widehat{\mathit{W}}+\tilde{\mathit{W}})-\frac{1}{\widehat{\theta }}\{E_{\mathrm{tot}}(\widehat{\mathit{W}}+\tilde{\mathit{W}})-E_{\mathrm{tot}}\left(\widehat{\mathit{W}}\right)\}],$$

(28)

where $\widehat{\theta }$ is the spatially homogeneous temperature in the equilibrium steady state (see Bulíček et al. [23] for details). The observation that Equation (28) can be used as a suitable Lyapunov type functional for stability analysis of homogeneous equilibrium steady states in a thermodynamically isolated systems or mechanically isolated systems immersed in a thermal bath is well known (see, for example, Šilhavý [48] or Grmela and Öttinger [31], Öttinger and Grmela [32]), and in the continuum thermodynamics setting it dates back to the works of Coleman [49], Gurtin [50], and Gurtin [51]. (The core idea can be found in earlier works, see especially Duhem [52].) Unfortunately, the same functional cannot be used in stability analysis of non-equilibrium spatially inhomogeneous steady states in thermodynamically open systems. This fact is clear from Equation (28) itself. If one works with a spatially inhomogeneous steady states, then $\widehat{\theta }$ in Equation (28) is a function, and Equation (28) does not define a functional at all.

4.2. Construction of Lyapunov Type Functional for Stability Analysis of a Spatially Inhomogeneous Steady State

Recently, Bulíček et al. [23] proposed a method for construction of Lyapunov type functionals for stability analysis of non-equilibrium spatially inhomogeneous steady states in thermodynamically open systems. In the ongoing analysis, we use the same ideas as in Bulíček et al. [23], but we restrict ourselves to the mechanical quantities only. This is a matter of convenience, since we are interested in mechanical quantities only, and the temperature evolution has no feedback on the mechanical part of the system of governing equations. Consequently, we do not need to work with the Lyapunov type functional for the full thermomechanical problem, and we can construct a simpler Lyapunov type functional solely for the mechanical quantities.

Using the net mechanical energy functional $E_{\mathrm{mech}}$ introduced in Equation (10), one can see that the net mechanical energy in a thermodynamically closed system must decay in time. Consequently, the functional

$${\mathcal{V}}_{\mathrm{eq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}}){=}_{\mathrm{def}}{E}_{\mathrm{mech}}(\widehat{\mathit{W}}+\tilde{\mathit{W}})-E_{\mathrm{mech}}\left(\widehat{\mathit{W}}\right)$$

(29)

can serve as a Lyapunov type functional for stability analysis of equilibrium spatially homogeneous state in Equation (19).

Following the methodology outlined in Bulíček et al. [23], we use the Lyapunov type functional for the equilibrium steady state in Equation (29), and we define the candidate for the Lyapunov type functional for the non-equilibrium steady state as

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}}){=}_{\mathrm{def}}{E}_{\mathrm{mech}}(\widehat{\mathit{W}}+\tilde{\mathit{W}})-E_{\mathrm{mech}}\left(\widehat{\mathit{W}}\right)-{{\mathrm{D}}_{\mathit{W}}{E}_{\mathrm{mech}}\left(\mathit{W}\right)|}_{\mathit{W}=\widehat{\mathit{W}}}\left[\tilde{\mathit{W}}\right],$$

(30)

where ${{\mathrm{D}}_{\mathit{W}}{E}_{\mathrm{mech}}\left(\mathit{W}\right)|}_{\mathit{W}=\widehat{\mathit{W}}}\left[\tilde{\mathit{W}}\right]$ denotes the Gâteaux derivative at point $\widehat{\mathit{W}}$ in the direction $\tilde{\mathit{W}}$. (This is essentially the affine correction trick introduced in a different context by Ericksen [53] (see also Dafermos [43]).) The (dimensionless) explicit formulae for the individual terms in Equation (30) read

$$\begin{array}{cc}\hfill E_{\mathrm{mech}}(\widehat{\mathit{W}}+\tilde{\mathit{W}})& ={\int }_{\Omega }[\frac{1}{2}{|\widehat{\mathit{v}}+\tilde{\mathit{v}}|}^2+\frac{\Xi }{2}(\mathrm{Tr}(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})-3-\ln \det (\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}))]\mathrm{d}\mathrm{v},\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill E_{\mathrm{mech}}\left(\widehat{\mathit{W}}\right)& ={\int }_{\Omega }[\frac{1}{2}{\left|\widehat{\mathit{v}}\right|}^2+\frac{\Xi }{2}(\mathrm{Tr}\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-3-\ln \det \widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})]\mathrm{d}\mathrm{v},\hfill \end{array}$$

(31b)

$$\begin{array}{cc}\hfill {{\mathrm{D}}_{\mathit{W}}{E}_{\mathrm{mech}}\left(\mathit{W}\right)|}_{\mathit{W}=\widehat{\mathit{W}}}\left[\tilde{\mathit{W}}\right]& ={\int }_{\Omega }\{\widehat{\mathit{v}}•\tilde{\mathit{v}}+\frac{\Xi }{2}[\mathrm{Tr}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-\mathrm{Tr}\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)]\}\mathrm{d}\mathrm{v}.\hfill \end{array}$$

(31c)

Using Equation (31) in Equation (30), we get, after some algebraic manipulations, the explicit formula for the proposed Lyapunov type functional

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})={\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}[-\ln \det (\mathbb{I}+{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})+\mathrm{Tr}\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)]\mathrm{d}\mathrm{v}.$$

(32)

It remains to show that the functional in Equation (32) has the properties in Equation (26) introduced in Section 4.1. First, we show that the condition in Equation (26a) holds for a neighbourhood of $\widehat{\mathit{W}}$, which means that we have to specify a norm on the state space.

The suitable norm is the norm introduced in Appendix A in Definition A3. This norm is a “standard” Lebesgue type norm. However, it turns out that it is convenient to use this norm for a characterisation of the evolution of a “shifted” state. The idea is the following. If we are given a constant-in-time tensor field $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$, which corresponds to the steady solution of Equation (9), and a state $\mathit{W}=[\mathit{v},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}]$, then we can introduce the shifted state $\mathit{Z}=[\mathit{v},{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}]$ that is defined as

$$\mathit{W}=[\mathit{v},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}]\mapsto \mathit{Z}=[\mathit{v},{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}],{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}{=}_{\mathrm{def}}{\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right)}^{\frac{1}{2}}.$$

(33)

This shifted state seems to be ideal for the investigation of the perturbations to the steady state $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$. Indeed, the steady state for Equation (9) is $\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$, which is in virtue of Equation (33) shifted to

$$\widehat{\mathit{Z}}=[\widehat{\mathit{v}},\mathbb{I}].$$

(34)

Now, instead of investigating the behaviour of the perturbation ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ with respect to $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$, we can investigate the behaviour of the shifted perturbation ${\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}$ with respect to the identity $\mathbb{I}$.

Lemma 1

(Relation between the proposed Lyapunov functional and a norm). Let $\widehat{\mathit{W}}$ and $\mathit{W}=\widehat{\mathit{W}}+\tilde{\mathit{W}}$ denote two states of the system governed by equations in Equation (9), and let $\widehat{\mathit{Z}}$ and $\mathit{Z}$ denote the corresponding shifted states. Furthermore, let ${\parallel ·\parallel }_{\mathrm{st}}$ denote the norm introduced in Definition A3. Then, there exists a positive constant $D(\Xi )$ such that

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\ge D{\parallel \widehat{\mathit{Z}}-\mathit{Z}\parallel }_{\mathrm{st}}^2.$$

(35)

Proof. 

We note that Equation (32) for the Lyapunov type functional can be rewritten as

$$\begin{array}{}& {\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}}=\frac{1}{2}{\parallel \widehat{\mathit{v}}-(\widehat{\mathit{v}}+\tilde{\mathit{v}})\parallel }_{L^2(\Omega )}^2\hfill & +{\int }_{\Omega }\frac{\Xi }{2}[\mathrm{Tr}\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}){\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right)-3-\ln \det \left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}){\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right)]\mathrm{d}\mathrm{v}\left(36\right)& \hfill =\frac{1}{2}{\parallel \widehat{\mathit{v}}-(\widehat{\mathit{v}}+\tilde{\mathit{v}})\parallel }_{L^2(\Omega )}^2+{\int }_{\Omega }\frac{\Xi }{2}[\mathrm{Tr}{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}^2-3-\ln \det {\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}^2]\mathrm{d}\mathrm{v},\end{array}$$

where we use the cyclic property of the trace and the properties of the determinant. Making use of Lemma A3 and the inequality in the integrand of the last term in Equation (36), we see that

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\ge \frac{1}{2}{\parallel \widehat{\mathit{v}}-(\widehat{\mathit{v}}+\tilde{\mathit{v}})\parallel }_{L^2(\Omega )}^2+\frac{\Xi }{2}{\int }_{\Omega }{|\mathbb{I}-{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}|}^2\mathrm{d}\mathrm{v}\ge \min \{\frac{1}{2},\frac{\Xi }{2}\}{\parallel \widehat{\mathit{Z}}-\mathit{Z}\parallel }_{\mathrm{st}}^2.$$

(37)

 □

The less straightforward part of the analysis of properties of proposed Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ is the evaluation of its time derivative $\frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}$. The formula for the time derivative is derived via a lengthy algebraic manipulation described in Appendix B, and it is given below in Lemma 2. Note that, although we are working with a thermodynamically open system, the formula for the time derivative does not contain boundary terms.

Lemma 2

(Explicit formula for the time derivative of the Lyapunov type functional). Let $\widehat{\mathit{W}}$ and $\mathit{W}=\widehat{\mathit{W}}+\tilde{\mathit{W}}$ denote two states of the system governed by equations in Equation (9), where the state $\widehat{\mathit{W}}$ is the steady state, that is it solves Equation (16). The formula for the time derivative of the functional ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ introduced in Equation (32) reads

$$\begin{array}{}& \frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})=-{\int }_{\Omega }\frac{2}{\mathrm{Re}}\tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}-{\int }_{\Omega }\Xi \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}-{\int }_{\Omega }\widehat{\mathbb{D}}\tilde{\mathit{v}}•\tilde{\mathit{v}}\mathrm{d}\mathrm{v}\hfill & -{\int }_{\Omega }\frac{\Xi }{2}\mathrm{Tr}\left[{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}(\tilde{\mathit{v}}•\nabla )\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right]\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}:(\tilde{\mathbb{L}}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\tilde{\mathbb{L}}}^{\top })\mathrm{d}\mathrm{v}& -{\int }_{\Omega }(1-\alpha )\frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}\left[{(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})}^{-1}\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right){\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right)}^{\top }\right]\mathrm{d}\mathrm{v}\left(38\right)& \hfill -{\int }_{\Omega }\alpha \frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}\left[{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}{\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2\right]\mathrm{d}\mathrm{v}.\end{array}$$

Proof. 

See Appendix B.  □

We note that the terms

$$\begin{array}{c}-{\int }_{\Omega }\frac{2}{\mathrm{Re}}\tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v},-{\int }_{\Omega }\alpha \frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}\left[{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}{\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2\right]\mathrm{d}\mathrm{v},\end{array}$$

(39a)

$$\begin{array}{c}-{\int }_{\Omega }(1-\alpha )\frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}\left[{(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}+\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})}^{-1}\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right){\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right)}^{\top }\right]\mathrm{d}\mathrm{v},\end{array}$$

(39b)

are negative provided that $\tilde{\mathit{v}}\ne \mathbf{0}$ and $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\ne \mathbb{O}$. If we were able to show that these damping terms are strong enough to balance all the remaining terms on the right-hand side of Equation (38), we would get the desired result (Equation (26b)) concerning the negativity of the time derivative. This should be possible at least for sufficiently small Reynolds number $\mathrm{Re}$ and Weissenberg number $\mathrm{Wi}$. The hypothesis follows from the observation that as $\mathrm{Re}$ and $\mathrm{Wi}$ tend to zero, then the magnitude of the damping terms increases, and it should outgrow the other terms in Equation (38) that do not depend on $\mathrm{Re}$ and $\mathrm{Wi}$. This observation is consistent with the expectation that low Reynolds number and low Weissenberg number flows are stable.

Now, the objective is to show that the hypothesis is true, and that the proposed functional indeed satisfies the condition in Equation (26b). In the quantification of the “sufficient smallness” of the Weissenberg number $\mathrm{Wi}$ and the Reynolds number $\mathrm{Re}$, we aim at a simple but very rough estimate based on the elementary use of Friedrichs–Poincaré, Cauchy–Schwarz, Young and Korn (in)equalities (see Nečas et al. [54] or Evans [55] or any other standard reference work on function spaces). A precise characterisation of the Reynolds number and the Weissenberg number that guarantee the negativity of the time derivative, and hence the stability, could be obtained by a variational technique known from the standard energy method (see Joseph [56] or Straughan [57]). This is however beyond the scope of the current contribution.

Lemma 3

(Estimate on the time derivative). Let $\widehat{\mathit{W}}$ and $\mathit{W}=\widehat{\mathit{W}}+\tilde{\mathit{W}}$ denote two states of the system governed by equations in Equation (9), where the state $\widehat{\mathit{W}}$ is the steady state, that is it solves Equation (16). Then, there exist constants $C_1(\widehat{\mathit{W}},\mathrm{Re},\mathrm{Wi},\Xi ,\Omega )$ and $C_2(\widehat{\mathit{W}},\mathrm{Re},\mathrm{Wi},\Xi )$ such that the time derivative of the functional ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ introduced in Equation (32) can be estimated from above as

$$\frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\le C_1{\parallel \nabla \tilde{\mathit{v}}\parallel }_{L^2(\Omega )}^2+C_2{\parallel \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\parallel }_{L^2(\Omega )}^2,$$

(40)

where we denote

$$C_1{=}_{\mathrm{def}}-\frac{1}{\mathrm{Re}}+C_P\underset{\mathit{x}\in \Omega }{\sup }\left|{\lambda }_{\min }(\widehat{\mathbb{D}})\right|+\frac{\Xi }{2}\underset{\mathit{x}\in \Omega }{\sup }|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}-\mathbb{I}|+C_P\frac{\Xi }{4}\underset{\mathit{x}\in \Omega }{\sup }{\left|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right|}^2\underset{\mathit{x}\in \Omega }{\sup }|\nabla \widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}|,$$

(41a)

$$C_2{=}_{\mathrm{def}}-\alpha \frac{\Xi }{2\mathrm{Wi}}\underset{\mathit{x}\in \Omega }{\inf }{\lambda }_{\min }({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1})+\frac{\Xi }{2}\underset{\mathit{x}\in \Omega }{\sup }|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}-\mathbb{I}|+\frac{\Xi }{4}\underset{\mathit{x}\in \Omega }{\sup }{\left|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\right|}^2\underset{\mathit{x}\in \Omega }{\sup }|\nabla \widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}|,$$

(41b)

and where ${\lambda }_{\min }(·)$ denotes the minimal eigenvalue of the corresponding matrix and $C_P$ denotes the domain dependent constant from Friedrichs–Poincaré inequality.

Proof. 

See Appendix C.  □

Lemma 4

(Estimate on the time derivative in terms of the norm on the shifted state space). Let $\widehat{\mathit{W}}$ and $\mathit{W}=\widehat{\mathit{W}}+\tilde{\mathit{W}}$ denote two states of the system governed by equations in Equation (9), where the state $\widehat{\mathit{W}}$ is the steady state, that is it solves Equation (16), and let $\widehat{\mathit{Z}}$ and $\mathit{Z}$ denote the corresponding shifted states (see Equation (33)). Let us further assume that the constants $C_1$ and $C_2$ in Lemma 3 are negative. Then, there exists a positive constant $C(C_1,C_2,\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})$ such that

$$\frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\le -C{\parallel \widehat{\mathit{Z}}-\mathit{Z}\parallel }_{\mathrm{st}}^2,$$

(42)

where ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ denotes the functional introduced in Equation (32) and ${\parallel ·\parallel }_{\mathrm{st}}$ is the norm introduced in Definition A3.

Proof. 

In virtue of Lemma 3, we already know the estimate in Equation (40). Making use of Friedrichs–Poincaré inequality ${\parallel \mathit{v}\parallel }_{L^2(\Omega )}^2\le C_P{\parallel \nabla \mathit{v}\parallel }_{L^2(\Omega )}^2,$ where $C_P$ is the domain dependent constant. We see that, if $C_1<0$ and $C_2<0$, then Equation (40) implies

$$\frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\le -\frac{\left|C_1\right|}{C_P}{\parallel \tilde{\mathit{v}}\parallel }_{L^2(\Omega )}^2-\left|C_2\right|{\parallel \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\parallel }_{L^2(\Omega )}^2.$$

(43)

Next, we use Lemma A4, which implies that

$$\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|=|\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}-{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}|\ge {\left|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right|}^{-2}|\mathbb{I}-{\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right)}^{\frac{1}{2}}|={\left|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right|}^{-2}|\mathbb{I}-{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}|,$$

(44)

where ${\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}$ denotes the shifted state introduced in Equation (33). Consequently, we see that

$$-|C_2|{\parallel \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\parallel }_{L^2(\Omega )}^2=-\left|C_2\right|{\int }_{\Omega }{\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2\mathrm{d}\mathrm{v}\le -\left|C_2\right|{\left(\underset{\mathit{x}\in \Omega }{\sup }\left|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\right|\right)}^{-4}{\int }_{\Omega }{|\mathbb{I}-{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}|}^2\mathrm{d}\mathrm{v},$$

(45)

which means that Equation (43) can be further manipulated to the form

$$\frac{\mathrm{d}{\mathcal{V}}_{\mathrm{neq}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})\le -\min \{\frac{\left|C_1\right|}{C_P},\left|C_2\right|{\left(\underset{\mathit{x}\in \Omega }{\sup }\right|{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\left|\right)}^{-4}\}({\parallel \tilde{\mathit{v}}\parallel }_{L^2(\Omega )}^2+{\parallel \mathbb{I}-{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}\parallel }_{L^2(\Omega )}^2).$$

(46)

The inequality in Equation (46) gives in virtue of the definition of the norm ${\parallel ·\parallel }_{\mathrm{st}}$ (see Definition A3 and the proposition in Equation (42)). (Recall that the transformation in Equation (33) implies that $\widehat{{\mathbb{Z}}_{{\kappa }_{p\left(t\right)}}}=\mathbb{I}$.)  □

5. Main Result

Using the estimate from Lemma 4 and the relation between the metric and the functional ${\mathcal{V}}_{\mathrm{neq}}$ (see Lemma 1), it is straightforward to prove the following theorem.

Theorem 1

(Sufficient conditions for unconditional asymptotic stability). Let the pair $\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$ solve the governing equations for the steady state in Equation (16) with the boundary conditions in Equation (17). If the Reynolds number $\mathrm{Re}$, the Weissenberg number $\mathrm{Wi}$ and the dimensionless shear modulus Ξ are such that the constants $C_1$ and $C_2$ introduced in Equation (41) satisfy

$$C_1<0,C_2<0,$$

(47)

then the spatially inhomogeneous non-equilibrium steady state $\widehat{\mathit{W}}$ is unconditionally asymptotically stable, that is

$$\mathrm{dist}(\widehat{\mathit{W}},\mathit{W})\stackrel{t\to +\infty }{\to }0,$$

(48)

holds for any initial perturbation $\mathit{W}$, where the metric $\mathrm{dist}(·,·)$ is the metric introduced in Equation (A7).

Proof. 

We first investigate the stability in the shifted state space (see Equation (33)), that is we investigate perturbation $\mathit{Z}$ with respect to $\widehat{\mathit{Z}}=[\widehat{\mathit{v}},\mathbb{I}]$. We introduce the functional ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ (see Equation (32) and the equivalent expression in Equation (36)). The functional satisfies the condition in Equation (26a) (see Lemma 1). Furthermore, if $C_1$ and $C_2$ are negative, then Lemma 4 implies that the functional ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ decreases along the trajectories in a desired manner, that is it satisfies Equation (26b). Consequently, the functional ${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})$ is a genuine Lyapunov type functional for any neighbourhood of the steady state $\widehat{\mathit{Z}}$, hence $\widehat{\mathit{Z}}$ is unconditionally asymptotically stable, ${\parallel \widehat{\mathit{Z}}-\mathit{Z}\parallel }_{\mathrm{st}}\stackrel{t\to +\infty }{\to }0$.

The convergence in the norm on the shifted space seems to be an obscure characterisation of the approach to the equilibrium. However, if we exploit the definition of the shifted state (see Equation (33), and the equality in Equation (A17) proved in Lemma A3), we see that

$${\parallel \widehat{\mathit{Z}}-\mathit{Z}\parallel }_{\mathrm{st}}={({\parallel \widehat{\mathit{v}}-\mathit{v}\parallel }_{L^2(\Omega )}^2+{\int }_{\Omega }{\left({\mathrm{dist}}_{\mathbb{P}\left(d\right),\mathrm{BW}}\right(\mathbb{I},{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-\frac{1}{2}}\left)\right)}^2\mathrm{d}\mathrm{v})}^{\frac{1}{2}}\ge E\mathrm{dist}(\widehat{\mathit{W}},\mathit{W}),$$

(49)

where the last inequality follows from the estimate in Equation (A10) in Lemma A2. Here, E is a positive constant that depends on $\widehat{\mathit{W}}$, and $\mathrm{dist}(·,·)$ is the metric introduced in Definition A2 (Equation (A7)). This metric is a natural one if we restrict ourselves to the set of positive definite tensor fields. The inequality in Equation (49) then implies Equation (48).  □

We note that if we want to investigate the spatially homogeneous steady state $\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]{=}_{\mathrm{def}}[\mathbf{0},\mathbb{I}]$, that is if we set the boundary condition $\mathit{V}=\mathbf{0}$, then

$$C_1=-\frac{1}{\mathrm{Re}},C_2=-\alpha \frac{\Xi }{2\mathrm{Wi}},$$

(50)

and the condition in Equation (47) is automatically satisfied without any restriction of the values of Reynolds number and the Weissenberg number. If the steady state is non-trivial, that is if $\widehat{\mathit{W}}=[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]\ne [\mathbf{0},\mathbb{I}]$, then the condition in Equation (47) must be evaluated. This is done in Section 6 for the Taylor–Couette flow. (Note that $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ and $\widehat{\mathit{v}}$ are solutions to Equation (16), hence they depend on the Weissenberg and Reynolds number.) Naturally, one can expect that the condition will hold for a sufficiently small Reynolds number and Weissenberg number.

Having the Lyapunov type functional given by Equation (32), it is interesting to see how the functional works in the case of close to the equilibrium setting, that is for $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\approx \mathbb{I}$, and for small perturbations, that is for small $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$. We see that if $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ is small and if $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ is close to the identity, then

$$-\ln \det (\mathbb{I}+{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})+\mathrm{Tr}\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)\approx \frac{1}{2}\mathrm{Tr}\left({\left({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right)}^2\right)\approx \frac{1}{2}{\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2,$$

(51)

and the proposed Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ can be approximated as ${\mathcal{V}}_{\mathrm{neq}}\approx {\mathcal{V}}_{\mathrm{naive}}$ where

$${\mathcal{V}}_{\mathrm{naive}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}}){=}_{\mathrm{def}}{\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{4}{\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2\mathrm{d}\mathrm{v}.$$

(52)

The functional ${\mathcal{V}}_{\mathrm{naive}}$ might be the first candidate for the Lyapunov type functional if the stability is investigated using the popular “energy method” (see, for example, Straughan [57]). The functional is clearly nonnegative, and it vanishes if and only if the perturbation vanishes. Moreover, the candidate ${\mathcal{V}}_{\mathrm{naive}}$ for the Lyapunov type functional is much simpler than ${\mathcal{V}}_{\mathrm{neq}}$. Indeed, the proximity of the perturbation to the non-equilibrium state $[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$ is measured using the standard Lebesgue space norms, and ${\mathcal{V}}_{\mathrm{naive}}$ does not depend on the value of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$. Therefore, it seems that ${\mathcal{V}}_{\mathrm{naive}}$ is a good candidate for the Lyapunov type functional for the analysis of arbitrary spatially inhomogeneous non-equilibrium steady state $[\widehat{\mathit{v}},\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}]$.

The inappropriateness of ${\mathcal{V}}_{\mathrm{naive}}$ for the stability analysis is in fact apparent even in a very trivial setting. Let us consider the spatially homogeneous equilibrium steady state $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\mathbb{I}$, $\widehat{\mathit{v}}=\mathbf{0}$ in a mechanically isolated container, that is we set $\mathit{V}=\mathbf{0}$ in the boundary condition in Equation (12). If we use the (exact) evolution equations for the perturbation velocity in Equation (A28b), and if we evaluate the time derivative of ${\mathcal{V}}_{\mathrm{naive}}$, then we get

$$\begin{array}{}& \frac{\mathrm{d}{\mathcal{V}}_{\mathrm{naive}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})={\int }_{\Omega }\tilde{\mathit{v}}•\frac{\partial \tilde{\mathit{v}}}{\partial t}\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}\mathrm{Tr}\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\frac{\partial \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}{\partial t}\right)\mathrm{d}\mathrm{v}\hfill \left(53\right)& \hfill =-{\int }_{\Omega }\frac{2}{\mathrm{Re}}\tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}-{\int }_{\Omega }\Xi \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}-{\int }_{\Omega }\widehat{\mathbb{D}}\tilde{\mathit{v}}•\tilde{\mathit{v}}\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}\mathrm{Tr}\left(\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\frac{\partial \widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}{\partial t}\right)\mathrm{d}\mathrm{v},\end{array}$$

(see also Equation (A34)). The last term on the right-hand side of Equation (53) can be evaluated using the (exact) evolution equation for $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ (see Equation (A29)). Substituting Equation (A29) into Equation (53) and using the fact that $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\mathbb{I}$ and $\widehat{\mathit{v}}=\mathbf{0}$, yields

$$\begin{array}{}& \frac{\mathrm{d}{\mathcal{V}}_{\mathrm{naive}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})=-{\int }_{\Omega }\frac{2}{\mathrm{Re}}\tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}-{\int }_{\Omega }\frac{\Xi }{2}\mathrm{Tr}\left[\right(\{\tilde{\mathit{v}}•\nabla \}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\left)\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right]\mathrm{d}\mathrm{v}+{\int }_{\Omega }\Xi \mathrm{Tr}\left(\tilde{\mathbb{D}}{\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2\right)\mathrm{d}\mathrm{v}\hfill \left(54\right)& \hfill -{\int }_{\Omega }\frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}[\alpha {\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^3+2\alpha {\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2+(1-2\alpha ){\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2]\mathrm{d}\mathrm{v}.\end{array}$$

Using the standard manipulation

$${\int }_{\Omega }\frac{\Xi }{2}\mathrm{Tr}\left[\right(\{\tilde{\mathit{v}}•\nabla \}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\left)\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right]\mathrm{d}\mathrm{v}={\int }_{\partial \Omega }\frac{\Xi }{4}{\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2(\tilde{\mathit{v}}•\mathit{n})\mathrm{d}\mathrm{s}-{\int }_{\Omega }\frac{\Xi }{4}(\mathrm{div}\tilde{\mathit{v}}){\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2\mathrm{d}\mathrm{v},$$

(55)

we see that the second term on the right-hand side of Equation (54) vanishes in virtue of the incompressibility constraint in Equation (A28a) and the boundary condition for $\tilde{\mathit{v}}$. Consequently, Equation (54) reduces to

$$\begin{array}{}& \frac{\mathrm{d}{\mathcal{V}}_{\mathrm{naive}}}{\mathrm{d}t}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})=-{\int }_{\Omega }\frac{2}{\mathrm{Re}}\tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{d}\mathrm{v}+{\int }_{\Omega }\Xi \mathrm{Tr}\left(\tilde{\mathbb{D}}{\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2\right)\mathrm{d}\mathrm{v}\hfill \left(56\right)& \hfill -{\int }_{\Omega }\frac{\Xi }{2\mathrm{Wi}}\mathrm{Tr}[\alpha {\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^3+2\alpha {\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2+(1-2\alpha ){\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2]\mathrm{d}\mathrm{v}.\end{array}$$

Let us now consider an initial perturbation that is chosen is such a way that ${\int }_{\Omega }\Xi \mathrm{Tr}(\tilde{\mathbb{D}}{\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^2)\mathrm{d}\mathrm{v}>0$, which can certainly be done. This positive value will dominate the right-hand side of Equation (56) provided that the Reynolds number and the Weissenberg number are large enough. Consequently, ${\mathcal{V}}_{\mathrm{naive}}$ will (initially) increase, and it would be useless as the Lyapunov type functional unless we a priori limit ourselves to small perturbations.

On the other hand, if we use the functional ${\mathcal{V}}_{\mathrm{neq}}$ in the case $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\mathbb{I}$ and $\widehat{\mathit{v}}=\mathbf{0}$, then we immediately see that the constants $C_1$ and $C_2$ in Equation (40) are negative, and that the equilibrium steady state is asymptotically stable with respect to any perturbation and any value of the Reynolds and the Weissenberg number! (Note that Guillopé and Saut [58] obtained only a conditional stability result in a Sobolev space norm for the equilibrium rest state $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\mathbb{I}$ and $\widehat{\mathit{v}}=\mathbf{0}$ (see their Corollary 3.5 and assumptions of Theorem 3.3).) Based on the analysis presented above, we can therefore claim that we have indeed benefited from a well constructed Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ and the choice of metric. Unlike the naive Lyapunov type functional ${\mathcal{V}}_{\mathrm{naive}}$, the proposed Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ seems to properly reflect the nonlinearity of the governing equations and the related energy storage mechanisms and the entropy production mechanisms.

6. Taylor–Couette Flow

Let us now consider a viscoelastic fluid described by the Giesekus model introduced in Section 2 with $\alpha =\frac{1}{2}$, and let us investigate the stability of steady flow in the standard Taylor–Couette flow geometry (see Figure 1). The objective is to show as how the theory introduced above works in a specific setting. The choice $\alpha =\frac{1}{2}$ is motivated by the simplicity of the expressions for the corresponding steady state.

The fluid is placed in between two infinite concentric cylinders of radii $R_1$ and $R_2$, with $R_1<R_2$. The cylinders are rotating with the angular velocities ${\Omega }_1$ (inner cylinder) and ${\Omega }_2$ (outer cylinder) along the common axis. The geometry naturally leads to the use of cylindrical coordinates $(r,\phi ,z)$; the normed basis vectors are denoted as ${\mathit{g}}_{\widehat{r}}$, ${\mathit{g}}_{\widehat{\phi }}$ and ${\mathit{g}}_{\widehat{z}}$ (see Figure 1). Since the domain is unbounded in the z-direction, we henceforth consider a periodic cell

$$\Omega {=}_{\mathrm{def}}\{(r,\phi ,z)\in {\mathbb{R}}^3|R_1<r<R_2,0\le \phi <2\pi ,\left|z\right|<h\}$$

(57)

where $h>0$ is arbitrary, and we use the notation ${\Gamma }_1{=}_{\mathrm{def}}\{(r,\phi ,z)\in {\mathbb{R}}^3|R_1<r<R_2,0\le \phi <2\pi ,\left|z\right|=h\}$ for the top and bottom base, and ${\Gamma }_2{=}_{\mathrm{def}}\{(r,\phi ,z)\in {\mathbb{R}}^3|r\in \{R_1,R_2\},0\le \phi <2\pi ,\left|z\right|<h\}$ for the cylindrical walls of the domain. The flow is driven by the rotation of the cylinders along the common axis.

6.1. Base Flow—Non-Equilibrium Steady State

The first task in the stability analysis is to find the steady solution to the governing equations. This solution is the spatially inhomogeneous non-equilibrium steady state $\widehat{\mathit{W}}$ as introduced in Section 3.2. The characteristic length and characteristic time have been chosen as $x_{\mathrm{char}}{=}_{\mathrm{def}}{R}_1$, $t_{\mathrm{char}}{=}_{\mathrm{def}}\frac{1}{{\Omega }_1}$. We use the periodic boundary condition on ${\Gamma }_1$ and the no-slip boundary condition for velocity field $\mathit{v}$ on ${\Gamma }_2$, that is ${\mathit{v}|}_{r=R_1}=R_1{\Omega }_1{\mathit{g}}_{\widehat{\phi }}$, ${\mathit{v}|}_{r=R_2}=R_2{\Omega }_2{\mathit{g}}_{\widehat{\phi }}$. These boundary conditions are consistent with the requirements on boundary conditions specified in Section 2.4. In their dimensionless form, the boundary conditions read

$${{\mathit{v}}^{\star }|}_{r^{\star }=1}={\mathit{g}}_{\widehat{\phi }},{{\mathit{v}}^{\star }|}_{r^{\star }=\frac{1}{\eta }}=\frac{\zeta }{\eta }{\mathit{g}}_{\widehat{\phi }},$$

(58)

where we have introduced two dimensionless parameters $\eta {=}_{\mathrm{def}}\frac{R_1}{R_2}$ and $\zeta {=}_{\mathrm{def}}\frac{{\Omega }_2}{{\Omega }_1}$. Hereafter, we work exclusively with the dimensionless variables and thus, for the sake of simplicity, we omit the star denoting them.

Since the problem has the rotational symmetry, we search for the steady non-equilibrium state in a special form. Namely, the solution to Equation (16) subject to boundary conditions in Equation (58) is sought in the form

$$\widehat{\mathit{v}}={\mathrm{v}}^{\widehat{\phi }}\left(r\right){\mathit{g}}_{\widehat{\phi }},\widehat{m}=\widehat{m}\left(r\right),\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\left[\begin{array}{ccc}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}\left(r\right)& {{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}\left(r\right)& 0\\ {{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{r}}\left(r\right)& {{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}\left(r\right)& 0\\ 0& 0& {{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}}\left(r\right)\end{array}\right]$$

(59)

Note that the chosen ansatz for the velocity field automatically satisfies the incompressibility condition. The assumptions lead to the following expressions for the velocity gradient, the symmetric part of the velocity gradient, the convective term, the divergence of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$, and the upper convected derivative of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$:

$$\begin{array}{}& \nabla \widehat{\mathit{v}}=\left[\begin{array}{ccc}0& -\omega & 0\\ r\frac{\mathrm{d}\omega }{\mathrm{d}r}+\omega & 0& 0\\ 0& 0& 0\end{array}\right],\widehat{\mathbb{D}}=\left[\begin{array}{ccc}0& \frac{r}{2}\frac{\mathrm{d}\omega }{\mathrm{d}r}& 0\\ \frac{r}{2}\frac{\mathrm{d}\omega }{\mathrm{d}r}& 0& 0\\ 0& 0& 0\end{array}\right],\frac{\mathrm{d}\widehat{\mathit{v}}}{\mathrm{d}t}=\left[\begin{array}{c}-r{\omega }^2\\ 0\\ 0\end{array}\right],\hfill \left(60\right)& \hfill \mathrm{div}\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}=\left[\begin{array}{c}\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}(r{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}})-\frac{{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}}{r}\\ \frac{\mathrm{d}{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{r}}}{\mathrm{d}r}+\frac{{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{r}}+{{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}}{r}\\ 0\end{array}\right],\stackrel{▿}{\overline{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}}=\left[\begin{array}{ccc}0& -r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}& 0\\ -r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}& -2r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{r}}& 0\\ 0& 0& 0\end{array}\right],\end{array}$$

where we introduce the angular velocity $\omega \left(r\right)$, ${\mathrm{v}}^{\widehat{\phi }}\left(r\right){=}_{\mathrm{def}}\omega \left(r\right)r$. Using Equation (60), we see that the governing equations for the velocity field in Equation (16b) reduce to

$$\left[\begin{array}{c}-r{\omega }^2\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}r}(\widehat{m}+\Xi ({{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}-\frac{1}{3}({{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}+{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}+{{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}})\left)\right)+\Xi \frac{{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}-{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}}{r}\\ \frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}(\frac{1}{\mathrm{Re}}{r}^3\frac{\mathrm{d}\omega }{\mathrm{d}r}+\Xi r^2\mathrm{B}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }})\\ 0\end{array}\right],$$

(61a)

while the governing equations in Equation (16c) for $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ read

$$\begin{array}{}& \left[\begin{array}{ccc}0& -r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}& 0\\ -r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}& -2r\frac{\mathrm{d}\omega }{\mathrm{d}r}{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{r}}& 0\\ 0& 0& 0\end{array}\right]=\hfill \left(61\mathrm{b}\right)& -\frac{1}{\mathrm{Wi}}\left[\begin{array}{ccc}\alpha ({({{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}})}^2+{({{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }})}^2)+(1-2\alpha ){{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}-(1-\alpha )& \alpha {{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}({{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}+{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }})+(1-2\alpha ){{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}& 0\\ \alpha {{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}({{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}+{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }})+(1-2\alpha ){{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}& \alpha ({({{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }})}^2+{({{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }})}^2)+(1-2\alpha ){{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}-(1-\alpha )& 0\\ 0& 0& \alpha {({{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}})}^2+(1-2\alpha ){{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}}-(1-\alpha )\end{array}\right].\hfill \end{array}$$

Assuming that $\frac{\mathrm{d}\omega }{\mathrm{d}r}\ne 0$ in $(R_1,R_2)$, Equation (61b) can be solved for ${{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}$, ${{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}$, ${{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}$ and ${{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}}$. However, for general $\alpha \in (0,1)$, the formulae for the aforementioned quantities are too complex to be written down here. Let us note however that for $\alpha =\frac{1}{2}$ the formulae simplify significantly; the solution to Equation (61b) which satisfies the condition of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ being positive definite in this case reads

$${{\mathrm{B}}^{\widehat{z}}}_{\widehat{z}}=1,{{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}=\frac{-1+\sqrt{1+c^2}}{c},{{\mathrm{B}}^{\widehat{r}}}_{\widehat{r}}=\sqrt{\frac{2}{c}{{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}},{{\mathrm{B}}^{\widehat{\phi }}}_{\widehat{\phi }}=\sqrt{2(c+\frac{1}{c}){{\mathrm{B}}^{\widehat{r}}}_{\widehat{\phi }}},$$

(62)

where we denote $c{=}_{\mathrm{def}}2\mathrm{Wi}r\frac{\mathrm{d}\omega }{\mathrm{d}r}$. Substituting Equation (62) into the second equation in Equation (61a) then yields an ordinary differential equation for the angular velocity $\omega$

$$0=\frac{\mathrm{d}}{\mathrm{d}r}(\frac{1}{\mathrm{Re}}{r}^3\frac{\mathrm{d}\omega }{\mathrm{d}r}+\Xi r^2\frac{-1+\sqrt{1+4{\mathrm{Wi}}^2{r}^2{\left(\frac{\mathrm{d}\omega }{\mathrm{d}r}\right)}^2}}{2\mathrm{Wi}r\frac{\mathrm{d}\omega }{\mathrm{d}r}}),$$

(63)

supplemented by the boundary conditions ${\omega |}_{r=1}=1$, ${\omega |}_{r=\frac{1}{\eta }}=\zeta$, which follow from Equation (58) and the fact that ${\mathrm{v}}^{\widehat{\phi }}\left(r\right)=\omega \left(r\right)r$. Equation (63) together with the boundary conditions constitute a boundary value problem which needs to be solved numerically.

6.2. Explicit Criterion for the Stability of Spatially Inhomogeneous Non-Equilibrium Steady State

Here, we explicitly compute constants $C_1$, $C_2$ defined by Equations (41a) and (41b) for the Taylor–Couette problem and for the specific values of the dimensionless numbers $\Xi$, $\mathrm{Re}$ and $\mathrm{Wi}$. Let us recall that, for the sake of simplicity, we set $\alpha =\frac{1}{2}$, and we consider the steady tensor field $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ given by Equation (62). We fix the values for the geometric parameter $\eta$ and angular velocities ratio $\zeta$ as $\eta =\frac{1}{2}$ and $\zeta =2$.

The angular velocity $\omega$ is obtained by solving Equation (63) which is a boundary-value problem for a second order nonlinear differential equation. The problem was solved numerically using the function solve_bvp from SciPy library version 1.0.0, which implements a fourth-order collocation algorithm with the control of residuals as described in Kierzenka and Shampine [59]. With the angular velocity $\omega$ in hand, we immediately get the steady velocity field $\widehat{\mathit{v}}=\omega \left(r\right)r{\mathit{g}}_{\widehat{\phi }}$, and the steady left Cauchy–Green tensor field $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ through Equation (62). The plots of the velocity field and the components of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ are shown in Figure 2.

Having computed the steady velocity field $\widehat{\mathit{v}}$ and the corresponding steady field $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$, we can evaluate the constants $C_1$ and $C_2$ in the estimate in Equation (40). The gradient of $\widehat{\mathit{v}}$ as well as the gradient of $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$ are again computed numerically from the obtained numerical solution. Finally, the Poincaré constant for the cylindrical annulus is determined via an explicit solution of the corresponding eigenvalue problem $-\Delta u=\lambda u$ for the Laplace operator with Dirichlet boundary condition, which leads, for the geometrical parameter $\eta =\frac{1}{2}$, to the value $C_P\approx 0.1025.$ The resulting stability regions in the $\mathrm{Re}$–$\mathrm{Wi}$ plane are shown in Figure 3 for a fixed value of the dimensionless shear modulus $\Xi$. As one might expect the spatially inhomogeneous steady state is indeed unconditionally asymptotically stable if the Weissenberg number $\mathrm{Wi}$ and the Reynolds number $\mathrm{Re}$ are small enough.

6.3. Numerical Experiments—Evolution of Various Initial Perturbations

Finally, we document the theoretically predicted behaviour by numerical experiments. The numerical experiments allow us to quantitatively track the evolution of key quantities such as the net kinetic energy, and also to quantitatively monitor the energy exchange between the fluid and its surroundings.

The governing equations were numerically solved using standard techniques. The weak forms of the governing equations were discretised in the space using the finite element method, while the time derivatives were approximated with the backward Euler method. The two-dimensional domain $\Omega$ was discretised by regular quadrilaterals. The mesh divided the annular region $\Omega$ into 80 pieces in the radial direction, and in 720 pieces in the azimuthal direction. The corresponding total number of degrees of freedom in all numerical experiments was over $1.3\times {10}^6$. The velocity field $\mathit{v}$ and the ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ field were approximated by biquadratic $Q2$ elements, and the pressure field m was approximated by the piecewise linear discontinuous elements ${P1}^d$ (see Korelc and Wriggers [60] for details). The finite element pair that was used for the velocity/pressure fields satisfied the Babuška-Brezzi condition, the finite element for ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ field was chosen to be the same as for the velocity in order to provide rich enough finite element space for the solution. The same finite elements were chosen for the two-dimensional simulation of other viscoelastic rate-type fluids (Oldroyd, Burgers and their various nonlinear versions) by Hron et al. [61] and Málek et al. [62]. In the three-dimensional case, low order elements can be used to decrease the overall cost of the calculation (see Tůma et al. [63]).

The numerical scheme was implemented in the AceGen/AceFEM system (see Korelc [64] and Korelc [65]). The main advantage of the system is that it provides automatic differentiation used for the computation of the exact tangent matrix needed by the Newton solver that treats all nonlinearities. The final set of linear equations was solved by the direct solver Intel MKL Pardiso. The stopping criterion for the Newton solver was set to ${10}^{-9}$.

Using the numerical scheme, we are ready to study the behaviour of various perturbations to the non-equilibrium steady state. In all scenarios described below, we use the dimensionless parameters

$$\Xi =0.1,\mathrm{Re}=50,\mathrm{Wi}=5,\alpha =\frac{1}{2}$$

(64)

and we fix the geometric parameter $\eta$ and angular velocities ratio $\zeta$ as in Section 6.2, that is $\eta =\frac{1}{2}$ and $\zeta =2$. The chosen values of $\eta$, $\zeta$ and $\Xi$ correspond to the stability diagram shown in Figure 3a. The values of Reynolds number and Weissenberg number are outside the region where we have proven the decay of the proposed Lyapunov type functional. Nevertheless, as we show below, the Lyapunov type functional is, in the cases being investigated below, still a decreasing function.

First, we start from the homogeneous steady state solution $[\mathit{v},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}},m]=[\mathbf{0},\mathbb{I},0]$, and we let the system to spontaneously evolve up to the time instant $t=1000$. (More precisely, the initial condition is $\mathit{v}=\mathbf{0}$ inside the domain $\Omega$, and Equation (17) holds on the boundary of $\Omega$. After the first computational time step, which is chosen as $\Delta t=0.05$, we get on the discrete level a divergence-free velocity field with the appropriate boundary condition. This discrete velocity field provides us a consistent initial condition for further computations. Therefore, we formally start the evolution not at $t=0$, but at $t=0.05$.) At this time instant, the system is almost relaxed and is close to the steady solution. The solution at $t=1000$ is used as a starting point for solving the steady governing equations (without the time derivatives) and the spatially inhomogeneous non-equilibrium steady state is obtained just in two Newton iterations. (The finite element solution coincides with the semi-analytical steady solution obtained in Section 6.1. This among others provides us a tool for the code verification.) Consequently, the finite element solution is in what follows used as the spatially inhomogeneous non-equilibrium steady state $\widehat{\mathit{W}}$.

Having obtained the numerical representation of the spatially inhomogeneous non-equilibrium steady state, we proceed with two scenarios concerning the specification of the initial perturbation.

6.3.1. Scenario A—Localised Perturbation of the Left Cauchy–Green Field

In the first scenario, we keep the initial velocity field perturbation equal to zero,

$${\tilde{\mathit{v}}|}_{t=0}=\mathbf{0},$$

(65)

while the initial perturbation in ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ is localised in space (see the first snapshot in Figure 4). Since the system is fully coupled, the perturbation in the ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ field triggers for $t>0$ a nontrivial evolution of the velocity perturbation $\tilde{\mathit{v}}$ (see Figure 5). This can be observed also in the plots showing the evolution of the net elastic stored energy and the net kinetic energy (see Figure 6).

Finally, we also investigate the time evolution of the proposed Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ and the naive Lyapunov type functional ${\mathcal{V}}_{\mathrm{naive}}$, and the net mechanical energy flux going through the boundary of $\Omega$ (see Figure 6). Although we work with the Reynolds number/Weissenberg number pair outside the guaranteed stability region, we see that the value of Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ still decreases in time, and that the perturbation vanishes for $t\to +\infty$. This indicates that the estimates on the time derivative of the proposed Lyapunov type functional are, at least for a class of perturbations, too strict and they might be improved. One should also note that the “net kinetic energy” of the perturbation, that is the functional ${\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}$, does not decrease for all $t>0$ (see Figure 6b). In fact, it experiences a transitional growth, and such a transient growth can be observed even for the Reynolds number/Weissenberg number values within the stability region. This is a natural observation. The elastic energy stored in the material can be released in the form of the kinetic energy. It is only the combination of the elastic energy and the kinetic energy that appears in the Lyapunov type functional that leads to a quantity that decays at any time.

Further, the net mechanical energy flux through the boundary fluctuates around the value that corresponds to the non-equilibrium steady state, and then it reaches the value that corresponds to the spatially inhomogeneous non-equilibrium steady state (see Figure 6d). This can again happen even if the Reynolds number/Weissenberg number take values within the stability region.

6.3.2. Scenario B—Global Perturbation of the Velocity Field

In the second scenario, we start with a nonzero velocity perturbation $\tilde{\mathit{v}}$, and the ${\mathbb{B}}_{{\kappa }_{p\left(t\right)}}$ field is kept unchanged,

$${\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}|}_{t=0}=\mathbb{O}.$$

(66)

The initial velocity $\mathit{v}$ is prescribed as

$${\mathit{v}|}_{t=0}=\Omega r{\mathit{g}}_{\widehat{\phi }},$$

(67)

where the angular velocity is the arithmetic mean of the two angular velocities $\Omega =\frac{{\Omega }_1+{\Omega }_2}{2}$. (Formally, we apply the same procedure as discussed in the previous section. The initial condition is $\mathit{v}=\Omega r{\mathit{g}}_{\widehat{\phi }}$ inside the domain $\Omega$, and Equation (17) holds on the boundary of $\Omega$. The actual computation starts after the first (formal) time step, when the discrete velocity field is divergence-free and it fulfills the boundary condition.) Again, as in the previous case, the non-zero perturbation in one unknown field ($\tilde{\mathit{v}}$) triggers for $t>0$ a nontrivial evolution of the other unknown field ($\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$) (see Figure 7 and Figure 8).

Moreover, this numerical experiment is instructive for yet another reason. In Figure 9c, we plot the time evolution of the values of the functionals ${\mathcal{V}}_{\mathrm{neq}}$ and ${\mathcal{V}}_{\mathrm{naive}}$. Clearly, the functional ${\mathcal{V}}_{\mathrm{naive}}$ (see Equation (52)), which is a naive candidate for the Lyapunov type functional, experiences a transitional growth. Interestingly, the proposed complex Lyapunov type functional ${\mathcal{V}}_{\mathrm{neq}}$ is still a decreasing function, although the Reynolds number/Weissenberg number values are outside the region, where we have actually proven the decay of the functional. This further indicates that the functional ${\mathcal{V}}_{\mathrm{naive}}$ is indeed not a good candidate for a Lyapunov type functional (see also Section 5 for further discussion).

7. Conclusions

We have investigated the stability of spatially inhomogeneous non-equilibrium steady states (flows) of viscoelastic fluids described by the Giesekus model. We have derived bounds on the values of the Reynolds number and the Weissenberg number that guarantee the flow stability subject to any finite perturbation. The stability has been investigated using the Lyapunov type functional given by the formula

$${\mathcal{V}}_{\mathrm{neq}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})={\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{2}[-\ln \det (\mathbb{I}+{\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}})+\mathrm{Tr}({\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}}^{-1}\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\left)\right]\mathrm{d}\mathrm{v}.$$

(68)

A few observations concerning the proposed Lyapunov type functional are at hand.

First, the proposed Lyapunov type functional has a relatively complicated form. In particular, it is not quadratic in the perturbation $\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$, and it depends on the spatially inhomogeneous non-equilibrium state $\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}$. This makes it remarkably different from a naive Lyapunov type functional of the form

$${\mathcal{V}}_{\mathrm{naive}}(\tilde{\mathit{W}}\parallel \widehat{\mathit{W}})={\int }_{\Omega }\frac{1}{2}{\left|\tilde{\mathit{v}}\right|}^2\mathrm{d}\mathrm{v}+{\int }_{\Omega }\frac{\Xi }{4}{\left|\widetilde{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}}\right|}^2\mathrm{d}\mathrm{v},$$

(69)

which might be a first try if the stability problem were analysed using the popular “energy method”. However, as we have shown, the complicated structure of the proposed functional ${\mathcal{V}}_{\mathrm{neq}}$ leads to a relatively simple and well structured formula for its time derivative, which in turn allows one to formulate conditions that guarantee the negativity of the time derivative. Furthermore, the complicated structure of the proposed functional ${\mathcal{V}}_{\mathrm{neq}}$ also leads to a simple relation between the functional and the metric

$$\mathrm{dist}(\widehat{\mathit{W}},\mathit{W}){=}_{\mathrm{def}}{({\parallel \widehat{\mathit{v}}-\mathit{v}\parallel }_{L^2(\Omega )}^2+{\left[{\mathrm{dist}}_{{\mathbb{P}}_{\Omega }\left(d\right),\mathrm{BW}}\right(\widehat{{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}},{\mathbb{B}}_{{\kappa }_{p\left(t\right)}}\left)\right]}^2)}^{\frac{1}{2}}$$

(70)

on the set of spatially distributed symmetric positive definite matrices.

Second, the Lyapunov type functional has been used in the investigation of stability of solution to the complete system of nonlinear governing equations. In particular, the evolution equations for the perturbation have been investigated without any simplification. This makes the present approach different from the “energy budget” analysis (see, for example, Joo and Shaqfeh [66], Byars et al. [67], Ganpule and Khomami [68], Smith et al. [69], Karapetsas and Tsamopoulos [70], Pettas et al. [71], and especially the work by Grillet et al. [19] who investigated the Giesekus model). The “energy budget” analysis, although valuable in the discussion of the nature of the instability mechanisms, is based on the linearised momentum equation for the perturbation and linearised constitutive equation for the “polymeric stress”. Consequently, the standard “energy budget” analysis is, unlike the present approach, a tool that cannot be used in the finite amplitude stability analysis of the complete system of nonlinear governing equations. One might also note that, despite the complexity of the proposed Lyapunov type functional, the formula for its time derivative is in fact quite simple compared to the formulae in the “energy budget” analysis. This happens even though the “energy budget” formulae paradoxically stem from various simplifications of the original system of governing equations.

Third, the Lyapunov type functional has been designed using thermodynamical arguments. In fact, the proposed Lyapunov type functional has been constructed from the net mechanical energy functional $E_{\mathrm{mech}}$ (see Equation (10), via Equation (30)). This makes the construction quite general, and one might speculate that a similar approach is very likely applicable to other popular viscoelastic rate-type models such as the PTT model (see Phan Thien and Tanner [72]) or the FENE-P model (see Bird et al. [73]), as well as complex viscoelastic rate-type models with, for example, stress diffusion terms (see Málek et al. [37] or Dostalík et al. [39]). Further, the construction of the Lyapunov type functional has been based on the method proposed by Bulíček et al. [23], and this method is speculated to work even for complex coupled thermo-mechanical systems. This naturally calls for the investigation of the applicability of the method in more complex settings such as flows of viscoelastic rate-type fluids with temperature dependent material parameters.

Fourth, thermodynamical type considerations such as the identification of the energy storage mechanisms and entropy producing mechanisms are known to play an important role in the rigorous mathematical theory of nonlinear partial differential equations governing the motion of viscoelastic fluids (see, for example, Hu and Lelièvre [74], Boyaval et al. [40], Barrett and Boyaval [75], Barrett and Boyaval [76], Barrett and Süli [47] or Bulíček et al. [77]). On the other hand, rigorous mathematical analysis of long-time behaviour of viscoelastic fluids is usually done without a direct appeal to thermodynamics, and the available results are quite limited especially if one considers thermodynamically open systems (see, for example, Guillopé and Saut [78], Nohel and Pego [79], Jourdain et al. [80] or Renardy [81]). (Usually, only stability of unidirectional steady flows in simple geometries is considered.) Consequently, the approach proposed in the current contribution might be of interest from the rigorous mathematical perspective as well. This means that one should deal with the weak solution to the governing equations, and that one should investigate the applicability of the presented arguments for a solution/perturbation that has only the smoothness that can be actually proven.