On diffusive variants of some classical viscoelastic rate-type models

TL;DR

This paper develops a thermodynamically consistent framework for adding stress diffusion to classical viscoelastic models, enabling more accurate modeling of complex fluid flows with stress diffusion effects.

Contribution

It introduces a method to incorporate stress diffusion into standard viscoelastic models while maintaining thermodynamic consistency, expanding their applicability.

Findings

Models with stress diffusion obey thermodynamic laws.

Framework applicable to multiple classical viscoelastic models.

Potential for improved flow simulations of complex fluids.

Abstract

We present a thermodynamically based approach to the design of models for viscoelastic fluids with stress diffusion effect. In particular, we show how to add a stress diffusion term to some standard viscoelastic rate-type models (Giesekus, FENE-P, Johnson--Segalman, Phan-Thien--Tanner and Bautista--Manero--Puig) so that the resulting models with the added stress diffusion term are thermodynamically consistent in the sense that they obey the first and the second law of thermodynamics. We point out the potential applications of the provided thermodynamical background in the study of flows of fluids described by the proposed models.

Figures

Tables

Table 1. Table 1: Mechanical part of Helmholtz free energy ψ m subscript 𝜓 𝑚 \psi_{m} and contribution of mechanical quantities to entropy production ξ m subscript 𝜉 𝑚 \xi_{m} for some viscoelastic rate-type models with a stress diffusion term. If the stress diffusion coefficient μ ~ ~ 𝜇 \tilde{\mu} is equal to zero, then the given Helmholtz free energy/entropy production pair leads to the corresponding classical model. The Helmholtz free energy ψ 𝜓 \psi and entropy production ξ 𝜉 \xi are obtained from the tabulated formulae via ψ = def − c V NSE ​ θ ​ ( ln ⁡ θ θ ref − 1 ) + ψ m subscript def 𝜓 superscript subscript 𝑐 V NSE 𝜃 𝜃 subscript 𝜃 ref 1 subscript 𝜓 𝑚 \psi=_{\mathrm{def}}-c_{\mathrm{V}}^{\mathrm{NSE}}\theta\left(\ln\frac{\theta}{\theta_{\mathrm{ref}}}-1\right)+\psi_{m} and ξ = ξ m + κ ​ | ∇ θ | 2 θ 2 𝜉 subscript 𝜉 𝑚 𝜅 superscript ∇ 𝜃 2 superscript 𝜃 2 \xi=\xi_{m}+\kappa\frac{\left|\nabla\theta\right|^{2}}{\theta^{2}} .

Helmholtz free energy ${\psi }_m$ (mechanical part)	Entropy production ${\xi }_m$ (contribution of mechanical quantities)
Giesekus
$\frac{\mu }{2\rho }\left[\mathrm{Tr}{𝔹}_2-3-\ln det{𝔹}_2\right]$	$\frac{1}{\theta }\{2\nu {𝔻}_{\delta }:{𝔻}_{\delta }+\frac{{\mu }^2}{2{\nu }_1}\mathrm{Tr}\left[\lambda {𝔹}_2^2+\left(1-3\lambda \right){𝔹}_2+\left(1+\lambda \right){𝔹}_2^{-1}+\left(3\lambda -2\right)𝕀\right]+\frac{\mu \tilde{\mu }}{2{\nu }_1}\nabla {𝔹}_2⋮\nabla {𝔹}_2\}$
FENE-P
$\frac{\mu }{2\rho }\left[-b\ln \left(1-\frac{1}{b}\mathrm{Tr}{𝔹}_2\right)-3-\ln det{𝔹}_2\right]$	$\frac{1}{\theta }\{2\nu {𝔻}_{\delta }:{𝔻}_{\delta }+\frac{{\mu }^2}{2{\nu }_1}\left[{\left(1-\frac{1}{b}\mathrm{Tr}{𝔹}_2\right)}^{-2}\mathrm{Tr}{𝔹}_2-6{\left(1-\frac{1}{b}\mathrm{Tr}{𝔹}_2\right)}^{-1}+\mathrm{Tr}{𝔹}_2^{-1}\right]+\frac{\mu \tilde{\mu }}{2{\nu }_1}\nabla \left[{\left(1-\frac{1}{b}\mathrm{Tr}{𝔹}_2\right)}^{-1}{𝔹}_2\right]⋮\nabla \left[{\left(1-\frac{1}{b}\mathrm{Tr}{𝔹}_2\right)}^{-1}{𝔹}_2\right]\}$
Johnson–Segalman
$\frac{\mu }{2\rho }\left[\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}-3-\ln det{𝔹}_{2,\mathrm{GS}}\right]$	$\frac{1}{\theta }\{2\nu {𝔻}_{\delta }:{𝔻}_{\delta }+\frac{{\mu }^2}{2{\nu }_1}\left[\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}+\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}^{-1}-6\right]+\frac{\mu \tilde{\mu }}{2{\nu }_1}\nabla {𝔹}_{2,\mathrm{GS}}⋮\nabla {𝔹}_{2,\mathrm{GS}}\}$
Phan–Thien–Tanner
$\frac{\mu }{2\rho }\left[\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}-3-\ln det{𝔹}_{2,\mathrm{GS}}\right]$	$\frac{1}{\theta }\{2\nu {𝔻}_{\delta }:{𝔻}_{\delta }+\frac{{\mu }^2}{2{\nu }_1}f\left({𝔹}_{2,\mathrm{GS}}\right)\left[\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}+\mathrm{Tr}{𝔹}_{2,\mathrm{GS}}^{-1}-6\right]+\frac{\mu \tilde{\mu }}{2{\nu }_1}\nabla {𝔹}_{2,\mathrm{GS}}⋮\nabla {𝔹}_{2,\mathrm{GS}}\}$
Bautista–Manero–Puig
$\frac{\mu }{2\rho }\left[\mathrm{Tr}{𝔹}_2-3-\ln det{𝔹}_2\right]+\frac{\chi }{2\rho }{\left({\phi }_0-\phi \right)}^2$	$\frac{1}{\theta }\{2\nu {𝔻}_{\delta }:{𝔻}_{\delta }+\frac{{\mu }^2\phi }{2}\left[\mathrm{Tr}{𝔹}_2+\mathrm{Tr}{𝔹}_2^{-1}-6\right]+\frac{\beta \mu \tilde{\mu }}{2}\nabla {𝔹}_2⋮\nabla {𝔹}_2+\chi \frac{{\left({\phi }_0-\phi \right)}^2}{\tau }\}$

Equations

$$\mathbb{F}={\mathbb{F}}_2{\mathbb{F}}_1,$$

$$\frac{\mathrm{d}\mathbb{F}}{\mathrm{d}t}=\mathbb{LF}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{=}_{\mathrm{def}}\frac{\partial }{\partial t}+\mathbf{v}\bullet \nabla ,$$

$$\mathbb{L}$$

$$\mathbb{D}$$

$$\mathbb{W}$$

$${\mathbb{C}}_2{=}_{\mathrm{def}}{\mathbb{F}}_2^{\top }{\mathbb{F}}_2,{\mathbb{B}}_2{=}_{\mathrm{def}}{\mathbb{F}}_2{\mathbb{F}}_2^{\top }.$$

$$\mathbb{L}={\mathbb{L}}_2+{\mathbb{F}}_2{\mathbb{L}}_1{\mathbb{F}}_2^{-1},$$

$$\text{accentset}▽\overline{{\mathbb{B}}_2}=-2{\mathbb{F}}_2{\mathbb{D}}_1{\mathbb{F}}_2^{\top },$$

$$\text{accentset}▽\mathbb{A}{=}_{\mathrm{def}}\frac{\mathrm{d}\mathbb{A}}{\mathrm{d}t}-\mathbb{L}\mathbb{A}-\mathbb{A}{\mathbb{L}}^{\top }.$$

$$\text{accentset}\text{medsquare}\mathbb{A}{=}_{\mathrm{def}}\frac{\mathrm{d}\mathbb{A}}{\mathrm{d}t}-a\left(\mathbb{D}\mathbb{A}+\mathbb{A}\mathbb{D}\right)-\left(\mathbb{W}\mathbb{A}+\mathbb{A}{\mathbb{W}}^{\top }\right),$$

$$\text{accentset}\text{medsquare}\mathbb{A}=\frac{\mathrm{d}\mathbb{A}}{\mathrm{d}t}-{\mathbb{L}}_{\mathrm{GS}}\mathbb{A}-\mathbb{A}{\mathbb{L}}_{\mathrm{GS}}^{\top },$$

$${\mathbb{L}}_{\mathrm{GS}}{=}_{\mathrm{def}}a\mathbb{D}+\mathbb{W}=\mathbb{L}+(a-1)\mathbb{D}.$$

$$\frac{\mathrm{d}{\mathbb{F}}_{\mathrm{GS}}}{\mathrm{d}t}={\mathbb{L}}_{\mathrm{GS}}{\mathbb{F}}_{\mathrm{GS}},{{\mathbb{F}}_{\mathrm{GS}}|}_{t=0}=\mathbb{I}.$$

$${\mathbb{L}}_{2,\mathrm{GS}}{=}_{\mathrm{def}}a{\mathbb{D}}_2+{\mathbb{W}}_2={\mathbb{L}}_2+(a-1){\mathbb{D}}_2,$$

$$\frac{\mathrm{d}{\mathbb{F}}_{2,\mathrm{GS}}}{\mathrm{d}t}={\mathbb{L}}_{2,\mathrm{GS}}{\mathbb{F}}_{2,\mathrm{GS}},{{\mathbb{F}}_{2,\mathrm{GS}}|}_{t=0}=\mathbb{I}.$$

$${\mathbb{C}}_{2,\mathrm{GS}}{=}_{\mathrm{def}}{\mathbb{F}}_{2,\mathrm{GS}}^{\top }{\mathbb{F}}_{2,\mathrm{GS}},{\mathbb{B}}_{2,\mathrm{GS}}{=}_{\mathrm{def}}{\mathbb{F}}_{2,\mathrm{GS}}{\mathbb{F}}_{2,\mathrm{GS}}^{\top }.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\det {\mathbb{F}}_{\mathrm{GS}}=\left(\det {\mathbb{F}}_{\mathrm{GS}}\right)\mathrm{Tr}\left(\frac{\mathrm{d}{\mathbb{F}}_{\mathrm{GS}}}{\mathrm{d}t}{\mathbb{F}}_{\mathrm{GS}}^{-1}\right)=\left(\det {\mathbb{F}}_{\mathrm{GS}}\right)\mathrm{Tr}{\mathbb{L}}_{\mathrm{GS}}=a\left(\det {\mathbb{F}}_{\mathrm{GS}}\right)\mathrm{Tr}\mathbb{L}.$$

$${\mathbb{F}}_{1,\mathrm{GS}}$$

$${\mathbb{L}}_{1,\mathrm{GS}}$$

$${\mathbb{F}}_{\mathrm{GS}}={\mathbb{F}}_{2,\mathrm{GS}}{\mathbb{F}}_{1,\mathrm{GS}},$$

$${\mathbb{L}}_{\mathrm{GS}}={\mathbb{L}}_{2,\mathrm{GS}}+{\mathbb{F}}_{2,\mathrm{GS}}{\mathbb{L}}_{1,\mathrm{GS}}{\mathbb{F}}_{2,\mathrm{GS}}^{-1},$$

$$\text{accentset}\text{medsquare}\overline{{\mathbb{B}}_{2,\mathrm{GS}}}=-2{\mathbb{F}}_{2,\mathrm{GS}}{\mathbb{D}}_{1,\mathrm{GS}}{\mathbb{F}}_{2,\mathrm{GS}}^{\top },$$

$${\mathbb{D}}_{1,\mathrm{GS}}{=}_{\mathrm{def}}\frac{1}{2}\left({\mathbb{L}}_{1,\mathrm{GS}}+{\mathbb{L}}_{1,\mathrm{GS}}^{\top }\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\ln \det {\mathbb{B}}_{2,\mathrm{GS}}$$

$${\mathbb{C}}_{2,\mathrm{GS}}:{\mathbb{D}}_{1,\mathrm{GS}}$$

$$\mathbb{I}:{\mathbb{D}}_{1,\mathrm{GS}}$$

$$\rho \frac{\mathrm{d}e}{\mathrm{d}t}=\mathbb{T}:\mathbb{D}-\mathrm{div}{\mathbf{j}}_e,$$

$$\frac{\mathrm{d}e}{\mathrm{d}t}=\theta \frac{\mathrm{d}\eta }{\mathrm{d}t}+\frac{\partial \psi }{\partial \mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}}\right)+\frac{\partial \psi }{\partial \ln \det {\mathbb{B}}_{2,\mathrm{GS}}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\ln \det {\mathbb{B}}_{2,\mathrm{GS}}\right)+\frac{\partial \psi }{\partial \phi }\frac{\mathrm{d}\phi }{\mathrm{d}t},$$

$$\psi ={\psi }_0(\theta )+{\psi }_1(\mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}},\ln \det {\mathbb{B}}_{2,\mathrm{GS}})+{\psi }_2(\phi ).$$

$${\psi }_0=-c_{\mathrm{V}}^{\mathrm{NSE}}\theta \left(\ln \frac{\theta }{{\theta }_{\mathrm{ref}}}-1\right),$$

$$\begin{gathered} \rho \frac{\mathrm{d}\eta }{\mathrm{d}t}+\mathrm{div}\left(\frac{{\mathbf{j}}_e}{\theta }\right)=\frac{1}{\theta }\{\left[{\mathbb{T}}_{\delta }-2\rho \frac{\partial \psi }{\partial \mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}}}a{\left({\mathbb{B}}_{2,\mathrm{GS}}\right)}_{\delta }\right]:{\mathbb{D}}_{\delta } \\ -\frac{1}{2}\mathrm{Tr}\left[\text{accentset}\text{medsquare}\overline{{\mathbb{B}}_{2,\mathrm{GS}}}\left(2\rho \frac{\partial \psi }{\partial \mathrm{Tr}{\mathbb{B}}_{2,\mathrm{GS}}}\mathbb{I}+2\rho \frac{\partial \psi }{\partial \ln \det {\mathbb{B}}_{2,\mathrm{GS}}}{\mathbb{B}}_{2,\mathrm{GS}}^{-1}\right)\right]-\rho \frac{\partial \psi }{\partial \phi }\frac{\mathrm{d}\phi }{\mathrm{d}t}-\frac{{\mathbf{j}}_e\bullet \nabla \theta }{\theta }\}, \end{gathered}$$

Full text

The manuscript has been published as Dostalík, M. and Průša, V. and Skřrivan, T.: On diffusive variants of some classical viscoelastic rate-type models, AIP Conference Proceedings, 2107, 1, 020002, 2019, doi:10.1063/1.5109493. Unfortunately, the published version of this manuscript contains several misprints.

•

Misprint in the specification of the additive constant in the Helmholtz free energy ansatz for the FENE-P model. The correct treatment of the FENE-P is discussed at the end of the manuscript is section “Helmholtz free energy for the FENE-P model”.

•

Misprint in the specification of the entropy production for the Giesekus model. The correct treatment of the entropy production for the Giesekus model is discussed at the end of the manuscript is section “Entropy production for the Giesekus model”.

We are fully responsible for the misprints.

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Kinematics

4. 2.1.1 Models with upper convected derivative

5. 2.1.2 Models with Gordon–Schowalter derivative

6. 2.2 Evolution equation for entropy

7. 2.3 Evolution equation for temperature

8. 3 Giesekus model

9. 3.1 Helmholtz free energy

10. 3.2 Entropy production

11. 3.3 Constitutive relations

12. 3.4 Governing equations

13. 4 FENE-P model

14. 4.1 Helmholtz free energy

15. 4.2 Entropy production

16. 4.3 Constitutive relations

17. 4.4 Governing equations

18. 5 Johnson–Segalman model

19. 5.1 Helmholtz free energy

20. 5.2 Entropy production

21. 5.3 Constitutive relations

22. 5.4 Governing equations

23. 6 Phan–Thien–Tanner model

24. 6.1 Helmholtz free energy

25. 6.2 Entropy production

26. 6.3 Constitutive relations

27. 6.4 Governing equations

28. 7 Bautista–Manero–Puig type model

29. 7.1 Helmholtz free energy

30. 7.2 Entropy production

31. 7.3 Constitutive relations

32. 7.4 Governing equations

33. 8 Implications of the thermodynamical background

34. 8.1 Boundary conditions – characterisation of thermodynamically isolated or open systems

35. 8.2 Evolution equation for temperature – complete system of governing equations for complex processes

36. 8.3 Qualitative behaviour – conserved quantities and stability analysis

37. 9 Conclusion

38. A Corrections

39. A.1 Helmholtz free energy for the FENE-P model

40. A.2 Entropy production for the Giesekus model

1 Introduction

The behaviour of complex fluids such as wormlike micellar solutions is on the macroscopic level often described using viscoelastic rate-type models with an added stress diffusion term, see Cates and Fielding (2006), Subbotin et al. (2011) or Fardin et al. (2015). The derivation of most of these phenomenological models usually proceeds in a rather ad hoc manner. One takes a well established viscoelastic rate-type model, and adds a stress diffusion term to the evolution equation for the extra stress tensor, while the stress diffusion term takes the form of Laplace operator acting on the extra stress tensor. In this way one obtains the diffusive Johnson–Segalman model, see for example Olmsted et al. (2000), out of the standard Johnoson–Segalman model, see Johnson and Segalman (1977). Similarly the diffusive Giesekus model, see for example Helgeson et al. (2009a, b) or Cheng et al. (2017), is built out of the standard Giesekus model, see Giesekus (1982), or the diffusive Rolie–Poly model, see for example Adams et al. (2011), Chung et al. (2011) or Carter et al. (2016), is built out of the standard Rolie–Poly model, see Likhtman and Graham (2003), and so forth. The presence of the stress diffusion term is motivated by several arguments.

First, the stress diffusion term naturally arises in the analysis of microscopic dumbbell models. If the presumably strong spatial inhomogeneity of the flow field is properly taken into account, see El-Kareh and Leal (1989), then a stress diffusion term appears in the averaged governing equations on the macroscopic level. (Here the term “diffusion” simply refers to the fact that the stress diffusion term is given by the Laplace operator.) Second, the addition of the stress diffusion term to the standard viscoelastic models has pleasing implications regarding the flow dynamics, in particular in the study of the shear banding phenomenon, see Fardin et al. (2015), Divoux et al. (2016) and Olmsted et al. (2000). Third, the presence of the stress diffusion term is convenient from a mathematical point of view since it provides a regularisation term in the equations, see for example Thomases (2011), Barrett and Boyaval (2011), Chupin and Martin (2015), Chupin et al. (2018) or Lukáčová-Medvid’ová et al. (2016).

This ad hoc approach might be acceptable provided that one is interested in the evolution of mechanical fields (pressure, velocity, extra stress) only, and if one is willing to completely ignore the temperature field. If the evolution of mechanical fields and thermal field is coupled, for example via the temperature dependence of the stress diffusion coefficient, which could be the case in practice, see Mohammadigoushki and Muller (2016), then a different approach to the model design must be taken. Namely thermodynamically consistent models describing the evolution of mechanical fields as well as the thermal field must be developed.

While a thermodynamical background is usually available for the standard viscoelastic models, see for example Wapperom and Hulsen (1998), Dressler et al. (1999), Rajagopal and Srinivasa (2000), Ellero et al. (2003) or Pavelka et al. (2018) for various approaches, much less is known about thermodynamically consistent models describing complex fluids such as wormlike micellar solutions. The existing thermodynamically consistent models are based on the modelling of internal structure of the fluid either by the means of a single internal parameter, see for example Manero et al. (2007), or by the means of multiple-species reactions in the spirit of the model introduced by Vasquez et al. (2007), see for example thermodynamical analysis by Grmela et al. (2010) or Germann et al. (2013).

On the other hand, thermodynamical analysis of models that contain the stress diffusion term in the form of Laplace operator has not been provided, to our best knowledge, until very recently, see Málek et al. (2018). In their work Málek et al. (2018) have followed the purely phenomenological thermodynamical approach to viscoelastic rate-type models introduced by Rajagopal and Srinivasa (2000), see also Málek et al. (2015a), Hron et al. (2017) or Málek and Průša (2017), and they have derived Maxwell/Oldroyd-B type models with the stress diffusion term both for compressible and incompressible fluids. In what follows we extend their analysis, and we derive diffusive variants of the Giesekus, FENE-P, Johnson–Segalman, Phan-Thien–Tanner and Bautista–Manero–Puig models.

The basic premise of the phenomenological thermodynamical approach by Rajagopal and Srinivasa (2000), see also Rajagopal and Srinivasa (2004), is that material of interest is characterised by the way it stores the energy and by the way it produces the entropy. Once these mechanisms are identified then the evolution equations simply follow from this identification and the underlying kinematical assumptions.

In their analysis Málek et al. (2018) have explored the possibility to interpret the stress diffusion effect in both ways. Namely, the stress diffusion term has been identified either as a symptom of a nonstandard energy storage mechanism or a nonstandard entropy production mechanism. In both cases one obtains a stress diffusion term in the evolution equation for the extra stress tensor, but the final models differ in the structure of the full Cauchy stress tensor and the specific heat capacity, see Málek et al. (2018) for details.

In the present analysis we focus on the second option, that is we interpret the stress diffusion as a consequence of a nonstandard entropy producing mechanism. This is consistent with the microscopic analysis by El-Kareh and Leal (1989), who have shown that the stress diffusion coefficient depends on the hydrodynamic resistance of one dumbbell bead. Consequently it is reasonable to interpret the stress diffusion as a dissipative (entropy producing) process. Concerning the characterisation of the entropy production, we in principle use the same ansatz as that used by Málek et al. (2018). In particular the corresponding term in the entropy production ansatz is quadratic in the gradient of the “extra stress” field.

Concerning the characterisation of the energy storage mechanism, which is in the given framework done via the specification of the Helmholtz free energy, we use formulae that are (Giesekus, Johnson–Segalman, Phan-Thien–Tanner models) structurally similar to that used by Málek et al. (2018), or that are known for the classical models without stress diffusion, see Wapperom and Hulsen (1998), Hu and Lelièvre (2007) and also Barrett and Boyaval (2017). However, contrary to the approach by Rajagopal and Srinivasa (2000) and Málek et al. (2018), we need to make different assumptions concerning the underlying kinematics. In particular, the kinematical assumptions must be designed in such a way that they allow one to replace the standard upper convected derivative by the Gordon–Schowalter derivative, see Gordon and Schowalter (1972). This requires substantial changes in the analysis provided by Málek et al. (2018).

Once the specification of the Helmholtz free energy (energy storage mechanism) and the entropy production (entropy production mechanism) is done, see Table 1, and once the underlying kinematical assumptions are known, then the governing equations for the mechanical variables as well as for the temperature follow immediately. This shows the flexibility of the adopted thermodynamical approach — the governing equations are indeed obtained by combining three simple pieces of information.

The paper is organised as follows. First, we briefly describe the thermodynamical framework introduced by Rajagopal and Srinivasa (2000), and we in detail describe the underlying kinematical assumptions that allow one to incorporate the Gordon–Schowalter derivative into the models. Once the preliminary work is done, we introduce the specific Helmholtz free energy and the entropy production that lead one after the other to the diffusive variants of the Giesekus, FENE-P, Johnson–Segalman, Phan-Thien–Tanner and Bautista–Manero–Puig models, and we explicitly write down the corresponding evolution equations for the temperature field. (The Helmholtz free energy/entropy production pairs for all the models are summarised in Table 1.) We also focus on the implied structure of the corresponding energy and entropy fluxes. This is essential piece of information provided that one wants to precisely characterize what is in the context of complex fluids meant by thermodynamically open/isolated systems. The paper is concluded by a brief discussion of the potential use of the provided thermodynamical background in the study of stability of flows of fluids described by the proposed models.

2 Preliminaries

The methodology introduced by Rajagopal and Srinivasa (2000) is purely phenomenological one in the sense that it does not rely on the knowledge of the internal microscopic structure of the given fluid. In particular, it does not operate with the notion of the conformation tensor. The additional tensorial variable arises as a consequence of the virtual decomposition of the total deformation to elastic and dissipative response, see Figure 2.1. This decomposition clearly mimics the standard spring–dashpot analogue for Maxwell type viscoelastic fluid, see for example Wineman and Rajagopal (2000).

If necessary, this decomposition can easily be generalised in order to design models that correspond to more involved spring-dashpot analogues such as Burgers model, see Karra and Rajagopal (2009a, b) or Málek et al. (2018), which is a model that is popular in the modelling of geomaterials and biological fluids, see for example Hron et al. (2014), Málek et al. (2015b); Málek et al. (2016) or Tůma et al. (2018). However, in order to get viscoelastic rate-type models with the Gordon–Schowalter derivative instead of the upper convected derivative, one needs to rethink the decomposition.

2.1 Kinematics

{ltxfigure}

[h] [Decomposition of deformation gradient to a dissipative and elastic part.]   [General decomposition of deformation gradient.] Viscoelastic fluid – kinematics.

2.1.1 Models with upper convected derivative

Let us now consider a general decomposition as shown in Figure 2.1. If the total deformation is seen as a composition of the two deformations, then the total deformation gradient $\mathbb{F}$ can be written as

[TABLE]

where $\mathbb{F}_{2}$ and $\mathbb{F}_{1}$ are the deformation gradients of the partial deformations. (Note that such a decomposition is frequent in other settings such as plasticity, see for example Gurtin et al. (2010).) Motivated by the standard relation

[TABLE]

between the spatial velocity gradient $\mathbb{L}=_{\mathrm{def}}\nabla\bm{v}$ and the deformation gradient $\mathbb{F}$, where

[TABLE]

denotes the material time derivative, we introduce new tensorial rate quantities $\mathbb{L}_{1}$ and $\mathbb{L}_{2}$, and also symmetric and skew-symmetric parts of these tensors. Their definitions along with the corresponding symmetric and asymmetric parts read

[TABLE]

Moreover, we define the right and left Cauchy–Green tensors associated to the partial deformation associated with the deformation gradient $\mathbb{F}_{2}$,

[TABLE]

By taking the material time derivative of Eq. 1, we obtain

[TABLE]

which then yields a kinematical identity

[TABLE]

where

[TABLE]

denotes the upper convected derivative. In the standard derivation by Rajagopal and Srinivasa (2000) one then uses thermodynamical arguments and writes the right hand side of Eq. 7 as a function of $\mathbb{B}_{\mathrm{2}}$.

2.1.2 Models with Gordon–Schowalter derivative

The evolution equations for both the Johnson–Segalman model, see Johnson and Segalman (1977), and the Phan–Thien–Tanner model, see Phan Thien and Tanner (1977), however contain the so-called Gordon–Schowalter derivative

[TABLE]

where $a\in[-1,1]$, see Gordon and Schowalter (1972). (Note that if $a=1$, then the Gordon–Schowalter derivative $\accentset{\medsquare}{{\mathbb{A}}}$ reduces to the standard upper-convected derivative $\accentset{\triangledown}{{\mathbb{A}}}$.) This type of derivative presumably takes into account “non-affine” molecular motion, see Johnson and Segalman (1977). Consequently, we would like to derive a kinematical formula similar to Eq. 7 for the Gordon–Schowalter derivative. The formula Eq. 9 can be rewritten as

[TABLE]

where we have defined

[TABLE]

Note that when $a=1$ the tensor $\mathbb{L}_{\mathrm{GS}}$ reduces to the standard velocity gradient $\mathbb{L}$. Motivated by Eq. 2 we then define tensor $\mathbb{F}_{\mathrm{GS}}$ as the solution of the differential equation

[TABLE]

Similarly to Eq. 11 we define

[TABLE]

and then define $\mathbb{F}_{2,\mathrm{GS}}$ as the solution of the differential equation

[TABLE]

As in the previous case we define the right and left Cauchy–Green tensors corresponding to the tensorial quantity $\mathbb{F}_{2,\mathrm{GS}}$

[TABLE]

We complete the definitions of the new tensorial quantities by defining $\mathbb{L}_{1,\mathrm{GS}}$ and $\mathbb{F}_{1,\mathrm{GS}}$ as111 The invertibility of the tensor $\mathbb{F}_{\mathrm{GS}}$ follows from the fact that

$\frac{\mathrm{d}{}}{\mathrm{d}{t}}\det\mathbb{F}_{\mathrm{GS}}=\left(\det\mathbb{F}_{\mathrm{GS}}\right)\operatorname{Tr}\left(\frac{\mathrm{d}{\mathbb{F}_{\mathrm{GS}}}}{\mathrm{d}{t}}\mathbb{F}_{\mathrm{GS}}^{-1}\right)=\left(\det\mathbb{F}_{\mathrm{GS}}\right)\operatorname{Tr}\mathbb{L}_{\mathrm{GS}}=a\left(\det\mathbb{F}_{\mathrm{GS}}\right)\operatorname{Tr}\mathbb{L}.$

(16)

We see that the solution of the differential equation Eq. 16, subject to the initial condition $\left.\det\mathbb{F}_{\mathrm{GS}}\right|_{t=0}=1$, reads $\det\mathbb{F}_{\mathrm{GS}}(t)=\exp\left(\int_{0}^{t}a\operatorname{Tr}\mathbb{L}(\tau)\,\mathrm{d}\tau\right)$. Hence, $\mathbb{F}_{\mathrm{GS}}$ remains invertible at all times. Similarly, we could show that the tensor $\mathbb{F}_{2,\mathrm{GS}}$ remains invertible at all times and thus Eq. 17a then yields that $\mathbb{F}_{1,\mathrm{GS}}$ is invertible as well.

[TABLE]

Definition Eq. 17a yields

[TABLE]

which is an analogue of the formula Eq. 1. By taking the material time of Eq. 18, we obtain

[TABLE]

which finally yields the kinematical identity

[TABLE]

where

[TABLE]

Formula Eq. 20 is an analogue of the kinematical identity Eq. 7. Further, one can easily show that

[TABLE]

where ${\mathbb{A}}:\mathbb{B}=_{\mathrm{def}}\operatorname{Tr}\left({\mathbb{A}}\mathbb{B}^{\top}\right)$ denotes the standard scalar product on the space of matrices. (Note that the first term on the right-hand-side of Eq. 22b vanishes for the isochoric deformation $\operatorname{div}\bm{v}=0$.) Moreover, using the kinematical identity Eq. 20, we can notice that

[TABLE]

which allows us to rewrite the right hand side of Eq. 22 completely in terms of $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ and $\mathbb{D}$.

2.2 Evolution equation for entropy

Once we specify the specific Helmholtz free energy $\psi$ (Helmholtz free energy per unit mass, $[\psi]=\nicefrac{\mathrm{J}}{\mathrm{kg}}$), we can use the general evolution equation for the specific internal energy $e$ of a continuous medium,

[TABLE]

see for example Gurtin et al. (2010), and we can derive the evolution equation for the specific entropy $\eta$. (Here $\rho$ denotes the density, $\mathbb{T}$ denotes the Cauchy stress tensor and $\bm{j}_{e}$ represents the non-mechanical contribution to the energy flux.) Indeed, using the chain rule, the standard set of thermodynamical identities such as $\psi=e-\theta\eta$, and the assumption that $\psi=\psi(\theta,\operatorname{Tr}\mathbb{B}_{\mathrm{2},\mathrm{GS}},\ln\det\mathbb{B}_{\mathrm{2},\mathrm{GS}},\varphi)$, one arrives at an expression for the time derivative of the internal energy in terms of the time derivative of the entropy $\eta$ and the time derivative of $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ and $\varphi$,

[TABLE]

see Hron et al. (2017) or Málek and Průša (2017) for details concerning a similar manipulation. In Eq. 25 the symbol $\theta$ denotes the temperature and $\varphi$ denotes a scalar quantity. (We will make use of $\varphi$ in the discussion of thermodynamics of Bautista–Manero–Puig type model.) Using Eq. 25 on the left-hand-side of Eq. 24 then yields the sought evolution equation for the entropy $\eta$.

In particular, for all the models studied in the present contribution we consider the following ansatz for the specific Helmholtz free energy

[TABLE]

For simplicity, the ansatz for the Helmholtz free energy splits the free energy into the thermal part $\psi_{0}$ and the mechanical part $\psi_{m}=\psi_{1}+\psi_{2}$. (See Hron et al. (2017) for more complicated formulae for the Helmholtz free energy.) The thermal part takes—in the simplest case of a fluid with a constant specific heat capacity—the form

[TABLE]

where $\theta_{\mathrm{ref}}$ is a reference temperature and $c_{\mathrm{V}}^{\mathrm{NSE}}$ is a positive material parameter (specific heat at constant volume). The postulated form of the Helmholtz free energy then yields an evolution equation for the specific entropy $\eta$ in the form

[TABLE]

where we have used kinematical identities Eq. 22 and Eq. 23, and where ${{\mathbb{A}}}_{\delta}=_{\mathrm{def}}{\mathbb{A}}-\frac{1}{3}\left(\operatorname{Tr}{\mathbb{A}}\right)\mathbb{I}$ denotes the traceless part of the corresponding tensor.

The derived entropy evolution equation Eq. 28 that is implied by the chosen ansatz for the Helmholtz free energy can then be “compared” with the general evolution equation for the entropy in a continuous medium that reads

[TABLE]

where $\xi$ denotes the required (known) entropy production and $\bm{j}_{\eta}$ denotes the entropy flux, see Málek and Průša (2017) for comments. This way one can identify the sought constitutive relations for $\mathbb{T}$ and $\bm{j}_{e}$ and the evolution equations for $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ and $\varphi$. Particular applications of this procedure in the derivation of specific viscoelastic rate-type models will be discussed in the following sections.

Concerning the stress diffusion term, the comparison of Eq. 28 and Eq. 29 will rely on using the identity

[TABLE]

where $c$ is a scalar quantity and ${\mathbb{A}}$ is a symmetric tensor, while the coordinate expression of the identity is

[TABLE]

where we have used the Einstein summation convention. Note that the triple product $\nabla{\mathbb{A}}\vdots\nabla{\mathbb{A}}$ on the left-hand-side is always a nonnegative quantity.

2.3 Evolution equation for temperature

Using the fact that $\eta=-\frac{\partial{\psi}}{\partial{\theta}}$, we can rewrite the evolution equation for the entropy Eq. 29 as

[TABLE]

Moreover appealing to the assumed simple structure of the Helmholtz free energy, namely to the splitting to the thermal and mechanical part, and the particular formula for the thermal part Eq. 27, we see that

[TABLE]

(Note that in the case of more complicated ansatz for the Helmholtz free energy one can get non-constant specific heat both with respect to the temperature and $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$, see for example Hron et al. (2017) or Málek et al. (2018).) Consequently, using the chain rule we see that Eq. 32 can be rewritten as an evolution equation for the temperature in the form

[TABLE]

This manipulation shows that once the Helmholtz free energy $\psi$ and entropy production $\xi$ are specified via a suitable ansatz, then the corresponding evolution equation for temperature is a simple consequence of the choice of $\psi$ and $\xi$.

3 Giesekus model

Let us now focus on specific viscoelastic rate-type models with a stress diffusion term. First, we identify the energy storage mechanisms and the entropy production mechanisms which yield the constitutive relations and in turn the evolution equations for a diffusive variant of the Giesekus model. Note that the system of evolution equations contains evolution equations for mechanical variables as well as for the temperature.

If we set $\tilde{\mu}=0$, that is if the stress diffusion term vanishes, see below, then we obtain, for mechanical variables, the original model proposed by Giesekus (1982). (No evolution equation for temperature was given in Giesekus (1982).) In this sense the analysis outlined below also provides a thermodynamical background for the original Giesekus model.

3.1 Helmholtz free energy

The ansatz for the specific Helmholtz free energy for the Giesekus model is chosen as

[TABLE]

where $\mu$ is a positive material constant and $\rho$ is a constant density. Notice that the Helmholtz free energy is, amongst others, specified in terms of the tensor $\mathbb{B}_{\mathrm{2}}$, which is a special case of the general tensorial quantity $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$. We thus set $a=1$, replace $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ with $\mathbb{B}_{\mathrm{2}}$ in Eq. 28, and we consider the upper convected derivative $\accentset{\triangledown}{{\mathbb{A}}}$ instead of the Gordon–Schowalter derivative $\accentset{\medsquare}{{\mathbb{A}}}$.

The postulated Helmholtz free energy yields the following formulae for the partial derivatives of the Helmholtz free energy that appear in Eq. 28

[TABLE]

By plugging Eq. 36 into Eq. 28 we arrive at

[TABLE]

This manipulation gives us the entropy evolution equation that is implied by the choice of the specific Helmholtz free energy.

3.2 Entropy production

Next, we postulate the entropy production for the Giesekus model with stress diffusion in the form

[TABLE]

where $\lambda\in[0,1]$ is a model parameter, $\nu$, $\nu_{1}$ are positive material constants, and $\tilde{\mu}$ is a material parameter that can be temperature dependent. The entropy production is chosen in such a way that it is a nonnegative function, which automatically guarantees the fulfillment of the second law of thermodynamics. The nonnegativity of the first, third and fourth term on the right-had side of Eq. 38 is obvious. Concerning the third term, we know that $\mathbb{B}_{\mathrm{2}}$ is a symmetric positive definite matrix, hence it is diagonalisable with the positive eigenvalues $\left\{\mu_{i}\right\}_{i=1}^{3}$ on the diagonal. Consequently, we see that

[TABLE]

where the function $\lambda\mu_{i}^{2}+(1-3\lambda)\mu_{i}+(1+\lambda)\frac{1}{\mu_{i}}+(3\lambda-2)$ is for $\lambda\in[0,1]$ and positive $\mu_{i}$ a nonnegative function that vanishes if and only if $\mu_{i}=1$. Finally, we can also observe that the entropy production vanishes in the equilibrium homogeneous rest state, that is if $\bm{v}=\bm{0}$, $\mathbb{B}_{\mathrm{2}}=\mathbb{I}$, and $\theta=\theta_{\mathrm{eq}}$, where $\theta_{\mathrm{eq}}$ is a spatially uniform temperature field.

The general evolution equation for entropy

[TABLE]

can then be transformed using the chosen entropy production $\xi$ and the identity Eq. 30, where we set $c=\tilde{\mu}$, ${\mathbb{A}}=\mathbb{B}_{\mathrm{2}}$, into the form

[TABLE]

3.3 Constitutive relations

Now we are in a position to read constitutive relations for $\mathbb{T}$, $\mathbb{B}_{\mathrm{2}}$, $\bm{j}_{e}$ and $\bm{j}_{\eta}$ out of Eq. 41 and Eq. 37. As discussed previously, all that needs to be done is to compare Eq. 41, which is the evolution equation for the entropy that is a consequence of the choice of the entropy production, with Eq. 37, which is the evolution equation for the entropy that is a consequence of the choice of the Helmholtz free energy. In particular, comparing the terms with the symmetric part of the velocity gradient, that is the term

[TABLE]

Making the comparison for the remaining terms in equations Eq. 37 and Eq. 41 then leads us to the constitutive relations

[TABLE]

Relations Eq. 43 are then used in the formulation of the evolution equations for the mechanical quantities as well as the temperature evolution equation. In particular, if we want to obtain the temperature evolution equation, we substitute for $\xi$ and $\bm{j}_{\eta}$ from Eq. 38 and Eq. 43d into the general temperature evolution equation Eq. 34.

3.4 Governing equations

Following the steps outlined above we see that the Helmholtz free energy ansatz Eq. 35 and the entropy production ansatz Eq. 38 lead to the following full system of governing equations for the diffusive Giesekus model

[TABLE]

and where $-m$ denotes the pressure. (See Málek and Průša (2017) and references therein concerning the notion of the pressure in incompressible fluids.) If body forces were present, the only change to the system would be the additional term $\rho\bm{b}$ in the balance of linear momentum Eq. 44b, that is Eq. 44b would read

[TABLE]

For the sake of simplicity, we however ignore this term in the remaining analysis.

One can note that if $\tilde{\mu}=0$ and $\lambda=0$, then the evolution equations for the mechanical variables, that is equations Eq. 44a, Eq. 44b and Eq. 44c, reduce to the governing equations of the standard Oldroyd-B fluid.

4 FENE-P model

Here, we identify the energy storage mechanisms and the entropy production mechanisms which yield the constitutive relations and in turn the evolution equations for a diffusive variant of the classical FENE-P, including the temperature evolution equation. The original model without the stress diffusion term was proposed by Bird et al. (1980), see also comments in Keunings (1997). The procedure for the derivation of the model is the same as in the case of the Giesekus model, hence we will proceed with less detail.

4.1 Helmholtz free energy

The ansatz for the specific Helmholtz free energy for the FENE-P model is chosen as

[TABLE]

where $\mu$ is a positive material constant and $b$ is a positive model constant which represents an upper bound on the trace of the tensor $\mathbb{B}_{\mathrm{2}}$, that is $\operatorname{Tr}\mathbb{B}_{\mathrm{2}}<b$. Notice that the Helmholtz free energy is, amongst others, specified in terms of the tensor $\mathbb{B}_{\mathrm{2}}$, the special case of the general tensorial quantity $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$. We thus set $a=1$, replace $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ with $\mathbb{B}_{\mathrm{2}}$ in Eq. 28, and we again consider the upper convected derivative instead of the Gordon–Schowalter derivative.

The postulated Helmholtz free energy yields the following formulae for the partial derivatives of the Helmholtz free energy that appear in Eq. 28

[TABLE]

By plugging Eq. 47 into Eq. 28 we get

[TABLE]

4.2 Entropy production

Now we postulate the entropy production for the FENE-P model with stress diffusion in the form

[TABLE]

where $\nu$, $\nu_{1}$ are positive material constants, and $\tilde{\mu}$ is a (possibly) temperature-dependent material parameter. Since $\mathbb{B}_{\mathrm{2}}$ is a symmetric positive definite matrix we can easily show that the second term of Eq. 49 is nonnegative by using the spectral decomposition of $\mathbb{B}_{\mathrm{2}}$, and by investigating the nonnegativity of the corresponding scalar function of eigenvalues of $\mathbb{B}_{\mathrm{2}}$. (See Eq. 39 for a detailed discussion concerning the same issue.) Consequently, the proposed entropy production is a nonnegative function.

The general evolution equation for entropy

[TABLE]

can then be transformed using the chosen entropy production $\xi$ and the identity Eq. 30, where we set $c=\tilde{\mu}$, ${\mathbb{A}}=\left(1-\frac{1}{b}\operatorname{Tr}\mathbb{B}_{\mathrm{2}}\right)^{-1}\mathbb{B}_{\mathrm{2}}$, into the form

[TABLE]

Equation Eq. 51 can now be compared with the derived entropy evolution equation Eq. 48.

4.3 Constitutive relations

In particular, equations Eq. 48 and Eq. 51 will coincide if we set

[TABLE]

Relations Eq. 52 can then be used for formulating the evolution equations for the mechanical quantities as well as the temperature evolution equation.

4.4 Governing equations

Following the steps outlined above we see that the Helmholtz free energy ansatz Eq. 46 and the entropy production ansatz Eq. 49 lead to the following full system of governing equations for the diffusive FENE-P model

[TABLE]

We again see that if $\tilde{\mu}=0$, that is if the stress diffusion term vanishes, then the governing equations for the mechanical variables coincide with the governing equations for the original FENE-P model without stress diffusion.

5 Johnson–Segalman model

Here, we identify the energy storage mechanisms and the entropy production mechanisms which yield the constitutive relations and in turn the evolution equations for a diffusive variant of the Johnson–Segalman model, including the temperature evolution equation. The original model without the stress diffusion term was proposed by Johnson and Segalman (1977). Unlike in the previous sections devoted to the analysis of the Giesekus and FENE-P models, we need to change the kinematical assumptions, since we want to obtain a model with the Gordon–Schowalter derivative instead of the upper convected derivative. The kinematical assumptions that lead to the Gordon–Schowalter derivative have been already discussed in the introductory sections, hence, regarding the discussion of the kinematics, we restrict ourselves to a simple statement that we are going to replace tensor $\mathbb{B}_{\mathrm{2}}$ by $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$.

5.1 Helmholtz free energy

The ansatz for the Helmholtz free energy of the Johnson–Segalman model is chosen as

[TABLE]

where $\mu$ is a positive material constant. Note that the chosen formula for the Helmholtz free energy is structurally identical to that for the Giesekus model, the only difference between Eq. 54 and Eq. 35 is that $\mathbb{B}_{\mathrm{2}}$ has been replaced by $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$. The postulated Helmholtz free energy yields the following formulae for the partial derivatives that appear in Eq. 28

[TABLE]

By plugging Eq. 55 into Eq. 28 we get

[TABLE]

5.2 Entropy production

We postulate the entropy production for the Johnson–Segalman model with stress diffusion in the form

[TABLE]

where $\nu$, $\nu_{1}$ are positive material constants, and $\tilde{\mu}$ is a (possibly) temperature-dependent material parameter. The chosen entropy production is a nonnegative function. The general evolution equation for entropy

[TABLE]

can then be transformed using the chosen entropy production $\xi$ and the identity Eq. 30, where we set $c=\tilde{\mu}$, ${\mathbb{A}}=\mathbb{B}_{\mathrm{2},\mathrm{GS}}$, into the form

[TABLE]

Equation Eq. 59 can now be compared with the derived entropy evolution equation Eq. 56 to yield constitutive relations.

5.3 Constitutive relations

In particular, equations Eq. 56 and Eq. 59 will coincide if we set

[TABLE]

Relations Eq. 60 can then be used for formulating the evolution equations for the mechanical quantities as well as the temperature evolution equation.

5.4 Governing equations

The full system of governing equations for the Johnson-Segalman model with stress diffusion reads

[TABLE]

We note that the governing equations are nearly the same as the governing equations for the diffusive Oldroyd-B model. (Recall that Oldroyd-B model is the Giesekus model with $\lambda=0$, see Eq. 44.) This is not surprising given the fact that the formulae for the Helmholtz free energy and the entropy production are structurally the same. The difference between Oldroyd-B model and the Johnson–Segalman models is in the choice of the kinematical quantity $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ versus $\mathbb{B}_{\mathrm{2}}$. This different choice is reflected in the final governing equations via the presence of the Gordon–Schowalter/upper convected derivative in the evolution equations, and by the presence of the factor $a$ in the formula for the Cauchy stress tensor, see Eq. 61e and Eq. 44e.

6 Phan–Thien–Tanner model

Here, we identify the energy storage mechanisms and the entropy production mechanisms which yield the constitutive relations and in turn the evolution equations for a diffusive variant of the Phan–Thien–Tanner model, including the temperature evolution equation. The Phan–Thien–Tanner model without a stress diffusion term has several variants. The original model was proposed by Phan Thien and Tanner (1977) and it corresponds to case Eq. 66a for the choice of the function appearing in the evolution equation for the extra stress tensor, see below. The exponential variant Eq. 66c was proposed by Phan Thien (1978). The quadratic variant Eq. 66b is far less seen in the literature, but is used for example by Ngamaramvaranggul and Webster (2002).

6.1 Helmholtz free energy

The ansatz for the Helmholtz free energy of the Phan–Thien–Tanner model is chosen as

[TABLE]

where $\mu$ is a positive material constant. The postulated Helmholtz free energy yields the following formulae for the partial derivatives that appear in Eq. 28

[TABLE]

By plugging Eq. 63 into Eq. 28 we get

[TABLE]

6.2 Entropy production

We postulate the entropy production for the Phan–Thien–Tanner model with stress diffusion in the form

[TABLE]

where $\nu$, $\nu_{1}$ are the positive material constants, $\tilde{\mu}$ is a (possibly) temperature-dependent material parameter, and the function $f$ is defined by one of the formulae

[TABLE]

where $p$ is a model parameter. Since the function $f$ appears in the entropy production Eq. 65, we see that the nonnegativity of the entropy production will be guaranteed provided that we will place some restrictions on $f$.

In particular, we want $f$ to be nonnegative for any symmetric positive definite matrix. The nonnegativity of $f$ is guaranteed for any $p$ provided that $f$ is given by Eq. 66c. On the other hand, if $f$ is given by Eq. 66a, then we need to restrict the values of $p$ to the interval $\left[0,\frac{1}{3}\right]$. This requirement follows from the spectral decomposition. If ${\mathbb{A}}$ is a symmetric positive definite matrix, then

[TABLE]

where $\left\{\lambda_{i}\right\}_{i=1}^{3}$ are positive eigenvalues of ${\mathbb{A}}$. Consequently $1+p\operatorname{Tr}\left({\mathbb{A}}-\mathbb{I}\right)>0$ provided that $p\in\left[0,\frac{1}{3}\right]$. Similarly, if $f$ is given by Eq. 66b, then we again need to restrict $p$ to the interval $\left[0,\frac{1}{3}\right]$.

Having the ansatz for the entropy production $\xi$, we see that the general evolution equation for entropy

[TABLE]

can then be transformed using the chosen entropy production $\xi$ and the identity Eq. 30, where we set $c=\tilde{\mu}$, ${\mathbb{A}}=\mathbb{B}_{\mathrm{2},\mathrm{GS}}$, into the following form

[TABLE]

Equation Eq. 68 can now be compared with the derived entropy evolution equation Eq. 64 to yield constitutive relations.

6.3 Constitutive relations

In particular, equations Eq. 64 and Eq. 68 will coincide if we set

[TABLE]

Relations Eq. 69 can then be used for formulating the evolution equations for the mechanical quantities as well as the temperature evolution equation.

6.4 Governing equations

The full system of governing equations for the diffusive variant of the Phan-Thien–Tanner model reads

[TABLE]

and the function $f$ is given by Eq. 66a, Eq. 66b or Eq. 66c. Adjusting the value of the parametr $a$ we can again select the convenient frame-indifferent derivative in the evolution equation for the extra stress tensor. Note also that if we have set $f=1$, then we would obtain the Johnson–Segalman model.

7 Bautista–Manero–Puig type model

Here, we identify the energy storage mechanisms and the entropy production mechanisms which yield the constitutive relations and in turn the evolution equations for a diffusive variant of Bautista–Manero–Puig model, including the temperature evolution equation.

The Bautista–Manero-Puig type models are based on the idea that the fluid can undergo microstructural changes that are triggered by the flow. In this respect, this class of models is a generalisation of the original idea by Fredrickson (1970). Fredrickson (1970) has considered the standard Navier–Stokes fluid, but he has assumed that the fluidity (the reciprocal value of viscosity) is given by an extra evolution equation. Bautista–Manero–Puig type models basically extend this idea to viscoelastic rate-type models, see Bautista et al. (1999) and the follow-up studies by Bautista et al. (2000, 2002); Manero et al. (2007); López-Aguilar et al. (2016).

The model presented below is based on the same idea that leads to the Bautista–Manero–Puig model. However, even if we consider the model without the stress diffusion term, then the present model is in certain aspects different from the original Bautista–Manero–Puig model. We will discuss this issue when we come to the point where the difference arises.

7.1 Helmholtz free energy

The ansatz for the Helmholtz free energy for the diffusive Bautista–Manero–Puig type model is chosen as

[TABLE]

where $\mu$, $\chi$, and $\varphi_{0}$ are positive material parameters. Notice that the Helmholtz free energy is a function of the tensor $\mathbb{B}_{\mathrm{2}}$, which is a special instance of the general tensorial quantity $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$. We thus set $a=1$, and we replace $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ with $\mathbb{B}_{\mathrm{2}}$ in Eq. 28. This will lead us to the use of the upper convected derivative instead of the Gordon–Schowalter derivative in the corresponding governing equation. (See also the derivation of the Giesekus model for a similar discussion.) The major difference between the Helmholtz free energy ansatz used in the previous sections and the Helmholtz free energy ansatz Eq. 71 used in the current section is the presence of a new positive scalar variable $\varphi$. This variable serves as a proxy to account for microstructural changes.

The postulated Helmholtz free energy yields the following formulae for the partial derivatives that appear in Eq. 28

[TABLE]

By plugging Eq. 72 into Eq. 28 we get

[TABLE]

7.2 Entropy production

The entropy production for the Bautista–Manero–Puig type model with stress diffusion is postulated in the form

[TABLE]

where $\nu$ and $\beta$ are positive material constants, and $\tilde{\mu}$ is a (possibly) temperature-dependent material parameter. We note that the chosen entropy production is nonnegative provided that we can guarantee that $\varphi$ is a positive quantity, which is the issue we will discuss later. The general evolution equation for entropy

[TABLE]

can then be transformed using the chosen entropy production $\xi$ and the identity Eq. 30, where we set $c=\tilde{\mu}$, ${\mathbb{A}}=\mathbb{B}_{\mathrm{2}}$, into the form

[TABLE]

Equation Eq. 76 can now be compared with the derived entropy evolution equation Eq. 73 to yield constitutive relations which guarantee the validity of the second law of thermodynamics.

7.3 Constitutive relations

In particular, equations Eq. 73 and Eq. 76 will coincide if we set

[TABLE]

where $\alpha$ and $\varphi_{\infty}$ are positive material parameters. Relations Eq. 77 can then be used for formulating the evolution equations for the mechanical quantities as well as the temperature evolution equation by appealing the standard procedure. Before we formulate the complete set of governing equations, we however need to discuss the positivity of $\varphi$, which is needed for nonnegativity of the entropy production ansatz Eq. 74.

Let us assume that $\varphi_{\infty}>\varphi_{0}>0$, and let us analyse the behaviour of $\varphi$ as it is implied by Eq. 77c. We see that if $\varphi<\varphi_{0}$, then the right-hand side of Eq. 77c is positive. It means that whenever we have $\varphi<\varphi_{0}$, then $\varphi$ is an increasing function, and, in particular the value of $\varphi$ can not decrease to zero. On the other hand if $\varphi>\varphi_{\infty}$, then the right-hand side of Eq. 77c is negative, hence $\varphi$ is a decreasing function. Consequently, if $\varphi\in\left[\varphi_{0},\varphi_{\infty}\right)$, then Eq. 77c predicts that $\varphi$ will stay in this interval forever. Moreover if ${\mathbb{D}}_{\delta}:{\mathbb{D}}_{\delta}$ is small then $\varphi$ has a tendency to go to the equilibrium value $\varphi_{0}$, and, on the other hand, if ${\mathbb{D}}_{\delta}:{\mathbb{D}}_{\delta}$ is large, then $\varphi$ is driven to $\varphi_{\infty}$. In this sense $\varphi$ indeed works as a variable monitoring the “flow-induced microstructural changes” from the equilibrium value $\varphi_{0}$ to the nonequilibrium (strong shear flow) value $\varphi_{\infty}$. A similar analysis can be done in the case $\varphi_{0}>\varphi_{\infty}>0$.

We note that Eq. 77c is the equation that is different from the equation used in the Bautista–Manero–Puig models, since in these models the evolution equation for $\varphi$ usually takes the form

[TABLE]

In the simple setting outlined above we however cannot guarantee the positivity of the product ${\left(\mathbb{B}_{\mathrm{2}}\right)}_{\delta}:{\mathbb{D}}_{\delta}$, and consequently the positivity of $\varphi$, hence we are forced to replace Eq. 78 by Eq. 77c.

7.4 Governing equations

The full system of governing equations for the diffusive Bautista–Manero–Puig type model reads

[TABLE]

8 Implications of the thermodynamical background

Having established the thermodynamical background for various viscoelastic rate-type models with stress diffusion, we can focus on its implications for the analysis of flows of fluids described by these models.

8.1 Boundary conditions – characterisation of thermodynamically isolated or open systems

First, as a by-product of the thermodynamical analysis, we have obtained explicit expressions for the energy flux and entropy flux. This means that we can say which choice of boundary conditions leads to thermodynamically isolated/open systems. Let us consider the fluid occupying the domain $\Omega$. In order to get a thermodynamically isolated system, we have to impose the no-penetration boundary condition on the velocity

[TABLE]

where $\bm{j}_{e}$ is given by a model specific formula. (See for example Eq. 43c for the non-mechanical energy flux in the diffusive Giesekus model.)

8.2 Evolution equation for temperature – complete system of governing equations for complex processes

Second, for all the discussed models we have an evolution equation for the temperature that properly takes into account all energy transfer mechanisms in the fluid described by the given model. This allows one to investigate not only the evolution of mechanical variables but also of the temperature field.

8.3 Qualitative behaviour – conserved quantities and stability analysis

Third, the outlined analysis also allows one to immediately identify the quantities that are conserved during the motion. Let us for example consider the system of governing equations for the diffusive Johnson–Segalman model in the absence of body forces, that is $\bm{b}=\bm{0}$. The system reads

[TABLE]

see Eq. 61, and where $\accentset{\medsquare}{\overline{\mathbb{B}_{\mathrm{2},\mathrm{GS}}}}$ denotes the Gordon–Schowalter derivative,

[TABLE]

Given the system of governing equations one can ask the question whether there is a quantity that is conserved in the motion.

We know that in isolated systems the quantity that is being conserved is the net total energy. (This means the total energy $\rho e+\frac{1}{2}\rho\left|\bm{v}\right|^{2}$ integrated over the domain $\Omega$.) Having studied the thermodynamical background of the models, we are able to find an explicit formula for the total energy. We know that the Helmholtz free energy is in the case of the diffusive Johnson–Segalman model given by the formula

[TABLE]

see Eq. 54 and Eq. 27. The internal energy is obtained from the Helmholtz free energy via the formula $e=\psi+\theta\eta$, where $\eta=-\frac{\partial{\psi}}{\partial{\theta}}$, which yields

[TABLE]

Consequently, for the diffusive Johnson–Segalman model, we expect that in the isolated system we get

[TABLE]

Let us see whether this piece of information can be directly derived using the evolution equations only.

We start with a standard manipulation, see for example Gurtin et al. (2010). We take the scalar product of the balance of liner momentum Eq. 81b with the velocity field $\bm{v}$, we integrate the equation over the domain $\Omega$, and we use the integration by parts formula, and the no-penetration boundary condition Eq. 80a, which yields

[TABLE]

If we take into account the incompressibility condition Eq. 81a and the structure of the Cauchy stress tensor, see Eq. 81e, and if we use the Stokes theorem, then we can rewrite Eq. 86 as

[TABLE]

where $\bm{n}$ denotes the unit outward normal to the boundary of $\Omega$.

Now take the trace of the evolution equation Eq. 81c for $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ which, upon using the definition of the Gordon–Schowalter derivative, yields

[TABLE]

Next we multiply the evolution equation Eq. 81c for $\mathbb{B}_{\mathrm{2},\mathrm{GS}}$ by $\mathbb{B}_{\mathrm{2},\mathrm{GS}}^{-1}$ and then we take the trace. This manipulation yields

[TABLE]

where we have used the incompressibility constraint $\operatorname{div}\bm{v}=\operatorname{Tr}\mathbb{D}=0$, and the well known formula $\frac{\mathrm{d}{}}{\mathrm{d}{t}}\ln\det{\mathbb{A}}=\operatorname{Tr}\left(\frac{\mathrm{d}{{\mathbb{A}}}}{\mathrm{d}{t}}{\mathbb{A}}^{-1}\right)$. Subtracting Eq. 89 from Eq. 88 and integrating over the domain $\Omega$ finally leads to the equation

[TABLE]

Having derived Eq. 90, we can multiply Eq. 90 by $\frac{\mu}{2\nu_{1}}$ and add it to Eq. 86, which yields

[TABLE]

(On the left hand side we have added constant $3$ to the term under the time derivative, which is an operation that does not change the resulting equation. The constant is added for the matter of convenience, it basically fixes the energy of the equilibrium state $\mathbb{B}_{\mathrm{2},\mathrm{GS}}=\mathbb{I}$ to zero.) Using Eq. 30 we rewrite the terms on the right hand side, and we get

[TABLE]

The next to the last integral can be converted to a surface integral using the Stokes theorem. Moreover, we see that the term under the divergence can be rewritten as $-\bm{j}_{e}-\kappa\nabla\theta$, see Eq. 60c, which yields

[TABLE]

If we are dealing with an isolated system (no work is done on the system, no heat flux into the system, no energy flux into the system), then all the surface integrals vanish a we see that the quantity

[TABLE]

decays in time. (All the remaining terms on the right hand side are negative.) This is not surprising provided the we know the thermodynamical background of the model.

In fact we could have seen it a priori without any manipulation of the governing equations. The fact that the quantity given by Eq. 94 decays in time is a consequence of the design of the model. In a thermodynamically isolated system the net mechanical energy, which is precisely the quantity given by the integral Eq. 94, must degrade to the thermal energy. Note that this simple consequence of the thermodynamical background of the model was rather difficult to obtain directly from the governing equations.

Let us further investigate Eq. 93. Integrating the evolution equation for the temperature, see Eq. 81d, and adding the result to Eq. 93 we get

[TABLE]

where we have again used the Stokes theorem. The equality tells us that the quantity on the left-hand side is conserved in thermodynamically isolated systems, that is in systems with boundary conditions specified via Eq. 80. If Eq. 80 holds, then Eq. 95 coincides with Eq. 85 which automatically follows from the thermodynamical background of the model. Manipulations of this type are essential in the mathematical analysis of the corresponding models, see for example Bulíček et al. (2018, 2018).

Besides the automatic identification of conserved quantities in thermodynamically isolated systems, the knowledge of thermodynamical background of the models can be also exploited in the stability analysis. For example, if we consider a thermodynamically isolated system such as a fluid in a thermodynamically isolated vessel, then we expect that the fluid will eventually reach a homogeneous steady state regardless of the initial state of the fluid. (Here homogeneous steady state means zero velocity field, $\mathbb{B}_{\mathrm{2},\mathrm{GS}}=\mathbb{I}$ and spatially homogeneous temperature field.) This natural tendency would be rather difficult to see directly from the governing equations.

However, since we know the thermodynamical background of the models, we know that “entropy rises” and “energy is conserved”. Consequently, we can design a functional that can serve as a Lyapunov-like functional that charaterises the approach to the homogeneous steady state, see Coleman (1970) and Gurtin (1973, 1975) for early studies in this direction. The Lyapunov-like functional is made up of the net entropy $S$ and the net total energy $E_{\mathrm{tot}}$

[TABLE]

and it is given by the formula

[TABLE]

where $\theta_{\mathrm{ref}}$ is the temperature at the steady state and $\widehat{E_{\mathrm{tot}}}$ is the energy at the steady state (a known constant), see Bulíček et al. (2017) for details. Fortunately, since we have completely explored the thermodynamical background of the models, we know the explicit formulae for the internal energy $e$ and the specific entropy $\eta$ for the given models, hence we can explicitly construct the functional.

Taking the time derivative of the functional $\mathcal{V}_{\theta_{\mathrm{ref}},\mathrm{eq}}$ yields

[TABLE]

where we have exploited the fact that the net total energy is conserved in a thermodynamically isolated system, $\frac{\mathrm{d}{E_{\mathrm{tot}}}}{\mathrm{d}{t}}=0$, and that the net entropy grows in a thermodynamically isolated system. (The growth of the net entropy follows by integration of the evolution equation Eq. 29 the domain $\Omega$, the positivity of the entropy production term $\xi$ and the no-flux boundary condition for the entropy flux.) This simple manipulation shows that the functional indeed decreases in time hence it is truly a Lyapunov-like functional.

The reader interested in the application of this procedure in the context of viscoelastic rate-type fluids is kindly referred to Málek et al. (2018). In this work the authors have analysed the stability of the equilibrium steady state of a fluid described by a diffusive variant of Oldroyd-B model. The successful identification of the energy storage and entropy producing mechanisms for other diffusive viscoelastic-type models, see above, opens up the possibility to redo the stability analysis by Málek et al. (2018) for these models as well.

Thermodynamical background can be however exploited in more complex settings than the stability analysis of the rest state. In fact, thermodynamical background is of use even in the stability analysis of steady states in thermodynamically open systems, see Bulíček et al. (2017). Bulíček et al. (2017) have conjectured, that the functional

[TABLE]

where

[TABLE]

could serve as a Lyapunov-like functional for the analysis of the steady state in a thermodynamically open system. Here $\widehat{\bm{W}}=\left[\widehat{\bm{v}},\widehat{\mathbb{B}_{\mathrm{2},\mathrm{GS}}},\widehat{\theta}\right]$ denotes the vector of state variables in the steady state, and $\widetilde{\bm{W}}=\left[\widetilde{\bm{v}},\widetilde{\mathbb{B}_{\mathrm{2},\mathrm{GS}}},\widetilde{\theta}\right]$ denotes the vector of the perturbations of the state variables with respect to the steady state, and $\left.\mathrm{D}_{\bm{W}}S_{\widehat{\theta}}\left(\bm{W}\right)\right|_{\bm{W}=\widehat{\bm{W}}}\left[\widetilde{\bm{W}}\right]$ and $\left.\mathrm{D}_{\bm{W}}E_{\mathrm{tot}}\left(\bm{W}\right)\right|_{\bm{W}=\widehat{\bm{W}}}\left[\widetilde{\bm{W}}\right]$ denote the Gâteaux derivative222 Let us recall that the Gâteaux derivative $\mathrm{D}\mathcal{M}(\bm{x})[\bm{y}]$ of a functional $\mathcal{M}$ at point $\bm{x}$ in the direction $\bm{y}$ is defined as $\mathrm{D}\mathcal{M}(\bm{x})[\bm{y}]=_{\mathrm{def}}\lim_{s\to 0}\frac{\mathcal{M}(\bm{x}+s\bm{y})-\mathcal{M}(\bm{x})}{s}$ which is tantamount to $\mathrm{D}\mathcal{M}(\bm{x})[\bm{y}]=_{\mathrm{def}}\left.\frac{\mathrm{d}{}}{\mathrm{d}{s}}\mathcal{M}(\bm{x}+s\bm{y})\right|_{s=0}$. If it is necessary to emphasize the variable against which we differentiate, we also write $\mathrm{D}_{\bm{x}}\mathcal{M}(\bm{x})[\bm{y}]$ instead of $\mathrm{D}\mathcal{M}(\bm{x})[\bm{y}]$. of the given functional with respect to $\bm{W}$ at point $\widehat{\bm{W}}$ in the direction $\widetilde{\bm{W}}$, see Bulíček et al. (2017) for details.

As in the previous case the critical piece of information necessary to build the Lyapunov-like functionals is the knowledge of the Helmholtz free energy, which allows one to obtain the explicit formulae for the entropy and the internal energy, and hence the explicit formulae for the functionals. A reader interested in the application of this procedure in the context of viscoelastic rate-type fluids is kindly referred to Dostalík et al. (2018), who have investigated the stability of a steady flow of a fluid described by the Giesekus model. The successful identification of the energy storage and entropy producing mechanisms for the diffusive viscoelastic-type models, see above, opens up the possibility to redo the stability analysis by Dostalík et al. (2018) for these models as well.

9 Conclusion

We have presented a purely phenomenological approach to the viscoelastic rate-type models with a stress diffusion term. In particular we have identified energy storage mechanisms and entropy production mechanisms that lead to the diffusive variants of the Giesekus, FENE-P, Johnson–Segalman, Phan-Thien–Tanner and Bautista–Manero–Puig models. In all these cases the stress diffusion term has been interpreted as a consequence of a nonstandard entropy production mechanism. Namely we have assumed that the entropy production contains a quadratic term in the gradient of the “extra stress tensor”, which in turn have lead to the presence of the stress diffusion term in the corresponding governing equation.

All the models have been based on the knowledge of the triad energy storage mechanisms/entropy production mechanisms/underlying kinematics. The governing equations—including the evolution equation for the temperature field—have been interpreted as a simple consequence of the specifications of the individual components of the triad. The nature of the outlined procedure allows one to freely choose individual components of the triad, and alternatively derive more complicated models. For example, the simple quadratic term in the gradient of the “extra stress tensor” which appears in the entropy production, can be replaced, if necessary, by a more complicated formula. The same holds for the thermal part of the Helmholtz free energy/entropy production that has been for simplicity chosen in the form that implies constant specific heat capacity and the standard Fourier’s law for the heat flux.

Finally, we have discussed the potential use of the provided thermodynamical background in the study of stability of flows of fluids described by the proposed models.

Mark Dostalík and Vít Průša acknowledge the support of the Czech Science Foundation project 18-12719S. Mark Dostalík has also been supported by Charles University Research programme No. UNCE/SCI/023.

Appendix A Corrections

A.1 Helmholtz free energy for the FENE-P model

The Helmholtz free energy ansatz for the FENE-P model should read

[TABLE]

instead of Eq. 46 which reads

[TABLE]

The change of the additive constant has no impact on the constitutive relations and the governing equations, since these depend on the derivatives of the Helmholtz free energy.

On the other hand, it is convenient to choose the Helmholtz free energy such that it vanishes at the equilibrium steady state. If the velocity field $\bm{v}$ vanishes, then Eq. 53c implies that the equilibrium steady state value of $\mathbb{B}_{\mathrm{2}}$ is

[TABLE]

Consequently if we want the (mechanical) part of the Helmholtz free energy to vanish for ${\mathbb{B}_{\mathrm{2}}}_{\mathrm{eq}}$, we have to replace Eq. 102 with Eq. 101. The fact that the mechanical part of Eq. 101 vanishes for ${\mathbb{B}_{\mathrm{2}}}_{\mathrm{eq}}$ given by Eq. 103 is straightforward to show by direct substitution. Moreover, if we take the formal limit $b\to+\infty$, we see that the Helmholtz free energy ansatz for the FENE-P model reduces to the standard ansatz for the Oldroyd-B model.

A.2 Entropy production for the Giesekus model

There is a misprint in the entropy production formula Eq. 38 for the Giesekus model. The wrong formula reads

[TABLE]

while the correct formula is

[TABLE]

The misprint is in the term $(1+\lambda)\mathbb{B}_{\mathrm{2}}^{-1}$ which should read $(1-\lambda)\mathbb{B}_{\mathrm{2}}^{-1}$. The same misprint appears in Eq. 39 and Eq. 44d and in Table 1.